Add convergence-rate CI for WENO5/3/1 and MUSCL2 (#1404)

* Add 2D isentropic vortex convergence test and hcid=283 GL IC

* Add convergence CI workflow for 2D isentropic vortex

* Use weak vortex (eps=0.01) to reveal WENO5 rate-5 convergence

With eps=5 the prim->cons covariance error O(eps^3 h^2) swamps WENO5's
O(eps^2 h^5) error at all practical resolutions, giving apparent rate 2.
Using eps=0.01 shifts the crossover to h~0.22 so that N=16..64 (h=0.63..0.16)
are firmly in the WENO5-dominated regime and show rate ~5, as in Johnsen &
Colonius (2006) who similarly used a weak vortex.

- hcid=283: vortex strength now read from patch_icpp(patch_id)%epsilon
  (defaults to 5 if not set, preserving backward compatibility)
- case.py: eps_vortex=0.01, c_max updated accordingly
- run_convergence.py: resolutions 16,32,64; WENO5 tolerance 4.0 (rate>=4)
- convergence.yml: updated to 16 32 64

* Add 1D advection convergence test; split CI into 1D and 2D jobs

CFL=0.02 eliminates RK3 temporal error so WENO5 shows rate 5.
WENO3 expected rate is 2 (not 3) on smooth-extremum problems per
Henrick et al. 2005 — JS smoothness indicators degrade at f'=0.
The 2D vortex job remains the authority on WENO3 rate 3.

* Extend CI resolutions: 1D to N=256, 2D to N=128

* Switch 1D convergence to single-fluid Euler; add muscl_lim=0 unlimited

1D runner now uses examples/1D_euler_convergence/case.py (single-fluid,
density sine wave, no non-conservative alpha equation) instead of the
two-fluid advection case.  The five-equation model non-conservative alpha
transport degrades unlimited MUSCL to ~1st order on two-fluid problems
even with identical EOS; the single-fluid Euler case gives clean rates.

Add muscl_lim=0 (unlimited central-difference) to MUSCL reconstruction:
- src/simulation/m_muscl.fpp: central slope = 0.5*(slopeL+slopeR)
- toolchain/mfc/params/definitions.py: add 0 to choices
- toolchain/mfc/case_validator.py: allow muscl_lim in {0..5}

1D and 2D convergence cases default to muscl_lim=0 for smooth problems
where TVD limiters would stall at smooth extrema.  Runner default CFL=0.02
so RK3 temporal error O(h^3) stays negligible for WENO5 rate-5 verification.

* Fix 2D convergence CI: start at N=32, lower WENO3 threshold to 1.8

N=16 (m=n=15) is below the minimum grid size for WENO5's 5-point stencil
and causes pre_process to abort.  Start all 2D tests at N=32.

WENO3 on the stationary isentropic vortex is pre-asymptotic at N=32-128:
fitted rate ~2.02, local rate 1.83->2.21 (converging toward 3). Lower
threshold from 2.2 to 1.8 (expected=3, tol=1.2).  The rate still clearly
exceeds the 1D result (1.77, Jiang-Shu smooth-extremum degradation), which
is the distinction this test is designed to demonstrate.

* 1D convergence: raise resolutions to 128-1024, tighten WENO1 threshold to 0.95

Default and CI resolutions: 128 256 512 1024.  WENO5 is capped at N=512
per-scheme (hits double-precision floor at N=1024 after 111k steps; error
~2.6e-12 vs 4.2e-12 at N=512, rate collapses to 0.69).  N=128-512 gives
fitted rate 4.99.

WENO1 threshold tightened from 0.6 to 0.95 (tol 0.05): pairwise rates
0.95->0.97->0.99 at N=128-1024, fitted 0.97.

MUSCL2 shows exact rate 2.00 at all doublings N=128-1024.

* Add per-scheme min_N/max_N bounds and tighten convergence thresholds

WENO5 capped at N=512 (machine-precision floor kills rate at N=1024,
error ~2.6e-12, rate collapses to 0.69). WENO3 starts at N=256 to skip
coarse pre-asymptotic points. WENO1 and MUSCL2 use full [128,1024] range.

Thresholds: WENO5>=4.8, WENO3>=1.8, WENO1>=0.95, MUSCL2>=1.9.

* Restrict convergence CI to workflow_dispatch only

Tests are too slow (~30-45 min each job) to run on every PR push.
Manual trigger via workflow_dispatch lets you run them deliberately
when verifying numerics, not on every commit.

* Trigger convergence CI on push to master and workflow_dispatch

* Trigger convergence CI on push to any branch

* Add TENO5, WENO7, TENO7 to convergence runners

Both 1D and 2D runners now test all 7 schemes: WENO5/3/1, MUSCL2, TENO5,
WENO7, TENO7.  Both case.py files gain --teno and --teno-ct flags.

Resolution bounds:
  1D: TENO5 [128,512], WENO7/TENO7 [128,256] (precision floor near N=512)
  2D: TENO5 [32,128],  WENO7/TENO7 [64,128]  (MFC stencil constraint: N>=35)

TENO CT values match the Shu-Osher examples: 1e-6 for order-5, 1e-9 for order-7.
weno_eps tightened to 1e-40 (was 1e-16) for all WENO/TENO schemes in case.py.

* Fix convergence CI: activate MFC venv before running runners

Scripts import numpy which lives in build/venv, not the system Python.

* Fix convergence runners: let MFC build case-specific binaries

Removed --no-build from run_case() in both runners. The generic
./mfc.sh build created a binary without the analytical IC (case.fpp),
so --no-build would silently fall back to that binary and crash
with SIGILL when the IC code wasn't found.

MFC's hash system rebuilds only when case.fpp changes, so each
scheme triggers one build then reuses the binary for all N values.

Also drop WENO7/TENO7 tolerance to 4.0 (threshold >=3.0): local
runs confirm rate ~3.7 due to RK3 temporal and spatial h^7 errors
being comparable at N=128-256.

* Use per-scheme CFL for WENO7/TENO7 to expose 7th-order spatial rate

With CFL=0.02 (1D) or CFL=0.4 (2D) and RK3 time integration, the temporal
error O(dt^3) is comparable to the O(h^7) spatial error at N=128-256, giving
a spurious fitted rate of ~3.7 instead of 7.

Fix: bake CFL=0.005 directly into WENO7/TENO7 extra_args.  This drops the
temporal error by (0.005/0.02)^3=1/64 in 1D and (0.005/0.4)^3=1/51200 in 2D,
making spatial error dominate.  All other schemes keep CFL=0.02 (1D) or 0.4 (2D).

Threshold updated from >=3.0 to >=6.5 (tol 0.5) in 1D and >=6.0 (tol 1.0) in 2D.
Also adds --cfl arg to 2D case.py so the runner can override it per scheme.

* Fix WENO7/TENO7 convergence test resolution windows

1D: use N=64,128 (not N=128,256).  At N=256 the WENO7 spatial error
(~1.7e-14) falls below the round-off accumulation floor (~2.5e-12 for
~28M cell-steps), giving a spurious rate of -0.21.  At N=64 the spatial
error is ~2.8e-10, well above the floor, and the rate between N=64 and
N=128 is 6.97.  Add N=64 to the CI resolution list accordingly.

2D: remove WENO7/TENO7 from the 2D test.  The isentropic vortex has a
primitive->conserved covariance floor O(eps^3*h^2).  For WENO7 scheme
error O(eps^2*h^7) to dominate requires h > eps^(1/5) = 0.40 (N < 25).
At N=64-128 the covariance floor gives rate ~2, masking the 7th-order
spatial rate regardless of CFL.  The 1D test cleanly verifies WENO7/TENO7
7th-order convergence on a pure advection problem without this floor.

* Add --num-ranks to convergence runners; use 4 MPI ranks in CI

Threads num_ranks through run_case() and test_scheme() in both runners.
GitHub ubuntu-latest runners have 4 vCPUs; using all 4 reduces simulation
time for schemes with large Nt (e.g. WENO7/TENO7 at 27k-56k steps).

* Fix readers to collect data from all MPI ranks

read_vf1_1d and read_cons_vf1 were reading only p_all/p0/, so with
num_ranks>1 only the first rank's N/num_ranks cells were included in
the L2 error. Now both readers iterate over p0..p{n-1} and concatenate.
For the L2 norm this is correct regardless of spatial ordering across
ranks. Validated: 4-rank WENO7 |vf0|=64 (not 16) and error matches
single-rank result exactly.

* Add RK3 temporal order verification runner and CI job

* Add RK1/RK2/RK3 temporal order verification; add --time-stepper to 1D case

* Add conservation checks (density, momentum, energy) to all convergence runners

* Fix CI: set OMPI_MCA_rmaps_base_oversubscribe=1 to allow 4 ranks on 2-core runners

* Add Sod shock tube L1 self-convergence runner across all MFC schemes

Tests WENO1/3/5/7, MUSCL-minmod/MC/VanLeer/SUPERBEE, TENO5/7 on the
1D Sod problem (shock + contact + rarefaction). Self-convergence: compare
N vs 2N via cell-averaging, no exact Riemann solver needed. Rate ~1 for
all schemes by Godunov theorem; WENO1 and SUPERBEE use wider tolerance
(0.5) due to contact-dominated O(h^0.5) L1 diffusion at these resolutions.
CI runs at 64/128/256/512 with 4 MPI ranks.

* Fix 2D convergence: set min_N=64 for WENO5/TENO5 and exclude momentum from conservation check

WENO5/TENO5 require (N/2) >= num_stcls_min*weno_order = 25 cells/rank
in a 2x2 decomposition; N=32 only gives 16 cells/rank which fails MFC's
decomposition check. Setting min_N=64 (32 cells/rank) fixes this.

The isentropic vortex has zero net linear momentum by symmetry, making
the relative conservation error formula ill-conditioned (divides near-zero
by 1e-300). Density and energy have large nonzero integrals and are the
meaningful conserved quantities to check.

* Fix Sod runner: add min_N per scheme, set min_N=128 for MUSCL-SUPERBEE

MUSCL-SUPERBEE's pre-asymptotic L1 rate at N=64 is ~0.40 (contact-dominated),
which pulls the fitted rate below the 0.5 threshold. Setting min_N=128 skips
that pre-asymptotic point; the asymptotic rates (N=128,256,512) are ~0.56.

* Fix convergence CI: only run on push to master and PRs

Previously triggered on every push to every branch, running expensive
1D/2D/temporal/Sod jobs (60-90 min total) on every dev push. Now
triggers on: push to master, any pull_request, or workflow_dispatch.

* Skip convergence example dirs in main regression test suite

* Improve error handling in convergence runners: print stderr on failure, add traceback, validate output sizes, fix silent PASS on zero data

* refactor: extract shared helpers from convergence test runners

* Delete docs/superpowers/plans/2026-05-05-convergence-ci.md

* test: integrate convergence runs into ./mfc.sh test

Convergence tests now register as TestCases (kind='convergence') with a
'Convergence -> {1D,2D,Sod,Temporal} -> <scheme>' trace prefix, dispatched
through _handle_case to the runner library in toolchain/mfc/test/convergence.py.
Skipped by default; opt in via --only Convergence.

Drops the four standalone run_*.py runner scripts and the _convergence_common
helper module — their content moved into toolchain/mfc/test/convergence.py.
The convergence.yml workflow shrinks from 4 hand-rolled jobs invoking python
runners to a single ./mfc.sh test --only Convergence call.

* test: add 2D advection WENO7/TENO7 convergence cases

The existing 2D vortex test can't measure 7th-order rate because the
primitive->conserved IC has an O(eps^3 h^2) covariance floor that dominates
WENO7's O(eps^2 h^7) scheme error at testable resolutions (h > eps^(1/5)
needed; eps=0.01 gives N < 25). Add a 2D smooth-advection case with rho =
1 + 0.2 sin(2pi(x+y)) and v=(1,1) on the unit square with periodic BCs;
linear primitives mean no covariance floor.

Cases registered: Convergence -> 2D Adv -> {WENO7, TENO7} at N=[80, 96],
CFL=0.005, 4 ranks, T=0.5 (one diagonal period). Both measure rate 6.86
(threshold 6.5) in ~5-8 min wall time per case.

* test: unify 2D convergence tests onto diagonal advection

Drop the 2D isentropic vortex case (eps=0.01 had a covariance floor that
forced wide tolerances 0.4-1.2 just to pass) and route all 2D convergence
schemes through the diagonal-advection case introduced for WENO7/TENO7.
Linear primitives → no covariance floor → tight 1D-like tolerances and
all 7 schemes (WENO1/3/5/7, MUSCL2, TENO5/7) on the same case.

Measured 2D rates (CFL=0.02 except WENO7/TENO7 at 0.005, T=0.5):
  WENO5  4.99 (need >= 4.7)    MUSCL2 2.01 (need >= 1.9)
  WENO3  1.74 (need >= 1.7)    TENO5  4.99 (need >= 4.7)
  WENO1  0.84 (need >= 0.8)    WENO7  6.86 (need >= 6.5)
                                TENO7  6.86 (need >= 6.5)

Suite wall time on 4 ranks: ~23 min (WENO7/TENO7 dominate at ~6/8 min).

Per-scheme resolution lists replace the min_N/max_N filter pattern (cleaner
when low-order schemes need different N than WENO7). Renames the runner
key '2d_vortex' → '2d_advection'.

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

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 < 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 >= 5.7)
  WENO3  1D 2.50  2D 2.49  3D 2.47  (need >= 2.2; effective 1.5 at extrema)
  WENO1  1D 2.00  2D 1.99  3D 1.97  (need >= 1.8)
  MUSCL2 1D 3.00  2D 3.01  3D 3.02  (need >= 2.7; unlimited central)
  TENO5  1D 6.00  2D 5.99  3D 5.99  (need >= 5.7)
  WENO7  1D 6.97  2D 6.86         (period mode; need >= 6.5)
  TENO7  1D 6.97  2D 6.86         (period mode; need >= 6.5)

Wall: 1D 2.7 min, 2D 15.5 min (WENO7/TENO7 period dominates), 3D 14.5 min.

* test: report scheme spatial order (rate-1) in cell-shift mode

In cell-shift mode T = K*h scales with h, so accumulated truncation error
T*h^p = h^(p+1) gives raw measured rate p+1 — one greater than the scheme's
spatial order. The +1 is an artifact of scaling T with h, not an extra order
of accuracy in the scheme. Subtract 1 from the displayed value so the runner
reports the scheme's canonical spatial order p in both modes; the spec's
expected_order field always means p directly.

Output now self-documents:
  cell-shift: 'Fitted spatial order: 5.00 (need >= 4.7)'
  period:     'Fitted rate: 6.97 (need >= 6.5)'

* test: refactor convergence module to match toolchain style

Apply the toolchain conventions (BenchCase pattern, MFCException, typing.List)
to convergence.py and trim newcomer-hostile bits:

- Replace dict-based spec with @dataclasses.dataclass ConvergenceSpec —
  fields and types are visible at the top of the module instead of scattered
  across .get(default) calls inside each runner.
- Drop the _RUNNERS string registry; spec.runner is now a Callable, so
  cases.py just stores the function reference directly (run_h_sweep,
  run_dt_sweep, run_sod_l1).
- Drop io.StringIO + redirect_stdout; runners build a list of lines and
  '\n'.join it. No stdlib magic.
- Rename runners to verb-noun (run_h_sweep / run_dt_sweep / run_sod_l1).
- Use common.MFCException for run failures (matches the rest of toolchain).
- Add a module docstring that names the three flavours up front.
- Rename time_mode='cell_shift' string flag → cell_shift int (0=period).
- Rename temporal-spec field 'N' → 'fixed_N' for clarity.
- Collapse cases.py registration into one loop driven by a small table of
  (schemes, dim_label, case_path, ndim, cons_vars, cell_shift, ppn) tuples.

Same numerical results pre/post (verified 1D suite + Sod WENO1 + Temporal RK1).

Author: sbryngelson

.github/workflows/convergence.yml

@@ -0,0 +1,36 @@
+name: Convergence
+
+on:
+  push:
+    branches: [master]
+  pull_request:
+  workflow_dispatch:
+
+env:
+  OMPI_MCA_rmaps_base_oversubscribe: 1
+
+jobs:
+  convergence:
+    name: "Convergence"
+    runs-on: ubuntu-latest
+    timeout-minutes: 240
+
+    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: Build MFC
+        run: ./mfc.sh build -j 4
+
+      - name: Run convergence tests
+        run: ./mfc.sh test --only Convergence --no-build -j 1

examples/1D_advection_convergence/case.py

@@ -0,0 +1,103 @@
+#!/usr/bin/env python3
+"""
+1D periodic advection convergence case.
+
+Two identical fluids (same gamma, same density = 1) with a sine-wave volume
+fraction.  Since both EOS are identical, alpha_1 advects passively at u = 1
+with no acoustic coupling.  After exactly one period (T = L/u = 1), the
+exact solution equals the IC, so L2(q_cons_vf1(T) - q_cons_vf1(0)) is
+purely the scheme's accumulated spatial truncation error.
+"""
+
+import argparse
+import json
+import math
+
+parser = argparse.ArgumentParser(description="1D advection convergence case")
+parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT")
+parser.add_argument("-N", type=int, default=64, help="Grid points (default: 64)")
+parser.add_argument("--order", type=int, default=5, help="WENO order: 1, 3, or 5")
+parser.add_argument("--muscl", action="store_true", help="Use MUSCL-2 instead of WENO")
+parser.add_argument("--cfl", type=float, default=0.4, help="CFL number (default: 0.4)")
+parser.add_argument("--mp-weno", action="store_true", help="Enable MP-WENO limiter")
+parser.add_argument("--muscl-lim", type=int, default=1, help="MUSCL limiter: 1=minmod 2=MC 3=VanAlbada 4=VanLeer 5=Superbee")
+args = parser.parse_args()
+
+gamma = 1.4
+N = args.N
+m = N - 1
+L = 1.0
+dx = L / N
+
+# Max wave speed: acoustic speed + convective speed
+# c_sound = sqrt(gamma * p / rho) = sqrt(gamma) ≈ 1.183 (for p=1, rho=1)
+c_max = math.sqrt(gamma) + 1.0
+dt = args.cfl * dx / c_max
+T_end = 1.0  # exactly one period: u=1, L=1
+Nt = max(4, math.ceil(T_end / dt))
+dt = T_end / Nt  # snap to land exactly on T_end
+
+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-16,
+        "weno_Re_flux": "F",
+        "weno_avg": "F",
+        "mapped_weno": "F" if args.order == 1 else "T",
+        "null_weights": "F",
+        "mp_weno": "T" if args.mp_weno else "F",
+    }
+
+print(
+    json.dumps(
+        {
+            "run_time_info": "F",
+            "x_domain%beg": 0.0,
+            "x_domain%end": L,
+            "m": m,
+            "n": 0,
+            "p": 0,
+            "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": 2,
+            "mpp_lim": "F",
+            "mixture_err": "F",
+            "time_stepper": 3,
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -1,
+            "bc_x%end": -1,
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "F",
+            "patch_icpp(1)%geometry": 1,
+            "patch_icpp(1)%x_centroid": 0.5,
+            "patch_icpp(1)%length_x": L,
+            "patch_icpp(1)%vel(1)": 1.0,
+            "patch_icpp(1)%pres": 1.0,
+            "patch_icpp(1)%alpha_rho(1)": "0.5 + 0.2 * sin(2.0 * pi * x / lx)",
+            "patch_icpp(1)%alpha_rho(2)": "0.5 - 0.2 * sin(2.0 * pi * x / lx)",
+            "patch_icpp(1)%alpha(1)": "0.5 + 0.2 * sin(2.0 * pi * x / lx)",
+            "patch_icpp(1)%alpha(2)": "0.5 - 0.2 * sin(2.0 * pi * x / lx)",
+            "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0),
+            "fluid_pp(1)%pi_inf": 0.0,
+            "fluid_pp(2)%gamma": 1.0 / (gamma - 1.0),
+            "fluid_pp(2)%pi_inf": 0.0,
+            **scheme_params,
+        }
+    )
+)

examples/1D_euler_convergence/case.py

@@ -0,0 +1,106 @@
+#!/usr/bin/env python3
+"""
+1D single-fluid Euler convergence case.
+
+Single fluid with a density sine wave: rho = 1 + 0.2*sin(2*pi*x).
+Constant velocity u=1 and pressure p=1.  For this IC, the Euler equations
+reduce to pure advection of all variables at speed u=1.  After exactly one
+period (T = L/u = 1), the exact solution equals the IC, so
+L2(rho(T) - rho(0)) measures the accumulated scheme spatial truncation error.
+
+No non-conservative alpha equation — clean benchmark for WENO/MUSCL rates.
+"""
+
+import argparse
+import json
+import math
+
+parser = argparse.ArgumentParser(description="1D Euler convergence case")
+parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT")
+parser.add_argument("-N", type=int, default=64, help="Grid points (default: 64)")
+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)")
+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
+N = args.N
+m = N - 1
+L = 1.0
+dx = L / N
+
+# 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 = 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:
+    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,
+            "m": m,
+            "n": 0,
+            "p": 0,
+            "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,
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "F",
+            "patch_icpp(1)%geometry": 1,
+            "patch_icpp(1)%x_centroid": 0.5,
+            "patch_icpp(1)%length_x": L,
+            "patch_icpp(1)%vel(1)": 1.0,
+            "patch_icpp(1)%pres": 1.0,
+            "patch_icpp(1)%alpha_rho(1)": "1.0 + 0.2 * sin(2.0 * pi * x / lx)",
+            "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,
+        }
+    )
+)

examples/1D_sod_convergence/case.py

@@ -0,0 +1,105 @@
+#!/usr/bin/env python3
+"""
+1D Sod shock tube convergence case.
+
+Standard Sod problem: rho_L=1, u_L=0, p_L=1; rho_R=0.125, u_R=0, p_R=0.1.
+Discontinuity at x=0.5, gamma=1.4, T_end=0.2 (shock, contact, rarefaction).
+Used for L1 self-convergence study; outflow BCs (-3) at both ends.
+"""
+
+import argparse
+import json
+import math
+
+parser = argparse.ArgumentParser(description="1D Sod shock tube convergence case")
+parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT")
+parser.add_argument("-N", type=int, default=128, help="Grid points (default: 128)")
+parser.add_argument("--order", type=int, default=5, help="WENO order: 1, 3, 5, or 7 (default: 5)")
+parser.add_argument("--muscl", action="store_true", help="Use MUSCL-2 instead of WENO")
+parser.add_argument("--muscl-lim", type=int, default=1, help="MUSCL limiter: 1=minmod 2=MC 4=VanLeer 5=SUPERBEE (default: 1)")
+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.3, help="CFL number (default: 0.3)")
+args = parser.parse_args()
+
+gamma = 1.4
+N = args.N
+m = N - 1
+L = 1.0
+dx = L / N
+
+c_max = math.sqrt(gamma) + 1.0  # conservative: sound speed + max velocity
+dt = args.cfl * dx / c_max
+T_end = 0.2
+Nt = max(4, 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.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,
+            "m": m,
+            "n": 0,
+            "p": 0,
+            "dt": dt,
+            "t_step_start": 0,
+            "t_step_stop": Nt,
+            "t_step_save": Nt,
+            "num_patches": 2,
+            "model_eqns": 2,
+            "alt_soundspeed": "F",
+            "num_fluids": 1,
+            "mpp_lim": "F",
+            "mixture_err": "F",
+            "time_stepper": 3,
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -3,
+            "bc_x%end": -3,
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "F",
+            "patch_icpp(1)%geometry": 1,
+            "patch_icpp(1)%x_centroid": 0.25,
+            "patch_icpp(1)%length_x": 0.5,
+            "patch_icpp(1)%vel(1)": 0.0,
+            "patch_icpp(1)%pres": 1.0,
+            "patch_icpp(1)%alpha_rho(1)": 1.0,
+            "patch_icpp(1)%alpha(1)": 1.0,
+            "patch_icpp(2)%geometry": 1,
+            "patch_icpp(2)%x_centroid": 0.75,
+            "patch_icpp(2)%length_x": 0.5,
+            "patch_icpp(2)%vel(1)": 0.0,
+            "patch_icpp(2)%pres": 0.1,
+            "patch_icpp(2)%alpha_rho(1)": 0.125,
+            "patch_icpp(2)%alpha(1)": 1.0,
+            "fluid_pp(1)%gamma": 1.0 / (gamma - 1.0),
+            "fluid_pp(1)%pi_inf": 0.0,
+            **scheme_params,
+        }
+    )
+)