Add a per-fluid equation-of-state selector

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

@@ -114,6 +114,13 @@ module m_constants
     integer, parameter :: BC_NO_SLIP_WALL = -16
     integer, parameter :: BC_DIRICHLET = -17
 
+    ! fluid_pp(:)%eos selector; values must match _EOS_NAMES in toolchain/mfc/params/definitions.py
+    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
+
     ! Synthetic turbulence array size limits
     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

@@ -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,27 @@ 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"
+
+        if num_fluids is None:
+            return
+
+        eos_stiffened_gas, eos_ideal_gas_mixture = 1, 2
+
+        for i in range(1, num_fluids + 1):
+            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 +2340,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": {
@@ -887,6 +895,7 @@ def _load():
         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