From 10787a2da6c4a154464c61df5e7e4823f686cb2b Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 23 Mar 2023 13:37:26 -0500 Subject: [PATCH 1/3] Account for coincident position in find_r_crossing --- src/mesh.cpp | 10 +-- tests/unit_tests/test_filter_mesh.py | 97 ++++++++++++++++++++++++++++ 2 files changed, 103 insertions(+), 4 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 02e5171bbb..306e560707 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1075,7 +1075,8 @@ double CylindricalMesh::find_r_crossing( const double inv_denominator = 1.0 / denominator; const double p = (u.x * r.x + u.y * r.y) * inv_denominator; - double D = p * p + (r0 * r0 - r.x * r.x - r.y * r.y) * inv_denominator; + double c = r.x * r.x + r.y * r.y - r0 * r0; + double D = p * p - c * inv_denominator; if (D < 0.0) return INFTY; @@ -1083,7 +1084,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 (-p - D > l) + if (-p - D > l && std::abs(c) > FP_COINCIDENT) return -p - D; if (-p + D > l) return -p + D; @@ -1303,12 +1304,13 @@ double SphericalMesh::find_r_crossing( // |r+s*u| = |r| + 2*s*r*u + s^2 (|u|==1 !) const double r0 = grid_[0][shell]; const double p = r.dot(u); - double D = p * p - r.dot(r) + r0 * r0; + double c = r.dot(r) - r0 * r0; + double D = p * p - c; 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) + if (-p - D > l && std::abs(c) > FP_COINCIDENT) return -p - D; if (-p + D > l) return -p + D; diff --git a/tests/unit_tests/test_filter_mesh.py b/tests/unit_tests/test_filter_mesh.py index 3bb836e469..5e5a84f7b0 100644 --- a/tests/unit_tests/test_filter_mesh.py +++ b/tests/unit_tests/test_filter_mesh.py @@ -1,6 +1,7 @@ import math import numpy as np +import pytest from uncertainties import unumpy import openmc @@ -112,6 +113,102 @@ def test_cylindrical_mesh_estimators(run_in_tmpdir): std_dev = unumpy.std_devs(delta) assert np.all(diff < 3*std_dev) + +def test_cylindrical_mesh_coincident(run_in_tmpdir): + """Test for cylindrical mesh boundary being coincident with a cell boundary""" + + fuel = openmc.Material() + fuel.add_nuclide('U235', 1.) + fuel.set_density('g/cm3', 4.5) + + zcyl = openmc.ZCylinder(r=1.25) + box = openmc.rectangular_prism(4.0, 4.0, boundary_type='reflective') + cell1 = openmc.Cell(fill=fuel, region=-zcyl) + cell2 = openmc.Cell(fill=None, region=+zcyl & box) + model = openmc.Model() + model.geometry = openmc.Geometry([cell1, cell2]) + + model.settings.particles = 100 + model.settings.batches = 10 + model.settings.inactive = 0 + + cyl_mesh = openmc.CylindricalMesh() + cyl_mesh.r_grid = [0., 1.25] + cyl_mesh.phi_grid = [0., 2*math.pi] + cyl_mesh.z_grid = [-1e10, 1e10] + cyl_mesh_filter = openmc.MeshFilter(cyl_mesh) + cell_filter = openmc.CellFilter([cell1]) + + tally1 = openmc.Tally() + tally1.filters = [cyl_mesh_filter] + tally1.scores = ['flux'] + tally2 = openmc.Tally() + tally2.filters = [cell_filter] + tally2.scores = ['flux'] + model.tallies = openmc.Tallies([tally1, tally2]) + + # Run OpenMC + sp_filename = model.run() + + # Get flux for each of the two tallies + with openmc.StatePoint(sp_filename) as sp: + t1 = sp.tallies[tally1.id] + t2 = sp.tallies[tally2.id] + mean1 = t1.mean.ravel()[0] + mean2 = t2.mean.ravel()[0] + + # The two tallies should be exactly the same + assert mean1 == pytest.approx(mean2) + + +def test_spherical_mesh_coincident(run_in_tmpdir): + """Test for spherical mesh boundary being coincident with a cell boundary""" + + fuel = openmc.Material() + fuel.add_nuclide('U235', 1.) + fuel.set_density('g/cm3', 4.5) + + sph = openmc.Sphere(r=1.25) + rcc = openmc.model.RectangularParallelepiped( + -2.0, 2.0, -2.0, 2.0, -2.0, 2.0, boundary_type='reflective') + cell1 = openmc.Cell(fill=fuel, region=-sph) + cell2 = openmc.Cell(fill=None, region=+sph & -rcc) + model = openmc.Model() + model.geometry = openmc.Geometry([cell1, cell2]) + + model.settings.particles = 100 + model.settings.batches = 10 + model.settings.inactive = 0 + + sph_mesh = openmc.SphericalMesh() + sph_mesh.r_grid = [0., 1.25] + sph_mesh.phi_grid = [0., 2*math.pi] + sph_mesh.theta_grid = [0., math.pi] + sph_mesh_filter = openmc.MeshFilter(sph_mesh) + cell_filter = openmc.CellFilter([cell1]) + + tally1 = openmc.Tally() + tally1.filters = [sph_mesh_filter] + tally1.scores = ['flux'] + tally2 = openmc.Tally() + tally2.filters = [cell_filter] + tally2.scores = ['flux'] + model.tallies = openmc.Tallies([tally1, tally2]) + + # Run OpenMC + sp_filename = model.run() + + # Get flux for each of the two tallies + with openmc.StatePoint(sp_filename) as sp: + t1 = sp.tallies[tally1.id] + t2 = sp.tallies[tally2.id] + mean1 = t1.mean.ravel()[0] + mean2 = t2.mean.ravel()[0] + + # The two tallies should be exactly the same + assert mean1 == pytest.approx(mean2) + + def test_get_reshaped_data(run_in_tmpdir): """Test that expanding MeshFilter dimensions works as expected""" From 2aef18c02ac82c7c90d59466b199a58883018bd1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 25 Mar 2023 11:55:53 -0500 Subject: [PATCH 2/3] Loosen tolerance on coincidence check for cyl/sph meshes and parametrize tests --- src/mesh.cpp | 7 ++++--- tests/unit_tests/test_filter_mesh.py | 19 +++++++++++-------- 2 files changed, 15 insertions(+), 11 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 306e560707..3ee489a68c 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -500,7 +500,8 @@ void StructuredMesh::raytrace_mesh( } // translate start and end positions, - // this needs to come after the get_indices call because it does its own translation + // this needs to come after the get_indices call because it does its own + // translation local_coords(r0); local_coords(r1); @@ -1084,7 +1085,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 (-p - D > l && std::abs(c) > FP_COINCIDENT) + if (-p - D > l && std::abs(c) > 1e-10) return -p - D; if (-p + D > l) return -p + D; @@ -1310,7 +1311,7 @@ double SphericalMesh::find_r_crossing( 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) > FP_COINCIDENT) + if (-p - D > l && std::abs(c) > 1e-10) return -p - D; if (-p + D > l) return -p + D; diff --git a/tests/unit_tests/test_filter_mesh.py b/tests/unit_tests/test_filter_mesh.py index 5e5a84f7b0..e3fae15be9 100644 --- a/tests/unit_tests/test_filter_mesh.py +++ b/tests/unit_tests/test_filter_mesh.py @@ -114,15 +114,16 @@ def test_cylindrical_mesh_estimators(run_in_tmpdir): assert np.all(diff < 3*std_dev) -def test_cylindrical_mesh_coincident(run_in_tmpdir): +@pytest.mark.parametrize("scale", [1.0, 1e2, 1e4, 1e5]) +def test_cylindrical_mesh_coincident(scale, run_in_tmpdir): """Test for cylindrical mesh boundary being coincident with a cell boundary""" fuel = openmc.Material() fuel.add_nuclide('U235', 1.) fuel.set_density('g/cm3', 4.5) - zcyl = openmc.ZCylinder(r=1.25) - box = openmc.rectangular_prism(4.0, 4.0, boundary_type='reflective') + zcyl = openmc.ZCylinder(r=1.25*scale) + box = openmc.rectangular_prism(4*scale, 4*scale, boundary_type='reflective') cell1 = openmc.Cell(fill=fuel, region=-zcyl) cell2 = openmc.Cell(fill=None, region=+zcyl & box) model = openmc.Model() @@ -133,7 +134,7 @@ def test_cylindrical_mesh_coincident(run_in_tmpdir): model.settings.inactive = 0 cyl_mesh = openmc.CylindricalMesh() - cyl_mesh.r_grid = [0., 1.25] + cyl_mesh.r_grid = [0., 1.25*scale] cyl_mesh.phi_grid = [0., 2*math.pi] cyl_mesh.z_grid = [-1e10, 1e10] cyl_mesh_filter = openmc.MeshFilter(cyl_mesh) @@ -161,16 +162,18 @@ def test_cylindrical_mesh_coincident(run_in_tmpdir): assert mean1 == pytest.approx(mean2) -def test_spherical_mesh_coincident(run_in_tmpdir): +@pytest.mark.parametrize("scale", [1.0, 1e2, 1e4, 1e5]) +def test_spherical_mesh_coincident(scale, run_in_tmpdir): """Test for spherical mesh boundary being coincident with a cell boundary""" fuel = openmc.Material() fuel.add_nuclide('U235', 1.) fuel.set_density('g/cm3', 4.5) - sph = openmc.Sphere(r=1.25) + sph = openmc.Sphere(r=1.25*scale) rcc = openmc.model.RectangularParallelepiped( - -2.0, 2.0, -2.0, 2.0, -2.0, 2.0, boundary_type='reflective') + -2*scale, 2*scale, -2*scale, 2*scale, -2*scale, 2*scale, + boundary_type='reflective') cell1 = openmc.Cell(fill=fuel, region=-sph) cell2 = openmc.Cell(fill=None, region=+sph & -rcc) model = openmc.Model() @@ -181,7 +184,7 @@ def test_spherical_mesh_coincident(run_in_tmpdir): model.settings.inactive = 0 sph_mesh = openmc.SphericalMesh() - sph_mesh.r_grid = [0., 1.25] + sph_mesh.r_grid = [0., 1.25*scale] sph_mesh.phi_grid = [0., 2*math.pi] sph_mesh.theta_grid = [0., math.pi] sph_mesh_filter = openmc.MeshFilter(sph_mesh) From 5b7fb427d3cb4c5ba9fa2a5555004c10e2e53db1 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 25 Mar 2023 15:36:11 -0500 Subject: [PATCH 3/3] Apply @pshriwise suggestion adding 0.1 to scale Co-authored-by: Patrick Shriwise --- tests/unit_tests/test_filter_mesh.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/unit_tests/test_filter_mesh.py b/tests/unit_tests/test_filter_mesh.py index e3fae15be9..c93bfa0377 100644 --- a/tests/unit_tests/test_filter_mesh.py +++ b/tests/unit_tests/test_filter_mesh.py @@ -114,7 +114,7 @@ def test_cylindrical_mesh_estimators(run_in_tmpdir): assert np.all(diff < 3*std_dev) -@pytest.mark.parametrize("scale", [1.0, 1e2, 1e4, 1e5]) +@pytest.mark.parametrize("scale", [0.1, 1.0, 1e2, 1e4, 1e5]) def test_cylindrical_mesh_coincident(scale, run_in_tmpdir): """Test for cylindrical mesh boundary being coincident with a cell boundary""" @@ -162,7 +162,7 @@ def test_cylindrical_mesh_coincident(scale, run_in_tmpdir): assert mean1 == pytest.approx(mean2) -@pytest.mark.parametrize("scale", [1.0, 1e2, 1e4, 1e5]) +@pytest.mark.parametrize("scale", [0.1, 1.0, 1e2, 1e4, 1e5]) def test_spherical_mesh_coincident(scale, run_in_tmpdir): """Test for spherical mesh boundary being coincident with a cell boundary"""