diff --git a/docs/source/io_formats/tallies.rst b/docs/source/io_formats/tallies.rst index dc57e577511..ea9621d3cb4 100644 --- a/docs/source/io_formats/tallies.rst +++ b/docs/source/io_formats/tallies.rst @@ -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 diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index b3911c8b3ba..023b96d13e6 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -169,6 +169,7 @@ Meshes openmc.CylindricalMesh openmc.SphericalMesh openmc.UnstructuredMesh + openmc.HexagonalMesh Geometry Plotting ----------------- diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 0d8189caa1d..326c16182b7 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -302,10 +302,10 @@ class StructuredMesh : public Mesh { 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 @@ class StructuredMesh : public Mesh { } }; + 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 @@ class StructuredMesh : public Mesh { //! \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 @@ class StructuredMesh : public Mesh { 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& 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 get_shape_tensor() const; @@ -452,6 +468,8 @@ class StructuredMesh : public Mesh { // Data members std::array shape_; //!< Number of mesh elements in each dimension + std::vector> correlated_axes_ = {{0}, {1}, {2}}; + std::vector axes_labels_ = {"x", "y", "z"}; protected: }; @@ -580,6 +598,11 @@ class CylindricalMesh : public PeriodicStructuredMesh { 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 @@ class SphericalMesh : public PeriodicStructuredMesh { 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,85 @@ class SphericalMesh : public PeriodicStructuredMesh { } }; +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& distances, double traveled_distance, + int& k_max) const override; + + std::pair, vector> plot( + Position plot_ll, Position plot_ur) const override; + + void to_hdf5_inner(hid_t group) const override; + + int set_grid(); + + // Data members + int num_rings_; + double pitch_; + vector grid_; + Direction q_; + Direction r_; + Direction q_dual_; + Direction r_dual_; + +private: + enum class Orientation { + y, //!< Flat side of lattice parallel to y-axis + x //!< Flat side of lattice parallel to x-axis + }; + + Orientation orientation_ {Orientation::y}; + + 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 { diff --git a/include/openmc/tallies/filter_meshsurface.h b/include/openmc/tallies/filter_meshsurface.h index 195995c699f..fc6915200cb 100644 --- a/include/openmc/tallies/filter_meshsurface.h +++ b/include/openmc/tallies/filter_meshsurface.h @@ -22,21 +22,6 @@ class MeshSurfaceFilter : public MeshFilter { // 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 diff --git a/openmc/filter.py b/openmc/filter.py index 87aeb70c3a7..754cd753d6f 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -30,13 +30,6 @@ '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 @@ def __new__(cls, name, bases, namespace, **kwargs): 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 @@ def get_pandas_dataframe(self, data_size, stride, **kwargs): # 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 @@ def get_pandas_dataframe(self, data_size, stride, **kwargs): 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. diff --git a/openmc/mesh.py b/openmc/mesh.py index 030a57218eb..7feac1d1954 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -300,6 +300,8 @@ def from_hdf5(cls, group: h5py.Group): 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: @@ -527,11 +529,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): @@ -542,11 +539,6 @@ def n_dimension(self): def _axis_labels(self): pass - @property - @abstractmethod - def _grids(self): - pass - @abstractmethod def get_indices_at_coords(self, coords: Sequence[float]) -> tuple: pass @@ -2533,6 +2525,273 @@ def get_indices_at_coords(self, coords: Sequence[float]) -> tuple: "get_indices_at_coords is not yet implemented for SphericalMesh" ) +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_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+1): + 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) diff --git a/src/mesh.cpp b/src/mesh.cpp index 5ab7ac3988b..5808d7c21e8 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -342,6 +342,8 @@ const std::unique_ptr& Mesh::create( model::meshes.push_back(make_unique(dataset)); } else if (mesh_type == SphericalMesh::mesh_type) { model::meshes.push_back(make_unique(dataset)); + } else if (mesh_type == HexagonalMesh::mesh_type) { + model::meshes.push_back(make_unique(dataset)); #ifdef OPENMC_DAGMC_ENABLED } else if (mesh_type == UnstructuredMesh::mesh_type && mesh_library == MOABMesh::mesh_lib_type) { @@ -767,6 +769,17 @@ std::string StructuredMesh::bin_label(int bin) const } } +std::string StructuredMesh::surface_label(int surface) const +{ + int axis = static_cast(std::floor(surface / 4)); + std::string label = axes_labels_[axis]; + const std::vector MINMAX = {"min", "max"}; + const std::vector OUTIN = {"Outgoing", "Incoming"}; + auto minmax = MINMAX[static_cast(std::floor((surface % 4) / 2))]; + auto outin = OUTIN[surface % 2]; + return fmt::format("{}, {}-{}", outin, minmax, label); +} + tensor::Tensor StructuredMesh::get_shape_tensor() const { return tensor::Tensor(shape_.data(), static_cast(n_dimension_)); @@ -1036,6 +1049,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_) { @@ -1183,14 +1201,13 @@ void StructuredMesh::raytrace_mesh( } // Calculate initial distances to next surfaces in all three dimensions - std::array distances; + std::array 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 @@ -1211,14 +1228,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 @@ -1226,23 +1246,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; @@ -1262,6 +1268,29 @@ void StructuredMesh::raytrace_mesh( } } +double StructuredMesh::distance_to_mesh(const MeshIndex& ijk, + const std::array& 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& bins, vector& lengths) const { @@ -1497,16 +1526,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]; } @@ -1670,16 +1698,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; @@ -1932,7 +1959,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) @@ -1940,10 +1966,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; @@ -1954,17 +1980,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]); @@ -1977,6 +2001,8 @@ int CylindricalMesh::set_grid() static_cast(grid_[1].size()) - 1, static_cast(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 " @@ -2062,7 +2088,6 @@ SphericalMesh::SphericalMesh(pugi::xml_node node) : PeriodicStructuredMesh {node} { n_dimension_ = 3; - grid_[0] = get_node_array(node, "r_grid"); grid_[1] = get_node_array(node, "theta_grid"); grid_[2] = get_node_array(node, "phi_grid"); @@ -2076,7 +2101,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]); @@ -2283,20 +2307,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))); } } @@ -2306,9 +2329,11 @@ int SphericalMesh::set_grid() static_cast(grid_[1].size()) - 1, static_cast(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; } @@ -2384,6 +2409,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(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(std::floor(surface / 4)); + int max = static_cast(std::floor((surface % 4) / 2)); + const std::vector MINMAX = {"min", "max"}; + const std::vector 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_ + 1) * num_rings_) * (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 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(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(q_r); + int r_i = static_cast(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(std::floor(bin / hexes)) + 1; + int sp_idx = bin % hexes; + if (sp_idx == 0) { + return ijk; + } + int rad = static_cast(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& 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> 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, 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 //============================================================================== @@ -2450,6 +2882,8 @@ extern "C" int openmc_extend_meshes( model::meshes.push_back(make_unique()); } else if (SphericalMesh::mesh_type == type) { model::meshes.push_back(make_unique()); + } else if (HexagonalMesh::mesh_type == type) { + model::meshes.push_back(make_unique()); } else { throw std::runtime_error {"Unknown mesh type: " + std::string(type)}; } diff --git a/src/tallies/filter_meshsurface.cpp b/src/tallies/filter_meshsurface.cpp index b26cd198b32..6012d441b44 100644 --- a/src/tallies/filter_meshsurface.cpp +++ b/src/tallies/filter_meshsurface.cpp @@ -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(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(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(model::meshes[mesh_].get())) + fatal_error("Only structured mesh is supported in MeshSurfaceFilter."); n_bins_ = model::meshes[mesh_]->n_surface_bins(); } diff --git a/tests/unit_tests/test_hexagonal_mesh.py b/tests/unit_tests/test_hexagonal_mesh.py new file mode 100644 index 00000000000..d76816e081e --- /dev/null +++ b/tests/unit_tests/test_hexagonal_mesh.py @@ -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)