Update of rates and addition of Photchem class - #48
Draft
jeonggyukim wants to merge 112 commits into
Draft
Conversation
photchem-full was opened 2025-02-25 (PR #48) and had not been updated since `Add PhotChem class` (58c6405). master moved 131 commits forward in the interim, including: - numpy 2 migration (numpy>=2, yt>=4.4) + removed-API fixes - SB99 dataset refactor + Zari+22 SFH support - ISRF API rewrite - find_files.py / load_sim.py docstring + behavior cleanups - README rewrites; conda/miniforge install docs - Various docs / autodoc improvements Dry-run on pyathena-ct-mvp showed zero conflicts (photchem-full's new files under data/microphysics/ + pyathena/microphysics/ touch disjoint paths from master's changes). Confirmed by re-running this merge with --no-commit -- same clean auto-merge here. The sibling sub-PR #67 (pyathena-ct-mvp) will rebase on this state before merging.
Static method `get_ct_ion_HII_OI_Draine11` was misnamed `get_ct_ion_OI_HII_Draine11` at the dispatch site (line 183), causing `AttributeError` when O I charge-exchange ionization rate is requested. Symptom: `get_ct_ion_rate(8, 8, T)` crashes. Caught by new charge-exchange rate test `test_ct_rate_balance.py::test_ct_ion_rate_finite_nonneg[8-8-O I]`.
Pin the current behavior of `ChargeTransferRate.get_ct_rec_rate` and `get_ct_ion_rate` before downstream code wires charge-exchange into the abundance sweep in `evolve_one_species`. Tests document: - API conventions (which (Z, N) labels reactant vs product per method); - Sign of the Boltzmann factor in `get_ct_ion_rate` via the N I endothermic pair (ct_ion must grow with T faster than ct_rec); - Near-resonance behavior for O+H (similar T-growth, Draine 2011 per-J-level detailed balance at high T); - Order-of-magnitude reference at T = 1e4 K; - Known-quirk pinning: S I CT-ion is a constant ~1e-14 placeholder in the Cloudy data table (`a=1e-5, b=c=d=dE=0`) -- documented so future rate upgrades flag this site. - High-q fallback behavior (q > 3 in `get_ct_ion_rate` returns 0; in `get_ct_rec_rate` returns Dalgarno generic estimate `1.92e-9*(q+1)`). These tests are the regression anchor for the charge-exchange completeness fixes described in `tigris-notes/docs-claude/tigris-ncr/pyathena_ct_fixes_plan.md` (MVP items 1-3); the wiring of CT into the abundance sweep follows in later commits on this branch. Pytest infrastructure (empty `tests/__init__.py` + this test file) is new to pyathena -- run with `pytest pyathena/microphysics/tests/`.
Relocate `pyathena/microphysics/tests/` -> `tests/microphysics/` so
new test suites (io, fields, obs, ...) can drop into sibling subdirs
without nesting inside the package. Adds:
- `tests/README.md` documenting layout, run commands, and the
convention for adding new tests + diagnostic plots.
- `tests/conftest.py` providing:
* `figures_dir` fixture -> `tests/figures/` (created on demand).
* `--no-figures` pytest CLI option to skip plot generation.
* `save_figures` fixture (True unless `--no-figures` passed).
* `ion_color(Z, q)` helper + `ion_colors` fixture matching the
per-ion color convention in
`pyathena/microphysics/photchem.py:_set_colors` (H=Greys,
He=Purples, C=Blues, N=Oranges, O=Greens, S=Reds; intensity
by ionization stage). Tests that produce diagnostic plots
use this so the same ion gets the same color across every
figure in the suite and in notebooks that use
`PhotChem.plt_rate_coeffs` / `plt_sed_sigma_pi`.
- `tests/figures/.gitignore` excluding all PNG/PDF output.
Also enriches the existing CT-rate balance test docstrings with a
provenance note clarifying that Draine 2011 and the Cloudy O+H
polynomial both fit Stancil et al. 1999 (+ Barragan et al. 2006),
so their close agreement is expected rather than independent
cross-validation. Cloudy poly is only valid for T > 200 K (see
`cloudy/source/atmdat_char_tran.cpp:103-134` reference comment).
Three new test files under `tests/microphysics/`, mirroring the existing `test_ct_rate_balance.py` style: smoke tests for finiteness and sign, reference-value spot-checks against published numbers, shape / monotonicity probes, and a diagnostic-plot test per file. `test_photx_sigma.py` (Verner+96 photoionization cross sections): - Threshold zero-cutoff, finite + positive above threshold for an ion catalog covering H, He, C, N, O, S. - Reference values: H I sigma at 13.6 eV ~ 6.3e-18, He I at 24.6 eV ~ 7.4e-18, He II at 54.4 eV ~ 1.6e-18 (Draine 2011 Tables 13.1-13.3; rtol 10-20%). - High-E decay shape (~E^-3 to E^-3.5 hydrogenic asymptote). - Round-trip `get_Eth` consistency with the sigma=0 boundary. - Diagnostic plot: sigma_pi(E) for the full catalog, per-ion colored by the `photchem.py` convention. `test_rec_rate.py` (Badnell radiative + dielectronic recomb): - Pins the (Z, N) = initial-ion convention used by `get_rr_rate`, `get_dr_rate`, `get_rec_rate`. The initial ion captures the electron: H II rec to H I is queried as (1, 0); He II rec is (2, 1); O III rec is (8, 6). - H II Case B at T = 1e4 K = 2.59e-13 cm^3/s (Draine 2011 Eq. 14.6). Case A > Case B verified. - DR temperature shape: O III DR rises from T=1e3 to T=1e5 K (Badnell 2006a Fig. 1). - Consistency: total = RR + DR exactly for non-H ions. - Diagnostic plot: RR (left panel) and DR (right panel) vs T for HII / HeII / HeIII / CII / NII / OII / OIII. `test_ci_rate.py` (Voronov 1997 collisional ionization): - (Z, N) = REACTANT-ion convention pinned. - Endothermic Boltzmann sign verified (rate grows with T for H I, He I, O I). - Cold cutoff: rate is exactly 0 at T = 100 K (the U > 80 branch). - H I CI at T = 2e5 K matches Voronov fit to 10%. - Isoelectronic ordering: O I -> O II -> O III ionization rate decreases with charge at fixed T = 1e5 K. - Diagnostic plot: beta_CI(T) for the followed-ion set. All plot tests use the shared `ion_colors` fixture from `tests/conftest.py` and the LaTeX `\,` thin-space convention for labels. Plot generation is on by default; pass `pytest --no-figures` to skip. 121 tests + 5 skipped (DR not defined for Z=1 or N=0; 3 plot tests skipped under `--no-figures`).
Replace the legend-on-the-side layout with inline rotated labels via `pyathena.plt_tools.line_annotation.line_annotate`, matched to each line's slope and underlined with a white-stroke patheffect so they stay legible over neighbours. Same-color-per-ion is preserved (via the `ion_colors` fixture), so visual identification across figures stays consistent. Additional polish: * `photx_sigma_overview.png`: clamp xlim to [3, 1e3] eV and ylim to [1e-20, 1e-16] cm^2 so the most relevant range (HII + soft X-ray) fills the panel. Add minor+major log grids. * `rec_rate_overview.png`: per-ion annotation x-coordinate chosen per-line to avoid overlapping labels in the dense Badnell-RR-decline / DR-peak regions. Grids added. * `ci_rate_overview.png`: same pattern -- per-line annotation x in the steep-rise portion. Grids added. * `ct_OH_rec_comparison.png`: keep the side legend (only 3 source lines + Draine per-J decomposition, dense overlap at low T makes inline labels overlap badly). Add grid. No assertion changes; 124 tests + 2 skipped continue to pass.
Round of figure refinements based on visual review:
* `rec_rate_overview.png`:
- Add S II, S III recombination -- completes the H/He/C/N/O/S
coverage (S is in the 10-followed-coolant set and a key
observational tracer; missing-from-figure was an omission).
- Split y-labels: left panel `\alpha_{rr}`, right panel
`\alpha_{dr}`.
- Tighten ylim top to 2e-10 (was 1e-9; the upper decade was
empty space).
* `ci_rate_overview.png`:
- Extend the ion catalog to high charge states (C IV/VI,
N IV/VII, O V/VIII, S VII/XII) since the T-range goes up to
1e8 K and at those temperatures the CIE balance is dominated
by highly-ionized species.
- `ion_color` now accepts an optional `num_ions` arg in
`conftest.py` so the per-ion color gradient spans
neutral->fully-stripped (Z+1) per element when high-q stages
appear. Low-q stages keep the photchem.py default.
- ylim top -> 1e-6 (was 1e-7) to show the full plateau region
at high T.
* `photx_sigma_overview.png`:
- Add S III (Eth=34.83 eV) since S III is in the followed-ion
set and a major observational tracer.
- Hand-tune each ion's inline-label x-position to avoid label
overlaps near the H I / O I (Eth ~ 13.6) and He I / C II
(Eth ~ 24) clusters.
* `ct_OH_rec_comparison.png`:
- ylim bottom raised 1e-12 -> 3e-11 to focus on the HII-region
regime where all three sources are within factor of 2.
* Whitespace normalization: collapse all " " (multiple spaces) after
commas and colons to single space across all four test files,
matching the surrounding pyathena code style.
127 tests + 2 skipped pass.
pyathena tests + microphysics regression suite + O I CT dispatch fix
`get_ct_rec_rate(Z, N, T)` now treats (Z, N) as the REACTANT ion
(the one that gains an electron from H I and recombines), matching
the convention already used by `get_ct_ion_rate`, `get_ci_rate`,
and `RecRate.get_rec_rate`. Previously (Z, N) labelled the product,
which was inconsistent with every other rate function in the
microphysics module and confusing in the abundance-sweep loop.
Examples in the new convention:
- O+ + H I -> O + H II is queried as get_ct_rec_rate(8, 7, T).
Previously this same physical rate was get_ct_rec_rate(8, 8, T).
- Fully-neutral reactants (q = 0) now return 0 (cannot recombine
further) instead of raising or returning a stale rate.
- q > 4 falls back to the Dalgarno-like generic estimate
`1.92e-9 * q`, identical in numerical value to the previous
`(q_product + 1)` multiplier (since q_reactant = q_product + 1).
Internal data lookup keeps the original product-indexed table layout
(matches the Cloudy ctrecombdata.dat row order); the API translates
by adding 1 to N before calling `_get_index2`.
Also fixed a bug in `photchem.py:plt_rate_coeffs`: the ionize-loop
axis was plotting `get_ct_rec_rate` instead of `get_ct_ion_rate` for
the CT-ionization curves. The new code calls the correct function;
under the reactant-indexed API the per-ion (Z, N) values stored in
`Ion` work directly for both calls.
…production
Test changes for the get_ct_rec_rate API refactor:
- Split the ion catalog into CT_ION_CATALOG (neutral reactant, queries
CT-ionization) and CT_REC_CATALOG (singly-ionized reactant, queries
CT-recombination). Previously a single catalog was shared, which only
worked under the inconsistent old convention.
- Per-test (Z, N) values shifted where they query CT-recombination
(NI ct_rec now (7, 6) for N+ reactant; OII ct_rec (8, 7); etc.).
- Updated Dalgarno-high-q pin to `1.92e-9 * q_reactant`.
- Updated module-level API docstring to describe the new convention.
New diagnostic figure: ct_O_equilibrium_Draine_fig14_5.png
Two-panel plot:
- Left panel reproduces the CORRECTED Draine 2011 Figure 14.5, which
plots the steady-state ratio
[n(O+)/n(O0)] / [n(H+)/n(H0)]
versus temperature in the low-density and high-density limits.
The version in the printed 2011 Princeton book is numerically
incorrect (acknowledged on Draine's book errata page,
https://www.astro.princeton.edu/~draine/book/errata_p1.pdf,
noted 2011.05.18 by E. B. Jenkins). The corrected figure on the
errata has y in [0, 1] with both curves asymptoting to about 8/9;
this test matches the corrected version.
- Right panel computes the self-consistent x_HII and x_OII versus
density in neutral gas at fixed chi_FUV = 1, xi_CR = 2e-16 s^-1,
no photoionization. Uses pyathena `get_xHII` with grain-assisted
recombination, then closes the oxygen balance with cosmic-ray
ionization (zeta_O = 2.7 * xi_CR per Draine 2011 sec 13.7),
CT-ion, CT-rec, and RR + DR. Two temperatures (50 K and 200 K)
are shown to illustrate the temperature dependence in the diffuse
/ cold neutral medium regime.
Bug in the realistic-equilibrium right panel of ct_O_equilibrium_Draine_fig14_5.png: the O+ source rate used `get_ct_ion_rate(8, 8, T)`, which returns the J-summed total `k0r + k1r + k2r`. This is the high-density (LTE-populated J levels) expression. In the diffuse / cold neutral medium (n_H << n_crit for [OI] 63 um), O0 sits in J=2 only, so only the J=2 channel of CT-ionization fires. Mixing the J-sum CT-ion rate with the J-sum CT-rec rate (which is correctly the full sink, since all final J products radiate back down to J=2) artificially boosted x_OII by a factor of about 3-4 above what the CT-equilibrium ratio K(T) predicts. The right panel then showed x_OII > x_HII at low T, inconsistent with K(T) < 1 in the left panel. Fix: call `get_ct_ion_HII_OI_Draine11(T, sum=False)` to get the per-J components and use only the J=2 entry (k0i, equivalent to Draine's k0r). Sink remains the J-summed `get_ct_rec_rate(8, 7, T)`. After the fix, x_OII / x_HII matches K(T) from the left panel to within a couple percent (0.011 at T = 50 K, 0.31 at T = 200 K). Also drops the buggy explanatory title text from the left panel.
…P populations
The CT-ionization rate of O0 depends on the fine-structure level
populations of O 3P (J = 2, 1, 0). In the low-density limit
(n_H << n_crit) all O0 sits in the J = 2 ground state and only the
k0r CT-ion channel fires. In the high-density LTE limit
(n_H >> n_crit) the J levels thermalize and the effective CT-ion
rate is the J-summed weighted Boltzmann average.
Add two helpers in tests/microphysics/test_ct_rate_balance.py:
- `_n_crit_OI(T, x_e)` -- Draine 2011 Eqs. 14.33-14.34, the
critical density for the [OI] 63 um J=1 -> J=2 fine-structure
transition. Sets the boundary between the two regimes.
- `_f_OI_J(T, n_H_over_n_crit)` -- partial-LTE populations of the
J = 2, 1, 0 levels, parametrized by n_H / n_crit. Returns
(1, 0, 0) at low density and the full Boltzmann fractions at high
density. Uses an effective two-level interpolation factor
beta = (n_H/n_crit) / (1 + n_H/n_crit).
Wire these into both the inline right-panel computation in
test_plot_oxygen_CT_equilibrium and the helper `_K_realistic_OH`
used elsewhere. Also extends the right-panel T grid from {50, 200}
K to {50, 200, 1000, 8000} K to cover the CNM through WIM regimes.
At n_H <~ 1e3 cm^-3 in the test plot range the interpolation makes
no visible difference (beta is small), but the smooth handling will
be needed in the multi-ion sweep itself where the chemistry traverses
n_H >~ 1e4 cm^-3 regions.
The previous commit added a `beta = n_H / (n_H + n_crit)` two-level
interpolation between low-density (J=2 only) and high-density LTE
limits for the O 3P fine-structure populations, with n_crit from
Draine 2011 Eqs. 14.33-14.34. This was a hack to side-step
implementing the full level-population solve.
Pyathena already has the exact solve in `cool.get_OI_lev(nH, T, xe,
xHI, xH2)` at `cool.py:1128`, which integrates the 3-level
statistical equilibrium with H I, H2 (para + ortho), and electron
collision rates from Draine 2011 Appendix F Table F.6. Returns the
J=2 / J=1 / J=0 population fractions directly.
Switch both the inline right-panel computation and the helper
`_K_realistic_OH` to call `get_OI_lev`. Delete the now-unused
`_n_crit_OI` and `_f_OI_J` interpolation helpers. The CT-ion rate
becomes
k_CT_ion_eff = f_J=2 * k0i + f_J=1 * k1i + f_J=0 * k2i
with f_J read from the proper 3-level solve. Behavior at the
test plot densities (n_H <~ 1e3 cm^-3, n_H << n_crit) is
indistinguishable from the previous beta interpolation, but the
new path is correct across the full density range and uses the
real collision rates rather than a two-level approximation.
In the realistic-equilibrium right panel of ct_O_equilibrium_Draine_fig14_5.png the code used `zeta_O = 2.7 * xi_CR` (later 2.84) and attributed the factor to "Draine 2011 sec 13.7" / "Indriolo & McCall 2012". Both citations were incorrect. Indriolo+12 does not use a 2.84 factor; Indriolo+15 ApJ 800 40 routes the O+ source through the H+ + O charge-transfer near-resonance rather than direct CR ionization of O, so the per-atom zeta_O factor is rarely the dominant term in observational inversions. Where the ~2-3 factor does come from: - Draine 2011 ch. 13.2 covers Auger ionization theory but does not tabulate per-element zeta values. - Tielens 2005, Physics+Chem of the ISM, sec 3.4.2 lists per-element CR ionization rates approximately 2-3x zeta_H for O. - Glassgold et al. 2012, ApJ 756 157 has the detailed Auger cascade calculation. In the test's parameter range (T = 50-200 K, n_H = 1-1e3 cm^-3) the direct CR-ion source term `zeta_O * n_O0` is about 4 orders of magnitude smaller than the CT-ion source `n_HII * k_CT_ion * n_O0`, so the precise factor barely matters. The code now documents this honestly and keeps 2.84 as a labeled placeholder. Also rename `num` / `den` to `source_per_OI` / `sink_per_OII` for clarity, since they are the per-O0 source rate and the per-O+ sink rate in the steady-state ratio `n(O+)/n(O0) = source / sink`.
Items 3 + 4 of the pyathena_ct_fixes_plan MVP. Add per-ion CT
contributions to drate / crate inside `evolve_one_species`, plus the
cross-element H I / H II coupling that accumulates the sum-over-
metals CT fluxes for the two H walks.
New private helper `_ct_rate_safe(kind, Z, N, T)` wraps
`ChargeTransferRate.get_ct_{ion,rec}_rate` so that ions with no
tabulated CT data (e.g., H itself, high-q metals beyond the
Kingdon-Ferland table) return 0 instead of raising. Cleaner than
sprinkling try/except into the sweep.
New pre-pass `_compute_metal_CT_fluxes(T)` builds four arrays:
HI_drate, HII_crate, HI_crate, HII_drate. These are the
sum-over-metals contributions per cell that the two H walks need.
Pre-computed once per timestep at the start-of-step state, in the
same spirit as the existing implicit-Euler treatment of each
species in isolation.
Inside `evolve_one_species`:
- Hydrogen (`element == 'H'`): the H I and H II walks read the
appropriate pair of pre-computed arrays from the `metal_CT` tuple.
- Metals (`element not in ('H', 'He')`): the four per-ion CT terms
are added directly. CT-rec when this ion has q >= 1 (it can lose
one charge). CT-ion of this ion when it is in the ionize branch.
Creation from X^(q-1) by CT-ion (recomb branch) and from X^(q+1)
by CT-rec (ionize branch).
- Helium: skipped both ways (no tabulated CT data; H-He self-CT is
small enough to ignore at the MVP level).
`evolve_all` now invokes the pre-pass once per timestep after
`calc_radiation_field`, then threads `metal_CT` into each
`evolve_one_species` call.
The MVP test suite (124 tests / 6 skipped under --no-figures)
continues to pass. The 1D OII<->HII resonance relaxation test and
the HIIWind setup will be added separately as items 5a and 5b.
Item 5a of pyathena_ct_fixes_plan MVP: `tests/microphysics/test_OII_HII_resonance.py` is a small regression test of the new CT wiring in `PhotChem`. It avoids the full PhotChem instantiation (which needs an SB99 SED file) by building a minimal stand-in object that carries only the attributes the CT methods read: `den`, `ions`, `ct`, `rc`, `ci`, and a stub `Fphot`. Two checks: 1. `test_metal_CT_fluxes_signs` -- verifies the four sum-over-metals arrays returned by `_compute_metal_CT_fluxes` against analytic expressions for a 1-zone partially ionized cell with O at 50/50. 2. `test_OII_HII_resonance_lock` -- mass-conserving implicit Euler on the (O0, O+) sub-system at T = 1e4 K, n_H = 100 cm^-3 (low-density limit so only the J=2 channel of CT-ion fires). Confirms x_OII relaxes to the analytic CT-equilibrium value, and pins K_low(T=1e4) ~ 8/9 (Draine 2011 Figure 14.5 corrected). The test surfaced a subtle point about pyathena's CT API: the J-summed rate `get_ct_ion_rate(8, 8, T)` returns the UNWEIGHTED sum k0i + k1i + k2i, which is not the physically right rate at either density limit. For the low-density regime it overcounts by ~4x; for high-density LTE it needs a Boltzmann-weighted (5, 3, 1) combination. The CT wiring in `evolve_one_species` uses this same unweighted sum and so currently does NOT include the J-population-weighted correction; this will need a follow-up that ties into `cool.get_OI_lev` (and equivalent solves for N, S, ... neutral metals which are not yet in pyathena). Also drops the direct CR ionization source term `zeta_O * n_O0` from the realistic-equilibrium right panel of `test_plot_oxygen_CT_equilibrium`. I had cited multiple sources (Indriolo+12, Tielens 2005 sec 3.4.2) for a `zeta_O ~ 2.84 zeta_p` factor; I was unable to verify any of these citations against the actual papers, so the unverifiable factor is dropped rather than carry a number I cannot back up. The term is ~1e4 times smaller than the CT-ion source in this parameter range anyway, so the figure is unchanged at plotting resolution.
Pyathena's `ct_rate.get_ct_ion_rate(Z=8, N=8, T)` dispatches to `get_ct_ion_HII_OI_Draine11(T)` which returns the UNWEIGHTED J-sum `k0i + k1i + k2i` of the three O 3P fine-structure channels. The unweighted sum is not the physically correct per-O0 rate at any density (overcounts by ~4x in either low-density or LTE limit). For the chemistry sweep in PhotChem (which targets HII regions operating well below n_crit ~ 2e4 cm^-3 for [OI] 63 um), the low-density limit applies, so the correct per-O0 CT-ion rate uses only the J = 2 ground-state channel k0i. Implement this in `_ct_rate_safe`: when (kind='ion', Z=8, N=8), call `get_ct_ion_HII_OI_Draine11(T, sum=False)` and return only the first component. Other neutrals (C0 / N0 / S0 / Mg0 / Si0) use scalar single-channel fits (Kingdon-Ferland generic or `get_ct_ion_MgI_HII`, `get_ct_ion_SiI_HII`), so there is no per-J data to mishandle for them and no special-case is needed. Update `test_metal_CT_fluxes_signs` to expect the J = 2 rate instead of the unweighted sum for the O0 leg. The proper upgrade -- weight by exact J populations from `cool.get_OI_lev` -- is deferred to a follow-up branch (see pyathena_ct_fixes_plan.md item 9 in tigris-notes). Impact in HII-region ionization state is small (the J=2 channel dominates either way at HII temperatures), so the low-density-limit fix is adequate for the MVP. PDR / CNM chemistry would feel the difference more, but pyathena's PhotChem is not the intended workhorse there.
Add docstring + inline comments explaining the two indexing schemes in play for the OI 3P system: - Level index (0/1/2) used by `cool.coolOI`, `cool.get_OI_lev`, and most modern radiation-transfer codes: numbered by EXCITATION ENERGY starting from ground. - Channel index (0/1/2) used by Draine 2011 sec 14.7 for the per-J CT-rec / CT-ion rates k0r, k1r, k2r and k0i, k1i, k2i: labeled by the destination J of the produced O atom. For OI 3P these happen to agree numerically (Hund's 3rd rule inverts the J ordering since 2p^4 is more than half-filled, making J=2 the ground state). For other ions with normal Hund ordering (e.g., CI 3P) the two conventions would disagree. The added comments + docstrings flag this so future readers don't get confused. Also expand the `get_ct_ion_HII_OI_Draine11` docstring to warn that the `sum=True` return value (unweighted k0i+k1i+k2i) is NOT a physically valid per-O0 rate at any density (overcounts by ~4x in either low-density or LTE limit). For correct per-O0 CT-ion rate, call with sum=False and weight by the J=2/1/0 populations from `cool.get_OI_lev`.
New top-level location for per-ion metal-line cooling work that
will support the multi-ion sweep + the R-type ionization-front
benchmark vs Cloudy. Existing `pyathena/microphysics/` code is
untouched -- this is purely additive.
Layout:
pyathena/photchem/
__init__.py -- re-exports solve_5level_steady_state +
cooling_from_populations
n_level.py -- 5-level statistical equilibrium solver
coolants/ -- empty for now; will house per-ion cooling
functions (OIII, NII, SII, NI, SI, SIII)
data/ -- empty for now; offline-built CHIANTI fit
coefficients per ion
The 5-level truncation is standard for nebular HII-region coolants:
np^2 ions (CI, NII, OIII, SIII), np^3 ions (NI, OII, SII), and
np^4 ions (OI, SI) all share a 5-level ground configuration. Higher
CHIANTI levels (typically 2s 2p^3 / 2p^4 / Rydberg) sit above
60000 K = exp(-6) ~ 0.002 Boltzmann factor at HII T -- negligible
contribution to cooling and to ionization-relevant J populations.
Two-level ions (CII 2P) are handled separately.
`solve_5level_steady_state(A, E, g, T, n_coll_times_q)` builds the
5x5 rate-balance matrix per cell (spontaneous decays + collisional
de-excitation downward; detailed-balance Boltzmann upward) and
solves with the closure sum_i f_i = 1. The colliders are pre-summed
by the caller into the `n_coll_times_q[i, j]` matrix; this keeps the
solver agnostic about the collider mix (HI, e-, H2, etc.).
Smoke-tested in two limits at T = 1e5 K:
- Zero collisions: f = (1, 0, 0, 0, 0) (all in ground).
- Collisions >> A: f matches the LTE Boltzmann distribution to
~5e-7 relative error.
Next: an offline CHIANTI-v11.0.2 table builder under data/ that
extracts (level energies, A coefficients, effective collision
strengths Y(T) on a T grid) for each ion and saves per-ion .npz
files. The runtime coolant code then reads only the .npz, so
ChiantiPy is a build-time-only dependency.
`pyathena/photchem/data/build_chianti_tables.py` is the offline
extractor. Reads CHIANTI v11.0.2 via ChiantiPy, takes the lowest
N levels for each followed coolant ion (5 for the np^2 / np^3 /
np^4 ions; 2 for CII), and writes a per-ion ASCII text file with
four sections: LEVELS, A_COEFFS, T_GRID, UPSILON_E.
CLI:
python -m pyathena.photchem.data.build_chianti_tables \
[--T-min 1e3] [--T-max 1e6] [--N-T 60] [--ions OIII NII ...]
Defaults to a 60-point log T grid from 1e3 to 1e6 K, covering the
HII-region target plus PDR / CNM (~1e3 K) and a margin to early
SN-shock onset (~1e5-1e6 K) for non-equilibrium conditions.
ASCII format chosen for:
- `cat` / git-diff inspection,
- C++ readability from tigris-NCR with no extra library
(`std::ifstream` + section-keyword loop),
- consistency with the existing `tigress_coolftn_ncr.txt` precedent.
Initial build commits the 10 followed-coolant tables in
`pyathena/photchem/data/`: CI, CII, NI, NII, OI, OII, OIII, SI,
SII, SIII. ChiantiPy is a build-time-only dependency; the runtime
coolant modules will parse the .txt directly.
Two ChiantiPy v0.16 quirks documented in the build script:
- CHIANTI v11 uses ecm = -1.0 for unmeasured experimental level
energies; we fall back to ecmth (theoretical) when ecm <= 0.
- Elvlc['conf'] is an internal integer index (always 0), NOT the
configuration string. The configuration string lives in
Elvlc['term']; the spin multiplicity 2S+1 is in Elvlc['spin'],
not Elvlc['mult'] (which is the level degeneracy 2J+1).
Wgfa / Scups convention: lvl1 = lower level, lvl2 = upper level.
Splups is the v10 name for what v11 calls Scups; the script
handles either.
Companion `read_ascii()` parser is included in the same module
for convenience; the runtime coolant files will import it.
When a channel cannot use analytic d_out (e.g. CIE-table-driven
cooling with no precomputed slope column), the substep-loop
bootstrap fallback is 1-point forward FD at dT_rel = 1e-3.
Decision lineage and trade-offs:
- Analytic d_out where the chain rule is tractable (Phase 4d-a, b1,
...): exact, free, bit-stable with the C++ port.
- 1-point forward FD at dT_rel = 1e-3 as the bootstrap: 1 extra
evaluate per cell per substep, rel error < 1% on the worst
channels (Lya / HRecomb), < 1e-9 on smooth channels (Dust /
FreeFree). Adequate for the substep loop's stiffness-damping role.
- Avoid central 2-point FD in production: 2 extra evals per cell
buys 50x precision the damping role does not need.
- Avoid the historical NCR dT_rel = 2e-2 (set at
`tigris-ncr/src/photchem/photchem.hpp:224`, dlntemp_ = 0.02):
same one-extra-eval cost at 1e-3 gives an order-of-magnitude
better accuracy. The 2e-2 plateau is fine for the original C++
runs but the modern Python rewrite gets a free precision upgrade.
- Avoid per-channel optimum dT_rel: the optimum sits on the
roundoff cliff. Platform-dependent (Apple silicon vs Intel FMA
fusion shifts it), numpy-version dependent, function-shape
dependent. Risks NaN under any perturbation. Conservative 1e-3
sits in the truncation-dominated plateau where error is platform-
independent and predictable.
Calibration test (tests/chemistry/test_fd_calibration.py):
- Sweeps dT_rel over logspace(-10, -1, 19) for both 2-point central
and 1-point forward FD on 3 representative channels:
Lya steep Boltzmann exp(-T_excite/T)
HRecomb steep power law alpha_B ~ T^-0.8
Dust smooth sqrt(T) * (T - T_dust)
- Asserts the conservative defaults (1e-3 central, 1e-1 forward)
give rel error below the target AND sit to the right of the
per-curve optimum (i.e. in the truncation regime).
- 6 parametrized cases (3 channels x 2 FD methods). Plain pytest
runtime under 2 s; left in CI without a slow marker.
Empirical numbers from the sweep:
Channel Method dT_rel=1e-3 dT_rel=2e-2 Optimum
Lya central 2-pt 2.3e-5 2.3e-3 3e-6 (5e-10)
Lya forward 1-pt 6.5e-3 6.7e-2 3e-8 (3e-7)
HRecomb central 2-pt 3.2e-6 3.2e-4 1e-5 (1.7e-9)
HRecomb forward 1-pt 3.1e-3 3.0e-2 1e-7 (5e-7)
Dust central 2-pt ~1e-12 ~5e-10 ~1e-7
Dust forward 1-pt ~1e-9 ~5e-7 ~1e-7
Documentation:
- The CoolingChannel ABC docstring (pyathena/chemistry/cooling/base.py)
now spells out the FD bootstrap convention inline so any future
channel author sees it at the point of use.
- Memory: feedback_fd_bootstrap_convention.md (Claude-visible).
- This commit adds the calibration test + the inline ABC comment +
the CHANGES log entry.
Validation: 535 passed, 4 skipped (up from 529 / 4; +6 net new
calibration cases, no regressions).
Two more channels gain a non-zero d_out via the project's documented 1-point forward FD bootstrap at dT_rel = 1e-3 rather than analytic. Both share the WD01 charge-parameter chain x = 1.7 * chi * sqrt(T) / (n_e * phi) + 50, whose analytic chain rule on the surrounding Gamma / Lambda formulas is mechanically tractable but ~40 ops through quotient rules on 4-term denominators and gives no meaningful accuracy gain for the substep damping role. Refactor pattern (per channel): - Lambda / Gamma computation moves into a private _compute_X(state, out) helper that reads state.T directly. - evaluate(state, out, d_out) calls _compute_X once for the value. When d_out is supplied, evaluate snapshots state.T into the T_orig scratch, scales state.T by (1 + 1e-3), calls _compute_X a second time into the out_tp scratch, restores state.T, and assembles d_out = mu * (out_tp - out) / (T_orig * 1e-3). Each channel gains two new scratch slots: heating:photoelectric:T_orig, heating:photoelectric:out_tp, cooling:grain_rec:T_orig, cooling:grain_rec:out_tp. Tests (tests/chemistry/test_phase4d_analytic_derivatives.py): - test_photoelectric_heating_d_out_matches_FD_bootstrap and test_grain_recombination_d_out_matches_FD_bootstrap compare each channel's own forward-FD d_out against the test harness's 5-pt central FD reference. Tolerance rtol = 5e-2 because forward FD has O(dT) truncation around 3-5% in the steepest regimes; this is the calibration-sweep empirical worst case. Long-term direction (Phase 7): - Mini-RAMSES uses precomputed cool_prime / heat_prime derivative columns alongside cool / heat tables; cubic Hermite at query returns Lambda and dLambda/dT in one pass with C^1 continuity. See cooling_module.f90:516-534. This is the right end state for tabulated cooling channels (CII/OI/CI/OII at Phase 4d-e + the Phase 7 metal-line tiers); the FD bootstrap path here is the interim while pyathena.chemistry.tables doesn't yet build derivative columns. Phase 4d running total after this commit: 9 channels with a non-zero d_out path -- 7 analytic (CR / H2Diss / Dust / FreeFreeH / Lya / HICollIon / HRecomb) and 2 FD bootstrap (PE / GrainRec). 12 remaining channels still write d_out[:] = 0. Validation: 537 passed, 4 skipped (up from 535 / 4; +2 net new analytic-derivative parity tests, no regressions).
Nebular, HISmith21, H2CollDiss, H2FormationHeating, H2PumpHeating
all gain a non-zero d_out via the project's 1-point forward FD
bootstrap at dT_rel = 1e-3. Same refactor pattern as the Phase 4d-b2
PE / GrainRec batch: Lambda / Gamma body moves into a private
_compute_lambda / _compute_gamma helper that reads state.T; the
public evaluate calls the helper once for the value, then for d_out
snapshots state.T, scales by (1 + 1e-3), calls the helper again,
restores state.T, and assembles
d_out = mu * (out_tp - out) / (T_orig * 1e-3).
Two new scratch slots per channel: <ns>:T_orig, <ns>:out_tp.
Why FD bootstrap and not analytic for these five:
- All five run through composite chain rules:
Nebular -- 6-term Horner polynomial * Boltzmann *
Hummer density-reduction factor
HISmith21 -- 4 Upsilon series each piecewise in T6
H2CollDiss -- log-log interpolated rate coefficients in a
gated temperature window
H2Form, H2Pump -- HM79 ncrit chain through saturation densities
and temperature-dependent grain rates
Each analytic derivative is mechanically tractable but ~40-80 ops
through quotient / product rules with no meaningful accuracy gain
for the substep damping role.
- None of these channels are tabulatable as Lambda(T) alone (they
depend on additional combinations of xe, n_H, xHI, xH2, Z_g);
Phase 7 tabulated Hermite path (mini-RAMSES style cool_prime
columns) does NOT cover them, so the FD bootstrap stays the
long-term production path.
Tests (tests/chemistry/test_phase4d_analytic_derivatives.py):
- One test per channel comparing the channel's own forward-FD d_out
against the 5-pt central FD reference at rtol = 5e-2 (the
bootstrap O(dT) truncation worst-case from the calibration sweep).
Per-channel subnormal-tail masks; Nebular's mask is at 1e-26
because the f_red density-reduction factor flips sign in the cold
tail and the comparison becomes precision noise there.
Phase 4d running total after this commit: 14 channels with non-zero
d_out -- 7 analytic and 7 FD bootstrap. 7 remaining: H2Moseley21,
H2Gong17 (4d-d) and CII / OI / CI / OII (4d-e).
Validation: 542 passed, 4 skipped (up from 537 / 4; +5 net new
FD-bootstrap parity tests, no regressions).
- H2Moseley21Cooling and H2Gong17Cooling gain FD-bootstrap d_out (forward FD at dT_rel = 1e-3). Lambda computation moves into a private `_compute_lambda(state, out)`; evaluate calls it twice. Two new scratch slots per channel (T_orig, out_tp). - pyathena.chemistry.equilibrium new module porting get_xHII, get_xCII, coeff_kcoll_H, coeff_alpha_rr_H, coeff_alpha_gr_H from pyathena.microphysics into the chemistry-rewrite namespace, plus eq_xHII_xe(T, nH, xH2, xi_CR, G_PE, G_CI, ...) wrapper that solves jointly for (xHII, xCII, xe) via damped fixed-point iteration on xe = xHII(xe) + xCII(xe). H conservation xHI = 1 - 2*xH2 - xHII enforced on return. - test_phase4d_analytic_derivatives: _build_state gains an H-conservation assert; new _build_state_eq helper derives (xHII, xCII, xe, xHI) from eq_xHII_xe. Prior CR / H2Dissociation / H2 formation / pump / colldiss tests get conservation-correct triplets (xe ~ 1.6e-4 from xCII in molecular regime). Two new tests for Phase 4d-d at chemical equilibrium across the molecular-regime grid T = logspace(1, 4, 30) x nH = logspace(0, 6, 14). Phase 4d running total: 16/20 channels (7 analytic, 9 FD bootstrap). 4 remaining (4d-e): CII, OI, CI, OII level-pop. Validation: tests/chemistry/ 404 passed, 2 skipped. Co-Authored-By: Claude Opus 4.7 <noreply@anthropic.com>
- Add `state.u_rad: Dict[band, ndarray]` as the primary radiation representation (energy density, cgs). Bands: PE, LW, FUV (alias), EUV. `chi_for(band)` returns the Draine-normalised dimensionless field; ISRF reference densities (`U_RAD_PE_ISRF_CGS = 7.613e-14`, `U_RAD_LW_ISRF_CGS = 1.335e-14`) match tigris-ncr/src/photchem/photchem.hpp. - Add `state.zeta_pi: Dict[species, ndarray]` (photoionisation) and `state.zeta_diss: Dict[species, ndarray]` (photodissociation), replacing the prior `state.xi_ph_HI` / `state.xi_diss_H2` flat attributes. Species without a key default to zero; arrays may be scalar or `(ncell,)`. - `xi_CR` (cosmic-ray ionisation) keeps its `xi_*` naming -- the `zeta_*` / `xi_*` split now matches PDR / HII literature convention used in the rest of the chemistry rewrite.
- networks/base.py: NetworkBase gains `closure_species` (tuple of
species names solved algebraically from conservation rather than
as Newton variables) and `slow_species` (per-network default
empty). `seed_equilibrium(state)` is the new hook for analytic /
fixed-point initialisation.
- networks/ncr3.py: NCRNetwork3 declares `closure_species =
('HI',)`; HI is derived from H conservation
`xHI = 1 - 2 xH2 - xHII` in `closure(state)`. `seed_equilibrium`
wraps `equilibrium_seed.eq_xHII_xe` to fill `xHII / xCII / xe / xHI`.
Rate accessors switch from `state.xi_ph_HI` / `xi_diss_H2` to
`state.zeta_pi['HI']` / `state.zeta_diss['H2']` via a new
`_get_rate(state, dict_attr, species)` helper that returns the
scalar / array value or 0 when missing.
- solvers/_stubs.py: RadiationStub builds the dict-layout
`state.zeta_pi` / `state.zeta_diss` (and populates
`state.u_rad['PE']` / `['LW']`) instead of the flat attributes.
- tests/chemistry/{test_driver_smoke, test_ncr3_network,
test_phase4d_analytic_derivatives, test_solvers_explicit_subcycling}:
fixtures migrate to the dict layout; assertions unchanged.
Reuse `_coeff_kcoll_H` / `_coeff_alpha_rr_H` / `_coeff_alpha_gr_H` from `networks/ncr3.py` instead of duplicating the formulas. The network is the single source of truth for the rate coefficients; the equilibrium seed and the time-dependent integrator now agree by construction. File renamed from `equilibrium.py` to `equilibrium_seed.py` (rename landed in the prior commit) to reflect its restricted role: produce an initial guess for the Newton solver in `equilibrium_solver_v2`, not solve the full equilibrium itself.
Preserve the first-pass equilibrium-solver experiments before replacing them. This module hosts three candidate routines: GaussSeidel (cheap but diverges in the cold molecular regime where C and D self-couple through xe), a naive linear-space Newton (over- shoots the cold-molecular fixed point by 4 decades), and a damped variant. Kept as reference for the redesign rationale; the v2 module in the next commit supersedes this file.
Replace the v1 Gauss-Seidel / linear-space Newton prototypes (deleted in this commit) with a Newton-Raphson solver that operates in log10 space (`y_i = log10(x_i)`) on the evolved-non-closure species. Key ingredients: - Closure species (e.g. `xHI` for NCR3) are removed from the Newton variable set and recomputed from `network.closure(state)` after each accepted step. F is `C - D * x` restricted to the Newton variables. - Log-space variables linearise the multi-decade transition from atomic to molecular regimes; the cold-molecular cell that needed a 4-decade jump in linear xH2 now needs a 4-unit jump in log xH2, which Newton handles in O(10) iterations. - Armijo backtracking line search guards against overshoot. A per-step trust region cap `max_log_step = 2.0` prevents single steps from displacing y by more than 100x. - `numpy.linalg.lstsq` fallback handles singular Jacobians (the HII region's H2 row is identically zero when no H2 is present). Returns a `SolveResult` with per-cell `converged` mask and inf-norm residual; the regression test in the next commit calls `max_iter = 60`, `tol = 1e-6` and reaches all 5 ISM-strip cells.
Time-marching alternative to ExplicitSubcyclingSolver for solving `dx/dt = f(x)` to t_target. Verbatim port of OrdinaryDiffEq.jl's `Rosenbrock23ConstantCache + perform_step!`: - Constants `d = 1/(2+sqrt 2)`, `c32 = 6+sqrt 2`. Stage 1-2 build the order-2 propagated solution, stage 3 produces the order-3 embedded estimate `utilde = (h/6) (k1 - 2 k2 + k3)`. - W = J - I/(h d) is LU-decomposed once per step via numpy.linalg.solve per cell, then reused for 3 RHS vectors. - Mixed-tolerance norm `sc = atol + rtol * max(|y|, |y_new|)`, inf-norm over species per cell. Step accepted when err_norm <= 1 or h <= h_min. Step-size update `factor = 0.9 * err_norm^(-1/3)` clipped to [0.1, 5]; exponent is `1/(p+1)` for order p = 2. - Per-cell independent h and t -- vectorisation comes from evaluating f / J as `(nvar, ncell)` and looping the per-cell solve, not from sharing a strip-wide dt. Done cells are masked so they neither advance nor reject. - Optional `y_lo` / `y_hi` clips after each accepted step for positivity / physical bounds. _rosenbrock_chemistry_adapter.make_chemistry_f_jac wraps a Network + ChemState into the (f, jac) closures the integrator needs. f evaluates `C - D * x` on the Newton variables; jac uses one-sided forward FD with `eps_J * max(|y|, floor_scale)` and reuses the analytic-vs-FD switching the v2 equilibrium solver already does.
- test_equilibrium_solver: regression test. Builds a 5-cell ISM strip (cold molecular nH=1e3 and nH=1e6, CNM, WNM, HII region) via `_build_strip`, seeds with `NCRNetwork3.seed_equilibrium`, and asserts `solve_equilibrium(max_iter=60, tol=1e-6)` converges on every cell. `_CELLS` is shared with the figure module. - test_equilibrium_solver_convergence_figure: per-cell time-marching comparison. Each cell is integrated independently by Rosenbrock23 and by ExplicitSubcyclingSolver (BE-marched, log- spaced sample times). The 2 x 5 figure plots species fractions (top) and signed relative imbalance `(C - D x) / (C + D x)` (bottom) per Newton species, with Rosenbrock solid and BE dashed, colour-coded `xHI = C2`, `xHII = C3`, `xH2 = C0`. Writes `tests/figures/equilibrium_solver_convergence.png` (regenerated by the test; not tracked). - Per-cell `tlim_yr` reflects the dominant relaxation timescale of the dominant evolved species (HII recomb in the HII / WNM / CNM cells, H2 formation in the cold molecular ones); the HII region tlim was bumped from 1e4 to 5e4 to clear the ~1.2e4 yr recomb timescale, and WNM trimmed from 1e8 to 5e7.
Below ~700 K all four collisional dissociation rate coefficients (`k9l`, `k9h`, `k10l`, `k10h`) underflow to 0 from `exp(-43900/T)` and similar, giving `log10(0) = -inf` runtime warnings. The result is discarded by the `T > temp_coll` mask on the next line, so the warnings are spurious; wrap the rate composition in `np.errstate(divide='ignore', invalid='ignore')` to suppress them locally without perturbing the numbers.
…aptive substep Port `UpdateChemistry` from `tigris-ncr/src/photchem/ncr_solver.hpp:513-571` and the post-step adaptive control from Zier+ 2024 (AREPO-RT, MNRAS 533, 268), Section 4.1.1. - Replace the three independent per-row implicit-Euler updates on (HI, HII, H2) with a 2x2 Cramer joint solve for (xH2, xHII) that substitutes `x_HI = 1 - 2 x_H2 - x_HII` into the C terms algebraically, followed by RAMSES-style clipping (renorm if `x_HII + 2 x_H2 > 1`, asymmetric clip branched on `x_H2 < 0.25`) and closure `x_HI = 1 - 2 x_H2 - x_HII`. The old per-row form did not preserve hydrogen conservation when the three species saw each other through stale entry values; the Cramer form is conservation-respecting by construction. - ChemistryConfig gains `f_chem_cap` (default 0.1) and `f_chem_target` (default 0.05) for the post-step adaptive controller. The metric is `f = max(|delta temp_mu|/temp_mu, |delta x_HI|, |delta x_HII|, |delta(2 x_H2)|)` -- relative on the substep-invariant `temp_mu = T / mu` (chemistry inside one substep changes mu without doing work; `T` would otherwise appear to drift while `temp_mu` is constant), absolute on the H species (since `x_HI + x_HII + 2 x_H2 = 1`, each term sits in `[0, 1]` and an absolute cap of 0.1 means "no more than 10 per cent of the H budget moves between species per substep" -- relative on x would over-constrain dt for trace species). - Substep flow: snapshot `T`, `temp_mu_old`, `(xHI, xHII, xH2)` at entry; cap dt at `min(_estimate_dt_sub, _dt_sub_next, dt_remaining)` where `_dt_sub_next` is the forward-controller carry-over from the previous substep. Try T update + chemistry at that dt; if observed `f > f_chem_cap`, restore the snapshot, halve dt, retry (existing T-cap path inside `_attempt_temp_mu_step` is unchanged). On accept update `_dt_sub_next = dt * min(2.0, f_chem_target / f)` (Kannan+ 2019 / Zier+ 2024 2x growth cap). - Add 11 ncell scratch buffers: `solver:cramer_a..f`, `solver:cramer_det`, `solver:xHI_pre`, `solver:xHII_pre`, `solver:xH2_pre`, `solver:f_denom`. - Rewrite the hand-rolled-BE regression test against the new Cramer expression: same idea (compare against an analytic closed-form at a tiny dt) but the closed form is now the 2x2 Cramer solution, not the per-row BE formula. Validation: `tests/chemistry/test_solvers_explicit_subcycling.py` all 8 pass.
- pyathena/chemistry/solvers/rosenbrock.py: `integrate_rosenbrock23` step-size bounds (`h_init`, `h_min`, `h_max`) default to fractions of `t_target` rather than absolute cgs-second scalars (`1e-3`, `1e-12`, `1e-20` previously). Defaults are `h_init = 1e-6 t_target`, `h_min = 1e-12 t_target`, `h_max = t_target` (natural ceiling: no single step covers more than one full window). Signature changes the three knobs from `float` to `Optional[float]` so explicit callers behave the same. - Drop `pyathena/chemistry/equilibrium_solver_v2.py`: the Newton-Armijo equilibrium solver is removed since the production path now reaches equilibrium by time-marching the BE substep (Cramer + Zier-style adaptive control); no root-finder is needed. - Drop `tests/chemistry/test_equilibrium_solver.py`: the Newton regression test went with the solver; the test had been holding the shared `_CELLS` / `_build_strip` fixtures only. - Rename `tests/chemistry/test_equilibrium_solver_convergence_figure.py` to `tests/chemistry/test_chem_solver_comparison.py` and pull the fixtures inline. The `convergence_figure` suffix referred to the Newton convergence test, which no longer exists. The renamed file hosts the time-marching Rosenbrock-vs-BE-vs-BE-reference figure.
Rebuild the time-marching comparison around the new BE substep (Cramer joint solve + Zier-style forward controller). Three solvers share a common dt scale and the same `_CELLS` strip: - Rosenbrock23 (`rtol = 1e-2`, `atol = 1e-12`) -- adaptive reference for embedded local-error control. - backward-Euler with the 10 percent rule (`f_chem_cap = 0.1`, `f_chem_target = 0.05`) -- production setting. - BE reference with the 0.1 percent rule (`f_chem_cap = 0.001`), seeded with `dt_init_frac = _DT_INIT_FRAC` so the early transient is resolved before the forward controller takes over -- treated as ground truth. The earlier file had three "smoke" tests asserting only `res.shape[0] >= 2`, which a solver returning all-NaN would pass. Replace them with real regression checks: - `test_state_finite_bounded_and_H_conserved`: every recorded trajectory point is finite, every species fraction sits in `[0, 1]`, and `x_HI + x_HII + 2 x_H2 == 1` to round-off. Catches NaN explosions, broken positivity clipping, and silent hydrogen leaks in the Cramer joint solve. - `test_solvers_agree_at_tlim`: at `t = tlim`, Rosenbrock and BE 10 percent both agree with BE reference on `(x_HI, x_HII, 2 x_H2)` within solver-appropriate tolerances (atol = 1e-2 for Rosenbrock, atol = 0.05 for BE). This is the main regression check -- catches Cramer or forward-controller drift, Rosenbrock atol / y_lo misconfiguration, or any future solver replacement that quietly converges to a different fixed point. Other changes: - BE is driven one substep at a time via `solver._do_one_substep` so the recorded trajectory reflects the solver's actual adaptive cadence; the old `solver.step(dt_call)` loop conflated sample cadence with substep count. - BE 10 percent shows as markers (no line) on the figure -- its substep cadence is meaningful to display per-step. Rosenbrock shows as solid line with dots, BE reference as a thin solid line as the trusted curve. - `_CELLS` gains `xHI_init`; `_build_strip` derives `x_HII` from hydrogen conservation. The HII region initial state now starts near its equilibrium (xHI ~ 1e-4) so the figure does not waste the budget on the initial ionisation ramp. - Initial dt fraction `_DT_INIT_FRAC = 1e-5` is a module-level constant shared by Rosenbrock `h_min`, BE-ref `dt_init_frac`, and the figure xlim so the three solvers and the plot start from the same dt scale. - Step count + wall-time summary printed to stdout (with `-s`) and written to `tests/figures/chem_solver_comparison_summary.txt`. - ncr3 collisional-dissociation `log10(0)` warning silenced (was noise; result is masked out by the `T > temp_coll` gate). - `_SEC_PER_YR` derived from `astropy.units` so the value stays consistent with the rest of the codebase.
Print the step-count / wall-time table to stdout only (visible with `pytest -s`); the same numbers are annotated on the figure, so the text file was redundant.
ltlancas
reviewed
Jul 17, 2026
Zero a band's mean photoionization cross section only when the whole band lies below the ionization threshold energy (band minimum wavelength past the threshold), instead of comparing the energy-weighted mean wavelength. The mean test wrongly zeroed bands that straddle the threshold, discarding real ionizing photons; the effect grows with band width.
The four level-population channels gain a temperature derivative via the forward finite-difference bootstrap used by the other 4d channels: the Lambda computation moves into _compute_lambda(state, out), evaluate calls it at T and T*(1+1e-3), and the difference quotient is multiplied by mu_at_entry. Two new scratch slots per channel (T_orig, out_tp). An analytic derivative would have to chain through the steady-state 2-/3-level solve. Tests compare against a 5-point central reference at rtol 5e-2 on the broad (T, nH) grid. New _slope_mask helper keeps cells whose log-slope |d ln Lambda / d ln T| exceeds 0.03 instead of applying an absolute |d_out| floor: the collisional d_out amplitude scales with nH, so an absolute floor masks entire low-density rows, while the log-slope cut drops only the dLambda/dT sign-change cells where the forward-vs-central curvature term dominates. CII and CI additionally mask a 0.02-dex neighborhood of their piecewise collision-partner boundaries (T = 500 K and T = 1000 K). Phase 4d complete: 20 of 20 channels have a d_out path. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
ltlancas
reviewed
Jul 22, 2026
| # Compute recombination rate of the ion with one more charge | ||
| n_p1 = den[idx+1] | ||
| rc_rate_p1 = n_p1*ne*self.rc.get_rec_rate(Z, N-1, T) | ||
| ci_rate = ne*self.ci.get_ci_rate(Z, N, T) |
Collaborator
There was a problem hiding this comment.
This is missing a factor of n no?
Owner
Author
There was a problem hiding this comment.
Confusing, but it's not. See:
pyathena/pyathena/microphysics/photchem.py
Lines 714 to 718 in 1370146
The convention (noted at L714-715) is that destruction rates are per-ion [1/time] while creation rates are volumetric [1/time/volume]. ci_rate is destruction of the current ion (X^q + e -> X^(q+1)), so it's nek only — the density factor comes in through the implicit-Euler denominator at L718, n_after = (n + cratedt)/(1 + drate*dt), which multiplies drate by n_after.
Collaborator
There was a problem hiding this comment.
Makes sense, I was going to re-write this a bit in order to use scipy solve_ivp and make things a little more explicitly conservative by constructing the time derivatives as a matrix of rates which are then multiplied by densities. It seems like you've set things up in this file to work with the HII wind problem though, so should I potentially put that somewhere else?
Owner
Author
There was a problem hiding this comment.
I would not restructure the ODE solver inside photchem.py. Backward-Euler is the most robust, inexpensive choice here (although only first-order accurate), and I think it will work well for our use case (NCR hydrogen + multi-ion chemistry).
Perhaps the right place is here:
https://github.com/jeonggyukim/pyathena/tree/photchem-full/pyathena/chemistry/solvers
where I have been testing a Rosenbrock implicit solver. In the future, it would be great if we can extend this to more complicated chemistry networks (e.g. GOW17) that may require more accurate, conservative ODE solver.
Adds tests/chemistry/perf/, skipped unless `pytest --bench`. The `--bench` flag is declared in tests/conftest.py because pytest only calls pytest_addoption on initial conftests; a perf-local declaration would leave the option unknown when the suite runs as `pytest tests/chemistry`. The result the phase turns on is not one the plan asked for. The solver takes a strip-MIN dt_sub, so every cell walks the substep count of the stiffest cell in its strip. A 256-cell strip spanning five ISM regimes charges 37,888 cell-substeps against an ideal of 11,161 (3.39x). Unsorted splitting recovers nothing, because a random chunk almost always contains the stiffest regime. Sorting by stiffness and splitting into four gives 15,232 (1.36x), saturating near 1.1x by eight chunks. The unit is cell-substeps rather than seconds, so the counts transfer to the C++ port unchanged. The plan's own SIMD criteria could not be evaluated with a Python proxy. The strip-size sweep measures NumPy per-call dispatch, not vectorisation: 256 cells cost 1.07x the wall time of one cell at equal substep count. The branch-free versus branched numba ratio came out at exactly the flop ratio (0.50x), which cannot separate "neither vectorises" from "both vectorise equally". Benchmark outputs under results/ are git-ignored, following the tests/figures/ precedent. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
CoolingChannels.update summed the per-hydrogen channel outputs and handed them to a solver expecting a volumetric rate, so solver:net_cool was short one factor of n_H. Channels emit n_H * Lambda_NCR under the per-H^2 convention of tables/chianti_v11/build_cool.py, where the volumetric rate is n_H^2 * Lambda_NCR. Because channel output scales as n_H and inv_heat_cap_per_temp_mu = (gamma-1)/(n_H mu_hyd k_B) scales as 1/n_H, the two cancelled exactly and d(temp_mu)/dt came out independent of density. The solver was correct only at n_H = 1 cm^-3: cooling ran 1/n_H too fast in diffuse gas and n_H too slow in dense gas. Measured over four decades at fixed T and composition, the pre-fix rate was constant to five significant figures where it must be linear in n_H. The C++ applies the same multiplication in the solver rather than the channel (tigris-ncr/src/photchem/ncr_solver.hpp:263, var.cool_rate = var.nH * cool * ...). Per-regime effect over a 1 Myr step, substep count before to after: hot (n_H 1e-2) 350 to 1, HII (1e1) 42 to 142, CNM (5e1) 27 to 40, cold molecular (1e3) 1 to 37. Gas at 1e6 K fell to 1.07e3 K and now falls to 9.51e5 K. test_aggregator_update_matches_channel_sum asserted the old contract and had to change; it codified the defect. New test_cooling_rate_scales_linearly_with_density asserts d(temp_mu)/dt is proportional to n_H across 1e-2 to 1e2 at rtol 1e-12. It runs against a stub channel with a density-independent Lambda_NCR, not the production set: level-population cooling saturates above its critical density, and photoelectric heating and grain recombination depend on psi = G sqrt(T)/n_e, so the production channels carry a genuine density dependence of about 18 percent over that range. Note the defect is invisible at n_H = 1, so a fixture pinned there would not catch it. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
ChargeTransferRate accepts the table interpolation modes; the constructor no longer raises NotImplementedError for them. This was the last of the three rate modules without a table path. Grid and encoding match ci_rate and rec_rate exactly (2000 points, 1 K to 1e9 K, shared log-T axis), so one temperature index can drive a lookup across all three. Two tables, one column per (Z, N) reactant: 114 ionization columns (92 identically zero) and 110 recombination. Columns are filled by calling the analytic getters, so the four hand-fitted ions and the per-ion Tmin/Tmax clamping are captured without the builder knowing about them. Temperature-independent charge branches stay outside the table. Recombination uses a clamp-aware path. Tabulating the analytic form directly cost 8.9e-2 worst-case relative error, entirely from the clamp kink: the analytic path clips T and then evaluates, leaving a slope discontinuity at each window edge, and interpolating across the straddling bin is what loses the accuracy. Refining the grid does not help (8.9e-2, 7.4e-3, 1.0e-2, 7.5e-4 across _TAB_N 2000 to 16000, non-monotonic, depending on whether a node lands near each kink). Recombination columns are therefore built from the unclamped fit with the clamp applied to T at lookup, which moves the kink out of the interpolant and keeps the clamp exact. Accuracy where the rate exceeds 1e-4 of each ion's peak: mode ionization recombination LogLog 7.5e-3 2.8e-4 Nqt2 1.1e-2 9.9e-3 Nqt1 1.5e-2 1.1e-2 Recombination under LogLog now holds the same rtol 1e-3 bar as ci_rate. Ionization keeps the clamped-form path at 7.5e-3: unclamping it would also need exp(-dE * 1e4 / T) at lookup, and a transcendental in the lookup works against tabulating. Only 8 of 120 rows carry a nonzero dE_rate1. The significance mask matters for C II, whose recombination fit has a root at 5397 K just below its Tmin2 of 5500 K. At the clamped floor it evaluates to 0.2 percent of its leading term, so relative error there is meaningless: unmasked it reads 3.3e-2 against 2.9e-3 inside the fit's own window. Exact remains the default, so the rtol 1e-12 parity test against microphysics is unaffected. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
TemperatureTable holds one copy of the build-and-look-up code that ci_rate, rec_rate and ct_rate each grew separately, so cooling channels can tabulate their temperature-only factors without a fourth. It uses the same grid as the three rate modules (2000 points, 1 K to 1e9 K), which is what lets one log-T index drive a lookup across every tabulated quantity in a cell. lookup(T, out, scratch) is allocation-free, which is the property that made this worth extracting: the substep loop is allocation-free by design (assert_no_alloc), so a lookup that allocated seven intermediate arrays per call could not be used there. Every step writes through out=, and np.take replaces fancy indexing. It also carries the clamp-aware pattern from the charge-transfer work. Pass an unclamped fn plus its (Tmin, Tmax) window and lookup clips T into it, keeping the slope discontinuity out of the interpolant. test_clamp_beats_tabulating_the_clamped_form pins that as a measured comparison rather than a claim in a comment. Accuracy on the free-free temperature factor: LogLog 4.8e-7, Nqt2 8.4e-4, Nqt1 1.1e-3. test_grid_endpoints_are_reproduced caught a real defect while this was being written. Clipping the float index to TAB_N - 2 before the floor discards the top node, since frac can then never reach 1, so the last grid point returned its neighbour's value. No channel uses this yet; wiring follows. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The helper as first committed only offered lookup(), which recomputed the axis position for every table. That does not deliver the benefit the module docstring claimed for it. Finding a cell's position costs a transcendental; interpolating at a known position does not. So tabulation only pays if a cell resolves the position once and interpolates many times. Lyman alpha is the case that exposed it: its two temperature-only factors replace three analytic transcendentals, but two independent LogLog lookups cost two transcendentals each. Four to replace three is a loss. One locate plus two interp calls costs two. locate(T, idx, frac) and interp(idx, frac, out, lo, hi) are now separate; lookup is the two back to back and remains the only option for a clamped table, whose clip changes T and so cannot share a position. Both stay allocation-free: interp reaches the upper node through a stored view onto _tab[1:] rather than forming idx + 1, and locate truncates with np.copyto(casting='unsafe') rather than astype. Scratch order is now (idx, frac, lo, hi). Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Adds CoolingChannel.T_FADE = (T_full, T_off) and ramps it in the aggregator with the C2 Hermite smootherstep. This is the first mechanism in the cooling path able to express that a channel does not apply at a given temperature; until now CoolingChannels.update called every channel's evaluate unconditionally at every temperature. The derivative uses the product rule, dL_eff/dT = w dL/dT + L dw/dT. Scaling d_out alone would drop the term carrying the switching itself and leave d_net_cool_d_temp_mu wrong exactly inside the ramp. C2 rather than C1 because the solver differentiates the cooling function. Fixes NaN in solver:net_cool above 1.13e7 K. The C I three-level solve overflows there and NaN * 0 is still NaN, so a channel whose abundance was 1e-20 poisoned the whole cooling rate. Channels wholly outside their window are now skipped rather than evaluated, and where the weight is zero the result is assigned zero rather than multiplied by it. Fades applied to the three dust-coupled channels over the C++ hot-gas window, on the thermal-sputtering argument, and to the four fine-structure channels a decade above where CIE puts their ion fraction below 1e-4, as margin against underionization. At nH = 1e-2, T = 1e6 K with metals zero, net Lambda_e goes from 3.70e-23 to 2.46e-24 against a metal-free CIE reference of 5.00e-24; the remaining factor of two is helium, which the network does not track. The NaN was found by the new CIE consistency test on its first run. With photoionization and cosmic rays off the network relaxes to collisional ionization equilibrium, so its cooling is directly comparable against an independent CIE table. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
The fade figure showed the test channel's invented (1e4, 1e5) window, which is not a value anything uses. It now plots every production T_FADE on one axis so the windows can be checked against where each ion or grain population is expected to survive, keeping the machinery verification on the test channel below it. Two things the figure makes visible. The derivative panel shows the analytic curve tracking the finite difference through the ramp, including the excursion to -0.95 that is the L dw/dT term; scaling d_out alone would leave it flat at +1 there. And the dust window (2e4, 3.5e4) is 0.24 dex against a full decade for the fine-structure channels, so its dw/dT is about sixty times larger. Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
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
2 participants
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.
Photochemistry
PhotChemclass. Calculates SED-averaged cross sections, mean photon energy, etc.. For a given radial density, temperature profile of ionized gas, evolve a set of ion abundances (SpeciesSet) to chemical equilibrium.TODO