diff --git a/include/openmc/capi.h b/include/openmc/capi.h index 605d4040ba..6bde7b6709 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -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(); diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 5a8367b8bb..b2e35819c6 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -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& bins, std::vector& 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 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& 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 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& bins) const override; @@ -239,7 +244,6 @@ public: void to_hdf5(hid_t group) const override; -private: std::vector> 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& bins, std::vector& lengths) const override; diff --git a/openmc/cmfd.py b/openmc/cmfd.py index a99b606193..91c2646337 100644 --- a/openmc/cmfd.py +++ b/openmc/cmfd.py @@ -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() diff --git a/openmc/lib/filter.py b/openmc/lib/filter.py index 1d16b68e07..a47b374e23 100644 --- a/openmc/lib/filter.py +++ b/openmc/lib/filter.py @@ -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): diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py index e96115049a..92684ec8db 100644 --- a/openmc/lib/mesh.py +++ b/openmc/lib/mesh.py @@ -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() diff --git a/src/mesh.cpp b/src/mesh.cpp index a5409f05a0..576ae3e772 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -162,6 +162,62 @@ int StructuredMesh::n_surface_bins() const return 4 * n_dimension_ * n_bins(); } +xt::xtensor +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 shape = {m}; + + // Create array of zeros + xt::xarray 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{}.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 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 -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 shape = {m}; - - // Create array of zeros - xt::xarray 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{}.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 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(node, "x_grid"); grid_[1] = get_node_array(node, "y_grid"); grid_[2] = get_node_array(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(model::meshes[index].get()); - if (!*mesh) { - set_errmsg("This function is only valid for regular meshes."); + StructuredMesh* mesh = dynamic_cast(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(model::meshes[index].get()); + std::strcpy(type, mesh->type().c_str()); + return 0; +} + RegularMesh* get_regular_mesh(int32_t index) { return dynamic_cast(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()); + if (std::strcmp(type, "regular") == 0) + { + model::meshes.push_back(std::make_unique()); + } else if (std::strcmp(type, "rectilinear") == 0) + { + model::meshes.push_back(std::make_unique()); + } 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(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(model::meshes[index].get()); // Copy dimension std::vector shape = {static_cast(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(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(model::meshes[index].get()); std::vector shape = {static_cast(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(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(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(m->grid_[0].size()) - 1, + static_cast(m->grid_[1].size()) - 1, + static_cast(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