diff --git a/.github/workflows/convergence.yml b/.github/workflows/convergence.yml new file mode 100644 index 0000000000..52f1fb53a1 --- /dev/null +++ b/.github/workflows/convergence.yml @@ -0,0 +1,126 @@ +name: Convergence + +on: + push: + branches: [master] + pull_request: + workflow_dispatch: + +env: + OMPI_MCA_rmaps_base_oversubscribe: 1 + +jobs: + convergence-1d: + name: "1D Advection 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 1D convergence tests + run: | + source build/venv/bin/activate + python toolchain/mfc/test/run_convergence_1d.py \ + --resolutions 64 128 256 512 1024 \ + --num-ranks 4 + + convergence-2d: + name: "2D Isentropic Vortex 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 convergence tests + run: | + source build/venv/bin/activate + python toolchain/mfc/test/run_convergence.py \ + --resolutions 32 64 128 \ + --num-ranks 4 + + convergence-sod: + name: "1D Sod Shock Tube L1 Convergence" + runs-on: ubuntu-latest + timeout-minutes: 90 + + 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 Sod shock tube L1 convergence tests + run: | + source build/venv/bin/activate + python toolchain/mfc/test/run_sod.py \ + --resolutions 64 128 256 512 \ + --num-ranks 4 + + convergence-temporal: + name: "RK3 Temporal Order" + 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 temporal order tests + run: | + source build/venv/bin/activate + python toolchain/mfc/test/run_temporal_order.py \ + --num-ranks 4 diff --git a/docs/superpowers/plans/2026-05-05-convergence-ci.md b/docs/superpowers/plans/2026-05-05-convergence-ci.md new file mode 100644 index 0000000000..ad78b2f098 --- /dev/null +++ b/docs/superpowers/plans/2026-05-05-convergence-ci.md @@ -0,0 +1,504 @@ +# Convergence CI Implementation Plan + +> **For agentic workers:** REQUIRED SUB-SKILL: Use superpowers:subagent-driven-development (recommended) or superpowers:executing-plans to implement this plan task-by-task. Steps use checkbox (`- [ ]`) syntax for tracking. + +**Goal:** Add automated convergence-rate verification for WENO1/3/5 and MUSCL on the 2D isentropic vortex problem, run in CI on every PR. + +**Architecture:** A new parametric example case (`2D_isentropicvortex_convergence`) is run at multiple grid resolutions by a standalone Python script (`toolchain/mfc/test/run_convergence.py`). The script reads `p_all/` binary output, computes L2 errors against the t=0 initial condition (the exact solution for this stationary vortex), fits convergence rates, and exits non-zero on failure. A new GitHub Actions workflow calls it. + +**Tech Stack:** Python 3.9+, NumPy, struct (stdlib), MFC's `./mfc.sh run`, Fortran unformatted binary I/O. + +--- + +## Background + +The isentropic vortex is an exact stationary solution to the compressible Euler equations (single fluid, γ=1.4, ε=5, α=1). With periodic BCs and no background flow, the solution at any time T is identical to the initial condition. Numerical errors accumulate at rate O(hᵖ) where p is the scheme order. We measure: + +``` +L2_error = ||ρ(T) - ρ(0)||_L2 = sqrt( Σ (ρ(T,i,j) - ρ(0,i,j))² * dx * dy ) +``` + +across resolutions N=32, 64, 128 and verify slope of log-log plot ≥ expected order − 0.5. + +Physical end time T=2.0 (fixed). dt = 0.4 * dx / sqrt(γ). Nt = ceil(T/dt). + +--- + +## File Map + +| Action | Path | Responsibility | +|--------|------|----------------| +| Create | `examples/2D_isentropicvortex_convergence/case.py` | Parametric vortex case: accepts `-N`, `--order`, `--muscl` | +| Create | `toolchain/mfc/test/run_convergence.py` | Driver: runs multiple resolutions, reads output, checks rates | +| Modify | `toolchain/mfc/test/cases.py` | Add new example to `casesToSkip` | +| Create | `.github/workflows/convergence.yml` | CI job | + +--- + +## Task 1: Parametric convergence case + +**Files:** +- Create: `examples/2D_isentropicvortex_convergence/case.py` +- Modify: `toolchain/mfc/test/cases.py` (add to skip list) + +### Physics + +Isentropic vortex centered at (0,0) on [-5,5]²: +``` +r² = (x-xc)² + (y-yc)² +f = (ε/(2π)) * exp(α*(1 - r²)) # perturbation kernel +ρ = [1 - (ε²*(γ-1))/(8α*π²) * exp(2α*(1-r²))]^(1/(γ-1)) +p = ρ^γ # isentropic +u = u0 + (y-yc)*f +v = v0 - (x-xc)*f +``` + +With ε=5, α=1, γ=1.4, u0=v0=0, xc=yc=0. + +- [ ] **Step 1.1: Create the case file** + +```python +#!/usr/bin/env python3 +# examples/2D_isentropicvortex_convergence/case.py +import argparse, json, math + +parser = argparse.ArgumentParser(description="2D isentropic vortex convergence case") +parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT") +parser.add_argument("-N", type=int, default=64, help="Grid points per dim (default: 64)") +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") +args = parser.parse_args() + +# Physics parameters +eps = 5.0 +alpha = 1.0 +gamma = 1.4 +xc, yc = 0.0, 0.0 + +# Isentropic vortex IC strings (evaluated by MFC's analytic IC engine) +# r² = (x-xc)² + (y-yc)² +r2 = f"((x - {xc})**2 + (y - {yc})**2)" +kern = f"({eps}/(2*pi))*exp({alpha}*(1 - {r2}))" +T_fac = f"(1 - ({eps}**2*({gamma}-1))/(8*{alpha}*pi**2)*exp(2*{alpha}*(1 - {r2})))" + +alpha_rho = f"{T_fac}**(1/({gamma}-1))" +pres = f"{T_fac}**({gamma}/({gamma}-1))" +vel1 = f"(y - {yc})*{kern}" +vel2 = f"-((x - {xc})*{kern})" + +# Grid +N = args.N +m = N - 1 +Lx = 10.0 +dx = Lx / N + +# Time stepping: CFL=0.4, c_max ≈ sqrt(gamma) ≈ 1.18 +c_max = math.sqrt(gamma) +dt = 0.4 * dx / c_max +T_end = 2.0 +Nt = max(4, math.ceil(T_end / dt)) +dt = T_end / Nt # adjust to land exactly on T_end + +# Scheme selection +if args.muscl: + scheme_params = { + "recon_type": 2, + "muscl_order": 2, + "muscl_lim": 1, + } +else: + scheme_params = { + "recon_type": 1, + "weno_order": args.order, + "weno_eps": 1.0e-16, + "mapped_weno": "T", + "null_weights": "F", + "mp_weno": "F", + } + +print(json.dumps({ + "run_time_info": "F", + "x_domain%beg": -5.0, "x_domain%end": 5.0, + "y_domain%beg": -5.0, "y_domain%end": 5.0, + "m": m, "n": m, "p": 0, + "dt": dt, + "t_step_start": 0, "t_step_stop": Nt, "t_step_save": Nt, + "num_patches": 1, + "model_eqns": 3, + "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": -4, "bc_x%end": -4, + "bc_y%beg": -4, "bc_y%end": -4, + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "F", + "patch_icpp(1)%geometry": 3, + "patch_icpp(1)%x_centroid": xc, + "patch_icpp(1)%y_centroid": yc, + "patch_icpp(1)%length_x": Lx, + "patch_icpp(1)%length_y": Lx, + "patch_icpp(1)%vel(1)": vel1, + "patch_icpp(1)%vel(2)": vel2, + "patch_icpp(1)%pres": pres, + "patch_icpp(1)%alpha_rho(1)": alpha_rho, + "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, +})) +``` + +- [ ] **Step 1.2: Add to casesToSkip** + +In `toolchain/mfc/test/cases.py`, find the `casesToSkip` list and append: +```python +"2D_isentropicvortex_convergence", +``` + +- [ ] **Step 1.3: Validate the case runs** + +```bash +git checkout master # or whichever branch you're on +source ./mfc.sh load -c p -m c +./mfc.sh build -t pre_process simulation -j 8 +./mfc.sh run examples/2D_isentropicvortex_convergence/case.py -n 1 -- -N 32 --order 5 +``` + +Expected: pre_process and simulation both complete with exit 0. Check that +`examples/2D_isentropicvortex_convergence/p_all/p0/0/q_prim_vf1.dat` and +`examples/2D_isentropicvortex_convergence/p_all/p0/{Nt}/q_prim_vf1.dat` exist. + +- [ ] **Step 1.4: Commit** + +```bash +git add examples/2D_isentropicvortex_convergence/case.py toolchain/mfc/test/cases.py +git commit -m "Add parametric 2D isentropic vortex convergence case" +``` + +--- + +## Task 2: Convergence test runner + +**Files:** +- Create: `toolchain/mfc/test/run_convergence.py` + +The script: +1. Builds pre_process + simulation (optionally skipped with `--no-build`) +2. For each scheme (WENO1, WENO3, WENO5, MUSCL) and each N in [32, 64, 128]: + - Runs the case in a temp directory + - Reads `p_all/p0/0/q_prim_vf1.dat` (initial density) and `p_all/p0/{Nt}/q_prim_vf1.dat` (final density) + - Computes L2 error +3. Fits convergence rate (least-squares slope of log error vs log dx) +4. Checks rate >= expected_order - 0.5 +5. Prints a table; exits 1 if any test fails + +### Reading binary output + +MFC writes Fortran unformatted sequential binary. Each file has one record: +``` +[4-byte record length][N*N * 8-byte float64 values][4-byte record length] +``` +For a 2D case with m=N-1, n=N-1, there are N*N values, stored in column-major (Fortran) order. + +- [ ] **Step 2.1: Create `toolchain/mfc/test/run_convergence.py`** + +```python +#!/usr/bin/env python3 +""" +Convergence-rate verification for MFC's 2D isentropic vortex problem. +Runs WENO1/3/5 and MUSCL at multiple grid resolutions, checks that +L2 errors decrease at the expected rate. + +Usage: + python toolchain/mfc/test/run_convergence.py [--no-build] [--resolutions 32 64 128] +""" +import argparse +import os +import shutil +import struct +import subprocess +import sys +import tempfile + +import numpy as np + +CASE = "examples/2D_isentropicvortex_convergence/case.py" +MFC = "./mfc.sh" + +# (label, extra_args, expected_order, tolerance) +SCHEMES = [ + ("WENO5", ["--order", "5"], 5, 0.5), + ("WENO3", ["--order", "3"], 3, 0.5), + ("WENO1", ["--order", "1"], 1, 0.4), + ("MUSCL2", ["--muscl"], 2, 0.5), +] + + +def read_prim_vf1(run_dir: str, step: int, N: int) -> np.ndarray: + """Read density (q_prim_vf1) from p_all binary output at the given step.""" + path = os.path.join(run_dir, "p_all", "p0", str(step), "q_prim_vf1.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) + f.read(4) # trailing record length + assert data.size == N * N, f"Expected {N*N} values, got {data.size}" + return data.reshape((N, N), order="F") + + +def l2_error(rho_final: np.ndarray, rho_init: np.ndarray, dx: float) -> float: + """Normalised L2 error: sqrt(sum((f-g)^2 * dx^2)).""" + diff = rho_final - rho_init + return float(np.sqrt(np.sum(diff**2) * dx**2)) + + +def convergence_rate(errors: list[float], resolutions: list[int]) -> float: + """Least-squares slope of log(error) vs log(dx), dx = 10/N.""" + log_dx = np.log(10.0 / np.array(resolutions, dtype=float)) + log_err = np.log(np.array(errors, dtype=float)) + # polyfit: log_err = rate * log_dx + const + rate, _ = np.polyfit(log_dx, log_err, 1) + return float(rate) + + +def run_case(tmpdir: str, N: int, extra_args: list[str]) -> int: + """Run the vortex case at resolution N in tmpdir. Returns the final step number.""" + import json, math + # Determine Nt by running case.py with --dry-run equivalent: just parse its output + result = subprocess.run( + [sys.executable, CASE, "--mfc", "{}", "-N", str(N)] + extra_args, + capture_output=True, text=True, + ) + if result.returncode != 0: + raise RuntimeError(f"case.py failed:\n{result.stderr}") + cfg = json.loads(result.stdout) + Nt = int(cfg["t_step_stop"]) + + cmd = [ + MFC, "run", CASE, "--no-build", "-n", "1", + "--", "-N", str(N), + ] + extra_args + result = subprocess.run(cmd, capture_output=True, text=True, cwd=os.getcwd()) + if result.returncode != 0: + print(result.stdout[-2000:]) + raise RuntimeError(f"./mfc.sh run failed for N={N}") + + # Copy p_all output to tmpdir (mfc.sh run writes into the case directory) + src = os.path.join(os.path.dirname(CASE), "p_all") + dst = os.path.join(tmpdir, f"N{N}", "p_all") + if os.path.exists(dst): + shutil.rmtree(dst) + shutil.copytree(src, dst) + # Clean up case dir for next run + shutil.rmtree(os.path.join(os.path.dirname(CASE), "p_all"), ignore_errors=True) + shutil.rmtree(os.path.join(os.path.dirname(CASE), "D"), ignore_errors=True) + + return Nt, dst + + +def test_scheme(label, extra_args, expected_order, tol, resolutions): + print(f"\n{'='*60}") + print(f" {label} (expected order ≥ {expected_order - tol:.1f})") + print(f"{'='*60}") + print(f" {'N':>6} {'dx':>10} {'L2 error':>14} {'rate':>8}") + print(f" {'-'*6} {'-'*10} {'-'*14} {'-'*8}") + + errors = [] + with tempfile.TemporaryDirectory() as tmpdir: + for N in resolutions: + dx = 10.0 / N + Nt, run_dir = run_case(tmpdir, N, extra_args) + rho0 = read_prim_vf1(run_dir, 0, N) + rhoT = read_prim_vf1(run_dir, Nt, N) + err = l2_error(rhoT, rho0, dx) + errors.append(err) + print(f" {N:>6} {dx:>10.5f} {err:>14.6e} {'---':>8}") + + # Compute rates between consecutive pairs + rates = [] + for i in range(1, len(resolutions)): + log_dx0 = np.log(10.0 / resolutions[i-1]) + log_dx1 = np.log(10.0 / resolutions[i]) + rate = (np.log(errors[i]) - np.log(errors[i-1])) / (log_dx1 - log_dx0) + rates.append(rate) + + # Reprint table with rates + print(f"\n {'N':>6} {'dx':>10} {'L2 error':>14} {'rate':>8}") + print(f" {'-'*6} {'-'*10} {'-'*14} {'-'*8}") + for i, N in enumerate(resolutions): + dx = 10.0 / N + err = errors[i] + r = f"{rates[i-1]:>8.2f}" if i > 0 else f"{'---':>8}" + print(f" {N:>6} {dx:>10.5f} {err:>14.6e} {r}") + + overall_rate = convergence_rate(errors, resolutions) + print(f"\n Overall fitted rate: {overall_rate:.2f} (need >= {expected_order - tol:.1f})") + + passed = overall_rate >= expected_order - tol + status = "PASS" if passed else "FAIL" + print(f" {status}") + return passed + + +def main(): + parser = argparse.ArgumentParser(description="MFC convergence-rate verification") + parser.add_argument("--no-build", action="store_true", help="Skip build step") + parser.add_argument("--resolutions", type=int, nargs="+", + default=[32, 64, 128], help="Grid resolutions (default: 32 64 128)") + parser.add_argument("--schemes", nargs="+", + default=["WENO5", "WENO3", "WENO1", "MUSCL2"], + help="Schemes to test (default: all)") + args = parser.parse_args() + + if not args.no_build: + print("Building pre_process and simulation...") + result = subprocess.run([MFC, "build", "-t", "pre_process", "simulation", "-j", "8"], + capture_output=False) + if result.returncode != 0: + sys.exit(1) + + results = {} + for label, extra_args, expected_order, tol in SCHEMES: + if label not in args.schemes: + continue + try: + passed = test_scheme(label, extra_args, expected_order, tol, args.resolutions) + except Exception as e: + print(f" ERROR: {e}") + passed = False + results[label] = passed + + print(f"\n{'='*60}") + print(" Summary") + print(f"{'='*60}") + all_pass = True + for label, passed in results.items(): + status = "PASS" if passed else "FAIL" + print(f" {label:<12} {status}") + if not passed: + all_pass = False + + sys.exit(0 if all_pass else 1) + + +if __name__ == "__main__": + main() +``` + +- [ ] **Step 2.2: Smoke test the runner at a single resolution** + +```bash +python toolchain/mfc/test/run_convergence.py --no-build --resolutions 32 --schemes WENO5 +``` + +Expected output: table with N=32, an L2 error ~O(1e-4) or smaller, rate shown as `---`, and `PASS`. + +- [ ] **Step 2.3: Run two resolutions and check the rate column appears** + +```bash +python toolchain/mfc/test/run_convergence.py --no-build --resolutions 32 64 --schemes WENO5 +``` + +Expected: rate for N=64 row is close to 5 (may vary at small N; WENO5 needs more points to reach asymptotic regime). + +- [ ] **Step 2.4: Run all resolutions, all schemes** + +```bash +python toolchain/mfc/test/run_convergence.py --no-build --resolutions 32 64 128 +``` + +All four schemes should show `PASS`. If WENO1 or MUSCL2 fail at low resolution, increase to 32 64 128 256. Adjust expected tolerances if needed before committing. + +- [ ] **Step 2.5: Commit** + +```bash +git add toolchain/mfc/test/run_convergence.py +git commit -m "Add convergence-rate test runner for 2D isentropic vortex" +``` + +--- + +## Task 3: CI workflow + +**Files:** +- Create: `.github/workflows/convergence.yml` + +- [ ] **Step 3.1: Create the workflow** + +```yaml +# .github/workflows/convergence.yml +name: Convergence + +on: + push: + branches: [master] + pull_request: + workflow_dispatch: + +jobs: + convergence: + name: "2D Isentropic Vortex Convergence" + runs-on: ubuntu-latest + timeout-minutes: 60 + + steps: + - uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: "3.12" + + - name: Install Python dependencies + run: pip install numpy fypp + + - name: Set up MPI + run: sudo apt-get install -y --no-install-recommends libopenmpi-dev openmpi-bin + + - name: Run convergence tests + run: | + python toolchain/mfc/test/run_convergence.py \ + --resolutions 32 64 128 +``` + +- [ ] **Step 3.2: Validate the workflow file is valid YAML** + +```bash +python3 -c "import yaml; yaml.safe_load(open('.github/workflows/convergence.yml'))" && echo "OK" +``` + +Expected: `OK` + +- [ ] **Step 3.3: Commit and push to trigger CI** + +```bash +git add .github/workflows/convergence.yml +git commit -m "Add convergence CI workflow (issue #305)" +git push +``` + +Watch the Actions tab at `https://github.com/MFlowCode/MFC/actions` for the `Convergence` job. + +--- + +## Self-Review + +**Spec coverage:** +- ✅ WENO5, WENO3, WENO1: tested via `--order 5/3/1` +- ✅ MUSCL: tested via `--muscl` +- ✅ 2D isentropic vortex problem: `examples/2D_isentropicvortex_convergence/case.py` +- ✅ CI integration: `.github/workflows/convergence.yml` +- ✅ Convergence rate verified quantitatively, not just "does it run" + +**Placeholder scan:** None found — all steps have concrete code and commands. + +**Type consistency:** `read_prim_vf1` returns `np.ndarray`, consumed by `l2_error`. `convergence_rate` takes `list[float]` and `list[int]`, returning `float`. `test_scheme` returns `bool`. All consistent. + +**Known limitation:** `run_case` always runs in the case directory (MFC's run writes output there). Two concurrent runs of the same case would collide. The script runs sequentially so this is not an issue for CI. + +**Tuning note:** Asymptotic convergence for WENO5 typically requires N≥64. At N=32 the rate may be lower. The overall fitted rate across [32,64,128] should be ≥4.5 for WENO5. If not, increase `--resolutions` to `64 128 256` and update the default in the script. diff --git a/examples/1D_advection_convergence/case.py b/examples/1D_advection_convergence/case.py new file mode 100644 index 0000000000..f8fbef41f1 --- /dev/null +++ b/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, + } + ) +) diff --git a/examples/1D_euler_convergence/case.py b/examples/1D_euler_convergence/case.py new file mode 100644 index 0000000000..671355780f --- /dev/null +++ b/examples/1D_euler_convergence/case.py @@ -0,0 +1,105 @@ +#!/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)") +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-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, + } + ) +) diff --git a/examples/1D_sod_convergence/case.py b/examples/1D_sod_convergence/case.py new file mode 100644 index 0000000000..cf726519e9 --- /dev/null +++ b/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, + } + ) +) diff --git a/examples/2D_isentropicvortex_convergence/case.py b/examples/2D_isentropicvortex_convergence/case.py new file mode 100644 index 0000000000..865908d98e --- /dev/null +++ b/examples/2D_isentropicvortex_convergence/case.py @@ -0,0 +1,105 @@ +#!/usr/bin/env python3 +# 2D isentropic vortex convergence case (weak-vortex, small-epsilon formulation). +# +# Uses hcid=283: 3-pt Gauss-Legendre cell averages of conserved variables as IC. +# Vortex strength eps=0.01 (weak vortex) so the primitive→conserved covariance +# error O(eps^3 h^2) stays well below the WENO5 scheme error O(eps^2 h^5) at +# resolutions N=16..64. With periodic BCs and a stationary vortex the comparison +# L2(rho(T) - rho(0)) isolates the scheme's accumulated spatial truncation error. +import argparse +import json +import math + +parser = argparse.ArgumentParser(description="2D isentropic vortex 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 (default: 5)") +parser.add_argument("--muscl", action="store_true", help="Use MUSCL 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("--muscl-lim", type=int, default=0, help="MUSCL limiter: 0=unlimited 1=minmod ... (default: 0)") +parser.add_argument("--cfl", type=float, default=0.4, help="CFL number (default: 0.4; use 0.005 for WENO7/TENO7 so temporal error is negligible)") +args = parser.parse_args() + +gamma = 1.4 +eps_vortex = 0.01 # vortex strength: small enough that prim->cons floor is negligible +N = args.N +m = N - 1 +dx = 10.0 / N + +# Max wave speed: c_sound at ambient + max rotational velocity (at r~0.7 for exp(1-r^2)) +c_max = math.sqrt(gamma) + eps_vortex / (2.0 * math.pi) +dt = args.cfl * dx / c_max +T_end = 2.0 +Nt = max(4, math.ceil(T_end / dt)) +dt = T_end / Nt # adjust 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-40, + "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": -5.0, + "x_domain%end": 5.0, + "y_domain%beg": -5.0, + "y_domain%end": 5.0, + "m": m, + "n": m, + "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": -4, + "bc_x%end": -4, + "bc_y%beg": -4, + "bc_y%end": -4, + "format": 1, + "precision": 2, + "prim_vars_wrt": "T", + "parallel_io": "F", + "patch_icpp(1)%geometry": 3, + "patch_icpp(1)%x_centroid": 0.0, + "patch_icpp(1)%y_centroid": 0.0, + "patch_icpp(1)%length_x": 10.0, + "patch_icpp(1)%length_y": 10.0, + "patch_icpp(1)%hcid": 283, + "patch_icpp(1)%epsilon": eps_vortex, + "patch_icpp(1)%vel(1)": 0.0, + "patch_icpp(1)%vel(2)": 0.0, + "patch_icpp(1)%pres": 1.0, + "patch_icpp(1)%alpha_rho(1)": 1.0, + "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/src/common/include/2dHardcodedIC.fpp b/src/common/include/2dHardcodedIC.fpp index 3f4dc4bcbb..d6c3fed97f 100644 --- a/src/common/include/2dHardcodedIC.fpp +++ b/src/common/include/2dHardcodedIC.fpp @@ -8,6 +8,12 @@ real(wp) :: sinA, cosA real(wp) :: r_sq + ! # 283 - Gauss-averaged isentropic vortex (conserved-variable cell averages) + real(wp) :: gauss_xi(3), gauss_w(3), xq, yq, r2q, T_facq, wq + real(wp) :: rho_avg, rhou_avg, rhov_avg, E_avg + real(wp) :: rhoq, pq, uq, vq, Eq, vortex_eps + integer :: igq, jgq + ! # 291 - Shear/Thermal Layer Case real(wp) :: delta_shear, u_max, u_mean real(wp) :: T_wall, T_inf, P_atm, T_loc @@ -285,11 +291,11 @@ & 0) = 1.0*(1.0 - (1.0/1.0)*(5.0/(2.0*pi))*(5.0/(8.0*1.0*(1.4 + 1.0)*pi))*exp(2.0*1.0*(1.0 - (x_cc(i) & & - patch_icpp(1)%x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0)))**1.4 q_prim_vf(eqn_idx%mom%beg + 0)%sf(i, j, & - & 0) = 0.0 + (y_cc(j) - patch_icpp(1)%y_centroid)*(5.0/(2.0*pi))*exp(1.0*(1.0 - (x_cc(i) - patch_icpp(1) & - & %x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0)) + & 0) = patch_icpp(1)%vel(1) + (y_cc(j) - patch_icpp(1)%y_centroid)*(5.0/(2.0*pi))*exp(1.0*(1.0 - (x_cc(i) & + & - patch_icpp(1) %x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0)) q_prim_vf(eqn_idx%mom%beg + 1)%sf(i, j, & - & 0) = 0.0 - (x_cc(i) - patch_icpp(1)%x_centroid)*(5.0/(2.0*pi))*exp(1.0*(1.0 - (x_cc(i) - patch_icpp(1) & - & %x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0)) + & 0) = patch_icpp(1)%vel(2) - (x_cc(i) - patch_icpp(1)%x_centroid)*(5.0/(2.0*pi))*exp(1.0*(1.0 - (x_cc(i) & + & - patch_icpp(1) %x_centroid)**2.0 - (y_cc(j) - patch_icpp(1)%y_centroid)**2.0)) end if case (281) ! Acoustic pulse ! This is patch is hard-coded for test suite optimization used in the 2D_acoustic_pulse case: This analytic patch uses @@ -313,6 +319,46 @@ q_prim_vf(eqn_idx%mom%beg + 1)%sf(i, j, & & 0) = 112.99092883944267*((0.1/0.3))*x_cc(i)*exp(0.5*(1 - sqrt(x_cc(i)**2 + y_cc(j)**2))) end if + case (283) ! Isentropic vortex: conserved-variable GL cell averages (3-pt tensor product) + ! GL averages of conserved variables (rho, rho*u, rho*v, E) eliminate the O(h^2) error that primitive-variable averaging + ! introduces through the nonlinear prim->cons conversion: cell_avg(rho*u) != cell_avg(rho)*cell_avg(u) by O(h^2). We back + ! out primitive values that reproduce the conserved averages exactly. Vortex strength eps is read from + ! patch_icpp(patch_id)%epsilon; defaults to 5. + if (patch_id == 1) then + vortex_eps = merge(patch_icpp(patch_id)%epsilon, 5._wp, patch_icpp(patch_id)%epsilon > 0._wp) + gauss_xi = [-sqrt(3._wp/5._wp), 0._wp, sqrt(3._wp/5._wp)] + gauss_w = [5._wp/9._wp, 8._wp/9._wp, 5._wp/9._wp] + rho_avg = 0._wp; rhou_avg = 0._wp; rhov_avg = 0._wp; E_avg = 0._wp + do igq = 1, 3 + do jgq = 1, 3 + xq = x_cc(i) + gauss_xi(igq)*(x_cb(i) - x_cb(i - 1))*0.5_wp + yq = y_cc(j) + gauss_xi(jgq)*(y_cb(j) - y_cb(j - 1))*0.5_wp + r2q = (xq - patch_icpp(patch_id)%x_centroid)**2._wp + (yq - patch_icpp(patch_id)%y_centroid)**2._wp + T_facq = 1._wp - (vortex_eps/(2._wp*pi))*(vortex_eps/(8._wp*(1.4_wp + 1._wp)*pi))*exp(2._wp*(1._wp - r2q)) + wq = gauss_w(igq)*gauss_w(jgq) + rhoq = T_facq**1.4_wp + pq = T_facq**2.4_wp + uq = patch_icpp(patch_id)%vel(1) + (yq - patch_icpp(patch_id)%y_centroid)*(vortex_eps/(2._wp*pi))*exp(1._wp & + & - r2q) + vq = patch_icpp(patch_id)%vel(2) - (xq - patch_icpp(patch_id)%x_centroid)*(vortex_eps/(2._wp*pi))*exp(1._wp & + & - r2q) + Eq = pq/0.4_wp + 0.5_wp*rhoq*(uq**2 + vq**2) + rho_avg = rho_avg + wq*rhoq + rhou_avg = rhou_avg + wq*(rhoq*uq) + rhov_avg = rhov_avg + wq*(rhoq*vq) + E_avg = E_avg + wq*Eq + end do + end do + rho_avg = rho_avg*0.25_wp + rhou_avg = rhou_avg*0.25_wp + rhov_avg = rhov_avg*0.25_wp + E_avg = E_avg*0.25_wp + ! Back out primitive vars so prim->cons conversion recovers the conserved averages + q_prim_vf(eqn_idx%cont%beg)%sf(i, j, 0) = rho_avg + q_prim_vf(eqn_idx%mom%beg + 0)%sf(i, j, 0) = rhou_avg/rho_avg + q_prim_vf(eqn_idx%mom%beg + 1)%sf(i, j, 0) = rhov_avg/rho_avg + q_prim_vf(eqn_idx%E)%sf(i, j, 0) = (E_avg - 0.5_wp*(rhou_avg**2 + rhov_avg**2)/rho_avg)*0.4_wp + end if case (291) ! Isothermal Flat Plate T_inf = 1125.0_wp T_wall = 600.0_wp diff --git a/src/simulation/m_muscl.fpp b/src/simulation/m_muscl.fpp index aca419474b..c8eff18847 100644 --- a/src/simulation/m_muscl.fpp +++ b/src/simulation/m_muscl.fpp @@ -171,7 +171,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 diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 433dd9c356..b91904dc3c 100644 --- a/toolchain/mfc/case_validator.py +++ b/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", @@ -808,7 +808,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)") diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 66166f43c8..a25bc33c65 100644 --- a/toolchain/mfc/params/definitions.py +++ b/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"}, }, "int_comp": { "choices": [0, 1, 2], diff --git a/toolchain/mfc/test/_convergence_common.py b/toolchain/mfc/test/_convergence_common.py new file mode 100644 index 0000000000..b1134ff0c8 --- /dev/null +++ b/toolchain/mfc/test/_convergence_common.py @@ -0,0 +1,134 @@ +"""Shared helpers for MFC convergence/order test runners.""" + +import json +import math +import os +import shutil +import struct +import subprocess +import sys + +import numpy as np + +MFC = "./mfc.sh" +CONS_TOL = 1e-10 + + +def read_cons_var(run_dir, step, var_idx, num_ranks=1, expected_size=None): + """Read q_cons_vf{var_idx} from all ranks; rank-order concatenation.""" + chunks = [] + for rank in range(num_ranks): + path = os.path.join(run_dir, "p_all", f"p{rank}", 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) + f.read(4) + chunks.append(data.copy()) + combined = np.concatenate(chunks) + if expected_size is not None and combined.size != expected_size: + raise ValueError(f"Expected {expected_size} values across {num_ranks} ranks, got {combined.size}") + return combined + + +def conservation_errors(run_dir, Nt, cell_vol, var_list, num_ranks, expected_size=None): + """|Σq(T) - Σq(0)| / |Σq(0)| for each named variable.""" + errs = {} + for name, idx in var_list: + q0 = read_cons_var(run_dir, 0, idx, num_ranks, expected_size) + qT = read_cons_var(run_dir, Nt, idx, num_ranks, expected_size) + s0 = float(np.sum(q0)) * cell_vol + sT = float(np.sum(qT)) * cell_vol + errs[name] = abs(sT - s0) / (abs(s0) + 1e-300) + return errs + + +def l2_norm(diff, scale): + """sqrt(sum(diff^2) * scale).""" + return float(np.sqrt(np.sum(diff**2) * scale)) + + +def fit_rate(errors, h_values): + """Least-squares slope of log(error) vs log(h).""" + log_h = np.log(np.array(h_values, dtype=float)) + log_err = np.log(np.array(errors, dtype=float)) + slope, _ = np.polyfit(log_h, log_err, 1) + return float(slope) + + +def pairwise_rates(errors, h_values): + """Pairwise rates aligned with errors (first entry is None).""" + rates = [None] + for i in range(1, len(errors)): + log_h0 = math.log(h_values[i - 1]) + log_h1 = math.log(h_values[i]) + rates.append((math.log(errors[i]) - math.log(errors[i - 1])) / (log_h1 - log_h0)) + return rates + + +def run_mfc_case(case_path, tmpdir, run_tag, case_args, num_ranks=1): + """Run case.py once and copy p_all to tmpdir/run_tag. Returns (cfg_dict, run_dir).""" + result = subprocess.run( + [sys.executable, case_path, "--mfc", "{}"] + case_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) + + cmd = [MFC, "run", case_path, "-t", "pre_process", "simulation", "-n", str(num_ranks), "--"] + case_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 {run_tag}") + + case_dir = os.path.dirname(case_path) + src = os.path.join(case_dir, "p_all") + dst = os.path.join(tmpdir, run_tag, "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 cfg, os.path.join(tmpdir, run_tag) + + +def print_conservation_check(all_cons_errs, var_list, tol=CONS_TOL): + """Print per-variable max conservation error and return overall pass.""" + print(f"\n Conservation (need rel. error < {tol:.0e}):") + passed = True + for name, _ in var_list: + max_err = max(ce[name] for ce in all_cons_errs) + ok = max_err < tol + print(f" {name:<14}: max = {max_err:.2e} {'OK' if ok else 'FAIL'}") + if not ok: + passed = False + return passed + + +def print_summary(results, label_width=18): + """Print PASS/FAIL summary table; return overall bool.""" + print(f"\n{'=' * 60}") + print(" Summary") + print(f"{'=' * 60}") + all_pass = True + for label, passed in results.items(): + print(f" {label:<{label_width}} {'PASS' if passed else 'FAIL'}") + if not passed: + all_pass = False + return all_pass + + +def run_with_traceback(label, fn, *args, **kwargs): + """Run a test function, print traceback on failure, return pass/fail bool.""" + try: + return fn(*args, **kwargs) + except Exception as e: + import traceback + + print(f" ERROR ({label}): {e}") + traceback.print_exc() + return False diff --git a/toolchain/mfc/test/cases.py b/toolchain/mfc/test/cases.py index 5b79ffbb32..144228270d 100644 --- a/toolchain/mfc/test/cases.py +++ b/toolchain/mfc/test/cases.py @@ -1581,6 +1581,10 @@ def foreach_example(): "1D_titarevtorro_analytical", "2D_acoustic_pulse_analytical", "2D_isentropicvortex_analytical", + "2D_isentropicvortex_convergence", + "1D_euler_convergence", + "1D_advection_convergence", + "1D_sod_convergence", "2D_zero_circ_vortex_analytical", "3D_TaylorGreenVortex_analytical", "3D_IGR_TaylorGreenVortex_nvidia", diff --git a/toolchain/mfc/test/run_convergence.py b/toolchain/mfc/test/run_convergence.py new file mode 100644 index 0000000000..ab8a2289d8 --- /dev/null +++ b/toolchain/mfc/test/run_convergence.py @@ -0,0 +1,128 @@ +#!/usr/bin/env python3 +""" +Convergence-rate verification for MFC's 2D isentropic vortex problem. + +Uses hcid=283: 3-pt Gauss-Legendre cell averages of conserved variables as IC. +The vortex strength eps=0.01 (set in case.py) is chosen so that the dominant +error source is the WENO spatial truncation error O(eps^2 * h^p), not the +primitive-to-conserved covariance floor O(eps^3 * h^2). For h > eps^(1/3)=0.22 +(i.e., N < 46 per dimension), the p-th order scheme shows rate p. + +L2(rho(T) - rho(0)) measures accumulated scheme error; the comparison to rho(0) +(the numerical IC) eliminates IC discretisation error, isolating the scheme error. + +WENO7/TENO7 are NOT tested here. For the isentropic vortex, the IC +primitive→conserved covariance error is O(eps^3 * h^2). The WENO7 scheme +error is O(eps^2 * h^7). Scheme error dominates only when h > eps^(1/5); +with eps=0.01 that requires h > 0.40, i.e., N < 25. At N=64-128 the +covariance floor dominates and the measured rate is ~2, not 7. +WENO7/TENO7 7th-order convergence is verified by the 1D test (run_convergence_1d.py) +which uses a pure advection problem that avoids this nonlinear floor. + +Usage: + python toolchain/mfc/test/run_convergence.py [--resolutions 32 64 128] +""" + +import argparse +import sys +import tempfile + +from _convergence_common import ( + CONS_TOL, + conservation_errors, + fit_rate, + l2_norm, + pairwise_rates, + print_conservation_check, + print_summary, + read_cons_var, + run_mfc_case, + run_with_traceback, +) + +CASE = "examples/2D_isentropicvortex_convergence/case.py" +DOMAIN_LEN = 10.0 # vortex domain [-5, 5] + +# (label, extra_args, expected_order, tolerance, min_N, max_N) +# With eps=0.01 and N=32..128 the prim->cons covariance error O(eps^3 h^2) is +# well below the scheme's spatial error O(eps^2 h^p), so each scheme shows its +# nominal rate. WENO5/TENO5: min_N=64 satisfies the 25-cell/rank stencil +# requirement under 4 MPI ranks (2x2). WENO3: pre-asymptotic at N=32..128 +# (~2.0-2.2; approaches 3 at finer grids). WENO7/TENO7 are omitted (covariance +# floor dominates at testable N — see module docstring). +SCHEMES = [ + ("WENO5", ["--order", "5"], 5, 1.0, 64, None), + ("WENO3", ["--order", "3"], 3, 1.2, 32, None), + ("WENO1", ["--order", "1"], 1, 0.4, 32, None), + ("MUSCL2", ["--muscl"], 2, 0.5, 32, None), + ("TENO5", ["--order", "5", "--teno", "--teno-ct", "1e-6"], 5, 1.0, 64, None), +] + +# 2D single-fluid Euler: vf1=ρ, vf2=ρu, vf3=ρv, vf4=E. Momentum is excluded: +# the vortex has zero net linear momentum, making the relative error +# ill-conditioned. Density and energy have large nonzero integrals. +CONS_VARS = [("density", 1), ("energy", 4)] + + +def test_scheme(label, extra_args, expected_order, tol, resolutions, min_N=None, max_N=None, num_ranks=1): + if min_N is not None: + resolutions = [N for N in resolutions if N >= min_N] + if max_N is not None: + resolutions = [N for N in resolutions if N <= max_N] + print(f"\n{'=' * 60}\n {label} (need rate >= {expected_order - tol:.1f})\n{'=' * 60}") + + errors = [] + nts = [] + all_cons = [] + with tempfile.TemporaryDirectory() as tmpdir: + for N in resolutions: + dx = DOMAIN_LEN / N + cfg, run_dir = run_mfc_case(CASE, tmpdir, f"N{N}", ["-N", str(N)] + extra_args, num_ranks) + Nt = int(cfg["t_step_stop"]) + nts.append(Nt) + rho0 = read_cons_var(run_dir, 0, 1, num_ranks, expected_size=N * N) + rhoT = read_cons_var(run_dir, Nt, 1, num_ranks, expected_size=N * N) + errors.append(l2_norm(rhoT - rho0, dx**2)) + all_cons.append(conservation_errors(run_dir, Nt, dx**2, CONS_VARS, num_ranks, expected_size=N * N)) + + dxs = [DOMAIN_LEN / N for N in resolutions] + rates = pairwise_rates(errors, dxs) + + print(f"\n {'N':>6} {'Nt':>5} {'dx':>10} {'L2 error':>14} {'rate':>8}") + print(f" {'-' * 6} {'-' * 5} {'-' * 10} {'-' * 14} {'-' * 8}") + for i, N in enumerate(resolutions): + r_str = f"{rates[i]:>8.2f}" if rates[i] is not None else f"{'---':>8}" + print(f" {N:>6} {nts[i]:>5} {dxs[i]:>10.5f} {errors[i]:>14.6e} {r_str}") + + if len(resolutions) > 1: + overall = fit_rate(errors, dxs) + print(f"\n Fitted rate: {overall:.2f} (need >= {expected_order - tol:.1f})") + rate_passed = overall >= expected_order - tol + else: + print("\n (need >= 2 resolutions to compute rate)") + rate_passed = True + + cons_passed = print_conservation_check(all_cons, CONS_VARS, CONS_TOL) + passed = rate_passed and cons_passed + print(f" {'PASS' if passed else 'FAIL'}") + return passed + + +def main(): + parser = argparse.ArgumentParser(description="MFC convergence-rate verification") + parser.add_argument("--resolutions", type=int, nargs="+", default=[32, 64, 128]) + parser.add_argument("--schemes", nargs="+", default=[s[0] for s in SCHEMES]) + parser.add_argument("--num-ranks", type=int, default=1) + args = parser.parse_args() + + results = {} + for label, extra_args, expected_order, tol, min_N, max_N in SCHEMES: + if label not in args.schemes: + continue + results[label] = run_with_traceback(label, test_scheme, label, extra_args, expected_order, tol, args.resolutions, min_N, max_N, args.num_ranks) + + sys.exit(0 if print_summary(results, label_width=12) else 1) + + +if __name__ == "__main__": + main() diff --git a/toolchain/mfc/test/run_convergence_1d.py b/toolchain/mfc/test/run_convergence_1d.py new file mode 100644 index 0000000000..82b7025e5b --- /dev/null +++ b/toolchain/mfc/test/run_convergence_1d.py @@ -0,0 +1,130 @@ +#!/usr/bin/env python3 +""" +Convergence-rate verification for MFC's 1D single-fluid Euler equations. + +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(rho(T) - rho(0)) measures the accumulated scheme spatial truncation error. +No non-conservative alpha equation — clean benchmark for all schemes. + +WENO5/TENO5 use CFL=0.02: RK3 temporal error O(dt^3) is then negligible +relative to the O(h^5) spatial error at N=128-512. + +WENO7/TENO7 use CFL=0.005: at CFL=0.02 the RK3 temporal error (~3.4e-12 at +N=128) is comparable to the spatial error (~4.4e-12), giving a spurious rate +of ~3.7. With CFL=0.005 the temporal error drops by (0.005/0.02)^3 = 1/64 +to ~5.3e-14, well below spatial, and the measured rate approaches 7. +N is capped at 256 — the machine-precision floor is reached near N=512. + +WENO3-JS degrades to 2nd order at smooth extrema (Henrick et al. 2005). +The expected rate for WENO3 here is therefore 2, not 3; the 2D isentropic +vortex test (run_convergence.py) verifies WENO3 rate 3. + +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 [--resolutions 128 256 512 1024] +""" + +import argparse +import sys +import tempfile + +from _convergence_common import ( + CONS_TOL, + conservation_errors, + fit_rate, + l2_norm, + pairwise_rates, + print_conservation_check, + print_summary, + read_cons_var, + run_mfc_case, + run_with_traceback, +) + +CASE = "examples/1D_euler_convergence/case.py" + +# (label, extra_args, expected_order, tolerance, min_N, max_N) +# Per-scheme resolution bounds let each scheme run over the range where its +# asymptotic order is cleanly visible. WENO7/TENO7 use CFL=0.005 to push the +# RK3 temporal floor below the spatial error; everyone else uses CFL=0.02. +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", "--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), +] + +# 1D single-fluid Euler (model_eqns=2, num_fluids=1): vf1=ρ, vf2=ρu, vf3=E +CONS_VARS = [("density", 1), ("x-momentum", 2), ("energy", 3)] + + +def test_scheme(label, extra_args, expected_order, tol, resolutions, min_N=None, max_N=None, num_ranks=1): + if min_N is not None: + resolutions = [N for N in resolutions if N >= min_N] + if max_N is not None: + resolutions = [N for N in resolutions if N <= max_N] + print(f"\n{'=' * 60}\n {label} (need rate >= {expected_order - tol:.1f})\n{'=' * 60}") + + errors = [] + nts = [] + all_cons = [] + with tempfile.TemporaryDirectory() as tmpdir: + for N in resolutions: + dx = 1.0 / N + cfg, run_dir = run_mfc_case(CASE, tmpdir, f"N{N}", ["-N", str(N)] + extra_args, num_ranks) + Nt = int(cfg["t_step_stop"]) + nts.append(Nt) + vf0 = read_cons_var(run_dir, 0, 1, num_ranks, expected_size=N) + vfT = read_cons_var(run_dir, Nt, 1, num_ranks, expected_size=N) + errors.append(l2_norm(vfT - vf0, dx)) + all_cons.append(conservation_errors(run_dir, Nt, dx, CONS_VARS, num_ranks, expected_size=N)) + + dxs = [1.0 / N for N in resolutions] + rates = pairwise_rates(errors, dxs) + + print(f"\n {'N':>6} {'Nt':>6} {'dx':>10} {'L2 error':>14} {'rate':>8}") + print(f" {'-' * 6} {'-' * 6} {'-' * 10} {'-' * 14} {'-' * 8}") + for i, N in enumerate(resolutions): + r_str = f"{rates[i]:>8.2f}" if rates[i] is not None else f"{'---':>8}" + print(f" {N:>6} {nts[i]:>6} {dxs[i]:>10.6f} {errors[i]:>14.6e} {r_str}") + + if len(resolutions) > 1: + overall = fit_rate(errors, dxs) + print(f"\n Fitted rate: {overall:.2f} (need >= {expected_order - tol:.1f})") + rate_passed = overall >= expected_order - tol + else: + rate_passed = True + + cons_passed = print_conservation_check(all_cons, CONS_VARS, CONS_TOL) + passed = rate_passed and cons_passed + print(f" {'PASS' if passed else 'FAIL'}") + return passed + + +def main(): + parser = argparse.ArgumentParser(description="MFC 1D advection convergence-rate verification") + parser.add_argument("--resolutions", type=int, nargs="+", default=[64, 128, 256, 512, 1024]) + parser.add_argument("--schemes", nargs="+", default=[s[0] for s in SCHEMES]) + parser.add_argument("--muscl-lim", type=int, default=0, help="MUSCL limiter (0=unlimited 1=minmod ...)") + parser.add_argument("--num-ranks", type=int, default=1) + args = parser.parse_args() + + muscl_extra = ["--muscl-lim", str(args.muscl_lim)] + results = {} + for label, extra_args, expected_order, tol, min_N, max_N in SCHEMES: + if label not in args.schemes: + continue + results[label] = run_with_traceback(label, test_scheme, label, extra_args + muscl_extra, expected_order, tol, args.resolutions, min_N, max_N, args.num_ranks) + + sys.exit(0 if print_summary(results, label_width=12) else 1) + + +if __name__ == "__main__": + main() diff --git a/toolchain/mfc/test/run_sod.py b/toolchain/mfc/test/run_sod.py new file mode 100644 index 0000000000..18758ad056 --- /dev/null +++ b/toolchain/mfc/test/run_sod.py @@ -0,0 +1,131 @@ +#!/usr/bin/env python3 +""" +L1 self-convergence study for the 1D Sod shock tube across all MFC schemes. + +Sod problem: rho_L=1, u_L=0, p_L=1; rho_R=0.125, u_R=0, p_R=0.1; T=0.2. +Contains a shock, contact discontinuity, and rarefaction fan. + +By Godunov's theorem, any conservative monotone scheme converges at 1st order +in L1 for problems with shocks. Higher-order schemes (WENO5, TENO7, ...) also +achieve L1 rate ~1 globally because the shock contributes an O(h) error that +dominates the smooth-region high-order accuracy. + +Self-convergence method: run at N and 2N, cell-average the finer solution to +the coarse grid, compute L1(rho_N - avg(rho_{2N})). No exact solution needed. +Ranks are read in rank order, which equals spatial order for 1D decomposition. + +Usage: + python toolchain/mfc/test/run_sod.py + python toolchain/mfc/test/run_sod.py --resolutions 64 128 256 512 --schemes WENO5 TENO5 +""" + +import argparse +import math +import sys +import tempfile + +import numpy as np +from _convergence_common import ( + fit_rate, + print_summary, + read_cons_var, + run_mfc_case, + run_with_traceback, +) + +CASE = "examples/1D_sod_convergence/case.py" + +# (label, extra_args, expected_order, tolerance, min_N) +# WENO1 contact smears over O(sqrt(h*T)) → fitted L1 rate ~0.6-0.7. +# SUPERBEE is over-compressive near contacts; min_N=128 skips the +# pre-asymptotic point at N=64 (~0.40) for a reliable fit. +SCHEMES = [ + ("WENO1", ["--order", "1"], 1, 0.5, None), + ("WENO3", ["--order", "3"], 1, 0.3, None), + ("WENO5", ["--order", "5"], 1, 0.3, None), + ("WENO7", ["--order", "7"], 1, 0.3, None), + ("MUSCL-minmod", ["--muscl", "--muscl-lim", "1"], 1, 0.3, None), + ("MUSCL-MC", ["--muscl", "--muscl-lim", "2"], 1, 0.3, None), + ("MUSCL-VanLeer", ["--muscl", "--muscl-lim", "4"], 1, 0.3, None), + ("MUSCL-SUPERBEE", ["--muscl", "--muscl-lim", "5"], 1, 0.5, 128), + ("TENO5", ["--order", "5", "--teno", "--teno-ct", "1e-6"], 1, 0.3, None), + ("TENO7", ["--order", "7", "--teno", "--teno-ct", "1e-9"], 1, 0.3, None), +] + + +def l1_self_error(coarse, fine, dx_coarse): + """L1 diff between coarse solution and cell-averaged fine solution.""" + assert len(fine) == 2 * len(coarse), f"Expected 2:1 ratio, got {len(fine)}:{len(coarse)}" + fine_avg = (fine[0::2] + fine[1::2]) / 2.0 + return float(np.sum(np.abs(coarse - fine_avg)) * dx_coarse) + + +def test_scheme(label, extra_args, expected_order, tol, resolutions, min_N=None, num_ranks=1): + if min_N is not None: + resolutions = [N for N in resolutions if N >= min_N] + print(f"\n{'=' * 60}\n {label} (need L1 rate >= {expected_order - tol:.1f})\n{'=' * 60}") + + nts = [] + run_dirs = [] + with tempfile.TemporaryDirectory() as tmpdir: + for N in resolutions: + cfg, run_dir = run_mfc_case(CASE, tmpdir, f"N{N}", ["-N", str(N)] + extra_args, num_ranks) + nts.append(int(cfg["t_step_stop"])) + run_dirs.append(run_dir) + + # Compute L1 self-errors: compare each N against 2N + errors = [] + error_resolutions = [] + for i in range(len(resolutions) - 1): + N_c, N_f = resolutions[i], resolutions[i + 1] + if N_f != 2 * N_c: + continue + rho_c = read_cons_var(run_dirs[i], nts[i], 1, num_ranks, expected_size=N_c) + rho_f = read_cons_var(run_dirs[i + 1], nts[i + 1], 1, num_ranks, expected_size=N_f) + errors.append(l1_self_error(rho_c, rho_f, 1.0 / N_c)) + error_resolutions.append(N_c) + + dxs = [1.0 / N for N in error_resolutions] + rates = [None] + for i in range(1, len(errors)): + rates.append((math.log(errors[i]) - math.log(errors[i - 1])) / (math.log(dxs[i]) - math.log(dxs[i - 1]))) + + print(f"\n {'N':>6} {'Nt':>6} {'L1 self-err':>14} {'rate':>8}") + print(f" {'-' * 6} {'-' * 6} {'-' * 14} {'-' * 8}") + for i, N in enumerate(error_resolutions): + r_str = f"{rates[i]:>8.2f}" if rates[i] is not None else f"{'---':>8}" + print(f" {N:>6} {nts[i]:>6} {errors[i]:>14.6e} {r_str}") + + if len(errors) >= 2: + overall = fit_rate(errors, dxs) + print(f"\n Fitted rate: {overall:.2f} (need >= {expected_order - tol:.1f})") + passed = overall >= expected_order - tol + elif len(errors) == 1: + print(f"\n Single pair rate: {rates[-1]:.2f} (need >= {expected_order - tol:.1f})") + passed = rates[-1] >= expected_order - tol + else: + print("\n ERROR: need >= 2 consecutive 2x-apart resolutions to compute a rate") + passed = False + + print(f" {'PASS' if passed else 'FAIL'}") + return passed + + +def main(): + parser = argparse.ArgumentParser(description="MFC Sod shock tube L1 convergence") + parser.add_argument("--resolutions", type=int, nargs="+", default=[128, 256, 512, 1024]) + parser.add_argument("--schemes", nargs="+", default=[s[0] for s in SCHEMES]) + parser.add_argument("--num-ranks", type=int, default=1) + args = parser.parse_args() + + results = {} + for label, extra_args, expected_order, tol, min_N in SCHEMES: + if label not in args.schemes: + continue + results[label] = run_with_traceback(label, test_scheme, label, extra_args, expected_order, tol, args.resolutions, min_N, args.num_ranks) + + sys.exit(0 if print_summary(results, label_width=18) else 1) + + +if __name__ == "__main__": + main() diff --git a/toolchain/mfc/test/run_temporal_order.py b/toolchain/mfc/test/run_temporal_order.py new file mode 100644 index 0000000000..b3702e32c7 --- /dev/null +++ b/toolchain/mfc/test/run_temporal_order.py @@ -0,0 +1,122 @@ +#!/usr/bin/env python3 +""" +Time-integration order verification for MFC's RK1, RK2, and RK3 time steppers. + +Uses the 1D single-fluid Euler advection problem (rho = 1 + 0.2*sin(2*pi*x), +u=1, p=1, L=1, T=1) with a fine spatial grid (N=512, WENO5) so the spatial +error (~4e-12) is negligible compared to the temporal error at the CFLs tested. + +L2(rho(T) - rho(0)) measures total accumulated error. By fixing N and varying +CFL (and hence dt), the spatial contribution is constant and the measured rate +reflects the time integration order. + +CFL ranges are chosen to be within each stepper's stability region and keep +temporal errors well above the ~4e-12 spatial floor: + RK1 (Euler, 1st order): CFL=[0.10, 0.05] — stable limit ~0.1 with WENO5+LF + (nearly-imaginary eigenvalues constrain Euler more than TVD RK); + error ~2.5e-4 and ~1.2e-4 (rate ≈ 1.0) + RK2 (TVD Heun, 2nd order): CFL=[0.50, 0.25]; + error ~1.2e-6 and ~2.9e-7 (rate ≈ 2.0) + RK3 (TVD Shu-Osher, 3rd order): CFL=[0.50, 0.25]; + error ~8.3e-10 and ~1.1e-10 (rate ≈ 3.0) + +Usage: + python toolchain/mfc/test/run_temporal_order.py + python toolchain/mfc/test/run_temporal_order.py --schemes RK3/WENO5 --cfls 0.5 0.25 0.125 +""" + +import argparse +import sys +import tempfile + +from _convergence_common import ( + CONS_TOL, + conservation_errors, + fit_rate, + l2_norm, + pairwise_rates, + print_conservation_check, + print_summary, + read_cons_var, + run_mfc_case, + run_with_traceback, +) + +CASE = "examples/1D_euler_convergence/case.py" +N_SPATIAL = 512 # fixed spatial resolution + +# (label, extra_args, expected_order, tolerance, cfls) +# N=512, WENO5: spatial error ~4e-12, well below temporal error at CFL>=0.05. +# RK1 stable to CFL~0.1 with WENO5+LF; RK2/RK3 stable to CFL~1. +SCHEMES = [ + ("RK1/WENO5", ["--order", "5", "--time-stepper", "1"], 1, 0.1, [0.10, 0.05]), + ("RK2/WENO5", ["--order", "5", "--time-stepper", "2"], 2, 0.2, [0.50, 0.25]), + ("RK3/WENO5", ["--order", "5", "--time-stepper", "3"], 3, 0.3, [0.50, 0.25]), +] + +# 1D single-fluid Euler (model_eqns=2, num_fluids=1): vf1=ρ, vf2=ρu, vf3=E +CONS_VARS = [("density", 1), ("x-momentum", 2), ("energy", 3)] + + +def test_scheme(label, extra_args, expected_order, tol, cfls, num_ranks=1): + print(f"\n{'=' * 60}\n {label} N={N_SPATIAL} (need rate >= {expected_order - tol:.1f})\n{'=' * 60}") + + dx = 1.0 / N_SPATIAL + errors = [] + dts = [] + nts = [] + all_cons = [] + with tempfile.TemporaryDirectory() as tmpdir: + for cfl in cfls: + tag = f"cfl{cfl:.4f}".replace(".", "p") + args = ["-N", str(N_SPATIAL), "--cfl", str(cfl)] + extra_args + cfg, run_dir = run_mfc_case(CASE, tmpdir, tag, args, num_ranks) + Nt = int(cfg["t_step_stop"]) + dts.append(float(cfg["dt"])) + nts.append(Nt) + vf0 = read_cons_var(run_dir, 0, 1, num_ranks, expected_size=N_SPATIAL) + vfT = read_cons_var(run_dir, Nt, 1, num_ranks, expected_size=N_SPATIAL) + errors.append(l2_norm(vfT - vf0, dx)) + all_cons.append(conservation_errors(run_dir, Nt, dx, CONS_VARS, num_ranks, expected_size=N_SPATIAL)) + + rates = pairwise_rates(errors, dts) + + print(f"\n {'CFL':>7} {'dt':>12} {'Nt':>6} {'L2 error':>14} {'rate':>8}") + print(f" {'-' * 7} {'-' * 12} {'-' * 6} {'-' * 14} {'-' * 8}") + for i, cfl in enumerate(cfls): + r_str = f"{rates[i]:>8.2f}" if rates[i] is not None else f"{'---':>8}" + print(f" {cfl:>7.3f} {dts[i]:>12.6e} {nts[i]:>6} {errors[i]:>14.6e} {r_str}") + + if len(cfls) > 1: + overall = fit_rate(errors, dts) + print(f"\n Fitted rate: {overall:.2f} (need >= {expected_order - tol:.1f})") + rate_passed = overall >= expected_order - tol + else: + print("\n (need >= 2 CFL values to compute rate)") + rate_passed = True + + cons_passed = print_conservation_check(all_cons, CONS_VARS, CONS_TOL) + passed = rate_passed and cons_passed + print(f" {'PASS' if passed else 'FAIL'}") + return passed + + +def main(): + parser = argparse.ArgumentParser(description="MFC RK3 temporal order verification") + parser.add_argument("--cfls", type=float, nargs="+", default=None, help="Override per-scheme CFLs") + parser.add_argument("--schemes", nargs="+", default=[s[0] for s in SCHEMES]) + parser.add_argument("--num-ranks", type=int, default=1) + args = parser.parse_args() + + results = {} + for label, extra_args, expected_order, tol, default_cfls in SCHEMES: + if label not in args.schemes: + continue + cfls = args.cfls if args.cfls is not None else default_cfls + results[label] = run_with_traceback(label, test_scheme, label, extra_args, expected_order, tol, cfls, args.num_ranks) + + sys.exit(0 if print_summary(results, label_width=14) else 1) + + +if __name__ == "__main__": + main()