diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index d34ad2087b..c6a870f8d5 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -134,6 +134,18 @@ public: StructuredMesh(pugi::xml_node node) : Mesh {node} {}; virtual ~StructuredMesh() = default; + using MeshIndex = std::array; + struct MeshDistance { + MeshDistance() {}; + MeshDistance(int _index, bool _maxSurface, double _distance): nextIndex{_index}, maxSurface(_maxSurface), distance{_distance} {}; + int nextIndex { -1 }; + double distance { INFTY }; + bool maxSurface { true }; + bool operator<(const MeshDistance& o) const { + return distance < o.distance; + } + }; + int get_bin(Position r) const override; int n_bins() const override; @@ -143,6 +155,18 @@ public: void bins_crossed(Position r0, Position r1, const Direction& u, vector& bins, vector& lengths) const override; + void surface_bins_crossed(Position r0, Position r1, const Direction& u, + vector& bins) const override; + + //! Determine the which cell or surface bins were crossed by a particle + // + //! \param[in] r0 Previous position of the particle + //! \param[in] r1 Current position of the particle + //! \param[in] u Particle direction + //! \param[in] tally Functor that eventually stores the tally data + template void raytrace_mesh(Position r0, Position r1, const Direction& u, + T tally) const; + //! Count number of bank sites in each mesh bin / energy bin // //! \param[in] Pointer to bank sites @@ -155,20 +179,20 @@ public: // //! \param[in] Array of mesh indices //! \return Mesh bin - virtual int get_bin_from_indices(const int* ijk) const; + virtual int get_bin_from_indices(const MeshIndex& ijk) const; //! Get mesh indices given a position // //! \param[in] r Position to get indices for - //! \param[out] ijk Array of mesh indices //! \param[out] in_mesh Whether position is in mesh - virtual void get_indices(Position r, int* ijk, bool* in_mesh) const; + //! \return ijk Array of mesh indices + virtual MeshIndex get_indices(Position r, bool& in_mesh) const; //! Get mesh indices corresponding to a mesh bin // //! \param[in] bin Mesh bin //! \param[out] ijk Mesh indices - virtual void get_indices_from_bin(int bin, int* ijk) const; + virtual MeshIndex get_indices_from_bin(int bin) const; //! Get mesh index in a particular direction //! @@ -176,25 +200,11 @@ public: //! \param[in] i Direction index virtual int get_index_in_direction(double r, int i) const = 0; - //! Check where a line segment intersects the mesh and if it intersects at all - // - //! \param[in,out] r0 In: starting position, out: intersection point - //! \param[in] r1 Ending position - //! \param[out] ijk Indices of the mesh bin containing the intersection point - //! \return Whether the line segment connecting r0 and r1 intersects mesh - virtual bool intersects(Position& r0, Position r1, int* ijk) const; - - //! Get the coordinate for the mesh grid boundary in the positive direction + //! Get the closest distance from to coordinate r to the grid for the mesh grid boundary in the negative direction //! //! \param[in] ijk Array of mesh indices - //! \param[in] i Direction index - virtual double positive_grid_boundary(int* ijk, int i) const = 0; - - //! Get the coordinate for the mesh grid boundary in the negative direction - //! - //! \param[in] ijk Array of mesh indices - //! \param[in] i Direction index - virtual double negative_grid_boundary(int* ijk, int i) const = 0; + //! \param[in] i Direction index + virtual MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const = 0; //! Get a label for the mesh bin std::string bin_label(int bin) const override; @@ -205,9 +215,7 @@ public: xt::xtensor shape_; //!< Number of mesh elements in each dimension protected: - virtual bool intersects_1d(Position& r0, Position r1, int* ijk) const; - virtual bool intersects_2d(Position& r0, Position r1, int* ijk) const; - virtual bool intersects_3d(Position& r0, Position r1, int* ijk) const; + }; //============================================================================== @@ -221,14 +229,11 @@ public: RegularMesh(pugi::xml_node node); // Overridden methods - void surface_bins_crossed(Position r0, Position r1, const Direction& u, - vector& bins) const override; - int get_index_in_direction(double r, int i) const override; - double positive_grid_boundary(int* ijk, int i) const override; - double negative_grid_boundary(int* ijk, int i) const override; + + MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override; std::pair, vector> plot( Position plot_ll, Position plot_ur) const override; @@ -236,6 +241,17 @@ public: void to_hdf5(hid_t group) const override; // New methods + //! Get the coordinate for the mesh grid boundary in the positive direction + //! + //! \param[in] ijk Array of mesh indices + //! \param[in] i Direction index + double positive_grid_boundary(const MeshIndex& ijk, int i) const; + + //! Get the coordinate for the mesh grid boundary in the negative direction + //! + //! \param[in] ijk Array of mesh indices + //! \param[in] i Direction index + double negative_grid_boundary(const MeshIndex& ijk, int i) const; //! Count number of bank sites in each mesh bin / energy bin // @@ -257,14 +273,46 @@ public: RectilinearMesh(pugi::xml_node node); // Overridden methods - void surface_bins_crossed(Position r0, Position r1, const Direction& u, - vector& bins) const override; + int get_index_in_direction(double r, int i) const override; + + MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override; + + std::pair, vector> plot( + Position plot_ll, Position plot_ur) const override; + + void to_hdf5(hid_t group) const override; + + // New methods + //! Get the coordinate for the mesh grid boundary in the positive direction + //! + //! \param[in] ijk Array of mesh indices + //! \param[in] i Direction index + double positive_grid_boundary(const MeshIndex& ijk, int i) const; + + //! Get the coordinate for the mesh grid boundary in the negative direction + //! + //! \param[in] ijk Array of mesh indices + //! \param[in] i Direction index + double negative_grid_boundary(const MeshIndex& ijk, int i) const; + + + vector> grid_; + + int set_grid(); +}; + +class CylindricalMesh : public StructuredMesh { +public: + // Constructors + CylindricalMesh() = default; + CylindricalMesh(pugi::xml_node node); + + // Overridden methods + virtual MeshIndex get_indices(Position r, bool& in_mesh) const override; int get_index_in_direction(double r, int i) const override; - double positive_grid_boundary(int* ijk, int i) const override; - - double negative_grid_boundary(int* ijk, int i) const override; + MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override; std::pair, vector> plot( Position plot_ll, Position plot_ur) const override; @@ -274,6 +322,78 @@ public: vector> grid_; int set_grid(); + +protected: + double find_r_crossing(const Position& r, const Direction& u, double l, int shell) const; + double find_phi_crossing(const Position& r, const Direction& u, double l, int shell) const; + StructuredMesh::MeshDistance find_z_crossing(const Position& r, const Direction& u, double l, int shell) const; + + bool full_phi { false }; + + constexpr inline int sanitize_angular_index(int idx, bool full, int N) const { + if ((idx > 0) and (idx <= N)) { + return idx; + } else if (full) { + return (idx + N - 1) % N + 1; + } else { + return 0; + } + } + + inline int sanitize_phi(int idx) const { + return sanitize_angular_index(idx, full_phi, shape_[1]); + } + +}; + + +class SphericalMesh : public StructuredMesh { +public: + // Constructors + SphericalMesh() = default; + SphericalMesh(pugi::xml_node node); + + // Overridden methods + virtual MeshIndex get_indices(Position r, bool& in_mesh) const override; + + int get_index_in_direction(double r, int i) const override; + + MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override; + + std::pair, vector> plot( + Position plot_ll, Position plot_ur) const override; + + void to_hdf5(hid_t group) const override; + + vector> grid_; + + int set_grid(); + +protected: + double find_r_crossing(const Position& r, const Direction& u, double l, int shell) const; + double find_theta_crossing(const Position& r, const Direction& u, double l, int shell) const; + double find_phi_crossing(const Position& r, const Direction& u, double l, int shell) const; + + bool full_theta { false }; + bool full_phi { false }; + + constexpr inline int sanitize_angular_index(int idx, bool full, int N) const { + if ((idx > 0) and (idx <= N)) { + return idx; + } else if (full) { + return (idx + N - 1) % N + 1; + } else { + return 0; + } + } + + inline int sanitize_theta(int idx) const { + return sanitize_angular_index(idx, full_theta, shape_[1]); + } + inline int sanitize_phi(int idx) const { + return sanitize_angular_index(idx, full_phi, shape_[2]); + } + }; // Abstract class for unstructured meshes diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py index dc528573ac..82545e2bf1 100644 --- a/openmc/lib/mesh.py +++ b/openmc/lib/mesh.py @@ -11,7 +11,7 @@ from . import _dll from .core import _FortranObjectWithID from .error import _error_handler -__all__ = ['RegularMesh', 'RectilinearMesh', 'meshes'] +__all__ = ['RegularMesh', 'RectilinearMesh', 'CylindricalMesh', 'SphericalMesh', 'meshes'] # Mesh functions _dll.openmc_extend_meshes.argtypes = [c_int32, c_char_p, POINTER(c_int32), @@ -56,6 +56,26 @@ _dll.openmc_regular_mesh_set_params.argtypes = [ _dll.openmc_regular_mesh_set_params.restype = c_int _dll.openmc_regular_mesh_set_params.errcheck = _error_handler +_dll.openmc_cylindrical_mesh_get_grid.argtypes = [c_int32, + POINTER(POINTER(c_double)), POINTER(c_int), POINTER(POINTER(c_double)), + POINTER(c_int), POINTER(POINTER(c_double)), POINTER(c_int)] +_dll.openmc_cylindrical_mesh_get_grid.restype = c_int +_dll.openmc_cylindrical_mesh_get_grid.errcheck = _error_handler +_dll.openmc_cylindrical_mesh_set_grid.argtypes = [c_int32, POINTER(c_double), + c_int, POINTER(c_double), c_int, POINTER(c_double), c_int] +_dll.openmc_cylindrical_mesh_set_grid.restype = c_int +_dll.openmc_cylindrical_mesh_set_grid.errcheck = _error_handler + +_dll.openmc_spherical_mesh_get_grid.argtypes = [c_int32, + POINTER(POINTER(c_double)), POINTER(c_int), POINTER(POINTER(c_double)), + POINTER(c_int), POINTER(POINTER(c_double)), POINTER(c_int)] +_dll.openmc_spherical_mesh_get_grid.restype = c_int +_dll.openmc_spherical_mesh_get_grid.errcheck = _error_handler +_dll.openmc_spherical_mesh_set_grid.argtypes = [c_int32, POINTER(c_double), + c_int, POINTER(c_double), c_int, POINTER(c_double), c_int] +_dll.openmc_spherical_mesh_set_grid.restype = c_int +_dll.openmc_spherical_mesh_set_grid.errcheck = _error_handler + class Mesh(_FortranObjectWithID): """Base class to represent mesh objects @@ -275,9 +295,184 @@ class RectilinearMesh(Mesh): ny, z_grid, nz) +class CylindricalMesh(Mesh): + """CylindricalMesh stored internally. + + This class exposes a mesh that is stored internally in the OpenMC + library. To obtain a view of a mesh with a given ID, use the + :data:`openmc.lib.meshes` mapping. + + Parameters + ---------- + index : int + Index in the `meshes` array. + + Attributes + ---------- + id : int + ID of the mesh + dimension : iterable of int + The number of mesh cells in each direction. + lower_left : numpy.ndarray + The lower-left corner of the structured mesh. + upper_right : numpy.ndarray + The upper-right corner of the structrued mesh. + width : numpy.ndarray + The width of mesh cells in each direction. + + """ + mesh_type = 'cylindrical' + + def __init__(self, uid=None, new=True, index=None): + super().__init__(uid, new, index) + + @property + def lower_left(self): + return self._get_parameters()[0] + + @property + def upper_right(self): + return self._get_parameters()[1] + + @property + def dimension(self): + return self._get_parameters()[2] + + @property + def width(self): + return self._get_parameters()[3] + + def _get_parameters(self): + gx = POINTER(c_double)() + nx = c_int() + gy = POINTER(c_double)() + ny = c_int() + gz = POINTER(c_double)() + nz = c_int() + # Call C API to get grid parameters + _dll.openmc_cylindrical_mesh_get_grid(self._index, gx, nx, gy, ny, gz, + nz) + + # Convert grid parameters to Numpy arrays + grid_x = as_array(gx, (nx.value,)) + grid_y = as_array(gy, (ny.value,)) + grid_z = as_array(gz, (nz.value,)) + + # Calculate lower_left, upper_right, width, and dimension from grid + lower_left = np.array((grid_x[0], grid_y[0], grid_z[0])) + upper_right = np.array((grid_x[-1], grid_y[-1], grid_z[-1])) + dimension = np.array((nx.value - 1, ny.value - 1, nz.value - 1)) + width = np.zeros(list(dimension) + [3]) + + for i, diff_x in enumerate(np.diff(grid_x)): + for j, diff_y in enumerate(np.diff(grid_y)): + for k, diff_z in enumerate(np.diff(grid_z)): + width[i, j, k, :] = diff_x, diff_y, diff_z + + return (lower_left, upper_right, dimension, width) + + def set_grid(self, x_grid, y_grid, z_grid): + nx = len(x_grid) + x_grid = (c_double*nx)(*x_grid) + ny = len(y_grid) + y_grid = (c_double*ny)(*y_grid) + nz = len(z_grid) + z_grid = (c_double*nz)(*z_grid) + _dll.openmc_cylindrical_mesh_set_grid(self._index, x_grid, nx, y_grid, + ny, z_grid, nz) + +class SphericalMesh(Mesh): + """SphericalMesh stored internally. + + This class exposes a mesh that is stored internally in the OpenMC + library. To obtain a view of a mesh with a given ID, use the + :data:`openmc.lib.meshes` mapping. + + Parameters + ---------- + index : int + Index in the `meshes` array. + + Attributes + ---------- + id : int + ID of the mesh + dimension : iterable of int + The number of mesh cells in each direction. + lower_left : numpy.ndarray + The lower-left corner of the structured mesh. + upper_right : numpy.ndarray + The upper-right corner of the structrued mesh. + width : numpy.ndarray + The width of mesh cells in each direction. + + """ + mesh_type = 'spherical' + + def __init__(self, uid=None, new=True, index=None): + super().__init__(uid, new, index) + + @property + def lower_left(self): + return self._get_parameters()[0] + + @property + def upper_right(self): + return self._get_parameters()[1] + + @property + def dimension(self): + return self._get_parameters()[2] + + @property + def width(self): + return self._get_parameters()[3] + + def _get_parameters(self): + gx = POINTER(c_double)() + nx = c_int() + gy = POINTER(c_double)() + ny = c_int() + gz = POINTER(c_double)() + nz = c_int() + # Call C API to get grid parameters + _dll.openmc_spherical_mesh_get_grid(self._index, gx, nx, gy, ny, gz, + nz) + + # Convert grid parameters to Numpy arrays + grid_x = as_array(gx, (nx.value,)) + grid_y = as_array(gy, (ny.value,)) + grid_z = as_array(gz, (nz.value,)) + + # Calculate lower_left, upper_right, width, and dimension from grid + lower_left = np.array((grid_x[0], grid_y[0], grid_z[0])) + upper_right = np.array((grid_x[-1], grid_y[-1], grid_z[-1])) + dimension = np.array((nx.value - 1, ny.value - 1, nz.value - 1)) + width = np.zeros(list(dimension) + [3]) + + for i, diff_x in enumerate(np.diff(grid_x)): + for j, diff_y in enumerate(np.diff(grid_y)): + for k, diff_z in enumerate(np.diff(grid_z)): + width[i, j, k, :] = diff_x, diff_y, diff_z + + return (lower_left, upper_right, dimension, width) + + def set_grid(self, x_grid, y_grid, z_grid): + nx = len(x_grid) + x_grid = (c_double*nx)(*x_grid) + ny = len(y_grid) + y_grid = (c_double*ny)(*y_grid) + nz = len(z_grid) + z_grid = (c_double*nz)(*z_grid) + _dll.openmc_spherical_mesh_set_grid(self._index, x_grid, nx, y_grid, + ny, z_grid, nz) + + _MESH_TYPE_MAP = { 'regular': RegularMesh, - 'rectilinear': RectilinearMesh + 'rectilinear': RectilinearMesh, + 'cylindrical': RectilinearMesh, + 'spherical': RectilinearMesh } diff --git a/openmc/mesh.py b/openmc/mesh.py index 9aacdac283..0d4b0bab91 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -79,6 +79,10 @@ class MeshBase(IDManagerMixin, ABC): return RegularMesh.from_hdf5(group) elif mesh_type == 'rectilinear': return RectilinearMesh.from_hdf5(group) + elif mesh_type == 'cylindrical': + return CylindricalMesh.from_hdf5(group) + elif mesh_type == 'spherical': + return SphericalMesh.from_hdf5(group) elif mesh_type == 'unstructured': return UnstructuredMesh.from_hdf5(group) else: @@ -624,6 +628,339 @@ class RectilinearMesh(MeshBase): return element +class CylindricalMesh(MeshBase): + """A 3D cylindrical mesh + + Parameters + ---------- + mesh_id : int + Unique identifier for the mesh + name : str + Name of the mesh + + Attributes + ---------- + id : int + Unique identifier for the mesh + name : str + Name of the mesh + dimension : Iterable of int + The number of mesh cells in each direction. + n_dimension : int + Number of mesh dimensions (always 3 for a CylindricalMesh). + r_grid : Iterable of float + Mesh boundary points along the r-axis. + phi_grid : Iterable of float + Mesh boundary points along the phi-axis. + z_grid : Iterable of float + Mesh boundary points along the z-axis. + indices : Iterable of tuple + An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1), + (2, 1, 1), ...] + + """ + + def __init__(self, mesh_id=None, name=''): + super().__init__(mesh_id, name) + + self._r_grid = None + self._phi_grid = [0, 360] + self._z_grid = None + + @property + def dimension(self): + return (len(self.r_grid) - 1, + len(self.phi_grid) - 1, + len(self.z_grid) - 1) + + @property + def n_dimension(self): + return 3 + + @property + def r_grid(self): + return self._r_grid + + @property + def phi_grid(self): + return self._phi_grid + + @property + def z_grid(self): + return self._z_grid + + @property + def indices(self): + nr = len(self.r_grid) - 1 + np = len(self.phi_grid) - 1 + nz = len(self.z_grid) - 1 + return ((r, p, z) + for z in range(1, nz + 1) + for p in range(1, np + 1) + for r in range(1, nr + 1)) + + @r_grid.setter + def r_grid(self, grid): + cv.check_type('mesh r_grid', grid, Iterable, Real) + self._r_grid = grid + + @phi_grid.setter + def phi_grid(self, grid): + cv.check_type('mesh phi_grid', grid, Iterable, Real) + self._phi_grid = np.array(grid) + + @z_grid.setter + def z_grid(self, grid): + cv.check_type('mesh z_grid', grid, Iterable, Real) + self._z_grid = np.array(grid) + + def __repr__(self): + fmt = '{0: <16}{1}{2}\n' + string = super().__repr__() + string += fmt.format('\tDimensions', '=\t', self.n_dimension) + r_grid_str = str(self._r_grid) if self._r_grid is None else len(self._r_grid) + string += fmt.format('\tN R pnts:', '=\t', r_grid_str) + if self._r_grid is not None: + string += fmt.format('\tR Min:', '=\t', self._r_grid[0]) + string += fmt.format('\tR Max:', '=\t', self._r_grid[-1]) + t_grid_str = str(self._phi_grid) if self._phi_grid is None else len(self._phi_grid) + string += fmt.format('\tN Phi pnts:', '=\t', t_grid_str) + if self._phi_grid is not None: + string += fmt.format('\tPhi Min:', '=\t', self._phi_grid[0]) + string += fmt.format('\tPhi Max:', '=\t', self._phi_grid[-1]) + p_grid_str = str(self._z_grid) if self._z_grid is None else len(self._z_grid) + string += fmt.format('\tN Z pnts:', '=\t', p_grid_str) + if self._z_grid is not None: + string += fmt.format('\tZ Min:', '=\t', self._z_grid[0]) + string += fmt.format('\tZ Max:', '=\t', self._z_grid[-1]) + return string + + @classmethod + def from_hdf5(cls, group): + mesh_id = int(group.name.split('/')[-1].lstrip('mesh ')) + + # Read and assign mesh properties + mesh = cls(mesh_id) + mesh.r_grid = group['r_grid'][()] + mesh.phi_grid = 180 / np.pi * group['p_grid'][()] + mesh.z_grid = group['z_grid'][()] + + return mesh + + def to_xml_element(self): + """Return XML representation of the mesh + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing mesh data + + """ + + element = ET.Element("mesh") + element.set("id", str(self._id)) + element.set("type", "cylindrical") + + subelement = ET.SubElement(element, "r_grid") + subelement.text = ' '.join(map(str, self.r_grid)) + + subelement = ET.SubElement(element, "p_grid") + subelement.text = ' '.join(map(str, self.phi_grid)) + + subelement = ET.SubElement(element, "z_grid") + subelement.text = ' '.join(map(str, self.z_grid)) + + return element + + def calc_mesh_volumes(self): + """Return Volumes for every mesh cell + + Returns + ------- + volumes : Iterable of float + Volumes + + """ + + # V = int_r int_phi int_phi r dr dphi dz + + V_r = np.array(self.r_grid)**2 / 3 + V_r = V_r[1:] - V_r[:-1] + V_p = -np.cos(np.pi * np.array(self.phi_grid) / 180.0) + V_p = V_p[1:] - V_p[:-1] + V_z = np.array(self.z_grid) * np.pi / 180 + V_z = V_z[1:] - V_z[:-1] + + return np.multiply.outer(np.outer(V_r, V_p), V_z) + + +class SphericalMesh(MeshBase): + """A 3D spherical mesh + + Parameters + ---------- + mesh_id : int + Unique identifier for the mesh + name : str + Name of the mesh + + Attributes + ---------- + id : int + Unique identifier for the mesh + name : str + Name of the mesh + dimension : Iterable of int + The number of mesh cells in each direction. + n_dimension : int + Number of mesh dimensions (always 3 for a RectilinearMesh). + r_grid : Iterable of float + Mesh boundary points along the r-axis. + theta_grid : Iterable of float + Mesh boundary points along the theta-axis. + phi_grid : Iterable of float + Mesh boundary points along the phi-axis. + indices : Iterable of tuple + An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1), + (2, 1, 1), ...] + + """ + + def __init__(self, mesh_id=None, name=''): + super().__init__(mesh_id, name) + + self._r_grid = None + self._theta_grid = [0, 180] + self._phi_grid = [0, 360] + + @property + def dimension(self): + return (len(self.r_grid) - 1, + len(self.theta_grid) - 1, + len(self.phi_grid) - 1) + + @property + def n_dimension(self): + return 3 + + @property + def r_grid(self): + return self._r_grid + + @property + def theta_grid(self): + return self._theta_grid + + @property + def phi_grid(self): + return self._phi_grid + + @property + def indices(self): + nr = len(self.r_grid) - 1 + nt = len(self.theta_grid) - 1 + np = len(self.phi_grid) - 1 + return ((r, t, p) + for p in range(1, np + 1) + for t in range(1, nt + 1) + for r in range(1, nr + 1)) + + @r_grid.setter + def r_grid(self, grid): + cv.check_type('mesh r_grid', grid, Iterable, Real) + self._r_grid = grid + + @theta_grid.setter + def theta_grid(self, grid): + cv.check_type('mesh theta_grid', grid, Iterable, Real) + self._theta_grid = np.array(grid) + + @phi_grid.setter + def phi_grid(self, grid): + cv.check_type('mesh phi_grid', grid, Iterable, Real) + self._phi_grid = np.array(grid) + + def __repr__(self): + fmt = '{0: <16}{1}{2}\n' + string = super().__repr__() + string += fmt.format('\tDimensions', '=\t', self.n_dimension) + r_grid_str = str(self._r_grid) if not self._r_grid else len(self._r_grid) + string += fmt.format('\tN R pnts:', '=\t', r_grid_str) + if self._r_grid: + string += fmt.format('\tR Min:', '=\t', self._r_grid[0]) + string += fmt.format('\tR Max:', '=\t', self._r_grid[-1]) + t_grid_str = str(self._theta_grid) if not self._theta_grid else len(self._theta_grid) + string += fmt.format('\tN Theta pnts:', '=\t', t_grid_str) + if self._theta_grid: + string += fmt.format('\tTheta Min:', '=\t', self._theta_grid[0]) + string += fmt.format('\tTheta Max:', '=\t', self._theta_grid[-1]) + p_grid_str = str(self._phi_grid) if not self._phi_grid else len(self._phi_grid) + string += fmt.format('\tN Phi pnts:', '=\t', p_grid_str) + if self._phi_grid: + string += fmt.format('\tPhi Min:', '=\t', self._phi_grid[0]) + string += fmt.format('\tPhi Max:', '=\t', self._phi_grid[-1]) + return string + + @classmethod + def from_hdf5(cls, group): + mesh_id = int(group.name.split('/')[-1].lstrip('mesh ')) + + # Read and assign mesh properties + mesh = cls(mesh_id) + mesh.r_grid = group['r_grid'][()] + mesh.theta_grid = 180 / np.pi * group['t_grid'][()] + mesh.phi_grid = 180 / np.pi * group['p_grid'][()] + + return mesh + + def to_xml_element(self): + """Return XML representation of the mesh + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing mesh data + + """ + + element = ET.Element("mesh") + element.set("id", str(self._id)) + element.set("type", "spherical") + + subelement = ET.SubElement(element, "r_grid") + subelement.text = ' '.join(map(str, self.r_grid)) + + subelement = ET.SubElement(element, "t_grid") + subelement.text = ' '.join(map(str, self.theta_grid)) + + subelement = ET.SubElement(element, "p_grid") + subelement.text = ' '.join(map(str, self.phi_grid)) + + return element + + def calc_mesh_volumes(self): + """Return Volumes for every mesh cell + + Returns + ------- + volumes : Iterable of float + Volumes + + """ + + # V = int_r int_theta int_phi r^2 dr sin(theta) dtheta dphi + + V_r = np.array(self.r_grid)**3 / 3 + V_r = V_r[1:] - V_r[:-1] + V_t = -np.cos(np.pi * np.array(self.theta_grid) / 180.0) + V_t = V_t[1:] - V_t[:-1] + V_p = np.array(self.phi_grid) * np.pi / 180 + V_p = V_p[1:] - V_p[:-1] + + return np.multiply.outer(np.outer(V_r, V_t), V_p) + + + class UnstructuredMesh(MeshBase): """A 3D unstructured mesh diff --git a/src/mesh.cpp b/src/mesh.cpp index cabdd6b3e3..d4948fb864 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -32,6 +32,8 @@ #include "openmc/tallies/tally.h" #include "openmc/xml_interface.h" +#include + #ifdef LIBMESH #include "libmesh/mesh_tools.h" #include "libmesh/numeric_vector.h" @@ -145,8 +147,7 @@ void Mesh::set_id(int32_t id) std::string StructuredMesh::bin_label(int bin) const { - vector ijk(n_dimension_); - get_indices_from_bin(bin, ijk.data()); + MeshIndex ijk = get_indices_from_bin(bin); if (n_dimension_ > 2) { return fmt::format("Mesh Index ({}, {}, {})", ijk[0], ijk[1], ijk[2]); @@ -240,18 +241,20 @@ void UnstructuredMesh::set_length_multiplier(double length_multiplier) specified_length_multiplier_ = true; } -void StructuredMesh::get_indices(Position r, int* ijk, bool* in_mesh) const +StructuredMesh::MeshIndex StructuredMesh::get_indices(Position r, bool& in_mesh) const { - *in_mesh = true; + MeshIndex ijk; + in_mesh = true; for (int i = 0; i < n_dimension_; ++i) { ijk[i] = get_index_in_direction(r[i], i); if (ijk[i] < 1 || ijk[i] > shape_[i]) - *in_mesh = false; + in_mesh = false; } + return ijk; } -int StructuredMesh::get_bin_from_indices(const int* ijk) const +int StructuredMesh::get_bin_from_indices(const MeshIndex& ijk) const { switch (n_dimension_) { case 1: @@ -265,8 +268,9 @@ int StructuredMesh::get_bin_from_indices(const int* ijk) const } } -void StructuredMesh::get_indices_from_bin(int bin, int* ijk) const +StructuredMesh::MeshIndex StructuredMesh::get_indices_from_bin(int bin) const { + MeshIndex ijk; if (n_dimension_ == 1) { ijk[0] = bin + 1; } else if (n_dimension_ == 2) { @@ -277,19 +281,19 @@ void StructuredMesh::get_indices_from_bin(int bin, int* ijk) const ijk[1] = (bin % (shape_[0] * shape_[1])) / shape_[0] + 1; ijk[2] = bin / (shape_[0] * shape_[1]) + 1; } + return ijk; } int StructuredMesh::get_bin(Position r) const { // Determine indices - vector ijk(n_dimension_); bool in_mesh; - get_indices(r, ijk.data(), &in_mesh); + MeshIndex ijk = get_indices(r, in_mesh); if (!in_mesh) return -1; // Convert indices to bin - return get_bin_from_indices(ijk.data()); + return get_bin_from_indices(ijk); } int StructuredMesh::n_bins() const @@ -357,360 +361,158 @@ xt::xtensor StructuredMesh::count_sites( return counts; } -bool StructuredMesh::intersects(Position& r0, Position r1, int* ijk) const +// raytrace through the mesh. The template class T will do the tallying. +// A modern optimizing compiler can recognize the noop method of T and eleminate that call entirely. +template +void StructuredMesh::raytrace_mesh(Position r0, Position r1, const Direction& u, + T tally) const { - switch (n_dimension_) { - case 1: - return intersects_1d(r0, r1, ijk); - case 2: - return intersects_2d(r0, r1, ijk); - case 3: - return intersects_3d(r0, r1, ijk); - default: - throw std::runtime_error {"Invalid number of mesh dimensions."}; - } -} -bool StructuredMesh::intersects_1d(Position& r0, Position r1, int* ijk) const -{ - // Copy coordinates of starting point - double x0 = r0.x; - double y0 = r0.y; - double z0 = r0.z; - - // Copy coordinates of ending point - double x1 = r1.x; - double y1 = r1.y; - double z1 = r1.z; - - // Copy coordinates of mesh lower_left and upper_right - double xm0 = lower_left_[0]; - double xm1 = upper_right_[0]; - - double min_dist = INFTY; - - // Check if line intersects left surface -- calculate the intersection point - // (y,z) - if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) { - double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0); - double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0); - if (check_intersection_point(xm0, x0, yi, yi, zi, zi, r0, min_dist)) { - ijk[0] = 1; - } - } - - // Check if line intersects right surface -- calculate the intersection point - // (y,z) - if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) { - double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0); - double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0); - if (check_intersection_point(xm1, x0, yi, yi, zi, zi, r0, min_dist)) { - ijk[0] = shape_[0]; - } - } - - return min_dist < INFTY; -} - -bool StructuredMesh::intersects_2d(Position& r0, Position r1, int* ijk) const -{ - // Copy coordinates of starting point - double x0 = r0.x; - double y0 = r0.y; - double z0 = r0.z; - - // Copy coordinates of ending point - double x1 = r1.x; - double y1 = r1.y; - double z1 = r1.z; - - // Copy coordinates of mesh lower_left - double xm0 = lower_left_[0]; - double ym0 = lower_left_[1]; - - // Copy coordinates of mesh upper_right - double xm1 = upper_right_[0]; - double ym1 = upper_right_[1]; - - double min_dist = INFTY; - - // Check if line intersects left surface -- calculate the intersection point - // (y,z) - if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) { - double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0); - double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0); - if (yi >= ym0 && yi < ym1) { - if (check_intersection_point(xm0, x0, yi, y0, zi, zi, r0, min_dist)) { - ijk[0] = 1; - ijk[1] = get_index_in_direction(yi, 1); - } - } - } - - // Check if line intersects back surface -- calculate the intersection point - // (x,z) - if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) { - double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0); - double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0); - if (xi >= xm0 && xi < xm1) { - if (check_intersection_point(xi, x0, ym0, y0, zi, zi, r0, min_dist)) { - ijk[0] = get_index_in_direction(xi, 0); - ijk[1] = 1; - } - } - } - - // Check if line intersects right surface -- calculate the intersection point - // (y,z) - if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) { - double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0); - double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0); - if (yi >= ym0 && yi < ym1) { - if (check_intersection_point(xm1, x0, yi, y0, zi, zi, r0, min_dist)) { - ijk[0] = shape_[0]; - ijk[1] = get_index_in_direction(yi, 1); - } - } - } - - // Check if line intersects front surface -- calculate the intersection point - // (x,z) - if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) { - double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0); - double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0); - if (xi >= xm0 && xi < xm1) { - if (check_intersection_point(xi, x0, ym1, y0, zi, zi, r0, min_dist)) { - ijk[0] = get_index_in_direction(xi, 0); - ijk[1] = shape_[1]; - } - } - } - - return min_dist < INFTY; -} - -bool StructuredMesh::intersects_3d(Position& r0, Position r1, int* ijk) const -{ - // Copy coordinates of starting point - double x0 = r0.x; - double y0 = r0.y; - double z0 = r0.z; - - // Copy coordinates of ending point - double x1 = r1.x; - double y1 = r1.y; - double z1 = r1.z; - - // Copy coordinates of mesh lower_left - double xm0 = lower_left_[0]; - double ym0 = lower_left_[1]; - double zm0 = lower_left_[2]; - - // Copy coordinates of mesh upper_right - double xm1 = upper_right_[0]; - double ym1 = upper_right_[1]; - double zm1 = upper_right_[2]; - - double min_dist = INFTY; - - // Check if line intersects left surface -- calculate the intersection point - // (y,z) - if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) { - double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0); - double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0); - if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) { - if (check_intersection_point(xm0, x0, yi, y0, zi, z0, r0, min_dist)) { - ijk[0] = 1; - ijk[1] = get_index_in_direction(yi, 1); - ijk[2] = get_index_in_direction(zi, 2); - } - } - } - - // Check if line intersects back surface -- calculate the intersection point - // (x,z) - if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) { - double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0); - double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0); - if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) { - if (check_intersection_point(xi, x0, ym0, y0, zi, z0, r0, min_dist)) { - ijk[0] = get_index_in_direction(xi, 0); - ijk[1] = 1; - ijk[2] = get_index_in_direction(zi, 2); - } - } - } - - // Check if line intersects bottom surface -- calculate the intersection - // point (x,y) - if ((z0 < zm0 && z1 > zm0) || (z0 > zm0 && z1 < zm0)) { - double xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0); - double yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0); - if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) { - if (check_intersection_point(xi, x0, yi, y0, zm0, z0, r0, min_dist)) { - ijk[0] = get_index_in_direction(xi, 0); - ijk[1] = get_index_in_direction(yi, 1); - ijk[2] = 1; - } - } - } - - // Check if line intersects right surface -- calculate the intersection point - // (y,z) - if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) { - double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0); - double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0); - if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) { - if (check_intersection_point(xm1, x0, yi, y0, zi, z0, r0, min_dist)) { - ijk[0] = shape_[0]; - ijk[1] = get_index_in_direction(yi, 1); - ijk[2] = get_index_in_direction(zi, 2); - } - } - } - - // Check if line intersects front surface -- calculate the intersection point - // (x,z) - if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) { - double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0); - double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0); - if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) { - if (check_intersection_point(xi, x0, ym1, y0, zi, z0, r0, min_dist)) { - ijk[0] = get_index_in_direction(xi, 0); - ijk[1] = shape_[1]; - ijk[2] = get_index_in_direction(zi, 2); - } - } - } - - // Check if line intersects top surface -- calculate the intersection point - // (x,y) - if ((z0 < zm1 && z1 > zm1) || (z0 > zm1 && z1 < zm1)) { - double xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0); - double yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0); - if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) { - if (check_intersection_point(xi, x0, yi, y0, zm1, z0, r0, min_dist)) { - ijk[0] = get_index_in_direction(xi, 0); - ijk[1] = get_index_in_direction(yi, 1); - ijk[2] = shape_[2]; - } - } - } - - return min_dist < INFTY; -} - -void StructuredMesh::bins_crossed(Position r0, Position r1, const Direction& u, - vector& bins, vector& lengths) const -{ - // ======================================================================== - // Determine where the track intersects the mesh and if it intersects at all. + // TODO: when c++-17 is available, use "if constexpr ()" to compile-time enable/disable tally calls + // for now, T template type needs to provide both surface and track methods, which might be empty. + // modern optimizing compilers will (hopefully) eleminate the complete code (including calculation of parameters) + // but for the future: be explicit // Compute the length of the entire track. double total_distance = (r1 - r0).norm(); if (total_distance == 0.0) return; - // While determining if this track intersects the mesh, offset the starting - // and ending coords by a bit. This avoid finite-precision errors that can - // occur when the mesh surfaces coincide with lattice or geometric surfaces. - Position last_r = r0 + TINY_BIT * u; - Position r = r1 - TINY_BIT * u; + const int n = n_dimension_; - // Determine the mesh indices for the starting and ending coords. Here, we - // use arrays for ijk0 and ijk1 instead of vector because we obtain a - // small performance improvement by forcing this data to live on the stack, - // rather than on the heap. We know the maximum length is 3, and by - // ensuring that all loops are only indexed up to n_dimension, we will not - // access any non-initialized values. The same concept is used throughout. - int n = n_dimension_; - int ijk0[3], ijk1[3]; - bool start_in_mesh; - get_indices(last_r, ijk0, &start_in_mesh); - bool end_in_mesh; - get_indices(r, ijk1, &end_in_mesh); + bool in_mesh; - // Reset coordinates and check for a mesh intersection if necessary. - if (start_in_mesh) { - // The initial coords lie in the mesh, use those coords for tallying. - last_r = r0; - } else { - // The initial coords do not lie in the mesh. Check to see if the particle - // eventually intersects the mesh and compute the relevant coords and - // indices. - if (!intersects(last_r, r, ijk0)) - return; - } - r = r1; + // Calculate index of current cell. Offset the position a tiny bit in direction of flight + MeshIndex ijk = get_indices(r0+TINY_BIT*u, in_mesh); + - // The TINY_BIT offsets above mean that the preceding logic cannot always find - // the correct ijk0 and ijk1 indices. For tracks shorter than 2*TINY_BIT, just - // assume the track lies in only one mesh bin. These tracks are very short so - // any error caused by this assumption will be small. It is important that - // ijk0 values are used rather than ijk1 because the previous logic guarantees - // ijk0 is a valid mesh bin. + // if track is very short, assume that it is completely inside one cell. + // Only the current cell will score and no surfaces if (total_distance < 2 * TINY_BIT) { - for (int i = 0; i < n; ++i) - ijk1[i] = ijk0[i]; + if (in_mesh) { + tally.track(ijk, 1.0); + } + return; } - // ======================================================================== - // Find which mesh cells are traversed and the length of each traversal. + // Calculate initial distances to next surfaces in all three dimensions + std::array distances; + for (int k = 0; k < n; ++k) { + distances[k] = distance_to_grid_boundary(ijk, k, r0, u, 0.0); + } + // Position is r = r0 + u * l, start at r0 + double l { 0.0 }; + + // Loop until r = r1 is eventually reached while (true) { - if (std::equal(ijk0, ijk0 + n, ijk1)) { - // The track ends in this cell. Use the particle end location rather - // than the mesh surface and stop iterating. - double distance = (last_r - r).norm(); - bins.push_back(get_bin_from_indices(ijk0)); - lengths.push_back(distance / total_distance); - break; - } - // The track exits this cell. Determine the distance to each mesh surface. - double d[3]; - for (int k = 0; k < n; ++k) { - if (std::fabs(u[k]) < FP_PRECISION) { - d[k] = INFTY; - } else if (u[k] > 0) { - double xyz_cross = positive_grid_boundary(ijk0, k); - d[k] = (xyz_cross - last_r[k]) / u[k]; - } else { - double xyz_cross = negative_grid_boundary(ijk0, k); - d[k] = (xyz_cross - last_r[k]) / u[k]; + if (in_mesh) { + + // find surface with minimal distance to current position + const auto k = std::min_element(distances.begin(), distances.end()) - distances.begin(); + + // Tally track length delta since last step + tally.track(ijk, (std::min(distances[k].distance, total_distance) - l) / total_distance); + + // update position and leave, if we have reached end position + l = distances[k].distance; + if (l >= total_distance) + return; + + // If we have not reached r1, we have hit a surface. Tally outward current + tally.surface(ijk, k, distances[k].maxSurface, 0); + + // Update cell and calculate distance to next surface in k-direction. + // The two other directions are still valid! + ijk[k] = distances[k].nextIndex; + distances[k] = distance_to_grid_boundary(ijk, k, r0, u, l); + + // Check if we have left the interiour of the mesh + in_mesh = ((ijk[k]>=1) and (ijk[k]<=shape_[k])); + + // If we are still inside the mesh, tally inward current for the next cell + if (in_mesh) tally.surface(ijk, k, !distances[k].maxSurface, 1); + + } 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 { 0 }; + for (int k = 0; k < n; ++k) { + if ((ijk[k] < 1 || ijk[k] > shape_[k]) and (distances[k].distance > l)) { + l = distances[k].distance; + k_max = k; + } } - } - // Pick the closest mesh surface and append this traversal to the output. - auto j = std::min_element(d, d + n) - d; - double distance = d[j]; - bins.push_back(get_bin_from_indices(ijk0)); - lengths.push_back(distance / total_distance); + // If r1 is not inside the mesh, exit here + if (l >= total_distance) + return; - // Translate to the oncoming mesh surface. - last_r += distance * u; - - // Increment the indices into the next mesh cell. - if (u[j] > 0.0) { - ++ijk0[j]; - } else { - --ijk0[j]; - } - - // If the next indices are invalid, then the track has left the mesh and - // we are done. - bool in_mesh = true; - for (int i = 0; i < n; ++i) { - if (ijk0[i] < 1 || ijk0[i] > shape_[i]) { - in_mesh = false; - break; + // Calculate the new cell index and update all distances to next surfaces. + ijk = get_indices(r0+(l+TINY_BIT)*u, in_mesh); + for (int k = 0; k < n; ++k) { + distances[k] = distance_to_grid_boundary(ijk, k, r0, u, l); } + + // If inside the mesh, Tally inward current + if (in_mesh) tally.surface(ijk, k_max, !distances[k_max].maxSurface, 1); + } - if (!in_mesh) - break; } } +void StructuredMesh::bins_crossed(Position r0, Position r1, const Direction& u, + vector& bins, vector& lengths) const { + + // Helper tally class. + // stores a pointer to the mesh class and references to bins and lengths parameters. + // Performs the actual tally through the track method. + struct TallyBins { + TallyBins(const StructuredMesh* _mesh, vector& _bins, vector& _lengths): + mesh(_mesh), bins(_bins), lengths(_lengths) {} + void surface(const MeshIndex& ijk, int k, bool max, int inwward) const {} + void track(const MeshIndex& ijk, double l) const { + bins.push_back(mesh->get_bin_from_indices(ijk)); + lengths.push_back(l); + } + + const StructuredMesh* mesh; + vector& bins; + vector& lengths; + }; + + // Perform the mesh raytrace with the helper class. + raytrace_mesh(r0, r1, u, TallyBins(this, bins, lengths)); +} + +void StructuredMesh::surface_bins_crossed( + Position r0, Position r1, const Direction& u, vector& bins) const +{ + + // Helper tally class. + // stores a pointer to the mesh class and a reference to the bins parameter. + // Performs the actual tally through the surface method. + struct TallyBins { + TallyBins(const StructuredMesh* _mesh, vector& _bins): + mesh(_mesh), bins(_bins) {} + void surface(const MeshIndex& ijk, int k, bool max, int inward) const { + int i_bin = 4 * mesh->n_dimension_ * mesh->get_bin_from_indices(ijk) + 4 * k + inward; + if (max) i_bin +=2; + bins.push_back(i_bin); + } + void track(const MeshIndex& idx, double l) const {} + + const StructuredMesh* mesh; + vector& bins; + }; + + // Perform the mesh raytrace with the helper class. + raytrace_mesh(r0, r1, u, TallyBins(this, bins)); +} + + //============================================================================== // RegularMesh implementation //============================================================================== @@ -802,160 +604,35 @@ int RegularMesh::get_index_in_direction(double r, int i) const return std::ceil((r - lower_left_[i]) / width_[i]); } -double RegularMesh::positive_grid_boundary(int* ijk, int i) const +double RegularMesh::positive_grid_boundary(const MeshIndex& ijk, int i) const { return lower_left_[i] + ijk[i] * width_[i]; } -double RegularMesh::negative_grid_boundary(int* ijk, int i) const +double RegularMesh::negative_grid_boundary(const MeshIndex& ijk, int i) const { return lower_left_[i] + (ijk[i] - 1) * width_[i]; } -void RegularMesh::surface_bins_crossed( - Position r0, Position r1, const Direction& u, vector& bins) const +StructuredMesh::MeshDistance RegularMesh::distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const { - // Determine indices for starting and ending location. - int n = n_dimension_; - vector ijk0(n), ijk1(n); - bool start_in_mesh; - get_indices(r0, ijk0.data(), &start_in_mesh); - bool end_in_mesh; - get_indices(r1, ijk1.data(), &end_in_mesh); + MeshDistance d; + d.nextIndex = ijk[i]; + if (std::fabs(u[i]) < FP_PRECISION) + return d; - // Check if the track intersects any part of the mesh. - if (!start_in_mesh) { - Position r0_copy = r0; - vector ijk0_copy(ijk0); - if (!intersects(r0_copy, r1, ijk0_copy.data())) - return; - } - - // ======================================================================== - // Find which mesh surfaces are crossed. - - // Calculate number of surface crossings - int n_cross = 0; - for (int i = 0; i < n; ++i) - n_cross += std::abs(ijk1[i] - ijk0[i]); - if (n_cross == 0) - return; - - // Bounding coordinates - Position xyz_cross; - for (int i = 0; i < n; ++i) { - if (u[i] > 0.0) { - xyz_cross[i] = positive_grid_boundary(ijk0.data(), i); - } else { - xyz_cross[i] = negative_grid_boundary(ijk0.data(), i); - } - } - - for (int j = 0; j < n_cross; ++j) { - // Set the distances to infinity - Position d {INFTY, INFTY, INFTY}; - - // Determine closest bounding surface. We need to treat - // special case where the cosine of the angle is zero since this would - // result in a divide-by-zero. - double distance = INFTY; - for (int i = 0; i < n; ++i) { - if (u[i] == 0) { - d[i] = INFTY; - } else { - d[i] = (xyz_cross[i] - r0[i]) / u[i]; - } - distance = std::min(distance, d[i]); - } - - // Loop over the dimensions - for (int i = 0; i < n; ++i) { - // Check whether distance is the shortest distance - if (distance == d[i]) { - - // Check whether the current indices are within the mesh bounds - bool in_mesh = true; - for (int j = 0; j < n; ++j) { - if (ijk0[j] < 1 || ijk0[j] > shape_[j]) { - in_mesh = false; - break; - } - } - - // Check whether particle is moving in positive i direction - if (u[i] > 0) { - - // Outward current on i max surface - if (in_mesh) { - int i_surf = 4 * i + 3; - int i_mesh = get_bin_from_indices(ijk0.data()); - int i_bin = 4 * n * i_mesh + i_surf - 1; - - bins.push_back(i_bin); - } - - // Advance position - ++ijk0[i]; - xyz_cross[i] += width_[i]; - in_mesh = true; - for (int j = 0; j < n; ++j) { - if (ijk0[j] < 1 || ijk0[j] > shape_[j]) { - in_mesh = false; - break; - } - } - - // If the particle crossed the surface, tally the inward current on - // i min surface - if (in_mesh) { - int i_surf = 4 * i + 2; - int i_mesh = get_bin_from_indices(ijk0.data()); - int i_bin = 4 * n * i_mesh + i_surf - 1; - - bins.push_back(i_bin); - } - - } else { - // The particle is moving in the negative i direction - - // Outward current on i min surface - if (in_mesh) { - int i_surf = 4 * i + 1; - int i_mesh = get_bin_from_indices(ijk0.data()); - int i_bin = 4 * n * i_mesh + i_surf - 1; - - bins.push_back(i_bin); - } - - // Advance position - --ijk0[i]; - xyz_cross[i] -= width_[i]; - in_mesh = true; - for (int j = 0; j < n; ++j) { - if (ijk0[j] < 1 || ijk0[j] > shape_[j]) { - in_mesh = false; - break; - } - } - - // If the particle crossed the surface, tally the inward current on - // i max surface - if (in_mesh) { - int i_surf = 4 * i + 4; - int i_mesh = get_bin_from_indices(ijk0.data()); - int i_bin = 4 * n * i_mesh + i_surf - 1; - - bins.push_back(i_bin); - } - } - } - } - - // Calculate new coordinates - r0 += distance * u; + d.maxSurface = (u[i] > 0); + if (d.maxSurface and (ijk[i] <= shape_[i])) { + d.nextIndex++; + d.distance = (positive_grid_boundary(ijk, i) - r0[i]) / u[i]; + } else if (not d.maxSurface and (ijk[i] >= 1)) { + d.nextIndex--; + d.distance = (negative_grid_boundary(ijk, i) - r0[i]) / u[i]; } + return d; } + std::pair, vector> RegularMesh::plot( Position plot_ll, Position plot_ur) const { @@ -1065,6 +742,8 @@ xt::xtensor RegularMesh::count_sites( return counts; } + + //============================================================================== // RectilinearMesh implementation //============================================================================== @@ -1083,16 +762,34 @@ RectilinearMesh::RectilinearMesh(pugi::xml_node node) : StructuredMesh {node} } } -double RectilinearMesh::positive_grid_boundary(int* ijk, int i) const +double RectilinearMesh::positive_grid_boundary(const MeshIndex& ijk, int i) const { return grid_[i][ijk[i]]; } -double RectilinearMesh::negative_grid_boundary(int* ijk, int i) const +double RectilinearMesh::negative_grid_boundary(const MeshIndex& ijk, int i) const { return grid_[i][ijk[i] - 1]; } +StructuredMesh::MeshDistance RectilinearMesh::distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const +{ + MeshDistance d; + d.nextIndex = ijk[i]; + if (std::fabs(u[i]) < FP_PRECISION) + return d; + + d.maxSurface = (u[i] > 0); + if (d.maxSurface and (ijk[i] <= shape_[i])) { + d.nextIndex++; + d.distance = (positive_grid_boundary(ijk, i) - r0[i]) / u[i]; + } else if (not d.maxSurface and (ijk[i] > 0)) { + d.nextIndex--; + d.distance = (negative_grid_boundary(ijk, i) - r0[i]) / u[i]; + } + return d; +} + int RectilinearMesh::set_grid() { shape_ = {static_cast(grid_[0].size()) - 1, @@ -1120,168 +817,7 @@ int RectilinearMesh::set_grid() return 0; } -void RectilinearMesh::surface_bins_crossed( - Position r0, Position r1, const Direction& u, vector& bins) const -{ - // Determine indices for starting and ending location. - int ijk0[3], ijk1[3]; - bool start_in_mesh; - get_indices(r0, ijk0, &start_in_mesh); - bool end_in_mesh; - get_indices(r1, ijk1, &end_in_mesh); - // If the starting coordinates do not lie in the mesh, compute the coords and - // mesh indices of the first intersection, and add the bin for this first - // intersection. Return if the particle does not intersect the mesh at all. - if (!start_in_mesh) { - // Compute the incoming intersection coordinates and indices. - if (!intersects(r0, r1, ijk0)) - return; - - // Determine which surface the particle entered. - double min_dist = INFTY; - int i_surf; - for (int i = 0; i < 3; ++i) { - if (u[i] > 0.0 && ijk0[i] == 1) { - double d = std::abs(r0[i] - grid_[i][0]); - if (d < min_dist) { - min_dist = d; - i_surf = 4 * i + 2; - } - } else if (u[i] < 0.0 && ijk0[i] == shape_[i]) { - double d = std::abs(r0[i] - grid_[i][shape_[i]]); - if (d < min_dist) { - min_dist = d; - i_surf = 4 * i + 4; - } - } // u[i] == 0 intentionally skipped - } - - // Add the incoming current bin. - int i_mesh = get_bin_from_indices(ijk0); - int i_bin = 4 * 3 * i_mesh + i_surf - 1; - bins.push_back(i_bin); - } - - // If the ending coordinates do not lie in the mesh, compute the coords and - // mesh indices of the last intersection, and add the bin for this last - // intersection. - if (!end_in_mesh) { - // Compute the outgoing intersection coordinates and indices. - intersects(r1, r0, ijk1); - - // Determine which surface the particle exited. - double min_dist = INFTY; - int i_surf; - for (int i = 0; i < 3; ++i) { - if (u[i] > 0.0 && ijk1[i] == shape_[i]) { - double d = std::abs(r1[i] - grid_[i][shape_[i]]); - if (d < min_dist) { - min_dist = d; - i_surf = 4 * i + 3; - } - } else if (u[i] < 0.0 && ijk1[i] == 1) { - double d = std::abs(r1[i] - grid_[i][0]); - if (d < min_dist) { - min_dist = d; - i_surf = 4 * i + 1; - } - } // u[i] == 0 intentionally skipped - } - - // Add the outgoing current bin. - int i_mesh = get_bin_from_indices(ijk1); - int i_bin = 4 * 3 * i_mesh + i_surf - 1; - bins.push_back(i_bin); - } - - // ======================================================================== - // Find which mesh surfaces are crossed. - - // Calculate number of surface crossings - int n_cross = 0; - for (int i = 0; i < 3; ++i) - n_cross += std::abs(ijk1[i] - ijk0[i]); - if (n_cross == 0) - return; - - // Bounding coordinates - Position xyz_cross; - for (int i = 0; i < 3; ++i) { - if (u[i] > 0.0) { - xyz_cross[i] = positive_grid_boundary(ijk0, i); - } else { - xyz_cross[i] = negative_grid_boundary(ijk0, i); - } - } - - for (int j = 0; j < n_cross; ++j) { - // Set the distances to infinity - Position d {INFTY, INFTY, INFTY}; - - // Determine closest bounding surface. We need to treat - // special case where the cosine of the angle is zero since this would - // result in a divide-by-zero. - double distance = INFTY; - for (int i = 0; i < 3; ++i) { - if (u[i] == 0) { - d[i] = INFTY; - } else { - d[i] = (xyz_cross[i] - r0[i]) / u[i]; - } - distance = std::min(distance, d[i]); - } - - // Loop over the dimensions - for (int i = 0; i < 3; ++i) { - // Check whether distance is the shortest distance - if (distance == d[i]) { - - // Check whether particle is moving in positive i direction - if (u[i] > 0) { - - // Outward current on i max surface - int i_surf = 4 * i + 3; - int i_mesh = get_bin_from_indices(ijk0); - int i_bin = 4 * 3 * i_mesh + i_surf - 1; - bins.push_back(i_bin); - - // Advance position - ++ijk0[i]; - xyz_cross[i] = grid_[i][ijk0[i]]; - - // Inward current on i min surface - i_surf = 4 * i + 2; - i_mesh = get_bin_from_indices(ijk0); - i_bin = 4 * 3 * i_mesh + i_surf - 1; - bins.push_back(i_bin); - - } else { - // The particle is moving in the negative i direction - - // Outward current on i min surface - int i_surf = 4 * i + 1; - int i_mesh = get_bin_from_indices(ijk0); - int i_bin = 4 * 3 * i_mesh + i_surf - 1; - bins.push_back(i_bin); - - // Advance position - --ijk0[i]; - xyz_cross[i] = grid_[i][ijk0[i] - 1]; - - // Inward current on i min surface - i_surf = 4 * i + 4; - i_mesh = get_bin_from_indices(ijk0); - i_bin = 4 * 3 * i_mesh + i_surf - 1; - bins.push_back(i_bin); - } - } - } - - // Calculate new coordinates - r0 += distance * u; - } -} int RectilinearMesh::get_index_in_direction(double r, int i) const { @@ -1330,6 +866,484 @@ void RectilinearMesh::to_hdf5(hid_t group) const close_group(mesh_group); } +//============================================================================== +// CylindricalMesh implementation +//============================================================================== + +CylindricalMesh::CylindricalMesh(pugi::xml_node node) : StructuredMesh {node} +{ + n_dimension_ = 3; + + grid_.resize(n_dimension_); + grid_[0] = get_node_array(node, "r_grid"); + grid_[1] = get_node_array(node, "p_grid"); + grid_[2] = get_node_array(node, "z_grid"); + + if (int err = set_grid()) { + fatal_error(openmc_err_msg); + } +} + +StructuredMesh::MeshIndex CylindricalMesh::get_indices(Position r, bool& in_mesh) const +{ + Position mapped_r; + + mapped_r[0] = std::hypot(r.x, r.y); + mapped_r[2] = r[2]; + if (mapped_r[0] < FP_PRECISION) { + mapped_r[1] = 0.0; + } else { + mapped_r[1] = std::atan2(r.y, r.x); + if (mapped_r[1] < 0) mapped_r[1] += 2*M_PI; + } + + MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh); + + mapped_r[1] = sanitize_phi(mapped_r[1]); + + return idx; +} + +double CylindricalMesh::find_r_crossing(const Position& r, const Direction& u, double l, int shell) const +{ + + if ((shell < 0) || (shell >= shape_[0])) return INFTY; + + //solve r.x^2 + r.y^2 == r0^2 + // x^2 + 2*s*u*x + s^2*u^2 + s^2*v^2+2*s*v*y + y^2 -r0^2 = 0 + // s^2 * (u^2 + v^2) + 2*s*(u*x+v*y) + x^2+y^2-r0^2 = 0 + + const double r0 = grid_[0][shell]; + const double denominator = u.x*u.x + u.y*u.y; + + // Direction of flight is in z-direction. Will never intersect r. + if (fabs(denominator) < FP_PRECISION) return INFTY; + + const double p = (u.x*r.x + u.y*r.y) / denominator; + double D = p*p + (r0*r0 - r.x*r.x - r.y*r.y) / denominator; + + if (D >= 0.0) { + D = std::sqrt(D); + + // the solution -p - D is always smaller as -p + D : Check this one first + if (-p - D > l) + return -p - D; + if ( -p + D > l) + return -p + D; + } + + return INFTY; +} + +double CylindricalMesh::find_phi_crossing(const Position& r, const Direction& u, double l, int shell) const +{ + // Phi grid is [0, 360], thus there is no real surface to cross + if (full_phi and (shape_[1]==1)) return INFTY; + + shell = sanitize_phi(shell); + + const double p0 = grid_[1][shell]; + + // solve y(s)/x(s) = tan(p0) = sin(p0)/cos(p0) + // => x(s) * cos(p0) = y(s) * sin(p0) + // => (y + s * v) * cos(p0) = (x + s * u) * sin(p0) + // = s * (v * cos(p0) - u * sin(p0)) = - (y * cos(p0) - x * sin(p0)) + + const double c0 = std::cos(p0); + const double s0 = std::sin(p0); + + const double denominator = (u.x * s0 - u.y * c0); + + // Check if direction of flight is not parallel to phi surface + if (std::fabs(denominator) > FP_PRECISION) { + const double s = - (r.x * s0 - r.y * c0) / denominator; + // Check if solution is in positive direction of flight and crosses the correct phi surface (not -phi) + if ((s>l) and ((c0*(r.x+s*u.x) + s0*(r.y+s*u.y))>0.0)) + return s; + } + + return INFTY; +} + +StructuredMesh::MeshDistance CylindricalMesh::find_z_crossing(const Position& r, const Direction& u, double l, int shell) const +{ + MeshDistance d; + d.nextIndex = shell; + + // Direction of flight is within xy-plane. Will never intersect z. + if (std::fabs(u.z) < FP_PRECISION) + return d; + + d.maxSurface = (u.z > 0.0); + if (d.maxSurface and (shell <= shape_[2])) { + d.nextIndex += 1; + d.distance = (grid_[2][shell] - r.z) / u.z; + } else if (not d.maxSurface and (shell > 0)) { + d.nextIndex -= 1; + d.distance = (grid_[2][shell-1] - r.z) / u.z; + } + return d; +} + +StructuredMesh::MeshDistance CylindricalMesh::distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, 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)) + ); + + } 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)) + ); + + } else { + return find_z_crossing(r0, u, l, ijk[i] ); + } + +} + +int CylindricalMesh::set_grid() +{ + shape_ = {static_cast(grid_[0].size()) - 1, + static_cast(grid_[1].size()) - 1, + static_cast(grid_[2].size()) - 1}; + + for (const auto& g : grid_) { + if (g.size() < 2) { + set_errmsg("r-, phi-, and z- grids for cylindrical meshes " + "must each have at least 2 points"); + return OPENMC_E_INVALID_ARGUMENT; + } + for (int i = 1; i < g.size(); ++i) { + if (g[i] <= g[i - 1]) { + set_errmsg("Values in for r-, phi-, and z- grids for " + "cylindrical meshes must be sorted and unique."); + return OPENMC_E_INVALID_ARGUMENT; + } + } + } + if (grid_[0].front() < 0.0) { + set_errmsg("r-grid for " + "cylindrical meshes must start at r >= 0."); + return OPENMC_E_INVALID_ARGUMENT; + } + if (grid_[1].front() < 0.0) { + set_errmsg("phi-grid for " + "cylindrical meshes must start at phi >= 0."); + return OPENMC_E_INVALID_ARGUMENT; + } + if (grid_[1].back() > 360) { + set_errmsg("phi-grids for " + "cylindrical meshes must end with theta <= 360 degree."); + + return OPENMC_E_INVALID_ARGUMENT; + } + + full_phi = (grid_[1].front() == 0.0) and (grid_[1].back()==360.0); + + // Transform phi-grid from degrees to radians + + std::transform(grid_[1].begin(), grid_[1].end(), grid_[1].begin(), [](double d){ return M_PI*d/180.0; }); + + lower_left_ = {grid_[0].front(), grid_[1].front(), grid_[2].front()}; + upper_right_ = {grid_[0].back(), grid_[1].back(), grid_[2].back()}; + + return 0; +} + + +int CylindricalMesh::get_index_in_direction(double r, int i) const +{ + return lower_bound_index(grid_[i].begin(), grid_[i].end(), r) + 1; +} + +std::pair, vector> CylindricalMesh::plot( + Position plot_ll, Position plot_ur) const +{ + fatal_error("Plot of cylindrical 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 CylindricalMesh::to_hdf5(hid_t group) const +{ + hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_)); + + write_dataset(mesh_group, "type", "cylindrical"); + write_dataset(mesh_group, "r_grid", grid_[0]); + write_dataset(mesh_group, "p_grid", grid_[1]); + write_dataset(mesh_group, "z_grid", grid_[2]); + + close_group(mesh_group); +} + +//============================================================================== +// SphericalMesh implementation +//============================================================================== + +SphericalMesh::SphericalMesh(pugi::xml_node node) : StructuredMesh {node} +{ + n_dimension_ = 3; + + grid_.resize(n_dimension_); + grid_[0] = get_node_array(node, "r_grid"); + grid_[1] = get_node_array(node, "t_grid"); + grid_[2] = get_node_array(node, "p_grid"); + + if (int err = set_grid()) { + fatal_error(openmc_err_msg); + } +} + +StructuredMesh::MeshIndex SphericalMesh::get_indices(Position r, bool& in_mesh) const +{ + Position mapped_r; + + mapped_r[0] = r.norm(); + if (mapped_r[0] < FP_PRECISION) { + mapped_r[1] = 0.0; + mapped_r[2] = 0.0; + } else { + mapped_r[1] = std::acos(r.z / mapped_r.x); + mapped_r[2] = std::atan2(r.y, r.x); + if (mapped_r[2] < 0) mapped_r[2] += 2*M_PI; + } + + MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh); + + mapped_r[1] = sanitize_theta(mapped_r[1]); + mapped_r[2] = sanitize_phi(mapped_r[2]); + + return idx; +} + +double SphericalMesh::find_r_crossing(const Position& r, const Direction& u, double l, int shell) const +{ + + ////std::cout << " r: " << shell << " " << grid_[0].size() << " " << shape_[0] << "\n"; + + if ((shell < 0) || (shell >= shape_[0])) return INFTY; + + // solve |r+s*u| = r0 + // |r+s*u| = |r| + 2*s*r*u + s^2 (|u|==1 !) + + const double r0 = grid_[0][shell]; + const double p = r.dot(u); + double D = p * p - r.dot(r) + r0 * r0; + + if (D >= 0.0) { + D = std::sqrt(D); + // the solution -p - D is always smaller as -p + D : Check this one first + if (-p - D > l) + return -p - D; + if ( -p + D > l) + return -p + D; + } + + return INFTY; +} + +double SphericalMesh::find_theta_crossing(const Position& r, const Direction& u, double l, int shell) const +{ + // Theta grid is [0, 180], thus there is no real surface to cross + if (full_theta and (shape_[1]==1)) return INFTY; + + shell = sanitize_theta(shell); + + + // solving z(s) = cos/theta) * r(s) with r(s) = r+s*u + // yields + // a*s^2 + 2*b*s + c == 0 with + // a = cos(theta)^2 - u.z * u.z + // b = r*u * cos(theta)^2 - u.z * r.z + // c = r*r * cos(theta)^2 - r.z^2 + + const double cos_t = std::cos(grid_[1][shell]); + const bool sgn = std::signbit(cos_t); + const double cos_t_2 = cos_t * cos_t; + + const double a = cos_t_2 - u.z * u.z; + const double b = (r.dot(u) * cos_t_2 - r.z * u.z); + const double c = (r.dot(r) * cos_t_2 - r.z * r.z); + + // if factor of s^2 is zero, direction of flight is parallel to theta surface + if (fabs(a) < FP_PRECISION) { + // if b vanishes, direction of flight is within theta surface and crossing is not possible + if (fabs(b) > FP_PRECISION) { + const double s = - 0.5 * c / b; + // Check if solution is in positive direction of flight and has correct sign + if ((s > l) and (std::signbit(r.z+s*u.z) == sgn)) return s; + } + return INFTY; + } + + const double p = b / a; + double D = p*p - c / a; + + if (D < 0.0) + return INFTY; + + D = std::sqrt(D); + + // the solution -p-D is always smaller as -p+D : Check this one first + double s = - p - D; + // Check if solution is in positive direction of flight and has correct sign + if ((s > l) and (std::signbit(r.z+s*u.z) == sgn)) return s; + + s = - p + D; + // Check if solution is in positive direction of flight and has correct sign + if ((s > l) and (std::signbit(r.z+s*u.z) == sgn)) return s; + + return INFTY; +} + +double SphericalMesh::find_phi_crossing(const Position& r, const Direction& u, double l, int shell) const +{ + // Phi grid is [0, 360], thus there is no real surface to cross + if (full_phi and (shape_[2]==1)) return INFTY; + + shell = sanitize_phi(shell); + + const double p0 = grid_[2][shell]; + + // solve y(s)/x(s) = tan(p0) = sin(p0)/cos(p0) + // => x(s) * cos(p0) = y(s) * sin(p0) + // => (y + s * v) * cos(p0) = (x + s * u) * sin(p0) + // = s * (v * cos(p0) - u * sin(p0)) = - (y * cos(p0) - x * sin(p0)) + + const double c0 = std::cos(p0); + const double s0 = std::sin(p0); + + const double denominator = (u.x * s0 - u.y * c0); + + // Check if direction of flight is not parallel to phi surface + if (std::fabs(denominator) > FP_PRECISION) { + const double s = - (r.x * s0 - r.y * c0) / denominator; + // Check if solution is in positive direction of flight and crosses the correct phi surface (not -phi) + if ((s>l) and ((c0*(r.x+s*u.x) + s0*(r.y+s*u.y))>0.0)) + return s; + } + + return INFTY; +} + + +StructuredMesh::MeshDistance SphericalMesh::distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, 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)) + ); + + } 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)) + ); + + } 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)) + ); + } + +} + +int SphericalMesh::set_grid() +{ + shape_ = {static_cast(grid_[0].size()) - 1, + static_cast(grid_[1].size()) - 1, + static_cast(grid_[2].size()) - 1}; + + for (const auto& g : grid_) { + if (g.size() < 2) { + set_errmsg("x-, y-, and z- grids for spherical meshes " + "must each have at least 2 points"); + return OPENMC_E_INVALID_ARGUMENT; + } + for (int i = 1; i < g.size(); ++i) { + if (g[i] <= g[i - 1]) { + set_errmsg("Values in for r-, theta-, and phi- grids for " + "spherical meshes must be sorted and unique."); + return OPENMC_E_INVALID_ARGUMENT; + } + } + if (g.front() < 0.0) { + set_errmsg("r-, theta-, and phi- grids for " + "spherical meshes must start at v >= 0."); + return OPENMC_E_INVALID_ARGUMENT; + } + } + if (grid_[1].back() > 180) { + set_errmsg("theta-grids for " + "spherical meshes must end with theta <= 180 degree."); + + return OPENMC_E_INVALID_ARGUMENT; + } + if (grid_[2].back() > 360) { + set_errmsg("phi-grids for " + "spherical meshes must end with phi <= 180 degree."); + return OPENMC_E_INVALID_ARGUMENT; + } + + full_theta = (grid_[1].front() == 0.0) and (grid_[1].back()==180.0); + full_phi = (grid_[2].front() == 0.0) and (grid_[2].back()==360.0); + + // Transform theta- and phi-grid from degrees to radians + for (int i=1; i<3; i++) + std::transform(grid_[i].begin(), grid_[i].end(), grid_[i].begin(), [](double d){ return M_PI*d/180.0; }); + + lower_left_ = {grid_[0].front(), grid_[1].front(), grid_[2].front()}; + upper_right_ = {grid_[0].back(), grid_[1].back(), grid_[2].back()}; + + return 0; +} + + +int SphericalMesh::get_index_in_direction(double r, int i) const +{ + return lower_bound_index(grid_[i].begin(), grid_[i].end(), r) + 1; +} + +std::pair, vector> SphericalMesh::plot( + Position plot_ll, Position plot_ur) const +{ + fatal_error("Plot of spherical 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 SphericalMesh::to_hdf5(hid_t group) const +{ + hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_)); + + write_dataset(mesh_group, "type", "spherical"); + write_dataset(mesh_group, "r_grid", grid_[0]); + write_dataset(mesh_group, "t_grid", grid_[1]); + write_dataset(mesh_group, "p_grid", grid_[2]); + + close_group(mesh_group); +} + + + + //============================================================================== // Helper functions for the C API //============================================================================== @@ -1357,6 +1371,13 @@ int check_mesh_type(int32_t index) return 0; } +template +bool is_mesh_type(int32_t index) +{ + T* mesh = dynamic_cast(model::meshes[index].get()); + return mesh; +} + //============================================================================== // C API functions //============================================================================== @@ -1367,15 +1388,14 @@ extern "C" int openmc_mesh_get_type(int32_t index, char* type) if (int err = check_mesh(index)) return err; - RegularMesh* mesh = dynamic_cast(model::meshes[index].get()); - if (mesh) { + if (is_mesh_type(index)) { std::strcpy(type, "regular"); - } else { - RectilinearMesh* mesh = - dynamic_cast(model::meshes[index].get()); - if (mesh) { - std::strcpy(type, "rectilinear"); - } + } else if (is_mesh_type(index)) { + std::strcpy(type, "rectilinear"); + } else if (is_mesh_type(index)) { + std::strcpy(type, "cylindrical"); + } else if (is_mesh_type(index)) { + std::strcpy(type, "spherical"); } return 0; } @@ -1393,6 +1413,10 @@ extern "C" int openmc_extend_meshes( model::meshes.push_back(make_unique()); } else if (std::strcmp(type, "rectilinear") == 0) { model::meshes.push_back(make_unique()); + } else if (std::strcmp(type, "cylindrical") == 0) { + model::meshes.push_back(make_unique()); + } else if (std::strcmp(type, "spherical") == 0) { + model::meshes.push_back(make_unique()); } else { throw std::runtime_error {"Unknown mesh type: " + std::string(type)}; } @@ -1597,6 +1621,108 @@ extern "C" int openmc_rectilinear_mesh_set_grid(int32_t index, return err; } +//! Get the cylindrical mesh grid +extern "C" int openmc_cylindrical_mesh_get_grid(int32_t index, double** grid_x, + int* nx, double** grid_y, int* ny, double** grid_z, int* nz) +{ + if (int err = check_mesh_type(index)) + return err; + CylindricalMesh* m = + dynamic_cast(model::meshes[index].get()); + + if (m->lower_left_.dimension() == 0) { + set_errmsg("Mesh parameters have not been set."); + return OPENMC_E_ALLOCATE; + } + + *grid_x = m->grid_[0].data(); + *nx = m->grid_[0].size(); + *grid_y = m->grid_[1].data(); + *ny = m->grid_[1].size(); + *grid_z = m->grid_[2].data(); + *nz = m->grid_[2].size(); + + return 0; +} + +//! Set the cylindrical mesh parameters +extern "C" int openmc_cylindrical_mesh_set_grid(int32_t index, + const double* grid_x, const int nx, const double* grid_y, const int ny, + const double* grid_z, const int nz) +{ + if (int err = check_mesh_type(index)) + return err; + CylindricalMesh* m = + dynamic_cast(model::meshes[index].get()); + + m->n_dimension_ = 3; + m->grid_.resize(m->n_dimension_); + + for (int i = 0; i < nx; i++) { + m->grid_[0].push_back(grid_x[i]); + } + for (int i = 0; i < ny; i++) { + m->grid_[1].push_back(grid_y[i]); + } + for (int i = 0; i < nz; i++) { + m->grid_[2].push_back(grid_z[i]); + } + + int err = m->set_grid(); + return err; +} + +//! Get the spherical mesh grid +extern "C" int openmc_spherical_mesh_get_grid(int32_t index, double** grid_x, + int* nx, double** grid_y, int* ny, double** grid_z, int* nz) +{ + if (int err = check_mesh_type(index)) + return err; + SphericalMesh* m = + dynamic_cast(model::meshes[index].get()); + + if (m->lower_left_.dimension() == 0) { + set_errmsg("Mesh parameters have not been set."); + return OPENMC_E_ALLOCATE; + } + + *grid_x = m->grid_[0].data(); + *nx = m->grid_[0].size(); + *grid_y = m->grid_[1].data(); + *ny = m->grid_[1].size(); + *grid_z = m->grid_[2].data(); + *nz = m->grid_[2].size(); + + return 0; +} + +//! Set the spherical mesh parameters +extern "C" int openmc_spherical_mesh_set_grid(int32_t index, + const double* grid_x, const int nx, const double* grid_y, const int ny, + const double* grid_z, const int nz) +{ + if (int err = check_mesh_type(index)) + return err; + SphericalMesh* m = + dynamic_cast(model::meshes[index].get()); + + m->n_dimension_ = 3; + m->grid_.resize(m->n_dimension_); + + for (int i = 0; i < nx; i++) { + m->grid_[0].push_back(grid_x[i]); + } + for (int i = 0; i < ny; i++) { + m->grid_[1].push_back(grid_y[i]); + } + for (int i = 0; i < nz; i++) { + m->grid_[2].push_back(grid_z[i]); + } + + int err = m->set_grid(); + return err; +} + #ifdef DAGMC MOABMesh::MOABMesh(pugi::xml_node node) : UnstructuredMesh(node) @@ -2442,6 +2568,10 @@ void read_meshes(pugi::xml_node root) model::meshes.push_back(make_unique(node)); } else if (mesh_type == "rectilinear") { model::meshes.push_back(make_unique(node)); + } else if (mesh_type == "cylindrical") { + model::meshes.push_back(make_unique(node)); + } else if (mesh_type == "spherical") { + model::meshes.push_back(make_unique(node)); #ifdef DAGMC } else if (mesh_type == "unstructured" && mesh_lib == "moab") { model::meshes.push_back(make_unique(node));