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.

Author: sbryngelson

examples/1D_advection_convergence/case.py

@@ -19,6 +19,8 @@
 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
@@ -39,7 +41,7 @@
     scheme_params = {
         "recon_type": 2,
         "muscl_order": 2,
-        "muscl_lim": 1,
+        "muscl_lim": args.muscl_lim,
     }
 else:
     scheme_params = {
@@ -50,7 +52,7 @@
         "weno_avg": "F",
         "mapped_weno": "F" if args.order == 1 else "T",
         "null_weights": "F",
-        "mp_weno": "F",
+        "mp_weno": "T" if args.mp_weno else "F",
     }
 
 print(

examples/1D_euler_convergence/case.py

@@ -0,0 +1,100 @@
+#!/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, 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("--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)")
+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 = 1.0
+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-16,
+        "weno_Re_flux": "F",
+        "weno_avg": "F",
+        "mapped_weno": "F" if (args.order == 1 or args.no_mapped) else "T",
+        "null_weights": "F",
+        "mp_weno": "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": 1,
+            "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)": "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/2D_isentropicvortex_convergence/case.py

@@ -15,6 +15,7 @@
 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, or 5 (default: 5)")
 parser.add_argument("--muscl", action="store_true", help="Use MUSCL instead of WENO")
+parser.add_argument("--muscl-lim", type=int, default=0, help="MUSCL limiter: 0=unlimited 1=minmod ... (default: 0)")
 args = parser.parse_args()
 
 gamma = 1.4
@@ -34,7 +35,7 @@
     scheme_params = {
         "recon_type": 2,
         "muscl_order": 2,
-        "muscl_lim": 1,
+        "muscl_lim": args.muscl_lim,
     }
 else:
     scheme_params = {

src/simulation/m_muscl.fpp

@@ -169,7 +169,9 @@ contains
                                     slopeR = v_rs_ws_${XYZ}$_muscl(j, k, l, i) - v_rs_ws_${XYZ}$_muscl(j - 1, k, l, i)
                                     slope = 0._wp
 
-                                    if (muscl_lim == 1) then  ! minmod
+                                    if (muscl_lim == 0) then  ! unlimited (central difference)
+                                        slope = 5e-1_wp*(slopeL + slopeR)
+                                    else if (muscl_lim == 1) then  ! minmod
                                         if (slopeL*slopeR > muscl_eps) then
                                             slope = min(abs(slopeL), abs(slopeR))
                                         end if

toolchain/mfc/case_validator.py

@@ -133,7 +133,7 @@
     "check_muscl": {
         "title": "MUSCL Reconstruction",
         "category": "Numerical Schemes",
-        "explanation": "muscl_order must be 1 or 2. Second order requires muscl_lim in {1,2,3,4,5}.",
+        "explanation": "muscl_order must be 1 or 2. Second order requires muscl_lim in {0,1,2,3,4,5}.",
     },
     "check_time_stepping": {
         "title": "Time Stepping",
@@ -790,7 +790,7 @@ def check_muscl_simulation(self):
             return
 
         self.prohibit(muscl_order == 2 and muscl_lim is None, "muscl_lim must be defined if using muscl_order = 2")
-        self.prohibit(muscl_lim is not None and (muscl_lim < 1 or muscl_lim > 5), "muscl_lim must be 1, 2, 3, 4, or 5")
+        self.prohibit(muscl_lim is not None and (muscl_lim < 0 or muscl_lim > 5), "muscl_lim must be 0 (unlimited), 1, 2, 3, 4, or 5")
         if muscl_eps is not None:
             self.prohibit(muscl_eps < 0, "muscl_eps must be >= 0 (use 0 for textbook limiter behavior)")
 

toolchain/mfc/params/definitions.py

@@ -609,8 +609,8 @@ def get_value_label(param_name: str, value: int) -> str:
         "value_labels": {1: "1st order", 2: "2nd order"},
     },
     "muscl_lim": {
-        "choices": [1, 2, 3, 4, 5],
-        "value_labels": {1: "minmod", 2: "MC", 3: "Van Albada", 4: "Van Leer", 5: "SUPERBEE"},
+        "choices": [0, 1, 2, 3, 4, 5],
+        "value_labels": {0: "unlimited", 1: "minmod", 2: "MC", 3: "Van Albada", 4: "Van Leer", 5: "SUPERBEE"},
     },
     # Time stepping
     "time_stepper": {

toolchain/mfc/test/run_convergence_1d.py

@@ -1,10 +1,11 @@
 #!/usr/bin/env python3
 """
-Convergence-rate verification for MFC's 1D periodic advection problem.
+Convergence-rate verification for MFC's 1D single-fluid Euler equations.
 
-Two identical fluids (same EOS, rho=1) with a sine-wave volume fraction.
+Single fluid with a density sine wave: rho = 1 + 0.2*sin(2*pi*x), u=1, p=1.
 After exactly one period (T=1, u=1, L=1), the exact solution equals the IC.
-L2(q_cons_vf1(T) - q_cons_vf1(0)) measures the accumulated scheme error.
+L2(rho(T) - rho(0)) measures the accumulated scheme spatial truncation error.
+No non-conservative alpha equation — clean benchmark for all schemes.
 
 CFL=0.02 by default so that RK3 temporal error O(dt^3)~O(h^3) is negligible
 relative to the spatial error at all tested resolutions, allowing WENO5 to
@@ -16,6 +17,10 @@
 the 2D isentropic vortex test (run_convergence.py) verifies WENO3 rate 3 on
 a problem without smooth-extremum contamination.
 
+MUSCL2 uses muscl_lim=0 (unlimited central-difference) by default.  TVD
+limiters clip slopes to zero at smooth extrema and stall at 1st order on the
+sine wave; the unlimited limiter preserves 2nd-order convergence everywhere.
+
 Usage:
     python toolchain/mfc/test/run_convergence_1d.py [--no-build] [--resolutions 32 64 128]
 """
@@ -32,7 +37,7 @@
 
 import numpy as np
 
-CASE = "examples/1D_advection_convergence/case.py"
+CASE = "examples/1D_euler_convergence/case.py"
 MFC = "./mfc.sh"
 
 SCHEMES = [
@@ -168,6 +173,7 @@ def main():
         help="Schemes to test",
     )
     parser.add_argument("--cfl", type=float, default=0.02, help="CFL number (default: 0.02; small so RK3 temporal error is negligible)")
+    parser.add_argument("--muscl-lim", type=int, default=0, help="MUSCL limiter (0=unlimited 1=minmod ...; default: 0)")
     args = parser.parse_args()
 
     if not args.no_build:
@@ -179,7 +185,7 @@ def main():
         if result.returncode != 0:
             sys.exit(1)
 
-    cfl_extra = ["--cfl", str(args.cfl)]
+    cfl_extra = ["--cfl", str(args.cfl), "--muscl-lim", str(args.muscl_lim)]
 
     results = {}
     for label, extra_args, expected_order, tol in SCHEMES: