From 4b0eddb282bcc2b7f3dca27d3e4e238d6e5fc258 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Tue, 5 May 2026 22:39:54 -0400 Subject: [PATCH 01/13] feat: add fp-stability command for Verrou-based FP instability testing MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Adds ./mfc.sh fp-stability — a persistent floating-point stability test suite using Verrou's random IEEE-754 rounding mode. For each registered test case the runner: 1. Generates initial conditions via pre_process 2. Runs simulation once with --rounding-mode=nearest (reference) 3. Runs simulation N times with --rounding-mode=random 4. Reports max L-inf deviation vs threshold (PASS/FAIL) Two cases probe known ill-conditioning in MFC: - sod_strong: 1-D Sod p_L/p_R=100,000 — HLLC xi-factor cancellation (s_L - vel_L)/(s_L - s_S) near sonic contact - water_stiffened: 1-D water shock pi_inf=4046 — pressure recovery p=(E-pi_inf)/gamma loses ~4 decimal digits on low-pressure side Requires a Verrou-enabled Valgrind at $VERROU_HOME/bin/valgrind (default: $HOME/.local/verrou). Silently skips if not found. Binaries are auto-discovered from build/install/ or passed explicitly. --- .gitignore | 1 + .../cases/sod_strong/pre_process.inp | 33 +++ .../cases/sod_strong/simulation.inp | 29 ++ .../cases/water_stiffened/pre_process.inp | 33 +++ .../cases/water_stiffened/simulation.inp | 29 ++ toolchain/main.py | 4 + toolchain/mfc/cli/commands.py | 60 +++++ toolchain/mfc/fp_stability.py | 253 ++++++++++++++++++ 8 files changed, 442 insertions(+) create mode 100644 tests/fp_stability/cases/sod_strong/pre_process.inp create mode 100644 tests/fp_stability/cases/sod_strong/simulation.inp create mode 100644 tests/fp_stability/cases/water_stiffened/pre_process.inp create mode 100644 tests/fp_stability/cases/water_stiffened/simulation.inp create mode 100644 toolchain/mfc/fp_stability.py diff --git a/.gitignore b/.gitignore index aba54411e1..b0fbea5382 100644 --- a/.gitignore +++ b/.gitignore @@ -37,6 +37,7 @@ docs/documentation/parameters.md /tests/*/** !/tests/*/golden.txt !/tests/*/golden-metadata.txt +!/tests/fp_stability/** # NVIDIA Nsight Compute *.nsys-rep diff --git a/tests/fp_stability/cases/sod_strong/pre_process.inp b/tests/fp_stability/cases/sod_strong/pre_process.inp new file mode 100644 index 0000000000..33693f3d58 --- /dev/null +++ b/tests/fp_stability/cases/sod_strong/pre_process.inp @@ -0,0 +1,33 @@ +&user_inputs +x_domain%beg = 0.0 +x_domain%end = 1.0 +m = 24 +n = 0 +p = 0 +t_step_start = 0 +num_patches = 2 +model_eqns = 2 +num_fluids = 1 +mpp_lim = F +weno_order = 5 +bc_x%beg = -3 +bc_x%end = -3 +precision = 2 +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 = 1000.0 +patch_icpp(1)%alpha_rho(1) = 10.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.01 +patch_icpp(2)%alpha_rho(1) = 0.01 +patch_icpp(2)%alpha(1) = 1.0 +fluid_pp(1)%gamma = 2.5000000000000004 +fluid_pp(1)%pi_inf = 0.0 +&end/ diff --git a/tests/fp_stability/cases/sod_strong/simulation.inp b/tests/fp_stability/cases/sod_strong/simulation.inp new file mode 100644 index 0000000000..5e3091e668 --- /dev/null +++ b/tests/fp_stability/cases/sod_strong/simulation.inp @@ -0,0 +1,29 @@ +&user_inputs +run_time_info = F +x_domain%beg = 0.0 +x_domain%end = 1.0 +m = 24 +n = 0 +p = 0 +dt = 0.0001 +t_step_start = 0 +t_step_stop = 5 +t_step_save = 5 +model_eqns = 2 +num_fluids = 1 +mpp_lim = F +mixture_err = F +time_stepper = 3 +weno_order = 5 +weno_eps = 1e-16 +riemann_solver = 2 +wave_speeds = 1 +avg_state = 2 +bc_x%beg = -3 +bc_x%end = -3 +precision = 2 +prim_vars_wrt = T +parallel_io = F +fluid_pp(1)%gamma = 2.5000000000000004 +fluid_pp(1)%pi_inf = 0.0 +&end/ diff --git a/tests/fp_stability/cases/water_stiffened/pre_process.inp b/tests/fp_stability/cases/water_stiffened/pre_process.inp new file mode 100644 index 0000000000..986eb0f81f --- /dev/null +++ b/tests/fp_stability/cases/water_stiffened/pre_process.inp @@ -0,0 +1,33 @@ +&user_inputs +x_domain%beg = 0.0 +x_domain%end = 1.0 +m = 24 +n = 0 +p = 0 +t_step_start = 0 +num_patches = 2 +model_eqns = 2 +num_fluids = 1 +mpp_lim = F +weno_order = 5 +bc_x%beg = -3 +bc_x%end = -3 +precision = 2 +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 = 100.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) = 1.0 +patch_icpp(2)%alpha(1) = 1.0 +fluid_pp(1)%gamma = 0.195312500 +fluid_pp(1)%pi_inf = 4046.31 +&end/ diff --git a/tests/fp_stability/cases/water_stiffened/simulation.inp b/tests/fp_stability/cases/water_stiffened/simulation.inp new file mode 100644 index 0000000000..06a8c554ea --- /dev/null +++ b/tests/fp_stability/cases/water_stiffened/simulation.inp @@ -0,0 +1,29 @@ +&user_inputs +run_time_info = F +x_domain%beg = 0.0 +x_domain%end = 1.0 +m = 24 +n = 0 +p = 0 +dt = 5e-5 +t_step_start = 0 +t_step_stop = 5 +t_step_save = 5 +model_eqns = 2 +num_fluids = 1 +mpp_lim = F +mixture_err = T +time_stepper = 3 +weno_order = 5 +weno_eps = 1e-16 +riemann_solver = 2 +wave_speeds = 1 +avg_state = 2 +bc_x%beg = -3 +bc_x%end = -3 +precision = 2 +prim_vars_wrt = T +parallel_io = F +fluid_pp(1)%gamma = 0.195312500 +fluid_pp(1)%pi_inf = 4046.31 +&end/ diff --git a/toolchain/main.py b/toolchain/main.py index fe3c0ed630..f7b1fa9e70 100644 --- a/toolchain/main.py +++ b/toolchain/main.py @@ -198,6 +198,10 @@ def __run(): from mfc import params_cmd params_cmd.params() + elif cmd == "fp-stability": + from mfc import fp_stability + + fp_stability.fp_stability() if __name__ == "__main__": diff --git a/toolchain/mfc/cli/commands.py b/toolchain/mfc/cli/commands.py index 207ab08dbf..5f240d620b 100644 --- a/toolchain/mfc/cli/commands.py +++ b/toolchain/mfc/cli/commands.py @@ -900,6 +900,65 @@ include_common=["targets", "mfc_config", "jobs", "verbose", "debug_log"], ) +FP_STABILITY_COMMAND = Command( + name="fp-stability", + help="Run floating-point stability tests using Verrou.", + description=( + "Runs each registered test case N times under Verrou's random IEEE-754 " + "rounding mode and compares against a nearest-rounding reference run. " + "Reports the max L∞ deviation and PASS/FAIL against per-case thresholds.\n\n" + "Requires a Verrou-enabled Valgrind at $VERROU_HOME/bin/valgrind " + "(defaults to $HOME/.local/verrou). The simulation and pre_process " + "binaries must be serial (no-MPI, no-GPU) debug builds.\n\n" + "Current test cases:\n" + " sod_strong 1-D Sod p_L/p_R=100,000 — HLLC xi-factor cancellation\n" + " water_stiffened 1-D water shock (pi_inf=4046) — pressure-recovery cancellation\n" + ), + include_common=["mfc_config", "verbose", "debug_log"], + arguments=[ + Argument( + name="sim-binary", + help="Path to a serial simulation binary (debug, no-MPI). Auto-discovered from build/install/ if omitted.", + default=None, + metavar="PATH", + ), + Argument( + name="pre-binary", + help="Path to a serial pre_process binary (no-MPI). Auto-discovered from build/install/ if omitted.", + default=None, + metavar="PATH", + ), + Argument( + name="verrou-binary", + help="Path to a Verrou-enabled valgrind binary. Defaults to $VERROU_HOME/bin/valgrind or $HOME/.local/verrou/bin/valgrind.", + default=None, + metavar="PATH", + ), + Argument( + name="samples", + short="N", + help="Number of random-rounding simulation runs per test case.", + type=int, + default=5, + metavar="N", + ), + ], + examples=[ + Example("./mfc.sh fp-stability", "Auto-discover binaries and run all cases"), + Example( + "./mfc.sh fp-stability --sim-binary build/install/abc123/bin/simulation", + "Specify simulation binary explicitly", + ), + Example("./mfc.sh fp-stability -N 10", "Run 10 random-rounding samples per case"), + ], + key_options=[ + ("--sim-binary PATH", "Serial simulation binary (debug, no-MPI)"), + ("--pre-binary PATH", "Serial pre_process binary"), + ("--verrou-binary PATH", "Verrou-enabled valgrind"), + ("-N, --samples N", "Random-rounding samples per case (default: 5)"), + ], +) + VIZ_COMMAND = Command( name="viz", help="Visualize post-processed MFC output.", @@ -1372,6 +1431,7 @@ BENCH_DIFF_COMMAND, COUNT_COMMAND, COUNT_DIFF_COMMAND, + FP_STABILITY_COMMAND, ], common_sets=[ COMMON_TARGETS, diff --git a/toolchain/mfc/fp_stability.py b/toolchain/mfc/fp_stability.py new file mode 100644 index 0000000000..3d907a4e81 --- /dev/null +++ b/toolchain/mfc/fp_stability.py @@ -0,0 +1,253 @@ +""" +Floating-point stability test suite using Verrou. + +Each test case runs the simulation N times under random IEEE-754 rounding +and measures the max deviation from a nearest-rounding reference. Cases +that exceed their threshold are reported as FAIL. + +Requires: + - Verrou-enabled Valgrind at $VERROU_HOME/bin/valgrind + (default: $HOME/.local/verrou) + - A serial (no-MPI, no-GPU) simulation binary + - A serial pre_process binary (to generate initial conditions) + +Usage: + ./mfc.sh fp-stability --sim-binary PATH --pre-binary PATH + ./mfc.sh fp-stability # auto-discovers newest installed binaries +""" + +import glob +import os +import shutil +import subprocess +import sys +import tempfile +import time + +from .common import MFC_ROOT_DIR, MFCException +from .printer import cons +from .state import ARG + +CASES_DIR = os.path.join(MFC_ROOT_DIR, "tests", "fp_stability", "cases") + +# Each case: +# name - subdirectory under CASES_DIR +# description - human-readable purpose +# compare - list of D/ filenames to compare (relative to save step) +# threshold - max L∞ deviation allowed (in conserved-variable units) +# ill_cond - description of the known ill-conditioning +CASES = [ + { + "name": "sod_strong", + "description": "1-D Sod, p_L/p_R=100,000, ideal gas", + "compare": ["cons.1.00.000005.dat", "cons.3.00.000005.dat"], + "threshold": 1e-10, + "ill_cond": "HLLC xi factor: (s_L - vel_L)/(s_L - s_S) cancels near sonic contact", + }, + { + "name": "water_stiffened", + "description": "1-D water shock, stiffened EOS (pi_inf=4046)", + "compare": ["cons.1.00.000005.dat", "prim.3.00.000005.dat"], + "threshold": 1e-10, + "ill_cond": "Pressure recovery: p = (E - pi_inf)/gamma loses ~4 digits (pi_inf/p_right ~ 40,000)", + }, +] + + +def _find_verrou() -> str: + verrou_home = os.environ.get("VERROU_HOME", os.path.join(os.path.expanduser("~"), ".local", "verrou")) + candidate = os.path.join(verrou_home, "bin", "valgrind") + if os.path.isfile(candidate) and os.access(candidate, os.X_OK): + return candidate + fallback = shutil.which("valgrind") + return fallback or "" + + +def _find_binary(name: str) -> str: + install_dir = os.path.join(os.getcwd(), "build", "install") + candidates = glob.glob(os.path.join(install_dir, "*", "bin", name)) + if not candidates: + return "" + return max(candidates, key=os.path.getmtime) + + +def _exclude_args(exclude_file: str) -> list: + """Return --exclude flag if a verrou exclusion file is provided.""" + if exclude_file and os.path.isfile(exclude_file): + return [f"--exclude={exclude_file}"] + return [] + + +def _run_preprocess(pp_bin: str, case_dir: str, work_dir: str): + """Generate initial conditions in work_dir.""" + shutil.copy2(os.path.join(case_dir, "pre_process.inp"), work_dir) + with open(os.path.join(work_dir, "pre.log"), "w") as f: + result = subprocess.run( + [pp_bin], + cwd=work_dir, + stdout=f, + stderr=subprocess.STDOUT, + check=False, + ) + if result.returncode != 0: + raise MFCException(f"pre_process failed (rc={result.returncode}). See {work_dir}/pre.log") + + +def _run_simulation_verrou( + verrou_bin: str, + sim_bin: str, + work_dir: str, + run_dir: str, + rounding_mode: str, + exclude_args: list, +): + """Copy IC files into a fresh tmpdir, run simulation under verrou, collect D/ output.""" + with tempfile.TemporaryDirectory(prefix="mfc-fps-") as tmpdir: + # Copy static inputs + for fname in ["simulation.inp", "indices.dat", "pre_time_data.dat", "io_time_data.dat"]: + src = os.path.join(work_dir, fname) + if os.path.exists(src): + shutil.copy2(src, tmpdir) + + # Copy ICs (p_all_ic → p_all) + shutil.copytree( + os.path.join(work_dir, "p_all"), + os.path.join(tmpdir, "p_all"), + ) + os.makedirs(os.path.join(tmpdir, "D")) + + log_path = os.path.join(run_dir, "verrou.log") + cmd = [ + verrou_bin, + "--tool=verrou", + f"--rounding-mode={rounding_mode}", + "--error-limit=no", + f"--log-file={log_path}", + *exclude_args, + sim_bin, + ] + + with open(os.path.join(run_dir, "sim.out"), "w") as f: + result = subprocess.run(cmd, cwd=tmpdir, stdout=f, stderr=subprocess.STDOUT, check=False) + + if result.returncode != 0: + raise MFCException(f"simulation (rounding={rounding_mode}) exited {result.returncode}. See {run_dir}/sim.out") + + # Collect D/ outputs + os.makedirs(run_dir, exist_ok=True) + d_src = os.path.join(tmpdir, "D") + for fn in os.listdir(d_src): + shutil.copy2(os.path.join(d_src, fn), run_dir) + + +def _max_diff(ref_dir: str, run_dir: str, compare_files: list) -> float: + try: + import numpy as np + except ImportError: + raise MFCException("numpy is required for fp-stability comparison") + + total = 0.0 + for fname in compare_files: + ref_path = os.path.join(ref_dir, fname) + run_path = os.path.join(run_dir, fname) + if not os.path.exists(ref_path): + raise MFCException(f"Reference output missing: {ref_path}") + if not os.path.exists(run_path): + raise MFCException(f"Run output missing: {run_path}") + ref = np.loadtxt(ref_path)[:, 1] + run = np.loadtxt(run_path)[:, 1] + total = max(total, float(np.max(np.abs(ref - run)))) + return total + + +def _run_case(case: dict, verrou_bin: str, sim_bin: str, pp_bin: str, n_samples: int) -> bool: + """Run one stability test case. Returns True if the case passes.""" + name = case["name"] + threshold = case["threshold"] + compare = case["compare"] + case_dir = os.path.join(CASES_DIR, name) + + cons.print(f"[bold]{name}[/bold]: {case['description']}") + cons.indent() + cons.print(f" ill-conditioning: {case['ill_cond']}") + cons.print(f" threshold: {threshold:.0e}") + + work_dir = tempfile.mkdtemp(prefix=f"mfc-fps-{name}-") + try: + # Generate ICs + cons.print(" [dim]running pre_process...[/dim]") + shutil.copy2(os.path.join(case_dir, "simulation.inp"), work_dir) + _run_preprocess(pp_bin, case_dir, work_dir) + + # Reference run (nearest rounding — deterministic baseline) + ref_dir = os.path.join(work_dir, "ref") + os.makedirs(ref_dir) + cons.print(" [dim]reference run (rounding=nearest)...[/dim]") + _run_simulation_verrou(verrou_bin, sim_bin, work_dir, ref_dir, "nearest", []) + + # Random-rounding samples + max_dev = 0.0 + cons.print(f" [dim]random-rounding runs (N={n_samples})...[/dim]") + for i in range(n_samples): + run_dir = os.path.join(work_dir, f"run_{i:02d}") + os.makedirs(run_dir) + _run_simulation_verrou(verrou_bin, sim_bin, work_dir, run_dir, "random", []) + dev = _max_diff(ref_dir, run_dir, compare) + max_dev = max(max_dev, dev) + + passed = max_dev <= threshold + + tag = "[bold green]PASS[/bold green]" if passed else "[bold red]FAIL[/bold red]" + cons.print(f" {tag} max_dev={max_dev:.3e} threshold={threshold:.0e}") + finally: + shutil.rmtree(work_dir, ignore_errors=True) + + cons.unindent() + cons.print() + return passed + + +def fp_stability(): + verrou_bin = ARG("verrou_binary") or _find_verrou() + if not verrou_bin or not os.path.isfile(verrou_bin): + cons.print("[bold yellow]SKIP[/bold yellow]: verrou not found. Install at $HOME/.local/verrou or set VERROU_HOME.") + sys.exit(0) + + sim_bin = ARG("sim_binary") or _find_binary("simulation") + if not sim_bin or not os.path.isfile(sim_bin): + raise MFCException("simulation binary not found. Build with --debug --no-mpi, or pass --sim-binary.") + + pp_bin = ARG("pre_binary") or _find_binary("pre_process") + if not pp_bin or not os.path.isfile(pp_bin): + raise MFCException("pre_process binary not found. Build with --no-mpi, or pass --pre-binary.") + + n_samples = ARG("samples") or 5 + + cons.print() + cons.print("[bold]MFC Floating-Point Stability Suite[/bold]") + cons.print(f" verrou: {verrou_bin}") + cons.print(f" simulation: {sim_bin}") + cons.print(f" pre_process: {pp_bin}") + cons.print(f" samples: {n_samples}") + cons.print() + + start = time.time() + results = [] + for case in CASES: + try: + ok = _run_case(case, verrou_bin, sim_bin, pp_bin, n_samples) + except MFCException as exc: + cons.print(f" [bold red]ERROR[/bold red]: {exc}") + ok = False + results.append((case["name"], ok)) + + elapsed = time.time() - start + n_pass = sum(1 for _, ok in results if ok) + n_fail = len(results) - n_pass + + cons.print(f"[bold]Results[/bold] ({elapsed:.0f}s): [green]{n_pass} passed[/green] [red]{n_fail} failed[/red]") + for name, ok in results: + mark = "[green]✓[/green]" if ok else "[red]✗[/red]" + cons.print(f" {mark} {name}") + + sys.exit(0 if n_fail == 0 else 1) From ae917e12a8e5fc2d4095931f47654f0c5e944fa5 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Tue, 5 May 2026 22:46:02 -0400 Subject: [PATCH 02/13] ci: add weekly FP stability workflow using Verrou --- .github/workflows/fp-stability.yml | 83 ++++++++++++++++++++++++++++++ 1 file changed, 83 insertions(+) create mode 100644 .github/workflows/fp-stability.yml diff --git a/.github/workflows/fp-stability.yml b/.github/workflows/fp-stability.yml new file mode 100644 index 0000000000..5b98f1f261 --- /dev/null +++ b/.github/workflows/fp-stability.yml @@ -0,0 +1,83 @@ +name: FP Stability + +# Runs the Verrou-based floating-point stability suite. +# +# What it tests (./mfc.sh fp-stability): +# sod_strong 1-D Sod shock, p_L/p_R=100,000, ideal gas +# 25 cells, 5 steps, WENO5 + HLLC +# Compares density & energy vs. threshold 1e-10 +# Probes: HLLC xi-factor cancellation near sonic contact +# (s_L - vel_L)/(s_L - s_S) when s_L ≈ s_S +# +# water_stiffened 1-D water shock, stiffened EOS (pi_inf=4046) +# 25 cells, 5 steps, WENO5 + HLLC +# Compares density & pressure vs. threshold 1e-10 +# Probes: pressure-recovery cancellation p=(E-pi_inf)/gamma +# loses ~4 decimal digits when pi_inf/p_right ~ 40,000 +# +# For each case: 1 nearest-rounding reference run + N random-rounding runs. +# PASS if max L∞ deviation across all N samples stays below threshold. +# +# Verrou (Valgrind 3.26.0 + edf-hpc/verrou@a58d434) is built once and cached. +# Build takes ~20 min uncached; cached runs restore in ~30 s. + +on: + schedule: + - cron: '0 5 * * 1' # Weekly Monday 5 AM UTC + workflow_dispatch: # Manual trigger + +concurrency: + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: true + +jobs: + fp-stability: + name: Floating-Point Stability (Verrou) + runs-on: ubuntu-latest + + steps: + - name: Clone + uses: actions/checkout@v4 + + - name: Cache Verrou + id: cache-verrou + uses: actions/cache@v4 + with: + path: ~/.local/verrou + key: verrou-a58d434-valgrind-3.26.0-${{ runner.os }} + + - name: Install system dependencies + run: | + sudo apt-get update -y + sudo apt-get install -y \ + build-essential automake python3 libc6-dbg \ + cmake gfortran + + - name: Build Verrou + if: steps.cache-verrou.outputs.cache-hit != 'true' + run: | + cd /tmp + wget -q https://sourceware.org/pub/valgrind/valgrind-3.26.0.tar.bz2 + tar xf valgrind-3.26.0.tar.bz2 + + git clone https://github.com/edf-hpc/verrou.git + git -C verrou checkout a58d434 + + # Merge Verrou into Valgrind source tree and patch + cp -r verrou valgrind-3.26.0/verrou + cd valgrind-3.26.0 + cat verrou/valgrind.*diff | patch -p1 + + ./autogen.sh + ./configure --enable-only64bit --prefix="$HOME/.local/verrou" + make -j"$(nproc)" + make install + + - name: Verify Verrou + run: ~/.local/verrou/bin/valgrind --version + + - name: Build MFC (debug, serial) + run: ./mfc.sh build -t pre_process simulation --no-mpi --debug -j"$(nproc)" + + - name: Run FP Stability Suite + run: ./mfc.sh fp-stability -N 5 From b0951f0060ead5b2b9426c49addc5b545c696777 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Tue, 5 May 2026 22:47:35 -0400 Subject: [PATCH 03/13] ci: trigger fp-stability on every push --- .github/workflows/fp-stability.yml | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/.github/workflows/fp-stability.yml b/.github/workflows/fp-stability.yml index 5b98f1f261..57caa40b32 100644 --- a/.github/workflows/fp-stability.yml +++ b/.github/workflows/fp-stability.yml @@ -22,9 +22,8 @@ name: FP Stability # Build takes ~20 min uncached; cached runs restore in ~30 s. on: - schedule: - - cron: '0 5 * * 1' # Weekly Monday 5 AM UTC - workflow_dispatch: # Manual trigger + push: + workflow_dispatch: concurrency: group: ${{ github.workflow }}-${{ github.ref }} From 5c8eaab38ce8eb3f89715a4eb47e4d7ff518e9be Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Tue, 5 May 2026 22:47:51 -0400 Subject: [PATCH 04/13] ci: remove concurrency group from fp-stability --- .github/workflows/fp-stability.yml | 4 ---- 1 file changed, 4 deletions(-) diff --git a/.github/workflows/fp-stability.yml b/.github/workflows/fp-stability.yml index 57caa40b32..1ead1d5cc8 100644 --- a/.github/workflows/fp-stability.yml +++ b/.github/workflows/fp-stability.yml @@ -25,10 +25,6 @@ on: push: workflow_dispatch: -concurrency: - group: ${{ github.workflow }}-${{ github.ref }} - cancel-in-progress: true - jobs: fp-stability: name: Floating-Point Stability (Verrou) From 13a6faf0bfacc0f84a1eb7682ce13ee25dd0fad7 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Wed, 6 May 2026 12:53:16 -0400 Subject: [PATCH 05/13] style: collapse multi-line raise ValueError to satisfy CI ruff formatter --- toolchain/mfc/params/registry.py | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/toolchain/mfc/params/registry.py b/toolchain/mfc/params/registry.py index 1d93eae624..82f895b2f4 100644 --- a/toolchain/mfc/params/registry.py +++ b/toolchain/mfc/params/registry.py @@ -233,11 +233,7 @@ def register(self, param: ParamDef) -> None: if param.name in self._params: existing = self._params[param.name] if existing.param_type != param.param_type: - raise ValueError( - f"Type mismatch for '{param.name}': " - f"existing type is {existing.param_type!r}, " - f"new type is {param.param_type!r}" - ) + raise ValueError(f"Type mismatch for '{param.name}': existing type is {existing.param_type!r}, new type is {param.param_type!r}") existing.tags.update(param.tags) for tag in param.tags: self._by_tag[tag].add(param.name) From 3e0a258f15ae89a0ec31ab962318a3f9c5a53388 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Wed, 6 May 2026 10:19:12 -0400 Subject: [PATCH 06/13] feat: add sod_standard baseline, enlarge cases, run dd_sym on failure MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Add sod_standard (p_L/p_R=10, well-conditioned, threshold=1e-13) as a canary that should always pass under any FP rounding - Enlarge sod_strong and water_stiffened from 25→50 cells and 5→10 steps to give Verrou more arithmetic to perturb and raise sensitivity - Run verrou_dd_sym automatically on any failing case: generates dd_run.sh and dd_cmp.py per case, sets PYTHONPATH for verrou Python libs, saves rddmin_summary and full logs to fp-stability-logs/ (uploaded as CI artifact) - Add python3-numpy to CI apt-get and upload-artifact step (if: always()) --- .github/workflows/fp-stability.yml | 24 ++- .../cases/sod_standard/pre_process.inp | 33 +++ .../cases/sod_standard/simulation.inp | 29 +++ .../cases/sod_strong/pre_process.inp | 2 +- .../cases/sod_strong/simulation.inp | 8 +- .../cases/water_stiffened/pre_process.inp | 2 +- .../cases/water_stiffened/simulation.inp | 8 +- toolchain/mfc/fp_stability.py | 203 +++++++++++++++--- 8 files changed, 270 insertions(+), 39 deletions(-) create mode 100644 tests/fp_stability/cases/sod_standard/pre_process.inp create mode 100644 tests/fp_stability/cases/sod_standard/simulation.inp diff --git a/.github/workflows/fp-stability.yml b/.github/workflows/fp-stability.yml index 1ead1d5cc8..c8f4c07ea9 100644 --- a/.github/workflows/fp-stability.yml +++ b/.github/workflows/fp-stability.yml @@ -3,20 +3,26 @@ name: FP Stability # Runs the Verrou-based floating-point stability suite. # # What it tests (./mfc.sh fp-stability): -# sod_strong 1-D Sod shock, p_L/p_R=100,000, ideal gas +# sod_standard 1-D standard Sod, p_L/p_R=10, ideal gas (well-conditioned baseline) # 25 cells, 5 steps, WENO5 + HLLC -# Compares density & energy vs. threshold 1e-10 +# Threshold 1e-13 — should always pass +# +# sod_strong 1-D Sod shock, p_L/p_R=100,000, ideal gas +# 50 cells, 10 steps, WENO5 + HLLC +# Threshold 1e-10 # Probes: HLLC xi-factor cancellation near sonic contact # (s_L - vel_L)/(s_L - s_S) when s_L ≈ s_S # # water_stiffened 1-D water shock, stiffened EOS (pi_inf=4046) -# 25 cells, 5 steps, WENO5 + HLLC -# Compares density & pressure vs. threshold 1e-10 +# 50 cells, 10 steps, WENO5 + HLLC +# Threshold 1e-10 # Probes: pressure-recovery cancellation p=(E-pi_inf)/gamma # loses ~4 decimal digits when pi_inf/p_right ~ 40,000 # # For each case: 1 nearest-rounding reference run + N random-rounding runs. # PASS if max L∞ deviation across all N samples stays below threshold. +# On FAIL: verrou_dd_sym runs to identify the responsible function symbols. +# Logs are uploaded as CI artifacts. # # Verrou (Valgrind 3.26.0 + edf-hpc/verrou@a58d434) is built once and cached. # Build takes ~20 min uncached; cached runs restore in ~30 s. @@ -45,7 +51,7 @@ jobs: run: | sudo apt-get update -y sudo apt-get install -y \ - build-essential automake python3 libc6-dbg \ + build-essential automake python3 python3-numpy libc6-dbg \ cmake gfortran - name: Build Verrou @@ -76,3 +82,11 @@ jobs: - name: Run FP Stability Suite run: ./mfc.sh fp-stability -N 5 + + - name: Upload FP stability logs + if: always() + uses: actions/upload-artifact@v4 + with: + name: fp-stability-logs + path: fp-stability-logs/ + if-no-files-found: ignore diff --git a/tests/fp_stability/cases/sod_standard/pre_process.inp b/tests/fp_stability/cases/sod_standard/pre_process.inp new file mode 100644 index 0000000000..8beed0a160 --- /dev/null +++ b/tests/fp_stability/cases/sod_standard/pre_process.inp @@ -0,0 +1,33 @@ +&user_inputs +x_domain%beg = 0.0 +x_domain%end = 1.0 +m = 24 +n = 0 +p = 0 +t_step_start = 0 +num_patches = 2 +model_eqns = 2 +num_fluids = 1 +mpp_lim = F +weno_order = 5 +bc_x%beg = -3 +bc_x%end = -3 +precision = 2 +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 = 2.5000000000000004 +fluid_pp(1)%pi_inf = 0.0 +&end/ diff --git a/tests/fp_stability/cases/sod_standard/simulation.inp b/tests/fp_stability/cases/sod_standard/simulation.inp new file mode 100644 index 0000000000..623dc58f1b --- /dev/null +++ b/tests/fp_stability/cases/sod_standard/simulation.inp @@ -0,0 +1,29 @@ +&user_inputs +run_time_info = F +x_domain%beg = 0.0 +x_domain%end = 1.0 +m = 24 +n = 0 +p = 0 +dt = 0.001 +t_step_start = 0 +t_step_stop = 5 +t_step_save = 5 +model_eqns = 2 +num_fluids = 1 +mpp_lim = F +mixture_err = F +time_stepper = 3 +weno_order = 5 +weno_eps = 1e-16 +riemann_solver = 2 +wave_speeds = 1 +avg_state = 2 +bc_x%beg = -3 +bc_x%end = -3 +precision = 2 +prim_vars_wrt = F +parallel_io = F +fluid_pp(1)%gamma = 2.5000000000000004 +fluid_pp(1)%pi_inf = 0.0 +&end/ diff --git a/tests/fp_stability/cases/sod_strong/pre_process.inp b/tests/fp_stability/cases/sod_strong/pre_process.inp index 33693f3d58..88fb152dd5 100644 --- a/tests/fp_stability/cases/sod_strong/pre_process.inp +++ b/tests/fp_stability/cases/sod_strong/pre_process.inp @@ -1,7 +1,7 @@ &user_inputs x_domain%beg = 0.0 x_domain%end = 1.0 -m = 24 +m = 49 n = 0 p = 0 t_step_start = 0 diff --git a/tests/fp_stability/cases/sod_strong/simulation.inp b/tests/fp_stability/cases/sod_strong/simulation.inp index 5e3091e668..ffa377ff33 100644 --- a/tests/fp_stability/cases/sod_strong/simulation.inp +++ b/tests/fp_stability/cases/sod_strong/simulation.inp @@ -2,13 +2,13 @@ run_time_info = F x_domain%beg = 0.0 x_domain%end = 1.0 -m = 24 +m = 49 n = 0 p = 0 -dt = 0.0001 +dt = 0.00005 t_step_start = 0 -t_step_stop = 5 -t_step_save = 5 +t_step_stop = 10 +t_step_save = 10 model_eqns = 2 num_fluids = 1 mpp_lim = F diff --git a/tests/fp_stability/cases/water_stiffened/pre_process.inp b/tests/fp_stability/cases/water_stiffened/pre_process.inp index 986eb0f81f..6fd8be999e 100644 --- a/tests/fp_stability/cases/water_stiffened/pre_process.inp +++ b/tests/fp_stability/cases/water_stiffened/pre_process.inp @@ -1,7 +1,7 @@ &user_inputs x_domain%beg = 0.0 x_domain%end = 1.0 -m = 24 +m = 49 n = 0 p = 0 t_step_start = 0 diff --git a/tests/fp_stability/cases/water_stiffened/simulation.inp b/tests/fp_stability/cases/water_stiffened/simulation.inp index 06a8c554ea..ea9e67f196 100644 --- a/tests/fp_stability/cases/water_stiffened/simulation.inp +++ b/tests/fp_stability/cases/water_stiffened/simulation.inp @@ -2,13 +2,13 @@ run_time_info = F x_domain%beg = 0.0 x_domain%end = 1.0 -m = 24 +m = 49 n = 0 p = 0 -dt = 5e-5 +dt = 2.5e-5 t_step_start = 0 -t_step_stop = 5 -t_step_save = 5 +t_step_stop = 10 +t_step_save = 10 model_eqns = 2 num_fluids = 1 mpp_lim = F diff --git a/toolchain/mfc/fp_stability.py b/toolchain/mfc/fp_stability.py index 3d907a4e81..14dcf0297b 100644 --- a/toolchain/mfc/fp_stability.py +++ b/toolchain/mfc/fp_stability.py @@ -5,6 +5,9 @@ and measures the max deviation from a nearest-rounding reference. Cases that exceed their threshold are reported as FAIL. +On failure, verrou_dd_sym is run to identify the minimal set of functions +responsible for the instability. Logs are saved to fp-stability-logs/. + Requires: - Verrou-enabled Valgrind at $VERROU_HOME/bin/valgrind (default: $HOME/.local/verrou) @@ -19,9 +22,11 @@ import glob import os import shutil +import stat import subprocess import sys import tempfile +import textwrap import time from .common import MFC_ROOT_DIR, MFCException @@ -33,21 +38,28 @@ # Each case: # name - subdirectory under CASES_DIR # description - human-readable purpose -# compare - list of D/ filenames to compare (relative to save step) +# compare - list of D/ filenames to compare # threshold - max L∞ deviation allowed (in conserved-variable units) -# ill_cond - description of the known ill-conditioning +# ill_cond - description of the known ill-conditioning (empty = none expected) CASES = [ + { + "name": "sod_standard", + "description": "1-D standard Sod, p_L/p_R=10, ideal gas (well-conditioned baseline)", + "compare": ["cons.1.00.000005.dat", "cons.3.00.000005.dat"], + "threshold": 1e-13, + "ill_cond": "", + }, { "name": "sod_strong", "description": "1-D Sod, p_L/p_R=100,000, ideal gas", - "compare": ["cons.1.00.000005.dat", "cons.3.00.000005.dat"], + "compare": ["cons.1.00.000010.dat", "cons.3.00.000010.dat"], "threshold": 1e-10, "ill_cond": "HLLC xi factor: (s_L - vel_L)/(s_L - s_S) cancels near sonic contact", }, { "name": "water_stiffened", "description": "1-D water shock, stiffened EOS (pi_inf=4046)", - "compare": ["cons.1.00.000005.dat", "prim.3.00.000005.dat"], + "compare": ["cons.1.00.000010.dat", "prim.3.00.000010.dat"], "threshold": 1e-10, "ill_cond": "Pressure recovery: p = (E - pi_inf)/gamma loses ~4 digits (pi_inf/p_right ~ 40,000)", }, @@ -71,11 +83,16 @@ def _find_binary(name: str) -> str: return max(candidates, key=os.path.getmtime) -def _exclude_args(exclude_file: str) -> list: - """Return --exclude flag if a verrou exclusion file is provided.""" - if exclude_file and os.path.isfile(exclude_file): - return [f"--exclude={exclude_file}"] - return [] +def _find_dd_sym(verrou_bin: str) -> str: + candidate = os.path.join(os.path.dirname(verrou_bin), "verrou_dd_sym") + return candidate if os.path.isfile(candidate) else "" + + +def _verrou_pythonpath(verrou_bin: str) -> str: + """Return the path that must be on PYTHONPATH for verrou_dd_sym's imports.""" + verrou_home = os.path.dirname(os.path.dirname(verrou_bin)) + matches = glob.glob(os.path.join(verrou_home, "lib", "python*", "site-packages", "valgrind")) + return matches[0] if matches else "" def _run_preprocess(pp_bin: str, case_dir: str, work_dir: str): @@ -99,17 +116,14 @@ def _run_simulation_verrou( work_dir: str, run_dir: str, rounding_mode: str, - exclude_args: list, ): """Copy IC files into a fresh tmpdir, run simulation under verrou, collect D/ output.""" with tempfile.TemporaryDirectory(prefix="mfc-fps-") as tmpdir: - # Copy static inputs for fname in ["simulation.inp", "indices.dat", "pre_time_data.dat", "io_time_data.dat"]: src = os.path.join(work_dir, fname) if os.path.exists(src): shutil.copy2(src, tmpdir) - # Copy ICs (p_all_ic → p_all) shutil.copytree( os.path.join(work_dir, "p_all"), os.path.join(tmpdir, "p_all"), @@ -123,7 +137,6 @@ def _run_simulation_verrou( f"--rounding-mode={rounding_mode}", "--error-limit=no", f"--log-file={log_path}", - *exclude_args, sim_bin, ] @@ -133,7 +146,6 @@ def _run_simulation_verrou( if result.returncode != 0: raise MFCException(f"simulation (rounding={rounding_mode}) exited {result.returncode}. See {run_dir}/sim.out") - # Collect D/ outputs os.makedirs(run_dir, exist_ok=True) d_src = os.path.join(tmpdir, "D") for fn in os.listdir(d_src): @@ -160,8 +172,140 @@ def _max_diff(ref_dir: str, run_dir: str, compare_files: list) -> float: return total -def _run_case(case: dict, verrou_bin: str, sim_bin: str, pp_bin: str, n_samples: int) -> bool: - """Run one stability test case. Returns True if the case passes.""" +def _write_dd_run_sh(path: str, verrou_bin: str, sim_bin: str, ic_dir: str): + """Write a dd_run.sh script for verrou_dd_sym. + + verrou_dd_sym calls: dd_run.sh RUNDIR + It sets VERROU_ROUNDING_MODE (reference) or VERROU_EXCLUDE (sample runs) + in the environment. The script must place D/ output files into RUNDIR. + """ + content = textwrap.dedent(f"""\ + #!/usr/bin/env bash + # Generated by mfc.sh fp-stability — do not edit by hand. + VERROU_BIN={verrou_bin!r} + SIM_BIN={sim_bin!r} + IC_DIR={ic_dir!r} + + RUNDIR="$1" + TMPDIR_RUN=$(mktemp -d) + trap 'rm -rf "$TMPDIR_RUN"' EXIT + + cp -r "$IC_DIR/p_all" "$TMPDIR_RUN/p_all" + cp "$IC_DIR/simulation.inp" "$TMPDIR_RUN/simulation.inp" + for fname in indices.dat pre_time_data.dat io_time_data.dat; do + [ -f "$IC_DIR/$fname" ] && cp "$IC_DIR/$fname" "$TMPDIR_RUN/" + done + mkdir -p "$TMPDIR_RUN/D" + + cd "$TMPDIR_RUN" + "$VERROU_BIN" --tool=verrou --error-limit=no "$SIM_BIN" + rc=$? + + [ -d "$TMPDIR_RUN/D" ] && cp -a "$TMPDIR_RUN/D/." "$RUNDIR/" + exit $rc + """) + with open(path, "w") as f: + f.write(content) + os.chmod(path, os.stat(path).st_mode | stat.S_IXUSR | stat.S_IXGRP | stat.S_IXOTH) + + +def _write_dd_cmp_py(path: str, compare_files: list, threshold: float): + """Write a dd_cmp.py script for verrou_dd_sym. + + verrou_dd_sym calls: dd_cmp.py REF_DIR RUN_DIR + Exits 0 if max L∞ deviation is within threshold, 1 otherwise. + """ + files_repr = repr(compare_files) + content = textwrap.dedent(f"""\ + #!/usr/bin/env python3 + # Generated by mfc.sh fp-stability — do not edit by hand. + import sys, os, numpy as np + + COMPARE_FILES = {files_repr} + THRESHOLD = {threshold!r} + + ref_dir, run_dir = sys.argv[1], sys.argv[2] + max_dev = 0.0 + for fname in COMPARE_FILES: + ref_p = os.path.join(ref_dir, fname) + run_p = os.path.join(run_dir, fname) + if not os.path.exists(ref_p) or not os.path.exists(run_p): + print(f"MISSING: {{fname}}") + sys.exit(1) + ref = np.loadtxt(ref_p)[:, 1] + run = np.loadtxt(run_p)[:, 1] + dev = float(np.max(np.abs(ref - run))) + max_dev = max(max_dev, dev) + + print(f"max_dev={{max_dev:.3e}} threshold={{THRESHOLD:.0e}}") + sys.exit(0 if max_dev <= THRESHOLD else 1) + """) + with open(path, "w") as f: + f.write(content) + os.chmod(path, os.stat(path).st_mode | stat.S_IXUSR | stat.S_IXGRP | stat.S_IXOTH) + + +def _run_dd_sym(case: dict, verrou_bin: str, sim_bin: str, work_dir: str, log_dir: str): + """Run verrou_dd_sym to find the minimal set of symbols causing instability.""" + name = case["name"] + dd_sym_bin = _find_dd_sym(verrou_bin) + if not dd_sym_bin: + cons.print(" [dim]verrou_dd_sym not found; skipping delta-debug[/dim]") + return + + dd_dir = os.path.join(log_dir, name) + os.makedirs(dd_dir, exist_ok=True) + + dd_run_sh = os.path.join(dd_dir, "dd_run.sh") + dd_cmp_py = os.path.join(dd_dir, "dd_cmp.py") + _write_dd_run_sh(dd_run_sh, verrou_bin, sim_bin, work_dir) + _write_dd_cmp_py(dd_cmp_py, case["compare"], case["threshold"]) + + # PYTHONPATH: verrou_dd_sym imports dd_config etc. from the valgrind/ subpackage + py_pkg = _verrou_pythonpath(verrou_bin) + env = os.environ.copy() + if py_pkg: + existing = env.get("PYTHONPATH", "") + env["PYTHONPATH"] = ":".join(filter(None, [py_pkg, existing])) + + log_file = os.path.join(dd_dir, "dd_sym.log") + cmd = [ + dd_sym_bin, + "--nruns=10", + "--rddmin=d", + "--reference-rounding=nearest", + dd_run_sh, + dd_cmp_py, + ] + + cons.print(" [dim]running verrou_dd_sym (--nruns=10 --rddmin=d)...[/dim]") + with open(log_file, "w") as f: + result = subprocess.run(cmd, cwd=dd_dir, env=env, stdout=f, stderr=subprocess.STDOUT, check=False) + + if result.returncode == 0: + summary_path = os.path.join(dd_dir, "dd.sym", "rddmin_summary") + if os.path.isfile(summary_path): + cons.print(" [bold yellow]dd_sym result[/bold yellow]:") + with open(summary_path) as f: + for line in f: + cons.print(f" {line.rstrip()}") + else: + cons.print(f" [dim]dd_sym done; see {log_file}[/dim]") + else: + cons.print(f" [bold yellow]dd_sym exited {result.returncode}[/bold yellow] (see {log_file})") + + cons.print(f" [dim]dd_sym logs: {dd_dir}[/dim]") + + +def _run_case( + case: dict, + verrou_bin: str, + sim_bin: str, + pp_bin: str, + n_samples: int, + log_dir: str, +) -> bool: + """Run one stability test case. Returns True if the case passes.""" name = case["name"] threshold = case["threshold"] compare = case["compare"] @@ -169,36 +313,40 @@ def _run_case(case: dict, verrou_bin: str, sim_bin: str, pp_bin: str, n_samples: cons.print(f"[bold]{name}[/bold]: {case['description']}") cons.indent() - cons.print(f" ill-conditioning: {case['ill_cond']}") + if case["ill_cond"]: + cons.print(f" ill-conditioning: {case['ill_cond']}") cons.print(f" threshold: {threshold:.0e}") work_dir = tempfile.mkdtemp(prefix=f"mfc-fps-{name}-") + passed = False try: - # Generate ICs cons.print(" [dim]running pre_process...[/dim]") shutil.copy2(os.path.join(case_dir, "simulation.inp"), work_dir) _run_preprocess(pp_bin, case_dir, work_dir) - # Reference run (nearest rounding — deterministic baseline) ref_dir = os.path.join(work_dir, "ref") os.makedirs(ref_dir) cons.print(" [dim]reference run (rounding=nearest)...[/dim]") - _run_simulation_verrou(verrou_bin, sim_bin, work_dir, ref_dir, "nearest", []) + _run_simulation_verrou(verrou_bin, sim_bin, work_dir, ref_dir, "nearest") - # Random-rounding samples max_dev = 0.0 cons.print(f" [dim]random-rounding runs (N={n_samples})...[/dim]") for i in range(n_samples): run_dir = os.path.join(work_dir, f"run_{i:02d}") os.makedirs(run_dir) - _run_simulation_verrou(verrou_bin, sim_bin, work_dir, run_dir, "random", []) + _run_simulation_verrou(verrou_bin, sim_bin, work_dir, run_dir, "random") dev = _max_diff(ref_dir, run_dir, compare) max_dev = max(max_dev, dev) passed = max_dev <= threshold - tag = "[bold green]PASS[/bold green]" if passed else "[bold red]FAIL[/bold red]" cons.print(f" {tag} max_dev={max_dev:.3e} threshold={threshold:.0e}") + + if not passed: + try: + _run_dd_sym(case, verrou_bin, sim_bin, work_dir, log_dir) + except Exception as exc: + cons.print(f" [bold yellow]dd_sym error[/bold yellow]: {exc}") finally: shutil.rmtree(work_dir, ignore_errors=True) @@ -223,19 +371,23 @@ def fp_stability(): n_samples = ARG("samples") or 5 + log_dir = os.path.join(os.getcwd(), "fp-stability-logs") + os.makedirs(log_dir, exist_ok=True) + cons.print() cons.print("[bold]MFC Floating-Point Stability Suite[/bold]") cons.print(f" verrou: {verrou_bin}") cons.print(f" simulation: {sim_bin}") cons.print(f" pre_process: {pp_bin}") cons.print(f" samples: {n_samples}") + cons.print(f" logs: {log_dir}") cons.print() start = time.time() results = [] for case in CASES: try: - ok = _run_case(case, verrou_bin, sim_bin, pp_bin, n_samples) + ok = _run_case(case, verrou_bin, sim_bin, pp_bin, n_samples, log_dir) except MFCException as exc: cons.print(f" [bold red]ERROR[/bold red]: {exc}") ok = False @@ -250,4 +402,7 @@ def fp_stability(): mark = "[green]✓[/green]" if ok else "[red]✗[/red]" cons.print(f" {mark} {name}") + if n_fail > 0: + cons.print(f"\n dd_sym logs in: {log_dir}") + sys.exit(0 if n_fail == 0 else 1) From d561934751c73ff2b08359d41775b1d7495941bb Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Wed, 6 May 2026 10:42:38 -0400 Subject: [PATCH 07/13] feat: add float proxy, VPREC sweep, dd_line, air_water_interface case - Add four new test cases total: sod_standard, sod_strong, water_stiffened, air_water_interface (two-fluid isobaric contact, stiffened EOS) - Feature B: float-proxy run (--rounding-mode=float) measures single-precision sensitivity without recompiling; skippable with --no-float-proxy - Feature C: VPREC mantissa sweep at [52,23,16,10] bits shows precision floor curve for each case; skippable with --no-vprec - Feature D: verrou_dd_sym symbol-level delta-debug on failure; --no-dd-sym - Feature E: verrou_dd_line source-line delta-debug on failure; --no-dd-line - Add --no-float-proxy/--no-vprec/--no-dd-sym/--no-dd-line CLI flags - Remove dead _max_diff() function (superseded by _max_diff_np()) - Update FP_STABILITY_COMMAND description to document all 4 cases and features --- toolchain/mfc/cli/commands.py | 50 ++++- toolchain/mfc/fp_stability.py | 342 +++++++++++++++++++++++----------- 2 files changed, 276 insertions(+), 116 deletions(-) diff --git a/toolchain/mfc/cli/commands.py b/toolchain/mfc/cli/commands.py index 5f240d620b..c0a0958934 100644 --- a/toolchain/mfc/cli/commands.py +++ b/toolchain/mfc/cli/commands.py @@ -905,14 +905,21 @@ help="Run floating-point stability tests using Verrou.", description=( "Runs each registered test case N times under Verrou's random IEEE-754 " - "rounding mode and compares against a nearest-rounding reference run. " + "rounding mode and compares against a nearest-rounding reference run. " "Reports the max L∞ deviation and PASS/FAIL against per-case thresholds.\n\n" "Requires a Verrou-enabled Valgrind at $VERROU_HOME/bin/valgrind " - "(defaults to $HOME/.local/verrou). The simulation and pre_process " + "(defaults to $HOME/.local/verrou). The simulation and pre_process " "binaries must be serial (no-MPI, no-GPU) debug builds.\n\n" - "Current test cases:\n" - " sod_strong 1-D Sod p_L/p_R=100,000 — HLLC xi-factor cancellation\n" - " water_stiffened 1-D water shock (pi_inf=4046) — pressure-recovery cancellation\n" + "Test cases:\n" + " sod_standard 1-D standard Sod, p_L/p_R=10 (well-conditioned baseline)\n" + " sod_strong 1-D Sod, p_L/p_R=100,000 — HLLC xi-factor cancellation\n" + " water_stiffened 1-D water shock (pi_inf=4046) — pressure-recovery cancellation\n" + " air_water_interface 1-D air/water contact (two-fluid) — mixed-cell cancellation\n\n" + "Additional features (skip with --no-* flags):\n" + " float proxy One run with --rounding-mode=float (single-precision sensitivity)\n" + " vprec sweep Runs at mantissa bits [52, 23, 16, 10] (precision floor curve)\n" + " dd_sym verrou_dd_sym bisection to responsible functions (on failure)\n" + " dd_line verrou_dd_line bisection to responsible source lines (on failure)\n" ), include_common=["mfc_config", "verbose", "debug_log"], arguments=[ @@ -942,6 +949,34 @@ default=5, metavar="N", ), + Argument( + name="no-float-proxy", + help="Skip the --rounding-mode=float single-precision sensitivity run.", + action=ArgAction.STORE_TRUE, + default=False, + dest="no_float_proxy", + ), + Argument( + name="no-vprec", + help="Skip the VPREC mantissa-bit precision sweep.", + action=ArgAction.STORE_TRUE, + default=False, + dest="no_vprec", + ), + Argument( + name="no-dd-sym", + help="Skip verrou_dd_sym function-level delta-debug on failure.", + action=ArgAction.STORE_TRUE, + default=False, + dest="no_dd_sym", + ), + Argument( + name="no-dd-line", + help="Skip verrou_dd_line source-line delta-debug on failure.", + action=ArgAction.STORE_TRUE, + default=False, + dest="no_dd_line", + ), ], examples=[ Example("./mfc.sh fp-stability", "Auto-discover binaries and run all cases"), @@ -950,12 +985,17 @@ "Specify simulation binary explicitly", ), Example("./mfc.sh fp-stability -N 10", "Run 10 random-rounding samples per case"), + Example("./mfc.sh fp-stability --no-vprec --no-dd-line", "Skip VPREC sweep and line debug"), ], key_options=[ ("--sim-binary PATH", "Serial simulation binary (debug, no-MPI)"), ("--pre-binary PATH", "Serial pre_process binary"), ("--verrou-binary PATH", "Verrou-enabled valgrind"), ("-N, --samples N", "Random-rounding samples per case (default: 5)"), + ("--no-float-proxy", "Skip float-rounding proxy run"), + ("--no-vprec", "Skip VPREC mantissa-bit sweep"), + ("--no-dd-sym", "Skip verrou_dd_sym on failure"), + ("--no-dd-line", "Skip verrou_dd_line on failure"), ], ) diff --git a/toolchain/mfc/fp_stability.py b/toolchain/mfc/fp_stability.py index 14dcf0297b..cf7b1da50c 100644 --- a/toolchain/mfc/fp_stability.py +++ b/toolchain/mfc/fp_stability.py @@ -1,12 +1,27 @@ """ Floating-point stability test suite using Verrou. -Each test case runs the simulation N times under random IEEE-754 rounding -and measures the max deviation from a nearest-rounding reference. Cases -that exceed their threshold are reported as FAIL. +Features +-------- +A. Stability suite (always) + N random-rounding samples per case, threshold-based PASS/FAIL. -On failure, verrou_dd_sym is run to identify the minimal set of functions -responsible for the instability. Logs are saved to fp-stability-logs/. +B. Float proxy (--no-float-proxy to skip) + One run with --rounding-mode=float — deterministic proxy for + single-precision sensitivity without recompiling. + +C. VPREC precision sweep (--no-vprec to skip) + One run per mantissa-bit level [52,23,16,10] with + --backend=vprec --vprec-mode=full; shows where each case breaks. + +D. verrou_dd_sym on failure (--no-dd-sym to skip) + Delta-debug bisection isolates the minimal set of *functions* causing + instability. + +E. verrou_dd_line on failure, after dd_sym (--no-dd-line to skip) + Further bisects to exact *source lines* within the responsible functions. + +Logs are saved to fp-stability-logs/ and uploaded as CI artifacts. Requires: - Verrou-enabled Valgrind at $VERROU_HOME/bin/valgrind @@ -15,8 +30,9 @@ - A serial pre_process binary (to generate initial conditions) Usage: + ./mfc.sh fp-stability + ./mfc.sh fp-stability --no-vprec --no-dd-line ./mfc.sh fp-stability --sim-binary PATH --pre-binary PATH - ./mfc.sh fp-stability # auto-discovers newest installed binaries """ import glob @@ -35,12 +51,16 @@ CASES_DIR = os.path.join(MFC_ROOT_DIR, "tests", "fp_stability", "cases") +# Mantissa-bit levels for the VPREC sweep (C). +# 52 = full double, 23 = single, 16 = half-ish, 10 = ultra-low. +VPREC_MANTISSA_BITS = [52, 23, 16, 10] + # Each case: # name - subdirectory under CASES_DIR # description - human-readable purpose # compare - list of D/ filenames to compare -# threshold - max L∞ deviation allowed (in conserved-variable units) -# ill_cond - description of the known ill-conditioning (empty = none expected) +# threshold - max L∞ deviation allowed (conserved-variable units) +# ill_cond - known ill-conditioning (empty string = none expected) CASES = [ { "name": "sod_standard", @@ -61,7 +81,14 @@ "description": "1-D water shock, stiffened EOS (pi_inf=4046)", "compare": ["cons.1.00.000010.dat", "prim.3.00.000010.dat"], "threshold": 1e-10, - "ill_cond": "Pressure recovery: p = (E - pi_inf)/gamma loses ~4 digits (pi_inf/p_right ~ 40,000)", + "ill_cond": "Pressure recovery: p=(E-pi_inf)/gamma loses ~4 digits (pi_inf/p_right~40,000)", + }, + { + "name": "air_water_interface", + "description": "1-D air/water isobaric contact (two-fluid, pi_inf=4046)", + "compare": ["cons.1.00.000005.dat", "cons.4.00.000005.dat", "cons.5.00.000005.dat"], + "threshold": 1e-10, + "ill_cond": "Mixed-cell pressure recovery: E-alpha_w*gamma_w*pi_inf cancels when alpha_w<<1", }, ] @@ -71,41 +98,36 @@ def _find_verrou() -> str: candidate = os.path.join(verrou_home, "bin", "valgrind") if os.path.isfile(candidate) and os.access(candidate, os.X_OK): return candidate - fallback = shutil.which("valgrind") - return fallback or "" + return shutil.which("valgrind") or "" def _find_binary(name: str) -> str: install_dir = os.path.join(os.getcwd(), "build", "install") candidates = glob.glob(os.path.join(install_dir, "*", "bin", name)) - if not candidates: - return "" - return max(candidates, key=os.path.getmtime) + return max(candidates, key=os.path.getmtime) if candidates else "" def _find_dd_sym(verrou_bin: str) -> str: - candidate = os.path.join(os.path.dirname(verrou_bin), "verrou_dd_sym") - return candidate if os.path.isfile(candidate) else "" + c = os.path.join(os.path.dirname(verrou_bin), "verrou_dd_sym") + return c if os.path.isfile(c) else "" + + +def _find_dd_line(verrou_bin: str) -> str: + c = os.path.join(os.path.dirname(verrou_bin), "verrou_dd_line") + return c if os.path.isfile(c) else "" def _verrou_pythonpath(verrou_bin: str) -> str: - """Return the path that must be on PYTHONPATH for verrou_dd_sym's imports.""" + """Path that must be on PYTHONPATH for verrou_dd_* imports (valgrind/ subdir).""" verrou_home = os.path.dirname(os.path.dirname(verrou_bin)) matches = glob.glob(os.path.join(verrou_home, "lib", "python*", "site-packages", "valgrind")) return matches[0] if matches else "" def _run_preprocess(pp_bin: str, case_dir: str, work_dir: str): - """Generate initial conditions in work_dir.""" shutil.copy2(os.path.join(case_dir, "pre_process.inp"), work_dir) with open(os.path.join(work_dir, "pre.log"), "w") as f: - result = subprocess.run( - [pp_bin], - cwd=work_dir, - stdout=f, - stderr=subprocess.STDOUT, - check=False, - ) + result = subprocess.run([pp_bin], cwd=work_dir, stdout=f, stderr=subprocess.STDOUT, check=False) if result.returncode != 0: raise MFCException(f"pre_process failed (rc={result.returncode}). See {work_dir}/pre.log") @@ -115,69 +137,90 @@ def _run_simulation_verrou( sim_bin: str, work_dir: str, run_dir: str, - rounding_mode: str, + rounding_mode: str = None, + extra_flags: list = None, ): - """Copy IC files into a fresh tmpdir, run simulation under verrou, collect D/ output.""" + """Copy ICs into a fresh tmpdir, run simulation under verrou, collect D/ output. + + rounding_mode is passed as --rounding-mode= when not None. + extra_flags are appended before the binary (e.g. --backend=vprec ...). + """ with tempfile.TemporaryDirectory(prefix="mfc-fps-") as tmpdir: for fname in ["simulation.inp", "indices.dat", "pre_time_data.dat", "io_time_data.dat"]: src = os.path.join(work_dir, fname) if os.path.exists(src): shutil.copy2(src, tmpdir) - - shutil.copytree( - os.path.join(work_dir, "p_all"), - os.path.join(tmpdir, "p_all"), - ) + shutil.copytree(os.path.join(work_dir, "p_all"), os.path.join(tmpdir, "p_all")) os.makedirs(os.path.join(tmpdir, "D")) log_path = os.path.join(run_dir, "verrou.log") - cmd = [ - verrou_bin, - "--tool=verrou", - f"--rounding-mode={rounding_mode}", - "--error-limit=no", - f"--log-file={log_path}", - sim_bin, - ] + cmd = [verrou_bin, "--tool=verrou", "--error-limit=no", f"--log-file={log_path}"] + if rounding_mode: + cmd.append(f"--rounding-mode={rounding_mode}") + cmd.extend(extra_flags or []) + cmd.append(sim_bin) with open(os.path.join(run_dir, "sim.out"), "w") as f: result = subprocess.run(cmd, cwd=tmpdir, stdout=f, stderr=subprocess.STDOUT, check=False) if result.returncode != 0: - raise MFCException(f"simulation (rounding={rounding_mode}) exited {result.returncode}. See {run_dir}/sim.out") + tag = rounding_mode or "vprec" + raise MFCException(f"simulation ({tag}) exited {result.returncode}. See {run_dir}/sim.out") os.makedirs(run_dir, exist_ok=True) - d_src = os.path.join(tmpdir, "D") - for fn in os.listdir(d_src): - shutil.copy2(os.path.join(d_src, fn), run_dir) + for fn in os.listdir(os.path.join(tmpdir, "D")): + shutil.copy2(os.path.join(tmpdir, "D", fn), run_dir) -def _max_diff(ref_dir: str, run_dir: str, compare_files: list) -> float: - try: - import numpy as np - except ImportError: - raise MFCException("numpy is required for fp-stability comparison") +def _max_diff_np(ref_dir: str, run_dir: str, compare_files: list) -> float: + import numpy as np total = 0.0 for fname in compare_files: - ref_path = os.path.join(ref_dir, fname) - run_path = os.path.join(run_dir, fname) - if not os.path.exists(ref_path): - raise MFCException(f"Reference output missing: {ref_path}") - if not os.path.exists(run_path): - raise MFCException(f"Run output missing: {run_path}") - ref = np.loadtxt(ref_path)[:, 1] - run = np.loadtxt(run_path)[:, 1] + ref_p, run_p = os.path.join(ref_dir, fname), os.path.join(run_dir, fname) + if not os.path.exists(ref_p) or not os.path.exists(run_p): + return float("inf") + ref = np.loadtxt(ref_p)[:, 1] + run = np.loadtxt(run_p)[:, 1] total = max(total, float(np.max(np.abs(ref - run)))) return total +def _run_float_proxy(case: dict, verrou_bin: str, sim_bin: str, work_dir: str, ref_dir: str) -> float: + """One run with --rounding-mode=float; returns L∞ deviation from nearest-ref.""" + run_dir = os.path.join(work_dir, "float_proxy") + os.makedirs(run_dir) + _run_simulation_verrou(verrou_bin, sim_bin, work_dir, run_dir, rounding_mode="float") + return _max_diff_np(ref_dir, run_dir, case["compare"]) + + +def _run_vprec_sweep(case: dict, verrou_bin: str, sim_bin: str, work_dir: str, ref_dir: str) -> list: + """Run at each mantissa-bit level. Returns [(bits, dev), ...].""" + results = [] + for bits in VPREC_MANTISSA_BITS: + run_dir = os.path.join(work_dir, f"vprec_{bits}") + os.makedirs(run_dir) + flags = [ + "--backend=vprec", + "--vprec-mode=full", + f"--vprec-precision-binary64={bits}", + "--vprec-range-binary64=11", + ] + try: + _run_simulation_verrou(verrou_bin, sim_bin, work_dir, run_dir, extra_flags=flags) + dev = _max_diff_np(ref_dir, run_dir, case["compare"]) + except MFCException: + dev = float("inf") + results.append((bits, dev)) + return results + + def _write_dd_run_sh(path: str, verrou_bin: str, sim_bin: str, ic_dir: str): - """Write a dd_run.sh script for verrou_dd_sym. + """Generate dd_run.sh for verrou_dd_sym / verrou_dd_line. - verrou_dd_sym calls: dd_run.sh RUNDIR - It sets VERROU_ROUNDING_MODE (reference) or VERROU_EXCLUDE (sample runs) - in the environment. The script must place D/ output files into RUNDIR. + verrou_dd_* calls: dd_run.sh RUNDIR + It sets VERROU_ROUNDING_MODE / VERROU_EXCLUDE / VERROU_SOURCE etc. in the + environment; the script just invokes verrou and copies D/ output to RUNDIR. """ content = textwrap.dedent(f"""\ #!/usr/bin/env bash @@ -210,18 +253,17 @@ def _write_dd_run_sh(path: str, verrou_bin: str, sim_bin: str, ic_dir: str): def _write_dd_cmp_py(path: str, compare_files: list, threshold: float): - """Write a dd_cmp.py script for verrou_dd_sym. + """Generate dd_cmp.py for verrou_dd_sym / verrou_dd_line. - verrou_dd_sym calls: dd_cmp.py REF_DIR RUN_DIR - Exits 0 if max L∞ deviation is within threshold, 1 otherwise. + verrou_dd_* calls: dd_cmp.py REF_DIR RUN_DIR + Exits 0 (stable) or 1 (unstable) based on threshold. """ - files_repr = repr(compare_files) content = textwrap.dedent(f"""\ #!/usr/bin/env python3 # Generated by mfc.sh fp-stability — do not edit by hand. import sys, os, numpy as np - COMPARE_FILES = {files_repr} + COMPARE_FILES = {compare_files!r} THRESHOLD = {threshold!r} ref_dir, run_dir = sys.argv[1], sys.argv[2] @@ -245,58 +287,78 @@ def _write_dd_cmp_py(path: str, compare_files: list, threshold: float): os.chmod(path, os.stat(path).st_mode | stat.S_IXUSR | stat.S_IXGRP | stat.S_IXOTH) -def _run_dd_sym(case: dict, verrou_bin: str, sim_bin: str, work_dir: str, log_dir: str): - """Run verrou_dd_sym to find the minimal set of symbols causing instability.""" - name = case["name"] - dd_sym_bin = _find_dd_sym(verrou_bin) - if not dd_sym_bin: - cons.print(" [dim]verrou_dd_sym not found; skipping delta-debug[/dim]") - return - - dd_dir = os.path.join(log_dir, name) - os.makedirs(dd_dir, exist_ok=True) - - dd_run_sh = os.path.join(dd_dir, "dd_run.sh") - dd_cmp_py = os.path.join(dd_dir, "dd_cmp.py") - _write_dd_run_sh(dd_run_sh, verrou_bin, sim_bin, work_dir) - _write_dd_cmp_py(dd_cmp_py, case["compare"], case["threshold"]) - - # PYTHONPATH: verrou_dd_sym imports dd_config etc. from the valgrind/ subpackage +def _dd_env(verrou_bin: str) -> dict: + """Environment with PYTHONPATH set for verrou_dd_* imports.""" py_pkg = _verrou_pythonpath(verrou_bin) env = os.environ.copy() if py_pkg: existing = env.get("PYTHONPATH", "") env["PYTHONPATH"] = ":".join(filter(None, [py_pkg, existing])) - - log_file = os.path.join(dd_dir, "dd_sym.log") - cmd = [ - dd_sym_bin, - "--nruns=10", - "--rddmin=d", - "--reference-rounding=nearest", - dd_run_sh, - dd_cmp_py, - ] - - cons.print(" [dim]running verrou_dd_sym (--nruns=10 --rddmin=d)...[/dim]") + return env + + +def _run_dd_tool( + dd_bin: str, + dd_dir: str, + dd_run_sh: str, + dd_cmp_py: str, + env: dict, + log_name: str, + summary_subdir: str, + label: str, +): + """Generic runner for verrou_dd_sym / verrou_dd_line.""" + log_file = os.path.join(dd_dir, log_name) + cmd = [dd_bin, "--nruns=10", "--rddmin=d", "--reference-rounding=nearest", dd_run_sh, dd_cmp_py] + cons.print(f" [dim]running {label} (--nruns=10 --rddmin=d)...[/dim]") with open(log_file, "w") as f: result = subprocess.run(cmd, cwd=dd_dir, env=env, stdout=f, stderr=subprocess.STDOUT, check=False) - if result.returncode == 0: - summary_path = os.path.join(dd_dir, "dd.sym", "rddmin_summary") + summary_path = os.path.join(dd_dir, summary_subdir, "rddmin_summary") if os.path.isfile(summary_path): - cons.print(" [bold yellow]dd_sym result[/bold yellow]:") + cons.print(f" [bold yellow]{label} result[/bold yellow]:") with open(summary_path) as f: for line in f: cons.print(f" {line.rstrip()}") else: - cons.print(f" [dim]dd_sym done; see {log_file}[/dim]") + cons.print(f" [dim]{label} done; see {log_file}[/dim]") else: - cons.print(f" [bold yellow]dd_sym exited {result.returncode}[/bold yellow] (see {log_file})") + cons.print(f" [bold yellow]{label} exited {result.returncode}[/bold yellow] (see {log_file})") + +def _run_dd_sym(case: dict, verrou_bin: str, sim_bin: str, work_dir: str, log_dir: str): + dd_bin = _find_dd_sym(verrou_bin) + if not dd_bin: + cons.print(" [dim]verrou_dd_sym not found; skipping delta-debug[/dim]") + return + + dd_dir = os.path.join(log_dir, case["name"]) + os.makedirs(dd_dir, exist_ok=True) + dd_run_sh = os.path.join(dd_dir, "dd_run.sh") + dd_cmp_py = os.path.join(dd_dir, "dd_cmp.py") + _write_dd_run_sh(dd_run_sh, verrou_bin, sim_bin, work_dir) + _write_dd_cmp_py(dd_cmp_py, case["compare"], case["threshold"]) + _run_dd_tool(dd_bin, dd_dir, dd_run_sh, dd_cmp_py, _dd_env(verrou_bin), "dd_sym.log", "dd.sym", "verrou_dd_sym") cons.print(f" [dim]dd_sym logs: {dd_dir}[/dim]") +def _run_dd_line(case: dict, verrou_bin: str, sim_bin: str, work_dir: str, log_dir: str): + dd_bin = _find_dd_line(verrou_bin) + if not dd_bin: + cons.print(" [dim]verrou_dd_line not found; skipping line-level debug[/dim]") + return + + # Reuse scripts written by dd_sym (or write them now if dd_sym was skipped). + dd_dir = os.path.join(log_dir, case["name"]) + os.makedirs(dd_dir, exist_ok=True) + dd_run_sh = os.path.join(dd_dir, "dd_run.sh") + dd_cmp_py = os.path.join(dd_dir, "dd_cmp.py") + if not os.path.isfile(dd_run_sh): + _write_dd_run_sh(dd_run_sh, verrou_bin, sim_bin, work_dir) + _write_dd_cmp_py(dd_cmp_py, case["compare"], case["threshold"]) + _run_dd_tool(dd_bin, dd_dir, dd_run_sh, dd_cmp_py, _dd_env(verrou_bin), "dd_line.log", "dd.line", "verrou_dd_line") + + def _run_case( case: dict, verrou_bin: str, @@ -304,8 +366,11 @@ def _run_case( pp_bin: str, n_samples: int, log_dir: str, + run_float: bool, + run_vprec: bool, + run_dd_sym: bool, + run_dd_line: bool, ) -> bool: - """Run one stability test case. Returns True if the case passes.""" name = case["name"] threshold = case["threshold"] compare = case["compare"] @@ -327,26 +392,56 @@ def _run_case( ref_dir = os.path.join(work_dir, "ref") os.makedirs(ref_dir) cons.print(" [dim]reference run (rounding=nearest)...[/dim]") - _run_simulation_verrou(verrou_bin, sim_bin, work_dir, ref_dir, "nearest") + _run_simulation_verrou(verrou_bin, sim_bin, work_dir, ref_dir, rounding_mode="nearest") + # --- A: random-rounding stability samples --- max_dev = 0.0 cons.print(f" [dim]random-rounding runs (N={n_samples})...[/dim]") for i in range(n_samples): run_dir = os.path.join(work_dir, f"run_{i:02d}") os.makedirs(run_dir) - _run_simulation_verrou(verrou_bin, sim_bin, work_dir, run_dir, "random") - dev = _max_diff(ref_dir, run_dir, compare) - max_dev = max(max_dev, dev) + _run_simulation_verrou(verrou_bin, sim_bin, work_dir, run_dir, rounding_mode="random") + max_dev = max(max_dev, _max_diff_np(ref_dir, run_dir, compare)) passed = max_dev <= threshold tag = "[bold green]PASS[/bold green]" if passed else "[bold red]FAIL[/bold red]" cons.print(f" {tag} max_dev={max_dev:.3e} threshold={threshold:.0e}") - if not passed: + # --- B: float proxy --- + if run_float: try: - _run_dd_sym(case, verrou_bin, sim_bin, work_dir, log_dir) - except Exception as exc: - cons.print(f" [bold yellow]dd_sym error[/bold yellow]: {exc}") + fdev = _run_float_proxy(case, verrou_bin, sim_bin, work_dir, ref_dir) + cons.print(f" float proxy: dev={fdev:.3e} (single-precision sensitivity)") + except MFCException as exc: + cons.print(f" [dim]float proxy error: {exc}[/dim]") + + # --- C: VPREC sweep --- + if run_vprec: + cons.print(" VPREC precision sweep:") + vprec_results = _run_vprec_sweep(case, verrou_bin, sim_bin, work_dir, ref_dir) + labels = {52: "double", 23: "single", 16: "~half", 10: "ultra-low"} + for bits, dev in vprec_results: + label = labels.get(bits, "") + label_str = f" ({label})" if label else "" + marker = "" + if dev == float("inf"): + marker = " [red]crashed[/red]" + elif dev > threshold: + marker = " [red]FAIL[/red]" + cons.print(f" {bits:2d} bits{label_str}: dev={dev:.3e}{marker}") + + # --- D/E: delta-debug on failure --- + if not passed: + if run_dd_sym: + try: + _run_dd_sym(case, verrou_bin, sim_bin, work_dir, log_dir) + except Exception as exc: + cons.print(f" [bold yellow]dd_sym error[/bold yellow]: {exc}") + if run_dd_line: + try: + _run_dd_line(case, verrou_bin, sim_bin, work_dir, log_dir) + except Exception as exc: + cons.print(f" [bold yellow]dd_line error[/bold yellow]: {exc}") finally: shutil.rmtree(work_dir, ignore_errors=True) @@ -370,6 +465,10 @@ def fp_stability(): raise MFCException("pre_process binary not found. Build with --no-mpi, or pass --pre-binary.") n_samples = ARG("samples") or 5 + run_float = not ARG("no_float_proxy") + run_vprec = not ARG("no_vprec") + run_dd_sym = not ARG("no_dd_sym") + run_dd_line = not ARG("no_dd_line") log_dir = os.path.join(os.getcwd(), "fp-stability-logs") os.makedirs(log_dir, exist_ok=True) @@ -380,6 +479,16 @@ def fp_stability(): cons.print(f" simulation: {sim_bin}") cons.print(f" pre_process: {pp_bin}") cons.print(f" samples: {n_samples}") + features = [] + if run_float: + features.append("float-proxy") + if run_vprec: + features.append("vprec-sweep") + if run_dd_sym: + features.append("dd_sym") + if run_dd_line: + features.append("dd_line") + cons.print(f" features: {', '.join(features) if features else 'stability only'}") cons.print(f" logs: {log_dir}") cons.print() @@ -387,7 +496,18 @@ def fp_stability(): results = [] for case in CASES: try: - ok = _run_case(case, verrou_bin, sim_bin, pp_bin, n_samples, log_dir) + ok = _run_case( + case, + verrou_bin, + sim_bin, + pp_bin, + n_samples, + log_dir, + run_float, + run_vprec, + run_dd_sym, + run_dd_line, + ) except MFCException as exc: cons.print(f" [bold red]ERROR[/bold red]: {exc}") ok = False @@ -403,6 +523,6 @@ def fp_stability(): cons.print(f" {mark} {name}") if n_fail > 0: - cons.print(f"\n dd_sym logs in: {log_dir}") + cons.print(f"\n dd_sym/dd_line logs in: {log_dir}") sys.exit(0 if n_fail == 0 else 1) From f932c4a2a765ffd6383de8d939be57239a292c99 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Wed, 6 May 2026 12:28:00 -0400 Subject: [PATCH 08/13] fix: add missing air_water_interface case input files --- .../cases/air_water_interface/pre_process.inp | 37 +++++++++++++++++++ .../cases/air_water_interface/simulation.inp | 31 ++++++++++++++++ 2 files changed, 68 insertions(+) create mode 100644 tests/fp_stability/cases/air_water_interface/pre_process.inp create mode 100644 tests/fp_stability/cases/air_water_interface/simulation.inp diff --git a/tests/fp_stability/cases/air_water_interface/pre_process.inp b/tests/fp_stability/cases/air_water_interface/pre_process.inp new file mode 100644 index 0000000000..01d9b01a6f --- /dev/null +++ b/tests/fp_stability/cases/air_water_interface/pre_process.inp @@ -0,0 +1,37 @@ +&user_inputs +x_domain%beg = 0.0 +x_domain%end = 1.0 +m = 24 +n = 0 +p = 0 +t_step_start = 0 +num_patches = 2 +model_eqns = 2 +num_fluids = 2 +mpp_lim = T +weno_order = 5 +bc_x%beg = -3 +bc_x%end = -3 +precision = 2 +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_rho(2) = 0.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 = 1.0 +patch_icpp(2)%alpha_rho(1) = 0.0 +patch_icpp(2)%alpha_rho(2) = 1.0 +patch_icpp(2)%alpha(1) = 0.0 +fluid_pp(1)%gamma = 2.5000000000000004 +fluid_pp(1)%pi_inf = 0.0 +fluid_pp(2)%gamma = 0.195312500 +fluid_pp(2)%pi_inf = 4046.31 +&end/ diff --git a/tests/fp_stability/cases/air_water_interface/simulation.inp b/tests/fp_stability/cases/air_water_interface/simulation.inp new file mode 100644 index 0000000000..1253b168f6 --- /dev/null +++ b/tests/fp_stability/cases/air_water_interface/simulation.inp @@ -0,0 +1,31 @@ +&user_inputs +run_time_info = F +x_domain%beg = 0.0 +x_domain%end = 1.0 +m = 24 +n = 0 +p = 0 +dt = 5e-5 +t_step_start = 0 +t_step_stop = 5 +t_step_save = 5 +model_eqns = 2 +num_fluids = 2 +mpp_lim = T +mixture_err = T +time_stepper = 3 +weno_order = 5 +weno_eps = 1e-16 +riemann_solver = 2 +wave_speeds = 1 +avg_state = 2 +bc_x%beg = -3 +bc_x%end = -3 +precision = 2 +prim_vars_wrt = F +parallel_io = F +fluid_pp(1)%gamma = 2.5000000000000004 +fluid_pp(1)%pi_inf = 0.0 +fluid_pp(2)%gamma = 0.195312500 +fluid_pp(2)%pi_inf = 4046.31 +&end/ From d363ca27fe624048a25e80b41e50932994c4d812 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Wed, 6 May 2026 12:41:41 -0400 Subject: [PATCH 09/13] fix: add alpha(2) for both patches in air_water_interface pre_process.inp --- tests/fp_stability/cases/air_water_interface/pre_process.inp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/tests/fp_stability/cases/air_water_interface/pre_process.inp b/tests/fp_stability/cases/air_water_interface/pre_process.inp index 01d9b01a6f..663be3d2e3 100644 --- a/tests/fp_stability/cases/air_water_interface/pre_process.inp +++ b/tests/fp_stability/cases/air_water_interface/pre_process.inp @@ -22,6 +22,7 @@ patch_icpp(1)%pres = 1.0 patch_icpp(1)%alpha_rho(1) = 1.0 patch_icpp(1)%alpha_rho(2) = 0.0 patch_icpp(1)%alpha(1) = 1.0 +patch_icpp(1)%alpha(2) = 0.0 patch_icpp(2)%geometry = 1 patch_icpp(2)%x_centroid = 0.75 patch_icpp(2)%length_x = 0.5 @@ -30,6 +31,7 @@ patch_icpp(2)%pres = 1.0 patch_icpp(2)%alpha_rho(1) = 0.0 patch_icpp(2)%alpha_rho(2) = 1.0 patch_icpp(2)%alpha(1) = 0.0 +patch_icpp(2)%alpha(2) = 1.0 fluid_pp(1)%gamma = 2.5000000000000004 fluid_pp(1)%pi_inf = 0.0 fluid_pp(2)%gamma = 0.195312500 From 0e7f365393ac97ee769d968df7410a04400c61ed Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Wed, 6 May 2026 14:52:22 -0400 Subject: [PATCH 10/13] fix: add reduced-energy (E-tilde) scheme for model_eqns=2 FP stability MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Store Ẽ = E - pi_inf_mix (reduced energy) in the conserved energy slot for the Allaire 5-equation model (model_eqns=2). This eliminates catastrophic cancellation when recovering pressure via p = (Ẽ - KE - qv)/gamma vs the old p = (E - pi_inf - KE - qv)/gamma where pi_inf >> p (e.g. water: pi_inf=4046 >> p=1). Key changes: - m_variables_conversion.fpp: scope Ẽ storage to model_eqns=2 only; add explicit model_eqns=1/3 branch (physical E) to avoid fallthrough to Tait EOS for model_eqns=3 (Saurel 6-equation) - m_riemann_solvers.fpp: HLL and LF non-MHD blocks use Ẽ as states with pi_inf added back for enthalpy H; HLLC Block 1 (model_eqns=3) and Block 4 (model_eqns=2) retain their existing correct handling; xi-factor denominators protected with sgm_eps (Fix 1) - m_rhs.fpp: add non-conservative Ẽ source S_Ẽ = -sum_i(pi_inf_i * rhs(alpha_i)) for x/y/z directions, gated on model_eqns==2 and .not.(bubbles_euler .or. mhd) - m_sim_helpers.fpp: igr branch drops pi_inf from pressure recovery (Ẽ stored) and restores it for physical H - m_cbc.fpp: energy flux drops dpi_inf_dt term (Ẽ not pi_inf in flux) --- src/common/m_variables_conversion.fpp | 11 ++- src/simulation/m_cbc.fpp | 7 +- src/simulation/m_rhs.fpp | 118 +++++++++++++++++++++++++- src/simulation/m_riemann_solvers.fpp | 48 ++++++----- src/simulation/m_sim_helpers.fpp | 4 +- 5 files changed, 158 insertions(+), 30 deletions(-) diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index 2417da1adf..c9c683f3b0 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -104,8 +104,12 @@ contains if (mhd) then ! MHD pressure: subtract magnetic pressure from total energy pres = (energy - dyn_p - pi_inf - qv - pres_mag)/gamma + else if (model_eqns == 2 .and. (bubbles_euler .neqv. .true.)) then + ! Allaire et al. (JCP 2002): store reduced energy Ẽ = E - pi_inf_mix so that p = (Ẽ - KE - qv)/gamma avoids + ! cancellation when pi_inf_mix >> p. + pres = (energy - dyn_p - qv)/gamma else if ((model_eqns /= 4) .and. (bubbles_euler .neqv. .true.)) then - ! Gamma/pi_inf model or five-equation model (Allaire et al. JCP 2002): p from mixture EOS + ! model_eqns=1 or 3: physical E stored, p = (E - pi_inf - KE - qv)/gamma pres = (energy - dyn_p - pi_inf - qv)/gamma else if ((model_eqns /= 4) .and. bubbles_euler) then ! Bubble-augmented pressure with void fraction correction @@ -919,8 +923,11 @@ contains ! MHD energy includes magnetic pressure contribution q_cons_vf(eqn_idx%E)%sf(j, k, l) = gamma*q_prim_vf(eqn_idx%E)%sf(j, k, & & l) + dyn_pres + pres_mag + pi_inf + qv + else if (model_eqns == 2 .and. (bubbles_euler .neqv. .true.)) then + ! Store reduced energy Ẽ = gamma*p + KE + qv (no pi_inf) for numerical stability. + q_cons_vf(eqn_idx%E)%sf(j, k, l) = gamma*q_prim_vf(eqn_idx%E)%sf(j, k, l) + dyn_pres + qv else if ((model_eqns /= 4) .and. (bubbles_euler .neqv. .true.)) then - ! Five-equation model (Allaire et al. JCP 2002): E = Gamma*p + 0.5*rho*|u|^2 + pi_inf + qv + ! model_eqns=1 or 3: store physical E = gamma*p + pi_inf + KE + qv q_cons_vf(eqn_idx%E)%sf(j, k, l) = gamma*q_prim_vf(eqn_idx%E)%sf(j, k, l) + dyn_pres + pi_inf + qv else if ((model_eqns /= 4) .and. (bubbles_euler)) then ! Bubble-augmented energy with void fraction correction diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp index d703a3c1e3..ba9b262b66 100644 --- a/src/simulation/m_cbc.fpp +++ b/src/simulation/m_cbc.fpp @@ -658,10 +658,11 @@ contains gamma = 1.0_wp/(Cp/Cv - 1.0_wp) end if else - E = gamma*pres + pi_inf + 5.e-1_wp*rho*vel_K_sum + ! Reduced energy Ẽ = E - pi_inf; add pi_inf back for physical H. + E = gamma*pres + 5.e-1_wp*rho*vel_K_sum end if - H = (E + pres)/rho + H = (E + pi_inf + pres)/rho ! Compute mixture sound speed call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv_local, vel_K_sum, 0._wp, c, qv) @@ -890,7 +891,7 @@ contains end do else flux_rs${XYZ}$_vf_l(-1, k, r, eqn_idx%E) = flux_rs${XYZ}$_vf_l(0, k, r, & - & eqn_idx%E) + ds(0)*(pres*dgamma_dt + gamma*dpres_dt + dpi_inf_dt + dqv_dt & + & eqn_idx%E) + ds(0)*(pres*dgamma_dt + gamma*dpres_dt + dqv_dt & & + rho*vel_dv_dt_sum + 5.e-1_wp*drho_dt*vel_K_sum) end if diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index 8589a7ff52..a8194ff103 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -1118,9 +1118,9 @@ contains type(vector_field), intent(in) :: flux_src_n_vf_arg ! CORRECTED DECLARATION FOR Kterm_arg: real(wp), allocatable, dimension(:,:,:), intent(in) :: Kterm_arg - integer :: j_adv, k_idx, l_idx, q_idx + integer :: j_adv, k_idx, l_idx, q_idx, i_fluid real(wp) :: local_inv_ds, local_term_coeff, local_flux1, local_flux2 - real(wp) :: local_q_cons_val, local_k_term_val + real(wp) :: local_q_cons_val, local_k_term_val, pi_inf_src logical :: use_standard_riemann select case (current_idir) @@ -1191,6 +1191,44 @@ contains $:END_GPU_PARALLEL_LOOP() end if end if + ! Ẽ non-conservative source for model_eqns=2: S_Ẽ = -\Sigma_i \pi_inf_i * rhs(\alpha_i) + if (model_eqns == 2 .and. (bubbles_euler .neqv. .true.) .and. .not. mhd) then + if (use_standard_riemann) then + $:GPU_PARALLEL_LOOP(collapse=3, private='[k_idx, l_idx, q_idx, local_inv_ds, local_term_coeff, & + & pi_inf_src, i_fluid]') + do q_idx = 0, p; do l_idx = 0, n; do k_idx = 0, m + local_inv_ds = 1._wp/dx(k_idx) + local_term_coeff = q_prim_vf_arg%vf(eqn_idx%cont%end + 1)%sf(k_idx, l_idx, q_idx) + pi_inf_src = 0._wp + do i_fluid = 1, num_fluids + pi_inf_src = pi_inf_src + pi_infs(i_fluid)*(flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & + & %sf(k_idx - 1, l_idx, & + & q_idx) - flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & + & %sf(k_idx, l_idx, q_idx)) + end do + rhs_vf_arg(eqn_idx%E)%sf(k_idx, l_idx, q_idx) = rhs_vf_arg(eqn_idx%E)%sf(k_idx, l_idx, & + & q_idx) - local_inv_ds*local_term_coeff*pi_inf_src + end do; end do; end do + $:END_GPU_PARALLEL_LOOP() + else + $:GPU_PARALLEL_LOOP(collapse=3, private='[k_idx, l_idx, q_idx, local_inv_ds, pi_inf_src, i_fluid]') + do q_idx = 0, p; do l_idx = 0, n; do k_idx = 0, m + local_inv_ds = 1._wp/dx(k_idx) + pi_inf_src = 0._wp + do i_fluid = 1, num_fluids + pi_inf_src = pi_inf_src + pi_infs(i_fluid)*q_cons_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & + & %sf(k_idx, l_idx, & + & q_idx)*(flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & + & %sf(k_idx, l_idx, & + & q_idx) - flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & + & %sf(k_idx - 1, l_idx, q_idx)) + end do + rhs_vf_arg(eqn_idx%E)%sf(k_idx, l_idx, q_idx) = rhs_vf_arg(eqn_idx%E)%sf(k_idx, l_idx, & + & q_idx) - local_inv_ds*pi_inf_src + end do; end do; end do + $:END_GPU_PARALLEL_LOOP() + end if + end if case (2) ! y-direction: loops q_idx (x), k_idx (y), l_idx (z); sf(q_idx, k_idx, l_idx); dy(k_idx); Kterm(q_idx,k_idx,l_idx) use_standard_riemann = (riemann_solver == 1 .or. riemann_solver == 4) @@ -1267,6 +1305,44 @@ contains $:END_GPU_PARALLEL_LOOP() end if end if + ! Ẽ non-conservative source for model_eqns=2: S_Ẽ = -\Sigma_i \pi_inf_i * rhs(\alpha_i) + if (model_eqns == 2 .and. (bubbles_euler .neqv. .true.) .and. .not. mhd) then + if (use_standard_riemann) then + $:GPU_PARALLEL_LOOP(collapse=3, private='[k_idx, l_idx, q_idx, local_inv_ds, local_term_coeff, & + & pi_inf_src, i_fluid]') + do l_idx = 0, p; do k_idx = 0, n; do q_idx = 0, m + local_inv_ds = 1._wp/dy(k_idx) + local_term_coeff = q_prim_vf_arg%vf(eqn_idx%cont%end + 2)%sf(q_idx, k_idx, l_idx) + pi_inf_src = 0._wp + do i_fluid = 1, num_fluids + pi_inf_src = pi_inf_src + pi_infs(i_fluid)*(flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & + & %sf(q_idx, k_idx - 1, & + & l_idx) - flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & + & %sf(q_idx, k_idx, l_idx)) + end do + rhs_vf_arg(eqn_idx%E)%sf(q_idx, k_idx, l_idx) = rhs_vf_arg(eqn_idx%E)%sf(q_idx, k_idx, & + & l_idx) - local_inv_ds*local_term_coeff*pi_inf_src + end do; end do; end do + $:END_GPU_PARALLEL_LOOP() + else + $:GPU_PARALLEL_LOOP(collapse=3, private='[k_idx, l_idx, q_idx, local_inv_ds, pi_inf_src, i_fluid]') + do l_idx = 0, p; do k_idx = 0, n; do q_idx = 0, m + local_inv_ds = 1._wp/dy(k_idx) + pi_inf_src = 0._wp + do i_fluid = 1, num_fluids + pi_inf_src = pi_inf_src + pi_infs(i_fluid)*q_cons_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & + & %sf(q_idx, k_idx, & + & l_idx)*(flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & + & %sf(q_idx, k_idx, & + & l_idx) - flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & + & %sf(q_idx, k_idx - 1, l_idx)) + end do + rhs_vf_arg(eqn_idx%E)%sf(q_idx, k_idx, l_idx) = rhs_vf_arg(eqn_idx%E)%sf(q_idx, k_idx, & + & l_idx) - local_inv_ds*pi_inf_src + end do; end do; end do + $:END_GPU_PARALLEL_LOOP() + end if + end if case (3) ! z-direction: loops l_idx (x), q_idx (y), k_idx (z); sf(l_idx, q_idx, k_idx); dz(k_idx); Kterm(l_idx,q_idx,k_idx) if (grid_geometry == 3) then @@ -1340,6 +1416,44 @@ contains $:END_GPU_PARALLEL_LOOP() end if end if + ! Ẽ non-conservative source for model_eqns=2: S_Ẽ = -\Sigma_i \pi_inf_i * rhs(\alpha_i) + if (model_eqns == 2 .and. (bubbles_euler .neqv. .true.) .and. .not. mhd) then + if (use_standard_riemann) then + $:GPU_PARALLEL_LOOP(collapse=3, private='[k_idx, l_idx, q_idx, local_inv_ds, local_term_coeff, & + & pi_inf_src, i_fluid]') + do k_idx = 0, p; do q_idx = 0, n; do l_idx = 0, m + local_inv_ds = 1._wp/dz(k_idx) + local_term_coeff = q_prim_vf_arg%vf(eqn_idx%cont%end + 3)%sf(l_idx, q_idx, k_idx) + pi_inf_src = 0._wp + do i_fluid = 1, num_fluids + pi_inf_src = pi_inf_src + pi_infs(i_fluid)*(flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & + & %sf(l_idx, q_idx, & + & k_idx - 1) - flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid & + & - 1)%sf(l_idx, q_idx, k_idx)) + end do + rhs_vf_arg(eqn_idx%E)%sf(l_idx, q_idx, k_idx) = rhs_vf_arg(eqn_idx%E)%sf(l_idx, q_idx, & + & k_idx) - local_inv_ds*local_term_coeff*pi_inf_src + end do; end do; end do + $:END_GPU_PARALLEL_LOOP() + else + $:GPU_PARALLEL_LOOP(collapse=3, private='[k_idx, l_idx, q_idx, local_inv_ds, pi_inf_src, i_fluid]') + do k_idx = 0, p; do q_idx = 0, n; do l_idx = 0, m + local_inv_ds = 1._wp/dz(k_idx) + pi_inf_src = 0._wp + do i_fluid = 1, num_fluids + pi_inf_src = pi_inf_src + pi_infs(i_fluid)*q_cons_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & + & %sf(l_idx, q_idx, & + & k_idx)*(flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & + & %sf(l_idx, q_idx, & + & k_idx) - flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & + & %sf(l_idx, q_idx, k_idx - 1)) + end do + rhs_vf_arg(eqn_idx%E)%sf(l_idx, q_idx, k_idx) = rhs_vf_arg(eqn_idx%E)%sf(l_idx, q_idx, & + & k_idx) - local_inv_ds*pi_inf_src + end do; end do; end do + $:END_GPU_PARALLEL_LOOP() + end if + end if end select end subroutine s_add_directional_advection_source_terms diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 33661b54d7..69b2fa8c72 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -429,10 +429,11 @@ contains H_R = (E_R + pres_R - pres_mag%R) & & /rho_R ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound) else - E_L = gamma_L*pres_L + pi_inf_L + 5.e-1*rho_L*vel_L_rms + qv_L - E_R = gamma_R*pres_R + pi_inf_R + 5.e-1*rho_R*vel_R_rms + qv_R - H_L = (E_L + pres_L)/rho_L - H_R = (E_R + pres_R)/rho_R + ! Reduced energy Ẽ = E - pi_inf; add pi_inf back for physical H. + E_L = gamma_L*pres_L + 5.e-1*rho_L*vel_L_rms + qv_L + E_R = gamma_R*pres_R + 5.e-1*rho_R*vel_R_rms + qv_R + H_L = (E_L + pi_inf_L + pres_L)/rho_L + H_R = (E_R + pi_inf_R + pres_R)/rho_R end if ! elastic energy update @@ -1119,10 +1120,11 @@ contains H_R = (E_R + pres_R - pres_mag%R) & & /rho_R ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound) else - E_L = gamma_L*pres_L + pi_inf_L + 5.e-1*rho_L*vel_L_rms + qv_L - E_R = gamma_R*pres_R + pi_inf_R + 5.e-1*rho_R*vel_R_rms + qv_R - H_L = (E_L + pres_L)/rho_L - H_R = (E_R + pres_R)/rho_R + ! Reduced energy Ẽ = E - pi_inf; add pi_inf back for physical H. + E_L = gamma_L*pres_L + 5.e-1*rho_L*vel_L_rms + qv_L + E_R = gamma_R*pres_R + 5.e-1*rho_R*vel_R_rms + qv_R + H_L = (E_L + pi_inf_L + pres_L)/rho_L + H_R = (E_R + pi_inf_R + pres_R)/rho_R end if ! elastic energy update @@ -1966,6 +1968,7 @@ contains end do end if + ! E_L/E_R are physical E (include pi_inf) for model_eqns=3 H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R @@ -2040,8 +2043,8 @@ contains s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) ! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) - xi_L = (s_L - vel_L(dir_idx(1)))/(s_L - s_S) - xi_R = (s_R - vel_R(dir_idx(1)))/(s_R - s_S) + xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps) + xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps) ! goes with numerical star velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star) xi_M = (5.e-1_wp + sign(0.5_wp, s_S)) @@ -2331,8 +2334,8 @@ contains s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) ! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) - xi_L = (s_L - vel_L(dir_idx(1)))/(s_L - s_S) - xi_R = (s_R - vel_R(dir_idx(1)))/(s_R - s_S) + xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps) + xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps) ! goes with numerical velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star) xi_M = (5.e-1_wp + sign(5.e-1_wp, s_S)) @@ -2688,8 +2691,8 @@ contains s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) ! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) - xi_L = (s_L - vel_L(dir_idx(1)))/(s_L - s_S) - xi_R = (s_R - vel_R(dir_idx(1)))/(s_R - s_S) + xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps) + xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps) ! goes with numerical velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star) xi_M = (5.e-1_wp + sign(5.e-1_wp, s_S)) @@ -2988,11 +2991,12 @@ contains H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R else - E_L = gamma_L*pres_L + pi_inf_L + 5.e-1*rho_L*vel_L_rms + qv_L - E_R = gamma_R*pres_R + pi_inf_R + 5.e-1*rho_R*vel_R_rms + qv_R + ! Reduced energy Ẽ = E - pi_inf_mix; add pi_inf back for physical enthalpy H. + E_L = gamma_L*pres_L + 5.e-1*rho_L*vel_L_rms + qv_L + E_R = gamma_R*pres_R + 5.e-1*rho_R*vel_R_rms + qv_R - H_L = (E_L + pres_L)/rho_L - H_R = (E_R + pres_R)/rho_R + H_L = (E_L + pi_inf_L + pres_L)/rho_L + H_R = (E_R + pi_inf_R + pres_R)/rho_R end if ! ENERGY ADJUSTMENTS FOR HYPOELASTIC ENERGY @@ -3051,8 +3055,8 @@ contains end do end if - H_L = (E_L + pres_L)/rho_L - H_R = (E_R + pres_R)/rho_R + H_L = (E_L + pi_inf_L + pres_L)/rho_L + H_R = (E_R + pi_inf_R + pres_R)/rho_R @:compute_average_state() @@ -3127,8 +3131,8 @@ contains s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) ! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) - xi_L = (s_L - vel_L(dir_idx(1)))/(s_L - s_S) - xi_R = (s_R - vel_R(dir_idx(1)))/(s_R - s_S) + xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps) + xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps) ! goes with numerical velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star) xi_M = (5.e-1_wp + sign(5.e-1_wp, s_S)) diff --git a/src/simulation/m_sim_helpers.fpp b/src/simulation/m_sim_helpers.fpp index 4a0978919e..abff6f32bf 100644 --- a/src/simulation/m_sim_helpers.fpp +++ b/src/simulation/m_sim_helpers.fpp @@ -121,8 +121,10 @@ contains end do if (igr) then + ! igr passes q_cons_vf which stores Ẽ = E - pi_inf_mix; drop pi_inf from recovery. E = q_prim_vf(eqn_idx%E)%sf(j, k, l) - pres = (E - pi_inf - qv - 5.e-1_wp*rho*vel_sum)/gamma + pres = (E - qv - 5.e-1_wp*rho*vel_sum)/gamma + E = E + pi_inf ! Ẽ → physical E for H below else pres = q_prim_vf(eqn_idx%E)%sf(j, k, l) E = gamma*pres + pi_inf + 5.e-1_wp*rho*vel_sum + qv From 7e80db0d575dd864faba8dfade8b73513d22c8dc Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Wed, 6 May 2026 15:08:33 -0400 Subject: [PATCH 11/13] ci: increase fp-stability case step counts to 50 for better sweep coverage --- .../fp_stability/cases/air_water_interface/simulation.inp | 4 ++-- tests/fp_stability/cases/sod_standard/simulation.inp | 4 ++-- tests/fp_stability/cases/sod_strong/simulation.inp | 4 ++-- tests/fp_stability/cases/water_stiffened/simulation.inp | 4 ++-- toolchain/mfc/fp_stability.py | 8 ++++---- 5 files changed, 12 insertions(+), 12 deletions(-) diff --git a/tests/fp_stability/cases/air_water_interface/simulation.inp b/tests/fp_stability/cases/air_water_interface/simulation.inp index 1253b168f6..09fbb69e2b 100644 --- a/tests/fp_stability/cases/air_water_interface/simulation.inp +++ b/tests/fp_stability/cases/air_water_interface/simulation.inp @@ -7,8 +7,8 @@ n = 0 p = 0 dt = 5e-5 t_step_start = 0 -t_step_stop = 5 -t_step_save = 5 +t_step_stop = 50 +t_step_save = 50 model_eqns = 2 num_fluids = 2 mpp_lim = T diff --git a/tests/fp_stability/cases/sod_standard/simulation.inp b/tests/fp_stability/cases/sod_standard/simulation.inp index 623dc58f1b..c75f9a2542 100644 --- a/tests/fp_stability/cases/sod_standard/simulation.inp +++ b/tests/fp_stability/cases/sod_standard/simulation.inp @@ -7,8 +7,8 @@ n = 0 p = 0 dt = 0.001 t_step_start = 0 -t_step_stop = 5 -t_step_save = 5 +t_step_stop = 50 +t_step_save = 50 model_eqns = 2 num_fluids = 1 mpp_lim = F diff --git a/tests/fp_stability/cases/sod_strong/simulation.inp b/tests/fp_stability/cases/sod_strong/simulation.inp index ffa377ff33..5c20727aef 100644 --- a/tests/fp_stability/cases/sod_strong/simulation.inp +++ b/tests/fp_stability/cases/sod_strong/simulation.inp @@ -7,8 +7,8 @@ n = 0 p = 0 dt = 0.00005 t_step_start = 0 -t_step_stop = 10 -t_step_save = 10 +t_step_stop = 50 +t_step_save = 50 model_eqns = 2 num_fluids = 1 mpp_lim = F diff --git a/tests/fp_stability/cases/water_stiffened/simulation.inp b/tests/fp_stability/cases/water_stiffened/simulation.inp index ea9e67f196..24b9bb5d06 100644 --- a/tests/fp_stability/cases/water_stiffened/simulation.inp +++ b/tests/fp_stability/cases/water_stiffened/simulation.inp @@ -7,8 +7,8 @@ n = 0 p = 0 dt = 2.5e-5 t_step_start = 0 -t_step_stop = 10 -t_step_save = 10 +t_step_stop = 50 +t_step_save = 50 model_eqns = 2 num_fluids = 1 mpp_lim = F diff --git a/toolchain/mfc/fp_stability.py b/toolchain/mfc/fp_stability.py index cf7b1da50c..23877669fb 100644 --- a/toolchain/mfc/fp_stability.py +++ b/toolchain/mfc/fp_stability.py @@ -65,28 +65,28 @@ { "name": "sod_standard", "description": "1-D standard Sod, p_L/p_R=10, ideal gas (well-conditioned baseline)", - "compare": ["cons.1.00.000005.dat", "cons.3.00.000005.dat"], + "compare": ["cons.1.00.000050.dat", "cons.3.00.000050.dat"], "threshold": 1e-13, "ill_cond": "", }, { "name": "sod_strong", "description": "1-D Sod, p_L/p_R=100,000, ideal gas", - "compare": ["cons.1.00.000010.dat", "cons.3.00.000010.dat"], + "compare": ["cons.1.00.000050.dat", "cons.3.00.000050.dat"], "threshold": 1e-10, "ill_cond": "HLLC xi factor: (s_L - vel_L)/(s_L - s_S) cancels near sonic contact", }, { "name": "water_stiffened", "description": "1-D water shock, stiffened EOS (pi_inf=4046)", - "compare": ["cons.1.00.000010.dat", "prim.3.00.000010.dat"], + "compare": ["cons.1.00.000050.dat", "prim.3.00.000050.dat"], "threshold": 1e-10, "ill_cond": "Pressure recovery: p=(E-pi_inf)/gamma loses ~4 digits (pi_inf/p_right~40,000)", }, { "name": "air_water_interface", "description": "1-D air/water isobaric contact (two-fluid, pi_inf=4046)", - "compare": ["cons.1.00.000005.dat", "cons.4.00.000005.dat", "cons.5.00.000005.dat"], + "compare": ["cons.1.00.000050.dat", "cons.4.00.000050.dat", "cons.5.00.000050.dat"], "threshold": 1e-10, "ill_cond": "Mixed-cell pressure recovery: E-alpha_w*gamma_w*pi_inf cancels when alpha_w<<1", }, From ab8ae6c519b939b2f994f1edd0934d0506e96fcb Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Wed, 6 May 2026 16:08:33 -0400 Subject: [PATCH 12/13] =?UTF-8?q?revert:=20remove=20reduced-energy=20(E-ti?= =?UTF-8?q?lde)=20scheme=20=E2=80=94=20moving=20to=20separate=20PR?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/common/m_variables_conversion.fpp | 11 +-- src/simulation/m_cbc.fpp | 7 +- src/simulation/m_rhs.fpp | 118 +------------------------- src/simulation/m_riemann_solvers.fpp | 48 +++++------ src/simulation/m_sim_helpers.fpp | 4 +- 5 files changed, 30 insertions(+), 158 deletions(-) diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index c9c683f3b0..2417da1adf 100644 --- a/src/common/m_variables_conversion.fpp +++ b/src/common/m_variables_conversion.fpp @@ -104,12 +104,8 @@ contains if (mhd) then ! MHD pressure: subtract magnetic pressure from total energy pres = (energy - dyn_p - pi_inf - qv - pres_mag)/gamma - else if (model_eqns == 2 .and. (bubbles_euler .neqv. .true.)) then - ! Allaire et al. (JCP 2002): store reduced energy Ẽ = E - pi_inf_mix so that p = (Ẽ - KE - qv)/gamma avoids - ! cancellation when pi_inf_mix >> p. - pres = (energy - dyn_p - qv)/gamma else if ((model_eqns /= 4) .and. (bubbles_euler .neqv. .true.)) then - ! model_eqns=1 or 3: physical E stored, p = (E - pi_inf - KE - qv)/gamma + ! Gamma/pi_inf model or five-equation model (Allaire et al. JCP 2002): p from mixture EOS pres = (energy - dyn_p - pi_inf - qv)/gamma else if ((model_eqns /= 4) .and. bubbles_euler) then ! Bubble-augmented pressure with void fraction correction @@ -923,11 +919,8 @@ contains ! MHD energy includes magnetic pressure contribution q_cons_vf(eqn_idx%E)%sf(j, k, l) = gamma*q_prim_vf(eqn_idx%E)%sf(j, k, & & l) + dyn_pres + pres_mag + pi_inf + qv - else if (model_eqns == 2 .and. (bubbles_euler .neqv. .true.)) then - ! Store reduced energy Ẽ = gamma*p + KE + qv (no pi_inf) for numerical stability. - q_cons_vf(eqn_idx%E)%sf(j, k, l) = gamma*q_prim_vf(eqn_idx%E)%sf(j, k, l) + dyn_pres + qv else if ((model_eqns /= 4) .and. (bubbles_euler .neqv. .true.)) then - ! model_eqns=1 or 3: store physical E = gamma*p + pi_inf + KE + qv + ! Five-equation model (Allaire et al. JCP 2002): E = Gamma*p + 0.5*rho*|u|^2 + pi_inf + qv q_cons_vf(eqn_idx%E)%sf(j, k, l) = gamma*q_prim_vf(eqn_idx%E)%sf(j, k, l) + dyn_pres + pi_inf + qv else if ((model_eqns /= 4) .and. (bubbles_euler)) then ! Bubble-augmented energy with void fraction correction diff --git a/src/simulation/m_cbc.fpp b/src/simulation/m_cbc.fpp index ba9b262b66..d703a3c1e3 100644 --- a/src/simulation/m_cbc.fpp +++ b/src/simulation/m_cbc.fpp @@ -658,11 +658,10 @@ contains gamma = 1.0_wp/(Cp/Cv - 1.0_wp) end if else - ! Reduced energy Ẽ = E - pi_inf; add pi_inf back for physical H. - E = gamma*pres + 5.e-1_wp*rho*vel_K_sum + E = gamma*pres + pi_inf + 5.e-1_wp*rho*vel_K_sum end if - H = (E + pi_inf + pres)/rho + H = (E + pres)/rho ! Compute mixture sound speed call s_compute_speed_of_sound(pres, rho, gamma, pi_inf, H, adv_local, vel_K_sum, 0._wp, c, qv) @@ -891,7 +890,7 @@ contains end do else flux_rs${XYZ}$_vf_l(-1, k, r, eqn_idx%E) = flux_rs${XYZ}$_vf_l(0, k, r, & - & eqn_idx%E) + ds(0)*(pres*dgamma_dt + gamma*dpres_dt + dqv_dt & + & eqn_idx%E) + ds(0)*(pres*dgamma_dt + gamma*dpres_dt + dpi_inf_dt + dqv_dt & & + rho*vel_dv_dt_sum + 5.e-1_wp*drho_dt*vel_K_sum) end if diff --git a/src/simulation/m_rhs.fpp b/src/simulation/m_rhs.fpp index a8194ff103..8589a7ff52 100644 --- a/src/simulation/m_rhs.fpp +++ b/src/simulation/m_rhs.fpp @@ -1118,9 +1118,9 @@ contains type(vector_field), intent(in) :: flux_src_n_vf_arg ! CORRECTED DECLARATION FOR Kterm_arg: real(wp), allocatable, dimension(:,:,:), intent(in) :: Kterm_arg - integer :: j_adv, k_idx, l_idx, q_idx, i_fluid + integer :: j_adv, k_idx, l_idx, q_idx real(wp) :: local_inv_ds, local_term_coeff, local_flux1, local_flux2 - real(wp) :: local_q_cons_val, local_k_term_val, pi_inf_src + real(wp) :: local_q_cons_val, local_k_term_val logical :: use_standard_riemann select case (current_idir) @@ -1191,44 +1191,6 @@ contains $:END_GPU_PARALLEL_LOOP() end if end if - ! Ẽ non-conservative source for model_eqns=2: S_Ẽ = -\Sigma_i \pi_inf_i * rhs(\alpha_i) - if (model_eqns == 2 .and. (bubbles_euler .neqv. .true.) .and. .not. mhd) then - if (use_standard_riemann) then - $:GPU_PARALLEL_LOOP(collapse=3, private='[k_idx, l_idx, q_idx, local_inv_ds, local_term_coeff, & - & pi_inf_src, i_fluid]') - do q_idx = 0, p; do l_idx = 0, n; do k_idx = 0, m - local_inv_ds = 1._wp/dx(k_idx) - local_term_coeff = q_prim_vf_arg%vf(eqn_idx%cont%end + 1)%sf(k_idx, l_idx, q_idx) - pi_inf_src = 0._wp - do i_fluid = 1, num_fluids - pi_inf_src = pi_inf_src + pi_infs(i_fluid)*(flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & - & %sf(k_idx - 1, l_idx, & - & q_idx) - flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & - & %sf(k_idx, l_idx, q_idx)) - end do - rhs_vf_arg(eqn_idx%E)%sf(k_idx, l_idx, q_idx) = rhs_vf_arg(eqn_idx%E)%sf(k_idx, l_idx, & - & q_idx) - local_inv_ds*local_term_coeff*pi_inf_src - end do; end do; end do - $:END_GPU_PARALLEL_LOOP() - else - $:GPU_PARALLEL_LOOP(collapse=3, private='[k_idx, l_idx, q_idx, local_inv_ds, pi_inf_src, i_fluid]') - do q_idx = 0, p; do l_idx = 0, n; do k_idx = 0, m - local_inv_ds = 1._wp/dx(k_idx) - pi_inf_src = 0._wp - do i_fluid = 1, num_fluids - pi_inf_src = pi_inf_src + pi_infs(i_fluid)*q_cons_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & - & %sf(k_idx, l_idx, & - & q_idx)*(flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & - & %sf(k_idx, l_idx, & - & q_idx) - flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & - & %sf(k_idx - 1, l_idx, q_idx)) - end do - rhs_vf_arg(eqn_idx%E)%sf(k_idx, l_idx, q_idx) = rhs_vf_arg(eqn_idx%E)%sf(k_idx, l_idx, & - & q_idx) - local_inv_ds*pi_inf_src - end do; end do; end do - $:END_GPU_PARALLEL_LOOP() - end if - end if case (2) ! y-direction: loops q_idx (x), k_idx (y), l_idx (z); sf(q_idx, k_idx, l_idx); dy(k_idx); Kterm(q_idx,k_idx,l_idx) use_standard_riemann = (riemann_solver == 1 .or. riemann_solver == 4) @@ -1305,44 +1267,6 @@ contains $:END_GPU_PARALLEL_LOOP() end if end if - ! Ẽ non-conservative source for model_eqns=2: S_Ẽ = -\Sigma_i \pi_inf_i * rhs(\alpha_i) - if (model_eqns == 2 .and. (bubbles_euler .neqv. .true.) .and. .not. mhd) then - if (use_standard_riemann) then - $:GPU_PARALLEL_LOOP(collapse=3, private='[k_idx, l_idx, q_idx, local_inv_ds, local_term_coeff, & - & pi_inf_src, i_fluid]') - do l_idx = 0, p; do k_idx = 0, n; do q_idx = 0, m - local_inv_ds = 1._wp/dy(k_idx) - local_term_coeff = q_prim_vf_arg%vf(eqn_idx%cont%end + 2)%sf(q_idx, k_idx, l_idx) - pi_inf_src = 0._wp - do i_fluid = 1, num_fluids - pi_inf_src = pi_inf_src + pi_infs(i_fluid)*(flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & - & %sf(q_idx, k_idx - 1, & - & l_idx) - flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & - & %sf(q_idx, k_idx, l_idx)) - end do - rhs_vf_arg(eqn_idx%E)%sf(q_idx, k_idx, l_idx) = rhs_vf_arg(eqn_idx%E)%sf(q_idx, k_idx, & - & l_idx) - local_inv_ds*local_term_coeff*pi_inf_src - end do; end do; end do - $:END_GPU_PARALLEL_LOOP() - else - $:GPU_PARALLEL_LOOP(collapse=3, private='[k_idx, l_idx, q_idx, local_inv_ds, pi_inf_src, i_fluid]') - do l_idx = 0, p; do k_idx = 0, n; do q_idx = 0, m - local_inv_ds = 1._wp/dy(k_idx) - pi_inf_src = 0._wp - do i_fluid = 1, num_fluids - pi_inf_src = pi_inf_src + pi_infs(i_fluid)*q_cons_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & - & %sf(q_idx, k_idx, & - & l_idx)*(flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & - & %sf(q_idx, k_idx, & - & l_idx) - flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & - & %sf(q_idx, k_idx - 1, l_idx)) - end do - rhs_vf_arg(eqn_idx%E)%sf(q_idx, k_idx, l_idx) = rhs_vf_arg(eqn_idx%E)%sf(q_idx, k_idx, & - & l_idx) - local_inv_ds*pi_inf_src - end do; end do; end do - $:END_GPU_PARALLEL_LOOP() - end if - end if case (3) ! z-direction: loops l_idx (x), q_idx (y), k_idx (z); sf(l_idx, q_idx, k_idx); dz(k_idx); Kterm(l_idx,q_idx,k_idx) if (grid_geometry == 3) then @@ -1416,44 +1340,6 @@ contains $:END_GPU_PARALLEL_LOOP() end if end if - ! Ẽ non-conservative source for model_eqns=2: S_Ẽ = -\Sigma_i \pi_inf_i * rhs(\alpha_i) - if (model_eqns == 2 .and. (bubbles_euler .neqv. .true.) .and. .not. mhd) then - if (use_standard_riemann) then - $:GPU_PARALLEL_LOOP(collapse=3, private='[k_idx, l_idx, q_idx, local_inv_ds, local_term_coeff, & - & pi_inf_src, i_fluid]') - do k_idx = 0, p; do q_idx = 0, n; do l_idx = 0, m - local_inv_ds = 1._wp/dz(k_idx) - local_term_coeff = q_prim_vf_arg%vf(eqn_idx%cont%end + 3)%sf(l_idx, q_idx, k_idx) - pi_inf_src = 0._wp - do i_fluid = 1, num_fluids - pi_inf_src = pi_inf_src + pi_infs(i_fluid)*(flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & - & %sf(l_idx, q_idx, & - & k_idx - 1) - flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid & - & - 1)%sf(l_idx, q_idx, k_idx)) - end do - rhs_vf_arg(eqn_idx%E)%sf(l_idx, q_idx, k_idx) = rhs_vf_arg(eqn_idx%E)%sf(l_idx, q_idx, & - & k_idx) - local_inv_ds*local_term_coeff*pi_inf_src - end do; end do; end do - $:END_GPU_PARALLEL_LOOP() - else - $:GPU_PARALLEL_LOOP(collapse=3, private='[k_idx, l_idx, q_idx, local_inv_ds, pi_inf_src, i_fluid]') - do k_idx = 0, p; do q_idx = 0, n; do l_idx = 0, m - local_inv_ds = 1._wp/dz(k_idx) - pi_inf_src = 0._wp - do i_fluid = 1, num_fluids - pi_inf_src = pi_inf_src + pi_infs(i_fluid)*q_cons_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & - & %sf(l_idx, q_idx, & - & k_idx)*(flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & - & %sf(l_idx, q_idx, & - & k_idx) - flux_src_n_vf_arg%vf(eqn_idx%adv%beg + i_fluid - 1) & - & %sf(l_idx, q_idx, k_idx - 1)) - end do - rhs_vf_arg(eqn_idx%E)%sf(l_idx, q_idx, k_idx) = rhs_vf_arg(eqn_idx%E)%sf(l_idx, q_idx, & - & k_idx) - local_inv_ds*pi_inf_src - end do; end do; end do - $:END_GPU_PARALLEL_LOOP() - end if - end if end select end subroutine s_add_directional_advection_source_terms diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 69b2fa8c72..33661b54d7 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -429,11 +429,10 @@ contains H_R = (E_R + pres_R - pres_mag%R) & & /rho_R ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound) else - ! Reduced energy Ẽ = E - pi_inf; add pi_inf back for physical H. - E_L = gamma_L*pres_L + 5.e-1*rho_L*vel_L_rms + qv_L - E_R = gamma_R*pres_R + 5.e-1*rho_R*vel_R_rms + qv_R - H_L = (E_L + pi_inf_L + pres_L)/rho_L - H_R = (E_R + pi_inf_R + pres_R)/rho_R + E_L = gamma_L*pres_L + pi_inf_L + 5.e-1*rho_L*vel_L_rms + qv_L + E_R = gamma_R*pres_R + pi_inf_R + 5.e-1*rho_R*vel_R_rms + qv_R + H_L = (E_L + pres_L)/rho_L + H_R = (E_R + pres_R)/rho_R end if ! elastic energy update @@ -1120,11 +1119,10 @@ contains H_R = (E_R + pres_R - pres_mag%R) & & /rho_R ! stagnation enthalpy here excludes magnetic energy (only used to find speed of sound) else - ! Reduced energy Ẽ = E - pi_inf; add pi_inf back for physical H. - E_L = gamma_L*pres_L + 5.e-1*rho_L*vel_L_rms + qv_L - E_R = gamma_R*pres_R + 5.e-1*rho_R*vel_R_rms + qv_R - H_L = (E_L + pi_inf_L + pres_L)/rho_L - H_R = (E_R + pi_inf_R + pres_R)/rho_R + E_L = gamma_L*pres_L + pi_inf_L + 5.e-1*rho_L*vel_L_rms + qv_L + E_R = gamma_R*pres_R + pi_inf_R + 5.e-1*rho_R*vel_R_rms + qv_R + H_L = (E_L + pres_L)/rho_L + H_R = (E_R + pres_R)/rho_R end if ! elastic energy update @@ -1968,7 +1966,6 @@ contains end do end if - ! E_L/E_R are physical E (include pi_inf) for model_eqns=3 H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R @@ -2043,8 +2040,8 @@ contains s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) ! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) - xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps) - xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps) + xi_L = (s_L - vel_L(dir_idx(1)))/(s_L - s_S) + xi_R = (s_R - vel_R(dir_idx(1)))/(s_R - s_S) ! goes with numerical star velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star) xi_M = (5.e-1_wp + sign(0.5_wp, s_S)) @@ -2334,8 +2331,8 @@ contains s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) ! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) - xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps) - xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps) + xi_L = (s_L - vel_L(dir_idx(1)))/(s_L - s_S) + xi_R = (s_R - vel_R(dir_idx(1)))/(s_R - s_S) ! goes with numerical velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star) xi_M = (5.e-1_wp + sign(5.e-1_wp, s_S)) @@ -2691,8 +2688,8 @@ contains s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) ! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) - xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps) - xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps) + xi_L = (s_L - vel_L(dir_idx(1)))/(s_L - s_S) + xi_R = (s_R - vel_R(dir_idx(1)))/(s_R - s_S) ! goes with numerical velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star) xi_M = (5.e-1_wp + sign(5.e-1_wp, s_S)) @@ -2991,12 +2988,11 @@ contains H_L = (E_L + pres_L)/rho_L H_R = (E_R + pres_R)/rho_R else - ! Reduced energy Ẽ = E - pi_inf_mix; add pi_inf back for physical enthalpy H. - E_L = gamma_L*pres_L + 5.e-1*rho_L*vel_L_rms + qv_L - E_R = gamma_R*pres_R + 5.e-1*rho_R*vel_R_rms + qv_R + E_L = gamma_L*pres_L + pi_inf_L + 5.e-1*rho_L*vel_L_rms + qv_L + E_R = gamma_R*pres_R + pi_inf_R + 5.e-1*rho_R*vel_R_rms + qv_R - H_L = (E_L + pi_inf_L + pres_L)/rho_L - H_R = (E_R + pi_inf_R + pres_R)/rho_R + H_L = (E_L + pres_L)/rho_L + H_R = (E_R + pres_R)/rho_R end if ! ENERGY ADJUSTMENTS FOR HYPOELASTIC ENERGY @@ -3055,8 +3051,8 @@ contains end do end if - H_L = (E_L + pi_inf_L + pres_L)/rho_L - H_R = (E_R + pi_inf_R + pres_R)/rho_R + H_L = (E_L + pres_L)/rho_L + H_R = (E_R + pres_R)/rho_R @:compute_average_state() @@ -3131,8 +3127,8 @@ contains s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) ! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) - xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps) - xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps) + xi_L = (s_L - vel_L(dir_idx(1)))/(s_L - s_S) + xi_R = (s_R - vel_R(dir_idx(1)))/(s_R - s_S) ! goes with numerical velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star) xi_M = (5.e-1_wp + sign(5.e-1_wp, s_S)) diff --git a/src/simulation/m_sim_helpers.fpp b/src/simulation/m_sim_helpers.fpp index abff6f32bf..4a0978919e 100644 --- a/src/simulation/m_sim_helpers.fpp +++ b/src/simulation/m_sim_helpers.fpp @@ -121,10 +121,8 @@ contains end do if (igr) then - ! igr passes q_cons_vf which stores Ẽ = E - pi_inf_mix; drop pi_inf from recovery. E = q_prim_vf(eqn_idx%E)%sf(j, k, l) - pres = (E - qv - 5.e-1_wp*rho*vel_sum)/gamma - E = E + pi_inf ! Ẽ → physical E for H below + pres = (E - pi_inf - qv - 5.e-1_wp*rho*vel_sum)/gamma else pres = q_prim_vf(eqn_idx%E)%sf(j, k, l) E = gamma*pres + pi_inf + 5.e-1_wp*rho*vel_sum + qv From ae16849f2582ebe40cecd63bcc613df4c3636093 Mon Sep 17 00:00:00 2001 From: Spencer Bryngelson Date: Wed, 6 May 2026 16:09:12 -0400 Subject: [PATCH 13/13] fix: protect HLLC xi-factor denominators from division by zero with sgm_eps --- src/simulation/m_riemann_solvers.fpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/simulation/m_riemann_solvers.fpp b/src/simulation/m_riemann_solvers.fpp index 33661b54d7..d035f32094 100644 --- a/src/simulation/m_riemann_solvers.fpp +++ b/src/simulation/m_riemann_solvers.fpp @@ -2040,8 +2040,8 @@ contains s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) ! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) - xi_L = (s_L - vel_L(dir_idx(1)))/(s_L - s_S) - xi_R = (s_R - vel_R(dir_idx(1)))/(s_R - s_S) + xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps) + xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps) ! goes with numerical star velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star) xi_M = (5.e-1_wp + sign(0.5_wp, s_S)) @@ -2331,8 +2331,8 @@ contains s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) ! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) - xi_L = (s_L - vel_L(dir_idx(1)))/(s_L - s_S) - xi_R = (s_R - vel_R(dir_idx(1)))/(s_R - s_S) + xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps) + xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps) ! goes with numerical velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star) xi_M = (5.e-1_wp + sign(5.e-1_wp, s_S)) @@ -2688,8 +2688,8 @@ contains s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) ! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) - xi_L = (s_L - vel_L(dir_idx(1)))/(s_L - s_S) - xi_R = (s_R - vel_R(dir_idx(1)))/(s_R - s_S) + xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps) + xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps) ! goes with numerical velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star) xi_M = (5.e-1_wp + sign(5.e-1_wp, s_S)) @@ -3127,8 +3127,8 @@ contains s_M = min(0._wp, s_L); s_P = max(0._wp, s_R) ! goes with q_star_L/R = xi_L/R * (variable) xi_L/R = ( ( s_L/R - u_L/R )/(s_L/R - s_star) ) - xi_L = (s_L - vel_L(dir_idx(1)))/(s_L - s_S) - xi_R = (s_R - vel_R(dir_idx(1)))/(s_R - s_S) + xi_L = (s_L - vel_L(dir_idx(1)))/min(s_L - s_S, -sgm_eps) + xi_R = (s_R - vel_R(dir_idx(1)))/max(s_R - s_S, sgm_eps) ! goes with numerical velocity in x/y/z directions xi_P/M = 0.5 +/m sgn(0.5,s_star) xi_M = (5.e-1_wp + sign(5.e-1_wp, s_S))