From e19b2910ddd5ad172d82ade01f04874c2ad42339 Mon Sep 17 00:00:00 2001 From: GuySten Date: Mon, 31 Aug 2026 18:42:29 +0300 Subject: [PATCH 1/6] claudes initial work --- CMakeLists.txt | 1 + include/openmc/particle_data.h | 13 ++ include/openmc/simulation.h | 6 + include/openmc/tallies/pulse_height.h | 71 ++++++++ include/openmc/tallies/tally_scoring.h | 9 +- src/finalize.cpp | 2 + src/initialize.cpp | 26 --- src/particle.cpp | 43 ++++- src/simulation.cpp | 43 +++++ src/tallies/pulse_height.cpp | 223 +++++++++++++++++++++++++ src/tallies/tally_scoring.cpp | 5 +- tests/unit_tests/test_pulse_height.py | 203 ++++++++++++++++++++++ 12 files changed, 613 insertions(+), 32 deletions(-) create mode 100644 include/openmc/tallies/pulse_height.h create mode 100644 src/tallies/pulse_height.cpp diff --git a/CMakeLists.txt b/CMakeLists.txt index 62c2ac8a151..fd931d04ff0 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -478,6 +478,7 @@ list(APPEND libopenmc_SOURCES src/tallies/filter_universe.cpp src/tallies/filter_weight.cpp src/tallies/filter_zernike.cpp + src/tallies/pulse_height.cpp src/tallies/tally.cpp src/tallies/tally_scoring.cpp src/tallies/trigger.cpp diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 24b7eb53c0a..20ee6f567c2 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -52,6 +52,12 @@ struct SourceSite { int parent_nuclide {-1}; int64_t parent_id {0}; int64_t progeny_id {0}; + // Dense global index, in [0, n_particles), of the primary at the root of + // this particle's tree. Propagated unchanged through every secondary + // generation so that per-history quantities (currently pulse height) can be + // reassembled after the tree has been transported across separate Particle + // objects and, under MPI, across ranks. + int64_t root_index {-1}; double wgt_born {1.0}; double wgt_ww_born {-1.0}; int64_t n_split {0}; @@ -556,6 +562,8 @@ class ParticleData : public GeometryState { vector pht_storage_; + int64_t root_index_ {-1}; + double keff_tally_absorption_ {0.0}; double keff_tally_collision_ {0.0}; double keff_tally_tracklength_ {0.0}; @@ -736,6 +744,11 @@ class ParticleData : public GeometryState { // Interim pulse height tally storage vector& pht_storage() { return pht_storage_; } + const vector& pht_storage() const { return pht_storage_; } + + // Index of the primary particle at the root of this particle's tree + int64_t& root_index() { return root_index_; } + int64_t root_index() const { return root_index_; } // Global tally accumulators double& keff_tally_absorption() { return keff_tally_absorption_; } diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 454752cd271..b34ed3af099 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -49,6 +49,12 @@ extern const RegularMesh* ufs_mesh; extern vector k_generation; extern vector work_index; +//! Snapshot of work_index taken during phase 1 of shared-secondary transport, +//! i.e. the partition of *primary* particles across MPI ranks. work_index +//! itself is overwritten by calculate_work() on every secondary generation, so +//! it cannot be used afterwards to map a root index back to its owning rank. +extern vector phase1_work_index; + extern int64_t simulation_tracks_completed; //!< Number of tracks completed on this rank diff --git a/include/openmc/tallies/pulse_height.h b/include/openmc/tallies/pulse_height.h new file mode 100644 index 00000000000..3e574406e31 --- /dev/null +++ b/include/openmc/tallies/pulse_height.h @@ -0,0 +1,71 @@ +//! \file pulse_height.h +//! \brief Deferred, per-history aggregation of pulse-height results +//! +//! A pulse-height tally scores one count per source history, in the bin +//! containing the total energy that history's entire particle tree deposited in +//! a given cell. In the default transport modes a whole tree is carried by a +//! single Particle object, so Particle::pht_storage() already holds the +//! per-history total by the time event_death() runs and can be scored directly. +//! +//! Under the shared secondary bank each secondary generation is transported as +//! a fresh set of Particle objects, redistributed across MPI ranks between +//! generations. A history's deposition is therefore spread over many Particle +//! objects on potentially many ranks. This module collects those fragments, +//! keyed by the root index carried on every SourceSite, and scores them once +//! per history after the generation loop has drained. + +#ifndef OPENMC_TALLIES_PULSE_HEIGHT_H +#define OPENMC_TALLIES_PULSE_HEIGHT_H + +#include + +#include "openmc/vector.h" + +namespace openmc { + +//============================================================================== +//! One flushed pulse-height fragment, tagged with the history it belongs to. +//============================================================================== + +struct PulseHeightContribution { + int64_t root_index; //!< index of the primary at the root of the tree + vector energy; //!< per-cell energy, indexed as pulse_height_cells +}; + +namespace simulation { + +//! Per-thread staging buffers, merged in finalize_pulse_height_tallies(). +extern vector> pht_thread_buffers; + +} // namespace simulation + +//! Allocate the per-thread staging buffers. Called from initialize_simulation() +//! when pulse-height tallies and the shared secondary bank are both active. +void init_pulse_height_buffers(); + +//! Release the staging buffers and the phase-1 partition snapshot. +void free_memory_pulse_height(); + +//! Stage one Particle's contribution to its history's pulse height. +// +//! Thread-safe by construction: each thread appends only to its own buffer. +//! Contributions that are identically zero in every cell are dropped; histories +//! that deposit nothing are recovered in finalize_pulse_height_tallies() by +//! iterating over the full root range rather than over staged entries. +// +//! \param root_index index of the primary at the root of this particle's tree +//! \param pht per-cell energy deposited by this particle alone +void stage_pulse_height(int64_t root_index, const vector& pht); + +//! Aggregate staged contributions by history and score them. +// +//! Sends each contribution to the rank that owns its root according to +//! simulation::phase1_work_index, sums per (history, cell), and scores every +//! owned history including those with no deposition. Must be called after the +//! last secondary generation has been transported and before tally results are +//! accumulated for the batch. +void finalize_pulse_height_tallies(); + +} // namespace openmc + +#endif // OPENMC_TALLIES_PULSE_HEIGHT_H diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index 4303ddb7a90..8305023e3c0 100644 --- a/include/openmc/tallies/tally_scoring.h +++ b/include/openmc/tallies/tally_scoring.h @@ -120,7 +120,14 @@ void score_surface_tally( // //! \param p The particle being tracked //! \param tallies A vector of the indices of the tallies to score to -void score_pulse_height_tally(Particle& p, const vector& tallies); +//! Score a completed per-history pulse-height result. +// +//! \param p particle used to drive filter matching; its cell and E_last are +//! temporarily overwritten and restored +//! \param pht per-cell deposited energy, indexed as model::pulse_height_cells +//! \param tallies indices of the pulse-height tallies to score into +void score_pulse_height_tally( + Particle& p, const vector& pht, const vector& tallies); } // namespace openmc diff --git a/src/finalize.cpp b/src/finalize.cpp index fd891d9dd84..79e17e2e120 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -25,6 +25,7 @@ #include "openmc/simulation.h" #include "openmc/source.h" #include "openmc/surface.h" +#include "openmc/tallies/pulse_height.h" #include "openmc/tallies/tally.h" #include "openmc/thermal.h" #include "openmc/timer.h" @@ -51,6 +52,7 @@ void free_memory() free_memory_source(); free_memory_mesh(); free_memory_tally(); + free_memory_pulse_height(); free_memory_bank(); free_memory_plot(); free_memory_weight_windows(); diff --git a/src/initialize.cpp b/src/initialize.cpp index 33dfeca0e51..3c726675978 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -389,28 +389,6 @@ int parse_command_line(int argc, char* argv[]) return 0; } -// TODO: Pulse-height tallies require per-history scoring across the full -// particle tree (parent + all descendants). The shared secondary bank -// transports each secondary as an independent Particle, breaking this -// assumption. A proper fix would defer pulse-height scoring: save -// (root_source_id, cell, pht_storage) per particle, then aggregate by -// root_source_id after all secondary generations complete before scoring -// into the histogram. For now, disable shared secondary when pulse-height -// tallies are present. -static void check_pulse_height_compatibility() -{ - if (settings::use_shared_secondary_bank) { - for (const auto& t : model::tallies) { - if (t->type_ == TallyType::PULSE_HEIGHT) { - settings::use_shared_secondary_bank = false; - warning("Pulse-height tallies are not yet compatible with the shared " - "secondary bank. Disabling shared secondary bank."); - break; - } - } - } -} - bool read_model_xml() { std::string model_filename = settings::path_input; @@ -505,8 +483,6 @@ bool read_model_xml() if (check_for_node(root, "tallies")) read_tallies_xml(root.child("tallies")); - check_pulse_height_compatibility(); - // Initialize distribcell_filters prepare_distribcell(); @@ -552,8 +528,6 @@ void read_separate_xml_files() read_tallies_xml(); - check_pulse_height_compatibility(); - // Initialize distribcell_filters prepare_distribcell(); diff --git a/src/particle.cpp b/src/particle.cpp index f7a0098acf9..28d60b17592 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -29,6 +29,7 @@ #include "openmc/source.h" #include "openmc/surface.h" #include "openmc/tallies/derivative.h" +#include "openmc/tallies/pulse_height.h" #include "openmc/tallies/tally.h" #include "openmc/tallies/tally_scoring.h" #include "openmc/track_output.h" @@ -109,10 +110,29 @@ bool Particle::create_secondary( if (settings::use_shared_secondary_bank) { bank.progeny_id = n_progeny()++; } + bank.root_index = root_index(); bank.wgt_born = wgt_born(); bank.wgt_ww_born = wgt_ww_born(); bank.n_split = n_split(); + // Remove the energy carried off by this secondary from the parent's interim + // pulse-height result for the cell the parent is currently in. In non-shared + // mode the equivalent subtraction is performed at revival by + // pht_secondary_particles(); doing it here instead is equivalent, because the + // secondary is born at the parent's position and therefore in the parent's + // current cell, and it avoids the exhaustive_find_cell() call needed there. + // Placing this after the energy-cutoff early return above means a secondary + // that is never created is never subtracted, matching the non-shared path. + if (settings::use_shared_secondary_bank && + !model::active_pulse_height_tallies.empty() && type.is_photon()) { + auto it = std::find(model::pulse_height_cells.begin(), + model::pulse_height_cells.end(), lowest_coord().cell()); + if (it != model::pulse_height_cells.end()) { + int index = std::distance(model::pulse_height_cells.begin(), it); + pht_storage()[index] -= bank.E; + } + } + local_secondary_bank().emplace_back(bank); return true; } @@ -143,6 +163,10 @@ void Particle::split(double wgt) if (settings::use_shared_secondary_bank) { bank.progeny_id = n_progeny()++; } + // A split clone belongs to the same history as its parent. No pulse-height + // subtraction is applied here: a split is a weight artifact, not a physical + // secondary, and its energy is not carried away from the parent. + bank.root_index = root_index(); local_secondary_bank().emplace_back(bank); } @@ -502,6 +526,11 @@ void Particle::event_revive_from_secondary(const SourceSite& site) from_source(&site); + // Inherit the root of the tree this secondary belongs to. from_source() does + // not copy this, because it is also used for primaries read from the source + // bank, whose root index is assigned in initialize_particle_track(). + root_index() = site.root_index; + n_event() = 0; if (!settings::use_shared_secondary_bank) { n_tracks()++; @@ -509,8 +538,8 @@ void Particle::event_revive_from_secondary(const SourceSite& site) bank_second_E() = 0.0; // Subtract secondary particle energy from interim pulse-height results. - // In shared secondary mode, this subtraction was already done on the parent - // particle during create_secondary(), so skip it here. + // In shared secondary mode this subtraction is performed on the parent in + // create_secondary(), so skip it here. if (!settings::use_shared_secondary_bank && !model::active_pulse_height_tallies.empty() && this->type().is_photon()) { // Since the birth cell of the particle has not been set we @@ -604,7 +633,15 @@ void Particle::event_death() keff_tally_leakage() = 0.0; if (!model::active_pulse_height_tallies.empty()) { - score_pulse_height_tally(*this, model::active_pulse_height_tallies); + if (settings::use_shared_secondary_bank) { + // This Particle carries only one fragment of its history's pulse. Stage + // it for aggregation by root index; scoring happens once per history in + // finalize_pulse_height_tallies() after all generations have drained. + stage_pulse_height(root_index(), pht_storage()); + } else { + score_pulse_height_tally( + *this, pht_storage(), model::active_pulse_height_tallies); + } } // Accumulate track count for this particle history diff --git a/src/simulation.cpp b/src/simulation.cpp index 03f40a726eb..27fa4b761b5 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -23,6 +23,7 @@ #include "openmc/state_point.h" #include "openmc/tallies/derivative.h" #include "openmc/tallies/filter.h" +#include "openmc/tallies/pulse_height.h" #include "openmc/tallies/tally.h" #include "openmc/tallies/trigger.h" #include "openmc/timer.h" @@ -351,6 +352,7 @@ const RegularMesh* ufs_mesh {nullptr}; vector k_generation; vector work_index; +vector phase1_work_index; int64_t simulation_tracks_completed {0}; @@ -740,6 +742,17 @@ void initialize_particle_track( // Reset pulse_height_storage std::fill(p.pht_storage().begin(), p.pht_storage().end(), 0); + // A primary is the root of its own tree. Secondaries overwrite this in + // Particle::event_revive_from_secondary() using the value carried on the + // bank site. Only meaningful in shared-secondary mode, where a history is + // spread over several Particle objects; harmless otherwise. + if (!is_secondary) { + p.root_index() = simulation::phase1_work_index.empty() + ? index_source - 1 + : simulation::phase1_work_index[mpi::rank] + + index_source - 1; + } + // set random number seed int64_t particle_seed = compute_transport_seed(p.id()); init_particle_seeds(particle_seed, p.seeds()); @@ -1001,6 +1014,15 @@ void transport_history_based_shared_secondary() simulation::shared_secondary_bank_read.clear(); simulation::shared_secondary_bank_write.clear(); + // Record the primary partition before calculate_work() starts rewriting + // work_index for each secondary generation. finalize_pulse_height_tallies() + // needs it to map a root index back to the rank that owns that history. + simulation::phase1_work_index = simulation::work_index; + + if (!model::active_pulse_height_tallies.empty()) { + init_pulse_height_buffers(); + } + if (mpi::master) { write_message(fmt::format(" Primary source particles: {}", settings::n_particles), @@ -1099,6 +1121,12 @@ void transport_history_based_shared_secondary() simulation::simulation_tracks_completed += alive_secondary; } // End of loop over secondary generations + // The full particle tree of every history is now complete, so per-history + // pulse-height results can be reassembled and scored. + if (!model::active_pulse_height_tallies.empty()) { + finalize_pulse_height_tallies(); + } + // Reset work so that fission bank etc works correctly calculate_work(settings::n_particles); } @@ -1135,6 +1163,15 @@ void transport_event_based_shared_secondary() simulation::shared_secondary_bank_read.clear(); simulation::shared_secondary_bank_write.clear(); + // Record the primary partition before calculate_work() starts rewriting + // work_index for each secondary generation. finalize_pulse_height_tallies() + // needs it to map a root index back to the rank that owns that history. + simulation::phase1_work_index = simulation::work_index; + + if (!model::active_pulse_height_tallies.empty()) { + init_pulse_height_buffers(); + } + if (mpi::master) { write_message(fmt::format(" Primary source particles: {}", settings::n_particles), @@ -1231,6 +1268,12 @@ void transport_event_based_shared_secondary() simulation::simulation_tracks_completed += alive_secondary; } // End of loop over secondary generations + // The full particle tree of every history is now complete, so per-history + // pulse-height results can be reassembled and scored. + if (!model::active_pulse_height_tallies.empty()) { + finalize_pulse_height_tallies(); + } + // Reset work so that fission bank etc works correctly calculate_work(settings::n_particles); } diff --git a/src/tallies/pulse_height.cpp b/src/tallies/pulse_height.cpp new file mode 100644 index 00000000000..aa5b27862a5 --- /dev/null +++ b/src/tallies/pulse_height.cpp @@ -0,0 +1,223 @@ +#include "openmc/tallies/pulse_height.h" + +#include // upper_bound +#include + +#include "openmc/message_passing.h" +#include "openmc/openmp_interface.h" +#include "openmc/particle.h" +#include "openmc/settings.h" +#include "openmc/simulation.h" +#include "openmc/tallies/tally.h" +#include "openmc/tallies/tally_scoring.h" + +namespace openmc { + +//============================================================================== +// Global variables +//============================================================================== + +namespace simulation { + +vector> pht_thread_buffers; + +} // namespace simulation + +//============================================================================== +// Non-member functions +//============================================================================== + +void init_pulse_height_buffers() +{ + simulation::pht_thread_buffers.resize(num_threads()); + for (auto& buffer : simulation::pht_thread_buffers) { + buffer.clear(); + } +} + +void free_memory_pulse_height() +{ + simulation::pht_thread_buffers.clear(); + simulation::pht_thread_buffers.shrink_to_fit(); + simulation::phase1_work_index.clear(); + simulation::phase1_work_index.shrink_to_fit(); +} + +void stage_pulse_height(int64_t root_index, const vector& pht) +{ + // A particle whose root was never assigned cannot be attributed to a + // history. This should not happen, but dropping the fragment is safer than + // adding it to an arbitrary history. + if (root_index < 0) + return; + + // Defensive: staging is only reachable from the shared-secondary drivers, + // which call init_pulse_height_buffers() before transporting anything. + if (simulation::pht_thread_buffers.empty()) + return; + + // Histories that deposit nothing still have to be scored, but they are + // recovered from the full root range in finalize_pulse_height_tallies() + // rather than from staged entries, so an all-zero fragment carries no + // information and is not worth moving between ranks. + bool nonzero = false; + for (double e : pht) { + if (e != 0.0) { + nonzero = true; + break; + } + } + if (!nonzero) + return; + + PulseHeightContribution contribution; + contribution.root_index = root_index; + contribution.energy = pht; + simulation::pht_thread_buffers[thread_num()].push_back( + std::move(contribution)); +} + +namespace { + +//! Rank that owns a given root index, from the phase-1 primary partition. +int owner_of_root(int64_t root_index) +{ + const auto& index = simulation::phase1_work_index; + auto it = std::upper_bound(index.begin(), index.end(), root_index); + return static_cast(std::distance(index.begin(), it)) - 1; +} + +} // namespace + +void finalize_pulse_height_tallies() +{ + int n_cells = model::pulse_height_cells.size(); + if (n_cells == 0) + return; + + // Range of root indices owned by this rank + int64_t first_root = simulation::phase1_work_index[mpi::rank]; + int64_t last_root = simulation::phase1_work_index[mpi::rank + 1]; + int64_t n_owned = last_root - first_root; + + // Per-history, per-cell deposited energy for the histories owned here. + // Entries left at zero correspond to histories whose tree deposited nothing + // in any pulse-height cell; those are still scored below, matching the + // behaviour of the non-shared path where every primary is scored at death. + vector totals(n_owned * n_cells, 0.0); + + // Flatten the per-thread staging buffers, folding in everything already + // destined for this rank and packing the rest by destination. +#ifdef OPENMC_MPI + vector send_counts(mpi::n_procs, 0); + vector send_roots; + vector send_energy; + vector> roots_by_rank(mpi::n_procs); + vector> energy_by_rank(mpi::n_procs); +#endif + + for (auto& buffer : simulation::pht_thread_buffers) { + for (auto& contribution : buffer) { + int64_t root = contribution.root_index; + int owner = owner_of_root(root); + if (owner == mpi::rank) { + int64_t offset = (root - first_root) * n_cells; + for (int c = 0; c < n_cells; ++c) { + totals[offset + c] += contribution.energy[c]; + } + } else { +#ifdef OPENMC_MPI + roots_by_rank[owner].push_back(root); + energy_by_rank[owner].insert(energy_by_rank[owner].end(), + contribution.energy.begin(), contribution.energy.end()); + send_counts[owner]++; +#endif + } + } + buffer.clear(); + } + +#ifdef OPENMC_MPI + if (mpi::n_procs > 1) { + // Concatenate the per-destination buffers into contiguous send buffers + vector send_displs(mpi::n_procs, 0); + int total_send = 0; + for (int r = 0; r < mpi::n_procs; ++r) { + send_displs[r] = total_send; + total_send += send_counts[r]; + } + send_roots.reserve(total_send); + send_energy.reserve(static_cast(total_send) * n_cells); + for (int r = 0; r < mpi::n_procs; ++r) { + send_roots.insert( + send_roots.end(), roots_by_rank[r].begin(), roots_by_rank[r].end()); + send_energy.insert( + send_energy.end(), energy_by_rank[r].begin(), energy_by_rank[r].end()); + roots_by_rank[r].clear(); + roots_by_rank[r].shrink_to_fit(); + energy_by_rank[r].clear(); + energy_by_rank[r].shrink_to_fit(); + } + + // Exchange how many contributions each rank is sending to each other rank + vector recv_counts(mpi::n_procs, 0); + MPI_Alltoall(send_counts.data(), 1, MPI_INT, recv_counts.data(), 1, MPI_INT, + mpi::intracomm); + + vector recv_displs(mpi::n_procs, 0); + int total_recv = 0; + for (int r = 0; r < mpi::n_procs; ++r) { + recv_displs[r] = total_recv; + total_recv += recv_counts[r]; + } + + // Root indices, one per contribution + vector recv_roots(total_recv); + MPI_Alltoallv(send_roots.data(), send_counts.data(), send_displs.data(), + MPI_INT64_T, recv_roots.data(), recv_counts.data(), recv_displs.data(), + MPI_INT64_T, mpi::intracomm); + + // Energies, n_cells per contribution + vector send_counts_e(mpi::n_procs); + vector send_displs_e(mpi::n_procs); + vector recv_counts_e(mpi::n_procs); + vector recv_displs_e(mpi::n_procs); + for (int r = 0; r < mpi::n_procs; ++r) { + send_counts_e[r] = send_counts[r] * n_cells; + send_displs_e[r] = send_displs[r] * n_cells; + recv_counts_e[r] = recv_counts[r] * n_cells; + recv_displs_e[r] = recv_displs[r] * n_cells; + } + vector recv_energy(static_cast(total_recv) * n_cells); + MPI_Alltoallv(send_energy.data(), send_counts_e.data(), + send_displs_e.data(), MPI_DOUBLE, recv_energy.data(), + recv_counts_e.data(), recv_displs_e.data(), MPI_DOUBLE, mpi::intracomm); + + for (int i = 0; i < total_recv; ++i) { + int64_t offset = (recv_roots[i] - first_root) * n_cells; + for (int c = 0; c < n_cells; ++c) { + totals[offset + c] += recv_energy[static_cast(i) * n_cells + c]; + } + } + } +#endif + + // Score one pulse per owned history. score_pulse_height_tally() drives filter + // matching off a Particle, so give each thread a default-constructed one; its + // cell and E_last are overwritten and restored inside the call. +#pragma omp parallel + { + Particle p; + vector pht(n_cells); + +#pragma omp for schedule(static) + for (int64_t i = 0; i < n_owned; ++i) { + for (int c = 0; c < n_cells; ++c) { + pht[c] = totals[i * n_cells + c]; + } + score_pulse_height_tally(p, pht, model::active_pulse_height_tallies); + } + } +} + +} // namespace openmc diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index d17a62dce14..6c0a684b556 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2720,7 +2720,8 @@ void score_surface_tally( match.bins_present_ = false; } -void score_pulse_height_tally(Particle& p, const vector& tallies) +void score_pulse_height_tally( + Particle& p, const vector& pht, const vector& tallies) { // The pulse height tally in OpenMC hijacks the logic of CellFilter and // EnergyFilter to score specific quantities related to particle pulse height. @@ -2756,7 +2757,7 @@ void score_pulse_height_tally(Particle& p, const vector& tallies) int index = std::distance(model::pulse_height_cells.begin(), it); // Temporarily change energy of particle to pulse-height value - p.E_last() = p.pht_storage()[index]; + p.E_last() = pht[index]; // Initialize an iterator over valid filter bin combinations. If // there are no valid combinations, use a continue statement to ensure diff --git a/tests/unit_tests/test_pulse_height.py b/tests/unit_tests/test_pulse_height.py index 1f27cc6f264..3bf1ef10a8e 100644 --- a/tests/unit_tests/test_pulse_height.py +++ b/tests/unit_tests/test_pulse_height.py @@ -57,3 +57,206 @@ def test_pulse_height(model, run_in_tmpdir): np.testing.assert_array_equal(t1, t2[::-1]) +# --------------------------------------------------------------------------- +# Shared secondary bank +# +# A pulse-height tally scores exactly one count 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 lives in one Particle +# object, so that total is available directly at particle death. Under the +# shared secondary bank each secondary generation is transported as a fresh set +# of Particle objects, redistributed across MPI ranks between generations, so a +# history's deposition is spread over many Particle objects and has to be +# reassembled before scoring. +# +# The comparison between modes cannot be exact. compute_particle_id() and +# compute_transport_seed() both take a different branch when the shared bank is +# active, so the two modes sample different random number streams and produce +# different realizations of the same distribution. Only count conservation is +# exact; everything else is compared statistically. +# +# These run single-rank. The cross-rank aggregation in +# finalize_pulse_height_tallies() is only exercised under MPI. +# +# One trap when adding tests here: never place an energy filter edge at a +# deposition value a history can hit exactly, such as the source energy of a +# monoenergetic source in a fully absorbing detector. Floating point in the +# accumulated sum then decides which side of the edge each history falls on, +# roughly 6% land above it and are dropped from the tally, and what looks like +# a physics discrepancy is only a difference in rounding. +# --------------------------------------------------------------------------- + + +def _detector_model(particle, radius, energy_bounds, particles=1000, + batches=10, shared_secondary=False): + """NaI sphere in a void, with a pulse-height tally on the detector cell.""" + openmc.reset_auto_ids() + model = openmc.Model() + + NaI = openmc.Material() + NaI.set_density('g/cm3', 3.7) + NaI.add_element('Na', 1.0) + NaI.add_element('I', 1.0) + + detector_surf = openmc.Sphere(r=radius) + outer_surf = openmc.Sphere(r=radius + 1.0, boundary_type='vacuum') + detector = openmc.Cell(name='detector', fill=NaI, region=-detector_surf) + outside = openmc.Cell(name='outside', region=+detector_surf & -outer_surf) + model.geometry = openmc.Geometry([detector, outside]) + + model.settings.run_mode = 'fixed source' + model.settings.batches = batches + model.settings.particles = particles + model.settings.photon_transport = True + model.settings.shared_secondary_bank = shared_secondary + model.settings.source = openmc.IndependentSource( + energy=openmc.stats.delta_function(1e6), + particle=particle + ) + + tally = openmc.Tally(name='pht') + tally.scores = ['pulse-height'] + tally.filters = [ + openmc.CellFilter(detector), + openmc.EnergyFilter(energy_bounds), + ] + model.tallies = openmc.Tallies([tally]) + + return model + + +def _pulse_height(statepoint_path): + """Return the pulse-height spectrum and its per-bin standard deviation.""" + with openmc.StatePoint(statepoint_path) as sp: + tally = sp.get_tally(name='pht') + return tally.mean.ravel().copy(), tally.std_dev.ravel().copy() + + +def _assert_spectra_consistent(a, a_err, b, b_err, n_sigma=5.0): + """Compare two spectra bin by bin against their combined standard errors. + + Only bins holding at least 1% of histories are compared. Bins in the tail + carry few counts per batch, so their batch-to-batch spread is a poor + estimate of their true uncertainty and would drive spurious failures. + """ + sigma = np.hypot(a_err, b_err) + significant = (0.5 * (a + b) > 0.01) & (sigma > 0.0) + assert significant.any(), "no bins carry enough counts to compare" + z = np.abs(a[significant] - b[significant]) / sigma[significant] + assert z.max() < n_sigma, f"largest per-bin discrepancy is {z.max():.1f} sigma" + + +def _mean_deposition(spectrum, std_dev, bounds): + """Mean deposited energy per history, and its standard error.""" + centers = 0.5 * (bounds[:-1] + bounds[1:]) + mean = float(centers @ spectrum) + err = float(np.sqrt(np.sum((centers * std_dev) ** 2))) + return mean, err + + +@pytest.mark.parametrize('shared_secondary', [False, True]) +@pytest.mark.parametrize('particle', ['photon', 'neutron']) +def test_pulse_height_count_conservation(particle, shared_secondary, + run_in_tmpdir): + """Every history scores exactly once, in both transport modes. + + Fixed-source results are normalized by the source strength divided by the + number of source particles (Tally::accumulate), so summing a pulse-height + tally over all its energy bins gives the number of scores per source + particle, which must be exactly one. The energy filter is wide enough that + no history can deposit outside it. + + This is the direct test for scoring per track rather than per history: if + each secondary were scored separately the sum would exceed one by the mean + number of tracks per history. A history whose deposition is dropped instead + makes the sum fall short. + """ + # Neutron capture in iodine releases several MeV of prompt gammas, so + # deposition is not bounded by the source energy in the neutron case. + upper = 1.1e6 if particle == 'photon' else 20.0e6 + model = _detector_model( + particle, radius=1.0, energy_bounds=np.linspace(0.0, upper, 101), + shared_secondary=shared_secondary, + ) + + spectrum, _ = _pulse_height(model.run()) + + assert spectrum.sum() == pytest.approx(1.0, abs=1e-9) + + +@pytest.mark.parametrize('particle', ['photon', 'neutron']) +def test_shared_secondary_matches_local(particle, run_in_tmpdir): + """Both modes give the same pulse-height distribution in a thin detector. + + A 1 cm NaI sphere at 1 MeV produces short cascades, so this mostly exercises + the aggregation bookkeeping rather than deep secondary trees. Compared are + the mean deposited energy per history and the spectrum shape bin by bin. + """ + upper = 1.1e6 if particle == 'photon' else 20.0e6 + bounds = np.linspace(0.0, upper, 51) + + local = _detector_model(particle, radius=1.0, energy_bounds=bounds, + particles=2000, shared_secondary=False) + local_spec, local_err = _pulse_height(local.run()) + + shared = _detector_model(particle, radius=1.0, energy_bounds=bounds, + particles=2000, shared_secondary=True) + shared_spec, shared_err = _pulse_height(shared.run()) + + # Count conservation must hold in both before the shapes are compared + assert local_spec.sum() == pytest.approx(1.0, abs=1e-9) + assert shared_spec.sum() == pytest.approx(1.0, abs=1e-9) + + local_mean, local_mean_err = _mean_deposition(local_spec, local_err, bounds) + shared_mean, shared_mean_err = _mean_deposition( + shared_spec, shared_err, bounds) + mean_sigma = np.hypot(local_mean_err, shared_mean_err) + assert mean_sigma > 0.0 + assert abs(local_mean - shared_mean) < 4.0 * mean_sigma + + _assert_spectra_consistent(local_spec, local_err, shared_spec, shared_err) + + +def test_shared_secondary_matches_local_thick_detector(run_in_tmpdir): + """Both modes agree when the secondary cascade is deep. + + A 50 cm NaI sphere is tens of mean free paths thick at 1 MeV, so nothing + escapes and each history spawns roughly 1.7 secondaries through Compton + scattering, fluorescence and bremsstrahlung. That makes this far more + sensitive than the thin-detector case to the subtraction performed on the + parent in create_secondary() being matched by the descendants' own + contributions once they are reassembled. + + The top bin deliberately extends past the source energy. Pulse height + accumulates as a telescoping sum of E_last() - E() over a history, which + equals the source energy in exact arithmetic but not in floating point: + measured on this geometry, about 94% of histories sum to exactly 1 MeV or + just below and about 6% land a few ULPs above it. EnergyFilter matches on + E >= bins.front() && E <= bins.back(), so a top edge sitting exactly at + 1 MeV drops that 6% entirely. They are then absent from the tally, count + conservation silently breaks, and the comparison below degenerates into a + comparison of rounding behaviour between two different random number + streams rather than of spectra. + """ + # Top bin straddles the full-energy peak so histories whose floating point + # sum lands marginally above 1 MeV are still binned + bounds = np.concatenate([np.linspace(0.0, 0.99e6, 20), [1.1e6]]) + + local = _detector_model('photon', radius=50.0, energy_bounds=bounds, + shared_secondary=False) + local_spec, local_err = _pulse_height(local.run()) + + shared = _detector_model('photon', radius=50.0, energy_bounds=bounds, + shared_secondary=True) + shared_spec, shared_err = _pulse_height(shared.run()) + + assert local_spec.sum() == pytest.approx(1.0, abs=1e-9) + assert shared_spec.sum() == pytest.approx(1.0, abs=1e-9) + + # The full-energy peak holds the great majority of histories, so its + # fraction is the sharpest single statistic available here. + peak_sigma = np.hypot(local_err[-1], shared_err[-1]) + assert peak_sigma > 0.0 + assert abs(local_spec[-1] - shared_spec[-1]) < 4.0 * peak_sigma + + _assert_spectra_consistent(local_spec, local_err, shared_spec, shared_err) From 9757fba2771dd8cc42ec4e94e417c30b4922b9db Mon Sep 17 00:00:00 2001 From: GuySten Date: Mon, 31 Aug 2026 18:45:38 +0300 Subject: [PATCH 2/6] updates --- include/openmc/particle_data.h | 5 ----- include/openmc/simulation.h | 4 ---- src/simulation.cpp | 8 ++++---- 3 files changed, 4 insertions(+), 13 deletions(-) diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 20ee6f567c2..a096fc8a17e 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -52,11 +52,6 @@ struct SourceSite { int parent_nuclide {-1}; int64_t parent_id {0}; int64_t progeny_id {0}; - // Dense global index, in [0, n_particles), of the primary at the root of - // this particle's tree. Propagated unchanged through every secondary - // generation so that per-history quantities (currently pulse height) can be - // reassembled after the tree has been transported across separate Particle - // objects and, under MPI, across ranks. int64_t root_index {-1}; double wgt_born {1.0}; double wgt_ww_born {-1.0}; diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index b34ed3af099..7a17835a9a0 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -49,10 +49,6 @@ extern const RegularMesh* ufs_mesh; extern vector k_generation; extern vector work_index; -//! Snapshot of work_index taken during phase 1 of shared-secondary transport, -//! i.e. the partition of *primary* particles across MPI ranks. work_index -//! itself is overwritten by calculate_work() on every secondary generation, so -//! it cannot be used afterwards to map a root index back to its owning rank. extern vector phase1_work_index; extern int64_t diff --git a/src/simulation.cpp b/src/simulation.cpp index 27fa4b761b5..c076ffef1f1 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -747,10 +747,10 @@ void initialize_particle_track( // bank site. Only meaningful in shared-secondary mode, where a history is // spread over several Particle objects; harmless otherwise. if (!is_secondary) { - p.root_index() = simulation::phase1_work_index.empty() - ? index_source - 1 - : simulation::phase1_work_index[mpi::rank] + - index_source - 1; + p.root_index() = + simulation::phase1_work_index.empty() + ? index_source - 1 + : simulation::phase1_work_index[mpi::rank] + index_source - 1; } // set random number seed From 42ae42c320090bf6f00ce96a6514e872dcb575f0 Mon Sep 17 00:00:00 2001 From: GuySten Date: Mon, 31 Aug 2026 18:52:53 +0300 Subject: [PATCH 3/6] update regression tests --- .../shared_neutron/results_true.dat | 20 +++---- .../shared_photon/results_true.dat | 54 +++++++++---------- 2 files changed, 37 insertions(+), 37 deletions(-) diff --git a/tests/regression_tests/pulse_height/shared_neutron/results_true.dat b/tests/regression_tests/pulse_height/shared_neutron/results_true.dat index 7d28d5d442e..9082081516d 100644 --- a/tests/regression_tests/pulse_height/shared_neutron/results_true.dat +++ b/tests/regression_tests/pulse_height/shared_neutron/results_true.dat @@ -1,6 +1,6 @@ tally 1: -4.890000E+00 -4.784900E+00 +4.880000E+00 +4.765400E+00 0.000000E+00 0.000000E+00 0.000000E+00 @@ -9,8 +9,8 @@ tally 1: 0.000000E+00 1.000000E-02 1.000000E-04 -6.000000E-02 -1.800000E-03 +7.000000E-02 +1.900000E-03 0.000000E+00 0.000000E+00 0.000000E+00 @@ -19,8 +19,6 @@ tally 1: 1.000000E-04 0.000000E+00 0.000000E+00 -1.000000E-02 -1.000000E-04 0.000000E+00 0.000000E+00 0.000000E+00 @@ -29,8 +27,6 @@ tally 1: 0.000000E+00 0.000000E+00 0.000000E+00 -1.000000E-02 -1.000000E-04 0.000000E+00 0.000000E+00 0.000000E+00 @@ -41,10 +37,10 @@ tally 1: 0.000000E+00 0.000000E+00 0.000000E+00 -1.000000E-02 -1.000000E-04 0.000000E+00 0.000000E+00 +1.000000E-02 +1.000000E-04 0.000000E+00 0.000000E+00 0.000000E+00 @@ -89,6 +85,8 @@ tally 1: 0.000000E+00 0.000000E+00 0.000000E+00 +1.000000E-02 +1.000000E-04 0.000000E+00 0.000000E+00 0.000000E+00 @@ -149,6 +147,8 @@ tally 1: 0.000000E+00 0.000000E+00 0.000000E+00 +1.000000E-02 +1.000000E-04 0.000000E+00 0.000000E+00 0.000000E+00 diff --git a/tests/regression_tests/pulse_height/shared_photon/results_true.dat b/tests/regression_tests/pulse_height/shared_photon/results_true.dat index 795f9b4c513..a6bc603b531 100644 --- a/tests/regression_tests/pulse_height/shared_photon/results_true.dat +++ b/tests/regression_tests/pulse_height/shared_photon/results_true.dat @@ -1,6 +1,6 @@ tally 1: -4.120000E+00 -3.409000E+00 +4.110000E+00 +3.393900E+00 3.000000E-02 5.000000E-04 1.000000E-02 @@ -14,7 +14,7 @@ tally 1: 1.000000E-02 1.000000E-04 2.000000E-02 -2.000000E-04 +4.000000E-04 0.000000E+00 0.000000E+00 1.000000E-02 @@ -23,12 +23,12 @@ tally 1: 1.000000E-04 0.000000E+00 0.000000E+00 +3.000000E-02 +3.000000E-04 2.000000E-02 2.000000E-04 -2.000000E-02 -2.000000E-04 -2.000000E-02 -4.000000E-04 +4.000000E-02 +8.000000E-04 1.000000E-02 1.000000E-04 0.000000E+00 @@ -43,8 +43,8 @@ tally 1: 0.000000E+00 0.000000E+00 0.000000E+00 -1.000000E-02 -1.000000E-04 +0.000000E+00 +0.000000E+00 2.000000E-02 2.000000E-04 0.000000E+00 @@ -55,8 +55,8 @@ tally 1: 1.000000E-04 1.000000E-02 1.000000E-04 -0.000000E+00 -0.000000E+00 +1.000000E-02 +1.000000E-04 0.000000E+00 0.000000E+00 0.000000E+00 @@ -67,12 +67,12 @@ tally 1: 0.000000E+00 0.000000E+00 0.000000E+00 -1.000000E-02 -1.000000E-04 +2.000000E-02 +2.000000E-04 3.000000E-02 3.000000E-04 2.000000E-02 -2.000000E-04 +4.000000E-04 0.000000E+00 0.000000E+00 2.000000E-02 @@ -103,20 +103,20 @@ tally 1: 0.000000E+00 1.000000E-02 1.000000E-04 -2.000000E-02 -2.000000E-04 -0.000000E+00 -0.000000E+00 -3.000000E-02 -5.000000E-04 1.000000E-02 1.000000E-04 0.000000E+00 0.000000E+00 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 +1.000000E-02 +1.000000E-04 3.000000E-02 5.000000E-04 3.000000E-02 @@ -131,8 +131,8 @@ tally 1: 1.000000E-04 1.000000E-02 1.000000E-04 -3.000000E-02 -5.000000E-04 +2.000000E-02 +2.000000E-04 0.000000E+00 0.000000E+00 1.000000E-02 @@ -156,7 +156,7 @@ tally 1: 0.000000E+00 0.000000E+00 3.000000E-02 -3.000000E-04 +5.000000E-04 0.000000E+00 0.000000E+00 0.000000E+00 @@ -193,9 +193,9 @@ tally 1: 0.000000E+00 0.000000E+00 0.000000E+00 +1.000000E-02 +1.000000E-04 0.000000E+00 0.000000E+00 -0.000000E+00 -0.000000E+00 -1.700000E-01 -5.900000E-03 +1.600000E-01 +6.400000E-03 From 8489a5b1ff65cdc8156aa4dda2862293908b755c Mon Sep 17 00:00:00 2001 From: GuySten Date: Mon, 31 Aug 2026 19:10:38 +0300 Subject: [PATCH 4/6] update structs --- openmc/lib/core.py | 1 + src/initialize.cpp | 18 ++++++++++-------- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 22580d52a46..94523d19377 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -32,6 +32,7 @@ class _SourceSite(Structure): ('parent_nuclide', c_int), ('parent_id', c_int64), ('progeny_id', c_int64), + ('root_index', c_int64), ('wgt_born', c_double), ('wgt_ww_born', c_double), ('n_split', c_int64), diff --git a/src/initialize.cpp b/src/initialize.cpp index 3c726675978..f9f21f743f1 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -161,7 +161,7 @@ void initialize_mpi(MPI_Comm intracomm) // Create bank datatype SourceSite b; - MPI_Aint disp[15]; + MPI_Aint disp[16]; MPI_Get_address(&b.r, &disp[0]); MPI_Get_address(&b.u, &disp[1]); MPI_Get_address(&b.E, &disp[2]); @@ -173,16 +173,17 @@ void initialize_mpi(MPI_Comm intracomm) MPI_Get_address(&b.parent_nuclide, &disp[8]); MPI_Get_address(&b.parent_id, &disp[9]); MPI_Get_address(&b.progeny_id, &disp[10]); - MPI_Get_address(&b.wgt_born, &disp[11]); - MPI_Get_address(&b.wgt_ww_born, &disp[12]); - MPI_Get_address(&b.n_split, &disp[13]); - MPI_Get_address(&b.n_collision, &disp[14]); - for (int i = 14; i >= 0; --i) { + MPI_Get_address(&b.root_index, &disp[11]); + MPI_Get_address(&b.wgt_born, &disp[12]); + MPI_Get_address(&b.wgt_ww_born, &disp[13]); + MPI_Get_address(&b.n_split, &disp[14]); + MPI_Get_address(&b.n_collision, &disp[15]); + for (int i = 15; i >= 0; --i) { disp[i] -= disp[0]; } // Block counts for each field - int blocks[] = {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; + int blocks[] = {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; // Types for each field MPI_Datatype types[] = { @@ -197,13 +198,14 @@ void initialize_mpi(MPI_Comm intracomm) MPI_INT, // parent_nuclide MPI_INT64_T, // parent_id MPI_INT64_T, // progeny_id + MPI_INT64_T, // root_index MPI_DOUBLE, // wgt_born MPI_DOUBLE, // wgt_ww_born MPI_INT64_T, // n_split MPI_INT // n_collision }; - MPI_Type_create_struct(15, blocks, disp, types, &mpi::source_site); + MPI_Type_create_struct(16, blocks, disp, types, &mpi::source_site); MPI_Type_commit(&mpi::source_site); CollisionTrackSite bc; From 4852d6363f49abc8bdb08795962d5d19df4d46cc Mon Sep 17 00:00:00 2001 From: GuySten Date: Mon, 31 Aug 2026 19:41:27 +0300 Subject: [PATCH 5/6] update regression test --- .../shared/results_true.dat | 54 +++++++++---------- 1 file changed, 27 insertions(+), 27 deletions(-) diff --git a/tests/regression_tests/weightwindows_pulse_height/shared/results_true.dat b/tests/regression_tests/weightwindows_pulse_height/shared/results_true.dat index 795f9b4c513..a6bc603b531 100644 --- a/tests/regression_tests/weightwindows_pulse_height/shared/results_true.dat +++ b/tests/regression_tests/weightwindows_pulse_height/shared/results_true.dat @@ -1,6 +1,6 @@ tally 1: -4.120000E+00 -3.409000E+00 +4.110000E+00 +3.393900E+00 3.000000E-02 5.000000E-04 1.000000E-02 @@ -14,7 +14,7 @@ tally 1: 1.000000E-02 1.000000E-04 2.000000E-02 -2.000000E-04 +4.000000E-04 0.000000E+00 0.000000E+00 1.000000E-02 @@ -23,12 +23,12 @@ tally 1: 1.000000E-04 0.000000E+00 0.000000E+00 +3.000000E-02 +3.000000E-04 2.000000E-02 2.000000E-04 -2.000000E-02 -2.000000E-04 -2.000000E-02 -4.000000E-04 +4.000000E-02 +8.000000E-04 1.000000E-02 1.000000E-04 0.000000E+00 @@ -43,8 +43,8 @@ tally 1: 0.000000E+00 0.000000E+00 0.000000E+00 -1.000000E-02 -1.000000E-04 +0.000000E+00 +0.000000E+00 2.000000E-02 2.000000E-04 0.000000E+00 @@ -55,8 +55,8 @@ tally 1: 1.000000E-04 1.000000E-02 1.000000E-04 -0.000000E+00 -0.000000E+00 +1.000000E-02 +1.000000E-04 0.000000E+00 0.000000E+00 0.000000E+00 @@ -67,12 +67,12 @@ tally 1: 0.000000E+00 0.000000E+00 0.000000E+00 -1.000000E-02 -1.000000E-04 +2.000000E-02 +2.000000E-04 3.000000E-02 3.000000E-04 2.000000E-02 -2.000000E-04 +4.000000E-04 0.000000E+00 0.000000E+00 2.000000E-02 @@ -103,20 +103,20 @@ tally 1: 0.000000E+00 1.000000E-02 1.000000E-04 -2.000000E-02 -2.000000E-04 -0.000000E+00 -0.000000E+00 -3.000000E-02 -5.000000E-04 1.000000E-02 1.000000E-04 0.000000E+00 0.000000E+00 +1.000000E-02 +1.000000E-04 +1.000000E-02 +1.000000E-04 0.000000E+00 0.000000E+00 0.000000E+00 0.000000E+00 +1.000000E-02 +1.000000E-04 3.000000E-02 5.000000E-04 3.000000E-02 @@ -131,8 +131,8 @@ tally 1: 1.000000E-04 1.000000E-02 1.000000E-04 -3.000000E-02 -5.000000E-04 +2.000000E-02 +2.000000E-04 0.000000E+00 0.000000E+00 1.000000E-02 @@ -156,7 +156,7 @@ tally 1: 0.000000E+00 0.000000E+00 3.000000E-02 -3.000000E-04 +5.000000E-04 0.000000E+00 0.000000E+00 0.000000E+00 @@ -193,9 +193,9 @@ tally 1: 0.000000E+00 0.000000E+00 0.000000E+00 +1.000000E-02 +1.000000E-04 0.000000E+00 0.000000E+00 -0.000000E+00 -0.000000E+00 -1.700000E-01 -5.900000E-03 +1.600000E-01 +6.400000E-03 From 1643bd7eaf0ec59d5402aa850ea087a724fb405c Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 8 Sep 2026 20:26:59 +0300 Subject: [PATCH 6/6] simplify design significantly --- include/openmc/particle_data.h | 24 ++++- include/openmc/simulation.h | 33 ++++++- include/openmc/tallies/pulse_height.h | 12 ++- openmc/lib/core.py | 3 +- src/bank.cpp | 12 +-- src/initialize.cpp | 22 ++--- src/particle.cpp | 30 +++--- src/physics.cpp | 2 +- src/physics_mg.cpp | 2 +- src/simulation.cpp | 106 ++++++++++++++++----- src/tallies/pulse_height.cpp | 132 ++++++++++++++------------ tests/unit_tests/test_pulse_height.py | 93 ++++++------------ 12 files changed, 286 insertions(+), 185 deletions(-) diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index a096fc8a17e..859a43911ef 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -50,13 +50,33 @@ struct SourceSite { // Extra attributes that don't show up in source written to file int parent_nuclide {-1}; - int64_t parent_id {0}; + + //! Index of an ancestor of this site, in one of two senses that never + //! overlap. Before the site is placed by collection it is the immediate + //! parent's slot in the current generation's work, which with progeny_id + //! gives the site its position in the collected bank and is used for nothing + //! else. Placement consumes that key, and resolve_root_indices() then + //! overwrites the field with the index of the primary at the root of the + //! site's history, which is what survives the sort and the MPI migration the + //! shared secondary bank performs between generations. Read it through + //! parent_slot() or root_index() so the sense is explicit at every use. + int64_t ancestor_index {0}; + int64_t progeny_id {0}; - int64_t root_index {-1}; double wgt_born {1.0}; double wgt_ww_born {-1.0}; int64_t n_split {0}; int n_collision {0}; + + //! Slot of the immediate parent within the current generation's work. Valid + //! from creation until the site is placed by collection. + int64_t& parent_slot() { return ancestor_index; } + int64_t parent_slot() const { return ancestor_index; } + + //! Index of the primary at the root of this site's history. Valid once the + //! site has been placed by collection. + int64_t& root_index() { return ancestor_index; } + int64_t root_index() const { return ancestor_index; } }; struct CollisionTrackSite { diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 7a17835a9a0..99df7c56889 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -6,6 +6,7 @@ #include "openmc/mesh.h" #include "openmc/particle.h" +#include "openmc/shared_array.h" #include "openmc/vector.h" #include @@ -49,8 +50,6 @@ extern const RegularMesh* ufs_mesh; extern vector k_generation; extern vector work_index; -extern vector phase1_work_index; - extern int64_t simulation_tracks_completed; //!< Number of tracks completed on this rank @@ -66,6 +65,36 @@ void allocate_banks(); //! Determine number of particles to transport per process void calculate_work(int64_t n_particles); +//! First primary index owned by a rank under the phase-1 partition +//! +//! Recomputes what calculate_work(settings::n_particles) would produce, so +//! that the primary partition is available after work_index has been +//! overwritten for a secondary generation. Valid for rank in [0, n_procs], +//! where n_procs returns the total primary count. +//! +//! \param rank MPI rank +//! \return Index of that rank's first primary +int64_t phase1_first_root(int rank); + +//! Rank owning a given root index under the phase-1 partition +//! +//! \param root Root index in [0, n_particles) +//! \return Rank that transported that primary +int phase1_owner_of_root(int64_t root); + +//! Replace the placement key of every site in a collected bank with the root +//! of its history +//! +//! Must be called once per generation, after the sites have been placed and +//! before any MPI migration, since the placement key indexes a bank that is +//! local to the rank that produced the sites. +//! +//! \param sites Bank whose sites were just collected +//! \param parents Bank the parents were transported from, or nullptr when the +//! parents were the primaries +void resolve_root_indices( + SharedArray& sites, const SharedArray* parents); + //! Initialize nuclear data before a simulation void initialize_data(); diff --git a/include/openmc/tallies/pulse_height.h b/include/openmc/tallies/pulse_height.h index 3e574406e31..6e771e0b3b5 100644 --- a/include/openmc/tallies/pulse_height.h +++ b/include/openmc/tallies/pulse_height.h @@ -29,6 +29,7 @@ namespace openmc { struct PulseHeightContribution { int64_t root_index; //!< index of the primary at the root of the tree + int64_t track_id; //!< id of the track that deposited this fragment vector energy; //!< per-cell energy, indexed as pulse_height_cells }; @@ -54,16 +55,23 @@ void free_memory_pulse_height(); //! iterating over the full root range rather than over staged entries. // //! \param root_index index of the primary at the root of this particle's tree +//! \param track_id this particle's id, used to give the fragments of a history +//! a canonical summation order in finalize_pulse_height_tallies() //! \param pht per-cell energy deposited by this particle alone -void stage_pulse_height(int64_t root_index, const vector& pht); +void stage_pulse_height( + int64_t root_index, int64_t track_id, const vector& pht); //! Aggregate staged contributions by history and score them. // //! Sends each contribution to the rank that owns its root according to -//! simulation::phase1_work_index, sums per (history, cell), and scores every +//! the phase-1 primary partition, sums per (history, cell), and scores every //! owned history including those with no deposition. Must be called after the //! last secondary generation has been transported and before tally results are //! accumulated for the batch. +// +//! Fragments are summed in order of track id rather than in arrival order, so +//! a history's total is bit-for-bit independent of thread scheduling and of the +//! rank a descendant happened to land on. void finalize_pulse_height_tallies(); } // namespace openmc diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 94523d19377..9b6af17555c 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -30,9 +30,8 @@ class _SourceSite(Structure): ('surf_id', c_int), ('particle', c_int32), ('parent_nuclide', c_int), - ('parent_id', c_int64), + ('ancestor_index', c_int64), ('progeny_id', c_int64), - ('root_index', c_int64), ('wgt_born', c_double), ('wgt_ww_born', c_double), ('n_split', c_int64), diff --git a/src/bank.cpp b/src/bank.cpp index 5b12b48fd0e..9034ff4cef3 100644 --- a/src/bank.cpp +++ b/src/bank.cpp @@ -79,8 +79,8 @@ void init_fission_bank(int64_t max) } // Performs an O(n) sort on a fission or secondary bank, by leveraging -// the parent_id and progeny_id fields of banked particles. See the following -// paper for more details: +// the ancestor_index and progeny_id fields of banked particles. See the +// following paper for more details: // "Reproducibility and Monte Carlo Eigenvalue Calculations," F.B. Brown and // T.M. Sutton, 1992 ANS Annual Meeting, Transactions of the American Nuclear // Society, Volume 65, Page 235. @@ -121,15 +121,15 @@ void sort_bank(SharedArray& bank, bool is_fission_bank) // Use parent and progeny indices to sort bank for (int64_t i = 0; i < bank.size(); i++) { const auto& site = bank[i]; - if (site.parent_id < 0 || - site.parent_id >= + if (site.parent_slot() < 0 || + site.parent_slot() >= static_cast(simulation::progeny_per_particle.size())) { fatal_error(fmt::format("Invalid parent_id {} for banked site (expected " "range [0, {})).", - site.parent_id, simulation::progeny_per_particle.size())); + site.parent_slot(), simulation::progeny_per_particle.size())); } int64_t idx = - simulation::progeny_per_particle[site.parent_id] + site.progeny_id; + simulation::progeny_per_particle[site.parent_slot()] + site.progeny_id; if (idx < 0 || idx >= bank.size()) { fatal_error("Mismatch detected between sum of all particle progeny and " "bank size during sorting."); diff --git a/src/initialize.cpp b/src/initialize.cpp index f9f21f743f1..41d66b02c0b 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -161,7 +161,7 @@ void initialize_mpi(MPI_Comm intracomm) // Create bank datatype SourceSite b; - MPI_Aint disp[16]; + MPI_Aint disp[15]; MPI_Get_address(&b.r, &disp[0]); MPI_Get_address(&b.u, &disp[1]); MPI_Get_address(&b.E, &disp[2]); @@ -171,19 +171,18 @@ void initialize_mpi(MPI_Comm intracomm) MPI_Get_address(&b.surf_id, &disp[6]); MPI_Get_address(&b.particle, &disp[7]); MPI_Get_address(&b.parent_nuclide, &disp[8]); - MPI_Get_address(&b.parent_id, &disp[9]); + MPI_Get_address(&b.ancestor_index, &disp[9]); MPI_Get_address(&b.progeny_id, &disp[10]); - MPI_Get_address(&b.root_index, &disp[11]); - MPI_Get_address(&b.wgt_born, &disp[12]); - MPI_Get_address(&b.wgt_ww_born, &disp[13]); - MPI_Get_address(&b.n_split, &disp[14]); - MPI_Get_address(&b.n_collision, &disp[15]); - for (int i = 15; i >= 0; --i) { + MPI_Get_address(&b.wgt_born, &disp[11]); + MPI_Get_address(&b.wgt_ww_born, &disp[12]); + MPI_Get_address(&b.n_split, &disp[13]); + MPI_Get_address(&b.n_collision, &disp[14]); + for (int i = 14; i >= 0; --i) { disp[i] -= disp[0]; } // Block counts for each field - int blocks[] = {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; + int blocks[] = {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; // Types for each field MPI_Datatype types[] = { @@ -196,16 +195,15 @@ void initialize_mpi(MPI_Comm intracomm) MPI_INT, // surf_id MPI_INT, // particle (enum) MPI_INT, // parent_nuclide - MPI_INT64_T, // parent_id + MPI_INT64_T, // ancestor_index MPI_INT64_T, // progeny_id - MPI_INT64_T, // root_index MPI_DOUBLE, // wgt_born MPI_DOUBLE, // wgt_ww_born MPI_INT64_T, // n_split MPI_INT // n_collision }; - MPI_Type_create_struct(16, blocks, disp, types, &mpi::source_site); + MPI_Type_create_struct(15, blocks, disp, types, &mpi::source_site); MPI_Type_commit(&mpi::source_site); CollisionTrackSite bc; diff --git a/src/particle.cpp b/src/particle.cpp index 28d60b17592..54b1c9e2575 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -106,11 +106,10 @@ bool Particle::create_secondary( bank.E = settings::run_CE ? E : g(); bank.time = time(); bank_second_E() += bank.E; - bank.parent_id = current_work(); + bank.parent_slot() = current_work(); if (settings::use_shared_secondary_bank) { bank.progeny_id = n_progeny()++; } - bank.root_index = root_index(); bank.wgt_born = wgt_born(); bank.wgt_ww_born = wgt_ww_born(); bank.n_split = n_split(); @@ -159,14 +158,14 @@ void Particle::split(double wgt) bank.wgt_ww_born = wgt_ww_born(); bank.n_split = n_split(); bank.n_collision = n_collision(); - bank.parent_id = current_work(); + bank.parent_slot() = current_work(); if (settings::use_shared_secondary_bank) { bank.progeny_id = n_progeny()++; } - // A split clone belongs to the same history as its parent. No pulse-height + // A split clone belongs to the same history as its parent, which it inherits + // through the ordinary root resolution at collection. No pulse-height // subtraction is applied here: a split is a weight artifact, not a physical // secondary, and its energy is not carried away from the parent. - bank.root_index = root_index(); local_secondary_bank().emplace_back(bank); } @@ -526,10 +525,15 @@ void Particle::event_revive_from_secondary(const SourceSite& site) from_source(&site); - // Inherit the root of the tree this secondary belongs to. from_source() does - // not copy this, because it is also used for primaries read from the source - // bank, whose root index is assigned in initialize_particle_track(). - root_index() = site.root_index; + // Inherit the root of the tree this secondary belongs to. Only in shared + // secondary mode: there the site comes from a collected bank, whose sites + // carry a resolved root index, and the history is spread over several + // Particle objects. On the local path the site is still carrying its + // placement key, and the Particle already holds the correct root from + // initialize_particle_track(), since the whole tree is transported here. + if (settings::use_shared_secondary_bank) { + root_index() = site.root_index(); + } n_event() = 0; if (!settings::use_shared_secondary_bank) { @@ -637,7 +641,7 @@ void Particle::event_death() // This Particle carries only one fragment of its history's pulse. Stage // it for aggregation by root index; scoring happens once per history in // finalize_pulse_height_tallies() after all generations have drained. - stage_pulse_height(root_index(), pht_storage()); + stage_pulse_height(root_index(), id(), pht_storage()); } else { score_pulse_height_tally( *this, pht_storage(), model::active_pulse_height_tallies); @@ -1107,8 +1111,10 @@ void add_surf_source_to_bank(Particle& p, const Surface& surf) site.delayed_group = p.delayed_group(); site.surf_id = surf.id_; site.particle = p.type(); - site.parent_id = p.id(); - site.progeny_id = p.n_progeny(); + // ancestor_index and progeny_id are deliberately left alone. They exist to + // give a site its place when a bank is collected or sorted; a surface source + // site is never placed by either, so neither field applies to it, and nothing + // reads them here since the surface source file format does not carry them. int64_t idx = simulation::surf_source_bank.thread_safe_append(site); } diff --git a/src/physics.cpp b/src/physics.cpp index 4bf459b458f..1dda7887fd2 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -229,7 +229,7 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) } // Set parent and progeny IDs - site.parent_id = p.current_work(); + site.parent_slot() = p.current_work(); site.progeny_id = p.n_progeny()++; // Store fission site in bank diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 212ba765ca0..20b44f1cc45 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -179,7 +179,7 @@ void create_fission_sites(Particle& p) } // Set parent and progeny ID - site.parent_id = p.current_work(); + site.parent_slot() = p.current_work(); site.progeny_id = p.n_progeny()++; // Store fission site in bank diff --git a/src/simulation.cpp b/src/simulation.cpp index c076ffef1f1..17bcf57595c 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -352,7 +352,6 @@ const RegularMesh* ufs_mesh {nullptr}; vector k_generation; vector work_index; -vector phase1_work_index; int64_t simulation_tracks_completed {0}; @@ -396,15 +395,15 @@ void collect_sorted_history_secondary_banks( // Place each secondary according to its parent and progeny identifiers for (const auto& bank : thread_banks) { for (const auto& site : bank) { - if (site.parent_id < 0 || - site.parent_id >= + if (site.parent_slot() < 0 || + site.parent_slot() >= static_cast(simulation::progeny_per_particle.size())) { - fatal_error(fmt::format("Invalid parent_id {} for banked site " + fatal_error(fmt::format("Invalid parent slot {} for banked site " "(expected range [0, {})).", - site.parent_id, simulation::progeny_per_particle.size())); + site.parent_slot(), simulation::progeny_per_particle.size())); } int64_t idx = - simulation::progeny_per_particle[site.parent_id] + site.progeny_id; + simulation::progeny_per_particle[site.parent_slot()] + site.progeny_id; if (idx < 0 || idx >= n_progeny) { fatal_error("Mismatch detected between sum of all particle progeny and " "secondary bank size during collection."); @@ -747,10 +746,7 @@ void initialize_particle_track( // bank site. Only meaningful in shared-secondary mode, where a history is // spread over several Particle objects; harmless otherwise. if (!is_secondary) { - p.root_index() = - simulation::phase1_work_index.empty() - ? index_source - 1 - : simulation::phase1_work_index[mpi::rank] + index_source - 1; + p.root_index() = phase1_first_root(mpi::rank) + index_source - 1; } // set random number seed @@ -821,6 +817,67 @@ int64_t compute_transport_seed(int64_t particle_id) } } +int64_t phase1_first_root(int rank) +{ + // Reproduces the partition calculate_work(settings::n_particles) produces, + // without depending on the current contents of simulation::work_index, which + // is overwritten for every secondary generation. The partition is a pure + // function of the primary count and the number of ranks, so it can be + // recomputed wherever it is needed instead of being snapshotted. + int64_t min_work = settings::n_particles / mpi::n_procs; + int64_t remainder = settings::n_particles % mpi::n_procs; + return rank < remainder + ? static_cast(rank) * (min_work + 1) + : remainder * (min_work + 1) + + (static_cast(rank) - remainder) * min_work; +} + +int phase1_owner_of_root(int64_t root) +{ + int64_t min_work = settings::n_particles / mpi::n_procs; + int64_t remainder = settings::n_particles % mpi::n_procs; + + // Ranks below the remainder carry one extra primary each. Roots below the + // boundary fall in that region; the rest divide evenly. When min_work is + // zero the boundary equals n_particles, so the second branch, which would + // divide by zero, is unreachable. + int64_t boundary = remainder * (min_work + 1); + if (root < boundary) { + return static_cast(root / (min_work + 1)); + } + return static_cast(remainder + (root - boundary) / min_work); +} + +void resolve_root_indices( + SharedArray& sites, const SharedArray* parents) +{ + // Every site has now been placed, so its placement key has been consumed and + // the field can be overwritten with the root of its history. A site's parent + // is a primary when parents is null, and an entry of the generation just + // transported otherwise; in the latter case that entry already carries its + // own root, so the value simply propagates down the tree. + int64_t n = sites.size(); + int64_t n_parents = parents ? parents->size() : 0; + int64_t first_root = phase1_first_root(mpi::rank); + +#pragma omp parallel for schedule(static) + for (int64_t i = 0; i < n; ++i) { + int64_t slot = sites[i].parent_slot(); + if (parents) { + if (slot < 0 || slot >= n_parents) { + // fatal_error aborts the process, so it is safe to call from inside a + // parallel region + fatal_error(fmt::format("Invalid parent slot {} while resolving root " + "index (expected range [0, {})).", + slot, n_parents)); + } + sites[i].root_index() = (*parents)[slot].root_index(); + } else { + sites[i].root_index() = first_root + slot; + } + } +} + void calculate_work(int64_t n_particles) { // Determine minimum amount of particles to simulate on each processor @@ -1014,11 +1071,6 @@ void transport_history_based_shared_secondary() simulation::shared_secondary_bank_read.clear(); simulation::shared_secondary_bank_write.clear(); - // Record the primary partition before calculate_work() starts rewriting - // work_index for each secondary generation. finalize_pulse_height_tallies() - // needs it to map a root index back to the rank that owns that history. - simulation::phase1_work_index = simulation::work_index; - if (!model::active_pulse_height_tallies.empty()) { init_pulse_height_buffers(); } @@ -1053,6 +1105,7 @@ void transport_history_based_shared_secondary() } } collect_sorted_history_secondary_banks(thread_banks); + resolve_root_indices(simulation::shared_secondary_bank_write, nullptr); thread_banks.clear(); simulation::simulation_tracks_completed += settings::n_particles; @@ -1112,10 +1165,13 @@ void transport_history_based_shared_secondary() p.local_secondary_bank().clear(); } } // End of transport loop over tracks in shared secondary bank - simulation::shared_secondary_bank_write = - std::move(simulation::shared_secondary_bank_read); - simulation::shared_secondary_bank_read = SharedArray(); + // The bank just transported is needed to resolve the roots of the sites it + // produced, so it is released after collection rather than recycled into + // the write bank beforehand. collect_sorted_history_secondary_banks(thread_banks); + resolve_root_indices(simulation::shared_secondary_bank_write, + &simulation::shared_secondary_bank_read); + simulation::shared_secondary_bank_read = SharedArray(); thread_banks.clear(); n_generation_depth++; simulation::simulation_tracks_completed += alive_secondary; @@ -1163,11 +1219,6 @@ void transport_event_based_shared_secondary() simulation::shared_secondary_bank_read.clear(); simulation::shared_secondary_bank_write.clear(); - // Record the primary partition before calculate_work() starts rewriting - // work_index for each secondary generation. finalize_pulse_height_tallies() - // needs it to map a root index back to the rank that owns that history. - simulation::phase1_work_index = simulation::work_index; - if (!model::active_pulse_height_tallies.empty()) { init_pulse_height_buffers(); } @@ -1213,6 +1264,15 @@ void transport_event_based_shared_secondary() // ensure reproducibility. sort_bank(simulation::shared_secondary_bank_write, false); + // Roots are resolved after the sort, which consumes the placement key, and + // before the migration below, which invalidates the parent slots by moving + // sites away from the rank whose bank they index. On the first pass the + // read bank is empty because the parents were the primaries. + resolve_root_indices(simulation::shared_secondary_bank_write, + simulation::shared_secondary_bank_read.size() > 0 + ? &simulation::shared_secondary_bank_read + : nullptr); + // Synchronize the shared secondary bank amongst all MPI ranks, such // that each MPI rank has an approximately equal number of secondary // tracks. diff --git a/src/tallies/pulse_height.cpp b/src/tallies/pulse_height.cpp index aa5b27862a5..e43a3f0db3e 100644 --- a/src/tallies/pulse_height.cpp +++ b/src/tallies/pulse_height.cpp @@ -1,7 +1,9 @@ #include "openmc/tallies/pulse_height.h" -#include // upper_bound +#include // sort #include +#include // iota +#include // make_pair #include "openmc/message_passing.h" #include "openmc/openmp_interface.h" @@ -39,11 +41,10 @@ void free_memory_pulse_height() { simulation::pht_thread_buffers.clear(); simulation::pht_thread_buffers.shrink_to_fit(); - simulation::phase1_work_index.clear(); - simulation::phase1_work_index.shrink_to_fit(); } -void stage_pulse_height(int64_t root_index, const vector& pht) +void stage_pulse_height( + int64_t root_index, int64_t track_id, const vector& pht) { // A particle whose root was never assigned cannot be attributed to a // history. This should not happen, but dropping the fragment is safer than @@ -72,23 +73,12 @@ void stage_pulse_height(int64_t root_index, const vector& pht) PulseHeightContribution contribution; contribution.root_index = root_index; + contribution.track_id = track_id; contribution.energy = pht; simulation::pht_thread_buffers[thread_num()].push_back( std::move(contribution)); } -namespace { - -//! Rank that owns a given root index, from the phase-1 primary partition. -int owner_of_root(int64_t root_index) -{ - const auto& index = simulation::phase1_work_index; - auto it = std::upper_bound(index.begin(), index.end(), root_index); - return static_cast(std::distance(index.begin(), it)) - 1; -} - -} // namespace - void finalize_pulse_height_tallies() { int n_cells = model::pulse_height_cells.size(); @@ -96,8 +86,8 @@ void finalize_pulse_height_tallies() return; // Range of root indices owned by this rank - int64_t first_root = simulation::phase1_work_index[mpi::rank]; - int64_t last_root = simulation::phase1_work_index[mpi::rank + 1]; + int64_t first_root = phase1_first_root(mpi::rank); + int64_t last_root = phase1_first_root(mpi::rank + 1); int64_t n_owned = last_root - first_root; // Per-history, per-cell deposited energy for the histories owned here. @@ -106,28 +96,39 @@ void finalize_pulse_height_tallies() // behaviour of the non-shared path where every primary is scored at death. vector totals(n_owned * n_cells, 0.0); - // Flatten the per-thread staging buffers, folding in everything already - // destined for this rank and packing the rest by destination. + // Fragments are collected rather than summed as they are found, so that the + // summation order below can be fixed by track id instead of left to arrival + // order, which depends on thread scheduling and on where a descendant + // landed. A track id is the particle's global slot within its generation plus + // the tracks completed in earlier generations, both global quantities, so the + // resulting order is the same for any thread or rank count. + // Root and track id travel together as a pair, so one exchange carries both + vector own_keys; + vector own_energy; + + // Flatten the per-thread staging buffers, keeping everything already destined + // for this rank and packing the rest by destination. #ifdef OPENMC_MPI vector send_counts(mpi::n_procs, 0); - vector send_roots; + vector send_keys; vector send_energy; - vector> roots_by_rank(mpi::n_procs); + vector> keys_by_rank(mpi::n_procs); vector> energy_by_rank(mpi::n_procs); #endif for (auto& buffer : simulation::pht_thread_buffers) { for (auto& contribution : buffer) { int64_t root = contribution.root_index; - int owner = owner_of_root(root); + int owner = phase1_owner_of_root(root); if (owner == mpi::rank) { - int64_t offset = (root - first_root) * n_cells; - for (int c = 0; c < n_cells; ++c) { - totals[offset + c] += contribution.energy[c]; - } + own_keys.push_back(root); + own_keys.push_back(contribution.track_id); + own_energy.insert(own_energy.end(), contribution.energy.begin(), + contribution.energy.end()); } else { #ifdef OPENMC_MPI - roots_by_rank[owner].push_back(root); + keys_by_rank[owner].push_back(root); + keys_by_rank[owner].push_back(contribution.track_id); energy_by_rank[owner].insert(energy_by_rank[owner].end(), contribution.energy.begin(), contribution.energy.end()); send_counts[owner]++; @@ -146,15 +147,15 @@ void finalize_pulse_height_tallies() send_displs[r] = total_send; total_send += send_counts[r]; } - send_roots.reserve(total_send); + send_keys.reserve(static_cast(total_send) * 2); send_energy.reserve(static_cast(total_send) * n_cells); for (int r = 0; r < mpi::n_procs; ++r) { - send_roots.insert( - send_roots.end(), roots_by_rank[r].begin(), roots_by_rank[r].end()); + send_keys.insert( + send_keys.end(), keys_by_rank[r].begin(), keys_by_rank[r].end()); send_energy.insert( send_energy.end(), energy_by_rank[r].begin(), energy_by_rank[r].end()); - roots_by_rank[r].clear(); - roots_by_rank[r].shrink_to_fit(); + keys_by_rank[r].clear(); + keys_by_rank[r].shrink_to_fit(); energy_by_rank[r].clear(); energy_by_rank[r].shrink_to_fit(); } @@ -171,36 +172,49 @@ void finalize_pulse_height_tallies() total_recv += recv_counts[r]; } - // Root indices, one per contribution - vector recv_roots(total_recv); - MPI_Alltoallv(send_roots.data(), send_counts.data(), send_displs.data(), - MPI_INT64_T, recv_roots.data(), recv_counts.data(), recv_displs.data(), - MPI_INT64_T, mpi::intracomm); - - // Energies, n_cells per contribution - vector send_counts_e(mpi::n_procs); - vector send_displs_e(mpi::n_procs); - vector recv_counts_e(mpi::n_procs); - vector recv_displs_e(mpi::n_procs); - for (int r = 0; r < mpi::n_procs; ++r) { - send_counts_e[r] = send_counts[r] * n_cells; - send_displs_e[r] = send_displs[r] * n_cells; - recv_counts_e[r] = recv_counts[r] * n_cells; - recv_displs_e[r] = recv_displs[r] * n_cells; - } + // Both payloads are a fixed number of items per contribution, so their + // counts and displacements are the contribution ones scaled + auto scaled = [&](const vector& v, int factor) { + vector out(v.size()); + for (size_t i = 0; i < v.size(); ++i) + out[i] = v[i] * factor; + return out; + }; + + vector recv_keys(static_cast(total_recv) * 2); + MPI_Alltoallv(send_keys.data(), scaled(send_counts, 2).data(), + scaled(send_displs, 2).data(), MPI_INT64_T, recv_keys.data(), + scaled(recv_counts, 2).data(), scaled(recv_displs, 2).data(), MPI_INT64_T, + mpi::intracomm); + vector recv_energy(static_cast(total_recv) * n_cells); - MPI_Alltoallv(send_energy.data(), send_counts_e.data(), - send_displs_e.data(), MPI_DOUBLE, recv_energy.data(), - recv_counts_e.data(), recv_displs_e.data(), MPI_DOUBLE, mpi::intracomm); + MPI_Alltoallv(send_energy.data(), scaled(send_counts, n_cells).data(), + scaled(send_displs, n_cells).data(), MPI_DOUBLE, recv_energy.data(), + scaled(recv_counts, n_cells).data(), scaled(recv_displs, n_cells).data(), + MPI_DOUBLE, mpi::intracomm); + + // Received fragments join the local ones, to be ordered together below + own_keys.insert(own_keys.end(), recv_keys.begin(), recv_keys.end()); + own_energy.insert(own_energy.end(), recv_energy.begin(), recv_energy.end()); + } +#endif - for (int i = 0; i < total_recv; ++i) { - int64_t offset = (recv_roots[i] - first_root) * n_cells; - for (int c = 0; c < n_cells; ++c) { - totals[offset + c] += recv_energy[static_cast(i) * n_cells + c]; - } + // Sum each history's fragments in track id order. Sorting by root first is + // not needed for the result, since fragments of different histories land in + // disjoint accumulators, but it keeps the accumulation local in memory. + vector order(own_keys.size() / 2); + std::iota(order.begin(), order.end(), int64_t {0}); + std::sort(order.begin(), order.end(), [&](int64_t a, int64_t b) { + return std::make_pair(own_keys[2 * a], own_keys[2 * a + 1]) < + std::make_pair(own_keys[2 * b], own_keys[2 * b + 1]); + }); + + for (int64_t i : order) { + int64_t offset = (own_keys[2 * i] - first_root) * n_cells; + for (int c = 0; c < n_cells; ++c) { + totals[offset + c] += own_energy[static_cast(i) * n_cells + c]; } } -#endif // Score one pulse per owned history. score_pulse_height_tally() drives filter // matching off a Particle, so give each thread a default-constructed one; its diff --git a/tests/unit_tests/test_pulse_height.py b/tests/unit_tests/test_pulse_height.py index 3bf1ef10a8e..6993e28d1e2 100644 --- a/tests/unit_tests/test_pulse_height.py +++ b/tests/unit_tests/test_pulse_height.py @@ -184,29 +184,41 @@ def test_pulse_height_count_conservation(particle, shared_secondary, assert spectrum.sum() == pytest.approx(1.0, abs=1e-9) -@pytest.mark.parametrize('particle', ['photon', 'neutron']) -def test_shared_secondary_matches_local(particle, run_in_tmpdir): - """Both modes give the same pulse-height distribution in a thin detector. - - A 1 cm NaI sphere at 1 MeV produces short cascades, so this mostly exercises - the aggregation bookkeeping rather than deep secondary trees. Compared are - the mean deposited energy per history and the spectrum shape bin by bin. - """ - upper = 1.1e6 if particle == 'photon' else 20.0e6 - bounds = np.linspace(0.0, upper, 51) - - local = _detector_model(particle, radius=1.0, energy_bounds=bounds, - particles=2000, shared_secondary=False) - local_spec, local_err = _pulse_height(local.run()) - - shared = _detector_model(particle, radius=1.0, energy_bounds=bounds, - particles=2000, shared_secondary=True) - shared_spec, shared_err = _pulse_height(shared.run()) +# Thin detectors exercise the aggregation bookkeeping on short cascades. The +# thick one is tens of mean free paths at 1 MeV, so nothing escapes and each +# history spawns roughly 1.7 secondaries through Compton scattering, +# fluorescence and bremsstrahlung, which makes it far more sensitive to the +# subtraction performed on the parent in create_secondary() being matched by +# the descendants' own contributions once they are reassembled. Its top bin +# extends past the source energy for the reason given above. +_COMPARISON_CASES = [ + ('photon', 1.0, np.linspace(0.0, 1.1e6, 51), 2000), + ('neutron', 1.0, np.linspace(0.0, 20.0e6, 51), 2000), + ('photon', 50.0, np.concatenate([np.linspace(0.0, 0.99e6, 20), [1.1e6]]), + 1000), +] + + +@pytest.mark.parametrize('particle,radius,bounds,particles', _COMPARISON_CASES, + ids=['thin-photon', 'thin-neutron', 'thick-photon']) +def test_shared_secondary_matches_local(particle, radius, bounds, particles, + run_in_tmpdir): + """Both transport modes give the same pulse-height distribution.""" + spectra = {} + for shared in (False, True): + model = _detector_model(particle, radius=radius, energy_bounds=bounds, + particles=particles, shared_secondary=shared) + spectra[shared] = _pulse_height(model.run()) + + (local_spec, local_err), (shared_spec, shared_err) = \ + spectra[False], spectra[True] # Count conservation must hold in both before the shapes are compared assert local_spec.sum() == pytest.approx(1.0, abs=1e-9) assert shared_spec.sum() == pytest.approx(1.0, abs=1e-9) + # Mean deposited energy is the sharpest single statistic available, and + # catches a systematic shift that a per-bin comparison can absorb local_mean, local_mean_err = _mean_deposition(local_spec, local_err, bounds) shared_mean, shared_mean_err = _mean_deposition( shared_spec, shared_err, bounds) @@ -215,48 +227,3 @@ def test_shared_secondary_matches_local(particle, run_in_tmpdir): assert abs(local_mean - shared_mean) < 4.0 * mean_sigma _assert_spectra_consistent(local_spec, local_err, shared_spec, shared_err) - - -def test_shared_secondary_matches_local_thick_detector(run_in_tmpdir): - """Both modes agree when the secondary cascade is deep. - - A 50 cm NaI sphere is tens of mean free paths thick at 1 MeV, so nothing - escapes and each history spawns roughly 1.7 secondaries through Compton - scattering, fluorescence and bremsstrahlung. That makes this far more - sensitive than the thin-detector case to the subtraction performed on the - parent in create_secondary() being matched by the descendants' own - contributions once they are reassembled. - - The top bin deliberately extends past the source energy. Pulse height - accumulates as a telescoping sum of E_last() - E() over a history, which - equals the source energy in exact arithmetic but not in floating point: - measured on this geometry, about 94% of histories sum to exactly 1 MeV or - just below and about 6% land a few ULPs above it. EnergyFilter matches on - E >= bins.front() && E <= bins.back(), so a top edge sitting exactly at - 1 MeV drops that 6% entirely. They are then absent from the tally, count - conservation silently breaks, and the comparison below degenerates into a - comparison of rounding behaviour between two different random number - streams rather than of spectra. - """ - # Top bin straddles the full-energy peak so histories whose floating point - # sum lands marginally above 1 MeV are still binned - bounds = np.concatenate([np.linspace(0.0, 0.99e6, 20), [1.1e6]]) - - local = _detector_model('photon', radius=50.0, energy_bounds=bounds, - shared_secondary=False) - local_spec, local_err = _pulse_height(local.run()) - - shared = _detector_model('photon', radius=50.0, energy_bounds=bounds, - shared_secondary=True) - shared_spec, shared_err = _pulse_height(shared.run()) - - assert local_spec.sum() == pytest.approx(1.0, abs=1e-9) - assert shared_spec.sum() == pytest.approx(1.0, abs=1e-9) - - # The full-energy peak holds the great majority of histories, so its - # fraction is the sharpest single statistic available here. - peak_sigma = np.hypot(local_err[-1], shared_err[-1]) - assert peak_sigma > 0.0 - assert abs(local_spec[-1] - shared_spec[-1]) < 4.0 * peak_sigma - - _assert_spectra_consistent(local_spec, local_err, shared_spec, shared_err)