diff --git a/README.md b/README.md index 796b4dd72..4f73fec9d 100644 --- a/README.md +++ b/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) diff --git a/docs/documentation/case.md b/docs/documentation/case.md index c37dbba0c..a537e5ed6 100644 --- a/docs/documentation/case.md +++ b/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 | diff --git a/docs/documentation/equations.md b/docs/documentation/equations.md index 867dc40b8..cebabf090 100644 --- a/docs/documentation/equations.md +++ b/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". --- diff --git a/docs/documentation/thermochemistry.md b/docs/documentation/thermochemistry.md new file mode 100644 index 000000000..7398d0593 --- /dev/null +++ b/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. diff --git a/examples/2D_reacting_mixing_layer/case.py b/examples/2D_reacting_mixing_layer/case.py index f6b322f04..5415ec33f 100644 --- a/examples/2D_reacting_mixing_layer/case.py +++ b/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) diff --git a/examples/2D_reacting_mixing_layer/flamelet_ic.py b/examples/2D_reacting_mixing_layer/flamelet_ic.py index f07b92d31..815632590 100644 --- a/examples/2D_reacting_mixing_layer/flamelet_ic.py +++ b/examples/2D_reacting_mixing_layer/flamelet_ic.py @@ -1,276 +1,31 @@ -"""Flamelet-based initial condition generator for a temporal reacting mixing layer. +"""Grid and MFC file layout for Cantera mixing-layer initial profiles. -Adapted from an external Cantera + Pyrometheus(flamelets) + JAX tool that used to write -MFC's internal binary restart format directly. Here the same 1-D flamelet solve instead -writes plain-text `prim..00..dat` files read by hcid=273 -(`src/common/include/2dHardcodedIC.fpp`), so `pre_process` — not this script — owns the -binary layout. - -AXIS CONVENTION (do not "fix" this to look more natural -- it is required by hcid=273): - MFC x -> cross-stream / flame coordinate (the flamelet profile varies here) - MFC y -> streamwise (periodic, statistically homogeneous; extruded by hcid=273) -This is the reverse of the more intuitive x=streamwise labeling, because -`ExtrusionHardcodedIC.fpp` mechanically requires the profile to vary with MFC's x and be -replicated across MFC's y. See docs/documentation/case.md's hcid=273 entry. - -The streamwise mean velocity (which must vary along MFC's x, i.e. the profile axis) is -written into the mom%beg (index 2) file slot -- the slot that would otherwise hold the -(always-zero, in this unperturbed base state) cross-stream velocity. `case(273)` in -2dHardcodedIC.fpp swaps it into its correct place (mom%end) at load time. +The common thermodynamic and flame solve lives in toolchain/mfc/flamelet.py. """ import contextlib -import json import os import sys -import time -from dataclasses import dataclass +from pathlib import Path -import jax.numpy as jnp import numpy as np -from pyrometheus.flamelets.domain import Domain, DomainConfig -from pyrometheus.flamelets.solver import FlameletSolver -from pyrometheus.flamelets.state import FlameletState -from pyrometheus.flamelets.utils import bell_profile, stoichiometric_mixture_fraction -from scipy.interpolate import interp1d - - -@dataclass -class SimulationFields: - """1-D flamelet profiles on the flame/cross-stream (MFC x) coordinate.""" - - mixture_fraction: jnp.ndarray - temperature: jnp.ndarray - pressure: jnp.ndarray - velocity: jnp.ndarray # streamwise (MFC y) velocity - mass_fractions: jnp.ndarray # shape (Ns, nx) - - -def diffusivity(pyro_gas, pressure, temperature, mass_fractions): - """Unity-Lewis diffusivity (scalar or array).""" - k = pyro_gas.get_mixture_thermal_conductivity_mixavg(temperature, mass_fractions) - rho = pyro_gas.get_density(pressure, temperature, mass_fractions) - cp = pyro_gas.get_mixture_specific_heat_cp_mass(temperature, mass_fractions) - return k / (rho * cp) - - -def new_dissipation_profile(z, z_st, val_st, bval): - """Bell-shaped scalar dissipation profile with value val_st at z_st, bval at ends.""" - profile = val_st * bell_profile(z) / bell_profile(z_st) - for ib in (0, -1): - profile = profile.at[ib].set(bval) - return profile - - -def configure_flamelet_solver(): - return { - "max_attempts": 10, - "bdf": { - "maxsteps": 10, - "time_step": 1e-5, - "newton": {"maxiter": 20, "tol": 1e-9}, - }, - "newton": {"maxiter": 10, "tol": 1e-8}, - "eos": { - "maxiter": 40, - "tol": 1e-8, - "update_size": 0.1, - "update_method": "gauss_newton", - }, - } - - -def streams(sol, fuel, pres, temp_ox, temp_fu, molefrac_ox, molefrac_fu, vort_thickness, mach_c): - """Thermodynamic state and velocities for the oxidizer and fuel streams - (temporal evolution: symmetric convective frame, u_ox = -u_fu).""" - molefrac_di = 1 - molefrac_ox - sol.TPX = temp_ox, pres, f"O2:{molefrac_ox}, N2:{molefrac_di}" - y_ox = sol.Y - h_ox = sol.enthalpy_mass - c_ox = np.sqrt((sol.cp_mass / sol.cv_mass) * sol.P / sol.density) - nu_ox = sol.viscosity / sol.density - rho_ox = sol.density - molefrac_di = 1 - molefrac_fu - sol.TPX = temp_fu, pres, f"{fuel}:{molefrac_fu}, N2:{molefrac_di}" - y_fu = sol.Y - h_fu = sol.enthalpy_mass - c_fu = np.sqrt((sol.cp_mass / sol.cv_mass) * sol.P / sol.density) - nu_fu = sol.viscosity / sol.density - rho_fu = sol.density - - z_st = stoichiometric_mixture_fraction(sol, y_ox, y_fu) - print(f"Stoichiometric mixture fraction: Z_st = {z_st:.3f}") - - u_ox = 0.5 * mach_c * (c_ox + c_fu) - u_fu = -0.5 * mach_c * (c_ox + c_fu) - delta_u = 2 * u_ox - - print(f"Convective Mach: Ma = {delta_u / (c_ox + c_fu)}") - print(f"Reynolds number: Re = {0.5 * delta_u * vort_thickness / nu_ox}") - return ( - (h_ox, y_ox, c_ox, rho_ox, u_ox, nu_ox), - (h_fu, y_fu, c_fu, rho_fu, u_fu, nu_fu), - z_st, - ) - - -def equilibrium_profiles(sol, pres, h_ox, h_fu, y_ox, y_fu, mixture_fraction): - """Equilibrium flamelet state at every mixture-fraction point.""" - num_z = len(mixture_fraction) - linear_enthalpy = h_ox + (h_fu - h_ox) * mixture_fraction - linear_mass_frac = (y_ox + (y_fu - y_ox) * mixture_fraction[:, None]).T - - state_equil = FlameletState( - enthalpy=linear_enthalpy, - mass_fractions=linear_mass_frac, - ) - temp_equil = jnp.zeros(num_z, dtype=jnp.float64) - for i in range(num_z): - sol.HPY = (state_equil.enthalpy[i], pres, state_equil.mass_fractions[:, i]) - sol.equilibrate("HP") - temp_equil = temp_equil.at[i].set(sol.T) - state_equil.mass_fractions = state_equil.mass_fractions.at[:, i].set(sol.Y) - - print(f"Equilibrium profiles: Max temperature T_max = {temp_equil.max():.3f} [K]") - return state_equil, temp_equil - - -def baseline_s_curve_state(z_st, pres, h_ox, h_fu, state_guess, temp_guess, mixture_fraction, solver_options, solver): - """Solve the flamelet at two high dissipation rates (also JIT-warms the solver).""" - - def _diss(val_st): - return new_dissipation_profile(mixture_fraction, z_st, val_st, 2 * solver.domain.jac[0] ** 2) - - def _visc(): - return new_dissipation_profile(mixture_fraction, z_st, 0, 0) - - t0 = time.time() - wmp = solver.warmup( - "flamelet_newton_step", - state_guess, - _diss(100), - _visc(), - temp_guess, - pres, - h_ox, - h_fu, - ) - wmp.enthalpy.block_until_ready() - print(f"Warm up time: {(time.time() - t0):.4e} s") - - t0 = time.time() - wmp, _, _ = solver.warmup( - "flamelet_time_step", - state_guess, - solver_options["bdf"]["newton"]["maxiter"], - solver_options["bdf"]["newton"]["tol"], - state_guess, - solver_options["bdf"]["time_step"], - _diss(100), - _visc(), - temp_guess, - pres, - h_ox, - h_fu, - ) - wmp.enthalpy.block_until_ready() - print(f"Warm up time: {(time.time() - t0):.4e} s") - - def _solve(val_st, try_newton, t_in, s_in): - return solver.solve( - True, - solver_options["newton"]["maxiter"], - solver_options["newton"]["tol"], - solver_options["bdf"]["newton"]["maxiter"], - solver_options["bdf"]["newton"]["tol"], - solver_options["bdf"]["time_step"], - solver_options["bdf"]["maxsteps"], - try_newton, - solver_options["max_attempts"], - _diss(val_st), - _visc(), - t_in, - pres, - h_ox, - h_fu, - s_in, - ) - - t0 = time.time() - state_base, temp_base = _solve(100, False, temp_guess, state_guess) - state_base.enthalpy.block_until_ready() - print(f"Run time: {(time.time() - t0):.4e} s") - - t0 = time.time() - state_base, temp_base = _solve(1000, True, temp_base, state_base) - state_base.enthalpy.block_until_ready() - print(f"Run time: {(time.time() - t0):.4e} s") - - return state_base, temp_base - - -def find_flame_dissipation_rate(pyro_gas, z_st, mixture_fraction, cross_coord, num_iter, sim_fields_cold, state_guess, temp_guess, pres, h_ox, h_fu, solver_options, solver): - """Iterate the scalar dissipation rate to match the cold mixture-fraction field.""" - cold_diff = diffusivity( - pyro_gas, - sim_fields_cold.pressure, - sim_fields_cold.temperature, - sim_fields_cold.mass_fractions, - ) - dz_dy = np.gradient(sim_fields_cold.mixture_fraction, cross_coord) - diss_rate = 2 * cold_diff * (dz_dy**2) - - sim_fields = SimulationFields( - mixture_fraction=sim_fields_cold.mixture_fraction, - velocity=sim_fields_cold.velocity, - pressure=sim_fields_cold.pressure, - temperature=sim_fields_cold.temperature, - mass_fractions=sim_fields_cold.mass_fractions, - ) - state_it = state_guess - temp_it = temp_guess - try_newton = False - - for j_chi in range(num_iter): - _diss = interp1d(sim_fields.mixture_fraction, diss_rate, fill_value="extrapolate")(mixture_fraction) - _diss[0] = 2 * solver.domain.jac[0] ** 2 - _diss[-1] = 2 * solver.domain.jac[0] ** 2 - - state_it, temp_it = solver.solve( - True, - solver_options["newton"]["maxiter"], - solver_options["newton"]["tol"], - solver_options["bdf"]["newton"]["maxiter"], - solver_options["bdf"]["newton"]["tol"], - solver_options["bdf"]["time_step"], - solver_options["bdf"]["maxsteps"], - try_newton, - solver_options["max_attempts"], - _diss, - jnp.zeros_like(_diss), - temp_it, - pres, - h_ox, - h_fu, - state_it, - ) - if not try_newton: - try_newton = True - - sim_fields.temperature = interp1d(mixture_fraction, temp_it)(sim_fields.mixture_fraction) - sim_fields.mass_fractions = interp1d(mixture_fraction, state_it.mass_fractions)(sim_fields.mixture_fraction) - - diff = diffusivity(pyro_gas, pres, sim_fields.temperature, sim_fields.mass_fractions) - diss_rate = 2 * diff * (dz_dy**2) - chi_st = interp1d(sim_fields.mixture_fraction, diss_rate)(z_st) - print(f"---> j = {j_chi}, chi_max = {diss_rate.max():.4f}, " f"chi_st = {chi_st:.4f}, T_max = {temp_it.max():.4f}") - - return sim_fields +# Also support direct `python case.py` from a source checkout. +sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "toolchain")) +from mfc.flamelet import ( + INITIALIZER_VERSION, + create_simulation_fields, + density, + ic_cache_valid, + mechanism_fingerprint, + reference_fluid_properties, + streams, + write_cache_key, +) def compute_grid(vort_thickness, cross_min, cross_max, points_per_cross, stream_min, stream_max, num_y): - """Pure grid arithmetic -- no Cantera/JAX. Always cheap to call, including on an + """Pure grid arithmetic -- no Cantera. Always cheap to call, including on an IC/ cache hit, so case.py never needs to run the flamelet solve just to learn its own domain size. @@ -301,111 +56,6 @@ def compute_grid(vort_thickness, cross_min, cross_max, points_per_cross, stream_ return cross_coord, grid -def reference_fluid_properties(sol, temperature_ox, pressure, mole_fraction_ox): - """Cheap (no equilibration/solve) reference gamma and viscosity for fluid_pp(1).""" - sol.TPX = temperature_ox, pressure, f"O2:{mole_fraction_ox}, N2:{1 - mole_fraction_ox}" - return {"gamma": float(sol.cp_mass / sol.cv_mass), "viscosity": float(sol.viscosity)} - - -def ic_cache_valid(ic_dir, file_extension, expected_lines, cache_key=None): - """True only if IC/ has prim.1..dat with exactly the current grid's expected - line count AND (if `cache_key` is given) a matching .cache_key.json. A bare "IC/ is - non-empty" check isn't enough: this same case.py is invoked with different --scale - values by different toolchain paths (e.g. `./mfc.sh validate` during precheck uses no - args/default scale, while a registered test passes its own --scale) that can share - this directory, so a cache populated by one grid size must not be silently reused by a - run expecting a different one -- reading it would desync the Fortran reader (hcid=273 - expects exactly len(cross_coord) lines) from the actual grid. The `cache_key` further - guards against silently reusing an IC generated with a different mode (--hot vs cold) - or different physical parameters that leave the line count unchanged.""" - path = os.path.join(ic_dir, f"prim.1.00.{file_extension}.dat") - if not os.path.isfile(path): - return False - with open(path) as fh: - if sum(1 for _ in fh) != expected_lines: - return False - if cache_key is not None: - key_path = os.path.join(ic_dir, ".cache_key.json") - if not os.path.isfile(key_path): - return False - try: - with open(key_path) as fh: - stored = json.load(fh) - except (OSError, ValueError): - return False - if stored != cache_key: - return False - return True - - -def write_cache_key(ic_dir, cache_key): - """Record the parameters an IC/ was generated with, so ic_cache_valid can detect a - stale cache (different --hot/cold mode or physical parameters at the same grid size).""" - with open(os.path.join(ic_dir, ".cache_key.json"), "w") as fh: - json.dump(cache_key, fh, sort_keys=True) - - -def create_simulation_fields(pyro_gas, sol, pres, temp_ox, temp_fu, cross_coord, vort_thickness, stream_ox, stream_fu, z_st, num_iter, cold): - """1-D flamelet profiles on the cross-stream (MFC x) coordinate `cross_coord`.""" - mollifier = 0.5 * (1 - np.tanh(2 * cross_coord / vort_thickness)) - - h_ox, y_ox, _, _, u_ox, nu_ox = stream_ox - h_fu, y_fu, _, _, u_fu, _ = stream_fu - - sim_fields_cold = SimulationFields( - mixture_fraction=mollifier, - temperature=pyro_gas.get_temperature_from_enthalpy( - h_ox + (h_fu - h_ox) * mollifier, - (y_ox + (y_fu - y_ox) * mollifier[:, None]).T, - temp_ox + (temp_fu - temp_ox) * mollifier, - ), - pressure=pres * jnp.ones_like(mollifier), - velocity=(u_ox + (u_fu - u_ox) * mollifier), - mass_fractions=(y_ox + (y_fu - y_ox) * mollifier[:, None]).T, - ) - - if cold: - return sim_fields_cold - - domain_config = DomainConfig(num_x=101, x_l=0, x_r=1) - domain = Domain(domain_config) - solver = FlameletSolver(domain, pyro_gas, (y_ox, y_fu)) - mixture_fraction = jnp.array(domain.x) - - state_equil, temp_equil = equilibrium_profiles(sol, pres, h_ox, h_fu, y_ox, y_fu, mixture_fraction) - solver_options = configure_flamelet_solver() - # baseline_s_curve_state's return value is unused -- its purpose is JIT warmup - # for the JAX solver calls that follow (matches the upstream tool's structure). - baseline_s_curve_state( - z_st, - pres, - h_ox, - h_fu, - state_equil, - temp_equil, - mixture_fraction, - solver_options, - solver, - ) - - sim_fields = find_flame_dissipation_rate( - pyro_gas, - z_st, - mixture_fraction, - cross_coord, - num_iter, - sim_fields_cold, - state_equil, - temp_equil, - pres, - h_ox, - h_fu, - solver_options, - solver, - ) - return sim_fields - - def write_hcid_ic(output_dir, cross_coord, density, streamwise_velocity, pressure, mass_fractions, file_extension="000000"): """Write hcid=273 IC text files: prim..00..dat, one `x value` pair per line. @@ -433,17 +83,17 @@ def write_hcid_ic(output_dir, cross_coord, density, streamwise_velocity, pressur def generate_ic_files( - *, output_dir, sol, pyro_gas, cross_coord, pressure, temperature_ox, temperature_fu, fuel, mole_fraction_ox, mole_fraction_fu, vort_thickness, mach_c, num_iter, cold, file_extension="000000" + *, output_dir, sol, cross_coord, pressure, temperature_ox, temperature_fu, fuel, mole_fraction_ox, mole_fraction_fu, vort_thickness, mach_c, strain_rate, cold, file_extension="000000" ): """Run the flamelet solve (expensive when cold=False) and write hcid=273 IC files on `cross_coord` (from `compute_grid`, so the file spacing exactly matches the grid case.py declares). - All Cantera/JAX/Pyrometheus stdout diagnostics are redirected to stderr: case.py's + All Cantera stdout diagnostics are redirected to stderr: case.py's contract requires its entire stdout to be exactly one JSON line. """ with contextlib.redirect_stdout(sys.stderr): - stream_ox, stream_fu, z_st = streams( + stream_ox, stream_fu, _ = streams( sol, fuel, pressure, @@ -456,7 +106,6 @@ def generate_ic_files( ) sim_fields = create_simulation_fields( - pyro_gas, sol, pressure, temperature_ox, @@ -465,8 +114,7 @@ def generate_ic_files( vort_thickness, stream_ox, stream_fu, - z_st, - num_iter, + strain_rate, cold, ) @@ -474,7 +122,7 @@ def generate_ic_files( pressure_1d = np.array(sim_fields.pressure) velocity_1d = np.array(sim_fields.velocity) mass_fractions_1d = np.array(sim_fields.mass_fractions) - density_1d = np.array(pyro_gas.get_density(pressure_1d, temperature_1d, mass_fractions_1d)) + density_1d = np.array(density(sol, pressure_1d, temperature_1d, mass_fractions_1d)) # Fail at generation time rather than writing a non-finite IC that would only # surface downstream as a cryptic VCFL=Inf crash (e.g. a diverged --hot solve). diff --git a/examples/2D_spatial_reacting_mixing_layer/case.py b/examples/2D_spatial_reacting_mixing_layer/case.py index 1743d4609..7c934f876 100644 --- a/examples/2D_spatial_reacting_mixing_layer/case.py +++ b/examples/2D_spatial_reacting_mixing_layer/case.py @@ -34,7 +34,7 @@ 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). See 2D_reacting_mixing_layer/case.py for why the default is cold.", + help="Run the Cantera counterflow flame solve for a reacting initial profile (slow; skipped by default). See 2D_reacting_mixing_layer/case.py for why the default is cold.", ) args = parser.parse_args() @@ -55,7 +55,7 @@ # delta_u = u_ox - u_fu (oxidizer coflow faster than the fuel jet). mach_ox = 0.9 mach_fu = 0.1 -num_iter = 5 +flame_strain_rate = 100.0 # Nominal counterflow strain rate [1/s] for --hot initialization. # Grid: x = streamwise (inflow at 0), y = cross-stream. stream_min, stream_max = 0.0, 15.0 @@ -91,6 +91,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(stream_coord) * len(cross_coord), "vort_thickness": vort_thickness, @@ -100,20 +103,14 @@ "mach_fu": mach_fu, "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(stream_coord) * 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_spatial( output_dir=ic_dir, sol=sol, - pyro_gas=pyro_gas, stream_coord=stream_coord, cross_coord=cross_coord, pressure=pressure, @@ -125,7 +122,7 @@ vort_thickness=vort_thickness, mach_ox=mach_ox, mach_fu=mach_fu, - num_iter=num_iter, + strain_rate=flame_strain_rate, cold=not args.hot, ) flamelet_ic.write_cache_key(ic_dir, cache_key) diff --git a/examples/2D_spatial_reacting_mixing_layer/flamelet_ic.py b/examples/2D_spatial_reacting_mixing_layer/flamelet_ic.py index fec1d55d0..d5c793571 100644 --- a/examples/2D_spatial_reacting_mixing_layer/flamelet_ic.py +++ b/examples/2D_spatial_reacting_mixing_layer/flamelet_ic.py @@ -1,384 +1,31 @@ -"""Flamelet-based initial/inflow condition generator for a spatially-evolving reacting -mixing layer. +"""Grid and MFC file layout for Cantera mixing-layer initial profiles. -The 1-D flamelet solver core here is identical to (and kept in sync by hand with) -examples/2D_reacting_mixing_layer/flamelet_ic.py -- duplicated rather than imported -across example directories, since each MFC example is a self-contained, copyable unit -(no example currently imports from a sibling directory). Only the writer and grid -functions differ, because this case needs a genuinely full 2D field, not an extrusion. - -AXIS CONVENTION (the reverse of the temporal example's x/y roles, and deliberately so): - MFC x -> streamwise (matches spatial_bf's hardcoded advecting-direction convention - in src/simulation/m_body_forces.fpp: theta_x includes a conv_vel*t term, - theta_y does not) - MFC y -> cross-stream (the flamelet profile varies here) -Uses hcid=274 (src/common/include/2dHardcodedIC.fpp): a full (x,y) field read with no -extrusion assumption, unlike hcid=273 (which mechanically requires cross-stream=x and -would conflict with spatial_bf's fixed axis convention -- see case.py's module -docstring). +The common thermodynamic and flame solve lives in toolchain/mfc/flamelet.py. """ import contextlib -import json import os import sys -import time -from dataclasses import dataclass +from pathlib import Path -import jax.numpy as jnp import numpy as np -from pyrometheus.flamelets.domain import Domain, DomainConfig -from pyrometheus.flamelets.solver import FlameletSolver -from pyrometheus.flamelets.state import FlameletState -from pyrometheus.flamelets.utils import bell_profile, stoichiometric_mixture_fraction -from scipy.interpolate import interp1d - - -@dataclass -class SimulationFields: - """1-D flamelet profiles on the flame/cross-stream (MFC y) coordinate.""" - - mixture_fraction: jnp.ndarray - temperature: jnp.ndarray - pressure: jnp.ndarray - velocity: jnp.ndarray # streamwise (MFC x) velocity - mass_fractions: jnp.ndarray # shape (Ns, ny) - - -def diffusivity(pyro_gas, pressure, temperature, mass_fractions): - """Unity-Lewis diffusivity (scalar or array).""" - k = pyro_gas.get_mixture_thermal_conductivity_mixavg(temperature, mass_fractions) - rho = pyro_gas.get_density(pressure, temperature, mass_fractions) - cp = pyro_gas.get_mixture_specific_heat_cp_mass(temperature, mass_fractions) - return k / (rho * cp) - - -def new_dissipation_profile(z, z_st, val_st, bval): - """Bell-shaped scalar dissipation profile with value val_st at z_st, bval at ends.""" - profile = val_st * bell_profile(z) / bell_profile(z_st) - for ib in (0, -1): - profile = profile.at[ib].set(bval) - return profile - - -def configure_flamelet_solver(): - return { - "max_attempts": 10, - "bdf": { - "maxsteps": 10, - "time_step": 1e-5, - "newton": {"maxiter": 20, "tol": 1e-9}, - }, - "newton": {"maxiter": 10, "tol": 1e-8}, - "eos": { - "maxiter": 40, - "tol": 1e-8, - "update_size": 0.1, - "update_method": "gauss_newton", - }, - } - - -def streams(sol, fuel, pres, temp_ox, temp_fu, molefrac_ox, molefrac_fu, vort_thickness, mach_ox, mach_fu): - """Thermodynamic state and velocities for the oxidizer and fuel streams, spatial - (lab-frame) evolution: u_ox = mach_ox*c_ox, u_fu = mach_fu*c_fu -- both streams - co-flow at their own absolute speed, giving a nonzero mean/convective velocity - (unlike the temporal case's symmetric convective frame, where u_ox = -u_fu by - construction and the mean is always zero -- unusable for spatial_bf%conv_vel, - which needs to be nonzero to drive the forcing).""" - molefrac_di = 1 - molefrac_ox - sol.TPX = temp_ox, pres, f"O2:{molefrac_ox}, N2:{molefrac_di}" - y_ox = sol.Y - h_ox = sol.enthalpy_mass - c_ox = np.sqrt((sol.cp_mass / sol.cv_mass) * sol.P / sol.density) - nu_ox = sol.viscosity / sol.density - rho_ox = sol.density - - molefrac_di = 1 - molefrac_fu - sol.TPX = temp_fu, pres, f"{fuel}:{molefrac_fu}, N2:{molefrac_di}" - y_fu = sol.Y - h_fu = sol.enthalpy_mass - c_fu = np.sqrt((sol.cp_mass / sol.cv_mass) * sol.P / sol.density) - nu_fu = sol.viscosity / sol.density - rho_fu = sol.density - - z_st = stoichiometric_mixture_fraction(sol, y_ox, y_fu) - print(f"Stoichiometric mixture fraction: Z_st = {z_st:.3f}") - - u_ox = mach_ox * c_ox - u_fu = mach_fu * c_fu - delta_u = u_ox - u_fu - - print(f"Convective Mach: Ma = {delta_u / (c_ox + c_fu)}") - print(f"Reynolds number: Re = {0.5 * delta_u * vort_thickness / nu_ox}") - return ( - (h_ox, y_ox, c_ox, rho_ox, u_ox, nu_ox), - (h_fu, y_fu, c_fu, rho_fu, u_fu, nu_fu), - z_st, - ) - - -def equilibrium_profiles(sol, pres, h_ox, h_fu, y_ox, y_fu, mixture_fraction): - """Equilibrium flamelet state at every mixture-fraction point.""" - num_z = len(mixture_fraction) - linear_enthalpy = h_ox + (h_fu - h_ox) * mixture_fraction - linear_mass_frac = (y_ox + (y_fu - y_ox) * mixture_fraction[:, None]).T - - state_equil = FlameletState( - enthalpy=linear_enthalpy, - mass_fractions=linear_mass_frac, - ) - temp_equil = jnp.zeros(num_z, dtype=jnp.float64) - for i in range(num_z): - sol.HPY = (state_equil.enthalpy[i], pres, state_equil.mass_fractions[:, i]) - sol.equilibrate("HP") - temp_equil = temp_equil.at[i].set(sol.T) - state_equil.mass_fractions = state_equil.mass_fractions.at[:, i].set(sol.Y) - - print(f"Equilibrium profiles: Max temperature T_max = {temp_equil.max():.3f} [K]") - return state_equil, temp_equil - - -def baseline_s_curve_state(z_st, pres, h_ox, h_fu, state_guess, temp_guess, mixture_fraction, solver_options, solver): - """Solve the flamelet at two high dissipation rates (also JIT-warms the solver).""" - - def _diss(val_st): - return new_dissipation_profile(mixture_fraction, z_st, val_st, 2 * solver.domain.jac[0] ** 2) - - def _visc(): - return new_dissipation_profile(mixture_fraction, z_st, 0, 0) - - t0 = time.time() - wmp = solver.warmup( - "flamelet_newton_step", - state_guess, - _diss(100), - _visc(), - temp_guess, - pres, - h_ox, - h_fu, - ) - wmp.enthalpy.block_until_ready() - print(f"Warm up time: {(time.time() - t0):.4e} s") - - t0 = time.time() - wmp, _, _ = solver.warmup( - "flamelet_time_step", - state_guess, - solver_options["bdf"]["newton"]["maxiter"], - solver_options["bdf"]["newton"]["tol"], - state_guess, - solver_options["bdf"]["time_step"], - _diss(100), - _visc(), - temp_guess, - pres, - h_ox, - h_fu, - ) - wmp.enthalpy.block_until_ready() - print(f"Warm up time: {(time.time() - t0):.4e} s") - - def _solve(val_st, try_newton, t_in, s_in): - return solver.solve( - True, - solver_options["newton"]["maxiter"], - solver_options["newton"]["tol"], - solver_options["bdf"]["newton"]["maxiter"], - solver_options["bdf"]["newton"]["tol"], - solver_options["bdf"]["time_step"], - solver_options["bdf"]["maxsteps"], - try_newton, - solver_options["max_attempts"], - _diss(val_st), - _visc(), - t_in, - pres, - h_ox, - h_fu, - s_in, - ) - - t0 = time.time() - state_base, temp_base = _solve(100, False, temp_guess, state_guess) - state_base.enthalpy.block_until_ready() - print(f"Run time: {(time.time() - t0):.4e} s") - - t0 = time.time() - state_base, temp_base = _solve(1000, True, temp_base, state_base) - state_base.enthalpy.block_until_ready() - print(f"Run time: {(time.time() - t0):.4e} s") - - return state_base, temp_base - - -def find_flame_dissipation_rate(pyro_gas, z_st, mixture_fraction, cross_coord, num_iter, sim_fields_cold, state_guess, temp_guess, pres, h_ox, h_fu, solver_options, solver): - """Iterate the scalar dissipation rate to match the cold mixture-fraction field.""" - cold_diff = diffusivity( - pyro_gas, - sim_fields_cold.pressure, - sim_fields_cold.temperature, - sim_fields_cold.mass_fractions, - ) - dz_dy = np.gradient(sim_fields_cold.mixture_fraction, cross_coord) - diss_rate = 2 * cold_diff * (dz_dy**2) - - sim_fields = SimulationFields( - mixture_fraction=sim_fields_cold.mixture_fraction, - velocity=sim_fields_cold.velocity, - pressure=sim_fields_cold.pressure, - temperature=sim_fields_cold.temperature, - mass_fractions=sim_fields_cold.mass_fractions, - ) - state_it = state_guess - temp_it = temp_guess - try_newton = False - - for j_chi in range(num_iter): - _diss = interp1d(sim_fields.mixture_fraction, diss_rate, fill_value="extrapolate")(mixture_fraction) - _diss[0] = 2 * solver.domain.jac[0] ** 2 - _diss[-1] = 2 * solver.domain.jac[0] ** 2 - - state_it, temp_it = solver.solve( - True, - solver_options["newton"]["maxiter"], - solver_options["newton"]["tol"], - solver_options["bdf"]["newton"]["maxiter"], - solver_options["bdf"]["newton"]["tol"], - solver_options["bdf"]["time_step"], - solver_options["bdf"]["maxsteps"], - try_newton, - solver_options["max_attempts"], - _diss, - jnp.zeros_like(_diss), - temp_it, - pres, - h_ox, - h_fu, - state_it, - ) - if not try_newton: - try_newton = True - - sim_fields.temperature = interp1d(mixture_fraction, temp_it)(sim_fields.mixture_fraction) - sim_fields.mass_fractions = interp1d(mixture_fraction, state_it.mass_fractions)(sim_fields.mixture_fraction) - - diff = diffusivity(pyro_gas, pres, sim_fields.temperature, sim_fields.mass_fractions) - diss_rate = 2 * diff * (dz_dy**2) - chi_st = interp1d(sim_fields.mixture_fraction, diss_rate)(z_st) - print(f"---> j = {j_chi}, chi_max = {diss_rate.max():.4f}, " f"chi_st = {chi_st:.4f}, T_max = {temp_it.max():.4f}") - - return sim_fields - - -def reference_fluid_properties(sol, temperature_ox, pressure, mole_fraction_ox): - """Cheap (no equilibration/solve) reference gamma and viscosity for fluid_pp(1).""" - sol.TPX = temperature_ox, pressure, f"O2:{mole_fraction_ox}, N2:{1 - mole_fraction_ox}" - return {"gamma": float(sol.cp_mass / sol.cv_mass), "viscosity": float(sol.viscosity)} - - -def ic_cache_valid(ic_dir, file_extension, expected_lines, cache_key=None): - """True only if IC/ has prim.1..dat with exactly the current grid's expected - line count AND (if `cache_key` is given) a matching .cache_key.json. A bare "IC/ is - non-empty" check isn't enough: this same case.py is invoked with different --scale - values by different toolchain paths (e.g. `./mfc.sh validate` during precheck uses no - args/default scale, while a registered test passes its own --scale) that can share - this directory, so a cache populated by one grid size must not be silently reused by a - run expecting a different one -- reading it would desync the Fortran reader (hcid=274 - expects exactly (m_glb+1)*(n_glb+1) lines) from the actual grid. The `cache_key` - further guards against silently reusing an IC generated with a different mode (--hot vs - cold) or different physical parameters that leave the line count unchanged.""" - path = os.path.join(ic_dir, f"prim.1.00.{file_extension}.dat") - if not os.path.isfile(path): - return False - with open(path) as fh: - if sum(1 for _ in fh) != expected_lines: - return False - if cache_key is not None: - key_path = os.path.join(ic_dir, ".cache_key.json") - if not os.path.isfile(key_path): - return False - try: - with open(key_path) as fh: - stored = json.load(fh) - except (OSError, ValueError): - return False - if stored != cache_key: - return False - return True - - -def write_cache_key(ic_dir, cache_key): - """Record the parameters an IC/ was generated with, so ic_cache_valid can detect a - stale cache (different --hot/cold mode or physical parameters at the same grid size).""" - with open(os.path.join(ic_dir, ".cache_key.json"), "w") as fh: - json.dump(cache_key, fh, sort_keys=True) - - -def create_simulation_fields(pyro_gas, sol, pres, temp_ox, temp_fu, cross_coord, vort_thickness, stream_ox, stream_fu, z_st, num_iter, cold): - """1-D flamelet profiles on the cross-stream (MFC y) coordinate `cross_coord`.""" - mollifier = 0.5 * (1 - np.tanh(2 * cross_coord / vort_thickness)) - - h_ox, y_ox, _, _, u_ox, nu_ox = stream_ox - h_fu, y_fu, _, _, u_fu, _ = stream_fu - - sim_fields_cold = SimulationFields( - mixture_fraction=mollifier, - temperature=pyro_gas.get_temperature_from_enthalpy( - h_ox + (h_fu - h_ox) * mollifier, - (y_ox + (y_fu - y_ox) * mollifier[:, None]).T, - temp_ox + (temp_fu - temp_ox) * mollifier, - ), - pressure=pres * jnp.ones_like(mollifier), - velocity=(u_ox + (u_fu - u_ox) * mollifier), - mass_fractions=(y_ox + (y_fu - y_ox) * mollifier[:, None]).T, - ) - - if cold: - return sim_fields_cold - - domain_config = DomainConfig(num_x=101, x_l=0, x_r=1) - domain = Domain(domain_config) - solver = FlameletSolver(domain, pyro_gas, (y_ox, y_fu)) - mixture_fraction = jnp.array(domain.x) - - state_equil, temp_equil = equilibrium_profiles(sol, pres, h_ox, h_fu, y_ox, y_fu, mixture_fraction) - solver_options = configure_flamelet_solver() - # baseline_s_curve_state's return value is unused -- its purpose is JIT warmup - # for the JAX solver calls that follow (matches the upstream tool's structure). - baseline_s_curve_state( - z_st, - pres, - h_ox, - h_fu, - state_equil, - temp_equil, - mixture_fraction, - solver_options, - solver, - ) - sim_fields = find_flame_dissipation_rate( - pyro_gas, - z_st, - mixture_fraction, - cross_coord, - num_iter, - sim_fields_cold, - state_equil, - temp_equil, - pres, - h_ox, - h_fu, - solver_options, - solver, - ) - return sim_fields +# Also support direct `python case.py` from a source checkout. +sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "toolchain")) +from mfc.flamelet import ( + INITIALIZER_VERSION, + create_simulation_fields, + density, + ic_cache_valid, + mechanism_fingerprint, + reference_fluid_properties, + streams, + write_cache_key, +) def compute_grid_spatial(vort_thickness, cross_min, cross_max, points_per_cross, stream_min, stream_max, points_per_stream): - """Pure grid arithmetic -- no Cantera/JAX. Always cheap to call, including on an + """Pure grid arithmetic -- no Cantera. Always cheap to call, including on an IC/ cache hit, so case.py never needs to run the flamelet solve just to learn its own domain size. @@ -448,7 +95,6 @@ def generate_ic_files_spatial( *, output_dir, sol, - pyro_gas, stream_coord, cross_coord, pressure, @@ -460,7 +106,7 @@ def generate_ic_files_spatial( vort_thickness, mach_ox, mach_fu, - num_iter, + strain_rate, cold, file_extension="000000", ): @@ -469,11 +115,11 @@ def generate_ic_files_spatial( axis -- it only develops streamwise variation once simulation starts evolving it via the inflow BC (bc_x%beg=-17) and spatial_bf forcing. - All Cantera/JAX/Pyrometheus stdout diagnostics are redirected to stderr: case.py's + All Cantera stdout diagnostics are redirected to stderr: case.py's contract requires its entire stdout to be exactly one JSON line. """ with contextlib.redirect_stdout(sys.stderr): - stream_ox, stream_fu, z_st = streams( + stream_ox, stream_fu, _ = streams( sol, fuel, pressure, @@ -487,7 +133,6 @@ def generate_ic_files_spatial( ) sim_fields = create_simulation_fields( - pyro_gas, sol, pressure, temperature_ox, @@ -496,8 +141,7 @@ def generate_ic_files_spatial( vort_thickness, stream_ox, stream_fu, - z_st, - num_iter, + strain_rate, cold, ) @@ -505,7 +149,7 @@ def generate_ic_files_spatial( pressure_1d = np.array(sim_fields.pressure) velocity_1d = np.array(sim_fields.velocity) mass_fractions_1d = np.array(sim_fields.mass_fractions) - density_1d = np.array(pyro_gas.get_density(pressure_1d, temperature_1d, mass_fractions_1d)) + density_1d = np.array(density(sol, pressure_1d, temperature_1d, mass_fractions_1d)) # Fail at generation time rather than writing a non-finite IC that would only # surface downstream as a cryptic VCFL=Inf crash (e.g. a diverged --hot solve). diff --git a/examples/3D_reacting_mixing_layer/README.md b/examples/3D_reacting_mixing_layer/README.md index 17f5c2c18..1abac1051 100644 --- a/examples/3D_reacting_mixing_layer/README.md +++ b/examples/3D_reacting_mixing_layer/README.md @@ -1,8 +1,8 @@ # 3D Temporal Reacting Mixing Layer (H2/N2 - air, Mc = 1.5) A temporally-evolving supersonic reacting shear layer between a hot air stream and an -N2-diluted hydrogen stream. The base state comes from a 1-D flamelet solve (Cantera + -Pyrometheus + JAX) extruded into 3D by `hcid=371`. This is the supersonic counterpart to +N2-diluted hydrogen stream. The base state comes from Cantera-only stream mixing +(or a counterflow flame with `--hot`) extruded into 3D by `hcid=371`. This is the supersonic counterpart to `examples/2D_reacting_mixing_layer`, which runs the same flamelet machinery at `Mc = 0.3`. ## Configuration @@ -52,7 +52,11 @@ The file spacing must match the run grid. A mismatch aborts in `pre_process`, so `--scale` shrinks the grid for cheap runs; `--scale 0.05` gives 32^3, which is what the `3D -> Chemistry -> Reacting Mixing Layer` regression test uses. `--hot` runs the full -flamelet Newton/BDF solve instead of the default cold mollified profile. +Cantera counterflow flame solve instead of the default cold mollified profile. +`flame_strain_rate` in `case.py` sets the nominal inlet strain rate (default 100/s). +The hot profile is mapped by mixture fraction onto the prescribed shear layer; it +replaces the former scalar-dissipation-matched flamelet initialization. See +@ref thermochemistry "Thermochemistry implementation" for the model and validation details. The mechanism ships alongside the case as `sandiego.yaml` (UC San Diego Combustion Research Group, ). diff --git a/examples/3D_reacting_mixing_layer/case.py b/examples/3D_reacting_mixing_layer/case.py index 82798de6e..555da450c 100644 --- a/examples/3D_reacting_mixing_layer/case.py +++ b/examples/3D_reacting_mixing_layer/case.py @@ -20,8 +20,8 @@ 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.") # See examples/2D_reacting_mixing_layer/case.py for why the default is the cold (mollified, -# non-reacting) profile and --hot runs the full flamelet Newton/BDF solve. -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).") +# non-reacting) profile and --hot runs the Cantera counterflow flame solve. +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() # The 2D temporal case's H2/air mixing layer at mach_c=1.5 instead of 0.3, to probe @@ -38,7 +38,7 @@ mole_fraction_fu = 0.5 vort_thickness = 1.0e-3 mach_c = 1.5 -num_iter = 5 +flame_strain_rate = 100.0 # Nominal counterflow strain rate [1/s] for --hot initialization. # Grid: x = streamwise (periodic), y = cross-stream (flamelet profile axis), z = spanwise # (periodic). Wang et al. (C&F 2024)'s temporal mixing-layer DNS at their resolution @@ -70,6 +70,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(x_coord) * len(cross_coord), "vort_thickness": vort_thickness, @@ -78,21 +81,15 @@ "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), "perturb_seed": perturb_seed, } if not flamelet_ic.ic_cache_valid(ic_dir, "000000", len(x_coord) * 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, x_coord=x_coord, pressure=pressure, @@ -103,7 +100,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, perturb_seed=perturb_seed, ) @@ -151,8 +148,7 @@ "chemistry": "T", "chem_params%diffusion": "T", "chem_params%reactions": "T", - # Unity-Lewis, matching the flamelet solve's own assumption (flamelet_ic.py's - # diffusivity() uses D_k = k/(rho*cp) for every species). + # Unity-Lewis, matching the Cantera counterflow initialization. "chem_params%transport_model": 2, "files_dir": ic_dir, "file_extension": "000000", diff --git a/examples/3D_reacting_mixing_layer/flamelet_ic.py b/examples/3D_reacting_mixing_layer/flamelet_ic.py index 461621bb7..add15d42c 100644 --- a/examples/3D_reacting_mixing_layer/flamelet_ic.py +++ b/examples/3D_reacting_mixing_layer/flamelet_ic.py @@ -1,291 +1,31 @@ -"""Flamelet-based initial condition generator for a 3D temporal reacting mixing layer. - -Shares its flamelet solve (Cantera + Pyrometheus(flamelets) + JAX) with -`examples/2D_reacting_mixing_layer/flamelet_ic.py`; only the grid and IC-writing -functions differ, to extrude the 1-D profile across two directions instead of one and to -match hcid=370's file format instead of hcid=273's. - -AXIS CONVENTION: - MFC x -> streamwise (periodic, statistically homogeneous) - MFC y -> cross-stream / flame coordinate (the flamelet profile varies here) - MFC z -> spanwise (periodic, statistically homogeneous) -This is the "intuitive" x=streamwise labeling -- unlike the 2D case, hcid=370 (3D -extrusion of a 2D profile from external data, `src/common/include/3dHardcodedIC.fpp`) -reads genuine (x, y) coordinates from file and only mechanically requires the *z* axis to -be the uniform/extruded one, so x and y are free to take their natural roles here. - -In-plane (x,y) velocity perturbations are baked directly into this file (`perturb_xy`, -ported from ../tmp-mixlyr/perturb.py's perturb_temporally_evolving_layer -- computes -wavenumbers relative to domain_length/vort_thickness, so it's dimensionally correct -regardless of the case's absolute length-scale choice, unlike MFC's own -`mixlayer_perturb` which hardcodes an absolute wavenumber range unsuited to this -millimeter-scale domain). The spanwise (z) component isn't -- and can't be, since hcid=370 -extrudes this (x,y) file uniformly across z -- so hcid=371 (`src/common/include/ -3dHardcodedIC.fpp`) additionally modulates the already-written (x,y) perturbation by a -closed-form z-dependent factor at IC-assignment time, using only the z domain's own extent -(no new case parameters needed). +"""Grid and MFC file layout for Cantera mixing-layer initial profiles. + +The common thermodynamic and flame solve lives in toolchain/mfc/flamelet.py. """ import contextlib -import json import os import sys -import time -from dataclasses import dataclass +from pathlib import Path -import jax.numpy as jnp import numpy as np -from pyrometheus.flamelets.domain import Domain, DomainConfig -from pyrometheus.flamelets.solver import FlameletSolver -from pyrometheus.flamelets.state import FlameletState -from pyrometheus.flamelets.utils import bell_profile, stoichiometric_mixture_fraction -from scipy.interpolate import interp1d - - -@dataclass -class SimulationFields: - """1-D flamelet profiles on the flame/cross-stream (MFC y) coordinate.""" - - mixture_fraction: jnp.ndarray - temperature: jnp.ndarray - pressure: jnp.ndarray - velocity: jnp.ndarray # streamwise (MFC x) velocity - mass_fractions: jnp.ndarray # shape (Ns, ny) - - -def diffusivity(pyro_gas, pressure, temperature, mass_fractions): - """Unity-Lewis diffusivity (scalar or array).""" - k = pyro_gas.get_mixture_thermal_conductivity_mixavg(temperature, mass_fractions) - rho = pyro_gas.get_density(pressure, temperature, mass_fractions) - cp = pyro_gas.get_mixture_specific_heat_cp_mass(temperature, mass_fractions) - return k / (rho * cp) - - -def new_dissipation_profile(z, z_st, val_st, bval): - """Bell-shaped scalar dissipation profile with value val_st at z_st, bval at ends.""" - profile = val_st * bell_profile(z) / bell_profile(z_st) - for ib in (0, -1): - profile = profile.at[ib].set(bval) - return profile - - -def configure_flamelet_solver(): - return { - "max_attempts": 10, - "bdf": { - "maxsteps": 10, - "time_step": 1e-5, - "newton": {"maxiter": 20, "tol": 1e-9}, - }, - "newton": {"maxiter": 10, "tol": 1e-8}, - "eos": { - "maxiter": 40, - "tol": 1e-8, - "update_size": 0.1, - "update_method": "gauss_newton", - }, - } - -def streams(sol, fuel, pres, temp_ox, temp_fu, molefrac_ox, molefrac_fu, vort_thickness, mach_c): - """Thermodynamic state and velocities for the oxidizer and fuel streams - (temporal evolution: symmetric convective frame, u_ox = -u_fu).""" - molefrac_di = 1 - molefrac_ox - sol.TPX = temp_ox, pres, f"O2:{molefrac_ox}, N2:{molefrac_di}" - y_ox = sol.Y - h_ox = sol.enthalpy_mass - c_ox = np.sqrt((sol.cp_mass / sol.cv_mass) * sol.P / sol.density) - nu_ox = sol.viscosity / sol.density - rho_ox = sol.density - - molefrac_di = 1 - molefrac_fu - sol.TPX = temp_fu, pres, f"{fuel}:{molefrac_fu}, N2:{molefrac_di}" - y_fu = sol.Y - h_fu = sol.enthalpy_mass - c_fu = np.sqrt((sol.cp_mass / sol.cv_mass) * sol.P / sol.density) - nu_fu = sol.viscosity / sol.density - rho_fu = sol.density - - z_st = stoichiometric_mixture_fraction(sol, y_ox, y_fu) - print(f"Stoichiometric mixture fraction: Z_st = {z_st:.3f}") - - u_ox = 0.5 * mach_c * (c_ox + c_fu) - u_fu = -0.5 * mach_c * (c_ox + c_fu) - delta_u = 2 * u_ox - - print(f"Convective Mach: Ma = {delta_u / (c_ox + c_fu)}") - print(f"Reynolds number: Re = {0.5 * delta_u * vort_thickness / nu_ox}") - return ( - (h_ox, y_ox, c_ox, rho_ox, u_ox, nu_ox), - (h_fu, y_fu, c_fu, rho_fu, u_fu, nu_fu), - z_st, - ) - - -def equilibrium_profiles(sol, pres, h_ox, h_fu, y_ox, y_fu, mixture_fraction): - """Equilibrium flamelet state at every mixture-fraction point.""" - num_z = len(mixture_fraction) - linear_enthalpy = h_ox + (h_fu - h_ox) * mixture_fraction - linear_mass_frac = (y_ox + (y_fu - y_ox) * mixture_fraction[:, None]).T - - state_equil = FlameletState( - enthalpy=linear_enthalpy, - mass_fractions=linear_mass_frac, - ) - temp_equil = jnp.zeros(num_z, dtype=jnp.float64) - for i in range(num_z): - sol.HPY = (state_equil.enthalpy[i], pres, state_equil.mass_fractions[:, i]) - sol.equilibrate("HP") - temp_equil = temp_equil.at[i].set(sol.T) - state_equil.mass_fractions = state_equil.mass_fractions.at[:, i].set(sol.Y) - - print(f"Equilibrium profiles: Max temperature T_max = {temp_equil.max():.3f} [K]") - return state_equil, temp_equil - - -def baseline_s_curve_state(z_st, pres, h_ox, h_fu, state_guess, temp_guess, mixture_fraction, solver_options, solver): - """Solve the flamelet at two high dissipation rates (also JIT-warms the solver).""" - - def _diss(val_st): - return new_dissipation_profile(mixture_fraction, z_st, val_st, 2 * solver.domain.jac[0] ** 2) - - def _visc(): - return new_dissipation_profile(mixture_fraction, z_st, 0, 0) - - t0 = time.time() - wmp = solver.warmup( - "flamelet_newton_step", - state_guess, - _diss(100), - _visc(), - temp_guess, - pres, - h_ox, - h_fu, - ) - wmp.enthalpy.block_until_ready() - print(f"Warm up time: {(time.time() - t0):.4e} s") - - t0 = time.time() - wmp, _, _ = solver.warmup( - "flamelet_time_step", - state_guess, - solver_options["bdf"]["newton"]["maxiter"], - solver_options["bdf"]["newton"]["tol"], - state_guess, - solver_options["bdf"]["time_step"], - _diss(100), - _visc(), - temp_guess, - pres, - h_ox, - h_fu, - ) - wmp.enthalpy.block_until_ready() - print(f"Warm up time: {(time.time() - t0):.4e} s") - - def _solve(val_st, try_newton, t_in, s_in): - return solver.solve( - True, - solver_options["newton"]["maxiter"], - solver_options["newton"]["tol"], - solver_options["bdf"]["newton"]["maxiter"], - solver_options["bdf"]["newton"]["tol"], - solver_options["bdf"]["time_step"], - solver_options["bdf"]["maxsteps"], - try_newton, - solver_options["max_attempts"], - _diss(val_st), - _visc(), - t_in, - pres, - h_ox, - h_fu, - s_in, - ) - - t0 = time.time() - state_base, temp_base = _solve(100, False, temp_guess, state_guess) - state_base.enthalpy.block_until_ready() - print(f"Run time: {(time.time() - t0):.4e} s") - - t0 = time.time() - # Was 1000 -- for this case's stream temperatures (T_ox=500K, T_fu=300K) and - # sandiego.yaml mechanism, chi=1000 1/s stalls the Newton solve (residual oscillates, - # never converges): it's past extinction for this fuel/oxidizer pair, not a numerics - # bug. This value is pure JIT-warmup scaffolding (state_base/temp_base's return here - # is discarded -- see this function's docstring), not a physical target, so any - # comfortably-converging value works; 150 is safely below both this case's actual - # chi_st (~286 1/s at vort_thickness=1e-3, from find_flame_dissipation_rate) and - # extinction. - state_base, temp_base = _solve(150, True, temp_base, state_base) - state_base.enthalpy.block_until_ready() - print(f"Run time: {(time.time() - t0):.4e} s") - - return state_base, temp_base - - -def find_flame_dissipation_rate(pyro_gas, z_st, mixture_fraction, cross_coord, num_iter, sim_fields_cold, state_guess, temp_guess, pres, h_ox, h_fu, solver_options, solver): - """Iterate the scalar dissipation rate to match the cold mixture-fraction field.""" - cold_diff = diffusivity( - pyro_gas, - sim_fields_cold.pressure, - sim_fields_cold.temperature, - sim_fields_cold.mass_fractions, - ) - dz_dy = np.gradient(sim_fields_cold.mixture_fraction, cross_coord) - diss_rate = 2 * cold_diff * (dz_dy**2) - - sim_fields = SimulationFields( - mixture_fraction=sim_fields_cold.mixture_fraction, - velocity=sim_fields_cold.velocity, - pressure=sim_fields_cold.pressure, - temperature=sim_fields_cold.temperature, - mass_fractions=sim_fields_cold.mass_fractions, - ) - state_it = state_guess - temp_it = temp_guess - try_newton = False - - for j_chi in range(num_iter): - _diss = interp1d(sim_fields.mixture_fraction, diss_rate, fill_value="extrapolate")(mixture_fraction) - _diss[0] = 2 * solver.domain.jac[0] ** 2 - _diss[-1] = 2 * solver.domain.jac[0] ** 2 - - state_it, temp_it = solver.solve( - True, - solver_options["newton"]["maxiter"], - solver_options["newton"]["tol"], - solver_options["bdf"]["newton"]["maxiter"], - solver_options["bdf"]["newton"]["tol"], - solver_options["bdf"]["time_step"], - solver_options["bdf"]["maxsteps"], - try_newton, - solver_options["max_attempts"], - _diss, - jnp.zeros_like(_diss), - temp_it, - pres, - h_ox, - h_fu, - state_it, - ) - if not try_newton: - try_newton = True - - sim_fields.temperature = interp1d(mixture_fraction, temp_it)(sim_fields.mixture_fraction) - sim_fields.mass_fractions = interp1d(mixture_fraction, state_it.mass_fractions)(sim_fields.mixture_fraction) - - diff = diffusivity(pyro_gas, pres, sim_fields.temperature, sim_fields.mass_fractions) - diss_rate = 2 * diff * (dz_dy**2) - chi_st = interp1d(sim_fields.mixture_fraction, diss_rate)(z_st) - print(f"---> j = {j_chi}, chi_max = {diss_rate.max():.4f}, " f"chi_st = {chi_st:.4f}, T_max = {temp_it.max():.4f}") - - return sim_fields +# Also support direct `python case.py` from a source checkout. +sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "toolchain")) +from mfc.flamelet import ( + INITIALIZER_VERSION, + create_simulation_fields, + density, + ic_cache_valid, + mechanism_fingerprint, + reference_fluid_properties, + streams, + write_cache_key, +) def compute_grid_3d(vort_thickness, cross_min, cross_max, points_per_cross, stream_min, stream_max, num_x, span_min, span_max, num_z): - """Pure grid arithmetic -- no Cantera/JAX. Always cheap to call, including on an + """Pure grid arithmetic -- no Cantera. Always cheap to call, including on an IC/ cache hit, so case.py never needs to run the flamelet solve just to learn its own domain size. @@ -328,105 +68,6 @@ def compute_grid_3d(vort_thickness, cross_min, cross_max, points_per_cross, stre return cross_coord, x_coord, grid -def reference_fluid_properties(sol, temperature_ox, pressure, mole_fraction_ox): - """Cheap (no equilibration/solve) reference gamma and viscosity for fluid_pp(1).""" - sol.TPX = temperature_ox, pressure, f"O2:{mole_fraction_ox}, N2:{1 - mole_fraction_ox}" - return {"gamma": float(sol.cp_mass / sol.cv_mass), "viscosity": float(sol.viscosity)} - - -def ic_cache_valid(ic_dir, file_extension, expected_lines, cache_key=None): - """True only if IC/ has prim.1..dat with exactly the current grid's expected - line count AND (if `cache_key` is given) a matching .cache_key.json. See - examples/2D_reacting_mixing_layer/flamelet_ic.py's docstring for why both checks - are needed.""" - path = os.path.join(ic_dir, f"prim.1.00.{file_extension}.dat") - if not os.path.isfile(path): - return False - with open(path) as fh: - if sum(1 for _ in fh) != expected_lines: - return False - if cache_key is not None: - key_path = os.path.join(ic_dir, ".cache_key.json") - if not os.path.isfile(key_path): - return False - try: - with open(key_path) as fh: - stored = json.load(fh) - except (OSError, ValueError): - return False - if stored != cache_key: - return False - return True - - -def write_cache_key(ic_dir, cache_key): - """Record the parameters an IC/ was generated with, so ic_cache_valid can detect a - stale cache (different --hot/cold mode or physical parameters at the same grid size).""" - with open(os.path.join(ic_dir, ".cache_key.json"), "w") as fh: - json.dump(cache_key, fh, sort_keys=True) - - -def create_simulation_fields(pyro_gas, sol, pres, temp_ox, temp_fu, cross_coord, vort_thickness, stream_ox, stream_fu, z_st, num_iter, cold): - """1-D flamelet profiles on the cross-stream (MFC y) coordinate `cross_coord`.""" - mollifier = 0.5 * (1 - np.tanh(2 * cross_coord / vort_thickness)) - - h_ox, y_ox, _, _, u_ox, nu_ox = stream_ox - h_fu, y_fu, _, _, u_fu, _ = stream_fu - - sim_fields_cold = SimulationFields( - mixture_fraction=mollifier, - temperature=pyro_gas.get_temperature_from_enthalpy( - h_ox + (h_fu - h_ox) * mollifier, - (y_ox + (y_fu - y_ox) * mollifier[:, None]).T, - temp_ox + (temp_fu - temp_ox) * mollifier, - ), - pressure=pres * jnp.ones_like(mollifier), - velocity=(u_ox + (u_fu - u_ox) * mollifier), - mass_fractions=(y_ox + (y_fu - y_ox) * mollifier[:, None]).T, - ) - - if cold: - return sim_fields_cold - - domain_config = DomainConfig(num_x=101, x_l=0, x_r=1) - domain = Domain(domain_config) - solver = FlameletSolver(domain, pyro_gas, (y_ox, y_fu)) - mixture_fraction = jnp.array(domain.x) - - state_equil, temp_equil = equilibrium_profiles(sol, pres, h_ox, h_fu, y_ox, y_fu, mixture_fraction) - solver_options = configure_flamelet_solver() - # baseline_s_curve_state's return value is unused -- its purpose is JIT warmup - # for the JAX solver calls that follow (matches the upstream tool's structure). - baseline_s_curve_state( - z_st, - pres, - h_ox, - h_fu, - state_equil, - temp_equil, - mixture_fraction, - solver_options, - solver, - ) - - sim_fields = find_flame_dissipation_rate( - pyro_gas, - z_st, - mixture_fraction, - cross_coord, - num_iter, - sim_fields_cold, - state_equil, - temp_equil, - pres, - h_ox, - h_fu, - solver_options, - solver, - ) - return sim_fields - - def perturb_xy(x_coord, cross_coord, vort_thickness, delta_u, seed, num_modes=10, num_blocks=5): """Solenoidal (x,y) velocity perturbation. Wavenumbers are computed relative to domain_length (x) and vort_thickness (y), not a hardcoded absolute range, so the @@ -534,7 +175,6 @@ def generate_ic_files( *, output_dir, sol, - pyro_gas, cross_coord, x_coord, pressure, @@ -545,7 +185,7 @@ def generate_ic_files( mole_fraction_fu, vort_thickness, mach_c, - num_iter, + strain_rate, cold, perturb_seed, file_extension="000000", @@ -554,11 +194,11 @@ def generate_ic_files( files on `cross_coord`/`x_coord` (from `compute_grid_3d`, so the file spacing exactly matches the grid case.py declares). - All Cantera/JAX/Pyrometheus stdout diagnostics are redirected to stderr: case.py's + All Cantera stdout diagnostics are redirected to stderr: case.py's contract requires its entire stdout to be exactly one JSON line. """ with contextlib.redirect_stdout(sys.stderr): - stream_ox, stream_fu, z_st = streams( + stream_ox, stream_fu, _ = streams( sol, fuel, pressure, @@ -571,7 +211,6 @@ def generate_ic_files( ) sim_fields = create_simulation_fields( - pyro_gas, sol, pressure, temperature_ox, @@ -580,8 +219,7 @@ def generate_ic_files( vort_thickness, stream_ox, stream_fu, - z_st, - num_iter, + strain_rate, cold, ) @@ -589,7 +227,7 @@ def generate_ic_files( pressure_1d = np.array(sim_fields.pressure) velocity_1d = np.array(sim_fields.velocity) mass_fractions_1d = np.array(sim_fields.mass_fractions) - density_1d = np.array(pyro_gas.get_density(pressure_1d, temperature_1d, mass_fractions_1d)) + density_1d = np.array(density(sol, pressure_1d, temperature_1d, mass_fractions_1d)) # Fail at generation time rather than writing a non-finite IC that would only # surface downstream as a cryptic VCFL=Inf crash (e.g. a diverged --hot solve). diff --git a/toolchain/mfc/build.py b/toolchain/mfc/build.py index eff5decaa..307cebf58 100644 --- a/toolchain/mfc/build.py +++ b/toolchain/mfc/build.py @@ -18,6 +18,7 @@ from .printer import cons from .run import input from .state import ARG, CFG, gpuConfigOptions +from .thermochem.fingerprint import generator_fingerprint, mechanism_fingerprint # Regex to parse build progress # Ninja format: [42/156] Building Fortran object ... @@ -305,9 +306,10 @@ def get_slug(self, case: Case) -> str: m.update(self.name.encode()) m.update(CFG().make_slug().encode()) m.update(case.get_fpp(self, False).encode()) + m.update(generator_fingerprint().encode()) if case.params.get("chemistry", "F") == "T": - m.update(case.get_cantera_solution().name.encode()) + m.update(mechanism_fingerprint(case.get_cantera_solution()).encode()) cfg = CFG() if cfg.gpu == gpuConfigOptions.ACC.value: diff --git a/toolchain/mfc/common.py b/toolchain/mfc/common.py index c59b5d086..b1a3d3fca 100644 --- a/toolchain/mfc/common.py +++ b/toolchain/mfc/common.py @@ -190,7 +190,7 @@ def get_py_program_output(filepath: str, arguments: typing.List[str] = None): """Run a case file and capture its stdout. sys.executable, not a bare "python3": a case file imports the same optional deps the - toolchain venv provides (cantera, pyrometheus, scipy, ...), so it must run under the + toolchain venv provides (cantera, numpy, ...), so it must run under the interpreter the toolchain itself is running under. Nearly every entry point activates the venv first, which makes PATH's python3 the venv's -- but a caller that invokes build/venv/bin/python3 DIRECTLY does not, and then case files silently ran under the diff --git a/toolchain/mfc/flamelet.py b/toolchain/mfc/flamelet.py new file mode 100644 index 000000000..0cec70c5f --- /dev/null +++ b/toolchain/mfc/flamelet.py @@ -0,0 +1,182 @@ +"""Cantera-only initial profiles shared by MFC's reacting mixing-layer examples. + +Hot profiles use a unity-Lewis counterflow diffusion flame tabulated against +Bilger mixture fraction, then mapped onto the prescribed tanh mixing layer. +This replaces the former self-consistent scalar-dissipation flamelet solve; +the counterflow's nominal strain rate is an explicit initialization parameter. +""" + +import json +import os +from dataclasses import dataclass + +import cantera as ct +import numpy as np + +from .thermochem.fingerprint import mechanism_fingerprint # noqa: F401 + +INITIALIZER_VERSION = "cantera-counterflow-v1" + + +@dataclass +class SimulationFields: + mixture_fraction: np.ndarray + temperature: np.ndarray + pressure: np.ndarray + velocity: np.ndarray + mass_fractions: np.ndarray + + +def stoichiometric_mixture_fraction(sol, y_ox, y_fu): + return 1.0 / (1.0 + sol.stoich_air_fuel_ratio(y_fu, y_ox, basis="mass")) + + +def density(sol, pressure, temperature, mass_fractions): + inv_weight = np.sum(mass_fractions / sol.molecular_weights[:, None], axis=0) + return pressure / (ct.gas_constant * temperature * inv_weight) + + +def _hot_mass_fractions(sol, pres, temp_ox, temp_fu, y_ox, y_fu, width, strain_rate, mixture_fraction): + if not np.isfinite(strain_rate) or strain_rate <= 0: + raise ValueError("Flame initialization requires a positive, finite strain rate") + original_state = sol.TPY + original_transport = sol.transport_model + try: + flame = ct.CounterflowDiffusionFlame(sol, width=width) + flame.P = pres + flame.fuel_inlet.T = temp_fu + flame.fuel_inlet.Y = y_fu + flame.oxidizer_inlet.T = temp_ox + flame.oxidizer_inlet.Y = y_ox + inlet_speed = strain_rate * width / 2 + sol.TPY = temp_fu, pres, y_fu + flame.fuel_inlet.mdot = sol.density * inlet_speed + sol.TPY = temp_ox, pres, y_ox + flame.oxidizer_inlet.mdot = sol.density * inlet_speed + flame.transport_model = "unity-Lewis-number" + flame.set_refine_criteria(ratio=3, slope=0.1, curve=0.2) + flame.solve(loglevel=0, auto=True) + if flame.extinct(): + raise ValueError("Hot initialization converged to an extinguished flame; reduce flame_strain_rate") + # Normalize against the inlet streams, not the flame-domain extrema: + # diffusion can make the boundary composition differ from the inlet. + z = np.empty(len(flame.grid)) + flame_y = flame.Y + for i, y in enumerate(flame_y.T): + sol.Y = y + z[i] = sol.mixture_fraction(y_fu, y_ox, basis="mass") + if np.any(np.diff(z) > 1e-7): + raise ValueError("Counterflow mixture fraction is not monotone; cannot map the flame profile") + # Supply exact stream endpoints and discard roundoff outside their range. + interior = (z > 1e-12) & (z < 1 - 1e-12) + z = np.concatenate(([0.0], z[interior], [1.0])) + y = np.column_stack((y_ox, flame_y[:, interior], y_fu)) + z, indices = np.unique(z, return_index=True) + y = y[:, indices] + mapped = np.array([np.interp(mixture_fraction, z, species) for species in y]) + if np.min(mapped) < -1e-8: + raise ValueError("Counterflow solve produced negative mass fractions") + mapped = np.maximum(mapped, 0) + return mapped / mapped.sum(axis=0) + finally: + sol.transport_model = original_transport + sol.TPY = original_state + + +def create_simulation_fields(sol, pres, temp_ox, temp_fu, cross_coord, vort_thickness, stream_ox, stream_fu, strain_rate, cold): + """Map cold enthalpy mixing or a hot counterflow flame onto the MFC grid.""" + if not np.isfinite(vort_thickness) or vort_thickness <= 0: + raise ValueError("Vorticity thickness must be positive and finite") + z = 0.5 * (1 - np.tanh(2 * cross_coord / vort_thickness)) + h_ox, y_ox, _, _, u_ox, _ = stream_ox + h_fu, y_fu, _, _, u_fu, _ = stream_fu + y = y_ox[:, None] + (y_fu - y_ox)[:, None] * z + if not cold: + y = _hot_mass_fractions(sol, pres, temp_ox, temp_fu, y_ox, y_fu, 10 * vort_thickness, strain_rate, z) + h = h_ox + (h_fu - h_ox) * z + temperature = np.empty_like(z) + for i in range(len(z)): + sol.HPY = h[i], pres, y[:, i] + temperature[i] = sol.T + if not np.all(np.isfinite(temperature)) or np.any(temperature <= 0): + raise ValueError("Invalid temperature in mixing-layer initialization") + return SimulationFields(z, temperature, np.full_like(z, pres), u_ox + (u_fu - u_ox) * z, y) + + +def streams(sol, fuel, pres, temp_ox, temp_fu, molefrac_ox, molefrac_fu, vort_thickness, mach_c, mach_fu=None): + """Thermodynamic state and velocities for the oxidizer and fuel streams + (temporal evolution: symmetric convective frame, u_ox = -u_fu).""" + molefrac_di = 1 - molefrac_ox + sol.TPX = temp_ox, pres, f"O2:{molefrac_ox}, N2:{molefrac_di}" + y_ox = sol.Y + h_ox = sol.enthalpy_mass + c_ox = np.sqrt((sol.cp_mass / sol.cv_mass) * sol.P / sol.density) + nu_ox = sol.viscosity / sol.density + rho_ox = sol.density + + molefrac_di = 1 - molefrac_fu + sol.TPX = temp_fu, pres, f"{fuel}:{molefrac_fu}, N2:{molefrac_di}" + y_fu = sol.Y + h_fu = sol.enthalpy_mass + c_fu = np.sqrt((sol.cp_mass / sol.cv_mass) * sol.P / sol.density) + nu_fu = sol.viscosity / sol.density + rho_fu = sol.density + + z_st = stoichiometric_mixture_fraction(sol, y_ox, y_fu) + print(f"Stoichiometric mixture fraction: Z_st = {z_st:.3f}") + + u_ox = 0.5 * mach_c * (c_ox + c_fu) if mach_fu is None else mach_c * c_ox + u_fu = -u_ox if mach_fu is None else mach_fu * c_fu + delta_u = u_ox - u_fu + + print(f"Convective Mach: Ma = {delta_u / (c_ox + c_fu)}") + print(f"Reynolds number: Re = {0.5 * delta_u * vort_thickness / nu_ox}") + return ( + (h_ox, y_ox, c_ox, rho_ox, u_ox, nu_ox), + (h_fu, y_fu, c_fu, rho_fu, u_fu, nu_fu), + z_st, + ) + + +def reference_fluid_properties(sol, temperature_ox, pressure, mole_fraction_ox): + """Cheap (no equilibration/solve) reference gamma and viscosity for fluid_pp(1).""" + sol.TPX = temperature_ox, pressure, f"O2:{mole_fraction_ox}, N2:{1 - mole_fraction_ox}" + return {"gamma": float(sol.cp_mass / sol.cv_mass), "viscosity": float(sol.viscosity)} + + +def ic_cache_valid(ic_dir, file_extension, expected_lines, cache_key=None): + """True only if IC/ has prim.1..dat with exactly the current grid's expected + line count AND (if `cache_key` is given) a matching .cache_key.json. A bare "IC/ is + non-empty" check isn't enough: this same case.py is invoked with different --scale + values by different toolchain paths (e.g. `./mfc.sh validate` during precheck uses no + args/default scale, while a registered test passes its own --scale) that can share + this directory, so a cache populated by one grid size must not be silently reused by a + run expecting a different one -- reading it would desync the Fortran reader (hcid=273 + expects exactly len(cross_coord) lines) from the actual grid. The `cache_key` further + guards against silently reusing an IC generated with a different mode (--hot vs cold) + or different physical parameters that leave the line count unchanged.""" + path = os.path.join(ic_dir, f"prim.1.00.{file_extension}.dat") + if not os.path.isfile(path): + return False + with open(path) as fh: + if sum(1 for _ in fh) != expected_lines: + return False + if cache_key is not None: + key_path = os.path.join(ic_dir, ".cache_key.json") + if not os.path.isfile(key_path): + return False + try: + with open(key_path) as fh: + stored = json.load(fh) + except (OSError, ValueError): + return False + if stored != cache_key: + return False + return True + + +def write_cache_key(ic_dir, cache_key): + """Record the parameters an IC/ was generated with, so ic_cache_valid can detect a + stale cache (different --hot/cold mode or physical parameters at the same grid size).""" + with open(os.path.join(ic_dir, ".cache_key.json"), "w") as fh: + json.dump(cache_key, fh, sort_keys=True) diff --git a/toolchain/mfc/run/input.py b/toolchain/mfc/run/input.py index dd3a122cf..a1fcb2677 100644 --- a/toolchain/mfc/run/input.py +++ b/toolchain/mfc/run/input.py @@ -7,7 +7,7 @@ from .. import case_validator, common from ..case import Case -# Note: pyrometheus and cantera are imported lazily in the methods that need them +# Note: thermochemistry generation and cantera are imported lazily where needed # to avoid slow startup times for commands that don't use chemistry features # Note: build is imported lazily to avoid circular import with build.py from ..printer import cons @@ -68,7 +68,7 @@ def get_cantera_solution(self): def generate_fpp(self, target) -> None: # Lazy import to avoid slow startup for commands that don't need chemistry - import pyrometheus as pyro + from ..thermochem import generate_fortran if target.isDependency: return @@ -83,8 +83,8 @@ def generate_fpp(self, target) -> None: modules_dir = os.path.join(target.get_staging_dirpath(self), "modules", target.name) common.create_directory(modules_dir) - # Determine the real type based on the single precision flag - real_type = "real(sp)" if (ARG("single") or ARG("mixed")) else "real(dp)" + # Match wp in m_precision_select; --mixed changes storage precision only. + real_type = "real(sp)" if ARG("single") else "real(dp)" if ARG("gpu") == gpuConfigOptions.MP.value: directive_str = "mp" @@ -96,7 +96,7 @@ def generate_fpp(self, target) -> None: # Write the generated Fortran code to the m_thermochem.f90 file with the chosen precision sol = self.get_cantera_solution() - thermochem_code = pyro.FortranCodeGenerator().generate("m_thermochem", sol, pyro.CodeGenerationOptions(scalar_type=real_type, directive_offload=directive_str)) + thermochem_code = generate_fortran(sol, scalar_type=real_type, offload=directive_str) common.file_write(os.path.join(modules_dir, "m_thermochem.f90"), thermochem_code, True) diff --git a/toolchain/mfc/test_flamelet.py b/toolchain/mfc/test_flamelet.py new file mode 100644 index 000000000..010076d45 --- /dev/null +++ b/toolchain/mfc/test_flamelet.py @@ -0,0 +1,54 @@ +"""Physical invariants of the Cantera-only mixing-layer initialization.""" + +import cantera as ct +import numpy as np +import pytest + +from mfc.flamelet import create_simulation_fields, density, mechanism_fingerprint, streams +from mfc.test_thermochem import ROOT + + +@pytest.mark.parametrize("cold", [True, False]) +@pytest.mark.parametrize("mechanism,fuel_fraction", [("h2o2.yaml", 1.0), (str(ROOT / "examples/3D_reacting_mixing_layer/sandiego.yaml"), 0.5)]) +def test_mixing_layer_profiles(mechanism, fuel_fraction, cold): + gas = ct.Solution(mechanism) + pressure = ct.one_atm + oxidizer, fuel, _ = streams(gas, "H2", pressure, 500, 300, 0.21, fuel_fraction, 0.001, 0.3) + coords = np.linspace(-0.01, 0.01, 401) + fields = create_simulation_fields(gas, pressure, 500, 300, coords, 0.001, oxidizer, fuel, 100, cold) + y = fields.mass_fractions + assert np.all(np.isfinite(y)) and np.min(y) >= 0 + np.testing.assert_allclose(y.sum(axis=0), 1, atol=1e-12) + np.testing.assert_allclose(y[:, 0], fuel[1], atol=1e-12) + np.testing.assert_allclose(y[:, -1], oxidizer[1], atol=1e-12) + np.testing.assert_allclose(fields.temperature[[0, -1]], [300, 500], atol=1e-6) + np.testing.assert_allclose(fields.velocity[[0, -1]], [fuel[4], oxidizer[4]], atol=1e-8) + rho = density(gas, fields.pressure, fields.temperature, y) + enthalpy = oxidizer[0] + (fuel[0] - oxidizer[0]) * fields.mixture_fraction + atoms = np.array([[gas.n_atoms(k, el) * gas.atomic_weight(el) / gas.molecular_weights[k] for k in range(gas.n_species)] for el in gas.element_names]) + mixed_y = oxidizer[1][:, None] + (fuel[1] - oxidizer[1])[:, None] * fields.mixture_fraction + # Unity-Lewis counterflow conserves the stream elemental mixture fractions. + np.testing.assert_allclose(atoms @ y, atoms @ mixed_y, atol=2e-5) + for i in range(len(coords)): + gas.TPY = fields.temperature[i], pressure, y[:, i] + assert gas.density == pytest.approx(rho[i], rel=1e-12) + # Cantera's HP inversion tolerance is scaled by cp*T, including near h=0. + assert gas.enthalpy_mass == pytest.approx(enthalpy[i], abs=1e-8 * gas.cp_mass * gas.T) + if cold: + np.testing.assert_allclose(y, mixed_y, atol=1e-14) + assert fields.temperature.max() <= 500 + 1e-6 + else: + assert fields.temperature.max() > 1500 + assert y[gas.species_index("H2O")].max() > 0.05 + + +def test_mechanism_fingerprint_ignores_state_but_tracks_rates(): + gas = ct.Solution("h2o2.yaml") + original = mechanism_fingerprint(gas) + gas.TPX = 1500, 2 * ct.one_atm, "H2:1" + assert mechanism_fingerprint(gas) == original + reaction = gas.reaction(2) + rate = reaction.rate + reaction.rate = ct.ArrheniusRate(2 * rate.pre_exponential_factor, rate.temperature_exponent, rate.activation_energy) + gas.modify_reaction(2, reaction) + assert mechanism_fingerprint(gas) != original diff --git a/toolchain/mfc/test_thermochem.py b/toolchain/mfc/test_thermochem.py new file mode 100644 index 000000000..17a0f5f34 --- /dev/null +++ b/toolchain/mfc/test_thermochem.py @@ -0,0 +1,225 @@ +"""Compile generated kernels and compare their numerical interface with Cantera.""" + +import shutil +import subprocess +from pathlib import Path +from types import SimpleNamespace + +import cantera as ct +import numpy as np +import pytest + +from mfc.thermochem import generate_fortran + +ROOT = Path(__file__).resolve().parents[2] +MECHANISMS = [ + "h2o2.yaml", + "gri30.yaml", + str(ROOT / "examples/3D_reacting_mixing_layer/sandiego.yaml"), + str(ROOT / "examples/2D_reactive_shock_bubble/h2o2_xe.yaml"), +] + +DRIVER = """ +program reference + use m_thermochem + implicit none + integer :: ierr + real(KIND) :: t, pressure, rho, cp, cv, h, e, mw, mu, lambda, te, th, pcheck + real(KIND) :: y(num_species), omega(num_species), creation(num_species), destruction(num_species) + real(KIND) :: diffusion(num_species), enthalpy(num_species) + do + read(*,*,iostat=ierr) t, pressure, y + if (ierr /= 0) exit + call get_density(pressure, t, y, rho) + call get_pressure(rho, t, y, pcheck) + call get_mixture_specific_heat_cp_mass(t, y, cp) + call get_mixture_specific_heat_cv_mass(t, y, cv) + call get_mixture_enthalpy_mass(t, y, h) + call get_mixture_energy_mass(t, y, e) + call get_mixture_molecular_weight(y, mw) + call get_temperature(e, 0.9_KIND*t, y, .true., te) + call get_temperature(h, 1.1_KIND*t, y, .false., th) + call get_mixture_viscosity_mixavg(t, y, mu) + call get_mixture_thermal_conductivity_mixavg(t, y, lambda) + call get_species_mass_diffusivities_mixavg(pressure, t, y, diffusion) + call get_species_enthalpies_rt(t, enthalpy) + call get_net_production_rates(rho, t, y, omega) + call get_creation_destruction_rates(rho, t, y, creation, destruction) + write(*,'(*(ES25.16E3,1X))') rho, pcheck, cp, cv, h, e, mw, te, th, mu, lambda, & + diffusion, enthalpy, omega, creation, destruction + end do +end program +""" + + +def compile_kernel(directory, gas, precision="dp", offload=None, *, source=None, driver_source=DRIVER, extra_flags=(), extra_sources=()): + compiler = shutil.which("gfortran") + if compiler is None: + pytest.skip("gfortran is required to validate generated Fortran") + module = directory / "m_thermochem.f90" + module.write_text(source if source is not None else generate_fortran(gas, scalar_type=f"real({precision})", offload=offload)) + driver = directory / "driver.f90" + driver.write_text(driver_source.replace("KIND", precision)) + executable = directory / "reference" + flags = {None: [], "acc": ["-fopenacc"], "mp": ["-fopenmp"]}[offload] + subprocess.run( + [compiler, "-cpp", "-O0", "-Wconversion", "-Werror=conversion", *flags, *extra_flags, *map(str, extra_sources), str(module), str(driver), "-o", str(executable)], + cwd=directory, + check=True, + capture_output=True, + text=True, + ) + return executable + + +def reference_states(gas): + rng = np.random.default_rng(23) + for temperature, pressure in [(300, ct.one_atm), (999, 0.1 * ct.one_atm), (1001, ct.one_atm), (1800, 10 * ct.one_atm), (2800, ct.one_atm)]: + y = rng.uniform(0.01, 1, gas.n_species) + gas.TPY = temperature, pressure, y / y.sum() + yield gas.T, gas.P, gas.Y + # Zero concentrations and the pure-species diffusion fallback. + gas.TPX = 800, ct.one_atm, "N2:1" + yield gas.T, gas.P, gas.Y + gas.TPX = 1200, ct.one_atm, "H2:2,O2:1,N2:3.76" + yield gas.T, gas.P, gas.Y + + +def compare_kernel(executable, gas, precision="dp"): + states = list(reference_states(gas)) + inputs = "\n".join(" ".join(map(str, [t, p, *y])) for t, p, y in states) + "\n" + result = subprocess.run([str(executable)], input=inputs, capture_output=True, text=True, check=True) + rows = np.array([np.fromstring(line, sep=" ") for line in result.stdout.splitlines()]) + assert len(rows) == len(states) + rtol = 3e-5 if precision == "sp" else 2e-11 + for actual, (t, p, y) in zip(rows, states): + gas.TPY = t, p, y + diffusion = gas.mix_diff_coeffs.copy() + # Preserve the old generator's pure-species self-diffusion convention. + # Cantera's mix_diff_coeffs may instead return zero in this degenerate limit. + pure = np.flatnonzero(y == 1) + for k in pure: + diffusion[k] = gas.binary_diff_coeffs[k, k] + expected = np.concatenate( + ( + [gas.density, gas.P, gas.cp_mass, gas.cv_mass, gas.enthalpy_mass, gas.int_energy_mass, gas.mean_molecular_weight, gas.T, gas.T, gas.viscosity, gas.thermal_conductivity], + diffusion, + gas.standard_enthalpies_RT, + gas.net_production_rates, + gas.creation_rates, + gas.destruction_rates, + ) + ) + np.testing.assert_allclose(actual, expected, rtol=rtol, atol=2e-8 if precision == "sp" else 1e-10) + ns = gas.n_species + net, creation, destruction = actual[11 + 2 * ns :].reshape(3, ns) + np.testing.assert_allclose(net, creation - destruction, rtol=rtol, atol=rtol * max(1, np.max(creation))) + elements = np.array([[gas.n_atoms(k, el) for k in range(ns)] for el in gas.element_names]) + residual = elements @ net + assert np.max(np.abs(residual)) < rtol * max(1, np.max(np.abs(net))) + + +@pytest.mark.parametrize("mechanism", MECHANISMS) +def test_generated_mechanisms(tmp_path, mechanism): + gas = ct.Solution(mechanism) + compare_kernel(compile_kernel(tmp_path, gas), gas) + + +@pytest.mark.parametrize("precision,offload", [("sp", None), ("dp", "acc"), ("dp", "mp")]) +def test_precision_and_directives(tmp_path, precision, offload): + gas = ct.Solution("h2o2.yaml") + compare_kernel(compile_kernel(tmp_path, gas, precision, offload), gas, precision) + + +def test_rejects_unsupported_thermo(): + species = ct.Species("N2", {"N": 2}) + species.thermo = ct.ConstantCp(200, 5000, ct.one_atm, [300, 0, 0, 30000]) + species.transport = ct.Solution("h2o2.yaml").species("N2").transport + gas = ct.Solution(thermo="ideal-gas", kinetics="gas", species=[species], transport_model="mixture-averaged") + with pytest.raises(ValueError, match="Species N2.*NASA7"): + generate_fortran(gas) + + +def test_rejects_custom_orders(): + gas = ct.Solution("h2o2.yaml") + reaction = ct.Reaction(equation="H2 + O => H + OH", rate=ct.ArrheniusRate(1e8, 0, 0)) + reaction.orders = {"H2": 0.5} + custom = ct.Solution(thermo="ideal-gas", kinetics="gas", species=gas.species(), reactions=[reaction], transport_model="mixture-averaged") + with pytest.raises(ValueError, match="Reaction 1.*custom reaction orders"): + generate_fortran(custom) + + +@pytest.mark.parametrize("mode", ["double", "single", "mixed"]) +def test_solver_working_precision(tmp_path, monkeypatch, mode): + from mfc.run import input as input_module + + monkeypatch.setattr(input_module, "ARG", lambda name: {"single": mode == "single", "mixed": mode == "mixed", "gpu": None}[name]) + case = input_module.MFCInputFile("case.py", str(tmp_path), {"chemistry": "T", "cantera_file": "h2o2.yaml"}) + monkeypatch.setattr(case, "get_fpp", lambda target: "") + target = SimpleNamespace(name="simulation", isDependency=False, get_staging_dirpath=lambda case: str(tmp_path)) + case.generate_fpp(target) + source = (tmp_path / "modules/simulation/m_thermochem.f90").read_text() + driver = DRIVER.replace("use m_thermochem", "use m_thermochem\n use m_precision_select, only: wp") + flags = [] if mode == "double" else [f"-DMFC_{mode.upper()}_PRECISION"] + gas = ct.Solution("h2o2.yaml") + executable = compile_kernel( + tmp_path, + gas, + "wp", + source=source, + driver_source=driver, + extra_flags=flags, + extra_sources=[ROOT / "src/common/m_precision_select.f90"], + ) + compare_kernel(executable, gas, "sp" if mode == "single" else "dp") + + +def test_long_species_names(tmp_path): + names = ["nitrogen_reference_one", "nitrogen_reference_two"] + nitrogen = ct.Solution("h2o2.yaml").species("N2") + species = [] + for name in names: + data = dict(nitrogen.input_data) + data["name"] = name + species.append(ct.Species.from_dict(data)) + gas = ct.Solution(thermo="ideal-gas", kinetics="gas", species=species, transport_model="mixture-averaged") + driver = """ +program names + use m_thermochem + implicit none + character(len=50) :: name + integer :: i, index + do i = 1, num_species + call get_species_name(i, name) + call get_species_index(name, index) + if (index /= i) stop 1 + print *, trim(name) + end do +end program +""" + executable = compile_kernel(tmp_path, gas, driver_source=driver) + result = subprocess.run([str(executable)], capture_output=True, text=True, check=True) + assert result.stdout.split() == names + + +@pytest.mark.parametrize("precision", ["sp", "dp"]) +def test_zero_concentration_falloff(tmp_path, precision): + gas = ct.Solution("h2o2.yaml") + driver = """ +program falloff + use m_thermochem + use, intrinsic :: ieee_arithmetic, only: ieee_is_finite + implicit none + real(KIND) :: concentrations(num_species), rates(num_reactions) + concentrations = 0.0_KIND + call get_fwd_rate_coefficients(1200.0_KIND, concentrations, rates) + if (.not. all(ieee_is_finite(rates))) stop 1 + print *, rates +end program +""" + executable = compile_kernel(tmp_path, gas, precision, driver_source=driver, extra_flags=["-ffpe-trap=invalid,zero,overflow"]) + result = subprocess.run([str(executable)], capture_output=True, text=True, check=True) + rates = np.fromstring(result.stdout, sep=" ") + for i, reaction in enumerate(gas.reactions()): + if isinstance(reaction.rate, ct.FalloffRate): + assert rates[i] == 0 diff --git a/toolchain/mfc/thermochem/LICENSE b/toolchain/mfc/thermochem/LICENSE new file mode 100644 index 000000000..99cdfa654 --- /dev/null +++ b/toolchain/mfc/thermochem/LICENSE @@ -0,0 +1,21 @@ +pyrometheus is licensed to you under the MIT/Expat/X-Consortium license. + +Copyright (c) 2020 University of Illinois Board of Trustees + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. diff --git a/toolchain/mfc/thermochem/__init__.py b/toolchain/mfc/thermochem/__init__.py new file mode 100644 index 000000000..cfa6d9095 --- /dev/null +++ b/toolchain/mfc/thermochem/__init__.py @@ -0,0 +1,10 @@ +"""MFC-owned, Fortran-only thermochemistry generation.""" + +__all__ = ["generate_fortran"] + + +def generate_fortran(solution, module_name="m_thermochem", scalar_type="real(dp)", offload=None): + """Load the expression machinery only when generating Fortran.""" + from .fortran import generate_fortran as generate + + return generate(solution, module_name, scalar_type, offload) diff --git a/toolchain/mfc/thermochem/expressions.py b/toolchain/mfc/thermochem/expressions.py new file mode 100644 index 000000000..06e5f4969 --- /dev/null +++ b/toolchain/mfc/thermochem/expressions.py @@ -0,0 +1,410 @@ +__copyright__ = """ +Copyright (C) 2020 University of Illinois Board of Trustees +""" + +__license__ = """ +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in +all copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN +THE SOFTWARE. +""" + +from functools import singledispatch + +import cantera as ct +import numpy as np +import pymbolic.primitives as p + +# {{{ polynomial processing + + +def nasa7_conditional(t, poly, part_gen): + # FIXME: Should check minTemp, maxTemp + return p.If( + p.Comparison(t, ">", poly.coeffs[0]), + part_gen(poly.coeffs[1:8], t), + part_gen(poly.coeffs[8:15], t), + ) + + +@singledispatch +def poly_to_expr(poly): + raise TypeError(f"unexpected argument type in poly_to_expr: {type(poly)}") + + +@poly_to_expr.register +def _(poly: ct.NasaPoly2, arg_name): + def gen(c, t): + assert len(c) == 7 + return c[0] + c[1] * t + c[2] * t**2 + c[3] * t**3 + c[4] * t**4 + + return nasa7_conditional(p.Variable(arg_name), poly, gen) + + +@singledispatch +def poly_to_enthalpy_expr(poly, arg_name): + raise TypeError("unexpected argument type in poly_to_enthalpy_expr: " f"{type(poly)}") + + +@poly_to_enthalpy_expr.register +def _(poly: ct.NasaPoly2, arg_name): + def gen(c, t): + assert len(c) == 7 + return c[0] + c[1] / 2 * t + c[2] / 3 * t**2 + c[3] / 4 * t**3 + c[4] / 5 * t**4 + c[5] / t + + return nasa7_conditional(p.Variable(arg_name), poly, gen) + + +@singledispatch +def poly_to_entropy_expr(poly, arg_name): + raise TypeError("unexpected argument type in poly_to_entropy_expr: " f"{type(poly)}") + + +@poly_to_entropy_expr.register +def _(poly: ct.NasaPoly2, arg_name): + log = p.Variable("log") + + def gen(c, t): + assert len(c) == 7 + return c[0] * log(t) + c[1] * t + c[2] / 2 * t**2 + c[3] / 3 * t**3 + c[4] / 4 * t**4 + c[6] + + return nasa7_conditional(p.Variable(arg_name), poly, gen) + + +# {{{ Data-handling helper + + +def _zeros_like(argument): + # FIXME: This mishandles NaNs. + return 0 * argument + + +# }}} + + +# {{{ Transport polynomials & mixture rules + + +def viscosity_polynomial_expr(c, t): + """Generate code for viscosity polynomials + + :returns: Viscosity polynomial expression with coefficients c in terms of + the temperature t as a :class:`pymbolic.primitives.Expression`. + """ + assert len(c) == 5 + return p.Variable("sqrt")(t) * (c[0] + c[1] * p.Variable("log")(t) + c[2] * p.Variable("log")(t) ** 2 + c[3] * p.Variable("log")(t) ** 3 + c[4] * p.Variable("log")(t) ** 4) ** 2 + + +def conductivity_polynomial_expr(c, t): + """Generate code for conductivity polynomials + + :returns: Conductivity polynomial expression with coefficients c in terms + of the temperature t as a :class:`pymbolic.primitives.Expression`. + """ + assert len(c) == 5 + return p.Variable("sqrt")(t) * (c[0] + c[1] * p.Variable("log")(t) + c[2] * p.Variable("log")(t) ** 2 + c[3] * p.Variable("log")(t) ** 3 + c[4] * p.Variable("log")(t) ** 4) + + +def diffusivity_polynomial_expr(c, t): + """Generate code for diffusivity polynomials + + :returns: Diffusivity polynomial expression with coefficients c in terms + of the temperature t as a :class:`pymbolic.primitives.Expression`. + """ + assert len(c) == 5 + return p.Variable("sqrt")(t) * t * (c[0] + c[1] * p.Variable("log")(t) + c[2] * p.Variable("log")(t) ** 2 + c[3] * p.Variable("log")(t) ** 3 + c[4] * p.Variable("log")(t) ** 4) + + +def viscosity_mixture_rule_wilke_expr(sol: ct.Solution, sp, x, mu): + """Generate code for species mixture rule. See [Kee_2003]_, chapter 12. + + :returns: Expression for the Wilke viscosity mixture rule + for species *sp* in terms of species mole fractions *w* + and viscosities *mu* as a :class:`pymbolic.primitives.Expression` + """ + w = sol.molecular_weights + sqrt = p.Variable("sqrt") + return sum([x[j] * (1 + sqrt((mu[sp] / mu[j]) * np.sqrt(w[j] / w[sp]))) ** 2 / np.sqrt(8 * (1 + (w[sp] / w[j]))) for j in range(sol.n_species)]) + + +def diffusivity_mixture_rule_denom_expr(sol: ct.Solution, j_sp, x, bdiff): + """See [Kee_2003]_, chapter 12 for details. + :returns: The denominator expression to the mixture rule + for mixture-averaged species diffusivities in terms + of the species mole fractions *x* and binary diffusivities *bdiff* as a + :class:`pymbolic.primitives.Expression` + """ + return sum(x[i_sp] / bdiff[i_sp][j_sp] for i_sp in range(sol.n_species)) + + +# }}} + + +# {{{ Equilibrium constants + + +def equilibrium_constants_expr(sol: ct.Solution, reaction_index, gibbs_rt): + """Generate code for equilibrium constants. + + :returns: Equilibrium constant expression for reaction with + index *reaction_index* in terms of the species Gibbs + functions *gibbs_rt* as a :class:`pymbolic.primitives.Expression` + """ + indices_reac = [sol.species_index(sp) for sp in sol.reaction(reaction_index).reactants] + indices_prod = [sol.species_index(sp) for sp in sol.reaction(reaction_index).products] + + # Stoichiometric coefficients + nu_reac = [sol.reactant_stoich_coeff(sol.species_index(sp), reaction_index) for sp in sol.reaction(reaction_index).reactants] + nu_prod = [sol.product_stoich_coeff(sol.species_index(sp), reaction_index) for sp in sol.reaction(reaction_index).products] + + sum_r = sum(nu_reac_i * gibbs_rt[indices_reac_i] for indices_reac_i, nu_reac_i in zip(indices_reac, nu_reac)) + sum_p = sum(nu_prod_i * gibbs_rt[indices_prod_i] for indices_prod_i, nu_prod_i in zip(indices_prod, nu_prod)) + + # Check if reaction is termolecular + sum_nu_net = sum(nu_prod) - sum(nu_reac) + if sum_nu_net != 0: + return sum_p - sum_r - sum_nu_net * p.Variable("c0") + else: + return sum_p - sum_r + + +# }}} + + +# {{{ Rate coefficients + + +def rate_coefficient_expr(rate_coeff: ct.Arrhenius, t): + """ + :returns: The rate coefficient expression for *rate_coeff* in terms + of the temperature *t* as a :class:`pymbolic.primitives.Expression` + """ + # Rate parameters + a = rate_coeff.pre_exponential_factor + b = rate_coeff.temperature_exponent + t_a = rate_coeff.activation_energy / ct.gas_constant + if t_a == 0: + # Weakly temperature-dependent rate + return a * t**b + else: + # Modified Arrhenius + return p.Variable("exp")(np.log(a) + b * p.Variable("log")(t) - t_a / t) + + +def third_body_efficiencies_expr(sol: ct.Solution, react: ct.Reaction, c): + """ + :returns: The third-body concentration expression for reaction *react* in + terms of the species concentrations *c* as a + :class:`pymbolic.primitives.Expression` + """ + + efficiencies = [react.third_body.efficiencies[sp] for sp in react.third_body.efficiencies] + indices_nondef = [sol.species_index(sp) for sp in react.third_body.efficiencies] + indices_default = [i for i in range(sol.n_species) if i not in indices_nondef] + sum_nondef = sum(eff_i * c[index_i] for eff_i, index_i in zip(np.array(efficiencies), indices_nondef)) + sum_default = react.third_body.default_efficiency * sum(c[i] for i in indices_default) + return sum_nondef + sum_default + + +def troe_falloff_center_expr(react: ct.Reaction, t): + """ + :returns: The Troe falloff center expression for reaction *react* in + terms of the temperature *t* as a + :class:`pymbolic.primitives.Expression` + """ + + if isinstance(react.rate, ct.TroeRate): + troe_params = react.rate.falloff_coeffs + elif isinstance(react.rate, ct.LindemannRate): + return 1 + else: + raise ValueError("Unexpected value of 'rate.type': " f" '{react.rate.type}'") + + troe_1 = (1.0 - troe_params[0]) * p.Variable("exp")(-t / troe_params[1]) + troe_2 = troe_params[0] * p.Variable("exp")(-t / troe_params[2]) + if len(troe_params) == 3: + return p.Variable("log10")(troe_1 + troe_2) + elif len(troe_params) == 4: + troe_3 = p.Variable("exp")(-troe_params[3] / t) + return p.Variable("log10")(troe_1 + troe_2 + troe_3) + else: + raise ValueError("Unexpected length of 'troe_params': " f" '{len(troe_params)}'") + return + + +def troe_falloff_factor_expr(react: ct.Reaction, i, red_pressure, falloff_center): + """ + :returns: The Troe falloff factor expression for reaction + *react* in terms of reduced pressure *red_pressure* and the + falloff center *falloff_center* as a + :class:`pymbolic.primitives.Expression` + + """ + if isinstance(react.rate, ct.TroeRate): + # Fortran MERGE need not short-circuit: guard the logarithm's argument. + positive = p.Comparison(red_pressure[i], ">", 0) + log_rp = p.Variable("log10")(p.If(positive, red_pressure[i], 1)) + c = -0.4 - 0.67 * falloff_center[i] + n = 0.75 - 1.27 * falloff_center[i] + return p.If(positive, (log_rp + c) / (n - 0.14 * (log_rp + c)), -1 / 0.14) + elif isinstance(react.rate, ct.LindemannRate): + return 0 + else: + raise ValueError("Unexpected value of 'rate.type': " f" '{react.rate.type}'") + + +def falloff_function_expr(react: ct.Reaction, i, falloff_factor, falloff_center): + """ + :returns: Falloff function expression for reaction *react* in + terms of the temperature *t*, falloff width factor + *falloff_factor*, and falloff center *falloff_center* as a + :class:`pymbolic.primitives.Expression` + + """ + + falloff_type = react.reaction_type.split("-")[1] + + if falloff_type == "Troe": + return 10 ** (falloff_center[i] / (1 + falloff_factor[i] ** 2)) + elif falloff_type == "Lindemann": + return 1 + else: + raise ValueError("Unexpected value of 'falloff_type': " f" '{falloff_type}'") + + +# }}} + + +# {{{ Rates of progress + + +def rate_of_progress_expr(sol: ct.Solution, reaction_index, c, k_fwd, log_k_eq): + """ + :returns: Rate of progress expression for reaction with + index *reaction_index* in terms of species concentrations *c* + with rate coefficients *k_fwd* and equilbrium constants *k_eq* + as a :class:`pymbolic.primitives.Expression` + """ + indices_reac = [sol.species_index(sp) for sp in sol.reaction(reaction_index).reactants] + indices_prod = [sol.species_index(sp) for sp in sol.reaction(reaction_index).products] + + if sol.reaction(reaction_index).orders: + nu_reac = [sol.reaction(reaction_index).orders[sp] for sp in sol.reaction(reaction_index).orders] + else: + nu_reac = [sol.reaction(reaction_index).reactants[sp] for sp in sol.reaction(reaction_index).reactants] + + r_fwd = np.prod([c[index] ** nu for index, nu in zip(indices_reac, nu_reac)]) + + if sol.reaction(reaction_index).reversible: + nu_prod = [sol.reaction(reaction_index).products[sp] for sp in sol.reaction(reaction_index).products] + r_rev = np.prod([c[index] ** nu for index, nu in zip(indices_prod, nu_prod)]) + return k_fwd[reaction_index] * (r_fwd - p.Variable("exp")(log_k_eq[reaction_index]) * r_rev) + else: + return k_fwd[reaction_index] * r_fwd + + +# }}} + + +# {{{ Species production rates + + +def production_rate_expr(sol: ct.Solution, species, r_net): + """ + :returns: Species production rate for species *species* in terms of + the net reaction rates of progress *r_net* as a + :class:`pymbolic.primitives.Expression` + """ + ones = _zeros_like(r_net[0]) + 1.0 + indices_fwd = [i for i, react in enumerate(sol.reactions()) if species in react.reactants] + indices_rev = [i for i, react in enumerate(sol.reactions()) if species in react.products] + nu_fwd = [sol.reactant_stoich_coeff(sol.species_index(species), react_index) for react_index in indices_fwd] + nu_rev = [sol.product_stoich_coeff(sol.species_index(species), prod_index) for prod_index in indices_rev] + sum_fwd = sum(nu * r_net[index] for nu, index in zip(nu_fwd, indices_fwd)) + sum_rev = sum(nu * r_net[index] for nu, index in zip(nu_rev, indices_rev)) + return (sum_rev - sum_fwd) * ones + + +def fwd_rate_of_progress_expr(sol: ct.Solution, reaction_index, c, k_fwd): + """ + :returns: Forward rate of progress (non-negative) for reaction + *reaction_index*, i.e. the forward part of + :func:`rate_of_progress_expr`. + """ + indices_reac = [sol.species_index(sp) for sp in sol.reaction(reaction_index).reactants] + if sol.reaction(reaction_index).orders: + nu_reac = [sol.reaction(reaction_index).orders[sp] for sp in sol.reaction(reaction_index).orders] + else: + nu_reac = [sol.reaction(reaction_index).reactants[sp] for sp in sol.reaction(reaction_index).reactants] + r_fwd = np.prod([c[index] ** nu for index, nu in zip(indices_reac, nu_reac)]) + return k_fwd[reaction_index] * r_fwd + + +def rev_rate_of_progress_expr(sol: ct.Solution, reaction_index, c, k_fwd, log_k_eq): + """ + :returns: Reverse rate of progress (non-negative) for reaction + *reaction_index*; zero for irreversible reactions. + """ + if not sol.reaction(reaction_index).reversible: + return _zeros_like(c[0]) + indices_prod = [sol.species_index(sp) for sp in sol.reaction(reaction_index).products] + nu_prod = [sol.reaction(reaction_index).products[sp] for sp in sol.reaction(reaction_index).products] + r_rev = np.prod([c[index] ** nu for index, nu in zip(indices_prod, nu_prod)]) + return k_fwd[reaction_index] * p.Variable("exp")(log_k_eq[reaction_index]) * r_rev + + +def _species_reaction_stoich(sol: ct.Solution, species): + """ + :returns: ``(idx_reactant, idx_product, nu_reactant, nu_product)`` -- the + reaction indices in which *species* appears as a reactant and as a + product, with the matching stoichiometric coefficients. Shared by the + creation/destruction splits, which differ only in how these are summed. + """ + si = sol.species_index(species) + idx_reactant = [i for i, react in enumerate(sol.reactions()) if species in react.reactants] + idx_product = [i for i, react in enumerate(sol.reactions()) if species in react.products] + nu_reactant = [sol.reactant_stoich_coeff(si, i) for i in idx_reactant] + nu_product = [sol.product_stoich_coeff(si, i) for i in idx_product] + return idx_reactant, idx_product, nu_reactant, nu_product + + +def creation_rate_expr(sol: ct.Solution, species, r_fwd, r_rev): + """ + :returns: Species creation rate for *species*: created as a product by + forward reactions and as a reactant by reverse reactions. Mirrors + Cantera's ``creation_rates``. + """ + ones = _zeros_like(r_fwd[0]) + 1.0 + idx_reactant, idx_product, nu_reactant, nu_product = _species_reaction_stoich(sol, species) + made = sum(nu * r_fwd[i] for nu, i in zip(nu_product, idx_product)) + sum(nu * r_rev[i] for nu, i in zip(nu_reactant, idx_reactant)) + return made * ones + + +def destruction_rate_expr(sol: ct.Solution, species, r_fwd, r_rev): + """ + :returns: Species destruction rate for *species*: consumed as a reactant by + forward reactions and as a product by reverse reactions. Mirrors + Cantera's ``destruction_rates``. creation - destruction == net. + """ + ones = _zeros_like(r_fwd[0]) + 1.0 + idx_reactant, idx_product, nu_reactant, nu_product = _species_reaction_stoich(sol, species) + lost = sum(nu * r_fwd[i] for nu, i in zip(nu_reactant, idx_reactant)) + sum(nu * r_rev[i] for nu, i in zip(nu_product, idx_product)) + return lost * ones + + +# }}} + +# vim diff --git a/toolchain/mfc/thermochem/fingerprint.py b/toolchain/mfc/thermochem/fingerprint.py new file mode 100644 index 000000000..bda8beea3 --- /dev/null +++ b/toolchain/mfc/thermochem/fingerprint.py @@ -0,0 +1,28 @@ +"""Content identities for mechanisms and the MFC-owned generator.""" + +import hashlib +import json +from functools import lru_cache +from pathlib import Path + + +def mechanism_fingerprint(solution): + """Hash the phase, species and reactions independently of the current state.""" + import cantera as ct + + data = dict(solution.input_data) + data.pop("state", None) + data["species"] = [species.input_data for species in solution.species()] + data["reactions"] = [reaction.input_data for reaction in solution.reactions()] + data["cantera-version"] = ct.__version__ + return hashlib.sha256(json.dumps(data, sort_keys=True).encode()).hexdigest() + + +@lru_cache(maxsize=1) +def generator_fingerprint(): + """Keep builds from different generator revisions in separate staging trees.""" + digest = hashlib.sha256() + root = Path(__file__).parent + for name in ("__init__.py", "fortran.py", "expressions.py", "module.f90.mako"): + digest.update((root / name).read_bytes()) + return digest.hexdigest() diff --git a/toolchain/mfc/thermochem/fortran.py b/toolchain/mfc/thermochem/fortran.py new file mode 100644 index 000000000..a88a56b38 --- /dev/null +++ b/toolchain/mfc/thermochem/fortran.py @@ -0,0 +1,247 @@ +"""Generate the mechanism-specific Fortran module used by MFC. See LICENSE. + +Adapted from the Pyrometheus 1.1.1 Fortran emitter. +""" + +import shlex +from functools import partial +from numbers import Integral +from pathlib import Path + +import cantera as ct +import pymbolic.primitives as p +from mako.template import Template +from pymbolic.mapper.stringifier import PREC_CALL, PREC_NONE, PREC_PRODUCT, StringifyMapper + +from . import expressions + +# {{{ code generation helpers + + +def pad_fortran(line, width): + line += " " * (width - 1 - len(line)) + line += "&" + return line + + +def wrap_line_base(line, level=0, width=80, indentation=" ", pad_func=lambda string, amount: string, lex_func=None): + """ + The input is a line of code at the given indentation level. Return the list + of lines that results from wrapping the line to the given width. Lines + subsequent to the first line in the returned list are padded with extra + indentation. The initial indentation level is not included in the input or + output lines. + + The `pad_func` argument is a function that adds line continuations. The + `lex_func` argument returns the list of tokens in the line. + """ + if lex_func is None: + lex_func = partial(shlex.split, posix=False) + + tokens = lex_func(line) + resulting_lines = [] + at_line_start = True + indentation_len = len(level * indentation) + current_line = "" + padding_width = width - indentation_len + for index, word in enumerate(tokens): + has_next_word = index < len(tokens) - 1 + word_len = len(word) + if not at_line_start: + next_len = indentation_len + len(current_line) + 1 + word_len + if next_len < width or (not has_next_word and next_len == width): + # The word goes on the same line. + current_line += " " + word + else: + # The word goes on the next line. + resulting_lines.append(pad_func(current_line, padding_width)) + at_line_start = True + current_line = indentation + if at_line_start: + current_line += word + at_line_start = False + resulting_lines.append(current_line) + return resulting_lines + + +def count_leading_spaces(s): + n = 0 + while n < len(s) and s[n] == " ": + n += 1 + return n + + +def wrap_code(s, indent=4): + lines = s.split("\n") + result_lines = [] + for ln in lines: + nspaces = count_leading_spaces(ln) + level, remainder = divmod(nspaces, indent) + + if remainder != 0: + raise ValueError(f"indentation of '{ln}' is not a multiple of " f"{indent}") + + result_lines.extend((level * indent) * " " + subln for subln in wrap_line_base(ln, level=level, indentation=" " * indent, pad_func=pad_fortran)) + + return "\n".join(result_lines) + + +def float_to_fortran(num, kind): + result = f"{num}" + if "." not in result and "e" not in result.lower(): + result += ".0" + result += f"_{kind}" + if num < 0: + result = "(%s)" % result + return result + + +def str_np(ary, kind): + return ", ".join(float_to_fortran(entry, kind) for entry in ary) + + +# }}} + + +def validate_mechanism(sol): + """Reject unsupported physics before emitting a partially valid module.""" + if sol.thermo_model != "ideal-gas": + raise ValueError(f"MFC thermochemistry requires ideal-gas thermodynamics, got {sol.thermo_model!r}") + if sol.transport_model not in ("mixture-averaged", "multicomponent", "unity-Lewis-number"): + raise ValueError(f"MFC thermochemistry requires gas transport data, got {sol.transport_model!r}") + for species in sol.species(): + if not isinstance(species.thermo, ct.NasaPoly2): + raise ValueError(f"Species {species.name}: MFC thermochemistry supports NASA7 polynomials only") + supported = {"Arrhenius", "three-body-Arrhenius", "falloff-Troe", "falloff-Lindemann"} + for i, reaction in enumerate(sol.reactions()): + label = f"Reaction {i + 1} ({reaction.equation})" + if reaction.reaction_type not in supported: + raise ValueError(f"{label}: unsupported rate type {reaction.reaction_type!r}") + if reaction.orders: + raise ValueError(f"{label}: custom reaction orders are not supported") + rates = (reaction.rate.low_rate, reaction.rate.high_rate) if reaction.reaction_type.startswith("falloff") else (reaction.rate,) + if any(rate.pre_exponential_factor <= 0 for rate in rates): + raise ValueError(f"{label}: Arrhenius pre-exponential factors must be positive") + + +def generate_fortran(solution, module_name="m_thermochem", scalar_type="real(dp)", offload=None): + """Emit MFC's thermodynamic, kinetics and transport interface from Cantera.""" + import re + + if not re.fullmatch(r"[A-Za-z][A-Za-z0-9_]{0,62}", module_name): + raise ValueError(f"Invalid Fortran module name: {module_name!r}") + if scalar_type not in ("real(sp)", "real(dp)"): + raise ValueError(f"Unsupported scalar type: {scalar_type!r}") + directives = {None: "! name", "acc": "!$acc routine seq", "mp": "!$omp declare target"} + if offload not in directives: + raise ValueError(f"Unsupported offload mode: {offload!r}") + validate_mechanism(solution) + kind = "sp" if scalar_type == "real(sp)" else "dp" + falloff = [(i, r) for i, r in enumerate(solution.reactions()) if r.reaction_type.startswith("falloff")] + three_body = [(i, r) for i, r in enumerate(solution.reactions()) if r.reaction_type == "three-body-Arrhenius"] + template = Template(filename=str(Path(__file__).with_name("module.f90.mako"))) + return wrap_code( + template.render( + ct=ct, + sol=solution, + str_np=partial(str_np, kind=kind), + cgm=FortranExpressionMapper(kind), + Variable=p.Variable, + float_to_fortran=partial(float_to_fortran, kind=kind), + real_type=scalar_type, + kind=kind, + species_name_length=max(map(len, solution.species_names)), + gpu_routine=f"#define GPU_ROUTINE(name) {directives[offload]}", + module_name=module_name, + ce=expressions, + falloff_reactions=falloff, + falloff_indices={i for i, _ in falloff}, + three_body_reactions=three_body, + ) + ) + + +# {{{ fortran expression generation + + +class FortranExpressionMapper(StringifyMapper): + """Converts expressions to Fortran code.""" + + def __init__(self, kind): + super().__init__() + self.kind = kind + + def map_constant(self, expr, enclosing_prec): + if isinstance(expr, bool): + if expr: + return ".true." + else: + return ".false." + else: + return float_to_fortran(expr, self.kind) + + def map_variable(self, expr, enclosing_prec): + return expr.name + + def map_lookup(self, expr, enclosing_prec): + return self.parenthesize_if_needed(self.format("%s%%%s", self.rec(expr.aggregate, PREC_CALL), expr.name), enclosing_prec, PREC_CALL) + + def map_subscript(self, expr, enclosing_prec): + def get_base_and_indices(expr): + if not hasattr(expr, "aggregate") or not hasattr(expr, "index"): + return expr, [] + + # Get current level indices + if isinstance(expr.index, tuple): + current_indices = list(expr.index) + else: + current_indices = [expr.index] + + # Only recurse if aggregate is another subscript + if hasattr(expr.aggregate, "aggregate") and hasattr(expr.aggregate, "index"): + base, prev_indices = get_base_and_indices(expr.aggregate) + return base, prev_indices + current_indices + else: + return expr.aggregate, current_indices + + # Get base array and all indices + base_array, all_indices = get_base_and_indices(expr) + + # Convert zero-based expression indices without applying real-literal kinds. + def convert_index(idx): + if isinstance(idx, Integral): + return str(idx + 1) + return f"({self.rec(idx, PREC_NONE)} + 1)" + + # Format indices, converting floats to integers and adding 1 + index_str = ", ".join(convert_index(idx) for idx in all_indices) + + # Format the final expression + return self.parenthesize_if_needed(self.format("%s(%s)", self.rec(base_array, PREC_CALL), index_str), enclosing_prec, PREC_CALL) + + def map_product(self, expr, enclosing_prec, *args, **kwargs): + # This differs from the superclass only by adding spaces + # around the operator, which provide an opportunity for + # line breaking. + return self.parenthesize_if_needed(self.join_rec(" * ", expr.children, PREC_PRODUCT, *args, **kwargs), enclosing_prec, PREC_PRODUCT) + + def map_logical_not(self, expr, enclosing_prec): + from pymbolic.mapper.stringifier import PREC_UNARY + + return self.parenthesize_if_needed(".not. " + self.rec(expr.child, PREC_UNARY), enclosing_prec, PREC_UNARY) + + def map_logical_or(self, expr, enclosing_prec): + from pymbolic.mapper.stringifier import PREC_LOGICAL_OR + + return self.parenthesize_if_needed(self.join_rec(" .or. ", expr.children, PREC_LOGICAL_OR), enclosing_prec, PREC_LOGICAL_OR) + + def map_logical_and(self, expr, enclosing_prec): + from pymbolic.mapper.stringifier import PREC_LOGICAL_AND + + return self.parenthesize_if_needed(self.join_rec(" .and. ", expr.children, PREC_LOGICAL_AND), enclosing_prec, PREC_LOGICAL_AND) + + def map_if(self, expr, enclosing_prec): + return self.format("merge(%s)" % self.join_rec(", ", [expr.then, expr.else_, expr.condition], PREC_NONE)) + + +# }}} diff --git a/toolchain/mfc/thermochem/module.f90.mako b/toolchain/mfc/thermochem/module.f90.mako new file mode 100644 index 000000000..0a364e4f6 --- /dev/null +++ b/toolchain/mfc/thermochem/module.f90.mako @@ -0,0 +1,776 @@ +! MFC-owned thermochemistry; derived from Pyrometheus 1.1.1 (MIT). + + +${gpu_routine} + +module ${module_name} + + implicit none + + integer, parameter :: sp = selected_real_kind(6,37) ! Single precision + integer, parameter :: dp = selected_real_kind(15,307) ! Double precision + + integer, parameter :: num_elements = ${sol.n_elements} + integer, parameter :: num_species = ${sol.n_species} + integer, parameter :: num_reactions = ${sol.n_reactions} + integer, parameter :: num_falloff = ${len(falloff_reactions)} + ${real_type}, parameter :: one_atm = ${float_to_fortran(ct.one_atm)} + ${real_type}, parameter :: gas_constant = ${float_to_fortran(ct.gas_constant)} + ${real_type}, parameter :: molecular_weights(${sol.n_species}) = & + (/ ${str_np(sol.molecular_weights)} /) + ${real_type}, parameter :: inv_molecular_weights(${sol.n_species}) = & + (/ ${str_np(1/sol.molecular_weights)} /) + + character(len=${species_name_length}), parameter :: species_names(${sol.n_species}) = & + (/ ${", ".join('"'+s.ljust(species_name_length)+'"' for s in sol.species_names)} /) + + character(len=4), parameter :: element_names(${sol.n_elements}) = & + (/ ${", ".join('"'+'{0: <4}'.format(e)+'"' for e in sol.element_names)} /) + +contains + + subroutine get_species_name(sp_index, sp_name) + + integer, intent(in) :: sp_index + character(len=*), intent(out) :: sp_name + + sp_name = species_names(sp_index) + + end subroutine get_species_name + + subroutine get_species_index(sp_name, sp_index) + + character(len=*), intent(in) :: sp_name + integer, intent(out) :: sp_index + + integer :: idx + + sp_index = 0 + loop:do idx = 1, num_species + if(trim(adjustl(sp_name)) .eq. trim(species_names(idx))) then + sp_index = idx + exit loop + end if + end do loop + + end subroutine get_species_index + + subroutine get_element_index(el_name, el_index) + + character(len=*), intent(in) :: el_name + integer, intent(out) :: el_index + + integer :: idx + + el_index = 0 + loop:do idx = 1, num_elements + if(trim(adjustl(el_name)) .eq. trim(element_names(idx))) then + el_index = idx + exit loop + end if + end do loop + + end subroutine get_element_index + + subroutine get_specific_gas_constant(mass_fractions, specific_gas_constant) + + GPU_ROUTINE(get_specific_gas_constant) + + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out) :: specific_gas_constant + + specific_gas_constant = gas_constant * ( & + %for i in range(sol.n_species): + + inv_molecular_weights(${i+1})*mass_fractions(${i+1}) & + %endfor + ) + + end subroutine get_specific_gas_constant + + subroutine get_density(pressure, temperature, mass_fractions, density) + + GPU_ROUTINE(get_density) + + ${real_type}, intent(in) :: pressure + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out) :: density + + ${real_type} :: mix_mol_weight + + call get_mixture_molecular_weight(mass_fractions, mix_mol_weight) + density = pressure * mix_mol_weight / (gas_constant * temperature) + + end subroutine get_density + + subroutine get_pressure(density, temperature, mass_fractions, pressure) + + GPU_ROUTINE(get_pressure) + + ${real_type}, intent(in) :: density + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out) :: pressure + + ${real_type} :: mix_mol_weight + + call get_mixture_molecular_weight(mass_fractions, mix_mol_weight) + pressure = density * gas_constant * temperature / mix_mol_weight + + end subroutine get_pressure + + subroutine get_mixture_molecular_weight(mass_fractions, mix_mol_weight) + + GPU_ROUTINE(get_mixture_molecular_weight) + + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out) :: mix_mol_weight + + mix_mol_weight = 1.0e0_${kind} / ( & + %for i in range(sol.n_species): + + inv_molecular_weights(${i+1})*mass_fractions(${i+1}) & + %endfor + ) + + end subroutine get_mixture_molecular_weight + + subroutine get_concentrations(density, mass_fractions, concentrations) + + GPU_ROUTINE(get_concentrations) + + ${real_type}, intent(in) :: density + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out), dimension(${sol.n_species}) :: concentrations + + %for i in range(sol.n_species): + concentrations(${i+1}) = density * & + inv_molecular_weights(${i+1}) * mass_fractions(${i+1}) + %endfor + + end subroutine get_concentrations + + subroutine get_mole_fractions(mix_mol_weight, mass_fractions, mole_fractions) + + GPU_ROUTINE(get_mole_fractions) + + ${real_type}, intent(in) :: mix_mol_weight + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out), dimension(${sol.n_species}) :: mole_fractions + + %for i in range(sol.n_species): + mole_fractions(${i+1}) = inv_molecular_weights(${i+1}) * & + mass_fractions(${i+1}) * mix_mol_weight + %endfor + + end subroutine get_mole_fractions + + subroutine get_mass_averaged_property(& + & mass_fractions, spec_property, mix_property) + + GPU_ROUTINE(get_mass_averaged_property) + + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(in), dimension(${sol.n_species}) :: spec_property + ${real_type}, intent(out) :: mix_property + + mix_property = ( & + %for i in range(sol.n_species): + + inv_molecular_weights(${i+1})*mass_fractions(${i+1}) & + *spec_property(${i+1}) & + %endfor + ) + + end subroutine get_mass_averaged_property + + subroutine get_mixture_specific_heat_cp_mass(temperature, mass_fractions, cp_mix) + + GPU_ROUTINE(get_mixture_specific_heat_cp_mass) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out) :: cp_mix + + ${real_type}, dimension(${sol.n_species}) :: cp0_r + + call get_species_specific_heats_r(temperature, cp0_r) + call get_mass_averaged_property(mass_fractions, cp0_r, cp_mix) + cp_mix = cp_mix * gas_constant + + end subroutine get_mixture_specific_heat_cp_mass + + subroutine get_mixture_specific_heat_cv_mass(temperature, mass_fractions, cv_mix) + + GPU_ROUTINE(get_mixture_specific_heat_cv_mass) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out) :: cv_mix + + ${real_type}, dimension(${sol.n_species}) :: cp0_r + + call get_species_specific_heats_r(temperature, cp0_r) + + %for i in range(sol.n_species): + cp0_r(${i+1}) = cp0_r(${i+1}) - 1.e0_${kind} + %endfor + + call get_mass_averaged_property(mass_fractions, cp0_r, cv_mix) + cv_mix = cv_mix * gas_constant + + end subroutine get_mixture_specific_heat_cv_mass + + subroutine get_mixture_enthalpy_mass(temperature, mass_fractions, h_mix) + + GPU_ROUTINE(get_mixture_enthalpy_mass) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out) :: h_mix + + ${real_type}, dimension(${sol.n_species}) :: h0_rt + + call get_species_enthalpies_rt(temperature, h0_rt) + call get_mass_averaged_property(mass_fractions, h0_rt, h_mix) + h_mix = h_mix * gas_constant * temperature + + end subroutine get_mixture_enthalpy_mass + + subroutine get_mixture_energy_mass(temperature, mass_fractions, e_mix) + + GPU_ROUTINE(get_mixture_energy_mass) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out) :: e_mix + + ${real_type}, dimension(${sol.n_species}) :: h0_rt + + call get_species_enthalpies_rt(temperature, h0_rt) + + %for i in range(sol.n_species): + h0_rt(${i+1}) = h0_rt(${i+1}) - 1.e0_${kind} + %endfor + + call get_mass_averaged_property(mass_fractions, h0_rt, e_mix) + e_mix = e_mix * gas_constant * temperature + + end subroutine get_mixture_energy_mass + + subroutine get_species_specific_heats_r(temperature, cp0_r) + + GPU_ROUTINE(get_species_specific_heats_r) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(out), dimension(${sol.n_species}) :: cp0_r + + %for i, sp in enumerate(sol.species()): + cp0_r(${i+1}) = ${cgm(ce.poly_to_expr(sp.thermo, "temperature"))} + %endfor + + end subroutine get_species_specific_heats_r + + subroutine get_species_enthalpies_rt(temperature, h0_rt) + + GPU_ROUTINE(get_species_enthalpies_rt) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(out), dimension(${sol.n_species}) :: h0_rt + + %for i, sp in enumerate(sol.species()): + h0_rt(${i+1}) = ${cgm(ce.poly_to_enthalpy_expr(sp.thermo, "temperature"))} + %endfor + + end subroutine get_species_enthalpies_rt + + subroutine get_species_entropies_r(temperature, s0_r) + + GPU_ROUTINE(get_species_entropies_r) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(out), dimension(${sol.n_species}) :: s0_r + + %for i, sp in enumerate(sol.species()): + s0_r(${i+1}) = ${cgm(ce.poly_to_entropy_expr(sp.thermo, "temperature"))} + %endfor + + end subroutine get_species_entropies_r + + subroutine get_species_gibbs_rt(temperature, g0_rt) + + GPU_ROUTINE(get_species_gibbs_rt) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(out), dimension(${sol.n_species}) :: g0_rt + + ${real_type}, dimension(${sol.n_species}) :: h0_rt + ${real_type}, dimension(${sol.n_species}) :: s0_r + + call get_species_enthalpies_rt(temperature, h0_rt) + call get_species_entropies_r(temperature, s0_r) + + %for i in range(sol.n_species): + g0_rt(${i+1}) = h0_rt(${i+1}) - s0_r(${i+1}) + %endfor + + end subroutine get_species_gibbs_rt + + subroutine get_equilibrium_constants(temperature, k_eq) + + GPU_ROUTINE(get_equilibrium_constants) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(out), dimension(${sol.n_reactions}) :: k_eq + + ${real_type} :: rt + ${real_type} :: c0 + + ${real_type}, dimension(${sol.n_species}) :: g0_rt + + rt = gas_constant * temperature + c0 = log(one_atm/rt) + + call get_species_gibbs_rt(temperature, g0_rt) + + %for i, react in enumerate(sol.reactions()): + %if react.reversible: + k_eq(${i+1}) = ${cgm( + ce.equilibrium_constants_expr(sol, i, Variable("g0_rt")))} + %else: + k_eq(${i+1}) = -0.1e0_${kind}*temperature + %endif + %endfor + + end subroutine get_equilibrium_constants + + subroutine get_temperature( & + & enthalpy_or_energy, t_guess, mass_fractions, do_energy, temperature) + + GPU_ROUTINE(get_temperature) + + logical, intent(in) :: do_energy + ${real_type}, intent(in) :: enthalpy_or_energy + ${real_type}, intent(in) :: t_guess + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out) :: temperature + + integer :: iter + integer, parameter :: num_iter = 500 + ${real_type}, parameter :: tol = 1.0e-06_${kind} + + ${real_type} :: iter_temp + ${real_type} :: iter_energy + ${real_type} :: iter_energy_deriv + ${real_type} :: iter_rhs + ${real_type} :: iter_deriv + + iter_rhs = 0.e0_${kind} + iter_deriv = 1.e0_${kind} + iter_temp = t_guess + + do iter = 1, num_iter + if(do_energy) then + call get_mixture_specific_heat_cv_mass(& + & iter_temp, mass_fractions, iter_energy_deriv) + call get_mixture_energy_mass(iter_temp, mass_fractions, iter_energy) + else + call get_mixture_specific_heat_cp_mass(& + & iter_temp, mass_fractions, iter_energy_deriv) + call get_mixture_enthalpy_mass(& + & iter_temp, mass_fractions, iter_energy) + endif + iter_rhs = enthalpy_or_energy - iter_energy + iter_deriv = (-1.e0_${kind})*iter_energy_deriv + iter_temp = iter_temp - iter_rhs / iter_deriv + if(abs(iter_rhs/iter_deriv) .lt. tol) exit + end do + + temperature = iter_temp + + end subroutine get_temperature + + %if falloff_reactions: + subroutine get_falloff_rates(temperature, concentrations, k_fwd) + + GPU_ROUTINE(get_falloff_rates) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: concentrations + ${real_type}, intent(out), dimension(${sol.n_reactions}) :: k_fwd + + ${real_type}, dimension(${len(falloff_reactions)}) :: k_high + ${real_type}, dimension(${len(falloff_reactions)}) :: k_low + ${real_type}, dimension(${len(falloff_reactions)}) :: reduced_pressure + ${real_type}, dimension(${len(falloff_reactions)}) :: falloff_center + ${real_type}, dimension(${len(falloff_reactions)}) :: falloff_factor + ${real_type}, dimension(${len(falloff_reactions)}) :: falloff_function + + %for i, (_, react) in enumerate(falloff_reactions): + k_high(${i+1}) = ${cgm(ce.rate_coefficient_expr( + react.rate.high_rate, + Variable("temperature")))} + %endfor + + %for i, (_, react) in enumerate(falloff_reactions): + k_low(${i+1}) = ${cgm(ce.rate_coefficient_expr( + react.rate.low_rate, + Variable("temperature")))} + %endfor + + %for i, (_, react) in enumerate(falloff_reactions): + reduced_pressure(${i+1}) = (${cgm( + ce.third_body_efficiencies_expr(sol, + react, + Variable("concentrations")))})*k_low(${i+1})/k_high(${i+1}) + %endfor + + %for i, (_, react) in enumerate(falloff_reactions): + falloff_center(${i+1}) = ${cgm(ce.troe_falloff_center_expr( + react, Variable("temperature")))} + %endfor + + %for i, (_, react) in enumerate(falloff_reactions): + falloff_factor(${i+1}) = ${cgm(ce.troe_falloff_factor_expr(react, i, + Variable("reduced_pressure"), Variable("falloff_center")))} + %endfor + + %for i, (_, react) in enumerate(falloff_reactions): + falloff_function(${i+1}) = ${cgm(ce.falloff_function_expr( + react, i, + Variable("falloff_factor"), + Variable("falloff_center")))} + %endfor + + %for i, (j, react) in enumerate(falloff_reactions): + k_fwd(${j+1}) = k_high(${i+1})*falloff_function(${i+1}) * & + reduced_pressure(${i+1})/(1.e0_${kind} + reduced_pressure(${i+1})) + %endfor + + end subroutine get_falloff_rates + + %endif + subroutine get_fwd_rate_coefficients(temperature, concentrations, k_fwd) + + GPU_ROUTINE(get_fwd_rate_coefficients) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: concentrations + ${real_type}, intent(out), dimension(${sol.n_reactions}) :: k_fwd + + %if falloff_reactions: + ${real_type}, dimension(${len(falloff_reactions)}) :: k_falloff + %endif + + %for i, react in enumerate(sol.reactions()): + %if i in falloff_indices: + k_fwd(${i+1}) = 0.e0_${kind} + %else: + k_fwd(${i+1}) = ${cgm(ce.rate_coefficient_expr(react.rate, + Variable("temperature")))} + %endif + %endfor + + %for j, react in three_body_reactions: + k_fwd(${j+1}) = k_fwd(${j+1}) * ( & + ${cgm(ce.third_body_efficiencies_expr( + sol, react, Variable("concentrations")))}) + %endfor + + %if falloff_reactions: + call get_falloff_rates(temperature, concentrations, k_fwd) + %endif + + end subroutine get_fwd_rate_coefficients + + subroutine get_net_rates_of_progress(temperature, concentrations, r_net) + + GPU_ROUTINE(get_net_rates_of_progress) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: concentrations + ${real_type}, intent(out), dimension(${sol.n_reactions}) :: r_net + + ${real_type}, dimension(${sol.n_reactions}) :: k_fwd + ${real_type}, dimension(${sol.n_reactions}) :: log_k_eq + + call get_fwd_rate_coefficients(temperature, concentrations, k_fwd) + call get_equilibrium_constants(temperature, log_k_eq) + %for i in range(sol.n_reactions): + r_net(${i+1}) = ${cgm(ce.rate_of_progress_expr(sol, i, + Variable("concentrations"), + Variable("k_fwd"), Variable("log_k_eq")))} + %endfor + + end subroutine get_net_rates_of_progress + + subroutine get_net_production_rates(density, temperature, mass_fractions, omega) + + GPU_ROUTINE(get_net_production_rates) + + ${real_type}, intent(in) :: density + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out), dimension(${sol.n_species}) :: omega + + ${real_type}, dimension(${sol.n_species}) :: concentrations + ${real_type}, dimension(${sol.n_reactions}) :: r_net + + call get_concentrations(density, mass_fractions, concentrations) + call get_net_rates_of_progress(temperature, concentrations, r_net) + + %for i, sp in enumerate(sol.species()): + omega(${i+1}) = ${cgm(ce.production_rate_expr(sol, + sp.name, Variable("r_net")))} + %endfor + + end subroutine get_net_production_rates + + subroutine get_fwd_rates_of_progress(temperature, concentrations, r_fwd) + + GPU_ROUTINE(get_fwd_rates_of_progress) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: concentrations + ${real_type}, intent(out), dimension(${sol.n_reactions}) :: r_fwd + + ${real_type}, dimension(${sol.n_reactions}) :: k_fwd + + call get_fwd_rate_coefficients(temperature, concentrations, k_fwd) + %for i in range(sol.n_reactions): + r_fwd(${i+1}) = ${cgm(ce.fwd_rate_of_progress_expr(sol, i, + Variable("concentrations"), Variable("k_fwd")))} + %endfor + + end subroutine get_fwd_rates_of_progress + + subroutine get_rev_rates_of_progress(temperature, concentrations, r_rev) + + GPU_ROUTINE(get_rev_rates_of_progress) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: concentrations + ${real_type}, intent(out), dimension(${sol.n_reactions}) :: r_rev + + ${real_type}, dimension(${sol.n_reactions}) :: k_fwd + ${real_type}, dimension(${sol.n_reactions}) :: log_k_eq + + call get_fwd_rate_coefficients(temperature, concentrations, k_fwd) + call get_equilibrium_constants(temperature, log_k_eq) + %for i in range(sol.n_reactions): + r_rev(${i+1}) = ${cgm(ce.rev_rate_of_progress_expr(sol, i, + Variable("concentrations"), + Variable("k_fwd"), Variable("log_k_eq")))} + %endfor + + end subroutine get_rev_rates_of_progress + + subroutine get_creation_rates(density, temperature, mass_fractions, cdot) + + GPU_ROUTINE(get_creation_rates) + + ${real_type}, intent(in) :: density + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out), dimension(${sol.n_species}) :: cdot + + ${real_type}, dimension(${sol.n_species}) :: concentrations + ${real_type}, dimension(${sol.n_reactions}) :: r_fwd, r_rev + + call get_concentrations(density, mass_fractions, concentrations) + call get_fwd_rates_of_progress(temperature, concentrations, r_fwd) + call get_rev_rates_of_progress(temperature, concentrations, r_rev) + + %for i, sp in enumerate(sol.species()): + cdot(${i+1}) = ${cgm(ce.creation_rate_expr(sol, sp.name, + Variable("r_fwd"), Variable("r_rev")))} + %endfor + + end subroutine get_creation_rates + + subroutine get_destruction_rates(density, temperature, mass_fractions, ddot) + + GPU_ROUTINE(get_destruction_rates) + + ${real_type}, intent(in) :: density + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out), dimension(${sol.n_species}) :: ddot + + ${real_type}, dimension(${sol.n_species}) :: concentrations + ${real_type}, dimension(${sol.n_reactions}) :: r_fwd, r_rev + + call get_concentrations(density, mass_fractions, concentrations) + call get_fwd_rates_of_progress(temperature, concentrations, r_fwd) + call get_rev_rates_of_progress(temperature, concentrations, r_rev) + + %for i, sp in enumerate(sol.species()): + ddot(${i+1}) = ${cgm(ce.destruction_rate_expr(sol, sp.name, + Variable("r_fwd"), Variable("r_rev")))} + %endfor + + end subroutine get_destruction_rates + + subroutine get_creation_destruction_rates(density, temperature, & + mass_fractions, cdot, ddot) + + GPU_ROUTINE(get_creation_destruction_rates) + + ${real_type}, intent(in) :: density + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out), dimension(${sol.n_species}) :: cdot, ddot + + ${real_type}, dimension(${sol.n_species}) :: concentrations + ${real_type}, dimension(${sol.n_reactions}) :: r_fwd, r_rev + + call get_concentrations(density, mass_fractions, concentrations) + call get_fwd_rates_of_progress(temperature, concentrations, r_fwd) + call get_rev_rates_of_progress(temperature, concentrations, r_rev) + + %for i, sp in enumerate(sol.species()): + cdot(${i+1}) = ${cgm(ce.creation_rate_expr(sol, sp.name, + Variable("r_fwd"), Variable("r_rev")))} + ddot(${i+1}) = ${cgm(ce.destruction_rate_expr(sol, sp.name, + Variable("r_fwd"), Variable("r_rev")))} + %endfor + + end subroutine get_creation_destruction_rates + + subroutine get_species_viscosities(temperature, viscosities) + + GPU_ROUTINE(get_species_viscosities) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(out), dimension(${sol.n_species}) :: viscosities + + %for sp in range(sol.n_species): + viscosities(${sp+1}) = ${cgm(ce.viscosity_polynomial_expr( + sol.get_viscosity_polynomial(sp), + Variable("temperature")))} + %endfor + + end subroutine get_species_viscosities + + subroutine get_species_thermal_conductivities(temperature, conductivities) + + GPU_ROUTINE(get_species_thermal_conductivities) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(out), dimension(${sol.n_species}) :: conductivities + + %for sp in range(sol.n_species): + conductivities(${sp+1}) = ${cgm(ce.conductivity_polynomial_expr( + sol.get_thermal_conductivity_polynomial(sp), + Variable("temperature")))} + %endfor + + end subroutine get_species_thermal_conductivities + + subroutine get_species_binary_mass_diffusivities(temperature, diffusivities) + + GPU_ROUTINE(get_species_binary_mass_diffusivities) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(out), dimension(${sol.n_species}, ${sol.n_species})& + :: diffusivities + + %for i in range(sol.n_species): + %for j in range(sol.n_species): + diffusivities(${i + 1}, ${j + 1}) = ${cgm(ce.diffusivity_polynomial_expr( + sol.get_binary_diff_coeffs_polynomial(i, j), + Variable("temperature")))} + %endfor + %endfor + + end subroutine get_species_binary_mass_diffusivities + + subroutine get_mixture_viscosity_mixavg(& + temperature, mass_fractions, mixture_viscosity_mixavg) + + GPU_ROUTINE(get_mixture_viscosity_mixavg) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out) :: mixture_viscosity_mixavg + + ${real_type} :: mix_mol_weight + ${real_type}, dimension(${sol.n_species}) :: & + mole_fractions, viscosities, mix_rule_f + + call get_mixture_molecular_weight(mass_fractions, mix_mol_weight) + call get_mole_fractions(mix_mol_weight, mass_fractions, mole_fractions) + call get_species_viscosities(temperature, viscosities) + + %for sp in range(sol.n_species): + mix_rule_f(${sp + 1}) = ${cgm(ce.viscosity_mixture_rule_wilke_expr(sol, sp, + Variable("mole_fractions"), Variable("viscosities")))} + %endfor + + mixture_viscosity_mixavg = sum(mole_fractions*viscosities/mix_rule_f) + + end subroutine get_mixture_viscosity_mixavg + + subroutine get_mixture_thermal_conductivity_mixavg(temperature, & + mass_fractions, mixture_thermal_conductivity_mixavg) + + GPU_ROUTINE(get_mixture_thermal_conductivity_mixavg) + + ${real_type}, intent(in) :: temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out) :: mixture_thermal_conductivity_mixavg + + ${real_type} :: mix_mol_weight + ${real_type}, dimension(${sol.n_species}) :: mole_fractions, conductivities + + call get_mixture_molecular_weight(mass_fractions, mix_mol_weight) + call get_mole_fractions(mix_mol_weight, mass_fractions, mole_fractions) + call get_species_thermal_conductivities(temperature, conductivities) + + mixture_thermal_conductivity_mixavg = 0.5_${kind}*(& + sum(mole_fractions*conductivities) + & + 1/sum(mole_fractions/conductivities)) + + end subroutine get_mixture_thermal_conductivity_mixavg + + subroutine get_species_mass_diffusivities_mixavg(& + pressure, temperature, mass_fractions, mass_diffusivities_mixavg) + + GPU_ROUTINE(get_species_mass_diffusivities_mixavg) + + ${real_type}, intent(in) :: pressure, temperature + ${real_type}, intent(in), dimension(${sol.n_species}) :: mass_fractions + ${real_type}, intent(out), dimension(${sol.n_species}) :: & + mass_diffusivities_mixavg + + ${real_type} :: mix_mol_weight + ${real_type}, dimension(${sol.n_species}) :: mole_fractions, x_sum, denom + ${real_type}, dimension(${sol.n_species}, ${sol.n_species}) :: bdiff_ij + + call get_mixture_molecular_weight(mass_fractions, mix_mol_weight) + call get_mole_fractions(mix_mol_weight, mass_fractions, mole_fractions) + call get_species_binary_mass_diffusivities(temperature, bdiff_ij) + + %for sp in range(sol.n_species): + x_sum(${sp + 1}) = ${cgm(ce.diffusivity_mixture_rule_denom_expr( + sol, sp, Variable("mole_fractions"), Variable("bdiff_ij")))} + %endfor + + %for sp in range(sol.n_species): + denom(${sp + 1}) = x_sum(${sp + 1}) - & + mole_fractions(${sp + 1})/bdiff_ij(${sp + 1}, ${sp + 1}) + %endfor + + %for sp in range(sol.n_species): + if (denom(${sp + 1}) .gt. 0e0_${kind}) then + mass_diffusivities_mixavg(${sp + 1}) = & + (mix_mol_weight - & + mole_fractions(${sp + 1})*molecular_weights(${sp + 1}))& + /(pressure * mix_mol_weight * denom(${sp + 1})) + else + mass_diffusivities_mixavg(${sp + 1}) = & + bdiff_ij(${sp + 1}, ${sp + 1}) / pressure + end if + %endfor + + end subroutine get_species_mass_diffusivities_mixavg + +end module ${module_name} diff --git a/toolchain/pyproject.toml b/toolchain/pyproject.toml index 893b114b4..e013f474a 100644 --- a/toolchain/pyproject.toml +++ b/toolchain/pyproject.toml @@ -35,9 +35,7 @@ dependencies = [ # Chemistry "cantera>=3.1.0", - "pyrometheus == 1.1.1", # ==1.1.1 for get_creation_destruction_rates; also provides pyrometheus.flamelets (mixing-layer examples) - "jax", # flamelet IC solve in the reacting mixing-layer examples (via pyrometheus.flamelets) - "scipy", # flamelet IC interpolation (scipy.interpolate) in the reacting mixing-layer examples + "pymbolic>=2022.2", # Fortran generation; allow the last Python 3.9-compatible releases # Frontier Profiling "astunparse==1.6.2",