diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index c3da79c04..ac30015cf 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -46,8 +46,6 @@ jobs: python-version: ["3.12"] mpi: [n, y] omp: [n, y] - dagmc: [n] - libmesh: [n] event: [n] include: @@ -60,33 +58,26 @@ jobs: - python-version: "3.14t" omp: n mpi: n - - dagmc: y - python-version: "3.12" - mpi: y - omp: y - - libmesh: y - python-version: "3.12" - mpi: y - omp: y - - libmesh: y - python-version: "3.12" + - python-version: "3.12" + omp: n mpi: n - omp: y + umesh_libs: y - event: y python-version: "3.12" omp: y mpi: n name: "Python ${{ matrix.python-version }} (omp=${{ matrix.omp }}, - mpi=${{ matrix.mpi }}, dagmc=${{ matrix.dagmc }}, - libmesh=${{ matrix.libmesh }}, event=${{ matrix.event }}" + mpi=${{ matrix.mpi }}, umesh_libs=${{ matrix.umesh_libs }}, + event=${{ matrix.event }}" env: MPI: ${{ matrix.mpi }} PHDF5: ${{ matrix.mpi }} OMP: ${{ matrix.omp }} - DAGMC: ${{ matrix.dagmc }} + DAGMC: ${{ matrix.umesh_libs }} + XDG: ${{ matrix.umesh_libs }} EVENT: ${{ matrix.event }} - LIBMESH: ${{ matrix.libmesh }} + LIBMESH: ${{ matrix.umesh_libs }} NPY_DISABLE_CPU_FEATURES: "AVX512F AVX512_SKX" OPENBLAS_NUM_THREADS: 1 PYTEST_ADDOPTS: --cov=openmc --cov-report=lcov:coverage-python.lcov @@ -146,6 +137,12 @@ jobs: sudo update-alternatives --set mpirun /usr/bin/mpirun.mpich sudo update-alternatives --set mpi-x86_64-linux-gnu /usr/include/x86_64-linux-gnu/mpich + - name: Optional apt dependencies for XDG + shell: bash + if: ${{ matrix.umesh_libs == 'y' }} + run: | + sudo apt install -y libembree-dev + - name: install shell: bash run: | diff --git a/CMakeLists.txt b/CMakeLists.txt index 61d2cc6c5..233b221b6 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -39,6 +39,7 @@ option(OPENMC_BUILD_TESTS "Build tests" option(OPENMC_ENABLE_PROFILE "Compile with profiling flags" OFF) option(OPENMC_ENABLE_COVERAGE "Compile with coverage analysis flags" OFF) option(OPENMC_USE_DAGMC "Enable support for DAGMC (CAD) geometry" OFF) +option(OPENMC_USE_XDG "Enable support for XDG discretized CAD geometry" OFF) option(OPENMC_USE_LIBMESH "Enable support for libMesh unstructured mesh tallies" OFF) option(OPENMC_USE_MPI "Enable MPI" OFF) option(OPENMC_USE_UWUW "Enable UWUW" OFF) @@ -156,6 +157,15 @@ if(OPENMC_USE_DAGMC) endif() endif() +#=============================================================================== +# XDG Mesh and Geometry Support +#=============================================================================== + +if(OPENMC_USE_XDG) + find_package(XDG REQUIRED PATH_SUFFIXES lib/cmake) + message(STATUS "Found XDG: ${XDG_DIR} (version ${XDG_VERSION})") +endif() + #=============================================================================== # libMesh Unstructured Mesh Support #=============================================================================== @@ -428,6 +438,7 @@ list(APPEND libopenmc_SOURCES src/string_utils.cpp src/summary.cpp src/surface.cpp + src/xdg.cpp src/tallies/derivative.cpp src/tallies/filter.cpp src/tallies/filter_azimuthal.cpp @@ -544,6 +555,11 @@ elseif(OPENMC_USE_UWUW) message(FATAL_ERROR "DAGMC must be enabled when UWUW is enabled.") endif() +if(OPENMC_USE_XDG) + target_compile_definitions(libopenmc PRIVATE OPENMC_XDG_ENABLED) + target_link_libraries(libopenmc xdg::xdg) +endif() + if(OPENMC_USE_LIBMESH) target_compile_definitions(libopenmc PRIVATE OPENMC_LIBMESH_ENABLED) target_link_libraries(libopenmc PkgConfig::LIBMESH) diff --git a/cmake/OpenMCConfig.cmake.in b/cmake/OpenMCConfig.cmake.in index 0e15060a3..44b2c30a7 100644 --- a/cmake/OpenMCConfig.cmake.in +++ b/cmake/OpenMCConfig.cmake.in @@ -26,6 +26,10 @@ if("@PNG_FOUND@") find_dependency(PNG) endif() +if(@OPENMC_USE_XDG@) + find_dependency(XDG REQUIRED HINTS @XDG_DIR@) +endif() + if(@OPENMC_USE_MPI@) find_dependency(MPI REQUIRED) endif() diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 0d8189caa..c3b84fd01 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -771,6 +771,9 @@ public: //! Get the library used for this unstructured mesh virtual std::string library() const = 0; + //! Get the mesh filename + virtual const std::string& filename() const { return filename_; } + // Data members bool output_ { true}; //!< Write tallies onto the unstructured mesh at the end of a run @@ -797,7 +800,36 @@ protected: //! \param[in] coords Coordinates of the tetrahedron //! \param[in] seed Random number generation seed //! \return Sampled position within the tetrahedron - Position sample_tet(std::array coords, uint64_t* seed) const; + template + Position sample_tet(span coords, uint64_t* seed) const + { + // Uniform distribution + double s = prn(seed); + double t = prn(seed); + double u = prn(seed); + + // From PyNE implementation of moab tet sampling C. Rocchini & P. Cignoni + // (2000) Generating Random Points in a Tetrahedron, Journal of Graphics + // Tools, 5:4, 9-12, DOI: 10.1080/10867651.2000.10487528 + if (s + t > 1) { + s = 1.0 - s; + t = 1.0 - t; + } + if (s + t + u > 1) { + if (t + u > 1) { + double old_t = t; + t = 1.0 - u; + u = 1.0 - s - old_t; + } else if (t + u <= 1) { + double old_s = s; + s = 1.0 - t - u; + u = old_s + t + u - 1; + } + } + V result = s * (coords[1] - coords[0]) + t * (coords[2] - coords[0]) + + u * (coords[3] - coords[0]) + coords[0]; + return {result[0], result[1], result[2]}; + } // Data members double length_multiplier_ { diff --git a/include/openmc/xdg.h b/include/openmc/xdg.h new file mode 100644 index 000000000..c7250bd53 --- /dev/null +++ b/include/openmc/xdg.h @@ -0,0 +1,125 @@ +#ifndef OPENMC_XDG_H +#define OPENMC_XDG_H + +namespace openmc { +extern "C" const bool XDG_ENABLED; +} + +#ifdef OPENMC_XDG_ENABLED + +#include "xdg/xdg.h" + +#include "openmc/mesh.h" +#include "openmc/position.h" +#include "openmc/xml_interface.h" + +namespace openmc { + +class XDGMesh : public UnstructuredMesh { + +public: + //---------------------------------------------------------------------------- + // Constructors + XDGMesh() = default; + XDGMesh(pugi::xml_node node); + XDGMesh(hid_t group); + XDGMesh(const std::string& filename, double length_multiplier = 1.0); + XDGMesh(std::shared_ptr external_xdg); + + static const std::string mesh_type; + + //---------------------------------------------------------------------------- + // Methods + + //! Get the underlying XDG instance + //! + //! \return Shared pointer to the XDG instance + const std::shared_ptr& xdg_instance() const { return xdg_; } + + //! Check whether a bin index is valid + //! + //! \param[in] bin Bin index to check + //! \return True if the bin index is in [0, n_bins()) + bool bin_is_valid(int bin) const { return bin >= 0 && bin < n_bins(); } + + std::string get_mesh_type() const override { return mesh_type; } + + //! Convert a mesh bin index to an XDG MeshID + //! + //! \param[in] bin Bin index to convert + //! \return XDG MeshID corresponding to the bin + xdg::MeshID bin_to_mesh_id(int bin) const; + + //! Convert an XDG MeshID to a mesh bin index + //! + //! \param[in] id XDG MeshID to convert + //! \return Bin index corresponding to the MeshID, or -1 if invalid + int mesh_id_to_bin(xdg::MeshID id) const; + + //! Get the name of the underlying mesh library + //! + //! \return Name of the mesh library being used in XDG + std::string library() const override; + + //---------------------------------------------------------------------------- + // Overridden Methods + + //! Perform any preparation needed to support use in mesh filters + void prepare_for_point_location() override; + + Position sample_element(int32_t bin, uint64_t* seed) const override; + + void bins_crossed(Position r0, Position r1, const Direction& u, + vector& bins, vector& lengths) const override; + + int get_bin(Position r) const override; + + int n_bins() const override; + + int n_surface_bins() const override; + + std::pair, vector> plot( + Position plot_ll, Position plot_ur) const override; + + //! Add a score to the mesh instance + void add_score(const std::string& score) override {}; + + //! Remove all scores from the mesh instance + void remove_scores() override {}; + + //! Set data for a score + void set_score_data(const std::string& score, const vector& values, + const vector& std_dev) override {}; + + //! Write the mesh with any current tally data + void write(const std::string& base_filename) const override; + + Position centroid(int bin) const override; + + int n_vertices() const override; + + Position vertex(int id) const override; + + std::vector connectivity(int id) const override; + + //! Get the volume of a mesh bin + // + //! \param[in] bin Bin to return the volume for + //! \return Volume of the bin + double volume(int bin) const override; + +private: + void initialize() override; + + //---------------------------------------------------------------------------- + // Private data members + std::shared_ptr xdg_; //!< XDG instance + xdg::MeshLibrary mesh_library_ { + xdg::MeshLibrary::LIBMESH}; //!< Mesh library type +}; + +} // namespace openmc + +#endif // OPENMC_XDG_ENABLED + +#endif // OPENMC_XDG_H diff --git a/openmc/__init__.py b/openmc/__init__.py index c204929c8..ce7b53b23 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -34,6 +34,7 @@ from openmc.plotter import * from openmc.search import * from openmc.polynomial import * from openmc.tracks import * +from openmc.xdg import * from .config import * # Import a few names from the model module diff --git a/openmc/filter.py b/openmc/filter.py index 87aeb70c3..54a497abb 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -920,11 +920,13 @@ class MeshFilter(Filter): def mesh(self, mesh): cv.check_type('filter mesh', mesh, openmc.MeshBase) self._mesh = mesh - if isinstance(mesh, openmc.UnstructuredMesh): + if isinstance(mesh, (openmc.UnstructuredMesh, openmc.XDGMesh)): if mesh.has_statepoint_data: self.bins = list(range(len(mesh.volumes))) else: self.bins = [] + elif isinstance(mesh, openmc.XDGMesh): + self.bins = [] else: self.bins = list(mesh.indices) diff --git a/openmc/lib/__init__.py b/openmc/lib/__init__.py index 9b135370f..6074611df 100644 --- a/openmc/lib/__init__.py +++ b/openmc/lib/__init__.py @@ -40,6 +40,9 @@ else: def _dagmc_enabled(): return c_bool.in_dll(_dll, "DAGMC_ENABLED").value +def _xdg_enabled(): + return c_bool.in_dll(_dll, "XDG_ENABLED").value + def _coord_levels(): return c_int.in_dll(_dll, "n_coord_levels").value diff --git a/openmc/mesh.py b/openmc/mesh.py index ff8e14577..3b37e0dad 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -301,6 +301,8 @@ class MeshBase(IDManagerMixin, ABC): return SphericalMesh.from_hdf5(group, mesh_id, mesh_name) elif mesh_type == 'unstructured': return UnstructuredMesh.from_hdf5(group, mesh_id, mesh_name) + elif mesh_type == 'xdg': + return openmc.XDGMesh.from_hdf5(group, mesh_id, mesh_name) else: raise ValueError('Unrecognized mesh type: "' + mesh_type + '"') @@ -346,6 +348,8 @@ class MeshBase(IDManagerMixin, ABC): mesh = CylindricalMesh.from_xml_element(elem) elif mesh_type == 'spherical': mesh = SphericalMesh.from_xml_element(elem) + elif mesh_type == 'xdg': + mesh = openmc.XDGMesh.from_xml_element(elem) elif mesh_type == 'unstructured': mesh = UnstructuredMesh.from_xml_element(elem) else: @@ -3262,18 +3266,17 @@ class UnstructuredMesh(MeshBase): def from_hdf5(cls, group: h5py.Group, mesh_id: int, name: str): filename = group["filename"][()].decode() library = group["library"][()].decode() - if "options" in group.attrs: - options = group.attrs['options'].decode() - else: - options = None - mesh = cls( - filename=filename, - library=library, - mesh_id=mesh_id, - name=name, - options=options, - ) + kwargs = {'filename': filename, + 'library': library, + 'mesh_id': mesh_id, + 'name': name} + + if "options" in group.attrs: + kwargs['options'] = group.attrs['options'].decode() + + mesh = cls(**kwargs) + mesh._has_statepoint_data = True vol_data = group["volumes"][()] mesh.volumes = np.reshape(vol_data, (vol_data.shape[0],)) diff --git a/openmc/xdg.py b/openmc/xdg.py new file mode 100644 index 000000000..60d6c63af --- /dev/null +++ b/openmc/xdg.py @@ -0,0 +1,107 @@ +from collections.abc import Iterable +from numbers import Integral +from functools import wraps + +import h5py +import lxml.etree as ET + +import openmc +import openmc.checkvalue as cv +from ._xml import get_text +from .bounding_box import BoundingBox +from .checkvalue import PathLike +from .utility_funcs import input_path + + +class XDGMesh(openmc.UnstructuredMesh): + """A 3D unstructured mesh supported by the XDG library. Various mesh + backends are supported in XDG, including MOAB and libMesh. For more + information on supported mesh types and libraries, see + https://xdg-org.github.io/xdg/. + + This class is meant to serve as a reference to an unstructured mesh to be + applied in an OpenMC model in one of the following ways: + + 1. As a volumetric mesh for an tally MeshFilter (see + :class:`openmc.MeshFilter`) + 2. As a surface mesh for an XDGUniverse (see :class:`openmc.XDGUniverse`) + 3. As a volumetric mesh for an XDGUniverse (see :class:`openmc.XDGUniverse`) + + + Parameters + ---------- + filename : path-like + Location of the unstructured mesh file. Supported files for 'moab' + library are .h5, .vtk and .e (exodus) depending on the MOAB build + configuration. Supported files for 'libmesh' library are exodus mesh + files .exo. + library : {'moab', 'libmesh'} + Mesh library used for the unstructured mesh tally + 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 + filename : str + Name of the file containing the unstructured mesh + library : {'moab', 'libmesh'} + Mesh library used for the unstructured mesh tally + """ + def __init__(self, filename: PathLike, library: str | None = None, mesh_id: int | None = None, + name: str = ''): + super().__init__(mesh_id=mesh_id, + name=name, + filename=filename, + library=library) + + # set library based on file extension if not provided by the user + if library is None: + if filename.suffix in ('.h5m', '.vtk', '.e'): + library = 'moab' + elif filename.suffix in ('.exo',): + library = 'libmesh' + + if library is None: + raise ValueError(f"Unable to determine mesh library from file extension {filename.suffix}. " + "Please specify the mesh library type explicitly.") + self.library = library + + def to_xml_element(self): + """Return XML representation of the mesh + + Returns + ------- + element : lxml.etree._Element + XML element containing mesh data + + """ + element = super().to_xml_element() + element.set("type", "xdg") + element.set("library", self._library) + + return element + + @classmethod + def from_xml_element(cls, elem: ET.Element): + """Generate unstructured mesh object from XML element + + Parameters + ---------- + elem : lxml.etree._Element + XML element + + Returns + ------- + openmc.UnstructuredMesh + UnstructuredMesh generated from an XML element + """ + mesh_id = int(get_text(elem, 'id')) + filename = get_text(elem, 'filename') + library = get_text(elem, 'library') + return cls(filename, library, mesh_id) diff --git a/src/initialize.cpp b/src/initialize.cpp index 78b414f78..d453777ad 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -42,6 +42,10 @@ #include "libmesh/libmesh.h" #endif +#ifdef OPENMC_XDG_ENABLED +#include "xdg/config.h" +#endif + int openmc_init(int argc, char* argv[], const void* intracomm) { using namespace openmc; @@ -85,6 +89,14 @@ int openmc_init(int argc, char* argv[], const void* intracomm) settings::libmesh_comm = &(settings::libmesh_init->comm()); } +#if defined(OPENMC_XDG_ENABLED) && defined(XDG_ENABLE_LIBMESH) + // Set the external libMesh initialization and communicator for XDG if + // libMesh was initialized externally. If libMesh was initialized internally, + // the XDG config will use the internal initialization and communicator. + xdg::config::external_libmesh_init = settings::libmesh_init.get(); + xdg::config::external_libmesh_comm = settings::libmesh_comm; +#endif + #endif // Start total and initialization timer diff --git a/src/mesh.cpp b/src/mesh.cpp index 89d1ff24b..30ac4cbab 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -56,6 +56,8 @@ #include "moab/FileOptions.hpp" #endif +#include "openmc/xdg.h" + namespace openmc { //============================================================================== @@ -348,6 +350,10 @@ const std::unique_ptr& Mesh::create( mesh_library == MOABMesh::mesh_lib_type) { model::meshes.push_back(make_unique(dataset)); #endif +#ifdef OPENMC_XDG_ENABLED + } else if (mesh_type == XDGMesh::mesh_type) { + model::meshes.push_back(make_unique(dataset)); +#endif #ifdef OPENMC_LIBMESH_ENABLED } else if (mesh_type == UnstructuredMesh::mesh_type && mesh_library == LibMesh::mesh_lib_type) { @@ -799,12 +805,12 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node) n_dimension_ = 3; // check the mesh type - if (check_for_node(node, "type")) { - auto temp = get_node_value(node, "type", true, true); - if (temp != mesh_type) { - fatal_error(fmt::format("Invalid mesh type: {}", temp)); - } - } + // if (check_for_node(node, "type")) { + // auto temp = get_node_value(node, "type", true, true); + // if (temp != mesh_type) { + // fatal_error(fmt::format("Invalid mesh type: {}", temp)); + // } + // } // check if a length unit multiplier was specified if (check_for_node(node, "length_multiplier")) { @@ -895,36 +901,6 @@ void UnstructuredMesh::determine_bounds() upper_right_ = {xmax, ymax, zmax}; } -Position UnstructuredMesh::sample_tet( - std::array coords, uint64_t* seed) const -{ - // Uniform distribution - double s = prn(seed); - double t = prn(seed); - double u = prn(seed); - - // From PyNE implementation of moab tet sampling C. Rocchini & P. Cignoni - // (2000) Generating Random Points in a Tetrahedron, Journal of Graphics - // Tools, 5:4, 9-12, DOI: 10.1080/10867651.2000.10487528 - if (s + t > 1) { - s = 1.0 - s; - t = 1.0 - t; - } - if (s + t + u > 1) { - if (t + u > 1) { - double old_t = t; - t = 1.0 - u; - u = 1.0 - s - old_t; - } else if (t + u <= 1) { - double old_s = s; - s = 1.0 - t - u; - u = old_s + t + u - 1; - } - } - return s * (coords[1] - coords[0]) + t * (coords[2] - coords[0]) + - u * (coords[3] - coords[0]) + coords[0]; -} - const std::string UnstructuredMesh::mesh_type = "unstructured"; std::string UnstructuredMesh::get_mesh_type() const @@ -992,7 +968,6 @@ void UnstructuredMesh::to_hdf5_inner(hid_t mesh_group) const connectivity.slice(i) = -1; } } - // warn users that some elements were skipped if (num_elem_skipped > 0) { warning(fmt::format("The connectivity of {} elements " @@ -3173,7 +3148,6 @@ std::string MOABMesh::library() const // Sample position within a tet for MOAB type tets Position MOABMesh::sample_element(int32_t bin, uint64_t* seed) const { - moab::EntityHandle tet_ent = get_ent_handle_from_bin(bin); // Get vertex coordinates for MOAB tet @@ -3191,12 +3165,8 @@ Position MOABMesh::sample_element(int32_t bin, uint64_t* seed) const fatal_error("Failed to get tet coords"); } - std::array tet_verts; - for (int i = 0; i < 4; i++) { - tet_verts[i] = {p[i][0], p[i][1], p[i][2]}; - } // Samples position within tet using Barycentric stuff - return this->sample_tet(tet_verts, seed); + return this->sample_tet({p, 4}, seed); } double MOABMesh::tet_volume(moab::EntityHandle tet) const @@ -3688,7 +3658,8 @@ Position LibMesh::sample_element(int32_t bin, uint64_t* seed) const tet_verts[i] = {node_ref(0), node_ref(1), node_ref(2)}; } // Samples position within tet using Barycentric coordinates - Position sampled_position = this->sample_tet(tet_verts, seed); + Position sampled_position = + this->sample_tet({tet_verts.begin(), tet_verts.end()}, seed); if (length_multiplier_ > 0.0) { return length_multiplier_ * sampled_position; } else { diff --git a/src/output.cpp b/src/output.cpp index 04a8caa60..6f427fb35 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -317,6 +317,7 @@ void print_build_info() std::string mpi(n); std::string phdf5(n); std::string dagmc(n); + std::string xdg(n); std::string libmesh(n); std::string png(n); std::string profiling(n); @@ -333,6 +334,9 @@ void print_build_info() #ifdef OPENMC_DAGMC_ENABLED dagmc = y; #endif +#ifdef OPENMC_XDG_ENABLED + xdg = y; +#endif #ifdef OPENMC_LIBMESH_ENABLED libmesh = y; #endif @@ -364,6 +368,7 @@ void print_build_info() fmt::print("Parallel HDF5 enabled: {}\n", phdf5); fmt::print("PNG support: {}\n", png); fmt::print("DAGMC support: {}\n", dagmc); + fmt::print("XDG support: {}\n", xdg); fmt::print("libMesh support: {}\n", libmesh); fmt::print("Coverage testing: {}\n", coverage); fmt::print("Profiling flags: {}\n", profiling); diff --git a/src/xdg.cpp b/src/xdg.cpp new file mode 100644 index 000000000..b2919a9b0 --- /dev/null +++ b/src/xdg.cpp @@ -0,0 +1,259 @@ +#include "openmc/xdg.h" + +#include // for remove_if +#include // for shared_ptr +#include // for string +#include // for pair +#include // for vector + +#include + +#ifdef OPENMC_XDG_ENABLED +#include "xdg/xdg.h" +#endif + +#include "openmc/constants.h" +#include "openmc/container_util.h" +#include "openmc/error.h" +#include "openmc/file_utils.h" +#include "openmc/geometry.h" +#include "openmc/geometry_aux.h" +#include "openmc/hdf5_interface.h" +#include "openmc/material.h" +#include "openmc/settings.h" +#include "openmc/string_utils.h" + +namespace openmc { + +#ifdef OPENMC_XDG_ENABLED +const bool XDG_ENABLED = true; +#else +const bool XDG_ENABLED = false; +#endif + +} // namespace openmc + +#ifdef OPENMC_XDG_ENABLED + +namespace openmc { + +//============================================================================== +// XDG Mesh implementation +//============================================================================== + +const std::string XDGMesh::mesh_type = "xdg"; + +XDGMesh::XDGMesh(pugi::xml_node node) : UnstructuredMesh(node) +{ + std::string mesh_lib = get_node_value(node, "library", true, true); + if (mesh_lib == "moab") { + mesh_library_ = xdg::MeshLibrary::MOAB; + } else if (mesh_lib == "libmesh") { + mesh_library_ = xdg::MeshLibrary::LIBMESH; + } else { + fatal_error(fmt::format("Invalid mesh library specified in XML: {}. " + "Valid options are 'moab' and 'libmesh'.", + mesh_lib)); + } + + initialize(); +} + +XDGMesh::XDGMesh(hid_t group) : UnstructuredMesh(group) +{ + std::string mesh_lib; + read_dataset(group, "library", mesh_lib); + if (mesh_lib == "moab") { + mesh_library_ = xdg::MeshLibrary::MOAB; + } else if (mesh_lib == "libmesh") { + mesh_library_ = xdg::MeshLibrary::LIBMESH; + } else { + fatal_error(fmt::format("Invalid mesh library specified in HDF5: {}. " + "Valid options are 'moab' and 'libmesh'.", + mesh_lib)); + } + initialize(); +} + +XDGMesh::XDGMesh(const std::string& filename, double length_multiplier) +{ + filename_ = filename; + set_length_multiplier(length_multiplier); + initialize(); +} + +XDGMesh::XDGMesh(std::shared_ptr external_xdg) +{ + xdg_ = external_xdg; + filename_ = "unknown (external file)"; + initialize(); +} + +void XDGMesh::initialize() +{ + if (xdg_) + return; + + // create XDGMesh instance + xdg_ = xdg::XDG::create(mesh_library_); + + // load XDGMesh file + if (!file_exists(filename_)) { + fatal_error(fmt::format("Mesh file \"{}\" does not exist", filename_)); + } + + xdg_->mesh_manager()->load_file(filename_); + xdg_->mesh_manager()->init(); + xdg_->mesh_manager()->parse_metadata(); + + auto global_bbox = xdg_->mesh_manager()->global_bounding_box(); + Position lower_left = + length_multiplier_ * + Position(global_bbox.min_x, global_bbox.min_y, global_bbox.min_z); + Position upper_right = + length_multiplier_ * + Position(global_bbox.max_x, global_bbox.max_y, global_bbox.max_z); + lower_left_ = {lower_left.x, lower_left.y, lower_left.z}; + upper_right_ = {upper_right.x, upper_right.y, upper_right.z}; +} + +void XDGMesh::prepare_for_point_location() +{ + xdg_->prepare_raytracer(); +} + +Position XDGMesh::sample_element(int32_t bin, uint64_t* seed) const +{ + auto vertices = xdg_->mesh_manager()->element_vertices(bin_to_mesh_id(bin)); + Position sampled_position = this->sample_tet(vertices, seed); + if (length_multiplier_ > 0.0) { + sampled_position *= length_multiplier_; + } + return sampled_position; +} + +void XDGMesh::bins_crossed(Position r0, Position r1, const Direction& u, + vector& bins, vector& lengths) const +{ + xdg::Position p0 {r0.x, r0.y, r0.z}; + xdg::Position p1 {r1.x, r1.y, r1.z}; + double inv_length = 1 / (p1 - p0).length(); + auto track_segments = xdg_->segments(p0, p1); + // remove elements with lengths of zero + track_segments.erase( + std::remove_if(track_segments.begin(), track_segments.end(), + [](const std::pair& p) { return p.second == 0.0; }), + track_segments.end()); + for (const auto& track_segment : track_segments) { + bins.push_back(mesh_id_to_bin(track_segment.first)); + lengths.push_back(track_segment.second * inv_length); + } +} + +int XDGMesh::get_bin(Position r) const +{ + xdg::Position p {r.x, r.y, r.z}; + if (length_multiplier_ > 0.0) { + p /= length_multiplier_; + } + return mesh_id_to_bin(xdg_->find_element(p)); +} + +int XDGMesh::n_bins() const +{ + return xdg_->mesh_manager()->num_volume_elements(); +} + +int XDGMesh::n_surface_bins() const +{ + return 4 * n_bins(); +} + +std::pair, vector> XDGMesh::plot( + Position plot_ll, Position plot_ur) const +{ + fatal_error("Plot of XDGMesh mesh not implemented"); + return {}; +} + +std::string XDGMesh::library() const +{ + if (mesh_library_ == xdg::MeshLibrary::LIBMESH) { + return "libmesh"; + } else if (mesh_library_ == xdg::MeshLibrary::MOAB) { + return "moab"; + } +} + +void XDGMesh::write(const std::string& base_filename) const +{ + warning("XDGMesh mesh write from C++ not implemented"); +} + +Position XDGMesh::centroid(int bin) const +{ + auto element_vertices = + xdg_->mesh_manager()->element_vertices(bin_to_mesh_id(bin)); + + xdg::Vertex center {0.0, 0.0, 0.0}; + for (const auto& v : element_vertices) { + center += v; + } + + center /= double(element_vertices.size()); + + if (length_multiplier_ > 0.0) { + center *= length_multiplier_; + } + + return {center[0], center[1], center[2]}; +} + +int XDGMesh::n_vertices() const +{ + return xdg_->mesh_manager()->num_vertices(); +} + +Position XDGMesh::vertex(int id) const +{ + xdg::MeshID mesh_id = xdg_->mesh_manager()->vertex_id(id); + auto v = xdg_->mesh_manager()->vertex_coordinates(mesh_id); + if (length_multiplier_ > 0.0) { + v *= length_multiplier_; + } + return {v[0], v[1], v[2]}; +} + +std::vector XDGMesh::connectivity(int id) const +{ + auto conn = xdg_->mesh_manager()->element_connectivity(bin_to_mesh_id(id)); + for (auto& c : conn) { + c = xdg_->mesh_manager()->vertex_index(c); + } + return conn; +} + +double XDGMesh::volume(int bin) const +{ + double vol = xdg_->mesh_manager()->element_volume(bin_to_mesh_id(bin)); + if (length_multiplier_ > 0.0) { + vol *= length_multiplier_ * length_multiplier_ * length_multiplier_; + } + return vol; +} + +xdg::MeshID XDGMesh::bin_to_mesh_id(int bin) const +{ + return xdg_->mesh_manager()->element_id(bin); +} + +int32_t XDGMesh::mesh_id_to_bin(xdg::MeshID id) const +{ + if (id < 0) + return -1; + return xdg_->mesh_manager()->element_index(id); +} + +} // namespace openmc + +#endif // XDG diff --git a/tests/regression_tests/unstructured_mesh/inputs_true.dat b/tests/regression_tests/unstructured_mesh/inputs_true.dat index e7a485b4d..eecb7af15 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true.dat @@ -48,7 +48,7 @@ 10 - + 1.0 1.0 diff --git a/tests/regression_tests/unstructured_mesh/inputs_true0.dat b/tests/regression_tests/unstructured_mesh/inputs_true0.dat index 2508c25b3..21cfb9a06 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true0.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true0.dat @@ -48,7 +48,7 @@ 10 - + 1.0 1.0 diff --git a/tests/regression_tests/unstructured_mesh/inputs_true1.dat b/tests/regression_tests/unstructured_mesh/inputs_true1.dat index 04bbdb4b2..26d8f95d8 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true1.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true1.dat @@ -48,7 +48,7 @@ 10 - + 1.0 1.0 diff --git a/tests/regression_tests/unstructured_mesh/inputs_true10.dat b/tests/regression_tests/unstructured_mesh/inputs_true10.dat index 169f0a6b4..343d7631d 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true10.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true10.dat @@ -48,7 +48,7 @@ 10 - + 1.0 1.0 diff --git a/tests/regression_tests/unstructured_mesh/inputs_true11.dat b/tests/regression_tests/unstructured_mesh/inputs_true11.dat index 8163750ca..5455d4db7 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true11.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true11.dat @@ -48,7 +48,7 @@ 10 - + 1.0 1.0 diff --git a/tests/regression_tests/unstructured_mesh/inputs_true12.dat b/tests/regression_tests/unstructured_mesh/inputs_true12.dat index be0655e98..7ebd8cff2 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true12.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true12.dat @@ -48,7 +48,7 @@ 10 - + 1.0 1.0 diff --git a/tests/regression_tests/unstructured_mesh/inputs_true13.dat b/tests/regression_tests/unstructured_mesh/inputs_true13.dat index 42c0cc122..d2dcce7b4 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true13.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true13.dat @@ -48,7 +48,7 @@ 10 - + 1.0 1.0 diff --git a/tests/regression_tests/unstructured_mesh/inputs_true14.dat b/tests/regression_tests/unstructured_mesh/inputs_true14.dat index a27c7445d..706ce231f 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true14.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true14.dat @@ -48,7 +48,7 @@ 10 - + 1.0 1.0 diff --git a/tests/regression_tests/unstructured_mesh/inputs_true15.dat b/tests/regression_tests/unstructured_mesh/inputs_true15.dat index c69bd15d0..101da0515 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true15.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true15.dat @@ -48,7 +48,7 @@ 10 - + 1.0 1.0 diff --git a/tests/regression_tests/unstructured_mesh/inputs_true2.dat b/tests/regression_tests/unstructured_mesh/inputs_true2.dat index ba9d477b5..1aae4efe5 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true2.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true2.dat @@ -48,7 +48,7 @@ 10 - + 1.0 1.0 diff --git a/tests/regression_tests/unstructured_mesh/inputs_true3.dat b/tests/regression_tests/unstructured_mesh/inputs_true3.dat index d840bfa09..aba6a10cb 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true3.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true3.dat @@ -48,7 +48,7 @@ 10 - + 1.0 1.0 diff --git a/tests/regression_tests/unstructured_mesh/inputs_true8.dat b/tests/regression_tests/unstructured_mesh/inputs_true8.dat index c04beb211..9c5ff141a 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true8.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true8.dat @@ -48,7 +48,7 @@ 10 - + 1.0 1.0 diff --git a/tests/regression_tests/unstructured_mesh/inputs_true9.dat b/tests/regression_tests/unstructured_mesh/inputs_true9.dat index 111898ff2..0a09889f9 100644 --- a/tests/regression_tests/unstructured_mesh/inputs_true9.dat +++ b/tests/regression_tests/unstructured_mesh/inputs_true9.dat @@ -48,7 +48,7 @@ 10 - + 1.0 1.0 diff --git a/tests/regression_tests/unstructured_mesh/test.py b/tests/regression_tests/unstructured_mesh/test.py index 7607531d8..3be8c1144 100644 --- a/tests/regression_tests/unstructured_mesh/test.py +++ b/tests/regression_tests/unstructured_mesh/test.py @@ -225,11 +225,11 @@ def model(): settings.batches = 10 # source setup - r = openmc.stats.Uniform(a=0.0, b=0.0) + r = openmc.stats.Uniform(a=0.0, b=9.0) cos_theta = openmc.stats.Discrete(x=[1.0], p=[1.0]) phi = openmc.stats.Discrete(x=[0.0], p=[1.0]) - space = openmc.stats.SphericalIndependent(r, cos_theta, phi) + space = openmc.stats.SphericalIndependent(r=r, cos_theta=cos_theta, phi=phi) energy = openmc.stats.Discrete(x=[15.e+06], p=[1.0]) source = openmc.IndependentSource(space=space, energy=energy) settings.source = source diff --git a/tests/regression_tests/xdg/__init__.py b/tests/regression_tests/xdg/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/xdg/tallies/__init__.py b/tests/regression_tests/xdg/tallies/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/xdg/tallies/inputs_true0_libmesh.dat b/tests/regression_tests/xdg/tallies/inputs_true0_libmesh.dat new file mode 100644 index 000000000..7cbf9a6f7 --- /dev/null +++ b/tests/regression_tests/xdg/tallies/inputs_true0_libmesh.dat @@ -0,0 +1,90 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 1000 + 10 + + + + + 1.0 1.0 + + + 0.0 1.0 + + + + 15000000.0 1.0 + + + + + + 10 10 10 + -10.0 -10.0 -10.0 + 10.0 10.0 10.0 + + + test_mesh_tets.e + + + 1 + + + 2 + + + 1 + flux + collision + + + 2 + flux + collision + + + diff --git a/tests/regression_tests/xdg/tallies/inputs_true0_moab.dat b/tests/regression_tests/xdg/tallies/inputs_true0_moab.dat new file mode 100644 index 000000000..b2b9925fc --- /dev/null +++ b/tests/regression_tests/xdg/tallies/inputs_true0_moab.dat @@ -0,0 +1,90 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 1000 + 10 + + + + + 1.0 1.0 + + + 0.0 1.0 + + + + 15000000.0 1.0 + + + + + + 10 10 10 + -10.0 -10.0 -10.0 + 10.0 10.0 10.0 + + + test_mesh_tets.e + + + 1 + + + 2 + + + 1 + flux + collision + + + 2 + flux + collision + + + diff --git a/tests/regression_tests/xdg/tallies/inputs_true1_libmesh.dat b/tests/regression_tests/xdg/tallies/inputs_true1_libmesh.dat new file mode 100644 index 000000000..5511a599c --- /dev/null +++ b/tests/regression_tests/xdg/tallies/inputs_true1_libmesh.dat @@ -0,0 +1,90 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 1000 + 10 + + + + + 1.0 1.0 + + + 0.0 1.0 + + + + 15000000.0 1.0 + + + + + + 10 10 10 + -10.0 -10.0 -10.0 + 10.0 10.0 10.0 + + + test_mesh_tets.exo + + + 1 + + + 2 + + + 1 + flux + collision + + + 2 + flux + collision + + + diff --git a/tests/regression_tests/xdg/tallies/inputs_true2_libmesh.dat b/tests/regression_tests/xdg/tallies/inputs_true2_libmesh.dat new file mode 100644 index 000000000..6d4eb6111 --- /dev/null +++ b/tests/regression_tests/xdg/tallies/inputs_true2_libmesh.dat @@ -0,0 +1,90 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 1000 + 10 + + + + + 1.0 1.0 + + + 0.0 1.0 + + + + 15000000.0 1.0 + + + + + + 10 10 10 + -10.0 -10.0 -10.0 + 10.0 10.0 10.0 + + + test_mesh_tets_w_holes.e + + + 1 + + + 2 + + + 1 + flux + collision + + + 2 + flux + collision + + + diff --git a/tests/regression_tests/xdg/tallies/inputs_true2_moab.dat b/tests/regression_tests/xdg/tallies/inputs_true2_moab.dat new file mode 100644 index 000000000..bb6546381 --- /dev/null +++ b/tests/regression_tests/xdg/tallies/inputs_true2_moab.dat @@ -0,0 +1,90 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 1000 + 10 + + + + + 1.0 1.0 + + + 0.0 1.0 + + + + 15000000.0 1.0 + + + + + + 10 10 10 + -10.0 -10.0 -10.0 + 10.0 10.0 10.0 + + + test_mesh_tets_w_holes.e + + + 1 + + + 2 + + + 1 + flux + collision + + + 2 + flux + collision + + + diff --git a/tests/regression_tests/xdg/tallies/inputs_true3_libmesh.dat b/tests/regression_tests/xdg/tallies/inputs_true3_libmesh.dat new file mode 100644 index 000000000..4fe742d99 --- /dev/null +++ b/tests/regression_tests/xdg/tallies/inputs_true3_libmesh.dat @@ -0,0 +1,90 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 1000 + 10 + + + + + 1.0 1.0 + + + 0.0 1.0 + + + + 15000000.0 1.0 + + + + + + 10 10 10 + -10.0 -10.0 -10.0 + 10.0 10.0 10.0 + + + test_mesh_tets_w_holes.exo + + + 1 + + + 2 + + + 1 + flux + collision + + + 2 + flux + collision + + + diff --git a/tests/regression_tests/xdg/tallies/inputs_true4_libmesh.dat b/tests/regression_tests/xdg/tallies/inputs_true4_libmesh.dat new file mode 100644 index 000000000..c09cb5aa3 --- /dev/null +++ b/tests/regression_tests/xdg/tallies/inputs_true4_libmesh.dat @@ -0,0 +1,90 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 1000 + 10 + + + + + 1.0 1.0 + + + 0.0 1.0 + + + + 15000000.0 1.0 + + + + + + 10 10 10 + -10.0 -10.0 -10.0 + 10.0 10.0 10.0 + + + test_mesh_tets.e + + + 1 + + + 2 + + + 1 + flux + collision + + + 2 + flux + collision + + + diff --git a/tests/regression_tests/xdg/tallies/inputs_true4_moab.dat b/tests/regression_tests/xdg/tallies/inputs_true4_moab.dat new file mode 100644 index 000000000..1f5366957 --- /dev/null +++ b/tests/regression_tests/xdg/tallies/inputs_true4_moab.dat @@ -0,0 +1,90 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 1000 + 10 + + + + + 1.0 1.0 + + + 0.0 1.0 + + + + 15000000.0 1.0 + + + + + + 10 10 10 + -10.0 -10.0 -10.0 + 10.0 10.0 10.0 + + + test_mesh_tets.e + + + 1 + + + 2 + + + 1 + flux + collision + + + 2 + flux + collision + + + diff --git a/tests/regression_tests/xdg/tallies/inputs_true5_libmesh.dat b/tests/regression_tests/xdg/tallies/inputs_true5_libmesh.dat new file mode 100644 index 000000000..945a16272 --- /dev/null +++ b/tests/regression_tests/xdg/tallies/inputs_true5_libmesh.dat @@ -0,0 +1,90 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 1000 + 10 + + + + + 1.0 1.0 + + + 0.0 1.0 + + + + 15000000.0 1.0 + + + + + + 10 10 10 + -10.0 -10.0 -10.0 + 10.0 10.0 10.0 + + + test_mesh_tets.exo + + + 1 + + + 2 + + + 1 + flux + collision + + + 2 + flux + collision + + + diff --git a/tests/regression_tests/xdg/tallies/inputs_true6_libmesh.dat b/tests/regression_tests/xdg/tallies/inputs_true6_libmesh.dat new file mode 100644 index 000000000..ab92115f6 --- /dev/null +++ b/tests/regression_tests/xdg/tallies/inputs_true6_libmesh.dat @@ -0,0 +1,90 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 1000 + 10 + + + + + 1.0 1.0 + + + 0.0 1.0 + + + + 15000000.0 1.0 + + + + + + 10 10 10 + -10.0 -10.0 -10.0 + 10.0 10.0 10.0 + + + test_mesh_tets_w_holes.e + + + 1 + + + 2 + + + 1 + flux + collision + + + 2 + flux + collision + + + diff --git a/tests/regression_tests/xdg/tallies/inputs_true6_moab.dat b/tests/regression_tests/xdg/tallies/inputs_true6_moab.dat new file mode 100644 index 000000000..a0fc43c9b --- /dev/null +++ b/tests/regression_tests/xdg/tallies/inputs_true6_moab.dat @@ -0,0 +1,90 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 1000 + 10 + + + + + 1.0 1.0 + + + 0.0 1.0 + + + + 15000000.0 1.0 + + + + + + 10 10 10 + -10.0 -10.0 -10.0 + 10.0 10.0 10.0 + + + test_mesh_tets_w_holes.e + + + 1 + + + 2 + + + 1 + flux + collision + + + 2 + flux + collision + + + diff --git a/tests/regression_tests/xdg/tallies/inputs_true7_libmesh.dat b/tests/regression_tests/xdg/tallies/inputs_true7_libmesh.dat new file mode 100644 index 000000000..569ab3da2 --- /dev/null +++ b/tests/regression_tests/xdg/tallies/inputs_true7_libmesh.dat @@ -0,0 +1,90 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + fixed source + 1000 + 10 + + + + + 1.0 1.0 + + + 0.0 1.0 + + + + 15000000.0 1.0 + + + + + + 10 10 10 + -10.0 -10.0 -10.0 + 10.0 10.0 10.0 + + + test_mesh_tets_w_holes.exo + + + 1 + + + 2 + + + 1 + flux + collision + + + 2 + flux + collision + + + diff --git a/tests/regression_tests/xdg/tallies/test.py b/tests/regression_tests/xdg/tallies/test.py new file mode 100644 index 000000000..283804f31 --- /dev/null +++ b/tests/regression_tests/xdg/tallies/test.py @@ -0,0 +1,311 @@ +import glob +from pathlib import Path + +import numpy as np +import pytest +import openmc +import openmc.lib +from openmc.xdg import XDGMesh + +from tests.testing_harness import PyAPITestHarness + + +class XDGMeshTallyTest(PyAPITestHarness): + + ELEM_PER_VOXEL = 12 + + def __init__( + self, + statepoint_name, + model, + inputs_true, + mesh_filename, + mesh_kind, + holes=None, + scale_factor=10.0, + ): + super().__init__(statepoint_name, model, inputs_true) + print(f"Running test with mesh file: {mesh_filename} and inputs file: {inputs_true}") + self.mesh_filename = mesh_filename + self.mesh_kind = mesh_kind + self.holes = holes + self.scale_bounding_cell(scale_factor) + + def scale_bounding_cell(self, scale_factor): + geometry = self._model.geometry + for surface in geometry.get_all_surfaces().values(): + if surface.boundary_type != 'vacuum': + continue + for coeff in surface._coefficients: + surface._coefficients[coeff] *= scale_factor + + def _compare_results(self): + with openmc.StatePoint(self._sp_name) as sp: + xdg_mesh = None + for mesh in sp.meshes.values(): + if isinstance(mesh, XDGMesh): + xdg_mesh = mesh + break + + assert xdg_mesh is not None + assert Path(xdg_mesh.filename).name == Path(self.mesh_filename).name + + if self.mesh_kind == "tet": + exp_vertex = (-10.0, -10.0, -10.0) + exp_centroid = (-8.75, -9.75, -9.25) + else: + exp_vertex = (-10.0, -10.0, 10.0) + exp_centroid = (-9.0, -9.0, 9.0) + + np.testing.assert_array_equal(xdg_mesh.vertices[0], exp_vertex) + np.testing.assert_array_equal(xdg_mesh.centroid(0), exp_centroid) + + reg_mesh_data = None + xdg_mesh_data = None + xdg_tally = None + for tally in sp.tallies.values(): + if tally.contains_filter(openmc.MeshFilter): + flt = tally.find_filter(openmc.MeshFilter) + if isinstance(flt.mesh, openmc.RegularMesh): + reg_mesh_data = self.get_mesh_tally_data(tally) + else: + xdg_tally = tally + xdg_mesh_data = self.get_mesh_tally_data(tally, structured=True) + + assert reg_mesh_data is not None + assert xdg_mesh_data is not None + + if self.holes: + reg_mesh_data = np.delete(reg_mesh_data, self.holes) + + decimals = 10 if xdg_tally.estimator == 'collision' else 8 + np.testing.assert_array_almost_equal( + np.sort(xdg_mesh_data), + np.sort(reg_mesh_data), + decimals + ) + + def get_mesh_tally_data(self, tally, structured=False): + data = tally.get_reshaped_data(value='mean') + if structured: + data = data.reshape((-1, self.ELEM_PER_VOXEL)) + else: + data.shape = (data.size, 1) + return np.sum(data, axis=1) + + def update_results(self): + """Update inputs_true.dat without storing results_true.dat.""" + try: + self._build_inputs() + inputs = self._get_inputs() + self._write_inputs(inputs) + self._overwrite_inputs() + self._run_openmc() + self._test_output_created() + finally: + self._cleanup() + + def _cleanup(self): + super()._cleanup() + output = glob.glob('tally*.vtk') + output += glob.glob('tally*.e') + for f in output: + Path(f).unlink(missing_ok=True) + + +@pytest.fixture +def model(): + openmc.reset_auto_ids() + + model = openmc.Model() + + ### Materials ### + materials = openmc.Materials() + + fuel_mat = openmc.Material(name="fuel") + fuel_mat.add_nuclide("U235", 1.0) + fuel_mat.set_density('g/cc', 4.5) + materials.append(fuel_mat) + + zirc_mat = openmc.Material(name="zircaloy") + zirc_mat.add_element("Zr", 1.0) + zirc_mat.set_density("g/cc", 5.77) + materials.append(zirc_mat) + + water_mat = openmc.Material(name="water") + water_mat.add_nuclide("H1", 2.0) + water_mat.add_nuclide("O16", 1.0) + water_mat.set_density("atom/b-cm", 0.07416) + materials.append(water_mat) + + model.materials = materials + + ### Geometry ### + fuel_boundary = openmc.model.RectangularParallelepiped(*[-5.0, 5.0]*3) + clad_boundary = openmc.model.RectangularParallelepiped(*[-6.0, 6.0]*3) + + fuel_cell = openmc.Cell(name="fuel") + fuel_cell.region = -fuel_boundary + fuel_cell.fill = fuel_mat + + clad_cell = openmc.Cell(name="clad") + clad_cell.region = +fuel_boundary & -clad_boundary + clad_cell.fill = zirc_mat + + # set bounding cell dimension to one + # this will be updated later according to the test case parameters + boundary = openmc.model.RectangularParallelepiped(*[-1.0, 1.0]*3, + boundary_type='vacuum') + + water_cell = openmc.Cell(name="water") + water_cell.region = +clad_boundary & -boundary + water_cell.fill = water_mat + + # create a containing universe + model.geometry = openmc.Geometry([fuel_cell, clad_cell, water_cell]) + + ### Reference Tally ### + + # create meshes and mesh filters + regular_mesh = openmc.RegularMesh() + regular_mesh.dimension = (10, 10, 10) + regular_mesh.lower_left = (-10.0, -10.0, -10.0) + regular_mesh.upper_right = (10.0, 10.0, 10.0) + + regular_mesh_filter = openmc.MeshFilter(mesh=regular_mesh) + regular_mesh_tally = openmc.Tally(name="regular mesh tally") + regular_mesh_tally.filters = [regular_mesh_filter] + regular_mesh_tally.scores = ['flux'] + + model.tallies = openmc.Tallies([regular_mesh_tally]) + + ### Settings ### + model.settings.run_mode = 'fixed source' + model.settings.particles = 1000 + model.settings.batches = 10 + + # source setup + r = openmc.stats.Uniform(a=0.0, b=9.0) + cos_theta = openmc.stats.delta_function(1.0) + phi = openmc.stats.delta_function(0.0) + + space = openmc.stats.SphericalIndependent(r, cos_theta, phi) + energy = openmc.stats.delta_function(15e6) + source = openmc.IndependentSource(space=space, energy=energy) + model.settings.source = source + + return model + + +MESH_CASES = ( + { + "mesh_filename": "test_mesh_tets.e", + "mesh_kind": "tet", + "holes": None, + "elem_per_voxel": 12, + "external_geom": False, + "libraries": ("moab", "libmesh"), + }, + { + "mesh_filename": "test_mesh_tets.exo", + "mesh_kind": "tet", + "holes": None, + "elem_per_voxel": 12, + "external_geom": False, + "libraries": ("libmesh",), + }, + { + "mesh_filename": "test_mesh_tets_w_holes.e", + "mesh_kind": "tet", + "holes": (333, 90, 77), + "elem_per_voxel": 12, + "external_geom": False, + "libraries": ("moab", "libmesh"), + }, + { + "mesh_filename": "test_mesh_tets_w_holes.exo", + "mesh_kind": "tet", + "holes": (333, 90, 77), + "elem_per_voxel": 12, + "external_geom": False, + "libraries": ("libmesh",), + }, + { + "mesh_filename": "test_mesh_tets.e", + "mesh_kind": "tet", + "holes": None, + "elem_per_voxel": 12, + "external_geom": True, + "libraries": ("moab", "libmesh"), + }, + { + "mesh_filename": "test_mesh_tets.exo", + "mesh_kind": "tet", + "holes": None, + "elem_per_voxel": 12, + "external_geom": True, + "libraries": ("libmesh",), + }, + { + "mesh_filename": "test_mesh_tets_w_holes.e", + "mesh_kind": "tet", + "holes": (333, 90, 77), + "elem_per_voxel": 12, + "external_geom": True, + "libraries": ("moab", "libmesh"), + }, + { + "mesh_filename": "test_mesh_tets_w_holes.exo", + "mesh_kind": "tet", + "holes": (333, 90, 77), + "elem_per_voxel": 12, + "external_geom": True, + "libraries": ("libmesh",), + } +) + +test_cases = [] +test_case_ids = [] +for i, mesh_case in enumerate(MESH_CASES): + for library in mesh_case["libraries"]: + test_cases.append({ + "inputs_true": f"inputs_true{i}_{library}.dat", + "library": library, + **{k: v for k, v in mesh_case.items() if k != "libraries"}, + }) + mesh_stem = Path(mesh_case["mesh_filename"]).stem + geom_mode = "external-geom" if mesh_case["external_geom"] else "bounded-geom" + holes = "holes" if mesh_case["holes"] else "solid" + test_case_ids.append(f"{library}-{mesh_stem}-{holes}-{geom_mode}") + + +@pytest.mark.parametrize("test_opts", test_cases, ids=test_case_ids) +def test_xdg_mesh_tallies(model, test_opts): + if not openmc.lib._xdg_enabled(): + pytest.skip("XDG is not enabled in this build.") + + # reference mesh tally + regular_mesh_tally = model.tallies[0] + regular_mesh_tally.estimator = 'collision' + + # add analogous XDG mesh tally + xdg_mesh = XDGMesh(test_opts["mesh_filename"], test_opts["library"]) + xdg_filter = openmc.MeshFilter(mesh=xdg_mesh) + + xdg_tally = openmc.Tally(name="xdg mesh tally") + xdg_tally.filters = [xdg_filter] + xdg_tally.scores = ['flux'] + xdg_tally.estimator = 'collision' + model.tallies.append(xdg_tally) + + harness = XDGMeshTallyTest( + 'statepoint.10.h5', + model, + test_opts["inputs_true"], + test_opts["mesh_filename"], + test_opts["mesh_kind"], + test_opts["holes"], + scale_factor=15.0 if test_opts["external_geom"] else 10.0, + ) + harness.ELEM_PER_VOXEL = test_opts["elem_per_voxel"] + harness.main() diff --git a/tests/regression_tests/xdg/tallies/test_mesh_hexes.e b/tests/regression_tests/xdg/tallies/test_mesh_hexes.e new file mode 120000 index 000000000..32a7bc68f --- /dev/null +++ b/tests/regression_tests/xdg/tallies/test_mesh_hexes.e @@ -0,0 +1 @@ +../../unstructured_mesh/test_mesh_hexes.e \ No newline at end of file diff --git a/tests/regression_tests/xdg/tallies/test_mesh_hexes.exo b/tests/regression_tests/xdg/tallies/test_mesh_hexes.exo new file mode 120000 index 000000000..723a379f7 --- /dev/null +++ b/tests/regression_tests/xdg/tallies/test_mesh_hexes.exo @@ -0,0 +1 @@ +../../unstructured_mesh/test_mesh_hexes.exo \ No newline at end of file diff --git a/tests/regression_tests/xdg/tallies/test_mesh_tets.e b/tests/regression_tests/xdg/tallies/test_mesh_tets.e new file mode 120000 index 000000000..65723ab81 --- /dev/null +++ b/tests/regression_tests/xdg/tallies/test_mesh_tets.e @@ -0,0 +1 @@ +../../unstructured_mesh/test_mesh_tets.e \ No newline at end of file diff --git a/tests/regression_tests/xdg/tallies/test_mesh_tets.exo b/tests/regression_tests/xdg/tallies/test_mesh_tets.exo new file mode 120000 index 000000000..676f4d799 --- /dev/null +++ b/tests/regression_tests/xdg/tallies/test_mesh_tets.exo @@ -0,0 +1 @@ +../../unstructured_mesh/test_mesh_tets.exo \ No newline at end of file diff --git a/tests/regression_tests/xdg/tallies/test_mesh_tets_w_holes.e b/tests/regression_tests/xdg/tallies/test_mesh_tets_w_holes.e new file mode 120000 index 000000000..2dae88efc --- /dev/null +++ b/tests/regression_tests/xdg/tallies/test_mesh_tets_w_holes.e @@ -0,0 +1 @@ +../../unstructured_mesh/test_mesh_tets_w_holes.e \ No newline at end of file diff --git a/tests/regression_tests/xdg/tallies/test_mesh_tets_w_holes.exo b/tests/regression_tests/xdg/tallies/test_mesh_tets_w_holes.exo new file mode 120000 index 000000000..fbd15ce88 --- /dev/null +++ b/tests/regression_tests/xdg/tallies/test_mesh_tets_w_holes.exo @@ -0,0 +1 @@ +../../unstructured_mesh/test_mesh_tets_w_holes.exo \ No newline at end of file diff --git a/tools/ci/gha-build-moab.sh b/tools/ci/gha-build-moab.sh new file mode 100755 index 000000000..e9ebd75a7 --- /dev/null +++ b/tools/ci/gha-build-moab.sh @@ -0,0 +1,29 @@ +#!/bin/bash +set -ex + +MOAB_BRANCH='5.5.1' +MOAB_REPO='https://bitbucket.org/fathomteam/moab/' +MOAB_INSTALL_DIR=$HOME/MOAB/ + +pushd $HOME + +if [[ -d $MOAB_INSTALL_DIR/lib ]]; then + popd + exit 0 +fi + +mkdir -p MOAB +cd MOAB +git clone -b $MOAB_BRANCH $MOAB_REPO moab +mkdir build +cd build +cmake ../moab \ + -DENABLE_HDF5=ON \ + -DENABLE_NETCDF=ON \ + -DBUILD_SHARED_LIBS=ON \ + -DCMAKE_INSTALL_PREFIX=$MOAB_INSTALL_DIR +make -j4 +make -j4 install +rm -rf $HOME/MOAB/moab $HOME/MOAB/build + +popd diff --git a/tools/ci/gha-install-dagmc.sh b/tools/ci/gha-install-dagmc.sh index 8d0648a50..21133c6d3 100755 --- a/tools/ci/gha-install-dagmc.sh +++ b/tools/ci/gha-install-dagmc.sh @@ -2,34 +2,23 @@ #!/bin/bash set -ex -# MOAB Variables -MOAB_BRANCH='5.5.1' -MOAB_REPO='https://bitbucket.org/fathomteam/moab/' -MOAB_INSTALL_DIR=$HOME/MOAB/ - # DAGMC Variables DAGMC_BRANCH='develop' DAGMC_REPO='https://github.com/svalinn/dagmc' DAGMC_INSTALL_DIR=$HOME/DAGMC/ -CURRENT_DIR=$(pwd) - -# MOAB Install -cd $HOME -mkdir MOAB && cd MOAB -git clone -b $MOAB_BRANCH $MOAB_REPO -mkdir build && cd build -cmake ../moab -DENABLE_HDF5=ON -DENABLE_NETCDF=ON -DBUILD_SHARED_LIBS=ON -DCMAKE_INSTALL_PREFIX=$MOAB_INSTALL_DIR -make -j && make -j install -rm -rf $HOME/MOAB/moab $HOME/MOAB/build - # DAGMC Install -cd $HOME -mkdir DAGMC && cd DAGMC -git clone -b $DAGMC_BRANCH $DAGMC_REPO -mkdir build && cd build -cmake ../dagmc -DBUILD_TALLY=ON -DCMAKE_INSTALL_PREFIX=$DAGMC_INSTALL_DIR -DBUILD_STATIC_LIBS=OFF -DMOAB_DIR=$MOAB_INSTALL_DIR -make -j install +pushd $HOME +mkdir -p DAGMC +cd DAGMC +git clone -b $DAGMC_BRANCH $DAGMC_REPO dagmc +mkdir build +cd build +cmake ../dagmc \ + -DBUILD_TALLY=ON \ + -DCMAKE_INSTALL_PREFIX=$DAGMC_INSTALL_DIR \ + -DBUILD_STATIC_LIBS=OFF \ + -DMOAB_DIR=$HOME/MOAB +make -j4 install rm -rf $HOME/DAGMC/dagmc $HOME/DAGMC/build - -cd $CURRENT_DIR +popd diff --git a/tools/ci/gha-install-xdg.sh b/tools/ci/gha-install-xdg.sh new file mode 100755 index 000000000..f960739ae --- /dev/null +++ b/tools/ci/gha-install-xdg.sh @@ -0,0 +1,43 @@ +#!/bin/bash +set -ex + +XDG_BRANCH='main' +XDG_REPO='https://github.com/xdg-org/xdg.git' +XDG_INSTALL_DIR=$HOME/XDG/ + +pushd $HOME + +mkdir -p XDG +cd XDG +git clone -b $XDG_BRANCH --recurse-submodules $XDG_REPO xdg +mkdir build +cd build + +cmake_args=( + ../xdg + -DCMAKE_INSTALL_PREFIX=$XDG_INSTALL_DIR + -DXDG_BUILD_TESTS=OFF + -DXDG_BUILD_TOOLS=ON + -DXDG_ENABLE_MOAB=ON + -DMOAB_DIR=$HOME/MOAB +) + +if [[ $LIBMESH == 'y' || $XDG == 'y' ]]; then + cmake_args+=( + -DXDG_ENABLE_LIBMESH=ON + -DCMAKE_PREFIX_PATH=$HOME/LIBMESH + ) +else + cmake_args+=(-DXDG_ENABLE_LIBMESH=OFF) +fi + +if [[ $MPI == 'y' ]]; then + cmake_args+=(-DXDG_LINK_MPI=ON) +fi + +cmake "${cmake_args[@]}" +make -j4 +make -j4 install +rm -rf $HOME/XDG/xdg $HOME/XDG/build + +popd diff --git a/tools/ci/gha-install.py b/tools/ci/gha-install.py index 5488b9547..e352ce3a0 100644 --- a/tools/ci/gha-install.py +++ b/tools/ci/gha-install.py @@ -3,7 +3,8 @@ import shutil import subprocess -def install(omp=False, mpi=False, phdf5=False, dagmc=False, libmesh=False): +def install(omp=False, mpi=False, phdf5=False, dagmc=False, xdg=False, + libmesh=False): # Create build directory and change to it shutil.rmtree('build', ignore_errors=True) os.mkdir('build') @@ -29,16 +30,26 @@ def install(omp=False, mpi=False, phdf5=False, dagmc=False, libmesh=False): else: cmake_cmd.append('-DHDF5_PREFER_PARALLEL=OFF') + prefix_paths = [] + if dagmc: cmake_cmd.append('-DOPENMC_USE_DAGMC=ON') cmake_cmd.append('-DOPENMC_USE_UWUW=ON') dagmc_path = os.environ.get('HOME') + '/DAGMC' - cmake_cmd.append('-DCMAKE_PREFIX_PATH=' + dagmc_path) + prefix_paths.append(dagmc_path) - if libmesh: + if xdg: + cmake_cmd.append('-DOPENMC_USE_XDG=ON') + xdg_path = os.environ.get('HOME') + '/XDG' + prefix_paths.append(xdg_path) + + if libmesh or xdg: cmake_cmd.append('-DOPENMC_USE_LIBMESH=ON') libmesh_path = os.environ.get('HOME') + '/LIBMESH' - cmake_cmd.append('-DCMAKE_PREFIX_PATH=' + libmesh_path) + prefix_paths.append(libmesh_path) + + if prefix_paths: + cmake_cmd.append('-DCMAKE_PREFIX_PATH=' + ';'.join(prefix_paths)) # Build in coverage mode for coverage testing cmake_cmd.append('-DOPENMC_ENABLE_COVERAGE=on') @@ -59,10 +70,11 @@ def main(): mpi = (os.environ.get('MPI') == 'y') phdf5 = (os.environ.get('PHDF5') == 'y') dagmc = (os.environ.get('DAGMC') == 'y') + xdg = (os.environ.get('XDG') == 'y') libmesh = (os.environ.get('LIBMESH') == 'y') # Build and install - install(omp, mpi, phdf5, dagmc, libmesh) + install(omp, mpi, phdf5, dagmc, xdg, libmesh) if __name__ == '__main__': main() diff --git a/tools/ci/gha-install.sh b/tools/ci/gha-install.sh index 4cf62afbb..8a125d908 100755 --- a/tools/ci/gha-install.sh +++ b/tools/ci/gha-install.sh @@ -9,20 +9,30 @@ pip install --upgrade numpy # Install NJOY 2016 ./tools/ci/gha-install-njoy.sh -# Install DAGMC if needed -if [[ $DAGMC = 'y' ]]; then - ./tools/ci/gha-install-dagmc.sh +# Build MOAB if needed +if [[ $DAGMC = 'y' || $XDG = 'y' ]]; then + ./tools/ci/gha-build-moab.sh fi # Install NCrystal and verify installation pip install 'ncrystal>=4.1.0' nctool --test +# Install DAGMC if needed +if [[ $DAGMC = 'y' ]]; then + ./tools/ci/gha-install-dagmc.sh +fi + # Install libMesh if needed -if [[ $LIBMESH = 'y' ]]; then +if [[ $LIBMESH = 'y' || $XDG = 'y' ]]; then ./tools/ci/gha-install-libmesh.sh fi +# Install XDG if needed +if [[ $XDG = 'y' ]]; then + ./tools/ci/gha-install-xdg.sh +fi + # Install MCPL pip install mcpl