diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index dfab96131df..8348d21312c 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 537930096d8..23d2af8eab4 100644 --- a/docs/source/methods/random_ray.rst +++ b/docs/source/methods/random_ray.rst @@ -1034,6 +1034,64 @@ The contents of this section, alongside the equations for the flat source and scalar flux, Equations :eq:`source_update` and :eq:`phi_sim` respectively, completes the set of equations for LS. +.. _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: ----------------------------- @@ -1196,6 +1254,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/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index c87ff48f301..eaed1af0207 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 f488aa293aa..3aca64ed126 100644 --- a/include/openmc/random_ray/source_region.h +++ b/include/openmc/random_ray/source_region.h @@ -167,6 +167,10 @@ 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). + Position* extent_min_; + Position* extent_max_; // 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,27 @@ class SourceRegionHandle { MomentMatrix& mom_matrix_t() { return *mom_matrix_t_; } const MomentMatrix mom_matrix_t() const { return *mom_matrix_t_; } + const Position extent_min() const { return *extent_min_; } + const Position extent_max() const { return *extent_max_; } + + // Grows the sampled bounding box to include a point + void expand_extent(const Position& p) + { + // Conditional stores: once the box has converged, no write is made + if (p.x < extent_min_->x) + extent_min_->x = p.x; + if (p.y < extent_min_->y) + extent_min_->y = p.y; + if (p.z < extent_min_->z) + extent_min_->z = p.z; + if (p.x > extent_max_->x) + extent_max_->x = p.x; + if (p.y > extent_max_->y) + extent_max_->y = p.y; + if (p.z > extent_max_->z) + extent_max_->z = p.z; + } + std::unordered_set& volume_task() { return *volume_task_; @@ -372,6 +397,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. The empty + // box has its minimum above its maximum. + Position extent_min_ {INFTY, INFTY, INFTY}; + Position extent_max_ {-INFTY, -INFTY, -INFTY}; + // 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 +444,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; @@ -498,6 +531,12 @@ class SourceRegionContainer { return mom_matrix_t_[sr]; } + Position& extent_min(int64_t sr) { return extent_min_[sr]; } + const Position extent_min(int64_t sr) const { return extent_min_[sr]; } + + Position& extent_max(int64_t sr) { return extent_max_[sr]; } + const Position extent_max(int64_t sr) const { return extent_max_[sr]; } + MomentArray& source_gradients(int64_t sr, int g) { return source_gradients_[index(sr, g)]; @@ -676,6 +715,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_; @@ -701,6 +743,8 @@ class SourceRegionContainer { vector centroid_t_; vector mom_matrix_; vector mom_matrix_t_; + vector extent_min_; + vector extent_max_; // 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 78bed32b8ff..87c75a16a70 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 @@ -1426,6 +1432,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': @@ -2519,6 +2527,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 f3567e6f9db..98ab8542e30 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 229752ca7f9..24800313b86 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,39 @@ 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. + Position lo = srh.extent_min() - srh.centroid(); + Position hi = srh.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..90a68733e1d 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.expand_extent(r); + srh.expand_extent(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 ab8126f1804..fc688f8941d 100644 --- a/src/random_ray/source_region.cpp +++ b/src/random_ray/source_region.cpp @@ -22,6 +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_), + extent_min_(&sr.extent_min_), extent_max_(&sr.extent_max_), 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()), @@ -103,6 +104,10 @@ 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_) { + extent_min_.push_back(sr.extent_min_); + extent_max_.push_back(sr.extent_max_); + } // Energy-dependent fields for (int g = 0; g < negroups_; ++g) { @@ -162,6 +167,8 @@ void SourceRegionContainer::assign( mom_matrix_.clear(); mom_matrix_t_.clear(); } + extent_min_.clear(); + extent_max_.clear(); scalar_flux_old_.clear(); scalar_flux_new_.clear(); @@ -236,6 +243,13 @@ 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); + if (track_extents_) { + handle.extent_min_ = &extent_min(sr); + handle.extent_max_ = &extent_max(sr); + } else { + handle.extent_min_ = nullptr; + handle.extent_max_ = 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); @@ -266,6 +280,11 @@ 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( + extent_min_.begin(), extent_min_.end(), Position {INFTY, INFTY, INFTY}); + std::fill( + extent_max_.begin(), extent_max_.end(), Position {-INFTY, -INFTY, -INFTY}); 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 85d1f428b2c..60f9ad192a9 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -342,6 +342,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/regression_tests/random_ray_linear_source_stability/__init__.py b/tests/regression_tests/random_ray_linear_source_stability/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/random_ray_linear_source_stability/inputs_true.dat b/tests/regression_tests/random_ray_linear_source_stability/inputs_true.dat new file mode 100644 index 00000000000..977889ab31c --- /dev/null +++ b/tests/regression_tests/random_ray_linear_source_stability/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_linear_source_stability/results_true.dat b/tests/regression_tests/random_ray_linear_source_stability/results_true.dat new file mode 100644 index 00000000000..d9e0f0ef4a8 --- /dev/null +++ b/tests/regression_tests/random_ray_linear_source_stability/results_true.dat @@ -0,0 +1,9 @@ +tally 1: +5.965230E+00 +1.205059E+00 +tally 2: +9.699577E-01 +3.233890E-02 +tally 3: +5.324099E-03 +9.829519E-07 diff --git a/tests/regression_tests/random_ray_linear_source_stability/test.py b/tests/regression_tests/random_ray_linear_source_stability/test.py new file mode 100644 index 00000000000..cfb0284fd7d --- /dev/null +++ b/tests/regression_tests/random_ray_linear_source_stability/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_linear_source_stability(): + # 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'