Add convergence-rate CI for WENO5/3/1 and MUSCL2

.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,
+        }
+    )
+)