From d723dd7f07337d7cca433951ceb1ca1a6a65a055 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Mon, 20 Dec 2021 19:41:52 +0100 Subject: [PATCH 01/17] Add Cylindrical & Spherical Mesh for Tally Filter New method for ray-tracing the mesh, used for both MeshFilter & MeshSurfaceFilter. For all Structured Meshes! meshes negative/positive_grid_boundary to non-virtual function for speed MeshIndex typedef now a std::array instead of int[3] -> Get rid of pointers (see e.g. core quidelines I.13) MeshDistance class to stroe all relevant information --- include/openmc/mesh.h | 188 +++++- openmc/lib/mesh.py | 199 +++++- openmc/mesh.py | 337 ++++++++++ src/mesh.cpp | 1432 ++++++++++++++++++++++------------------- 4 files changed, 1469 insertions(+), 687 deletions(-) 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)); From 1c8a94f41eec3e77edd2088afa595cd15ef9f68c Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Mon, 20 Dec 2021 19:42:56 +0100 Subject: [PATCH 02/17] Update Documentation to reflect new meshes --- docs/source/io_formats/tallies.rst | 15 ++++++++++++--- src/relaxng/tallies.rnc | 18 ++++++++++++++++++ 2 files changed, 30 insertions(+), 3 deletions(-) diff --git a/docs/source/io_formats/tallies.rst b/docs/source/io_formats/tallies.rst index 7f95f29850..50ba7b3d5a 100644 --- a/docs/source/io_formats/tallies.rst +++ b/docs/source/io_formats/tallies.rst @@ -312,8 +312,8 @@ element with the tag name ````. This element has the following attributes/sub-elements: :type: - The type of mesh. This can be either "regular", "rectilinear", or - "unstructured". + The type of mesh. This can be either "regular", "rectilinear", + "cylindrical", "spherical", or "unstructured". :dimension: The number of mesh cells in each direction. (For regular mesh only.) @@ -340,7 +340,16 @@ attributes/sub-elements: The mesh divisions along the y-axis. (For rectilinear mesh only.) :z_grid: - The mesh divisions along the z-axis. (For rectilinear mesh only.) + The mesh divisions along the z-axis. (For rectilinear and cylindrical meshes only.) + + :r_grid: + The mesh divisions along the r-axis. (For cylindrical and spherical meshes only.) + + :p_grid: + The mesh divisions along the phi-axis. (For cylindrical and spherical meshes only.) + + :t_grid: + The mesh divisions along the theta-axis. (For spherical mesh only.) :library: The mesh library used to represent an unstructured mesh. This can be either diff --git a/src/relaxng/tallies.rnc b/src/relaxng/tallies.rnc index 3d504a66c8..583b274c60 100644 --- a/src/relaxng/tallies.rnc +++ b/src/relaxng/tallies.rnc @@ -24,6 +24,24 @@ element tallies { attribute y_grid { list { xsd:double+ } }) & (element z_grid { list { xsd:double+ } } | attribute z_grid { list { xsd:double+ } }) + ) | ( + (element type { ( "cylindrical" ) } | + attribute type { ( "cylindrical" ) }) & + (element r_grid { list { xsd:double+ } } | + attribute r_grid { list { xsd:double+ } }) & + (element p_grid { list { xsd:double+ } } | + attribute p_grid { list { xsd:double+ } }) & + (element z_grid { list { xsd:double+ } } | + attribute z_grid { list { xsd:double+ } }) + ) | ( + (element type { ( "spherical" ) } | + attribute type { ( "spherical" ) }) & + (element r_grid { list { xsd:double+ } } | + attribute r_grid { list { xsd:double+ } }) & + (element t_grid { list { xsd:double+ } } | + attribute p_grid { list { xsd:double+ } }) & + (element p_grid { list { xsd:double+ } } | + attribute z_grid { list { xsd:double+ } }) ) ) }* & From efe36357a87c9a960d3d01f2eff1c6ff1ebd9f41 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Mon, 20 Dec 2021 19:43:40 +0100 Subject: [PATCH 03/17] Update tests to check new meshes --- .../filter_mesh/inputs_true.dat | 54 +++++++++++++--- .../filter_mesh/results_true.dat | 2 +- tests/regression_tests/filter_mesh/test.py | 18 +++++- tests/unit_tests/test_lib.py | 62 +++++++++++++++++++ 4 files changed, 125 insertions(+), 11 deletions(-) diff --git a/tests/regression_tests/filter_mesh/inputs_true.dat b/tests/regression_tests/filter_mesh/inputs_true.dat index 37716c02e4..38b20b7e33 100644 --- a/tests/regression_tests/filter_mesh/inputs_true.dat +++ b/tests/regression_tests/filter_mesh/inputs_true.dat @@ -53,36 +53,58 @@ -7.5 -6.617647058823529 -5.735294117647059 -4.852941176470589 -3.9705882352941178 -3.0882352941176467 -2.2058823529411766 -1.3235294117647065 -0.4411764705882355 0.4411764705882355 1.3235294117647065 2.2058823529411757 3.0882352941176467 3.9705882352941178 4.852941176470587 5.735294117647058 6.617647058823529 7.5 1.0 1.223224374241637 1.4962778697388448 1.8302835609029084 2.2388474634702153 2.7386127875258306 3.3499379133114306 4.09772570775871 5.012437964687018 6.131336292779302 7.500000000000001 + + 0.0 0.4411764705882353 0.8823529411764706 1.3235294117647058 1.7647058823529411 2.2058823529411766 2.6470588235294117 3.0882352941176467 3.5294117647058822 3.9705882352941178 4.411764705882353 4.852941176470588 5.294117647058823 5.735294117647059 6.1764705882352935 6.617647058823529 7.0588235294117645 7.5 + 0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0 180.0 200.0 220.0 240.0 260.0 280.0 300.0 320.0 340.0 360.0 + -7.5 -6.5625 -5.625 -4.6875 -3.75 -2.8125 -1.875 -0.9375 0.0 0.9375 1.875 2.8125 3.75 4.6875 5.625 6.5625 7.5 + + + 0.0 0.4411764705882353 0.8823529411764706 1.3235294117647058 1.7647058823529411 2.2058823529411766 2.6470588235294117 3.0882352941176467 3.5294117647058822 3.9705882352941178 4.411764705882353 4.852941176470588 5.294117647058823 5.735294117647059 6.1764705882352935 6.617647058823529 7.0588235294117645 7.5 + 0.0 22.5 45.0 67.5 90.0 112.5 135.0 157.5 180.0 + 0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0 180.0 200.0 220.0 240.0 260.0 280.0 300.0 320.0 340.0 360.0 + 1 - + 1 2 - + 2 3 - + 3 4 - + 4 + + 5 + + + 5 + + + 6 + + + 6 + 1 total - 5 + 7 current @@ -90,7 +112,7 @@ total - 6 + 8 current @@ -98,7 +120,7 @@ total - 7 + 9 current @@ -106,7 +128,23 @@ total - 8 + 10 + current + + + 5 + total + + + 11 + current + + + 6 + total + + + 12 current diff --git a/tests/regression_tests/filter_mesh/results_true.dat b/tests/regression_tests/filter_mesh/results_true.dat index 332c8b36cf..ec9fe47a3c 100644 --- a/tests/regression_tests/filter_mesh/results_true.dat +++ b/tests/regression_tests/filter_mesh/results_true.dat @@ -1 +1 @@ -97cf099a24099d2af794752e145841a9ef9f24980c53ee71e3fdd400d4249a7e6f3f52f63132fe2fe8d1973604a72d87910ab13bc3721f7d0bec533b87756847 \ No newline at end of file +cacc0cedfa243b31f6c44e22bedea400744af14ad139631ec7141c0812ade18e1b53f59e65865ef90d296716a555718cc06108db9c2d434d9d04ec9bb878a81d \ No newline at end of file diff --git a/tests/regression_tests/filter_mesh/test.py b/tests/regression_tests/filter_mesh/test.py index c8aa871a83..eddf0a9913 100644 --- a/tests/regression_tests/filter_mesh/test.py +++ b/tests/regression_tests/filter_mesh/test.py @@ -51,18 +51,32 @@ def model(): recti_mesh.y_grid = np.linspace(-7.5, 7.5, 18) recti_mesh.z_grid = np.logspace(0, np.log10(7.5), 11) + cyl_mesh = openmc.CylindricalMesh() + cyl_mesh.r_grid = np.linspace(0, 7.5, 18) + cyl_mesh.phi_grid = np.linspace(0, 360, 19) + cyl_mesh.z_grid = np.linspace(-7.5, 7.5, 17) + + sph_mesh = openmc.SphericalMesh() + sph_mesh.r_grid = np.linspace(0, 7.5, 18) + sph_mesh.theta_grid = np.linspace(0, 180, 9) + sph_mesh.phi_grid = np.linspace(0, 360, 19) + # Create filters reg_filters = [ openmc.MeshFilter(mesh_1d), openmc.MeshFilter(mesh_2d), openmc.MeshFilter(mesh_3d), - openmc.MeshFilter(recti_mesh) + openmc.MeshFilter(recti_mesh), + openmc.MeshFilter(cyl_mesh), + openmc.MeshFilter(sph_mesh) ] surf_filters = [ openmc.MeshSurfaceFilter(mesh_1d), openmc.MeshSurfaceFilter(mesh_2d), openmc.MeshSurfaceFilter(mesh_3d), - openmc.MeshSurfaceFilter(recti_mesh) + openmc.MeshSurfaceFilter(recti_mesh), + openmc.MeshSurfaceFilter(cyl_mesh), + openmc.MeshSurfaceFilter(sph_mesh) ] # Create tallies diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index ccff33ac39..3572e9863a 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -580,6 +580,68 @@ def test_rectilinear_mesh(lib_init): msf = openmc.lib.MeshSurfaceFilter(mesh) assert msf.mesh == mesh +def test_cylindrical_mesh(lib_init): + deg2rad = lambda deg: deg*np.pi/180 + mesh = openmc.lib.CylindricalMesh() + x_grid = [0., 5., 10.] + y_grid = [0., 10., 20.] + z_grid = [10., 20., 30.] + mesh.set_grid(x_grid, y_grid, z_grid) + assert np.all(mesh.lower_left == (0., 0., 10.)) + assert np.all(mesh.upper_right == (10., deg2rad(20.), 30.)) + assert np.all(mesh.dimension == (2, 2, 2)) + for i, diff_x in enumerate(np.diff(x_grid)): + for j, diff_y in enumerate(np.diff(y_grid)): + for k, diff_z in enumerate(np.diff(z_grid)): + assert np.all(mesh.width[i, j, k, :] == (5, deg2rad(10), 10)) + + with pytest.raises(exc.AllocationError): + mesh2 = openmc.lib.CylindricalMesh(mesh.id) + + meshes = openmc.lib.meshes + assert isinstance(meshes, Mapping) + assert len(meshes) == 3 + + mesh = meshes[mesh.id] + assert isinstance(mesh, openmc.lib.CylindricalMesh) + + mf = openmc.lib.MeshFilter(mesh) + assert mf.mesh == mesh + + msf = openmc.lib.MeshSurfaceFilter(mesh) + assert msf.mesh == mesh + +def test_spherical_mesh(lib_init): + deg2rad = lambda deg: deg*np.pi/180 + mesh = openmc.lib.SphericalMesh() + x_grid = [0., 5., 10.] + y_grid = [0., 10., 20.] + z_grid = [10., 20., 30.] + mesh.set_grid(x_grid, y_grid, z_grid) + assert np.all(mesh.lower_left == (0., 0., deg2rad(10.))) + assert np.all(mesh.upper_right == (10., deg2rad(20.), deg2rad(30.))) + assert np.all(mesh.dimension == (2, 2, 2)) + for i, diff_x in enumerate(np.diff(x_grid)): + for j, diff_y in enumerate(np.diff(y_grid)): + for k, diff_z in enumerate(np.diff(z_grid)): + assert np.all(abs(mesh.width[i, j, k, :] - (5, deg2rad(10), deg2rad(10))) < 1e-16) + + with pytest.raises(exc.AllocationError): + mesh2 = openmc.lib.SphericalMesh(mesh.id) + + meshes = openmc.lib.meshes + assert isinstance(meshes, Mapping) + assert len(meshes) == 4 + + mesh = meshes[mesh.id] + assert isinstance(mesh, openmc.lib.SphericalMesh) + + mf = openmc.lib.MeshFilter(mesh) + assert mf.mesh == mesh + + msf = openmc.lib.MeshSurfaceFilter(mesh) + assert msf.mesh == mesh + def test_restart(lib_init, mpi_intracomm): # Finalize and re-init to make internal state consistent with XML. From f92993d4433606d3fc456909c30dc45855837b0e Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Wed, 22 Dec 2021 09:24:13 +0100 Subject: [PATCH 04/17] Update tests with changed results.results_true --- .../filter_mesh/results_true.dat | 2 +- .../mgxs_library_mesh/results_true.dat | 16 ++++++++-------- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/tests/regression_tests/filter_mesh/results_true.dat b/tests/regression_tests/filter_mesh/results_true.dat index ec9fe47a3c..60b3f1928b 100644 --- a/tests/regression_tests/filter_mesh/results_true.dat +++ b/tests/regression_tests/filter_mesh/results_true.dat @@ -1 +1 @@ -cacc0cedfa243b31f6c44e22bedea400744af14ad139631ec7141c0812ade18e1b53f59e65865ef90d296716a555718cc06108db9c2d434d9d04ec9bb878a81d \ No newline at end of file +79f82f287ff427b896ab8aedc81f4671949571316b90e89806f0753ed0628dd5ac32f8d0e4a3f5f5f5c65a9e8a06a6869eb5fc896b4bb423870877afca2a7d4f \ No newline at end of file diff --git a/tests/regression_tests/mgxs_library_mesh/results_true.dat b/tests/regression_tests/mgxs_library_mesh/results_true.dat index 0acc2f407c..440c75a213 100644 --- a/tests/regression_tests/mgxs_library_mesh/results_true.dat +++ b/tests/regression_tests/mgxs_library_mesh/results_true.dat @@ -181,11 +181,11 @@ mesh 1 group in nuclide mean std. dev. x y surf 3 1 1 x-max in 1 total 0.1816 0.006129 -2 1 1 x-max out 1 total 0.2718 0.095077 +2 1 1 x-max out 1 total 0.1730 0.004382 1 1 1 x-min in 1 total 0.0000 0.000000 0 1 1 x-min out 1 total 0.0000 0.000000 7 1 1 y-max in 1 total 0.1790 0.011921 -6 1 1 y-max out 1 total 0.2290 0.034821 +6 1 1 y-max out 1 total 0.1802 0.016599 5 1 1 y-min in 1 total 0.0000 0.000000 4 1 1 y-min out 1 total 0.0000 0.000000 19 1 2 x-max in 1 total 0.1872 0.012447 @@ -194,23 +194,23 @@ 16 1 2 x-min out 1 total 0.0000 0.000000 23 1 2 y-max in 1 total 0.0000 0.000000 22 1 2 y-max out 1 total 0.0000 0.000000 -21 1 2 y-min in 1 total 0.2290 0.034821 +21 1 2 y-min in 1 total 0.1802 0.016599 20 1 2 y-min out 1 total 0.1790 0.011921 11 2 1 x-max in 1 total 0.0000 0.000000 10 2 1 x-max out 1 total 0.0000 0.000000 -9 2 1 x-min in 1 total 0.2718 0.095077 +9 2 1 x-min in 1 total 0.1730 0.004382 8 2 1 x-min out 1 total 0.1816 0.006129 15 2 1 y-max in 1 total 0.1778 0.009484 -14 2 1 y-max out 1 total 0.2326 0.042782 +14 2 1 y-max out 1 total 0.1822 0.010077 13 2 1 y-min in 1 total 0.0000 0.000000 12 2 1 y-min out 1 total 0.0000 0.000000 27 2 2 x-max in 1 total 0.0000 0.000000 -26 2 2 x-max out 1 total 0.0260 0.026000 +26 2 2 x-max out 1 total 0.0000 0.000000 25 2 2 x-min in 1 total 0.1952 0.015948 24 2 2 x-min out 1 total 0.1872 0.012447 31 2 2 y-max in 1 total 0.0000 0.000000 -30 2 2 y-max out 1 total 0.0244 0.024400 -29 2 2 y-min in 1 total 0.2326 0.042782 +30 2 2 y-max out 1 total 0.0000 0.000000 +29 2 2 y-min in 1 total 0.1822 0.010077 28 2 2 y-min out 1 total 0.1778 0.009484 mesh 1 group in nuclide mean std. dev. x y z From 776c6baf1fa81cf3845df0d22d38c8bd00536c64 Mon Sep 17 00:00:00 2001 From: Olaf Jochen Schumann Date: Wed, 22 Dec 2021 16:19:58 +0100 Subject: [PATCH 05/17] Changes request by @gridley --- include/openmc/mesh.h | 57 +++- openmc/lib/mesh.py | 4 +- src/eigenvalue.cpp | 2 +- src/mesh.cpp | 288 ++++++++---------- .../filter_mesh/results_true.dat | 2 +- 5 files changed, 177 insertions(+), 176 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index c6a870f8d5..f6908be432 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -123,6 +123,9 @@ public: //! \param[in] bin Mesh bin to generate a label for virtual std::string bin_label(int bin) const = 0; + //! Return the mesh type + virtual std::string get_mesh_type() const = 0; + // Data members int id_ {-1}; //!< User-specified ID int n_dimension_; //!< Number of dimensions @@ -135,9 +138,12 @@ public: virtual ~StructuredMesh() = default; using MeshIndex = std::array; + struct MeshDistance { - MeshDistance() {}; - MeshDistance(int _index, bool _maxSurface, double _distance): nextIndex{_index}, maxSurface(_maxSurface), distance{_distance} {}; + MeshDistance() = default; + MeshDistance(int _index, bool _maxSurface, double _distance): + nextIndex{_index}, maxSurface{_maxSurface}, distance{_distance} + { } int nextIndex { -1 }; double distance { INFTY }; bool maxSurface { true }; @@ -185,13 +191,13 @@ public: // //! \param[in] r Position to get indices for //! \param[out] in_mesh Whether position is in mesh - //! \return ijk Array of mesh indices + //! \return Array of mesh indices virtual MeshIndex get_indices(Position r, bool& in_mesh) const; //! Get mesh indices corresponding to a mesh bin // //! \param[in] bin Mesh bin - //! \param[out] ijk Mesh indices + //! \return ijk Mesh indices virtual MeshIndex get_indices_from_bin(int bin) const; //! Get mesh index in a particular direction @@ -200,19 +206,30 @@ public: //! \param[in] i Direction index virtual int get_index_in_direction(double r, int i) const = 0; - //! Get the closest distance from to coordinate r to the grid for the mesh grid boundary in the negative direction + //! Get the closest distance from the coordinate r to the grid surface + //! in i direction that bounds mesh cell ijk and that is larger than l + //! The coordinate r does not have to be inside the mesh cell ijk. In + //! curved coordinates, multiple crossings of the same surface can happen, + //! these are selected by the parameter l //! //! \param[in] ijk Array of mesh indices - //! \param[in] i Direction index + //! \param[in] i direction index of grid surface + //! \param[in] r0 position, from where to calculate the distance + //! \param[in] u direction of flight. actual position is r0 + l * u + //! \param[in] l actual chord length + //! \return MeshDistance struct with closest distance, next cell index in i-direction and min/max surface indicator virtual MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const = 0; //! Get a label for the mesh bin std::string bin_label(int bin) const override; + //! Get shape as xt::xtensor + xt::xtensor get_x_shape() const; + // Data members xt::xtensor lower_left_; //!< Lower-left coordinates of mesh xt::xtensor upper_right_; //!< Upper-right coordinates of mesh - xt::xtensor shape_; //!< Number of mesh elements in each dimension + std::array shape_; //!< Number of mesh elements in each dimension protected: @@ -231,7 +248,9 @@ public: // Overridden methods int get_index_in_direction(double r, int i) const override; + virtual std::string get_mesh_type() const override; + static const std::string mesh_type; MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override; @@ -274,6 +293,10 @@ public: // Overridden methods int get_index_in_direction(double r, int i) const override; + + virtual std::string get_mesh_type() const override; + + static const std::string mesh_type; MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override; @@ -296,7 +319,7 @@ public: double negative_grid_boundary(const MeshIndex& ijk, int i) const; - vector> grid_; + array, 3> grid_; int set_grid(); }; @@ -312,6 +335,10 @@ public: int get_index_in_direction(double r, int i) const override; + virtual std::string get_mesh_type() const override; + + static const std::string mesh_type; + MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override; std::pair, vector> plot( @@ -319,7 +346,7 @@ public: void to_hdf5(hid_t group) const override; - vector> grid_; + array, 3> grid_; int set_grid(); @@ -358,6 +385,10 @@ public: int get_index_in_direction(double r, int i) const override; + virtual std::string get_mesh_type() const override; + + static const std::string mesh_type; + MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override; std::pair, vector> plot( @@ -365,7 +396,7 @@ public: void to_hdf5(hid_t group) const override; - vector> grid_; + array, 3> grid_; int set_grid(); @@ -405,6 +436,8 @@ public: UnstructuredMesh(pugi::xml_node node); UnstructuredMesh(const std::string& filename); + static const std::string mesh_type; + // Overridden Methods void surface_bins_crossed(Position r0, Position r1, const Direction& u, vector& bins) const override; @@ -476,6 +509,8 @@ public: MOABMesh(const std::string& filename, double length_multiplier = 1.0); MOABMesh(std::shared_ptr external_mbi); + static const std::string mesh_lib_type; + // Overridden Methods void bins_crossed(Position r0, Position r1, const Direction& u, @@ -623,6 +658,8 @@ public: LibMesh(pugi::xml_node node); LibMesh(const std::string& filename, double length_multiplier = 1.0); + static const std::string mesh_lib_type; + // Overridden Methods void bins_crossed(Position r0, Position r1, const Direction& u, vector& bins, vector& lengths) const override; diff --git a/openmc/lib/mesh.py b/openmc/lib/mesh.py index 82545e2bf1..43f66d7a12 100644 --- a/openmc/lib/mesh.py +++ b/openmc/lib/mesh.py @@ -471,8 +471,8 @@ class SphericalMesh(Mesh): _MESH_TYPE_MAP = { 'regular': RegularMesh, 'rectilinear': RectilinearMesh, - 'cylindrical': RectilinearMesh, - 'spherical': RectilinearMesh + 'cylindrical': CylindricalMesh, + 'spherical': SphericalMesh } diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index b35281f4ad..5584210a1f 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -577,7 +577,7 @@ void ufs_count_sites() #ifdef OPENMC_MPI // Send source fraction to all processors - int n_bins = xt::prod(simulation::ufs_mesh->shape_)(); + int n_bins = simulation::ufs_mesh->n_bins(); MPI_Bcast( simulation::source_frac.data(), n_bins, MPI_DOUBLE, 0, mpi::intracomm); #endif diff --git a/src/mesh.cpp b/src/mesh.cpp index d4948fb864..bd139e94be 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -32,8 +32,6 @@ #include "openmc/tallies/tally.h" #include "openmc/xml_interface.h" -#include - #ifdef LIBMESH #include "libmesh/mesh_tools.h" #include "libmesh/numeric_vector.h" @@ -158,6 +156,10 @@ std::string StructuredMesh::bin_label(int bin) const } } +xt::xtensor StructuredMesh::get_x_shape() const { + return xt::adapt(shape_, { n_dimension_}); +} + //============================================================================== // Unstructured Mesh implementation //============================================================================== @@ -168,7 +170,7 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node) // check the mesh type if (check_for_node(node, "type")) { auto temp = get_node_value(node, "type", true, true); - if (temp != "unstructured") { + if (temp != mesh_type) { fatal_error(fmt::format("Invalid mesh type: {}", temp)); } } @@ -197,6 +199,8 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node) } } +const std::string UnstructuredMesh::mesh_type = "unstructured"; + void UnstructuredMesh::surface_bins_crossed( Position r0, Position r1, const Direction& u, vector& bins) const { @@ -212,7 +216,7 @@ void UnstructuredMesh::to_hdf5(hid_t group) const { hid_t mesh_group = create_group(group, fmt::format("mesh {}", id_)); - write_dataset(mesh_group, "type", "unstructured"); + write_dataset(mesh_group, "type", mesh_type); write_dataset(mesh_group, "filename", filename_); write_dataset(mesh_group, "library", this->library()); // write volume of each element @@ -298,7 +302,7 @@ int StructuredMesh::get_bin(Position r) const int StructuredMesh::n_bins() const { - return xt::prod(shape_)(); + return std::accumulate(shape_.begin(), shape_.begin() + n_dimension_, 1, std::multiplies<>()); } int StructuredMesh::n_surface_bins() const @@ -428,7 +432,7 @@ void StructuredMesh::raytrace_mesh(Position r0, Position r1, const Direction& u, 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 + // Check if we have left the interior 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 @@ -524,14 +528,15 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node} fatal_error("Must specify on a regular mesh."); } - shape_ = get_node_xarray(node, "dimension"); - int n = n_dimension_ = shape_.size(); + xt::xtensor shape = get_node_xarray(node, "dimension"); + int n = n_dimension_ = shape.size(); if (n != 1 && n != 2 && n != 3) { fatal_error("Mesh must be one, two, or three dimensions."); } + std::copy(shape.begin(), shape.end(), shape_.begin()); // Check that dimensions are all greater than zero - if (xt::any(shape_ <= 0)) { + if (xt::any(shape <= 0)) { fatal_error("All entries on the element for a tally " "mesh must be positive."); } @@ -545,7 +550,7 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node} } // Make sure lower_left and dimension match - if (shape_.size() != lower_left_.size()) { + if (shape.size() != lower_left_.size()) { fatal_error("Number of entries on must be the same " "as the number of entries on ."); } @@ -571,7 +576,7 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node} } // Set width and upper right coordinate - upper_right_ = xt::eval(lower_left_ + shape_ * width_); + upper_right_ = xt::eval(lower_left_ + shape * width_); } else if (check_for_node(node, "upper_right")) { upper_right_ = get_node_xarray(node, "upper_right"); @@ -590,13 +595,13 @@ RegularMesh::RegularMesh(pugi::xml_node node) : StructuredMesh {node} } // Set width - width_ = xt::eval((upper_right_ - lower_left_) / shape_); + width_ = xt::eval((upper_right_ - lower_left_) / shape); } else { fatal_error("Must specify either or on a mesh."); } // Set volume fraction - volume_frac_ = 1.0 / xt::prod(shape_)(); + volume_frac_ = 1.0 / xt::prod(shape)(); } int RegularMesh::get_index_in_direction(double r, int i) const @@ -604,6 +609,12 @@ int RegularMesh::get_index_in_direction(double r, int i) const return std::ceil((r - lower_left_[i]) / width_[i]); } +const std::string RegularMesh::mesh_type = "regular"; + +std::string RegularMesh::get_mesh_type() const { + return mesh_type; +} + double RegularMesh::positive_grid_boundary(const MeshIndex& ijk, int i) const { return lower_left_[i] + ijk[i] * width_[i]; @@ -679,7 +690,7 @@ void RegularMesh::to_hdf5(hid_t group) const hid_t mesh_group = create_group(group, "mesh " + std::to_string(id_)); write_dataset(mesh_group, "type", "regular"); - write_dataset(mesh_group, "dimension", shape_); + write_dataset(mesh_group, "dimension", get_x_shape()); write_dataset(mesh_group, "lower_left", lower_left_); write_dataset(mesh_group, "upper_right", upper_right_); write_dataset(mesh_group, "width", width_); @@ -752,7 +763,6 @@ RectilinearMesh::RectilinearMesh(pugi::xml_node node) : StructuredMesh {node} { n_dimension_ = 3; - grid_.resize(n_dimension_); grid_[0] = get_node_array(node, "x_grid"); grid_[1] = get_node_array(node, "y_grid"); grid_[2] = get_node_array(node, "z_grid"); @@ -762,6 +772,12 @@ RectilinearMesh::RectilinearMesh(pugi::xml_node node) : StructuredMesh {node} } } +const std::string RectilinearMesh::mesh_type = "rectilinear"; + +std::string RectilinearMesh::get_mesh_type() const { + return mesh_type; +} + double RectilinearMesh::positive_grid_boundary(const MeshIndex& ijk, int i) const { return grid_[i][ijk[i]]; @@ -802,12 +818,10 @@ int RectilinearMesh::set_grid() "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 x-, y-, and z- grids for " - "rectilinear meshes must be sorted and unique."); - return OPENMC_E_INVALID_ARGUMENT; - } + if (not std::is_sorted(g.begin(), g.end(), std::less_equal())) { + set_errmsg("Values in for x-, y-, and z- grids for " + "rectilinear meshes must be sorted and unique."); + return OPENMC_E_INVALID_ARGUMENT; } } @@ -817,8 +831,6 @@ int RectilinearMesh::set_grid() return 0; } - - int RectilinearMesh::get_index_in_direction(double r, int i) const { return lower_bound_index(grid_[i].begin(), grid_[i].end(), r) + 1; @@ -874,7 +886,6 @@ 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"); @@ -884,6 +895,12 @@ CylindricalMesh::CylindricalMesh(pugi::xml_node node) : StructuredMesh {node} } } +const std::string CylindricalMesh::mesh_type = "cylindrical"; + +std::string CylindricalMesh::get_mesh_type() const { + return mesh_type; +} + StructuredMesh::MeshIndex CylindricalMesh::get_indices(Position r, bool& in_mesh) const { Position mapped_r; @@ -1020,12 +1037,10 @@ int CylindricalMesh::set_grid() "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 (not std::is_sorted(g.begin(), g.end(), std::less_equal())) { + 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) { @@ -1093,7 +1108,6 @@ 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"); @@ -1103,6 +1117,12 @@ SphericalMesh::SphericalMesh(pugi::xml_node node) : StructuredMesh {node} } } +const std::string SphericalMesh::mesh_type = "spherical"; + +std::string SphericalMesh::get_mesh_type() const { + return mesh_type; +} + StructuredMesh::MeshIndex SphericalMesh::get_indices(Position r, bool& in_mesh) const { Position mapped_r; @@ -1128,8 +1148,6 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices(Position r, bool& in_mesh) 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 @@ -1275,12 +1293,10 @@ int SphericalMesh::set_grid() "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 (not std::is_sorted(g.begin(), g.end(), std::less_equal())) { + 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 " @@ -1333,7 +1349,7 @@ 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, "type", SphericalMesh::mesh_type); write_dataset(mesh_group, "r_grid", grid_[0]); write_dataset(mesh_group, "t_grid", grid_[1]); write_dataset(mesh_group, "p_grid", grid_[2]); @@ -1388,15 +1404,8 @@ extern "C" int openmc_mesh_get_type(int32_t index, char* type) if (int err = check_mesh(index)) return err; - if (is_mesh_type(index)) { - std::strcpy(type, "regular"); - } 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"); - } + std::strcpy(type, model::meshes[index].get()->get_mesh_type().c_str()); + return 0; } @@ -1409,13 +1418,14 @@ extern "C" int openmc_extend_meshes( std::string mesh_type; for (int i = 0; i < n; ++i) { - if (std::strcmp(type, "regular") == 0) { + //if (std::strcmp(type, RegularMesh::mesh_type.c_str()) == 0) { + if (RegularMesh::mesh_type == type) { model::meshes.push_back(make_unique()); - } else if (std::strcmp(type, "rectilinear") == 0) { + } else if (RectilinearMesh::mesh_type == type) { model::meshes.push_back(make_unique()); - } else if (std::strcmp(type, "cylindrical") == 0) { + } else if (CylindricalMesh::mesh_type == type) { model::meshes.push_back(make_unique()); - } else if (std::strcmp(type, "spherical") == 0) { + } else if (SphericalMesh::mesh_type == type) { model::meshes.push_back(make_unique()); } else { throw std::runtime_error {"Unknown mesh type: " + std::string(type)}; @@ -1436,14 +1446,14 @@ extern "C" int openmc_add_unstructured_mesh( bool valid_lib = false; #ifdef DAGMC - if (lib_name == "moab") { + if (lib_name == MOABMesh::mesh_lib_type) { model::meshes.push_back(std::move(make_unique(mesh_file))); valid_lib = true; } #endif #ifdef LIBMESH - if (lib_name == "libmesh") { + if (lib_name == LibMesh::mesh_lib_type) { model::meshes.push_back(std::move(make_unique(mesh_file))); valid_lib = true; } @@ -1515,9 +1525,8 @@ extern "C" int openmc_regular_mesh_set_dimension( RegularMesh* mesh = dynamic_cast(model::meshes[index].get()); // Copy dimension - vector shape = {static_cast(n)}; - mesh->shape_ = xt::adapt(dims, n, xt::no_ownership(), shape); - mesh->n_dimension_ = mesh->shape_.size(); + mesh->n_dimension_ = n; + std::copy(dims, dims + n, mesh->shape_.begin()); return 0; } @@ -1553,15 +1562,15 @@ extern "C" int openmc_regular_mesh_set_params( if (ll && ur) { m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape); m->upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape); - m->width_ = (m->upper_right_ - m->lower_left_) / m->shape_; + m->width_ = (m->upper_right_ - m->lower_left_) / m->get_x_shape(); } else if (ll && width) { m->lower_left_ = xt::adapt(ll, n, xt::no_ownership(), shape); m->width_ = xt::adapt(width, n, xt::no_ownership(), shape); - m->upper_right_ = m->lower_left_ + m->shape_ * m->width_; + m->upper_right_ = m->lower_left_ + m->get_x_shape() * m->width_; } else if (ur && width) { m->upper_right_ = xt::adapt(ur, n, xt::no_ownership(), shape); m->width_ = xt::adapt(width, n, xt::no_ownership(), shape); - m->lower_left_ = m->upper_right_ - m->shape_ * m->width_; + m->lower_left_ = m->upper_right_ - m->get_x_shape() * m->width_; } else { set_errmsg("At least two parameters must be specified."); return OPENMC_E_INVALID_ARGUMENT; @@ -1570,42 +1579,23 @@ extern "C" int openmc_regular_mesh_set_params( return 0; } -//! Get the rectilinear mesh grid -extern "C" int openmc_rectilinear_mesh_get_grid(int32_t index, double** grid_x, - int* nx, double** grid_y, int* ny, double** grid_z, int* nz) -{ - if (int err = check_mesh_type(index)) - return err; - RectilinearMesh* m = - dynamic_cast(model::meshes[index].get()); - if (m->lower_left_.dimension() == 0) { - set_errmsg("Mesh parameters have not been set."); - return OPENMC_E_ALLOCATE; - } - - *grid_x = m->grid_[0].data(); - *nx = m->grid_[0].size(); - *grid_y = m->grid_[1].data(); - *ny = m->grid_[1].size(); - *grid_z = m->grid_[2].data(); - *nz = m->grid_[2].size(); - - return 0; -} - -//! Set the rectilienar mesh parameters -extern "C" int openmc_rectilinear_mesh_set_grid(int32_t index, +//! Set the mesh parameters for rectilinear, cylindrical and spharical meshes +template +int openmc_structured_mesh_set_grid_impl(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)) + if (int err = check_mesh_type(index)) return err; - RectilinearMesh* m = - dynamic_cast(model::meshes[index].get()); + + C* m = dynamic_cast(model::meshes[index].get()); m->n_dimension_ = 3; - m->grid_.resize(m->n_dimension_); + + m->grid_[0].reserve(nx); + m->grid_[1].reserve(ny); + m->grid_[2].reserve(nz); for (int i = 0; i < nx; i++) { m->grid_[0].push_back(grid_x[i]); @@ -1618,17 +1608,17 @@ extern "C" int openmc_rectilinear_mesh_set_grid(int32_t index, } int err = m->set_grid(); - return err; + return err; } -//! Get the cylindrical mesh grid -extern "C" int openmc_cylindrical_mesh_get_grid(int32_t index, double** grid_x, +//! Get the mesh parameters for rectilinear, cylindrical and spharical meshes +template +int openmc_structured_mesh_get_grid_impl(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)) + if (int err = check_mesh_type(index)) return err; - CylindricalMesh* m = - dynamic_cast(model::meshes[index].get()); + C* m = dynamic_cast(model::meshes[index].get()); if (m->lower_left_.dimension() == 0) { set_errmsg("Mesh parameters have not been set."); @@ -1645,55 +1635,44 @@ extern "C" int openmc_cylindrical_mesh_get_grid(int32_t index, double** grid_x, return 0; } + +//! Get the rectilinear mesh grid +extern "C" int openmc_rectilinear_mesh_get_grid(int32_t index, double** grid_x, + int* nx, double** grid_y, int* ny, double** grid_z, int* nz) +{ + return openmc_structured_mesh_get_grid_impl(index, grid_x, nx, grid_y, ny, grid_z, nz); +} + + +//! Set the rectilienar mesh parameters +extern "C" int openmc_rectilinear_mesh_set_grid(int32_t index, + const double* grid_x, const int nx, const double* grid_y, const int ny, + const double* grid_z, const int nz) +{ + return openmc_structured_mesh_set_grid_impl(index, grid_x, nx, grid_y, ny, grid_z, nz); +} + +//! 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) +{ + return openmc_structured_mesh_get_grid_impl(index, grid_x, nx, grid_y, ny, grid_z, nz); +} + //! 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; + return openmc_structured_mesh_set_grid_impl(index, grid_x, nx, grid_y, ny, grid_z, nz); } //! 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; + return openmc_structured_mesh_get_grid_impl(index, grid_x, nx, grid_y, ny, grid_z, nz);; } //! Set the spherical mesh parameters @@ -1701,30 +1680,13 @@ 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; + return openmc_structured_mesh_set_grid_impl(index, grid_x, nx, grid_y, ny, grid_z, nz); } #ifdef DAGMC +const std::string MOABMesh::mesh_lib_type = "moab"; + MOABMesh::MOABMesh(pugi::xml_node node) : UnstructuredMesh(node) { initialize(); @@ -1994,7 +1956,7 @@ double MOABMesh::volume(int bin) const std::string MOABMesh::library() const { - return "moab"; + return MOABMesh::mesh_lib_type; } double MOABMesh::tet_volume(moab::EntityHandle tet) const @@ -2313,6 +2275,8 @@ void MOABMesh::write(const std::string& base_filename) const #ifdef LIBMESH +const std::string LibMesh::mesh_lib_type = "libmesh"; + LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node) { initialize(); @@ -2387,7 +2351,7 @@ Position LibMesh::centroid(int bin) const std::string LibMesh::library() const { - return "libmesh"; + return LibMesh::mesh_lib_type; } int LibMesh::n_bins() const @@ -2564,23 +2528,23 @@ void read_meshes(pugi::xml_node root) } // Read mesh and add to vector - if (mesh_type == "regular") { + if (mesh_type == RegularMesh::mesh_type) { model::meshes.push_back(make_unique(node)); - } else if (mesh_type == "rectilinear") { + } else if (mesh_type == RectilinearMesh::mesh_type) { model::meshes.push_back(make_unique(node)); - } else if (mesh_type == "cylindrical") { + } else if (mesh_type == CylindricalMesh::mesh_type) { model::meshes.push_back(make_unique(node)); - } else if (mesh_type == "spherical") { + } else if (mesh_type == SphericalMesh::mesh_type) { model::meshes.push_back(make_unique(node)); #ifdef DAGMC - } else if (mesh_type == "unstructured" && mesh_lib == "moab") { + } else if (mesh_type == UnstructuredMesh::mesh_type && mesh_lib == MOABMesh::mesh_lib_type) { model::meshes.push_back(make_unique(node)); #endif #ifdef LIBMESH - } else if (mesh_type == "unstructured" && mesh_lib == "libmesh") { + } else if (mesh_type == UnstructuredMesh::mesh_type && mesh_lib == LibMesh::mesh_lib_type) { model::meshes.push_back(make_unique(node)); #endif - } else if (mesh_type == "unstructured") { + } else if (mesh_type == UnstructuredMesh::mesh_type) { fatal_error("Unstructured mesh support is not enabled or the mesh " "library is invalid."); } else { diff --git a/tests/regression_tests/filter_mesh/results_true.dat b/tests/regression_tests/filter_mesh/results_true.dat index 60b3f1928b..4247f6ed1c 100644 --- a/tests/regression_tests/filter_mesh/results_true.dat +++ b/tests/regression_tests/filter_mesh/results_true.dat @@ -1 +1 @@ -79f82f287ff427b896ab8aedc81f4671949571316b90e89806f0753ed0628dd5ac32f8d0e4a3f5f5f5c65a9e8a06a6869eb5fc896b4bb423870877afca2a7d4f \ No newline at end of file +6e02d01fc115c54e4c60477a9f963149bc83c82e6a7bdb2d29125b131546150c59de5ac3aae5b9469b89b343ecf2574f6c949b918ccb43231e8666ee2b7c1b7c From 8b0b11b6c5a880ece27127e5fd1930e539dcafad Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Wed, 22 Dec 2021 18:37:55 +0100 Subject: [PATCH 06/17] Bug: Forget to declare and define UnstructuredMesh::get_mesh_type() --- include/openmc/mesh.h | 3 +- src/mesh.cpp | 8 +- .../cpp_driver/CMakeLists.txt | 6 + .../.nfs000000002b2100e400000105 | 356 ++++++++++++++++++ .../source_dlopen/CMakeLists.txt | 6 + .../CMakeLists.txt | 6 + 6 files changed, 382 insertions(+), 3 deletions(-) create mode 100644 tests/regression_tests/cpp_driver/CMakeLists.txt create mode 100644 tests/regression_tests/mgxs_library_mesh/.nfs000000002b2100e400000105 create mode 100644 tests/regression_tests/source_dlopen/CMakeLists.txt create mode 100644 tests/regression_tests/source_parameterized_dlopen/CMakeLists.txt diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index f6908be432..990479ea8a 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -437,7 +437,8 @@ public: UnstructuredMesh(const std::string& filename); static const std::string mesh_type; - + virtual std::string get_mesh_type() const override; + // Overridden Methods void surface_bins_crossed(Position r0, Position r1, const Direction& u, vector& bins) const override; diff --git a/src/mesh.cpp b/src/mesh.cpp index bd139e94be..9d4c14b185 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -201,6 +201,10 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node) const std::string UnstructuredMesh::mesh_type = "unstructured"; +virtual std::string UnstructuredMesh::get_mesh_type() const { + return mesh_type; +} + void UnstructuredMesh::surface_bins_crossed( Position r0, Position r1, const Direction& u, vector& bins) const { @@ -1956,7 +1960,7 @@ double MOABMesh::volume(int bin) const std::string MOABMesh::library() const { - return MOABMesh::mesh_lib_type; + return mesh_lib_type; } double MOABMesh::tet_volume(moab::EntityHandle tet) const @@ -2351,7 +2355,7 @@ Position LibMesh::centroid(int bin) const std::string LibMesh::library() const { - return LibMesh::mesh_lib_type; + return mesh_lib_type; } int LibMesh::n_bins() const diff --git a/tests/regression_tests/cpp_driver/CMakeLists.txt b/tests/regression_tests/cpp_driver/CMakeLists.txt new file mode 100644 index 0000000000..e6aa3181b1 --- /dev/null +++ b/tests/regression_tests/cpp_driver/CMakeLists.txt @@ -0,0 +1,6 @@ + +cmake_minimum_required(VERSION 3.3 FATAL_ERROR) +project(openmc_cpp_driver CXX) +add_executable(cpp_driver driver.cpp) +find_package(OpenMC REQUIRED HINTS /home/intern/schumann/src/openmc/openmc-ojs/build) +target_link_libraries(cpp_driver OpenMC::libopenmc) diff --git a/tests/regression_tests/mgxs_library_mesh/.nfs000000002b2100e400000105 b/tests/regression_tests/mgxs_library_mesh/.nfs000000002b2100e400000105 new file mode 100644 index 0000000000..0acc2f407c --- /dev/null +++ b/tests/regression_tests/mgxs_library_mesh/.nfs000000002b2100e400000105 @@ -0,0 +1,356 @@ + mesh 1 group in nuclide mean std. dev. + x y z +0 1 1 1 1 total 0.102319 0.005483 +2 1 2 1 1 total 0.104659 0.002878 +1 2 1 1 1 total 0.107122 0.005105 +3 2 2 1 1 total 0.103856 0.003459 + mesh 1 group in nuclide mean std. dev. + x y z +0 1 1 1 1 total 0.072455 0.005799 +2 1 2 1 1 total 0.072899 0.003254 +1 2 1 1 1 total 0.074187 0.005441 +3 2 2 1 1 total 0.074241 0.003779 + mesh 1 group in nuclide mean std. dev. + x y z +0 1 1 1 1 total 0.072455 0.005799 +2 1 2 1 1 total 0.072912 0.003251 +1 2 1 1 1 total 0.074140 0.005442 +3 2 2 1 1 total 0.074192 0.003783 + mesh 1 group in nuclide mean std. dev. + x y z +0 1 1 1 1 total 0.012943 0.000836 +2 1 2 1 1 total 0.013363 0.000553 +1 2 1 1 1 total 0.013980 0.000715 +3 2 2 1 1 total 0.013286 0.000621 + mesh 1 group in nuclide mean std. dev. + x y z +0 1 1 1 1 total 0.001231 0.000979 +2 1 2 1 1 total 0.001332 0.000691 +1 2 1 1 1 total 0.001346 0.000770 +3 2 2 1 1 total 0.001258 0.000815 + mesh 1 group in nuclide mean std. dev. + x y z +0 1 1 1 1 total 0.011712 0.000762 +2 1 2 1 1 total 0.012031 0.000468 +1 2 1 1 1 total 0.012635 0.000646 +3 2 2 1 1 total 0.012028 0.000556 + mesh 1 group in nuclide mean std. dev. + x y z +0 1 1 1 1 total 0.030549 0.001995 +2 1 2 1 1 total 0.031338 0.001160 +1 2 1 1 1 total 0.032944 0.001703 +3 2 2 1 1 total 0.031480 0.001423 + mesh 1 group in nuclide mean std. dev. + x y z +0 1 1 1 1 total 2.265259e+06 147469.851100 +2 1 2 1 1 total 2.326873e+06 90605.035909 +1 2 1 1 1 total 2.443628e+06 124917.160682 +3 2 2 1 1 total 2.326322e+06 107513.428248 + mesh 1 group in nuclide mean std. dev. + x y z +0 1 1 1 1 total 0.089376 0.004670 +2 1 2 1 1 total 0.091296 0.002370 +1 2 1 1 1 total 0.093142 0.004401 +3 2 2 1 1 total 0.090570 0.002842 + mesh 1 group in nuclide mean std. dev. + x y z +0 1 1 1 1 total 0.089670 0.005421 +2 1 2 1 1 total 0.094557 0.003843 +1 2 1 1 1 total 0.094972 0.005963 +3 2 2 1 1 total 0.088538 0.002530 + mesh 1 group in group out legendre nuclide mean std. dev. + x y z +0 1 1 1 1 1 P0 total 0.089670 0.005421 +1 1 1 1 1 1 P1 total 0.029864 0.001887 +2 1 1 1 1 1 P2 total 0.015945 0.001296 +3 1 1 1 1 1 P3 total 0.009511 0.000946 +8 1 2 1 1 1 P0 total 0.094524 0.003835 +9 1 2 1 1 1 P1 total 0.031760 0.001520 +10 1 2 1 1 1 P2 total 0.017210 0.000821 +11 1 2 1 1 1 P3 total 0.009315 0.000631 +4 2 1 1 1 1 P0 total 0.094835 0.005973 +5 2 1 1 1 1 P1 total 0.032936 0.001884 +6 2 1 1 1 1 P2 total 0.017196 0.001626 +7 2 1 1 1 1 P3 total 0.010278 0.001345 +12 2 2 1 1 1 P0 total 0.088412 0.002491 +13 2 2 1 1 1 P1 total 0.029615 0.001524 +14 2 2 1 1 1 P2 total 0.016925 0.000613 +15 2 2 1 1 1 P3 total 0.009759 0.000466 + mesh 1 group in group out legendre nuclide mean std. dev. + x y z +0 1 1 1 1 1 P0 total 0.089670 0.005421 +1 1 1 1 1 1 P1 total 0.029864 0.001887 +2 1 1 1 1 1 P2 total 0.015945 0.001296 +3 1 1 1 1 1 P3 total 0.009511 0.000946 +8 1 2 1 1 1 P0 total 0.094557 0.003843 +9 1 2 1 1 1 P1 total 0.031747 0.001512 +10 1 2 1 1 1 P2 total 0.017201 0.000817 +11 1 2 1 1 1 P3 total 0.009329 0.000635 +4 2 1 1 1 1 P0 total 0.094972 0.005963 +5 2 1 1 1 1 P1 total 0.032983 0.001886 +6 2 1 1 1 1 P2 total 0.017162 0.001631 +7 2 1 1 1 1 P3 total 0.010240 0.001336 +12 2 2 1 1 1 P0 total 0.088538 0.002530 +13 2 2 1 1 1 P1 total 0.029663 0.001532 +14 2 2 1 1 1 P2 total 0.016931 0.000627 +15 2 2 1 1 1 P3 total 0.009762 0.000466 + mesh 1 group in group out nuclide mean std. dev. + x y z +0 1 1 1 1 1 total 1.000000 0.056046 +2 1 2 1 1 1 total 1.000346 0.042442 +1 2 1 1 1 1 total 1.001447 0.056345 +3 2 2 1 1 1 total 1.001422 0.026106 + mesh 1 group in group out nuclide mean std. dev. + x y z +0 1 1 1 1 1 total 0.031401 0.002755 +2 1 2 1 1 1 total 0.033268 0.001692 +1 2 1 1 1 1 total 0.033756 0.002609 +3 2 2 1 1 1 total 0.030234 0.001308 + mesh 1 group in group out nuclide mean std. dev. + x y z +0 1 1 1 1 1 total 1.0 0.056046 +2 1 2 1 1 1 total 1.0 0.042360 +1 2 1 1 1 1 total 1.0 0.056488 +3 2 2 1 1 1 total 1.0 0.025630 + mesh 1 group in group out legendre nuclide mean std. dev. + x y z +0 1 1 1 1 1 P0 total 0.089376 0.006848 +1 1 1 1 1 1 P1 total 0.029766 0.002346 +2 1 1 1 1 1 P2 total 0.015893 0.001493 +3 1 1 1 1 1 P3 total 0.009480 0.001043 +8 1 2 1 1 1 P0 total 0.091296 0.004536 +9 1 2 1 1 1 P1 total 0.030675 0.001712 +10 1 2 1 1 1 P2 total 0.016622 0.000925 +11 1 2 1 1 1 P3 total 0.008997 0.000662 +4 2 1 1 1 1 P0 total 0.093142 0.006860 +5 2 1 1 1 1 P1 total 0.032348 0.002224 +6 2 1 1 1 1 P2 total 0.016889 0.001722 +7 2 1 1 1 1 P3 total 0.010094 0.001376 +12 2 2 1 1 1 P0 total 0.090570 0.003669 +13 2 2 1 1 1 P1 total 0.030338 0.001794 +14 2 2 1 1 1 P2 total 0.017338 0.000806 +15 2 2 1 1 1 P3 total 0.009997 0.000559 + mesh 1 group in group out legendre nuclide mean std. dev. + x y z +0 1 1 1 1 1 P0 total 0.089376 0.008485 +1 1 1 1 1 1 P1 total 0.029766 0.002878 +2 1 1 1 1 1 P2 total 0.015893 0.001738 +3 1 1 1 1 1 P3 total 0.009480 0.001171 +8 1 2 1 1 1 P0 total 0.091328 0.005967 +9 1 2 1 1 1 P1 total 0.030686 0.002151 +10 1 2 1 1 1 P2 total 0.016628 0.001164 +11 1 2 1 1 1 P3 total 0.009000 0.000764 +4 2 1 1 1 1 P0 total 0.093277 0.008645 +5 2 1 1 1 1 P1 total 0.032394 0.002878 +6 2 1 1 1 1 P2 total 0.016913 0.001969 +7 2 1 1 1 1 P3 total 0.010109 0.001491 +12 2 2 1 1 1 P0 total 0.090699 0.004369 +13 2 2 1 1 1 P1 total 0.030381 0.001963 +14 2 2 1 1 1 P2 total 0.017362 0.000925 +15 2 2 1 1 1 P3 total 0.010011 0.000618 + mesh 1 group out nuclide mean std. dev. + x y z +0 1 1 1 1 total 1.0 0.105972 +2 1 2 1 1 total 1.0 0.060624 +1 2 1 1 1 total 1.0 0.084855 +3 2 2 1 1 total 1.0 0.053024 + mesh 1 group out nuclide mean std. dev. + x y z +0 1 1 1 1 total 1.0 0.108202 +2 1 2 1 1 total 1.0 0.058908 +1 2 1 1 1 total 1.0 0.085583 +3 2 2 1 1 total 1.0 0.052287 + mesh 1 group in nuclide mean std. dev. + x y z +0 1 1 1 1 total 8.523138e-10 3.481466e-11 +2 1 2 1 1 total 9.044552e-10 2.851430e-11 +1 2 1 1 1 total 8.643697e-10 4.046801e-11 +3 2 2 1 1 total 8.667214e-10 1.753703e-11 + mesh 1 group in nuclide mean std. dev. + x y z +0 1 1 1 1 total 0.030356 0.001982 +2 1 2 1 1 total 0.031140 0.001152 +1 2 1 1 1 total 0.032736 0.001693 +3 2 2 1 1 total 0.031283 0.001414 + mesh 1 group in group out nuclide mean std. dev. + x y z +0 1 1 1 1 1 total 0.031233 0.002783 +2 1 2 1 1 1 total 0.033040 0.001647 +1 2 1 1 1 1 total 0.033581 0.002609 +3 2 2 1 1 1 total 0.030070 0.001287 + mesh 1 group in nuclide mean std. dev. + x y surf +3 1 1 x-max in 1 total 0.1816 0.006129 +2 1 1 x-max out 1 total 0.2718 0.095077 +1 1 1 x-min in 1 total 0.0000 0.000000 +0 1 1 x-min out 1 total 0.0000 0.000000 +7 1 1 y-max in 1 total 0.1790 0.011921 +6 1 1 y-max out 1 total 0.2290 0.034821 +5 1 1 y-min in 1 total 0.0000 0.000000 +4 1 1 y-min out 1 total 0.0000 0.000000 +19 1 2 x-max in 1 total 0.1872 0.012447 +18 1 2 x-max out 1 total 0.1952 0.015948 +17 1 2 x-min in 1 total 0.0000 0.000000 +16 1 2 x-min out 1 total 0.0000 0.000000 +23 1 2 y-max in 1 total 0.0000 0.000000 +22 1 2 y-max out 1 total 0.0000 0.000000 +21 1 2 y-min in 1 total 0.2290 0.034821 +20 1 2 y-min out 1 total 0.1790 0.011921 +11 2 1 x-max in 1 total 0.0000 0.000000 +10 2 1 x-max out 1 total 0.0000 0.000000 +9 2 1 x-min in 1 total 0.2718 0.095077 +8 2 1 x-min out 1 total 0.1816 0.006129 +15 2 1 y-max in 1 total 0.1778 0.009484 +14 2 1 y-max out 1 total 0.2326 0.042782 +13 2 1 y-min in 1 total 0.0000 0.000000 +12 2 1 y-min out 1 total 0.0000 0.000000 +27 2 2 x-max in 1 total 0.0000 0.000000 +26 2 2 x-max out 1 total 0.0260 0.026000 +25 2 2 x-min in 1 total 0.1952 0.015948 +24 2 2 x-min out 1 total 0.1872 0.012447 +31 2 2 y-max in 1 total 0.0000 0.000000 +30 2 2 y-max out 1 total 0.0244 0.024400 +29 2 2 y-min in 1 total 0.2326 0.042782 +28 2 2 y-min out 1 total 0.1778 0.009484 + mesh 1 group in nuclide mean std. dev. + x y z +0 1 1 1 1 total 4.600555 0.368212 +2 1 2 1 1 total 4.572518 0.204129 +1 2 1 1 1 total 4.493168 0.329560 +3 2 2 1 1 total 4.489908 0.228574 + mesh 1 group in nuclide mean std. dev. + x y z +0 1 1 1 1 total 4.600555 0.368212 +2 1 2 1 1 total 4.571696 0.203824 +1 2 1 1 1 total 4.496020 0.330019 +3 2 2 1 1 total 4.492825 0.229077 + mesh 1 delayedgroup group in nuclide mean std. dev. + x y z +0 1 1 1 1 1 total 0.000007 4.371033e-07 +1 1 1 1 2 1 total 0.000035 2.256193e-06 +2 1 1 1 3 1 total 0.000033 2.153958e-06 +3 1 1 1 4 1 total 0.000075 4.829351e-06 +4 1 1 1 5 1 total 0.000031 1.979966e-06 +5 1 1 1 6 1 total 0.000013 8.294008e-07 +12 1 2 1 1 1 total 0.000007 2.763697e-07 +13 1 2 1 2 1 total 0.000036 1.426535e-06 +14 1 2 1 3 1 total 0.000034 1.361895e-06 +15 1 2 1 4 1 total 0.000077 3.053480e-06 +16 1 2 1 5 1 total 0.000031 1.251884e-06 +17 1 2 1 6 1 total 0.000013 5.244097e-07 +6 2 1 1 1 1 total 0.000007 3.644103e-07 +7 2 1 1 2 1 total 0.000038 1.880974e-06 +8 2 1 1 3 1 total 0.000036 1.795741e-06 +9 2 1 1 4 1 total 0.000080 4.026199e-06 +10 2 1 1 5 1 total 0.000033 1.650685e-06 +11 2 1 1 6 1 total 0.000014 6.914661e-07 +18 2 2 1 1 1 total 0.000007 3.263844e-07 +19 2 2 1 2 1 total 0.000036 1.684696e-06 +20 2 2 1 3 1 total 0.000034 1.608357e-06 +21 2 2 1 4 1 total 0.000076 3.606069e-06 +22 2 2 1 5 1 total 0.000031 1.478437e-06 +23 2 2 1 6 1 total 0.000013 6.193123e-07 + mesh 1 delayedgroup group out nuclide mean std. dev. + x y z +0 1 1 1 1 1 total 0.0 0.000000 +1 1 1 1 2 1 total 1.0 1.414214 +2 1 1 1 3 1 total 1.0 1.414214 +3 1 1 1 4 1 total 1.0 0.579241 +4 1 1 1 5 1 total 1.0 1.414214 +5 1 1 1 6 1 total 0.0 0.000000 +12 1 2 1 1 1 total 1.0 1.414214 +13 1 2 1 2 1 total 0.0 0.000000 +14 1 2 1 3 1 total 1.0 0.866166 +15 1 2 1 4 1 total 1.0 0.868547 +16 1 2 1 5 1 total 1.0 0.873899 +17 1 2 1 6 1 total 1.0 1.414214 +6 2 1 1 1 1 total 0.0 0.000000 +7 2 1 1 2 1 total 1.0 0.654642 +8 2 1 1 3 1 total 1.0 1.414214 +9 2 1 1 4 1 total 1.0 1.414214 +10 2 1 1 5 1 total 0.0 0.000000 +11 2 1 1 6 1 total 0.0 0.000000 +18 2 2 1 1 1 total 1.0 0.867501 +19 2 2 1 2 1 total 0.0 0.000000 +20 2 2 1 3 1 total 0.0 0.000000 +21 2 2 1 4 1 total 1.0 0.867501 +22 2 2 1 5 1 total 1.0 1.414214 +23 2 2 1 6 1 total 1.0 1.414214 + mesh 1 delayedgroup group in nuclide mean std. dev. + x y z +0 1 1 1 1 1 total 0.000221 0.000018 +1 1 1 1 2 1 total 0.001141 0.000091 +2 1 1 1 3 1 total 0.001090 0.000087 +3 1 1 1 4 1 total 0.002443 0.000194 +4 1 1 1 5 1 total 0.001002 0.000080 +5 1 1 1 6 1 total 0.000420 0.000033 +12 1 2 1 1 1 total 0.000221 0.000011 +13 1 2 1 2 1 total 0.001142 0.000059 +14 1 2 1 3 1 total 0.001090 0.000056 +15 1 2 1 4 1 total 0.002444 0.000126 +16 1 2 1 5 1 total 0.001002 0.000052 +17 1 2 1 6 1 total 0.000420 0.000022 +6 2 1 1 1 1 total 0.000221 0.000013 +7 2 1 1 2 1 total 0.001140 0.000065 +8 2 1 1 3 1 total 0.001088 0.000062 +9 2 1 1 4 1 total 0.002440 0.000140 +10 2 1 1 5 1 total 0.001000 0.000057 +11 2 1 1 6 1 total 0.000419 0.000024 +18 2 2 1 1 1 total 0.000219 0.000014 +19 2 2 1 2 1 total 0.001132 0.000072 +20 2 2 1 3 1 total 0.001081 0.000069 +21 2 2 1 4 1 total 0.002424 0.000155 +22 2 2 1 5 1 total 0.000994 0.000064 +23 2 2 1 6 1 total 0.000416 0.000027 + mesh 1 delayedgroup nuclide mean std. dev. + x y z +0 1 1 1 1 total 0.013336 0.001054 +1 1 1 1 2 total 0.032739 0.002588 +2 1 1 1 3 total 0.120780 0.009548 +3 1 1 1 4 total 0.302780 0.023936 +4 1 1 1 5 total 0.849490 0.067157 +5 1 1 1 6 total 2.853000 0.225544 +12 1 2 1 1 total 0.013336 0.000716 +13 1 2 1 2 total 0.032739 0.001758 +14 1 2 1 3 total 0.120780 0.006485 +15 1 2 1 4 total 0.302780 0.016257 +16 1 2 1 5 total 0.849490 0.045611 +17 1 2 1 6 total 2.853000 0.153186 +6 2 1 1 1 total 0.013336 0.000744 +7 2 1 1 2 total 0.032739 0.001827 +8 2 1 1 3 total 0.120780 0.006740 +9 2 1 1 4 total 0.302780 0.016897 +10 2 1 1 5 total 0.849490 0.047407 +11 2 1 1 6 total 2.853000 0.159216 +18 2 2 1 1 total 0.013336 0.000872 +19 2 2 1 2 total 0.032739 0.002141 +20 2 2 1 3 total 0.120780 0.007899 +21 2 2 1 4 total 0.302780 0.019803 +22 2 2 1 5 total 0.849490 0.055560 +23 2 2 1 6 total 2.853000 0.186598 + mesh 1 delayedgroup group in group out nuclide mean std. dev. + x y z +0 1 1 1 1 1 1 total 0.000000 0.000000 +1 1 1 1 2 1 1 total 0.000029 0.000030 +2 1 1 1 3 1 1 total 0.000028 0.000028 +3 1 1 1 4 1 1 total 0.000083 0.000034 +4 1 1 1 5 1 1 total 0.000028 0.000028 +5 1 1 1 6 1 1 total 0.000000 0.000000 +12 1 2 1 1 1 1 total 0.000029 0.000029 +13 1 2 1 2 1 1 total 0.000000 0.000000 +14 1 2 1 3 1 1 total 0.000052 0.000032 +15 1 2 1 4 1 1 total 0.000055 0.000034 +16 1 2 1 5 1 1 total 0.000059 0.000037 +17 1 2 1 6 1 1 total 0.000033 0.000033 +6 2 1 1 1 1 1 total 0.000000 0.000000 +7 2 1 1 2 1 1 total 0.000113 0.000053 +8 2 1 1 3 1 1 total 0.000034 0.000034 +9 2 1 1 4 1 1 total 0.000029 0.000029 +10 2 1 1 5 1 1 total 0.000000 0.000000 +11 2 1 1 6 1 1 total 0.000000 0.000000 +18 2 2 1 1 1 1 total 0.000053 0.000033 +19 2 2 1 2 1 1 total 0.000000 0.000000 +20 2 2 1 3 1 1 total 0.000000 0.000000 +21 2 2 1 4 1 1 total 0.000053 0.000033 +22 2 2 1 5 1 1 total 0.000031 0.000031 +23 2 2 1 6 1 1 total 0.000025 0.000025 diff --git a/tests/regression_tests/source_dlopen/CMakeLists.txt b/tests/regression_tests/source_dlopen/CMakeLists.txt new file mode 100644 index 0000000000..afb9b64176 --- /dev/null +++ b/tests/regression_tests/source_dlopen/CMakeLists.txt @@ -0,0 +1,6 @@ + +cmake_minimum_required(VERSION 3.3 FATAL_ERROR) +project(openmc_sources CXX) +add_library(source SHARED source_sampling.cpp) +find_package(OpenMC REQUIRED HINTS /home/intern/schumann/src/openmc/openmc-ojs/build) +target_link_libraries(source OpenMC::libopenmc) diff --git a/tests/regression_tests/source_parameterized_dlopen/CMakeLists.txt b/tests/regression_tests/source_parameterized_dlopen/CMakeLists.txt new file mode 100644 index 0000000000..6b2e24073d --- /dev/null +++ b/tests/regression_tests/source_parameterized_dlopen/CMakeLists.txt @@ -0,0 +1,6 @@ + +cmake_minimum_required(VERSION 3.3 FATAL_ERROR) +project(openmc_sources CXX) +add_library(source SHARED parameterized_source_sampling.cpp) +find_package(OpenMC REQUIRED HINTS /home/intern/schumann/src/openmc/openmc-ojs/build) +target_link_libraries(source OpenMC::libopenmc) From c44249815f2f2702d19ef5ea7f160ac0255a95c5 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Wed, 22 Dec 2021 18:58:09 +0100 Subject: [PATCH 07/17] typo, foget virtual --- src/mesh.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 9d4c14b185..01483e549c 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -201,7 +201,7 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node) const std::string UnstructuredMesh::mesh_type = "unstructured"; -virtual std::string UnstructuredMesh::get_mesh_type() const { +std::string UnstructuredMesh::get_mesh_type() const { return mesh_type; } From 5bdc633fba1e0d9fcd03fe704ef3a5fc89f78064 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Wed, 22 Dec 2021 19:56:13 +0100 Subject: [PATCH 08/17] CR at end of line removed --- tests/regression_tests/filter_mesh/results_true.dat | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/regression_tests/filter_mesh/results_true.dat b/tests/regression_tests/filter_mesh/results_true.dat index 4247f6ed1c..97735130e1 100644 --- a/tests/regression_tests/filter_mesh/results_true.dat +++ b/tests/regression_tests/filter_mesh/results_true.dat @@ -1 +1 @@ -6e02d01fc115c54e4c60477a9f963149bc83c82e6a7bdb2d29125b131546150c59de5ac3aae5b9469b89b343ecf2574f6c949b918ccb43231e8666ee2b7c1b7c +6e02d01fc115c54e4c60477a9f963149bc83c82e6a7bdb2d29125b131546150c59de5ac3aae5b9469b89b343ecf2574f6c949b918ccb43231e8666ee2b7c1b7c \ No newline at end of file From b4c19a7950d0faa5fa309eb45a37b82263d76c1e Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Mon, 27 Dec 2021 11:34:43 +0100 Subject: [PATCH 09/17] Update to calc_mesh_volumes Update of calc_mesh_volumes method on CylindricalMesh and SphericalMesh class * use np.diff * Bug in CylindricalMesh for calculation of delta-phi (copy-paste error from SpherialMesh) --- openmc/mesh.py | 18 ++++++------------ 1 file changed, 6 insertions(+), 12 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 0d4b0bab91..d5753be5c8 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -784,12 +784,9 @@ class CylindricalMesh(MeshBase): # 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] + V_r = np.diff(np.array(self.r_grid)**2 / 2) + V_p = np.diff(np.array(self.phi_grid) * np.pi / 180.0) + V_z = np.diff(np.array(self.z_grid)) return np.multiply.outer(np.outer(V_r, V_p), V_z) @@ -950,12 +947,9 @@ class SphericalMesh(MeshBase): # 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] + V_r = np.diff(np.array(self.r_grid)**3 / 3) + V_t = np.diff(-np.cos(np.pi * np.array(self.theta_grid) / 180.0)) + V_p = np.diff(np.array(self.phi_grid) * np.pi / 180) return np.multiply.outer(np.outer(V_r, V_t), V_p) From d1a1cae2786266e5e8b71618bf1693acc6a882c9 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Tue, 4 Jan 2022 13:38:59 +0100 Subject: [PATCH 10/17] Update include/openmc/mesh.h Co-authored-by: Patrick Shriwise --- include/openmc/mesh.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 990479ea8a..48d84d288b 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -164,7 +164,7 @@ public: 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 + //! Determine which cell or surface bins were crossed by a particle // //! \param[in] r0 Previous position of the particle //! \param[in] r1 Current position of the particle From 5384734f627f19efdd04bca60b90eaf8cbae8c4b Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Tue, 4 Jan 2022 15:17:11 +0100 Subject: [PATCH 11/17] Apply suggestions from code review added suggestions from pshriwise Co-authored-by: Patrick Shriwise --- include/openmc/mesh.h | 4 ++-- openmc/mesh.py | 8 ++------ src/mesh.cpp | 25 +++++++++++-------------- 3 files changed, 15 insertions(+), 22 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 48d84d288b..74b899db2a 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -146,7 +146,7 @@ public: { } int nextIndex { -1 }; double distance { INFTY }; - bool maxSurface { true }; + bool max_surface { true }; bool operator<(const MeshDistance& o) const { return distance < o.distance; } @@ -355,7 +355,7 @@ protected: 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 }; + bool full_phi_ { false }; constexpr inline int sanitize_angular_index(int idx, bool full, int N) const { if ((idx > 0) and (idx <= N)) { diff --git a/openmc/mesh.py b/openmc/mesh.py index d5753be5c8..11ad8989e1 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -691,7 +691,7 @@ class CylindricalMesh(MeshBase): @property def indices(self): - nr = len(self.r_grid) - 1 + nr, np, nz = self.dimension np = len(self.phi_grid) - 1 nz = len(self.z_grid) - 1 return ((r, p, z) @@ -782,8 +782,6 @@ class CylindricalMesh(MeshBase): """ - # V = int_r int_phi int_phi r dr dphi dz - V_r = np.diff(np.array(self.r_grid)**2 / 2) V_p = np.diff(np.array(self.phi_grid) * np.pi / 180.0) V_z = np.diff(np.array(self.z_grid)) @@ -854,7 +852,7 @@ class SphericalMesh(MeshBase): @property def indices(self): - nr = len(self.r_grid) - 1 + nr, nt, np = self.dimension nt = len(self.theta_grid) - 1 np = len(self.phi_grid) - 1 return ((r, t, p) @@ -945,8 +943,6 @@ class SphericalMesh(MeshBase): """ - # V = int_r int_theta int_phi r^2 dr sin(theta) dtheta dphi - V_r = np.diff(np.array(self.r_grid)**3 / 3) V_t = np.diff(-np.cos(np.pi * np.array(self.theta_grid) / 180.0)) V_p = np.diff(np.array(self.phi_grid) * np.pi / 180) diff --git a/src/mesh.cpp b/src/mesh.cpp index 01483e549c..899349e350 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -391,7 +391,7 @@ void StructuredMesh::raytrace_mesh(Position r0, Position r1, const Direction& u, bool in_mesh; // 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); + MeshIndex ijk = get_indices(r0 + TINY_BIT*u, in_mesh); // if track is very short, assume that it is completely inside one cell. @@ -480,7 +480,7 @@ void StructuredMesh::bins_crossed(Position r0, Position r1, const Direction& u, 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 surface(const MeshIndex& ijk, int k, bool max, int inward) const {} void track(const MeshIndex& ijk, double l) const { bins.push_back(mesh->get_bin_from_indices(ijk)); lengths.push_back(l); @@ -949,7 +949,7 @@ double CylindricalMesh::find_r_crossing(const Position& r, const Direction& u, 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) + if (-p + D > l) return -p + D; } @@ -979,7 +979,7 @@ double CylindricalMesh::find_phi_crossing(const Position& r, const Direction& u, 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)) + if ((s > l) and ((c0*(r.x + s*u.x) + s0*(r.y + s*u.y)) > 0.0)) return s; } @@ -1166,7 +1166,7 @@ double SphericalMesh::find_r_crossing(const Position& r, const Direction& u, dou // 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) + if (-p + D > l) return -p + D; } @@ -1193,14 +1193,14 @@ double SphericalMesh::find_theta_crossing(const Position& r, const Direction& u, 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); + 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; + 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; } @@ -1216,11 +1216,11 @@ double SphericalMesh::find_theta_crossing(const Position& r, const Direction& u, D = std::sqrt(D); // the solution -p-D is always smaller as -p+D : Check this one first - double s = - p - D; + 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; + 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; @@ -1361,9 +1361,6 @@ void SphericalMesh::to_hdf5(hid_t group) const close_group(mesh_group); } - - - //============================================================================== // Helper functions for the C API //============================================================================== @@ -1615,7 +1612,7 @@ int openmc_structured_mesh_set_grid_impl(int32_t index, return err; } -//! Get the mesh parameters for rectilinear, cylindrical and spharical meshes +//! Get the mesh parameters for rectilinear, cylindrical and spherical meshes template int openmc_structured_mesh_get_grid_impl(int32_t index, double** grid_x, int* nx, double** grid_y, int* ny, double** grid_z, int* nz) From 32d33f25cbea73f256670dd1038e219d7c3ea824 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Tue, 4 Jan 2022 15:21:36 +0100 Subject: [PATCH 12/17] Remove acidentially included files --- .../cpp_driver/CMakeLists.txt | 6 - .../.nfs000000002b2100e400000105 | 356 ------------------ .../source_dlopen/CMakeLists.txt | 6 - .../CMakeLists.txt | 6 - 4 files changed, 374 deletions(-) delete mode 100644 tests/regression_tests/cpp_driver/CMakeLists.txt delete mode 100644 tests/regression_tests/mgxs_library_mesh/.nfs000000002b2100e400000105 delete mode 100644 tests/regression_tests/source_dlopen/CMakeLists.txt delete mode 100644 tests/regression_tests/source_parameterized_dlopen/CMakeLists.txt diff --git a/tests/regression_tests/cpp_driver/CMakeLists.txt b/tests/regression_tests/cpp_driver/CMakeLists.txt deleted file mode 100644 index e6aa3181b1..0000000000 --- a/tests/regression_tests/cpp_driver/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ - -cmake_minimum_required(VERSION 3.3 FATAL_ERROR) -project(openmc_cpp_driver CXX) -add_executable(cpp_driver driver.cpp) -find_package(OpenMC REQUIRED HINTS /home/intern/schumann/src/openmc/openmc-ojs/build) -target_link_libraries(cpp_driver OpenMC::libopenmc) diff --git a/tests/regression_tests/mgxs_library_mesh/.nfs000000002b2100e400000105 b/tests/regression_tests/mgxs_library_mesh/.nfs000000002b2100e400000105 deleted file mode 100644 index 0acc2f407c..0000000000 --- a/tests/regression_tests/mgxs_library_mesh/.nfs000000002b2100e400000105 +++ /dev/null @@ -1,356 +0,0 @@ - mesh 1 group in nuclide mean std. dev. - x y z -0 1 1 1 1 total 0.102319 0.005483 -2 1 2 1 1 total 0.104659 0.002878 -1 2 1 1 1 total 0.107122 0.005105 -3 2 2 1 1 total 0.103856 0.003459 - mesh 1 group in nuclide mean std. dev. - x y z -0 1 1 1 1 total 0.072455 0.005799 -2 1 2 1 1 total 0.072899 0.003254 -1 2 1 1 1 total 0.074187 0.005441 -3 2 2 1 1 total 0.074241 0.003779 - mesh 1 group in nuclide mean std. dev. - x y z -0 1 1 1 1 total 0.072455 0.005799 -2 1 2 1 1 total 0.072912 0.003251 -1 2 1 1 1 total 0.074140 0.005442 -3 2 2 1 1 total 0.074192 0.003783 - mesh 1 group in nuclide mean std. dev. - x y z -0 1 1 1 1 total 0.012943 0.000836 -2 1 2 1 1 total 0.013363 0.000553 -1 2 1 1 1 total 0.013980 0.000715 -3 2 2 1 1 total 0.013286 0.000621 - mesh 1 group in nuclide mean std. dev. - x y z -0 1 1 1 1 total 0.001231 0.000979 -2 1 2 1 1 total 0.001332 0.000691 -1 2 1 1 1 total 0.001346 0.000770 -3 2 2 1 1 total 0.001258 0.000815 - mesh 1 group in nuclide mean std. dev. - x y z -0 1 1 1 1 total 0.011712 0.000762 -2 1 2 1 1 total 0.012031 0.000468 -1 2 1 1 1 total 0.012635 0.000646 -3 2 2 1 1 total 0.012028 0.000556 - mesh 1 group in nuclide mean std. dev. - x y z -0 1 1 1 1 total 0.030549 0.001995 -2 1 2 1 1 total 0.031338 0.001160 -1 2 1 1 1 total 0.032944 0.001703 -3 2 2 1 1 total 0.031480 0.001423 - mesh 1 group in nuclide mean std. dev. - x y z -0 1 1 1 1 total 2.265259e+06 147469.851100 -2 1 2 1 1 total 2.326873e+06 90605.035909 -1 2 1 1 1 total 2.443628e+06 124917.160682 -3 2 2 1 1 total 2.326322e+06 107513.428248 - mesh 1 group in nuclide mean std. dev. - x y z -0 1 1 1 1 total 0.089376 0.004670 -2 1 2 1 1 total 0.091296 0.002370 -1 2 1 1 1 total 0.093142 0.004401 -3 2 2 1 1 total 0.090570 0.002842 - mesh 1 group in nuclide mean std. dev. - x y z -0 1 1 1 1 total 0.089670 0.005421 -2 1 2 1 1 total 0.094557 0.003843 -1 2 1 1 1 total 0.094972 0.005963 -3 2 2 1 1 total 0.088538 0.002530 - mesh 1 group in group out legendre nuclide mean std. dev. - x y z -0 1 1 1 1 1 P0 total 0.089670 0.005421 -1 1 1 1 1 1 P1 total 0.029864 0.001887 -2 1 1 1 1 1 P2 total 0.015945 0.001296 -3 1 1 1 1 1 P3 total 0.009511 0.000946 -8 1 2 1 1 1 P0 total 0.094524 0.003835 -9 1 2 1 1 1 P1 total 0.031760 0.001520 -10 1 2 1 1 1 P2 total 0.017210 0.000821 -11 1 2 1 1 1 P3 total 0.009315 0.000631 -4 2 1 1 1 1 P0 total 0.094835 0.005973 -5 2 1 1 1 1 P1 total 0.032936 0.001884 -6 2 1 1 1 1 P2 total 0.017196 0.001626 -7 2 1 1 1 1 P3 total 0.010278 0.001345 -12 2 2 1 1 1 P0 total 0.088412 0.002491 -13 2 2 1 1 1 P1 total 0.029615 0.001524 -14 2 2 1 1 1 P2 total 0.016925 0.000613 -15 2 2 1 1 1 P3 total 0.009759 0.000466 - mesh 1 group in group out legendre nuclide mean std. dev. - x y z -0 1 1 1 1 1 P0 total 0.089670 0.005421 -1 1 1 1 1 1 P1 total 0.029864 0.001887 -2 1 1 1 1 1 P2 total 0.015945 0.001296 -3 1 1 1 1 1 P3 total 0.009511 0.000946 -8 1 2 1 1 1 P0 total 0.094557 0.003843 -9 1 2 1 1 1 P1 total 0.031747 0.001512 -10 1 2 1 1 1 P2 total 0.017201 0.000817 -11 1 2 1 1 1 P3 total 0.009329 0.000635 -4 2 1 1 1 1 P0 total 0.094972 0.005963 -5 2 1 1 1 1 P1 total 0.032983 0.001886 -6 2 1 1 1 1 P2 total 0.017162 0.001631 -7 2 1 1 1 1 P3 total 0.010240 0.001336 -12 2 2 1 1 1 P0 total 0.088538 0.002530 -13 2 2 1 1 1 P1 total 0.029663 0.001532 -14 2 2 1 1 1 P2 total 0.016931 0.000627 -15 2 2 1 1 1 P3 total 0.009762 0.000466 - mesh 1 group in group out nuclide mean std. dev. - x y z -0 1 1 1 1 1 total 1.000000 0.056046 -2 1 2 1 1 1 total 1.000346 0.042442 -1 2 1 1 1 1 total 1.001447 0.056345 -3 2 2 1 1 1 total 1.001422 0.026106 - mesh 1 group in group out nuclide mean std. dev. - x y z -0 1 1 1 1 1 total 0.031401 0.002755 -2 1 2 1 1 1 total 0.033268 0.001692 -1 2 1 1 1 1 total 0.033756 0.002609 -3 2 2 1 1 1 total 0.030234 0.001308 - mesh 1 group in group out nuclide mean std. dev. - x y z -0 1 1 1 1 1 total 1.0 0.056046 -2 1 2 1 1 1 total 1.0 0.042360 -1 2 1 1 1 1 total 1.0 0.056488 -3 2 2 1 1 1 total 1.0 0.025630 - mesh 1 group in group out legendre nuclide mean std. dev. - x y z -0 1 1 1 1 1 P0 total 0.089376 0.006848 -1 1 1 1 1 1 P1 total 0.029766 0.002346 -2 1 1 1 1 1 P2 total 0.015893 0.001493 -3 1 1 1 1 1 P3 total 0.009480 0.001043 -8 1 2 1 1 1 P0 total 0.091296 0.004536 -9 1 2 1 1 1 P1 total 0.030675 0.001712 -10 1 2 1 1 1 P2 total 0.016622 0.000925 -11 1 2 1 1 1 P3 total 0.008997 0.000662 -4 2 1 1 1 1 P0 total 0.093142 0.006860 -5 2 1 1 1 1 P1 total 0.032348 0.002224 -6 2 1 1 1 1 P2 total 0.016889 0.001722 -7 2 1 1 1 1 P3 total 0.010094 0.001376 -12 2 2 1 1 1 P0 total 0.090570 0.003669 -13 2 2 1 1 1 P1 total 0.030338 0.001794 -14 2 2 1 1 1 P2 total 0.017338 0.000806 -15 2 2 1 1 1 P3 total 0.009997 0.000559 - mesh 1 group in group out legendre nuclide mean std. dev. - x y z -0 1 1 1 1 1 P0 total 0.089376 0.008485 -1 1 1 1 1 1 P1 total 0.029766 0.002878 -2 1 1 1 1 1 P2 total 0.015893 0.001738 -3 1 1 1 1 1 P3 total 0.009480 0.001171 -8 1 2 1 1 1 P0 total 0.091328 0.005967 -9 1 2 1 1 1 P1 total 0.030686 0.002151 -10 1 2 1 1 1 P2 total 0.016628 0.001164 -11 1 2 1 1 1 P3 total 0.009000 0.000764 -4 2 1 1 1 1 P0 total 0.093277 0.008645 -5 2 1 1 1 1 P1 total 0.032394 0.002878 -6 2 1 1 1 1 P2 total 0.016913 0.001969 -7 2 1 1 1 1 P3 total 0.010109 0.001491 -12 2 2 1 1 1 P0 total 0.090699 0.004369 -13 2 2 1 1 1 P1 total 0.030381 0.001963 -14 2 2 1 1 1 P2 total 0.017362 0.000925 -15 2 2 1 1 1 P3 total 0.010011 0.000618 - mesh 1 group out nuclide mean std. dev. - x y z -0 1 1 1 1 total 1.0 0.105972 -2 1 2 1 1 total 1.0 0.060624 -1 2 1 1 1 total 1.0 0.084855 -3 2 2 1 1 total 1.0 0.053024 - mesh 1 group out nuclide mean std. dev. - x y z -0 1 1 1 1 total 1.0 0.108202 -2 1 2 1 1 total 1.0 0.058908 -1 2 1 1 1 total 1.0 0.085583 -3 2 2 1 1 total 1.0 0.052287 - mesh 1 group in nuclide mean std. dev. - x y z -0 1 1 1 1 total 8.523138e-10 3.481466e-11 -2 1 2 1 1 total 9.044552e-10 2.851430e-11 -1 2 1 1 1 total 8.643697e-10 4.046801e-11 -3 2 2 1 1 total 8.667214e-10 1.753703e-11 - mesh 1 group in nuclide mean std. dev. - x y z -0 1 1 1 1 total 0.030356 0.001982 -2 1 2 1 1 total 0.031140 0.001152 -1 2 1 1 1 total 0.032736 0.001693 -3 2 2 1 1 total 0.031283 0.001414 - mesh 1 group in group out nuclide mean std. dev. - x y z -0 1 1 1 1 1 total 0.031233 0.002783 -2 1 2 1 1 1 total 0.033040 0.001647 -1 2 1 1 1 1 total 0.033581 0.002609 -3 2 2 1 1 1 total 0.030070 0.001287 - mesh 1 group in nuclide mean std. dev. - x y surf -3 1 1 x-max in 1 total 0.1816 0.006129 -2 1 1 x-max out 1 total 0.2718 0.095077 -1 1 1 x-min in 1 total 0.0000 0.000000 -0 1 1 x-min out 1 total 0.0000 0.000000 -7 1 1 y-max in 1 total 0.1790 0.011921 -6 1 1 y-max out 1 total 0.2290 0.034821 -5 1 1 y-min in 1 total 0.0000 0.000000 -4 1 1 y-min out 1 total 0.0000 0.000000 -19 1 2 x-max in 1 total 0.1872 0.012447 -18 1 2 x-max out 1 total 0.1952 0.015948 -17 1 2 x-min in 1 total 0.0000 0.000000 -16 1 2 x-min out 1 total 0.0000 0.000000 -23 1 2 y-max in 1 total 0.0000 0.000000 -22 1 2 y-max out 1 total 0.0000 0.000000 -21 1 2 y-min in 1 total 0.2290 0.034821 -20 1 2 y-min out 1 total 0.1790 0.011921 -11 2 1 x-max in 1 total 0.0000 0.000000 -10 2 1 x-max out 1 total 0.0000 0.000000 -9 2 1 x-min in 1 total 0.2718 0.095077 -8 2 1 x-min out 1 total 0.1816 0.006129 -15 2 1 y-max in 1 total 0.1778 0.009484 -14 2 1 y-max out 1 total 0.2326 0.042782 -13 2 1 y-min in 1 total 0.0000 0.000000 -12 2 1 y-min out 1 total 0.0000 0.000000 -27 2 2 x-max in 1 total 0.0000 0.000000 -26 2 2 x-max out 1 total 0.0260 0.026000 -25 2 2 x-min in 1 total 0.1952 0.015948 -24 2 2 x-min out 1 total 0.1872 0.012447 -31 2 2 y-max in 1 total 0.0000 0.000000 -30 2 2 y-max out 1 total 0.0244 0.024400 -29 2 2 y-min in 1 total 0.2326 0.042782 -28 2 2 y-min out 1 total 0.1778 0.009484 - mesh 1 group in nuclide mean std. dev. - x y z -0 1 1 1 1 total 4.600555 0.368212 -2 1 2 1 1 total 4.572518 0.204129 -1 2 1 1 1 total 4.493168 0.329560 -3 2 2 1 1 total 4.489908 0.228574 - mesh 1 group in nuclide mean std. dev. - x y z -0 1 1 1 1 total 4.600555 0.368212 -2 1 2 1 1 total 4.571696 0.203824 -1 2 1 1 1 total 4.496020 0.330019 -3 2 2 1 1 total 4.492825 0.229077 - mesh 1 delayedgroup group in nuclide mean std. dev. - x y z -0 1 1 1 1 1 total 0.000007 4.371033e-07 -1 1 1 1 2 1 total 0.000035 2.256193e-06 -2 1 1 1 3 1 total 0.000033 2.153958e-06 -3 1 1 1 4 1 total 0.000075 4.829351e-06 -4 1 1 1 5 1 total 0.000031 1.979966e-06 -5 1 1 1 6 1 total 0.000013 8.294008e-07 -12 1 2 1 1 1 total 0.000007 2.763697e-07 -13 1 2 1 2 1 total 0.000036 1.426535e-06 -14 1 2 1 3 1 total 0.000034 1.361895e-06 -15 1 2 1 4 1 total 0.000077 3.053480e-06 -16 1 2 1 5 1 total 0.000031 1.251884e-06 -17 1 2 1 6 1 total 0.000013 5.244097e-07 -6 2 1 1 1 1 total 0.000007 3.644103e-07 -7 2 1 1 2 1 total 0.000038 1.880974e-06 -8 2 1 1 3 1 total 0.000036 1.795741e-06 -9 2 1 1 4 1 total 0.000080 4.026199e-06 -10 2 1 1 5 1 total 0.000033 1.650685e-06 -11 2 1 1 6 1 total 0.000014 6.914661e-07 -18 2 2 1 1 1 total 0.000007 3.263844e-07 -19 2 2 1 2 1 total 0.000036 1.684696e-06 -20 2 2 1 3 1 total 0.000034 1.608357e-06 -21 2 2 1 4 1 total 0.000076 3.606069e-06 -22 2 2 1 5 1 total 0.000031 1.478437e-06 -23 2 2 1 6 1 total 0.000013 6.193123e-07 - mesh 1 delayedgroup group out nuclide mean std. dev. - x y z -0 1 1 1 1 1 total 0.0 0.000000 -1 1 1 1 2 1 total 1.0 1.414214 -2 1 1 1 3 1 total 1.0 1.414214 -3 1 1 1 4 1 total 1.0 0.579241 -4 1 1 1 5 1 total 1.0 1.414214 -5 1 1 1 6 1 total 0.0 0.000000 -12 1 2 1 1 1 total 1.0 1.414214 -13 1 2 1 2 1 total 0.0 0.000000 -14 1 2 1 3 1 total 1.0 0.866166 -15 1 2 1 4 1 total 1.0 0.868547 -16 1 2 1 5 1 total 1.0 0.873899 -17 1 2 1 6 1 total 1.0 1.414214 -6 2 1 1 1 1 total 0.0 0.000000 -7 2 1 1 2 1 total 1.0 0.654642 -8 2 1 1 3 1 total 1.0 1.414214 -9 2 1 1 4 1 total 1.0 1.414214 -10 2 1 1 5 1 total 0.0 0.000000 -11 2 1 1 6 1 total 0.0 0.000000 -18 2 2 1 1 1 total 1.0 0.867501 -19 2 2 1 2 1 total 0.0 0.000000 -20 2 2 1 3 1 total 0.0 0.000000 -21 2 2 1 4 1 total 1.0 0.867501 -22 2 2 1 5 1 total 1.0 1.414214 -23 2 2 1 6 1 total 1.0 1.414214 - mesh 1 delayedgroup group in nuclide mean std. dev. - x y z -0 1 1 1 1 1 total 0.000221 0.000018 -1 1 1 1 2 1 total 0.001141 0.000091 -2 1 1 1 3 1 total 0.001090 0.000087 -3 1 1 1 4 1 total 0.002443 0.000194 -4 1 1 1 5 1 total 0.001002 0.000080 -5 1 1 1 6 1 total 0.000420 0.000033 -12 1 2 1 1 1 total 0.000221 0.000011 -13 1 2 1 2 1 total 0.001142 0.000059 -14 1 2 1 3 1 total 0.001090 0.000056 -15 1 2 1 4 1 total 0.002444 0.000126 -16 1 2 1 5 1 total 0.001002 0.000052 -17 1 2 1 6 1 total 0.000420 0.000022 -6 2 1 1 1 1 total 0.000221 0.000013 -7 2 1 1 2 1 total 0.001140 0.000065 -8 2 1 1 3 1 total 0.001088 0.000062 -9 2 1 1 4 1 total 0.002440 0.000140 -10 2 1 1 5 1 total 0.001000 0.000057 -11 2 1 1 6 1 total 0.000419 0.000024 -18 2 2 1 1 1 total 0.000219 0.000014 -19 2 2 1 2 1 total 0.001132 0.000072 -20 2 2 1 3 1 total 0.001081 0.000069 -21 2 2 1 4 1 total 0.002424 0.000155 -22 2 2 1 5 1 total 0.000994 0.000064 -23 2 2 1 6 1 total 0.000416 0.000027 - mesh 1 delayedgroup nuclide mean std. dev. - x y z -0 1 1 1 1 total 0.013336 0.001054 -1 1 1 1 2 total 0.032739 0.002588 -2 1 1 1 3 total 0.120780 0.009548 -3 1 1 1 4 total 0.302780 0.023936 -4 1 1 1 5 total 0.849490 0.067157 -5 1 1 1 6 total 2.853000 0.225544 -12 1 2 1 1 total 0.013336 0.000716 -13 1 2 1 2 total 0.032739 0.001758 -14 1 2 1 3 total 0.120780 0.006485 -15 1 2 1 4 total 0.302780 0.016257 -16 1 2 1 5 total 0.849490 0.045611 -17 1 2 1 6 total 2.853000 0.153186 -6 2 1 1 1 total 0.013336 0.000744 -7 2 1 1 2 total 0.032739 0.001827 -8 2 1 1 3 total 0.120780 0.006740 -9 2 1 1 4 total 0.302780 0.016897 -10 2 1 1 5 total 0.849490 0.047407 -11 2 1 1 6 total 2.853000 0.159216 -18 2 2 1 1 total 0.013336 0.000872 -19 2 2 1 2 total 0.032739 0.002141 -20 2 2 1 3 total 0.120780 0.007899 -21 2 2 1 4 total 0.302780 0.019803 -22 2 2 1 5 total 0.849490 0.055560 -23 2 2 1 6 total 2.853000 0.186598 - mesh 1 delayedgroup group in group out nuclide mean std. dev. - x y z -0 1 1 1 1 1 1 total 0.000000 0.000000 -1 1 1 1 2 1 1 total 0.000029 0.000030 -2 1 1 1 3 1 1 total 0.000028 0.000028 -3 1 1 1 4 1 1 total 0.000083 0.000034 -4 1 1 1 5 1 1 total 0.000028 0.000028 -5 1 1 1 6 1 1 total 0.000000 0.000000 -12 1 2 1 1 1 1 total 0.000029 0.000029 -13 1 2 1 2 1 1 total 0.000000 0.000000 -14 1 2 1 3 1 1 total 0.000052 0.000032 -15 1 2 1 4 1 1 total 0.000055 0.000034 -16 1 2 1 5 1 1 total 0.000059 0.000037 -17 1 2 1 6 1 1 total 0.000033 0.000033 -6 2 1 1 1 1 1 total 0.000000 0.000000 -7 2 1 1 2 1 1 total 0.000113 0.000053 -8 2 1 1 3 1 1 total 0.000034 0.000034 -9 2 1 1 4 1 1 total 0.000029 0.000029 -10 2 1 1 5 1 1 total 0.000000 0.000000 -11 2 1 1 6 1 1 total 0.000000 0.000000 -18 2 2 1 1 1 1 total 0.000053 0.000033 -19 2 2 1 2 1 1 total 0.000000 0.000000 -20 2 2 1 3 1 1 total 0.000000 0.000000 -21 2 2 1 4 1 1 total 0.000053 0.000033 -22 2 2 1 5 1 1 total 0.000031 0.000031 -23 2 2 1 6 1 1 total 0.000025 0.000025 diff --git a/tests/regression_tests/source_dlopen/CMakeLists.txt b/tests/regression_tests/source_dlopen/CMakeLists.txt deleted file mode 100644 index afb9b64176..0000000000 --- a/tests/regression_tests/source_dlopen/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ - -cmake_minimum_required(VERSION 3.3 FATAL_ERROR) -project(openmc_sources CXX) -add_library(source SHARED source_sampling.cpp) -find_package(OpenMC REQUIRED HINTS /home/intern/schumann/src/openmc/openmc-ojs/build) -target_link_libraries(source OpenMC::libopenmc) diff --git a/tests/regression_tests/source_parameterized_dlopen/CMakeLists.txt b/tests/regression_tests/source_parameterized_dlopen/CMakeLists.txt deleted file mode 100644 index 6b2e24073d..0000000000 --- a/tests/regression_tests/source_parameterized_dlopen/CMakeLists.txt +++ /dev/null @@ -1,6 +0,0 @@ - -cmake_minimum_required(VERSION 3.3 FATAL_ERROR) -project(openmc_sources CXX) -add_library(source SHARED parameterized_source_sampling.cpp) -find_package(OpenMC REQUIRED HINTS /home/intern/schumann/src/openmc/openmc-ojs/build) -target_link_libraries(source OpenMC::libopenmc) From 53c15b9e7b7d3021e44f9fd32ba2ed12dc97655d Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Tue, 4 Jan 2022 16:57:30 +0100 Subject: [PATCH 13/17] More updates after suggestions from @pshriwise --- include/openmc/mesh.h | 18 +++--- openmc/mesh.py | 15 +++-- src/mesh.cpp | 129 ++++++++++++++++++++++-------------------- 3 files changed, 87 insertions(+), 75 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 74b899db2a..65f6334774 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -141,8 +141,8 @@ public: struct MeshDistance { MeshDistance() = default; - MeshDistance(int _index, bool _maxSurface, double _distance): - nextIndex{_index}, maxSurface{_maxSurface}, distance{_distance} + MeshDistance(int _index, bool _max_surface, double _distance): + nextIndex{_index}, max_surface{_max_surface}, distance{_distance} { } int nextIndex { -1 }; double distance { INFTY }; @@ -350,7 +350,7 @@ public: int set_grid(); -protected: +private: double find_r_crossing(const Position& r, const Direction& u, double l, int shell) const; double find_phi_crossing(const Position& r, const Direction& u, double l, int shell) const; StructuredMesh::MeshDistance find_z_crossing(const Position& r, const Direction& u, double l, int shell) const; @@ -368,7 +368,7 @@ protected: } inline int sanitize_phi(int idx) const { - return sanitize_angular_index(idx, full_phi, shape_[1]); + return sanitize_angular_index(idx, full_phi_, shape_[1]); } }; @@ -400,13 +400,13 @@ public: int set_grid(); -protected: +private: double find_r_crossing(const Position& r, const Direction& u, double l, int shell) const; double find_theta_crossing(const Position& r, const Direction& u, double l, int shell) const; double find_phi_crossing(const Position& r, const Direction& u, double l, int shell) const; - bool full_theta { false }; - bool full_phi { false }; + 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)) { @@ -419,10 +419,10 @@ protected: } inline int sanitize_theta(int idx) const { - return sanitize_angular_index(idx, full_theta, shape_[1]); + 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]); + return sanitize_angular_index(idx, full_phi_, shape_[2]); } }; diff --git a/openmc/mesh.py b/openmc/mesh.py index 11ad8989e1..5c57468af2 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -649,9 +649,11 @@ class CylindricalMesh(MeshBase): n_dimension : int Number of mesh dimensions (always 3 for a CylindricalMesh). r_grid : Iterable of float - Mesh boundary points along the r-axis. + Mesh boundary points along the r-axis. + Requirement is r >= 0. phi_grid : Iterable of float - Mesh boundary points along the phi-axis. + Mesh boundary points along the phi-axis. + The default value is [0, 360], i.e. the full phi range. z_grid : Iterable of float Mesh boundary points along the z-axis. indices : Iterable of tuple @@ -808,13 +810,16 @@ class SphericalMesh(MeshBase): dimension : Iterable of int The number of mesh cells in each direction. n_dimension : int - Number of mesh dimensions (always 3 for a RectilinearMesh). + Number of mesh dimensions (always 3 for a SphericalMesh). r_grid : Iterable of float Mesh boundary points along the r-axis. + Requirement is r >= 0. theta_grid : Iterable of float - Mesh boundary points along the theta-axis. + Mesh boundary points along the theta-axis in degrees. + The default value is [0, 180], i.e. the full theta range. phi_grid : Iterable of float - Mesh boundary points along the phi-axis. + Mesh boundary points along the phi-axis in degrees. + The default value is [0, 360], i.e. the full phi range. indices : Iterable of tuple An iterable of mesh indices for each mesh element, e.g. [(1, 1, 1), (2, 1, 1), ...] diff --git a/src/mesh.cpp b/src/mesh.cpp index 899349e350..f35717dd00 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -388,10 +388,14 @@ void StructuredMesh::raytrace_mesh(Position r0, Position r1, const Direction& u, const int n = n_dimension_; + // Flag if position is inside the mesh bool in_mesh; + // Position is r = r0 + u * traveled_distance, start at r0 + double traveled_distance { 0.0 }; + // 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); + MeshIndex ijk = get_indices(r0 + TINY_BIT * u, in_mesh); // if track is very short, assume that it is completely inside one cell. @@ -409,8 +413,6 @@ void StructuredMesh::raytrace_mesh(Position r0, Position r1, const Direction& u, 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) { @@ -421,26 +423,26 @@ void StructuredMesh::raytrace_mesh(Position r0, Position r1, const Direction& u, 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); + tally.track(ijk, (std::min(distances[k].distance, total_distance) - traveled_distance) / total_distance); // update position and leave, if we have reached end position - l = distances[k].distance; - if (l >= total_distance) + traveled_distance = distances[k].distance; + if (traveled_distance >= 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); + 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].nextIndex; - distances[k] = distance_to_grid_boundary(ijk, k, r0, u, l); + distances[k] = distance_to_grid_boundary(ijk, k, r0, u, traveled_distance); // Check if we have left the interior of the mesh - in_mesh = ((ijk[k]>=1) and (ijk[k]<=shape_[k])); + 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); + if (in_mesh) tally.surface(ijk, k, !distances[k].max_surface, true); } else { // not inside mesh @@ -448,24 +450,24 @@ void StructuredMesh::raytrace_mesh(Position r0, Position r1, const Direction& u, // 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; + if ((ijk[k] < 1 || ijk[k] > shape_[k]) and (distances[k].distance > traveled_distance)) { + traveled_distance = distances[k].distance; k_max = k; } } // If r1 is not inside the mesh, exit here - if (l >= total_distance) + if (traveled_distance >= total_distance) return; // Calculate the new cell index and update all distances to next surfaces. - ijk = get_indices(r0+(l+TINY_BIT)*u, in_mesh); + ijk = get_indices(r0 + (traveled_distance + TINY_BIT) * u, in_mesh); for (int k = 0; k < n; ++k) { - distances[k] = distance_to_grid_boundary(ijk, k, r0, u, l); + distances[k] = distance_to_grid_boundary(ijk, k, r0, u, traveled_distance); } // If inside the mesh, Tally inward current - if (in_mesh) tally.surface(ijk, k_max, !distances[k_max].maxSurface, 1); + if (in_mesh) tally.surface(ijk, k_max, !distances[k_max].max_surface, true); } } @@ -477,10 +479,10 @@ void StructuredMesh::bins_crossed(Position r0, Position r1, const Direction& u, // 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): + struct TrackAggregator { + TrackAggregator(const StructuredMesh* _mesh, vector& _bins, vector& _lengths): mesh(_mesh), bins(_bins), lengths(_lengths) {} - void surface(const MeshIndex& ijk, int k, bool max, int inward) const {} + void surface(const MeshIndex& ijk, int k, bool max, bool inward) const {} void track(const MeshIndex& ijk, double l) const { bins.push_back(mesh->get_bin_from_indices(ijk)); lengths.push_back(l); @@ -492,7 +494,7 @@ void StructuredMesh::bins_crossed(Position r0, Position r1, const Direction& u, }; // Perform the mesh raytrace with the helper class. - raytrace_mesh(r0, r1, u, TallyBins(this, bins, lengths)); + raytrace_mesh(r0, r1, u, TrackAggregator(this, bins, lengths)); } void StructuredMesh::surface_bins_crossed( @@ -502,12 +504,13 @@ void StructuredMesh::surface_bins_crossed( // 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): + struct SurfaceAggregator { + SurfaceAggregator(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; + void surface(const MeshIndex& ijk, int k, bool max, bool inward) const { + int i_bin = 4 * mesh->n_dimension_ * mesh->get_bin_from_indices(ijk) + 4 * k; + if (max) i_bin += 2; + if (inward) i_bin += 1; bins.push_back(i_bin); } void track(const MeshIndex& idx, double l) const {} @@ -517,7 +520,7 @@ void StructuredMesh::surface_bins_crossed( }; // Perform the mesh raytrace with the helper class. - raytrace_mesh(r0, r1, u, TallyBins(this, bins)); + raytrace_mesh(r0, r1, u, SurfaceAggregator(this, bins)); } @@ -636,11 +639,11 @@ StructuredMesh::MeshDistance RegularMesh::distance_to_grid_boundary(const MeshIn if (std::fabs(u[i]) < FP_PRECISION) return d; - d.maxSurface = (u[i] > 0); - if (d.maxSurface and (ijk[i] <= shape_[i])) { + d.max_surface = (u[i] > 0); + if (d.max_surface 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)) { + } else if (not d.max_surface and (ijk[i] >= 1)) { d.nextIndex--; d.distance = (negative_grid_boundary(ijk, i) - r0[i]) / u[i]; } @@ -799,11 +802,11 @@ StructuredMesh::MeshDistance RectilinearMesh::distance_to_grid_boundary(const Me if (std::fabs(u[i]) < FP_PRECISION) return d; - d.maxSurface = (u[i] > 0); - if (d.maxSurface and (ijk[i] <= shape_[i])) { + d.max_surface = (u[i] > 0); + if (d.max_surface 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)) { + } else if (not d.max_surface and (ijk[i] > 0)) { d.nextIndex--; d.distance = (negative_grid_boundary(ijk, i) - r0[i]) / u[i]; } @@ -920,7 +923,7 @@ StructuredMesh::MeshIndex CylindricalMesh::get_indices(Position r, bool& in_mesh MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh); - mapped_r[1] = sanitize_phi(mapped_r[1]); + idx[1] = sanitize_phi(idx[1]); return idx; } @@ -940,18 +943,21 @@ double CylindricalMesh::find_r_crossing(const Position& r, const Direction& u, d // 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; + // inverse of dominator to help the compiler to speed things up + const double inv_denominator = 1.0 / denominator; - if (D >= 0.0) { - D = std::sqrt(D); + const double p = (u.x*r.x + u.y*r.y) * inv_denominator; + double D = p*p + (r0*r0 - r.x*r.x - r.y*r.y) * inv_denominator; + + if (D < 0.0) return INFTY; + + 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; - } + // 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; } @@ -959,7 +965,7 @@ double CylindricalMesh::find_r_crossing(const Position& r, const Direction& u, d 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; + if (full_phi_ and (shape_[1]==1)) return INFTY; shell = sanitize_phi(shell); @@ -995,11 +1001,11 @@ StructuredMesh::MeshDistance CylindricalMesh::find_z_crossing(const Position& r, if (std::fabs(u.z) < FP_PRECISION) return d; - d.maxSurface = (u.z > 0.0); - if (d.maxSurface and (shell <= shape_[2])) { + d.max_surface = (u.z > 0.0); + if (d.max_surface and (shell <= shape_[2])) { d.nextIndex += 1; d.distance = (grid_[2][shell] - r.z) / u.z; - } else if (not d.maxSurface and (shell > 0)) { + } else if (not d.max_surface and (shell > 0)) { d.nextIndex -= 1; d.distance = (grid_[2][shell-1] - r.z) / u.z; } @@ -1064,7 +1070,7 @@ int CylindricalMesh::set_grid() return OPENMC_E_INVALID_ARGUMENT; } - full_phi = (grid_[1].front() == 0.0) and (grid_[1].back()==360.0); + full_phi_ = (grid_[1].front() == 0.0) and (grid_[1].back()==360.0); // Transform phi-grid from degrees to radians @@ -1143,8 +1149,8 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices(Position r, bool& in_mesh) 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]); + idx[1] = sanitize_theta(idx[1]); + idx[2] = sanitize_phi(idx[2]); return idx; } @@ -1176,7 +1182,7 @@ double SphericalMesh::find_r_crossing(const Position& r, const Direction& u, dou 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; + if (full_theta_ and (shape_[1]==1)) return INFTY; shell = sanitize_theta(shell); @@ -1199,16 +1205,18 @@ double SphericalMesh::find_theta_crossing(const Position& r, const Direction& u, // 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; - } + if (fabs(b) < FP_PRECISION) return INFTY; + + 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; + + // no crossing is possible return INFTY; } const double p = b / a; - double D = p*p - c / a; + double D = p * p - c / a; if (D < 0.0) return INFTY; @@ -1230,7 +1238,7 @@ double SphericalMesh::find_theta_crossing(const Position& r, const Direction& u, 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; + if (full_phi_ and (shape_[2]==1)) return INFTY; shell = sanitize_phi(shell); @@ -1320,8 +1328,8 @@ int SphericalMesh::set_grid() 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); + 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++) @@ -1419,7 +1427,6 @@ extern "C" int openmc_extend_meshes( std::string mesh_type; for (int i = 0; i < n; ++i) { - //if (std::strcmp(type, RegularMesh::mesh_type.c_str()) == 0) { if (RegularMesh::mesh_type == type) { model::meshes.push_back(make_unique()); } else if (RectilinearMesh::mesh_type == type) { From 674f006aaa6bdee6a69fbe364dcde9488e44daba Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Tue, 4 Jan 2022 18:09:00 +0100 Subject: [PATCH 14/17] changed is_sorted to adjacent_find std::is_sorted with less_equal is a standard violation --- src/mesh.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index f35717dd00..415d018e87 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -825,7 +825,7 @@ int RectilinearMesh::set_grid() "must each have at least 2 points"); return OPENMC_E_INVALID_ARGUMENT; } - if (not std::is_sorted(g.begin(), g.end(), std::less_equal())) { + if (std::adjacent_find(g.begin(), g.end(), std::greater_equal<>()) != g.end()) { set_errmsg("Values in for x-, y-, and z- grids for " "rectilinear meshes must be sorted and unique."); return OPENMC_E_INVALID_ARGUMENT; @@ -1047,7 +1047,7 @@ int CylindricalMesh::set_grid() "must each have at least 2 points"); return OPENMC_E_INVALID_ARGUMENT; } - if (not std::is_sorted(g.begin(), g.end(), std::less_equal())) { + if (std::adjacent_find(g.begin(), g.end(), std::greater_equal<>()) != g.end()) { set_errmsg("Values in for r-, phi-, and z- grids for " "cylindrical meshes must be sorted and unique."); return OPENMC_E_INVALID_ARGUMENT; @@ -1305,7 +1305,7 @@ int SphericalMesh::set_grid() "must each have at least 2 points"); return OPENMC_E_INVALID_ARGUMENT; } - if (not std::is_sorted(g.begin(), g.end(), std::less_equal())) { + if (std::adjacent_find(g.begin(), g.end(), std::greater_equal<>()) != g.end()) { set_errmsg("Values in for r-, theta-, and phi- grids for " "spherical meshes must be sorted and unique."); return OPENMC_E_INVALID_ARGUMENT; From fbc2c68eb7e65e9f71e419b9addb147e087db393 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Tue, 11 Jan 2022 15:13:19 +0100 Subject: [PATCH 15/17] Apply suggestions from code review Co-authored-by: Patrick Shriwise --- src/mesh.cpp | 24 ++++++++++-------------- 1 file changed, 10 insertions(+), 14 deletions(-) diff --git a/src/mesh.cpp b/src/mesh.cpp index 415d018e87..a86994f539 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1015,22 +1015,22 @@ StructuredMesh::MeshDistance CylindricalMesh::find_z_crossing(const Position& r, StructuredMesh::MeshDistance CylindricalMesh::distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const { - if (i==0) { + if (i == 0) { return std::min( - MeshDistance(ijk[i]+1, true, find_r_crossing(r0, u, l, ijk[i] )), + 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) { + } 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), 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] ); + return find_z_crossing(r0, u, l, ijk[i] ); } } @@ -1073,7 +1073,6 @@ int CylindricalMesh::set_grid() 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()}; @@ -1157,12 +1156,10 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices(Position r, bool& in_mesh) double SphericalMesh::find_r_crossing(const Position& r, const Direction& u, double l, int shell) const { - 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; @@ -1185,7 +1182,6 @@ double SphericalMesh::find_theta_crossing(const Position& r, const Direction& u, 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 @@ -1269,24 +1265,24 @@ double SphericalMesh::find_phi_crossing(const Position& r, const Direction& u, d StructuredMesh::MeshDistance SphericalMesh::distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const { - if (i==0) { + if (i == 0) { return std::min( - MeshDistance(ijk[i]+1, true, find_r_crossing(r0, u, l, ijk[i] )), + 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) { + } 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), 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), 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)) ); } From 3b3294662e1a0abfa987c37b4a8a15287d419b47 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Wed, 12 Jan 2022 19:43:48 +0000 Subject: [PATCH 16/17] Addressing comments and suggestions of @paulromano manually run clang-format on mesh.cpp and mesh.h --- docs/source/io_formats/tallies.rst | 4 +- include/openmc/mesh.h | 94 +++--- openmc/mesh.py | 28 +- src/mesh.cpp | 510 ++++++++++++++++------------- src/relaxng/tallies.rnc | 12 +- 5 files changed, 361 insertions(+), 287 deletions(-) diff --git a/docs/source/io_formats/tallies.rst b/docs/source/io_formats/tallies.rst index 50ba7b3d5a..3cde6fc854 100644 --- a/docs/source/io_formats/tallies.rst +++ b/docs/source/io_formats/tallies.rst @@ -345,10 +345,10 @@ attributes/sub-elements: :r_grid: The mesh divisions along the r-axis. (For cylindrical and spherical meshes only.) - :p_grid: + :phi_grid: The mesh divisions along the phi-axis. (For cylindrical and spherical meshes only.) - :t_grid: + :theta_grid: The mesh divisions along the theta-axis. (For spherical mesh only.) :library: diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 65f6334774..ae1e648b38 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -141,13 +141,14 @@ public: struct MeshDistance { MeshDistance() = default; - MeshDistance(int _index, bool _max_surface, double _distance): - nextIndex{_index}, max_surface{_max_surface}, distance{_distance} - { } - int nextIndex { -1 }; - double distance { INFTY }; - bool max_surface { true }; - bool operator<(const MeshDistance& o) const { + MeshDistance(int _index, bool _max_surface, double _distance) + : next_index {_index}, max_surface {_max_surface}, distance {_distance} + {} + int next_index {-1}; + double distance {INFTY}; + bool max_surface {true}; + bool operator<(const MeshDistance& o) const + { return distance < o.distance; } }; @@ -169,9 +170,10 @@ public: //! \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; + //! \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 // @@ -206,19 +208,21 @@ public: //! \param[in] i Direction index virtual int get_index_in_direction(double r, int i) const = 0; - //! Get the closest distance from the coordinate r to the grid surface + //! Get the closest distance from the coordinate r to the grid surface //! in i direction that bounds mesh cell ijk and that is larger than l - //! The coordinate r does not have to be inside the mesh cell ijk. In + //! The coordinate r does not have to be inside the mesh cell ijk. In //! curved coordinates, multiple crossings of the same surface can happen, - //! these are selected by the parameter l + //! these are selected by the parameter l //! //! \param[in] ijk Array of mesh indices //! \param[in] i direction index of grid surface //! \param[in] r0 position, from where to calculate the distance //! \param[in] u direction of flight. actual position is r0 + l * u //! \param[in] l actual chord length - //! \return MeshDistance struct with closest distance, next cell index in i-direction and min/max surface indicator - virtual MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const = 0; + //! \return MeshDistance struct with closest distance, next cell index in + //! i-direction and min/max surface indicator + virtual MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, + const Position& r0, const Direction& u, double l) const = 0; //! Get a label for the mesh bin std::string bin_label(int bin) const override; @@ -232,7 +236,6 @@ public: std::array shape_; //!< Number of mesh elements in each dimension protected: - }; //============================================================================== @@ -252,7 +255,8 @@ public: static const std::string mesh_type; - MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override; + 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; @@ -297,8 +301,9 @@ public: virtual std::string get_mesh_type() const override; static const std::string mesh_type; - - MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override; + + 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; @@ -318,7 +323,6 @@ public: //! \param[in] i Direction index double negative_grid_boundary(const MeshIndex& ijk, int i) const; - array, 3> grid_; int set_grid(); @@ -339,7 +343,8 @@ public: static const std::string mesh_type; - MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override; + 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; @@ -351,13 +356,17 @@ public: int set_grid(); private: - double find_r_crossing(const Position& r, const Direction& u, double l, int shell) const; - double find_phi_crossing(const Position& r, const Direction& u, double l, int shell) const; - StructuredMesh::MeshDistance find_z_crossing(const Position& r, const Direction& u, double l, int shell) const; + 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 }; + bool full_phi_ {false}; - constexpr inline int sanitize_angular_index(int idx, bool full, int N) const { + constexpr inline int sanitize_angular_index(int idx, bool full, int N) const + { if ((idx > 0) and (idx <= N)) { return idx; } else if (full) { @@ -367,13 +376,12 @@ private: } } - inline int sanitize_phi(int idx) const { + inline int sanitize_phi(int idx) const + { return sanitize_angular_index(idx, full_phi_, shape_[1]); } - }; - class SphericalMesh : public StructuredMesh { public: // Constructors @@ -389,7 +397,8 @@ public: static const std::string mesh_type; - MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override; + 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; @@ -401,14 +410,18 @@ public: int set_grid(); private: - double find_r_crossing(const Position& r, const Direction& u, double l, int shell) const; - double find_theta_crossing(const Position& r, const Direction& u, double l, int shell) const; - double find_phi_crossing(const Position& r, const Direction& u, double l, int shell) const; + 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 }; + bool full_theta_ {false}; + bool full_phi_ {false}; - constexpr inline int sanitize_angular_index(int idx, bool full, int N) const { + constexpr inline int sanitize_angular_index(int idx, bool full, int N) const + { if ((idx > 0) and (idx <= N)) { return idx; } else if (full) { @@ -418,13 +431,14 @@ private: } } - inline int sanitize_theta(int idx) const { + inline int sanitize_theta(int idx) const + { return sanitize_angular_index(idx, full_theta_, shape_[1]); } - inline int sanitize_phi(int idx) const { + inline int sanitize_phi(int idx) const + { return sanitize_angular_index(idx, full_phi_, shape_[2]); } - }; // Abstract class for unstructured meshes @@ -438,7 +452,7 @@ public: static const std::string mesh_type; virtual std::string get_mesh_type() const override; - + // Overridden Methods void surface_bins_crossed(Position r0, Position r1, const Direction& u, vector& bins) const override; diff --git a/openmc/mesh.py b/openmc/mesh.py index 5c57468af2..3cd208e89d 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -725,13 +725,13 @@ class CylindricalMesh(MeshBase): 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) + phi_grid_str = str(self._phi_grid) if self._phi_grid is None else len(self._phi_grid) + string += fmt.format('\tN Phi pnts:', '=\t', phi_grid_str) if self._phi_grid is not None: string += fmt.format('\tPhi Min:', '=\t', self._phi_grid[0]) string += fmt.format('\tPhi Max:', '=\t', self._phi_grid[-1]) - 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) + 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]) @@ -744,7 +744,7 @@ class CylindricalMesh(MeshBase): # 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.phi_grid = 180 / np.pi * group['phi_grid'][()] mesh.z_grid = group['z_grid'][()] return mesh @@ -766,7 +766,7 @@ class CylindricalMesh(MeshBase): subelement = ET.SubElement(element, "r_grid") subelement.text = ' '.join(map(str, self.r_grid)) - subelement = ET.SubElement(element, "p_grid") + subelement = ET.SubElement(element, "phi_grid") subelement.text = ' '.join(map(str, self.phi_grid)) subelement = ET.SubElement(element, "z_grid") @@ -889,13 +889,13 @@ class SphericalMesh(MeshBase): 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) + theta_grid_str = str(self._theta_grid) if not self._theta_grid else len(self._theta_grid) + string += fmt.format('\tN Theta pnts:', '=\t', theta_grid_str) if self._theta_grid: string += fmt.format('\tTheta Min:', '=\t', self._theta_grid[0]) string += fmt.format('\tTheta Max:', '=\t', self._theta_grid[-1]) - 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) + phi_grid_str = str(self._phi_grid) if not self._phi_grid else len(self._phi_grid) + string += fmt.format('\tN Phi pnts:', '=\t', phi_grid_str) if self._phi_grid: string += fmt.format('\tPhi Min:', '=\t', self._phi_grid[0]) string += fmt.format('\tPhi Max:', '=\t', self._phi_grid[-1]) @@ -908,8 +908,8 @@ class SphericalMesh(MeshBase): # 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'][()] + mesh.theta_grid = 180 / np.pi * group['theta_grid'][()] + mesh.phi_grid = 180 / np.pi * group['phi_grid'][()] return mesh @@ -930,10 +930,10 @@ class SphericalMesh(MeshBase): subelement = ET.SubElement(element, "r_grid") subelement.text = ' '.join(map(str, self.r_grid)) - subelement = ET.SubElement(element, "t_grid") + subelement = ET.SubElement(element, "theta_grid") subelement.text = ' '.join(map(str, self.theta_grid)) - subelement = ET.SubElement(element, "p_grid") + subelement = ET.SubElement(element, "phi_grid") subelement.text = ' '.join(map(str, self.phi_grid)) return element diff --git a/src/mesh.cpp b/src/mesh.cpp index a86994f539..6654684acf 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -33,9 +33,9 @@ #include "openmc/xml_interface.h" #ifdef LIBMESH +#include "libmesh/mesh_modification.h" #include "libmesh/mesh_tools.h" #include "libmesh/numeric_vector.h" -#include "libmesh/mesh_modification.h" #endif namespace openmc { @@ -156,8 +156,9 @@ std::string StructuredMesh::bin_label(int bin) const } } -xt::xtensor StructuredMesh::get_x_shape() const { - return xt::adapt(shape_, { n_dimension_}); +xt::xtensor StructuredMesh::get_x_shape() const +{ + return xt::adapt(shape_, {n_dimension_}); } //============================================================================== @@ -201,7 +202,8 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node) const std::string UnstructuredMesh::mesh_type = "unstructured"; -std::string UnstructuredMesh::get_mesh_type() const { +std::string UnstructuredMesh::get_mesh_type() const +{ return mesh_type; } @@ -249,7 +251,8 @@ void UnstructuredMesh::set_length_multiplier(double length_multiplier) specified_length_multiplier_ = true; } -StructuredMesh::MeshIndex StructuredMesh::get_indices(Position r, bool& in_mesh) const +StructuredMesh::MeshIndex StructuredMesh::get_indices( + Position r, bool& in_mesh) const { MeshIndex ijk; in_mesh = true; @@ -306,7 +309,8 @@ int StructuredMesh::get_bin(Position r) const int StructuredMesh::n_bins() const { - return std::accumulate(shape_.begin(), shape_.begin() + n_dimension_, 1, std::multiplies<>()); + return std::accumulate( + shape_.begin(), shape_.begin() + n_dimension_, 1, std::multiplies<>()); } int StructuredMesh::n_surface_bins() const @@ -369,17 +373,19 @@ xt::xtensor StructuredMesh::count_sites( return counts; } -// 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 +// 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 { - // 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 + // 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(); @@ -392,11 +398,11 @@ void StructuredMesh::raytrace_mesh(Position r0, Position r1, const Direction& u, bool in_mesh; // Position is r = r0 + u * traveled_distance, start at r0 - double traveled_distance { 0.0 }; + double traveled_distance {0.0}; - // Calculate index of current cell. Offset the position a tiny bit in direction of flight + // 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); - // if track is very short, assume that it is completely inside one cell. // Only the current cell will score and no surfaces @@ -413,21 +419,23 @@ void StructuredMesh::raytrace_mesh(Position r0, Position r1, const Direction& u, distances[k] = distance_to_grid_boundary(ijk, k, r0, u, 0.0); } - // Loop until r = r1 is eventually reached while (true) { if (in_mesh) { // find surface with minimal distance to current position - const auto k = std::min_element(distances.begin(), distances.end()) - distances.begin(); + 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) - traveled_distance) / total_distance); + // Tally track length delta since last step + tally.track(ijk, + (std::min(distances[k].distance, total_distance) - traveled_distance) / + total_distance); // update position and leave, if we have reached end position traveled_distance = distances[k].distance; - if (traveled_distance >= total_distance) + if (traveled_distance >= total_distance) return; // If we have not reached r1, we have hit a surface. Tally outward current @@ -435,55 +443,64 @@ void StructuredMesh::raytrace_mesh(Position r0, Position r1, const Direction& u, // 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, traveled_distance); + ijk[k] = distances[k].next_index; + distances[k] = + distance_to_grid_boundary(ijk, k, r0, u, traveled_distance); - // Check if we have left the interior of the mesh - in_mesh = ((ijk[k] >= 1) and (ijk[k] <= shape_[k])); + // Check if we have left the interior of the mesh + in_mesh = ((ijk[k] >= 1) && (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].max_surface, true); + if (in_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 { 0 }; + // 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 > traveled_distance)) { + if ((ijk[k] < 1 || ijk[k] > shape_[k]) && + (distances[k].distance > traveled_distance)) { traveled_distance = distances[k].distance; k_max = k; } } // If r1 is not inside the mesh, exit here - if (traveled_distance >= total_distance) + if (traveled_distance >= total_distance) return; // Calculate the new cell index and update all distances to next surfaces. ijk = get_indices(r0 + (traveled_distance + TINY_BIT) * u, in_mesh); for (int k = 0; k < n; ++k) { - distances[k] = distance_to_grid_boundary(ijk, k, r0, u, traveled_distance); + distances[k] = + distance_to_grid_boundary(ijk, k, r0, u, traveled_distance); } // If inside the mesh, Tally inward current - if (in_mesh) tally.surface(ijk, k_max, !distances[k_max].max_surface, true); - + if (in_mesh) + tally.surface(ijk, k_max, !distances[k_max].max_surface, true); } } } void StructuredMesh::bins_crossed(Position r0, Position r1, const Direction& u, - vector& bins, vector& lengths) const { + 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. + // stores a pointer to the mesh class and references to bins and lengths + // parameters. Performs the actual tally through the track method. struct TrackAggregator { - TrackAggregator(const StructuredMesh* _mesh, vector& _bins, vector& _lengths): - mesh(_mesh), bins(_bins), lengths(_lengths) {} + TrackAggregator( + const StructuredMesh* _mesh, vector& _bins, vector& _lengths) + : mesh(_mesh), bins(_bins), lengths(_lengths) + {} void surface(const MeshIndex& ijk, int k, bool max, bool inward) const {} - void track(const MeshIndex& ijk, double l) const { + void track(const MeshIndex& ijk, double l) const + { bins.push_back(mesh->get_bin_from_indices(ijk)); lengths.push_back(l); } @@ -505,12 +522,17 @@ void StructuredMesh::surface_bins_crossed( // stores a pointer to the mesh class and a reference to the bins parameter. // Performs the actual tally through the surface method. struct SurfaceAggregator { - SurfaceAggregator(const StructuredMesh* _mesh, vector& _bins): - mesh(_mesh), bins(_bins) {} - void surface(const MeshIndex& ijk, int k, bool max, bool inward) const { - int i_bin = 4 * mesh->n_dimension_ * mesh->get_bin_from_indices(ijk) + 4 * k; - if (max) i_bin += 2; - if (inward) i_bin += 1; + SurfaceAggregator(const StructuredMesh* _mesh, vector& _bins) + : mesh(_mesh), bins(_bins) + {} + void surface(const MeshIndex& ijk, int k, bool max, bool inward) const + { + int i_bin = + 4 * mesh->n_dimension_ * mesh->get_bin_from_indices(ijk) + 4 * k; + if (max) + i_bin += 2; + if (inward) + i_bin += 1; bins.push_back(i_bin); } void track(const MeshIndex& idx, double l) const {} @@ -518,12 +540,11 @@ void StructuredMesh::surface_bins_crossed( const StructuredMesh* mesh; vector& bins; }; - + // Perform the mesh raytrace with the helper class. raytrace_mesh(r0, r1, u, SurfaceAggregator(this, bins)); } - //============================================================================== // RegularMesh implementation //============================================================================== @@ -618,7 +639,8 @@ int RegularMesh::get_index_in_direction(double r, int i) const const std::string RegularMesh::mesh_type = "regular"; -std::string RegularMesh::get_mesh_type() const { +std::string RegularMesh::get_mesh_type() const +{ return mesh_type; } @@ -632,25 +654,26 @@ double RegularMesh::negative_grid_boundary(const MeshIndex& ijk, int i) const return lower_left_[i] + (ijk[i] - 1) * width_[i]; } -StructuredMesh::MeshDistance RegularMesh::distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const +StructuredMesh::MeshDistance RegularMesh::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) + d.next_index = ijk[i]; + if (std::abs(u[i]) < FP_PRECISION) return d; d.max_surface = (u[i] > 0); - if (d.max_surface and (ijk[i] <= shape_[i])) { - d.nextIndex++; + if (d.max_surface && (ijk[i] <= shape_[i])) { + d.next_index++; d.distance = (positive_grid_boundary(ijk, i) - r0[i]) / u[i]; - } else if (not d.max_surface and (ijk[i] >= 1)) { - d.nextIndex--; + } else if (!d.max_surface && (ijk[i] >= 1)) { + d.next_index--; 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 { @@ -760,8 +783,6 @@ xt::xtensor RegularMesh::count_sites( return counts; } - - //============================================================================== // RectilinearMesh implementation //============================================================================== @@ -781,33 +802,38 @@ RectilinearMesh::RectilinearMesh(pugi::xml_node node) : StructuredMesh {node} const std::string RectilinearMesh::mesh_type = "rectilinear"; -std::string RectilinearMesh::get_mesh_type() const { +std::string RectilinearMesh::get_mesh_type() const +{ return mesh_type; } -double RectilinearMesh::positive_grid_boundary(const MeshIndex& 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(const MeshIndex& 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 +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) + d.next_index = ijk[i]; + if (std::abs(u[i]) < FP_PRECISION) return d; d.max_surface = (u[i] > 0); - if (d.max_surface and (ijk[i] <= shape_[i])) { - d.nextIndex++; + if (d.max_surface && (ijk[i] <= shape_[i])) { + d.next_index++; d.distance = (positive_grid_boundary(ijk, i) - r0[i]) / u[i]; - } else if (not d.max_surface and (ijk[i] > 0)) { - d.nextIndex--; + } else if (!d.max_surface && (ijk[i] > 0)) { + d.next_index--; d.distance = (negative_grid_boundary(ijk, i) - r0[i]) / u[i]; } return d; @@ -825,7 +851,8 @@ int RectilinearMesh::set_grid() "must each have at least 2 points"); return OPENMC_E_INVALID_ARGUMENT; } - if (std::adjacent_find(g.begin(), g.end(), std::greater_equal<>()) != g.end()) { + if (std::adjacent_find(g.begin(), g.end(), std::greater_equal<>()) != + g.end()) { set_errmsg("Values in for x-, y-, and z- grids for " "rectilinear meshes must be sorted and unique."); return OPENMC_E_INVALID_ARGUMENT; @@ -894,7 +921,7 @@ CylindricalMesh::CylindricalMesh(pugi::xml_node node) : StructuredMesh {node} n_dimension_ = 3; grid_[0] = get_node_array(node, "r_grid"); - grid_[1] = get_node_array(node, "p_grid"); + grid_[1] = get_node_array(node, "phi_grid"); grid_[2] = get_node_array(node, "z_grid"); if (int err = set_grid()) { @@ -904,11 +931,13 @@ CylindricalMesh::CylindricalMesh(pugi::xml_node node) : StructuredMesh {node} const std::string CylindricalMesh::mesh_type = "cylindrical"; -std::string CylindricalMesh::get_mesh_type() const { +std::string CylindricalMesh::get_mesh_type() const +{ return mesh_type; } -StructuredMesh::MeshIndex CylindricalMesh::get_indices(Position r, bool& in_mesh) const +StructuredMesh::MeshIndex CylindricalMesh::get_indices( + Position r, bool& in_mesh) const { Position mapped_r; @@ -918,9 +947,10 @@ StructuredMesh::MeshIndex CylindricalMesh::get_indices(Position r, bool& in_mesh 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; + if (mapped_r[1] < 0) + mapped_r[1] += 2 * M_PI; } - + MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh); idx[1] = sanitize_phi(idx[1]); @@ -928,31 +958,35 @@ StructuredMesh::MeshIndex CylindricalMesh::get_indices(Position r, bool& in_mesh return idx; } -double CylindricalMesh::find_r_crossing(const Position& r, const Direction& u, double l, int shell) const -{ +double CylindricalMesh::find_r_crossing( + const Position& r, const Direction& u, double l, int shell) const +{ - if ((shell < 0) || (shell >= shape_[0])) return INFTY; + if ((shell < 0) || (shell >= shape_[0])) + return INFTY; - //solve r.x^2 + r.y^2 == r0^2 + // 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; + 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; + if (std::abs(denominator) < FP_PRECISION) + return INFTY; // inverse of dominator to help the compiler to speed things up const double inv_denominator = 1.0 / denominator; - const double p = (u.x*r.x + u.y*r.y) * inv_denominator; - double D = p*p + (r0*r0 - r.x*r.x - r.y*r.y) * inv_denominator; + const double p = (u.x * r.x + u.y * r.y) * inv_denominator; + double D = p * p + (r0 * r0 - r.x * r.x - r.y * r.y) * inv_denominator; + + if (D < 0.0) + return INFTY; - if (D < 0.0) return INFTY; - 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; @@ -962,77 +996,81 @@ double CylindricalMesh::find_r_crossing(const Position& r, const Direction& u, d return INFTY; } -double CylindricalMesh::find_phi_crossing(const Position& r, const Direction& u, double l, int shell) const -{ +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; + if (full_phi_ && (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) + // 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)) + if (std::abs(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) && ((c0 * (r.x + s * u.x) + s0 * (r.y + s * u.y)) > 0.0)) return s; } - - return INFTY; + + return INFTY; } -StructuredMesh::MeshDistance CylindricalMesh::find_z_crossing(const Position& r, const Direction& u, double l, int shell) const -{ +StructuredMesh::MeshDistance CylindricalMesh::find_z_crossing( + const Position& r, const Direction& u, double l, int shell) const +{ MeshDistance d; - d.nextIndex = shell; + d.next_index = shell; // Direction of flight is within xy-plane. Will never intersect z. - if (std::fabs(u.z) < FP_PRECISION) + if (std::abs(u.z) < FP_PRECISION) return d; d.max_surface = (u.z > 0.0); - if (d.max_surface and (shell <= shape_[2])) { - d.nextIndex += 1; + if (d.max_surface && (shell <= shape_[2])) { + d.next_index += 1; d.distance = (grid_[2][shell] - r.z) / u.z; - } else if (not d.max_surface and (shell > 0)) { - d.nextIndex -= 1; - d.distance = (grid_[2][shell-1] - r.z) / u.z; + } else if (!d.max_surface && (shell > 0)) { + d.next_index -= 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 +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)) - ); + 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)) - ); + 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] ); - } - + return find_z_crossing(r0, u, l, ijk[i]); + } } int CylindricalMesh::set_grid() @@ -1047,33 +1085,35 @@ int CylindricalMesh::set_grid() "must each have at least 2 points"); return OPENMC_E_INVALID_ARGUMENT; } - if (std::adjacent_find(g.begin(), g.end(), std::greater_equal<>()) != g.end()) { + if (std::adjacent_find(g.begin(), g.end(), std::greater_equal<>()) != + g.end()) { 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; + 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; + 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."); - + "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); + full_phi_ = (grid_[1].front() == 0.0) && (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; }); + 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()}; @@ -1081,7 +1121,6 @@ int CylindricalMesh::set_grid() 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; @@ -1091,7 +1130,7 @@ 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]}; @@ -1103,7 +1142,7 @@ void CylindricalMesh::to_hdf5(hid_t group) const 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, "phi_grid", grid_[1]); write_dataset(mesh_group, "z_grid", grid_[2]); close_group(mesh_group); @@ -1118,8 +1157,8 @@ SphericalMesh::SphericalMesh(pugi::xml_node node) : StructuredMesh {node} n_dimension_ = 3; grid_[0] = get_node_array(node, "r_grid"); - grid_[1] = get_node_array(node, "t_grid"); - grid_[2] = get_node_array(node, "p_grid"); + grid_[1] = get_node_array(node, "theta_grid"); + grid_[2] = get_node_array(node, "phi_grid"); if (int err = set_grid()) { fatal_error(openmc_err_msg); @@ -1128,11 +1167,13 @@ SphericalMesh::SphericalMesh(pugi::xml_node node) : StructuredMesh {node} const std::string SphericalMesh::mesh_type = "spherical"; -std::string SphericalMesh::get_mesh_type() const { +std::string SphericalMesh::get_mesh_type() const +{ return mesh_type; } -StructuredMesh::MeshIndex SphericalMesh::get_indices(Position r, bool& in_mesh) const +StructuredMesh::MeshIndex SphericalMesh::get_indices( + Position r, bool& in_mesh) const { Position mapped_r; @@ -1143,9 +1184,10 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices(Position r, bool& in_mesh) } 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; + if (mapped_r[2] < 0) + mapped_r[2] += 2 * M_PI; } - + MeshIndex idx = StructuredMesh::get_indices(mapped_r, in_mesh); idx[1] = sanitize_theta(idx[1]); @@ -1154,16 +1196,18 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices(Position r, bool& in_mesh) return idx; } -double SphericalMesh::find_r_crossing(const Position& r, const Direction& u, double l, int shell) const +double SphericalMesh::find_r_crossing( + const Position& r, const Direction& u, double l, int shell) const { - if ((shell < 0) || (shell >= shape_[0])) return INFTY; - + 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); + 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 @@ -1176,21 +1220,23 @@ double SphericalMesh::find_r_crossing(const Position& r, const Direction& u, dou return INFTY; } -double SphericalMesh::find_theta_crossing(const Position& r, const Direction& u, double l, int shell) const +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; - + if (full_theta_ && (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 + // 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 double cos_t = std::cos(grid_[1][shell]); const bool sgn = std::signbit(cos_t); const double cos_t_2 = cos_t * cos_t; @@ -1199,21 +1245,24 @@ double SphericalMesh::find_theta_crossing(const Position& r, const Direction& u, 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) return INFTY; - + if (std::abs(a) < FP_PRECISION) { + // if b vanishes, direction of flight is within theta surface and crossing + // is not possible + if (std::abs(b) < FP_PRECISION) + return INFTY; + 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; - + // Check if solution is in positive direction of flight and has correct sign + if ((s > l) && (std::signbit(r.z + s * u.z) == sgn)) + return s; + // no crossing is possible return INFTY; } const double p = b / a; double D = p * p - c / a; - + if (D < 0.0) return INFTY; @@ -1221,72 +1270,76 @@ double SphericalMesh::find_theta_crossing(const Position& r, const Direction& u, // 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; + // Check if solution is in positive direction of flight and has correct sign + if ((s > l) && (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; + // Check if solution is in positive direction of flight and has correct sign + if ((s > l) && (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 -{ +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; + if (full_phi_ && (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) + // 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)) + if (std::abs(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) && ((c0 * (r.x + s * u.x) + s0 * (r.y + s * u.y)) > 0.0)) return s; } - - return INFTY; + + return INFTY; } - -StructuredMesh::MeshDistance SphericalMesh::distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const +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)) - ); + 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)) - ); + 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)) - ); - } + 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() @@ -1301,35 +1354,37 @@ int SphericalMesh::set_grid() "must each have at least 2 points"); return OPENMC_E_INVALID_ARGUMENT; } - if (std::adjacent_find(g.begin(), g.end(), std::greater_equal<>()) != g.end()) { + if (std::adjacent_find(g.begin(), g.end(), std::greater_equal<>()) != + g.end()) { 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; + 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."); - + "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."); + "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); + full_theta_ = (grid_[1].front() == 0.0) && (grid_[1].back() == 180.0); + full_phi_ = (grid_[2].front() == 0.0) && (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; }); + // 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()}; @@ -1337,7 +1392,6 @@ int SphericalMesh::set_grid() 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; @@ -1347,7 +1401,7 @@ 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]}; @@ -1359,8 +1413,8 @@ void SphericalMesh::to_hdf5(hid_t group) const write_dataset(mesh_group, "type", SphericalMesh::mesh_type); write_dataset(mesh_group, "r_grid", grid_[0]); - write_dataset(mesh_group, "t_grid", grid_[1]); - write_dataset(mesh_group, "p_grid", grid_[2]); + write_dataset(mesh_group, "theta_grid", grid_[1]); + write_dataset(mesh_group, "phi_grid", grid_[2]); close_group(mesh_group); } @@ -1423,7 +1477,7 @@ extern "C" int openmc_extend_meshes( std::string mesh_type; for (int i = 0; i < n; ++i) { - if (RegularMesh::mesh_type == type) { + if (RegularMesh::mesh_type == type) { model::meshes.push_back(make_unique()); } else if (RectilinearMesh::mesh_type == type) { model::meshes.push_back(make_unique()); @@ -1583,12 +1637,11 @@ extern "C" int openmc_regular_mesh_set_params( return 0; } - //! Set the mesh parameters for rectilinear, cylindrical and spharical meshes -template -int openmc_structured_mesh_set_grid_impl(int32_t index, - const double* grid_x, const int nx, const double* grid_y, const int ny, - const double* grid_z, const int nz) +template +int openmc_structured_mesh_set_grid_impl(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; @@ -1612,11 +1665,11 @@ int openmc_structured_mesh_set_grid_impl(int32_t index, } int err = m->set_grid(); - return err; + return err; } //! Get the mesh parameters for rectilinear, cylindrical and spherical meshes -template +template int openmc_structured_mesh_get_grid_impl(int32_t index, double** grid_x, int* nx, double** grid_y, int* ny, double** grid_z, int* nz) { @@ -1639,28 +1692,29 @@ int openmc_structured_mesh_get_grid_impl(int32_t index, double** grid_x, return 0; } - //! Get the rectilinear mesh grid extern "C" int openmc_rectilinear_mesh_get_grid(int32_t index, double** grid_x, int* nx, double** grid_y, int* ny, double** grid_z, int* nz) { - return openmc_structured_mesh_get_grid_impl(index, grid_x, nx, grid_y, ny, grid_z, nz); + return openmc_structured_mesh_get_grid_impl( + index, grid_x, nx, grid_y, ny, grid_z, nz); } - //! Set the rectilienar mesh parameters extern "C" int openmc_rectilinear_mesh_set_grid(int32_t index, const double* grid_x, const int nx, const double* grid_y, const int ny, const double* grid_z, const int nz) { - return openmc_structured_mesh_set_grid_impl(index, grid_x, nx, grid_y, ny, grid_z, nz); + return openmc_structured_mesh_set_grid_impl( + index, grid_x, nx, grid_y, ny, grid_z, nz); } //! 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) { - return openmc_structured_mesh_get_grid_impl(index, grid_x, nx, grid_y, ny, grid_z, nz); + return openmc_structured_mesh_get_grid_impl( + index, grid_x, nx, grid_y, ny, grid_z, nz); } //! Set the cylindrical mesh parameters @@ -1668,7 +1722,8 @@ 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) { - return openmc_structured_mesh_set_grid_impl(index, grid_x, nx, grid_y, ny, grid_z, nz); + return openmc_structured_mesh_set_grid_impl( + index, grid_x, nx, grid_y, ny, grid_z, nz); } //! Get the spherical mesh grid @@ -1676,7 +1731,9 @@ 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) { - return openmc_structured_mesh_get_grid_impl(index, grid_x, nx, grid_y, ny, grid_z, nz);; + return openmc_structured_mesh_get_grid_impl( + index, grid_x, nx, grid_y, ny, grid_z, nz); + ; } //! Set the spherical mesh parameters @@ -1684,7 +1741,8 @@ 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) { - return openmc_structured_mesh_set_grid_impl(index, grid_x, nx, grid_y, ny, grid_z, nz); + return openmc_structured_mesh_set_grid_impl( + index, grid_x, nx, grid_y, ny, grid_z, nz); } #ifdef DAGMC @@ -2541,11 +2599,13 @@ void read_meshes(pugi::xml_node root) } else if (mesh_type == SphericalMesh::mesh_type) { model::meshes.push_back(make_unique(node)); #ifdef DAGMC - } else if (mesh_type == UnstructuredMesh::mesh_type && mesh_lib == MOABMesh::mesh_lib_type) { + } else if (mesh_type == UnstructuredMesh::mesh_type && + mesh_lib == MOABMesh::mesh_lib_type) { model::meshes.push_back(make_unique(node)); #endif #ifdef LIBMESH - } else if (mesh_type == UnstructuredMesh::mesh_type && mesh_lib == LibMesh::mesh_lib_type) { + } else if (mesh_type == UnstructuredMesh::mesh_type && + mesh_lib == LibMesh::mesh_lib_type) { model::meshes.push_back(make_unique(node)); #endif } else if (mesh_type == UnstructuredMesh::mesh_type) { diff --git a/src/relaxng/tallies.rnc b/src/relaxng/tallies.rnc index 583b274c60..4e46c45a17 100644 --- a/src/relaxng/tallies.rnc +++ b/src/relaxng/tallies.rnc @@ -29,8 +29,8 @@ element tallies { attribute type { ( "cylindrical" ) }) & (element r_grid { list { xsd:double+ } } | attribute r_grid { list { xsd:double+ } }) & - (element p_grid { list { xsd:double+ } } | - attribute p_grid { list { xsd:double+ } }) & + (element phi_grid { list { xsd:double+ } } | + attribute phi_grid { list { xsd:double+ } }) & (element z_grid { list { xsd:double+ } } | attribute z_grid { list { xsd:double+ } }) ) | ( @@ -38,10 +38,10 @@ element tallies { attribute type { ( "spherical" ) }) & (element r_grid { list { xsd:double+ } } | attribute r_grid { list { xsd:double+ } }) & - (element t_grid { list { xsd:double+ } } | - attribute p_grid { list { xsd:double+ } }) & - (element p_grid { list { xsd:double+ } } | - attribute z_grid { list { xsd:double+ } }) + (element theta_grid { list { xsd:double+ } } | + attribute theta_grid { list { xsd:double+ } }) & + (element phi_grid { list { xsd:double+ } } | + attribute phi_grid { list { xsd:double+ } }) ) ) }* & From f21c52abb85a6eaa3c181b0c2ecc158d688973a1 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Wed, 12 Jan 2022 21:05:09 +0000 Subject: [PATCH 17/17] update test input --- tests/regression_tests/filter_mesh/inputs_true.dat | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/regression_tests/filter_mesh/inputs_true.dat b/tests/regression_tests/filter_mesh/inputs_true.dat index 38b20b7e33..e188768517 100644 --- a/tests/regression_tests/filter_mesh/inputs_true.dat +++ b/tests/regression_tests/filter_mesh/inputs_true.dat @@ -55,13 +55,13 @@ 0.0 0.4411764705882353 0.8823529411764706 1.3235294117647058 1.7647058823529411 2.2058823529411766 2.6470588235294117 3.0882352941176467 3.5294117647058822 3.9705882352941178 4.411764705882353 4.852941176470588 5.294117647058823 5.735294117647059 6.1764705882352935 6.617647058823529 7.0588235294117645 7.5 - 0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0 180.0 200.0 220.0 240.0 260.0 280.0 300.0 320.0 340.0 360.0 + 0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0 180.0 200.0 220.0 240.0 260.0 280.0 300.0 320.0 340.0 360.0 -7.5 -6.5625 -5.625 -4.6875 -3.75 -2.8125 -1.875 -0.9375 0.0 0.9375 1.875 2.8125 3.75 4.6875 5.625 6.5625 7.5 0.0 0.4411764705882353 0.8823529411764706 1.3235294117647058 1.7647058823529411 2.2058823529411766 2.6470588235294117 3.0882352941176467 3.5294117647058822 3.9705882352941178 4.411764705882353 4.852941176470588 5.294117647058823 5.735294117647059 6.1764705882352935 6.617647058823529 7.0588235294117645 7.5 - 0.0 22.5 45.0 67.5 90.0 112.5 135.0 157.5 180.0 - 0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0 180.0 200.0 220.0 240.0 260.0 280.0 300.0 320.0 340.0 360.0 + 0.0 22.5 45.0 67.5 90.0 112.5 135.0 157.5 180.0 + 0.0 20.0 40.0 60.0 80.0 100.0 120.0 140.0 160.0 180.0 200.0 220.0 240.0 260.0 280.0 300.0 320.0 340.0 360.0 1