fix(photochem): exclude flux from VODE convergence; hard-code atol for DTypeFront; add docs by chongchonghe · Pull Request #1956 · quokka-astro/quokka

chongchonghe · 2026-06-15T12:24:50Z

Description

Add hard-coded VODE absolute tolerances for the DTypeFront photoionization problem, and exclude radiation flux from VODE convergence checks. docs/markdown/photoionization.md now has complete documentation for the photoionization module: governing equations (M1 radiative transfer, hydrogen thermochemistry, on-the-spot approximation), numerical scheme (operator splitting with VODE), and the full VODE tolerance design rationale. docs/markdown/parameters.md has the complete integrator.* parameter reference.

The SetAtolFromPhysics machinery (automatic atol derivation from physical scales) has been moved to PR #1980.

Summary of changes

File	Change
`src/problems/DTypeFront/testDTypeFront.cpp`	`Erad_floor = a_rad * (0.01 K)^4` with documentation of the relationship to `atol_rad_num`
`inputs/DTypeFront.toml`	Hard-coded `integrator.atol_*` values with documentation
`src/networks/photoionization/actual_rhs.H`	Clamp `N_gamma >= 0` in RHS and Jacobian to prevent sign reversal when photons are fully absorbed
`docs/markdown/photoionization.md`	New: full photoionization module documentation — governing equations, numerical scheme, VODE tolerance design
`docs/markdown/parameters.md`	New Integrator section documenting all `integrator.*` parameters
`extern/Microphysics/.../vode_dvode.H`	Zero flux error weight — excludes it from all convergence norms
`extern/Microphysics/.../vode_dvnlsd.H`	Skip flux indices in corrector DEL and accumulated error ACNRM
`extern/Microphysics/.../vode_dvstep.H`	Remove flux change-factor step rejection check
`extern/Microphysics/.../vode_dvhin.H`	Do not constrain initial step by flux tolerance

Key design decisions

Hard-coded tolerances for DTypeFront — atol_spec = 1.4e-3 cm^-3, atol_enuc = 1.24e8 erg/g, atol_rad_num = 5e-2 cm^-3, radiation_failure_tolerance = 0.5 cm^-3. atol_spec and atol_enuc correspond to a negligibility floor of 1e-5 × n_H and a 1 K temperature accuracy at n_H = 139.9 cm^-3. atol_rad_num = 5e-2 is the performance plateau for this problem (see Performance below).
atol_rad_num and radiation_failure_tolerance are coupled — radiation_failure_tolerance must exceed the maximum negative excursion that the BDF predictor-corrector produces at the ionization front before converging (~0.05–0.2 cm^-3 for DTypeFront). Setting radiation_failure_tolerance = 0.5 provides a safe margin. Note: Microphysics commit 47696e49 introduced a 1.5× final-state factor for species (vode_final_state_species_failure_tolerance_factor), which allows species_failure_tolerance = atol_spec without manual inflation — but no equivalent factor exists for radiation, so radiation_failure_tolerance requires its own margin.
Erad_floor is separate from atol_rad_num — the former is a compile-time constexpr in RadSystem_Traits controlling the M1 hyperbolic floor; the latter is the VODE tolerance. They must satisfy N_gamma_floor << atol_rad_num so dark cells return in one BDF step.
Flux is excluded from VODE convergence checks — flux is a passive scalar (one-way coupled: dF/dt depends on n_HI but flux does not feed back into species, energy, or N_gamma). N_gamma has an additional isotropic source term (stellar emission) that flux lacks, so flux is not analytically determined by N_gamma — but it remains a diagnostic quantity whose convergence should not drive the timestep.

Performance

All benchmarks: ROCm GPU (64^3, AMD GPU on moth.anu.edu.au).

vs. development branch (64^3, stop_time = 2.895e12 ≈ 3 t_s):

Branch	Total (s)	PhotoChem (s)	PhotoChem %
`development`	439.8	436.1	99.1 %
this PR	14.2	10.5	73.9 %
speedup	31.0×	41.6×

Effect of atol_rad_num tuning within this PR (64^3, stop_time = 1.93e13 = 20 t_s):

atol_rad_num was swept from 1.25e-6 to 0.2 cm^-3. Performance plateaus at ~5e-2 cm^-3: below this VODE spends extra steps in dark cells; above it the step-rejection loop (N_gamma overshooting near −radiation_failure_tolerance at the ionization front) erodes the gain.

`atol_rad_num`	Total (s)	PhotoChem (s)	PhotoChem %
1.25e-6 (initial)	218.7	194.3	88.8 %
5e-2 (final)	98.3	73.9	75.2 %
speedup	2.22×	2.63×

Related issues

N/A

Checklist

I have added a description (see above).
I have added a link to any related issues (if applicable; see above).
I have read the Contributing Guide.
I have added tests for any new physics that this PR adds to the code.
(For quokka-astro org members) I have manually triggered the GPU tests with the magic comment /azp run.

Replace 1e-10 atol defaults with values grounded in the physical scales of the problem: atol_spec = 1e-3 cm^-3 (1e-4 × n_H negligibility floor) and atol_rad_num = 1e-2 cm^-3 (below the photoionization-rate significance threshold of ~0.8 cm^-3). These prevent the solver from pursuing spurious precision on trace species and negligible photon densities. All three physics checks continue to pass. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

atol_rad_num = 1e-2 allowed VODE to step the photon number density into negative values (N_gamma starts at ~1e-89 cm^-3, far below the tolerance floor), reversing the photoionization term sign and causing Newton-iteration failures. Reduce to 1e-5, which is safely below the empirical failure threshold (~3e-4) while remaining five decades more relaxed than the original default of 1e-10. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

- photochemistry.hpp: pass N_gamma_eff = max(Erad - Erad_floor, 0) / E_photon to VODE so dark cells (with only floor radiation) see N_gamma = 0 and take one trivial step; add Erad_floor_ back in write-back for energy conservation - actual_rhs.H: clamp n_gamma = max(rn[0], 0) in RHS and Jacobian to prevent sign reversal if VODE steps N_gamma slightly below zero within tolerance - testDTypeFront.cpp: set Erad_floor = a_rad * (10 K)^4 as a physically motivated radiation floor (vs the former 1e-99) - DTypeFront.toml: relax atol_rad_num from 1e-5 to 1e-2 (floor subtraction eliminates the sign-reversal numerical stability constraint); set radiation_failure_tolerance = 0.1 so VODE dipping N_gamma to -atol_rad_num does not abort the burn Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…floor subtraction Erad_floor = a_rad * (1 K)^4 gives N_gamma_floor ~ 3.5e-4 cm^-3, which is: - below atol_rad_num = 0.01, so VODE returns in one step (absolute control) - small enough that cumulative ionization over 10 steps is ~3.5e-3 cm^-3 per cell, leaving x_HI > 99.99% in dark cells (neutral temperature check passes) - physically motivated vs the arbitrary 1e-99 Floor subtraction in photochemistry.hpp is reverted — not needed when the floor is chosen correctly. The Erad_floor value itself ensures dark cells see negligible N_gamma without any special handling in the burner. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

for more information, see https://pre-commit.ci

gemini-code-assist

Code Review

This pull request adjusts resolution and integrator tolerance parameters in DTypeFront.toml, clamps n_gamma to zero in actual_rhs.H to prevent negative values from VODE steps, and updates Erad_floor in testDTypeFront.cpp. The review feedback suggests ensuring mathematical consistency in the Jacobian by updating derivatives with respect to state.rn[0] to account for the clamped n_gamma, and updating the hardcoded initial radiation energy density in testDTypeFront.cpp to use the new Erad_floor value.

Important

The consumer version of Gemini Code Assist on GitHub is being sunset. Starting June 18, 2026, new organization installations will be blocked, and all code review activity will officially cease on July 17, 2026.
For more details on the timeline and next steps, please review the Help Documentation.

Erad_floor = a_rad*(0.1K)^4 gives N_gamma_floor = 3.5e-8 cm^-3, six orders of magnitude below atol_rad_num = 0.01, so VODE returns in one BDF step in dark cells regardless of the number of timesteps. The equilibrium ionisation from floor photons is x_HII_eq ~ 0.002%, well below the 0.01% neutral-temperature test threshold, making the strategy robust for arbitrarily long (e.g. 1M step) production runs without needing floor subtraction in the photochemistry solver. Remove the 10-step hard cap; stop_time now governs (~15 CFL steps). Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…ical parameters Users now specify integrator.typical_n_H, integrator.desired_accuracy_on_T_at_typical_n_H, and integrator.typical_minimal_radiation_T instead of hand-tuning integrator.atol_spec / atol_enuc / atol_rad_num / atol_rad_flux. SetAtolFromPhysics<problem_t>() injects derived values into ParmParse before init_extern_parameters() reads them, leaving the Microphysics submodule untouched. The new keys and existing integrator.atol_* keys are mutually exclusive. Neither present is a backward-compatible no-op. Erad_floor (the M1 radiation floor) remains a compile-time constexpr in each problem's RadSystem_Traits; typical_minimal_radiation_T only sets the VODE tolerance, not Erad_floor. DTypeFront.toml updated to demonstrate the interface. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

…d_floor to 0.01K radiation_failure_tolerance is a physical guard (max allowed negative photon count per cm^3), not a numerical tolerance — defaulting to 0.01 cm^-3 avoids false burn failures in near-vacuum cells regardless of how tight atol_rad_num becomes. Erad_floor lowered from a_rad*(0.1K)^4 to a_rad*(0.01K)^4 so that typical_minimal_radiation_T=10K still gives atol_rad_num/N_gamma_floor ~ 1e4, ensuring dark cells converge in one BDF step. Add docs/markdown/photoionization.md documenting atol/rtol parameter choices. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

VODE's built-in defaults (~1e-10) are unusably tight for photochemistry and will cause the integrator to stall. Require either typical_* (recommended) or atol_* parameters to be present. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

E_photon cancels, so the user only needs T_min and T_floor: 1e-6 × (T_min / T_floor)⁴ ≥ 10⁴ ↔ T_floor ≤ T_min / 316 Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

- Add Physical constants section defining a_rad, k_B, m_p, c_v, E_photon - Introduce T_min, T_floor, N_gamma_floor with explicit definitions - Fix incorrect ratio in Relationship example (was 10^4, now 10^6) - Add missing rtol_rad_flux to the rtol description - Make step 4 actionable with "if this fails" guidance - Tighten prose throughout Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

for more information, see https://pre-commit.ci

Co-Authored-By: Claude <noreply@anthropic.com>

Temporary change for GPU validation — revert before merging. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

chongchonghe · 2026-06-15T16:16:59Z

/azp run

chongchonghe · 2026-06-22T17:27:26Z

@BenWibking Thanks for the detailed review! I've addressed your comments:

Microphysics (new branch chong/claude/radiation-failure-tolerance-factor on fork):

Added vode_final_state_radiation_failure_tolerance_factor = 1.5 in vode_type.H, mirroring the species_failure_tolerance design — internal VODE step checks use radiation_failure_tolerance directly; the final interpolated state check uses 1.5x to account for non-monotonic interpolation.
Updated actual_integrator.H, actual_integrator_sdc.H, and both integrator_cleanup functions.

parameters.md:

Replaced typical_* sections with a single tolerance parameter table (the SetAtolFromPhysics machinery has moved to feat(photochem): add SetAtolFromPhysics to derive VODE atol from physical parameters #1980).
Updated radiation_failure_tolerance description: it's a numerical tolerance at internal nodes, should equal atol_rad_num; final state relaxed to 1.5x.
Added notes that integrate_energy, scale_system, and X_reject_buffer don't work with number densities.
Noted subtract_internal_energy and use_number_densities must be set correctly for Quokka.

photoionization.md:

Added citations [@Krumholz2007], [@Koyama2002], [@Osterbrock1989].
Removed all SetAtolFromPhysics / typical_* content (moved to feat(photochem): add SetAtolFromPhysics to derive VODE atol from physical parameters #1980).
Clarified radiation_failure_tolerance and species_failure_tolerance design.

Could you take another look when you have a moment?

Generated by Claude Code

chongchonghe · 2026-06-22T17:29:28Z

/azp run quick

chongchonghe · 2026-06-22T17:29:35Z

/azp run rocm-quick

azure-pipelines · 2026-06-22T17:29:37Z

Azure Pipelines could not run because the pipeline triggers exclude this branch/path.

azure-pipelines · 2026-06-22T17:29:44Z

Azure Pipelines could not run because the pipeline triggers exclude this branch/path.

…view atol_rad_num = 1.0e-2 cm^-3 matches atol_spec = 1.0e-2 cm^-3, making the photon number density tolerance comparable to the electron number density tolerance. radiation_failure_tolerance = 0.1 (10x atol_rad_num) follows the same pattern as DTypeFront. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

Per Ben's review: both tolerances now describe primary (internal substep rejection) and secondary (final interpolated state at 1.5x) uses. Removed "physical guard" / "maximum allowed" language in favor of precise VODE mechanism descriptions. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

chongchonghe · 2026-06-23T01:08:40Z

/azp run quick

azure-pipelines · 2026-06-23T01:08:53Z

Azure Pipelines successfully started running 1 pipeline(s).

chongchonghe · 2026-06-23T01:12:45Z

/azp run rocm-quick

azure-pipelines · 2026-06-23T01:12:55Z

Azure Pipelines successfully started running 1 pipeline(s).

BenWibking

I think your description of the cooling is still grammatically confusing. I have suggested two rewordings below.

chongchonghe · 2026-06-23T02:09:27Z

OK, comments addressed. Ready for review again. @BenWibking

BenWibking

LGTM.

chongchonghe and others added 8 commits June 13, 2026 15:44

setup toml for testing

59cad2e

use ben's integrator params

9376b3f

nx=8

a50228c

[pre-commit.ci] auto fixes from pre-commit.com hooks

6a2019a

for more information, see https://pre-commit.ci

gemini-code-assist Bot reviewed Jun 15, 2026

View reviewed changes

Comment thread src/networks/photoionization/actual_rhs.H

Comment thread src/problems/DTypeFront/testDTypeFront.cpp Outdated

chongchonghe changed the title ~~perf(photochem): set physically grounded VODE tolerances for DTypeFront (7.8× speedup)~~ perf(photochem): set physically grounded VODE tolerances for DTypeFront (8.9× overhead reduction) Jun 15, 2026

chongchonghe and others added 10 commits June 16, 2026 00:28

update plan

7ad29fe

docs: show derived atol values as comments in DTypeFront.toml

e07527b

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

docs: document integrator parameters in parameters.md

eb18f0b

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

docs: clarify when typical_* parameters are required in parameters.md

3168f80

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

docs: simplify atol/N_floor check to temperature-only formula

1e64722

E_photon cancels, so the user only needs T_min and T_floor: 1e-6 × (T_min / T_floor)⁴ ≥ 10⁴ ↔ T_floor ≤ T_min / 316 Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

docs: replace "1 ppm" with "one millionth" for clarity

52a23a2

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

chongchonghe changed the title ~~perf(photochem): set physically grounded VODE tolerances for DTypeFront (8.9× overhead reduction)~~ feat(photochem): derive VODE atol from physical parameters via SetAtolFromPhysics Jun 15, 2026

pre-commit-ci Bot and others added 7 commits June 15, 2026 16:00

[pre-commit.ci] auto fixes from pre-commit.com hooks

2608401

for more information, see https://pre-commit.ci

fix(docs): add photoionization.md to mdBook SUMMARY.md

3ae4d47

Co-Authored-By: Claude <noreply@anthropic.com>

cleanup

2faacb8

revert nx

84f5ac9

ci: scope AMD GPU pipeline to DTypeFront only

7eee4c6

Temporary change for GPU validation — revert before merging. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

ci: scope CUDA pipeline to DTypeFront only

6045d2c

Temporary change for GPU validation — revert before merging. Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

doc

1325b15

Merge branch 'development' into chong/photo-opt2

ebf8002

BenWibking requested changes Jun 22, 2026

View reviewed changes

Comment thread docs/markdown/photoionization.md Outdated

Comment thread docs/markdown/photoionization.md

Comment thread docs/markdown/photoionization.md Outdated

Comment thread docs/markdown/photoionization.md Outdated

Comment thread docs/markdown/photoionization.md Outdated

chongchonghe and others added 3 commits June 23, 2026 10:34

kurmholz :

d7a2c07

Merge branch 'development' into chong/photo-opt2

ad12f38

Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>

BenWibking requested changes Jun 23, 2026

View reviewed changes

Comment thread docs/markdown/photoionization.md Outdated

Comment thread docs/markdown/photoionization.md Outdated

rewrite

cef8805

BenWibking previously approved these changes Jun 23, 2026

View reviewed changes

dosubot Bot added the lgtm This PR has been approved by a maintainer label Jun 23, 2026

chongchonghe enabled auto-merge June 23, 2026 02:49

chongchonghe mentioned this pull request Jun 23, 2026

Add v/c work term in photoionization #1988

Open

revert CI

ca45302

chongchonghe dismissed BenWibking’s stale review via ca45302 June 23, 2026 03:20

chongchonghe requested a review from BenWibking June 23, 2026 03:20

BenWibking approved these changes Jun 23, 2026

View reviewed changes

chongchonghe added this pull request to the merge queue Jun 23, 2026

Merged via the queue into development with commit b47bb38 Jun 23, 2026
51 checks passed

chongchonghe deleted the chong/photo-opt2 branch June 24, 2026 01:23

Uh oh!

Conversation

chongchonghe commented Jun 15, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Description

Summary of changes

Key design decisions

Performance

Related issues

Checklist

Uh oh!

gemini-code-assist Bot left a comment

Choose a reason for hiding this comment

Code Review

Uh oh!

Uh oh!

Uh oh!

chongchonghe commented Jun 15, 2026

Uh oh!

chongchonghe commented Jun 22, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

chongchonghe commented Jun 22, 2026

Uh oh!

chongchonghe commented Jun 22, 2026

Uh oh!

azure-pipelines Bot commented Jun 22, 2026

Uh oh!

azure-pipelines Bot commented Jun 22, 2026

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

chongchonghe commented Jun 23, 2026

Uh oh!

azure-pipelines Bot commented Jun 23, 2026

Uh oh!

chongchonghe commented Jun 23, 2026

Uh oh!

azure-pipelines Bot commented Jun 23, 2026

Uh oh!

BenWibking left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

chongchonghe commented Jun 23, 2026

Uh oh!

BenWibking left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

chongchonghe commented Jun 15, 2026 •

edited

Loading

chongchonghe commented Jun 22, 2026 •

edited

Loading