Skip to content

Fix sign guard on the cliq liquid extrapolation in thermr#400

Open
ramic-k wants to merge 1 commit into
njoy:developfrom
ramic-k:fix/thermr-cliq-sign-guard
Open

Fix sign guard on the cliq liquid extrapolation in thermr#400
ramic-k wants to merge 1 commit into
njoy:developfrom
ramic-k:fix/thermr-cliq-sign-guard

Conversation

@ramic-k

@ramic-k ramic-k commented Jun 7, 2026

Copy link
Copy Markdown

Fixes #399.

What

The small-α extrapolation in sig(),

s = sab(1,1) + log(alpha(1)/a)/2 - cliq*b**2/a

assumes cliq >= 0, which holds when S(α₁,β) decreases from β₁ to β₂.
The activation guard tests decay along the α axis instead
(sab(1,1) > sab(2,1)). A scattering law whose β=0 row decays in α
while S(α₁,β) rises from β=0 to β₂ — e.g. a coherent one-phonon law
for a crystalline powder, where the acoustic intensity at small Q sits
in the first nonzero β bin — makes cliq negative, so the
extrapolation grows like exp(+|cliq|·b²/a) with a clamped at 1e-6
and s clamped only from below. The result is unphysical thermal
inelastic cross sections (up to ~1e150+ barns with the reproducer in
the linked issue; ~1e91 b with the real graphite evaluation we first
hit this on) or runaway adaptive refinement.

This PR requires decay along both axes before activating cliq, at
both computation sites:

if (sab(1,1).gt.sab(2,1).and.sab(1,1).gt.sab(1,2))&
  cliq=(sab(1,1)-sab(1,2))*alpha(1)/beta(2)**2

Why existing results are unaffected

For liquid evaluations — the intended target of the extrapolation —
both decays hold (checked explicitly for ENDF/B-VIII.0 H-in-H₂O:
S(α₁,β₁)=6.17e5 > S(α₁,β₂)=6.03e3), so cliq keeps its current
value. Materials that never activated the guard keep cliq=0. Only the
pathological negative-cliq case changes, from garbage to the SCT/
tabulated treatment of the same region.

On macOS/arm64 the repo test suite gives identical results with and
without this change (same pass/fail sets; the 14 failures present on
this platform are present on unpatched develop as well).

The small-alpha extrapolation in sig(),

    s = sab(1,1) + log(alpha(1)/a)/2 - cliq*b**2/a

assumes cliq >= 0, which holds when S(alpha1,beta) decreases from
beta1 to beta2. The activation guard, however, tests decay along the
alpha axis (sab(1,1) > sab(2,1)). For a scattering law whose beta=0
row decays in alpha while S(alpha1,beta) rises from beta=0 to beta2
-- e.g. coherent one-phonon laws for crystalline powders, where the
acoustic intensity at alpha1 sits in the first nonzero beta bin --
cliq goes negative and the extrapolation grows like
exp(+|cliq|*b**2/a) with a clamped at 1e-6 and no upper bound on s,
producing unphysical inelastic cross sections (up to ~1e91 b) above
a few tenths of an eV, or runaway adaptive refinement.

Require decay along both axes before activating cliq, at both
computation sites. Liquid laws satisfy both conditions already
(verified against the ENDF/B-VIII.0 H-in-H2O evaluation), so existing
evaluations are unaffected.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant