crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Add openmc_mesh_get_volumes C API function (#2869)

Browse source
This commit is contained in:
Paul Romano 2024-02-13 06:44:05 -06:00 • committed by GitHub
parent f14fc55e60
commit bf33a9e11b
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
5 changed files with 54 additions and 2 deletions
Download patch file Download diff file Expand all files Collapse all files

1
include/openmc/capi.h
View file

@ -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);

28
openmc/lib/mesh.py
View file

@ -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'

2
openmc/material.py
View file

@ -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:

11
src/mesh.cpp
View file

@ -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)
{

14
tests/unit_tests/test_lib.py
View file

@ -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)