State-dependent equations of state: Mie-Gruneisen (linear and cubic Hugoniot), JWL and Vinet, with temperature, per-phase evaluation on the 5- and 6-equation models

.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

docs/documentation/case.md

@@ -732,6 +732,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
 

src/common/m_derived_types.fpp

@@ -386,22 +386,42 @@ module m_derived_types
     !> Derived type annexing the physical parameters (PP) of the fluids. These include the specific heat ratio function and liquid
     !! stiffness function.
     type physical_parameters
-        real(wp)               :: gamma          !< Sp. heat ratio
-        real(wp)               :: pi_inf         !< Liquid stiffness
-        real(wp), dimension(2) :: Re             !< Reynolds number
-        real(wp)               :: cv             !< heat capacity
-        real(wp)               :: qv             !< reference energy per unit mass for SGEOS, q (see Le Metayer (2004))
-        real(wp)               :: qvp            !< reference entropy per unit mass for SGEOS, q' (see Le Metayer (2004))
+        real(wp)               :: gamma              !< Sp. heat ratio
+        real(wp)               :: pi_inf             !< Liquid stiffness
+        real(wp), dimension(2) :: Re                 !< Reynolds number
+        real(wp)               :: cv                 !< heat capacity
+        real(wp)               :: qv                 !< reference energy per unit mass for SGEOS, q (see Le Metayer (2004))
+        real(wp)               :: qvp                !< reference entropy per unit mass for SGEOS, q' (see Le Metayer (2004))
         real(wp)               :: G
-        integer                :: eos            !< Equation of state selector (eos_* in m_constants)
-        logical                :: non_newtonian  !< Enable Herschel-Bulkley non-Newtonian viscosity
-        real(wp)               :: K              !< HB consistency index
-        real(wp)               :: nn             !< HB flow behavior index
-        real(wp)               :: tau0           !< HB yield stress (0 => power-law)
-        real(wp)               :: hb_m           !< Papanastasiou regularization parameter
-        real(wp)               :: mu_min         !< Lower viscosity clamp (inactive sentinel = dflt_real)
-        real(wp)               :: mu_max         !< Upper viscosity clamp (required when non_newtonian)
-        real(wp)               :: mu_bulk        !< Bulk viscosity for NN (inactive sentinel = dflt_real)
+        integer                :: eos                !< Equation of state selector (eos_* in m_constants)
+        real(wp)               :: mg_rho0            !< Mie-Gruneisen reference density
+        real(wp)               :: mg_c0              !< Mie-Gruneisen bulk sound speed at mg_rho0
+        real(wp)               :: mg_s               !< Mie-Gruneisen linear Hugoniot slope, u_s = c0 + s u_p
+        real(wp)               :: mg_gruneisen       !< Gruneisen coefficient Gamma_G (not the shear modulus G)
+        real(wp)               :: mg_gruneisen_a     !< d(Gamma_G)/d(mu): Gamma_G = Gamma_0 + a mu, zero keeps it constant
+        real(wp)               :: mg_t0              !< temperature at the reference density (for T output)
+        real(wp)               :: mg_s2, mg_s3       !< u_s = c0 + s u_p + s2 u_p^2 + s3 u_p^3; zero keeps the fit linear
+        real(wp)               :: jwl_a              !< JWL A
+        real(wp)               :: jwl_b              !< JWL B
+        real(wp)               :: jwl_r1             !< JWL R1
+        real(wp)               :: jwl_r2             !< JWL R2
+        real(wp)               :: jwl_omega          !< JWL omega (its Gruneisen coefficient)
+        real(wp)               :: jwl_rho0           !< JWL reference density
+        real(wp)               :: jwl_t0             !< temperature at the reference density (for T output)
+        real(wp)               :: vinet_k0           !< Vinet bulk modulus at rho0
+        real(wp)               :: vinet_k0p          !< Vinet pressure derivative of the bulk modulus
+        real(wp)               :: vinet_rho0         !< Vinet reference density
+        real(wp)               :: vinet_gruneisen    !< Gruneisen coefficient at rho0
+        real(wp)               :: vinet_gruneisen_a  !< d(Gamma_G)/d(mu)
+        real(wp)               :: vinet_t0           !< temperature at the reference density (for T output)
+        logical                :: non_newtonian      !< Enable Herschel-Bulkley non-Newtonian viscosity
+        real(wp)               :: K                  !< HB consistency index
+        real(wp)               :: nn                 !< HB flow behavior index
+        real(wp)               :: tau0               !< HB yield stress (0 => power-law)
+        real(wp)               :: hb_m               !< Papanastasiou regularization parameter
+        real(wp)               :: mu_min             !< Lower viscosity clamp (inactive sentinel = dflt_real)
+        real(wp)               :: mu_max             !< Upper viscosity clamp (required when non_newtonian)
+        real(wp)               :: mu_bulk            !< Bulk viscosity for NN (inactive sentinel = dflt_real)
     end type physical_parameters
 
     !> Derived type annexing the physical parameters required for sub-grid bubble models

src/common/m_global_parameters_common.fpp

@@ -54,6 +54,17 @@ module m_global_parameters_common
     !! written as p + B = const*rho**n.
     real(wp), allocatable, dimension(:) :: gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps
     $:GPU_DECLARE(create='[gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps]')
+    !> Per-fluid EOS selector and Mie-Gruneisen reference curve, resolved once at init like the arrays above
+    integer, allocatable, dimension(:) :: eoss
+    !> Reference state and Gruneisen closure Gamma_G = Gamma_0 + a mu of a state-dependent EOS, whatever its family
+    real(wp), allocatable, dimension(:) :: rho0s, t0s, gruneisen0s, gruneisen_as
+    real(wp), allocatable, dimension(:) :: mg_c0s, mg_ss, mg_s2s, mg_s3s
+    real(wp), allocatable, dimension(:) :: mg_mu_maxs  !< Compression at which a cubic Hugoniot fit turns over
+    real(wp), allocatable, dimension(:) :: jwl_as, jwl_bs, jwl_r1s, jwl_r2s
+    real(wp), allocatable, dimension(:) :: vinet_k0s, vinet_k0ps
+    logical :: any_state_dependent_eos                 !< True when some fluid's coefficients vary with density; set at init
+    $:GPU_DECLARE(create='[eoss, rho0s, t0s, gruneisen0s, gruneisen_as, mg_c0s, mg_ss, mg_s2s, mg_s3s, mg_mu_maxs, jwl_as, &
+                  & jwl_bs, jwl_r1s, jwl_r2s, vinet_k0s, vinet_k0ps, any_state_dependent_eos]')
     !> @}
 
     !> @name Fluids participating in shear and bulk viscosity