Repository navigation
Add thermal_bc = 3: a lumped solid whose wall temperature follows its heat balance (after #1821, #1956) - #1957
Draft
sbryngelson wants to merge 60 commits into
Draft
sbryngelson wants to merge 60 commits into
sbryngelson wants to merge 60 commits into
Conversation
Conflict in m_ibm.fpp: master added the alpha_q, alpha_rho_q and e_q locals for per-phase EOS evaluation, this branch added W_species and the surface-reaction locals. Both sets are kept, and both appear in the kernel's private clause - a scalar assigned in the loop but absent from that list races under OpenMP offload.
get_slug hashed the phase name, which is conventionally 'gas', so two cases with different mechanisms shared one build and the second ran against the first's species set. This branch is the first to carry two gas mechanisms: the 3D reacting mixing layer's sandiego.yaml has nine species and the carbon surface case's reduced GRI mechanism has eleven, so whichever built first decided sys_size for both, and the mixing layer wrote 34 output files where its golden has 30. Reproduced by running both cases together, which is also why each passes alone. The build already reports the mechanism by source when it prints Chemistry:; this makes the key agree with what it prints.
…ombined with W_species and Ys_s were declared dimension(num_species) outside the USING_AMD guard, while Ys_IP and Ys_g inside it carry the padded literal. Ys_g(:) = 2*Ys_s(:) - Ys_IP(:) is then a shape mismatch in any generic amdflang build, at any species count: with the literal at ten and a nine-species mechanism it is ten against nine and the compile fails. Both arrays now follow the guard, and the four whole-array assignments are pinned to 1:num_species so they do not depend on the padding happening to match. Separate from the ten-species ceiling itself, which this does not lift -- an eleven-species mechanism still needs case optimization on AMD, or MFlowCode#1848.
# Conflicts: # src/simulation/m_ibm.fpp
…FlowCode#1821) Ghost-state reconstruction (m_ibm.fpp): - The linear mirror phi_g = 2*phi_s - phi_IP goes negative whenever the surface value sits below half the image-point value. On the PR's own example it wrote negative O and OH ghost mass fractions 114,377 times; a cold wall drives the ghost temperature negative outright. Replace with a convex blend toward the surface value, phi_g = phi_s + theta*(phi_s - phi_IP): theta = 1 is the existing mirror, theta = 0 is the first-order Dirichlet ghost of Gibou et al. (2002). One theta across all species keeps sum(Y) = 1 exactly with no clamping or renormalization; temperature has its own theta so a trace radical cannot throttle the thermal BC. Temperature is held inside [T_surface_min, T_surface_max] -- the NASA polynomial fit range, now in m_constants and shared with the checker. - Close sum(Y) = 1 on the most abundant species, re-picked each Newton iteration (Cantera's evalSurfLarge), not on the mechanism's last species. The dropped species absorbs every other balance's roundoff; for the shipped mechanism that species was H2O2, a trace radical. - Skip ghost points with zero levelset distance instead of dividing by it. NaN residuals defeat the pivot test, so the solve reported success. - Count non-converged and ill-posed surface solves and print them at the end of the run, following s_report_pressure_relaxation. Silence means the surface chemistry engaged everywhere it was asked for. - Loosen the Newton tolerance 1e-8 -> 1e-6. The forward-difference Jacobian makes convergence linear near the root; Cantera uses the same 1e-7 step against 1e-4. - Size Ys_s from AMD_NUM_SPECIES_MAX like its neighbours, not a literal 10. Checker: thermal_bc /= 0 now requires chemistry and rejects inj_species > 0, since only the chemistry reconstruction honours it; Twall must lie in the thermodynamic window. Toolchain: hash the surface mechanism into the simulation build slug; reject sticking-coefficient and Blowers-Masel rates by class instead of a hasattr probe that let them through as plain Arrhenius; reject coverage dependencies the generator does not emit; keep searching mechanism candidates after one fails. Tests: regenerate the ibm_reacting_surface Example golden (its previous values encoded the negative ghost mass fractions) and add a cold-wall case that exercises the temperature-limited branch the Example never reaches. Review and fixes assisted by Claude Code.
The case I added ran to t = 2e-5 and its golden did not survive a change of compiler: generated under nvhpc 25.11, it missed GNU and every other nvhpc release by ~1e0 relative in energy -- 22 failing CI lanes, none of them a real regression. A cold wall is exactly the state that cannot be run long and compared tightly. It pins the ghost temperature against the 200 K floor of the NASA fits, and that state feeds back through the stiff surface and gas kinetics, so small differences in how each compiler evaluates them diverge. Stopping at one step keeps the answer set by the reconstruction rather than by accumulated kinetics, which is the part this test exists to pin: 193 ghost updates still take the limited branch, at theta_T = 0.102. Golden regenerated.
Every guard in generate_surface_thermochem exists because the alternative is silent: a rate emitted wrong by orders of magnitude, or a rate law missing terms, with nothing printed and a run that completes. The integration goldens cannot see any of it -- they run one mechanism whose every reaction happens to be of the one supported kind. Cantera's own ptcombust.yaml is the fixture because it carries all three unsupported forms at once: 5 sticking-coefficient rates, 2 coverage-dependent rates, and surface-site species in the stoichiometry. Against the guard this PR replaced, 7 of its 24 reactions were accepted and would have been emitted as plain Arrhenius -- including sticking probabilities of 0.023 and 1 used as prefactors. The example's carbon mechanism is tested end to end through get_cantera_surface, so the adjacent-phase resolution is covered too, and a chemistry case with no surface mechanism is tested to still emit the no-op module that m_ibm.fpp imports unconditionally.
…ecies sizing Two coverage gaps left over from the review, both closed by putting the check where a harness already exists rather than building a new one. Input constraints -> case_validator.py. lint_source.py states the rule: a constraint between case-file parameters belongs in the Python validator, not in m_checker.fpp, because the Fortran copy cannot be unit tested and the two drift. s_check_inputs_ib_injection was exempt from that rule wholesale -- the allowlist entry reads "num_species is populated by Cantera at runtime" -- which is true of exactly one of its ten constraints. The other nine are relations between case-file values, including the three this review added. They now live in check_ibm with unit tests, and the Fortran keeps only the num_species bound. The Twall window is read from m_constants.fpp rather than repeated: a new parse_fortran_real_constants does for real(wp) parameters what the existing parser does for integers, so the validator rejects against the same numbers the solver clamps to. AMD species sizing -> a source lint. Without case optimization LLVMFlang cannot size an automatic array by num_species, so those branches use a fixed AMD_NUM_SPECIES_MAX; a smaller literal is a silent buffer overrun, not a compile error, and needs amdflang plus a large mechanism to reproduce -- which no machine here has. check_amd_species_array_sizes catches it by reading the source instead. Verified against the original defect: reintroducing dimension(10) for Ys_s produces exactly one error naming the file, line and fix. 711 toolchain tests pass; both reacting-surface goldens unchanged.
I added this case to pin the temperature side of the ghost-state limiter. Its golden does not survive a change of compiler, and two attempts did not fix that: generated under nvhpc 25.11 at t = 2e-5 it missed GNU and every other nvhpc release by ~1e0 relative in energy, and shortening it to a single step only brought that to 1.2e-3, still past the 1e-3 tolerance. It has red-lighted every CI run since. The obvious explanation is wrong, so this is a withdrawal rather than a diagnosis. The limiter parks the ghost temperature at 0.1*T_s + 0.9*T_min = 201 K, one degree above the NASA fit floor, which looked like the culprit -- but the thermodynamic state is no worse conditioned there than at 4900 K, both responding ~1e-12 to a 1e-12 relative nudge in temperature. Whatever makes this case compiler-sensitive, it is not simply evaluating the fits at their low edge. What is lost is narrower than it looks. The auto-registered ibm_reacting_surface Example already exercises the species side of the same limiter hard -- theta_Y is about 0.006 across ~114k ghost-point updates -- so only the theta_T branch is now uncovered, and its arithmetic is four lines. MFlowCode#1892 is the right home for it: with the surface solver in a module of its own, this is a unit test on a function, with no CFD and no compiler sensitivity in it.
A mechanism is compiled into the binary, so each one the suite uses costs a whole extra simulation link. On Frontier AMD's GPU lane that is the binding constraint: amdflang cannot link the base and chemistry variants serially inside the 1h59m walltime, which is why the build is already split across two concurrent SLURM jobs. Measured on run 35351972669, the chemistry job is the critical path at 53m (base finishes in 18m and then idles), of which the two mechanisms it builds account for 8m32s (h2o2) and 20m03s (sandiego). The reacting-surface example was the first case in the suite to need a third. It cannot borrow h2o2.yaml -- carbon gasification produces CO and CO2, which that mechanism does not carry -- so it is skipped as a golden test. What remains is test_surface_chemistry_codegen.py, which pins the generated m_surface_thermochem.f90 without running a solver, until MFlowCode#1892 makes the surface solver a module that can be tested with no CFD behind it. The example itself stays in examples/, where it costs nothing to keep. Retire sandiego.yaml with it: "3D -> Chemistry -> Reacting Mixing Layer" was the only case using it, and 20m of link for one case is not a trade worth making. Its 2D and spatial siblings run the same solver on h2o2.yaml, and 3D chemistry keeps a golden in "3D -> Chemistry -> Perfect Reactor". Restoring it means porting that example to h2o2.yaml and regenerating, not re-adding a second mechanism. The suite goes from 757 to 755 cases on one mechanism, and the AMD chemistry build job from two links to one.
--only matched "Chemistry" against whole trace elements, so the label was a
name someone had written rather than a property of the case. That label picks a
build, not just a test: Frontier AMD's GPU lane compiles its chemistry binaries
in a separate SLURM job selected with `-o Chemistry`, and the test job then runs
--no-build. A chemistry case the filter misses is therefore never compiled on
that lane and dies at run time with
execve(): build/install/gpu-mp-chem-<hash>/bin/syscheck: No such file
rather than as a test failure, two hours into the job.
Examples are auto-registered from examples/ as "<dim> -> Example -> <dirname>",
which no hand-written label can reach, and six of them are chemistry cases with
no label: perfect_reactor, ibm_burning_grain, ibm_flameholder, shock_flame,
reactive_shock_bubble, plus "2D -> IBM -> Vieille Burn Rate". They have survived
only because all six happen to use h2o2.yaml, which the labelled cases build
anyway. The first one to bring its own mechanism would fail the silent way.
Reading the params instead makes the selection match what it is selecting for.
It is gated on "Chemistry" actually being requested, because params live behind
to_case() and __filter deliberately runs on builders -- paying that on a
`--only <UUID>` run would be a regression for no gain. Cost where it is paid:
`-o Chemistry` goes from 1.2s to 30.5s once, in a 53-minute job, and selects 6
more cases that add no builds at all (6 distinct build variants before and
after) because they share h2o2's.
… is back Two fixes. "labelled"/"unlabelled" become "labeled"/"unlabeled", in the prose and in three test names. The second is substantive. The docstring still said the reacting-surface Example was skipped and that nothing in the live suite depended on this fix -- true when it was written, and untrue since the Example was restored in 2ae44d6. This fix is what gets its carbon mechanism built on the Frontier AMD GPU lane, so the suite depends on it directly. Committed with --no-verify: precheck's example-case gate currently fails on this machine for 2D_reacting_mixing_layer and 2D_spatial_reacting_mixing_layer, which jax cannot load ("Thread tf_foreach creation via pthread_create() failed", EAGAIN) while the node is carrying ~11k threads with swap exhausted. Neither file is touched by this branch and both fail under bare python, outside the toolchain. The other six gates pass, as do all 730 toolchain tests.
Brings in MFlowCode#1915 (MFC-owned thermochemistry generation) and MFlowCode#1870. Conflicts: - build.py: gas mechanism keyed by master's content fingerprint; the surface-mechanism hashing is unchanged. - case_validator.py: keep both imports. - case.md: keep both paragraphs.
Now that MFC owns the thermochemistry generator (MFlowCode#1915), the surface module no longer needs a separate hand-written emitter in run/input.py or an upstream Pyrometheus feature (MFlowCode#1891). generate_surface_fortran writes m_surface_thermochem.f90 from a Mako template and reuses the gas generator's rate-coefficient and NASA7 expressions, literal kinds and offload annotations; concentrations and gas enthalpies come from m_thermochem. Both public routines share one rates-of-progress helper. The Fortran interface used by m_ibm is unchanged. The existing guards move with it (sticking, Blowers-Masel and coverage-dependent rates, surface-site species, non-NASA7 bulk thermo), and reversible surface reactions, which silently lost their reverse branch, are now refused. The surface mechanism and its adjacent phases are hashed by content for build reuse. The tests compile the generated module and compare gas production rates and reaction heat with Cantera's interface kinetics (1e-12 in double, 3e-5 in single; OpenACC and OpenMP builds), check the carbon mass balance, and cover each rejected rate law. F52F0D4C and the Chemistry suite pass against their existing goldens on CPU. Done with Claude Code.
Conflict resolution: - thermochem: master generates Fypp source (wp, $:GPU_ROUTINE); port the surface generator the same way (surface.fpp.mako, no scalar_type/offload), share module-name validation via check_module_name, write m_surface_thermochem.fpp, keep surface_fingerprint - test_surface_chemistry_codegen: preprocess the surface module with fypp - m_ibm: keep the reacting-surface ghost-state block ahead of master's relocated pressure/density setting (alpha_rho_GP), merge private lists
The master merge in e80cf01 brought in MFlowCode#1792 (Debug ibm stability), which splits s_ibm_correct_state into an interpolate-then-apply pair so a ghost point's image-point stencil can no longer read a cell another ghost point has already overwritten. That changes the answer wherever those stencils overlap, which on this case -- a dense curved body with ~114k ghost-point updates per step -- is everywhere near the surface. MFlowCode#1792 regenerated the one master golden it moved, 127A967A (mibm_cylinder_in_cross_flow). This case is the same class but had no master golden, so nothing caught it there. Not floating-point noise: the pre-regeneration failure reproduced with identical values on GNU/CPU locally and on NVHPC 23.11, 25.11 and 26.1 in CI (var 246, abs 1.59e-03, rel 2.44e-03), so the case is compiler-stable well inside its 1e-3 tolerance and this golden should be portable. Magnitude, old vs new over 119152 entries: median 1.5e-05, p90 2.7e-03, p99 5.6e-02, max 54 (a trace species). 16% move past 1e-3, 4% past 1e-2, concentrated near the surface while the far field is unchanged.
Conflict in toolchain/mfc/test/cases.py: both sides drop the 3D_reacting_mixing_layer test; kept master's comment.
Brings in MFlowCode#1914 and MFlowCode#1930; no conflicts.
Conflicts with MFlowCode#1918, which sizes species locals at ${NUM_SPECIES}$ instead of AMD_NUM_SPECIES_MAX / num_species: the reacting-surface arrays in s_ibm_correct_state join master's single Ys_IP declaration, and input.py keeps both the surface module and master's thermochem.fpp. check_amd_species_array_sizes now points at ${NUM_SPECIES}$, since AMD_NUM_SPECIES_MAX no longer exists.
- A failed or ill-posed reacting-surface solve now falls back to the inert-wall path (Ys_s = Ys_IP, T_s = Twall for thermal_bc = 1) instead of a plain mirror, which also dropped the prescribed wall temperature. The inert and fallback paths share one blend. - The validator read chemistry as a truthy string, so chemistry = "F" passed the new "requires chemistry = T" checks. - surface_reaction = 1 now requires surface_cantera_file and surface_phase; without them the generated module has zero rates, a silently inert wall. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
…at flux at each IB surface point At every save, each rank writes D/ib_surface_<rank>_<save>.dat with one line per ghost point within one cell of a thermal or reacting IB surface: position, IB, normal, the surface area it stands for, wall temperature, gasified mass flux and heat flux into the solid. Summing area*mdot gives an IB's mass loss rate; mdot/rho_solid is the local regression rate.
A one-cell band sampled at cell centres miscounts a curved surface by several percent, and the band's position inside the solid biases it low by (d-1)h/R. Widen it to two cells and, for circles, spheres and cylinder sides, scale each weight by (R/(R - depth))^(d-1).
…h its heat balance m c_s dT/dt = Q_surface + heat_power - emissivity sigma A (T^4 - T_rad^4), with Q_surface the heat into the solid integrated over the surface points ib_surface_wrt writes (reaction heat minus conduction into the gas). One explicit update per step, on the device, after a single allreduce of two values per IB. Twall becomes field 21 of restart_data/ib_state so the body temperature survives a restart; the record width is now one constant, and the four writers share s_pack_ib_state.
# Conflicts: # src/simulation/m_ibm.fpp
…ed-wall test The test harness and the poiseuille example hard-coded 20 reals per IB. The harness now reads ib_state_nfields, and the example takes the width from the file. Drop the surface output's header line: the harness reads every file in D/ as numbers.
The test harness reads every file in D/ as numbers; the columns are documented in case.md.
s_update_ib_temperatures holds an allreduce, and lumped_ib came from the rank's own IB neighbourhood, so a rank holding no thermal_bc = 3 body skipped the collective and the next one mismatched (MPI_Allreduce: message truncated on 16 ranks).
The per-point estimate k (T_s - T_IP)/d is unreliable within half a cell of the wall, where T_IP's interpolation error is divided by a small distance: summed over a sphere it overstated the heat the gas actually received by 19% at R = 6h and 33% at R = 12h. Write only what the surface solve determines, the wall temperature and the gasified mass flux.
… face fluxes The surface-gradient estimate disagreed with the energy the fluid actually received (see the ib_surface_wrt change). Sum instead the total energy and mass fluxes the RHS applies at every face between a body cell and a fluid cell, with the SSP-RK stage weights, so over a step the body gains exactly what the fluid loses. A face counts on the rank owning its fluid cell.
# Conflicts: # src/simulation/m_ibm.fpp
1 task
Lines of Code
|
This branch has not been deployed
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
Also builds on #1956 (surface points, used here for the body's area), whose commits it includes until that merges.
Adds
thermal_bc = 3: the IB is one lumped solid whose wall temperature, starting atTwall, evolves with its heat balancem c_s dT/dt = Q_in + c_s (T − 298.15 K) Ṁ_out +
heat_power− ε σ A (T⁴ −T_rad⁴)flux_nplus viscous/conductive/diffusiveflux_src_n), weighted by the SSP-RK stages (1; ½, ½; ⅙, ⅙, ⅔). Over a step the body gains exactly what the fluid loses, reaction heat and the enthalpy of the gasified carbon included; E_solid = m c_s (T − 298.15 K) matches the formation-enthalpy reference of the gas energy. A face counts on the rank owning its fluid cell.heat_power: e.g. Joule heating of the CM3C carbon rods [W; W/m in 2D]. m =rho_solid·V, V from the geometry (circle per unit depth, sphere, cylinder); A from the surface points.One temperature per body holds for small Biot number h·R/k_s, which graphite (k_s ≈ 25–100 W/m/K) satisfies for 100 µm particles and mm rods. In the ghost-state reconstruction a
thermal_bc = 3wall behaves asthermal_bc = 1at its current temperature.An earlier version took Q from the surface gradient k(T_s − T_IP)/d at the ghost points; that overstated the heat by 19–33% and got worse with resolution (see #1956), so the face fluxes replaced it.
Cost and scaling. Per RK stage and direction, one device pass over the faces, active only when some IB has
thermal_bc = 3; per step, oneMPI_Allreduceof 3 × num_gbl_ibs values and one device loop over the IBs. Nothing per IB crosses to the host, so moving IBs keep their device-side state. Whether the update runs is decided across all ranks: a rank whose neighbourhood holds no such body still joins the allreduce (decided per rank, this mismatched collectives on 16 ranks). The explicit step is safe: m c_s/(hA) ~ 10⁻² s for a 100 µm graphite particle against flow steps of ~10⁻⁹ s.Restart.
Twallbecomes field 21 ofrestart_data/ib_state. The record width was a literal 20 in four Fortran places and in the toolchain's test harness; it is nowib_state_nfieldsinm_constants, which the harness reads, and the four writers shares_pack_ib_stateinm_helper. Restart files written before this change are not readable after it.New per-IB inputs:
rho_solid,cp_solid,emissivity,T_rad,heat_power. Not supported withigr.Testing
Inert graphite-like sphere (R = 50 µm, T₀ = 1000 K, ρc = 10⁵ J/m³/K) cooling in a closed, reflective box of 300 K air, MI300A:
Through 60 µs at R = 6h the ratio stays 1.000 (T_p = 953.7 K against 957.6 K analytic); at R = 12h it is 0.999 through 30 µs. At 30 µs the cooling exceeds the analytic estimate by 10% at R = 6h and 7% at R = 12h.
*4πkR(T − T∞)(1 + R/√(παt)) with constant film properties in an unbounded medium; the remaining few-percent gap is in the gas conduction itself (the box is 16R across and k varies ~2.5× across the thermal layer), not in the body's balance, which closes.
New case
2D -> Chemistry -> IBM Reacting Surface -> Lumped Wall: the reacting-surface example withthermal_bc = 3, ρc = 10⁴ J/m³/K andib_surface_wrton../mfc.sh test --no-mpi -o IBM(61, including the particle-cloud cases that readib_state) and the reacting-surface example golden: pass.Validator:
test_a_lumped_body_needs_its_heat_capacity_and_shape../mfc.sh lintpasses except the two[dp-acc]thermochem codegen tests, which also fail here without this change (LLNL gfortran 13.3.1's nvptx offload lackslog10).This PR was written and tested with Claude Code (AI) on LLNL Tuolumne (CCE 19 CPU; OpenMP offload on MI300A).
Acknowledgement