Add 2D axisymmetric WENO5 convergence test and fix three bugs

sbryngelson · sbryngelson · commit 743016a35b30 · 2026-05-08T12:47:37.000-04:00
- Add examples/2D_axisym_convergence/case.py: 2D cylindrical density
  sine wave convergence case with HLLC (riemann_solver=2) and periodic
  x-BCs so the exact solution is valid at every cell.

- Add toolchain/mfc/test/run_convergence_axisym.py: convergence runner
  that refines Nz=[64,128,256] with Nr=32 fixed, reads alpha_rho from
  p_all binary, compares r-averaged density to the analytic solution,
  and asserts fitted WENO5 rate &gt;= 4.8 (observed: 5.08).

- Fix HLL and LF Riemann solvers (m_riemann_solvers.fpp): void-fraction
  geometric source was incorrectly set to flux_rs instead of 0.

- Fix single-fluid cyl_coord alpha drift (m_time_steppers.fpp): after
  each RK3 stage, reset alpha=1 for num_fluids=1 + cyl_coord cases to
  prevent pressure NaN from contact-wave-speed variation across faces.

- Add convergence-2d-axisym CI job to .github/workflows/convergence.yml.
diff --git a/.github/workflows/convergence.yml b/.github/workflows/convergence.yml
@@ -97,6 +97,34 @@ jobs:
             --resolutions 64 128 256 512 \
             --num-ranks 4
 
+  convergence-2d-axisym:
+    name: "2D Axisymmetric WENO5 Convergence"
+    runs-on: ubuntu-latest
+    timeout-minutes: 60
+
+    steps:
+      - uses: actions/checkout@v4
+
+      - name: Setup Ubuntu
+        run: |
+          sudo apt update -y
+          sudo apt install -y cmake gcc g++ python3 python3-dev \
+                   openmpi-bin libopenmpi-dev
+
+      - name: Setup Python
+        uses: actions/setup-python@v5
+        with:
+          python-version: "3.12"
+
+      - name: Initialize MFC
+        run: ./mfc.sh init
+
+      - name: Run 2D axisymmetric convergence test
+        run: |
+          source build/venv/bin/activate
+          python toolchain/mfc/test/run_convergence_axisym.py \
+            --resolutions 64 128 256
+
   convergence-temporal:
     name: "RK3 Temporal Order"
     runs-on: ubuntu-latest
diff --git a/examples/2D_axisym_convergence/case.py b/examples/2D_axisym_convergence/case.py
@@ -0,0 +1,85 @@
+#!/usr/bin/env python3
+"""
+2D axisymmetric convergence case: density sine wave in x (axial/z), cyl_coord=T.
+
+In MFC 2D cylindrical: x=axial(z), y=radial(r).
+rho = 1 + 0.2*sin(2*pi*x), u_x=1, p=1, u_y=0.
+Exact solution at time T: rho(x,T) = 1 + 0.2*sin(2*pi*(x-T)).
+Domain [0,5] = 5 full periods of the sine IC; periodic x-BCs make the exact
+solution valid everywhere (no upstream contamination from ghost cells).
+"""
+
+import argparse
+import json
+import math
+import sys
+
+parser = argparse.ArgumentParser()
+parser.add_argument("--mfc", type=str, default=None)
+parser.add_argument("-N", type=int, default=64, help="Grid cells in x (axial) direction")
+parser.add_argument("--order", type=int, default=5)
+parser.add_argument("--cfl", type=float, default=0.02)
+parser.add_argument("--nr", type=int, default=6, help="Radial cells (fixed)")
+parser.add_argument("--Tfinal", type=float, default=0.1, help="Final simulation time")
+args = parser.parse_args()
+
+N = args.N  # axial cells (x-direction, refined)
+Nr = args.nr  # radial cells (y-direction, fixed)
+# Domain [0,5] = 5 periods of sin(2*pi*x); periodic x-BCs make it exact everywhere
+Lx = 5.0
+Lr = 0.5  # radial domain length
+
+dx = Lx / N
+dt = args.cfl * dx
+Nt = math.ceil(args.Tfinal / dt)
+dt = args.Tfinal / Nt  # adjust to hit Tfinal exactly
+
+case = {
+    # x=axial, y=radial
+    "m": N - 1,
+    "n": Nr - 1,
+    "p": 0,
+    "x_domain%beg": 0.0,
+    "x_domain%end": Lx,
+    "y_domain%beg": 0.0,
+    "y_domain%end": Lr,
+    "cyl_coord": "T",
+    "dt": dt,
+    "t_step_start": 0,
+    "t_step_stop": Nt,
+    "t_step_save": Nt,
+    "weno_order": args.order,
+    "weno_eps": 1e-16,
+    "mapped_weno": "F",
+    "wenoz": "F",
+    "teno": "F",
+    "mp_weno": "F",
+    "time_stepper": 3,
+    "riemann_solver": 2,
+    "wave_speeds": 1,
+    "avg_state": 2,
+    "model_eqns": 2,
+    "num_fluids": 1,
+    "fluid_pp(1)%gamma": 1.4,
+    "fluid_pp(1)%pi_inf": 0.0,
+    # x: periodic (domain [0,5] = 5 full periods of the sine IC; wave advects cleanly)
+    # y: symmetry axis at r=0, extrapolation at outer r (same as CI tests)
+    "bc_x%beg": -1,
+    "bc_x%end": -1,
+    "bc_y%beg": -2,
+    "bc_y%end": -3,
+    "num_patches": 1,
+    "patch_icpp(1)%geometry": 3,
+    "patch_icpp(1)%x_centroid": Lx / 2,
+    "patch_icpp(1)%y_centroid": Lr / 2,
+    "patch_icpp(1)%length_x": Lx,
+    "patch_icpp(1)%length_y": Lr,
+    "patch_icpp(1)%vel(1)": 1.0,
+    "patch_icpp(1)%vel(2)": 0.0,
+    "patch_icpp(1)%pres": 1.0,
+    "patch_icpp(1)%alpha_rho(1)": "1.0 + 0.2*sin(2.0*acos(-1.0_wp)*x)",
+    "patch_icpp(1)%alpha(1)": 1.0,
+}
+
+if args.mfc is not None:
+    print(json.dumps(case))
diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp
@@ -775,7 +775,7 @@ contains
                                     ! Geometrical source of the void fraction(s) is zero
                                     $:GPU_LOOP(parallelism='[seq]')
                                     do i = eqn_idx%adv%beg, eqn_idx%adv%end
-                                        flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i)
+                                        flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = 0._wp
                                     end do
                                 end if
 
@@ -1385,7 +1385,7 @@ contains
                                     ! Geometrical source of the void fraction(s) is zero
                                     $:GPU_LOOP(parallelism='[seq]')
                                     do i = eqn_idx%adv%beg, eqn_idx%adv%end
-                                        flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = flux_rs${XYZ}$_vf(j, k, l, i)
+                                        flux_gsrc_rs${XYZ}$_vf(j, k, l, i) = 0._wp
                                     end do
                                 end if
 
diff --git a/src/simulation/m_time_steppers.fpp b/src/simulation/m_time_steppers.fpp
@@ -548,6 +548,20 @@ contains
 
             if (adv_n) call s_comp_alpha_from_n(q_cons_ts(1)%vf)
 
+            ! Single-fluid cyl_coord: alpha is trivially 1 but drifts due to varying HLLC contact velocity across faces. Reset to
+            ! prevent pressure NaN.
+            if (num_fluids == 1 .and. cyl_coord .and. .not. bubbles_euler) then
+                $:GPU_PARALLEL_LOOP(collapse=3)
+                do l = 0, p
+                    do k = 0, n
+                        do j = 0, m
+                            q_cons_ts(1)%vf(eqn_idx%adv%beg)%sf(j, k, l) = 1._wp
+                        end do
+                    end do
+                end do
+                $:END_GPU_PARALLEL_LOOP()
+            end if
+
             if (ib) then
                 ! check if any IBMS are moving, and if so, update the markers, ghost points, levelsets, and levelset norms
                 if (moving_immersed_boundary_flag) then
diff --git a/toolchain/mfc/test/run_convergence_axisym.py b/toolchain/mfc/test/run_convergence_axisym.py
@@ -0,0 +1,134 @@
+#!/usr/bin/env python3
+"""
+Convergence-rate verification for WENO5 on a 2D axisymmetric (cyl_coord=T) grid.
+
+Density sine wave in z: rho = 1 + 0.2*sin(2*pi*z), u_z=1, p=1, u_r=0.
+Exact solution at time T: rho_exact(z,T) = 1 + 0.2*sin(2*pi*(z-T)).
+Nr is held fixed; Nz is refined.
+L2 error is computed by averaging density over r then comparing to exact solution.
+"""
+
+import argparse
+import json
+import math
+import os
+import shutil
+import struct
+import subprocess
+import sys
+import tempfile
+
+import numpy as np
+
+CASE = "examples/2D_axisym_convergence/case.py"
+MFC = "./mfc.sh"
+NR = 32  # fixed radial cells (n=31; WENO5 needs n+1 >= num_stcls_min*5 = 25)
+TFINAL = 0.1  # short final time avoids long-time cylindrical instability
+
+
+def read_field(run_dir: str, step: int, var_idx: int, nz: int, nr: int) -> np.ndarray:
+    """Read 2D field. In MFC: x=axial(Nz), y=radial(Nr). Fortran col-major -> shape (Nz, Nr)."""
+    path = os.path.join(run_dir, "p_all", "p0", str(step), f"q_cons_vf{var_idx}.dat")
+    with open(path, "rb") as f:
+        rec_len = struct.unpack("i", f.read(4))[0]
+        data = np.frombuffer(f.read(rec_len), dtype=np.float64).copy()
+        f.read(4)
+    if data.size != nr * nz:
+        raise ValueError(f"Expected {nr * nz} values, got {data.size}")
+    # Fortran stores (x, y) = (axial, radial) in column-major: first index varies fastest
+    return data.reshape(nz, nr, order="F")
+
+
+def run_case(tmpdir: str, nz: int, extra_args: list) -> tuple:
+    """Run the 2D axisym case at resolution nz. Returns (Nt, run_dir)."""
+    result = subprocess.run(
+        [sys.executable, CASE, "--mfc", "{}", "-N", str(nz), "--nr", str(NR), "--Tfinal", str(TFINAL)] + extra_args,
+        capture_output=True,
+        text=True,
+        check=False,
+    )
+    if result.returncode != 0:
+        raise RuntimeError(f"case.py failed:\n{result.stderr}")
+    cfg = json.loads(result.stdout)
+    Nt = int(cfg["t_step_stop"])
+    dt = float(cfg["dt"])
+
+    cmd = [MFC, "run", CASE, "-t", "pre_process", "simulation", "-n", "1", "--", "-N", str(nz), "--nr", str(NR), "--Tfinal", str(TFINAL)] + extra_args
+    result = subprocess.run(cmd, capture_output=True, text=True, cwd=os.getcwd(), check=False)
+    if result.returncode != 0:
+        print(result.stdout[-3000:])
+        print(result.stderr)
+        raise RuntimeError(f"./mfc.sh run failed for Nz={nz}")
+
+    case_dir = os.path.dirname(CASE)
+    src = os.path.join(case_dir, "p_all")
+    dst = os.path.join(tmpdir, f"N{nz}", "p_all")
+    if os.path.exists(dst):
+        shutil.rmtree(dst)
+    shutil.copytree(src, dst)
+    shutil.rmtree(src, ignore_errors=True)
+    shutil.rmtree(os.path.join(case_dir, "D"), ignore_errors=True)
+    return Nt, dt, os.path.join(tmpdir, f"N{nz}")
+
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--resolutions", type=int, nargs="+", default=[64, 128, 256])
+    parser.add_argument("--cfl", type=float, default=0.02)
+    parser.add_argument("--order", type=int, default=5)
+    args = parser.parse_args()
+
+    extra = ["--cfl", str(args.cfl), "--order", str(args.order)]
+
+    print(f"\n{'=' * 60}")
+    print(f"  WENO{args.order} on 2D axisymmetric (cyl_coord=T) grid")
+    print(f"  Sine wave in z, Nr={NR} fixed, Nz refined, T={TFINAL}")
+    print(f"{'=' * 60}")
+
+    errors, nts = [], []
+    with tempfile.TemporaryDirectory() as tmpdir:
+        for nz in args.resolutions:
+            Lx = 5.0  # must match case.py
+            dz = Lx / nz
+            Nt, dt, run_dir = run_case(tmpdir, nz, extra)
+            nts.append(Nt)
+            T_actual = Nt * dt  # actual final time
+
+            # Read final density field (Nz x Nr), average over r
+            rhoT = read_field(run_dir, Nt, 1, nz, NR)
+            rhoT_z = rhoT.mean(axis=1)  # shape (nz,)
+
+            # Exact solution: rho(z, T) = 1 + 0.2*sin(2*pi*(z - T))
+            x_cc = (np.arange(nz) + 0.5) * dz
+            rho_exact = 1.0 + 0.2 * np.sin(2.0 * np.pi * (x_cc - T_actual))
+
+            err = float(np.sqrt(np.sum((rhoT_z - rho_exact) ** 2) * dz))
+            errors.append(err)
+            print(f"  Nz={nz:4d}: Nt={Nt}, err={err:.4e}")
+
+    Lx = 5.0
+    rates = [None]
+    for i in range(1, len(args.resolutions)):
+        nz0, nz1 = args.resolutions[i - 1], args.resolutions[i]
+        rates.append((math.log(errors[i]) - math.log(errors[i - 1])) / (math.log(Lx / nz1) - math.log(Lx / nz0)))
+
+    Lx = 5.0
+    print(f"\n  {'Nz':>6}  {'Nt':>6}  {'dz':>10}  {'L2 error':>14}  {'rate':>8}")
+    print(f"  {'-' * 6}  {'-' * 6}  {'-' * 10}  {'-' * 14}  {'-' * 8}")
+    for i, nz in enumerate(args.resolutions):
+        r = f"{rates[i]:>8.2f}" if rates[i] is not None else f"{'---':>8}"
+        print(f"  {nz:>6}  {nts[i]:>6}  {Lx / nz:>10.6f}  {errors[i]:>14.6e}  {r}")
+
+    if len(args.resolutions) > 1:
+        log_dz = np.log(Lx / np.array(args.resolutions, dtype=float))
+        log_err = np.log(np.array(errors, dtype=float))
+        rate, _ = np.polyfit(log_dz, log_err, 1)
+        expected = args.order - 0.2
+        passed = rate >= expected
+        print(f"\n  Fitted rate: {rate:.2f}  (need >= {expected:.1f})")
+        print(f"  {'PASS' if passed else 'FAIL'}")
+        sys.exit(0 if passed else 1)
+
+
+if __name__ == "__main__":
+    main()
