diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 114531bf754..47f52af6ca6 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -8,6 +8,7 @@ #include "openmc/tensor.h" +#include "openmc/constants.h" #include "openmc/error.h" #include "openmc/file_utils.h" #include "openmc/hdf5_interface.h" @@ -289,9 +290,13 @@ std::pair WeightWindows::get_weight_window( if (E < energy_bounds_.front() || E > energy_bounds_.back()) return {false, {}}; - // Get mesh index for particle's position + // Get mesh index for particle's position. The position is nudged along the + // direction of travel so that a particle sitting exactly on a mesh boundary + // (e.g., at a surface checkpoint where a geometry surface coincides with a + // mesh plane) is assigned to the element it is entering rather than always + // to the element on the lower-coordinate side. const auto& mesh = this->mesh(); - int mesh_bin = mesh->get_bin(p.r()); + int mesh_bin = mesh->get_bin(p.r() + TINY_BIT * p.u()); // particle is outside the weight window mesh if (mesh_bin < 0) diff --git a/tests/regression_tests/random_ray_auto_convert_bootstrap/results_true.dat b/tests/regression_tests/random_ray_auto_convert_bootstrap/results_true.dat index c32200aeff1..92dd0f70bce 100644 --- a/tests/regression_tests/random_ray_auto_convert_bootstrap/results_true.dat +++ b/tests/regression_tests/random_ray_auto_convert_bootstrap/results_true.dat @@ -1,9 +1,9 @@ tally 1: -8.094106E+01 -6.974487E+02 -3.181906E-04 -1.455701E-08 -9.840422E+02 -1.110377E+05 -8.515372E+01 -8.078508E+02 +8.087496E+01 +7.339588E+02 +2.446860E-04 +8.405655E-09 +9.761237E+02 +1.018073E+05 +9.815520E+01 +1.065938E+03 diff --git a/tests/regression_tests/weightwindows/survival_biasing/local/results_true.dat b/tests/regression_tests/weightwindows/survival_biasing/local/results_true.dat index 8a4d3368cc4..b8f32dc2f9f 100644 --- a/tests/regression_tests/weightwindows/survival_biasing/local/results_true.dat +++ b/tests/regression_tests/weightwindows/survival_biasing/local/results_true.dat @@ -1 +1 @@ -0d7b17d4e364bda2be21feb2e5b2aae4be69fc8ed634c4f45062e8efc61b7bd24dcf448837692fd603afe3756501fe8821d134f2d6289179618f85178ddb8793 \ No newline at end of file +7a9d895b430f0c7696fc36f73dc68179d484664a29bd168df3b28a2b0549552d6216e4e9b32de4449d8134a42a84883d7e7e7364f0b9e4cc75174e95ecaf6b1c \ No newline at end of file diff --git a/tests/regression_tests/weightwindows/survival_biasing/shared/results_true.dat b/tests/regression_tests/weightwindows/survival_biasing/shared/results_true.dat index 2e88e2f0e81..f99e3c14d80 100644 --- a/tests/regression_tests/weightwindows/survival_biasing/shared/results_true.dat +++ b/tests/regression_tests/weightwindows/survival_biasing/shared/results_true.dat @@ -1 +1 @@ -2f1c8efbf6ab0b0c7319fadd7ee8f27b545c1ebab94ebe15cf1133b54227487542966b6d798970c8b37bc19e7041e687687d5285803241bb83f4c8dc40a4d276 \ No newline at end of file +5ee91c9d9285829528fb2e53a1c5bf069318259e0d4e8768852309a921bc6dd77852de251bd6e91245f68bc33b0fdc2321cfc3b315432d822603976f74079396 \ No newline at end of file diff --git a/tests/unit_tests/weightwindows/test_ww_surface_checkpoint.py b/tests/unit_tests/weightwindows/test_ww_surface_checkpoint.py new file mode 100644 index 00000000000..911de6ade2a --- /dev/null +++ b/tests/unit_tests/weightwindows/test_ww_surface_checkpoint.py @@ -0,0 +1,91 @@ +import openmc +import pytest + + +@pytest.mark.parametrize("direction", [1.0, -1.0]) +def test_ww_surface_on_mesh_boundary(run_in_tmpdir, direction): + """Weight windows applied at a surface that coincides with a weight window + mesh boundary must use the mesh element the particle is entering, + regardless of the direction of travel. + + At the surface checkpoint the particle sits exactly on the plane x = 0. + Structured mesh index lookups resolve a position exactly on an element + boundary to the element on the lower-coordinate side, so when the window + was looked up at the particle's position alone, a particle moving in +x + was given the window of the element it was leaving. It was then not split + and the +x case failed, while the -x case passed only because there the + lower-coordinate element happens to be the one being entered. The lookup + now nudges the position along the direction of travel. + """ + + # Void-like one-group material so the only weight window checks happen at + # surface crossings + groups = openmc.mgxs.EnergyGroups([0.0, 20.0e6]) + xsdata = openmc.XSdata('void', groups) + xsdata.order = 0 + xsdata.set_total([0.0]) + xsdata.set_absorption([0.0]) + xsdata.set_scatter_matrix([[[0.0]]]) + mg_library = openmc.MGXSLibrary(groups) + mg_library.add_xsdata(xsdata) + mg_library.export_to_hdf5('mgxs.h5') + + mat = openmc.Material() + mat.add_macroscopic('void') + materials = openmc.Materials([mat]) + materials.cross_sections = 'mgxs.h5' + + # Two cells separated by the plane x = 0, which is also a mesh boundary + x_min = openmc.XPlane(-1.0, boundary_type='vacuum') + x_mid = openmc.XPlane(0.0) + x_max = openmc.XPlane(1.0, boundary_type='vacuum') + yz = openmc.model.RectangularPrism(2.0, 2.0, axis='x', + boundary_type='vacuum') + left = openmc.Cell(fill=mat, region=+x_min & -x_mid & -yz) + right = openmc.Cell(fill=mat, region=+x_mid & -x_max & -yz) + geometry = openmc.Geometry([left, right]) + + # Particle is born in the window of its birth element and must be split + # five ways when it enters the element on the other side of x = 0 + birth_cell, target_cell = (left, right) if direction > 0 else (right, left) + mesh = openmc.RectilinearMesh() + mesh.x_grid = [-1.0, 0.0, 1.0] + mesh.y_grid = [-1.0, 1.0] + mesh.z_grid = [-1.0, 1.0] + birth_bounds = (0.5, 1.5) + target_bounds = (0.1, 0.2) + if direction > 0: + lower = [birth_bounds[0], target_bounds[0]] + upper = [birth_bounds[1], target_bounds[1]] + else: + lower = [target_bounds[0], birth_bounds[0]] + upper = [target_bounds[1], birth_bounds[1]] + ww = openmc.WeightWindows(mesh, lower_ww_bounds=lower, + upper_ww_bounds=upper) + + settings = openmc.Settings() + settings.energy_mode = 'multi-group' + settings.run_mode = 'fixed source' + settings.particles = 10 + settings.batches = 1 + settings.source = openmc.IndependentSource( + space=openmc.stats.Point((-0.5*direction, 0.0, 0.0)), + angle=openmc.stats.Monodirectional((direction, 0.0, 0.0)), + ) + settings.weight_windows = ww + settings.weight_window_checkpoints = {'collision': False, 'surface': True} + + tally = openmc.Tally() + tally.filters = [ + openmc.CellFilter([target_cell]), + openmc.WeightFilter([0.0, 0.5, 2.0]) + ] + tally.scores = ['flux'] + + model = openmc.Model(geometry, materials, settings, openmc.Tallies([tally])) + model.run(apply_tally_results=True) + + # All flux in the target cell should be carried by split particles with + # weight 0.2; total flux is conserved (1 cm path length per source particle) + flux = tally.mean.squeeze() + assert flux == pytest.approx([1.0, 0.0])