From b25cf8635db75d73d56990b570e6bbeb132ed4a3 Mon Sep 17 00:00:00 2001 From: GuySten Date: Thu, 1 Jan 2026 21:38:10 +0200 Subject: [PATCH 1/4] simplify and fix RectLattice distance calculation lower mesh resolution --- src/lattice.cpp | 47 ++-- .../lattice_corner_crossing/__init__.py | 0 .../lattice_corner_crossing/inputs_true.dat | 61 ++++++ .../lattice_corner_crossing/results_true.dat | 201 ++++++++++++++++++ .../lattice_corner_crossing/test.py | 111 ++++++++++ 5 files changed, 385 insertions(+), 35 deletions(-) create mode 100644 tests/regression_tests/lattice_corner_crossing/__init__.py create mode 100644 tests/regression_tests/lattice_corner_crossing/inputs_true.dat create mode 100644 tests/regression_tests/lattice_corner_crossing/results_true.dat create mode 100644 tests/regression_tests/lattice_corner_crossing/test.py diff --git a/src/lattice.cpp b/src/lattice.cpp index efbfcb2163a..60bceb01862 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -260,46 +260,23 @@ std::pair> RectLattice::distance( // Determine the oncoming edge. double x0 {copysign(0.5 * pitch_[0], u.x)}; double y0 {copysign(0.5 * pitch_[1], u.y)}; + double z0; - // Left and right sides - double d {INFTY}; - array lattice_trans; - if ((std::abs(x - x0) > FP_PRECISION) && u.x != 0) { - d = (x0 - x) / u.x; - if (u.x > 0) { - lattice_trans = {1, 0, 0}; - } else { - lattice_trans = {-1, 0, 0}; - } + double d = std::min((x0 - x) / u.x, (y0 - y) / u.y); + if (is_3d_) { + z0 = copysign(0.5 * pitch_[2], u.z); + d = std::min(d, (z0 - z) / u.z); } - // Front and back sides - if ((std::abs(y - y0) > FP_PRECISION) && u.y != 0) { - double this_d = (y0 - y) / u.y; - if (this_d < d) { - d = this_d; - if (u.y > 0) { - lattice_trans = {0, 1, 0}; - } else { - lattice_trans = {0, -1, 0}; - } - } - } + array lattice_trans = {0, 0, 0}; - // Top and bottom sides + if (std::abs(x + u.x * d - x0) < FP_PRECISION) + lattice_trans[0] = copysign(1, u.x); + if (std::abs(y + u.y * d - y0) < FP_PRECISION) + lattice_trans[1] = copysign(1, u.y); if (is_3d_) { - double z0 {copysign(0.5 * pitch_[2], u.z)}; - if ((std::abs(z - z0) > FP_PRECISION) && u.z != 0) { - double this_d = (z0 - z) / u.z; - if (this_d < d) { - d = this_d; - if (u.z > 0) { - lattice_trans = {0, 0, 1}; - } else { - lattice_trans = {0, 0, -1}; - } - } - } + if (std::abs(z + u.z * d - z0) < FP_PRECISION) + lattice_trans[2] = copysign(1, u.z); } return {d, lattice_trans}; diff --git a/tests/regression_tests/lattice_corner_crossing/__init__.py b/tests/regression_tests/lattice_corner_crossing/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/lattice_corner_crossing/inputs_true.dat b/tests/regression_tests/lattice_corner_crossing/inputs_true.dat new file mode 100644 index 00000000000..1156d77f0e1 --- /dev/null +++ b/tests/regression_tests/lattice_corner_crossing/inputs_true.dat @@ -0,0 +1,61 @@ + + + + + + + + + + + + + + + + + + + 10.0 10.0 + 2 2 + -10.0 -10.0 + +2 3 +3 2 + + + + + + + + + + + fixed source + 1000 + 10 + 100000000 + + + -0.7071067811865476 -0.7071067811865475 0.0 + + + + + + + 10 10 1 + -20.0 -20.0 -1000000000.0 + 20.0 20.0 1000000000.0 + + + 1 + + + 1 + flux + tracklength + + + diff --git a/tests/regression_tests/lattice_corner_crossing/results_true.dat b/tests/regression_tests/lattice_corner_crossing/results_true.dat new file mode 100644 index 00000000000..bd81cf53534 --- /dev/null +++ b/tests/regression_tests/lattice_corner_crossing/results_true.dat @@ -0,0 +1,201 @@ +tally 1: +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.648775E-03 +7.016009E-06 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +2.981847E-03 +8.891414E-06 +3.331677E-03 +1.110007E-05 +6.742237E-03 +4.545776E-05 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +4.388206E-03 +1.925636E-05 +2.032104E-03 +4.129446E-06 +3.014946E-03 +9.089898E-06 +7.058968E-03 +4.982903E-05 +9.013188E-04 +8.123755E-07 +7.501461E-04 +5.627191E-07 +0.000000E+00 +0.000000E+00 +5.134898E-03 +2.636718E-05 +1.207302E-03 +1.457579E-06 +3.045925E-03 +9.277657E-06 +4.132399E-03 +1.707672E-05 +1.055271E-02 +5.829711E-05 +4.014909E-03 +1.611949E-05 +2.698215E-03 +7.280363E-06 +1.751215E-02 +1.073504E-04 +1.356908E-02 +9.589157E-05 +5.451466E-03 +2.971848E-05 +3.165676E-04 +1.002151E-07 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +3.724938E-03 +1.387516E-05 +2.812863E-03 +7.278837E-06 +1.000700E+01 +1.001402E+01 +2.345904E-02 +2.127641E-04 +2.451653E-02 +1.750725E-04 +8.695904E-03 +5.553981E-05 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +4.947565E-03 +2.447840E-05 +1.725721E-02 +8.878803E-05 +5.656444E+01 +3.199538E+02 +2.159961E-02 +1.034364E-04 +3.091063E-02 +5.543317E-04 +4.232209E-03 +1.791159E-05 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +4.716170E-03 +2.224226E-05 +1.321855E-02 +5.288265E-05 +5.309576E-02 +7.984865E-04 +5.656179E+01 +3.199238E+02 +2.399311E-02 +2.936670E-04 +4.680753E-02 +9.526509E-04 +1.074699E-02 +1.154978E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +4.526425E-03 +2.048853E-05 +3.983003E-03 +9.553596E-06 +1.726789E-02 +9.579848E-05 +3.355605E-02 +3.331438E-04 +3.578422E-02 +3.560242E-04 +5.649209E+01 +3.191358E+02 +2.591495E-02 +4.169753E-04 +9.484744E-03 +8.065715E-05 +0.000000E+00 +0.000000E+00 +4.232587E-03 +9.113263E-06 +5.787378E-03 +2.194940E-05 +6.193626E-03 +2.357423E-05 +4.552787E-03 +6.987360E-06 +9.641308E-03 +3.506339E-05 +1.303892E-02 +4.378408E-05 +1.831744E-02 +8.836683E-05 +3.024725E+01 +9.148969E+01 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +1.488204E-03 +2.214752E-06 +5.067487E-03 +1.520987E-05 +7.550927E-03 +2.993991E-05 +3.215559E-03 +1.033982E-05 +4.258525E-03 +1.205064E-05 +2.115487E-03 +4.475286E-06 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 diff --git a/tests/regression_tests/lattice_corner_crossing/test.py b/tests/regression_tests/lattice_corner_crossing/test.py new file mode 100644 index 00000000000..74b063dfbb9 --- /dev/null +++ b/tests/regression_tests/lattice_corner_crossing/test.py @@ -0,0 +1,111 @@ +import openmc +import pytest + +import numpy as np +from tests.testing_harness import PyAPITestHarness + + +@pytest.fixture +def model(): + model = openmc.model.Model() + + + # length of lattice on each side + lat_size = 20.0 + + # angle that we're crossing the corner at + phi = np.pi/4.0 + + # number of lattice elements across + numlat = 2 + + air = openmc.Material(name='air') + air.set_density('g/cc', 0.001) + air.add_nuclide('N14', 1.) + + metal = openmc.Material(name='metal') + metal.set_density('g/cc', 7) + metal.add_nuclide('Fe56', 1.) + + model.materials.extend([air, metal]) + + + left = openmc.XPlane(x0=-lat_size/2.0, name='left') + right = openmc.XPlane(x0=lat_size/2.0, name='right') + bottom = openmc.YPlane(y0=-lat_size/2.0, name='bottom') + top = openmc.YPlane(y0=lat_size/2.0, name='top') + cyl = openmc.ZCylinder(x0=0, y0=0, r=lat_size) + cyl.boundary_type = 'vacuum' + + in_lattice_region = +left & -right & +bottom & -top + + outside_lattice = openmc.Cell(name='outside_lattice') + outside_lattice.region = -cyl & ~in_lattice_region + outside_lattice.fill = air + + inside_lattice = openmc.Cell(name='inside_lattice') + inside_lattice.region = in_lattice_region + + root = openmc.Universe(name='root universe') + + # can we straightforwardly define infinite universes? + metal_uni = openmc.Universe() + metal_cell = openmc.Cell() + metal_cell.fill = metal + metal_cell.region = -openmc.Sphere(r=1e90) + metal_uni.add_cell(metal_cell) + + air_uni = openmc.Universe() + air_cell = openmc.Cell() + air_cell.fill = air + air_cell.region = -openmc.Sphere(r=1e90) + air_uni.add_cell(air_cell) + + # define a checkerboard lattice + lattice = openmc.RectLattice() + lattice.lower_left = [-lat_size/2.0, -lat_size/2.0] + lattice.pitch = [lat_size/numlat, lat_size/numlat] + lattice.universes = [[metal_uni if (i%2)==(j%2) else air_uni for i in range(numlat)] + for j in range(numlat)] + inside_lattice.fill = lattice + root.add_cells([outside_lattice, inside_lattice]) + + model.geometry = openmc.Geometry(root) + + + # Instantiate a Settings object, set all runtime parameters + model.settings = openmc.Settings() + model.settings.run_mode = 'fixed source' + model.settings.batches = 10 + model.settings.particles = int(1e3) + model.settings.max_lost_particles = int(1e8) + + # define a source located outside the lattice and pointing straight into its + # corner at 45 degrees + angular_distr = openmc.stats.Monodirectional(reference_uvw=[np.cos(phi), np.sin(phi), 0.0]) + offset = 0.0 #1e-6 + spatial_distr = openmc.stats.Point(xyz=[-np.cos(phi)+offset, -np.sin(phi), 0.0]) + model.settings.source = openmc.source.Source(space=spatial_distr, angle=angular_distr) + + # Instantiate a tally mesh + mesh = openmc.RegularMesh() + mesh.dimension = [10, 10, 1] + mesh.lower_left = [-lat_size, -lat_size, -1e9] + mesh.upper_right = [lat_size, lat_size, 1e9] + mesh_filter = openmc.MeshFilter(mesh) + + # Instantiate the Tally + tally = openmc.Tally(tally_id=1) + tally.filters = [mesh_filter] + tally.scores = ['flux'] + tally.estimator = 'tracklength' + model.tallies = openmc.Tallies([tally]) + return model + + + +def test_lattice_corner_crossing(model): + # Ensure we account for potential corner crossings + # in floating point precision. + harness = PyAPITestHarness('statepoint.10.h5', model) + harness.main() From ed3217db94779b624355bbd3e9d78f488665a2fd Mon Sep 17 00:00:00 2001 From: GuySten Date: Tue, 6 Jan 2026 00:29:00 +0200 Subject: [PATCH 2/4] add guards against divide by zero --- src/lattice.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/lattice.cpp b/src/lattice.cpp index 60bceb01862..a0f2ea5f402 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -262,10 +262,11 @@ std::pair> RectLattice::distance( double y0 {copysign(0.5 * pitch_[1], u.y)}; double z0; - double d = std::min((x0 - x) / u.x, (y0 - y) / u.y); + double d = std::min( + u.x != 0.0 ? (x0 - x) / u.x : INFTY, u.y != 0.0 ? (y0 - y) / u.y : INFTY); if (is_3d_) { z0 = copysign(0.5 * pitch_[2], u.z); - d = std::min(d, (z0 - z) / u.z); + d = std::min(d, u.z != 0.0 ? (z0 - z) / u.z : INFTY); } array lattice_trans = {0, 0, 0}; From 568b4e15d53c8a6279e7cf59238dba5c6213f0be Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 5 Jan 2026 21:12:21 -0600 Subject: [PATCH 3/4] Guard against zero in checks for lattice_trans --- src/lattice.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/lattice.cpp b/src/lattice.cpp index a0f2ea5f402..92d451f61fd 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -269,14 +269,14 @@ std::pair> RectLattice::distance( d = std::min(d, u.z != 0.0 ? (z0 - z) / u.z : INFTY); } + // Determine which lattice boundaries are being crossed array lattice_trans = {0, 0, 0}; - - if (std::abs(x + u.x * d - x0) < FP_PRECISION) + if (u.x != 0.0 && std::abs(x + u.x * d - x0) < FP_PRECISION) lattice_trans[0] = copysign(1, u.x); - if (std::abs(y + u.y * d - y0) < FP_PRECISION) + if (u.y != 0.0 && std::abs(y + u.y * d - y0) < FP_PRECISION) lattice_trans[1] = copysign(1, u.y); if (is_3d_) { - if (std::abs(z + u.z * d - z0) < FP_PRECISION) + if (u.z != 0.0 && std::abs(z + u.z * d - z0) < FP_PRECISION) lattice_trans[2] = copysign(1, u.z); } From b82f8d85a02f47c3c7f18efa16a80fc462657509 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 5 Jan 2026 21:33:49 -0600 Subject: [PATCH 4/4] Update lattice corner test --- .../lattice_corner_crossing/inputs_true.dat | 41 +++--- .../lattice_corner_crossing/test.py | 128 +++++++----------- 2 files changed, 69 insertions(+), 100 deletions(-) diff --git a/tests/regression_tests/lattice_corner_crossing/inputs_true.dat b/tests/regression_tests/lattice_corner_crossing/inputs_true.dat index 1156d77f0e1..4b49b14031a 100644 --- a/tests/regression_tests/lattice_corner_crossing/inputs_true.dat +++ b/tests/regression_tests/lattice_corner_crossing/inputs_true.dat @@ -1,41 +1,38 @@ - - + + - - + + - - - - - + + + + + 10.0 10.0 2 2 -10.0 -10.0 -2 3 -3 2 +1 2 +2 1 - - - - - - - + + + + + fixed source 1000 10 - 100000000 -0.7071067811865476 -0.7071067811865475 0.0 @@ -45,9 +42,9 @@ - 10 10 1 - -20.0 -20.0 -1000000000.0 - 20.0 20.0 1000000000.0 + 10 10 + -20.0 -20.0 + 20.0 20.0 1 diff --git a/tests/regression_tests/lattice_corner_crossing/test.py b/tests/regression_tests/lattice_corner_crossing/test.py index 74b063dfbb9..4684d144e1e 100644 --- a/tests/regression_tests/lattice_corner_crossing/test.py +++ b/tests/regression_tests/lattice_corner_crossing/test.py @@ -1,111 +1,83 @@ +""" +This test is designed to ensure that we account for potential corner crossings +in floating point precision. + +""" + +from math import pi, cos, sin + import openmc import pytest -import numpy as np from tests.testing_harness import PyAPITestHarness @pytest.fixture def model(): - model = openmc.model.Model() - - - # length of lattice on each side - lat_size = 20.0 - - # angle that we're crossing the corner at - phi = np.pi/4.0 - - # number of lattice elements across - numlat = 2 - - air = openmc.Material(name='air') - air.set_density('g/cc', 0.001) - air.add_nuclide('N14', 1.) - - metal = openmc.Material(name='metal') - metal.set_density('g/cc', 7) - metal.add_nuclide('Fe56', 1.) - - model.materials.extend([air, metal]) + model = openmc.Model() + # Length of lattice on each side + lat_size = 20.0 - left = openmc.XPlane(x0=-lat_size/2.0, name='left') - right = openmc.XPlane(x0=lat_size/2.0, name='right') - bottom = openmc.YPlane(y0=-lat_size/2.0, name='bottom') - top = openmc.YPlane(y0=lat_size/2.0, name='top') - cyl = openmc.ZCylinder(x0=0, y0=0, r=lat_size) - cyl.boundary_type = 'vacuum' + # Angle that we're crossing the corner at + phi = pi/4.0 - in_lattice_region = +left & -right & +bottom & -top + air = openmc.Material() + air.set_density('g/cm3', 0.001) + air.add_nuclide('N14', 1.0) + metal = openmc.Material() + metal.set_density('g/cm3', 7.0) + metal.add_nuclide('Fe56', 1.0) - outside_lattice = openmc.Cell(name='outside_lattice') - outside_lattice.region = -cyl & ~in_lattice_region - outside_lattice.fill = air + metal_cell = openmc.Cell(fill=metal) + metal_uni = openmc.Universe(cells=[metal_cell]) - inside_lattice = openmc.Cell(name='inside_lattice') - inside_lattice.region = in_lattice_region + air_cell = openmc.Cell(fill=air) + air_uni = openmc.Universe(cells=[air_cell]) - root = openmc.Universe(name='root universe') + # Define a checkerboard lattice + lattice = openmc.RectLattice() + lattice.lower_left = (-lat_size/2.0, -lat_size/2.0) + lattice.pitch = (lat_size/2, lat_size/2) + lattice.universes = [ + [metal_uni, air_uni], + [air_uni, metal_uni] + ] - # can we straightforwardly define infinite universes? - metal_uni = openmc.Universe() - metal_cell = openmc.Cell() - metal_cell.fill = metal - metal_cell.region = -openmc.Sphere(r=1e90) - metal_uni.add_cell(metal_cell) + box = openmc.model.RectangularPrism(lat_size, lat_size) + cyl = openmc.ZCylinder(r=lat_size, boundary_type='vacuum') + outside_lattice = openmc.Cell(region=-cyl & +box, fill=air) + inside_lattice = openmc.Cell(region=-box, fill=lattice) - air_uni = openmc.Universe() - air_cell = openmc.Cell() - air_cell.fill = air - air_cell.region = -openmc.Sphere(r=1e90) - air_uni.add_cell(air_cell) + model.geometry = openmc.Geometry([outside_lattice, inside_lattice]) - # define a checkerboard lattice - lattice = openmc.RectLattice() - lattice.lower_left = [-lat_size/2.0, -lat_size/2.0] - lattice.pitch = [lat_size/numlat, lat_size/numlat] - lattice.universes = [[metal_uni if (i%2)==(j%2) else air_uni for i in range(numlat)] - for j in range(numlat)] - inside_lattice.fill = lattice - root.add_cells([outside_lattice, inside_lattice]) - - model.geometry = openmc.Geometry(root) - - - # Instantiate a Settings object, set all runtime parameters - model.settings = openmc.Settings() + # Set all runtime parameters model.settings.run_mode = 'fixed source' model.settings.batches = 10 - model.settings.particles = int(1e3) - model.settings.max_lost_particles = int(1e8) + model.settings.particles = 1000 - # define a source located outside the lattice and pointing straight into its + # Define a source located outside the lattice and pointing straight into its # corner at 45 degrees - angular_distr = openmc.stats.Monodirectional(reference_uvw=[np.cos(phi), np.sin(phi), 0.0]) - offset = 0.0 #1e-6 - spatial_distr = openmc.stats.Point(xyz=[-np.cos(phi)+offset, -np.sin(phi), 0.0]) - model.settings.source = openmc.source.Source(space=spatial_distr, angle=angular_distr) + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Point((-cos(phi), -sin(phi), 0.0)), + angle=openmc.stats.Monodirectional((cos(phi), sin(phi), 0.0)) + ) - # Instantiate a tally mesh + # Create a mesh tally mesh = openmc.RegularMesh() - mesh.dimension = [10, 10, 1] - mesh.lower_left = [-lat_size, -lat_size, -1e9] - mesh.upper_right = [lat_size, lat_size, 1e9] + mesh.dimension = (10, 10) + mesh.lower_left = (-lat_size, -lat_size) + mesh.upper_right = (lat_size, lat_size) mesh_filter = openmc.MeshFilter(mesh) - - # Instantiate the Tally tally = openmc.Tally(tally_id=1) tally.filters = [mesh_filter] tally.scores = ['flux'] tally.estimator = 'tracklength' - model.tallies = openmc.Tallies([tally]) + model.tallies = [tally] + return model - - + def test_lattice_corner_crossing(model): - # Ensure we account for potential corner crossings - # in floating point precision. harness = PyAPITestHarness('statepoint.10.h5', model) harness.main()