Replace Pyrometheus with MFC-owned thermochemistry generation

README.md

@@ -156,7 +156,7 @@ MFC is an SPEChpc benchmark candidate, part of the JSC JUPITER Early Access Prog
 * Ideal and stiffened gas equations of state
 * Body forces
 * Sound wave generation (one- and two-way sources)
-* Chemistry and multi-species transport via [Pyrometheus](https://github.com/pyrometheus/pyrometheus)
+* Chemistry and multi-species transport via MFC-owned Fortran generation from Cantera mechanisms
 	* Stiff time-integration support 
 * Magnetohydrodynamics (MHD)
 	* Includes non- and standard-relativistic (RMHD)

docs/documentation/case.md

@@ -1253,6 +1253,8 @@ When ``cyl_coord = 'T'`` is set in 2D the following constraints must be met:
 
 - `cantera_file` specifies the chemical mechanism file. If the file is part of the standard Cantera library, only the filename is required. Otherwise, the file must be located in the same directory as your `case.py` file
 
+MFC generates and compiles the mechanism's Fortran routines itself. Supported mechanism features and the Cantera-only mixing-layer initialization are described in @ref thermochemistry "Thermochemistry implementation".
+
 ### 18. Chemistry-Specific Boundary Conditions
 
 | Parameter          | Type    | Description                                                                 |

docs/documentation/equations.md

@@ -678,7 +678,7 @@ Enthalpy flux with diffusion:
 
 \f[q_\text{diff} = \lambda\,\frac{\partial T}{\partial x} + \sum_k h_k\,\dot{m}_k\f]
 
-Reaction mechanisms are code-generated via Pyrometheus (\cite Cisneros26), which provides symbolic abstractions for thermochemistry that enable portable GPU computation and automatic differentiation of chemical source terms.
+Reaction mechanisms are compiled into Fortran by MFC's own thermochemistry generator, using Cantera to load mechanism and transport data. The generator derives from the Fortran path of Pyrometheus (\cite Cisneros26), with its MIT attribution retained. The generated routines support CPU, OpenACC, and OpenMP execution; MFC does not require Pyrometheus, JAX, or automatic differentiation. See @ref thermochemistry "Thermochemistry implementation".
 
 ---
 

docs/documentation/thermochemistry.md

@@ -0,0 +1,98 @@
+@page thermochemistry Thermochemistry implementation
+
+# Mechanism generation
+
+MFC owns the Fortran thermochemistry generator in `toolchain/mfc/thermochem/`.
+Cantera loads a mechanism, and `generate_fortran` produces `m_thermochem.f90`
+in the target's build staging directory. The existing CMake build compiles that
+module into MFC. Neither Cantera nor Python is called inside the flow solver.
+Pyrometheus and JAX are not installation or runtime requirements.
+
+The implementation was adapted from the MIT-licensed Pyrometheus 1.1.1 Fortran
+generator. Its license is retained in `toolchain/mfc/thermochem/LICENSE`.
+Changes are made directly in MFC; no external fork or upstream generator release
+is required. Pymbolic constructs scalar expressions during generation and Mako
+renders the Fortran template. Neither is a runtime solver dependency.
+
+| File | Responsibility |
+|---|---|
+| `thermochem/expressions.py` | NASA7, reaction-rate, equilibrium and transport expressions |
+| `thermochem/fortran.py` | Supported-feature checks, expression formatting and module generation |
+| `thermochem/module.f90.mako` | Fortran interface and numerical routines |
+| `thermochem/fingerprint.py` | Mechanism and generator content identities for build reuse |
+| `run/input.py` | Mechanism resolution and generation for each target |
+
+The generator supports ideal-gas mixtures with NASA7 thermodynamics, elementary
+and third-body Arrhenius reactions, and Troe or Lindemann falloff reactions.
+It requires gas transport data and positive Arrhenius pre-exponential factors.
+Other rate types, custom reaction orders and other thermodynamic models are
+rejected before emission with a species or reaction identifier. Cantera's ability
+to parse a mechanism does not imply that the generator supports every feature in it.
+
+The generated module provides species metadata, caloric and ideal-gas properties,
+temperature inversion, net production rates, fused creation/destruction rates,
+mixture viscosity and thermal conductivity, and species diffusivities. MFC still
+owns reaction time integration, including alpha-QSS, and spatial transport
+discretization. The generator emits single- or double-precision routines and the
+selected CPU, OpenACC or OpenMP annotations. Nonchemistry builds retain the existing
+dummy `h2o2.yaml` module to satisfy the shared Fortran interfaces.
+
+The initial ownership change preserves the previous numerical formulas. In
+particular, the pure-species diffusion limit retains the self-diffusion coefficient;
+Cantera's mixture diffusion query may return zero for that degenerate state.
+Temperature inversion also retains the existing Newton iteration and tolerance.
+Real constants use the selected working precision: single for single-precision
+builds and double for default and mixed-storage builds. Single-precision results
+can therefore differ in roundoff from the former double-literal expressions.
+Troe falloff guards the logarithm at zero reduced pressure, including when a
+compiler evaluates both arguments of a Fortran `merge` expression.
+
+# Mixing-layer initial conditions
+
+The temporal 2D, spatial 2D and temporal 3D reacting mixing-layer examples share
+`toolchain/mfc/flamelet.py`. Their local `flamelet_ic.py` files retain grid generation,
+coordinate conventions, perturbations and MFC input-file layouts.
+
+The default cold profile uses the prescribed hyperbolic-tangent mixture fraction,
+linearly mixed stream mass fractions and specific enthalpy, and Cantera's HP
+temperature recovery. Velocity is interpolated between the prescribed stream
+velocities. Density follows the ideal-gas mixture equation of state.
+
+With `--hot`, Cantera solves a counterflow diffusion flame using unity-Lewis transport.
+The counterflow width is ten vorticity thicknesses. Each inlet velocity is
+`flame_strain_rate * width / 2`, with mass flux determined by its stream density;
+`flame_strain_rate` is a nominal inlet strain parameter in inverse seconds, not the
+computed local strain. Its default is 100/s and is set in each example's `case.py`.
+The solution is tabulated against Bilger mixture fraction normalized to the inlet
+compositions, then mapped to the mixing layer's prescribed tanh profile. Temperature
+is recovered at the linearly mixed stream enthalpy. Extinguished flames or invalid
+profiles fail explicitly instead of silently producing a cold initialization.
+
+This hot initialization is a counterflow-based seed for the evolving MFC flow. It
+replaces the former JAX flamelet solver and its iterative scalar-dissipation matching;
+it is not numerically identical to that initialization or a steady solution of the
+mixing-layer equations. Changes to the initialization method and mechanism invalidate
+the example caches. No precomputed flame data or optional Pyrometheus installation
+is required.
+
+# Validation
+
+After bootstrapping the toolchain, run:
+
+```sh
+PYTHONPATH=toolchain build/venv/bin/pytest -q \
+    toolchain/mfc/test_thermochem.py toolchain/mfc/test_flamelet.py
+./mfc.sh test --no-mpi -j 8 --only Chemistry
+```
+
+The kernel tests compile generated Fortran with GNU Fortran and compare with
+Cantera for `h2o2.yaml`, `gri30.yaml`, the bundled San Diego mechanism, and the
+hydrogen/xenon mechanism. They cover thermodynamics, energy/enthalpy inversion,
+reaction production and destruction, elemental conservation, transport, single
+precision, mixed-storage working-precision compatibility, long species names,
+zero-concentration falloff with floating-point exception traps, and compilation
+with OpenACC and OpenMP directives. The directive tests
+execute on the host; they do not validate GPU offload on accelerator hardware.
+
+Initialization tests check stream limits, normalization, elemental composition,
+enthalpy, density and a burning hot profile for both mixing-layer mechanisms.

examples/2D_reacting_mixing_layer/case.py

@@ -17,13 +17,8 @@
 parser = argparse.ArgumentParser(prog="2D_reacting_mixing_layer", formatter_class=argparse.ArgumentDefaultsHelpFormatter)
 parser.add_argument("--mfc", type=json.loads, default="{}", metavar="DICT", help="MFC's toolchain's internal state.")
 parser.add_argument("--scale", type=float, default=1.0, help="Scales cross-stream grid resolution; use <1 for cheap runs.")
-# Default is the cold (non-reacting, mollified) profile: the flamelet Newton/BDF solve
-# runs on a fixed 101-point mixture-fraction grid independent of --scale, so it is not
-# cheap to skip via a smaller physical grid -- and this file must load fast every time,
-# since it's invoked by `./mfc.sh validate`/precheck (batched over all examples) and
-# twice per `./mfc.sh run` (see toolchain/mfc/run/run.py + build.py). Pass --hot for the
-# real, physically-converged reacting profile.
-parser.add_argument("--hot", action="store_true", help="Run the full flamelet Newton/BDF solve for a physically-converged reacting profile (slow; skipped by default).")
+# Cold enthalpy-mixed streams by default; --hot maps a Cantera counterflow flame.
+parser.add_argument("--hot", action="store_true", help="Run the Cantera counterflow flame solve for a reacting initial profile (slow; skipped by default).")
 args = parser.parse_args()
 
 # Physical parameters: representative temporal H2/air mixing layer.
@@ -35,7 +30,7 @@
 mole_fraction_fu = 1.0
 vort_thickness = 1.0e-3
 mach_c = 0.3
-num_iter = 5
+flame_strain_rate = 100.0  # Nominal counterflow strain rate [1/s] for --hot initialization.
 
 # Grid: x = cross-stream (flamelet profile axis), y = streamwise (periodic, extruded).
 cross_min, cross_max = -10.0, 10.0
@@ -54,6 +49,9 @@
 # Key the cache on grid size + mode + physics so a cached IC isn't silently reused across
 # a --hot/cold switch or a physical-parameter change that leaves the line count unchanged.
 cache_key = {
+    "grid": grid,
+    "pressure": pressure,
+    "fuel": fuel,
     "cold": not args.hot,
     "lines": len(cross_coord),
     "vort_thickness": vort_thickness,
@@ -62,20 +60,14 @@
     "mach_c": mach_c,
     "mole_fraction_ox": mole_fraction_ox,
     "mole_fraction_fu": mole_fraction_fu,
-    "num_iter": num_iter,
+    "flame_strain_rate": flame_strain_rate,
+    "initializer": flamelet_ic.INITIALIZER_VERSION,
+    "mechanism": flamelet_ic.mechanism_fingerprint(sol),
 }
 if not flamelet_ic.ic_cache_valid(ic_dir, "000000", len(cross_coord), cache_key):
-    import jax.numpy as jnp
-    from pyrometheus.codegen.python import PythonCodeGenerator
-    from pyrometheus.flamelets.make_pyro import make_pyro_object
-
-    pyro_cls = PythonCodeGenerator.get_thermochem_class(sol)
-    pyro_gas = make_pyro_object(pyro_cls, jnp)
-
     flamelet_ic.generate_ic_files(
         output_dir=ic_dir,
         sol=sol,
-        pyro_gas=pyro_gas,
         cross_coord=cross_coord,
         pressure=pressure,
         temperature_ox=temperature_ox,
@@ -85,7 +77,7 @@
         mole_fraction_fu=mole_fraction_fu,
         vort_thickness=vort_thickness,
         mach_c=mach_c,
-        num_iter=num_iter,
+        strain_rate=flame_strain_rate,
         cold=not args.hot,
     )
     flamelet_ic.write_cache_key(ic_dir, cache_key)