diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index bbc816d7eda..01021e8d22d 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -634,6 +634,20 @@ found in the :ref:`random ray user guide `. *Default*: None + :source_shape: + Specifies the assumed shape of the source distribution within each + source region. Options are "flat", "linear", or "linear_xy". + + *Default*: flat + + :source_gradient_limiter: + Specifies whether to rescale linear source gradients as needed so that + the source shape modeled within each source region remains non-negative + over the region's bounding box, as sampled by the rays that have crossed + it (bool). Only used when the source shape is "linear" or "linear_xy". + + *Default*: false + :volume_normalized_flux_tallies: Specifies whether to normalize flux tallies by volume (bool). The default is 'False'. When enabled, flux tallies will be reported in units diff --git a/docs/source/methods/random_ray.rst b/docs/source/methods/random_ray.rst index c139f05d29a..8e2d848ac3a 100644 --- a/docs/source/methods/random_ray.rst +++ b/docs/source/methods/random_ray.rst @@ -1067,6 +1067,64 @@ accumulated centroid it feeds, so the omitted term has no persistent sign and its contribution to the accumulated flux shrinks with the number of batches, while the original form carries less variance there. +.. _methods_random_ray_gradient_limiter: + +~~~~~~~~~~~~~~~~~~~~~~~~ +Source Gradient Limiting +~~~~~~~~~~~~~~~~~~~~~~~~ + +The fitted source gradient :math:`\boldsymbol{\vec{Q}}_{i,g} = +\mathbf{M}_i^{-1} \boldsymbol{\vec{q}}_{i,g}` amplifies noise in the fitted +moments along any thin extent of a region, so a poorly sampled region can +carry a spuriously steep gradient and emit a negative source over part of +its extent. Rays crossing that part can carry negative angular flux +downstream, which optically thin media with scattering ratios near one can +amplify. + +When the source gradient limiter is enabled, each group's gradient is +rescaled so that the modeled source stays non-negative over the region's +axis-aligned bounding box. The box is accumulated from the endpoints of +every ray segment that has crossed the region past the ray's inactive +length. These lie on the region's boundary except where a ray starts or +ends inside it. The linear term is lowest at a corner of the box, where it +reaches + +.. math:: + :label: gradient-limiter-bound + + \sum_{d \in \{x, y, z\}} \; \min_{x_d \in \{x^{\min}_{i,d},\, + x^{\max}_{i,d}\}} \left(\boldsymbol{\vec{Q}}_{i,g}\right)_d \left(x_d - + r_{\mathrm{c},i,d}\right), + +where :math:`x^{\min}_{i}` and :math:`x^{\max}_{i}` are the box bounds, +:math:`\mathbf{r}_{\mathrm{c},i}` is the centroid, and :math:`d` indexes +their components. Whenever the flat source :math:`Q_{i,g}` plus this +minimum is negative, the gradient is scaled by the ratio of the flat source +to the magnitude of the minimum, so that the modeled source reaches zero at +that corner. Because the linear term integrates to zero over the region, +the rescaling preserves the region's mean emission, and gradients that pass +the test are left untouched. A group whose flat source is not positive has +its gradient zeroed. Once the region's extreme points along each axis have +been sampled, the box contains the region and the modeled source is +non-negative throughout it. The bound is exact for axis-aligned box regions +and conservative for others: a sphere is limited by up to a factor of +:math:`\sqrt{3}` more than necessary, and a thin region lying diagonally to +the axes by much more, as its bounding box is far larger than the region. + +This is the treatment `MPACT `_ applies in its limited linear +source approximation, with the same mean-preserving factor. MPACT finds +the minimum source exactly, over the entrance and exit points of every +segment crossing the region, which requires the fixed set of tracks that +deterministic MOC lays down once. Random ray samples new rays every batch, +so no such segment set exists when the source is built, and the sampled +bounding box takes its place. + +The limiter is off by default because a steep fit can also be physical, as +in the optically thick regions of deep-penetration problems, where +limiting discards real shape information and alters the solution at +depth. It is best reserved for simulations that negative sources +destabilize. + .. _methods-shannon-entropy-random-ray: ----------------------------- @@ -1229,6 +1287,7 @@ in random ray particle transport are: .. _Tramm-2020: https://doi.org/10.1051/EPJCONF/202124703021 .. _Cosgrove-2023: https://doi.org/10.1080/00295639.2023.2270618 .. _Ferrer-2016: https://doi.org/10.13182/NSE15-6 +.. _Choi-2024: https://doi.org/10.1080/00295639.2023.2224234 .. _Gunow-2018: https://dspace.mit.edu/handle/1721.1/119030 .. only:: html diff --git a/docs/source/usersguide/random_ray.rst b/docs/source/usersguide/random_ray.rst index 88b5c8fab1c..fabddcd3d95 100644 --- a/docs/source/usersguide/random_ray.rst +++ b/docs/source/usersguide/random_ray.rst @@ -979,6 +979,21 @@ in the :attr:`openmc.Settings.random_ray` dictionary to ``'linear'`` as:: LS enables the use of coarser mesh discretizations and lower ray populations, offsetting the increased computation per ray. +In poorly sampled source regions, fitted gradients can become spuriously +steep, producing negative sources that may destabilize optically thin, +scattering-dominated problems. If this occurs, a gradient limiter can be +enabled as:: + + settings.random_ray['source_gradient_limiter'] = True + +The limiter rescales a region's gradient as needed so that the modeled +source stays non-negative over the region's bounding box, as sampled by +the rays that have crossed it, preserving the region's mean emission. The +limiter is off by default, as limiting also clips physically steep source +shapes such as those found in optically thick regions of deep-penetration +problems; see the :ref:`methods documentation +` for details. + While OpenMC has no specific mode for 2D simulations, such simulations can be performed implicitly by leaving one of the dimensions of the geometry unbounded or by imposing reflective boundary conditions with no variation in between them diff --git a/include/openmc/bounding_box.h b/include/openmc/bounding_box.h index 4fabe1b7093..653b0ca184f 100644 --- a/include/openmc/bounding_box.h +++ b/include/openmc/bounding_box.h @@ -27,6 +27,23 @@ struct BoundingBox { return {{INFTY, INFTY, INFTY}, {-INFTY, -INFTY, -INFTY}}; } + //! Expand the bounding box to include a point. + inline void expand(const Position& p) + { + if (p.x < min.x) + min.x = p.x; + if (p.y < min.y) + min.y = p.y; + if (p.z < min.z) + min.z = p.z; + if (p.x > max.x) + max.x = p.x; + if (p.y > max.y) + max.y = p.y; + if (p.z > max.z) + max.z = p.z; + } + inline BoundingBox operator&(const BoundingBox& other) { BoundingBox result = *this; diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index c1f2fdf39a0..5b66e8e099f 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -122,6 +122,9 @@ class FlatSourceDomain { static bool volume_normalized_flux_tallies_; // If the user wants outputs based on the adjoint flux static bool adjoint_requested_; + // If the user wants linear source gradients rescaled so the modeled source + // stays non-negative over each source region + static bool source_gradient_limiter_; // The solve currently being executed static RandomRaySolve solve_; static bool fw_cadis_local_; diff --git a/include/openmc/random_ray/source_region.h b/include/openmc/random_ray/source_region.h index 82244c695e6..f9177536d8e 100644 --- a/include/openmc/random_ray/source_region.h +++ b/include/openmc/random_ray/source_region.h @@ -1,6 +1,7 @@ #ifndef OPENMC_RANDOM_RAY_SOURCE_REGION_H #define OPENMC_RANDOM_RAY_SOURCE_REGION_H +#include "openmc/bounding_box.h" #include "openmc/openmp_interface.h" #include "openmc/position.h" #include "openmc/random_ray/moment_matrix.h" @@ -167,6 +168,9 @@ class SourceRegionHandle { Position* centroid_t_; MomentMatrix* mom_matrix_; MomentMatrix* mom_matrix_t_; + // Bounding box of the ray segment endpoints sampled in this region, kept + // only when the source gradient limiter is enabled (see SourceRegion). + BoundingBox* extent_; // A set of volume tally tasks. This more complicated data structure is // convenient for ensuring that volumes are only tallied once per source // region, regardless of how many energy groups are used for tallying. @@ -259,6 +263,9 @@ class SourceRegionHandle { MomentMatrix& mom_matrix_t() { return *mom_matrix_t_; } const MomentMatrix mom_matrix_t() const { return *mom_matrix_t_; } + BoundingBox& extent() { return *extent_; } + const BoundingBox& extent() const { return *extent_; } + std::unordered_set& volume_task() { return *volume_task_; @@ -372,6 +379,14 @@ class SourceRegion { MomentMatrix mom_matrix_t_ {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}; //!< The spatial moment matrix accumulated over all iterations + // Bounding box of the ray segment endpoints sampled in this region. Segment + // endpoints lie on the region boundary, so the box converges to the + // region's true extent. It is accumulated only when the source gradient + // limiter is enabled, which bounds the linear source over it. It starts + // inverted (minimum above maximum), the empty box, rather than at the + // default infinite box. + BoundingBox extent_ {BoundingBox::inverted()}; + // A set of volume tally tasks. This more complicated data structure is // convenient for ensuring that volumes are only tallied once per source // region, regardless of how many energy groups are used for tallying. @@ -411,10 +426,10 @@ class SourceRegionContainer { public: //---------------------------------------------------------------------------- // Constructors - SourceRegionContainer( - int negroups, bool is_linear, bool is_adaptive, bool is_strict_adaptive) + SourceRegionContainer(int negroups, bool is_linear, bool is_adaptive, + bool is_strict_adaptive, bool track_extents) : negroups_(negroups), is_linear_(is_linear), is_adaptive_(is_adaptive), - is_strict_adaptive_(is_strict_adaptive) + is_strict_adaptive_(is_strict_adaptive), track_extents_(track_extents) {} SourceRegionContainer() = default; @@ -504,6 +519,9 @@ class SourceRegionContainer { return mom_matrix_t_[sr]; } + BoundingBox& extent(int64_t sr) { return extents_[sr]; } + const BoundingBox& extent(int64_t sr) const { return extents_[sr]; } + MomentArray& source_gradients(int64_t sr, int g) { return source_gradients_[index(sr, g)]; @@ -682,6 +700,9 @@ class SourceRegionContainer { bool is_linear_ {false}; bool is_adaptive_ {false}; bool is_strict_adaptive_ {false}; + // Whether the sampled bounding boxes are stored (linear source with the + // source gradient limiter enabled) + bool track_extents_ {false}; // SoA storage for scalar fields (one item per source region) vector material_; @@ -712,6 +733,9 @@ class SourceRegionContainer { vector centroid_offset_; vector mom_matrix_; vector mom_matrix_t_; + // One box per region rather than separate minimum and maximum arrays: the + // two corners are always read, grown, and reset together. + vector extents_; // A set of volume tally tasks. This more complicated data structure is // convenient for ensuring that volumes are only tallied once per source // region, regardless of how many energy groups are used for tallying. diff --git a/openmc/settings.py b/openmc/settings.py index 100c283873e..5d71df091c5 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -210,6 +210,12 @@ class Settings: :source_shape: Assumed shape of the source distribution within each source region. Options are 'flat' (default), 'linear', or 'linear_xy'. + :source_gradient_limiter: + Whether to rescale linear source gradients as needed so that the + source shape modeled within each source region remains + non-negative over the region's bounding box, as sampled by the + rays that have crossed it (bool). The default is 'False'. Only + used when the source shape is 'linear' or 'linear_xy'. :volume_normalized_flux_tallies: Whether to normalize flux tallies by volume (bool). The default is 'False'. When enabled, flux tallies will be reported in units of @@ -1436,6 +1442,8 @@ def random_ray(self, random_ray: dict): ('flat', 'linear', 'linear_xy')) elif key == 'volume_normalized_flux_tallies': cv.check_type('volume normalized flux tallies', value, bool) + elif key == 'source_gradient_limiter': + cv.check_type('source gradient limiter', value, bool) elif key == 'adjoint': cv.check_type('adjoint', value, bool) elif key == 'source_region_meshes': @@ -2529,6 +2537,10 @@ def _random_ray_from_xml_element(self, root, meshes=None): self.random_ray['adjoint'] = ( child.text in ('true', '1') ) + elif child.tag == 'source_gradient_limiter': + self.random_ray['source_gradient_limiter'] = ( + child.text in ('true', '1') + ) elif child.tag == 'adjoint_source': self.random_ray['adjoint_source'] = [] for subelem in child.findall('source'): diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index d2f13f6b7d9..fc3320dcaff 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -33,6 +33,7 @@ RandomRayVolumeEstimator FlatSourceDomain::resolved_volume_estimator_ { RandomRayVolumeEstimator::AUTO}; bool FlatSourceDomain::volume_normalized_flux_tallies_ {false}; bool FlatSourceDomain::adjoint_requested_ {false}; +bool FlatSourceDomain::source_gradient_limiter_ {false}; RandomRaySolve FlatSourceDomain::solve_ {RandomRaySolve::FORWARD}; bool FlatSourceDomain::fw_cadis_local_ {false}; double FlatSourceDomain::diagonal_stabilization_rho_ {1.0}; @@ -61,8 +62,9 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_) bool is_adaptive = is_adaptive_family(resolved_volume_estimator_); bool is_strict_adaptive = resolved_volume_estimator_ == RandomRayVolumeEstimator::STRICT_ADAPTIVE; - source_regions_ = SourceRegionContainer( - negroups_, is_linear, is_adaptive, is_strict_adaptive); + // The sampled bounding boxes exist only for the source gradient limiter + source_regions_ = SourceRegionContainer(negroups_, is_linear, is_adaptive, + is_strict_adaptive, is_linear && source_gradient_limiter_); // Initialize tally volumes if (volume_normalized_flux_tallies_) { diff --git a/src/random_ray/linear_source_domain.cpp b/src/random_ray/linear_source_domain.cpp index 3d045dca85c..e61634e96c5 100644 --- a/src/random_ray/linear_source_domain.cpp +++ b/src/random_ray/linear_source_domain.cpp @@ -1,5 +1,7 @@ #include "openmc/random_ray/linear_source_domain.h" +#include + #include "openmc/cell.h" #include "openmc/geometry.h" #include "openmc/material.h" @@ -133,6 +135,40 @@ void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh) srh.source_gradients(g) = {0.0, 0.0, 0.0}; } } + + // If enabled by the user, limit the source gradients so the modeled local + // source q(r) = q_flat + (r - centroid) . q_gradient stays non-negative + // over the region's bounding box as sampled by the ray segment endpoints. + // The linear term is lowest at the box corner each gradient component + // points away from, so its minimum is the sum, over the three axes, of the + // gradient component times the offset from the centroid to that face. + // Once the region's extreme points along each axis have been sampled the + // box contains the region, and the modeled source is non-negative + // throughout it whenever the flat source covers the dip. When it does not, + // the gradient is scaled by their ratio, which preserves the region's mean + // emission, since the linear term integrates to zero over the region; + // gradients that pass are left untouched. A non-positive flat source + // leaves no shape to keep, so its cap is zero and its gradient is scaled + // away. A region with no sampled box yet carries no gradient to limit. + if (source_gradient_limiter_ && material != MATERIAL_VOID && + srh.extent().min.x <= srh.extent().max.x) { + // Offsets from the centroid to the box faces. The centroid is the + // length-weighted mean of segment midpoints, all of which lie in the + // box, so lo <= 0 <= hi and the dip below is non-negative. + const BoundingBox& extent = srh.extent(); + Position lo = extent.min - srh.centroid(); + Position hi = extent.max - srh.centroid(); + for (int g = 0; g < negroups_; g++) { + MomentArray& gradient = srh.source_gradients(g); + double cap = std::max(srh.source(g), 0.0); + double dip = std::max(-gradient.x * lo.x, -gradient.x * hi.x) + + std::max(-gradient.y * lo.y, -gradient.y * hi.y) + + std::max(-gradient.z * lo.z, -gradient.z * hi.z); + if (dip > cap) { + gradient *= cap / dip; + } + } + } } void LinearSourceDomain::normalize_scalar_flux_and_volumes( diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index dde5023e44f..b5a4e9dd933 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -642,6 +642,14 @@ void RandomRay::attenuate_flux_linear_source( moment_matrix_estimate *= distance; srh.mom_matrix() += moment_matrix_estimate; + // With the source gradient limiter enabled, grow the region's sampled + // bounding box with this segment's endpoints, which lie on the region + // boundary (or inside it, where the ray starts or ends). + if (FlatSourceDomain::source_gradient_limiter_) { + srh.extent().expand(r); + srh.extent().expand(r + distance * u()); + } + srh.n_hits() += 1; } diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 29d98bf1dbd..bc9322e5d45 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -291,6 +291,7 @@ void openmc_finalize_random_ray() FlatSourceDomain::resolved_volume_estimator_ = RandomRayVolumeEstimator::AUTO; FlatSourceDomain::volume_normalized_flux_tallies_ = false; FlatSourceDomain::adjoint_requested_ = false; + FlatSourceDomain::source_gradient_limiter_ = false; FlatSourceDomain::solve_ = RandomRaySolve::FORWARD; FlatSourceDomain::fw_cadis_local_ = false; FlatSourceDomain::fw_cadis_local_targets_.clear(); @@ -693,6 +694,10 @@ void RandomRaySimulation::print_results_random_ray( fatal_error("Invalid random ray source shape"); } fmt::print(" Source Shape = {}\n", shape); + if (RandomRay::source_shape_ != RandomRaySourceShape::FLAT) { + fmt::print(" Source Gradient Limiter = {}\n", + FlatSourceDomain::source_gradient_limiter_ ? "ON" : "OFF"); + } std::string sample_method; switch (RandomRay::sample_method_) { case RandomRaySampleMethod::PRNG: diff --git a/src/random_ray/source_region.cpp b/src/random_ray/source_region.cpp index 16598fd4b3a..70b23172505 100644 --- a/src/random_ray/source_region.cpp +++ b/src/random_ray/source_region.cpp @@ -22,7 +22,7 @@ SourceRegionHandle::SourceRegionHandle(SourceRegion& sr) position_(&sr.position_), centroid_(&sr.centroid_), centroid_iteration_(&sr.centroid_iteration_), centroid_t_(&sr.centroid_t_), mom_matrix_(&sr.mom_matrix_), mom_matrix_t_(&sr.mom_matrix_t_), - volume_task_(&sr.volume_task_), mesh_(&sr.mesh_), + extent_(&sr.extent_), volume_task_(&sr.volume_task_), mesh_(&sr.mesh_), parent_sr_(&sr.parent_sr_), scalar_flux_old_(sr.scalar_flux_old_.data()), scalar_flux_new_(sr.scalar_flux_new_.data()), source_(sr.source_.data()), external_source_(sr.external_source_.data()), @@ -104,6 +104,9 @@ void SourceRegionContainer::push_back(const SourceRegion& sr) mom_matrix_.push_back(sr.mom_matrix_); mom_matrix_t_.push_back(sr.mom_matrix_t_); } + if (track_extents_) { + extents_.push_back(sr.extent_); + } // Energy-dependent fields for (int g = 0; g < negroups_; ++g) { @@ -164,6 +167,7 @@ void SourceRegionContainer::assign( mom_matrix_.clear(); mom_matrix_t_.clear(); } + extents_.clear(); scalar_flux_old_.clear(); scalar_flux_new_.clear(); @@ -238,6 +242,7 @@ SourceRegionHandle SourceRegionContainer::get_source_region_handle(int64_t sr) handle.centroid_t_ = ¢roid_t(sr); handle.mom_matrix_ = &mom_matrix(sr); handle.mom_matrix_t_ = &mom_matrix_t(sr); + handle.extent_ = track_extents_ ? &extent(sr) : nullptr; handle.source_gradients_ = &source_gradients(sr, 0); handle.flux_moments_old_ = &flux_moments_old(sr, 0); handle.flux_moments_new_ = &flux_moments_new(sr, 0); @@ -270,6 +275,8 @@ void SourceRegionContainer::adjoint_reset() MomentMatrix {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}); std::fill(mom_matrix_t_.begin(), mom_matrix_t_.end(), MomentMatrix {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}); + // The sampled bounding boxes are re-accumulated alongside the centroids + std::fill(extents_.begin(), extents_.end(), BoundingBox::inverted()); if (settings::run_mode == RunMode::FIXED_SOURCE) { std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 0.0); } else { diff --git a/src/settings.cpp b/src/settings.cpp index 8ffec4a5355..5e2da3a52e4 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -341,6 +341,10 @@ void get_run_parameters(pugi::xml_node node_base) FlatSourceDomain::adjoint_requested_ = get_node_value_bool(random_ray_node, "adjoint"); } + if (check_for_node(random_ray_node, "source_gradient_limiter")) { + FlatSourceDomain::source_gradient_limiter_ = + get_node_value_bool(random_ray_node, "source_gradient_limiter"); + } if (check_for_node(random_ray_node, "sample_method")) { std::string temp_str = get_node_value(random_ray_node, "sample_method", true, true); diff --git a/tests/cpp_unit_tests/test_surface.cpp b/tests/cpp_unit_tests/test_surface.cpp index 1529cf21031..db5e7261a73 100644 --- a/tests/cpp_unit_tests/test_surface.cpp +++ b/tests/cpp_unit_tests/test_surface.cpp @@ -26,6 +26,27 @@ std::unique_ptr make_surface( } // anonymous namespace +TEST_CASE("Expand bounding box to include points") +{ + BoundingBox bbox = BoundingBox::inverted(); + + bbox.expand({1.0, -2.0, 3.0}); + CHECK(bbox.min.x == 1.0); + CHECK(bbox.min.y == -2.0); + CHECK(bbox.min.z == 3.0); + CHECK(bbox.max.x == 1.0); + CHECK(bbox.max.y == -2.0); + CHECK(bbox.max.z == 3.0); + + bbox.expand({-4.0, 0.0, 2.0}); + CHECK(bbox.min.x == -4.0); + CHECK(bbox.min.y == -2.0); + CHECK(bbox.min.z == 2.0); + CHECK(bbox.max.x == 1.0); + CHECK(bbox.max.y == 0.0); + CHECK(bbox.max.z == 3.0); +} + TEST_CASE("General plane bounding box") { pugi::xml_document doc; diff --git a/tests/regression_tests/random_ray_source_gradient_limiter/__init__.py b/tests/regression_tests/random_ray_source_gradient_limiter/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/random_ray_source_gradient_limiter/inputs_true.dat b/tests/regression_tests/random_ray_source_gradient_limiter/inputs_true.dat new file mode 100644 index 00000000000..977889ab31c --- /dev/null +++ b/tests/regression_tests/random_ray_source_gradient_limiter/inputs_true.dat @@ -0,0 +1,98 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 90 + 60 + 30 + + + 100.0 1.0 + + + material + 1 + + + multi-group + + 500.0 + 100.0 + + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + + true + linear + true + naive + + + + + + + + 12 12 12 + 0.0 0.0 0.0 + 30.0 30.0 30.0 + + + + + 1 + + + 2 + + + 3 + + + 3 + flux + tracklength + + + 2 + flux + tracklength + + + 1 + flux + tracklength + + + diff --git a/tests/regression_tests/random_ray_source_gradient_limiter/results_true.dat b/tests/regression_tests/random_ray_source_gradient_limiter/results_true.dat new file mode 100644 index 00000000000..598775f4873 --- /dev/null +++ b/tests/regression_tests/random_ray_source_gradient_limiter/results_true.dat @@ -0,0 +1,9 @@ +tally 1: +5.862734E+00 +1.147901E+00 +tally 2: +9.701948E-01 +3.234962E-02 +tally 3: +5.323626E-03 +9.818594E-07 diff --git a/tests/regression_tests/random_ray_source_gradient_limiter/test.py b/tests/regression_tests/random_ray_source_gradient_limiter/test.py new file mode 100644 index 00000000000..92968b2a2a1 --- /dev/null +++ b/tests/regression_tests/random_ray_source_gradient_limiter/test.py @@ -0,0 +1,59 @@ +import os + +import openmc +from openmc.examples import random_ray_three_region_cube + +from tests.testing_harness import TolerantPyAPITestHarness + + +class MGXSTestHarness(TolerantPyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = 'mgxs.h5' + if os.path.exists(f): + os.remove(f) + + +def test_random_ray_source_gradient_limiter(): + # A linear source run with the gradient limiter enabled and firing in + # both of its regimes: the naive volume estimator and an overlay + # source-region mesh leave the example's optically thin interior with + # noisy fitted gradients, and the absorber's steep attenuation over + # regions a few mean free paths thick gives physically steep ones, which + # the limiter clips as well. The example's three cubic regions are + # replaced by spherical ones so that the curved boundaries cut the mesh + # cells into pieces whose centroids sit off-center in their bounding + # boxes, which is where the limiter's bound differs from a symmetric one. + openmc.reset_auto_ids() + model = random_ray_three_region_cube() + source_mat, void_mat, absorber_mat = model.materials + width = 30.0 + x0 = openmc.XPlane(0.0, boundary_type='reflective') + y0 = openmc.YPlane(0.0, boundary_type='reflective') + z0 = openmc.ZPlane(0.0, boundary_type='reflective') + x1 = openmc.XPlane(width, boundary_type='vacuum') + y1 = openmc.YPlane(width, boundary_type='vacuum') + z1 = openmc.ZPlane(width, boundary_type='vacuum') + domain = +x0 & -x1 & +y0 & -y1 & +z0 & -z1 + source_sphere = openmc.Sphere(r=5.0) + void_sphere = openmc.Sphere(r=12.5) + model.geometry = openmc.Geometry([ + openmc.Cell(fill=source_mat, region=-source_sphere & domain), + openmc.Cell(fill=void_mat, + region=+source_sphere & -void_sphere & domain), + openmc.Cell(fill=absorber_mat, region=+void_sphere & domain), + ]) + model.settings.source[0].constraints = {'domains': [source_mat]} + model.settings.random_ray['source_shape'] = 'linear' + model.settings.random_ray['source_gradient_limiter'] = True + model.settings.random_ray['volume_estimator'] = 'naive' + mesh = openmc.RegularMesh() + mesh.lower_left = (0.0, 0.0, 0.0) + mesh.upper_right = (width, width, width) + mesh.dimension = (12, 12, 12) + model.settings.random_ray['source_region_meshes'] = [ + (mesh, [model.geometry.root_universe])] + model.settings.inactive = 30 + model.settings.batches = 60 + harness = MGXSTestHarness('statepoint.60.h5', model) + harness.main() diff --git a/tests/unit_tests/test_settings.py b/tests/unit_tests/test_settings.py index bdb3ea8fe9f..12195b876a1 100644 --- a/tests/unit_tests/test_settings.py +++ b/tests/unit_tests/test_settings.py @@ -86,6 +86,7 @@ def test_export_to_xml(run_in_tmpdir): 'source_region_meshes': [(source_region_mesh, [root_universe])], 'volume_estimator': 'hybrid', 'source_shape': 'linear', + 'source_gradient_limiter': True, 'volume_normalized_flux_tallies': True, 'adjoint': False, 'sample_method': 'halton' @@ -184,6 +185,7 @@ def test_export_to_xml(run_in_tmpdir): assert recovered_mesh.upper_right == [2., 2., 2.] assert s.random_ray['volume_estimator'] == 'hybrid' assert s.random_ray['source_shape'] == 'linear' + assert s.random_ray['source_gradient_limiter'] assert s.random_ray['volume_normalized_flux_tallies'] assert not s.random_ray['adjoint'] assert s.random_ray['sample_method'] == 'halton'