Skip to content

Support pulse-height tallies in shared secondary mode - #4089

Open
GuySten wants to merge 5 commits into
openmc-dev:developfrom
GuySten:shared-pht
Open

Support pulse-height tallies in shared secondary mode#4089
GuySten wants to merge 5 commits into
openmc-dev:developfrom
GuySten:shared-pht

Conversation

@GuySten

@GuySten GuySten commented Aug 31, 2026

Copy link
Copy Markdown
Contributor

Support pulse-height tallies with the shared secondary bank

Removes the restriction that forced shared_secondary_bank off whenever a
pulse-height tally was present, and implements the deferred per-history scoring
described in the TODO that guard carried.

Problem

A pulse-height tally scores once per source history, in the bin containing the
total energy that history's entire particle tree deposited in a cell. In the
default transport modes the whole tree is carried by a single Particle, so
pht_storage_ already holds the per-history total when event_death() runs.

Under the shared secondary bank each generation is transported as a fresh set of
Particle objects, load balanced across MPI ranks between generations. Three
things broke:

  1. pht_storage_ is a per-Particle member reset in initialize_particle_track,
    so every descendant's deposition was dropped from its root's total.
  2. event_death() scored unconditionally, so each fragment of one physical pulse
    was scored as a separate pulse, shredding the spectrum toward low energy.
  3. The subtraction removing a secondary's energy from its parent was never
    implemented on this path. event_revive_from_secondary() skipped it on the
    grounds that create_secondary() had already done it; create_secondary()
    contained no pulse-height code. The compensating pre-add in
    sample_secondary_photons() still ran unconditionally.

Approach

Every particle carries a root index: the dense global index of its primary in
the source bank, in [0, n_particles). It is propagated through SourceSite, so
it survives bank sorting and MPI migration.

Each Particle still accumulates its own fragment as before. At death the
fragment is staged in a per-thread buffer keyed by root index rather than scored.
Once the generation loop drains, contributions are routed by MPI_Alltoallv to
the rank owning each root, summed, and scored once per history.

Root index rather than parent_id/progeny_id because those identify the
immediate parent and would require walking a chain through generations that have
already been discarded. It is also dense, so it indexes an array directly and its
owner rank is a binary search in the phase-1 partition.

simulation::phase1_work_index snapshots that partition, because
calculate_work() overwrites work_index on every secondary generation.

Changes

  • particle_data.hSourceSite::root_index, ParticleData::root_index_ and accessors
  • particle.cpp — propagate the root index; implement the parent-side subtraction in
    create_secondary(); route event_death() to staging in shared mode
  • tallies/pulse_height.{h,cpp} (new) — staging buffers, owner exchange, deferred scoring
  • simulation.cpp / simulation.hphase1_work_index; init and finalize hooks in both
    shared drivers
  • tally_scoring.cppscore_pulse_height_tally() takes the energy vector as an argument
    instead of reading p.pht_storage()
  • initialize.cpp — remove check_pulse_height_compatibility()

Testing

New tests in tests/unit_tests/test_pulse_height.py:

  • Count conservation (photon/neutron × shared/local). Fixed-source results are
    normalized per source particle, so the tally summed over all bins must be exactly
    1.0. Per-track scoring pushes this above one; a dropped history pushes it below.
  • Thin and thick detector mode comparison. Exact comparison is not possible:
    compute_particle_id() and compute_transport_seed() both branch on the shared
    bank, so the two modes sample different random number streams. Mean deposited
    energy and per-bin shape are compared against combined standard errors.

Notes for reviewers

  • Regression results need regenerating. tests/regression_tests/pulse_height
    already parametrizes over shared/local, but its stored shared_* results are
    byte-identical to the local ones because the guard was silently disabling the
    feature. Run with --update and check the regenerated spectra against the local
    ones before committing them.
  • MPI coverage. The unit tests run single-rank. The cross-rank aggregation is
    only exercised with more than one process.
  • Reproducibility. Summation order into a given history's accumulator depends
    on thread scheduling and on which rank a descendant landed on, so results differ
    at round-off between thread counts. Given that sort_bank() exists specifically
    for reproducibility, this is a regression in kind, though only in the last bits.
    A deterministic variant (stage per generation, sort by root index, fold in order)
    costs a sort per generation and can be added behind a flag if wanted.
  • Memory. Staging holds one record per track that actually deposited, not per
    track; in a typical detector problem most histories never reach a pulse-height
    cell. The per-rank totals array is 8 * n_owned_roots * n_pulse_height_cells
    bytes.

Checklist

  • I have performed a self-review of my own code
  • I have run clang-format (version 18) on any C++ source files (if applicable)
  • I have followed the style guidelines for Python source files (if applicable)
  • I have made corresponding changes to the documentation (if applicable)
  • I have added tests that prove my fix is effective or that my feature works (if applicable)

@GuySten
GuySten marked this pull request as ready for review August 31, 2026 18:02
@GuySten
GuySten requested a review from paulromano as a code owner August 31, 2026 18:02
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant