Skip to content

AngleDistribution.forward_fraction returns uninitialized memory for ACE-derived distributions, so removal_xs is nondeterministic there #21

Description

@shimwell

AngleDistribution.forward_fraction allocates its result with np.empty and then fills only the entries whose distribution is a Legendre or a Tabulated1D:

fractions = np.empty(len(self.energy))

for i, mu_i in enumerate(self.mu):
    if isinstance(mu_i, Legendre):
        ...
    elif isinstance(mu_i, Tabulated1D):
        ...

return fractions

An AngleDistribution built by AngleDistribution.from_ace contains neither type. Its entries are Tabular (for a tabulated cosine distribution or the 32 equiprobable bins) and Uniform (for isotropic), so every entry falls through both branches and is returned as whatever np.empty happened to hand back.

That propagates: IncidentNeutron.removal_xs calls forward_fraction on the elastic angular distribution, so for any nuclide read from an ACE table it returns garbage — and garbage that changes between calls, because the value depends on what was last in that block of memory:

>>> import endf, endf.ace
>>> t = endf.ace.get_tables('tests/Li6.ace')[0]
>>> n = endf.IncidentNeutron.from_ace(t)
>>> n.removal_xs('294K', 0.0).y[0]
np.float64(10.40895)
>>> for cutoff in (-1.0, 0.0, 0.5):
...     print(cutoff, n.removal_xs('294K', cutoff).y[0])
-1.0 10.40895
0.0 -97.93729010250001
0.5 -97.93729010250001

The same call gives a different answer depending on what ran before it. The first line is the total cross section unchanged, i.e. a forward fraction that happened to be zero; the others are negative removal cross sections, which is not a physical answer.

The ENDF path is unaffected — AngleDistribution.from_dict produces Legendre and Tabulated1D — which is presumably why the tests in tests/test_removal_xs.py (all ENDF) pass.

Two things would fix it:

  1. fractions = np.zeros(len(self.energy)) instead of np.empty, so an unhandled type is at least deterministic.
  2. Handle the two ACE shapes, which both have a well-defined answer. For a Tabular, the forward fraction is 1 - C(mu_cutoff) where C is the CDF the class already computes; for a Uniform over [a, b] it is (b - max(a, mu_cutoff)) / (b - a), clamped to [0, 1].

With (2), removal_xs becomes meaningful for ACE data rather than merely deterministic.

Found while porting the reader to Rust: there is no behaviour to reproduce here, so the Rust forward_fraction implements (2), and its removal cross section is compared against the Python one only on the ENDF path, where the Python answer is well-defined.

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions