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

Merge pull request #1931 from ojschumann/mesh

Browse source 
Add Cylindrical and Spherical Mesh for MeshFilter and MeshSurfaceFilter
This commit is contained in:
Paul Romano 2022-01-12 21:12:10 -06:00 • committed by GitHub
commit 0cf7d92837
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
12 changed files with 1773 additions and 775 deletions
Download patch file Download diff file Expand all files Collapse all files

15
docs/source/io_formats/tallies.rst
View file

@ -312,8 +312,8 @@ element with the tag name ``<mesh>``. This element has the following
attributes/sub-elements:
:type:
The type of mesh. This can be either "regular", "rectilinear", or
"unstructured".
The type of mesh. This can be either "regular", "rectilinear",
"cylindrical", "spherical", or "unstructured".
:dimension:
The number of mesh cells in each direction. (For regular mesh only.)
@ -340,7 +340,16 @@ attributes/sub-elements:
The mesh divisions along the y-axis. (For rectilinear mesh only.)
:z_grid:
The mesh divisions along the z-axis. (For rectilinear mesh only.)
The mesh divisions along the z-axis. (For rectilinear and cylindrical meshes only.)
:r_grid:
The mesh divisions along the r-axis. (For cylindrical and spherical meshes only.)
:phi_grid:
The mesh divisions along the phi-axis. (For cylindrical and spherical meshes only.)
:theta_grid:
The mesh divisions along the theta-axis. (For spherical mesh only.)
:library:
The mesh library used to represent an unstructured mesh. This can be either

246
include/openmc/mesh.h
View file

@ -123,6 +123,9 @@ public:
//! \param[in] bin Mesh bin to generate a label for
virtual std::string bin_label(int bin) const = 0;
//! Return the mesh type
virtual std::string get_mesh_type() const = 0;
// Data members
int id_ {-1}; //!< User-specified ID
int n_dimension_; //!< Number of dimensions
@ -134,6 +137,22 @@ public:
StructuredMesh(pugi::xml_node node) : Mesh {node} {};
virtual ~StructuredMesh() = default;
using MeshIndex = std::array<int, 3>;
struct MeshDistance {
MeshDistance() = default;
MeshDistance(int _index, bool _max_surface, double _distance)
: next_index {_index}, max_surface {_max_surface}, distance {_distance}
{}
int next_index {-1};
double distance {INFTY};
bool max_surface {true};
bool operator<(const MeshDistance& o) const
{
return distance < o.distance;
}
};
int get_bin(Position r) const override;
int n_bins() const override;
@ -143,6 +162,19 @@ public:
void bins_crossed(Position r0, Position r1, const Direction& u,
vector<int>& bins, vector<double>& lengths) const override;
void surface_bins_crossed(Position r0, Position r1, const Direction& u,
vector<int>& bins) const override;
//! Determine which cell or surface bins were crossed by a particle
//
//! \param[in] r0 Previous position of the particle
//! \param[in] r1 Current position of the particle
//! \param[in] u Particle direction
//! \param[in] tally Functor that eventually stores the tally data
template<class T>
void raytrace_mesh(
Position r0, Position r1, const Direction& u, T tally) const;
//! Count number of bank sites in each mesh bin / energy bin
//
//! \param[in] Pointer to bank sites
@ -155,20 +187,20 @@ public:
//
//! \param[in] Array of mesh indices
//! \return Mesh bin
virtual int get_bin_from_indices(const int* ijk) const;
virtual int get_bin_from_indices(const MeshIndex& ijk) const;
//! Get mesh indices given a position
//
//! \param[in] r Position to get indices for
//! \param[out] ijk Array of mesh indices
//! \param[out] in_mesh Whether position is in mesh
virtual void get_indices(Position r, int* ijk, bool* in_mesh) const;
//! \return Array of mesh indices
virtual MeshIndex get_indices(Position r, bool& in_mesh) const;
//! Get mesh indices corresponding to a mesh bin
//
//! \param[in] bin Mesh bin
//! \param[out] ijk Mesh indices
virtual void get_indices_from_bin(int bin, int* ijk) const;
//! \return ijk Mesh indices
virtual MeshIndex get_indices_from_bin(int bin) const;
//! Get mesh index in a particular direction
//!
@ -176,38 +208,34 @@ public:
//! \param[in] i Direction index
virtual int get_index_in_direction(double r, int i) const = 0;
//! Check where a line segment intersects the mesh and if it intersects at all
//
//! \param[in,out] r0 In: starting position, out: intersection point
//! \param[in] r1 Ending position
//! \param[out] ijk Indices of the mesh bin containing the intersection point
//! \return Whether the line segment connecting r0 and r1 intersects mesh
virtual bool intersects(Position& r0, Position r1, int* ijk) const;
//! Get the coordinate for the mesh grid boundary in the positive direction
//! Get the closest distance from the coordinate r to the grid surface
//! in i direction that bounds mesh cell ijk and that is larger than l
//! The coordinate r does not have to be inside the mesh cell ijk. In
//! curved coordinates, multiple crossings of the same surface can happen,
//! these are selected by the parameter l
//!
//! \param[in] ijk Array of mesh indices
//! \param[in] i Direction index
virtual double positive_grid_boundary(int* ijk, int i) const = 0;
//! Get the coordinate for the mesh grid boundary in the negative direction
//!
//! \param[in] ijk Array of mesh indices
//! \param[in] i Direction index
virtual double negative_grid_boundary(int* ijk, int i) const = 0;
//! \param[in] i direction index of grid surface
//! \param[in] r0 position, from where to calculate the distance
//! \param[in] u direction of flight. actual position is r0 + l * u
//! \param[in] l actual chord length
//! \return MeshDistance struct with closest distance, next cell index in
//! i-direction and min/max surface indicator
virtual MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
const Position& r0, const Direction& u, double l) const = 0;
//! Get a label for the mesh bin
std::string bin_label(int bin) const override;
//! Get shape as xt::xtensor
xt::xtensor<int, 1> get_x_shape() const;
// Data members
xt::xtensor<double, 1> lower_left_; //!< Lower-left coordinates of mesh
xt::xtensor<double, 1> upper_right_; //!< Upper-right coordinates of mesh
xt::xtensor<int, 1> shape_; //!< Number of mesh elements in each dimension
std::array<int, 3> shape_; //!< Number of mesh elements in each dimension
protected:
virtual bool intersects_1d(Position& r0, Position r1, int* ijk) const;
virtual bool intersects_2d(Position& r0, Position r1, int* ijk) const;
virtual bool intersects_3d(Position& r0, Position r1, int* ijk) const;
};
//==============================================================================
@ -221,14 +249,14 @@ public:
RegularMesh(pugi::xml_node node);
// Overridden methods
void surface_bins_crossed(Position r0, Position r1, const Direction& u,
vector<int>& bins) const override;
int get_index_in_direction(double r, int i) const override;
double positive_grid_boundary(int* ijk, int i) const override;
virtual std::string get_mesh_type() const override;
double negative_grid_boundary(int* ijk, int i) const override;
static const std::string mesh_type;
MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
const Position& r0, const Direction& u, double l) const override;
std::pair<vector<double>, vector<double>> plot(
Position plot_ll, Position plot_ur) const override;
@ -236,6 +264,17 @@ public:
void to_hdf5(hid_t group) const override;
// New methods
//! Get the coordinate for the mesh grid boundary in the positive direction
//!
//! \param[in] ijk Array of mesh indices
//! \param[in] i Direction index
double positive_grid_boundary(const MeshIndex& ijk, int i) const;
//! Get the coordinate for the mesh grid boundary in the negative direction
//!
//! \param[in] ijk Array of mesh indices
//! \param[in] i Direction index
double negative_grid_boundary(const MeshIndex& ijk, int i) const;
//! Count number of bank sites in each mesh bin / energy bin
//
@ -257,25 +296,151 @@ public:
RectilinearMesh(pugi::xml_node node);
// Overridden methods
void surface_bins_crossed(Position r0, Position r1, const Direction& u,
vector<int>& bins) const override;
int get_index_in_direction(double r, int i) const override;
double positive_grid_boundary(int* ijk, int i) const override;
virtual std::string get_mesh_type() const override;
double negative_grid_boundary(int* ijk, int i) const override;
static const std::string mesh_type;
MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
const Position& r0, const Direction& u, double l) const override;
std::pair<vector<double>, vector<double>> plot(
Position plot_ll, Position plot_ur) const override;
void to_hdf5(hid_t group) const override;
vector<vector<double>> grid_;
// New methods
//! Get the coordinate for the mesh grid boundary in the positive direction
//!
//! \param[in] ijk Array of mesh indices
//! \param[in] i Direction index
double positive_grid_boundary(const MeshIndex& ijk, int i) const;
//! Get the coordinate for the mesh grid boundary in the negative direction
//!
//! \param[in] ijk Array of mesh indices
//! \param[in] i Direction index
double negative_grid_boundary(const MeshIndex& ijk, int i) const;
array<vector<double>, 3> grid_;
int set_grid();
};
class CylindricalMesh : public StructuredMesh {
public:
// Constructors
CylindricalMesh() = default;
CylindricalMesh(pugi::xml_node node);
// Overridden methods
virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
int get_index_in_direction(double r, int i) const override;
virtual std::string get_mesh_type() const override;
static const std::string mesh_type;
MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
const Position& r0, const Direction& u, double l) const override;
std::pair<vector<double>, vector<double>> plot(
Position plot_ll, Position plot_ur) const override;
void to_hdf5(hid_t group) const override;
array<vector<double>, 3> grid_;
int set_grid();
private:
double find_r_crossing(
const Position& r, const Direction& u, double l, int shell) const;
double find_phi_crossing(
const Position& r, const Direction& u, double l, int shell) const;
StructuredMesh::MeshDistance find_z_crossing(
const Position& r, const Direction& u, double l, int shell) const;
bool full_phi_ {false};
constexpr inline int sanitize_angular_index(int idx, bool full, int N) const
{
if ((idx > 0) and (idx <= N)) {
return idx;
} else if (full) {
return (idx + N - 1) % N + 1;
} else {
return 0;
}
}
inline int sanitize_phi(int idx) const
{
return sanitize_angular_index(idx, full_phi_, shape_[1]);
}
};
class SphericalMesh : public StructuredMesh {
public:
// Constructors
SphericalMesh() = default;
SphericalMesh(pugi::xml_node node);
// Overridden methods
virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
int get_index_in_direction(double r, int i) const override;
virtual std::string get_mesh_type() const override;
static const std::string mesh_type;
MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
const Position& r0, const Direction& u, double l) const override;
std::pair<vector<double>, vector<double>> plot(
Position plot_ll, Position plot_ur) const override;
void to_hdf5(hid_t group) const override;
array<vector<double>, 3> grid_;
int set_grid();
private:
double find_r_crossing(
const Position& r, const Direction& u, double l, int shell) const;
double find_theta_crossing(
const Position& r, const Direction& u, double l, int shell) const;
double find_phi_crossing(
const Position& r, const Direction& u, double l, int shell) const;
bool full_theta_ {false};
bool full_phi_ {false};
constexpr inline int sanitize_angular_index(int idx, bool full, int N) const
{
if ((idx > 0) and (idx <= N)) {
return idx;
} else if (full) {
return (idx + N - 1) % N + 1;
} else {
return 0;
}
}
inline int sanitize_theta(int idx) const
{
return sanitize_angular_index(idx, full_theta_, shape_[1]);
}
inline int sanitize_phi(int idx) const
{
return sanitize_angular_index(idx, full_phi_, shape_[2]);
}
};
// Abstract class for unstructured meshes
class UnstructuredMesh : public Mesh {
@ -285,6 +450,9 @@ public:
UnstructuredMesh(pugi::xml_node node);
UnstructuredMesh(const std::string& filename);
static const std::string mesh_type;
virtual std::string get_mesh_type() const override;
// Overridden Methods
void surface_bins_crossed(Position r0, Position r1, const Direction& u,
vector<int>& bins) const override;
@ -356,6 +524,8 @@ public:
MOABMesh(const std::string& filename, double length_multiplier = 1.0);
MOABMesh(std::shared_ptr<moab::Interface> external_mbi);
static const std::string mesh_lib_type;
// Overridden Methods
void bins_crossed(Position r0, Position r1, const Direction& u,
@ -503,6 +673,8 @@ public:
LibMesh(pugi::xml_node node);
LibMesh(const std::string& filename, double length_multiplier = 1.0);
static const std::string mesh_lib_type;
// Overridden Methods
void bins_crossed(Position r0, Position r1, const Direction& u,
vector<int>& bins, vector<double>& lengths) const override;

199
openmc/lib/mesh.py
View file

@ -11,7 +11,7 @@ from . import _dll
from .core import _FortranObjectWithID
from .error import _error_handler
__all__ = ['RegularMesh', 'RectilinearMesh', 'meshes']
__all__ = ['RegularMesh', 'RectilinearMesh', 'CylindricalMesh', 'SphericalMesh', 'meshes']
# Mesh functions
_dll.openmc_extend_meshes.argtypes = [c_int32, c_char_p, POINTER(c_int32),
@ -56,6 +56,26 @@ _dll.openmc_regular_mesh_set_params.argtypes = [
_dll.openmc_regular_mesh_set_params.restype = c_int
_dll.openmc_regular_mesh_set_params.errcheck = _error_handler
_dll.openmc_cylindrical_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_cylindrical_mesh_get_grid.restype = c_int
_dll.openmc_cylindrical_mesh_get_grid.errcheck = _error_handler
_dll.openmc_cylindrical_mesh_set_grid.argtypes = [c_int32, POINTER(c_double),
c_int, POINTER(c_double), c_int, POINTER(c_double), c_int]
_dll.openmc_cylindrical_mesh_set_grid.restype = c_int
_dll.openmc_cylindrical_mesh_set_grid.errcheck = _error_handler
_dll.openmc_spherical_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_spherical_mesh_get_grid.restype = c_int
_dll.openmc_spherical_mesh_get_grid.errcheck = _error_handler
_dll.openmc_spherical_mesh_set_grid.argtypes = [c_int32, POINTER(c_double),
c_int, POINTER(c_double), c_int, POINTER(c_double), c_int]
_dll.openmc_spherical_mesh_set_grid.restype = c_int
_dll.openmc_spherical_mesh_set_grid.errcheck = _error_handler
class Mesh(_FortranObjectWithID):
"""Base class to represent mesh objects
@ -275,9 +295,184 @@ class RectilinearMesh(Mesh):
ny, z_grid, nz)
class CylindricalMesh(Mesh):
"""CylindricalMesh 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 = 'cylindrical'
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_cylindrical_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_cylindrical_mesh_set_grid(self._index, x_grid, nx, y_grid,
ny, z_grid, nz)
class SphericalMesh(Mesh):
"""SphericalMesh 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 = 'spherical'
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_spherical_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_spherical_mesh_set_grid(self._index, x_grid, nx, y_grid,
ny, z_grid, nz)
_MESH_TYPE_MAP = {
'regular': RegularMesh,
'rectilinear': RectilinearMesh
'rectilinear': RectilinearMesh,
'cylindrical': CylindricalMesh,
'spherical': SphericalMesh
}

332
openmc/mesh.py
View file

@ -79,6 +79,10 @@ class MeshBase(IDManagerMixin, ABC):
return RegularMesh.from_hdf5(group)
elif mesh_type == 'rectilinear':
return RectilinearMesh.from_hdf5(group)
elif mesh_type == 'cylindrical':
return CylindricalMesh.from_hdf5(group)
elif mesh_type == 'spherical':
return SphericalMesh.from_hdf5(group)
elif mesh_type == 'unstructured':
return UnstructuredMesh.from_hdf5(group)
else:
@ -624,6 +628,334 @@ class RectilinearMesh(MeshBase):
return element
class CylindricalMesh(MeshBase):
"""A 3D cylindrical mesh
Parameters
----------
mesh_id : int
Unique identifier for the mesh
name : str
Name of the mesh
Attributes
----------
id : int
Unique identifier for the mesh
name : str
Name of the mesh
dimension : Iterable of int
The number of mesh cells in each direction.
n_dimension : int
Number of mesh dimensions (always 3 for a CylindricalMesh).
r_grid : Iterable of float
Mesh boundary points along the r-axis.
Requirement is r >= 0.
phi_grid : Iterable of float
Mesh boundary points along the phi-axis.
The default value is [0, 360], i.e. the full phi range.
z_grid : Iterable of float
Mesh boundary points along the z-axis.
indices : Iterable of tuple
An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1),
(2, 1, 1), ...]
"""
def __init__(self, mesh_id=None, name=''):
super().__init__(mesh_id, name)
self._r_grid = None
self._phi_grid = [0, 360]
self._z_grid = None
@property
def dimension(self):
return (len(self.r_grid) - 1,
len(self.phi_grid) - 1,
len(self.z_grid) - 1)
@property
def n_dimension(self):
return 3
@property
def r_grid(self):
return self._r_grid
@property
def phi_grid(self):
return self._phi_grid
@property
def z_grid(self):
return self._z_grid
@property
def indices(self):
nr, np, nz = self.dimension
np = len(self.phi_grid) - 1
nz = len(self.z_grid) - 1
return ((r, p, z)
for z in range(1, nz + 1)
for p in range(1, np + 1)
for r in range(1, nr + 1))
@r_grid.setter
def r_grid(self, grid):
cv.check_type('mesh r_grid', grid, Iterable, Real)
self._r_grid = grid
@phi_grid.setter
def phi_grid(self, grid):
cv.check_type('mesh phi_grid', grid, Iterable, Real)
self._phi_grid = np.array(grid)
@z_grid.setter
def z_grid(self, grid):
cv.check_type('mesh z_grid', grid, Iterable, Real)
self._z_grid = np.array(grid)
def __repr__(self):
fmt = '{0: <16}{1}{2}\n'
string = super().__repr__()
string += fmt.format('\tDimensions', '=\t', self.n_dimension)
r_grid_str = str(self._r_grid) if self._r_grid is None else len(self._r_grid)
string += fmt.format('\tN R pnts:', '=\t', r_grid_str)
if self._r_grid is not None:
string += fmt.format('\tR Min:', '=\t', self._r_grid[0])
string += fmt.format('\tR Max:', '=\t', self._r_grid[-1])
phi_grid_str = str(self._phi_grid) if self._phi_grid is None else len(self._phi_grid)
string += fmt.format('\tN Phi pnts:', '=\t', phi_grid_str)
if self._phi_grid is not None:
string += fmt.format('\tPhi Min:', '=\t', self._phi_grid[0])
string += fmt.format('\tPhi Max:', '=\t', self._phi_grid[-1])
z_grid_str = str(self._z_grid) if self._z_grid is None else len(self._z_grid)
string += fmt.format('\tN Z pnts:', '=\t', z_grid_str)
if self._z_grid is not None:
string += fmt.format('\tZ Min:', '=\t', self._z_grid[0])
string += fmt.format('\tZ Max:', '=\t', self._z_grid[-1])
return string
@classmethod
def from_hdf5(cls, group):
mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
# Read and assign mesh properties
mesh = cls(mesh_id)
mesh.r_grid = group['r_grid'][()]
mesh.phi_grid = 180 / np.pi * group['phi_grid'][()]
mesh.z_grid = group['z_grid'][()]
return mesh
def to_xml_element(self):
"""Return XML representation of the mesh
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing mesh data
"""
element = ET.Element("mesh")
element.set("id", str(self._id))
element.set("type", "cylindrical")
subelement = ET.SubElement(element, "r_grid")
subelement.text = ' '.join(map(str, self.r_grid))
subelement = ET.SubElement(element, "phi_grid")
subelement.text = ' '.join(map(str, self.phi_grid))
subelement = ET.SubElement(element, "z_grid")
subelement.text = ' '.join(map(str, self.z_grid))
return element
def calc_mesh_volumes(self):
"""Return Volumes for every mesh cell
Returns
-------
volumes : Iterable of float
Volumes
"""
V_r = np.diff(np.array(self.r_grid)**2 / 2)
V_p = np.diff(np.array(self.phi_grid) * np.pi / 180.0)
V_z = np.diff(np.array(self.z_grid))
return np.multiply.outer(np.outer(V_r, V_p), V_z)
class SphericalMesh(MeshBase):
"""A 3D spherical mesh
Parameters
----------
mesh_id : int
Unique identifier for the mesh
name : str
Name of the mesh
Attributes
----------
id : int
Unique identifier for the mesh
name : str
Name of the mesh
dimension : Iterable of int
The number of mesh cells in each direction.
n_dimension : int
Number of mesh dimensions (always 3 for a SphericalMesh).
r_grid : Iterable of float
Mesh boundary points along the r-axis.
Requirement is r >= 0.
theta_grid : Iterable of float
Mesh boundary points along the theta-axis in degrees.
The default value is [0, 180], i.e. the full theta range.
phi_grid : Iterable of float
Mesh boundary points along the phi-axis in degrees.
The default value is [0, 360], i.e. the full phi range.
indices : Iterable of tuple
An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1),
(2, 1, 1), ...]
"""
def __init__(self, mesh_id=None, name=''):
super().__init__(mesh_id, name)
self._r_grid = None
self._theta_grid = [0, 180]
self._phi_grid = [0, 360]
@property
def dimension(self):
return (len(self.r_grid) - 1,
len(self.theta_grid) - 1,
len(self.phi_grid) - 1)
@property
def n_dimension(self):
return 3
@property
def r_grid(self):
return self._r_grid
@property
def theta_grid(self):
return self._theta_grid
@property
def phi_grid(self):
return self._phi_grid
@property
def indices(self):
nr, nt, np = self.dimension
nt = len(self.theta_grid) - 1
np = len(self.phi_grid) - 1
return ((r, t, p)
for p in range(1, np + 1)
for t in range(1, nt + 1)
for r in range(1, nr + 1))
@r_grid.setter
def r_grid(self, grid):
cv.check_type('mesh r_grid', grid, Iterable, Real)
self._r_grid = grid
@theta_grid.setter
def theta_grid(self, grid):
cv.check_type('mesh theta_grid', grid, Iterable, Real)
self._theta_grid = np.array(grid)
@phi_grid.setter
def phi_grid(self, grid):
cv.check_type('mesh phi_grid', grid, Iterable, Real)
self._phi_grid = np.array(grid)
def __repr__(self):
fmt = '{0: <16}{1}{2}\n'
string = super().__repr__()
string += fmt.format('\tDimensions', '=\t', self.n_dimension)
r_grid_str = str(self._r_grid) if not self._r_grid else len(self._r_grid)
string += fmt.format('\tN R pnts:', '=\t', r_grid_str)
if self._r_grid:
string += fmt.format('\tR Min:', '=\t', self._r_grid[0])
string += fmt.format('\tR Max:', '=\t', self._r_grid[-1])
theta_grid_str = str(self._theta_grid) if not self._theta_grid else len(self._theta_grid)
string += fmt.format('\tN Theta pnts:', '=\t', theta_grid_str)
if self._theta_grid:
string += fmt.format('\tTheta Min:', '=\t', self._theta_grid[0])
string += fmt.format('\tTheta Max:', '=\t', self._theta_grid[-1])
phi_grid_str = str(self._phi_grid) if not self._phi_grid else len(self._phi_grid)
string += fmt.format('\tN Phi pnts:', '=\t', phi_grid_str)
if self._phi_grid:
string += fmt.format('\tPhi Min:', '=\t', self._phi_grid[0])
string += fmt.format('\tPhi Max:', '=\t', self._phi_grid[-1])
return string
@classmethod
def from_hdf5(cls, group):
mesh_id = int(group.name.split('/')[-1].lstrip('mesh '))
# Read and assign mesh properties
mesh = cls(mesh_id)
mesh.r_grid = group['r_grid'][()]
mesh.theta_grid = 180 / np.pi * group['theta_grid'][()]
mesh.phi_grid = 180 / np.pi * group['phi_grid'][()]
return mesh
def to_xml_element(self):
"""Return XML representation of the mesh
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing mesh data
"""
element = ET.Element("mesh")
element.set("id", str(self._id))
element.set("type", "spherical")
subelement = ET.SubElement(element, "r_grid")
subelement.text = ' '.join(map(str, self.r_grid))
subelement = ET.SubElement(element, "theta_grid")
subelement.text = ' '.join(map(str, self.theta_grid))
subelement = ET.SubElement(element, "phi_grid")
subelement.text = ' '.join(map(str, self.phi_grid))
return element
def calc_mesh_volumes(self):
"""Return Volumes for every mesh cell
Returns
-------
volumes : Iterable of float
Volumes
"""
V_r = np.diff(np.array(self.r_grid)**3 / 3)
V_t = np.diff(-np.cos(np.pi * np.array(self.theta_grid) / 180.0))
V_p = np.diff(np.array(self.phi_grid) * np.pi / 180)
return np.multiply.outer(np.outer(V_r, V_t), V_p)
class UnstructuredMesh(MeshBase):
"""A 3D unstructured mesh

2
src/eigenvalue.cpp
View file

@ -577,7 +577,7 @@ void ufs_count_sites()
#ifdef OPENMC_MPI
// Send source fraction to all processors
int n_bins = xt::prod(simulation::ufs_mesh->shape_)();
int n_bins = simulation::ufs_mesh->n_bins();
MPI_Bcast(
simulation::source_frac.data(), n_bins, MPI_DOUBLE, 0, mpi::intracomm);
#endif

1584
src/mesh.cpp
View file

File diff suppressed because it is too large Load diff

18
src/relaxng/tallies.rnc
View file

@ -24,6 +24,24 @@ element tallies {
attribute y_grid { list { xsd:double+ } }) &
(element z_grid { list { xsd:double+ } } |
attribute z_grid { list { xsd:double+ } })
) | (
(element type { ( "cylindrical" ) } |
attribute type { ( "cylindrical" ) }) &
(element r_grid { list { xsd:double+ } } |
attribute r_grid { list { xsd:double+ } }) &
(element phi_grid { list { xsd:double+ } } |
attribute phi_grid { list { xsd:double+ } }) &
(element z_grid { list { xsd:double+ } } |
attribute z_grid { list { xsd:double+ } })
) | (
(element type { ( "spherical" ) } |
attribute type { ( "spherical" ) }) &
(element r_grid { list { xsd:double+ } } |
attribute r_grid { list { xsd:double+ } }) &
(element theta_grid { list { xsd:double+ } } |
attribute theta_grid { list { xsd:double+ } }) &
(element phi_grid { list { xsd:double+ } } |
attribute phi_grid { list { xsd:double+ } })
)
)
}* &

54
tests/regression_tests/filter_mesh/inputs_true.dat
View file

@ -53,36 +53,58 @@
<y_grid>-7.5 -6.617647058823529 -5.735294117647059 -4.852941176470589 -3.9705882352941178 -3.0882352941176467 -2.2058823529411766 -1.3235294117647065 -0.4411764705882355 0.4411764705882355 1.3235294117647065 2.2058823529411757 3.0882352941176467 3.9705882352941178 4.852941176470587 5.735294117647058 6.617647058823529 7.5</y_grid>
<z_grid>1.0 1.223224374241637 1.4962778697388448 1.8302835609029084 2.2388474634702153 2.7386127875258306 3.3499379133114306 4.09772570775871 5.012437964687018 6.131336292779302 7.500000000000001</z_grid>
</mesh>
<mesh id="5" type="cylindrical">
<r_grid>0.0 0.4411764705882353 0.8823529411764706 1.3235294117647058 1.7647058823529411 2.2058823529411766 2.6470588235294117 3.0882352941176467 3.5294117647058822 3.9705882352941178 4.411764705882353 4.852941176470588 5.294117647058823 5.735294117647059 6.1764705882352935 6.617647058823529 7.0588235294117645 7.5</r_grid>
<phi_grid>0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0 180.0 200.0 220.0 240.0 260.0 280.0 300.0 320.0 340.0 360.0</phi_grid>
<z_grid>-7.5 -6.5625 -5.625 -4.6875 -3.75 -2.8125 -1.875 -0.9375 0.0 0.9375 1.875 2.8125 3.75 4.6875 5.625 6.5625 7.5</z_grid>
</mesh>
<mesh id="6" type="spherical">
<r_grid>0.0 0.4411764705882353 0.8823529411764706 1.3235294117647058 1.7647058823529411 2.2058823529411766 2.6470588235294117 3.0882352941176467 3.5294117647058822 3.9705882352941178 4.411764705882353 4.852941176470588 5.294117647058823 5.735294117647059 6.1764705882352935 6.617647058823529 7.0588235294117645 7.5</r_grid>
<theta_grid>0.0 22.5 45.0 67.5 90.0 112.5 135.0 157.5 180.0</theta_grid>
<phi_grid>0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0 180.0 200.0 220.0 240.0 260.0 280.0 300.0 320.0 340.0 360.0</phi_grid>
</mesh>
<filter id="1" type="mesh">
<bins>1</bins>
</filter>
<filter id="5" type="meshsurface">
<filter id="7" type="meshsurface">
<bins>1</bins>
</filter>
<filter id="2" type="mesh">
<bins>2</bins>
</filter>
<filter id="6" type="meshsurface">
<filter id="8" type="meshsurface">
<bins>2</bins>
</filter>
<filter id="3" type="mesh">
<bins>3</bins>
</filter>
<filter id="7" type="meshsurface">
<filter id="9" type="meshsurface">
<bins>3</bins>
</filter>
<filter id="4" type="mesh">
<bins>4</bins>
</filter>
<filter id="8" type="meshsurface">
<filter id="10" type="meshsurface">
<bins>4</bins>
</filter>
<filter id="5" type="mesh">
<bins>5</bins>
</filter>
<filter id="11" type="meshsurface">
<bins>5</bins>
</filter>
<filter id="6" type="mesh">
<bins>6</bins>
</filter>
<filter id="12" type="meshsurface">
<bins>6</bins>
</filter>
<tally id="1">
<filters>1</filters>
<scores>total</scores>
</tally>
<tally id="2">
<filters>5</filters>
<filters>7</filters>
<scores>current</scores>
</tally>
<tally id="3">
@ -90,7 +112,7 @@
<scores>total</scores>
</tally>
<tally id="4">
<filters>6</filters>
<filters>8</filters>
<scores>current</scores>
</tally>
<tally id="5">
@ -98,7 +120,7 @@
<scores>total</scores>
</tally>
<tally id="6">
<filters>7</filters>
<filters>9</filters>
<scores>current</scores>
</tally>
<tally id="7">
@ -106,7 +128,23 @@
<scores>total</scores>
</tally>
<tally id="8">
<filters>8</filters>
<filters>10</filters>
<scores>current</scores>
</tally>
<tally id="9">
<filters>5</filters>
<scores>total</scores>
</tally>
<tally id="10">
<filters>11</filters>
<scores>current</scores>
</tally>
<tally id="11">
<filters>6</filters>
<scores>total</scores>
</tally>
<tally id="12">
<filters>12</filters>
<scores>current</scores>
</tally>
</tallies>

2
tests/regression_tests/filter_mesh/results_true.dat
View file

@ -1 +1 @@
97cf099a24099d2af794752e145841a9ef9f24980c53ee71e3fdd400d4249a7e6f3f52f63132fe2fe8d1973604a72d87910ab13bc3721f7d0bec533b87756847
6e02d01fc115c54e4c60477a9f963149bc83c82e6a7bdb2d29125b131546150c59de5ac3aae5b9469b89b343ecf2574f6c949b918ccb43231e8666ee2b7c1b7c

18
tests/regression_tests/filter_mesh/test.py
View file

@ -51,18 +51,32 @@ def model():
recti_mesh.y_grid = np.linspace(-7.5, 7.5, 18)
recti_mesh.z_grid = np.logspace(0, np.log10(7.5), 11)
cyl_mesh = openmc.CylindricalMesh()
cyl_mesh.r_grid = np.linspace(0, 7.5, 18)
cyl_mesh.phi_grid = np.linspace(0, 360, 19)
cyl_mesh.z_grid = np.linspace(-7.5, 7.5, 17)
sph_mesh = openmc.SphericalMesh()
sph_mesh.r_grid = np.linspace(0, 7.5, 18)
sph_mesh.theta_grid = np.linspace(0, 180, 9)
sph_mesh.phi_grid = np.linspace(0, 360, 19)
# Create filters
reg_filters = [
openmc.MeshFilter(mesh_1d),
openmc.MeshFilter(mesh_2d),
openmc.MeshFilter(mesh_3d),
openmc.MeshFilter(recti_mesh)
openmc.MeshFilter(recti_mesh),
openmc.MeshFilter(cyl_mesh),
openmc.MeshFilter(sph_mesh)
]
surf_filters = [
openmc.MeshSurfaceFilter(mesh_1d),
openmc.MeshSurfaceFilter(mesh_2d),
openmc.MeshSurfaceFilter(mesh_3d),
openmc.MeshSurfaceFilter(recti_mesh)
openmc.MeshSurfaceFilter(recti_mesh),
openmc.MeshSurfaceFilter(cyl_mesh),
openmc.MeshSurfaceFilter(sph_mesh)
]
# Create tallies

16
tests/regression_tests/mgxs_library_mesh/results_true.dat
View file

@ -181,11 +181,11 @@
mesh 1 group in nuclide mean std. dev.
x y surf
3 1 1 x-max in 1 total 0.1816 0.006129
2 1 1 x-max out 1 total 0.2718 0.095077
2 1 1 x-max out 1 total 0.1730 0.004382
1 1 1 x-min in 1 total 0.0000 0.000000
0 1 1 x-min out 1 total 0.0000 0.000000
7 1 1 y-max in 1 total 0.1790 0.011921
6 1 1 y-max out 1 total 0.2290 0.034821
6 1 1 y-max out 1 total 0.1802 0.016599
5 1 1 y-min in 1 total 0.0000 0.000000
4 1 1 y-min out 1 total 0.0000 0.000000
19 1 2 x-max in 1 total 0.1872 0.012447
@ -194,23 +194,23 @@
16 1 2 x-min out 1 total 0.0000 0.000000
23 1 2 y-max in 1 total 0.0000 0.000000
22 1 2 y-max out 1 total 0.0000 0.000000
21 1 2 y-min in 1 total 0.2290 0.034821
21 1 2 y-min in 1 total 0.1802 0.016599
20 1 2 y-min out 1 total 0.1790 0.011921
11 2 1 x-max in 1 total 0.0000 0.000000
10 2 1 x-max out 1 total 0.0000 0.000000
9 2 1 x-min in 1 total 0.2718 0.095077
9 2 1 x-min in 1 total 0.1730 0.004382
8 2 1 x-min out 1 total 0.1816 0.006129
15 2 1 y-max in 1 total 0.1778 0.009484
14 2 1 y-max out 1 total 0.2326 0.042782
14 2 1 y-max out 1 total 0.1822 0.010077
13 2 1 y-min in 1 total 0.0000 0.000000
12 2 1 y-min out 1 total 0.0000 0.000000
27 2 2 x-max in 1 total 0.0000 0.000000
26 2 2 x-max out 1 total 0.0260 0.026000
26 2 2 x-max out 1 total 0.0000 0.000000
25 2 2 x-min in 1 total 0.1952 0.015948
24 2 2 x-min out 1 total 0.1872 0.012447
31 2 2 y-max in 1 total 0.0000 0.000000
30 2 2 y-max out 1 total 0.0244 0.024400
29 2 2 y-min in 1 total 0.2326 0.042782
30 2 2 y-max out 1 total 0.0000 0.000000
29 2 2 y-min in 1 total 0.1822 0.010077
28 2 2 y-min out 1 total 0.1778 0.009484
mesh 1 group in nuclide mean std. dev.
x y z

62
tests/unit_tests/test_lib.py
View file

@ -580,6 +580,68 @@ def test_rectilinear_mesh(lib_init):
msf = openmc.lib.MeshSurfaceFilter(mesh)
assert msf.mesh == mesh
def test_cylindrical_mesh(lib_init):
deg2rad = lambda deg: deg*np.pi/180
mesh = openmc.lib.CylindricalMesh()
x_grid = [0., 5., 10.]
y_grid = [0., 10., 20.]
z_grid = [10., 20., 30.]
mesh.set_grid(x_grid, y_grid, z_grid)
assert np.all(mesh.lower_left == (0., 0., 10.))
assert np.all(mesh.upper_right == (10., deg2rad(20.), 30.))
assert np.all(mesh.dimension == (2, 2, 2))
for i, diff_x in enumerate(np.diff(x_grid)):
for j, diff_y in enumerate(np.diff(y_grid)):
for k, diff_z in enumerate(np.diff(z_grid)):
assert np.all(mesh.width[i, j, k, :] == (5, deg2rad(10), 10))
with pytest.raises(exc.AllocationError):
mesh2 = openmc.lib.CylindricalMesh(mesh.id)
meshes = openmc.lib.meshes
assert isinstance(meshes, Mapping)
assert len(meshes) == 3
mesh = meshes[mesh.id]
assert isinstance(mesh, openmc.lib.CylindricalMesh)
mf = openmc.lib.MeshFilter(mesh)
assert mf.mesh == mesh
msf = openmc.lib.MeshSurfaceFilter(mesh)
assert msf.mesh == mesh
def test_spherical_mesh(lib_init):
deg2rad = lambda deg: deg*np.pi/180
mesh = openmc.lib.SphericalMesh()
x_grid = [0., 5., 10.]
y_grid = [0., 10., 20.]
z_grid = [10., 20., 30.]
mesh.set_grid(x_grid, y_grid, z_grid)
assert np.all(mesh.lower_left == (0., 0., deg2rad(10.)))
assert np.all(mesh.upper_right == (10., deg2rad(20.), deg2rad(30.)))
assert np.all(mesh.dimension == (2, 2, 2))
for i, diff_x in enumerate(np.diff(x_grid)):
for j, diff_y in enumerate(np.diff(y_grid)):
for k, diff_z in enumerate(np.diff(z_grid)):
assert np.all(abs(mesh.width[i, j, k, :] - (5, deg2rad(10), deg2rad(10))) < 1e-16)
with pytest.raises(exc.AllocationError):
mesh2 = openmc.lib.SphericalMesh(mesh.id)
meshes = openmc.lib.meshes
assert isinstance(meshes, Mapping)
assert len(meshes) == 4
mesh = meshes[mesh.id]
assert isinstance(mesh, openmc.lib.SphericalMesh)
mf = openmc.lib.MeshFilter(mesh)
assert mf.mesh == mesh
msf = openmc.lib.MeshSurfaceFilter(mesh)
assert msf.mesh == mesh
def test_restart(lib_init, mpi_intracomm):
# Finalize and re-init to make internal state consistent with XML.