diff --git a/CMakeLists.txt b/CMakeLists.txt index d79746209..4dff35418 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -425,6 +425,7 @@ list(APPEND libopenmc_SOURCES src/tallies/filter_surface.cpp src/tallies/filter_time.cpp src/tallies/filter_universe.cpp + src/tallies/filter_weight.cpp src/tallies/filter_zernike.cpp src/tallies/tally.cpp src/tallies/tally_scoring.cpp diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index f3eb7ce83..1d8a93d47 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -143,6 +143,7 @@ Constructing Tallies openmc.SpatialLegendreFilter openmc.SphericalHarmonicsFilter openmc.TimeFilter + openmc.WeightFilter openmc.ZernikeFilter openmc.ZernikeRadialFilter openmc.ParentNuclideFilter diff --git a/docs/source/pythonapi/capi.rst b/docs/source/pythonapi/capi.rst index c8e0e874d..63fc5d11b 100644 --- a/docs/source/pythonapi/capi.rst +++ b/docs/source/pythonapi/capi.rst @@ -91,6 +91,7 @@ Classes Tally UniverseFilter UnstructuredMesh + WeightFilter WeightWindows ZernikeFilter ZernikeRadialFilter diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h index 3e982d0cf..ee635b183 100644 --- a/include/openmc/tallies/filter.h +++ b/include/openmc/tallies/filter.h @@ -44,6 +44,7 @@ enum class FilterType { SURFACE, TIME, UNIVERSE, + WEIGHT, ZERNIKE, ZERNIKE_RADIAL }; diff --git a/include/openmc/tallies/filter_weight.h b/include/openmc/tallies/filter_weight.h new file mode 100644 index 000000000..1fe9d75d3 --- /dev/null +++ b/include/openmc/tallies/filter_weight.h @@ -0,0 +1,51 @@ +#ifndef OPENMC_TALLIES_FILTER_WEIGHT_H +#define OPENMC_TALLIES_FILTER_WEIGHT_H + +#include + +#include "openmc/span.h" +#include "openmc/tallies/filter.h" +#include "openmc/vector.h" + +namespace openmc { + +//============================================================================== +//! Bins the weights of the particles. +//============================================================================== + +class WeightFilter : public Filter { +public: + //---------------------------------------------------------------------------- + // Constructors, destructors + + ~WeightFilter() = default; + + //---------------------------------------------------------------------------- + // Methods + + std::string type_str() const override { return "weight"; } + FilterType type() const override { return FilterType::WEIGHT; } + + 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& bins() const { return bins_; } + void set_bins(span bins); + +protected: + //---------------------------------------------------------------------------- + // Data members + vector bins_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_WEIGHT_H diff --git a/openmc/filter.py b/openmc/filter.py index 755777c5e..4d37f1055 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -25,7 +25,7 @@ _FILTER_TYPES = ( 'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre', 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance', - 'collision', 'time', 'parentnuclide' + 'collision', 'time', 'parentnuclide', 'weight' ) _CURRENT_NAMES = ( @@ -2328,3 +2328,26 @@ class EnergyFunctionFilter(Filter): {self.short_name.lower(): filter_bins})]) return df + + +class WeightFilter(RealFilter): + """Bins tally events based on the incoming particle weight. + + Parameters + ---------- + Values : Iterable of float + A list or iterable of the weight boundaries, as float values. + filter_id : int + Unique identifier for the filter + + Attributes + ---------- + id : int + Unique identifier for the filter + bins : numpy.ndarray + An array of integer values representing the weights by which to filter + num_bins : int + The number of filter bins + values : numpy.ndarray + Array of weight boundaries + """ diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py index b6086c567..3a76b52a1 100644 --- a/openmc/lib/filter.py +++ b/openmc/lib/filter.py @@ -610,6 +610,10 @@ class UniverseFilter(Filter): filter_type = 'universe' +class WeightFilter(Filter): + filter_type = 'weight' + + class ZernikeFilter(Filter): filter_type = 'zernike' diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp index 27b750b93..e90aa7490 100644 --- a/src/hdf5_interface.cpp +++ b/src/hdf5_interface.cpp @@ -92,7 +92,7 @@ void get_shape_attr(hid_t obj_id, const char* name, hsize_t* dims) H5Aclose(attr); } -hid_t create_group(hid_t parent_id, char const* name) +hid_t create_group(hid_t parent_id, const char* name) { hid_t out = H5Gcreate(parent_id, name, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); if (out < 0) { diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index 6430d9ae1..17e57a987 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -36,6 +36,7 @@ #include "openmc/tallies/filter_surface.h" #include "openmc/tallies/filter_time.h" #include "openmc/tallies/filter_universe.h" +#include "openmc/tallies/filter_weight.h" #include "openmc/tallies/filter_zernike.h" #include "openmc/xml_interface.h" @@ -154,6 +155,8 @@ Filter* Filter::create(const std::string& type, int32_t id) return Filter::create(id); } else if (type == "universe") { return Filter::create(id); + } else if (type == "weight") { + return Filter::create(id); } else if (type == "zernike") { return Filter::create(id); } else if (type == "zernikeradial") { diff --git a/src/tallies/filter_weight.cpp b/src/tallies/filter_weight.cpp new file mode 100644 index 000000000..31f4bd1bf --- /dev/null +++ b/src/tallies/filter_weight.cpp @@ -0,0 +1,63 @@ +#include "openmc/tallies/filter_weight.h" + +#include // for is_sorted +#include // for runtime_error + +#include + +#include "openmc/search.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +//============================================================================== +// WeightFilter implementation +//============================================================================== + +void WeightFilter::from_xml(pugi::xml_node node) +{ + auto bins = get_node_array(node, "bins"); + this->set_bins(bins); +} + +void WeightFilter::set_bins(span bins) +{ + if (!std::is_sorted(bins.begin(), bins.end())) { + throw std::runtime_error {"Weight bins must be monotonically increasing."}; + } + + // Clear existing bins + bins_.clear(); + bins_.reserve(bins.size()); + + // Copy bins + bins_.insert(bins_.end(), bins.begin(), bins.end()); + n_bins_ = bins_.size() - 1; +} + +void WeightFilter::get_all_bins( + const Particle& p, TallyEstimator estimator, FilterMatch& match) const +{ + // Get particle weight + double wgt = p.wgt_last(); + + // Bin the weight + if (wgt >= bins_.front() && wgt <= bins_.back()) { + auto bin = lower_bound_index(bins_.begin(), bins_.end(), wgt); + match.bins_.push_back(bin); + match.weights_.push_back(1.0); + } +} + +void WeightFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + write_dataset(filter_group, "bins", bins_); +} + +std::string WeightFilter::text_label(int bin) const +{ + return fmt::format("Weight [{}, {}]", bins_[bin], bins_[bin + 1]); +} + +} // namespace openmc diff --git a/tests/unit_tests/test_filter_weight.py b/tests/unit_tests/test_filter_weight.py new file mode 100644 index 000000000..878929ee0 --- /dev/null +++ b/tests/unit_tests/test_filter_weight.py @@ -0,0 +1,44 @@ +import openmc +import numpy as np + + +def test_weightfilter(run_in_tmpdir): + steel = openmc.Material(name='Stainless Steel') + steel.set_density('g/cm3', 8.00) + steel.add_nuclide('Fe56', 1.0) + + sphere = openmc.Sphere(r=50.0, boundary_type='vacuum') + cell = openmc.Cell(region=-sphere, fill=steel) + model = openmc.Model() + model.geometry = openmc.Geometry([cell]) + model.settings.particles = 100 + model.settings.batches = 10 + + model.settings.source = openmc.IndependentSource( + energy=openmc.stats.delta_function(14e6), + ) + model.settings.run_mode = "fixed source" + + radius = list(range(1, 50)) + sphere_mesh = openmc.SphericalMesh(radius) + mesh_filter = openmc.MeshFilter(sphere_mesh) + weight_filter = openmc.WeightFilter( + [0.999, 0.9999, 0.99999, 0.999999, 1.0, 1.000001 ,1.00001, 1.0001, 1.001] + ) + + tally = openmc.Tally() + tally.filters = [mesh_filter, weight_filter] + tally.estimator = 'analog' + tally.scores = ['flux'] + model.tallies = openmc.Tallies([tally]) + + # Run OpenMC + model.run(apply_tally_results=True) + + # Get current binned by mu + neutron_flux = tally.mean.reshape(48, 8) + + # All contributions should show up in the fourth bin + assert np.all(neutron_flux[:, 3] != 0.0) + neutron_flux[:, 3] = 0.0 + assert np.all(neutron_flux == 0.0) diff --git a/tests/unit_tests/test_filters.py b/tests/unit_tests/test_filters.py index c3b737192..60a60e3f5 100644 --- a/tests/unit_tests/test_filters.py +++ b/tests/unit_tests/test_filters.py @@ -312,3 +312,24 @@ def test_energy_filter(): msg = 'Unable to set "filter value" to "-1.2" since it is less than "0.0"' with raises(ValueError, match=msg): openmc.EnergyFilter([-1.2, 0.25, 0.5]) + + +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]] + + assert np.allclose(f.bins, expected_bins) + assert len(f.bins) == 3 + + # Make sure __repr__ works + repr(f) + + # to_xml_element() + elem = f.to_xml_element() + assert elem.tag == 'filter' + assert elem.attrib['type'] == 'weight' + + # from_xml_element() + new_f = openmc.Filter.from_xml_element(elem) + assert new_f.id == f.id + assert np.allclose(new_f.bins, f.bins)