Repository navigation
Add patch_ib%thermal_layer: start a hot IB in its conduction solution (after #1821) - #1955
Draft
sbryngelson wants to merge 44 commits into
Draft
sbryngelson wants to merge 44 commits into
sbryngelson wants to merge 44 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.
Reverts the surface half of c23572f. Skipping that Example left nothing in CI exercising the surface boundary condition: no remaining case set surface_cantera_file, surface_phase or surface_reaction, so ~559 lines of m_ibm.fpp shipped with only the codegen unit tests behind them, and those check the generated Fortran without ever running a solver. The branch has six commits fixing compiler-specific GPU offload bugs in exactly that code, which is the wrong place to be running blind. It was also an unnecessary trade. The constraint on the Frontier AMD GPU lane is that the chemistry build job fits its 1h59m walltime while being the critical path, not a literal count of mechanisms -- and retiring sandiego.yaml freed more of that budget than the carbon mechanism needs. Measured on run 35351972669: sandiego's simulation link was 20m03s and h2o2's 8m32s, of a 53m job. Carbon is h2o2's size class (11 species / 33 reactions against 10 / 29), so h2o2 + carbon should come in under the h2o2 + sandiego pair it replaces. Net against the tip of this branch: the suite keeps two gas mechanisms, the same count as master, and trades a 3D mixing-layer golden that its 2D and spatial siblings already cover for the only end-to-end test of the feature this branch exists to add.
Conflicts: m_particle_cloud.fpp - master moved cloud generation to pre_process, where a cloud IB now carries only position, kinematics and radius. The thermal_bc/Twall/surface_reaction defaults this branch set there move to s_assign_particle_cloud_ib_defaults in simulation/m_start_up.fpp, which is where master now fills the rest of a cloud patch. m_ibm.fpp - master extracted s_compute_ghost_point_pressure/_velocity out of s_ibm_correct_state, taking v_blow and the slip/rotation handling with them. The reacting-surface block stays in the loop; it keeps its own norm/buf for the Stefan-flow superposition, and the GPU private list is master's plus this branch's surface variables and convergence reductions. lint_test_suite.py - master's new gate asked whether a case's trace carries a "Chemistry" segment. This branch had already made the --only Chemistry label derive from the params (an auto-registered Example can never say "Chemistry" in its trace), so the gate now asks case_filter_labels rather than re-deriving the rule, and the ibm_reacting_surface golden is covered by the chem pre-build. test_case_validator.py - both sides appended a test class; both kept.
… 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>
A thermal_bc = 1 circle or sphere otherwise meets ambient gas as a temperature step at t = 0. thermal_layer = delta > 0 initializes the surrounding gas in the conduction solution for a wall brought to Twall a time delta^2/(4 alpha) earlier, T = T_inf + (Twall - T_inf) (R/r)^((d-1)/2) erfc((r - R)/delta), at fixed pressure and composition (exact for a sphere, leading order for a cylinder). Off by default.
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
Adds
patch_ib(i)%thermal_layer= δ [m] (default 0, off) forthermal_bc = 1circles (geometry 2) and spheres (geometry 8).A hot isothermal IB otherwise meets the ambient gas as a temperature step at t = 0, and the mirrored ghost state doubles that step (T_ghost = 2 T_wall − T_gas ≈ 3300 K for an 1800 K wall in 300 K air). On fine cells the gas next to the wall heats faster than it can expand: in a 3D bed of 191 carbon spheres (d = 100 µm, 12 cells/d), near-wall cells reached 1.87 bar within 4 steps and the run went NaN on step 5. That failure did not depend on dt (CFL 0.5, 0.15 and 0.05 all failed at t ≈ 3–6e-8 s) or on surface chemistry (identical with
surface_reaction = 0).With
thermal_layer = δ, pre_process starts the gas in the conduction solution for a wall impulsively brought toTwalla time t₀ = δ²/(4α) earlier:T = T∞ + (T_wall − T∞) (R/r)^((d−1)/2) erfc((r − R)/δ)
at fixed pressure and composition, so ρ scales as T∞/T. This is exact for a sphere and the leading-order (short-time) term for a cylinder (Carslaw & Jaeger). Overlapping layers take the IB with the largest temperature change; periodic directions use the nearest image. It acts only on the initial condition.
Testing
2D -> Chemistry -> IBM Reacting Surface -> Thermal Layer(A74B7381): the reacting-surface example at the Example grid cap with δ = D/2. The initial density goes from 0.29 kg/m³ (1200 K) at the wall to the 1.18 kg/m³ freestream.thermal_layermust be ≥ 0 and, when > 0, needsthermal_bc = 1on geometry 2 or 8.Twall= 1800 K: without the layer it aborts at step 4; with δ = 0.2 d (2.4 cells) it runs 835 steps to t·U/d = 5 cleanly. At 24 cells/d the bed also runs without the layer../mfc.sh lintpasses except the two[dp-acc]thermochem codegen tests, which also fail without this change here (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