test: add cell-shift mode + 3D convergence cases

sbryngelson · sbryngelson · commit d8d604206b75 · 2026-05-10T14:25:09.000-04:00
Cell-shift mode in convergence runner: T = K*h per resolution, compare
q(T) against np.roll(q(0), +K) in each dim (analytical solution for periodic
linear advection with v=1). Cost is O(1) in N (Nt = K*c/CFL independent of
resolution) instead of period mode's O(N) Nt. Forces num_ranks=1 for
unambiguous Fortran-order data layout.

Migrate 1D/2D advection schemes (WENO1/3/5, MUSCL2, TENO5) to cell-shift mode
with K=1, expected_order = scheme_order + 1. WENO7/TENO7 stay in period mode
since cell-shift signal hits machine precision at typical N (h^8 &lt; 1e-15 at
N=64). Add 3D advection case + 5 cell-shift cases (skip WENO7/TENO7).

Empirical rates land cleanly across dimensionalities:
  WENO5  1D 6.00  2D 6.00  3D 6.00  (need &gt;= 5.7)
  WENO3  1D 2.50  2D 2.49  3D 2.47  (need &gt;= 2.2; effective 1.5 at extrema)
  WENO1  1D 2.00  2D 1.99  3D 1.97  (need &gt;= 1.8)
  MUSCL2 1D 3.00  2D 3.01  3D 3.02  (need &gt;= 2.7; unlimited central)
  TENO5  1D 6.00  2D 5.99  3D 5.99  (need &gt;= 5.7)
  WENO7  1D 6.97  2D 6.86         (period mode; need &gt;= 6.5)
  TENO7  1D 6.97  2D 6.86         (period mode; need &gt;= 6.5)

Wall: 1D 2.7 min, 2D 15.5 min (WENO7/TENO7 period dominates), 3D 14.5 min.
diff --git a/examples/1D_euler_convergence/case.py b/examples/1D_euler_convergence/case.py
@@ -26,6 +26,7 @@
 parser.add_argument("--no-mapped", action="store_true", help="Disable mapped WENO")
 parser.add_argument("--muscl-lim", type=int, default=0, help="MUSCL limiter: 0=unlimited 1=minmod ... (default: 0)")
 parser.add_argument("--time-stepper", type=int, default=3, help="Time stepper: 1=Euler 2=RK2 3=RK3 (default: 3)")
+parser.add_argument("--t-end", type=float, default=None, help="Override total simulation time (default: 1.0 = one period)")
 args = parser.parse_args()
 
 gamma = 1.4
@@ -37,8 +38,8 @@
 # c_sound = sqrt(gamma * p / rho) = sqrt(gamma) for p=1, rho=1
 c_max = math.sqrt(gamma) + 1.0  # acoustic + convective
 dt = args.cfl * dx / c_max
-T_end = 1.0
-Nt = max(4, math.ceil(T_end / dt))
+T_end = args.t_end if args.t_end is not None else 1.0
+Nt = max(1, math.ceil(T_end / dt))
 dt = T_end / Nt
 
 if args.muscl:
diff --git a/examples/2D_advection_convergence/case.py b/examples/2D_advection_convergence/case.py
@@ -30,6 +30,7 @@
 parser.add_argument("--no-mapped", action="store_true", help="Disable mapped WENO")
 parser.add_argument("--muscl-lim", type=int, default=0, help="MUSCL limiter: 0=unlimited 1=minmod ... (default: 0)")
 parser.add_argument("--time-stepper", type=int, default=3, help="Time stepper: 1=Euler 2=RK2 3=RK3 (default: 3)")
+parser.add_argument("--t-end", type=float, default=None, help="Override total simulation time (default: 0.5 = one diagonal period)")
 args = parser.parse_args()
 
 gamma = 1.4
@@ -41,8 +42,8 @@
 # Max wave speed: |u| + c with c = sqrt(gamma * p / rho_min). rho_min = 1 - 0.2.
 c_max = math.sqrt(gamma / 0.8) + 1.0
 dt = args.cfl * dx / c_max
-T_end = 0.5  # one full diagonal period: phase shift = 2*T_end = 1.0
-Nt = max(4, math.ceil(T_end / dt))
+T_end = args.t_end if args.t_end is not None else 0.5
+Nt = max(1, math.ceil(T_end / dt))
 dt = T_end / Nt
 
 if args.muscl:
diff --git a/examples/3D_advection_convergence/case.py b/examples/3D_advection_convergence/case.py
@@ -0,0 +1,119 @@
+#!/usr/bin/env python3
+"""
+3D single-fluid Euler convergence case (smooth diagonal advection).
+
+Density wave rho = 1 + 0.2*sin(2*pi*(x+y+z)) with uniform velocity (1, 1, 1)
+on the unit cube with periodic BCs. Pressure p = 1 throughout. The wave
+phase x+y+z - 3t returns to x+y+z after T = 1/3.
+
+Usually invoked from the convergence test framework with --t-end set so
+T = K / N for some integer K (cell-shift mode), allowing analytical
+comparison via integer-cell np.roll of the IC.
+"""
+
+import argparse
+import json
+import math
+
+parser = argparse.ArgumentParser(description="3D Euler diagonal-advection convergence case")
+parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT")
+parser.add_argument("-N", type=int, default=32, help="Grid points per dim (default: 32)")
+parser.add_argument("--order", type=int, default=5, help="WENO order: 1, 3, 5, or 7")
+parser.add_argument("--muscl", action="store_true", help="Use MUSCL-2 instead of WENO")
+parser.add_argument("--teno", action="store_true", help="Use TENO instead of WENO")
+parser.add_argument("--teno-ct", type=float, default=1e-6, help="TENO CT threshold (default: 1e-6)")
+parser.add_argument("--cfl", type=float, default=0.4, help="CFL number (default: 0.4; use 0.005 for WENO7/TENO7)")
+parser.add_argument("--no-mapped", action="store_true", help="Disable mapped WENO")
+parser.add_argument("--muscl-lim", type=int, default=0, help="MUSCL limiter: 0=unlimited 1=minmod ... (default: 0)")
+parser.add_argument("--time-stepper", type=int, default=3, help="Time stepper: 1=Euler 2=RK2 3=RK3 (default: 3)")
+parser.add_argument("--t-end", type=float, default=None, help="Override total simulation time (default: 1/3 = one diagonal period)")
+args = parser.parse_args()
+
+gamma = 1.4
+N = args.N
+m = N - 1
+L = 1.0
+dx = L / N
+
+c_max = math.sqrt(gamma / 0.8) + 1.0
+dt = args.cfl * dx / c_max
+T_end = args.t_end if args.t_end is not None else 1.0 / 3.0
+Nt = max(1, math.ceil(T_end / dt))
+dt = T_end / Nt
+
+if args.muscl:
+    scheme_params = {
+        "recon_type": 2,
+        "muscl_order": 2,
+        "muscl_lim": args.muscl_lim,
+    }
+else:
+    scheme_params = {
+        "recon_type": 1,
+        "weno_order": args.order,
+        "weno_eps": 1.0e-40,
+        "weno_Re_flux": "F",
+        "weno_avg": "F",
+        "mapped_weno": "F" if (args.order == 1 or args.no_mapped or args.teno) else "T",
+        "null_weights": "F",
+        "mp_weno": "F",
+        "teno": "T" if args.teno else "F",
+        **({"teno_CT": args.teno_ct} if args.teno else {}),
+    }
+
+print(
+    json.dumps(
+        {
+            "run_time_info": "F",
+            "x_domain%beg": 0.0,
+            "x_domain%end": L,
+            "y_domain%beg": 0.0,
+            "y_domain%end": L,
+            "z_domain%beg": 0.0,
+            "z_domain%end": L,
+            "m": m,
+            "n": m,
+            "p": m,
+            "dt": dt,
+            "t_step_start": 0,
+            "t_step_stop": Nt,
+            "t_step_save": Nt,
+            "num_patches": 1,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 1,
+            "mpp_lim": "F",
+            "mixture_err": "F",
+            "time_stepper": args.time_stepper,
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -1,
+            "bc_x%end": -1,
+            "bc_y%beg": -1,
+            "bc_y%end": -1,
+            "bc_z%beg": -1,
+            "bc_z%end": -1,
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "F",
+            "patch_icpp(1)%geometry": 9,
+            "patch_icpp(1)%x_centroid": 0.5,
+            "patch_icpp(1)%y_centroid": 0.5,
+            "patch_icpp(1)%z_centroid": 0.5,
+            "patch_icpp(1)%length_x": L,
+            "patch_icpp(1)%length_y": L,
+            "patch_icpp(1)%length_z": L,
+            "patch_icpp(1)%vel(1)": 1.0,
+            "patch_icpp(1)%vel(2)": 1.0,
+            "patch_icpp(1)%vel(3)": 1.0,
+            "patch_icpp(1)%pres": 1.0,
+            "patch_icpp(1)%alpha_rho(1)": "1.0 + 0.2 * sin(2.0 * pi * (x / lx + y / ly + z / lz))",
+            "patch_icpp(1)%alpha(1)": 1.0,
+            "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0),
+            "fluid_pp(1)%pi_inf": 0.0,
+            **scheme_params,
+        }
+    )
+)
diff --git a/toolchain/mfc/test/cases.py b/toolchain/mfc/test/cases.py
@@ -12,41 +12,56 @@
 # the filter handle (`./mfc.sh test --only Convergence`); convergence cases
 # are skipped by default.
 
-# 1D Euler advection (rho = 1 + 0.2*sin(2*pi*x), u=1, p=1, T=1):
-# WENO5/TENO5 use CFL=0.02 so RK3 temporal floor is below O(h^5) at N>=128.
-# WENO7/TENO7 cap at N=128 and use CFL=0.005 (machine-precision floor near N=512).
-# WENO3-JS degrades to 2nd order at smooth extrema (Henrick et al. 2005).
-# MUSCL2 uses unlimited slope (limiters clip to 1st order at smooth extrema).
+# Advection convergence cases use cell-shift mode by default: T = K*h per
+# resolution, compare q(T) to np.roll(q(0), -K) per dim. Spatial error scales
+# as T*h^p = h^(p+1), so measured rate = p+1 (p = scheme order). Wins ~10-100x
+# vs running a full period since Nt = K*c/CFL is independent of N.
+# WENO7/TENO7 stay in period mode: at typical N their cell-shift error hits
+# machine precision (h^8 < 1e-15 at N=64) before any rate signal develops.
 _CONS_VARS_1D = [("density", 1), ("x-momentum", 2), ("energy", 3)]
-_RES_1D_DEFAULT = [64, 128, 256, 512, 1024]
+_CONS_VARS_2D = [("density", 1), ("energy", 4)]
+_CONS_VARS_3D = [("density", 1), ("energy", 5)]
+
+# (label, extra_args, expected_order, tol, resolutions)
+# expected_order = scheme order p + 1 in cell-shift mode (T scales with h).
+# WENO3-JS reduces to 2 at smooth extrema → cell-shift expected = 3.
+# MUSCL2 unlimited central → effective order 2 → cell-shift expected = 3.
 _CONVERGENCE_1D_SCHEMES = [
-    ("WENO5", ["--order", "5", "--cfl", "0.02"], 5, 0.2, 128, 512),
-    ("WENO3", ["--order", "3", "--cfl", "0.02"], 2, 0.2, 256, None),
-    ("WENO1", ["--order", "1", "--cfl", "0.02"], 1, 0.05, 128, None),
-    ("MUSCL2", ["--muscl", "--muscl-lim", "0", "--cfl", "0.02"], 2, 0.1, 128, None),
-    ("TENO5", ["--order", "5", "--teno", "--teno-ct", "1e-6", "--cfl", "0.02"], 5, 0.2, 128, 512),
-    ("WENO7", ["--order", "7", "--cfl", "0.005"], 7, 0.5, 64, 128),
-    ("TENO7", ["--order", "7", "--teno", "--teno-ct", "1e-9", "--cfl", "0.005"], 7, 0.5, 64, 128),
+    ("WENO5", ["--order", "5", "--cfl", "0.02"], 6, 0.3, [32, 64, 128]),
+    ("WENO3", ["--order", "3", "--cfl", "0.02"], 2.5, 0.3, [64, 128, 256]),
+    ("WENO1", ["--order", "1", "--cfl", "0.02"], 2, 0.2, [64, 128, 256]),
+    ("MUSCL2", ["--muscl", "--muscl-lim", "0", "--cfl", "0.02"], 3, 0.3, [64, 128, 256]),
+    ("TENO5", ["--order", "5", "--teno", "--teno-ct", "1e-6", "--cfl", "0.02"], 6, 0.3, [32, 64, 128]),
+]
+# WENO7/TENO7 in 1D: period mode (full period T=1.0, see 1D case.py).
+_CONVERGENCE_1D_PERIOD_SCHEMES = [
+    ("WENO7", ["--order", "7", "--cfl", "0.005"], 7, 0.5, [64, 128]),
+    ("TENO7", ["--order", "7", "--teno", "--teno-ct", "1e-9", "--cfl", "0.005"], 7, 0.5, [64, 128]),
 ]
 
-# 2D diagonal advection (rho = 1 + 0.2 sin(2pi(x+y)), v=(1,1), period T=0.5):
-# linear primitives mean no covariance floor — clean asymptotic rates for all
-# WENO/MUSCL/TENO orders. WENO3 reduces to 2nd at smooth extrema (Henrick 2005).
-# 4 ranks (2x2 decomp): per-rank stencil constraint n+1 >= 5*weno_order/2.
-# WENO5/TENO5 use [64, 96, 128] (33>=25 at N=64). WENO7/TENO7 use [80, 96]
-# (41>=35 at N=80) and CFL=0.005 to push RK3 temporal error below O(h^7).
-# (label, extra_args, expected_order, tol, resolutions)
-_CONS_VARS_2D = [("density", 1), ("energy", 4)]
 _CONVERGENCE_2D_SCHEMES = [
-    ("WENO5", ["--order", "5", "--cfl", "0.02"], 5, 0.3, [64, 96, 128]),
-    ("WENO3", ["--order", "3", "--cfl", "0.02"], 2, 0.3, [32, 64, 128]),
-    ("WENO1", ["--order", "1", "--cfl", "0.02"], 1, 0.2, [32, 64, 128]),
-    ("MUSCL2", ["--muscl", "--muscl-lim", "0", "--cfl", "0.02"], 2, 0.1, [32, 64, 128]),
-    ("TENO5", ["--order", "5", "--teno", "--teno-ct", "1e-6", "--cfl", "0.02"], 5, 0.3, [64, 96, 128]),
+    ("WENO5", ["--order", "5", "--cfl", "0.02"], 6, 0.3, [32, 64, 96]),
+    ("WENO3", ["--order", "3", "--cfl", "0.02"], 2.5, 0.3, [32, 64, 128]),
+    ("WENO1", ["--order", "1", "--cfl", "0.02"], 2, 0.2, [32, 64, 128]),
+    ("MUSCL2", ["--muscl", "--muscl-lim", "0", "--cfl", "0.02"], 3, 0.3, [32, 64, 128]),
+    ("TENO5", ["--order", "5", "--teno", "--teno-ct", "1e-6", "--cfl", "0.02"], 6, 0.3, [32, 64, 96]),
+]
+_CONVERGENCE_2D_PERIOD_SCHEMES = [
     ("WENO7", ["--order", "7", "--cfl", "0.005"], 7, 0.5, [80, 96]),
     ("TENO7", ["--order", "7", "--teno", "--teno-ct", "1e-9", "--cfl", "0.005"], 7, 0.5, [80, 96]),
 ]
 
+# 3D diagonal advection: only cell-shift mode (period T=1/3 with N^3 cells
+# would dominate CI even at N=64). WENO7/TENO7 skipped — at N=64 with K=1
+# the spatial error signal is below machine precision.
+_CONVERGENCE_3D_SCHEMES = [
+    ("WENO5", ["--order", "5", "--cfl", "0.02"], 6, 0.3, [32, 64]),
+    ("WENO3", ["--order", "3", "--cfl", "0.02"], 2.5, 0.3, [32, 64]),
+    ("WENO1", ["--order", "1", "--cfl", "0.02"], 2, 0.2, [32, 64]),
+    ("MUSCL2", ["--muscl", "--muscl-lim", "0", "--cfl", "0.02"], 3, 0.3, [32, 64]),
+    ("TENO5", ["--order", "5", "--teno", "--teno-ct", "1e-6", "--cfl", "0.02"], 6, 0.3, [32, 64]),
+]
+
 # Sod L1 self-convergence: any conservative monotone scheme converges at L1
 # rate ~1 (Godunov). SUPERBEE is over-compressive; min_N=128 skips its
 # pre-asymptotic point.
@@ -71,41 +86,41 @@
     ("RK3", ["--order", "5", "--time-stepper", "3"], 3, 0.3, [0.50, 0.25]),
 ]
 
-# 2D diagonal advection (smooth density wave; no covariance floor): exercises
-# WENO7 / TENO7 in 2D. CFL=0.005 keeps RK3 temporal error below O(h^7).
-# With 4 ranks (2x2 decomp) the per-rank stencil constraint n+1 >= 5*weno_order=35
-# sets min_N=70 (rounded up to 80). N capped at 96 to keep CI cost reasonable:
-# 7th-order error scaling between N=80 and N=96 is ~3.6x, well above the
-# machine-precision floor and clean for rate measurement.
-_RES_2D_ADV_DEFAULT = [80, 96]
-_CONVERGENCE_2D_ADV_SCHEMES = [
-    ("WENO7", ["--order", "7", "--cfl", "0.005"], 7, 0.5, 80, 96),
-    ("TENO7", ["--order", "7", "--teno", "--teno-ct", "1e-9", "--cfl", "0.005"], 7, 0.5, 80, 96),
-]
-
 
 def add_convergence_cases(cases):
     num_ranks = 4
 
-    for label, extra_args, expected, tol, min_N, max_N in _CONVERGENCE_1D_SCHEMES:
+    def _adv_spec(case_path, ndim, cons_vars, extra_args, expected, tol, resolutions, time_mode):
+        # cell_shift forces num_ranks=1 inside the runner; period mode keeps the suite default.
+        return {
+            "runner": f"{ndim}d_advection",
+            "case_path": case_path,
+            "extra_args": extra_args,
+            "expected_order": expected,
+            "tol": tol,
+            "resolutions": resolutions,
+            "ndim": ndim,
+            "domain_len": 1.0,
+            "cons_vars": cons_vars,
+            "primary_idx": 1,
+            "num_ranks": num_ranks,
+            "time_mode": time_mode,
+            "cell_shift": 1,
+        }
+
+    for label, extra_args, expected, tol, resolutions in _CONVERGENCE_1D_SCHEMES:
         cases.append(
             define_convergence_case(
                 f"Convergence -> 1D -> {label}",
-                spec={
-                    "runner": "1d_advection",
-                    "case_path": "examples/1D_euler_convergence/case.py",
-                    "extra_args": extra_args,
-                    "expected_order": expected,
-                    "tol": tol,
-                    "resolutions": _RES_1D_DEFAULT,
-                    "min_N": min_N,
-                    "max_N": max_N,
-                    "ndim": 1,
-                    "domain_len": 1.0,
-                    "cons_vars": _CONS_VARS_1D,
-                    "primary_idx": 1,
-                    "num_ranks": num_ranks,
-                },
+                spec=_adv_spec("examples/1D_euler_convergence/case.py", 1, _CONS_VARS_1D, extra_args, expected, tol, resolutions, "cell_shift"),
+                ppn=1,
+            )
+        )
+    for label, extra_args, expected, tol, resolutions in _CONVERGENCE_1D_PERIOD_SCHEMES:
+        cases.append(
+            define_convergence_case(
+                f"Convergence -> 1D -> {label}",
+                spec=_adv_spec("examples/1D_euler_convergence/case.py", 1, _CONS_VARS_1D, extra_args, expected, tol, resolutions, "period"),
                 ppn=num_ranks,
             )
         )
@@ -114,23 +129,28 @@ def add_convergence_cases(cases):
         cases.append(
             define_convergence_case(
                 f"Convergence -> 2D -> {label}",
-                spec={
-                    "runner": "2d_advection",
-                    "case_path": "examples/2D_advection_convergence/case.py",
-                    "extra_args": extra_args,
-                    "expected_order": expected,
-                    "tol": tol,
-                    "resolutions": resolutions,
-                    "ndim": 2,
-                    "domain_len": 1.0,
-                    "cons_vars": _CONS_VARS_2D,
-                    "primary_idx": 1,
-                    "num_ranks": num_ranks,
-                },
+                spec=_adv_spec("examples/2D_advection_convergence/case.py", 2, _CONS_VARS_2D, extra_args, expected, tol, resolutions, "cell_shift"),
+                ppn=1,
+            )
+        )
+    for label, extra_args, expected, tol, resolutions in _CONVERGENCE_2D_PERIOD_SCHEMES:
+        cases.append(
+            define_convergence_case(
+                f"Convergence -> 2D -> {label}",
+                spec=_adv_spec("examples/2D_advection_convergence/case.py", 2, _CONS_VARS_2D, extra_args, expected, tol, resolutions, "period"),
                 ppn=num_ranks,
             )
         )
 
+    for label, extra_args, expected, tol, resolutions in _CONVERGENCE_3D_SCHEMES:
+        cases.append(
+            define_convergence_case(
+                f"Convergence -> 3D -> {label}",
+                spec=_adv_spec("examples/3D_advection_convergence/case.py", 3, _CONS_VARS_3D, extra_args, expected, tol, resolutions, "cell_shift"),
+                ppn=1,
+            )
+        )
+
     for label, extra_args, expected, tol, min_N in _CONVERGENCE_SOD_SCHEMES:
         cases.append(
             define_convergence_case(
@@ -1747,6 +1767,7 @@ def foreach_example():
                 "1D_advection_convergence",
                 "1D_sod_convergence",
                 "2D_advection_convergence",
+                "3D_advection_convergence",
                 "2D_zero_circ_vortex_analytical",
                 "3D_TaylorGreenVortex_analytical",
                 "3D_IGR_TaylorGreenVortex_nvidia",
diff --git a/toolchain/mfc/test/convergence.py b/toolchain/mfc/test/convergence.py
