diff --git a/CMakeLists.txt b/CMakeLists.txt index 0f4cc1b52..b4011434e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -416,6 +416,7 @@ list(APPEND libopenmc_SOURCES src/tallies/filter_meshborn.cpp src/tallies/filter_meshsurface.cpp src/tallies/filter_mu.cpp + src/tallies/filter_musurface.cpp src/tallies/filter_particle.cpp src/tallies/filter_polar.cpp src/tallies/filter_sph_harm.cpp diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index 5ae9f20ed..611d2c216 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -133,6 +133,7 @@ Constructing Tallies openmc.EnergyFilter openmc.EnergyoutFilter openmc.MuFilter + openmc.MuSurfaceFilter openmc.PolarFilter openmc.AzimuthalFilter openmc.DistribcellFilter diff --git a/include/openmc/tallies/filter.h b/include/openmc/tallies/filter.h index 1166c0eee..210ab284b 100644 --- a/include/openmc/tallies/filter.h +++ b/include/openmc/tallies/filter.h @@ -36,6 +36,7 @@ enum class FilterType { MESHBORN, MESH_SURFACE, MU, + MUSURFACE, PARTICLE, POLAR, SPHERICAL_HARMONICS, diff --git a/include/openmc/tallies/filter_mu.h b/include/openmc/tallies/filter_mu.h index 942ee60c2..b0ee40eac 100644 --- a/include/openmc/tallies/filter_mu.h +++ b/include/openmc/tallies/filter_mu.h @@ -40,7 +40,7 @@ public: void set_bins(gsl::span bins); -private: +protected: //---------------------------------------------------------------------------- // Data members diff --git a/include/openmc/tallies/filter_musurface.h b/include/openmc/tallies/filter_musurface.h new file mode 100644 index 000000000..2ca19e3a2 --- /dev/null +++ b/include/openmc/tallies/filter_musurface.h @@ -0,0 +1,34 @@ +#ifndef OPENMC_TALLIES_FILTER_MU_SURFACE_H +#define OPENMC_TALLIES_FILTER_MU_SURFACE_H + +#include + +#include "openmc/tallies/filter_mu.h" +#include "openmc/vector.h" + +namespace openmc { + +//============================================================================== +//! Bins the incoming-outgoing direction cosine. This is only used for surface +//! crossings. +//============================================================================== + +class MuSurfaceFilter : public MuFilter { +public: + //---------------------------------------------------------------------------- + // Constructors, destructors + + ~MuSurfaceFilter() = default; + + //---------------------------------------------------------------------------- + // Methods + + std::string type_str() const override { return "musurface"; } + FilterType type() const override { return FilterType::MUSURFACE; } + + void get_all_bins(const Particle& p, TallyEstimator estimator, + FilterMatch& match) const override; +}; + +} // namespace openmc +#endif // OPENMC_TALLIES_FILTER_MU_SURFACE_H diff --git a/openmc/filter.py b/openmc/filter.py index 0f884cfcd..b5926af5a 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -22,7 +22,7 @@ from ._xml import get_text _FILTER_TYPES = ( 'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', - 'energyout', 'mu', 'polar', 'azimuthal', 'distribcell', 'delayedgroup', + 'energyout', 'mu', 'musurface', 'polar', 'azimuthal', 'distribcell', 'delayedgroup', 'energyfunction', 'cellfrom', 'materialfrom', 'legendre', 'spatiallegendre', 'sphericalharmonics', 'zernike', 'zernikeradial', 'particle', 'cellinstance', 'collision', 'time' @@ -765,7 +765,7 @@ class ParticleFilter(Filter): @Filter.bins.setter def bins(self, bins): cv.check_type('bins', bins, Sequence, str) - bins = np.atleast_1d(bins) + bins = np.atleast_1d(bins) for edge in bins: cv.check_value('filter bin', edge, _PARTICLES) self._bins = bins @@ -1850,6 +1850,44 @@ class MuFilter(RealFilter): cv.check_less_than('filter value', x, 1., equality=True) +class MuSurfaceFilter(MuFilter): + """Bins tally events based on the angle of surface crossing. + + This filter bins events based on the cosine of the angle between the + direction of the particle and the normal to the surface at the point it + crosses. Only used in conjunction with a SurfaceFilter and current score. + + .. versionadded:: 0.15.1 + + Parameters + ---------- + values : int or Iterable of Real + A grid of surface crossing angles which the events will be divided into. + Values represent the cosine of the angle between the direction of the + particle and the normal to the surface at the point it crosses. If an + iterable is given, the values will be used explicitly as grid points. If + a single int is given, the range [-1, 1] will be divided equally into + that number of bins. + filter_id : int + Unique identifier for the filter + + Attributes + ---------- + values : numpy.ndarray + An array of values for which each successive pair constitutes a range of + surface crossing angle cosines for a single bin. + id : int + Unique identifier for the filter + bins : numpy.ndarray + An array of shape (N, 2) where each row is a pair of cosines of surface + crossing angle for a single filter + num_bins : Integral + The number of filter bins + + """ + # Note: inherits implementation from MuFilter + + class PolarFilter(RealFilter): """Bins tally events based on the incident particle's direction. diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py index 340c2fa34..b30f5e662 100644 --- a/openmc/lib/filter.py +++ b/openmc/lib/filter.py @@ -21,9 +21,9 @@ __all__ = [ 'CellInstanceFilter', 'CollisionFilter', 'DistribcellFilter', 'DelayedGroupFilter', 'EnergyFilter', 'EnergyoutFilter', 'EnergyFunctionFilter', 'LegendreFilter', 'MaterialFilter', 'MaterialFromFilter', 'MeshFilter', 'MeshBornFilter', - 'MeshSurfaceFilter', 'MuFilter', 'ParticleFilter', - 'PolarFilter', 'SphericalHarmonicsFilter', 'SpatialLegendreFilter', 'SurfaceFilter', - 'UniverseFilter', 'ZernikeFilter', 'ZernikeRadialFilter', 'filters' + 'MeshSurfaceFilter', 'MuFilter', 'MuSurfaceFilter', 'ParticleFilter', 'PolarFilter', + 'SphericalHarmonicsFilter', 'SpatialLegendreFilter', 'SurfaceFilter', 'UniverseFilter', + 'ZernikeFilter', 'ZernikeRadialFilter', 'filters' ] # Tally functions @@ -540,6 +540,10 @@ class MuFilter(Filter): filter_type = 'mu' +class MuSurfaceFilter(Filter): + filter_type = 'musurface' + + class ParticleFilter(Filter): filter_type = 'particle' @@ -642,6 +646,7 @@ _FILTER_TYPE_MAP = { 'meshborn': MeshBornFilter, 'meshsurface': MeshSurfaceFilter, 'mu': MuFilter, + 'musurface': MuSurfaceFilter, 'particle': ParticleFilter, 'polar': PolarFilter, 'sphericalharmonics': SphericalHarmonicsFilter, diff --git a/src/tallies/filter.cpp b/src/tallies/filter.cpp index ff7a3416b..074212db4 100644 --- a/src/tallies/filter.cpp +++ b/src/tallies/filter.cpp @@ -26,6 +26,7 @@ #include "openmc/tallies/filter_meshborn.h" #include "openmc/tallies/filter_meshsurface.h" #include "openmc/tallies/filter_mu.h" +#include "openmc/tallies/filter_musurface.h" #include "openmc/tallies/filter_particle.h" #include "openmc/tallies/filter_polar.h" #include "openmc/tallies/filter_sph_harm.h" @@ -133,6 +134,8 @@ Filter* Filter::create(const std::string& type, int32_t id) return Filter::create(id); } else if (type == "mu") { return Filter::create(id); + } else if (type == "musurface") { + return Filter::create(id); } else if (type == "particle") { return Filter::create(id); } else if (type == "polar") { diff --git a/src/tallies/filter_musurface.cpp b/src/tallies/filter_musurface.cpp new file mode 100644 index 000000000..58fe39b91 --- /dev/null +++ b/src/tallies/filter_musurface.cpp @@ -0,0 +1,36 @@ +#include "openmc/tallies/filter_musurface.h" + +#include // for abs, copysign + +#include "openmc/search.h" +#include "openmc/surface.h" +#include "openmc/tallies/tally_scoring.h" + +namespace openmc { + +void MuSurfaceFilter::get_all_bins( + const Particle& p, TallyEstimator estimator, FilterMatch& match) const +{ + // Get surface normal (and make sure it is a unit vector) + const auto surf {model::surfaces[std::abs(p.surface()) - 1].get()}; + auto n = surf->normal(p.r()); + n /= n.norm(); + + // Determine whether normal should be pointing in or out + if (p.surface() < 0) + n *= -1; + + // Determine cosine of angle between normal and particle direction + double mu = p.u().dot(n); + if (std::abs(mu) > 1.0) + mu = std::copysign(1.0, mu); + + // Find matching bin + if (mu >= bins_.front() && mu <= bins_.back()) { + auto bin = lower_bound_index(bins_.begin(), bins_.end(), mu); + match.bins_.push_back(bin); + match.weights_.push_back(1.0); + } +} + +} // namespace openmc diff --git a/tests/regression_tests/filter_musurface/__init__.py b/tests/regression_tests/filter_musurface/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/filter_musurface/inputs_true.dat b/tests/regression_tests/filter_musurface/inputs_true.dat new file mode 100644 index 000000000..031f62159 --- /dev/null +++ b/tests/regression_tests/filter_musurface/inputs_true.dat @@ -0,0 +1,37 @@ + + + + + + + + + + + + + + + + + + + + eigenvalue + 1000 + 5 + 0 + + + + 1 + + + -1.0 -0.5 0.0 0.5 1.0 + + + 1 2 + current + + + diff --git a/tests/regression_tests/filter_musurface/results_true.dat b/tests/regression_tests/filter_musurface/results_true.dat new file mode 100644 index 000000000..39c9f2b92 --- /dev/null +++ b/tests/regression_tests/filter_musurface/results_true.dat @@ -0,0 +1,11 @@ +k-combined: +1.157005E-01 7.587090E-03 +tally 1: +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +8.770000E-01 +1.608710E-01 +3.909000E+00 +3.063035E+00 diff --git a/tests/regression_tests/filter_musurface/test.py b/tests/regression_tests/filter_musurface/test.py new file mode 100644 index 000000000..f2ec96b49 --- /dev/null +++ b/tests/regression_tests/filter_musurface/test.py @@ -0,0 +1,43 @@ +import numpy as np +from math import pi + +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def model(): + model = openmc.Model() + fuel = openmc.Material() + fuel.set_density('g/cm3', 10.0) + fuel.add_nuclide('U235', 1.0) + zr = openmc.Material() + zr.set_density('g/cm3', 1.0) + zr.add_nuclide('Zr90', 1.0) + + cyl1 = openmc.ZCylinder(r=1.0) + cyl2 = openmc.ZCylinder(r=3.0, boundary_type='vacuum') + cell1 = openmc.Cell(fill=fuel, region=-cyl1) + cell2 = openmc.Cell(fill=zr, region=+cyl1 & -cyl2) + model.geometry = openmc.Geometry([cell1, cell2]) + + model.settings.batches = 5 + model.settings.inactive = 0 + model.settings.particles = 1000 + + # Create a tally for current through the first surface binned by mu + surf_filter = openmc.SurfaceFilter([cyl1]) + mu_filter = openmc.MuSurfaceFilter([-1.0, -0.5, 0.0, 0.5, 1.0]) + tally = openmc.Tally() + tally.filters = [surf_filter, mu_filter] + tally.scores = ['current'] + model.tallies.append(tally) + + return model + + +def test_filter_musurface(model): + harness = PyAPITestHarness('statepoint.5.h5', model) + harness.main() diff --git a/tests/unit_tests/test_filter_musurface.py b/tests/unit_tests/test_filter_musurface.py new file mode 100644 index 000000000..ca0db71f0 --- /dev/null +++ b/tests/unit_tests/test_filter_musurface.py @@ -0,0 +1,38 @@ +import openmc + + +def test_musurface(run_in_tmpdir): + sphere = openmc.Sphere(r=1.0, boundary_type='vacuum') + cell = openmc.Cell(region=-sphere) + model = openmc.Model() + model.geometry = openmc.Geometry([cell]) + model.settings.particles = 1000 + model.settings.batches = 10 + E = 1.0 + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Point(), + angle=openmc.stats.Isotropic(), + energy=openmc.stats.delta_function(E), + ) + model.settings.run_mode = "fixed source" + + filter1 = openmc.MuSurfaceFilter(200) + filter2 = openmc.SurfaceFilter(sphere) + tally = openmc.Tally() + tally.filters = [filter1, filter2] + tally.scores = ['current'] + model.tallies = openmc.Tallies([tally]) + + # Run OpenMC + sp_filename = model.run() + + # Get current binned by mu + with openmc.StatePoint(sp_filename) as sp: + current_mu = sp.tallies[tally.id].mean.ravel() + + # All contributions should show up in last bin + assert current_mu[-1] == 1.0 + for element in current_mu[:-1]: + assert element == 0.0 + +