From c487090b3f6c19b9de1adade6c6ccb261c1e1e08 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 6 Mar 2023 23:38:00 -0600 Subject: [PATCH 1/4] Correcting spherical to cartesian coordinate conversion --- openmc/mesh.py | 8 +- tests/unit_tests/test_mesh_to_vtk.py | 156 ++++++++++++++++++++++++--- 2 files changed, 147 insertions(+), 17 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 5294b2dd6..a0f86314d 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -446,7 +446,7 @@ class RegularMesh(StructuredMesh): if self._width is not None: self._width = None warnings.warn("Unsetting width attribute.") - + if self.lower_left is not None and any(np.isclose(self.lower_left, upper_right)): raise ValueError("Mesh cannot have zero thickness in any dimension") @@ -1533,9 +1533,9 @@ class SphericalMesh(StructuredMesh): pts_cartesian = np.copy(pts_spherical) r, theta, phi = pts_spherical[:, 0], pts_spherical[:, 1], pts_spherical[:, 2] - pts_cartesian[:, 0] = r * np.sin(phi) * np.cos(theta) - pts_cartesian[:, 1] = r * np.sin(phi) * np.sin(theta) - pts_cartesian[:, 2] = r * np.cos(phi) + pts_cartesian[:, 0] = r * np.sin(theta) * np.cos(phi) + pts_cartesian[:, 1] = r * np.sin(theta) * np.sin(phi) + pts_cartesian[:, 2] = r * np.cos(theta) return super().write_data_to_vtk( points=pts_cartesian, diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index a9dcd7243..06a5606e7 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -1,3 +1,5 @@ +from itertools import product + import numpy as np from pathlib import Path import pytest @@ -7,29 +9,77 @@ from vtk.util import numpy_support as nps import openmc +@pytest.fixture +def model(): + openmc.reset_auto_ids() + + surf1 = openmc.Sphere(r=10, boundary_type='vacuum') + surf2 = openmc.XPlane(x0=-0.001, boundary_type='vacuum') + + cell = openmc.Cell(region=-surf1 & -surf2) + + geometry = openmc.Geometry([cell]) + + settings = openmc.Settings() + settings.batches = 2 + settings.particles = 100 + settings.run_mode = 'fixed source' + + source = openmc.Source() + source.angle = openmc.stats.Isotropic() + source.energy = openmc.stats.Discrete([1.0e6], [1.0]) + source.space = openmc.stats.Point((-0.01, -0.01, -0.01)) + + settings.source = source + + model = openmc.Model(geometry=geometry, settings=settings) + + return model + regular_mesh = openmc.RegularMesh() -regular_mesh.lower_left = (0, 0, 0) -regular_mesh.upper_right = (1, 1, 1) +regular_mesh.lower_left = (-10, -10, -10) +regular_mesh.upper_right = (10, 10, 10) regular_mesh.dimension = [30, 20, 10] rectilinear_mesh = openmc.RectilinearMesh() -rectilinear_mesh.x_grid = np.linspace(1, 2, num=30) -rectilinear_mesh.y_grid = np.linspace(1, 2, num=30) -rectilinear_mesh.z_grid = np.linspace(1, 2, num=30) +rectilinear_mesh.x_grid = np.linspace(-10, 10, 6) +rectilinear_mesh.y_grid = np.logspace(0, 1, 7) +rectilinear_mesh.y_grid = \ + np.concatenate((-rectilinear_mesh.y_grid[::-1], rectilinear_mesh.y_grid)) +rectilinear_mesh.z_grid = np.linspace(-10, 10, 11) cylinder_mesh = openmc.CylindricalMesh() -cylinder_mesh.r_grid = np.linspace(1, 2, num=30) -cylinder_mesh.phi_grid = np.linspace(0, np.pi, num=50) -cylinder_mesh.z_grid = np.linspace(0, 1, num=30) +cylinder_mesh.r_grid = np.linspace(0, 10, 23) +cylinder_mesh.phi_grid = np.linspace(0, np.pi, 21) +cylinder_mesh.z_grid = np.linspace(0, 1, 15) spherical_mesh = openmc.SphericalMesh() -spherical_mesh.r_grid = np.linspace(1, 2, num=30) -spherical_mesh.phi_grid = np.linspace(0, np.pi, num=50) -spherical_mesh.theta_grid = np.linspace(0, np.pi / 2, num=30) +spherical_mesh.r_grid = np.linspace(1, 10, 30) +spherical_mesh.phi_grid = np.linspace(0, np.pi, 25) +spherical_mesh.theta_grid = np.linspace(0, np.pi / 2, 15) + +MESHES = [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh] + +x_plane = openmc.XPlane(x0=-0.001, boundary_type='vacuum') +y_plane = openmc.YPlane(y0=-0.001, boundary_type='vacuum') +z_plane = openmc.ZPlane(z0=-0.001, boundary_type='vacuum') + +SURFS = [x_plane, y_plane, z_plane] -@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh]) +def ids(mesh): + if isinstance(mesh, openmc.CylindricalMesh): + return 'cylindrical_mesh' + elif isinstance(mesh, openmc.RegularMesh): + return 'regular_mesh' + elif isinstance(mesh, openmc.RectilinearMesh): + return 'rectilinear_mesh' + elif isinstance(mesh, openmc.SphericalMesh): + return 'spherical_mesh' + + +@pytest.mark.parametrize("mesh", MESHES, ids=ids) def test_write_data_to_vtk(mesh, tmpdir): # BUILD filename = Path(tmpdir) / "out.vtk" @@ -69,7 +119,7 @@ def test_write_data_to_vtk(mesh, tmpdir): assert all(data1 == 1.0) -@pytest.mark.parametrize("mesh", [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh]) +@pytest.mark.parametrize("mesh", MESHES, ids=ids) def test_write_data_to_vtk_size_mismatch(mesh): """Checks that an error is raised when the size of the dataset doesn't match the mesh number of cells @@ -92,6 +142,7 @@ def test_write_data_to_vtk_size_mismatch(mesh): with pytest.raises(ValueError, match=expected_error_msg): mesh.write_data_to_vtk(filename="out.vtk", datasets={"label": data}) +<<<<<<< HEAD def test_write_data_to_vtk_round_trip(run_in_tmpdir): cmesh = openmc.CylindricalMesh() cmesh.r_grid = (0.0, 1.0, 2.0) @@ -149,3 +200,82 @@ def test_write_data_to_vtk_round_trip(run_in_tmpdir): # checks that the vtk cell values are equal to the data assert np.array_equal(vtk_values, data) +======= + +def mesh_surf_id(param): + + if isinstance(param, openmc.MeshBase): + return ids(param) + + if isinstance(param, openmc.XPlane): + return 'XPlane' + elif isinstance(param, openmc.YPlane): + return 'YPlane' + elif isinstance(param, openmc.ZPlane): + return 'ZPlane' + + +@pytest.mark.parametrize("mesh,surface", product(MESHES, SURFS), ids=mesh_surf_id) +def test_vtk_write_ordering(model, mesh, surface): + + tally = openmc.Tally() + tally.scores = ['flux'] + # use the mesh on the specified tally + mesh_filter = openmc.MeshFilter(mesh) + tally.filters = [mesh_filter] + + model.tallies = openmc.Tallies([tally]) + + # run the problem + sp_filename = model.run() + + with openmc.StatePoint(sp_filename) as sp: + mean = sp.tallies[tally.id].mean + + # write the data to a VTK file + vtk_filename = 'test.vtk' + mesh.write_data_to_vtk(vtk_filename, datasets={'mean': mean}) + + # read file + reader = vtk.vtkStructuredGridReader() + reader.SetFileName(str(vtk_filename)) + reader.Update() + + # check name of datasets + vtk_grid = reader.GetOutput() + array = vtk_grid.GetCellData().GetArray(0) + vtk_data = nps.vtk_to_numpy(array) + + # convenience function for determining if a mesh + # element has vertices in the geometry. This + # particular geometry allows us to assume that tally results + # in the element should be zero if none of its vertices lie in the geometry + def in_geom(cell): + point_ids = cell.GetPointIds() + + for i in range(point_ids.GetNumberOfIds()): + p = vtk_grid.GetPoint(point_ids.GetId(i)) + if model.geometry.find(p): + return True + + return False + + # reshape mean according to mesh dimensions + mean = mean.reshape(mesh.dimension[::-1]).T + centroid = [0.0, 0.0, 0.0] + + # check that tally and vtk array results are zero where expected + for ijk in mesh.indices: + ijk = tuple(n - 1 for n in ijk) + # get the cell from the stuctured mesh object + cell = vtk_grid.GetCell(*ijk) + if not in_geom(cell): + cell.GetCentroid(centroid) + err_msg = f'IJK: {ijk} should be zero but is not. Centroid: {centroid}' + assert mean[ijk] == 0.0, err_msg + + # need to get flat index with axes reversed due to ordering passed into the VTK file + flat_idx = np.ravel_multi_index(tuple(ijk[::-1]), mesh.dimension[::-1]) + assert vtk_data[flat_idx] == 0.0, err_msg + +>>>>>>> 641319115 (Correcting spherical to cartesian coordinate conversion) From ae57ad0e5be45619f260a3cd3d3df911467ec628 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Mon, 6 Mar 2023 23:39:15 -0600 Subject: [PATCH 2/4] Improving test coverage for mesh to vtk capability --- tests/unit_tests/test_mesh_to_vtk.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index 06a5606e7..22cdf1c1e 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -56,7 +56,7 @@ cylinder_mesh.z_grid = np.linspace(0, 1, 15) spherical_mesh = openmc.SphericalMesh() spherical_mesh.r_grid = np.linspace(1, 10, 30) -spherical_mesh.phi_grid = np.linspace(0, np.pi, 25) +spherical_mesh.phi_grid = np.linspace(0, 0.8*np.pi, 25) spherical_mesh.theta_grid = np.linspace(0, np.pi / 2, 15) MESHES = [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh] From 285edc77b982e7abf9f89511ee2e82a59b686e53 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Wed, 8 Mar 2023 22:47:40 -0600 Subject: [PATCH 3/4] Corrections after rebase --- tests/unit_tests/test_mesh_to_vtk.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index 22cdf1c1e..80da85d97 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -142,7 +142,6 @@ def test_write_data_to_vtk_size_mismatch(mesh): with pytest.raises(ValueError, match=expected_error_msg): mesh.write_data_to_vtk(filename="out.vtk", datasets={"label": data}) -<<<<<<< HEAD def test_write_data_to_vtk_round_trip(run_in_tmpdir): cmesh = openmc.CylindricalMesh() cmesh.r_grid = (0.0, 1.0, 2.0) @@ -200,7 +199,6 @@ def test_write_data_to_vtk_round_trip(run_in_tmpdir): # checks that the vtk cell values are equal to the data assert np.array_equal(vtk_values, data) -======= def mesh_surf_id(param): @@ -277,5 +275,3 @@ def test_vtk_write_ordering(model, mesh, surface): # need to get flat index with axes reversed due to ordering passed into the VTK file flat_idx = np.ravel_multi_index(tuple(ijk[::-1]), mesh.dimension[::-1]) assert vtk_data[flat_idx] == 0.0, err_msg - ->>>>>>> 641319115 (Correcting spherical to cartesian coordinate conversion) From 9b708e7ad47292a73455e26dda134ea0c6086625 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Tue, 21 Mar 2023 14:55:42 -0500 Subject: [PATCH 4/4] Apply suggestions from @paulromano Co-authored-by: Paul Romano --- tests/unit_tests/test_mesh_to_vtk.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index 80da85d97..4c0a43877 100644 --- a/tests/unit_tests/test_mesh_to_vtk.py +++ b/tests/unit_tests/test_mesh_to_vtk.py @@ -59,7 +59,7 @@ spherical_mesh.r_grid = np.linspace(1, 10, 30) spherical_mesh.phi_grid = np.linspace(0, 0.8*np.pi, 25) spherical_mesh.theta_grid = np.linspace(0, np.pi / 2, 15) -MESHES = [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh] +MESHES = [cylinder_mesh, regular_mesh, rectilinear_mesh, spherical_mesh] x_plane = openmc.XPlane(x0=-0.001, boundary_type='vacuum') y_plane = openmc.YPlane(y0=-0.001, boundary_type='vacuum') @@ -201,11 +201,9 @@ def test_write_data_to_vtk_round_trip(run_in_tmpdir): assert np.array_equal(vtk_values, data) def mesh_surf_id(param): - if isinstance(param, openmc.MeshBase): return ids(param) - - if isinstance(param, openmc.XPlane): + elif isinstance(param, openmc.XPlane): return 'XPlane' elif isinstance(param, openmc.YPlane): return 'YPlane'