From 6910e3657257ca1dcd6a00370614655afe6b3201 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sun, 26 Mar 2023 23:22:25 -0500 Subject: [PATCH 1/9] Adjustment to coincidence checking and case for r_grid[0] == 0 --- src/mesh.cpp | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 3ee489a68..5eb400a4b 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1066,6 +1066,9 @@ double CylindricalMesh::find_r_crossing( // s^2 * (u^2 + v^2) + 2*s*(u*x+v*y) + x^2+y^2-r0^2 = 0 const double r0 = grid_[0][shell]; + if (r0 == 0.0) + return INFTY; + const double denominator = u.x * u.x + u.y * u.y; // Direction of flight is in z-direction. Will never intersect r. @@ -1085,6 +1088,9 @@ double CylindricalMesh::find_r_crossing( D = std::sqrt(D); // the solution -p - D is always smaller as -p + D : Check this one first + if (std::abs(c) <= 1e-10) + return INFTY; + if (-p - D > l && std::abs(c) > 1e-10) return -p - D; if (-p + D > l) @@ -1304,14 +1310,19 @@ double SphericalMesh::find_r_crossing( // solve |r+s*u| = r0 // |r+s*u| = |r| + 2*s*r*u + s^2 (|u|==1 !) const double r0 = grid_[0][shell]; + if (r0 == 0.0) + return INFTY; const double p = r.dot(u); double c = r.dot(r) - r0 * r0; double D = p * p - c; + if (std::abs(c) <= 1e-10) + return INFTY; + if (D >= 0.0) { D = std::sqrt(D); // the solution -p - D is always smaller as -p + D : Check this one first - if (-p - D > l && std::abs(c) > 1e-10) + if (-p - D > l) return -p - D; if (-p + D > l) return -p + D; From 160360a2a6d2cf1f9bdbbdefa921c2bbca8f394e Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 27 Mar 2023 15:09:15 -0500 Subject: [PATCH 2/9] Remove old condition for cylindrical mesh --- src/mesh.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 5eb400a4b..3d94d21a0 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1091,7 +1091,7 @@ double CylindricalMesh::find_r_crossing( if (std::abs(c) <= 1e-10) return INFTY; - if (-p - D > l && std::abs(c) > 1e-10) + if (-p - D > l) return -p - D; if (-p + D > l) return -p + D; From 13a9e8fae4567f2a19156f087f6e4f038f032b20 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 27 Mar 2023 15:17:15 -0500 Subject: [PATCH 3/9] Adding tests with void geometry for cylindrical and spherical mesh --- tests/unit_tests/test_cylindrical_mesh.py | 89 +++++++++++++++++++ tests/unit_tests/test_spherical_mesh.py | 100 ++++++++++++++++++++++ 2 files changed, 189 insertions(+) diff --git a/tests/unit_tests/test_cylindrical_mesh.py b/tests/unit_tests/test_cylindrical_mesh.py index 83ba09cc5..c8346bb4a 100644 --- a/tests/unit_tests/test_cylindrical_mesh.py +++ b/tests/unit_tests/test_cylindrical_mesh.py @@ -103,3 +103,92 @@ def test_offset_mesh(model, estimator, origin): mean[i, j, k] == 0.0 else: mean[i, j, k] != 0.0 + +@pytest.fixture() +def void_coincident_geom_model(): + """A model with many geometric boundaries coincident with mesh boundaries + across many scales + """ + openmc.reset_auto_ids() + model = openmc.model.Model() + + model.materials = openmc.Materials() + radii = [0.1,1, 5, 50, 100, 150, 250] + spheres = [openmc.Sphere(r=ri) for ri in radii] + spheres[-1].boundary_type = 'vacuum' + + regions = openmc.model.subdivide(spheres)[:-1] + cells = [openmc.Cell(region=r, fill=None) for r in regions] + geom = openmc.Geometry(cells) + + model.geometry = geom + + settings = openmc.Settings(run_mode='fixed source') + settings.batches = 2 + settings.particles = 1000 + model.settings = settings + + mesh = openmc.CylindricalMesh() + mesh.r_grid = np.linspace(0, 250, 501) + mesh.z_grid = [-250, 250] + mesh.phi_grid = np.linspace(0, 2*np.pi, 2) + mesh_filter = openmc.MeshFilter(mesh) + + tally = openmc.Tally() + tally.scores = ['flux'] + tally.filters = [mesh_filter] + + model.tallies = openmc.Tallies([tally]) + + return model + + +# convenience function for checking tally results +# in the following tests +def _check_void_cylindrical_tally(statepoint_filename): + with openmc.StatePoint(statepoint_filename) as sp: + flux_tally = sp.tallies[1] + mesh = flux_tally.find_filter(openmc.MeshFilter).mesh + neutron_flux = flux_tally.get_reshaped_data().squeeze() / mesh.volumes.flatten() + flux_diff = np.diff(neutron_flux) + # ensure flux values are monotonically decreasing due to + # geometric attenuation + assert (flux_diff < 0.0).all() + + +def test_void_geom_pnt_src(run_in_tmpdir, void_coincident_geom_model): + src = openmc.Source() + src.space = openmc.stats.Point() + src.angle = openmc.stats.PolarAzimuthal(mu=openmc.stats.Discrete([0.0], [1.0])) + src.energy = openmc.stats.Discrete([14.06e6], [1]) + void_coincident_geom_model.settings.source = src + + sp_filename = void_coincident_geom_model.run() + _check_void_cylindrical_tally(sp_filename) + + +def test_void_geom_boundary_src(run_in_tmpdir, void_coincident_geom_model): + bbox = void_coincident_geom_model.geometry.bounding_box + + outer_r = bbox[1][0] - 0.0001 + + radial_vals = np.linspace(0.0, 2.0*np.pi, 100) + + sources = [] + + energy = openmc.stats.Discrete([14.06e6], [1]) + for val in radial_vals: + src = openmc.Source() + src.energy = energy + + pnt = np.array([np.cos(val), np.sin(val), 0.0]) + u = -pnt + src.space = openmc.stats.Point(outer_r*pnt) + src.angle = openmc.stats.Monodirectional(u) + + sources.append(src) + + void_coincident_geom_model.settings.source = sources + sp_filename = void_coincident_geom_model.run() + + _check_void_cylindrical_tally(sp_filename) \ No newline at end of file diff --git a/tests/unit_tests/test_spherical_mesh.py b/tests/unit_tests/test_spherical_mesh.py index 73153f5cc..860156671 100644 --- a/tests/unit_tests/test_spherical_mesh.py +++ b/tests/unit_tests/test_spherical_mesh.py @@ -107,3 +107,103 @@ def test_offset_mesh(model, estimator, origin): mean[i, j, k] == 0.0 else: mean[i, j, k] != 0.0 + +# Some void geometry tests to check our radial intersection methods on +# spherical and cylindrical meshes + +@pytest.fixture() +def void_coincident_geom_model(): + """A model with many geometric boundaries coincident with mesh boundaries + across many scales + """ + openmc.reset_auto_ids() + + model = openmc.Model() + + model.materials = openmc.Materials() + radii = [0.1, 1, 5, 50, 100, 150, 250] + spheres = [openmc.Sphere(r=ri) for ri in radii] + spheres[-1].boundary_type = 'vacuum' + + regions = openmc.model.subdivide(spheres)[:-1] + cells = [openmc.Cell(region=r, fill=None) for r in regions] + geom = openmc.Geometry(cells) + + model.geometry = geom + + settings = openmc.Settings(run_mode='fixed source') + settings.batches = 2 + settings.particles = 5000 + model.settings = settings + + mesh = openmc.SphericalMesh() + mesh.r_grid = np.linspace(0, 250, 501) + mesh_filter = openmc.MeshFilter(mesh) + + tally = openmc.Tally() + tally.scores = ['flux'] + tally.filters = [mesh_filter] + + model.tallies = openmc.Tallies([tally]) + + return model + + +# convenience function for checking tally results +# in the following tests +def _check_void_spherical_tally(statepoint_filename): + with openmc.StatePoint(statepoint_filename) as sp: + flux_tally = sp.tallies[1] + mesh = flux_tally.find_filter(openmc.MeshFilter).mesh + neutron_flux = flux_tally.get_reshaped_data().squeeze() / mesh.volumes.flatten() + flux_diff = np.diff(neutron_flux) + # ensure flux values are monotonically decreasing due to + # geometric attenuation + assert (flux_diff < 0.0).all() + + +def test_void_geom_pnt_src(run_in_tmpdir, void_coincident_geom_model): + # add isotropic point source + src = openmc.Source() + src.space = openmc.stats.Point() + src.energy = openmc.stats.Discrete([14.06e6], [1]) + void_coincident_geom_model.settings.source = src + + sp_filename = void_coincident_geom_model.run() + + _check_void_spherical_tally(sp_filename) + + +def test_void_geom_boundary_src(run_in_tmpdir, void_coincident_geom_model): + # update source to a number of points on the outside of the sphere + # with directions pointing toward the origin + + # sources + phi_vals = np.linspace(0, np.pi, 20) + theta_vals = np.linspace(0, 2.0*np.pi, 20) + + bbox = void_coincident_geom_model.geometry.bounding_box + # can't source particles directly on the geometry boundary + outer_r = bbox[1][0] - 0.00001 + + sources = [] + + energy = openmc.stats.Discrete([14.06e6], [1]) + + for phi, theta in zip(phi_vals, theta_vals): + + src = openmc.Source() + src.energy = energy + + pnt = np.array([np.sin(phi)*np.cos(theta), np.sin(phi)*np.sin(theta), np.cos(phi)]) + u = -pnt + src.space = openmc.stats.Point(outer_r*pnt) + src.angle = openmc.stats.Monodirectional(u) + + sources.append(src) + + void_coincident_geom_model.settings.source = sources + + sp_filename = void_coincident_geom_model.run() + + _check_void_spherical_tally(sp_filename) From ab403b2fc8cf0aed15065c7caeafdfa79f4c07af Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 27 Mar 2023 16:15:38 -0500 Subject: [PATCH 4/9] Improvements to tests --- tests/unit_tests/test_cylindrical_mesh.py | 22 ++++++++++++-------- tests/unit_tests/test_spherical_mesh.py | 25 ++++++++++++----------- 2 files changed, 26 insertions(+), 21 deletions(-) diff --git a/tests/unit_tests/test_cylindrical_mesh.py b/tests/unit_tests/test_cylindrical_mesh.py index c8346bb4a..27c4d5a1f 100644 --- a/tests/unit_tests/test_cylindrical_mesh.py +++ b/tests/unit_tests/test_cylindrical_mesh.py @@ -149,12 +149,12 @@ def _check_void_cylindrical_tally(statepoint_filename): with openmc.StatePoint(statepoint_filename) as sp: flux_tally = sp.tallies[1] mesh = flux_tally.find_filter(openmc.MeshFilter).mesh - neutron_flux = flux_tally.get_reshaped_data().squeeze() / mesh.volumes.flatten() - flux_diff = np.diff(neutron_flux) - # ensure flux values are monotonically decreasing due to - # geometric attenuation - assert (flux_diff < 0.0).all() - + neutron_flux = flux_tally.get_reshaped_data().squeeze() + # we expect the tally results to be the same as the mesh grid width + # for these cases + print(neutron_flux) + d_r = mesh.r_grid[1] - mesh.r_grid[0] + assert neutron_flux == pytest.approx(d_r) def test_void_geom_pnt_src(run_in_tmpdir, void_coincident_geom_model): src = openmc.Source() @@ -168,11 +168,15 @@ def test_void_geom_pnt_src(run_in_tmpdir, void_coincident_geom_model): def test_void_geom_boundary_src(run_in_tmpdir, void_coincident_geom_model): + # update source to a number of points on the outside of the cylinder + # with directions pointing toward the origin bbox = void_coincident_geom_model.geometry.bounding_box - outer_r = bbox[1][0] - 0.0001 + # can't source particle directly on the geometry boundary + outer_r = bbox[1][0] - 1e-08 - radial_vals = np.linspace(0.0, 2.0*np.pi, 100) + n_sources = 100 + radial_vals = np.linspace(0.0, 2.0*np.pi, n_sources) sources = [] @@ -185,7 +189,7 @@ def test_void_geom_boundary_src(run_in_tmpdir, void_coincident_geom_model): u = -pnt src.space = openmc.stats.Point(outer_r*pnt) src.angle = openmc.stats.Monodirectional(u) - + src.strength = 0.5/n_sources sources.append(src) void_coincident_geom_model.settings.source = sources diff --git a/tests/unit_tests/test_spherical_mesh.py b/tests/unit_tests/test_spherical_mesh.py index 860156671..0f02e9699 100644 --- a/tests/unit_tests/test_spherical_mesh.py +++ b/tests/unit_tests/test_spherical_mesh.py @@ -155,11 +155,11 @@ def _check_void_spherical_tally(statepoint_filename): with openmc.StatePoint(statepoint_filename) as sp: flux_tally = sp.tallies[1] mesh = flux_tally.find_filter(openmc.MeshFilter).mesh - neutron_flux = flux_tally.get_reshaped_data().squeeze() / mesh.volumes.flatten() - flux_diff = np.diff(neutron_flux) - # ensure flux values are monotonically decreasing due to - # geometric attenuation - assert (flux_diff < 0.0).all() + neutron_flux = flux_tally.get_reshaped_data().squeeze() + # the flux values for each bin should equal the width + # width of the mesh bins + d_r = mesh.r_grid[1] - mesh.r_grid[0] + assert neutron_flux == pytest.approx(d_r) def test_void_geom_pnt_src(run_in_tmpdir, void_coincident_geom_model): @@ -169,22 +169,21 @@ def test_void_geom_pnt_src(run_in_tmpdir, void_coincident_geom_model): src.energy = openmc.stats.Discrete([14.06e6], [1]) void_coincident_geom_model.settings.source = src + # run model and check tally results sp_filename = void_coincident_geom_model.run() - _check_void_spherical_tally(sp_filename) def test_void_geom_boundary_src(run_in_tmpdir, void_coincident_geom_model): # update source to a number of points on the outside of the sphere # with directions pointing toward the origin - - # sources - phi_vals = np.linspace(0, np.pi, 20) - theta_vals = np.linspace(0, 2.0*np.pi, 20) + n_sources = 20 + phi_vals = np.linspace(0, np.pi, n_sources) + theta_vals = np.linspace(0, 2.0*np.pi, n_sources) bbox = void_coincident_geom_model.geometry.bounding_box # can't source particles directly on the geometry boundary - outer_r = bbox[1][0] - 0.00001 + outer_r = bbox[1][0] - 1e-08 sources = [] @@ -199,11 +198,13 @@ def test_void_geom_boundary_src(run_in_tmpdir, void_coincident_geom_model): u = -pnt src.space = openmc.stats.Point(outer_r*pnt) src.angle = openmc.stats.Monodirectional(u) + # set source strengths so that we can still expect + # a tally value of 0.5 + src.strength = 0.5/n_sources sources.append(src) void_coincident_geom_model.settings.source = sources sp_filename = void_coincident_geom_model.run() - _check_void_spherical_tally(sp_filename) From 63fa42843b04966d75bb4eb0c3a91e12e531d2b9 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 27 Mar 2023 16:50:18 -0500 Subject: [PATCH 5/9] Removing unused print statements --- tests/unit_tests/test_cylindrical_mesh.py | 2 -- tests/unit_tests/test_spherical_mesh.py | 1 - 2 files changed, 3 deletions(-) diff --git a/tests/unit_tests/test_cylindrical_mesh.py b/tests/unit_tests/test_cylindrical_mesh.py index 27c4d5a1f..cbdb1d315 100644 --- a/tests/unit_tests/test_cylindrical_mesh.py +++ b/tests/unit_tests/test_cylindrical_mesh.py @@ -98,7 +98,6 @@ def test_offset_mesh(model, estimator, origin): centroids = mesh.centroids for ijk in mesh.indices: i, j, k = np.array(ijk) - 1 - print(centroids[:, i, j, k]) if model.geometry.find(centroids[:, i, j, k]): mean[i, j, k] == 0.0 else: @@ -152,7 +151,6 @@ def _check_void_cylindrical_tally(statepoint_filename): neutron_flux = flux_tally.get_reshaped_data().squeeze() # we expect the tally results to be the same as the mesh grid width # for these cases - print(neutron_flux) d_r = mesh.r_grid[1] - mesh.r_grid[0] assert neutron_flux == pytest.approx(d_r) diff --git a/tests/unit_tests/test_spherical_mesh.py b/tests/unit_tests/test_spherical_mesh.py index 0f02e9699..ba3167d6f 100644 --- a/tests/unit_tests/test_spherical_mesh.py +++ b/tests/unit_tests/test_spherical_mesh.py @@ -102,7 +102,6 @@ def test_offset_mesh(model, estimator, origin): centroids = mesh.centroids for ijk in mesh.indices: i, j, k = np.array(ijk) - 1 - print(centroids[:, i, j, k]) if model.geometry.find(centroids[:, i, j, k]): mean[i, j, k] == 0.0 else: From 26c031a041e110b1ac944044c5d011e4b874b20d Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 27 Mar 2023 16:54:42 -0500 Subject: [PATCH 6/9] Run tests in tmpdir to keep repo clean --- tests/unit_tests/test_cylindrical_mesh.py | 2 +- tests/unit_tests/test_spherical_mesh.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_cylindrical_mesh.py b/tests/unit_tests/test_cylindrical_mesh.py index cbdb1d315..2be0f4962 100644 --- a/tests/unit_tests/test_cylindrical_mesh.py +++ b/tests/unit_tests/test_cylindrical_mesh.py @@ -75,7 +75,7 @@ def label(p): return f'estimator:{p}' @pytest.mark.parametrize('estimator,origin', test_cases, ids=label) -def test_offset_mesh(model, estimator, origin): +def test_offset_mesh(run_in_tmpdir, model, estimator, origin): """Tests that the mesh has been moved based on tally results """ mesh = model.tallies[0].filters[0].mesh diff --git a/tests/unit_tests/test_spherical_mesh.py b/tests/unit_tests/test_spherical_mesh.py index ba3167d6f..ac2781378 100644 --- a/tests/unit_tests/test_spherical_mesh.py +++ b/tests/unit_tests/test_spherical_mesh.py @@ -79,7 +79,7 @@ def label(p): return f'estimator:{p}' @pytest.mark.parametrize('estimator,origin', test_cases, ids=label) -def test_offset_mesh(model, estimator, origin): +def test_offset_mesh(run_in_tmpdir, model, estimator, origin): """Tests that the mesh has been moved based on tally results """ mesh = model.tallies[0].filters[0].mesh From e5863eec1b80a24ffb471eb8287125c1e18d00db Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 28 Mar 2023 07:35:22 -0500 Subject: [PATCH 7/9] Correcting variable name --- tests/unit_tests/test_cylindrical_mesh.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/unit_tests/test_cylindrical_mesh.py b/tests/unit_tests/test_cylindrical_mesh.py index 2be0f4962..caf933360 100644 --- a/tests/unit_tests/test_cylindrical_mesh.py +++ b/tests/unit_tests/test_cylindrical_mesh.py @@ -113,10 +113,10 @@ def void_coincident_geom_model(): model.materials = openmc.Materials() radii = [0.1,1, 5, 50, 100, 150, 250] - spheres = [openmc.Sphere(r=ri) for ri in radii] - spheres[-1].boundary_type = 'vacuum' + cylinders = [openmc.Sphere(r=ri) for ri in radii] + cylinders[-1].boundary_type = 'vacuum' - regions = openmc.model.subdivide(spheres)[:-1] + regions = openmc.model.subdivide(cylinders)[:-1] cells = [openmc.Cell(region=r, fill=None) for r in regions] geom = openmc.Geometry(cells) From 86e9a4bd4b1c619e3f39b761a3a81658665ded19 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 28 Mar 2023 08:47:42 -0500 Subject: [PATCH 8/9] Adding constants for coincidence checking --- include/openmc/constants.h | 4 ++++ src/mesh.cpp | 4 ++-- src/surface.cpp | 2 +- 3 files changed, 7 insertions(+), 3 deletions(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 2a867c0ca..b0031fae3 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -51,6 +51,10 @@ constexpr double FP_PRECISION {1e-14}; constexpr double FP_REL_PRECISION {1e-5}; constexpr double FP_COINCIDENT {1e-12}; +// Coincidence tolerances +constexpr double TORUS_TOL {1e-10}; +constexpr double RADIAL_MESH_TOL {1e-10}; + // Maximum number of random samples per history constexpr int MAX_SAMPLE {100000}; diff --git a/src/mesh.cpp b/src/mesh.cpp index 3d94d21a0..e72a5cf7d 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1088,7 +1088,7 @@ double CylindricalMesh::find_r_crossing( D = std::sqrt(D); // the solution -p - D is always smaller as -p + D : Check this one first - if (std::abs(c) <= 1e-10) + if (std::abs(c) <= RADIAL_MESH_TOL) return INFTY; if (-p - D > l) @@ -1316,7 +1316,7 @@ double SphericalMesh::find_r_crossing( double c = r.dot(r) - r0 * r0; double D = p * p - c; - if (std::abs(c) <= 1e-10) + if (std::abs(c) <= RADIAL_MESH_TOL) return INFTY; if (D >= 0.0) { diff --git a/src/surface.cpp b/src/surface.cpp index c0b3f0a59..2403d0c6a 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1044,7 +1044,7 @@ double torus_distance(double x1, double x2, double x3, double u1, double u2, // zero but possibly small and positive. A tolerance is set to discard that // zero. double distance = INFTY; - double cutoff = coincident ? 1e-10 : 0.0; + double cutoff = coincident ? TORUS_TOL : 0.0; for (int i = 0; i < 4; ++i) { if (roots[i].imag() == 0) { double root = roots[i].real(); From 6cb4401230e9189169a4754180ec3d6a2b3e20b1 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 28 Mar 2023 10:11:14 -0500 Subject: [PATCH 9/9] Correcting surface type in cylinder tests --- tests/unit_tests/test_cylindrical_mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit_tests/test_cylindrical_mesh.py b/tests/unit_tests/test_cylindrical_mesh.py index caf933360..4382b49ba 100644 --- a/tests/unit_tests/test_cylindrical_mesh.py +++ b/tests/unit_tests/test_cylindrical_mesh.py @@ -113,7 +113,7 @@ def void_coincident_geom_model(): model.materials = openmc.Materials() radii = [0.1,1, 5, 50, 100, 150, 250] - cylinders = [openmc.Sphere(r=ri) for ri in radii] + cylinders = [openmc.ZCylinder(r=ri) for ri in radii] cylinders[-1].boundary_type = 'vacuum' regions = openmc.model.subdivide(cylinders)[:-1]