diff --git a/CMakeLists.txt b/CMakeLists.txt index 9fe133a22e3..75c0224b106 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -455,6 +455,8 @@ list(APPEND libopenmc_SOURCES src/tallies/filter_polar.cpp src/tallies/filter_reaction.cpp src/tallies/filter_sph_harm.cpp + src/tallies/filter_sptl_expansion.cpp + src/tallies/filter_sptl_fourier.cpp src/tallies/filter_sptl_legendre.cpp src/tallies/filter_surface.cpp src/tallies/filter_time.cpp diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index b3911c8b3ba..a9a2dc88b5c 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -142,6 +142,7 @@ Constructing Tallies openmc.EnergyFunctionFilter openmc.LegendreFilter openmc.SpatialLegendreFilter + openmc.SpatialFourierFilter openmc.SphericalHarmonicsFilter openmc.TimeFilter openmc.WeightFilter diff --git a/docs/source/pythonapi/capi.rst b/docs/source/pythonapi/capi.rst index dab45481dd9..a51796121c9 100644 --- a/docs/source/pythonapi/capi.rst +++ b/docs/source/pythonapi/capi.rst @@ -87,6 +87,7 @@ Classes ReactionFilter RectilinearMesh RegularMesh + SpatialFourierFilter SpatialLegendreFilter SphericalHarmonicsFilter SphericalMesh diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 911654d318f..f72519fc22d 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -192,6 +192,12 @@ int openmc_set_n_batches( int openmc_simulation_finalize(); int openmc_simulation_init(); int openmc_source_bank(void** ptr, int64_t* n); +int openmc_spatial_fourier_filter_get_order(int32_t index, int* order); +int openmc_spatial_fourier_filter_get_params( + int32_t index, int* axis, double* min, double* max); +int openmc_spatial_fourier_filter_set_order(int32_t index, int order); +int openmc_spatial_fourier_filter_set_params( + int32_t index, const int* axis, const double* min, const double* max); int openmc_spatial_legendre_filter_get_order(int32_t index, int* order); int openmc_spatial_legendre_filter_get_params( int32_t index, int* axis, double* min, double* max); diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h index 77b0d9f420d..ab8db12da90 100644 --- a/include/openmc/tallies/filter.h +++ b/include/openmc/tallies/filter.h @@ -43,6 +43,7 @@ enum class FilterType { POLAR, REACTION, SPHERICAL_HARMONICS, + SPATIAL_FOURIER, SPATIAL_LEGENDRE, SURFACE, TIME, diff --git a/include/openmc/tallies/filter_sptl_expansion.h b/include/openmc/tallies/filter_sptl_expansion.h new file mode 100644 index 00000000000..03e7f78ebc1 --- /dev/null +++ b/include/openmc/tallies/filter_sptl_expansion.h @@ -0,0 +1,70 @@ +#ifndef OPENMC_TALLIES_FILTER_SPTL_EXPANSION_H +#define OPENMC_TALLIES_FILTER_SPTL_EXPANSION_H + +#include "openmc/tallies/filter.h" + +namespace openmc { + +//============================================================================== +//! Abstract filter for functional expansions of the particle's normalized +//! position along a Cartesian axis +//============================================================================== + +class SpatialExpansionFilter : public Filter { +public: + enum class Axis { x, y, z }; + + //---------------------------------------------------------------------------- + // Constructors, destructors + + ~SpatialExpansionFilter() = default; + + //---------------------------------------------------------------------------- + // Methods + + void from_xml(pugi::xml_node node) override; + + void to_statepoint(hid_t filter_group) const override; + + //---------------------------------------------------------------------------- + // Accessors + + int order() const { return order_; } + + //! Set the expansion order and the corresponding number of bins + virtual void set_order(int order) = 0; + + Axis axis() const { return axis_; } + void set_axis(Axis axis); + + double min() const { return min_; } + double max() const { return max_; } + void set_minmax(double min, double max); + +protected: + //---------------------------------------------------------------------------- + // Methods + + //! Get the particle's coordinate along the expansion axis + double position(const Particle& p) const; + + //! Get the name of the expansion axis ("x", "y", or "z") + const char* axis_label() const; + + //---------------------------------------------------------------------------- + // Data members + + int order_; + + //! The Cartesian coordinate axis that the expansion is applied to. + Axis axis_; + + //! The minimum coordinate along the reference axis that the expansion covers. + double min_; + + //! The maximum coordinate along the reference axis that the expansion covers. + double max_; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_SPTL_EXPANSION_H diff --git a/include/openmc/tallies/filter_sptl_fourier.h b/include/openmc/tallies/filter_sptl_fourier.h new file mode 100644 index 00000000000..164799d00c3 --- /dev/null +++ b/include/openmc/tallies/filter_sptl_fourier.h @@ -0,0 +1,39 @@ +#ifndef OPENMC_TALLIES_FILTER_SPTL_FOURIER_H +#define OPENMC_TALLIES_FILTER_SPTL_FOURIER_H + +#include + +#include "openmc/tallies/filter_sptl_expansion.h" + +namespace openmc { + +//============================================================================== +//! Gives Fourier moments of the particle's normalized position along an axis +//============================================================================== + +class SpatialFourierFilter : public SpatialExpansionFilter { +public: + //---------------------------------------------------------------------------- + // Constructors, destructors + + ~SpatialFourierFilter() = default; + + //---------------------------------------------------------------------------- + // Methods + + std::string type_str() const override { return "spatialfourier"; } + FilterType type() const override { return FilterType::SPATIAL_FOURIER; } + + void get_all_bins(const Particle& p, TallyEstimator estimator, + FilterMatch& match) const override; + + std::string text_label(int bin) const override; + + //---------------------------------------------------------------------------- + // Accessors + + void set_order(int order) override; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_SPTL_FOURIER_H diff --git a/include/openmc/tallies/filter_sptl_legendre.h b/include/openmc/tallies/filter_sptl_legendre.h index b6c380e9b8d..5b7f2fd5aff 100644 --- a/include/openmc/tallies/filter_sptl_legendre.h +++ b/include/openmc/tallies/filter_sptl_legendre.h @@ -3,17 +3,15 @@ #include -#include "openmc/tallies/filter.h" +#include "openmc/tallies/filter_sptl_expansion.h" namespace openmc { -enum class LegendreAxis { x, y, z }; - //============================================================================== //! Gives Legendre moments of the particle's normalized position along an axis //============================================================================== -class SpatialLegendreFilter : public Filter { +class SpatialLegendreFilter : public SpatialExpansionFilter { public: //---------------------------------------------------------------------------- // Constructors, destructors @@ -26,42 +24,15 @@ class SpatialLegendreFilter : public Filter { std::string type_str() const override { return "spatiallegendre"; } FilterType type() const override { return FilterType::SPATIAL_LEGENDRE; } - 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 - int order() const { return order_; } - void set_order(int order); - - LegendreAxis axis() const { return axis_; } - void set_axis(LegendreAxis axis); - - double min() const { return min_; } - double max() const { return max_; } - void set_minmax(double min, double max); - -private: - //---------------------------------------------------------------------------- - // Data members - - int order_; - - //! The Cartesian coordinate axis that the Legendre expansion is applied to. - LegendreAxis axis_; - - //! The minimum coordinate along the reference axis that the expansion covers. - double min_; - - //! The maximum coordinate along the reference axis that the expansion covers. - double max_; + void set_order(int order) override; }; } // namespace openmc diff --git a/openmc/filter.py b/openmc/filter.py index 87aeb70c3a7..4626a10c789 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -25,7 +25,8 @@ 'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', 'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom', 'materialfrom', 'legendre', - 'spatiallegendre', 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', + 'spatialfourier', 'spatiallegendre', 'sphericalharmonics', + 'zernike', 'zernikeradial', 'particle', 'particleproduction', 'cellinstance', 'collision', 'time', 'parentnuclide', 'weight', 'meshborn', 'meshsurface', 'meshmaterial', 'reaction', ) diff --git a/openmc/filter_expansion.py b/openmc/filter_expansion.py index b79c8fc79e7..550533a1583 100644 --- a/openmc/filter_expansion.py +++ b/openmc/filter_expansion.py @@ -137,17 +137,18 @@ def from_hdf5(cls, group, **kwargs): return out -class SpatialLegendreFilter(ExpansionFilter): - r"""Score Legendre expansion moments in space up to specified order. +class SpatialExpansionFilter(ExpansionFilter): + """Abstract base class for spatial functional expansion filters. - This filter allows scores to be multiplied by Legendre polynomials of the - the particle's position along a particular axis, normalized to a given - range, up to a user-specified order. + This class provides common functionality for filters that expand + tally data along a spatial axis (x, y, or z) within a bounded region. + Subclasses must implement the order setter to define their specific + bin structure. Parameters ---------- order : int - Maximum Legendre polynomial order + Maximum expansion order axis : {'x', 'y', 'z'} Axis along which to take the expansion minimum : float @@ -160,7 +161,7 @@ class SpatialLegendreFilter(ExpansionFilter): Attributes ---------- order : int - Maximum Legendre polynomial order + Maximum expansion order axis : {'x', 'y', 'z'} Axis along which to take the expansion minimum : float @@ -197,11 +198,6 @@ def __repr__(self): string += '{: <16}=\t{}\n'.format('\tID', self.id) return string - @ExpansionFilter.order.setter - def order(self, order): - ExpansionFilter.order.__set__(self, order) - self.bins = [f'P{i}' for i in range(order + 1)] - @property def axis(self): return self._axis @@ -229,27 +225,13 @@ def maximum(self, maximum): cv.check_type('maximum', maximum, Real) self._maximum = maximum - @classmethod - def from_hdf5(cls, group, **kwargs): - if group['type'][()].decode() != cls.short_name.lower(): - raise ValueError("Expected HDF5 data for filter type '" - + cls.short_name.lower() + "' but got '" - + group['type'][()].decode() + " instead") - - filter_id = int(group.name.split('/')[-1].lstrip('filter ')) - order = group['order'][()] - axis = group['axis'][()].decode() - min_, max_ = group['min'][()], group['max'][()] - - return cls(order, axis, min_, max_, filter_id) - def to_xml_element(self): """Return XML Element representing the filter. Returns ------- element : lxml.etree._Element - XML element containing Legendre filter data + XML element containing spatial expansion filter data """ element = super().to_xml_element() @@ -259,7 +241,6 @@ def to_xml_element(self): subelement.text = str(self.minimum) subelement = ET.SubElement(element, 'max') subelement.text = str(self.maximum) - return element @classmethod @@ -271,6 +252,111 @@ def from_xml_element(cls, elem, **kwargs): maximum = float(get_text(elem, "max")) return cls(order, axis, minimum, maximum, filter_id=filter_id) + @classmethod + def from_hdf5(cls, group, **kwargs): + if group['type'][()].decode() != cls.short_name.lower(): + raise ValueError("Expected HDF5 data for filter type '" + + cls.short_name.lower() + "' but got '" + + group['type'][()].decode() + " instead") + + filter_id = int(group.name.split('/')[-1].lstrip('filter ')) + order = group['order'][()] + axis = group['axis'][()].decode() + min_, max_ = group['min'][()], group['max'][()] + + return cls(order, axis, min_, max_, filter_id) + + +class SpatialFourierFilter(SpatialExpansionFilter): + r"""Score Fourier expansion moments in space up to specified order. + + This filter allows scores to be multiplied by Fourier basis functions of + the particle's position along a particular axis, normalized to a given + range, up to a user-specified order. + + Parameters + ---------- + order : int + Maximum Fourier expansion order + axis : {'x', 'y', 'z'} + Axis along which to take the expansion + minimum : float + Minimum value along selected axis + maximum : float + Maximum value along selected axis + filter_id : int or None + Unique identifier for the filter + + Attributes + ---------- + order : int + Maximum Fourier expansion order + axis : {'x', 'y', 'z'} + Axis along which to take the expansion + minimum : float + Minimum value along selected axis + maximum : float + Maximum value along selected axis + id : int + Unique identifier for the filter + num_bins : int + The number of filter bins (2*order + 1) + + """ + + @ExpansionFilter.order.setter + def order(self, order): + ExpansionFilter.order.__set__(self, order) + self.bins = ['a0 (constant)'] + [None]*2*order + for i in range(1, order + 1): + a = 2*i - 1 + b = 2*i + self.bins[a] = f'a{i} (cos)' + self.bins[b] = f'b{i} (sin)' + + +class SpatialLegendreFilter(SpatialExpansionFilter): + r"""Score Legendre expansion moments in space up to specified order. + + This filter allows scores to be multiplied by Legendre polynomials of the + the particle's position along a particular axis, normalized to a given + range, up to a user-specified order. + + Parameters + ---------- + order : int + Maximum Legendre polynomial order + axis : {'x', 'y', 'z'} + Axis along which to take the expansion + minimum : float + Minimum value along selected axis + maximum : float + Maximum value along selected axis + filter_id : int or None + Unique identifier for the filter + + Attributes + ---------- + order : int + Maximum Legendre polynomial order + axis : {'x', 'y', 'z'} + Axis along which to take the expansion + minimum : float + Minimum value along selected axis + maximum : float + Maximum value along selected axis + id : int + Unique identifier for the filter + num_bins : int + The number of filter bins + + """ + + @ExpansionFilter.order.setter + def order(self, order): + ExpansionFilter.order.__set__(self, order) + self.bins = [f'P{i}' for i in range(order + 1)] + class SphericalHarmonicsFilter(ExpansionFilter): r"""Score spherical harmonic expansion moments up to specified order. diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py index 574a37443ae..50b8d80d8a0 100644 --- a/openmc/lib/filter.py +++ b/openmc/lib/filter.py @@ -23,7 +23,7 @@ 'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshBornFilter', 'MeshMaterialFilter', 'MeshSurfaceFilter', 'MuFilter', 'MuSurfaceFilter', 'ParentNuclideFilter', 'ParticleFilter', 'ParticleProductionFilter', 'PolarFilter', - 'ReactionFilter', 'SphericalHarmonicsFilter', 'SpatialLegendreFilter', + 'ReactionFilter', 'SphericalHarmonicsFilter', 'SpatialFourierFilter', 'SpatialLegendreFilter', 'SurfaceFilter', 'TimeFilter', 'UniverseFilter', 'WeightFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters' ] @@ -136,6 +136,12 @@ _dll.openmc_particle_filter_get_bins.argtypes = [c_int32, POINTER(c_int32)] _dll.openmc_particle_filter_get_bins.restype = c_int _dll.openmc_particle_filter_get_bins.errcheck = _error_handler +_dll.openmc_spatial_fourier_filter_get_order.argtypes = [c_int32, POINTER(c_int)] +_dll.openmc_spatial_fourier_filter_get_order.restype = c_int +_dll.openmc_spatial_fourier_filter_get_order.errcheck = _error_handler +_dll.openmc_spatial_fourier_filter_set_order.argtypes = [c_int32, c_int] +_dll.openmc_spatial_fourier_filter_set_order.restype = c_int +_dll.openmc_spatial_fourier_filter_set_order.errcheck = _error_handler _dll.openmc_spatial_legendre_filter_get_order.argtypes = [c_int32, POINTER(c_int)] _dll.openmc_spatial_legendre_filter_get_order.restype = c_int _dll.openmc_spatial_legendre_filter_get_order.errcheck = _error_handler @@ -639,6 +645,25 @@ def order(self, order): _dll.openmc_sphharm_filter_set_order(self._index, order) +class SpatialFourierFilter(Filter): + filter_type = 'spatialfourier' + + def __init__(self, order=None, uid=None, new=True, index=None): + super().__init__(uid, new, index) + if order is not None: + self.order = order + + @property + def order(self): + temp_order = c_int() + _dll.openmc_spatial_fourier_filter_get_order(self._index, temp_order) + return temp_order.value + + @order.setter + def order(self, order): + _dll.openmc_spatial_fourier_filter_set_order(self._index, order) + + class SpatialLegendreFilter(Filter): filter_type = 'spatiallegendre' @@ -724,6 +749,7 @@ class ZernikeRadialFilter(ZernikeFilter): 'polar': PolarFilter, 'reaction': ReactionFilter, 'sphericalharmonics': SphericalHarmonicsFilter, + 'spatialfourier': SpatialFourierFilter, 'spatiallegendre': SpatialLegendreFilter, 'surface': SurfaceFilter, 'time': TimeFilter, diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index badb9107733..956925754d7 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -35,6 +35,7 @@ #include "openmc/tallies/filter_polar.h" #include "openmc/tallies/filter_reaction.h" #include "openmc/tallies/filter_sph_harm.h" +#include "openmc/tallies/filter_sptl_fourier.h" #include "openmc/tallies/filter_sptl_legendre.h" #include "openmc/tallies/filter_surface.h" #include "openmc/tallies/filter_time.h" @@ -156,6 +157,8 @@ Filter* Filter::create(const std::string& type, int32_t id) return Filter::create(id); } else if (type == "surface") { return Filter::create(id); + } else if (type == "spatialfourier") { + return Filter::create(id); } else if (type == "spatiallegendre") { return Filter::create(id); } else if (type == "sphericalharmonics") { diff --git a/src/tallies/filter_sptl_expansion.cpp b/src/tallies/filter_sptl_expansion.cpp new file mode 100644 index 00000000000..44d83ee3ef2 --- /dev/null +++ b/src/tallies/filter_sptl_expansion.cpp @@ -0,0 +1,78 @@ +#include "openmc/tallies/filter_sptl_expansion.h" + +#include "openmc/xml_interface.h" + +namespace openmc { + +void SpatialExpansionFilter::from_xml(pugi::xml_node node) +{ + this->set_order(std::stoi(get_node_value(node, "order"))); + + auto axis = get_node_value(node, "axis"); + switch (axis[0]) { + case 'x': + this->set_axis(Axis::x); + break; + case 'y': + this->set_axis(Axis::y); + break; + case 'z': + this->set_axis(Axis::z); + break; + default: + throw std::runtime_error { + "Axis for spatial expansion filters must be 'x', 'y', or 'z'"}; + } + + double min = std::stod(get_node_value(node, "min")); + double max = std::stod(get_node_value(node, "max")); + this->set_minmax(min, max); +} + +void SpatialExpansionFilter::set_axis(Axis axis) +{ + axis_ = axis; +} + +void SpatialExpansionFilter::set_minmax(double min, double max) +{ + if (max <= min) { + throw std::invalid_argument { + "Maximum value must be greater than minimum value"}; + } + min_ = min; + max_ = max; +} + +double SpatialExpansionFilter::position(const Particle& p) const +{ + if (axis_ == Axis::x) { + return p.r().x; + } else if (axis_ == Axis::y) { + return p.r().y; + } else { + return p.r().z; + } +} + +const char* SpatialExpansionFilter::axis_label() const +{ + if (axis_ == Axis::x) { + return "x"; + } else if (axis_ == Axis::y) { + return "y"; + } else { + return "z"; + } +} + +void SpatialExpansionFilter::to_statepoint(hid_t filter_group) const +{ + Filter::to_statepoint(filter_group); + write_dataset(filter_group, "order", order_); + write_dataset(filter_group, "axis", axis_label()); + write_dataset(filter_group, "min", min_); + write_dataset(filter_group, "max", max_); +} + +} // namespace openmc diff --git a/src/tallies/filter_sptl_fourier.cpp b/src/tallies/filter_sptl_fourier.cpp new file mode 100644 index 00000000000..c1529e0b157 --- /dev/null +++ b/src/tallies/filter_sptl_fourier.cpp @@ -0,0 +1,151 @@ +#include "openmc/tallies/filter_sptl_fourier.h" + +#include // For pair + +#include + +#include "openmc/capi.h" +#include "openmc/constants.h" +#include "openmc/error.h" + +namespace openmc { + +void SpatialFourierFilter::set_order(int order) +{ + if (order < 0) { + throw std::invalid_argument {"Fourier order must be non-negative."}; + } + order_ = order; + n_bins_ = 2 * order_ + 1; +} + +void SpatialFourierFilter::get_all_bins( + const Particle& p, TallyEstimator estimator, FilterMatch& match) const +{ + // Get the coordinate along the axis of interest. + double x = this->position(p); + + if (x >= min_ && x <= max_) { + // Compute the normalized coordinate value on [0, 1] + double x_norm = (x - min_) / (max_ - min_); + + // Compute and return the Fourier weights. + vector wgt(n_bins_); + wgt[0] = 1.0; // a_0: constant term + for (int n = 1; n <= order_; ++n) { + double arg = 2.0 * PI * n * x_norm; + wgt[2 * n - 1] = std::cos(arg); + wgt[2 * n] = std::sin(arg); + } + for (int i = 0; i < n_bins_; ++i) { + match.bins_.push_back(i); + match.weights_.push_back(wgt[i]); + } + } +} + +std::string SpatialFourierFilter::text_label(int bin) const +{ + std::string func_str; + if (bin == 0) { + func_str = "a0 (constant)"; + } else if (bin % 2 == 1) { + int n = (bin + 1) / 2; + func_str = fmt::format("a{} (cos)", n); + } else { + int n = bin / 2; + func_str = fmt::format("b{} (sin)", n); + } + return fmt::format( + "Fourier expansion, {} axis, {}", this->axis_label(), func_str); +} + +//============================================================================== +// C-API functions +//============================================================================== + +std::pair check_sptl_fourier_filter(int32_t index) +{ + // Make sure this is a valid index to an allocated filter. + int err = verify_filter(index); + if (err) { + return {err, nullptr}; + } + + // Get a pointer to the filter and downcast. + const auto& filt_base = model::tally_filters[index].get(); + auto* filt = dynamic_cast(filt_base); + + // Check the filter type. + if (!filt) { + set_errmsg("Not a spatial Fourier filter."); + err = OPENMC_E_INVALID_TYPE; + } + return {err, filt}; +} + +extern "C" int openmc_spatial_fourier_filter_get_order( + int32_t index, int* order) +{ + // Check the filter. + auto check_result = check_sptl_fourier_filter(index); + auto err = check_result.first; + auto filt = check_result.second; + if (err) + return err; + + // Output the order. + *order = filt->order(); + return 0; +} + +extern "C" int openmc_spatial_fourier_filter_get_params( + int32_t index, int* axis, double* min, double* max) +{ + // Check the filter. + auto check_result = check_sptl_fourier_filter(index); + auto err = check_result.first; + auto filt = check_result.second; + if (err) + return err; + + // Output the params. + *axis = static_cast(filt->axis()); + *min = filt->min(); + *max = filt->max(); + return 0; +} + +extern "C" int openmc_spatial_fourier_filter_set_order(int32_t index, int order) +{ + // Check the filter. + auto check_result = check_sptl_fourier_filter(index); + auto err = check_result.first; + auto filt = check_result.second; + if (err) + return err; + + // Update the filter. + filt->set_order(order); + return 0; +} + +extern "C" int openmc_spatial_fourier_filter_set_params( + int32_t index, const int* axis, const double* min, const double* max) +{ + // Check the filter. + auto check_result = check_sptl_fourier_filter(index); + auto err = check_result.first; + auto filt = check_result.second; + if (err) + return err; + + // Update the filter. + if (axis) + filt->set_axis(static_cast(*axis)); + if (min && max) + filt->set_minmax(*min, *max); + return 0; +} + +} // namespace openmc diff --git a/src/tallies/filter_sptl_legendre.cpp b/src/tallies/filter_sptl_legendre.cpp index cf5ef2aed2d..ed0bcf972ec 100644 --- a/src/tallies/filter_sptl_legendre.cpp +++ b/src/tallies/filter_sptl_legendre.cpp @@ -7,35 +7,9 @@ #include "openmc/capi.h" #include "openmc/error.h" #include "openmc/math_functions.h" -#include "openmc/xml_interface.h" namespace openmc { -void SpatialLegendreFilter::from_xml(pugi::xml_node node) -{ - this->set_order(std::stoi(get_node_value(node, "order"))); - - auto axis = get_node_value(node, "axis"); - switch (axis[0]) { - case 'x': - this->set_axis(LegendreAxis::x); - break; - case 'y': - this->set_axis(LegendreAxis::y); - break; - case 'z': - this->set_axis(LegendreAxis::z); - break; - default: - throw std::runtime_error { - "Axis for SpatialLegendreFilter must be 'x', 'y', or 'z'"}; - } - - double min = std::stod(get_node_value(node, "min")); - double max = std::stod(get_node_value(node, "max")); - this->set_minmax(min, max); -} - void SpatialLegendreFilter::set_order(int order) { if (order < 0) { @@ -45,33 +19,11 @@ void SpatialLegendreFilter::set_order(int order) n_bins_ = order_ + 1; } -void SpatialLegendreFilter::set_axis(LegendreAxis axis) -{ - axis_ = axis; -} - -void SpatialLegendreFilter::set_minmax(double min, double max) -{ - if (max <= min) { - throw std::invalid_argument { - "Maximum value must be greater than minimum value"}; - } - min_ = min; - max_ = max; -} - void SpatialLegendreFilter::get_all_bins( const Particle& p, TallyEstimator estimator, FilterMatch& match) const { // Get the coordinate along the axis of interest. - double x; - if (axis_ == LegendreAxis::x) { - x = p.r().x; - } else if (axis_ == LegendreAxis::y) { - x = p.r().y; - } else { - x = p.r().z; - } + double x = this->position(p); if (x >= min_ && x <= max_) { // Compute the normalized coordinate value. @@ -87,30 +39,10 @@ void SpatialLegendreFilter::get_all_bins( } } -void SpatialLegendreFilter::to_statepoint(hid_t filter_group) const -{ - Filter::to_statepoint(filter_group); - write_dataset(filter_group, "order", order_); - if (axis_ == LegendreAxis::x) { - write_dataset(filter_group, "axis", "x"); - } else if (axis_ == LegendreAxis::y) { - write_dataset(filter_group, "axis", "y"); - } else { - write_dataset(filter_group, "axis", "z"); - } - write_dataset(filter_group, "min", min_); - write_dataset(filter_group, "max", max_); -} - std::string SpatialLegendreFilter::text_label(int bin) const { - if (axis_ == LegendreAxis::x) { - return fmt::format("Legendre expansion, x axis, P{}", bin); - } else if (axis_ == LegendreAxis::y) { - return fmt::format("Legendre expansion, y axis, P{}", bin); - } else { - return fmt::format("Legendre expansion, z axis, P{}", bin); - } + return fmt::format( + "Legendre expansion, {} axis, P{}", this->axis_label(), bin); } //============================================================================== @@ -196,7 +128,7 @@ extern "C" int openmc_spatial_legendre_filter_set_params( // Update the filter. if (axis) - filt->set_axis(static_cast(*axis)); + filt->set_axis(static_cast(*axis)); if (min && max) filt->set_minmax(*min, *max); return 0; diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 3fe48c1b021..a315a4444de 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -166,7 +166,8 @@ Tally::Tally(pugi::xml_node node) if (sf->cosine() == SphericalHarmonicsCosine::scatter) { estimator_ = TallyEstimator::ANALOG; } - } else if (filt_type == FilterType::SPATIAL_LEGENDRE || + } else if (filt_type == FilterType::SPATIAL_FOURIER || + filt_type == FilterType::SPATIAL_LEGENDRE || filt_type == FilterType::ZERNIKE || filt_type == FilterType::ZERNIKE_RADIAL) { estimator_ = TallyEstimator::COLLISION; diff --git a/tests/regression_tests/tallies/inputs_true.dat b/tests/regression_tests/tallies/inputs_true.dat index 40829f865a1..463f1d2f5e0 100644 --- a/tests/regression_tests/tallies/inputs_true.dat +++ b/tests/regression_tests/tallies/inputs_true.dat @@ -342,19 +342,31 @@ 4 - + 4 + x + -182.07 + 182.07 - + + 5 + x + -182.07 + 182.07 + + + 4 + + 1 2 3 4 6 8 - + 1 2 5 3 6 - + 10 21 22 23 60 - + 21 22 23 27 28 29 60 @@ -427,76 +439,86 @@ 11 - scatter nu-scatter flux total + flux analog - 11 + 12 + flux + analog + + + 13 + scatter nu-scatter flux total + analog + + + 13 flux total collision - - 11 + + 13 flux total tracklength - - 12 + + 14 total - - 15 + + 17 scatter - - 13 + + 15 absorption delayed-nu-fission events fission inverse-velocity kappa-fission (n,2n) (n,n1) (n,gamma) nu-fission scatter elastic total prompt-nu-fission fission-q-prompt fission-q-recoverable decay-rate tracklength - - 13 + + 15 U235 O16 total absorption delayed-nu-fission events fission inverse-velocity kappa-fission (n,2n) (n,n1) (n,gamma) nu-fission scatter elastic total prompt-nu-fission fission-q-prompt fission-q-recoverable decay-rate tracklength - - 13 + + 15 absorption delayed-nu-fission events fission inverse-velocity kappa-fission (n,2n) (n,n1) (n,gamma) nu-fission scatter elastic total prompt-nu-fission fission-q-prompt fission-q-recoverable decay-rate analog - - 13 + + 15 U235 O16 total absorption delayed-nu-fission events fission inverse-velocity kappa-fission (n,2n) (n,n1) (n,gamma) nu-fission scatter elastic total prompt-nu-fission fission-q-prompt fission-q-recoverable decay-rate analog - - 13 + + 15 absorption delayed-nu-fission events fission inverse-velocity kappa-fission (n,2n) (n,n1) (n,gamma) nu-fission scatter elastic total prompt-nu-fission fission-q-prompt fission-q-recoverable decay-rate collision - - 13 + + 15 U235 O16 total absorption delayed-nu-fission events fission inverse-velocity kappa-fission (n,2n) (n,n1) (n,gamma) nu-fission scatter elastic total prompt-nu-fission fission-q-prompt fission-q-recoverable decay-rate collision - - 14 + + 16 flux tracklength - - 14 + + 16 flux analog - - 14 + + 16 flux collision - + H1-production H2-production H3-production He3-production He4-production heating damage-energy diff --git a/tests/regression_tests/tallies/results_true.dat b/tests/regression_tests/tallies/results_true.dat index 1d3aca4b079..bf15a5303c7 100644 --- a/tests/regression_tests/tallies/results_true.dat +++ b/tests/regression_tests/tallies/results_true.dat @@ -1 +1 @@ -d01c3accd5b4de2aa166a77df28cfe42f5738a44c2480752fcfae7564a507362fff006b6dffb7b1dfe248e14bacef0070cabacea5d75c5996653e5605f7c7384 \ No newline at end of file +7cb0047ce4c66c9dfc563dafc95c88781fbfc7135a41ce22ea5028cfe0b1853d949e87d9b1698309dea8562495dd16334e1ea8b32f12524259fdd2d9219aff97 \ No newline at end of file diff --git a/tests/regression_tests/tallies/test.py b/tests/regression_tests/tallies/test.py index d20067ed33f..6e7b58d991a 100644 --- a/tests/regression_tests/tallies/test.py +++ b/tests/regression_tests/tallies/test.py @@ -105,6 +105,20 @@ def test_tallies(): legendre_tally.scores = ['scatter', 'nu-scatter'] legendre_tally.estimator = 'analog' + spatial_legendre_filter = SpatialLegendreFilter( + order=4, axis='x', minimum=-182.07, maximum=182.07) + spatial_legendre_tally = Tally() + spatial_legendre_tally.filters = [spatial_legendre_filter] + spatial_legendre_tally.scores = ['flux'] + spatial_legendre_tally.estimator = 'analog' + + spatial_fourier_filter = SpatialFourierFilter( + order=5, axis='x', minimum=-182.07, maximum=182.07) + spatial_fourier_tally = Tally() + spatial_fourier_tally.filters = [spatial_fourier_filter] + spatial_fourier_tally.scores = ['flux'] + spatial_fourier_tally.estimator = 'analog' + harmonics_filter = SphericalHarmonicsFilter(order=4) harmonics_tally = Tally() harmonics_tally.filters = [harmonics_filter] @@ -168,7 +182,8 @@ def test_tallies(): azimuthal_tally1, azimuthal_tally2, azimuthal_tally3, cellborn_tally, dg_tally, energy_tally, energyout_tally, transfer_tally, material_tally, mu_tally1, mu_tally2, - polar_tally1, polar_tally2, polar_tally3, legendre_tally, + polar_tally1, polar_tally2, polar_tally3, + legendre_tally, spatial_legendre_tally, spatial_fourier_tally, harmonics_tally, harmonics_tally2, harmonics_tally3, universe_tally, collision_tally] model.tallies += score_tallies diff --git a/tests/unit_tests/test_filters.py b/tests/unit_tests/test_filters.py index 2e7481c8767..d16f328f239 100644 --- a/tests/unit_tests/test_filters.py +++ b/tests/unit_tests/test_filters.py @@ -125,6 +125,36 @@ def test_spatial_legendre(): assert new_f.axis == f.axis +def test_spatial_fourier(): + n = 5 + axis = 'x' + f = openmc.SpatialFourierFilter(n, axis, -10., 10.) + assert f.order == n + assert f.axis == axis + assert f.minimum == -10. + assert f.maximum == 10. + assert f.bins[0] == 'a0 (constant)' + assert f.bins[-1] == 'b5 (sin)' + assert f.bins[-2] == 'a5 (cos)' + assert len(f.bins) == 2*n + 1 + + # Make sure __repr__ works + repr(f) + + # to_xml_element() + elem = f.to_xml_element() + assert elem.tag == 'filter' + assert elem.attrib['type'] == 'spatialfourier' + assert elem.find('order').text == str(n) + assert elem.find('axis').text == str(axis) + + # from_xml_element() + new_f = openmc.Filter.from_xml_element(elem) + assert new_f.id == f.id + assert new_f.order == f.order + assert new_f.axis == f.axis + + def test_spherical_harmonics(): n = 3 f = openmc.SphericalHarmonicsFilter(n)