Add a Mie-Gruneisen equation-of-state backend

src/common/m_constants.fpp

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

@@ -394,6 +394,10 @@ module m_derived_types
         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)
+        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)
         logical                :: non_newtonian  !< Enable Herschel-Bulkley non-Newtonian viscosity
         real(wp)               :: K              !< HB consistency index
         real(wp)               :: nn             !< HB flow behavior index

src/common/m_global_parameters_common.fpp

@@ -54,6 +54,11 @@ 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
+    real(wp), allocatable, dimension(:) :: mg_rho0s, mg_c0s, mg_ss, mg_gruneisens
+    logical :: any_state_dependent_eos  !< True when some fluid's coefficients vary with density; set at init
+    $:GPU_DECLARE(create='[eoss, mg_rho0s, mg_c0s, mg_ss, mg_gruneisens, any_state_dependent_eos]')
     !> @}
 
     !> @name Fluids participating in shear and bulk viscosity

src/common/m_variables_conversion.fpp

@@ -28,8 +28,8 @@ module m_variables_conversion
         & s_compute_species_fraction, s_compute_mixture_coefficients, s_compute_energy, s_compute_speed_of_sound, f_bulk_modulus, &
         & f_pressure, f_phase_internal_energy, f_isentrope_exponent, f_isentrope_pressure, f_sg_thermal, f_pressure_on_isentrope, &
         & s_compute_mixture_coefficients_dt, s_compute_speed_of_sound_avg, s_compute_fast_magnetosonic_speed, f_elastic_energy, &
-        & f_hypoelastic_energy, f_relativistic_enthalpy, s_finalize_variables_conversion_module, gammas, isentrope_n, pi_infs, &
-        & isentrope_B, cvs, qvs, qvps
+        & f_hypoelastic_energy, f_relativistic_enthalpy, s_eos_coefficients, s_finalize_variables_conversion_module, gammas, &
+        & isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps
 
     real(wp), allocatable, dimension(:)   :: Gs_vc
     integer, allocatable, dimension(:)    :: bubrs_vc
@@ -267,6 +267,8 @@ contains
         $:GPU_UPDATE(device='[enforce_density_floor_vc, preserve_qbmm_number_vc, lagrange_beta_index_vc]')
 
         @:ALLOCATE(gammas (1:num_fluids))
+        @:ALLOCATE(eoss (1:num_fluids), mg_rho0s (1:num_fluids), mg_c0s (1:num_fluids), mg_ss (1:num_fluids), &
+                   & mg_gruneisens (1:num_fluids))
         @:ALLOCATE(isentrope_n (1:num_fluids))
         @:ALLOCATE(pi_infs(1:num_fluids))
         @:ALLOCATE(isentrope_B(1:num_fluids))
@@ -275,6 +277,7 @@ contains
         @:ALLOCATE(qvps    (1:num_fluids))
         @:ALLOCATE(Gs_vc     (1:num_fluids))
 
+        any_state_dependent_eos = .false.
         do i = 1, num_fluids
             gammas(i) = fluid_pp(i)%gamma
             isentrope_n(i) = f_isentrope_exponent(gammas(i))
@@ -292,8 +295,15 @@ contains
             cvs(i) = fluid_pp(i)%cv
             qvs(i) = fluid_pp(i)%qv
             qvps(i) = fluid_pp(i)%qvp
+            eoss(i) = fluid_pp(i)%eos
+            mg_rho0s(i) = fluid_pp(i)%mg_rho0
+            mg_c0s(i) = fluid_pp(i)%mg_c0
+            mg_ss(i) = fluid_pp(i)%mg_s
+            mg_gruneisens(i) = fluid_pp(i)%mg_gruneisen
+            if (fluid_pp(i)%eos == eos_mie_gruneisen) any_state_dependent_eos = .true.
         end do
-        $:GPU_UPDATE(device='[gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc]')
+        $:GPU_UPDATE(device='[gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc, eoss, mg_rho0s, mg_c0s, mg_ss, &
+                     & mg_gruneisens, any_state_dependent_eos]')
 
         @:ALLOCATE(Res_vc(1:2, 1:max(1, Re_size_max)))
         Res_vc = dflt_real
@@ -1173,7 +1183,7 @@ contains
 
         if (allocated(rho_sf)) deallocate (rho_sf, gamma_sf, pi_inf_sf)
 
-        @:DEALLOCATE(gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc)
+        @:DEALLOCATE(gammas, isentrope_n, pi_infs, isentrope_B, cvs, qvs, qvps, Gs_vc, eoss, mg_rho0s, mg_c0s, mg_ss, mg_gruneisens)
         if (allocated(bubrs_vc)) then
             @:DEALLOCATE(bubrs_vc)
         end if
@@ -1284,6 +1294,49 @@ contains
 
     end subroutine s_compute_energy
 
+    !> Coefficients of fluid i at density rho in the form rho e = Gamma p + Pi that every operator here consumes, with dPi/drho. Any
+    !! Mie-Gruneisen EOS p = p_ref + rho Gamma_G (e - e_ref) is this form with Gamma = 1/Gamma_G and Pi = rho e_ref - p_ref/Gamma_G;
+    !! a new family adds one case supplying its reference curve. Gamma_G is constant, so dGamma/drho = 0. Stiffened and ideal gas
+    !! keep the constants resolved at init, bit for bit.
+    subroutine s_eos_coefficients(rho, i, gamma, pi_inf, dpi)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+
+        real(wp), intent(in)  :: rho
+        integer, intent(in)   :: i
+        real(wp), intent(out) :: gamma, pi_inf, dpi
+        real(wp)              :: mu, d, p_ref, e_ref, dp_dmu, de_dmu, G0
+
+        select case (eoss(i))
+        case (eos_mie_gruneisen)
+            ! Linear-Hugoniot reference curve, u_s = c0 + s u_p: p_H = rho0 c0^2 mu (1 + mu)/(1 - (s - 1) mu)^2 on
+            ! compression, extended linearly on release, with the Hugoniot energy e_H = p_H mu/(2 rho0 (1 + mu)).
+            ! Pole at mu = 1/(s - 1); the validator warns when the initial state is near it.
+            mu = rho/mg_rho0s(i) - 1._wp
+            if (mu >= 0._wp) then
+                d = 1._wp - (mg_ss(i) - 1._wp)*mu
+                p_ref = mg_rho0s(i)*mg_c0s(i)**2*mu*(1._wp + mu)/(d*d)
+                dp_dmu = mg_rho0s(i)*mg_c0s(i)**2*((1._wp + 2._wp*mu)*d + 2._wp*(mg_ss(i) - 1._wp)*mu*(1._wp + mu))/(d*d*d)
+            else
+                p_ref = mg_rho0s(i)*mg_c0s(i)**2*mu
+                dp_dmu = mg_rho0s(i)*mg_c0s(i)**2
+            end if
+            e_ref = p_ref*mu/(2._wp*mg_rho0s(i)*(1._wp + mu))
+            de_dmu = (dp_dmu*mu*(1._wp + mu) + p_ref)/(2._wp*mg_rho0s(i)*(1._wp + mu)**2)
+            G0 = mg_gruneisens(i)
+        case default
+            gamma = gammas(i)
+            pi_inf = pi_infs(i)
+            dpi = 0._wp
+            return
+        end select
+
+        gamma = 1._wp/G0
+        pi_inf = rho*e_ref - p_ref/G0
+        dpi = e_ref + (rho*de_dmu - dp_dmu/G0)/mg_rho0s(i)  ! d/drho = (1/rho0) d/dmu
+
+    end subroutine s_eos_coefficients
+
     !> Exponent of the stiffened-gas isentrope p + B = const rho**n. Precomputed per fluid as isentrope_n.
     function f_isentrope_exponent(gamma) result(n)
 

src/post_process/m_global_parameters.fpp

@@ -203,6 +203,10 @@ contains
         ! Fluids physical parameters (post-specific; G = dflt_real differs from pre/sim)
         do i = 1, num_fluids_max
             fluid_pp(i)%eos = eos_stiffened_gas
+            fluid_pp(i)%mg_rho0 = dflt_real
+            fluid_pp(i)%mg_c0 = dflt_real
+            fluid_pp(i)%mg_s = dflt_real
+            fluid_pp(i)%mg_gruneisen = dflt_real
             fluid_pp(i)%gamma = dflt_real
             fluid_pp(i)%pi_inf = dflt_real
             fluid_pp(i)%cv = 0._wp

src/pre_process/m_global_parameters.fpp

@@ -387,6 +387,10 @@ contains
         ! Fluids physical parameters
         do i = 1, num_fluids_max
             fluid_pp(i)%eos = eos_stiffened_gas
+            fluid_pp(i)%mg_rho0 = dflt_real
+            fluid_pp(i)%mg_c0 = dflt_real
+            fluid_pp(i)%mg_s = dflt_real
+            fluid_pp(i)%mg_gruneisen = dflt_real
             fluid_pp(i)%gamma = dflt_real
             fluid_pp(i)%pi_inf = dflt_real
             fluid_pp(i)%cv = 0._wp

src/simulation/m_global_parameters.fpp

@@ -444,6 +444,10 @@ contains
         ! Fluids physical parameters (sim-specific; Re(:) and G=0._wp differ from post)
         do i = 1, num_fluids_max
             fluid_pp(i)%eos = eos_stiffened_gas
+            fluid_pp(i)%mg_rho0 = dflt_real
+            fluid_pp(i)%mg_c0 = dflt_real
+            fluid_pp(i)%mg_s = dflt_real
+            fluid_pp(i)%mg_gruneisen = dflt_real
             fluid_pp(i)%gamma = dflt_real
             fluid_pp(i)%pi_inf = dflt_real
             fluid_pp(i)%cv = 0._wp

toolchain/mfc/case_validator.py

@@ -42,10 +42,12 @@
     "check_eos_selector": {
         "title": "Equation of State Selector",
         "category": "Thermodynamic Constraints",
-        "math": r"\Pi_\infty = 0 \;\; \text{for an ideal gas}",
+        "math": r"\rho e = \Gamma\,p + \Pi(\rho), \quad \Gamma = 1/\Gamma_G, \quad \Pi(\rho) = \rho\, e_{\mathrm{ref}}(\rho) - p_{\mathrm{ref}}(\rho)/\Gamma_G",
         "explanation": (
-            "An ideal gas is the stiffened-gas equation of state with no stiffness, so selecting it and also supplying a "
-            "nonzero pi_inf is contradictory. The selector determines the stiffness, not the input."
+            "Every backend supplies the same two coefficients. An ideal gas is the stiffened-gas form with no stiffness, so "
+            "selecting it and also supplying a nonzero pi_inf is contradictory. Mie-Gruneisen supplies a linear-Hugoniot "
+            "reference curve (mg_rho0, mg_c0, mg_s, mg_gruneisen) whose reference energy already carries the formation energy, so "
+            "qv must be zero; its parameters are read only when that backend is selected."
         ),
         "references": ["Wilfong26"],
     },
@@ -951,18 +953,52 @@ def check_eos_selector(self):
             return
         eos_names = CONSTRAINTS["fluid_pp(1)%eos"]["names"]
         eos_ideal_gas = eos_names["ideal_gas"]
+        eos_mg = eos_names["mie_gruneisen"]
         eos_values = set(eos_names.values())
         bub_fac = 1 if self.get("bubbles_euler", "F") == "T" else 0
         for i in range(1, num_fluids + 1 + bub_fac):
             eos = self.get(f"fluid_pp({i})%eos")
+            mg = {k: self.get(f"fluid_pp({i})%mg_{k}") for k in ("rho0", "c0", "s", "gruneisen")}
+            # An unset selector is stiffened gas, so stray mg_* parameters must be caught before the early return.
+            self.prohibit(
+                (eos if eos is not None else eos_names["stiffened_gas"]) != eos_mg and any(v is not None for v in mg.values()),
+                f"fluid_pp({i})%mg_* are only read when fluid_pp({i})%eos = 'mie_gruneisen'",
+            )
             if eos is None:
                 continue
-            self.prohibit(eos not in eos_values, f"fluid_pp({i})%eos must be 'stiffened_gas' or 'ideal_gas'")
+            self.prohibit(eos not in eos_values, f"fluid_pp({i})%eos must be 'stiffened_gas', 'ideal_gas' or 'mie_gruneisen'")
             pi_inf = self.get(f"fluid_pp({i})%pi_inf")
             self.prohibit(
                 eos == eos_ideal_gas and pi_inf is not None and pi_inf != 0,
                 f"fluid_pp({i})%eos = 'ideal_gas' requires fluid_pp({i})%pi_inf = 0 (an ideal gas has no stiffness)",
             )
+            self.prohibit(
+                eos == eos_mg and any(v is None for v in mg.values()),
+                f"fluid_pp({i})%eos = 'mie_gruneisen' requires fluid_pp({i})%mg_rho0, mg_c0, mg_s and mg_gruneisen",
+            )
+            if eos == eos_mg and all(v is not None for v in mg.values()):
+                self.prohibit(mg["rho0"] <= 0 or mg["c0"] <= 0 or mg["gruneisen"] <= 0, f"fluid_pp({i})%mg_rho0, mg_c0 and mg_gruneisen must be positive")
+                self.prohibit(mg["s"] < 1, f"fluid_pp({i})%mg_s must be >= 1 (u_s = c0 + s u_p; s < 1 gives no shock)")
+                qv = self.get(f"fluid_pp({i})%qv")
+                self.prohibit(
+                    qv is not None and qv != 0,
+                    f"fluid_pp({i})%qv must be 0 with eos = 'mie_gruneisen'; the reference energy e_ref(rho) carries it",
+                )
+                # The linear Hugoniot has a pole at mu = 1/(s - 1), its maximum compression. Only the initial
+                # state can be checked here; the solver does not guard the runtime density.
+                if mg["s"] > 1:
+                    rho_pole = mg["rho0"] * (1.0 + 1.0 / (mg["s"] - 1.0))
+                    num_patches = self.get("num_patches", 0) or 0
+                    for j in range(1, num_patches + 1):
+                        ar = self.get(f"patch_icpp({j})%alpha_rho({i})")
+                        a = self.get(f"patch_icpp({j})%alpha({i})")
+                        if not all(isinstance(v, (int, float)) for v in (ar, a)) or a <= 0:
+                            continue
+                        self.warn(
+                            ar / a > 0.8 * rho_pole,
+                            f"patch_icpp({j}) starts fluid {i} at rho = {ar/a:.4g}, within 20% of the Mie-Gruneisen "
+                            f"Hugoniot pole rho0*s/(s-1) = {rho_pole:.4g}; the reference curve is unphysical beyond it",
+                        )
 
     def check_stiffened_eos(self):
         """Checks constraints on stiffened equation of state fluids parameters"""

toolchain/mfc/params/definitions.py

@@ -896,19 +896,21 @@ def _load():
                 _r(f"{px}sph_har_coeff({ll},{mm})", REAL)
 
     # Values must match the hand-written eos_* constants in src/common/m_constants.fpp.
-    _EOS_NAMES = {"stiffened_gas": 1, "ideal_gas": 2}
-    _EOS_VALUE_LABELS = {1: "stiffened-gas", 2: "ideal-gas"}
+    _EOS_NAMES = {"stiffened_gas": 1, "ideal_gas": 2, "mie_gruneisen": 3}
+    _EOS_VALUE_LABELS = {1: "stiffened-gas", 2: "ideal-gas", 3: "Mie-Gruneisen"}
 
     # fluid_pp (10 fluids)
     # Members present in physical_parameters: gamma, pi_inf, Re, cv, qv, qvp, G.
     # mul0/ss/pv/gamma_v/M_v/mu_v/k_v/cp_v/D_v were removed from the Fortran type
     # by upstream #1085/#1093 — they must NOT be registered (namelist read would crash).
     for f in range(1, NF + 1):
         px = f"fluid_pp({f})%"
-        CONSTRAINTS[f"fluid_pp({f})%eos"] = {"choices": [1, 2], "value_labels": _EOS_VALUE_LABELS, "names": _EOS_NAMES}
+        CONSTRAINTS[f"fluid_pp({f})%eos"] = {"choices": [1, 2, 3], "value_labels": _EOS_VALUE_LABELS, "names": _EOS_NAMES}
         for a, sym in [("gamma", r"\f$\gamma_k\f$"), ("pi_inf", r"\f$\pi_{\infty,k}\f$"), ("cv", r"\f$c_{v,k}\f$"), ("qv", r"\f$q_{v,k}\f$"), ("qvp", r"\f$q'_{v,k}\f$")]:
             _r(f"{px}{a}", REAL, math=sym)
         _r(f"{px}eos", INT, math=r"\f$\mathrm{EOS}_k\f$")
+        for a, sym in [("mg_rho0", r"\f$\rho_{0,k}\f$"), ("mg_c0", r"\f$c_{0,k}\f$"), ("mg_s", r"\f$s_k\f$"), ("mg_gruneisen", r"\f$\Gamma_{G,k}\f$")]:
+            _r(f"{px}{a}", REAL, math=sym)
         _r(f"{px}G", REAL, {"hypoelasticity"}, math=r"\f$G_k\f$")
         _r(f"{px}Re(1)", REAL, {"viscosity"}, math=r"\f$\mathrm{Re}_k\f$ (shear)")
         _r(f"{px}Re(2)", REAL, {"viscosity"}, math=r"\f$\mathrm{Re}_k\f$ (bulk)")

toolchain/mfc/test_case_validator.py

@@ -389,3 +389,40 @@ def test_not_tripped_without_alt_soundspeed(self):
 
 if __name__ == "__main__":
     unittest.main()
+
+
+class TestMieGruneisenSelector(ConstraintTestCase):
+    """fluid_pp(i)%eos = 'mie_gruneisen' requires its reference curve and forbids what it makes redundant."""
+
+    MG = {"fluid_pp(1)%eos": 3, "fluid_pp(1)%mg_rho0": 8930.0, "fluid_pp(1)%mg_c0": 3940.0, "fluid_pp(1)%mg_s": 1.49, "fluid_pp(1)%mg_gruneisen": 2.0}
+
+    def warnings_for(self, params):
+        validator = CaseValidator(dict(params))
+        validator.validate("simulation")
+        return "\n".join(validator.warnings)
+
+    def test_accepts_complete_curve(self):
+        self.assertAccepts({**BASE, **self.MG, "fluid_pp(1)%qv": 0.0})
+
+    def test_requires_all_four_parameters(self):
+        p = {**BASE, **self.MG}
+        del p["fluid_pp(1)%mg_s"]
+        self.assertRejects(p, "requires fluid_pp(1)%mg_rho0, mg_c0, mg_s and mg_gruneisen")
+
+    def test_rejects_parameters_under_stiffened_gas(self):
+        self.assertRejects({**BASE, "fluid_pp(1)%mg_c0": 3940.0}, "only read when fluid_pp(1)%eos = 'mie_gruneisen'")
+
+    def test_rejects_formation_energy(self):
+        self.assertRejects({**BASE, **self.MG, "fluid_pp(1)%qv": 1.0e5}, "qv must be 0 with eos = 'mie_gruneisen'")
+
+    def test_rejects_slope_below_one(self):
+        self.assertRejects({**BASE, **self.MG, "fluid_pp(1)%mg_s": 0.9}, "mg_s must be >= 1")
+
+    def test_warns_near_the_hugoniot_pole(self):
+        # rho_pole = rho0 * s/(s-1) = 8930 * 1.49/0.49 ~ 27157; start at 90% of it
+        p = {**BASE, **self.MG, "patch_icpp(1)%alpha(1)": 1.0, "patch_icpp(1)%alpha_rho(1)": 0.9 * 8930.0 * 1.49 / 0.49}
+        self.assertIn("Hugoniot pole", self.warnings_for(p))
+
+    def test_no_pole_warning_at_reference_density(self):
+        p = {**BASE, **self.MG, "patch_icpp(1)%alpha(1)": 1.0, "patch_icpp(1)%alpha_rho(1)": 8930.0}
+        self.assertNotIn("Hugoniot pole", self.warnings_for(p))