Add a per-fluid equation-of-state selector

.claude/rules/common-pitfalls.md

@@ -48,7 +48,12 @@ covered in `docs/documentation/contributing.md`.
   `CASE_OPT_EXTRA_LINES` literal in `toolchain/mfc/params/generators/fortran_gen.py` (covers `num_dims`,
   `num_vels`, `weno_polyn`, `muscl_polyn`, `weno_num_stencils`, `wenojs`);
   multi-variable declaration lines (`bc_x/y/z`, `x/y/z_domain`, `x/y/z_output`, post's
-  `G`); and the MPI broadcast residue in `m_mpi_proxy` (computed variables that are not
+  `G`); enum constants for compound registry keys (a `CONSTRAINTS` key with a `names` dict whose
+  key contains `%` or `(`, e.g. `fluid_pp(:)%eos`): `generate_constants_fpp` **silently** skips
+  these (the `{param}_{name}` form, `fluid_pp(1)%eos_stiffened_gas`, is not a valid Fortran
+  identifier), so their constants must be hand-written in `m_constants.fpp` (as `eos_*` is) or
+  they simply never exist; and the
+  MPI broadcast residue in `m_mpi_proxy` (computed variables that are not
   namelist-bound: `m_glb`/`n_glb`/`p_glb`, `cfl_dt`, `bc_io`, and complex struct-member
   array loops — these cannot be auto-generated and stay hand-listed). Everything else — scalar declarations, plain arrays (`FORTRAN_ARRAY_DIMS`
   table in `definitions.py`), derived-type namelist declarations including `GPU_DECLARE`

docs/documentation/case.md

@@ -461,6 +461,8 @@ The parameters define material's property of compressible fluids that are used i
 
 - `fluid_pp(i)%%gamma` and `fluid_pp(i)%%pi_inf` define \f$\Gamma\f$ and \f$\Pi\f$ as parameters of $i$-th fluid that are used in stiffened gas equation of state.
 
+- `fluid_pp(i)%%eos` selects the equation of state of the $i$-th fluid. The accepted values are `stiffened_gas` (the default) and `ideal_gas_mixture` (requires a chemistry build, backed by Pyrometheus). Every fluid in a run must use the same family. For a non-chemistry ideal gas, use `stiffened_gas` with `pi_inf = 0`; `ideal_gas_mixture` is the Pyrometheus mixture backend and is only valid in a chemistry build.
+
 - `fluid_pp(i)%%Re(1)` and `fluid_pp(i)%%Re(2)` define the shear and volume viscosities of $i$-th fluid, respectively.
 
 When these parameters are undefined, fluids are treated as inviscid.

src/common/m_checker_common.fpp

@@ -26,12 +26,32 @@ contains
         integer(kind=8), intent(in) :: n_global
 
         if (check_total_cells) call s_check_total_cells(n_global)
+        call s_check_eos
         #:if USING_AMD
             call s_check_amd
         #:endif
 
     end subroutine s_check_inputs_common
 
+    !> Reject unsupported EOS selectors and intra-cell mixing; only stiffened_gas (non-chemistry) and ideal_gas_mixture (chemistry)
+    !! have a backend, and every fluid in a run must share one family.
+    impure subroutine s_check_eos
+
+        integer :: i
+
+        ! Every slot is default-assigned and broadcast up to num_fluids_max, so check all of
+        ! them: a selector left on an unused slot still reaches the solver. Input-time only.
+
+        do i = 1, num_fluids_max
+            @:PROHIBIT(chemistry .and. fluid_pp(i)%eos /= eos_ideal_gas_mixture, &
+                       & "fluid_pp(:)%eos must be 'ideal_gas_mixture' for every fluid when chemistry is enabled")
+            @:PROHIBIT(.not. chemistry .and. fluid_pp(i)%eos /= eos_stiffened_gas, &
+                       & "fluid_pp(:)%eos selector is not supported; only 'stiffened_gas' is available " &
+                       & // "(or 'ideal_gas_mixture' with a chemistry build)")
+        end do
+
+    end subroutine s_check_eos
+
     !> Verify that the total number of grid cells meets the minimum required by the number of dimensions and MPI ranks.
     impure subroutine s_check_total_cells(n_global)
 

src/common/m_constants.fpp

@@ -118,7 +118,12 @@ module m_constants
     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
 
-    ! Named values for enumerated case parameters (e.g. riemann_solver_hllc).
-    ! AUTO-GENERATED from "names" in toolchain/mfc/params/definitions.py.
+    ! Enum values are auto-generated from "names" in definitions.py by the include below, except
+    ! compound keys ("%" or "("), which generate_constants_fpp silently skips. So eos_* is
+    ! hand-written here and must match _EOS_NAMES in definitions.py (see common-pitfalls.md).
+    ! test_eos_selector.py::test_fortran_and_python_enums_agree guards the two against drift.
+    ! Only backends with a thermodynamics adapter belong here; add a value when its backend lands.
+    integer, parameter :: eos_stiffened_gas = 1
+    integer, parameter :: eos_ideal_gas_mixture = 2
     #:include 'generated_constants.fpp'
 end module m_constants

src/common/m_derived_types.fpp

@@ -382,6 +382,7 @@ 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
+        integer                :: eos            !< Equation of state selector (eos_* in m_constants)
         real(wp)               :: gamma          !< Sp. heat ratio
         real(wp)               :: pi_inf         !< Liquid stiffness
         real(wp), dimension(2) :: Re             !< Reynolds number

src/common/m_global_parameters_common.fpp

@@ -16,7 +16,7 @@ module m_global_parameters_common
     use m_derived_types
     use m_thermochem, only: num_species
     use m_constants, only: model_eqns_gamma_law, model_eqns_5eq, model_eqns_6eq, recon_type_weno, recon_type_muscl, name_len, &
-        & dflt_int, dflt_real
+        & dflt_int, dflt_real, eos_stiffened_gas, eos_ideal_gas_mixture
 
     implicit none
 

src/post_process/m_global_parameters.fpp

@@ -207,6 +207,7 @@ contains
             fluid_pp(i)%cv = 0._wp
             fluid_pp(i)%qv = 0._wp
             fluid_pp(i)%qvp = 0._wp
+            fluid_pp(i)%eos = merge(eos_ideal_gas_mixture, eos_stiffened_gas, chemistry)
             fluid_pp(i)%G = dflt_real
             fluid_pp(i)%non_newtonian = .false.
             fluid_pp(i)%K = dflt_real

src/pre_process/m_global_parameters.fpp

@@ -391,6 +391,7 @@ contains
             fluid_pp(i)%cv = 0._wp
             fluid_pp(i)%qv = 0._wp
             fluid_pp(i)%qvp = 0._wp
+            fluid_pp(i)%eos = merge(eos_ideal_gas_mixture, eos_stiffened_gas, chemistry)
             fluid_pp(i)%G = 0._wp
             fluid_pp(i)%non_newtonian = .false.
             fluid_pp(i)%K = dflt_real

src/simulation/m_global_parameters.fpp

@@ -450,6 +450,7 @@ contains
             fluid_pp(i)%cv = 0._wp
             fluid_pp(i)%qv = 0._wp
             fluid_pp(i)%qvp = 0._wp
+            fluid_pp(i)%eos = merge(eos_ideal_gas_mixture, eos_stiffened_gas, chemistry)
             fluid_pp(i)%Re(:) = dflt_real
             fluid_pp(i)%G = 0._wp
             fluid_pp(i)%non_newtonian = .false.

toolchain/mfc/case_validator.py

@@ -41,6 +41,17 @@
         "explanation": ("MFC uses the transformed stiffened gas parameter. A common mistake is entering the physical gamma (e.g., 1.4 for air) instead of the transformed value 1/(gamma-1) = 2.5."),
         "references": ["Wilfong26", "Allaire02"],
     },
+    "check_eos": {
+        "title": "Equation of State Selection",
+        "category": "Thermodynamic Constraints",
+        "explanation": (
+            "The per-fluid eos selector exposes only the backends with a thermodynamics adapter: "
+            "'stiffened_gas' (default) and 'ideal_gas_mixture' (chemistry, Pyrometheus). A single run "
+            "uses one family for every fluid, so intra-cell EOS mixing is rejected, as are values "
+            "outside the enumeration."
+        ),
+        "references": ["Wilfong26"],
+    },
     "check_patch_physics": {
         "title": "Patch Initial Condition Constraints",
         "category": "Thermodynamic Constraints",
@@ -782,6 +793,28 @@ def check_stiffened_eos(self):
                 self.prohibit(gamma is not None, f"model_eqns = 1 does not support fluid_pp({i})%gamma")
                 self.prohibit(pi_inf is not None, f"model_eqns = 1 does not support fluid_pp({i})%pi_inf")
 
+    def check_eos(self):
+        """Restricts the per-fluid EOS selector to the currently supported adapters"""
+        chemistry = self.get("chemistry", "F") == "T"
+
+        eos_names = CONSTRAINTS["fluid_pp(1)%eos"]["names"]
+        eos_ideal_gas_mixture = eos_names["ideal_gas_mixture"]
+        eos_values = set(eos_names.values())
+
+        # Every fluid_pp slot is default-assigned and MPI-broadcast up to num_fluids_max, so
+        # validate all of them rather than stopping at num_fluids: a selector left on an unused
+        # slot still reaches the solver. This is input-time only, so the wider loop costs nothing.
+        num_fluids_max = get_fortran_constants().get("num_fluids_max", 10)
+        for i in range(1, num_fluids_max + 1):
+            eos = self.get(f"fluid_pp({i})%eos")
+            if eos is None:
+                continue
+            # The "choices" constraint is enforced by validate_constraints, a separate layer from
+            # CaseValidator, so membership is re-checked here rather than assumed.
+            self.prohibit(eos not in eos_values, f"fluid_pp({i})%eos must be 'stiffened_gas' or 'ideal_gas_mixture'")
+            self.prohibit(chemistry and eos != eos_ideal_gas_mixture, f"fluid_pp({i})%eos must be 'ideal_gas_mixture' when chemistry is enabled")
+            self.prohibit(not chemistry and eos == eos_ideal_gas_mixture, f"fluid_pp({i})%eos = 'ideal_gas_mixture' requires a chemistry build")
+
     def check_surface_tension(self):
         """Checks constraints on surface tension"""
         surface_tension = self.get("surface_tension", "F") == "T"
@@ -2356,6 +2389,7 @@ def validate_common(self):
         self.check_phase_change()
         self.check_ibm()
         self.check_stiffened_eos()
+        self.check_eos()
         self.check_eos_parameter_sanity()
         self.check_surface_tension()
         self.check_mhd()

toolchain/mfc/params/definitions.py

@@ -112,6 +112,7 @@ def _fc(name: str, default: int) -> int:
         "cv": "Specific heat at constant volume",
         "qv": "Heat of formation",
         "qvp": "Heat of formation derivative",
+        "eos": "Equation of state selector",
     },
 }
 
@@ -386,6 +387,13 @@ def get_value_label(param_name: str, value: int) -> str:
     "p": {"min": 0},
 }
 
+# Values must match the hand-written eos_* constants in src/common/m_constants.fpp;
+# generate_constants_fpp skips compound keys, so these entries only drive name resolution and validation.
+_EOS_VALUE_LABELS = {1: "stiffened-gas", 2: "ideal-gas mixture"}
+_EOS_NAMES = {"stiffened_gas": 1, "ideal_gas_mixture": 2}
+for _f in range(1, NF + 1):
+    CONSTRAINTS[f"fluid_pp({_f})%eos"] = {"choices": [1, 2], "value_labels": _EOS_VALUE_LABELS, "names": _EOS_NAMES}
+
 # Parameter dependencies (requires, recommends)
 DEPENDENCIES = {
     "bubbles_euler": {
@@ -904,13 +912,15 @@ def _load():
                 _r(f"{px}sph_har_coeff({ll},{mm})", REAL)
 
     # fluid_pp (10 fluids)
-    # Members present in physical_parameters: gamma, pi_inf, Re, cv, qv, qvp, G.
+    # Members present in physical_parameters: gamma, pi_inf, Re, cv, qv, qvp, eos, G,
+    # non_newtonian, K, nn, tau0, hb_m, mu_min, mu_max, mu_bulk.
     # 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})%"
         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$")
         _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/params/descriptions.py

@@ -347,6 +347,7 @@
     (r"fluid_pp\((\d+)\)%cv", "Specific heat at constant volume for fluid {0}"),
     (r"fluid_pp\((\d+)\)%qv", "Heat of formation for fluid {0}"),
     (r"fluid_pp\((\d+)\)%qvp", "Heat of formation prime for fluid {0}"),
+    (r"fluid_pp\((\d+)\)%eos", "Equation of state selector for fluid {0}"),
     (r"fluid_pp\((\d+)\)%Re\((\d+)\)", "Reynolds number component {1} for fluid {0}"),
     (r"fluid_pp\((\d+)\)%non_newtonian", "Enable Herschel-Bulkley non-Newtonian viscosity for fluid {0}"),
     (r"fluid_pp\((\d+)\)%K", "HB consistency index for fluid {0}"),

toolchain/mfc/params/generators/fortran_gen.py

@@ -260,6 +260,9 @@ def generate_constants_fpp() -> str:
 
     lines = [_HEADER.rstrip()]
     for param in sorted(CONSTRAINTS):
+        # Compound keys (e.g. fluid_pp(1)%eos) are not valid Fortran identifiers; hand-written in m_constants.fpp
+        if "%" in param or "(" in param:
+            continue
         names = CONSTRAINTS[param].get("names")
         if not names:
             continue

toolchain/mfc/params_tests/test_eos_selector.py

@@ -0,0 +1,99 @@
+"""
+Tests for the per-fluid equation-of-state selector, fluid_pp(i)%eos.
+
+Covers the enum itself (Fortran/Python agreement), the readable-name to integer
+resolution done by Case, and the check_eos constraints in case_validator.
+"""
+
+import unittest
+
+from ..case import Case
+from ..case_validator import CaseConstraintError, CaseValidator
+from ..common import MFCException
+from ..params.definitions import _EOS_NAMES
+from ..params.namelist_parser import get_fortran_constants
+from .negative_tests import BASE_CASE
+
+
+def _eos_errors(overrides):
+    """Validate BASE_CASE plus overrides, returning only the eos-related messages."""
+    params = dict(BASE_CASE)
+    params.update(overrides)
+    try:
+        CaseValidator(Case(params).params).validate("pre_process")
+    except CaseConstraintError as exc:
+        return [line for line in str(exc).splitlines() if "%eos" in line]
+    return []
+
+
+class TestEosEnum(unittest.TestCase):
+    """The enum is hand-written in m_constants.fpp and restated in definitions.py."""
+
+    def test_fortran_and_python_enums_agree(self):
+        """_EOS_NAMES must match the eos_* parameters in m_constants.fpp.
+
+        generate_constants_fpp skips compound registry keys, so these constants are
+        hand-written on the Fortran side and nothing else forces the two to agree.
+        """
+        fortran = {name[len("eos_") :]: value for name, value in get_fortran_constants().items() if name.startswith("eos_")}
+        self.assertEqual(fortran, _EOS_NAMES)
+
+    def test_only_implemented_backends_are_exposed(self):
+        """Reserved values must not appear until they have a backend and a check_eos branch."""
+        self.assertEqual(set(_EOS_NAMES), {"stiffened_gas", "ideal_gas_mixture"})
+
+
+class TestEosNameResolution(unittest.TestCase):
+    """Case converts the readable name in a case file to the integer the namelist carries."""
+
+    def test_name_resolves_to_integer(self):
+        case = Case({"fluid_pp(1)%eos": "stiffened_gas"})
+        self.assertEqual(case.params["fluid_pp(1)%eos"], _EOS_NAMES["stiffened_gas"])
+
+    def test_integer_passes_through(self):
+        case = Case({"fluid_pp(1)%eos": _EOS_NAMES["ideal_gas_mixture"]})
+        self.assertEqual(case.params["fluid_pp(1)%eos"], _EOS_NAMES["ideal_gas_mixture"])
+
+    def test_unknown_name_rejected(self):
+        with self.assertRaises(MFCException) as ctx:
+            Case({"fluid_pp(1)%eos": "jwl"})
+        self.assertIn("stiffened_gas", str(ctx.exception))
+
+    def test_resolution_applies_to_every_fluid_slot(self):
+        """CONSTRAINTS is registered per slot, so slot 10 must resolve like slot 1."""
+        case = Case({"fluid_pp(10)%eos": "stiffened_gas"})
+        self.assertEqual(case.params["fluid_pp(10)%eos"], _EOS_NAMES["stiffened_gas"])
+
+
+class TestCheckEos(unittest.TestCase):
+    """check_eos constraints, on a non-chemistry build (BASE_CASE sets no chemistry)."""
+
+    def test_base_case_has_no_eos_errors(self):
+        self.assertEqual(_eos_errors({}), [])
+
+    def test_stiffened_gas_accepted(self):
+        self.assertEqual(_eos_errors({"fluid_pp(1)%eos": "stiffened_gas"}), [])
+
+    def test_ideal_gas_mixture_requires_chemistry(self):
+        errors = _eos_errors({"fluid_pp(1)%eos": "ideal_gas_mixture"})
+        self.assertTrue(errors)
+        self.assertIn("requires a chemistry build", " ".join(errors))
+
+    def test_value_outside_enum_rejected(self):
+        """The choices constraint covers integers the enum does not define."""
+        errors = _eos_errors({"fluid_pp(1)%eos": 99})
+        self.assertTrue(errors)
+
+    def test_unused_slot_is_validated(self):
+        """Slots above num_fluids are default-assigned and broadcast, so they are checked too.
+
+        BASE_CASE sets num_fluids = 1; fluid_pp(3) is an unused slot whose value still
+        reaches every rank through the fluid_pp member-loop broadcast.
+        """
+        errors = _eos_errors({"fluid_pp(3)%eos": "ideal_gas_mixture"})
+        self.assertTrue(errors)
+        self.assertIn("fluid_pp(3)%eos", " ".join(errors))
+
+
+if __name__ == "__main__":
+    unittest.main()