diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 6dff76838..0431c95ab 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -504,7 +504,10 @@ attributes/sub-elements: independent distributions of r-, cos_theta-, and phi-coordinates where cos_theta is the cosine of the angle with respect to the z-axis, phi is the azimuthal angle, and the sphere is centered on the coordinate - (x0,y0,z0). + (x0,y0,z0). A "mesh" spatial distribution samples source sites from a mesh element + based on the relative strengths provided in the node. Source locations + within an element are sampled isotropically. If no strengths are provided, + the space within the mesh is uniformly sampled. *Default*: None diff --git a/docs/source/pythonapi/stats.rst b/docs/source/pythonapi/stats.rst index fb6383fc7..ffb4fcda2 100644 --- a/docs/source/pythonapi/stats.rst +++ b/docs/source/pythonapi/stats.rst @@ -57,6 +57,7 @@ Spatial Distributions openmc.stats.SphericalIndependent openmc.stats.Box openmc.stats.Point + openmc.stats.MeshSpatial .. autosummary:: :toctree: generated diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 760a08e90..58dbe7f0e 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -187,7 +187,9 @@ The spatial distribution can be set equal to a sub-class of :class:`openmc.stats.CartesianIndependent`. To independently specify distributions using spherical or cylindrical coordinates, you can use :class:`openmc.stats.SphericalIndependent` or -:class:`openmc.stats.CylindricalIndependent`, respectively. +:class:`openmc.stats.CylindricalIndependent`, respectively. Meshes can also be +used to represent spatial distributions with :class:`openmc.stats.MeshSpatial` +by specifying a mesh and source strengths for each mesh element. The angular distribution can be set equal to a sub-class of :class:`openmc.stats.UnitSphere` such as :class:`openmc.stats.Isotropic`, diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index 5e426b878..dda23927e 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -4,6 +4,7 @@ #include "pugixml.hpp" #include "openmc/distribution.h" +#include "openmc/mesh.h" #include "openmc/position.h" namespace openmc { @@ -94,6 +95,32 @@ private: Position origin_; //!< Cartesian coordinates of the sphere center }; +//============================================================================== +//! Distribution of points within a mesh +//============================================================================== + +class MeshSpatial : public SpatialDistribution { +public: + explicit MeshSpatial(pugi::xml_node node); + + //! Sample a position from the distribution + //! \param seed Pseudorandom number seed pointer + //! \return Sampled position + Position sample(uint64_t* seed) const; + + const Mesh* mesh() const { return model::meshes.at(mesh_idx_).get(); } + + int32_t n_sources() const { return this->mesh()->n_bins(); } + +private: + int32_t mesh_idx_ {C_NONE}; + double total_strength_ {0.0}; + // TODO: move to an independent class in the future that's similar + // to a discrete distribution without outcomes + std::vector mesh_CDF_; + std::vector mesh_strengths_; +}; + //============================================================================== //! Uniform distribution of points over a box //============================================================================== diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index eece0f1ff..7abc8a832 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -7,8 +7,8 @@ #include #include "hdf5.h" -#include "xtensor/xtensor.hpp" #include "pugixml.hpp" +#include "xtensor/xtensor.hpp" #include "openmc/memory.h" // for unique_ptr #include "openmc/particle.h" @@ -41,7 +41,7 @@ namespace openmc { // Constants //============================================================================== -enum class ElementType { UNSUPPORTED=-1, LINEAR_TET, LINEAR_HEX }; +enum class ElementType { UNSUPPORTED = -1, LINEAR_TET, LINEAR_HEX }; //============================================================================== // Global variables @@ -75,6 +75,19 @@ public: // Methods + //! Sample a mesh volume using a certain seed + // + //! \param[in] seed Seed to use for random sampling + //! \param[in] bin Bin value of the tet sampled + //! \return sampled position within tet + virtual Position sample(uint64_t* seed, int32_t bin) const = 0; + + //! Get the volume of a mesh bin + // + //! \param[in] bin Bin to return the volume for + //! \return Volume of the bin + virtual double volume(int bin) const = 0; + //! Determine which bins were crossed by a particle // //! \param[in] r0 Previous position of the particle @@ -160,6 +173,10 @@ public: } }; + Position sample(uint64_t* seed, int32_t bin) const override; + + double volume(int bin) const override; + int get_bin(Position r) const override; int n_bins() const override; @@ -461,6 +478,7 @@ public: 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; @@ -509,10 +527,6 @@ public: //! \return element connectivity as IDs of the vertices virtual std::vector connectivity(int id) const = 0; - //! Get the volume of a mesh bin - // - //! \param[in] bin Bin to return the volume for - //! \return Volume of the bin virtual double volume(int bin) const = 0; //! Get the library used for this unstructured mesh @@ -523,6 +537,8 @@ public: true}; //!< Write tallies onto the unstructured mesh at the end of a run std::string filename_; //!< Path to unstructured mesh file + ElementType element_type(int bin) const; + protected: //! Set the length multiplier to apply to each point in the mesh void set_length_multiplier(const double length_multiplier); @@ -532,6 +548,14 @@ protected: 1.0}; //!< Constant multiplication factor to apply to mesh coordinates bool specified_length_multiplier_ {false}; + //! Sample barycentric coordinates given a seed and the vertex positions and + //! return the sampled position + // + //! \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; + private: //! Setup method for the mesh. Builds data structures, //! sets up element mapping, creates bounding boxes, etc. @@ -552,6 +576,8 @@ public: // Overridden Methods + Position sample(uint64_t* seed, int32_t bin) const override; + void bins_crossed(Position r0, Position r1, const Direction& u, vector& bins, vector& lengths) const override; @@ -687,7 +713,7 @@ private: std::pair get_score_tags(std::string score) const; // Data members - moab::Range ehs_; //!< Range of tetrahedra EntityHandle's in the mesh + moab::Range ehs_; //!< Range of tetrahedra EntityHandle's in the mesh moab::Range verts_; //!< Range of vertex EntityHandle's in the mesh moab::EntityHandle tetset_; //!< EntitySet containing all tetrahedra moab::EntityHandle kdtree_root_; //!< Root of the MOAB KDTree @@ -705,8 +731,8 @@ class LibMesh : public UnstructuredMesh { public: // Constructors LibMesh(pugi::xml_node node); - LibMesh(const std::string & filename, double length_multiplier = 1.0); - LibMesh(libMesh::MeshBase & input_mesh, double length_multiplier = 1.0); + LibMesh(const std::string& filename, double length_multiplier = 1.0); + LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier = 1.0); static const std::string mesh_lib_type; @@ -714,6 +740,8 @@ public: void bins_crossed(Position r0, Position r1, const Direction& u, vector& bins, vector& lengths) const override; + Position sample(uint64_t* seed, int32_t bin) const override; + int get_bin(Position r) const override; int n_bins() const override; @@ -759,8 +787,11 @@ private: int get_bin_from_element(const libMesh::Elem* elem) const; // Data members - unique_ptr unique_m_ = nullptr; //!< pointer to the libMesh MeshBase instance, only used if mesh is created inside OpenMC - libMesh::MeshBase* m_; //!< pointer to libMesh MeshBase instance, always set during intialization + unique_ptr unique_m_ = + nullptr; //!< pointer to the libMesh MeshBase instance, only used if mesh is + //!< created inside OpenMC + libMesh::MeshBase* m_; //!< pointer to libMesh MeshBase instance, always set + //!< during intialization vector> pl_; //!< per-thread point locators unique_ptr diff --git a/openmc/mesh.py b/openmc/mesh.py index dda7257f4..904d4a2ad 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -1,5 +1,6 @@ from abc import ABC, abstractmethod from collections.abc import Iterable +from collections import OrderedDict from math import pi from numbers import Real, Integral from pathlib import Path @@ -188,7 +189,7 @@ class StructuredMesh(MeshBase): Returns a numpy.ndarray representing the mesh element centroid coordinates with a shape equal to (ndim, dim1, ..., dimn). Can be unpacked along the first dimension with xx, yy, zz = mesh.centroids. - + """ ndim = self.n_dimension @@ -1943,3 +1944,25 @@ class UnstructuredMesh(MeshBase): length_multiplier = float(get_text(elem, 'length_multiplier', 1.0)) return cls(filename, library, mesh_id, '', length_multiplier) + + +def _read_meshes(elem): + """Generate dictionary of meshes from a given XML node + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + dict + A dictionary with mesh IDs as keys and openmc.MeshBase + instanaces as values + """ + out = dict() + for mesh_elem in elem.findall('mesh'): + mesh = MeshBase.from_xml_element(mesh_elem) + out[mesh.id] = mesh + + return out \ No newline at end of file diff --git a/openmc/settings.py b/openmc/settings.py index 22a83509f..516cbbc4a 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -11,10 +11,12 @@ from typing import Optional from xml.etree import ElementTree as ET import openmc.checkvalue as cv +from openmc.stats.multivariate import MeshSpatial from . import RegularMesh, Source, VolumeCalculation, WeightWindows from ._xml import clean_indentation, get_text, reorder_attributes from openmc.checkvalue import PathLike +from .mesh import _read_meshes class RunMode(Enum): @@ -990,6 +992,10 @@ class Settings: def _create_source_subelement(self, root): for source in self.source: root.append(source.to_xml_element()) + if isinstance(source.space, MeshSpatial): + path = f"./mesh[@id='{source.space.mesh.id}']" + if root.find(path) is None: + root.append(source.space.mesh.to_xml_element()) def _create_volume_calcs_subelement(self, root): for calc in self.volume_calculations: @@ -1340,9 +1346,11 @@ class Settings: threshold = float(get_text(elem, 'threshold')) self.keff_trigger = {'type': trigger, 'threshold': threshold} - def _source_from_xml_element(self, root): + def _source_from_xml_element(self, root, meshes=None): for elem in root.findall('source'): - self.source.append(Source.from_xml_element(elem)) + src = Source.from_xml_element(elem, meshes) + # add newly constructed source object to the list + self.source.append(src) def _volume_calcs_from_xml_element(self, root): volume_elems = root.findall("volume_calc") @@ -1721,7 +1729,7 @@ class Settings: settings._rel_max_lost_particles_from_xml_element(elem) settings._generations_per_batch_from_xml_element(elem) settings._keff_trigger_from_xml_element(elem) - settings._source_from_xml_element(elem) + settings._source_from_xml_element(elem, meshes) settings._volume_calcs_from_xml_element(elem) settings._output_from_xml_element(elem) settings._statepoint_from_xml_element(elem) @@ -1781,4 +1789,5 @@ class Settings: """ tree = ET.parse(path) root = tree.getroot() - return cls.from_xml_element(root) + meshes = _read_meshes(root) + return cls.from_xml_element(root, meshes) diff --git a/openmc/source.py b/openmc/source.py index 9293a5953..7f88bcd52 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -258,13 +258,16 @@ class Source: return element @classmethod - def from_xml_element(cls, elem: ET.Element) -> 'openmc.Source': + def from_xml_element(cls, elem: ET.Element, meshes=None) -> 'openmc.Source': """Generate source from an XML element Parameters ---------- elem : xml.etree.ElementTree.Element XML element + meshes : dict + Dictionary with mesh IDs as keys and openmc.MeshBase instaces as + values Returns ------- @@ -313,7 +316,7 @@ class Source: space = elem.find('space') if space is not None: - source.space = Spatial.from_xml_element(space) + source.space = Spatial.from_xml_element(space, meshes) angle = elem.find('angle') if angle is not None: diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 8bb319aa4..7372020b0 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -9,6 +9,7 @@ import numpy as np import openmc.checkvalue as cv from .._xml import get_text from .univariate import Univariate, Uniform, PowerLaw +from ..mesh import MeshBase class UnitSphere(ABC): @@ -261,7 +262,7 @@ class Spatial(ABC): @classmethod @abstractmethod - def from_xml_element(cls, elem): + def from_xml_element(cls, elem, meshes=None): distribution = get_text(elem, 'type') if distribution == 'cartesian': return CartesianIndependent.from_xml_element(elem) @@ -273,6 +274,8 @@ class Spatial(ABC): return Box.from_xml_element(elem) elif distribution == 'point': return Point.from_xml_element(elem) + elif distribution == 'mesh': + return MeshSpatial.from_xml_element(elem, meshes) class CartesianIndependent(Spatial): @@ -617,6 +620,133 @@ class CylindricalIndependent(Spatial): return cls(r, phi, z, origin=origin) +class MeshSpatial(Spatial): + """Spatial distribution for a mesh. + + This distribution specifies a mesh to sample over with source strengths + specified for each mesh element. + + .. versionadded:: 0.13.3 + + Parameters + ---------- + mesh : openmc.MeshBase + The mesh instance used for sampling + strengths : iterable of float, optional + An iterable of values that represents the weights of each element. If no + source strengths are specified, they will be equal for all mesh + elements. + volume_normalized : bool, optional + Whether or not the strengths will be multiplied by element volumes at + runtime. Default is True. + + Attributes + ---------- + mesh : openmc.MeshBase + The mesh instance used for sampling + strengths : numpy.ndarray or None + An array of source strengths for each mesh element + volume_normalized : bool + Whether or not the strengths will be multiplied by element volumes at + runtime. + """ + + def __init__(self, mesh, strengths=None, volume_normalized=True): + self.mesh = mesh + self.strengths = strengths + self.volume_normalized = volume_normalized + + @property + def mesh(self): + return self._mesh + + @mesh.setter + def mesh(self, mesh): + if mesh is not None: + cv.check_type('mesh instance', mesh, MeshBase) + self._mesh = mesh + + @property + def volume_normalized(self): + return self._volume_normalized + + @volume_normalized.setter + def volume_normalized(self, volume_normalized): + cv.check_type('Multiply strengths by element volumes', volume_normalized, bool) + self._volume_normalized = volume_normalized + + @property + def strengths(self): + return self._strengths + + @strengths.setter + def strengths(self, given_strengths): + if given_strengths is not None: + cv.check_type('strengths array passed in', given_strengths, Iterable, Real) + self._strengths = np.asarray(given_strengths, dtype=float).flatten() + else: + self._strengths = None + + @property + def num_strength_bins(self): + if self.strengths is None: + raise ValueError('Strengths are not set') + return self.strengths.size + + def to_xml_element(self): + """Return XML representation of the spatial distribution + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing spatial distribution data + + """ + element = ET.Element('space') + element.set('type', 'mesh') + element.set("mesh_id", str(self.mesh.id)) + element.set("volume_normalized", str(self.volume_normalized)) + + if self.strengths is not None: + subelement = ET.SubElement(element, 'strengths') + subelement.text = ' '.join(str(e) for e in self.strengths) + + return element + + @classmethod + def from_xml_element(cls, elem, meshes): + """Generate spatial distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + meshes : dict + A dictionary with mesh IDs as keys and openmc.MeshBase instances as + values + + Returns + ------- + openmc.stats.MeshSpatial + Spatial distribution generated from XML element + + """ + + mesh_id = int(elem.get('mesh_id')) + + # check if this mesh has been read in from another location already + if mesh_id not in meshes: + raise RuntimeError(f'Could not locate mesh with ID "{mesh_id}"') + + volume_normalized = elem.get("volume_normalized") + volume_normalized = get_text(elem, 'volume_normalized').lower() == 'true' + strengths = get_text(elem, 'strengths') + if strengths is not None: + strengths = [float(b) for b in get_text(elem, 'strengths').split()] + + return cls(meshes[mesh_id], strengths, volume_normalized) + + class Box(Spatial): """Uniform distribution of coordinates in a rectangular cuboid. diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index a06f84d80..fb1fbf099 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -1,7 +1,9 @@ #include "openmc/distribution_spatial.h" #include "openmc/error.h" +#include "openmc/mesh.h" #include "openmc/random_lcg.h" +#include "openmc/search.h" #include "openmc/xml_interface.h" namespace openmc { @@ -137,7 +139,8 @@ SphericalIndependent::SphericalIndependent(pugi::xml_node node) pugi::xml_node node_dist = node.child("cos_theta"); cos_theta_ = distribution_from_xml(node_dist); } else { - // If no distribution was specified, default to a single point at cos_theta=0 + // If no distribution was specified, default to a single point at + // cos_theta=0 double x[] {0.0}; double p[] {1.0}; cos_theta_ = make_unique(x, p, 1); @@ -180,6 +183,84 @@ Position SphericalIndependent::sample(uint64_t* seed) const return {x, y, z}; } +//============================================================================== +// MeshSpatial implementation +//============================================================================== + +MeshSpatial::MeshSpatial(pugi::xml_node node) +{ + // No in-tet distributions implemented, could include distributions for the + // barycentric coords Read in unstructured mesh from mesh_id value + int32_t mesh_id = std::stoi(get_node_value(node, "mesh_id")); + // Get pointer to spatial distribution + mesh_idx_ = model::mesh_map.at(mesh_id); + + auto mesh_ptr = + dynamic_cast(model::meshes.at(mesh_idx_).get()); + if (!mesh_ptr) { + fatal_error("Only unstructured mesh is supported for source sampling."); + } + + // ensure that the unstructured mesh contains only linear tets + for (int bin = 0; bin < mesh_ptr->n_bins(); bin++) { + if (mesh_ptr->element_type(bin) != ElementType::LINEAR_TET) { + fatal_error( + "Mesh specified for source must contain only linear tetrahedra."); + } + } + + int32_t n_bins = this->n_sources(); + std::vector strengths(n_bins, 0.0); + + mesh_CDF_.resize(n_bins + 1); + mesh_CDF_[0] = {0.0}; + total_strength_ = 0.0; + + // Create cdfs for sampling for an element over a mesh + // Volume scheme is weighted by the volume of each tet + // File scheme is weighted by an array given in the xml file + mesh_strengths_ = std::vector(n_bins, 1.0); + if (check_for_node(node, "strengths")) { + strengths = get_node_array(node, "strengths"); + if (strengths.size() != n_bins) { + fatal_error( + fmt::format("Number of entries in the source strengths array {} does " + "not match the number of entities in mesh {} ({}).", + strengths.size(), mesh_id, n_bins)); + } + mesh_strengths_ = std::move(strengths); + } + + if (get_node_value_bool(node, "volume_normalized")) { + for (int i = 0; i < n_bins; i++) { + mesh_strengths_[i] *= mesh()->volume(i); + } + } + + total_strength_ = + std::accumulate(mesh_strengths_.begin(), mesh_strengths_.end(), 0.0); + + for (int i = 0; i < n_bins; i++) { + mesh_CDF_[i + 1] = mesh_CDF_[i] + mesh_strengths_[i] / total_strength_; + } + + if (fabs(mesh_CDF_.back() - 1.0) > FP_COINCIDENT) { + fatal_error( + fmt::format("Mesh sampling CDF is incorrectly formed. Final value is: {}", + mesh_CDF_.back())); + } + mesh_CDF_.back() = 1.0; +} + +Position MeshSpatial::sample(uint64_t* seed) const +{ + // Create random variable for sampling element from mesh + double eta = prn(seed); + // Sample over the CDF defined in initialization above + int32_t elem_idx = lower_bound_index(mesh_CDF_.begin(), mesh_CDF_.end(), eta); + return mesh()->sample(seed, elem_idx); +} + //============================================================================== // SpatialBox implementation //============================================================================== diff --git a/src/initialize.cpp b/src/initialize.cpp index 1f9e40548..6383af64e 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -410,7 +410,6 @@ void read_separate_xml_files() // Finalize cross sections having assigned temperatures finalize_cross_sections(); - read_tallies_xml(); // Initialize distribcell_filters diff --git a/src/mesh.cpp b/src/mesh.cpp index d9ad0ce8a..121ca1f8a 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -195,6 +195,36 @@ UnstructuredMesh::UnstructuredMesh(pugi::xml_node node) : Mesh(node) } } +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 @@ -236,8 +266,8 @@ void UnstructuredMesh::to_hdf5(hid_t group) const // write element types and connectivity vector volumes; - xt::xtensor connectivity ({static_cast(this->n_bins()), 8}); - xt::xtensor elem_types ({static_cast(this->n_bins()), 1}); + xt::xtensor connectivity({static_cast(this->n_bins()), 8}); + xt::xtensor elem_types({static_cast(this->n_bins()), 1}); for (int i = 0; i < this->n_bins(); i++) { auto conn = this->connectivity(i); @@ -245,17 +275,20 @@ void UnstructuredMesh::to_hdf5(hid_t group) const // write linear tet element if (conn.size() == 4) { - xt::view(elem_types, i, xt::all()) = static_cast(ElementType::LINEAR_TET); - xt::view(connectivity, i, xt::all()) = xt::xarray({conn[0], conn[1], conn[2], conn[3], - -1, -1, -1, -1}); - // write linear hex element + xt::view(elem_types, i, xt::all()) = + static_cast(ElementType::LINEAR_TET); + xt::view(connectivity, i, xt::all()) = + xt::xarray({conn[0], conn[1], conn[2], conn[3], -1, -1, -1, -1}); + // write linear hex element } else if (conn.size() == 8) { - xt::view(elem_types, i, xt::all()) = static_cast(ElementType::LINEAR_HEX); - xt::view(connectivity, i, xt::all()) = xt::xarray({conn[0], conn[1], conn[2], conn[3], - conn[4], conn[5], conn[6], conn[7]}); + xt::view(elem_types, i, xt::all()) = + static_cast(ElementType::LINEAR_HEX); + xt::view(connectivity, i, xt::all()) = xt::xarray({conn[0], conn[1], + conn[2], conn[3], conn[4], conn[5], conn[6], conn[7]}); } else { num_elem_skipped++; - xt::view(elem_types, i, xt::all()) = static_cast(ElementType::UNSUPPORTED); + xt::view(elem_types, i, xt::all()) = + static_cast(ElementType::UNSUPPORTED); xt::view(connectivity, i, xt::all()) = -1; } } @@ -283,6 +316,18 @@ void UnstructuredMesh::set_length_multiplier(double length_multiplier) specified_length_multiplier_ = true; } +ElementType UnstructuredMesh::element_type(int bin) const +{ + auto conn = connectivity(bin); + + if (conn.size() == 4) + return ElementType::LINEAR_TET; + else if (conn.size() == 8) + return ElementType::LINEAR_HEX; + else + return ElementType::UNSUPPORTED; +} + StructuredMesh::MeshIndex StructuredMesh::get_indices( Position r, bool& in_mesh) const { @@ -327,6 +372,16 @@ StructuredMesh::MeshIndex StructuredMesh::get_indices_from_bin(int bin) const return ijk; } +Position StructuredMesh::sample(uint64_t* seed, int32_t bin) const +{ + fatal_error("Position sampling on structured meshes is not yet implemented"); +} + +double StructuredMesh::volume(int bin) const +{ + fatal_error("Unable to get volume of structured mesh, not yet implemented"); +} + int StructuredMesh::get_bin(Position r) const { // Determine indices @@ -1819,7 +1874,6 @@ void MOABMesh::initialize() fatal_error("Failed to get all vertex handles"); } - // make an entity set for all tetrahedra // this is used for convenience later in output rval = mbi_->create_meshset(moab::MESHSET_SET, tetset_); @@ -2049,6 +2103,35 @@ std::string MOABMesh::library() const return mesh_lib_type; } +// Sample position within a tet for MOAB type tets +Position MOABMesh::sample(uint64_t* seed, int32_t bin) const +{ + + moab::EntityHandle tet_ent = get_ent_handle_from_bin(bin); + + // Get vertex coordinates for MOAB tet + const moab::EntityHandle* conn1; + int conn1_size; + moab::ErrorCode rval = mbi_->get_connectivity(tet_ent, conn1, conn1_size); + if (rval != moab::MB_SUCCESS || conn1_size != 4) { + fatal_error(fmt::format( + "Failed to get tet connectivity or connectivity size ({}) is invalid.", + conn1_size)); + } + moab::CartVect p[4]; + rval = mbi_->get_coords(conn1, conn1_size, p[0].array()); + if (rval != moab::MB_SUCCESS) { + 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); +} + double MOABMesh::tet_volume(moab::EntityHandle tet) const { vector conn; @@ -2169,10 +2252,12 @@ std::pair, vector> MOABMesh::plot( return {}; } -int MOABMesh::get_vert_idx_from_handle(moab::EntityHandle vert) const { +int MOABMesh::get_vert_idx_from_handle(moab::EntityHandle vert) const +{ int idx = vert - verts_[0]; if (idx >= n_vertices()) { - fatal_error(fmt::format("Invalid vertex idx {} (# vertices {})", idx, n_vertices())); + fatal_error( + fmt::format("Invalid vertex idx {} (# vertices {})", idx, n_vertices())); } return idx; } @@ -2244,11 +2329,13 @@ Position MOABMesh::centroid(int bin) const return {centroid[0], centroid[1], centroid[2]}; } -int MOABMesh::n_vertices() const { +int MOABMesh::n_vertices() const +{ return verts_.size(); } -Position MOABMesh::vertex(int id) const { +Position MOABMesh::vertex(int id) const +{ moab::ErrorCode rval; @@ -2263,7 +2350,8 @@ Position MOABMesh::vertex(int id) const { return {coords[0], coords[1], coords[2]}; } -std::vector MOABMesh::connectivity(int bin) const { +std::vector MOABMesh::connectivity(int bin) const +{ moab::ErrorCode rval; auto tet = get_ent_handle_from_bin(bin); @@ -2417,14 +2505,15 @@ const std::string LibMesh::mesh_lib_type = "libmesh"; LibMesh::LibMesh(pugi::xml_node node) : UnstructuredMesh(node) { - // filename_ and length_multiplier_ will already be set by the UnstructuredMesh constructor + // filename_ and length_multiplier_ will already be set by the + // UnstructuredMesh constructor set_mesh_pointer_from_filename(filename_); set_length_multiplier(length_multiplier_); initialize(); } // create the mesh from a pointer to a libMesh Mesh -LibMesh::LibMesh(libMesh::MeshBase & input_mesh, double length_multiplier) +LibMesh::LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier) { m_ = &input_mesh; set_length_multiplier(length_multiplier); @@ -2501,6 +2590,20 @@ void LibMesh::initialize() bbox_ = libMesh::MeshTools::create_bounding_box(*m_); } +// Sample position within a tet for LibMesh type tets +Position LibMesh::sample(uint64_t* seed, int32_t bin) const +{ + const auto& elem = get_element_from_bin(bin); + // Get tet vertex coordinates from LibMesh + std::array tet_verts; + for (int i = 0; i < elem.n_nodes(); i++) { + auto node_ref = elem.node_ref(i); + tet_verts[i] = {node_ref(0), node_ref(1), node_ref(2)}; + } + // Samples position within tet using Barycentric coordinates + return this->sample_tet(tet_verts, seed); +} + Position LibMesh::centroid(int bin) const { const auto& elem = this->get_element_from_bin(bin); @@ -2682,7 +2785,7 @@ const libMesh::Elem& LibMesh::get_element_from_bin(int bin) const double LibMesh::volume(int bin) const { - return m_->elem_ref(bin).volume(); + return this->get_element_from_bin(bin).volume(); } #endif // LIBMESH diff --git a/src/settings.cpp b/src/settings.cpp index feff4fbd3..2f8fbd1d0 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -281,6 +281,9 @@ void read_settings_xml(pugi::xml_node root) } } + // Check for user meshes and allocate + read_meshes(root); + // Look for deprecated cross_sections.xml file in settings.xml if (check_for_node(root, "cross_sections")) { warning( @@ -566,9 +569,6 @@ void read_settings_xml(pugi::xml_node root) } } - // Read meshes - read_meshes(root); - // Shannon Entropy mesh if (check_for_node(root, "entropy_mesh")) { int temp = std::stoi(get_node_value(root, "entropy_mesh")); diff --git a/src/source.cpp b/src/source.cpp index 3a80e1156..5862a41e2 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -94,6 +94,8 @@ IndependentSource::IndependentSource(pugi::xml_node node) space_ = UPtrSpace {new CylindricalIndependent(node_space)}; } else if (type == "spherical") { space_ = UPtrSpace {new SphericalIndependent(node_space)}; + } else if (type == "mesh") { + space_ = UPtrSpace {new MeshSpatial(node_space)}; } else if (type == "box") { space_ = UPtrSpace {new SpatialBox(node_space)}; } else if (type == "fission") { diff --git a/tests/unit_tests/test_mesh_tets.e b/tests/unit_tests/test_mesh_tets.e new file mode 120000 index 000000000..8a6287b4d --- /dev/null +++ b/tests/unit_tests/test_mesh_tets.e @@ -0,0 +1 @@ +../regression_tests/unstructured_mesh/test_mesh_tets.e \ No newline at end of file diff --git a/tests/unit_tests/test_source_mesh.py b/tests/unit_tests/test_source_mesh.py new file mode 100644 index 000000000..596b22f9a --- /dev/null +++ b/tests/unit_tests/test_source_mesh.py @@ -0,0 +1,203 @@ +from itertools import product +from pathlib import Path +from subprocess import call + +import pytest +import numpy as np +import openmc +import openmc.lib + +from tests import cdtemp +from tests.regression_tests import config + + +TETS_PER_VOXEL = 12 + +# This test uses a geometry file with cells that match a regular mesh. Each cell +# in the geometry corresponds to 12 tetrahedra in the unstructured mesh file. +@pytest.fixture +def model(): + openmc.reset_auto_ids() + + ### Materials ### + materials = openmc.Materials() + + 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) + + ### Geometry ### + # This test uses a geometry file that resembles a regular mesh. + # 12 tets are used to match each voxel in the geometry. + + # create a regular mesh that matches the superimposed mesh + regular_mesh = openmc.RegularMesh(mesh_id=10) + regular_mesh.lower_left = (-10, -10, -10) + regular_mesh.dimension = (10, 10, 10) + regular_mesh.width = (2, 2, 2) + + root_cell, _ = regular_mesh.build_cells(bc=['vacuum']*6) + + geometry = openmc.Geometry(root=[root_cell]) + + ### Settings ### + settings = openmc.Settings() + settings.run_mode = 'fixed source' + settings.particles = 100 + settings.batches = 2 + + return openmc.model.Model(geometry=geometry, + materials=materials, + settings=settings) + +### Setup test cases ### +param_values = (['libmesh', 'moab'], # mesh libraries + ['uniform', 'manual']) # Element weighting schemes + +test_cases = [] +for i, (lib, schemes) in enumerate(product(*param_values)): + test_cases.append({'library' : lib, + 'source_strengths' : schemes}) + +def ids(params): + """Test naming function for clarity""" + return f"{params['library']}-{params['source_strengths']}" + +@pytest.mark.parametrize("test_cases", test_cases, ids=ids) +def test_unstructured_mesh_sampling(model, request, test_cases): + # skip the test if the library is not enabled + if test_cases['library'] == 'moab' and not openmc.lib._dagmc_enabled(): + pytest.skip("DAGMC (and MOAB) mesh not enabled in this build.") + + if test_cases['library'] == 'libmesh' and not openmc.lib._libmesh_enabled(): + pytest.skip("LibMesh is not enabled in this build.") + + # setup mesh source ### + mesh_filename = Path(request.fspath).parent / "test_mesh_tets.e" + uscd_mesh = openmc.UnstructuredMesh(mesh_filename, test_cases['library']) + + # subtract one to account for root cell produced by RegularMesh.build_cells + n_cells = len(model.geometry.get_all_cells()) - 1 + + # set source weights according to test case + if test_cases['source_strengths'] == 'uniform': + vol_norm = True + strengths = None + elif test_cases['source_strengths'] == 'manual': + vol_norm = False + # assign random weights + strengths = np.random.rand(n_cells*TETS_PER_VOXEL) + + # create the spatial distribution based on the mesh + space = openmc.stats.MeshSpatial(uscd_mesh, strengths, vol_norm) + + energy = openmc.stats.Discrete(x=[15.e+06], p=[1.0]) + source = openmc.Source(space=space, energy=energy) + model.settings.source = source + + with cdtemp([mesh_filename]): + model.export_to_xml() + + n_measurements = 100 + n_samples = 1000 + + cell_counts = np.zeros((n_cells, n_measurements)) + + # This model contains 1000 geometry cells. Each cell is a hex + # corresponding to 12 of the tets. This test runs 1000 samples. This + # results in the following average for each cell + openmc.lib.init([]) + + # perform many sets of samples and track counts for each cell + for m in range(n_measurements): + sites = openmc.lib.sample_external_source(n_samples) + cells = [openmc.lib.find_cell(s.r) for s in sites] + + for c in cells: + # subtract one from index to account for root cell + cell_counts[c[0]._index - 1, m] += 1 + + # make sure particle transport is successful + openmc.lib.run() + openmc.lib.finalize() + + # normalize cell counts to get sampling frequency per particle + cell_counts /= n_samples + + # get the mean and std. dev. of the cell counts + mean = cell_counts.mean(axis=1) + std_dev = cell_counts.std(axis=1) + + if test_cases['source_strengths'] == 'uniform': + exp_vals = np.ones(n_cells) / n_cells + else: + # sum up the source strengths for each tet, these are the expected true mean + # of the sampling frequency for that cell + exp_vals = strengths.reshape(-1, 12).sum(axis=1) / sum(strengths) + + diff = np.abs(mean - exp_vals) + assert((diff < 2*std_dev).sum() / diff.size >= 0.95) + assert((diff < 6*std_dev).sum() / diff.size >= 0.997) + + +def test_strengths_size_failure(request, model): + # setup mesh source ### + mesh_filename = Path(request.fspath).parent / "test_mesh_tets.e" + uscd_mesh = openmc.UnstructuredMesh(mesh_filename, 'libmesh') + + # intentionally incorrectly sized to trigger an error + n_cells = len(model.geometry.get_all_cells()) + strengths = np.random.rand(n_cells*TETS_PER_VOXEL) + + # create the spatial distribution based on the mesh + space = openmc.stats.MeshSpatial(uscd_mesh, strengths) + + energy = openmc.stats.Discrete(x=[15.e+06], p=[1.0]) + source = openmc.Source(space=space, energy=energy) + model.settings.source = source + + # skip the test if unstructured mesh is not available + if not openmc.lib._libmesh_enabled(): + if openmc.lib._dagmc_enabled(): + source.space.mesh.library = 'moab' + else: + pytest.skip("Unstructured mesh support unavailable.") + + # make sure that an incorrrectly sized strengths array causes a failure + source.space.strengths = source.space.strengths[:-1] + + mesh_filename = Path(request.fspath).parent / source.space.mesh.filename + + with pytest.raises(RuntimeError, match=r'strengths array'), cdtemp([mesh_filename]): + model.export_to_xml() + openmc.run() + +def test_roundtrip(run_in_tmpdir, model, request): + if not openmc.lib._libmesh_enabled() and not openmc.lib._dagmc_enabled(): + pytest.skip("Unstructured mesh is not enabled in this build.") + + mesh_filename = Path(request.fspath).parent / 'test_mesh_tets.e' + ucd_mesh = openmc.UnstructuredMesh(mesh_filename, library='libmesh') + + if not openmc.lib._libmesh_enabled(): + ucd_mesh.library = 'moab' + + n_cells = len(model.geometry.get_all_cells()) + + space_out = openmc.MeshSpatial(ucd_mesh) + space_out.strengths = np.random.rand(n_cells*TETS_PER_VOXEL) + model.settings.source = openmc.Source(space=space_out) + + # write out the model + model.export_to_xml() + + model_in = openmc.Model.from_xml() + + space_in = model_in.settings.source[0].space + + np.testing.assert_equal(space_out.strengths, space_in.strengths) + + assert space_in.mesh.id == space_out.mesh.id + assert space_in.volume_normalized == space_out.volume_normalized diff --git a/tests/unit_tests/weightwindows/test.py b/tests/unit_tests/weightwindows/test.py index d73a52f45..1dfad56aa 100644 --- a/tests/unit_tests/weightwindows/test.py +++ b/tests/unit_tests/weightwindows/test.py @@ -225,7 +225,7 @@ def test_lower_ww_bounds_shape(): assert ww.lower_ww_bounds.shape == (2, 3, 4, 1) -def test_roundtrip(model, wws): +def test_roundtrip(run_in_tmpdir, model, wws): model.settings.weight_windows = wws # write the model with weight windows to XML