diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 063ed8d9fa..5f98152cd6 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -106,6 +106,7 @@ int openmc_mesh_filter_set_translation(int32_t index, double translation[3]); int openmc_mesh_get_id(int32_t index, int32_t* id); int openmc_mesh_set_id(int32_t index, int32_t id); int openmc_mesh_get_n_elements(int32_t index, size_t* n); +int openmc_mesh_get_volumes(int32_t index, double* volumes); int openmc_mesh_material_volumes(int32_t index, int n_sample, int bin, int result_size, void* result, int* hits, uint64_t* seed); int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh); diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py index 599b750fb0..7f17dd81ff 100644 --- a/openmc/lib/mesh.py +++ b/openmc/lib/mesh.py @@ -39,6 +39,9 @@ _dll.openmc_mesh_set_id.errcheck = _error_handler _dll.openmc_mesh_get_n_elements.argtypes = [c_int32, POINTER(c_size_t)] _dll.openmc_mesh_get_n_elements.restype = c_int _dll.openmc_mesh_get_n_elements.errcheck = _error_handler +_dll.openmc_mesh_get_volumes.argtypes = [c_int32, POINTER(c_double)] +_dll.openmc_mesh_get_volumes.restype = c_int +_dll.openmc_mesh_get_volumes.errcheck = _error_handler _dll.openmc_mesh_material_volumes.argtypes = [ c_int32, c_int, c_int, c_int, POINTER(_MaterialVolume), POINTER(c_int), POINTER(c_uint64)] @@ -149,11 +152,18 @@ class Mesh(_FortranObjectWithID): _dll.openmc_mesh_set_id(self._index, mesh_id) @property - def n_elements(self): + def n_elements(self) -> int: n = c_size_t() _dll.openmc_mesh_get_n_elements(self._index, n) return n.value + @property + def volumes(self) -> np.ndarray: + volumes = np.empty((self.n_elements,)) + _dll.openmc_mesh_get_volumes( + self._index, volumes.ctypes.data_as(POINTER(c_double))) + return volumes + def material_volumes( self, n_samples: int = 10_000, @@ -276,6 +286,10 @@ class RegularMesh(Mesh): are given, it is assumed that the mesh is an x-y mesh. width : numpy.ndarray The width of mesh cells in each direction. + n_elements : int + Total number of mesh elements. + volumes : numpy.ndarray + Volume of each mesh element in [cm^3] """ mesh_type = 'regular' @@ -358,6 +372,10 @@ class RectilinearMesh(Mesh): The upper-right corner of the structrued mesh. width : numpy.ndarray The width of mesh cells in each direction. + n_elements : int + Total number of mesh elements. + volumes : numpy.ndarray + Volume of each mesh element in [cm^3] """ mesh_type = 'rectilinear' @@ -457,6 +475,10 @@ class CylindricalMesh(Mesh): The upper-right corner of the structrued mesh. width : numpy.ndarray The width of mesh cells in each direction. + n_elements : int + Total number of mesh elements. + volumes : numpy.ndarray + Volume of each mesh element in [cm^3] """ mesh_type = 'cylindrical' @@ -555,6 +577,10 @@ class SphericalMesh(Mesh): The upper-right corner of the structrued mesh. width : numpy.ndarray The width of mesh cells in each direction. + n_elements : int + Total number of mesh elements. + volumes : numpy.ndarray + Volume of each mesh element in [cm^3] """ mesh_type = 'spherical' diff --git a/openmc/material.py b/openmc/material.py index d42c489aee..c31862a18a 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -571,7 +571,7 @@ class Material(IDManagerMixin): for component, params in components.items(): cv.check_type('component', component, str) - if isinstance(params, float): + if isinstance(params, Real): params = {'percent': params} else: diff --git a/src/mesh.cpp b/src/mesh.cpp index f343b27e33..bc6781ad31 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1881,6 +1881,17 @@ extern "C" int openmc_mesh_get_n_elements(int32_t index, size_t* n) return 0; } +//! Get the volume of each element in the mesh +extern "C" int openmc_mesh_get_volumes(int32_t index, double* volumes) +{ + if (int err = check_mesh(index)) + return err; + for (int i = 0; i < model::meshes[index]->n_bins(); ++i) { + volumes[i] = model::meshes[index]->volume(i); + } + return 0; +} + extern "C" int openmc_mesh_material_volumes(int32_t index, int n_sample, int bin, int result_size, void* result, int* hits, uint64_t* seed) { diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index 6ae13ef76e..1310a6a7fa 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -566,6 +566,8 @@ def test_regular_mesh(lib_init): assert mesh.upper_right == pytest.approx(ur) assert mesh.width == pytest.approx(width) + np.testing.assert_allclose(mesh.volumes, 1.0) + meshes = openmc.lib.meshes assert isinstance(meshes, Mapping) assert len(meshes) == 1 @@ -644,6 +646,8 @@ def test_rectilinear_mesh(lib_init): for k, diff_z in enumerate(np.diff(z_grid)): assert np.all(mesh.width[i, j, k, :] == (10, 10, 10)) + np.testing.assert_allclose(mesh.volumes, 1000.0) + with pytest.raises(exc.AllocationError): mesh2 = openmc.lib.RectilinearMesh(mesh.id) @@ -688,6 +692,9 @@ def test_cylindrical_mesh(lib_init): for k, _ in enumerate(np.diff(z_grid)): assert np.allclose(mesh.width[i, j, k, :], (5, deg2rad(10), 10)) + np.testing.assert_allclose(mesh.volumes[::2], 10/360 * pi * 5**2 * 10) + np.testing.assert_allclose(mesh.volumes[1::2], 10/360 * pi * (10**2 - 5**2) * 10) + with pytest.raises(exc.AllocationError): mesh2 = openmc.lib.CylindricalMesh(mesh.id) @@ -734,6 +741,13 @@ def test_spherical_mesh(lib_init): for k, _ in enumerate(np.diff(phi_grid)): assert np.allclose(mesh.width[i, j, k, :], (5, deg2rad(10), deg2rad(10))) + dtheta = lambda d1, d2: np.cos(deg2rad(d1)) - np.cos(deg2rad(d2)) + f = 1/3 * deg2rad(10.) + np.testing.assert_allclose(mesh.volumes[::4], f * 5**3 * dtheta(0., 10.)) + np.testing.assert_allclose(mesh.volumes[1::4], f * (10**3 - 5**3) * dtheta(0., 10.)) + np.testing.assert_allclose(mesh.volumes[2::4], f * 5**3 * dtheta(10., 20.)) + np.testing.assert_allclose(mesh.volumes[3::4], f * (10**3 - 5**3) * dtheta(10., 20.)) + with pytest.raises(exc.AllocationError): mesh2 = openmc.lib.SphericalMesh(mesh.id)