Repository navigation
Add chem_params%transport_dt: limit dt by the local chemistry transport - #1953
sbryngelson wants to merge 12 commits into
Conversation
The viscous dt limit used the constant fluid_pp%Re, while chemistry integrates temperature-dependent viscosity, conductivity and species diffusion explicitly. Near hot walls and flames those grow by an order of magnitude or more. With transport_dt = T the viscous limit uses 1/(rho*D) for the largest of mu/rho, lambda/(rho*cp) and the mixture-averaged species diffusivities. Off by default.
There was a problem hiding this comment.
🟡 Changes recommended
IGR can access unallocated fields, and active non-Newtonian or bulk-viscosity limits can be discarded.
2 open findings
What changed in this PR
Adds chemistry-aware transport limits to CFL time-step calculation.
Changes:
- Introduces and validates
chem_params%transport_dt. - Computes local viscous, thermal, and species diffusivity limits.
- Adds documentation, validator coverage, and a golden regression case.
| File | Description |
|---|---|
toolchain/mfc/test/cases.py |
Adds transport-limited chemistry case. |
toolchain/mfc/test_case_validator.py |
Tests parameter constraints. |
toolchain/mfc/params/descriptions.py |
Describes the parameter. |
toolchain/mfc/params/definitions.py |
Registers the parameter. |
toolchain/mfc/case_validator.py |
Validates required configuration. |
tests/88B73ED4/golden.txt |
Adds regression results. |
tests/88B73ED4/golden-metadata.txt |
Records golden generation metadata. |
src/simulation/m_time_steppers.fpp |
Applies local transport limits. |
src/simulation/m_sim_helpers.fpp |
Enables the VCFL path. |
src/simulation/m_global_parameters.fpp |
Adds simulation default. |
src/pre_process/m_global_parameters.fpp |
Adds preprocessing default. |
src/post_process/m_global_parameters.fpp |
Adds postprocessing default. |
src/common/m_derived_types.fpp |
Extends chemistry parameters. |
src/common/m_chemistry.fpp |
Computes effective transport Reynolds numbers. |
docs/documentation/case.md |
Documents usage and constraints. |
🧠 Review effort: Balanced
Give feedback about Copilot approvals in this survey to enter a drawing for a $150 gift card.
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## master #1953 +/- ##
==========================================
+ Coverage 61.75% 61.81% +0.05%
==========================================
Files 86 86
Lines 22787 22802 +15
Branches 3357 3355 -2
==========================================
+ Hits 14073 14095 +22
+ Misses 6225 6218 -7
Partials 2489 2489 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Explicit diffusion is stable while the per-direction diffusion numbers D dt/dx^2 sum below the limit, so in d dimensions dt <= cfl dx^2/(d D). Without the factor a 2D shock-driven curtain of hot, freely moving reacting particles diverged at 5 us; with it the run completes.
… rejecting them The transport limit replaced Re outright, dropping the Herschel-Bulkley mu_max (Re(1)) and bulk viscosity (Re(2)) that the fluxes still apply. Take the tighter of each; Cantera replaces only the constant shear Re in the fluxes, so Re(1) is kept only for non-Newtonian fluids. Without viscous, Re holds dflt_real and is not combined.
# Conflicts: # src/simulation/m_time_steppers.fpp
|
Validation from a particle-resolved IB campaign on Tuolumne (MI300A); build = master a6add01 + #1955, #1956, #1957, #1962, #1969 and this PR (revision before f72eb74, whose HB/bulk-viscosity bounds are inactive here). On/off spot check: 2D cylinder, Re 40, M = 0.2, 20 cells/D, CFL 0.5, inert moist air,
The steady result is unchanged to 5–6 digits; the transport limit took 1.5× the steps. (The corrected Nu also uses the Tw = T∞ twin, whose on/off runs were averaged over different windows, hence the 0.15 %.) Without transport_dt, the default limit is not enough for these IB chemistry runs. Same case, CFL 0.5, transport_dt off:
Cost observation. The limit takes max_k 𝒟_k over all species, so in air at 300 K it is set by atomic-H diffusivity (≈ 7.8 ν) although Y_H = 0 everywhere in these runs (gas reactions off). This roughly doubles the step count at 80 cells/D compared with a limit from the species present. In a Mach 2, Re 5000 adiabatic sphere, the low-pressure wake drove dt to about 4e-11 s and made the run about 10× costlier, with diffusion otherwise negligible; we ran that case with transport_dt off. Weighting by the species actually present (or skipping species below a mass-fraction floor) might keep the stability benefit at lower cost. That is a suggestion, not a blocker: an absent species can appear once reactions are on. Made with Claude Code. |
The explicit energy equation diffuses temperature at lambda/(rho*cv) on the grid scale, so lambda/(rho*cp) let dt exceed the stability limit by gamma (unstable at the default cfl_target with transport_model = 2). Made with Claude Code.
|
@sbryngelson I pushed 3faed31, which fixes a stability bug in Bug. The heat-conduction term of the limit used λ/(ρ c_p). The explicit energy equation conducts the grid-scale temperature mode at λ/(ρ c_v), because a cell can't expand within one step. So the limit allowed a dt too large by a factor of γ. Who it affected.
The non-chemistry thermal limit in Evidence. 1D N₂ temperature step with h2o2.yaml, mixture-averaged quantities from Cantera, scanning
The RK3 diffusion limit is 0.628. Model 1 with H present is unchanged: stable to 0.62, unstable at 0.64. Change. Tests. The Chemistry suite passes 23/23, including the Found during verification work with Claude Code. |
Lines of Code
|
Summary
Adds
chem_params%transport_dt(logical, defaultF).With chemistry, the viscous fluxes use temperature-dependent μ and the diffusion step integrates λ and the mixture-averaged species diffusivities explicitly, but the adaptive dt only checks the viscous limit with the constant
fluid_pp%Re(the thermal TCFL applies only to the non-chemistryheat_conductionmodel). Mixture transport grows steeply with temperature: in air with CO near an 1800 K wall the thermal diffusivity is ~30× and the largest species diffusivity (H) ~160× the 300 K kinematic viscosity that the current limit uses.With
transport_dt = T,s_compute_dtreplacesRein the viscous limit by1/(d ρ D), D = max(μ/ρ [viscous], λ/(ρ c_p) and max_k D_k [chem diffusion]) and d the number of dimensions, from a news_compute_transport_dt_reinm_chemistry. The existing VCFL branch is reused unchanged. The factor d is needed because explicit diffusion is stable only while the per-direction diffusion numbers D Δt/Δx² sum below the limit; an earlier version of this PR omitted it.Not changed:
run_time.inf's VCFL column still usesfluid_pp%Re; the routine that writes it has no temperature field.Testing
1D -> Chemistry -> Dual Isothermal Wall Gradient -> Transport-Limited dt(88B73ED4): the existing isothermal-wall case on a 10 µm domain, where the H₂ diffusion limit (dt ≈ 2.2e-11 s) is tighter than the acoustic one (≈ 4.4e-11 s). All 40 steps report(VCFL)../mfc.sh test --no-mpi -o Chemistry: 16 pass, unchanged (option off by default).TestChemistryTransportDt(requires chemistry,cfl_adap_dt/cfl_const_dt, andviscousorchem_params%diffusion).s_compute_dtkernel goes from 0.87 ms to 9.4 ms per call with the option on (mixture-averaged diffusivities in every cell), 0.16% → 1.6% of GPU time; time per step 0.950 s either way. No communication is added: still the one 5-entry allreduce per step. Off, the cost is a branch../mfc.sh lintpasses excepttest_thermochem.py::test_precision_and_directives[dp-acc], which also fails without this change here (LLNL gfortran 13.3.1's nvptx offload lackslog10).Tested on LLNL Tuolumne (CCE 19 CPU; OpenMP offload on MI300A).
Why d matters (2D): a Mach 2 shock driving a curtain of 38 hot, freely moving, reacting graphite particles diverged at t = 5.0 µs without
transport_dt. Withtransport_dtand the limit Δx²/(2D) the same case runs to the end (11.5 µs), its step set by the transport limit (VCFL). The blowup is a checkerboard, odd–even in both space and time, in the gas next to a particle: the explicit-diffusion mode at its stability limit. In the hot, low-pressure wake (1300–1700 K, 0.3–0.5 bar) the largest diffusivity is H-atom diffusion, D ≈ 5e-3 m²/s. With Δx = 3.1 µm and the acoustic step dt = 8.9e-10 s, D·dt/Δx² ≈ 0.46 per direction (stable) but ≈ 0.93 summed over x and y (unstable). Without d the transport limit, ≈ 9.8e-10 s, would sit above the acoustic step and never activate; with d it is ≈ 4.9e-10 s.The 1D test case and its golden are unaffected (d = 1);
./mfc.sh test -o Chemistrypasses (16).