feat: add 3D reacting mixing-layer example (H2/N2-air, Mach 1.5) + hcid=371

.lychee.toml

@@ -33,4 +33,5 @@ exclude = [
     "https://code\\.visualstudio\\.com/?$",           # Root page returns 403 to automated requests
     "https://stackoverflow\\.com",                    # Returns 403 to automated requests
     "https://marketplace\\.visualstudio\\.com",        # Returns 503 to automated requests
+    "https://web\\.eng\\.ucsd\\.edu",                    # San Diego mechanism page has an untrusted SSL cert
 ]

docs/documentation/case.md

@@ -282,6 +282,13 @@ The code provides three pre-built patches for dimensional extrusion of initial c
   `bf_spatial_support` body force) so no single extrusion axis applies. The file's line
   count, origin, and (uniform) cell spacing must match the run grid; a mismatched file is
   rejected with a fatal error, so regenerate the IC whenever the grid changes.
+- `case(371)`: `case(370)` plus a closed-form spanwise (z) modulation, so the IC has
+  genuine 3D content from step 0. The cross-stream (mom%%beg+1) velocity read from the
+  file is scaled by `1 + 0.5*cos(k_z z)` and the spanwise (mom%%end) component is set
+  from that result; the streamwise component is left as read. `k_z = 2*pi/L_z` uses the
+  global z extent, so the IC does not depend on the MPI decomposition and is continuous
+  across a periodic `bc_z`. Assumes uniform z spacing. Used by
+  `examples/3D_reacting_mixing_layer`.
 
 Setup: Only requires specifying `files_dir` and filename pattern via `file_extension`. The files are located, for example, at `examples/1D_flamelet/IC`, and their format is `prim.XX.YY.file_extension.dat`.
 Implementation: All variables and file handling are managed in the `case.py` file of the simulation.

examples/3D_reacting_mixing_layer/README.md

@@ -0,0 +1,58 @@
+# 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
+`examples/2D_reacting_mixing_layer`, which runs the same flamelet machinery at `Mc = 0.3`.
+
+## Configuration
+
+| Parameter | Value |
+|---|---|
+| Oxidizer stream | air, `X_O2 = 0.21`, 500 K |
+| Fuel stream | `X_H2 = 0.5`, balance N2, 300 K |
+| Pressure | 101325 Pa |
+| Vorticity thickness `delta_omega` | 1.0e-3 m |
+| Convective Mach number `Mc` | 1.5 |
+| Domain | `[15, 20, 10] delta_omega` in `(x, y, z)` |
+| Grid | 14 pts/`delta_omega` in x and z, 28 in y, so 210 x 560 x 140 |
+| Time step | 1e-9 s, RK3 |
+| Numerics | WENO5 (mapped, monotonicity-preserving), HLLC |
+| Boundaries | periodic in x and z, ghost-cell extrapolation in y |
+| Chemistry | San Diego mechanism, 9 species, unity-Lewis transport |
+
+x is streamwise, y is cross-stream (the flamelet profile axis), z is spanwise. The
+resolution density follows the temporal mixing-layer DNS of Wang et al. (*Combustion and
+Flame*, 2024), with the box reduced to a quarter of theirs in each direction.
+
+The fuel stream is diluted because pure H2's sound speed is over 4x the oxidizer's, so
+reaching `Mc = 1.5` would demand a velocity split of roughly 2650 m/s. `X_H2 = 0.5` brings
+that to about 1394 m/s and moves the stoichiometric mixture fraction off the domain edge.
+
+## Initial condition
+
+`case.py` calls `flamelet_ic.py` to solve a 1-D flamelet on the cross-stream grid and write
+it as `prim.<n>.00.000000.dat` under `IC/`. `hcid=370` would extrude those files uniformly
+across z and leave the flow z-invariant, so this case uses `hcid=371`: it scales the file's
+cross-stream velocity by `1 + 0.5*cos(k_z z)` and sets the spanwise component from the
+result, giving the IC 3D content at step 0. `k_z = 2*pi/L_z` is taken over the global
+domain, so the IC does not depend on the MPI decomposition.
+
+The in-plane `(x, y)` perturbation is baked into the files by `perturb_xy`, seeded by
+`perturb_seed` in `case.py`. `IC/` is gitignored and regenerated on a fresh checkout, so
+the fixed seed is what keeps the case reproducible. Regeneration is skipped when `IC/`
+already matches the grid and the physical parameters, tracked in `.cache_key.json`.
+
+The file spacing must match the run grid. A mismatch aborts in `pre_process`, so delete
+`IC/` and let it regenerate after changing the grid or `--scale`.
+
+## Running
+
+    ./mfc.sh run examples/3D_reacting_mixing_layer/case.py -n 8
+
+`--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.
+
+The mechanism ships alongside the case as `sandiego.yaml` (UC San Diego Combustion
+Research Group, <https://web.eng.ucsd.edu/mae/groups/combustion/mechanism.html>).

examples/3D_reacting_mixing_layer/case.py

@@ -0,0 +1,188 @@
+#!/usr/bin/env python3
+"""3D temporal reacting (H2/air) mixing layer, initialized from a 1-D flamelet solve
+extruded via hcid=371. See flamelet_ic.py's module docstring for the axis convention
+(x = streamwise, y = cross-stream/profile, z = spanwise) and for how the 3D velocity
+perturbation that seeds the cascade is baked in: in-plane (x,y) via flamelet_ic.py's
+perturb_xy, spanwise (z) via a closed-form modulation added directly in Fortran.
+"""
+
+import argparse
+import json
+import os
+
+import cantera as ct
+import flamelet_ic
+
+current_dir = os.path.dirname(os.path.abspath(__file__))
+ctfile = os.path.join(current_dir, "sandiego.yaml")
+
+parser = argparse.ArgumentParser(prog="3D_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.")
+# 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).")
+args = parser.parse_args()
+
+# The 2D temporal case's H2/air mixing layer at mach_c=1.5 instead of 0.3, to probe
+# compressibility effects the subsonic case can't. Tuning below (dilution, resolution,
+# flamelet chi) is calibrated against sandiego.yaml; swapping mechanisms invalidates it.
+pressure = 101_325.0
+temperature_ox = 500.0
+temperature_fu = 300.0
+fuel = "H2"
+mole_fraction_ox = 0.21
+# Diluted with N2 rather than pure H2: pure H2's sound speed (~1318 m/s at 300 K) is over
+# 4x c_ox, so reaching mach_c=1.5 demanded delta_u ~ 2650 m/s and crashed on VCFL. X_H2=0.5
+# gives delta_u ~ 1394 m/s and moves Z_st off the domain edge to ~0.30.
+mole_fraction_fu = 0.5
+vort_thickness = 1.0e-3
+mach_c = 1.5
+num_iter = 5
+
+# 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
+# density (14 pts/delta_omega) but a quarter of their box: [15,20,10] delta_omega.
+points_per_delta = 14.0 * args.scale
+cross_min, cross_max = -10.0, 10.0  # y: 20 delta_omega
+# y is doubled relative to x/z: it carries the flame gradients and the Mach-1.5 braid
+# shocklets that crashed the uniform-14 run. dt stays at 1e-9 s.
+points_per_cross = 2.0 * points_per_delta
+stream_min, stream_max = -7.5, 7.5  # x: 15 delta_omega
+num_x = round((stream_max - stream_min) * points_per_delta)
+span_min, span_max = -5.0, 5.0  # z: 10 delta_omega
+num_z = round((span_max - span_min) * points_per_delta)
+
+t_step_stop = 20000000
+# 5000, not 10000: at the measured 1-GPU rate (3.82 s/step) a 10000-step checkpoint
+# interval is >10 h, so a job that hits walltime banks nothing.
+t_step_save = 5000
+
+# Fixed: IC/ is gitignored and regenerated on every checkout, so an unseeded perturbation
+# would make the case irreproducible.
+perturb_seed = 20260824
+
+sol = ct.Solution(ctfile)
+cross_coord, x_coord, grid = flamelet_ic.compute_grid_3d(vort_thickness, cross_min, cross_max, points_per_cross, stream_min, stream_max, num_x, span_min, span_max, num_z)
+fluid = flamelet_ic.reference_fluid_properties(sol, temperature_ox, pressure, mole_fraction_ox)
+
+ic_dir = os.path.join(current_dir, "IC")
+# 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 = {
+    "cold": not args.hot,
+    "lines": len(x_coord) * len(cross_coord),
+    "vort_thickness": vort_thickness,
+    "temperature_ox": temperature_ox,
+    "temperature_fu": temperature_fu,
+    "mach_c": mach_c,
+    "mole_fraction_ox": mole_fraction_ox,
+    "mole_fraction_fu": mole_fraction_fu,
+    "num_iter": num_iter,
+    "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,
+        temperature_ox=temperature_ox,
+        temperature_fu=temperature_fu,
+        fuel=fuel,
+        mole_fraction_ox=mole_fraction_ox,
+        mole_fraction_fu=mole_fraction_fu,
+        vort_thickness=vort_thickness,
+        mach_c=mach_c,
+        num_iter=num_iter,
+        cold=not args.hot,
+        perturb_seed=perturb_seed,
+    )
+    flamelet_ic.write_cache_key(ic_dir, cache_key)
+
+case = {
+    "run_time_info": "T",
+    "x_domain%beg": grid["x_domain_beg"],
+    "x_domain%end": grid["x_domain_end"],
+    "y_domain%beg": grid["y_domain_beg"],
+    "y_domain%end": grid["y_domain_end"],
+    "z_domain%beg": grid["z_domain_beg"],
+    "z_domain%end": grid["z_domain_end"],
+    "m": grid["m"],
+    "n": grid["n"],
+    "p": grid["p"],
+    "cyl_coord": "F",
+    "dt": 1.0e-9,
+    "t_step_start": 0,
+    "t_step_stop": t_step_stop,
+    "t_step_save": t_step_save,
+    "model_eqns": 2,
+    "alt_soundspeed": "F",
+    "mixture_err": "F",
+    "mpp_lim": "F",
+    "time_stepper": 3,
+    "avg_state": 1,
+    "weno_order": 5,
+    "weno_eps": 1e-16,
+    "mapped_weno": "T",
+    "null_weights": "F",
+    "mp_weno": "T",
+    "weno_Re_flux": "F",
+    "riemann_solver": 2,
+    "wave_speeds": 1,
+    "bc_x%beg": -1,
+    "bc_x%end": -1,
+    "bc_y%beg": -3,
+    "bc_y%end": -3,
+    "bc_z%beg": -1,
+    "bc_z%end": -1,
+    "num_patches": 1,
+    "num_fluids": 1,
+    "viscous": "T",
+    "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).
+    "chem_params%transport_model": 2,
+    "files_dir": ic_dir,
+    "file_extension": "000000",
+    "format": 1,
+    "precision": 2,
+    "prim_vars_wrt": "T",
+    "parallel_io": "T",
+    "fluid_pp(1)%gamma": 1.0 / (fluid["gamma"] - 1.0),
+    "fluid_pp(1)%pi_inf": 0.0,
+    "fluid_pp(1)%Re(1)": 1.0 / fluid["viscosity"],
+    "patch_icpp(1)%geometry": 9,
+    # hcid=371: hcid=370's file read (flamelet base state + in-plane perturbation from
+    # flamelet_ic.py's perturb_xy) plus a closed-form spanwise modulation in Fortran
+    # (src/common/include/3dHardcodedIC.fpp), so the IC is 3D from step 0. Not
+    # mixlayer_perturb: its wavenumber range is fixed in absolute units (~6 m wavelengths),
+    # meaningless for this millimeter-scale domain.
+    "patch_icpp(1)%hcid": 371,
+    "patch_icpp(1)%x_centroid": 0.5 * (grid["x_domain_beg"] + grid["x_domain_end"]),
+    "patch_icpp(1)%y_centroid": 0.5 * (grid["y_domain_beg"] + grid["y_domain_end"]),
+    "patch_icpp(1)%z_centroid": 0.5 * (grid["z_domain_beg"] + grid["z_domain_end"]),
+    "patch_icpp(1)%length_x": grid["x_domain_end"] - grid["x_domain_beg"],
+    "patch_icpp(1)%length_y": grid["y_domain_end"] - grid["y_domain_beg"],
+    "patch_icpp(1)%length_z": grid["z_domain_end"] - grid["z_domain_beg"],
+    "patch_icpp(1)%vel(1)": 0.0,
+    "patch_icpp(1)%vel(2)": 0.0,
+    "patch_icpp(1)%vel(3)": 0.0,
+    "patch_icpp(1)%pres": pressure,
+    "patch_icpp(1)%alpha_rho(1)": 1,
+    "patch_icpp(1)%alpha(1)": 1,
+    "cantera_file": ctfile,
+}
+
+print(json.dumps(case))