From 7797626163e9bf18701eda2503a32231b6f2a854 Mon Sep 17 00:00:00 2001
From: Fahad Nabid <55717325+fahnab666@users.noreply.github.com>
Date: Tue, 21 Jul 2026 02:51:51 -0400
Subject: [PATCH 1/3] Add generic Mie-Gruneisen equation-of-state core

---
 docs/module_categories.json        |   3 +-
 misc/mg_eos_verification.py        | 119 +++++++++++++++++++++++++++++
 src/common/m_eos_mie_gruneisen.fpp |  88 +++++++++++++++++++++
 3 files changed, 209 insertions(+), 1 deletion(-)
 create mode 100644 misc/mg_eos_verification.py
 create mode 100644 src/common/m_eos_mie_gruneisen.fpp

diff --git a/docs/module_categories.json b/docs/module_categories.json
index 1391908cad..af20960047 100644
--- a/docs/module_categories.json
+++ b/docs/module_categories.json
@@ -34,7 +34,8 @@
             "m_acoustic_src",
             "m_body_forces",
             "m_pressure_relaxation",
-            "m_collisions"
+            "m_collisions",
+            "m_eos_mie_gruneisen"
         ]
     },
     {
diff --git a/misc/mg_eos_verification.py b/misc/mg_eos_verification.py
new file mode 100644
index 0000000000..227bbd1759
--- /dev/null
+++ b/misc/mg_eos_verification.py
@@ -0,0 +1,119 @@
+#!/usr/bin/env python3
+"""Manufactured verification for the generic Mie-Gruneisen EOS core in
+src/common/m_eos_mie_gruneisen.fpp.
+
+The Mie-Gruneisen backend is not yet selectable at runtime, so it has no golden
+regression case; this standalone check is its verification instead. It mirrors the
+Fortran formulas line for line and exercises the two properties that fix the generic
+assembler: exact reduction to the existing stiffened-gas EOS, and agreement of the
+analytic sound speed with a finite-difference of the frozen identity for a density
+dependent reference curve.
+
+Reference: Arienti, Morano and Shepherd, GALCIT FM99-8 (2004); Menikoff and Plohr,
+Rev. Mod. Phys. 61, 75 (1989). Run with python3; exits nonzero on any failure.
+"""
+
+import sys
+
+
+# --- generic Mie-Gruneisen assembler (mirrors m_eos_mie_gruneisen.fpp) --------------
+def mg_pressure(gamma_mg, rho, e_int, p_ref, e_ref):
+    return p_ref + gamma_mg * rho * (e_int - e_ref)
+
+
+def mg_internal_energy(gamma_mg, rho, pres, p_ref, e_ref):
+    return e_ref + (pres - p_ref) / (gamma_mg * rho)
+
+
+def mg_sound_speed_sq(gamma_mg, rho, pres, p_ref, dp_ref, de_ref):
+    return dp_ref - gamma_mg * rho * de_ref + ((1.0 + gamma_mg) * pres - p_ref) / rho
+
+
+def stiffened_reference(gamma_s, pi_inf, rho):
+    """Stiffened gas as a Mie-Gruneisen reference curve."""
+    return gamma_s - 1.0, -pi_inf, 0.0, pi_inf / rho, -pi_inf / rho**2
+
+
+# --- analytic stiffened-gas oracle --------------------------------------------------
+def sg_pressure(gamma_s, rho, e_int, pi_inf):
+    return (gamma_s - 1.0) * rho * e_int - gamma_s * pi_inf
+
+
+def sg_sound_speed_sq(gamma_s, rho, pres, pi_inf):
+    return gamma_s * (pres + pi_inf) / rho
+
+
+def rel(a, b):
+    return abs(a - b) / max(1.0, abs(b))
+
+
+# --- Test A: exact reduction to stiffened gas (ideal gas is the pi_inf = 0 sub-case) -
+def test_reduction():
+    # (gamma_s, pi_inf, rho, e): ideal air, water (stiffened), a stiff solid, monatomic
+    cases = [
+        (1.4, 0.0, 1.2, 2.5e5),
+        (4.4, 6.0e8, 1000.0, 3.0e5),
+        (3.0, 1.0e5, 50.0, 2.0e4),
+        (1.667, 0.0, 0.1, 2.0e6),
+    ]
+    worst = 0.0
+    ok = True
+    for gamma_s, pi_inf, rho, e in cases:
+        gamma_mg, p_ref, dp_ref, e_ref, de_ref = stiffened_reference(gamma_s, pi_inf, rho)
+        p_mg = mg_pressure(gamma_mg, rho, e, p_ref, e_ref)
+        e_back = mg_internal_energy(gamma_mg, rho, p_mg, p_ref, e_ref)
+        c2_mg = mg_sound_speed_sq(gamma_mg, rho, p_mg, p_ref, dp_ref, de_ref)
+        errs = (
+            rel(p_mg, sg_pressure(gamma_s, rho, e, pi_inf)),
+            rel(e_back, e),
+            rel(c2_mg, sg_sound_speed_sq(gamma_s, rho, p_mg, pi_inf)),
+        )
+        worst = max(worst, *errs)
+        ok = ok and all(x < 1.0e-13 for x in errs)
+    print(f"[A] stiffened-gas reduction: worst rel err = {worst:.2e}  ->  {'PASS' if ok else 'FAIL'}")
+    return ok
+
+
+# --- Test B: analytic c^2 vs finite difference for a JWL-form isentrope reference ----
+# Synthetic (non-calibrated) products isentrope, exercising a density-dependent p_ref.
+def jwl_reference(rho, A, B, R1, R2, omega, rho0):
+    V = rho0 / rho
+    e1, e2 = A * pow(2.718281828459045, -R1 * V), B * pow(2.718281828459045, -R2 * V)
+    p_ref = e1 + e2
+    e_ref = e1 / (rho0 * R1) + e2 / (rho0 * R2)
+    dp_ref = (V / rho) * (R1 * e1 + R2 * e2)
+    de_ref = p_ref / rho**2
+    return omega, p_ref, dp_ref, e_ref, de_ref
+
+
+def test_soundspeed_fd():
+    A, B, R1, R2, omega, rho0 = 5.0e11, 8.0e9, 4.5, 1.2, 0.3, 1600.0
+    worst = 0.0
+    ok = True
+    for rho in (900.0, 1600.0, 2400.0):
+        for e in (1.0e6, 4.0e6):
+
+            def p_of(r, en):
+                g, p_r, _, e_r, _ = jwl_reference(r, A, B, R1, R2, omega, rho0)
+                return mg_pressure(g, r, en, p_r, e_r)
+
+            g, p_r, dp_r, e_r, de_r = jwl_reference(rho, A, B, R1, R2, omega, rho0)
+            pres = mg_pressure(g, rho, e, p_r, e_r)
+            c2_analytic = mg_sound_speed_sq(g, rho, pres, p_r, dp_r, de_r)
+
+            # frozen identity c^2 = dp/drho|_e + (p/rho^2) dp/de|_rho, dp/de|_rho = omega*rho
+            h = rho * 1.0e-6
+            dpdrho_e = (p_of(rho + h, e) - p_of(rho - h, e)) / (2.0 * h)
+            c2_fd = dpdrho_e + (pres / rho**2) * (omega * rho)
+
+            err = rel(c2_analytic, c2_fd)
+            worst = max(worst, err)
+            ok = ok and err < 1.0e-6
+    print(f"[B] sound speed vs finite difference: worst rel err = {worst:.2e}  ->  {'PASS' if ok else 'FAIL'}")
+    return ok
+
+
+if __name__ == "__main__":
+    passed = test_reduction() & test_soundspeed_fd()
+    print("ALL PASS" if passed else "FAILURE")
+    sys.exit(0 if passed else 1)
diff --git a/src/common/m_eos_mie_gruneisen.fpp b/src/common/m_eos_mie_gruneisen.fpp
new file mode 100644
index 0000000000..30865d7ced
--- /dev/null
+++ b/src/common/m_eos_mie_gruneisen.fpp
@@ -0,0 +1,88 @@
+!>
+!! @file m_eos_mie_gruneisen.fpp
+!! @brief Contains module m_eos_mie_gruneisen
+
+#:include 'macros.fpp'
+
+!> @brief Generic Mie-Gruneisen equation of state with a constant Gruneisen coefficient, structured as a reference-curve provider
+!! plus a curve-agnostic assembler. A Mie-Gruneisen material obeys p = p_ref(rho) + Gamma*rho*(e - e_ref(rho)); the reference pair
+!! (p_ref, e_ref) and its density derivatives are the only thing that distinguishes one Mie-Gruneisen material from another.
+!! Following the production hydrocode structure (Arienti, Morano and Shepherd, GALCIT FM99-8, 2004; LLNL singularity-eos), one
+!! provider returns the reference curve for a chosen material and one assembler consumes it, so stiffened gas, ideal gas, and JWL
+!! differ only in the provider. This PR ships the stiffened-gas reference curve, which recovers the existing stiffened-gas pressure
+!! and sound speed exactly (Gamma = gamma_s - 1, p_ref = -pi_inf, e_ref = pi_inf/rho); the JWL principal isentrope is the next
+!! provider, added with its consumer. This is a leaf module (it uses only the precision kinds) and no solver path calls it yet: it
+!! is the verified thermodynamic core for the Mie-Gruneisen backend, landed ahead of the reference curve that consumes it.
+module m_eos_mie_gruneisen
+
+    use m_precision_select
+
+    implicit none
+
+    private
+    public :: s_mg_stiffened_reference, f_mg_pressure, f_mg_internal_energy, f_mg_sound_speed_sq
+
+contains
+
+    !> Stiffened gas expressed as a Mie-Gruneisen reference curve. Returns the Gruneisen coefficient, the reference pressure and
+    !! energy, and their density derivatives, which together make the generic assembler reproduce stiffened gas. This is the first
+    !! concrete reference curve and the template the JWL principal isentrope follows. The reference pressure is density-independent
+    !! for stiffened gas, so its derivative is zero; ideal gas is the pi_inf = 0 sub-case.
+    pure subroutine s_mg_stiffened_reference(gamma_s, pi_inf, rho, gamma_mg, p_ref, dp_ref, e_ref, de_ref)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+
+        real(wp), intent(in)  :: gamma_s, pi_inf, rho
+        real(wp), intent(out) :: gamma_mg, p_ref, dp_ref, e_ref, de_ref
+
+        gamma_mg = gamma_s - 1._wp
+        p_ref = -pi_inf
+        dp_ref = 0._wp
+        e_ref = pi_inf/rho
+        de_ref = -pi_inf/rho**2
+
+    end subroutine s_mg_stiffened_reference
+
+    !> Mie-Gruneisen pressure from density and specific internal energy against a reference curve: p = p_ref + Gamma*rho*(e -
+    !! e_ref).
+    pure function f_mg_pressure(gamma_mg, rho, e_int, p_ref, e_ref) result(pres)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+
+        real(wp), intent(in) :: gamma_mg, rho, e_int, p_ref, e_ref
+        real(wp)             :: pres
+
+        pres = p_ref + gamma_mg*rho*(e_int - e_ref)
+
+    end function f_mg_pressure
+
+    !> Specific internal energy from pressure, the inverse of f_mg_pressure: e = e_ref + (p - p_ref)/(Gamma*rho). The Gamma*rho
+    !! divisor is positive for every admissible state (Gamma > 0, rho > 0).
+    pure function f_mg_internal_energy(gamma_mg, rho, pres, p_ref, e_ref) result(e_int)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+
+        real(wp), intent(in) :: gamma_mg, rho, pres, p_ref, e_ref
+        real(wp)             :: e_int
+
+        e_int = e_ref + (pres - p_ref)/(gamma_mg*rho)
+
+    end function f_mg_internal_energy
+
+    !> Squared sound speed c^2 = dp/drho at constant entropy for a constant Gruneisen coefficient, from the frozen identity c^2 =
+    !! dp/drho|_e + (p/rho^2) dp/de|_rho with dp/de|_rho = Gamma*rho (Arienti et al. 2004, Eq. 12; Menikoff and Plohr 1989). Written
+    !! in terms of pressure so no separate internal-energy argument is needed: c^2 = dp_ref - Gamma*rho*de_ref + ((1 + Gamma)*p -
+    !! p_ref)/rho, where dp_ref and de_ref are the reference-curve density derivatives dp_ref/drho and de_ref/drho. The result can
+    !! be negative outside the admissible domain, so the caller tests it before taking a square root.
+    pure function f_mg_sound_speed_sq(gamma_mg, rho, pres, p_ref, dp_ref, de_ref) result(c_sq)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+
+        real(wp), intent(in) :: gamma_mg, rho, pres, p_ref, dp_ref, de_ref
+        real(wp)             :: c_sq
+
+        c_sq = dp_ref - gamma_mg*rho*de_ref + ((1._wp + gamma_mg)*pres - p_ref)/rho
+
+    end function f_mg_sound_speed_sq
+
+end module m_eos_mie_gruneisen

From 0da82484966cfc88150e57bd93ebd9f36d56164f Mon Sep 17 00:00:00 2001
From: Fahad Nabid <55717325+fahnab666@users.noreply.github.com>
Date: Tue, 21 Jul 2026 03:02:54 -0400
Subject: [PATCH 2/3] Condense doc comments to the why-only essentials

---
 misc/mg_eos_verification.py        | 18 ++++++----------
 src/common/m_eos_mie_gruneisen.fpp | 33 ++++++++++--------------------
 2 files changed, 17 insertions(+), 34 deletions(-)

diff --git a/misc/mg_eos_verification.py b/misc/mg_eos_verification.py
index 227bbd1759..8e28637a53 100644
--- a/misc/mg_eos_verification.py
+++ b/misc/mg_eos_verification.py
@@ -1,16 +1,10 @@
 #!/usr/bin/env python3
-"""Manufactured verification for the generic Mie-Gruneisen EOS core in
-src/common/m_eos_mie_gruneisen.fpp.
-
-The Mie-Gruneisen backend is not yet selectable at runtime, so it has no golden
-regression case; this standalone check is its verification instead. It mirrors the
-Fortran formulas line for line and exercises the two properties that fix the generic
-assembler: exact reduction to the existing stiffened-gas EOS, and agreement of the
-analytic sound speed with a finite-difference of the frozen identity for a density
-dependent reference curve.
-
-Reference: Arienti, Morano and Shepherd, GALCIT FM99-8 (2004); Menikoff and Plohr,
-Rev. Mod. Phys. 61, 75 (1989). Run with python3; exits nonzero on any failure.
+"""Manufactured verification for src/common/m_eos_mie_gruneisen.fpp.
+
+The Mie-Gruneisen backend is not yet runtime-selectable, so it has no golden case; this
+script mirrors the Fortran formulas and checks exact reduction to stiffened gas and the
+analytic sound speed against a finite difference of the frozen identity (Arienti et al.
+2004; Menikoff and Plohr 1989). Run with python3; exits nonzero on failure.
 """
 
 import sys
diff --git a/src/common/m_eos_mie_gruneisen.fpp b/src/common/m_eos_mie_gruneisen.fpp
index 30865d7ced..70f0203b83 100644
--- a/src/common/m_eos_mie_gruneisen.fpp
+++ b/src/common/m_eos_mie_gruneisen.fpp
@@ -4,15 +4,10 @@
 
 #:include 'macros.fpp'
 
-!> @brief Generic Mie-Gruneisen equation of state with a constant Gruneisen coefficient, structured as a reference-curve provider
-!! plus a curve-agnostic assembler. A Mie-Gruneisen material obeys p = p_ref(rho) + Gamma*rho*(e - e_ref(rho)); the reference pair
-!! (p_ref, e_ref) and its density derivatives are the only thing that distinguishes one Mie-Gruneisen material from another.
-!! Following the production hydrocode structure (Arienti, Morano and Shepherd, GALCIT FM99-8, 2004; LLNL singularity-eos), one
-!! provider returns the reference curve for a chosen material and one assembler consumes it, so stiffened gas, ideal gas, and JWL
-!! differ only in the provider. This PR ships the stiffened-gas reference curve, which recovers the existing stiffened-gas pressure
-!! and sound speed exactly (Gamma = gamma_s - 1, p_ref = -pi_inf, e_ref = pi_inf/rho); the JWL principal isentrope is the next
-!! provider, added with its consumer. This is a leaf module (it uses only the precision kinds) and no solver path calls it yet: it
-!! is the verified thermodynamic core for the Mie-Gruneisen backend, landed ahead of the reference curve that consumes it.
+!> @brief Generic Mie-Gruneisen equation of state, p = p_ref(rho) + Gamma*rho*(e - e_ref(rho)) with constant Gamma, split into a
+!! reference-curve provider and a curve-agnostic assembler (Arienti et al., GALCIT FM99-8, 2004). Materials differ only in the
+!! provider: stiffened gas ships here and the JWL principal isentrope follows with its consumer. Leaf module; no solver path calls
+!! it yet.
 module m_eos_mie_gruneisen
 
     use m_precision_select
@@ -24,10 +19,8 @@ module m_eos_mie_gruneisen
 
 contains
 
-    !> Stiffened gas expressed as a Mie-Gruneisen reference curve. Returns the Gruneisen coefficient, the reference pressure and
-    !! energy, and their density derivatives, which together make the generic assembler reproduce stiffened gas. This is the first
-    !! concrete reference curve and the template the JWL principal isentrope follows. The reference pressure is density-independent
-    !! for stiffened gas, so its derivative is zero; ideal gas is the pi_inf = 0 sub-case.
+    !> Stiffened gas as a Mie-Gruneisen reference curve: Gamma = gamma_s - 1, p_ref = -pi_inf, e_ref = pi_inf/rho, plus the density
+    !! derivatives. Template for the JWL isentrope; ideal gas is the pi_inf = 0 sub-case.
     pure subroutine s_mg_stiffened_reference(gamma_s, pi_inf, rho, gamma_mg, p_ref, dp_ref, e_ref, de_ref)
 
         $:GPU_ROUTINE(parallelism='[seq]')
@@ -43,8 +36,7 @@ contains
 
     end subroutine s_mg_stiffened_reference
 
-    !> Mie-Gruneisen pressure from density and specific internal energy against a reference curve: p = p_ref + Gamma*rho*(e -
-    !! e_ref).
+    !> p = p_ref + Gamma*rho*(e - e_ref).
     pure function f_mg_pressure(gamma_mg, rho, e_int, p_ref, e_ref) result(pres)
 
         $:GPU_ROUTINE(parallelism='[seq]')
@@ -56,8 +48,7 @@ contains
 
     end function f_mg_pressure
 
-    !> Specific internal energy from pressure, the inverse of f_mg_pressure: e = e_ref + (p - p_ref)/(Gamma*rho). The Gamma*rho
-    !! divisor is positive for every admissible state (Gamma > 0, rho > 0).
+    !> Inverse of f_mg_pressure: e = e_ref + (p - p_ref)/(Gamma*rho).
     pure function f_mg_internal_energy(gamma_mg, rho, pres, p_ref, e_ref) result(e_int)
 
         $:GPU_ROUTINE(parallelism='[seq]')
@@ -69,11 +60,9 @@ contains
 
     end function f_mg_internal_energy
 
-    !> Squared sound speed c^2 = dp/drho at constant entropy for a constant Gruneisen coefficient, from the frozen identity c^2 =
-    !! dp/drho|_e + (p/rho^2) dp/de|_rho with dp/de|_rho = Gamma*rho (Arienti et al. 2004, Eq. 12; Menikoff and Plohr 1989). Written
-    !! in terms of pressure so no separate internal-energy argument is needed: c^2 = dp_ref - Gamma*rho*de_ref + ((1 + Gamma)*p -
-    !! p_ref)/rho, where dp_ref and de_ref are the reference-curve density derivatives dp_ref/drho and de_ref/drho. The result can
-    !! be negative outside the admissible domain, so the caller tests it before taking a square root.
+    !> Squared sound speed from the frozen identity c^2 = dp/drho|_e + (p/rho^2)*dp/de|_rho (Menikoff and Plohr 1989): c^2 = dp_ref
+    !! - Gamma*rho*de_ref + ((1 + Gamma)*p - p_ref)/rho. Can be negative outside the admissible domain, so the caller tests it
+    !! before taking a square root.
     pure function f_mg_sound_speed_sq(gamma_mg, rho, pres, p_ref, dp_ref, de_ref) result(c_sq)
 
         $:GPU_ROUTINE(parallelism='[seq]')

From c6a5882233770f0e4000052ecb6e2c8155f6a6d0 Mon Sep 17 00:00:00 2001
From: Fahad Nabid <55717325+fahnab666@users.noreply.github.com>
Date: Tue, 21 Jul 2026 03:48:04 -0400
Subject: [PATCH 3/3] Expose thermodynamic derivatives and move manufactured
 tests into the CI-run toolchain suite

---
 misc/mg_eos_verification.py        | 113 ------------------------
 src/common/m_eos_mie_gruneisen.fpp |  35 ++++++--
 toolchain/mfc/test_mg_eos.py       | 134 +++++++++++++++++++++++++++++
 3 files changed, 164 insertions(+), 118 deletions(-)
 delete mode 100644 misc/mg_eos_verification.py
 create mode 100644 toolchain/mfc/test_mg_eos.py

diff --git a/misc/mg_eos_verification.py b/misc/mg_eos_verification.py
deleted file mode 100644
index 8e28637a53..0000000000
--- a/misc/mg_eos_verification.py
+++ /dev/null
@@ -1,113 +0,0 @@
-#!/usr/bin/env python3
-"""Manufactured verification for src/common/m_eos_mie_gruneisen.fpp.
-
-The Mie-Gruneisen backend is not yet runtime-selectable, so it has no golden case; this
-script mirrors the Fortran formulas and checks exact reduction to stiffened gas and the
-analytic sound speed against a finite difference of the frozen identity (Arienti et al.
-2004; Menikoff and Plohr 1989). Run with python3; exits nonzero on failure.
-"""
-
-import sys
-
-
-# --- generic Mie-Gruneisen assembler (mirrors m_eos_mie_gruneisen.fpp) --------------
-def mg_pressure(gamma_mg, rho, e_int, p_ref, e_ref):
-    return p_ref + gamma_mg * rho * (e_int - e_ref)
-
-
-def mg_internal_energy(gamma_mg, rho, pres, p_ref, e_ref):
-    return e_ref + (pres - p_ref) / (gamma_mg * rho)
-
-
-def mg_sound_speed_sq(gamma_mg, rho, pres, p_ref, dp_ref, de_ref):
-    return dp_ref - gamma_mg * rho * de_ref + ((1.0 + gamma_mg) * pres - p_ref) / rho
-
-
-def stiffened_reference(gamma_s, pi_inf, rho):
-    """Stiffened gas as a Mie-Gruneisen reference curve."""
-    return gamma_s - 1.0, -pi_inf, 0.0, pi_inf / rho, -pi_inf / rho**2
-
-
-# --- analytic stiffened-gas oracle --------------------------------------------------
-def sg_pressure(gamma_s, rho, e_int, pi_inf):
-    return (gamma_s - 1.0) * rho * e_int - gamma_s * pi_inf
-
-
-def sg_sound_speed_sq(gamma_s, rho, pres, pi_inf):
-    return gamma_s * (pres + pi_inf) / rho
-
-
-def rel(a, b):
-    return abs(a - b) / max(1.0, abs(b))
-
-
-# --- Test A: exact reduction to stiffened gas (ideal gas is the pi_inf = 0 sub-case) -
-def test_reduction():
-    # (gamma_s, pi_inf, rho, e): ideal air, water (stiffened), a stiff solid, monatomic
-    cases = [
-        (1.4, 0.0, 1.2, 2.5e5),
-        (4.4, 6.0e8, 1000.0, 3.0e5),
-        (3.0, 1.0e5, 50.0, 2.0e4),
-        (1.667, 0.0, 0.1, 2.0e6),
-    ]
-    worst = 0.0
-    ok = True
-    for gamma_s, pi_inf, rho, e in cases:
-        gamma_mg, p_ref, dp_ref, e_ref, de_ref = stiffened_reference(gamma_s, pi_inf, rho)
-        p_mg = mg_pressure(gamma_mg, rho, e, p_ref, e_ref)
-        e_back = mg_internal_energy(gamma_mg, rho, p_mg, p_ref, e_ref)
-        c2_mg = mg_sound_speed_sq(gamma_mg, rho, p_mg, p_ref, dp_ref, de_ref)
-        errs = (
-            rel(p_mg, sg_pressure(gamma_s, rho, e, pi_inf)),
-            rel(e_back, e),
-            rel(c2_mg, sg_sound_speed_sq(gamma_s, rho, p_mg, pi_inf)),
-        )
-        worst = max(worst, *errs)
-        ok = ok and all(x < 1.0e-13 for x in errs)
-    print(f"[A] stiffened-gas reduction: worst rel err = {worst:.2e}  ->  {'PASS' if ok else 'FAIL'}")
-    return ok
-
-
-# --- Test B: analytic c^2 vs finite difference for a JWL-form isentrope reference ----
-# Synthetic (non-calibrated) products isentrope, exercising a density-dependent p_ref.
-def jwl_reference(rho, A, B, R1, R2, omega, rho0):
-    V = rho0 / rho
-    e1, e2 = A * pow(2.718281828459045, -R1 * V), B * pow(2.718281828459045, -R2 * V)
-    p_ref = e1 + e2
-    e_ref = e1 / (rho0 * R1) + e2 / (rho0 * R2)
-    dp_ref = (V / rho) * (R1 * e1 + R2 * e2)
-    de_ref = p_ref / rho**2
-    return omega, p_ref, dp_ref, e_ref, de_ref
-
-
-def test_soundspeed_fd():
-    A, B, R1, R2, omega, rho0 = 5.0e11, 8.0e9, 4.5, 1.2, 0.3, 1600.0
-    worst = 0.0
-    ok = True
-    for rho in (900.0, 1600.0, 2400.0):
-        for e in (1.0e6, 4.0e6):
-
-            def p_of(r, en):
-                g, p_r, _, e_r, _ = jwl_reference(r, A, B, R1, R2, omega, rho0)
-                return mg_pressure(g, r, en, p_r, e_r)
-
-            g, p_r, dp_r, e_r, de_r = jwl_reference(rho, A, B, R1, R2, omega, rho0)
-            pres = mg_pressure(g, rho, e, p_r, e_r)
-            c2_analytic = mg_sound_speed_sq(g, rho, pres, p_r, dp_r, de_r)
-
-            # frozen identity c^2 = dp/drho|_e + (p/rho^2) dp/de|_rho, dp/de|_rho = omega*rho
-            h = rho * 1.0e-6
-            dpdrho_e = (p_of(rho + h, e) - p_of(rho - h, e)) / (2.0 * h)
-            c2_fd = dpdrho_e + (pres / rho**2) * (omega * rho)
-
-            err = rel(c2_analytic, c2_fd)
-            worst = max(worst, err)
-            ok = ok and err < 1.0e-6
-    print(f"[B] sound speed vs finite difference: worst rel err = {worst:.2e}  ->  {'PASS' if ok else 'FAIL'}")
-    return ok
-
-
-if __name__ == "__main__":
-    passed = test_reduction() & test_soundspeed_fd()
-    print("ALL PASS" if passed else "FAILURE")
-    sys.exit(0 if passed else 1)
diff --git a/src/common/m_eos_mie_gruneisen.fpp b/src/common/m_eos_mie_gruneisen.fpp
index 70f0203b83..24963c2cb8 100644
--- a/src/common/m_eos_mie_gruneisen.fpp
+++ b/src/common/m_eos_mie_gruneisen.fpp
@@ -15,7 +15,7 @@ module m_eos_mie_gruneisen
     implicit none
 
     private
-    public :: s_mg_stiffened_reference, f_mg_pressure, f_mg_internal_energy, f_mg_sound_speed_sq
+    public :: s_mg_stiffened_reference, f_mg_pressure, f_mg_internal_energy, f_mg_dp_drho, f_mg_dp_de, f_mg_sound_speed_sq
 
 contains
 
@@ -60,9 +60,34 @@ contains
 
     end function f_mg_internal_energy
 
-    !> Squared sound speed from the frozen identity c^2 = dp/drho|_e + (p/rho^2)*dp/de|_rho (Menikoff and Plohr 1989): c^2 = dp_ref
-    !! - Gamma*rho*de_ref + ((1 + Gamma)*p - p_ref)/rho. Can be negative outside the admissible domain, so the caller tests it
-    !! before taking a square root.
+    !> Thermodynamic derivative dp/drho at constant specific internal energy, in terms of pressure: dp_ref + (p - p_ref)/rho -
+    !! Gamma*rho*de_ref, where dp_ref and de_ref are the reference-curve density derivatives.
+    pure function f_mg_dp_drho(gamma_mg, rho, pres, p_ref, dp_ref, de_ref) result(dp_drho)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+
+        real(wp), intent(in) :: gamma_mg, rho, pres, p_ref, dp_ref, de_ref
+        real(wp)             :: dp_drho
+
+        dp_drho = dp_ref + (pres - p_ref)/rho - gamma_mg*rho*de_ref
+
+    end function f_mg_dp_drho
+
+    !> Thermodynamic derivative dp/de at constant density: Gamma*rho, the defining property of the Gruneisen coefficient.
+    pure function f_mg_dp_de(gamma_mg, rho) result(dp_de)
+
+        $:GPU_ROUTINE(parallelism='[seq]')
+
+        real(wp), intent(in) :: gamma_mg, rho
+        real(wp)             :: dp_de
+
+        dp_de = gamma_mg*rho
+
+    end function f_mg_dp_de
+
+    !> Squared sound speed, assembled from the exposed derivatives per the frozen identity c^2 = dp/drho|_e + (p/rho^2)*dp/de|_rho
+    !! (Arienti et al. 2004; Menikoff and Plohr 1989). Can be negative outside the admissible domain, so the caller tests it before
+    !! taking a square root.
     pure function f_mg_sound_speed_sq(gamma_mg, rho, pres, p_ref, dp_ref, de_ref) result(c_sq)
 
         $:GPU_ROUTINE(parallelism='[seq]')
@@ -70,7 +95,7 @@ contains
         real(wp), intent(in) :: gamma_mg, rho, pres, p_ref, dp_ref, de_ref
         real(wp)             :: c_sq
 
-        c_sq = dp_ref - gamma_mg*rho*de_ref + ((1._wp + gamma_mg)*pres - p_ref)/rho
+        c_sq = f_mg_dp_drho(gamma_mg, rho, pres, p_ref, dp_ref, de_ref) + pres/rho**2*f_mg_dp_de(gamma_mg, rho)
 
     end function f_mg_sound_speed_sq
 
diff --git a/toolchain/mfc/test_mg_eos.py b/toolchain/mfc/test_mg_eos.py
new file mode 100644
index 0000000000..a87521a90a
--- /dev/null
+++ b/toolchain/mfc/test_mg_eos.py
@@ -0,0 +1,134 @@
+"""Manufactured verification for src/common/m_eos_mie_gruneisen.fpp.
+
+The Mie-Gruneisen backend is not yet runtime-selectable, so it has no golden case; these
+tests mirror the Fortran formulas line for line and run in CI with the toolchain unit
+tests. They check exact reduction to stiffened gas, the analytic sound speed and
+thermodynamic derivatives against finite differences of the frozen identity, and the
+isentrope self-consistency of the reference curves (Arienti et al., GALCIT FM99-8, 2004;
+Menikoff and Plohr, Rev. Mod. Phys. 61, 75, 1989).
+"""
+
+from math import exp
+
+
+# Generic Mie-Gruneisen assembler, mirroring m_eos_mie_gruneisen.fpp
+def mg_pressure(gamma_mg, rho, e_int, p_ref, e_ref):
+    return p_ref + gamma_mg * rho * (e_int - e_ref)
+
+
+def mg_internal_energy(gamma_mg, rho, pres, p_ref, e_ref):
+    return e_ref + (pres - p_ref) / (gamma_mg * rho)
+
+
+def mg_dp_drho(gamma_mg, rho, pres, p_ref, dp_ref, de_ref):
+    return dp_ref + (pres - p_ref) / rho - gamma_mg * rho * de_ref
+
+
+def mg_dp_de(gamma_mg, rho):
+    return gamma_mg * rho
+
+
+def mg_sound_speed_sq(gamma_mg, rho, pres, p_ref, dp_ref, de_ref):
+    return mg_dp_drho(gamma_mg, rho, pres, p_ref, dp_ref, de_ref) + pres / rho**2 * mg_dp_de(gamma_mg, rho)
+
+
+def stiffened_reference(gamma_s, pi_inf, rho):
+    """Stiffened gas as a Mie-Gruneisen reference curve."""
+    return gamma_s - 1.0, -pi_inf, 0.0, pi_inf / rho, -pi_inf / rho**2
+
+
+def jwl_reference(rho, A, B, R1, R2, omega, rho0):
+    """Synthetic (non-calibrated) JWL-form principal isentrope, exercising a density
+    dependent reference curve."""
+    V = rho0 / rho
+    e1, e2 = A * exp(-R1 * V), B * exp(-R2 * V)
+    p_ref = e1 + e2
+    e_ref = e1 / (rho0 * R1) + e2 / (rho0 * R2)
+    dp_ref = (V / rho) * (R1 * e1 + R2 * e2)
+    de_ref = p_ref / rho**2
+    return omega, p_ref, dp_ref, e_ref, de_ref
+
+
+# ideal air, water (stiffened), a stiff solid, monatomic: (gamma_s, pi_inf, rho, e)
+SG_CASES = [
+    (1.4, 0.0, 1.2, 2.5e5),
+    (4.4, 6.0e8, 1000.0, 3.0e5),
+    (3.0, 1.0e5, 50.0, 2.0e4),
+    (1.667, 0.0, 0.1, 2.0e6),
+]
+
+JWL_PARAMS = (5.0e11, 8.0e9, 4.5, 1.2, 0.3, 1600.0)
+JWL_STATES = [(rho, e) for rho in (900.0, 1600.0, 2400.0) for e in (1.0e6, 4.0e6)]
+
+
+def rel(a, b):
+    return abs(a - b) / max(1.0, abs(b))
+
+
+def test_stiffened_gas_reduction():
+    """The generic assembler fed the stiffened reference curve must reproduce the
+    analytic stiffened-gas pressure, energy round-trip, and sound speed exactly."""
+    for gamma_s, pi_inf, rho, e in SG_CASES:
+        gamma_mg, p_ref, dp_ref, e_ref, de_ref = stiffened_reference(gamma_s, pi_inf, rho)
+        p_mg = mg_pressure(gamma_mg, rho, e, p_ref, e_ref)
+        assert rel(p_mg, (gamma_s - 1.0) * rho * e - gamma_s * pi_inf) < 1.0e-13
+        assert rel(mg_internal_energy(gamma_mg, rho, p_mg, p_ref, e_ref), e) < 1.0e-13
+        c2 = mg_sound_speed_sq(gamma_mg, rho, p_mg, p_ref, dp_ref, de_ref)
+        assert rel(c2, gamma_s * (p_mg + pi_inf) / rho) < 1.0e-13
+
+
+def test_derivatives_and_sound_speed_vs_finite_difference():
+    """Analytic dp/drho|_e, dp/de|_rho, and c^2 against central finite differences of
+    p(rho, e) for the density-dependent JWL-form reference curve."""
+    A, B, R1, R2, omega, rho0 = JWL_PARAMS
+
+    def p_of(r, en):
+        g, p_r, _, e_r, _ = jwl_reference(r, A, B, R1, R2, omega, rho0)
+        return mg_pressure(g, r, en, p_r, e_r)
+
+    for rho, e in JWL_STATES:
+        g, p_ref, dp_ref, e_ref, de_ref = jwl_reference(rho, A, B, R1, R2, omega, rho0)
+        pres = mg_pressure(g, rho, e, p_ref, e_ref)
+
+        h = rho * 1.0e-6
+        fd_dp_drho = (p_of(rho + h, e) - p_of(rho - h, e)) / (2.0 * h)
+        assert rel(mg_dp_drho(g, rho, pres, p_ref, dp_ref, de_ref), fd_dp_drho) < 1.0e-6
+
+        he = e * 1.0e-6
+        fd_dp_de = (p_of(rho, e + he) - p_of(rho, e - he)) / (2.0 * he)
+        assert rel(mg_dp_de(g, rho), fd_dp_de) < 1.0e-9
+
+        c2_fd = fd_dp_drho + (pres / rho**2) * mg_dp_de(g, rho)
+        assert rel(mg_sound_speed_sq(g, rho, pres, p_ref, dp_ref, de_ref), c2_fd) < 1.0e-6
+
+
+def test_isentrope_self_consistency():
+    """Both reference curves must satisfy de_ref/drho = p_ref/rho^2, the density form
+    of de = -p dv along an isentrope; this is what makes them valid MG references and
+    catches per-mass vs per-initial-volume scaling errors in e_ref."""
+    A, B, R1, R2, omega, rho0 = JWL_PARAMS
+    for rho in (900.0, 1600.0, 2400.0):
+        _, p_ref, _, e_ref, de_ref = jwl_reference(rho, A, B, R1, R2, omega, rho0)
+        assert rel(de_ref, p_ref / rho**2) < 1.0e-13
+        h = rho * 1.0e-6
+        _, _, _, e_hi, _ = jwl_reference(rho + h, A, B, R1, R2, omega, rho0)
+        _, _, _, e_lo, _ = jwl_reference(rho - h, A, B, R1, R2, omega, rho0)
+        assert rel((e_hi - e_lo) / (2.0 * h), de_ref) < 1.0e-6
+    for gamma_s, pi_inf, rho, _ in SG_CASES:
+        _, p_ref, _, _, de_ref = stiffened_reference(gamma_s, pi_inf, rho)
+        assert rel(de_ref, p_ref / rho**2) < 1.0e-13
+
+
+def test_admissibility_scan():
+    """c^2 > 0 over the synthetic envelope for both reference curves."""
+    A, B, R1, R2, omega, rho0 = JWL_PARAMS
+    for rho, e in JWL_STATES:
+        g, p_ref, dp_ref, e_ref, de_ref = jwl_reference(rho, A, B, R1, R2, omega, rho0)
+        pres = mg_pressure(g, rho, e, p_ref, e_ref)
+        assert mg_sound_speed_sq(g, rho, pres, p_ref, dp_ref, de_ref) > 0.0
+    for gamma_s, pi_inf, rho, e in SG_CASES:
+        g, p_ref, dp_ref, e_ref, de_ref = stiffened_reference(gamma_s, pi_inf, rho)
+        # admissible states carry thermal energy on top of the reference energy;
+        # below e_ref the stiffened gas is in tension past -pi_inf and c^2 < 0 is correct
+        pres = mg_pressure(g, rho, e_ref + e, p_ref, e_ref)
+        assert mg_sound_speed_sq(g, rho, pres, p_ref, dp_ref, de_ref) > 0.0