Add spatially supported body force + reacting mixing-layer examples

.gitignore

@@ -55,6 +55,7 @@ docs/post_process/readme.md
 docs/api/readme.md
 
 examples/**/IC/*.dat
+examples/**/IC/.cache_key.json
 examples/*batch/*/
 examples/**/D/*
 examples/**/p*

.typos.toml

@@ -22,6 +22,7 @@ TKE = "TKE"
 HSA = "HSA"
 infp = "infp"
 Sur = "Sur"
+equil = "equil"       # abbreviation for "equilibrium" (flamelet chemistry)
 chioces = "chioces"   # typo for "choices" - tests constraint key validation
 reqires = "reqires"   # typo for "requires" - tests dependency key validation
 choises = "choises"   # appears in comment explaining validation purpose

docs/documentation/case.md

@@ -271,6 +271,19 @@ The code provides three pre-built patches for dimensional extrusion of initial c
 - `case(170)`: Load 1D profile from data files
 - `case(270)`: Extrude 1D data to 2D domain
 - `case(370)`: Extrude 2D data to 3D domain
+- `case(273)`: Extrude 1D data to 2D domain, with the mom%%beg data column carrying the
+  extruded-axis (mom%%end) velocity profile instead of its own (zeroed) component. Used by
+  `examples/2D_reacting_mixing_layer` to give a temporally-evolving mixing layer a nonzero
+  streamwise velocity profile along the extruded axis, which `case(270)` cannot represent since
+  it always zeros that component.
+- `case(274)`: Load a full 2D `(x, y)` field with no extrusion -- one data file per variable,
+  `(m_glb+1)*(n_glb+1)` lines each in x-major order, covering all primitive variables
+  directly (unlike `case(270)`/`case(273)`, no component is zeroed or repurposed). Used by
+  `examples/2D_spatial_reacting_mixing_layer` for a spatially-evolving mixing layer, where
+  the cross-stream profile must vary with the streamwise coordinate too (via the
+  `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.
 
 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.
@@ -1055,20 +1068,27 @@ This parameter enables the use of true `pi_\infty` in bubble dynamics models whe
 
 ### 13. Body Forces
 
-| Parameter         | Type    | Description                                |
-| ---:              | :---:   | :---                                       |
-| `bf_x[y,z]`       | Logical | Enable body forces in the x[y,z] direction |
-| `k_x[y,y]`        | Real    | Magnitude of oscillating acceleration      |
-| `w_x[y,z]`        | Real    | Frequency of oscillating acceleration      |
-| `p_x[y,z]`        | Real    | Phase shift of oscillating acceleration    |
-| `g_x[y,z]`        | Real    | Magnitude of background acceleration       |
+| Parameter            | Type    | Description                                                              |
+| ---:                 | :---:   | :---                                                                     |
+| `bf_x[y,z]`          | Logical | Enable body forces in the x[y,z] direction                               |
+| `k_x[y,z]`           | Real    | Magnitude of oscillating acceleration                                    |
+| `w_x[y,z]`           | Real    | Frequency of oscillating acceleration                                    |
+| `p_x[y,z]`           | Real    | Phase shift of oscillating acceleration                                  |
+| `g_x[y,z]`           | Real    | Magnitude of background acceleration                                     |
+| `bf_spatial_support` | Logical | Enable spatially supported body force (Wei & Freund, JFM 2005)           |
 
 `k_x[y,z]`, `w_x[y,z]`, `p_x[y,z]`, and `g_x[y,z]` define an oscillating acceleration in the `x[y,z]` direction with the form
 
 \f[ a_{x[y,z]} = g_{x[y,z]} + k_{x[y,z]}\sin\left(w_{x[y,z]}t + p_{x[y,z]}\right). \f]
 
 By convention, positive accelerations in the `x[y,z]` direction are in the positive `x[y,z]` direction.
 
+When `bf_spatial_support` is enabled, the body force is given a localized spatial envelope following Wei & Freund (JFM, 2005), parameterized through the `spatial_bf` derived type with attributes `amp`, `x_centroid`, `y_centroid`, `conv_vel`, `sigma`, and per-mode arrays `freq(1:8)`, `phase(1:8)`.
+Unlike the domain-wide `k/w/p/g` forcing above, it seeds a downstream instability at a fixed streamwise location with a ladder of forcing frequencies.
+`x` is the streamwise direction by convention for this forcing: the `x` component advects with the mean flow via a `conv_vel*t` term, while the cross-stream (`y`) component does not.
+The source is constructed from a streamfunction, so it is divergence-free, and the corresponding `u*f` work term is added to the energy equation.
+This forcing is implemented for **2D cases only**; enabling it in 1D or 3D is rejected at startup.
+
 ### 14. Magnetohydrodynamics (MHD)
 
 | Parameter                | Type    | Description                                              |

examples/2D_reacting_mixing_layer/case.py

@@ -0,0 +1,154 @@
+#!/usr/bin/env python3
+"""Temporal reacting (H2/air) mixing layer, initialized from a 1-D flamelet solve
+extruded via hcid=273. See flamelet_ic.py's module docstring for the required MFC
+axis convention (x = cross-stream/profile, y = streamwise/periodic/extruded).
+"""
+
+import argparse
+import json
+import os
+
+import cantera as ct
+import flamelet_ic
+
+current_dir = os.path.dirname(os.path.abspath(__file__))
+ctfile = "h2o2.yaml"
+
+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).")
+args = parser.parse_args()
+
+# Physical parameters: representative temporal H2/air mixing layer.
+pressure = 101_325.0
+temperature_ox = 500.0
+temperature_fu = 300.0
+fuel = "H2"
+mole_fraction_ox = 0.21
+mole_fraction_fu = 1.0
+vort_thickness = 1.0e-3
+mach_c = 0.3
+num_iter = 5
+
+# Grid: x = cross-stream (flamelet profile axis), y = streamwise (periodic, extruded).
+cross_min, cross_max = -10.0, 10.0
+points_per_cross = 8.0 * args.scale
+stream_min, stream_max = -5.0, 5.0
+num_y = 32  # weno_order=5 needs n+1 >= 25 (num_stcls_min*weno_order)
+
+t_step_stop = 200
+t_step_save = 100
+
+sol = ct.Solution(ctfile)
+cross_coord, grid = flamelet_ic.compute_grid(vort_thickness, cross_min, cross_max, points_per_cross, stream_min, stream_max, num_y)
+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(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,
+}
+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,
+        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,
+    )
+    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"],
+    "m": grid["m"],
+    "n": grid["n"],
+    "p": 0,
+    "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": -3,
+    "bc_x%end": -3,
+    "bc_y%beg": -1,
+    "bc_y%end": -1,
+    "num_patches": 1,
+    "num_fluids": 1,
+    "viscous": "T",
+    "chemistry": "T",
+    "chem_params%diffusion": "T",
+    "chem_params%reactions": "T",
+    "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": 3,
+    "patch_icpp(1)%hcid": 273,
+    "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)%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)%vel(1)": 0.0,
+    "patch_icpp(1)%vel(2)": 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))