Tier-2 GPU-AMR-CFD

The flow engine for dispersion (Phase 2) and fire (Phase 3). Variable-density low-Mach formulation; passive scalar dispersion is handled via the auxiliary-variable transport hooks.

Time integration

SDC (Spectral Deferred Correction). Default 2 corrector iterations. Used in V&V; production may drop to 1 corrector once accuracy parity holds.

Pressure projection

Approximate-projection MAC (face-centred) + nodal projection; geometric MLMG bottom solver, optional Hypre.

Advection

AMReX-Hydro Godunov PPM (default).

Diffusion

Implicit Crank-Nicolson, MLMG.

Reactions

Stiff CVODE (PelePhysics ReactorCvode) when do_react = 1.

LES

Smagorinsky default; aux-LES wired by patch 0001.

Subcycling

Level-by-level subcycling (PeleLMeX default).

Non-default options

• peleLM.sdc.nIter = 2
• peleLM.do_react = 0 — 1 for fire/explosion
• peleLM.do_temporals = 1 — time-series probes
• peleLM.steady_state_atol — early-stop when Δstate plateaus
• peleLM.advection_scheme · peleLM.cfl — see Convection
• peleLM.aux_les_contrib = 1 patch 0001
• peleLM.D_x_ramp · D_x_u_cal patch 0002
• peleLM.turbinflow_clip_vel patch 0003
• peleLM.turbinflow_lowpass_cutoff patch 0005

Refinement

2× per level (configurable); n_error_buf = 2 base cells.

Tagging

Concentration-gradient tagger (default); obstacle-aware tagger added with patch 0005 helper.

Subcycling

Level-by-level (PeleLMeX default).

Block factor

blocking_factor = 8 — every box dim must be divisible by 8.

MLMG tolerances

32K cells: atol=1e-8, maxiter=200. 262K + turb inflow: atol=1e-1, maxiter=20 (post-projection divu ≈ 1e-4 acceptable).

Auto Hypre

Build flag detects when to fall through to Hypre on native MLMG stall.

Non-default options

• amr.max_level = 0 — 0 uniform; 2 standard for AMR
• amr.n_cell = 128 64 32 — base grid
• amr.blocking_factor = 8
• amr.max_grid_size = 32 — maxgrid32 in α-ceiling workaround
• amr.regrid_int = 2
• amr.n_error_buf = 2
• amr.max_dt = 0.12 s
• amr.max_step — guardrail; for 1500 s @ dt=0.12 set ≥ 12500

AMReX EB cut cells with conservative state redistribution. Sharp interface. Required for canyon-flow and surface-pressure benchmarks (CODASC, EB-001 → EB-004).

Numerical method

State redistribution + small-cell combinatorial merging; AMReX-Hydro Godunov MOL fallback near cut faces; flux redistribution for conservation.

Wall model

no-slip default; upcoming wall-function for coarse-grid accuracy without cell-merging penalty.

Geometry kinds

box · sphere · cylinder · plane · STL signed-distance.

Non-default options

• eb2.geom_type = "all_regular" → "box" / "sphere" / "stl" / …
• eb2.stl_file — STL path when geom_type=stl
• eb2.small_volfrac = 1e-3 — small-cell merge threshold
• mac_proj.bottom_solver = "native" — set "hypre" for stiff EB cases (Hypre v2.28 in tree)
• nodal_proj.bottom_solver = "native" — same as above for nodal projection

Volume-of-Solid (VOS) penalization. Runtime-configurable via JSON → ParmParse — drop boxes / cylinders / spheres without a mesh-generation step. Augmented with an immersed-boundary wall drag (Item 1, 2026-05-04), a wake-spike softener (Item 3) and accelerated relaxation by substepping (Item 4).

Penalty form

F = −α · u inside obstacle (NOT α/dt — the inverse-dt form creates a stiffness feedback loop).

α default

α = 100 [1/s] (relaxation time ≈ 0.05 s for u = 4 m/s wind).

α ceiling

α ≥ 1000 collapses MLMG at sim time ≈ 280 s — cap α ≤ 500. Higher α = sharper wall but slower projection (research-mode trade).

Numerical method

Explicit pointwise penalty + optional semi-implicit Brinkman update (5× speed factor on obstacle-heavy cases). Substepping = M sub-cycles within one outer dt.

IBM wall drag

Item 1: velocity correction in a 2-cell band outside the obstacle, closes the Brinkman ramp under-prediction (≈ 9× concentration gap vs EB at x = 200 m → within 20%).

Initial condition

Mandatory: u = 0 inside obstacles. Without this, log-law IC creates u ≈ 4 m/s inside solid; penalty pulse generates a 1600 m/s shock at corners → CFL collapse to dt ≈ 3e-6 s.

CFL ceiling

≈ 0.024 s on 8 m grid (limited by flow acceleration around corners, not by Brinkman itself).

Non-default options

• prob.use_brinkman = 0 → 1 enables
• prob.brinkman_alpha = 100.0 — relaxation rate [1/s]; cap at 500
• prob.brinkman_substeps = 1 — Item 4; raise to 4–8 for accelerated transient
• prob.brinkman_semi_implicit = 0 — 1 enables semi-implicit update (5× faster)
• prob.wall_drag = 0 — 1 enables Item 1 IBM correction
• prob.wall_drag_band_cells = 2
• prob.wake_soft = 0 — Item 3 wake-spike softener
• prob.brinkman_n_obstacles · brinkman_obstacle_type[i] · brinkman_lo[i] · brinkman_hi[i] · brinkman_centre[i] · brinkman_radius[i]

Model

Smagorinsky · static. Cs = 0.18 neutral / 0.10 stable F (workaround for class F/A while gravity-projection limitation persists).

Aux-var coupling

patch 0001 adds μt/Sct,n to diff_aux_cc[n] — without it, the tracer carries only laminar viscosity and LES is defeated for dispersion.

Sc_t

0.7.

Numerical method

Cell-centred strain-rate from finite-difference velocity gradients; μ_t averaged to faces for diffusion solve.

Non-default options

• peleLM.use_les = 1
• peleLM.les_model = "Smagorinsky" — only model in tree today
• peleLM.smagorinsky_constant = 0.18 — 0.10 for class F
• peleLM.turb_Sc = 0.7
• peleLM.aux_les_contrib = 1 patch 0001

Pre-computed homogeneous-isotropic-turbulence box, time-superposed on the inflow plane (Klein 2003 family). Drives the realistic LES boundary that the dispersion benchmarks require. Two safety patches in tree.

Numerical method

Box read at start; bilinear interpolation in time; two-mode urms scaling; Sc_t-coupled mass continuity at inflow plane.

ParmParse pitfall

module uses bare turbinflows = ..., NOT peleLM.turbinflows. Wrong namespace silently runs laminar.

Velocity clamp

patch 0003 — clip per-cell turbulent perturbation to clip_vel × urms. Required at urms ≥ 0.6 to avoid step-1028 collapse even with Hypre.

Spectral filter

patch 0005 — low-pass on the HIT box for the AMR + turb path. Default off.

Non-default options

• turbinflows = inflow_lo_x — bare key (no peleLM. prefix)
• turbinflow.{name}.turb_file — path to HIT box
• turbinflow.{name}.iso_scale = 1.0
• turbinflow.{name}.urms0 = 0.4 m/s — 10% TI at 4 m/s wind
• turbinflow.{name}.k0 = 4
• turbinflow.{name}.N = 64
• peleLM.turbinflow_clip_vel = 0.5 patch 0003
• peleLM.turbinflow_lowpass_cutoff = 0.0 patch 0005

Method

Godunov finite-volume with PPM (Piecewise-Parabolic Method, 4th-order) slope limiter — default for momentum, species, aux scalars.

Alternates

PLM (Piecewise-Linear, 2nd-order, more diffusive but very robust); MOL (Method-of-Lines) used near EB cut cells.

Conservation near EB

Flux redistribution + state redistribution merge small cells, preserving global mass/momentum.

Reflux

Level-by-level subcycling refluxes coarse-fine fluxes at the end of each fine sub-step.

Non-default options

• peleLM.advection_scheme = "Godunov_PPM" — "Godunov_PLM" / "MOL"
• peleLM.cfl = 0.9 — bumped from 0.5 post-sweep, verified stable
• peleLM.max_dt · peleLM.init_dt · peleLM.dt_change_max — dt growth governor
• amr.max_dt = 0.12 — required for VC-003 (no-AMR 8 m) to avoid MLMG divergence

Time scheme

Crank-Nicolson (CN, 2nd-order) for momentum / species / thermal / aux diffusion.

Linear solver

AMReX MLMG (geometric multigrid). V-cycle default; W-cycle and semi-coarsening selectable.

Bottom solver

Native by default; Hypre v2.28 available (no MPI) when MLMG stalls.

Transport model

Constant (not Fuego). All values in CGS for Pele compatibility. fixed_Le = 0 mandatory — otherwise Constant-transport diffusivity is silently ignored in favour of Sc-derived D.

D(x) ramp

patch 0002 — distance-dependent eddy diffusivity, scaled by u_ref / u_cal (calibration wind 4.0 m/s).

Non-default options

• transport.const_viscosity = 30.0 [Poise = 3 Pa·s] · ν = 2.5 m²/s before LES
• transport.const_diffusivity = 25000.0 [cm²/s = 2.5 m²/s, Pasquill-D eddy D]
• transport.const_conductivity = 4.3e8 [erg/cm/s/K]
• peleLM.fixed_Le = 0 — mandatory
• peleLM.D_x_ramp = 1 · peleLM.D_x_u_cal = 4.0 — D(x) ramp; winner D_far=20, transition x=200 m (sweep 2026-05-03)
• mlmg.atol = 1e-8 · mlmg.maxiter = 200 — relax to atol=1e-1 / maxiter=20 for 262K + turb inflow
• mlmg.bottom_solver = "native" → "hypre" for stiff cases
• mlmg.semicoarsen_dir · mlmg.use_W_cycle — escape hatches for anisotropic stiffness

Unbiased ray tracing through participating media. Default for all VC-R cases — pelerad.mc_wall_flux = 1 auto-writes mc_wall_flux.csv. Higher ray counts are research mode.

Mode 1

Full-field incident radiation everywhere; G(x) MultiFab.

Mode 2

Monitor probes (directional or omni). Position + optional surface normal → q_inc per probe (analytical κ at probe).

Mode 3

Hemisphere wall flux — exact integration over wall hemisphere. Default-on per case.

Mode 4

q_inc(x) field + iso-contour — directional q_inc field + iso-line export (built 2026-04-29).

Numerical method

Woodcock tracking default (κ_max auto-computed at start); analog ray-marching available. Anisotropic scattering supported with Henyey-Greenstein phase function.

Wall flux formula

MC uses the exact hemispheric integral; P1 uses Marshak (not Lambert — Lambert was the bug).

Non-default options

• pelerad.mc_n_rays — rays per sample point / probe; raise for research-mode V&V
• pelerad.mc_woodcock = 1 — 0 = analog ray-march
• pelerad.mc_wall_flux = 1 — Mode 3, default-on per case
• pelerad.mc_skip_wall_flux = 0 — opt-out
• pelerad.mc_include_net_flux = 0 — opt-in
• pelerad.mc_monitor_positions · mc_monitor_normals · mc_monitor_T_surf · mc_monitor_eps_surf · mc_monitor_samples — Mode 2 probes
• pelerad.mc_compute_q_inc_field = 0 — Mode 4 toggle
• pelerad.mc_q_inc_field_normals — per-direction iso-contour
• pelerad.mc_use_sobol = 1 · mc_use_hot_ray = "auto" — Sobol on by default; hot-ray auto-fires for narrow sources (see 2.5 variance reduction)

Gas absorption

WSGGM from RADCAL. Composition selectable via pelerad.wsggm_composition (e.g. H₂O+CO₂ at AR Confluence ratios from the Bressloff thesis).

Auto-rule

If spectral_model not set: gray for κ-uniform / single-band cases, WSGGM otherwise. pelerad.spectral_model = gray|wsggm always wins.

Soot

Rayleigh small-particle approximation; volume-fraction field consumed by κ.

Droplet Mie

upcoming — Mie tables drafted; integration with cryogenic spray pending.

WSGGM accumulation bug

Fixed 2026-04-28 (Phase D).

Both samplers are on by default; the user only sees them as a toggle when they want to disable variance reduction for a baseline reference run.

Sobol primary-ray

Default for every MC run. 2-D Sobol sequence on (azimuth, polar) for Modes 1/2/3/4. 44× variance reduction at N = 1024 on smooth integrands. Replaces pseudo-random primary directions with no per-case tuning required.

Hot-ray adaptive

Auto-selected for narrow / localised emitters. Cumber (2009) two-pass importance sampling on top of Sobol; slope −1.4 observed on localised-source cases. Decision rule: enabled when the emitter footprint covers fewer than ~1% of domain cells (configurable via pelerad.hot_ray_auto_threshold) — fires automatically for jet-fire flames, vent-stack ports, and small-radius pool fires.

Override

Auto-detection can be forced on or off per-case if the user disagrees with the heuristic.

Non-default options

• pelerad.mc_use_sobol = 1 — default; set 0 only for baseline pseudo-random reference runs
• pelerad.mc_use_hot_ray = "auto" — "1" force on, "0" force off
• pelerad.hot_ray_auto_threshold = 0.01 — emitter cell-fraction below which hot-ray auto-enables
• pelerad.hot_ray_pilot_rays · pelerad.hot_ray_passes — Cumber two-pass tuning

Explicit Boussinesq buoyancy for dense-gas (Cl₂, NH₃, LPG vapour) and light-gas (H₂, CH₄) effects. Multi-component (MW) coupling fixed in 21 unit tests. Beyond |Δρ/ρ| ≈ 0.03 the low-Mach explicit-Boussinesq stability ceiling is hit; anelastic reformulation deferred.

Form

Boussinesq density perturbation from concentration + temperature; gravity term added to momentum.

MW coupling

Δρ/ρ = release_conc · (MW/MW_air − 1) + thermal term.

Workaround

keep release_conc · |MW/MW_air − 1| < 0.03.

Non-default options

• prob.boussinesq = 0 — 1 enables
• prob.MW_release = 29.0 g/mol — air-equivalent default
• prob.MW_air = 28.97 g/mol
• peleLM.gravity = 0 0 0 — class F/A workaround keeps gravity off; stability is encoded via Cs/D

Implemented in pelelmex_prob.H::initdata() + time-resolved injection in the SDC RHS. Q is read by Tier-1 source-term library and matched at postprocess via virtual-source σ shift.

Continuous

Steady Q at finite-radius origin; matches Tier-1 Q with virtual-source σ shift.

Instantaneous puff

Finite-time top-hat (Patch 0002 part A — superseded after virtual-source fix landed in postprocess).

Volumetric

Scattered injection across cylinder/box for elevated stacks.

Numerical method

Mass + species + scalar source vector built per cell from analytic mask; integrated by SDC alongside the rest of the RHS.

Non-default options

• prob.release_x / y / z — centre of source
• prob.release_radius = 16.0 m — PG default; reduce for stack
• prob.release_conc = 1.0 — mass fraction at source
• prob.release_duration = -1.0 — negative = continuous
• prob.wind_speed = 4.0 m/s — consumed by D(x) ramp
• prob.hydrostatic_ic = 0 — class F/A workaround

Solver

PeleLMeX low-Mach + AgniKawach pelelmex_prob.H.

Tracer

Single passive aux variable; multi-species mixtures iterate Tier-2 in QRA loop.

Turbulence

Smagorinsky LES (Cs = 0.18) + HIT TurbInflow (urms = 0.4 m/s default).

Geometry

EB · Brinkman+IBM · PDR — composable; CODASC uses EB, industrial sites use Brinkman+IBM + PDR.

Buoyancy

Boussinesq+MW (within |Δρ/ρ| ≲ 0.03 ceiling).

Source

Volumetric Gaussian-radius release. Virtual-source σ shift in postprocess.

Validated config

VC-005 inputs.pg-turb-8m-optimal: 8 m uniform grid, 262K cells, 1024×512×256 m domain. Centerline ratio at x=500m: 0.84. σ_y at x=104m: 16.5 vs PG-D 16.7.