diff --git a/src/common/m_constants.fpp b/src/common/m_constants.fpp index 5f9088961..f317ef9bf 100644 --- a/src/common/m_constants.fpp +++ b/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 diff --git a/src/common/m_derived_types.fpp b/src/common/m_derived_types.fpp index d80c885b4..bca9ca5ee 100644 --- a/src/common/m_derived_types.fpp +++ b/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 diff --git a/src/common/m_global_parameters_common.fpp b/src/common/m_global_parameters_common.fpp index 3c230ff66..7c46f1458 100644 --- a/src/common/m_global_parameters_common.fpp +++ b/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 diff --git a/src/common/m_variables_conversion.fpp b/src/common/m_variables_conversion.fpp index f8bbd640f..feedc6894 100644 --- a/src/common/m_variables_conversion.fpp +++ b/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) diff --git a/src/post_process/m_global_parameters.fpp b/src/post_process/m_global_parameters.fpp index 34effbcb9..725c92490 100644 --- a/src/post_process/m_global_parameters.fpp +++ b/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 diff --git a/src/pre_process/m_global_parameters.fpp b/src/pre_process/m_global_parameters.fpp index d4a5c559b..2d75411c7 100644 --- a/src/pre_process/m_global_parameters.fpp +++ b/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 diff --git a/src/simulation/m_global_parameters.fpp b/src/simulation/m_global_parameters.fpp index e1a258766..618196f07 100644 --- a/src/simulation/m_global_parameters.fpp +++ b/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 diff --git a/toolchain/mfc/case_validator.py b/toolchain/mfc/case_validator.py index 9d7351245..853dfcc02 100644 --- a/toolchain/mfc/case_validator.py +++ b/toolchain/mfc/case_validator.py @@ -42,9 +42,13 @@ "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 a case that selects it may not " "set pi_inf at all: the selector determines the stiffness, not the input." + "An ideal gas is the stiffened-gas equation of state with no stiffness, so a case that selects it may not " + "set pi_inf at all: the selector determines the stiffness, not the input. 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"], }, @@ -950,17 +954,51 @@ 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'") self.prohibit( eos == eos_ideal_gas and self.get(f"fluid_pp({i})%pi_inf") is not None, f"fluid_pp({i})%eos = 'ideal_gas' has no stiffness; do not set fluid_pp({i})%pi_inf", ) + 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""" diff --git a/toolchain/mfc/params/definitions.py b/toolchain/mfc/params/definitions.py index 15f496211..48a3a0805 100644 --- a/toolchain/mfc/params/definitions.py +++ b/toolchain/mfc/params/definitions.py @@ -896,8 +896,8 @@ 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. @@ -905,10 +905,12 @@ def _load(): # 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)") diff --git a/toolchain/mfc/test_case_validator.py b/toolchain/mfc/test_case_validator.py index 4cb7fe6c9..ce0304189 100644 --- a/toolchain/mfc/test_case_validator.py +++ b/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)) diff --git a/toolchain/mfc/test_eos_mie_gruneisen.py b/toolchain/mfc/test_eos_mie_gruneisen.py new file mode 100644 index 000000000..274987b06 --- /dev/null +++ b/toolchain/mfc/test_eos_mie_gruneisen.py @@ -0,0 +1,100 @@ +"""Manufactured checks for the Mie-Gruneisen backend in m_variables_conversion. + +The Fortran is s_eos_coefficients: a linear-Hugoniot reference curve mapped onto MFC's +rho e = Gamma p + Pi form. These tests pin the maths that mapping rests on, without a solver run. +""" + +import pytest + + +def mg_reference(rho, rho0, c0, s): + """p_ref, e_ref and their d/drho for the linear-Hugoniot curve, as the Fortran computes them.""" + mu = rho / rho0 - 1.0 + if mu >= 0.0: + d = 1.0 - (s - 1.0) * mu + p = rho0 * c0**2 * mu * (1.0 + mu) / d**2 + dp_dmu = rho0 * c0**2 * ((1.0 + 2.0 * mu) * d + 2.0 * (s - 1.0) * mu * (1.0 + mu)) / d**3 + else: + p = rho0 * c0**2 * mu + dp_dmu = rho0 * c0**2 + e = p * mu / (2.0 * rho0 * (1.0 + mu)) + de_dmu = (dp_dmu * mu * (1.0 + mu) + p) / (2.0 * rho0 * (1.0 + mu) ** 2) + return p, e, dp_dmu / rho0, de_dmu / rho0 + + +def eos_coefficients(rho, rho0, c0, s, G0): + """Gamma, Pi, dPi/drho - the same three numbers s_eos_coefficients returns.""" + p, e, dp, de = mg_reference(rho, rho0, c0, s) + return 1.0 / G0, rho * e - p / G0, e + rho * de - dp / G0 + + +# Copper-like, no calibrated material: order-of-magnitude values only. +RHO0, C0, S, G0 = 8930.0, 3940.0, 1.49, 2.0 +STATES = [(r, p) for r in (7500.0, 8930.0, 9800.0, 11000.0) for p in (1e8, 5e9, 2e10)] + + +def _fd(f, x, h): + return (f(x + h) - f(x - h)) / (2.0 * h) + + +@pytest.mark.parametrize("rho", [7500.0, 8931.0, 9800.0, 11000.0]) +def test_reference_curve_derivatives(rho): + """The analytic dp_ref/drho and de_ref/drho match central differences of the curve itself.""" + _, _, dp, de = mg_reference(rho, RHO0, C0, S) + h = rho * 1e-6 + dp_fd = _fd(lambda r: mg_reference(r, RHO0, C0, S)[0], rho, h) + de_fd = _fd(lambda r: mg_reference(r, RHO0, C0, S)[1], rho, h) + assert dp == pytest.approx(dp_fd, rel=1e-8) + assert de == pytest.approx(de_fd, rel=1e-8) + + +@pytest.mark.parametrize("rho,pres", STATES) +def test_gamma_pi_form_is_the_mie_gruneisen_pressure(rho, pres): + """rho e = Gamma p + Pi inverts to exactly p = p_ref + rho Gamma_G (e - e_ref).""" + gamma, pi, _ = eos_coefficients(rho, RHO0, C0, S, G0) + e = (gamma * pres + pi) / rho # the energy MFC stores for this p + p_ref, e_ref, _, _ = mg_reference(rho, RHO0, C0, S) + p_mg = p_ref + rho * G0 * (e - e_ref) + assert p_mg == pytest.approx(pres, rel=1e-12) + + +@pytest.mark.parametrize("rho,pres", STATES) +def test_sound_speed_matches_isentrope(rho, pres): + """c^2 = [((Gamma+1)p + Pi)/rho - dPi/drho - p dGamma/drho]/Gamma equals (dp/drho)_s, integrated numerically.""" + gamma, pi, dpi = eos_coefficients(rho, RHO0, C0, S, G0) + c2 = (((gamma + 1.0) * pres + pi) / rho - dpi) / gamma # dGamma/drho = 0 for constant Gamma_G + + # Walk the isentrope: de = p/rho^2 drho, then re-evaluate p at the new state. + e = (gamma * pres + pi) / rho + h = rho * 1e-6 + + def p_at(r, e_): + g, pi_, _ = eos_coefficients(r, RHO0, C0, S, G0) + return (r * e_ - pi_) / g + + c2_fd = (p_at(rho + h, e + pres / rho**2 * h) - p_at(rho - h, e - pres / rho**2 * h)) / (2.0 * h) + assert c2 > 0.0 + assert c2 == pytest.approx(c2_fd, rel=1e-6) + + +def test_stiffened_gas_is_the_degenerate_member(): + """p_ref = -gamma pi_inf, e_ref = 0, Gamma_G = gamma - 1 reproduces MFC's stored gammas and pi_infs.""" + gam, pinf = 4.4, 6.0e8 # water-like stiffened gas + Gamma_G = gam - 1.0 + p_ref, e_ref = -gam * pinf, 0.0 + rho = 1000.0 + Gamma = 1.0 / Gamma_G + Pi = rho * e_ref - p_ref / Gamma_G + assert Gamma == pytest.approx(1.0 / (gam - 1.0), rel=1e-15) # what MFC stores as gammas(i) + assert Pi == pytest.approx(gam * pinf / (gam - 1.0), rel=1e-15) # what MFC stores as pi_infs(i) + + +def test_release_branch_is_c1_at_rho0(): + """Compression and release branches meet with matching p_ref and dp_ref at mu = 0.""" + eps = 1e-9 * RHO0 + p_up, _, dp_up, _ = mg_reference(RHO0 + eps, RHO0, C0, S) + p_dn, _, dp_dn, _ = mg_reference(RHO0 - eps, RHO0, C0, S) + assert p_up == pytest.approx(C0**2 * eps, rel=1e-6) # linear through zero from above ... + assert p_dn == pytest.approx(-(C0**2) * eps, rel=1e-6) # ... and from below, same slope + assert dp_up == pytest.approx(dp_dn, rel=1e-6) + assert dp_up == pytest.approx(C0**2, rel=1e-6) # bulk modulus rho0 c0^2 over rho0