From f2c936cf5b5187c8bd575fee8689e5fcee1a8166 Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Thu, 5 Feb 2026 18:00:03 -0600 Subject: [PATCH] Implement filter for secondary particle production binned by energy (#3453) Co-authored-by: Paul Romano --- docs/source/pythonapi/base.rst | 1 + include/openmc/constants.h | 2 +- include/openmc/particle_data.h | 18 +++ include/openmc/tallies/filter.h | 1 + include/openmc/tallies/filter_energy.h | 32 +++++ openmc/filter.py | 115 +++++++++++++++--- src/particle.cpp | 9 ++ src/physics.cpp | 2 + src/physics_mg.cpp | 2 + src/tallies/filter.cpp | 2 + src/tallies/filter_delayedgroup.cpp | 6 +- src/tallies/filter_energy.cpp | 44 +++++++ src/tallies/tally.cpp | 2 + .../photon_production_inputs_true.dat | 15 +++ .../photon_production/inputs_true.dat | 15 +++ .../photon_production/results_true.dat | 11 ++ .../photon_production/test.py | 30 ++++- tests/unit_tests/test_filters.py | 47 ++++++- 18 files changed, 328 insertions(+), 26 deletions(-) diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index dea8c4427..7bd03f379 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -132,6 +132,7 @@ Constructing Tallies openmc.MeshSurfaceFilter openmc.EnergyFilter openmc.EnergyoutFilter + openmc.ParticleProductionFilter openmc.MuFilter openmc.MuSurfaceFilter openmc.PolarFilter diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 980cff447..8185214fc 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -301,7 +301,7 @@ enum class TallyEstimator { ANALOG, TRACKLENGTH, COLLISION }; enum class TallyEvent { SURFACE, LATTICE, KILL, SCATTER, ABSORB }; // Tally score type -- if you change these, make sure you also update the -// _SCORES dictionary in openmc/capi/tally.py +// _SCORES dictionary in openmc/lib/tally.py // // These are kept as a normal enum and made negative, since variables which // store one of these enum values usually also may be responsible for storing diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 75166b36a..5b632bbce 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -535,6 +535,11 @@ private: vector secondary_bank_; + // Keep track of how many secondary particles were created in the collision + // and what the starting index is in the secondary bank for this particle + int n_secondaries_ {0}; + int secondary_bank_index_ {0}; + int64_t current_work_; vector flux_derivs_; @@ -689,7 +694,20 @@ public: // secondary particle bank SourceSite& secondary_bank(int i) { return secondary_bank_[i]; } + const SourceSite& secondary_bank(int i) const { return secondary_bank_[i]; } decltype(secondary_bank_)& secondary_bank() { return secondary_bank_; } + decltype(secondary_bank_) const& secondary_bank() const + { + return secondary_bank_; + } + + // Number of secondaries created in a collision + int& n_secondaries() { return n_secondaries_; } + const int& n_secondaries() const { return n_secondaries_; } + + // Starting index in secondary bank for this collision + int& secondary_bank_index() { return secondary_bank_index_; } + const int& secondary_bank_index() const { return secondary_bank_index_; } // Current simulation work index int64_t& current_work() { return current_work_; } diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h index 65098597a..0fef97e50 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_PRODUCTION, POLAR, SPHERICAL_HARMONICS, SPATIAL_LEGENDRE, diff --git a/include/openmc/tallies/filter_energy.h b/include/openmc/tallies/filter_energy.h index cf8a8aa0f..6bc9ff6b1 100644 --- a/include/openmc/tallies/filter_energy.h +++ b/include/openmc/tallies/filter_energy.h @@ -1,6 +1,7 @@ #ifndef OPENMC_TALLIES_FILTER_ENERGY_H #define OPENMC_TALLIES_FILTER_ENERGY_H +#include "openmc/particle.h" #include "openmc/span.h" #include "openmc/tallies/filter.h" #include "openmc/vector.h" @@ -72,5 +73,36 @@ public: std::string text_label(int bin) const override; }; +//============================================================================== +//! Bins the outgoing energy of secondary particles +//! +//! This filter can be used to get the photon production matrix for multigroup +//! photon transport, the energy distribution of secondary neutrons, etc. Unlike +//! other energy filters, the weight that is applied is equal to the weight of +//! the secondary particle. Thus, to get secondary production it should be used +//! in conjunction with the "events" score. +//============================================================================== + +class ParticleProductionFilter : public EnergyFilter { +public: + //---------------------------------------------------------------------------- + // Methods + + std::string type_str() const override { return "particleproduction"; } + FilterType type() const override { return FilterType::PARTICLE_PRODUCTION; } + + void get_all_bins(const Particle& p, TallyEstimator estimator, + FilterMatch& match) const override; + + std::string text_label(int bin) const override; + + void from_xml(pugi::xml_node node) override; + + void to_statepoint(hid_t filter_group) const override; + +protected: + ParticleType secondary_type_; //!< Type of secondary particle to filter +}; + } // namespace openmc #endif // OPENMC_TALLIES_FILTER_ENERGY_H diff --git a/openmc/filter.py b/openmc/filter.py index 53ec93e21..7561f6be6 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -126,9 +126,9 @@ class Filter(IDManagerMixin, metaclass=FilterMeta): def __gt__(self, other): if type(self) is not type(other): if self.short_name in _FILTER_TYPES and \ - other.short_name in _FILTER_TYPES: + other.short_name in _FILTER_TYPES: delta = _FILTER_TYPES.index(self.short_name) - \ - _FILTER_TYPES.index(other.short_name) + _FILTER_TYPES.index(other.short_name) return delta > 0 else: return False @@ -275,7 +275,6 @@ class Filter(IDManagerMixin, metaclass=FilterMeta): if filter_type == subclass.short_name.lower(): return subclass.from_xml_element(elem, **kwargs) - def can_merge(self, other): """Determine if filter can be merged with another. @@ -426,6 +425,7 @@ class Filter(IDManagerMixin, metaclass=FilterMeta): class WithIDFilter(Filter): """Abstract parent for filters of types with IDs (Cell, Material, etc.).""" + def __init__(self, bins, filter_id=None): bins = np.atleast_1d(bins) @@ -629,6 +629,7 @@ class CellInstanceFilter(Filter): DistribcellFilter """ + def __init__(self, bins, filter_id=None): self.bins = bins self.id = filter_id @@ -749,6 +750,7 @@ class ParticleFilter(Filter): The number of filter bins """ + def __eq__(self, other): if type(self) is not type(other): return False @@ -1103,6 +1105,7 @@ class MeshMaterialFilter(MeshFilter): Unique identifier for the filter """ + def __init__(self, mesh: openmc.MeshBase, bins, filter_id=None): self.mesh = mesh self.bins = bins @@ -1437,6 +1440,7 @@ class RealFilter(Filter): The number of filter bins """ + def __init__(self, values, filter_id=None): self.values = np.asarray(values) self.bins = np.vstack((self.values[:-1], self.values[1:])).T @@ -1712,7 +1716,8 @@ class EnergyFilter(RealFilter): """ - cv.check_value('group_structure', group_structure, openmc.mgxs.GROUP_STRUCTURES.keys()) + cv.check_value('group_structure', group_structure, + openmc.mgxs.GROUP_STRUCTURES.keys()) return cls(openmc.mgxs.GROUP_STRUCTURES[group_structure.upper()]) @@ -1742,6 +1747,82 @@ class EnergyoutFilter(EnergyFilter): """ + +class ParticleProductionFilter(EnergyFilter): + """Bins tally events based on energy of secondary particles. + + This filter bins the energies of secondary particles (e.g., photons, + electrons, or recoils) produced in a reaction. This is useful for + constructing production matrices or analyzing secondary particle spectra. + Note that unlike other energy filters, the weight that is applied is equal + to the weight of the secondary particle. Thus, to obtain secondary particle + production, it should be used in conjunction with the "events" score. + + The incident particle type can be filtered using :class:`ParticleFilter`. + + .. versionadded:: 0.15.4 + + Parameters + ---------- + particle : str, int, openmc.ParticleType + Type of secondary particle to tally ('photon', 'neutron', etc.) + values : Iterable of Real + A list of energy boundaries in [eV]; each successive pair defines a bin. + filter_id : int, optional + Unique identifier for the filter + + Attributes + ---------- + values : numpy.ndarray + Energy boundaries in [eV] + bins : numpy.ndarray + Array of (low, high) energy bin pairs + num_bins : int + Number of filter bins + particle : str + The secondary particle type this filter applies to + """ + + def __init__(self, particle, values, filter_id=None): + super().__init__(values, filter_id) + self.particle = particle + + def __repr__(self): + string = type(self).__name__ + '\n' + string += '{: <16}=\t{}\n'.format('\tParticle', self.particle) + string += '{: <16}=\t{}\n'.format('\tValues', self.values) + string += '{: <16}=\t{}\n'.format('\tID', self.id) + return string + + @property + def particle(self) -> openmc.ParticleType: + return self._particle + + @particle.setter + def particle(self, particle): + self._particle = openmc.ParticleType(particle) + + def to_xml_element(self): + element = super().to_xml_element() + subelement = ET.SubElement(element, 'particle') + subelement.text = str(self.particle) + return element + + @classmethod + def from_xml_element(cls, elem, **kwargs): + filter_id = int(elem.get('id')) + values = [float(x) for x in get_text(elem, 'bins').split()] + particle = get_text(elem, 'particle') + return cls(particle, values, filter_id=filter_id) + + @classmethod + def from_hdf5(cls, group, **kwargs): + filter_id = int(group.name.split('/')[-1].lstrip('filter ')) + bins = group['bins'][()] + particle = group['particle'][()].decode() + return cls(particle, bins, filter_id=filter_id) + + class TimeFilter(RealFilter): """Bins tally events based on the particle's time. @@ -1909,7 +1990,7 @@ class DistribcellFilter(Filter): # Make sure there is only 1 bin. if not len(bins) == 1: msg = (f'Unable to add bins "{bins}" to a DistribcellFilter since ' - 'only a single distribcell can be used per tally') + 'only a single distribcell can be used per tally') raise ValueError(msg) # Check the type and extract the id, if necessary. @@ -2062,7 +2143,7 @@ class DistribcellFilter(Filter): # requests Summary geometric information filter_bins = _repeat_and_tile( np.arange(self.num_bins), stride, data_size) - df = pd.DataFrame({self.short_name.lower() : filter_bins}) + df = pd.DataFrame({self.short_name.lower(): filter_bins}) # Concatenate with DataFrame of distribcell instance IDs if level_df is not None: @@ -2101,6 +2182,7 @@ class MuFilter(RealFilter): The number of filter bins """ + def __init__(self, values, filter_id=None): if isinstance(values, Integral): values = np.linspace(-1., 1., values + 1) @@ -2266,6 +2348,7 @@ class DelayedGroupFilter(Filter): The number of filter bins """ + def check_bins(self, bins): # Check the bin values. for g in bins: @@ -2334,9 +2417,9 @@ class EnergyFunctionFilter(Filter): def __gt__(self, other): if type(self) is not type(other): if self.short_name in _FILTER_TYPES and \ - other.short_name in _FILTER_TYPES: + other.short_name in _FILTER_TYPES: delta = _FILTER_TYPES.index(self.short_name) - \ - _FILTER_TYPES.index(other.short_name) + _FILTER_TYPES.index(other.short_name) return delta > 0 else: return False @@ -2346,9 +2429,9 @@ class EnergyFunctionFilter(Filter): def __lt__(self, other): if type(self) is not type(other): if self.short_name in _FILTER_TYPES and \ - other.short_name in _FILTER_TYPES: + other.short_name in _FILTER_TYPES: delta = _FILTER_TYPES.index(self.short_name) - \ - _FILTER_TYPES.index(other.short_name) + _FILTER_TYPES.index(other.short_name) return delta < 0 else: return False @@ -2359,14 +2442,16 @@ class EnergyFunctionFilter(Filter): string = type(self).__name__ + '\n' string += '{: <16}=\t{}\n'.format('\tEnergy', self.energy) string += '{: <16}=\t{}\n'.format('\tInterpolant', self.y) - string += '{: <16}=\t{}\n'.format('\tInterpolation', self.interpolation) + string += '{: <16}=\t{}\n'.format('\tInterpolation', + self.interpolation) return hash(string) def __repr__(self): string = type(self).__name__ + '\n' string += '{: <16}=\t{}\n'.format('\tEnergy', self.energy) string += '{: <16}=\t{}\n'.format('\tInterpolant', self.y) - string += '{: <16}=\t{}\n'.format('\tInterpolation', self.interpolation) + string += '{: <16}=\t{}\n'.format('\tInterpolation', + self.interpolation) string += '{: <16}=\t{}\n'.format('\tID', self.id) return string @@ -2452,10 +2537,12 @@ class EnergyFunctionFilter(Filter): @interpolation.setter def interpolation(self, val): cv.check_type('interpolation', val, str) - cv.check_value('interpolation', val, self.INTERPOLATION_SCHEMES.values()) + cv.check_value('interpolation', val, + self.INTERPOLATION_SCHEMES.values()) if val == 'quadratic' and len(self.energy) < 3: - raise ValueError('Quadratic interpolation requires 3 or more values.') + raise ValueError( + 'Quadratic interpolation requires 3 or more values.') if val == 'cubic' and len(self.energy) < 4: raise ValueError('Cubic interpolation requires 3 or more values.') diff --git a/src/particle.cpp b/src/particle.cpp index a1176abc7..a035e76b9 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -87,6 +87,9 @@ bool Particle::create_secondary( return false; } + // Increment number of secondaries created (for ParticleProductionFilter) + n_secondaries()++; + auto& bank = secondary_bank().emplace_back(); bank.particle = type; bank.wgt = wgt; @@ -336,6 +339,7 @@ void Particle::event_cross_surface() void Particle::event_collide() { + // Score collision estimate of keff if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) { keff_tally_collision() += wgt() * macro_xs().nu_fission / macro_xs().total; @@ -383,6 +387,11 @@ void Particle::event_collide() n_bank() = 0; bank_second_E() = 0.0; wgt_bank() = 0.0; + + // Clear number of secondaries in this collision. This is + // distinct from the number of created neutrons n_bank() above! + n_secondaries() = 0; + zero_delayed_bank(); // Reset fission logical diff --git a/src/physics.cpp b/src/physics.cpp index 05a9e59d9..1f72e4fac 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -44,6 +44,7 @@ void collision(Particle& p) { // Add to collision counter for particle ++(p.n_collision()); + p.secondary_bank_index() = p.secondary_bank().size(); // Sample reaction for the material the particle is in switch (p.type().pdg_number()) { @@ -253,6 +254,7 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) } } else { p.secondary_bank().push_back(site); + p.n_secondaries()++; } // Increment the number of neutrons born delayed diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 866d4d728..58a21a678 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -27,6 +27,7 @@ void collision_mg(Particle& p) { // Add to the collision counter for the particle p.n_collision()++; + p.secondary_bank_index() = p.secondary_bank().size(); // Sample the reaction type sample_reaction(p); @@ -200,6 +201,7 @@ void create_fission_sites(Particle& p) } } else { p.secondary_bank().push_back(site); + p.n_secondaries()++; } // Set the delayed group on the particle as well diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index 79817981d..16e4c0550 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 == "particleproduction") { + return Filter::create(id); } else if (type == "polar") { return Filter::create(id); } else if (type == "surface") { diff --git a/src/tallies/filter_delayedgroup.cpp b/src/tallies/filter_delayedgroup.cpp index 01e39e554..46adf529b 100644 --- a/src/tallies/filter_delayedgroup.cpp +++ b/src/tallies/filter_delayedgroup.cpp @@ -27,7 +27,7 @@ void DelayedGroupFilter::set_groups(span groups) } else if (group > MAX_DELAYED_GROUPS) { throw std::invalid_argument { "Encountered delayedgroup bin with index " + std::to_string(group) + - " which is greater than MAX_DELATED_GROUPS (" + + " which is greater than MAX_DELAYED_GROUPS (" + std::to_string(MAX_DELAYED_GROUPS) + ")"}; } groups_.push_back(group); @@ -39,6 +39,10 @@ void DelayedGroupFilter::set_groups(span groups) void DelayedGroupFilter::get_all_bins( const Particle& p, TallyEstimator estimator, FilterMatch& match) const { + // Note that the bin is set to zero here, but bins outside zero are + // tallied to regardless. This is because that logic has to be handled + // in the scoring code instead where looping over the delayed + // group takes place (tally_scoring.cpp). match.bins_.push_back(0); match.weights_.push_back(1.0); } diff --git a/src/tallies/filter_energy.cpp b/src/tallies/filter_energy.cpp index 0b954cce3..0b039440a 100644 --- a/src/tallies/filter_energy.cpp +++ b/src/tallies/filter_energy.cpp @@ -117,6 +117,50 @@ std::string EnergyoutFilter::text_label(int bin) const "Outgoing Energy [{}, {})", bins_.at(bin), bins_.at(bin + 1)); } +//============================================================================== +// ParticleProductionFilter implementation +//============================================================================== + +void ParticleProductionFilter::get_all_bins( + const Particle& p, TallyEstimator estimator, FilterMatch& match) const +{ + int start_idx = p.secondary_bank_index(); + int end_idx = start_idx + p.n_secondaries(); + + // Loop over secondary bank entries + for (int bank_idx = start_idx; bank_idx < end_idx; bank_idx++) { + // Check if this is the correct type of secondary, then match its energy if + // it's the right type + const auto& site = p.secondary_bank(bank_idx); + if (site.particle == secondary_type_) { + if (site.E >= bins_.front() && site.E <= bins_.back()) { + auto bin = lower_bound_index(bins_.begin(), bins_.end(), site.E); + match.bins_.push_back(bin); + match.weights_.push_back(site.wgt); + } + } + } +} + +std::string ParticleProductionFilter::text_label(int bin) const +{ + return fmt::format("Secondary {}, Energy [{}, {})", secondary_type_.str(), + bins_.at(bin), bins_.at(bin + 1)); +} + +void ParticleProductionFilter::from_xml(pugi::xml_node node) +{ + EnergyFilter::from_xml(node); + std::string p = get_node_value(node, "particle"); + secondary_type_ = ParticleType {p}; +} + +void ParticleProductionFilter::to_statepoint(hid_t filter_group) const +{ + EnergyFilter::to_statepoint(filter_group); + write_dataset(filter_group, "particle", secondary_type_.str()); +} + //============================================================================== // C-API functions //============================================================================== diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 2432f5c2f..8cdf3a905 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -171,6 +171,8 @@ Tally::Tally(pugi::xml_node node) filt_type == FilterType::ZERNIKE || filt_type == FilterType::ZERNIKE_RADIAL) { estimator_ = TallyEstimator::COLLISION; + } else if (filt_type == FilterType::PARTICLE_PRODUCTION) { + estimator_ = TallyEstimator::ANALOG; } } diff --git a/tests/regression_tests/model_xml/photon_production_inputs_true.dat b/tests/regression_tests/model_xml/photon_production_inputs_true.dat index 07eebaa3e..cc618eaff 100644 --- a/tests/regression_tests/model_xml/photon_production_inputs_true.dat +++ b/tests/regression_tests/model_xml/photon_production_inputs_true.dat @@ -41,6 +41,16 @@ neutron photon electron positron + + neutron + + + 0.0 20000000.0 + + + 0.0 100000.0 300000.0 500000.0 2000000.0 20000000.0 + photon + 1 2 current @@ -63,5 +73,10 @@ total heating (n,gamma) analog + + 5 3 4 + events + analog + diff --git a/tests/regression_tests/photon_production/inputs_true.dat b/tests/regression_tests/photon_production/inputs_true.dat index 07eebaa3e..cc618eaff 100644 --- a/tests/regression_tests/photon_production/inputs_true.dat +++ b/tests/regression_tests/photon_production/inputs_true.dat @@ -41,6 +41,16 @@ neutron photon electron positron + + neutron + + + 0.0 20000000.0 + + + 0.0 100000.0 300000.0 500000.0 2000000.0 20000000.0 + photon + 1 2 current @@ -63,5 +73,10 @@ total heating (n,gamma) analog + + 5 3 4 + events + analog + diff --git a/tests/regression_tests/photon_production/results_true.dat b/tests/regression_tests/photon_production/results_true.dat index 413f6f0ca..ec319aa0e 100644 --- a/tests/regression_tests/photon_production/results_true.dat +++ b/tests/regression_tests/photon_production/results_true.dat @@ -138,3 +138,14 @@ tally 4: 2.173936E+08 0.000000E+00 0.000000E+00 +tally 5: +2.000000E-02 +4.000000E-04 +8.200000E-03 +6.724000E-05 +5.280000E-02 +2.787840E-03 +3.546000E-01 +1.257412E-01 +5.063000E-01 +2.563397E-01 diff --git a/tests/regression_tests/photon_production/test.py b/tests/regression_tests/photon_production/test.py index 150448a12..65d0688e9 100644 --- a/tests/regression_tests/photon_production/test.py +++ b/tests/regression_tests/photon_production/test.py @@ -25,7 +25,8 @@ def model(): inner_cyl_right.region = -cyl & +x_plane_center & -x_plane_right outer_cyl.region = ~(-cyl & +x_plane_left & -x_plane_right) inner_cyl_right.fill = mat - model.geometry = openmc.Geometry([inner_cyl_left, inner_cyl_right, outer_cyl]) + model.geometry = openmc.Geometry( + [inner_cyl_left, inner_cyl_right, outer_cyl]) source = openmc.IndependentSource() source.space = openmc.stats.Point((0, 0, 0)) @@ -38,17 +39,19 @@ def model(): model.settings.batches = 1 model.settings.photon_transport = True model.settings.electron_treatment = 'ttb' - model.settings.cutoff = {'energy_photon' : 1000.0} + model.settings.cutoff = {'energy_photon': 1000.0} model.settings.source = source surface_filter = openmc.SurfaceFilter(cyl) - particle_filter = openmc.ParticleFilter(['neutron', 'photon', 'electron', 'positron']) + particle_filter = openmc.ParticleFilter( + ['neutron', 'photon', 'electron', 'positron']) current_tally = openmc.Tally() current_tally.filters = [surface_filter, particle_filter] current_tally.scores = ['current'] tally_tracklength = openmc.Tally() tally_tracklength.filters = [particle_filter] - tally_tracklength.scores = ['total', '(n,gamma)'] # heating doesn't work with tracklength + # heating doesn't work with tracklength + tally_tracklength.scores = ['total', '(n,gamma)'] tally_tracklength.nuclides = ['Al27', 'total'] tally_tracklength.estimator = 'tracklength' tally_collision = openmc.Tally() @@ -61,8 +64,25 @@ def model(): tally_analog.scores = ['total', 'heating', '(n,gamma)'] tally_analog.nuclides = ['Al27', 'total'] tally_analog.estimator = 'analog' + + # This is an analog tally tracking the energy distribution of photons + # generated by neutrons. The sum of the tally should give the total + # number of photons generated per source neutron. + ene_filter = openmc.EnergyFilter([0.0, 20e6]) # incident neutron energy + + # Track source energies of secondary gammas + ene2_filter = openmc.ParticleProductionFilter( + 'photon', [0.0, 100e3, 300e3, 500e3, 2e6, 20e6]) + + neutron_only = openmc.ParticleFilter(['neutron']) + tally_gam_ene = openmc.Tally() + tally_gam_ene.filters = [neutron_only, ene_filter, ene2_filter] + tally_gam_ene.scores = ['events'] + tally_gam_ene.estimator = 'analog' + model.tallies.extend([current_tally, tally_tracklength, - tally_collision, tally_analog]) + tally_collision, tally_analog, + tally_gam_ene]) return model diff --git a/tests/unit_tests/test_filters.py b/tests/unit_tests/test_filters.py index 8c56a310e..ba3226e3a 100644 --- a/tests/unit_tests/test_filters.py +++ b/tests/unit_tests/test_filters.py @@ -17,14 +17,16 @@ def box_model(): model.settings.particles = 100 model.settings.batches = 10 model.settings.inactive = 0 - model.settings.source = openmc.IndependentSource(space=openmc.stats.Point()) + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Point()) return model def test_cell_instance(): c1 = openmc.Cell() c2 = openmc.Cell() - f = openmc.CellInstanceFilter([(c1, 0), (c1, 1), (c1, 2), (c2, 0), (c2, 1)]) + f = openmc.CellInstanceFilter( + [(c1, 0), (c1, 1), (c1, 2), (c2, 0), (c2, 1)]) # Make sure __repr__ works repr(f) @@ -234,7 +236,8 @@ def test_first_moment(run_in_tmpdir, box_model): flux, scatter = sp.tallies[plain_tally.id].mean.ravel() # Check that first moment matches - first_score = lambda t: sp.tallies[t.id].mean.ravel()[0] + def first_score(t): + return sp.tallies[t.id].mean.ravel()[0] assert first_score(leg_tally) == scatter assert first_score(leg_sptl_tally) == scatter assert first_score(sph_scat_tally) == scatter @@ -258,7 +261,8 @@ def test_lethargy_bin_width(): assert len(f.lethargy_bin_width) == 175 energy_bins = openmc.mgxs.GROUP_STRUCTURES['VITAMIN-J-175'] assert f.lethargy_bin_width[0] == np.log10(energy_bins[1]/energy_bins[0]) - assert f.lethargy_bin_width[-1] == np.log10(energy_bins[-1]/energy_bins[-2]) + assert f.lethargy_bin_width[-1] == np.log10( + energy_bins[-1]/energy_bins[-2]) def test_energyfunc(): @@ -292,7 +296,8 @@ def test_tabular_from_energyfilter(): # 'histogram' is the default assert tab.interpolation == 'histogram' - tab = efilter.get_tabular(values=np.array([10, 10, 5]), interpolation='linear-linear') + tab = efilter.get_tabular(values=np.array( + [10, 10, 5]), interpolation='linear-linear') assert tab.interpolation == 'linear-linear' @@ -314,6 +319,38 @@ def test_energy_filter(): openmc.EnergyFilter([-1.2, 0.25, 0.5]) +def test_particle_production_filter(): + energy_bins = [1e3, 1e4, 1e5, 1e6] + f = openmc.ParticleProductionFilter('photon', energy_bins) + + assert f.particle == openmc.ParticleType.PHOTON + assert f.num_bins == 3 + assert f.bins.shape == (3, 2) + assert np.allclose(f.values, energy_bins) + + # __repr__ check + repr(f) + + # to_xml_element() + elem = f.to_xml_element() + assert elem.tag == 'filter' + assert elem.attrib['type'] == 'particleproduction' + assert elem.find('particle').text == 'photon' + assert elem.find('bins').text.split()[0] == str(energy_bins[0]) + + # from_xml_element() + new_f = openmc.Filter.from_xml_element(elem) + assert new_f.id == f.id + assert new_f.particle == f.particle + assert np.allclose(new_f.bins, f.bins) + + # pandas output + df = f.get_pandas_dataframe(data_size=3, stride=1) + assert df.shape[0] == 3 + assert "particleproduction low [eV]" in df.columns + assert "particleproduction high [eV]" in df.columns + + def test_weight(): f = openmc.WeightFilter([0.01, 0.1, 1.0, 10.0]) expected_bins = [[0.01, 0.1], [0.1, 1.0], [1.0, 10.0]]