crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 14:15:42 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

First pass at exposing rectilinear mesh through C API

Browse source
This commit is contained in:
Shikhar Kumar 2021-02-08 22:25:17 -05:00
parent 86e7403a21
commit f83ab39b99
6 changed files with 537 additions and 219 deletions
Download patch file Download diff file Expand all files Collapse all files

16
include/openmc/capi.h
View file

@ -28,7 +28,8 @@ extern "C" {
int openmc_extend_cells(int32_t n, int32_t* index_start, int32_t* index_end);
int openmc_extend_filters(int32_t n, int32_t* index_start, int32_t* index_end);
int openmc_extend_materials(int32_t n, int32_t* index_start, int32_t* index_end);
int openmc_extend_meshes(int32_t n, int32_t* index_start, int32_t* index_end);
int openmc_extend_meshes(int32_t n, const char* type, int32_t* index_start,
int32_t* index_end);
int openmc_extend_tallies(int32_t n, int32_t* index_start, int32_t* index_end);
int openmc_filter_get_id(int32_t index, int32_t* id);
int openmc_filter_get_type(int32_t index, char* type);
@ -73,11 +74,7 @@ extern "C" {
int openmc_mesh_filter_get_mesh(int32_t index, int32_t* index_mesh);
int openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh);
int openmc_mesh_get_id(int32_t index, int32_t* id);
int openmc_mesh_get_dimension(int32_t index, int** id, int* n);
int openmc_mesh_get_params(int32_t index, double** ll, double** ur, double** width, int* n);
int openmc_mesh_set_id(int32_t index, int32_t id);
int openmc_mesh_set_dimension(int32_t index, int n, const int* dims);
int openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur, const double* width);
int openmc_meshsurface_filter_get_mesh(int32_t index, int32_t* index_mesh);
int openmc_meshsurface_filter_set_mesh(int32_t index, int32_t index_mesh);
int openmc_new_filter(const char* type, int32_t* index);
@ -86,6 +83,15 @@ extern "C" {
int openmc_plot_geometry();
int openmc_id_map(const void* slice, int32_t* data_out);
int openmc_property_map(const void* slice, double* data_out);
int openmc_rectilinear_mesh_get_grid(int32_t index, double** grid_x, int* nx,
double** grid_y, int * ny, double** grid_z, int* nz);
int openmc_rectilinear_mesh_set_grid(int32_t index, const double* grid_x,
const int nx, const double* grid_y, const int ny,
const double* grid_z, const int nz);
int openmc_regular_mesh_get_dimension(int32_t index, int** id, int* n);
int openmc_regular_mesh_get_params(int32_t index, double** ll, double** ur, double** width, int* n);
int openmc_regular_mesh_set_dimension(int32_t index, int n, const int* dims);
int openmc_regular_mesh_set_params(int32_t index, int n, const double* ll, const double* ur, const double* width);
int openmc_reset();
int openmc_reset_timers();
int openmc_run();

28
include/openmc/mesh.h
View file

@ -46,6 +46,8 @@ public:
Mesh(pugi::xml_node node);
virtual ~Mesh() = default;
virtual std::string type() const = 0;
// Methods
//! Determine which bins were crossed by a particle
@ -114,6 +116,14 @@ public:
void bins_crossed(const Particle& p, std::vector<int>& bins,
std::vector<double>& lengths) const override;
//! Count number of bank sites in each mesh bin / energy bin
//
//! \param[in] Pointer to bank sites
//! \param[in] Number of bank sites
//! \param[out] Whether any bank sites are outside the mesh
xt::xtensor<double, 1> count_sites(const Particle::Bank* bank,
int64_t length, bool* outside) const;
//! Get bin given mesh indices
//
//! \param[in] Array of mesh indices
@ -186,6 +196,8 @@ public:
// Overriden methods
std::string type() const override {return "regular";}
void surface_bins_crossed(const Particle& p, std::vector<int>& bins)
const override;
@ -200,16 +212,6 @@ public:
void to_hdf5(hid_t group) const override;
// New methods
//! Count number of bank sites in each mesh bin / energy bin
//
//! \param[in] bank Array of bank sites
//! \param[out] Whether any bank sites are outside the mesh
//! \return Array indicating number of sites in each mesh/energy bin
xt::xtensor<double, 1> count_sites(const Particle::Bank* bank, int64_t length,
bool* outside) const;
// Data members
double volume_frac_; //!< Volume fraction of each mesh element
@ -221,10 +223,13 @@ class RectilinearMesh : public StructuredMesh
{
public:
// Constructors
RectilinearMesh() = default;
RectilinearMesh(pugi::xml_node node);
// Overriden methods
std::string type() const override {return "rectilinear";}
void surface_bins_crossed(const Particle& p, std::vector<int>& bins)
const override;
@ -239,7 +244,6 @@ public:
void to_hdf5(hid_t group) const override;
private:
std::vector<std::vector<double>> grid_;
};
@ -252,6 +256,8 @@ public:
UnstructuredMesh(pugi::xml_node);
~UnstructuredMesh() = default;
std::string type() const override {return "unstructured";}
void bins_crossed(const Particle& p,
std::vector<int>& bins,
std::vector<double>& lengths) const override;

201
openmc/cmfd.py
View file

@ -60,15 +60,17 @@ class CMFDMesh:
Attributes
----------
lower_left : Iterable of float
The lower-left corner of the structured mesh. If only two coordinates
are given, it is assumed that the mesh is an x-y mesh.
The lower-left corner of a regular structured mesh. If only two
coordinates are given, it is assumed that the mesh is an x-y mesh.
upper_right : Iterable of float
The upper-right corner of the structured mesh. If only two coordinates
are given, it is assumed that the mesh is an x-y mesh.
The upper-right corner of a regular structured mesh. If only two
coordinates are given, it is assumed that the mesh is an x-y mesh.
dimension : Iterable of int
The number of mesh cells in each direction.
The number of mesh cells in each direction for a regular structured
mesh.
width : Iterable of float
The width of mesh cells in each direction.
The width of mesh cells in each direction for a regular structured
mesh
energy : Iterable of float
Energy bins in eV, listed in ascending order (e.g. [0.0, 0.625e-1,
20.0e6]) for CMFD tallies and acceleration. If no energy bins are
@ -96,6 +98,14 @@ class CMFDMesh:
is extremely important to use in reflectors as neutrons will not
contribute to any tallies far away from fission source neutron regions.
A ``1`` must be used to identify any fission source region.
mesh_type : str
Type of structured mesh to use. Acceptable values are:
* "regular" - Use RegularMesh to define CMFD mesh
* "rectilinear" - Use RectilinearMesh to define CMFD
grid : Iterable of Iterable of float
Grid used to define RectilinearMesh. First dimension must have length
3 where grid[0], grid[1], and grid[2] correspond to the x-, y-, and
z-grids respectively
"""
@ -107,6 +117,28 @@ class CMFDMesh:
self._energy = None
self._albedo = None
self._map = None
self._mesh_type = 'regular'
self._grid = None
def __repr__(self):
outstr = type(self).__name__ + '\n'
if self._mesh_type == 'regular':
outstr += (self._get_repr(self._lower_left, "Lower left") + "\n" +
self._get_repr(self._upper_right, "Upper right") + "\n" +
self._get_repr(self._dimension, "Dimension") + "\n" +
self._get_repr(self._width, "Width") + "\n" +
self._get_repr(self._albedo, "Albedo"))
elif self._mesh_type == 'rectilinear':
outstr += (self._get_repr(self._grid[0], "X-grid") + "\n" +
self._get_repr(self._grid[1], "Y-grid") + "\n" +
self._get_repr(self._grid[2], "Z-grid"))
return outstr
def _get_repr(self, list_var, label):
outstr = "\t{:<11} = ".format(label)
if list(list_var):
outstr += ", ".join(str(i) for i in list_var)
return outstr
@property
def lower_left(self):
@ -136,17 +168,27 @@ class CMFDMesh:
def map(self):
return self._map
@property
def mesh_type(self):
return self._mesh_type
@property
def grid(self):
return self._grid
@lower_left.setter
def lower_left(self, lower_left):
check_type('CMFD mesh lower_left', lower_left, Iterable, Real)
check_length('CMFD mesh lower_left', lower_left, 2, 3)
self._lower_left = lower_left
self._display_mesh_warning('regular', 'CMFD mesh lower_left')
@upper_right.setter
def upper_right(self, upper_right):
check_type('CMFD mesh upper_right', upper_right, Iterable, Real)
check_length('CMFD mesh upper_right', upper_right, 2, 3)
self._upper_right = upper_right
self._display_mesh_warning('regular', 'CMFD mesh upper_right')
@dimension.setter
def dimension(self, dimension):
@ -163,6 +205,7 @@ class CMFDMesh:
for w in width:
check_greater_than('CMFD mesh width', w, 0)
self._width = width
self._display_mesh_warning('regular', 'CMFD mesh width')
@energy.setter
def energy(self, energy):
@ -187,6 +230,31 @@ class CMFDMesh:
check_value('CMFD mesh map', m, [0, 1])
self._map = mesh_map
@mesh_type.setter
def mesh_type(self, mesh_type):
check_value('CMFD mesh type', mesh_type, ['regular', 'rectilinear'])
self._mesh_type = mesh_type
@grid.setter
def grid(self, grid):
grid_length = 3
dims = ['x', 'y', 'z']
check_length('CMFD mesh grid', grid, grid_length)
for i in range(grid_length):
check_type('CMFD mesh {}-grid'.format(dims[i]), grid[i], Iterable,
Real)
check_greater_than('CMFD mesh {}-grid length'.format(dims[i]),
len(grid[i]), 1)
self._grid = np.array(grid)
self._display_mesh_warning('rectilinear', 'CMFD mesh grid')
def _display_mesh_warning(self, mesh_type, variable_label):
if self._mesh_type != mesh_type:
warn_msg = 'Setting {} if mesh type is not set to {} ' \
'will have no effect'.format(variable_label, mesh_type)
warnings.warn(warn_msg, RuntimeWarning)
class CMFDRun:
r"""Class for running CMFD acceleration through the C API.
@ -566,48 +634,56 @@ class CMFDRun:
def mesh(self, cmfd_mesh):
check_type('CMFD mesh', cmfd_mesh, CMFDMesh)
# Check dimension defined
if cmfd_mesh.dimension is None:
raise ValueError('CMFD mesh requires spatial '
'dimensions to be specified')
if cmfd_mesh.mesh_type == 'regular':
# Check dimension defined
if cmfd_mesh.dimension is None:
raise ValueError('CMFD regular mesh requires spatial '
'dimensions to be specified')
# Check lower left defined
if cmfd_mesh.lower_left is None:
raise ValueError('CMFD mesh requires lower left coordinates '
'to be specified')
# Check lower left defined
if cmfd_mesh.lower_left is None:
raise ValueError('CMFD regular mesh requires lower left '
'coordinates to be specified')
# Check that both upper right and width both not defined
if cmfd_mesh.upper_right is not None and cmfd_mesh.width is not None:
raise ValueError('Both upper right coordinates and width '
'cannot be specified for CMFD mesh')
# Check that both upper right and width both not defined
if cmfd_mesh.upper_right is not None and cmfd_mesh.width is not None:
raise ValueError('Both upper right coordinates and width '
'cannot be specified for CMFD regular mesh')
# Check that at least one of width or upper right is defined
if cmfd_mesh.upper_right is None and cmfd_mesh.width is None:
raise ValueError('CMFD mesh requires either upper right '
'coordinates or width to be specified')
# Check that at least one of width or upper right is defined
if cmfd_mesh.upper_right is None and cmfd_mesh.width is None:
raise ValueError('CMFD regular mesh requires either upper right '
'coordinates or width to be specified')
# Check width and lower length are same dimension and define
# upper_right
if cmfd_mesh.width is not None:
check_length('CMFD mesh width', cmfd_mesh.width,
len(cmfd_mesh.lower_left))
cmfd_mesh.upper_right = np.array(cmfd_mesh.lower_left) + \
np.array(cmfd_mesh.width) * np.array(cmfd_mesh.dimension)
# Check width and lower length are same dimension and define
# upper_right
if cmfd_mesh.width is not None:
check_length('CMFD mesh width', cmfd_mesh.width,
len(cmfd_mesh.lower_left))
cmfd_mesh.upper_right = np.array(cmfd_mesh.lower_left) + \
np.array(cmfd_mesh.width) * np.array(cmfd_mesh.dimension)
# Check upper_right and lower length are same dimension and define
# width
elif cmfd_mesh.upper_right is not None:
check_length('CMFD mesh upper right', cmfd_mesh.upper_right,
len(cmfd_mesh.lower_left))
# Check upper right coordinates are greater than lower left
if np.any(np.array(cmfd_mesh.upper_right) <=
np.array(cmfd_mesh.lower_left)):
raise ValueError('CMFD regular mesh requires upper right '
'coordinates to be greater than lower '
'left coordinates')
cmfd_mesh.width = np.true_divide((np.array(cmfd_mesh.upper_right) -
np.array(cmfd_mesh.lower_left)),
np.array(cmfd_mesh.dimension))
elif cmfd_mesh.mesh_type == 'rectilinear':
# Check dimension defined
if cmfd_mesh.grid is None:
raise ValueError('CMFD rectilinear mesh requires spatial '
'grid to be specified')
cmfd_mesh.dimension = [len(cmfd_mesh.grid[i]) - 1 for i in range(3)]
# Check upper_right and lower length are same dimension and define
# width
elif cmfd_mesh.upper_right is not None:
check_length('CMFD mesh upper right', cmfd_mesh.upper_right,
len(cmfd_mesh.lower_left))
# Check upper right coordinates are greater than lower left
if np.any(np.array(cmfd_mesh.upper_right) <=
np.array(cmfd_mesh.lower_left)):
raise ValueError('CMFD mesh requires upper right '
'coordinates to be greater than lower '
'left coordinates')
cmfd_mesh.width = np.true_divide((np.array(cmfd_mesh.upper_right) -
np.array(cmfd_mesh.lower_left)),
np.array(cmfd_mesh.dimension))
self._mesh = cmfd_mesh
@norm.setter
@ -1088,10 +1164,17 @@ class CMFDRun:
# Overwrite CMFD mesh properties
cmfd_mesh_name = 'mesh ' + str(cmfd_group.attrs['mesh_id'])
cmfd_mesh = f['tallies']['meshes'][cmfd_mesh_name]
self._mesh.dimension = cmfd_mesh['dimension'][()]
self._mesh.lower_left = cmfd_mesh['lower_left'][()]
self._mesh.upper_right = cmfd_mesh['upper_right'][()]
self._mesh.width = cmfd_mesh['width'][()]
self._mesh.mesh_type = cmfd_mesh['type'][()].decode()
if self._mesh.mesh_type == 'regular':
self._mesh.dimension = cmfd_mesh['dimension'][()]
self._mesh.lower_left = cmfd_mesh['lower_left'][()]
self._mesh.upper_right = cmfd_mesh['upper_right'][()]
self._mesh.width = cmfd_mesh['width'][()]
elif self._mesh.mesh_type == 'rectilinear':
x_grid = cmfd_mesh['x_grid'][()]
y_grid = cmfd_mesh['y_grid'][()]
z_grid = cmfd_mesh['z_grid'][()]
self._mesh.grid = [x_grid, y_grid, z_grid]
# Define variables that exist only on master process
if openmc.lib.master():
@ -1370,8 +1453,7 @@ class CMFDRun:
# Compute fission source
cmfd_src = (np.sum(self._nfissxs[:,:,:,:,:] *
cmfd_flux[:,:,:,:,np.newaxis], axis=3) *
vol[:,:,:,np.newaxis])
cmfd_flux[:,:,:,:,np.newaxis], axis=3))
# Normalize source such that it sums to 1.0
self._cmfd_src = cmfd_src / np.sum(cmfd_src)
@ -2804,15 +2886,22 @@ class CMFDRun:
def _create_cmfd_tally(self):
"""Creates all tallies in-memory that are used to solve CMFD problem"""
# Create Mesh object based on CMFDMesh, stored internally
cmfd_mesh = openmc.lib.RegularMesh()
# Create Mesh object based on CMFDMesh mesh_type, stored internally
if self._mesh.mesh_type == 'regular':
cmfd_mesh = openmc.lib.RegularMesh()
# Set dimension and parameters of mesh object
cmfd_mesh.dimension = self._mesh.dimension
cmfd_mesh.set_parameters(lower_left=self._mesh.lower_left,
upper_right=self._mesh.upper_right,
width=self._mesh.width)
elif self._mesh.mesh_type == 'rectilinear':
cmfd_mesh = openmc.lib.RectilinearMesh()
# Set grid of mesh object
x_grid, y_grid, z_grid = self._mesh.grid
cmfd_mesh.set_grid(x_grid, y_grid, z_grid)
# Store id of mesh object
self._mesh_id = cmfd_mesh.id
# Set dimension and parameters of mesh object
cmfd_mesh.dimension = self._mesh.dimension
cmfd_mesh.set_parameters(lower_left=self._mesh.lower_left,
upper_right=self._mesh.upper_right,
width=self._mesh.width)
# Create mesh Filter object, stored internally
mesh_filter = openmc.lib.MeshFilter()

6
openmc/lib/filter.py
View file

@ -11,7 +11,7 @@ from . import _dll
from .core import _FortranObjectWithID
from .error import _error_handler
from .material import Material
from .mesh import RegularMesh
from .mesh import _get_mesh
__all__ = [
@ -315,7 +315,7 @@ class MeshFilter(Filter):
def mesh(self):
index_mesh = c_int32()
_dll.openmc_mesh_filter_get_mesh(self._index, index_mesh)
return RegularMesh(index=index_mesh.value)
return _get_mesh(index_mesh.value)
@mesh.setter
def mesh(self, mesh):
@ -334,7 +334,7 @@ class MeshSurfaceFilter(Filter):
def mesh(self):
index_mesh = c_int32()
_dll.openmc_meshsurface_filter_get_mesh(self._index, index_mesh)
return RegularMesh(index=index_mesh.value)
return _get_mesh(index_mesh.value)
@mesh.setter
def mesh(self, mesh):

246
openmc/lib/mesh.py
View file

@ -1,5 +1,6 @@
from collections.abc import Mapping
from ctypes import c_int, c_int32, c_double, POINTER
from ctypes import (c_int, c_int32, c_char_p, c_double, POINTER,
create_string_buffer)
from weakref import WeakValueDictionary
from numpy.ctypeslib import as_array
@ -9,38 +10,97 @@ from . import _dll
from .core import _FortranObjectWithID
from .error import _error_handler
__all__ = ['RegularMesh', 'meshes']
__all__ = ['RegularMesh', 'RectilinearMesh', 'meshes']
# Mesh functions
_dll.openmc_extend_meshes.argtypes = [c_int32, POINTER(c_int32), POINTER(c_int32)]
_dll.openmc_extend_meshes.argtypes = [c_int32, c_char_p, POINTER(c_int32),
POINTER(c_int32)]
_dll.openmc_extend_meshes.restype = c_int
_dll.openmc_extend_meshes.errcheck = _error_handler
_dll.openmc_mesh_get_id.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_mesh_get_id.restype = c_int
_dll.openmc_mesh_get_id.errcheck = _error_handler
_dll.openmc_mesh_get_dimension.argtypes = [c_int32, POINTER(POINTER(c_int)), POINTER(c_int)]
_dll.openmc_mesh_get_dimension.restype = c_int
_dll.openmc_mesh_get_dimension.errcheck = _error_handler
_dll.openmc_mesh_get_params.argtypes = [
c_int32, POINTER(POINTER(c_double)), POINTER(POINTER(c_double)),
POINTER(POINTER(c_double)), POINTER(c_int)]
_dll.openmc_mesh_get_params.restype = c_int
_dll.openmc_mesh_get_params.errcheck = _error_handler
_dll.openmc_mesh_set_id.argtypes = [c_int32, c_int32]
_dll.openmc_mesh_set_id.restype = c_int
_dll.openmc_mesh_set_id.errcheck = _error_handler
_dll.openmc_mesh_set_dimension.argtypes = [c_int32, c_int, POINTER(c_int)]
_dll.openmc_mesh_set_dimension.restype = c_int
_dll.openmc_mesh_set_dimension.errcheck = _error_handler
_dll.openmc_mesh_set_params.argtypes = [
c_int32, c_int, POINTER(c_double), POINTER(c_double), POINTER(c_double)]
_dll.openmc_mesh_set_params.restype = c_int
_dll.openmc_mesh_set_params.errcheck = _error_handler
_dll.openmc_get_mesh_index.argtypes = [c_int32, POINTER(c_int32)]
_dll.openmc_get_mesh_index.restype = c_int
_dll.openmc_get_mesh_index.errcheck = _error_handler
_dll.n_meshes.argtypes = []
_dll.n_meshes.restype = c_int
_dll.openmc_rectilinear_mesh_get_grid.argtypes = [c_int32,
POINTER(POINTER(c_double)), POINTER(c_int), POINTER(POINTER(c_double)),
POINTER(c_int), POINTER(POINTER(c_double)), POINTER(c_int)]
_dll.openmc_rectilinear_mesh_get_grid.restype = c_int
_dll.openmc_rectilinear_mesh_get_grid.errcheck = _error_handler
_dll.openmc_rectilinear_mesh_set_grid.argtypes = [c_int32, POINTER(c_double),
c_int, POINTER(c_double), c_int, POINTER(c_double), c_int]
_dll.openmc_rectilinear_mesh_set_grid.restype = c_int
_dll.openmc_rectilinear_mesh_set_grid.errcheck = _error_handler
_dll.openmc_regular_mesh_get_dimension.argtypes = [c_int32,
POINTER(POINTER(c_int)), POINTER(c_int)]
_dll.openmc_regular_mesh_get_dimension.restype = c_int
_dll.openmc_regular_mesh_get_dimension.errcheck = _error_handler
_dll.openmc_regular_mesh_get_params.argtypes = [
c_int32, POINTER(POINTER(c_double)), POINTER(POINTER(c_double)),
POINTER(POINTER(c_double)), POINTER(c_int)]
_dll.openmc_regular_mesh_get_params.restype = c_int
_dll.openmc_regular_mesh_get_params.errcheck = _error_handler
_dll.openmc_regular_mesh_set_dimension.argtypes = [c_int32, c_int,
POINTER(c_int)]
_dll.openmc_regular_mesh_set_dimension.restype = c_int
_dll.openmc_regular_mesh_set_dimension.errcheck = _error_handler
_dll.openmc_regular_mesh_set_params.argtypes = [
c_int32, c_int, POINTER(c_double), POINTER(c_double), POINTER(c_double)]
_dll.openmc_regular_mesh_set_params.restype = c_int
_dll.openmc_regular_mesh_set_params.errcheck = _error_handler
class Mesh(_FortranObjectWithID):
"""Base class to represent mesh objects
"""
__instances = WeakValueDictionary()
def __new__(cls, uid=None, new=True, index=None):
mapping = meshes
if index is None:
if new:
# Determine ID to assign
if uid is None:
uid = max(mapping, default=0) + 1
else:
if uid in mapping:
raise AllocationError('A mesh with ID={} has already '
'been allocated.'.format(uid))
# Set the mesh type -- note that mesh type attribute only
# exists on subclasses!
index = c_int32()
_dll.openmc_extend_meshes(1, cls.mesh_type.encode(), index,
None)
index = index.value
else:
index = mapping[uid]._index
if index not in cls.__instances:
instance = super().__new__(cls)
instance._index = index
if uid is not None:
instance.id = uid
cls.__instances[index] = instance
return cls.__instances[index]
@property
def id(self):
mesh_id = c_int32()
_dll.openmc_mesh_get_id(self._index, mesh_id)
return mesh_id.value
@id.setter
def id(self, mesh_id):
_dll.openmc_mesh_set_id(self._index, mesh_id)
class RegularMesh(_FortranObjectWithID):
@ -71,57 +131,24 @@ class RegularMesh(_FortranObjectWithID):
The width of mesh cells in each direction.
"""
__instances = WeakValueDictionary()
mesh_type = 'regular'
def __new__(cls, uid=None, new=True, index=None):
mapping = meshes
if index is None:
if new:
# Determine ID to assign
if uid is None:
uid = max(mapping, default=0) + 1
else:
if uid in mapping:
raise AllocationError('A mesh with ID={} has already '
'been allocated.'.format(uid))
index = c_int32()
_dll.openmc_extend_meshes(1, index, None)
index = index.value
else:
index = mapping[uid]._index
if index not in cls.__instances:
instance = super().__new__(cls)
instance._index = index
if uid is not None:
instance.id = uid
cls.__instances[index] = instance
return cls.__instances[index]
@property
def id(self):
mesh_id = c_int32()
_dll.openmc_mesh_get_id(self._index, mesh_id)
return mesh_id.value
@id.setter
def id(self, mesh_id):
_dll.openmc_mesh_set_id(self._index, mesh_id)
def __init__(self, uid=None, new=True, index=None):
super().__init__(uid, new, index)
@property
def dimension(self):
dims = POINTER(c_int)()
n = c_int()
_dll.openmc_mesh_get_dimension(self._index, dims, n)
_dll.openmc_regular_mesh_get_dimension(self._index, dims, n)
return tuple(as_array(dims, (n.value,)))
@dimension.setter
def dimension(self, dimension):
n = len(dimension)
dimension = (c_int*n)(*dimension)
_dll.openmc_mesh_set_dimension(self._index, n, dimension)
_dll.openmc_regular_mesh_set_dimension(self._index, n, dimension)
@property
def lower_left(self):
@ -140,7 +167,7 @@ class RegularMesh(_FortranObjectWithID):
ur = POINTER(c_double)()
w = POINTER(c_double)()
n = c_int()
_dll.openmc_mesh_get_params(self._index, ll, ur, w, n)
_dll.openmc_regular_mesh_get_params(self._index, ll, ur, w, n)
return (
as_array(ll, (n.value,)),
as_array(ur, (n.value,)),
@ -157,7 +184,103 @@ class RegularMesh(_FortranObjectWithID):
if width is not None:
n = len(width)
width = (c_double*n)(*width)
_dll.openmc_mesh_set_params(self._index, n, lower_left, upper_right, width)
_dll.openmc_regular_mesh_set_params(self._index, n, lower_left, upper_right, width)
class RectilinearMesh(Mesh):
"""RectilinearMesh stored internally.
This class exposes a mesh that is stored internally in the OpenMC
library. To obtain a view of a mesh with a given ID, use the
:data:`openmc.lib.meshes` mapping.
Parameters
----------
index : int
Index in the `meshes` array.
Attributes
----------
id : int
ID of the mesh
dimension : iterable of int
The number of mesh cells in each direction.
lower_left : numpy.ndarray
The lower-left corner of the structured mesh.
upper_right : numpy.ndarray
The upper-right corner of the structrued mesh.
width : numpy.ndarray
The width of mesh cells in each direction.
"""
mesh_type = 'rectilinear'
def __init__(self, uid=None, new=True, index=None):
super().__init__(uid, new, index)
@property
def lower_left(self):
return self._get_parameters()[0]
@property
def upper_right(self):
return self._get_parameters()[1]
@property
def dimension(self):
return self._get_parameters()[2]
@property
def width(self):
return self._get_parameters()[3]
def _get_parameters(self):
gx = POINTER(c_double)()
nx = c_int()
gy = POINTER(c_double)()
ny = c_int()
gz = POINTER(c_double)()
nz = c_int()
# Call C API to get grid parameters
_dll.openmc_rectilinear_mesh_get_grid(self._index, gx, nx, gy, ny, gz,
nz)
# Convert grid parameters to Numpy arrays
grid_x = as_array(gx, (nx.value,))
grid_y = as_array(gy, (ny.value,))
grid_z = as_array(gz, (nz.value,))
# Calculate lower_left, upper_right, width, and dimension from grid
lower_left = np.array((grid_x[0], grid_y[0], grid_z[0]))
upper_right = np.array((grid_x[-1], grid_y[-1], grid_z[-1]))
dimension = np.array((nx.value - 1, ny.value - 1, nz.value - 1))
width = np.zeros(list(dimension) + [3])
for i, diff_x in enumerate(np.diff(grid_x)):
for j, diff_y in enumerate(np.diff(grid_y)):
for k, diff_z in enumerate(np.diff(grid_z)):
width[i, j, k, :] = diff_x, diff_y, diff_z
return (lower_left, upper_right, dimension, width)
def set_grid(self, x_grid, y_grid, z_grid):
nx = len(x_grid)
x_grid = (c_double*nx)(*x_grid)
ny = len(y_grid)
y_grid = (c_double*ny)(*y_grid)
nz = len(z_grid)
z_grid = (c_double*nz)(*z_grid)
_dll.openmc_rectilinear_mesh_set_grid(self._index, x_grid, nx, y_grid,
ny, z_grid, nz)
_MESH_TYPE_MAP = {
'regular': RegularMesh,
'rectilinear': RectilinearMesh
}
def _get_mesh(index):
mesh_type = create_string_buffer(20)
_dll.openmc_mesh_get_type(index, mesh_type)
mesh_type = mesh_type.value.decode()
return _MESH_TYPE_MAP[mesh_type](index=index)
class _MeshMapping(Mapping):
@ -168,11 +291,12 @@ class _MeshMapping(Mapping):
except (AllocationError, InvalidIDError) as e:
# __contains__ expects a KeyError to work correctly
raise KeyError(str(e))
return RegularMesh(index=index.value)
return _get_mesh(index.value)
def __iter__(self):
for i in range(len(self)):
yield RegularMesh(index=i).id
yield _get_mesh(i).id
def __len__(self):
return _dll.n_meshes()

259
src/mesh.cpp
View file

@ -162,6 +162,62 @@ int StructuredMesh::n_surface_bins() const
return 4 * n_dimension_ * n_bins();
}
xt::xtensor<double, 1>
StructuredMesh::count_sites(const Particle::Bank* bank,
int64_t length,
bool* outside) const
{
// Determine shape of array for counts
std::size_t m = this->n_bins();
std::vector<std::size_t> shape = {m};
// Create array of zeros
xt::xarray<double> cnt {shape, 0.0};
bool outside_ = false;
for (int64_t i = 0; i < length; i++) {
const auto& site = bank[i];
// determine scoring bin for entropy mesh
int mesh_bin = get_bin(site.r);
// if outside mesh, skip particle
if (mesh_bin < 0) {
outside_ = true;
continue;
}
// Add to appropriate bin
cnt(mesh_bin) += site.wgt;
}
// Create copy of count data. Since ownership will be acquired by xtensor,
// std::allocator must be used to avoid Valgrind mismatched free() / delete
// warnings.
int total = cnt.size();
double* cnt_reduced = std::allocator<double>{}.allocate(total);
#ifdef OPENMC_MPI
// collect values from all processors
MPI_Reduce(cnt.data(), cnt_reduced, total, MPI_DOUBLE, MPI_SUM, 0,
mpi::intracomm);
// Check if there were sites outside the mesh for any processor
if (outside) {
MPI_Reduce(&outside_, outside, 1, MPI_C_BOOL, MPI_LOR, 0, mpi::intracomm);
}
#else
std::copy(cnt.data(), cnt.data() + total, cnt_reduced);
if (outside) *outside = outside_;
#endif
// Adapt reduced values in array back into an xarray
auto arr = xt::adapt(cnt_reduced, total, xt::acquire_ownership(), shape);
xt::xarray<double> counts = arr;
return counts;
}
bool StructuredMesh::intersects(Position& r0, Position r1, int* ijk) const
{
switch(n_dimension_) {
@ -809,7 +865,7 @@ void RegularMesh::to_hdf5(hid_t group) const
{
hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_));
write_dataset(mesh_group, "type", "regular");
write_dataset(mesh_group, "type", type());
write_dataset(mesh_group, "dimension", shape_);
write_dataset(mesh_group, "lower_left", lower_left_);
write_dataset(mesh_group, "upper_right", upper_right_);
@ -818,62 +874,6 @@ void RegularMesh::to_hdf5(hid_t group) const
close_group(mesh_group);
}
xt::xtensor<double, 1>
RegularMesh::count_sites(const Particle::Bank* bank,
int64_t length,
bool* outside) const
{
// Determine shape of array for counts
std::size_t m = this->n_bins();
std::vector<std::size_t> shape = {m};
// Create array of zeros
xt::xarray<double> cnt {shape, 0.0};
bool outside_ = false;
for (int64_t i = 0; i < length; i++) {
const auto& site = bank[i];
// determine scoring bin for entropy mesh
int mesh_bin = get_bin(site.r);
// if outside mesh, skip particle
if (mesh_bin < 0) {
outside_ = true;
continue;
}
// Add to appropriate bin
cnt(mesh_bin) += site.wgt;
}
// Create copy of count data. Since ownership will be acquired by xtensor,
// std::allocator must be used to avoid Valgrind mismatched free() / delete
// warnings.
int total = cnt.size();
double* cnt_reduced = std::allocator<double>{}.allocate(total);
#ifdef OPENMC_MPI
// collect values from all processors
MPI_Reduce(cnt.data(), cnt_reduced, total, MPI_DOUBLE, MPI_SUM, 0,
mpi::intracomm);
// Check if there were sites outside the mesh for any processor
if (outside) {
MPI_Reduce(&outside_, outside, 1, MPI_C_BOOL, MPI_LOR, 0, mpi::intracomm);
}
#else
std::copy(cnt.data(), cnt.data() + total, cnt_reduced);
if (outside) *outside = outside_;
#endif
// Adapt reduced values in array back into an xarray
auto arr = xt::adapt(cnt_reduced, total, xt::acquire_ownership(), shape);
xt::xarray<double> counts = arr;
return counts;
}
//==============================================================================
// RectilinearMesh implementation
//==============================================================================
@ -883,7 +883,7 @@ RectilinearMesh::RectilinearMesh(pugi::xml_node node)
{
n_dimension_ = 3;
grid_.resize(3);
grid_.resize(n_dimension_);
grid_[0] = get_node_array<double>(node, "x_grid");
grid_[1] = get_node_array<double>(node, "y_grid");
grid_[2] = get_node_array<double>(node, "z_grid");
@ -1122,7 +1122,7 @@ void RectilinearMesh::to_hdf5(hid_t group) const
{
hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_));
write_dataset(mesh_group, "type", "rectilinear");
write_dataset(mesh_group, "type", type());
write_dataset(mesh_group, "x_grid", grid_[0]);
write_dataset(mesh_group, "y_grid", grid_[1]);
write_dataset(mesh_group, "z_grid", grid_[2]);
@ -1145,12 +1145,13 @@ check_mesh(int32_t index)
}
int
check_regular_mesh(int32_t index, RegularMesh** mesh)
check_mesh_type(int32_t index, const std::string& mesh_compare_type)
{
if (int err = check_mesh(index)) return err;
*mesh = dynamic_cast<RegularMesh*>(model::meshes[index].get());
if (!*mesh) {
set_errmsg("This function is only valid for regular meshes.");
StructuredMesh* mesh = dynamic_cast<StructuredMesh*>(model::meshes[index].get());
auto mesh_type = mesh->type();
if (mesh_compare_type != mesh_type) {
set_errmsg("This function is only valid for " + mesh_type + " meshes.");
return OPENMC_E_INVALID_TYPE;
}
return 0;
@ -1160,17 +1161,37 @@ check_regular_mesh(int32_t index, RegularMesh** mesh)
// C API functions
//==============================================================================
// Return the type of mesh as a C string
extern "C" int
openmc_mesh_get_type(int32_t index, char* type)
{
if (int err = check_mesh(index)) return err;
StructuredMesh* mesh = dynamic_cast<StructuredMesh*>(model::meshes[index].get());
std::strcpy(type, mesh->type().c_str());
return 0;
}
RegularMesh* get_regular_mesh(int32_t index) {
return dynamic_cast<RegularMesh*>(model::meshes[index].get());
}
//! Extend the meshes array by n elements
extern "C" int
openmc_extend_meshes(int32_t n, int32_t* index_start, int32_t* index_end)
openmc_extend_meshes(int32_t n, const char* type, int32_t* index_start,
int32_t* index_end)
{
if (index_start) *index_start = model::meshes.size();
for (int i = 0; i < n; ++i) {
model::meshes.push_back(std::make_unique<RegularMesh>());
if (std::strcmp(type, "regular") == 0)
{
model::meshes.push_back(std::make_unique<RegularMesh>());
} else if (std::strcmp(type, "rectilinear") == 0)
{
model::meshes.push_back(std::make_unique<RectilinearMesh>());
} else {
throw std::runtime_error{"Unknown mesh type: " + std::string(type)};
}
}
if (index_end) *index_end = model::meshes.size() - 1;
@ -1209,23 +1230,23 @@ openmc_mesh_set_id(int32_t index, int32_t id)
return 0;
}
//! Get the dimension of a mesh
//! Get the dimension of a regular mesh
extern "C" int
openmc_mesh_get_dimension(int32_t index, int** dims, int* n)
openmc_regular_mesh_get_dimension(int32_t index, int** dims, int* n)
{
RegularMesh* mesh;
if (int err = check_regular_mesh(index, &mesh)) return err;
if (int err = check_mesh_type(index, "regular")) return err;
RegularMesh* mesh = dynamic_cast<RegularMesh*>(model::meshes[index].get());
*dims = mesh->shape_.data();
*n = mesh->n_dimension_;
return 0;
}
//! Set the dimension of a mesh
//! Set the dimension of a regular mesh
extern "C" int
openmc_mesh_set_dimension(int32_t index, int n, const int* dims)
openmc_regular_mesh_set_dimension(int32_t index, int n, const int* dims)
{
RegularMesh* mesh;
if (int err = check_regular_mesh(index, &mesh)) return err;
if (int err = check_mesh_type(index, "regular")) return err;
RegularMesh* mesh = dynamic_cast<RegularMesh*>(model::meshes[index].get());
// Copy dimension
std::vector<std::size_t> shape = {static_cast<std::size_t>(n)};
@ -1234,12 +1255,13 @@ openmc_mesh_set_dimension(int32_t index, int n, const int* dims)
return 0;
}
//! Get the mesh parameters
//! Get the regular mesh parameters
extern "C" int
openmc_mesh_get_params(int32_t index, double** ll, double** ur, double** width, int* n)
openmc_regular_mesh_get_params(int32_t index, double** ll, double** ur,
double** width, int* n)
{
RegularMesh* m;
if (int err = check_regular_mesh(index, &m)) return err;
if (int err = check_mesh_type(index, "regular")) return err;
RegularMesh* m = dynamic_cast<RegularMesh*>(model::meshes[index].get());
if (m->lower_left_.dimension() == 0) {
set_errmsg("Mesh parameters have not been set.");
@ -1253,13 +1275,13 @@ openmc_mesh_get_params(int32_t index, double** ll, double** ur, double** width,
return 0;
}
//! Set the mesh parameters
//! Set the regular mesh parameters
extern "C" int
openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur,
const double* width)
openmc_regular_mesh_set_params(int32_t index, int n, const double* ll,
const double* ur, const double* width)
{
RegularMesh* m;
if (int err = check_regular_mesh(index, &m)) return err;
if (int err = check_mesh_type(index, "regular")) return err;
RegularMesh* m = dynamic_cast<RegularMesh*>(model::meshes[index].get());
std::vector<std::size_t> shape = {static_cast<std::size_t>(n)};
if (ll && ur) {
@ -1282,6 +1304,77 @@ openmc_mesh_set_params(int32_t index, int n, const double* ll, const double* ur,
return 0;
}
//! Get the rectilinear mesh grid
extern "C" int
openmc_rectilinear_mesh_get_grid(int32_t index, double** grid_x, int* nx,
double** grid_y, int * ny, double** grid_z, int* nz)
{
if (int err = check_mesh_type(index, "rectilinear")) return err;
RectilinearMesh* m = dynamic_cast<RectilinearMesh*>(model::meshes[index].get());
if (m->lower_left_.dimension() == 0) {
set_errmsg("Mesh parameters have not been set.");
return OPENMC_E_ALLOCATE;
}
*grid_x = m->grid_[0].data();
*nx = m->grid_[0].size();
*grid_y = m->grid_[1].data();
*ny = m->grid_[1].size();
*grid_z = m->grid_[2].data();
*nz = m->grid_[2].size();
return 0;
}
//! Set the regular mesh parameters
extern "C" int
openmc_rectilinear_mesh_set_grid(int32_t index, const double* grid_x,
const int nx, const double* grid_y, const int ny,
const double* grid_z, const int nz)
{
if (int err = check_mesh_type(index, "rectilinear")) return err;
RectilinearMesh* m = dynamic_cast<RectilinearMesh*>(model::meshes[index].get());
m->n_dimension_ = 3;
m->grid_.resize(m->n_dimension_);
for (int i = 0; i < nx; i++) {
m->grid_[0].push_back(grid_x[i]);
}
for (int i = 0; i < ny; i++) {
m->grid_[1].push_back(grid_y[i]);
}
for (int i = 0; i < nz; i++) {
m->grid_[2].push_back(grid_z[i]);
}
m->shape_ = {static_cast<int>(m->grid_[0].size()) - 1,
static_cast<int>(m->grid_[1].size()) - 1,
static_cast<int>(m->grid_[2].size()) - 1};
for (const auto& g : m->grid_) {
if (g.size() < 2) {
set_errmsg("x-, y-, and z- grids for rectilinear meshes "
"must each have at least 2 points");
return OPENMC_E_INVALID_ARGUMENT;
}
for (int i = 1; i < g.size(); ++i) {
if (g[i] <= g[i-1]) {
std::cout << g[i] << " " << g[i-1] << "\n";
set_errmsg("Values in for x-, y-, and z- grids for "
"rectilinear meshes must be sorted and unique.");
return OPENMC_E_INVALID_ARGUMENT;
}
}
}
m->lower_left_ = {m->grid_[0].front(), m->grid_[1].front(), m->grid_[2].front()};
m->upper_right_ = {m->grid_[0].back(), m->grid_[1].back(), m->grid_[2].back()};
return 0;
}
#ifdef DAGMC
UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node)
@ -1577,7 +1670,7 @@ UnstructuredMesh::to_hdf5(hid_t group) const
{
hid_t mesh_group = create_group(group, fmt::format("mesh {}", id_));
write_dataset(mesh_group, "type", "unstructured");
write_dataset(mesh_group, "type", type());
write_dataset(mesh_group, "filename", filename_);
// write volume and centroid of each tet