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..4c0a43877 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, 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] + +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 @@ -149,3 +199,77 @@ 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) + elif 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