Merge master into debug-ibm-stability

Conflict in m_data_output.fpp: this branch guards the stability scan so immersed-boundary cells cannot trip a violation, master added the alpha_rho argument for per-phase EOS evaluation and the Mie-Gruneisen Hugoniot-limit reporting. Kept the guard and moved master's work inside it. s_report_icfl_violation, new on this branch, merged without conflict but still called s_compute_cell_state and s_compute_speed_of_sound on the old signature; updated those three calls and declared alpha_rho for them. Also corrected a comment typo that fails the spell-check gate.

Author: sbryngelson

.claude/rules/common-pitfalls.md

@@ -60,6 +60,19 @@ covered in `docs/documentation/contributing.md`.
   QBMM/viscous and MHD HLLD, while both Lagrange bubble cases *complete* with out-of-tolerance
   answers. Measured 2026-08-29 on MI210. A compile-only check returns green, so any future attempt to
   drop these must run the tests, not just build.
+- **CCE OpenACC (19.0.0 through 21.0.2, `-O2`; `-O0`/`-O1` correct; OpenMP offload unaffected): a device
+  routine that contains any `GPU_LOOP` (itself or in anything it calls) must be called with scalars,
+  never with an array element as an actual argument.** Every `routine` level is affected, including the
+  conforming `loop vector` inside `routine vector`. With both ingredients present the
+  element is misaddressed: an `intent(in)` element reads as garbage, an `intent(out)` element is
+  never written. Either ingredient alone is fine, which is why master's `s_compute_pressure(q%sf(j,k,l),...)`
+  works (no loop) and `s_compute_mixture_coefficients` works (scalar actuals). PR #1811 added the
+  Newton and RK4 loops to the EOS helpers and every call that passed `%sf(j,k,l)` or `blkmod1(k,l,q)`
+  ended in `NaN(s) in timestep output` on the Frontier CCE OpenACC lanes only, bit-identical on every
+  other backend. Fix: copy elements to locals before the call, receive into a local. 37-line
+  reproducer and the bisection: sbryngelson/compiler-bugs `cce/acc-routine-element-by-reference`,
+  MFC #1815. Do not "fix" it by deleting the `seq` directives instead: they are the idiom master
+  uses in every device routine.
 - The same "call it from the loop body" rule covers `m_thermochem`: calling `get_species_*` from
   inside a `GPU_ROUTINE` rather than from the kernel gave CCE OpenMP a runtime
   `Memory access fault by GPU node-N ... Reason: Unknown` on the first step (exit 134), while every

.github/scripts/submit-slurm-job.sh

@@ -301,6 +301,17 @@ while :; do
         # another node. Note bench-pair.sh probes only after building both trees, so
         # a fault there discards those builds and the resubmit repeats them.
         faulted_node=$(bash "$SCRIPT_DIR/node-exclude.sh" node-from "$output_file")
+        # Fall back to SLURM's own record when the MFC_FAULT_NODE marker is
+        # unreadable. A job that dies before its .out is flushed (or before NFS
+        # makes it visible) leaves no marker, so node-from returns empty; the
+        # merge below then adds nothing and SLURM re-draws the SAME bad node. A
+        # dead-GPU V100 (atl1-1-02-006-34-0, cuInit 999) ate both attempts of run
+        # 34183404644 exactly this way. sacct knows the node whether or not the
+        # .out exists, so identification no longer depends on the marker.
+        if [ -z "$faulted_node" ]; then
+            faulted_node=$(sacct -j "$job_id" -X -n -o NodeList 2>/dev/null | head -n1 | tr -d ' ')
+            case "$faulted_node" in ""|None*|*[,\[]*) faulted_node="" ;; esac
+        fi
         if [ "$node_attempt" -lt "$MFC_MAX_NODE_RESUBMITS" ]; then
             node_attempt=$((node_attempt + 1))
             node_exclude=$(bash "$SCRIPT_DIR/node-exclude.sh" merge "$node_exclude" "$faulted_node")

.github/workflows/common/build.sh

@@ -41,9 +41,16 @@ source .github/scripts/retry-build.sh
 # aborts in MPI_Init ("OPAL ERROR: Unreachable in file ext3x_client.c"). mpirun
 # is unaffected, and it is how MFC launches every binary anyway. Output is left
 # on stdout so a future failure is diagnosable from the CI log.
+#
+# --bind-to none is required, not cosmetic: on a busy GPU node SLURM hands the job
+# an offset/partial cpuset (e.g. cores 32-47), and Open MPI's default --bind-to
+# core then dies with "hwloc_set_cpubind Error for bitmap" before the binary even
+# starts. The run templates, coverage_build.py and preflight.sh all already pass
+# it; this smoke-test was the lone launch that did not, and it red-crossed jobs
+# that backfilled onto shared GPU nodes.
 validate_cmd=""
 if [ "$job_cluster" = "phoenix" ]; then
-    validate_cmd='syscheck_bin=$(find build/install -name syscheck -type f 2>/dev/null | head -1); [ -z "$syscheck_bin" ] || mpirun -np 1 "$syscheck_bin"'
+    validate_cmd='syscheck_bin=$(find build/install -name syscheck -type f 2>/dev/null | head -1); [ -z "$syscheck_bin" ] || mpirun --bind-to none -np 1 "$syscheck_bin"'
 fi
 
 # --- Variant selection ---

docs/documentation/case.md

@@ -733,6 +733,7 @@ To restart the simulation from $k$-th time step, see @ref running "Restarting Ca
 | `pi_inf_wrt`            | Logical | Add the liquid stiffness function to the database |
 | `pres_inf_wrt`          | Logical | Add the liquid stiffness to the formatted database	 |
 | `c_wrt`                 | Logical | Add the sound speed to the database	 |
+| `T_wrt`                 | Logical | Add each fluid's temperature to the database (needs `cv` > 0) |
 | `omega_wrt(i)`          | Logical | Add the $i$-direction vorticity to the database	 |
 | `schlieren_wrt`         | Logical | Add the numerical schlieren to the database|
 | `qm_wrt`                | Logical | Add the Q-criterion to the database|

docs/documentation/contributing.md

@@ -464,8 +464,8 @@ means supplying these, not grepping for `gammas`:
 | `f_pressure` / `s_compute_energy` | \f$p(e)\f$ and \f$E(p)\f$ |
 | `f_bulk_modulus` | \f$K(p)\f$ - every sound speed in MFC is \f$K/\rho\f$, differing only in how phases are mixed |
 | `s_compute_speed_of_sound` / `_avg` | that mixing: Wood's law, 6-equation, bubble-diluted |
-| `f_phase_internal_energy` | per-phase internal energy (6-equation model) |
-| `f_isentrope_exponent` / `f_isentrope_pressure` / `f_pressure_on_isentrope` | the isentrope \f$p + B = \textrm{const}\,\rho^n\f$ |
+| `s_phase_internal_energy` | per-phase internal energy (6-equation model) |
+| `f_isentrope_exponent` / `f_isentrope_pressure` | the isentrope \f$p + B = \textrm{const}\,\rho^n\f$ |
 | `f_sg_thermal` | the thermal law \f$p + B = (n-1)c_v\rho T\f$ |
 
 The first six are *mechanical* - they need only \f$p, \rho, e, c\f$. The last two are *caloric* and

examples/1D_isentropic_release/case.py

@@ -0,0 +1,72 @@
+"""
+Isentropic release of a fluid whose reference curve is an isentrope (JWL or Vinet): a Riemann problem
+between two states of the same fluid. Everything left of the contact keeps the left state's entropy,
+so the fan and the left star state must lie on the closed-form isentrope through (rho0, p0).
+"""
+
+import argparse
+import json
+import math
+
+parser = argparse.ArgumentParser(description="1D JWL isentropic release")
+parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT")
+parser.add_argument("-N", type=int, default=400)
+parser.add_argument("--cfl", type=float, default=0.4)
+parser.add_argument("--eos", choices=["jwl", "vinet"], default="jwl", help="reference curve: both are isentropes")
+args = parser.parse_args()
+
+rho0, p0, rho_r, p_r = 1.0, 1.0, 0.3, 0.1
+if args.eos == "jwl":
+    fluid = {f"fluid_pp(1)%jwl_{k}": v for k, v in {"a": 6.0, "b": 0.15, "r1": 4.0, "r2": 1.0, "omega": 0.3, "rho0": rho0}.items()}
+else:
+    fluid = {f"fluid_pp(1)%vinet_{k}": v for k, v in {"k0": 2.0, "k0p": 4.0, "gruneisen": 0.3, "rho0": rho0}.items()}
+N, L, T_end = args.N, 1.0, 0.15
+c_max = math.sqrt((1.3 * p0 + 6.15) / rho0)  # generous bound on c + |u| for either fit
+dt = args.cfl * (L / N) / c_max
+Nt = math.ceil(T_end / dt)
+dt = T_end / Nt
+
+case = {
+    "run_time_info": "F",
+    "x_domain%beg": 0.0,
+    "x_domain%end": L,
+    "m": N - 1,
+    "n": 0,
+    "p": 0,
+    "dt": dt,
+    "t_step_start": 0,
+    "t_step_stop": Nt,
+    "t_step_save": Nt,
+    "num_patches": 2,
+    "model_eqns": 2,
+    "num_fluids": 1,
+    "time_stepper": 3,
+    "recon_type": "weno",
+    "weno_order": 5,
+    "weno_eps": 1.0e-16,
+    "mapped_weno": "T",
+    "riemann_solver": 2,
+    "wave_speeds": 1,
+    "avg_state": 2,
+    "bc_x%beg": -3,
+    "bc_x%end": -3,
+    "format": 1,
+    "precision": 2,
+    "prim_vars_wrt": "T",
+    "parallel_io": "F",
+    "fluid_pp(1)%eos": args.eos,
+    **fluid,
+}
+for pid, (x_c, rho, pres) in enumerate([(0.25, rho0, p0), (0.75, rho_r, p_r)], start=1):
+    case.update(
+        {
+            f"patch_icpp({pid})%geometry": 1,
+            f"patch_icpp({pid})%x_centroid": x_c,
+            f"patch_icpp({pid})%length_x": 0.5 * L,
+            f"patch_icpp({pid})%alpha_rho(1)": rho,
+            f"patch_icpp({pid})%alpha(1)": 1.0,
+            f"patch_icpp({pid})%vel(1)": 0.0,
+            f"patch_icpp({pid})%pres": pres,
+        }
+    )
+print(json.dumps(case))

examples/1D_mg_acoustic/case.py

@@ -0,0 +1,81 @@
+"""
+Right-moving acoustic pulse in a single Mie-Gruneisen fluid at its reference state.
+A rectangle patch carries a simple-wave perturbation (drho, dp = c^2 drho, du = c drho/rho0), so
+only the right-going characteristic is excited; the harness tracks the centroid of drho against
+the general analytic c. Patches rather than an analytic IC: an IC expression is compiled in, and
+every distinct one costs the test suite a full rebuild.
+"""
+
+import argparse
+import json
+import math
+
+parser = argparse.ArgumentParser(description="1D Mie-Gruneisen acoustic pulse")
+parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT")
+parser.add_argument("-N", type=int, default=200)
+parser.add_argument("--cfl", type=float, default=0.4)
+parser.add_argument("--a", type=float, default=0.0, help="Gruneisen slope: Gamma_G = Gamma_0 + a mu")
+args = parser.parse_args()
+
+rho0, p0, c0, s, gruneisen = 1.0, 1.0, 1.0, 1.5, 0.4
+c = math.sqrt(c0**2 + (1.0 + gruneisen) * p0 / rho0 + args.a * p0 / (rho0 * gruneisen))  # the frozen speed at the reference state
+amp, x0, width = 1.0e-4, 0.25, 0.1
+N, L, T_end = args.N, 1.0, 0.3
+dt = args.cfl * (L / N) / c
+Nt = math.ceil(T_end / dt)
+dt = T_end / Nt
+
+print(
+    json.dumps(
+        {
+            "run_time_info": "F",
+            "x_domain%beg": 0.0,
+            "x_domain%end": L,
+            "m": N - 1,
+            "n": 0,
+            "p": 0,
+            "dt": dt,
+            "t_step_start": 0,
+            "t_step_stop": Nt,
+            "t_step_save": Nt,
+            "model_eqns": 2,
+            "num_fluids": 1,
+            "time_stepper": 3,
+            "recon_type": "weno",
+            "weno_order": 5,
+            "weno_eps": 1.0e-16,
+            "mapped_weno": "T",
+            "riemann_solver": 2,
+            "wave_speeds": 1,
+            "avg_state": 2,
+            "bc_x%beg": -3,
+            "bc_x%end": -3,
+            "format": 1,
+            "precision": 2,
+            "prim_vars_wrt": "T",
+            "parallel_io": "F",
+            "num_patches": 2,
+            "patch_icpp(1)%geometry": 1,
+            "patch_icpp(1)%x_centroid": 0.5,
+            "patch_icpp(1)%length_x": L,
+            "patch_icpp(1)%alpha_rho(1)": rho0,
+            "patch_icpp(1)%alpha(1)": 1.0,
+            "patch_icpp(1)%vel(1)": 0.0,
+            "patch_icpp(1)%pres": p0,
+            "patch_icpp(2)%geometry": 1,
+            "patch_icpp(2)%alter_patch(1)": "T",
+            "patch_icpp(2)%x_centroid": x0,
+            "patch_icpp(2)%length_x": width,
+            "patch_icpp(2)%alpha_rho(1)": rho0 + amp,
+            "patch_icpp(2)%alpha(1)": 1.0,
+            "patch_icpp(2)%vel(1)": c / rho0 * amp,
+            "patch_icpp(2)%pres": p0 + c**2 * amp,
+            "fluid_pp(1)%eos": "mie_gruneisen",
+            "fluid_pp(1)%mg_rho0": rho0,
+            "fluid_pp(1)%mg_c0": c0,
+            "fluid_pp(1)%mg_s": s,
+            "fluid_pp(1)%mg_gruneisen": gruneisen,
+            **({"fluid_pp(1)%mg_gruneisen_a": args.a} if args.a else {}),
+        }
+    )
+)

examples/1D_mg_impact/case.py

@@ -0,0 +1,75 @@
+"""
+Symmetric impact of two Mie-Gruneisen slabs approaching at relative speed U.
+Each slab is brought to rest by a shock with particle-velocity jump U/2, so the shock state
+lies exactly on the Hugoniot u_s = c0 + s u_p; the harness checks the shock speed and the
+plateau density against that relation.
+"""
+
+import argparse
+import json
+import math
+
+parser = argparse.ArgumentParser(description="1D Mie-Gruneisen symmetric impact")
+parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT")
+parser.add_argument("-N", type=int, default=800)
+parser.add_argument("--U", type=float, default=1.0, help="closing speed of the two slabs")
+parser.add_argument("--cfl", type=float, default=0.4)
+parser.add_argument("--s2", type=float, default=0.0, help="quadratic Hugoniot coefficient u_s = c0 + s u_p + s2 u_p^2 + s3 u_p^3")
+parser.add_argument("--s3", type=float, default=0.0, help="cubic Hugoniot coefficient")
+args = parser.parse_args()
+
+rho0, p0, c0, s, gruneisen = 1.0, 1.0e-3, 1.0, 1.5, 0.4
+N, L, T_end = args.N, 1.0, 0.2
+dt = args.cfl * (L / N) / (c0 + (s + 1.0 + (args.s2 + args.s3 * args.U) * args.U) * args.U)
+Nt = math.ceil(T_end / dt)
+dt = T_end / Nt
+
+case = {
+    "run_time_info": "F",
+    "x_domain%beg": 0.0,
+    "x_domain%end": L,
+    "m": N - 1,
+    "n": 0,
+    "p": 0,
+    "dt": dt,
+    "t_step_start": 0,
+    "t_step_stop": Nt,
+    "t_step_save": Nt,
+    "num_patches": 2,
+    "model_eqns": 2,
+    "num_fluids": 1,
+    "time_stepper": 3,
+    "recon_type": "weno",
+    "weno_order": 5,
+    "weno_eps": 1.0e-16,
+    "mapped_weno": "T",
+    "riemann_solver": 2,
+    "wave_speeds": 1,
+    "avg_state": 2,
+    "bc_x%beg": -3,
+    "bc_x%end": -3,
+    "format": 1,
+    "precision": 2,
+    "prim_vars_wrt": "T",
+    "parallel_io": "F",
+    "fluid_pp(1)%eos": "mie_gruneisen",
+    "fluid_pp(1)%mg_rho0": rho0,
+    "fluid_pp(1)%mg_c0": c0,
+    "fluid_pp(1)%mg_s": s,
+    "fluid_pp(1)%mg_gruneisen": gruneisen,
+    **({"fluid_pp(1)%mg_s2": args.s2} if args.s2 else {}),
+    **({"fluid_pp(1)%mg_s3": args.s3} if args.s3 else {}),
+}
+for pid, (x_c, vel) in enumerate([(0.25, 0.5 * args.U), (0.75, -0.5 * args.U)], start=1):
+    case.update(
+        {
+            f"patch_icpp({pid})%geometry": 1,
+            f"patch_icpp({pid})%x_centroid": x_c,
+            f"patch_icpp({pid})%length_x": 0.5 * L,
+            f"patch_icpp({pid})%alpha_rho(1)": rho0,
+            f"patch_icpp({pid})%alpha(1)": 1.0,
+            f"patch_icpp({pid})%vel(1)": vel,
+            f"patch_icpp({pid})%pres": p0,
+        }
+    )
+print(json.dumps(case))

src/common/m_constants.fpp

@@ -118,6 +118,11 @@ module m_constants
     !! cannot be auto-generated, so these are hand-written.
     integer, parameter :: eos_stiffened_gas = 1
     integer, parameter :: eos_ideal_gas = 2
+    integer, parameter :: eos_mie_gruneisen = 3
+    integer, parameter :: eos_jwl = 4
+    integer, parameter :: eos_vinet = 5
+    integer, parameter :: eos_rk4_steps = 8  !< fixed-step RK4 along a phasic isentrope or a reference temperature
+    integer, parameter :: ode_isentrope = 1, ode_reference_temperature = 2  !< the two ODEs s_rk4 integrates
     integer, parameter :: num_synth_shells_max = 50  !< Max energy shells for synthetic turbulence
     integer, parameter :: num_turb_sources_max = 10  !< Max Gaussian forcing zones for synthetic turbulence