diff --git a/CMakeLists.txt b/CMakeLists.txt index 11ae71ac65..0dcaa6fb37 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -156,9 +156,9 @@ endif() find_package(HDF5 REQUIRED COMPONENTS C HL) if(HDF5_IS_PARALLEL) if(NOT OPENMC_USE_MPI) - message(FATAL_ERROR "Parallel HDF5 was detected, but the detected compiler,\ - ${CMAKE_CXX_COMPILER}, does not support MPI. An MPI-capable compiler must \ - be used with parallel HDF5.") + message(FATAL_ERROR "Parallel HDF5 was detected, but MPI was not enabled.\ + To use parallel HDF5, OpenMC needs to be built with MPI support by passing\ + -DOPENMC_USE_MPI=ON when calling cmake.") endif() message(STATUS "Using parallel HDF5") endif() diff --git a/openmc/material.py b/openmc/material.py index a1b3367a90..a2c6138d28 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -1029,7 +1029,8 @@ class Material(IDManagerMixin): return nuclides - def get_activity(self, units: str = 'Bq/cm3', by_nuclide: bool = False): + def get_activity(self, units: str = 'Bq/cm3', by_nuclide: bool = False, + volume: Optional[float] = None): """Returns the activity of the material or for each nuclide in the material in units of [Bq], [Bq/g] or [Bq/cm3]. @@ -1044,6 +1045,11 @@ class Material(IDManagerMixin): by_nuclide : bool Specifies if the activity should be returned for the material as a whole or per nuclide. Default is False. + volume : float, optional + Volume of the material. If not passed, defaults to using the + :attr:`Material.volume` attribute. + + .. versionadded:: 0.13.3 Returns ------- @@ -1057,7 +1063,7 @@ class Material(IDManagerMixin): cv.check_type('by_nuclide', by_nuclide, bool) if units == 'Bq': - multiplier = self.volume + multiplier = volume if volume is not None else self.volume elif units == 'Bq/cm3': multiplier = 1 elif units == 'Bq/g': @@ -1070,7 +1076,8 @@ class Material(IDManagerMixin): return activity if by_nuclide else sum(activity.values()) - def get_decay_heat(self, units: str = 'W', by_nuclide: bool = False): + def get_decay_heat(self, units: str = 'W', by_nuclide: bool = False, + volume: Optional[float] = None): """Returns the decay heat of the material or for each nuclide in the material in units of [W], [W/g] or [W/cm3]. @@ -1085,6 +1092,11 @@ class Material(IDManagerMixin): by_nuclide : bool Specifies if the decay heat should be returned for the material as a whole or per nuclide. Default is False. + volume : float, optional + Volume of the material. If not passed, defaults to using the + :attr:`Material.volume` attribute. + + .. versionadded:: 0.13.3 Returns ------- @@ -1098,7 +1110,7 @@ class Material(IDManagerMixin): cv.check_type('by_nuclide', by_nuclide, bool) if units == 'W': - multiplier = self.volume + multiplier = volume if volume is not None else self.volume elif units == 'W/cm3': multiplier = 1 elif units == 'W/g': @@ -1113,11 +1125,19 @@ class Material(IDManagerMixin): return decayheat if by_nuclide else sum(decayheat.values()) - def get_nuclide_atoms(self): + def get_nuclide_atoms(self, volume: Optional[float] = None): """Return number of atoms of each nuclide in the material .. versionadded:: 0.13.1 + Parameters + ---------- + volume : float, optional + Volume of the material. If not passed, defaults to using the + :attr:`Material.volume` attribute. + + .. versionadded:: 0.13.3 + Returns ------- dict @@ -1125,11 +1145,13 @@ class Material(IDManagerMixin): atoms present in the material. """ - if self.volume is None: + if volume is None: + volume = self.volume + if volume is None: raise ValueError("Volume must be set in order to determine atoms.") atoms = {} for nuclide, atom_per_bcm in self.get_nuclide_atom_densities().items(): - atoms[nuclide] = 1.0e24 * atom_per_bcm * self.volume + atoms[nuclide] = 1.0e24 * atom_per_bcm * volume return atoms def get_mass_density(self, nuclide: Optional[str] = None): @@ -1154,7 +1176,7 @@ class Material(IDManagerMixin): mass_density += density_i return mass_density - def get_mass(self, nuclide: Optional[str] = None): + def get_mass(self, nuclide: Optional[str] = None, volume: Optional[float] = None): """Return mass of one or all nuclides. Note that this method requires that the :attr:`Material.volume` has @@ -1165,6 +1187,12 @@ class Material(IDManagerMixin): nuclides : str, optional Nuclide for which mass is desired. If not specified, the density for the entire material is given. + volume : float, optional + Volume of the material. If not passed, defaults to using the + :attr:`Material.volume` attribute. + + .. versionadded:: 0.13.3 + Returns ------- @@ -1172,9 +1200,11 @@ class Material(IDManagerMixin): Mass of the nuclide/material in [g] """ - if self.volume is None: + if volume is None: + volume = self.volume + if volume is None: raise ValueError("Volume must be set in order to determine mass.") - return self.volume*self.get_mass_density(nuclide) + return volume*self.get_mass_density(nuclide) def clone(self, memo: Optional[dict] = None): """Create a copy of this material with a new unique ID. diff --git a/openmc/mesh.py b/openmc/mesh.py index 7b69de8a03..f606b85261 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") @@ -1578,6 +1578,7 @@ 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) + self.origin[0] pts_cartesian[:, 1] = r * np.sin(phi) * np.sin(theta) + self.origin[1] pts_cartesian[:, 2] = r * np.cos(phi) + self.origin[2] diff --git a/openmc/model/model.py b/openmc/model/model.py index 79708d62b3..3e0422eb1e 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -244,7 +244,7 @@ class Model: def from_model_xml(cls, path='model.xml'): """Create model from single XML file - .. vesionadded:: 0.13.3 + .. versionadded:: 0.13.3 Parameters ---------- @@ -740,7 +740,7 @@ class Model: if len(self.settings.volume_calculations) == 0: # Then there is no volume calculation specified - raise ValueError("The Settings.volume_calculation attribute must" + raise ValueError("The Settings.volume_calculations attribute must" " be specified before executing this method!") with _change_directory(Path(cwd)): diff --git a/tests/unit_tests/test_material.py b/tests/unit_tests/test_material.py index e435c7493a..3b93b6879a 100644 --- a/tests/unit_tests/test_material.py +++ b/tests/unit_tests/test_material.py @@ -377,6 +377,9 @@ def test_get_nuclide_atoms(): atoms = mat.get_nuclide_atoms() assert atoms['Li6'] == pytest.approx(mat.density * mat.volume) + atoms = mat.get_nuclide_atoms(volume=10.0) + assert atoms['Li6'] == pytest.approx(mat.density * 10.0) + def test_mass(): m = openmc.Material() @@ -394,6 +397,9 @@ def test_mass(): assert m.get_mass() == pytest.approx(20.0) assert m.fissionable_mass == pytest.approx(10.0) + # Test with volume specified as argument + assert m.get_mass('Zr90', volume=1.0) == pytest.approx(1.0) + def test_materials(run_in_tmpdir): m1 = openmc.Material() @@ -544,11 +550,14 @@ def test_get_activity(): m4.volume = 10. assert pytest.approx(m4.get_activity(units='Bq')) == 355978108155965.94*3/2*10 # [Bq] + # Test with volume specified as argument + assert pytest.approx(m4.get_activity(units='Bq', volume=1.0)) == 355978108155965.94*3/2 + def test_get_decay_heat(): # Set chain file for testing openmc.config['chain_file'] = Path(__file__).parents[1] / 'chain_simple.xml' - + """Tests the decay heat of stable, metastable and active materials""" m1 = openmc.Material() m1.add_nuclide("U235", 0.2) @@ -589,7 +598,10 @@ def test_get_decay_heat(): # volume is required to calculate total decay heat m4.volume = 10. assert pytest.approx(m4.get_decay_heat(units='W')) == 40175.15720273193*3/2*10 # [W] - + + # Test with volume specified as argument + assert pytest.approx(m4.get_decay_heat(units='W', volume=1.0)) == 40175.15720273193*3/2 + def test_decay_photon_energy(): # Set chain file for testing diff --git a/tests/unit_tests/test_mesh_to_vtk.py b/tests/unit_tests/test_mesh_to_vtk.py index a9dcd7243f..4c0a438774 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