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 3ee489a68..e72a5cf7d 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,7 +1088,10 @@ 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) > 1e-10) + if (std::abs(c) <= RADIAL_MESH_TOL) + return INFTY; + + if (-p - D > l) return -p - D; if (-p + D > l) return -p + D; @@ -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) <= RADIAL_MESH_TOL) + 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; 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(); diff --git a/tests/unit_tests/test_cylindrical_mesh.py b/tests/unit_tests/test_cylindrical_mesh.py index 83ba09cc5..4382b49ba 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 @@ -98,8 +98,99 @@ 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: 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] + cylinders = [openmc.ZCylinder(r=ri) for ri in radii] + cylinders[-1].boundary_type = 'vacuum' + + regions = openmc.model.subdivide(cylinders)[:-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() + # we expect the tally results to be the same as the mesh grid width + # for these cases + 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() + 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): + # 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 + + # can't source particle directly on the geometry boundary + outer_r = bbox[1][0] - 1e-08 + + n_sources = 100 + radial_vals = np.linspace(0.0, 2.0*np.pi, n_sources) + + 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) + 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_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..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 @@ -102,8 +102,108 @@ 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: 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() + # 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): + # 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 + + # 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 + 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] - 1e-08 + + 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) + # 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)