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

@@ -463,6 +463,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. Only `stiffened_gas` (the default) and `ideal_gas_mixture` (requires a chemistry build, backed by Pyrometheus) are currently supported; the enumeration reserves `mie_gruneisen`, `jwl`, and `table` for future backends. 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,29 @@ 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 unimplemented 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
+
+        do i = 1, min(num_fluids + merge(1, 0, bubbles_euler), 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,13 @@ 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).
+    integer, parameter :: eos_stiffened_gas = 1
+    integer, parameter :: eos_ideal_gas_mixture = 2
+    integer, parameter :: eos_mie_gruneisen = 3
+    integer, parameter :: eos_jwl = 4
+    integer, parameter :: eos_table = 5
     #: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, model_eqns_4eq, recon_type_weno, &
-        & recon_type_muscl, name_len, dflt_int, dflt_real
+        & recon_type_muscl, name_len, 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

@@ -427,6 +427,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,16 @@
         "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; unimplemented values and intra-cell EOS mixing are rejected."
+        ),
+        "references": ["Wilfong26"],
+    },
     "check_patch_physics": {
         "title": "Patch Initial Condition Constraints",
         "category": "Thermodynamic Constraints",
@@ -747,6 +757,32 @@ 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"""
+        num_fluids = self.get("num_fluids")
+        chemistry = self.get("chemistry", "F") == "T"
+        bubbles_euler = self.get("bubbles_euler", "F") == "T"
+
+        if num_fluids is None:
+            return
+
+        eos_names = CONSTRAINTS["fluid_pp(1)%eos"]["names"]
+        eos_stiffened_gas, eos_ideal_gas_mixture = eos_names["stiffened_gas"], eos_names["ideal_gas_mixture"]
+
+        # Allow one extra fluid property slot when using bubbles_euler
+        bub_fac = 1 if (bubbles_euler) else 0
+
+        for i in range(1, num_fluids + 1 + bub_fac):
+            eos = self.get(f"fluid_pp({i})%eos")
+            if eos is None:
+                continue
+            self.prohibit(
+                eos not in (eos_stiffened_gas, eos_ideal_gas_mixture),
+                f"fluid_pp({i})%eos selects an equation of state that is not yet implemented; " "only 'stiffened_gas' and 'ideal_gas_mixture' are available",
+            )
+            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"
@@ -2309,6 +2345,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", 3: "Mie-Grueneisen", 4: "JWL", 5: "tabulated"}
+_EOS_NAMES = {"stiffened_gas": 1, "ideal_gas_mixture": 2, "mie_gruneisen": 3, "jwl": 4, "table": 5}
+for _f in range(1, NF + 1):
+    CONSTRAINTS[f"fluid_pp({_f})%eos"] = {"choices": [1, 2, 3, 4, 5], "value_labels": _EOS_VALUE_LABELS, "names": _EOS_NAMES}
+
 # Parameter dependencies (requires, recommends)
 DEPENDENCIES = {
     "bubbles_euler": {
@@ -880,13 +888,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

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

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