From fd474019699b420c533d87fb4c04114db78f1941 Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 26 Aug 2025 03:44:15 -0700 Subject: [PATCH 1/5] initial implementation wip --- include/openmc/particle.h | 3 +- include/openmc/particle_data.h | 3 ++ include/openmc/tallies/filter.h | 1 + include/openmc/tallies/filter_particle.h | 40 ++++++++++++++ include/openmc/tallies/tally.h | 1 + include/openmc/tallies/tally_scoring.h | 6 ++- openmc/filter.py | 29 ++++++++-- openmc/lib/filter.py | 10 ++-- src/mesh.cpp | 2 +- src/particle.cpp | 21 ++++++-- src/particle_restart.cpp | 1 + src/random_ray/random_ray.cpp | 2 +- src/simulation.cpp | 5 +- src/tallies/filter.cpp | 2 + src/tallies/filter_particle.cpp | 69 +++++++++++++++++++++++- src/tallies/tally.cpp | 11 +++- src/tallies/tally_scoring.cpp | 8 +-- 17 files changed, 190 insertions(+), 24 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index 0f37719b94f..cdfbfbb6b82 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -62,7 +62,8 @@ class Particle : public ParticleData { //! site may have been produced from an external source, from fission, or //! simply as a secondary particle. //! \param src Source site data - void from_source(const SourceSite* src); + //! \param particle ParticleType incident particle type + void from_source(const SourceSite* src, ParticleType particle); // Coarse-grained particle events void event_calculate_xs(); diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index ec393845f81..235fbf59ac7 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -469,6 +469,7 @@ class ParticleData : public GeometryState { CacheDataMG mg_xs_cache_; ParticleType type_ {ParticleType::neutron}; + ParticleType type_last_ {ParticleType::neutron}; double E_; double E_last_; @@ -566,6 +567,8 @@ class ParticleData : public GeometryState { // Particle type (n, p, e, gamma, etc) ParticleType& type() { return type_; } const ParticleType& type() const { return type_; } + ParticleType& type_last() { return type_last_; } + const ParticleType& type_last() const { return type_last_; } // Current particle energy, energy before collision, // and corresponding multigroup group indices. Energy diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h index 65098597a5d..6fd9d3bf8d5 100644 --- a/include/openmc/tallies/filter.h +++ b/include/openmc/tallies/filter.h @@ -39,6 +39,7 @@ enum class FilterType { MUSURFACE, PARENT_NUCLIDE, PARTICLE, + PARTICLE_OUT, POLAR, SPHERICAL_HARMONICS, SPATIAL_LEGENDRE, diff --git a/include/openmc/tallies/filter_particle.h b/include/openmc/tallies/filter_particle.h index 863a6d282fe..fff955ee5f2 100644 --- a/include/openmc/tallies/filter_particle.h +++ b/include/openmc/tallies/filter_particle.h @@ -48,5 +48,45 @@ class ParticleFilter : public Filter { vector particles_; }; +//============================================================================== +//! Bins by type of outgoing particle (e.g. neutron, photon). +//============================================================================== + +class ParticleOutFilter : public Filter { +public: + //---------------------------------------------------------------------------- + // Constructors, destructors + + ~ParticleOutFilter() = default; + + //---------------------------------------------------------------------------- + // Methods + + std::string type_str() const override { return "particleout"; } + FilterType type() const override { return FilterType::PARTICLE_OUT; } + + void from_xml(pugi::xml_node node) override; + + void get_all_bins(const Particle& p, TallyEstimator estimator, + FilterMatch& match) const override; + + void to_statepoint(hid_t filter_group) const override; + + std::string text_label(int bin) const override; + + //---------------------------------------------------------------------------- + // Accessors + + const vector& particles() const { return particles_; } + + void set_particles(span particles); + +private: + //---------------------------------------------------------------------------- + // Data members + + vector particles_; +}; + } // namespace openmc #endif // OPENMC_TALLIES_FILTER_PARTICLE_H diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 8c088c46089..4d4ed4f4735 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -202,6 +202,7 @@ extern std::unordered_map tally_map; extern vector> tallies; extern vector active_tallies; extern vector active_analog_tallies; +extern vector active_particleout_analog_tallies; extern vector active_tracklength_tallies; extern vector active_collision_tallies; extern vector active_meshsurf_tallies; diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index 28f1f16223d..043ea6d03d6 100644 --- a/include/openmc/tallies/tally_scoring.h +++ b/include/openmc/tallies/tally_scoring.h @@ -72,14 +72,16 @@ void score_collision_tally(Particle& p); //! Analog tallies are triggered at every collision, not every event. // //! \param p The particle being tracked -void score_analog_tally_ce(Particle& p); +//! \param tallies A vector of the indices of the tallies to score to +void score_analog_tally_ce(Particle& p, const vector& tallies); //! Score tallies based on a simple count of events (for multigroup). // //! Analog tallies are triggered at every collision, not every event. // //! \param p The particle being tracked -void score_analog_tally_mg(Particle& p); +//! \param tallies A vector of the indices of the tallies to score to +void score_analog_tally_mg(Particle& p, const vector& tallies); //! Score tallies using a tracklength estimate of the flux. // diff --git a/openmc/filter.py b/openmc/filter.py index 6a666d2a02c..ae039cbcef0 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -24,9 +24,9 @@ 'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', 'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre', - 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance', - 'collision', 'time', 'parentnuclide', 'weight', 'meshborn', 'meshsurface', - 'meshmaterial', + 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'particleout', + 'cellinstance', 'collision', 'time', 'parentnuclide', 'weight', 'meshborn', + 'meshsurface', 'meshmaterial', ) _CURRENT_NAMES = ( @@ -785,6 +785,29 @@ def from_xml_element(cls, elem, **kwargs): filter_id = int(get_text(elem, "id")) bins = get_elem_list(elem, "bins", str) or [] return cls(bins, filter_id=filter_id) + + +class ParticleoutFilter(ParticleFilter): + """Bins tally events based on the outgoing particle type. + + Parameters + ---------- + bins : str, or sequence of str + The particles to tally represented as strings ('neutron', 'photon', + 'electron', 'positron'). + filter_id : int + Unique identifier for the filter + + Attributes + ---------- + bins : sequence of str + The particles to tally + id : int + Unique identifier for the filter + num_bins : Integral + The number of filter bins + + """ class ParentNuclideFilter(ParticleFilter): diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py index 55b681d89ad..3ed03ec590a 100644 --- a/openmc/lib/filter.py +++ b/openmc/lib/filter.py @@ -22,9 +22,9 @@ 'EnergyFilter', 'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter', 'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshBornFilter', 'MeshMaterialFilter', 'MeshSurfaceFilter', 'MuFilter', 'MuSurfaceFilter', - 'ParentNuclideFilter', 'ParticleFilter', 'PolarFilter', 'SphericalHarmonicsFilter', - 'SpatialLegendreFilter', 'SurfaceFilter', 'TimeFilter', 'UniverseFilter', - 'WeightFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters' + 'ParentNuclideFilter', 'ParticleFilter', 'ParticleoutFilter', 'PolarFilter', + 'SphericalHarmonicsFilter', 'SpatialLegendreFilter', 'SurfaceFilter', 'TimeFilter', + 'UniverseFilter', 'WeightFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters' ] # Tally functions @@ -564,6 +564,10 @@ def bins(self): return [ParticleType(i) for i in particle_i] +class ParticleoutFilter(ParticleFilter): + filter_type = 'particleout' + + class PolarFilter(Filter): filter_type = 'polar' diff --git a/src/mesh.cpp b/src/mesh.cpp index 3e4ff1a3eca..97243e0d7d4 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -358,7 +358,7 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, site.r[ax1] = min1 + (i1 + 0.5) * d1; site.r[ax2] = min2 + (i2 + 0.5) * d2; - p.from_source(&site); + p.from_source(&site, site.particle); // Determine particle's location if (!exhaustive_find_cell(p)) { diff --git a/src/particle.cpp b/src/particle.cpp index fb41d82e950..6dfb6c5ff95 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -110,7 +110,7 @@ void Particle::split(double wgt) } } -void Particle::from_source(const SourceSite* src) +void Particle::from_source(const SourceSite* src, ParticleType particle) { // Reset some attributes clear(); @@ -141,6 +141,7 @@ void Particle::from_source(const SourceSite* src) E() = data::mg.energy_bin_avg_[g()]; } E_last() = E(); + type_last() = particle; time() = src->time; time_last() = src->time; parent_nuclide() = src->parent_nuclide; @@ -160,6 +161,7 @@ void Particle::event_calculate_xs() // Store pre-collision particle properties wgt_last() = wgt(); E_last() = E(); + type_last() = type(); u_last() = u(); r_last() = r(); time_last() = time(); @@ -348,9 +350,9 @@ void Particle::event_collide() score_collision_tally(*this); if (!model::active_analog_tallies.empty()) { if (settings::run_CE) { - score_analog_tally_ce(*this); + score_analog_tally_ce(*this, model::active_analog_tallies); } else { - score_analog_tally_mg(*this); + score_analog_tally_mg(*this, model::active_analog_tallies); } } @@ -408,6 +410,8 @@ void Particle::event_revive_from_secondary() wgt() = 0.0; } + auto type = this->type(); + // Check for secondary particles if this particle is dead if (!alive()) { // Write final position for this particle @@ -419,11 +423,20 @@ void Particle::event_revive_from_secondary() if (secondary_bank().empty()) return; - from_source(&secondary_bank().back()); + from_source(&secondary_bank().back(), type); secondary_bank().pop_back(); n_event() = 0; bank_second_E() = 0.0; + // Score tallies affected by secondary particles + if (!model::active_particleout_analog_tallies.empty()) { + if (settings::run_CE) { + score_analog_tally_ce(*this, model::active_particleout_analog_tallies); + } else { + score_analog_tally_mg(*this, model::active_particleout_analog_tallies); + } + } + // Subtract secondary particle energy from interim pulse-height results if (!model::active_pulse_height_tallies.empty() && this->type() == ParticleType::photon) { diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp index 6b7778211af..ea4070b951c 100644 --- a/src/particle_restart.cpp +++ b/src/particle_restart.cpp @@ -64,6 +64,7 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode) p.r_last() = p.r(); p.u_last() = p.u(); p.E_last() = p.E(); + p.type_last() = p.type(); p.g_last() = p.g(); p.time_last() = p.time(); diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index 27f674c4278..5e59acb3691 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -796,7 +796,7 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain) } site.E = 0.0; - this->from_source(&site); + this->from_source(&site, site.particle); // Locate ray if (lowest_coord().cell() == C_NONE) { diff --git a/src/simulation.cpp b/src/simulation.cpp index b9986f47375..a0adbaac9f9 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -567,7 +567,8 @@ void initialize_history(Particle& p, int64_t index_source) // set defaults if (settings::run_mode == RunMode::EIGENVALUE) { // set defaults for eigenvalue simulations from primary bank - p.from_source(&simulation::source_bank[index_source - 1]); + p.from_source( + &simulation::source_bank[index_source - 1], ParticleType::neutron); } else if (settings::run_mode == RunMode::FIXED_SOURCE) { // initialize random number seed int64_t id = (simulation::total_gen + overall_generation() - 1) * @@ -576,7 +577,7 @@ void initialize_history(Particle& p, int64_t index_source) uint64_t seed = init_seed(id, STREAM_SOURCE); // sample from external source distribution or custom library then set auto site = sample_external_source(&seed); - p.from_source(&site); + p.from_source(&site, site.particle); } p.current_work() = index_source; diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index 79817981db0..af23718b24b 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -146,6 +146,8 @@ Filter* Filter::create(const std::string& type, int32_t id) return Filter::create(id); } else if (type == "particle") { return Filter::create(id); + } else if (type == "particleout") { + return Filter::create(id); } else if (type == "polar") { return Filter::create(id); } else if (type == "surface") { diff --git a/src/tallies/filter_particle.cpp b/src/tallies/filter_particle.cpp index eef1d1e63c5..80295cded8c 100644 --- a/src/tallies/filter_particle.cpp +++ b/src/tallies/filter_particle.cpp @@ -6,6 +6,10 @@ namespace openmc { +//============================================================================== +// ParticleFilter implementation +//============================================================================== + void ParticleFilter::from_xml(pugi::xml_node node) { auto particles = get_node_array(node, "bins"); @@ -34,8 +38,11 @@ void ParticleFilter::set_particles(span particles) void ParticleFilter::get_all_bins( const Particle& p, TallyEstimator estimator, FilterMatch& match) const { + // Get the pre-collision type of the particle. + auto type = p.type_last(); + for (auto i = 0; i < particles_.size(); i++) { - if (particles_[i] == p.type()) { + if (particles_[i] == type) { match.bins_.push_back(i); match.weights_.push_back(1.0); } @@ -58,6 +65,66 @@ std::string ParticleFilter::text_label(int bin) const return fmt::format("Particle: {}", particle_type_to_str(p)); } +//============================================================================== +// ParticleOutFilter implementation +//============================================================================== + +void ParticleOutFilter::from_xml(pugi::xml_node node) +{ + auto particles = get_node_array(node, "bins"); + + // Convert to vector of ParticleType + vector types; + for (auto& p : particles) { + types.push_back(str_to_particle_type(p)); + } + this->set_particles(types); +} + +void ParticleOutFilter::set_particles(span particles) +{ + // Clear existing particles + particles_.clear(); + particles_.reserve(particles.size()); + + // Set particles and number of bins + for (auto p : particles) { + particles_.push_back(p); + } + n_bins_ = particles_.size(); +} + +void ParticleOutFilter::get_all_bins( + const Particle& p, TallyEstimator estimator, FilterMatch& match) const +{ + for (auto i = 0; i < particles_.size(); i++) { + if (particles_[i] == p.type()) { + match.bins_.push_back(i); + match.weights_.push_back(1.0); + } + } +} + +void ParticleOutFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + vector particles; + for (auto p : particles_) { + particles.push_back(particle_type_to_str(p)); + } + write_dataset(filter_group, "bins", particles); +} + +std::string ParticleOutFilter::text_label(int bin) const +{ + const auto& p = particles_.at(bin); + return fmt::format("ParticleOut: {}", particle_type_to_str(p)); +} + +//============================================================================== +// C-API functions +//============================================================================== + extern "C" int openmc_particle_filter_get_bins(int32_t idx, int bins[]) { if (int err = verify_filter(idx)) diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index ba183860269..cf6760bda19 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -55,6 +55,7 @@ std::unordered_map tally_map; vector> tallies; vector active_tallies; vector active_analog_tallies; +vector active_particleout_analog_tallies; vector active_tracklength_tallies; vector active_collision_tallies; vector active_meshsurf_tallies; @@ -152,7 +153,8 @@ Tally::Tally(pugi::xml_node node) // Change the tally estimator if a filter demands it FilterType filt_type = f->type(); if (filt_type == FilterType::ENERGY_OUT || - filt_type == FilterType::LEGENDRE) { + filt_type == FilterType::LEGENDRE || + filt_type == FilterType::PARTICLE_OUT) { estimator_ = TallyEstimator::ANALOG; } else if (filt_type == FilterType::SPHERICAL_HARMONICS) { auto sf = dynamic_cast(f); @@ -569,7 +571,6 @@ void Tally::set_scores(const vector& scores) fatal_error("Cannot tally flux with an outgoing energy filter."); break; - case SCORE_TOTAL: case SCORE_ABSORPTION: case SCORE_FISSION: if (energyout_present) @@ -578,6 +579,7 @@ void Tally::set_scores(const vector& scores) "outgoing energy filter"); break; + case SCORE_TOTAL: case SCORE_SCATTER: if (legendre_present) estimator_ = TallyEstimator::ANALOG; @@ -1068,6 +1070,7 @@ void setup_active_tallies() { model::active_tallies.clear(); model::active_analog_tallies.clear(); + model::active_particleout_analog_tallies.clear(); model::active_tracklength_tallies.clear(); model::active_collision_tallies.clear(); model::active_meshsurf_tallies.clear(); @@ -1079,12 +1082,15 @@ void setup_active_tallies() if (tally.active_) { model::active_tallies.push_back(i); + bool particleout_present = tally.has_filter(FilterType::PARTICLE_OUT); switch (tally.type_) { case TallyType::VOLUME: switch (tally.estimator_) { case TallyEstimator::ANALOG: model::active_analog_tallies.push_back(i); + if (particleout_present) + model::active_particleout_analog_tallies.push_back(i); break; case TallyEstimator::TRACKLENGTH: model::active_tracklength_tallies.push_back(i); @@ -1122,6 +1128,7 @@ void free_memory_tally() model::active_tallies.clear(); model::active_analog_tallies.clear(); + model::active_particleout_analog_tallies.clear(); model::active_tracklength_tallies.clear(); model::active_collision_tallies.clear(); model::active_meshsurf_tallies.clear(); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index e73fb90f311..a1220f2139a 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2303,7 +2303,7 @@ void score_general_mg(Particle& p, int i_tally, int start_index, } } -void score_analog_tally_ce(Particle& p) +void score_analog_tally_ce(Particle& p, const vector& tallies) { // Since electrons/positrons are not transported, we assign a flux of zero. // Note that the heating score does NOT use the flux and will be non-zero for @@ -2313,7 +2313,7 @@ void score_analog_tally_ce(Particle& p) ? 1.0 : 0.0; - for (auto i_tally : model::active_analog_tallies) { + for (auto i_tally : tallies) { const Tally& tally {*model::tallies[i_tally]}; // Initialize an iterator over valid filter bin combinations. If there are @@ -2355,9 +2355,9 @@ void score_analog_tally_ce(Particle& p) match.bins_present_ = false; } -void score_analog_tally_mg(Particle& p) +void score_analog_tally_mg(Particle& p, const vector& tallies) { - for (auto i_tally : model::active_analog_tallies) { + for (auto i_tally : tallies) { const Tally& tally {*model::tallies[i_tally]}; // Initialize an iterator over valid filter bin combinations. If there are From e60d78dd92e1cc6fc4bec981f102541f070b477f Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 26 Aug 2025 21:21:51 +0300 Subject: [PATCH 2/5] store progeny energy in SourceSite E_last --- include/openmc/particle_data.h | 1 + openmc/lib/core.py | 1 + src/initialize.cpp | 28 +++++++++++++++------------- src/mesh.cpp | 1 + src/particle.cpp | 2 +- src/physics.cpp | 1 + src/physics_mg.cpp | 1 + src/random_ray/random_ray.cpp | 1 + 8 files changed, 22 insertions(+), 14 deletions(-) diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 235fbf59ac7..109631a7b5a 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -44,6 +44,7 @@ struct SourceSite { Position r; Direction u; double E; + double E_last; double time {0.0}; double wgt {1.0}; int delayed_group {0}; diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 9f8db69d577..f65742d8d63 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -22,6 +22,7 @@ class _SourceSite(Structure): _fields_ = [('r', c_double*3), ('u', c_double*3), ('E', c_double), + ('E_last', c_double), ('time', c_double), ('wgt', c_double), ('delayed_group', c_int), diff --git a/src/initialize.cpp b/src/initialize.cpp index 2da78137d10..bab50aa8c09 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -157,26 +157,28 @@ void initialize_mpi(MPI_Comm intracomm) // Create bank datatype SourceSite b; - MPI_Aint disp[11]; + MPI_Aint disp[12]; MPI_Get_address(&b.r, &disp[0]); MPI_Get_address(&b.u, &disp[1]); MPI_Get_address(&b.E, &disp[2]); - MPI_Get_address(&b.time, &disp[3]); - MPI_Get_address(&b.wgt, &disp[4]); - MPI_Get_address(&b.delayed_group, &disp[5]); - 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.progeny_id, &disp[10]); - for (int i = 10; i >= 0; --i) { + MPI_Get_address(&b.E_last, &disp[3]); + MPI_Get_address(&b.time, &disp[4]); + MPI_Get_address(&b.wgt, &disp[5]); + MPI_Get_address(&b.delayed_group, &disp[6]); + MPI_Get_address(&b.surf_id, &disp[7]); + MPI_Get_address(&b.particle, &disp[8]); + MPI_Get_address(&b.parent_nuclide, &disp[9]); + MPI_Get_address(&b.parent_id, &disp[10]); + MPI_Get_address(&b.progeny_id, &disp[11]); + for (int i = 11; i >= 0; --i) { disp[i] -= disp[0]; } - int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1}; + int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, - MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG}; - MPI_Type_create_struct(11, blocks, disp, types, &mpi::source_site); + MPI_DOUBLE, MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, + MPI_LONG}; + MPI_Type_create_struct(12, blocks, disp, types, &mpi::source_site); MPI_Type_commit(&mpi::source_site); } #endif // OPENMC_MPI diff --git a/src/mesh.cpp b/src/mesh.cpp index 97243e0d7d4..4586177f739 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -321,6 +321,7 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, SourceSite site; site.E = 1.0; + site.E_last = 1.0; site.particle = ParticleType::neutron; for (int axis = 0; axis < 3; ++axis) { diff --git a/src/particle.cpp b/src/particle.cpp index 6dfb6c5ff95..1ebfbf2fc40 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -132,6 +132,7 @@ void Particle::from_source(const SourceSite* src, ParticleType particle) r_last_current() = src->r; r_last() = src->r; u_last() = src->u; + E_last() = src->E_last; if (settings::run_CE) { E() = src->E; g() = 0; @@ -140,7 +141,6 @@ void Particle::from_source(const SourceSite* src, ParticleType particle) g_last() = static_cast(src->E); E() = data::mg.energy_bin_avg_[g()]; } - E_last() = E(); type_last() = particle; time() = src->time; time_last() = src->time; diff --git a/src/physics.cpp b/src/physics.cpp index 3c06e543de1..1d463174988 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -209,6 +209,7 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) site.r = p.r(); site.particle = ParticleType::neutron; site.time = p.time(); + site.E_last = p.E(); site.wgt = 1. / weight; site.parent_id = p.id(); site.progeny_id = p.n_progeny()++; diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 5ebef9c1417..6445ccea356 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -138,6 +138,7 @@ void create_fission_sites(Particle& p) site.r = p.r(); site.particle = ParticleType::neutron; site.time = p.time(); + site.E_last = p.E(); site.wgt = 1. / weight; site.parent_id = p.id(); site.progeny_id = p.n_progeny()++; diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index 5e59acb3691..a3f88345622 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -796,6 +796,7 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain) } site.E = 0.0; + site.E_last = 0.0; this->from_source(&site, site.particle); // Locate ray From fac1c90d84c27eecceefb0e09f3d721d7019012b Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 26 Aug 2025 22:06:13 +0300 Subject: [PATCH 3/5] move E_last to be later in SourceSite structs --- include/openmc/particle_data.h | 2 +- openmc/lib/core.py | 2 +- src/initialize.cpp | 14 +++++++------- 3 files changed, 9 insertions(+), 9 deletions(-) diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 109631a7b5a..c45edadf700 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -44,7 +44,6 @@ struct SourceSite { Position r; Direction u; double E; - double E_last; double time {0.0}; double wgt {1.0}; int delayed_group {0}; @@ -52,6 +51,7 @@ struct SourceSite { ParticleType particle; // Extra attributes that don't show up in source written to file + double E_last; int parent_nuclide {-1}; int64_t parent_id; int64_t progeny_id; diff --git a/openmc/lib/core.py b/openmc/lib/core.py index f65742d8d63..68ab289aa1e 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -22,12 +22,12 @@ class _SourceSite(Structure): _fields_ = [('r', c_double*3), ('u', c_double*3), ('E', c_double), - ('E_last', c_double), ('time', c_double), ('wgt', c_double), ('delayed_group', c_int), ('surf_id', c_int), ('particle', c_int), + ('E_last', c_double), ('parent_nuclide', c_int), ('parent_id', c_int64), ('progeny_id', c_int64)] diff --git a/src/initialize.cpp b/src/initialize.cpp index bab50aa8c09..6d28f0bf84b 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -161,12 +161,12 @@ void initialize_mpi(MPI_Comm intracomm) MPI_Get_address(&b.r, &disp[0]); MPI_Get_address(&b.u, &disp[1]); MPI_Get_address(&b.E, &disp[2]); - MPI_Get_address(&b.E_last, &disp[3]); - MPI_Get_address(&b.time, &disp[4]); - MPI_Get_address(&b.wgt, &disp[5]); - MPI_Get_address(&b.delayed_group, &disp[6]); - MPI_Get_address(&b.surf_id, &disp[7]); - MPI_Get_address(&b.particle, &disp[8]); + MPI_Get_address(&b.time, &disp[3]); + MPI_Get_address(&b.wgt, &disp[4]); + MPI_Get_address(&b.delayed_group, &disp[5]); + MPI_Get_address(&b.surf_id, &disp[6]); + MPI_Get_address(&b.particle, &disp[7]); + MPI_Get_address(&b.E_last, &disp[8]); MPI_Get_address(&b.parent_nuclide, &disp[9]); MPI_Get_address(&b.parent_id, &disp[10]); MPI_Get_address(&b.progeny_id, &disp[11]); @@ -176,7 +176,7 @@ void initialize_mpi(MPI_Comm intracomm) int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, - MPI_DOUBLE, MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, + MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_DOUBLE, MPI_INT, MPI_LONG, MPI_LONG}; MPI_Type_create_struct(12, blocks, disp, types, &mpi::source_site); MPI_Type_commit(&mpi::source_site); From b3020d1a30692fe97d0da8aec498af260b0d69dd Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 26 Aug 2025 22:46:09 +0300 Subject: [PATCH 4/5] try new approach --- include/openmc/particle.h | 3 +-- include/openmc/particle_data.h | 1 - openmc/lib/core.py | 1 - src/initialize.cpp | 18 ++++++++--------- src/mesh.cpp | 3 +-- src/particle.cpp | 36 ++++++++++++++++++++-------------- src/physics.cpp | 1 - src/physics_mg.cpp | 1 - src/random_ray/random_ray.cpp | 3 +-- src/simulation.cpp | 5 ++--- 10 files changed, 34 insertions(+), 38 deletions(-) diff --git a/include/openmc/particle.h b/include/openmc/particle.h index cdfbfbb6b82..0f37719b94f 100644 --- a/include/openmc/particle.h +++ b/include/openmc/particle.h @@ -62,8 +62,7 @@ class Particle : public ParticleData { //! site may have been produced from an external source, from fission, or //! simply as a secondary particle. //! \param src Source site data - //! \param particle ParticleType incident particle type - void from_source(const SourceSite* src, ParticleType particle); + void from_source(const SourceSite* src); // Coarse-grained particle events void event_calculate_xs(); diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index c45edadf700..235fbf59ac7 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -51,7 +51,6 @@ struct SourceSite { ParticleType particle; // Extra attributes that don't show up in source written to file - double E_last; int parent_nuclide {-1}; int64_t parent_id; int64_t progeny_id; diff --git a/openmc/lib/core.py b/openmc/lib/core.py index 68ab289aa1e..9f8db69d577 100644 --- a/openmc/lib/core.py +++ b/openmc/lib/core.py @@ -27,7 +27,6 @@ class _SourceSite(Structure): ('delayed_group', c_int), ('surf_id', c_int), ('particle', c_int), - ('E_last', c_double), ('parent_nuclide', c_int), ('parent_id', c_int64), ('progeny_id', c_int64)] diff --git a/src/initialize.cpp b/src/initialize.cpp index 6d28f0bf84b..2da78137d10 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -157,7 +157,7 @@ void initialize_mpi(MPI_Comm intracomm) // Create bank datatype SourceSite b; - MPI_Aint disp[12]; + MPI_Aint disp[11]; MPI_Get_address(&b.r, &disp[0]); MPI_Get_address(&b.u, &disp[1]); MPI_Get_address(&b.E, &disp[2]); @@ -166,19 +166,17 @@ void initialize_mpi(MPI_Comm intracomm) MPI_Get_address(&b.delayed_group, &disp[5]); MPI_Get_address(&b.surf_id, &disp[6]); MPI_Get_address(&b.particle, &disp[7]); - MPI_Get_address(&b.E_last, &disp[8]); - MPI_Get_address(&b.parent_nuclide, &disp[9]); - MPI_Get_address(&b.parent_id, &disp[10]); - MPI_Get_address(&b.progeny_id, &disp[11]); - for (int i = 11; i >= 0; --i) { + MPI_Get_address(&b.parent_nuclide, &disp[8]); + MPI_Get_address(&b.parent_id, &disp[9]); + MPI_Get_address(&b.progeny_id, &disp[10]); + for (int i = 10; i >= 0; --i) { disp[i] -= disp[0]; } - int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}; + int blocks[] {3, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1}; MPI_Datatype types[] {MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, MPI_DOUBLE, - MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_DOUBLE, MPI_INT, MPI_LONG, - MPI_LONG}; - MPI_Type_create_struct(12, blocks, disp, types, &mpi::source_site); + MPI_DOUBLE, MPI_INT, MPI_INT, MPI_INT, MPI_INT, MPI_LONG, MPI_LONG}; + MPI_Type_create_struct(11, blocks, disp, types, &mpi::source_site); MPI_Type_commit(&mpi::source_site); } #endif // OPENMC_MPI diff --git a/src/mesh.cpp b/src/mesh.cpp index 4586177f739..3e4ff1a3eca 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -321,7 +321,6 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, SourceSite site; site.E = 1.0; - site.E_last = 1.0; site.particle = ParticleType::neutron; for (int axis = 0; axis < 3; ++axis) { @@ -359,7 +358,7 @@ void Mesh::material_volumes(int nx, int ny, int nz, int table_size, site.r[ax1] = min1 + (i1 + 0.5) * d1; site.r[ax2] = min2 + (i2 + 0.5) * d2; - p.from_source(&site, site.particle); + p.from_source(&site); // Determine particle's location if (!exhaustive_find_cell(p)) { diff --git a/src/particle.cpp b/src/particle.cpp index 1ebfbf2fc40..f6002bfd23f 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -88,6 +88,23 @@ bool Particle::create_secondary( bank.E = settings::run_CE ? E : g(); bank.time = time(); bank_second_E() += bank.E; + + // Score tallies affected by secondary particles + if (!model::active_particleout_analog_tallies.empty()) { + // Create secondary particle for tallying purposes only + Particle p; + p.from_source(&bank); + p.u_last() = this->u(); + p.r_last() = this->r(); + p.E_last() = this->E(); + p.type_last() = this->type(); + + if (settings::run_CE) { + score_analog_tally_ce(p, model::active_particleout_analog_tallies); + } else { + score_analog_tally_mg(p, model::active_particleout_analog_tallies); + } + } return true; } @@ -110,7 +127,7 @@ void Particle::split(double wgt) } } -void Particle::from_source(const SourceSite* src, ParticleType particle) +void Particle::from_source(const SourceSite* src) { // Reset some attributes clear(); @@ -124,6 +141,7 @@ void Particle::from_source(const SourceSite* src, ParticleType particle) // Copy attributes from source bank site type() = src->particle; + type_last() = src->particle; wgt() = src->wgt; wgt_last() = src->wgt; r() = src->r; @@ -132,7 +150,6 @@ void Particle::from_source(const SourceSite* src, ParticleType particle) r_last_current() = src->r; r_last() = src->r; u_last() = src->u; - E_last() = src->E_last; if (settings::run_CE) { E() = src->E; g() = 0; @@ -141,7 +158,7 @@ void Particle::from_source(const SourceSite* src, ParticleType particle) g_last() = static_cast(src->E); E() = data::mg.energy_bin_avg_[g()]; } - type_last() = particle; + E_last() = E(); time() = src->time; time_last() = src->time; parent_nuclide() = src->parent_nuclide; @@ -410,8 +427,6 @@ void Particle::event_revive_from_secondary() wgt() = 0.0; } - auto type = this->type(); - // Check for secondary particles if this particle is dead if (!alive()) { // Write final position for this particle @@ -423,20 +438,11 @@ void Particle::event_revive_from_secondary() if (secondary_bank().empty()) return; - from_source(&secondary_bank().back(), type); + from_source(&secondary_bank().back()); secondary_bank().pop_back(); n_event() = 0; bank_second_E() = 0.0; - // Score tallies affected by secondary particles - if (!model::active_particleout_analog_tallies.empty()) { - if (settings::run_CE) { - score_analog_tally_ce(*this, model::active_particleout_analog_tallies); - } else { - score_analog_tally_mg(*this, model::active_particleout_analog_tallies); - } - } - // Subtract secondary particle energy from interim pulse-height results if (!model::active_pulse_height_tallies.empty() && this->type() == ParticleType::photon) { diff --git a/src/physics.cpp b/src/physics.cpp index 1d463174988..3c06e543de1 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -209,7 +209,6 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) site.r = p.r(); site.particle = ParticleType::neutron; site.time = p.time(); - site.E_last = p.E(); site.wgt = 1. / weight; site.parent_id = p.id(); site.progeny_id = p.n_progeny()++; diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 6445ccea356..5ebef9c1417 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -138,7 +138,6 @@ void create_fission_sites(Particle& p) site.r = p.r(); site.particle = ParticleType::neutron; site.time = p.time(); - site.E_last = p.E(); site.wgt = 1. / weight; site.parent_id = p.id(); site.progeny_id = p.n_progeny()++; diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index a3f88345622..27f674c4278 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -796,8 +796,7 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain) } site.E = 0.0; - site.E_last = 0.0; - this->from_source(&site, site.particle); + this->from_source(&site); // Locate ray if (lowest_coord().cell() == C_NONE) { diff --git a/src/simulation.cpp b/src/simulation.cpp index a0adbaac9f9..b9986f47375 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -567,8 +567,7 @@ void initialize_history(Particle& p, int64_t index_source) // set defaults if (settings::run_mode == RunMode::EIGENVALUE) { // set defaults for eigenvalue simulations from primary bank - p.from_source( - &simulation::source_bank[index_source - 1], ParticleType::neutron); + p.from_source(&simulation::source_bank[index_source - 1]); } else if (settings::run_mode == RunMode::FIXED_SOURCE) { // initialize random number seed int64_t id = (simulation::total_gen + overall_generation() - 1) * @@ -577,7 +576,7 @@ void initialize_history(Particle& p, int64_t index_source) uint64_t seed = init_seed(id, STREAM_SOURCE); // sample from external source distribution or custom library then set auto site = sample_external_source(&seed); - p.from_source(&site, site.particle); + p.from_source(&site); } p.current_work() = index_source; From 4e47f041b7968f572418a0121e7613af7f499049 Mon Sep 17 00:00:00 2001 From: GuySten Date: Wed, 27 Aug 2025 18:30:10 +0300 Subject: [PATCH 5/5] pre allocate tmp parricle --- src/particle.cpp | 21 +++++++++++++-------- 1 file changed, 13 insertions(+), 8 deletions(-) diff --git a/src/particle.cpp b/src/particle.cpp index f6002bfd23f..9d723bba858 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -38,6 +38,10 @@ namespace openmc { +namespace simulation { +thread_local Particle tmp_particle; +} + //============================================================================== // Particle implementation //============================================================================== @@ -92,17 +96,18 @@ bool Particle::create_secondary( // Score tallies affected by secondary particles if (!model::active_particleout_analog_tallies.empty()) { // Create secondary particle for tallying purposes only - Particle p; - p.from_source(&bank); - p.u_last() = this->u(); - p.r_last() = this->r(); - p.E_last() = this->E(); - p.type_last() = this->type(); + simulation::tmp_particle.from_source(&bank); + simulation::tmp_particle.u_last() = this->u(); + simulation::tmp_particle.r_last() = this->r(); + simulation::tmp_particle.E_last() = this->E(); + simulation::tmp_particle.type_last() = this->type(); if (settings::run_CE) { - score_analog_tally_ce(p, model::active_particleout_analog_tallies); + score_analog_tally_ce( + simulation::tmp_particle, model::active_particleout_analog_tallies); } else { - score_analog_tally_mg(p, model::active_particleout_analog_tallies); + score_analog_tally_mg( + simulation::tmp_particle, model::active_particleout_analog_tallies); } } return true;