mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-21 14:35:27 -04:00
Merge a83e50173b into f1fb6721f0
This commit is contained in:
GuySten
GitHub
commit
696b5d4bcf
10 changed files with 1181 additions and 188 deletions
|
|
@ -362,7 +362,7 @@ attributes/sub-elements:
|
|||
|
||||
:type:
|
||||
The type of mesh. This can be either "regular", "rectilinear",
|
||||
"cylindrical", "spherical", or "unstructured".
|
||||
"cylindrical", "spherical", "hexagonal", or "unstructured".
|
||||
|
||||
:dimension:
|
||||
The number of mesh cells in each direction. (For regular mesh only.)
|
||||
|
|
@ -379,6 +379,9 @@ attributes/sub-elements:
|
|||
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. (For regular mesh only.)
|
||||
|
||||
:pitch:
|
||||
The mesh radial pitch. (For hexagonal mesh only.)
|
||||
|
||||
:width:
|
||||
The width of mesh cells in each direction. (For regular mesh only.)
|
||||
|
||||
|
|
@ -389,7 +392,7 @@ 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 and cylindrical meshes only.)
|
||||
The mesh divisions along the z-axis. (For rectilinear, cylindrical and hexagonal meshes only.)
|
||||
|
||||
:r_grid:
|
||||
The mesh divisions along the r-axis. (For cylindrical and spherical meshes only.)
|
||||
|
|
@ -400,8 +403,14 @@ attributes/sub-elements:
|
|||
:theta_grid:
|
||||
The mesh divisions along the theta-axis. (For spherical mesh only.)
|
||||
|
||||
:num_rings:
|
||||
The number of hexagonal rings. (For hexagonal mesh only.)
|
||||
|
||||
:orientation:
|
||||
The orientation of the hexagonal mesh, either "x" or "y". (For hexagonal mesh only.)
|
||||
|
||||
:origin:
|
||||
The origin in cartesian coordinates. (For cylindrical and spherical meshes only.)
|
||||
The origin in cartesian coordinates. (For cylindrical, spherical and hexagonal meshes only.)
|
||||
|
||||
:library:
|
||||
The mesh library used to represent an unstructured mesh. This can be either
|
||||
|
|
|
|||
|
|
@ -170,6 +170,7 @@ Meshes
|
|||
openmc.CylindricalMesh
|
||||
openmc.SphericalMesh
|
||||
openmc.UnstructuredMesh
|
||||
openmc.HexagonalMesh
|
||||
|
||||
Geometry Plotting
|
||||
-----------------
|
||||
|
|
|
|||
|
|
@ -302,10 +302,10 @@ public:
|
|||
|
||||
struct MeshDistance {
|
||||
MeshDistance() = default;
|
||||
MeshDistance(int _index, bool _max_surface, double _distance)
|
||||
: next_index {_index}, max_surface {_max_surface}, distance {_distance}
|
||||
MeshDistance(MeshIndex _offset, bool _max_surface, double _distance)
|
||||
: offset {_offset}, max_surface {_max_surface}, distance {_distance}
|
||||
{}
|
||||
int next_index {-1};
|
||||
MeshIndex offset {0, 0, 0};
|
||||
bool max_surface {true};
|
||||
double distance {INFTY};
|
||||
bool operator<(const MeshDistance& o) const
|
||||
|
|
@ -314,6 +314,8 @@ public:
|
|||
}
|
||||
};
|
||||
|
||||
virtual void sanitize_index(MeshIndex& idx) const {};
|
||||
|
||||
Position sample_element(int32_t bin, uint64_t* seed) const override
|
||||
{
|
||||
return sample_element(get_indices_from_bin(bin), seed);
|
||||
|
|
@ -364,6 +366,13 @@ public:
|
|||
//! \return Array of mesh indices
|
||||
virtual MeshIndex get_indices(Position r, bool& in_mesh) const;
|
||||
|
||||
//! Check if mesh indices are inside mesh
|
||||
//
|
||||
//! \param[in] Array of mesh indices
|
||||
//! \param[in] k Suspect axis
|
||||
//! \return are indices inside mesh?
|
||||
virtual bool valid_index(const MeshIndex& ijk, int k) const;
|
||||
|
||||
//! Get mesh indices corresponding to a mesh bin
|
||||
//
|
||||
//! \param[in] bin Mesh bin
|
||||
|
|
@ -416,9 +425,16 @@ public:
|
|||
virtual MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
|
||||
const Position& r0, const Direction& u, double l) const = 0;
|
||||
|
||||
virtual double distance_to_mesh(const MeshIndex& ijk,
|
||||
const std::array<MeshDistance, 4>& distances, double traveled_distance,
|
||||
int& k_max) const;
|
||||
|
||||
//! Get a label for the mesh bin
|
||||
std::string bin_label(int bin) const override;
|
||||
|
||||
//! Get a label for the surface
|
||||
virtual std::string surface_label(int surface) const;
|
||||
|
||||
//! Get mesh dimensions as a tensor
|
||||
tensor::Tensor<int> get_shape_tensor() const;
|
||||
|
||||
|
|
@ -452,6 +468,8 @@ public:
|
|||
|
||||
// Data members
|
||||
std::array<int, 3> shape_; //!< Number of mesh elements in each dimension
|
||||
std::vector<std::vector<int>> correlated_axes_ = {{0}, {1}, {2}};
|
||||
std::vector<std::string> axes_labels_ = {"x", "y", "z"};
|
||||
|
||||
protected:
|
||||
};
|
||||
|
|
@ -580,6 +598,11 @@ public:
|
|||
CylindricalMesh(hid_t group);
|
||||
|
||||
// Overridden methods
|
||||
void sanitize_index(MeshIndex& idx) const override
|
||||
{
|
||||
idx[1] = sanitize_phi(idx[1]);
|
||||
}
|
||||
|
||||
virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
|
||||
|
||||
int get_index_in_direction(double r, int i) const override;
|
||||
|
|
@ -645,6 +668,12 @@ public:
|
|||
SphericalMesh(hid_t group);
|
||||
|
||||
// Overridden methods
|
||||
void sanitize_index(MeshIndex& idx) const override
|
||||
{
|
||||
idx[1] = sanitize_theta(idx[1]);
|
||||
idx[2] = sanitize_phi(idx[2]);
|
||||
}
|
||||
|
||||
virtual MeshIndex get_indices(Position r, bool& in_mesh) const override;
|
||||
|
||||
int get_index_in_direction(double r, int i) const override;
|
||||
|
|
@ -706,6 +735,84 @@ private:
|
|||
}
|
||||
};
|
||||
|
||||
class HexagonalMesh : public PeriodicStructuredMesh {
|
||||
public:
|
||||
// Constructors
|
||||
HexagonalMesh() = default;
|
||||
HexagonalMesh(pugi::xml_node node);
|
||||
HexagonalMesh(hid_t group);
|
||||
|
||||
// Overridden methods
|
||||
|
||||
//! Get a label for the mesh bin
|
||||
std::string bin_label(int bin) const override;
|
||||
|
||||
std::string surface_label(int surface) const override;
|
||||
|
||||
int n_bins() const override;
|
||||
|
||||
MeshIndex get_indices(Position r, bool& in_mesh) const override;
|
||||
|
||||
//! Get mesh indices corresponding to a mesh bin
|
||||
//
|
||||
//! \param[in] bin Mesh bin
|
||||
//! \return ijk Mesh indices
|
||||
MeshIndex get_indices_from_bin(int bin) const override;
|
||||
|
||||
bool valid_index(const MeshIndex& ijk, int k) const override;
|
||||
|
||||
int get_bin_from_indices(const MeshIndex& ijk) const override;
|
||||
|
||||
std::string get_mesh_type() const override;
|
||||
|
||||
static const std::string mesh_type;
|
||||
|
||||
Position sample_element(const MeshIndex& ijk, uint64_t* seed) const override;
|
||||
|
||||
MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i,
|
||||
const Position& r0, const Direction& u, double l) const override;
|
||||
|
||||
double distance_to_mesh(const MeshIndex& ijk,
|
||||
const std::array<MeshDistance, 4>& distances, double traveled_distance,
|
||||
int& k_max) const override;
|
||||
|
||||
std::pair<vector<double>, vector<double>> plot(
|
||||
Position plot_ll, Position plot_ur) const override;
|
||||
|
||||
void to_hdf5_inner(hid_t group) const override;
|
||||
|
||||
int set_grid();
|
||||
|
||||
enum class Orientation {
|
||||
y, //!< Flat side of lattice parallel to y-axis
|
||||
x //!< Flat side of lattice parallel to x-axis
|
||||
};
|
||||
|
||||
// Data members
|
||||
int num_rings_;
|
||||
double pitch_;
|
||||
vector<double> grid_;
|
||||
Direction q_;
|
||||
Direction r_;
|
||||
Direction q_dual_;
|
||||
Direction r_dual_;
|
||||
Orientation orientation_ {Orientation::y};
|
||||
|
||||
private:
|
||||
StructuredMesh::MeshDistance find_z_crossing(
|
||||
const Position& r, const Direction& u, double l, int shell) const;
|
||||
|
||||
double volume(const MeshIndex& ijk) const override;
|
||||
|
||||
int get_index_in_direction(double r, int i) const override
|
||||
{
|
||||
fatal_error("This function is not implemented");
|
||||
return -1;
|
||||
}
|
||||
|
||||
int get_index_in_z_direction(double z) const;
|
||||
};
|
||||
|
||||
// Abstract class for unstructured meshes
|
||||
class UnstructuredMesh : public Mesh {
|
||||
|
||||
|
|
|
|||
|
|
@ -22,21 +22,6 @@ public:
|
|||
// Accessors
|
||||
|
||||
void set_mesh(int32_t mesh) override;
|
||||
|
||||
enum class MeshDir {
|
||||
OUT_LEFT, // x min
|
||||
IN_LEFT, // x min
|
||||
OUT_RIGHT, // x max
|
||||
IN_RIGHT, // x max
|
||||
OUT_BACK, // y min
|
||||
IN_BACK, // y min
|
||||
OUT_FRONT, // y max
|
||||
IN_FRONT, // y max
|
||||
OUT_BOTTOM, // z min
|
||||
IN_BOTTOM, // z min
|
||||
OUT_TOP, // z max
|
||||
IN_TOP // z max
|
||||
};
|
||||
};
|
||||
|
||||
} // namespace openmc
|
||||
|
|
|
|||
|
|
@ -30,13 +30,6 @@ _FILTER_TYPES = (
|
|||
'weight', 'meshborn', 'meshsurface', 'meshmaterial', 'reaction',
|
||||
)
|
||||
|
||||
_CURRENT_NAMES = (
|
||||
'x-min out', 'x-min in', 'x-max out', 'x-max in',
|
||||
'y-min out', 'y-min in', 'y-max out', 'y-max in',
|
||||
'z-min out', 'z-min in', 'z-max out', 'z-max in'
|
||||
)
|
||||
|
||||
|
||||
|
||||
class FilterMeta(ABCMeta):
|
||||
"""Metaclass for filters that ensures class names are appropriate."""
|
||||
|
|
@ -79,7 +72,7 @@ class FilterMeta(ABCMeta):
|
|||
|
||||
|
||||
def _repeat_and_tile(bins, repeat_factor, data_size):
|
||||
filter_bins = np.repeat(bins, repeat_factor)
|
||||
filter_bins = np.repeat(bins, repeat_factor, axis=0)
|
||||
tile_factor = data_size // len(filter_bins)
|
||||
return np.tile(filter_bins, tile_factor)
|
||||
|
||||
|
|
@ -992,16 +985,9 @@ class MeshFilter(Filter):
|
|||
# Append mesh ID as outermost index of multi-index
|
||||
mesh_key = f'mesh {self.mesh.id}'
|
||||
|
||||
# Determine index base (0-based for unstructured, 1-based otherwise)
|
||||
idx_start = 0 if isinstance(self.mesh, openmc.UnstructuredMesh) else 1
|
||||
|
||||
# Generate a multi-index sub-column for each axis
|
||||
for label, dim_size in zip(self.mesh._axis_labels, self.mesh.dimension):
|
||||
filter_dict[mesh_key, label] = _repeat_and_tile(
|
||||
np.arange(idx_start, idx_start + dim_size), stride, data_size)
|
||||
stride *= dim_size
|
||||
|
||||
return pd.DataFrame(filter_dict)
|
||||
columns = [(mesh_key, label) for label in self.mesh._axis_labels]
|
||||
indices = _repeat_and_tile(list(self.mesh.indices), stride, data_size)
|
||||
return pd.DataFrame(indices, columns=columns)
|
||||
|
||||
def to_xml_element(self):
|
||||
"""Return XML Element representing the Filter.
|
||||
|
|
@ -1282,15 +1268,34 @@ class MeshSurfaceFilter(MeshFilter):
|
|||
def shape(self):
|
||||
return (self.num_bins,)
|
||||
|
||||
@staticmethod
|
||||
def _current_names(mesh):
|
||||
if isinstance(mesh, openmc.HexagonalMesh):
|
||||
names = []
|
||||
ax0, ax1, ax2, ax3 = mesh._axis_labels
|
||||
pairs = [(ax0, ax1), (ax1, ax2), (ax2, ax0)]
|
||||
for ax0, ax1 in pairs:
|
||||
for minmax0, minmax1 in [('max','min'),('min','max')]:
|
||||
for inout in ('out','in'):
|
||||
names.append(f"{ax0}-{minmax0} {ax1}-{minmax1} {inout}")
|
||||
for minmax in ('min', 'max'):
|
||||
for inout in ('out','in'):
|
||||
names.append(f"{ax3}-{minmax} {inout}")
|
||||
return names
|
||||
|
||||
names = []
|
||||
for ax in mesh._axis_labels:
|
||||
for minmax in ('min', 'max'):
|
||||
for inout in ('out','in'):
|
||||
names.append(f"{ax}-{minmax} {inout}")
|
||||
return names
|
||||
|
||||
@MeshFilter.mesh.setter
|
||||
def mesh(self, mesh):
|
||||
cv.check_type('filter mesh', mesh, openmc.MeshBase)
|
||||
self._mesh = mesh
|
||||
|
||||
# Take the product of mesh indices and current names
|
||||
n_dim = mesh.n_dimension
|
||||
self.bins = [mesh_tuple + (surf,) for mesh_tuple, surf in
|
||||
product(mesh.indices, _CURRENT_NAMES[:4*n_dim])]
|
||||
product(mesh.indices, self._current_names(mesh))]
|
||||
|
||||
def get_pandas_dataframe(self, data_size, stride, **kwargs):
|
||||
"""Builds a Pandas DataFrame for the Filter's bins.
|
||||
|
|
@ -1320,47 +1325,19 @@ class MeshSurfaceFilter(MeshFilter):
|
|||
Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
|
||||
|
||||
"""
|
||||
# Initialize Pandas DataFrame
|
||||
df = pd.DataFrame()
|
||||
|
||||
# Initialize dictionary to build Pandas Multi-index column
|
||||
filter_dict = {}
|
||||
|
||||
# Append mesh ID as outermost index of multi-index
|
||||
mesh_key = f'mesh {self.mesh.id}'
|
||||
|
||||
# Find mesh dimensions - use 3D indices for simplicity
|
||||
n_surfs = 4 * len(self.mesh.dimension)
|
||||
if len(self.mesh.dimension) == 3:
|
||||
nx, ny, nz = self.mesh.dimension
|
||||
elif len(self.mesh.dimension) == 2:
|
||||
nx, ny = self.mesh.dimension
|
||||
nz = 1
|
||||
else:
|
||||
nx = self.mesh.dimension
|
||||
ny = nz = 1
|
||||
n_surfs = len(self._current_names(self._mesh))
|
||||
|
||||
# Generate multi-index sub-column for x-axis
|
||||
filter_dict[mesh_key, 'x'] = _repeat_and_tile(
|
||||
np.arange(1, nx + 1), n_surfs * stride, data_size)
|
||||
|
||||
# Generate multi-index sub-column for y-axis
|
||||
if len(self.mesh.dimension) > 1:
|
||||
filter_dict[mesh_key, 'y'] = _repeat_and_tile(
|
||||
np.arange(1, ny + 1), n_surfs * nx * stride, data_size)
|
||||
|
||||
# Generate multi-index sub-column for z-axis
|
||||
if len(self.mesh.dimension) > 2:
|
||||
filter_dict[mesh_key, 'z'] = _repeat_and_tile(
|
||||
np.arange(1, nz + 1), n_surfs * nx * ny * stride, data_size)
|
||||
|
||||
# Generate multi-index sub-column for surface
|
||||
filter_dict[mesh_key, 'surf'] = _repeat_and_tile(
|
||||
_CURRENT_NAMES[:n_surfs], stride, data_size)
|
||||
|
||||
# Initialize a Pandas DataFrame from the mesh dictionary
|
||||
return pd.concat([df, pd.DataFrame(filter_dict)])
|
||||
columns = [(mesh_key, label) for label in self.mesh._axis_labels]
|
||||
indices = _repeat_and_tile(list(self.mesh.indices), stride*n_surfs, data_size)
|
||||
filter_dict = dict(zip(columns,indices.T))
|
||||
surfs = _repeat_and_tile(self._current_names(self._mesh), stride, data_size)
|
||||
filter_dict[mesh_key, 'surf'] = surfs
|
||||
|
||||
return pd.DataFrame(filter_dict)
|
||||
|
||||
class CollisionFilter(Filter):
|
||||
"""Bins tally events based on the number of collisions.
|
||||
|
|
|
|||
|
|
@ -18,7 +18,7 @@ from ..mesh import MeshMaterialVolumes
|
|||
|
||||
__all__ = [
|
||||
'Mesh', 'RegularMesh', 'RectilinearMesh', 'CylindricalMesh',
|
||||
'SphericalMesh', 'UnstructuredMesh', 'meshes', 'MeshMaterialVolumes'
|
||||
'SphericalMesh', 'HexagonalMesh', 'UnstructuredMesh', 'meshes', 'MeshMaterialVolumes'
|
||||
]
|
||||
|
||||
|
||||
|
|
@ -108,6 +108,11 @@ _dll.openmc_spherical_mesh_set_grid.argtypes = [c_int32, POINTER(c_double),
|
|||
_dll.openmc_spherical_mesh_set_grid.restype = c_int
|
||||
_dll.openmc_spherical_mesh_set_grid.errcheck = _error_handler
|
||||
|
||||
_dll.openmc_hexagonal_mesh_get_grid.argtypes = [c_int32,
|
||||
POINTER(POINTER(c_double)), POINTER(c_int), POINTER(c_double),
|
||||
POINTER(POINTER(c_double)), POINTER(c_double), POINTER(c_char_p)]
|
||||
_dll.openmc_hexagonal_mesh_get_grid.restype = c_int
|
||||
_dll.openmc_hexagonal_mesh_get_grid.errcheck = _error_handler
|
||||
|
||||
class Mesh(_FortranObjectWithID):
|
||||
"""Base class to represent mesh objects
|
||||
|
|
@ -727,6 +732,91 @@ class SphericalMesh(Mesh):
|
|||
ntheta, phi_grid, nphi)
|
||||
|
||||
|
||||
class HexagonalMesh(Mesh):
|
||||
"""HexagonalMesh 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
|
||||
z_grid : numpy.ndarray
|
||||
1-D array of mesh boundary points along the z-axis.
|
||||
pitch : float
|
||||
Radial pitch of the hexagonal mesh in cm.
|
||||
num_rings : int
|
||||
Number of radial ring positions in the xy-plane
|
||||
orientation : {'x', 'y'}
|
||||
The orientation of the lattice. The 'x' orientation means that each
|
||||
lattice element has two faces that are perpendicular to the x-axis,
|
||||
while the 'y' orientation means that each lattice element has two faces
|
||||
that are perpendicular to the y-axis. By default, the orientation is
|
||||
'y'.
|
||||
origin : numpy.ndarray
|
||||
1-D array of length 3 the (x,y,z) origin of the mesh in
|
||||
cartesian coordinates
|
||||
n_elements : int
|
||||
Total number of mesh elements.
|
||||
volumes : numpy.ndarray
|
||||
Volume of each mesh element in [cm^3]
|
||||
bounding_box : openmc.BoundingBox
|
||||
Axis-aligned bounding box of the mesh
|
||||
|
||||
"""
|
||||
mesh_type = 'hexagonal'
|
||||
|
||||
def __init__(self, uid=None, new=True, index=None):
|
||||
super().__init__(uid, new, index)
|
||||
|
||||
@property
|
||||
def n_elements(self):
|
||||
z_grid, nr, *_ = self._get_parameters()
|
||||
return (z_grid.size-1)*(3*nr*(nr-1)+1)
|
||||
|
||||
@property
|
||||
def num_rings(self):
|
||||
return self._get_parameters()[1]
|
||||
|
||||
@property
|
||||
def pitch(self):
|
||||
return self._get_parameters()[2]
|
||||
|
||||
@property
|
||||
def orientation(self):
|
||||
return self._get_parameters()[4]
|
||||
|
||||
@property
|
||||
def origin(self):
|
||||
return self._get_parameters()[3]
|
||||
|
||||
@property
|
||||
def z_grid(self):
|
||||
return self._get_parameters()[0]
|
||||
|
||||
def _get_parameters(self):
|
||||
gz = POINTER(c_double)()
|
||||
nz = c_int()
|
||||
nr = c_int()
|
||||
orig = POINTER(c_double)()
|
||||
orient = c_char()
|
||||
p = c_double()
|
||||
# Call C API to get grid parameters
|
||||
_dll.openmc_hexagonal_mesh_get_grid(self._index, gz, nz, nr, orig, p, orient)
|
||||
|
||||
# Convert grid parameters to Numpy arrays
|
||||
grid_z = as_array(gz, (nz.value,))
|
||||
origin = as_array(orig, (3,))
|
||||
|
||||
return (grid_z, nr, pitch, origin, orientation.decode())
|
||||
|
||||
class UnstructuredMesh(Mesh):
|
||||
pass
|
||||
|
||||
|
|
@ -736,6 +826,7 @@ _MESH_TYPE_MAP = {
|
|||
'rectilinear': RectilinearMesh,
|
||||
'cylindrical': CylindricalMesh,
|
||||
'spherical': SphericalMesh,
|
||||
'hexagonal': HexagonalMesh,
|
||||
'unstructured': UnstructuredMesh
|
||||
}
|
||||
|
||||
|
|
|
|||
287
openmc/mesh.py
287
openmc/mesh.py
|
|
@ -299,6 +299,8 @@ class MeshBase(IDManagerMixin, ABC):
|
|||
return CylindricalMesh.from_hdf5(group, mesh_id, mesh_name)
|
||||
elif mesh_type == 'spherical':
|
||||
return SphericalMesh.from_hdf5(group, mesh_id, mesh_name)
|
||||
elif mesh_type == 'hexagonal':
|
||||
return HexagonalMesh.from_hdf5(group, mesh_id, mesh_name)
|
||||
elif mesh_type == 'unstructured':
|
||||
return UnstructuredMesh.from_hdf5(group, mesh_id, mesh_name)
|
||||
else:
|
||||
|
|
@ -526,11 +528,6 @@ class StructuredMesh(MeshBase):
|
|||
def __init__(self, *args, **kwargs):
|
||||
super().__init__(*args, **kwargs)
|
||||
|
||||
@property
|
||||
@abstractmethod
|
||||
def dimension(self):
|
||||
pass
|
||||
|
||||
@property
|
||||
@abstractmethod
|
||||
def n_dimension(self):
|
||||
|
|
@ -541,11 +538,6 @@ class StructuredMesh(MeshBase):
|
|||
def _axis_labels(self):
|
||||
pass
|
||||
|
||||
@property
|
||||
@abstractmethod
|
||||
def _grids(self):
|
||||
pass
|
||||
|
||||
@abstractmethod
|
||||
def get_indices_at_coords(self, coords: Sequence[float]) -> tuple:
|
||||
pass
|
||||
|
|
@ -2720,6 +2712,281 @@ class SphericalMesh(StructuredMesh):
|
|||
|
||||
return (r_index, theta_index, phi_index)
|
||||
|
||||
class HexagonalMesh(StructuredMesh):
|
||||
"""A 3D hexagonal mesh
|
||||
|
||||
Parameters
|
||||
----------
|
||||
z_grid : numpy.ndarray
|
||||
1-D array of mesh boundary points along the z-axis.
|
||||
pitch : float
|
||||
Radial pitch of the hexagonal mesh in cm.
|
||||
num_rings : int
|
||||
Number of radial ring positions in the xy-plane
|
||||
orientation : {'x', 'y'}
|
||||
The orientation of the lattice. The 'x' orientation means that each
|
||||
lattice element has two faces that are perpendicular to the x-axis,
|
||||
while the 'y' orientation means that each lattice element has two faces
|
||||
that are perpendicular to the y-axis. By default, the orientation is
|
||||
'y'.
|
||||
origin : numpy.ndarray
|
||||
1-D array of length 3 the (x,y,z) origin of the mesh in
|
||||
cartesian coordinates
|
||||
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
|
||||
z_grid : numpy.ndarray
|
||||
1-D array of mesh boundary points along the z-axis.
|
||||
pitch : float
|
||||
Radial pitch of the hexagonal mesh in cm.
|
||||
num_rings : int
|
||||
Number of radial ring positions in the xy-plane
|
||||
orientation : {'x', 'y'}
|
||||
The orientation of the lattice. The 'x' orientation means that each
|
||||
lattice element has two faces that are perpendicular to the x-axis,
|
||||
while the 'y' orientation means that each lattice element has two faces
|
||||
that are perpendicular to the y-axis. By default, the orientation is
|
||||
'y'.
|
||||
origin : numpy.ndarray
|
||||
1-D array of length 3 the (x,y,z) origin of the mesh in
|
||||
cartesian coordinates
|
||||
indices : Iterable of tuple
|
||||
An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1),
|
||||
(2, 1, 1), ...]
|
||||
lower_left : numpy.ndarray
|
||||
The lower-left corner of the structured mesh. If only two coordinate
|
||||
are given, it is assumed that the mesh is an x-y mesh.
|
||||
upper_right : numpy.ndarray
|
||||
The upper-right corner of the structured mesh. If only two coordinate
|
||||
are given, it is assumed that the mesh is an x-y mesh.
|
||||
bounding_box : openmc.BoundingBox
|
||||
Axis-aligned bounding box of the mesh as defined by the upper-right and
|
||||
lower-left coordinates.
|
||||
|
||||
"""
|
||||
|
||||
def __init__(
|
||||
self,
|
||||
z_grid: Sequence[float],
|
||||
pitch: float,
|
||||
num_rings: int,
|
||||
orientation: str = 'y',
|
||||
origin: Sequence[float] = (0., 0., 0.),
|
||||
mesh_id: int | None = None,
|
||||
name: str = '',
|
||||
):
|
||||
super().__init__(mesh_id, name)
|
||||
|
||||
self.z_grid = z_grid
|
||||
self.pitch = pitch
|
||||
self.num_rings = num_rings
|
||||
self.orientation = orientation
|
||||
self.origin = origin
|
||||
|
||||
@property
|
||||
def _axis_labels(self):
|
||||
return ('r', 'q', 's', 'z')
|
||||
|
||||
@property
|
||||
def origin(self):
|
||||
return self._origin
|
||||
|
||||
@origin.setter
|
||||
def origin(self, coords):
|
||||
cv.check_type('mesh origin', coords, Iterable, Real)
|
||||
cv.check_length("mesh origin", coords, 3)
|
||||
self._origin = np.asarray(coords, dtype=float)
|
||||
|
||||
@property
|
||||
def z_grid(self):
|
||||
return self._z_grid
|
||||
|
||||
@z_grid.setter
|
||||
def z_grid(self, grid):
|
||||
cv.check_type('mesh z_grid', grid, Iterable, Real)
|
||||
cv.check_length('mesh z_grid', grid, 2)
|
||||
cv.check_increasing('mesh z_grid', grid)
|
||||
self._z_grid = np.asarray(grid, dtype=float)
|
||||
|
||||
@property
|
||||
def pitch(self):
|
||||
return self._pitch
|
||||
|
||||
@pitch.setter
|
||||
def pitch(self, pitch):
|
||||
cv.check_type('mesh radial pitch', pitch, Real)
|
||||
cv.check_greater_than('mesh radial pitch', pitch, 0.0)
|
||||
self._pitch = pitch
|
||||
|
||||
@property
|
||||
def n_elements(self):
|
||||
return (self.z_grid.size - 1)*(3*self.num_rings*(self.num_rings - 1) + 1)
|
||||
|
||||
@property
|
||||
def dimension(self):
|
||||
return (self.n_elements,)
|
||||
|
||||
@property
|
||||
def n_dimension(self):
|
||||
return 4
|
||||
|
||||
@property
|
||||
def lower_left(self):
|
||||
p = self.pitch
|
||||
r = self.pitch/np.sqrt(3.0)
|
||||
n = self.num_rings
|
||||
offset = np.array([0,0,self.z_grid[0]])
|
||||
if self.orientation == 'x':
|
||||
offset[0] = -p*(n-0.5)
|
||||
offset[1] = -(1+1.5*(n-1))*r
|
||||
else:
|
||||
offset[0] = -(1+1.5*(n-1))*r
|
||||
offset[1] = -p*(n-0.5)
|
||||
return self.origin + offset
|
||||
|
||||
@property
|
||||
def upper_right(self):
|
||||
p = self.pitch
|
||||
r = self.pitch/np.sqrt(3.0)
|
||||
n = self.num_rings
|
||||
offset = np.array([0,0,self.z_grid[-1]])
|
||||
if self.orientation == 'x':
|
||||
offset[0] = p*(n-0.5)
|
||||
offset[1] = (1+1.5*(n-1))*r
|
||||
else:
|
||||
offset[0] = (1+1.5*(n-1))*r
|
||||
offset[1] = p*(n-0.5)
|
||||
return self.origin + offset
|
||||
|
||||
@property
|
||||
def num_rings(self):
|
||||
return self._num_rings
|
||||
|
||||
@num_rings.setter
|
||||
def num_rings(self, num_rings):
|
||||
cv.check_type('mesh num_rings', num_rings, Integral)
|
||||
cv.check_greater_than('mesh num_rings', num_rings, 0)
|
||||
self._num_rings = num_rings
|
||||
|
||||
@property
|
||||
def orientation(self):
|
||||
return self._orientation
|
||||
|
||||
@orientation.setter
|
||||
def orientation(self, orientation):
|
||||
cv.check_value('orientation', orientation.lower(), ('x', 'y'))
|
||||
self._orientation = orientation.lower()
|
||||
|
||||
@property
|
||||
def indices(self):
|
||||
idx = []
|
||||
for j in range(self.z_grid.size-1):
|
||||
idx.append((0, 0, 0, j))
|
||||
for rad in range(1, self.num_rings):
|
||||
for i in range(rad):
|
||||
idx.append((i, rad-i, -rad, j))
|
||||
idx.append((rad, -i, i-rad, j))
|
||||
idx.append((rad-i, -rad, i, j))
|
||||
idx.append((-i, i-rad, rad, j))
|
||||
idx.append((-rad, i, rad-i, j))
|
||||
idx.append((i-rad, rad, -i, j))
|
||||
return idx
|
||||
|
||||
def __repr__(self):
|
||||
fmt = '{0: <16}{1}{2}\n'
|
||||
string = super().__repr__()
|
||||
string += fmt.format('\tOrigin', '=\t', self.origin)
|
||||
string += fmt.format('\tPitch', '=\t', self.pitch)
|
||||
string += fmt.format('\tN rings', '=\t', self.num_rings)
|
||||
string += fmt.format('\tOrientation', '=\t', self.orientation)
|
||||
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: h5py.Group, mesh_id: int, name: str):
|
||||
# Read and assign mesh properties
|
||||
mesh = cls(
|
||||
z_grid = group['z_grid'][()],
|
||||
pitch = group.attrs['pitch'],
|
||||
num_rings = group.attrs['num_rings'],
|
||||
mesh_id=mesh_id,
|
||||
name=name
|
||||
)
|
||||
if 'origin' in group:
|
||||
mesh.origin = group['origin'][()]
|
||||
if 'orientation' in group.attrs:
|
||||
mesh.orientation = group.attrs['orientation']
|
||||
|
||||
return mesh
|
||||
|
||||
def to_xml_element(self):
|
||||
"""Return XML representation of the mesh
|
||||
|
||||
Returns
|
||||
-------
|
||||
element : lxml.etree._Element
|
||||
XML element containing mesh data
|
||||
|
||||
"""
|
||||
|
||||
element = super().to_xml_element()
|
||||
element.set("type", "hexagonal")
|
||||
element.set("pitch", str(self.pitch))
|
||||
element.set("num_rings", str(self.num_rings))
|
||||
element.set("orientation", str(self.orientation))
|
||||
|
||||
subelement = ET.SubElement(element, "z_grid")
|
||||
subelement.text = ' '.join(map(str, self.z_grid))
|
||||
|
||||
subelement = ET.SubElement(element, "origin")
|
||||
subelement.text = ' '.join(map(str, self.origin))
|
||||
|
||||
return element
|
||||
|
||||
@classmethod
|
||||
def from_xml_element(cls, elem: ET.Element):
|
||||
"""Generate a spherical mesh from an XML element
|
||||
|
||||
Parameters
|
||||
----------
|
||||
elem : lxml.etree._Element
|
||||
XML element
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.SphericalMesh
|
||||
Spherical mesh object
|
||||
|
||||
"""
|
||||
mesh_id = int(get_text(elem, 'id'))
|
||||
mesh = cls(
|
||||
mesh_id=mesh_id,
|
||||
z_grid = get_elem_list(elem, "z_grid", float),
|
||||
pitch = float(get_text(elem, "pitch")),
|
||||
num_rings = int(get_text(elem, "num_rings")),
|
||||
orientation = get_text(elem, "orientation"),
|
||||
origin = get_elem_list(elem, "origin", float) or [0., 0., 0.],
|
||||
)
|
||||
|
||||
return mesh
|
||||
|
||||
def get_indices_at_coords(self, coords: Sequence[float]) -> tuple:
|
||||
raise NotImplementedError(
|
||||
"get_indices_at_coords is not yet implemented for SphericalMesh"
|
||||
)
|
||||
|
||||
|
||||
def require_statepoint_data(func):
|
||||
@wraps(func)
|
||||
|
|
|
|||
564
src/mesh.cpp
564
src/mesh.cpp
|
|
@ -343,6 +343,8 @@ const std::unique_ptr<Mesh>& Mesh::create(
|
|||
model::meshes.push_back(make_unique<CylindricalMesh>(dataset));
|
||||
} else if (mesh_type == SphericalMesh::mesh_type) {
|
||||
model::meshes.push_back(make_unique<SphericalMesh>(dataset));
|
||||
} else if (mesh_type == HexagonalMesh::mesh_type) {
|
||||
model::meshes.push_back(make_unique<HexagonalMesh>(dataset));
|
||||
#ifdef OPENMC_DAGMC_ENABLED
|
||||
} else if (mesh_type == UnstructuredMesh::mesh_type &&
|
||||
mesh_library == MOABMesh::mesh_lib_type) {
|
||||
|
|
@ -768,6 +770,17 @@ std::string StructuredMesh::bin_label(int bin) const
|
|||
}
|
||||
}
|
||||
|
||||
std::string StructuredMesh::surface_label(int surface) const
|
||||
{
|
||||
int axis = static_cast<int>(std::floor(surface / 4));
|
||||
std::string label = axes_labels_[axis];
|
||||
const std::vector<std::string> MINMAX = {"min", "max"};
|
||||
const std::vector<std::string> OUTIN = {"Outgoing", "Incoming"};
|
||||
auto minmax = MINMAX[static_cast<int>(std::floor((surface % 4) / 2))];
|
||||
auto outin = OUTIN[surface % 2];
|
||||
return fmt::format("{}, {}-{}", outin, minmax, label);
|
||||
}
|
||||
|
||||
tensor::Tensor<int> StructuredMesh::get_shape_tensor() const
|
||||
{
|
||||
return tensor::Tensor<int>(shape_.data(), static_cast<size_t>(n_dimension_));
|
||||
|
|
@ -1037,6 +1050,11 @@ StructuredMesh::MeshIndex StructuredMesh::get_indices(
|
|||
return ijk;
|
||||
}
|
||||
|
||||
bool StructuredMesh::valid_index(const MeshIndex& ijk, int k) const
|
||||
{
|
||||
return ((ijk[k] >= 1) && (ijk[k] <= shape_[k]));
|
||||
}
|
||||
|
||||
int StructuredMesh::get_bin_from_indices(const MeshIndex& ijk) const
|
||||
{
|
||||
switch (n_dimension_) {
|
||||
|
|
@ -1197,14 +1215,13 @@ void StructuredMesh::raytrace_mesh(
|
|||
}
|
||||
|
||||
// Calculate initial distances to next surfaces in all three dimensions
|
||||
std::array<MeshDistance, 3> distances;
|
||||
std::array<MeshDistance, 4> distances;
|
||||
for (int k = 0; k < n; ++k) {
|
||||
distances[k] = distance_to_grid_boundary(ijk, k, local_r, u, 0.0);
|
||||
}
|
||||
|
||||
// Loop until r = r1 is eventually reached
|
||||
while (true) {
|
||||
|
||||
if (in_mesh) {
|
||||
|
||||
// find surface with minimal distance to current position
|
||||
|
|
@ -1225,14 +1242,17 @@ void StructuredMesh::raytrace_mesh(
|
|||
// current
|
||||
tally.surface(ijk, k, distances[k].max_surface, false);
|
||||
|
||||
// Update cell and calculate distance to next surface in k-direction.
|
||||
// The two other directions are still valid!
|
||||
ijk[k] = distances[k].next_index;
|
||||
distances[k] =
|
||||
distance_to_grid_boundary(ijk, k, local_r, u, traveled_distance);
|
||||
// Update cell and calculate distance to next surface in correlated
|
||||
// directions.
|
||||
for (int j = 0; j < 3; ++j)
|
||||
ijk[j] += distances[k].offset[j];
|
||||
sanitize_index(ijk);
|
||||
for (auto j : correlated_axes_[k])
|
||||
distances[j] =
|
||||
distance_to_grid_boundary(ijk, j, local_r, u, traveled_distance);
|
||||
|
||||
// Check if we have left the interior of the mesh
|
||||
in_mesh = ((ijk[k] >= 1) && (ijk[k] <= shape_[k]));
|
||||
in_mesh = valid_index(ijk, k);
|
||||
|
||||
// If we are still inside the mesh, tally inward current for the next
|
||||
// cell
|
||||
|
|
@ -1240,23 +1260,9 @@ void StructuredMesh::raytrace_mesh(
|
|||
tally.surface(ijk, k, !distances[k].max_surface, true);
|
||||
|
||||
} else { // not inside mesh
|
||||
|
||||
// For all directions outside the mesh, find the distance that we need
|
||||
// to travel to reach the next surface. Use the largest distance, as
|
||||
// only this will cross all outer surfaces.
|
||||
int k_max {-1};
|
||||
for (int k = 0; k < n; ++k) {
|
||||
if ((ijk[k] < 1 || ijk[k] > shape_[k]) &&
|
||||
(distances[k].distance > traveled_distance)) {
|
||||
traveled_distance = distances[k].distance;
|
||||
k_max = k;
|
||||
}
|
||||
}
|
||||
// Assure some distance is traveled
|
||||
if (k_max == -1) {
|
||||
traveled_distance += TINY_BIT;
|
||||
}
|
||||
|
||||
int k_max;
|
||||
traveled_distance =
|
||||
distance_to_mesh(ijk, distances, traveled_distance, k_max);
|
||||
// If r1 is not inside the mesh, exit here
|
||||
if (traveled_distance >= total_distance)
|
||||
return;
|
||||
|
|
@ -1276,6 +1282,29 @@ void StructuredMesh::raytrace_mesh(
|
|||
}
|
||||
}
|
||||
|
||||
double StructuredMesh::distance_to_mesh(const MeshIndex& ijk,
|
||||
const std::array<MeshDistance, 4>& distances, double traveled_distance,
|
||||
int& k_max) const
|
||||
{
|
||||
// For all directions outside the mesh, find the distance that we need
|
||||
// to travel to reach the next surface. Use the largest distance, as
|
||||
// only this will cross all outer surfaces.
|
||||
const int n = n_dimension_;
|
||||
k_max = -1;
|
||||
for (int k = 0; k < n; ++k) {
|
||||
if ((ijk[k] < 1 || ijk[k] > shape_[k]) &&
|
||||
(distances[k].distance > traveled_distance)) {
|
||||
traveled_distance = distances[k].distance;
|
||||
k_max = k;
|
||||
}
|
||||
}
|
||||
// Assure some distance is traveled
|
||||
if (k_max == -1) {
|
||||
traveled_distance += TINY_BIT;
|
||||
}
|
||||
return traveled_distance;
|
||||
}
|
||||
|
||||
void StructuredMesh::bins_crossed(Position r0, Position r1, const Direction& u,
|
||||
vector<int>& bins, vector<double>& lengths) const
|
||||
{
|
||||
|
|
@ -1511,16 +1540,15 @@ StructuredMesh::MeshDistance RegularMesh::distance_to_grid_boundary(
|
|||
double l) const
|
||||
{
|
||||
MeshDistance d;
|
||||
d.next_index = ijk[i];
|
||||
if (std::abs(u[i]) < FP_PRECISION)
|
||||
return d;
|
||||
|
||||
d.max_surface = (u[i] > 0);
|
||||
if (d.max_surface && (ijk[i] <= shape_[i])) {
|
||||
d.next_index++;
|
||||
++d.offset[i];
|
||||
d.distance = (positive_grid_boundary(ijk, i) - r0[i]) / u[i];
|
||||
} else if (!d.max_surface && (ijk[i] >= 1)) {
|
||||
d.next_index--;
|
||||
--d.offset[i];
|
||||
d.distance = (negative_grid_boundary(ijk, i) - r0[i]) / u[i];
|
||||
}
|
||||
|
||||
|
|
@ -1684,16 +1712,15 @@ StructuredMesh::MeshDistance RectilinearMesh::distance_to_grid_boundary(
|
|||
double l) const
|
||||
{
|
||||
MeshDistance d;
|
||||
d.next_index = ijk[i];
|
||||
if (std::abs(u[i]) < FP_PRECISION)
|
||||
return d;
|
||||
|
||||
d.max_surface = (u[i] > 0);
|
||||
if (d.max_surface && (ijk[i] <= shape_[i])) {
|
||||
d.next_index++;
|
||||
++d.offset[i];
|
||||
d.distance = (positive_grid_boundary(ijk, i) - r0[i]) / u[i];
|
||||
} else if (!d.max_surface && (ijk[i] > 0)) {
|
||||
d.next_index--;
|
||||
--d.offset[i];
|
||||
d.distance = (negative_grid_boundary(ijk, i) - r0[i]) / u[i];
|
||||
}
|
||||
return d;
|
||||
|
|
@ -1946,7 +1973,6 @@ StructuredMesh::MeshDistance CylindricalMesh::find_z_crossing(
|
|||
const Position& r, const Direction& u, double l, int shell) const
|
||||
{
|
||||
MeshDistance d;
|
||||
d.next_index = shell;
|
||||
|
||||
// Direction of flight is within xy-plane. Will never intersect z.
|
||||
if (std::abs(u.z) < FP_PRECISION)
|
||||
|
|
@ -1954,10 +1980,10 @@ StructuredMesh::MeshDistance CylindricalMesh::find_z_crossing(
|
|||
|
||||
d.max_surface = (u.z > 0.0);
|
||||
if (d.max_surface && (shell <= shape_[2])) {
|
||||
d.next_index += 1;
|
||||
++d.offset[2];
|
||||
d.distance = (grid_[2][shell] - r.z) / u.z;
|
||||
} else if (!d.max_surface && (shell > 0)) {
|
||||
d.next_index -= 1;
|
||||
--d.offset[2];
|
||||
d.distance = (grid_[2][shell - 1] - r.z) / u.z;
|
||||
}
|
||||
return d;
|
||||
|
|
@ -1968,17 +1994,15 @@ StructuredMesh::MeshDistance CylindricalMesh::distance_to_grid_boundary(
|
|||
double l) const
|
||||
{
|
||||
if (i == 0) {
|
||||
|
||||
return std::min(
|
||||
MeshDistance(ijk[i] + 1, true, find_r_crossing(r0, u, l, ijk[i])),
|
||||
MeshDistance(ijk[i] - 1, false, find_r_crossing(r0, u, l, ijk[i] - 1)));
|
||||
MeshDistance({1, 0, 0}, true, find_r_crossing(r0, u, l, ijk[i])),
|
||||
MeshDistance({-1, 0, 0}, false, find_r_crossing(r0, u, l, ijk[i] - 1)));
|
||||
|
||||
} else if (i == 1) {
|
||||
|
||||
return std::min(MeshDistance(sanitize_phi(ijk[i] + 1), true,
|
||||
find_phi_crossing(r0, u, l, ijk[i])),
|
||||
MeshDistance(sanitize_phi(ijk[i] - 1), false,
|
||||
find_phi_crossing(r0, u, l, ijk[i] - 1)));
|
||||
return std::min(
|
||||
MeshDistance({0, 1, 0}, true, find_phi_crossing(r0, u, l, ijk[i])),
|
||||
MeshDistance({0, -1, 0}, false, find_phi_crossing(r0, u, l, ijk[i] - 1)));
|
||||
|
||||
} else {
|
||||
return find_z_crossing(r0, u, l, ijk[i]);
|
||||
|
|
@ -1991,6 +2015,8 @@ int CylindricalMesh::set_grid()
|
|||
static_cast<int>(grid_[1].size()) - 1,
|
||||
static_cast<int>(grid_[2].size()) - 1};
|
||||
|
||||
axes_labels_ = {"r", "phi", "z"};
|
||||
|
||||
for (const auto& g : grid_) {
|
||||
if (g.size() < 2) {
|
||||
set_errmsg("r-, phi-, and z- grids for cylindrical meshes "
|
||||
|
|
@ -2076,7 +2102,6 @@ SphericalMesh::SphericalMesh(pugi::xml_node node)
|
|||
: PeriodicStructuredMesh {node}
|
||||
{
|
||||
n_dimension_ = 3;
|
||||
|
||||
grid_[0] = get_node_array<double>(node, "r_grid");
|
||||
grid_[1] = get_node_array<double>(node, "theta_grid");
|
||||
grid_[2] = get_node_array<double>(node, "phi_grid");
|
||||
|
|
@ -2090,7 +2115,6 @@ SphericalMesh::SphericalMesh(pugi::xml_node node)
|
|||
SphericalMesh::SphericalMesh(hid_t group) : PeriodicStructuredMesh {group}
|
||||
{
|
||||
n_dimension_ = 3;
|
||||
|
||||
read_dataset(group, "r_grid", grid_[0]);
|
||||
read_dataset(group, "theta_grid", grid_[1]);
|
||||
read_dataset(group, "phi_grid", grid_[2]);
|
||||
|
|
@ -2297,20 +2321,19 @@ StructuredMesh::MeshDistance SphericalMesh::distance_to_grid_boundary(
|
|||
|
||||
if (i == 0) {
|
||||
return std::min(
|
||||
MeshDistance(ijk[i] + 1, true, find_r_crossing(r0, u, l, ijk[i])),
|
||||
MeshDistance(ijk[i] - 1, false, find_r_crossing(r0, u, l, ijk[i] - 1)));
|
||||
MeshDistance({1, 0, 0}, true, find_r_crossing(r0, u, l, ijk[i])),
|
||||
MeshDistance({-1, 0, 0}, false, find_r_crossing(r0, u, l, ijk[i] - 1)));
|
||||
|
||||
} else if (i == 1) {
|
||||
return std::min(MeshDistance(sanitize_theta(ijk[i] + 1), true,
|
||||
find_theta_crossing(r0, u, l, ijk[i])),
|
||||
MeshDistance(sanitize_theta(ijk[i] - 1), false,
|
||||
find_theta_crossing(r0, u, l, ijk[i] - 1)));
|
||||
return std::min(
|
||||
MeshDistance({0, 1, 0}, true, find_theta_crossing(r0, u, l, ijk[i])),
|
||||
MeshDistance(
|
||||
{0, -1, 0}, false, find_theta_crossing(r0, u, l, ijk[i] - 1)));
|
||||
|
||||
} else {
|
||||
return std::min(MeshDistance(sanitize_phi(ijk[i] + 1), true,
|
||||
find_phi_crossing(r0, u, l, ijk[i])),
|
||||
MeshDistance(sanitize_phi(ijk[i] - 1), false,
|
||||
find_phi_crossing(r0, u, l, ijk[i] - 1)));
|
||||
return std::min(
|
||||
MeshDistance({0, 0, 1}, true, find_phi_crossing(r0, u, l, ijk[i])),
|
||||
MeshDistance({0, 0, -1}, false, find_phi_crossing(r0, u, l, ijk[i] - 1)));
|
||||
}
|
||||
}
|
||||
|
||||
|
|
@ -2320,9 +2343,11 @@ int SphericalMesh::set_grid()
|
|||
static_cast<int>(grid_[1].size()) - 1,
|
||||
static_cast<int>(grid_[2].size()) - 1};
|
||||
|
||||
axes_labels_ = {"r", "theta", "phi"};
|
||||
|
||||
for (const auto& g : grid_) {
|
||||
if (g.size() < 2) {
|
||||
set_errmsg("x-, y-, and z- grids for spherical meshes "
|
||||
set_errmsg("r-, theta-, and phi- grids for spherical meshes "
|
||||
"must each have at least 2 points");
|
||||
return OPENMC_E_INVALID_ARGUMENT;
|
||||
}
|
||||
|
|
@ -2398,6 +2423,413 @@ double SphericalMesh::volume(const MeshIndex& ijk) const
|
|||
(std::cos(theta_i) - std::cos(theta_o)) * (phi_o - phi_i);
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// HexagonalMesh implementation
|
||||
//==============================================================================
|
||||
|
||||
HexagonalMesh::HexagonalMesh(pugi::xml_node node)
|
||||
: PeriodicStructuredMesh {node}
|
||||
{
|
||||
grid_ = get_node_array<double>(node, "z_grid");
|
||||
num_rings_ = std::stoi(get_node_value(node, "num_rings"));
|
||||
pitch_ = std::stod(get_node_value(node, "pitch"));
|
||||
origin_ = get_node_position(node, "origin");
|
||||
|
||||
// Read the orientation. Default to 'y'.
|
||||
if (check_for_node(node, "orientation")) {
|
||||
std::string orientation = get_node_value(node, "orientation");
|
||||
if (orientation == "x") {
|
||||
orientation_ = Orientation::x;
|
||||
} else if (orientation != "y") {
|
||||
fatal_error("Unrecognized orientation '" + orientation + "'");
|
||||
}
|
||||
}
|
||||
|
||||
if (int err = set_grid()) {
|
||||
fatal_error(openmc_err_msg);
|
||||
}
|
||||
}
|
||||
|
||||
HexagonalMesh::HexagonalMesh(hid_t group) : PeriodicStructuredMesh {group}
|
||||
{
|
||||
read_dataset(group, "z_grid", grid_);
|
||||
read_attribute(group, "num_rings", num_rings_);
|
||||
read_attribute(group, "pitch", pitch_);
|
||||
read_dataset(group, "origin", origin_);
|
||||
|
||||
if (attribute_exists(group, "orientation")) {
|
||||
std::string orientation;
|
||||
read_attribute(group, "orientation", orientation);
|
||||
if (orientation == "x") {
|
||||
orientation_ = Orientation::x;
|
||||
} else if (orientation != "y") {
|
||||
fatal_error("Unrecognized orientation '" + orientation + "'");
|
||||
}
|
||||
}
|
||||
if (int err = set_grid()) {
|
||||
fatal_error(openmc_err_msg);
|
||||
}
|
||||
}
|
||||
|
||||
const std::string HexagonalMesh::mesh_type = "hexagonal";
|
||||
|
||||
std::string HexagonalMesh::get_mesh_type() const
|
||||
{
|
||||
return mesh_type;
|
||||
}
|
||||
|
||||
std::string HexagonalMesh::bin_label(int bin) const
|
||||
{
|
||||
MeshIndex ijk = get_indices_from_bin(bin);
|
||||
return fmt::format(
|
||||
"Mesh Index ({}, {}, {}, {})", ijk[0], ijk[1], -ijk[0] - ijk[1], ijk[2]);
|
||||
}
|
||||
|
||||
std::string HexagonalMesh::surface_label(int surface) const
|
||||
{
|
||||
int axis = static_cast<int>(std::floor(surface / 4));
|
||||
int max = static_cast<int>(std::floor((surface % 4) / 2));
|
||||
const std::vector<std::string> MINMAX = {"min", "max"};
|
||||
const std::vector<std::string> OUTIN = {"Outgoing", "Incoming"};
|
||||
auto minmax = MINMAX[max];
|
||||
auto outin = OUTIN[surface % 2];
|
||||
std::string label = axes_labels_[axis];
|
||||
if (axis == 3) {
|
||||
return fmt::format("{}, {}-{}", outin, label, minmax);
|
||||
}
|
||||
auto minmax2 = MINMAX[1 - max];
|
||||
std::string label2 = axes_labels_[(axis + 1) % 2];
|
||||
return fmt::format("{}, {}-{} {}-{}", outin, label, minmax, label2, minmax2);
|
||||
}
|
||||
|
||||
int HexagonalMesh::n_bins() const
|
||||
{
|
||||
return (1 + 3 * num_rings_ * (num_rings_ - 1)) * (grid_.size() - 1);
|
||||
}
|
||||
|
||||
bool HexagonalMesh::valid_index(const MeshIndex& ijk, int k) const
|
||||
{
|
||||
if (k == 3)
|
||||
return ((ijk[2] >= 1) && (ijk[2] < grid_.size()));
|
||||
int ring =
|
||||
std::max({std::abs(ijk[0]), std::abs(ijk[1]), std::abs(ijk[0] + ijk[1])});
|
||||
return (ring < num_rings_);
|
||||
}
|
||||
|
||||
int HexagonalMesh::get_bin_from_indices(const MeshIndex& ijk) const
|
||||
{
|
||||
int q = ijk[0];
|
||||
int r = ijk[1];
|
||||
int k = ijk[2];
|
||||
int hexes = 3 * num_rings_ * (num_rings_ - 1) + 1;
|
||||
int bin = (k - 1) * (grid_.size() - 1) * hexes;
|
||||
int rad = std::max({std::abs(r), std::abs(q), std::abs(r + q)});
|
||||
if (rad == 0)
|
||||
return bin;
|
||||
|
||||
bin += 3 * rad * (rad - 1) + 1;
|
||||
|
||||
if ((q >= 0) and (r >= 0)) {
|
||||
bin += q;
|
||||
} else if ((r < 0) && ((-r - q) < 0)) {
|
||||
bin += (rad - r);
|
||||
} else if ((q >= 0) && ((-q - r) >= 0)) {
|
||||
bin += (2 * rad - q - r);
|
||||
} else if ((q < 0) && (r < 0)) {
|
||||
bin += (3 * rad - q);
|
||||
} else if ((r >= 0) && ((-q - r) >= 0)) {
|
||||
bin += (4 * rad + r);
|
||||
} else {
|
||||
bin += (5 * rad + q + r);
|
||||
}
|
||||
return bin;
|
||||
}
|
||||
|
||||
StructuredMesh::MeshIndex HexagonalMesh::get_indices(
|
||||
Position r, bool& in_mesh) const
|
||||
{
|
||||
r = local_coords(r);
|
||||
|
||||
double q = q_dual_.dot(r);
|
||||
double r0 = r_dual_.dot(r);
|
||||
double s = -q - r0;
|
||||
|
||||
double q_r = std::round(q);
|
||||
double r_r = std::round(r0);
|
||||
double s_r = std::round(s);
|
||||
std::array<double, 3> diff = {
|
||||
std::abs(q - q_r), std::abs(r0 - r_r), std::abs(s - s_r)};
|
||||
auto max_it = std::max_element(diff.begin(), diff.end());
|
||||
int i = static_cast<int>(std::distance(diff.begin(), max_it));
|
||||
|
||||
if (i == 0)
|
||||
q_r = -s_r - r_r;
|
||||
else if (i == 1)
|
||||
r_r = -q_r - s_r;
|
||||
|
||||
int z_i = get_index_in_z_direction(r.z);
|
||||
int q_i = static_cast<int>(q_r);
|
||||
int r_i = static_cast<int>(r_r);
|
||||
int s_i = -q_i - r_i;
|
||||
|
||||
MeshIndex ijk;
|
||||
in_mesh = true;
|
||||
|
||||
ijk[0] = q_i;
|
||||
ijk[1] = r_i;
|
||||
ijk[2] = z_i;
|
||||
|
||||
if (std::max({std::abs(q_i), std::abs(r_i), std::abs(s_i)}) >= num_rings_)
|
||||
in_mesh = false;
|
||||
if (ijk[2] < 1 || ijk[2] > grid_.size() - 1)
|
||||
in_mesh = false;
|
||||
|
||||
return ijk;
|
||||
}
|
||||
|
||||
StructuredMesh::MeshIndex HexagonalMesh::get_indices_from_bin(int bin) const
|
||||
{
|
||||
MeshIndex ijk = {0, 0, 0};
|
||||
int hexes = 3 * num_rings_ * (num_rings_ - 1) + 1;
|
||||
ijk[2] = static_cast<int>(std::floor(bin / hexes)) + 1;
|
||||
int sp_idx = bin % hexes;
|
||||
if (sp_idx == 0) {
|
||||
return ijk;
|
||||
}
|
||||
int rad = static_cast<int>(std::floor((3 + std::sqrt(12 * sp_idx - 3)) / 6));
|
||||
int offset = sp_idx - (3 * rad * (rad - 1) + 1);
|
||||
int side = (offset / rad) % 6;
|
||||
int dist = offset % rad;
|
||||
switch (side) {
|
||||
case 0:
|
||||
ijk[0] = dist;
|
||||
ijk[1] = rad - dist;
|
||||
break;
|
||||
case 1:
|
||||
ijk[0] = rad;
|
||||
ijk[1] = -dist;
|
||||
break;
|
||||
case 2:
|
||||
ijk[0] = rad - dist;
|
||||
ijk[1] = -rad;
|
||||
break;
|
||||
case 3:
|
||||
ijk[0] = -dist;
|
||||
ijk[1] = dist - rad;
|
||||
break;
|
||||
case 4:
|
||||
ijk[0] = -rad;
|
||||
ijk[1] = dist;
|
||||
break;
|
||||
case 5:
|
||||
ijk[0] = dist - rad;
|
||||
ijk[1] = rad;
|
||||
break;
|
||||
}
|
||||
return ijk;
|
||||
}
|
||||
|
||||
Position HexagonalMesh::sample_element(
|
||||
const MeshIndex& ijk, uint64_t* seed) const
|
||||
{
|
||||
double radius = pitch_ / std::sqrt(3.0);
|
||||
double q0 = ijk[0];
|
||||
double r0 = ijk[1];
|
||||
double z = uniform_distribution(grid_[ijk[2] - 1], grid_[ijk[2]], seed);
|
||||
double q, r, s;
|
||||
do {
|
||||
double rad = std::cbrt(prn(seed)) * radius;
|
||||
double phi = uniform_distribution(0.0, 2 * PI, seed);
|
||||
double x = rad * std::cos(phi);
|
||||
double y = rad * std::sin(phi);
|
||||
q = q_dual_[0] * x + q_dual_[1] * y;
|
||||
r = r_dual_[0] * x + r_dual_[1] * y;
|
||||
s = -q - r;
|
||||
} while (std::max({std::abs(q), std::abs(r), std::abs(s)}) > 1);
|
||||
|
||||
return origin_ + radius * ((q + q0) * q_ + (r + r0) * r_) + z;
|
||||
}
|
||||
|
||||
StructuredMesh::MeshDistance HexagonalMesh::find_z_crossing(
|
||||
const Position& r, const Direction& u, double l, int shell) const
|
||||
{
|
||||
MeshDistance d;
|
||||
|
||||
// Direction of flight is within xy-plane. Will never intersect z.
|
||||
if (std::abs(u.z) < FP_PRECISION)
|
||||
return d;
|
||||
|
||||
d.max_surface = (u.z > 0.0);
|
||||
if (d.max_surface && (shell < grid_.size())) {
|
||||
++d.offset[2];
|
||||
d.distance = (grid_[shell] - r.z) / u.z;
|
||||
} else if (!d.max_surface && (shell > 0)) {
|
||||
--d.offset[2];
|
||||
d.distance = (grid_[shell - 1] - r.z) / u.z;
|
||||
}
|
||||
return d;
|
||||
}
|
||||
|
||||
StructuredMesh::MeshDistance HexagonalMesh::distance_to_grid_boundary(
|
||||
const MeshIndex& ijk, int i, const Position& r0, const Direction& u,
|
||||
double l) const
|
||||
{
|
||||
auto r = local_coords(r0);
|
||||
if (i == 3)
|
||||
return find_z_crossing(r, u, l, ijk[2]);
|
||||
r.z = 0.0;
|
||||
|
||||
MeshDistance d;
|
||||
Direction s_ = -q_ - r_;
|
||||
Direction dir;
|
||||
MeshIndex idx = {ijk[0], ijk[1], -ijk[0] - ijk[1]};
|
||||
MeshIndex offset = {0, 0, 0};
|
||||
double distance;
|
||||
if (i == 0) {
|
||||
dir = 0.5 * (q_ - r_);
|
||||
++offset[0];
|
||||
--offset[1];
|
||||
} else if (i == 1) {
|
||||
dir = 0.5 * (r_ - s_);
|
||||
++offset[1];
|
||||
--offset[2];
|
||||
} else {
|
||||
dir = 0.5 * (s_ - q_);
|
||||
--offset[0];
|
||||
++offset[2];
|
||||
}
|
||||
dir /= dir.norm();
|
||||
if (std::abs(u.dot(dir)) < FP_PRECISION)
|
||||
return d;
|
||||
|
||||
d.max_surface = (u.dot(dir) > 0.0);
|
||||
if (d.max_surface) {
|
||||
distance = (idx[0] * q_ + idx[1] * r_ + dir - r).norm() / u.dot(dir);
|
||||
idx[0] += offset[0];
|
||||
idx[1] += offset[1];
|
||||
idx[2] += offset[2];
|
||||
} else {
|
||||
distance = -(idx[0] * q_ + idx[1] * r_ - dir - r).norm() / u.dot(dir);
|
||||
idx[0] -= offset[0];
|
||||
idx[1] -= offset[1];
|
||||
idx[2] -= offset[2];
|
||||
}
|
||||
int radius = std::max({std::abs(idx[0]), std::abs(idx[1]), std::abs(idx[2])});
|
||||
if (d.max_surface && (radius < num_rings_)) {
|
||||
d.offset[0] += offset[0];
|
||||
d.offset[1] += offset[1];
|
||||
d.distance = distance;
|
||||
} else if (!d.max_surface && (radius < num_rings_)) {
|
||||
d.offset[0] -= offset[0];
|
||||
d.offset[1] -= offset[1];
|
||||
d.distance = distance;
|
||||
}
|
||||
return d;
|
||||
}
|
||||
|
||||
double HexagonalMesh::distance_to_mesh(const MeshIndex& ijk,
|
||||
const std::array<MeshDistance, 4>& distances, double traveled_distance,
|
||||
int& k_max) const
|
||||
{
|
||||
// For all directions outside the mesh, find the distance that we need
|
||||
// to travel to reach the next surface. Use the largest distance, as
|
||||
// only this will cross all outer surfaces.
|
||||
const int n = n_dimension_;
|
||||
k_max = -1;
|
||||
|
||||
for (int k = 0; k < n; ++k) {
|
||||
if (distances[k].distance <= traveled_distance)
|
||||
continue;
|
||||
if ((k < 2) && (std::abs(ijk[k]) < num_rings_))
|
||||
continue;
|
||||
if ((k == 2) && (std::abs(ijk[0] + ijk[1]) < num_rings_))
|
||||
continue;
|
||||
if ((k == 3) && (ijk[2] >= 1) && (ijk[2] < grid_.size()))
|
||||
continue;
|
||||
traveled_distance = distances[k].distance;
|
||||
k_max = k;
|
||||
}
|
||||
// Assure some distance is traveled
|
||||
if (k_max == -1) {
|
||||
traveled_distance += TINY_BIT;
|
||||
}
|
||||
return traveled_distance;
|
||||
}
|
||||
|
||||
int HexagonalMesh::set_grid()
|
||||
{
|
||||
n_dimension_ = 4;
|
||||
correlated_axes_ = {{0, 1, 2}, {0, 1, 2}, {0, 1, 2}, {3}};
|
||||
|
||||
axes_labels_ = {"r", "q", "s", "z"};
|
||||
|
||||
if (orientation_ == Orientation::x) {
|
||||
q_ = {std::sqrt(3.0), 0.0, 0.0};
|
||||
r_ = {0.5 * std::sqrt(3), -1.5, 0.0};
|
||||
q_dual_ = {std::sqrt(3.0) / 3.0, 1.0 / 3.0, 0.0};
|
||||
r_dual_ = {0.0, -2.0 / 3.0, 0.0};
|
||||
} else {
|
||||
q_ = {1.5, -0.5 * std::sqrt(3.0), 0.0};
|
||||
r_ = {0.0, -std::sqrt(3.0), 0.0};
|
||||
q_dual_ = {2.0 / 3.0, 0.0, 0.0};
|
||||
r_dual_ = {-1.0 / 3.0, -std::sqrt(3.0) / 3.0, 0.0};
|
||||
}
|
||||
|
||||
double size = pitch_ / std::sqrt(3.0);
|
||||
|
||||
q_ *= size;
|
||||
r_ *= size;
|
||||
q_dual_ /= size;
|
||||
r_dual_ /= size;
|
||||
|
||||
if (grid_.size() < 2) {
|
||||
set_errmsg("z- grid for hexagonal meshes "
|
||||
"must have at least 2 points");
|
||||
return OPENMC_E_INVALID_ARGUMENT;
|
||||
}
|
||||
if (std::adjacent_find(grid_.begin(), grid_.end(), std::greater_equal<>()) !=
|
||||
grid_.end()) {
|
||||
set_errmsg("Values in z- grid for "
|
||||
"hexagonal meshes must be sorted and unique.");
|
||||
return OPENMC_E_INVALID_ARGUMENT;
|
||||
}
|
||||
return 0;
|
||||
}
|
||||
|
||||
int HexagonalMesh::get_index_in_z_direction(double z) const
|
||||
{
|
||||
return lower_bound_index(grid_.begin(), grid_.end(), z) + 1;
|
||||
}
|
||||
|
||||
std::pair<vector<double>, vector<double>> HexagonalMesh::plot(
|
||||
Position plot_ll, Position plot_ur) const
|
||||
{
|
||||
fatal_error("Plot of hexagonal Mesh not implemented");
|
||||
|
||||
// Figure out which axes lie in the plane of the plot.
|
||||
array<vector<double>, 2> axis_lines;
|
||||
return {axis_lines[0], axis_lines[1]};
|
||||
}
|
||||
|
||||
void HexagonalMesh::to_hdf5_inner(hid_t mesh_group) const
|
||||
{
|
||||
write_dataset(mesh_group, "z_grid", grid_);
|
||||
if (orientation_ == Orientation::x)
|
||||
write_attribute(mesh_group, "orientation", "x");
|
||||
write_attribute(mesh_group, "num_rings", num_rings_);
|
||||
write_attribute(mesh_group, "pitch", pitch_);
|
||||
write_dataset(mesh_group, "origin", origin_);
|
||||
}
|
||||
|
||||
double HexagonalMesh::volume(const MeshIndex& ijk) const
|
||||
{
|
||||
double f = 1.5 * std::sqrt(3.0);
|
||||
double size = pitch_ / std::sqrt(3.0);
|
||||
double z_i = grid_[ijk[2] - 1];
|
||||
double z_o = grid_[ijk[2]];
|
||||
return f * (z_o - z_i) * size * size;
|
||||
}
|
||||
|
||||
//==============================================================================
|
||||
// Helper functions for the C API
|
||||
//==============================================================================
|
||||
|
|
@ -2464,6 +2896,8 @@ extern "C" int openmc_extend_meshes(
|
|||
model::meshes.push_back(make_unique<CylindricalMesh>());
|
||||
} else if (SphericalMesh::mesh_type == type) {
|
||||
model::meshes.push_back(make_unique<SphericalMesh>());
|
||||
} else if (HexagonalMesh::mesh_type == type) {
|
||||
model::meshes.push_back(make_unique<HexagonalMesh>());
|
||||
} else {
|
||||
throw std::runtime_error {"Unknown mesh type: " + std::string(type)};
|
||||
}
|
||||
|
|
@ -2857,6 +3291,28 @@ extern "C" int openmc_spherical_mesh_set_grid(int32_t index,
|
|||
index, grid_x, nx, grid_y, ny, grid_z, nz);
|
||||
}
|
||||
|
||||
//! Get the hexagonal mesh grid
|
||||
extern "C" int openmc_hexagonal_mesh_get_grid(int32_t index, double** grid_z,
|
||||
int* nz, int* nr, double** origin, double* pitch, const char** orient)
|
||||
{
|
||||
if (int err = check_mesh_type<HexagonalMesh>(index))
|
||||
return err;
|
||||
HexagonalMesh* m = dynamic_cast<HexagonalMesh*>(model::meshes[index].get());
|
||||
|
||||
if (m->grid_.empty()) {
|
||||
set_errmsg("Mesh parameters have not been set.");
|
||||
return OPENMC_E_ALLOCATE;
|
||||
}
|
||||
|
||||
*grid_z = m->grid_.data();
|
||||
*nz = m->grid_.size();
|
||||
*nr = m->num_rings_;
|
||||
*origin = &m->origin_.x;
|
||||
*pitch = m->pitch_;
|
||||
*orient = (m->orientation_ == HexagonalMesh::Orientation::y) ? "y" : "x";
|
||||
return 0;
|
||||
}
|
||||
|
||||
#ifdef OPENMC_DAGMC_ENABLED
|
||||
|
||||
const std::string MOABMesh::mesh_lib_type = "moab";
|
||||
|
|
|
|||
|
|
@ -25,62 +25,24 @@ void MeshSurfaceFilter::get_all_bins(
|
|||
|
||||
std::string MeshSurfaceFilter::text_label(int bin) const
|
||||
{
|
||||
auto& mesh = *model::meshes[mesh_];
|
||||
int n_dim = mesh.n_dimension_;
|
||||
auto mesh = dynamic_cast<StructuredMesh*>(model::meshes[mesh_].get());
|
||||
int n_dim = mesh->n_dimension_;
|
||||
|
||||
// Get flattend mesh index and surface index.
|
||||
int i_mesh = bin / (4 * n_dim);
|
||||
MeshDir surf_dir = static_cast<MeshDir>(bin % (4 * n_dim));
|
||||
int i_surf = bin % (4 * n_dim);
|
||||
|
||||
// Get mesh index part of label.
|
||||
std::string out = MeshFilter::text_label(i_mesh);
|
||||
|
||||
// Get surface part of label.
|
||||
switch (surf_dir) {
|
||||
case MeshDir::OUT_LEFT:
|
||||
out += " Outgoing, x-min";
|
||||
break;
|
||||
case MeshDir::IN_LEFT:
|
||||
out += " Incoming, x-min";
|
||||
break;
|
||||
case MeshDir::OUT_RIGHT:
|
||||
out += " Outgoing, x-max";
|
||||
break;
|
||||
case MeshDir::IN_RIGHT:
|
||||
out += " Incoming, x-max";
|
||||
break;
|
||||
case MeshDir::OUT_BACK:
|
||||
out += " Outgoing, y-min";
|
||||
break;
|
||||
case MeshDir::IN_BACK:
|
||||
out += " Incoming, y-min";
|
||||
break;
|
||||
case MeshDir::OUT_FRONT:
|
||||
out += " Outgoing, y-max";
|
||||
break;
|
||||
case MeshDir::IN_FRONT:
|
||||
out += " Incoming, y-max";
|
||||
break;
|
||||
case MeshDir::OUT_BOTTOM:
|
||||
out += " Outgoing, z-min";
|
||||
break;
|
||||
case MeshDir::IN_BOTTOM:
|
||||
out += " Incoming, z-min";
|
||||
break;
|
||||
case MeshDir::OUT_TOP:
|
||||
out += " Outgoing, z-max";
|
||||
break;
|
||||
case MeshDir::IN_TOP:
|
||||
out += " Incoming, z-max";
|
||||
break;
|
||||
}
|
||||
|
||||
out += " ";
|
||||
out += mesh->surface_label(i_surf);
|
||||
return out;
|
||||
}
|
||||
|
||||
void MeshSurfaceFilter::set_mesh(int32_t mesh)
|
||||
{
|
||||
mesh_ = mesh;
|
||||
if (!dynamic_cast<StructuredMesh*>(model::meshes[mesh_].get()))
|
||||
fatal_error("Only structured mesh is supported in MeshSurfaceFilter.");
|
||||
n_bins_ = model::meshes[mesh_]->n_surface_bins();
|
||||
}
|
||||
|
||||
|
|
|
|||
138
tests/unit_tests/test_hexagonal_mesh.py
Normal file
138
tests/unit_tests/test_hexagonal_mesh.py
Normal file
|
|
@ -0,0 +1,138 @@
|
|||
from itertools import product, permutations
|
||||
|
||||
import openmc
|
||||
import numpy as np
|
||||
|
||||
import pytest
|
||||
|
||||
pitch = 1.25
|
||||
|
||||
@pytest.fixture(params=['x','y'])
|
||||
def model(request):
|
||||
openmc.reset_auto_ids()
|
||||
|
||||
orientation = request.param
|
||||
|
||||
water = openmc.Material(name='water')
|
||||
water.add_element('H', 2.0)
|
||||
water.add_element('O', 1.0)
|
||||
water.set_density('g/cc', 1.0)
|
||||
|
||||
outer = openmc.Cell(fill=water, cell_id=100)
|
||||
cell10 = openmc.Cell(fill=water, cell_id=10)
|
||||
cell00 = openmc.Cell(fill=water, cell_id=0)
|
||||
cell01 = openmc.Cell(fill=water, cell_id=1)
|
||||
cell02 = openmc.Cell(fill=water, cell_id=2)
|
||||
cell03 = openmc.Cell(fill=water, cell_id=3)
|
||||
cell04 = openmc.Cell(fill=water, cell_id=4)
|
||||
cell05 = openmc.Cell(fill=water, cell_id=5)
|
||||
|
||||
univ10 = openmc.Universe(cells=[cell10])
|
||||
univ00 = openmc.Universe(cells=[cell00])
|
||||
univ01 = openmc.Universe(cells=[cell01])
|
||||
univ02 = openmc.Universe(cells=[cell02])
|
||||
univ03 = openmc.Universe(cells=[cell03])
|
||||
univ04 = openmc.Universe(cells=[cell04])
|
||||
univ05 = openmc.Universe(cells=[cell05])
|
||||
|
||||
plane1 = openmc.ZPlane(-10.0, boundary_type='vacuum')
|
||||
plane2 = openmc.ZPlane(10.0, boundary_type='vacuum')
|
||||
|
||||
lat = openmc.HexLattice()
|
||||
lat.center = (0., 0.)
|
||||
lat.pitch = (pitch,)
|
||||
lat.universes = [[univ00, univ01, univ02, univ03, univ04, univ05], [univ10]]
|
||||
lat.outer = openmc.Universe(cells=[outer])
|
||||
lat.orientation = orientation
|
||||
|
||||
hex_prism = openmc.model.HexagonalPrism(
|
||||
edge_length=2*pitch,
|
||||
orientation=orientation,
|
||||
boundary_type='vacuum'
|
||||
)
|
||||
cell = openmc.Cell(region=-hex_prism & +plane1 & -plane2, fill=lat)
|
||||
|
||||
geom = openmc.Geometry([cell])
|
||||
|
||||
source = openmc.IndependentSource()
|
||||
source.space = openmc.stats.Point()
|
||||
source.energy = openmc.stats.Discrete([10000], [1.0])
|
||||
|
||||
settings = openmc.Settings()
|
||||
settings.particles = 2000
|
||||
settings.batches = 10
|
||||
settings.run_mode = 'fixed source'
|
||||
|
||||
# build
|
||||
mesh = openmc.HexagonalMesh(
|
||||
z_grid=[-10.0, 10.0],
|
||||
num_rings = 2,
|
||||
pitch = pitch,
|
||||
orientation = orientation,
|
||||
)
|
||||
tally = openmc.Tally()
|
||||
|
||||
mesh_filter = openmc.MeshFilter(mesh)
|
||||
cell_filter = openmc.CellFilter([cell10, cell00, cell01, cell02, cell03, cell04, cell05])
|
||||
tally.filters.append(mesh_filter)
|
||||
tally.filters.append(cell_filter)
|
||||
|
||||
tally.scores.append("total")
|
||||
|
||||
tallies = openmc.Tallies([tally])
|
||||
|
||||
return openmc.Model(geometry=geom, settings=settings, tallies=tallies)
|
||||
|
||||
|
||||
@pytest.mark.parametrize("estimator", ["collision", "tracklength"])
|
||||
def test_correct_locations(model, run_in_tmpdir, estimator):
|
||||
tally, = model.tallies
|
||||
tally.estimator = estimator
|
||||
model.run(apply_tally_results=True)
|
||||
df = tally.get_pandas_dataframe()
|
||||
df = df[df['mean']>0]
|
||||
df = df.set_index("cell")["mesh 1"]
|
||||
assert len(df) == 7
|
||||
|
||||
if tally.filters[0].mesh.orientation == "y":
|
||||
assert tuple(df.loc[10]) == (0,0,0,0)
|
||||
assert tuple(df.loc[0]) == (0,-1,1,0)
|
||||
assert tuple(df.loc[1]) == (1,-1,0,0)
|
||||
assert tuple(df.loc[2]) == (1,0,-1,0)
|
||||
assert tuple(df.loc[3]) == (0,1,-1,0)
|
||||
assert tuple(df.loc[4]) == (-1,1,0,0)
|
||||
assert tuple(df.loc[5]) == (-1,0,1,0)
|
||||
else:
|
||||
assert tuple(df.loc[10]) == (0,0,0,0)
|
||||
assert tuple(df.loc[0]) == (1,0,-1,0)
|
||||
assert tuple(df.loc[1]) == (0,1,-1,0)
|
||||
assert tuple(df.loc[2]) == (-1,1,0,0)
|
||||
assert tuple(df.loc[3]) == (-1,0,1,0)
|
||||
assert tuple(df.loc[4]) == (0,-1,1,0)
|
||||
assert tuple(df.loc[5]) == (1,-1,0,0)
|
||||
|
||||
def test_meshsurface(model, run_in_tmpdir):
|
||||
tally, = model.tallies
|
||||
orientation = tally.filters[0].mesh.orientation
|
||||
|
||||
tally = openmc.Tally()
|
||||
mesh = openmc.HexagonalMesh(
|
||||
z_grid=[-10.0, 10.0],
|
||||
num_rings = 2,
|
||||
pitch = pitch,
|
||||
orientation = orientation,
|
||||
)
|
||||
|
||||
tally.filters = [openmc.MeshSurfaceFilter(mesh)]
|
||||
tally.scores = ['current']
|
||||
|
||||
model.tallies = [tally]
|
||||
|
||||
model.run(apply_tally_results=True)
|
||||
|
||||
df = tally.get_pandas_dataframe()
|
||||
|
||||
z_max_in = df['mesh 2', 'surf']=='z-max in'
|
||||
z_min_in = df['mesh 2', 'surf']=='z-min in'
|
||||
|
||||
assert np.all(df[z_max_in & z_min_in]['mean']==0.0)
|
||||
Loading…
Add table
Add a link
Reference in a new issue