diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 5383487d4fa..174ef071373 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -478,6 +478,7 @@ class ParticleData : public GeometryState { CacheDataMG mg_xs_cache_; ParticleType type_ {ParticleType::neutron}; + ParticleType type_last_ {ParticleType::neutron}; double E_; double E_last_; @@ -575,6 +576,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 3beeb9d5ac5..8fc29aaf9e3 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_timed_tracklength_tallies; extern vector active_collision_tallies; diff --git a/include/openmc/tallies/tally_scoring.h b/include/openmc/tallies/tally_scoring.h index c3ab779e6a1..2efbfd94c42 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/particle.cpp b/src/particle.cpp index 6b4c332a8b2..30af7e9730a 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -38,6 +38,10 @@ namespace openmc { +namespace simulation { +thread_local Particle tmp_particle; +} + //============================================================================== // Particle implementation //============================================================================== @@ -88,6 +92,24 @@ 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 + 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( + simulation::tmp_particle, model::active_particleout_analog_tallies); + } else { + score_analog_tally_mg( + simulation::tmp_particle, model::active_particleout_analog_tallies); + } + } return true; } @@ -124,6 +146,7 @@ void Particle::from_source(const SourceSite* src) // Copy attributes from source bank site type() = src->particle; + type_last() = src->particle; wgt() = src->wgt; wgt_last() = src->wgt; r() = src->r; @@ -161,6 +184,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(); @@ -355,9 +379,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); } } 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/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 b9c615ecb5b..60d0ec00b66 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -58,6 +58,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_timed_tracklength_tallies; vector active_collision_tallies; @@ -157,7 +158,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); @@ -323,8 +325,7 @@ Tally::Tally(pugi::xml_node node) if (has_energyout && i_nuc == -1) { fatal_error(fmt::format( "Error on tally {}: Cannot use a " - "'nuclide_density' or 'temperature' derivative on a tally with " - "an " + "'nuclide_density' or 'temperature' derivative on a tally with an " "outgoing energy filter and 'total' nuclide rate. Instead, tally " "each nuclide in the material individually.", id_)); @@ -577,7 +578,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) @@ -586,6 +586,7 @@ void Tally::set_scores(const vector& scores) "outgoing energy filter"); break; + case SCORE_TOTAL: case SCORE_SCATTER: if (legendre_present) estimator_ = TallyEstimator::ANALOG; @@ -1109,6 +1110,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_timed_tracklength_tallies.clear(); model::active_collision_tallies.clear(); @@ -1125,12 +1127,15 @@ void setup_active_tallies() bool mesh_present = (tally.get_filter() || tally.get_filter()); auto time_filter = tally.get_filter(); + 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: if (time_filter && mesh_present) { @@ -1173,6 +1178,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_timed_tracklength_tallies.clear(); model::active_collision_tallies.clear(); diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 67e851644a9..1a02760dc1c 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -2311,7 +2311,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 @@ -2321,7 +2321,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 @@ -2363,9 +2363,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