From 05ea1d8a03a5b89fe33a8db269b7fb0a6191986e Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Sat, 2 Dec 2023 11:35:23 -0600 Subject: [PATCH] Mesh Source Class (#2759) Co-authored-by: Paul Romano Co-authored-by: Jonathan Shimwell --- .github/workflows/format-check.yml | 9 +- docs/source/io_formats/settings.rst | 13 +- docs/source/pythonapi/base.rst | 1 + include/openmc/distribution.h | 7 +- include/openmc/distribution_multi.h | 2 + include/openmc/distribution_spatial.h | 15 ++- include/openmc/mesh.h | 75 +++++++++--- include/openmc/source.h | 42 ++++++- openmc/bounding_box.py | 4 + openmc/settings.py | 20 ++- openmc/source.py | 167 ++++++++++++++++++++++++-- openmc/stats/multivariate.py | 1 + src/distribution.cpp | 12 +- src/distribution_multi.cpp | 19 +++ src/distribution_spatial.cpp | 88 +++++++++++--- src/mesh.cpp | 78 ++++++++++-- src/settings.cpp | 23 +--- src/source.cpp | 155 ++++++++++++++++-------- tests/unit_tests/test_source_mesh.py | 167 +++++++++++++++++++++++++- 19 files changed, 752 insertions(+), 146 deletions(-) diff --git a/.github/workflows/format-check.yml b/.github/workflows/format-check.yml index 4914d5e5ea9..36ccece3cc5 100644 --- a/.github/workflows/format-check.yml +++ b/.github/workflows/format-check.yml @@ -1,6 +1,13 @@ name: C++ Format Check -on: pull_request +on: + # allow workflow to be run manually + workflow_dispatch: + + pull_request: + branches: + - develop + - master jobs: cpp-linter: diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 03abe9a3c4f..9b855b25029 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -473,6 +473,8 @@ pseudo-random number generator. *Default*: 1 +.. _source_element: + -------------------- ```` Element -------------------- @@ -491,7 +493,8 @@ attributes/sub-elements: *Default*: 1.0 :type: - Indicator of source type. One of ``independent``, ``file``, or ``compiled``. + Indicator of source type. One of ``independent``, ``file``, ``compiled``, or ``mesh``. + The type of the source will be determined by this attribute if it is present. :particle: The source particle type, either ``neutron`` or ``photon``. @@ -664,6 +667,14 @@ attributes/sub-elements: *Default*: false + :mesh: + For mesh sources, this indicates the ID of the corresponding mesh. + + :source: + For mesh sources, this sub-element specifies the source for an individual + mesh element and follows the format for :ref:`source_element`. The number of + ```` sub-elements should correspond to the number of mesh elements. + .. _univariate: Univariate Probability Distributions diff --git a/docs/source/pythonapi/base.rst b/docs/source/pythonapi/base.rst index ba86db117bb..705aecee9f1 100644 --- a/docs/source/pythonapi/base.rst +++ b/docs/source/pythonapi/base.rst @@ -25,6 +25,7 @@ Simulation Settings openmc.IndependentSource openmc.FileSource openmc.CompiledSource + openmc.MeshSource openmc.SourceParticle openmc.VolumeCalculation openmc.Settings diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index 65c68ad6457..e70c803de2f 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -7,6 +7,7 @@ #include // for size_t #include "pugixml.hpp" +#include #include "openmc/constants.h" #include "openmc/memory.h" // for unique_ptr @@ -43,9 +44,9 @@ class DiscreteIndex { public: DiscreteIndex() {}; DiscreteIndex(pugi::xml_node node); - DiscreteIndex(const double* p, int n); + DiscreteIndex(gsl::span p); - void assign(const double* p, int n); + void assign(gsl::span p); //! Sample a value from the distribution //! \param seed Pseudorandom number seed pointer @@ -77,7 +78,7 @@ class DiscreteIndex { class Discrete : public Distribution { public: explicit Discrete(pugi::xml_node node); - Discrete(const double* x, const double* p, int n); + Discrete(const double* x, const double* p, size_t n); //! Sample a value from the distribution //! \param seed Pseudorandom number seed pointer diff --git a/include/openmc/distribution_multi.h b/include/openmc/distribution_multi.h index e171d58ab61..9e84d03d57f 100644 --- a/include/openmc/distribution_multi.h +++ b/include/openmc/distribution_multi.h @@ -22,6 +22,8 @@ class UnitSphereDistribution { explicit UnitSphereDistribution(pugi::xml_node node); virtual ~UnitSphereDistribution() = default; + static unique_ptr create(pugi::xml_node node); + //! Sample a direction from the distribution //! \param seed Pseudorandom number seed pointer //! \return Direction sampled diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index 9fba9506025..bd10a299605 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -19,6 +19,8 @@ class SpatialDistribution { //! Sample a position from the distribution virtual Position sample(uint64_t* seed) const = 0; + + static unique_ptr create(pugi::xml_node node); }; //============================================================================== @@ -102,16 +104,27 @@ class SphericalIndependent : public SpatialDistribution { class MeshSpatial : public SpatialDistribution { public: explicit MeshSpatial(pugi::xml_node node); + explicit MeshSpatial(int32_t mesh_id, gsl::span strengths); //! Sample a position from the distribution //! \param seed Pseudorandom number seed pointer //! \return Sampled position Position sample(uint64_t* seed) const override; - const Mesh* mesh() const { return model::meshes.at(mesh_idx_).get(); } + //! Sample the mesh for an element and position within that element + //! \param seed Pseudorandom number seed pointer + //! \return Sampled element index and position within that element + std::pair sample_mesh(uint64_t* seed) const; + //! For unstructured meshes, ensure that elements are all linear tetrahedra + void check_element_types() const; + + // Accessors + const Mesh* mesh() const { return model::meshes.at(mesh_idx_).get(); } int32_t n_sources() const { return this->mesh()->n_bins(); } + double total_strength() { return this->elem_idx_dist_.integral(); } + private: int32_t mesh_idx_ {C_NONE}; DiscreteIndex elem_idx_dist_; //!< Distribution of diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 47f54c023fc..de556321f22 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -10,6 +10,7 @@ #include "pugixml.hpp" #include "xtensor/xtensor.hpp" +#include "openmc/error.h" #include "openmc/memory.h" // for unique_ptr #include "openmc/particle.h" #include "openmc/position.h" @@ -81,12 +82,12 @@ class Mesh { //! Return a position in the local coordinates of the mesh virtual Position local_coords(const Position& r) const { return r; }; - //! Sample a mesh volume using a certain seed + //! Sample a position within a mesh element // - //! \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; + //! \param[in] bin Bin value of the mesh element sampled + //! \param[inout] seed Seed to use for random sampling + //! \return sampled position within mesh element + virtual Position sample_element(int32_t bin, uint64_t* seed) const = 0; //! Determine which bins were crossed by a particle // @@ -183,7 +184,12 @@ class StructuredMesh : public Mesh { } }; - Position sample(uint64_t* seed, int32_t bin) const override; + Position sample_element(int32_t bin, uint64_t* seed) const override + { + return sample_element(get_indices_from_bin(bin), seed); + }; + + virtual Position sample_element(const MeshIndex& ijk, uint64_t* seed) const; int get_bin(Position r) const override; @@ -240,6 +246,30 @@ class StructuredMesh : public Mesh { //! \param[in] i Direction index virtual int get_index_in_direction(double r, int i) const = 0; + //! Get the coordinate for the mesh grid boundary in the positive direction + //! + //! \param[in] ijk Array of mesh indices + //! \param[in] i Direction index + virtual double positive_grid_boundary(const MeshIndex& ijk, int i) const + { + auto msg = + fmt::format("Attempting to call positive_grid_boundary on a {} mesh.", + get_mesh_type()); + fatal_error(msg); + }; + + //! Get the coordinate for the mesh grid boundary in the negative direction + //! + //! \param[in] ijk Array of mesh indices + //! \param[in] i Direction index + virtual double negative_grid_boundary(const MeshIndex& ijk, int i) const + { + auto msg = + fmt::format("Attempting to call negative_grid_boundary on a {} mesh.", + get_mesh_type()); + fatal_error(msg); + }; + //! Get the closest distance from the coordinate r to the grid surface //! in i direction that bounds mesh cell ijk and that is larger than l //! The coordinate r does not have to be inside the mesh cell ijk. In @@ -322,18 +352,17 @@ class RegularMesh : public StructuredMesh { void to_hdf5(hid_t group) const override; - // New methods //! Get the coordinate for the mesh grid boundary in the positive direction //! //! \param[in] ijk Array of mesh indices //! \param[in] i Direction index - double positive_grid_boundary(const MeshIndex& ijk, int i) const; + double positive_grid_boundary(const MeshIndex& ijk, int i) const override; //! Get the coordinate for the mesh grid boundary in the negative direction //! //! \param[in] ijk Array of mesh indices //! \param[in] i Direction index - double negative_grid_boundary(const MeshIndex& ijk, int i) const; + double negative_grid_boundary(const MeshIndex& ijk, int i) const override; //! Count number of bank sites in each mesh bin / energy bin // @@ -373,18 +402,17 @@ class RectilinearMesh : public StructuredMesh { void to_hdf5(hid_t group) const override; - // New methods //! Get the coordinate for the mesh grid boundary in the positive direction //! //! \param[in] ijk Array of mesh indices //! \param[in] i Direction index - double positive_grid_boundary(const MeshIndex& ijk, int i) const; + double positive_grid_boundary(const MeshIndex& ijk, int i) const override; //! Get the coordinate for the mesh grid boundary in the negative direction //! //! \param[in] ijk Array of mesh indices //! \param[in] i Direction index - double negative_grid_boundary(const MeshIndex& ijk, int i) const; + double negative_grid_boundary(const MeshIndex& ijk, int i) const override; //! Return the volume for a given mesh index double volume(const MeshIndex& ijk) const override; @@ -410,6 +438,8 @@ class CylindricalMesh : public PeriodicStructuredMesh { static const std::string mesh_type; + Position sample_element(const MeshIndex& ijk, uint64_t* seed) const override; + MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override; @@ -420,10 +450,16 @@ class CylindricalMesh : public PeriodicStructuredMesh { double volume(const MeshIndex& ijk) const override; - array, 3> grid_; + // grid accessors + double r(int i) const { return grid_[0][i]; } + double phi(int i) const { return grid_[1][i]; } + double z(int i) const { return grid_[2][i]; } int set_grid(); + // Data members + array, 3> grid_; + private: double find_r_crossing( const Position& r, const Direction& u, double l, int shell) const; @@ -466,6 +502,8 @@ class SphericalMesh : public PeriodicStructuredMesh { static const std::string mesh_type; + Position sample_element(const MeshIndex& ijk, uint64_t* seed) const override; + MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const override; @@ -474,10 +512,15 @@ class SphericalMesh : public PeriodicStructuredMesh { void to_hdf5(hid_t group) const override; - array, 3> grid_; + double r(int i) const { return grid_[0][i]; } + double theta(int i) const { return grid_[1][i]; } + double phi(int i) const { return grid_[2][i]; } int set_grid(); + // Data members + array, 3> grid_; + private: double find_r_crossing( const Position& r, const Direction& u, double l, int shell) const; @@ -621,7 +664,7 @@ class MOABMesh : public UnstructuredMesh { // Overridden Methods - Position sample(uint64_t* seed, int32_t bin) const 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; @@ -789,7 +832,7 @@ class LibMesh : public UnstructuredMesh { 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; + Position sample_element(int32_t bin, uint64_t* seed) const override; int get_bin(Position r) const override; diff --git a/include/openmc/source.h b/include/openmc/source.h index ad2aae562bb..d5ea9bd1d56 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -50,6 +50,8 @@ class Source { // Methods that can be overridden virtual double strength() const { return 1.0; } + + static unique_ptr create(pugi::xml_node node); }; //============================================================================== @@ -101,12 +103,13 @@ class IndependentSource : public Source { class FileSource : public Source { public: // Constructors - explicit FileSource(std::string path); - explicit FileSource(const vector& sites) : sites_ {sites} {} + explicit FileSource(pugi::xml_node node); + explicit FileSource(const std::string& path); // Methods SourceSite sample(uint64_t* seed) const override; - + void load_sites_from_file( + const std::string& path); //!< Load source sites from file private: vector sites_; //!< Source sites from a file }; @@ -118,7 +121,7 @@ class FileSource : public Source { class CompiledSourceWrapper : public Source { public: // Constructors, destructors - CompiledSourceWrapper(std::string path, std::string parameters); + CompiledSourceWrapper(pugi::xml_node node); ~CompiledSourceWrapper(); // Defer implementation to custom source library @@ -129,6 +132,8 @@ class CompiledSourceWrapper : public Source { double strength() const override { return compiled_source_->strength(); } + void setup(const std::string& path, const std::string& parameters); + private: void* shared_library_; //!< library from dlopen unique_ptr compiled_source_; @@ -136,6 +141,35 @@ class CompiledSourceWrapper : public Source { typedef unique_ptr create_compiled_source_t(std::string parameters); +//============================================================================== +//! Mesh-based source with different distributions for each element +//============================================================================== + +class MeshSource : public Source { +public: + // Constructors + explicit MeshSource(pugi::xml_node node); + + //! Sample from the external source distribution + //! \param[inout] seed Pseudorandom seed pointer + //! \return Sampled site + SourceSite sample(uint64_t* seed) const override; + + // Properties + double strength() const override { return space_->total_strength(); } + + // Accessors + const std::unique_ptr& source(int32_t i) const + { + return sources_.size() == 1 ? sources_[0] : sources_[i]; + } + +private: + // Data members + unique_ptr space_; //!< Mesh spatial + vector> sources_; //!< Source distributions +}; + //============================================================================== // Functions //============================================================================== diff --git a/openmc/bounding_box.py b/openmc/bounding_box.py index b74931d73ae..e5655cf8cae 100644 --- a/openmc/bounding_box.py +++ b/openmc/bounding_box.py @@ -95,6 +95,10 @@ def __or__(self, other: BoundingBox) -> BoundingBox: new |= other return new + def __contains__(self, point): + """Check whether or not a point is in the bounding box""" + return all(point > self.lower_left) and all(point < self.upper_right) + @property def center(self) -> np.ndarray: return (self[0] + self[1]) / 2 diff --git a/openmc/settings.py b/openmc/settings.py index d8f74dd1179..5d0800ac13d 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -12,7 +12,7 @@ import openmc.checkvalue as cv from openmc.stats.multivariate import MeshSpatial -from . import (RegularMesh, SourceBase, IndependentSource, +from . import (RegularMesh, SourceBase, MeshSource, IndependentSource, VolumeCalculation, WeightWindows, WeightWindowGenerator) from ._xml import clean_indentation, get_text, reorder_attributes from openmc.checkvalue import PathLike @@ -1072,13 +1072,19 @@ def _create_max_order_subelement(self, root): element = ET.SubElement(root, "max_order") element.text = str(self._max_order) - def _create_source_subelement(self, root): + def _create_source_subelement(self, root, mesh_memo=None): for source in self.source: root.append(source.to_xml_element()) if isinstance(source, IndependentSource) and 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()) + if isinstance(source, MeshSource): + path = f"./mesh[@id='{source.mesh.id}']" + if root.find(path) is None: + root.append(source.mesh.to_xml_element()) + if mesh_memo is not None: + mesh_memo.add(source.mesh.id) def _create_volume_calcs_subelement(self, root): for calc in self.volume_calculations: @@ -1365,7 +1371,8 @@ def _create_weight_windows_subelement(self, root, mesh_memo=None): path = f"./mesh[@id='{ww.mesh.id}']" if root.find(path) is None: root.append(ww.mesh.to_xml_element()) - if mesh_memo is not None: mesh_memo.add(ww.mesh.id) + if mesh_memo is not None: + mesh_memo.add(ww.mesh.id) if self._weight_windows_on is not None: elem = ET.SubElement(root, "weight_windows_on") @@ -1787,7 +1794,7 @@ def to_xml_element(self, mesh_memo=None): self._create_max_write_lost_particles_subelement(element) self._create_generations_per_batch_subelement(element) self._create_keff_trigger_subelement(element) - self._create_source_subelement(element) + self._create_source_subelement(element, mesh_memo) self._create_output_subelement(element) self._create_statepoint_subelement(element) self._create_sourcepoint_subelement(element) @@ -1874,6 +1881,11 @@ def from_xml_element(cls, elem, meshes=None): Settings object """ + # read all meshes under the settings node and update + settings_meshes = _read_meshes(elem) + meshes = {} if meshes is None else meshes + meshes.update(settings_meshes) + settings = cls() settings._eigenvalue_from_xml_element(elem) settings._run_mode_from_xml_element(elem) diff --git a/openmc/source.py b/openmc/source.py index 98b7673b176..c339d6ed806 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -18,7 +18,7 @@ from openmc.stats.multivariate import UnitSphere, Spatial from openmc.stats.univariate import Univariate from ._xml import get_text -from .mesh import MeshBase +from .mesh import MeshBase, StructuredMesh, UnstructuredMesh class SourceBase(ABC): @@ -31,7 +31,7 @@ class SourceBase(ABC): Attributes ---------- - type : {'independent', 'file', 'compiled'} + type : {'independent', 'file', 'compiled', 'mesh'} Indicator of source type. strength : float Strength of the source @@ -61,7 +61,6 @@ def populate_xml_element(self, element): XML element containing source data """ - pass def to_xml_element(self) -> ET.Element: """Return XML representation of the source @@ -79,7 +78,7 @@ def to_xml_element(self) -> ET.Element: return element @classmethod - def from_xml_element(cls, elem: ET.Element, meshes=None) -> openmc.SourceBase: + def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase: """Generate source from an XML element Parameters @@ -114,6 +113,8 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> openmc.SourceBase: return CompiledSource.from_xml_element(elem) elif source_type == 'file': return FileSource.from_xml_element(elem) + elif source_type == 'mesh': + return MeshSource.from_xml_element(elem, meshes) else: raise ValueError(f'Source type {source_type} is not recognized') @@ -292,13 +293,13 @@ def domain_type(self, domain_type): def populate_xml_element(self, element): """Add necessary source information to an XML element + Returns ------- element : lxml.etree._Element XML element containing source data """ - super().populate_xml_element(element) element.set("particle", self.particle) if self.space is not None: element.append(self.space.to_xml_element()) @@ -315,7 +316,7 @@ def populate_xml_element(self, element): id_elem.text = ' '.join(str(uid) for uid in self.domain_ids) @classmethod - def from_xml_element(cls, elem: ET.Element, meshes=None) -> 'openmc.SourceBase': + def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase: """Generate source from an XML element Parameters @@ -382,6 +383,154 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> 'openmc.SourceBase': return source +class MeshSource(SourceBase): + """A source with a spatial distribution over mesh elements + + This class represents a mesh-based source in which random positions are + uniformly sampled within mesh elements and each element can have independent + angle, energy, and time distributions. The element sampled is chosen based + on the relative strengths of the sources applied to the elements. The + strength of the mesh source as a whole is the sum of all source strengths + applied to the elements. + + .. versionadded:: 0.14.1 + + Parameters + ---------- + mesh : openmc.MeshBase + The mesh over which source sites will be generated. + sources : iterable of openmc.SourceBase + Sources for each element in the mesh. If spatial distributions are set + on any of the source objects, they will be ignored during source site + sampling. + + Attributes + ---------- + mesh : openmc.MeshBase + The mesh over which source sites will be generated. + sources : numpy.ndarray or iterable of openmc.SourceBase + The set of sources to apply to each element. The shape of this array + must match the shape of the mesh with and exception in the case of + unstructured mesh, which allows for application of 1-D array or + iterable. + strength : float + Strength of the source + type : str + Indicator of source type: 'mesh' + + """ + def __init__(self, mesh: MeshBase, sources: Sequence[SourceBase]): + self.mesh = mesh + self.sources = sources + + @property + def type(self) -> str: + return "mesh" + + @property + def mesh(self) -> MeshBase: + return self._mesh + + @property + def strength(self) -> float: + return sum(s.strength for s in self.sources.flat) + + @property + def sources(self) -> np.ndarray: + return self._sources + + @mesh.setter + def mesh(self, m): + cv.check_type('source mesh', m, MeshBase) + self._mesh = m + + @sources.setter + def sources(self, s): + cv.check_iterable_type('mesh sources', s, SourceBase, max_depth=3) + + s = np.asarray(s) + + if isinstance(self.mesh, StructuredMesh) and s.shape != self.mesh.dimension: + raise ValueError('The shape of the source array' + f'({s.shape}) does not match the ' + f'dimensions of the structured mesh ({self.mesh.dimension})') + elif isinstance(self.mesh, UnstructuredMesh): + if len(s.shape) > 1: + raise ValueError('Sources must be a 1-D array for unstructured mesh') + + self._sources = s + for src in self._sources.flat: + if isinstance(src, IndependentSource) and src.space is not None: + warnings.warn('Some sources on the mesh have spatial ' + 'distributions that will be ignored at runtime.') + break + + @strength.setter + def strength(self, val): + cv.check_type('mesh source strength', val, Real) + self.set_total_strength(val) + + def set_total_strength(self, strength: float): + """Scales the element source strengths based on a desired total strength. + + Parameters + ---------- + strength : float + Total source strength + + """ + current_strength = self.strength if self.strength != 0.0 else 1.0 + + for s in self.sources.flat: + s.strength *= strength / current_strength + + def normalize_source_strengths(self): + """Update all element source strengths such that they sum to 1.0.""" + self.set_total_strength(1.0) + + def populate_xml_element(self, elem: ET.Element): + """Add necessary source information to an XML element + + Returns + ------- + element : lxml.etree._Element + XML element containing source data + + """ + elem.set("mesh", str(self.mesh.id)) + + # write in the order of mesh indices + for idx in self.mesh.indices: + idx = tuple(i - 1 for i in idx) + elem.append(self.sources[idx].to_xml_element()) + + @classmethod + def from_xml_element(cls, elem: ET.Element, meshes) -> openmc.MeshSource: + """ + Generate MeshSource from an XML element + + Parameters + ---------- + elem : lxml.etree._Element + XML element + meshes : dict + A dictionary with mesh IDs as keys and openmc.MeshBase instances as + values + + Returns + ------- + openmc.MeshSource + MeshSource generated from the XML element + """ + mesh_id = int(get_text(elem, 'mesh')) + + mesh = meshes[mesh_id] + + sources = [SourceBase.from_xml_element(e) for e in elem.iterchildren('source')] + sources = np.asarray(sources).reshape(mesh.dimension, order='F') + return cls(mesh, sources) + + def Source(*args, **kwargs): """ A function for backward compatibility of sources. Will be removed in the @@ -459,15 +608,13 @@ def populate_xml_element(self, element): XML element containing source data """ - super().populate_xml_element(element) - element.set("library", self.library) if self.parameters is not None: element.set("parameters", self.parameters) @classmethod - def from_xml_element(cls, elem: ET.Element, meshes=None) -> openmc.CompiledSource: + def from_xml_element(cls, elem: ET.Element) -> openmc.CompiledSource: """Generate a compiled source from an XML element Parameters @@ -551,8 +698,6 @@ def populate_xml_element(self, element): XML element containing source data """ - super().populate_xml_element(element) - if self.path is not None: element.set("file", self.path) diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 09a2f8752ff..22c3d7ea86b 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -710,6 +710,7 @@ def to_xml_element(self): """ element = ET.Element('space') + element.set('type', 'mesh') element.set("mesh_id", str(self.mesh.id)) element.set("volume_normalized", str(self.volume_normalized)) diff --git a/src/distribution.cpp b/src/distribution.cpp index 9c76147ee27..3026630b335 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -27,17 +27,17 @@ DiscreteIndex::DiscreteIndex(pugi::xml_node node) auto params = get_node_array(node, "parameters"); std::size_t n = params.size() / 2; - assign(params.data() + n, n); + assign({params.data() + n, n}); } -DiscreteIndex::DiscreteIndex(const double* p, int n) +DiscreteIndex::DiscreteIndex(gsl::span p) { - assign(p, n); + assign(p); } -void DiscreteIndex::assign(const double* p, int n) +void DiscreteIndex::assign(gsl::span p) { - prob_.assign(p, p + n); + prob_.assign(p.begin(), p.end()); this->init_alias(); } @@ -126,7 +126,7 @@ Discrete::Discrete(pugi::xml_node node) : di_(node) x_.assign(params.begin(), params.begin() + n); } -Discrete::Discrete(const double* x, const double* p, int n) : di_(p, n) +Discrete::Discrete(const double* x, const double* p, size_t n) : di_({p, n}) { x_.assign(x, x + n); diff --git a/src/distribution_multi.cpp b/src/distribution_multi.cpp index 0735f0994ed..b7b3efe5268 100644 --- a/src/distribution_multi.cpp +++ b/src/distribution_multi.cpp @@ -12,6 +12,25 @@ namespace openmc { +unique_ptr UnitSphereDistribution::create( + pugi::xml_node node) +{ + // Check for type of angular distribution + std::string type; + if (check_for_node(node, "type")) + type = get_node_value(node, "type", true, true); + if (type == "isotropic") { + return UPtrAngle {new Isotropic()}; + } else if (type == "monodirectional") { + return UPtrAngle {new Monodirectional(node)}; + } else if (type == "mu-phi") { + return UPtrAngle {new PolarAzimuthal(node)}; + } else { + fatal_error(fmt::format( + "Invalid angular distribution for external source: {}", type)); + } +} + //============================================================================== // UnitSphereDistribution implementation //============================================================================== diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index e1a4c20c21c..8f34ea6b936 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -8,6 +8,36 @@ namespace openmc { +//============================================================================== +// SpatialDistribution implementation +//============================================================================== + +unique_ptr SpatialDistribution::create(pugi::xml_node node) +{ + // Check for type of spatial distribution and read + std::string type; + if (check_for_node(node, "type")) + type = get_node_value(node, "type", true, true); + if (type == "cartesian") { + return UPtrSpace {new CartesianIndependent(node)}; + } else if (type == "cylindrical") { + return UPtrSpace {new CylindricalIndependent(node)}; + } else if (type == "spherical") { + return UPtrSpace {new SphericalIndependent(node)}; + } else if (type == "mesh") { + return UPtrSpace {new MeshSpatial(node)}; + } else if (type == "box") { + return UPtrSpace {new SpatialBox(node)}; + } else if (type == "fission") { + return UPtrSpace {new SpatialBox(node, true)}; + } else if (type == "point") { + return UPtrSpace {new SpatialPoint(node)}; + } else { + fatal_error(fmt::format( + "Invalid spatial distribution for external source: {}", type)); + } +} + //============================================================================== // CartesianIndependent implementation //============================================================================== @@ -189,27 +219,23 @@ Position SphericalIndependent::sample(uint64_t* seed) const MeshSpatial::MeshSpatial(pugi::xml_node node) { + + if (get_node_value(node, "type", true, true) != "mesh") { + fatal_error(fmt::format( + "Incorrect spatial type '{}' for a MeshSpatial distribution")); + } + // 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."); - } + const auto mesh_ptr = model::meshes.at(mesh_idx_).get(); - // 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."); - } - } + check_element_types(); - int32_t n_bins = this->n_sources(); + size_t n_bins = this->n_sources(); std::vector strengths(n_bins, 1.0); // Create cdfs for sampling for an element over a mesh @@ -227,18 +253,44 @@ MeshSpatial::MeshSpatial(pugi::xml_node node) if (get_node_value_bool(node, "volume_normalized")) { for (int i = 0; i < n_bins; i++) { - strengths[i] *= mesh()->volume(i); + strengths[i] *= this->mesh()->volume(i); } } - elem_idx_dist_.assign(strengths.data(), n_bins); + elem_idx_dist_.assign(strengths); } -Position MeshSpatial::sample(uint64_t* seed) const +MeshSpatial::MeshSpatial(int32_t mesh_idx, gsl::span strengths) + : mesh_idx_(mesh_idx) +{ + check_element_types(); + elem_idx_dist_.assign(strengths); +} + +void MeshSpatial::check_element_types() const +{ + const auto umesh_ptr = dynamic_cast(this->mesh()); + if (umesh_ptr) { + // ensure that the unstructured mesh contains only linear tets + for (int bin = 0; bin < umesh_ptr->n_bins(); bin++) { + if (umesh_ptr->element_type(bin) != ElementType::LINEAR_TET) { + fatal_error( + "Mesh specified for source must contain only linear tetrahedra."); + } + } + } +} + +std::pair MeshSpatial::sample_mesh(uint64_t* seed) const { - // Sample over the CDF defined in initialization above + // Sample the CDF defined in initialization above int32_t elem_idx = elem_idx_dist_.sample(seed); - return mesh()->sample(seed, elem_idx); + return {elem_idx, mesh()->sample_element(elem_idx, seed)}; +} + +Position MeshSpatial::sample(uint64_t* seed) const +{ + return this->sample_mesh(seed).second; } //============================================================================== diff --git a/src/mesh.cpp b/src/mesh.cpp index 12052bed62f..ea66552de40 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -27,6 +27,7 @@ #include "openmc/hdf5_interface.h" #include "openmc/memory.h" #include "openmc/message_passing.h" +#include "openmc/random_dist.h" #include "openmc/search.h" #include "openmc/settings.h" #include "openmc/tallies/filter.h" @@ -165,6 +166,23 @@ xt::xtensor StructuredMesh::get_x_shape() const return xt::adapt(tmp_shape, {n_dimension_}); } +Position StructuredMesh::sample_element( + const MeshIndex& ijk, uint64_t* seed) const +{ + // lookup the lower/upper bounds for the mesh element + double x_min = negative_grid_boundary(ijk, 0); + double x_max = positive_grid_boundary(ijk, 0); + + double y_min = (n_dimension_ >= 2) ? negative_grid_boundary(ijk, 1) : 0.0; + double y_max = (n_dimension_ >= 2) ? positive_grid_boundary(ijk, 1) : 0.0; + + double z_min = (n_dimension_ == 3) ? negative_grid_boundary(ijk, 2) : 0.0; + double z_max = (n_dimension_ == 3) ? positive_grid_boundary(ijk, 2) : 0.0; + + return {x_min + (x_max - x_min) * prn(seed), + y_min + (y_max - y_min) * prn(seed), z_min + (z_max - z_min) * prn(seed)}; +} + //============================================================================== // Unstructured Mesh implementation //============================================================================== @@ -381,11 +399,6 @@ 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"); -} - int StructuredMesh::get_bin(Position r) const { // Determine indices @@ -1077,6 +1090,30 @@ StructuredMesh::MeshIndex CylindricalMesh::get_indices( return idx; } +Position CylindricalMesh::sample_element( + const MeshIndex& ijk, uint64_t* seed) const +{ + double r_min = this->r(ijk[0] - 1); + double r_max = this->r(ijk[0]); + + double phi_min = this->phi(ijk[1] - 1); + double phi_max = this->phi(ijk[1]); + + double z_min = this->z(ijk[2] - 1); + double z_max = this->z(ijk[2]); + + double r_min_sq = r_min * r_min; + double r_max_sq = r_max * r_max; + double r = std::sqrt(uniform_distribution(r_min_sq, r_max_sq, seed)); + double phi = uniform_distribution(phi_min, phi_max, seed); + double z = uniform_distribution(z_min, z_max, seed); + + double x = r * std::cos(phi); + double y = r * std::sin(phi); + + return origin_ + Position(x, y, z); +} + double CylindricalMesh::find_r_crossing( const Position& r, const Direction& u, double l, int shell) const { @@ -1338,6 +1375,33 @@ StructuredMesh::MeshIndex SphericalMesh::get_indices( return idx; } +Position SphericalMesh::sample_element( + const MeshIndex& ijk, uint64_t* seed) const +{ + double r_min = this->r(ijk[0] - 1); + double r_max = this->r(ijk[0]); + + double theta_min = this->theta(ijk[1] - 1); + double theta_max = this->theta(ijk[1]); + + double phi_min = this->phi(ijk[2] - 1); + double phi_max = this->phi(ijk[2]); + + double cos_theta = uniform_distribution(theta_min, theta_max, seed); + double sin_theta = std::sin(std::acos(cos_theta)); + double phi = uniform_distribution(phi_min, phi_max, seed); + double r_min_cub = std::pow(r_min, 3); + double r_max_cub = std::pow(r_max, 3); + // might be faster to do rejection here? + double r = std::cbrt(uniform_distribution(r_min_cub, r_max_cub, seed)); + + double x = r * std::cos(phi) * sin_theta; + double y = r * std::sin(phi) * sin_theta; + double z = r * cos_theta; + + return origin_ + Position(x, y, z); +} + double SphericalMesh::find_r_crossing( const Position& r, const Direction& u, double l, int shell) const { @@ -2185,7 +2249,7 @@ std::string MOABMesh::library() const } // Sample position within a tet for MOAB type tets -Position MOABMesh::sample(uint64_t* seed, int32_t bin) const +Position MOABMesh::sample_element(int32_t bin, uint64_t* seed) const { moab::EntityHandle tet_ent = get_ent_handle_from_bin(bin); @@ -2672,7 +2736,7 @@ void LibMesh::initialize() } // Sample position within a tet for LibMesh type tets -Position LibMesh::sample(uint64_t* seed, int32_t bin) const +Position LibMesh::sample_element(int32_t bin, uint64_t* seed) const { const auto& elem = get_element_from_bin(bin); // Get tet vertex coordinates from LibMesh diff --git a/src/settings.cpp b/src/settings.cpp index c0a3aacf3cc..a5256c9c26d 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -457,28 +457,7 @@ void read_settings_xml(pugi::xml_node root) // Get point to list of elements and make sure there is at least one for (pugi::xml_node node : root.children("source")) { - if (check_for_node(node, "file")) { - auto path = get_node_value(node, "file", false, true); - if (ends_with(path, ".mcpl") || ends_with(path, ".mcpl.gz")) { - auto sites = mcpl_source_sites(path); - model::external_sources.push_back(make_unique(sites)); - } else { - model::external_sources.push_back(make_unique(path)); - } - } else if (check_for_node(node, "library")) { - // Get shared library path and parameters - auto path = get_node_value(node, "library", false, true); - std::string parameters; - if (check_for_node(node, "parameters")) { - parameters = get_node_value(node, "parameters", false, true); - } - - // Create custom source - model::external_sources.push_back( - make_unique(path, parameters)); - } else { - model::external_sources.push_back(make_unique(node)); - } + model::external_sources.push_back(Source::create(node)); } // Check if the user has specified to read surface source diff --git a/src/source.cpp b/src/source.cpp index 3ddd903bf4c..bf52b568a26 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -22,6 +22,7 @@ #include "openmc/geometry.h" #include "openmc/hdf5_interface.h" #include "openmc/material.h" +#include "openmc/mcpl_interface.h" #include "openmc/memory.h" #include "openmc/message_passing.h" #include "openmc/mgxs_interface.h" @@ -31,6 +32,7 @@ #include "openmc/settings.h" #include "openmc/simulation.h" #include "openmc/state_point.h" +#include "openmc/string_utils.h" #include "openmc/xml_interface.h" namespace openmc { @@ -44,6 +46,39 @@ namespace model { vector> external_sources; } +//============================================================================== +// Source create implementation +//============================================================================== + +unique_ptr Source::create(pugi::xml_node node) +{ + // if the source type is present, use it to determine the type + // of object to create + if (check_for_node(node, "type")) { + std::string source_type = get_node_value(node, "type"); + if (source_type == "independent") { + return make_unique(node); + } else if (source_type == "file") { + return make_unique(node); + } else if (source_type == "compiled") { + return make_unique(node); + } else if (source_type == "mesh") { + return make_unique(node); + } else { + fatal_error(fmt::format("Invalid source type '{}' found.", source_type)); + } + } else { + // support legacy source format + if (check_for_node(node, "file")) { + return make_unique(node); + } else if (check_for_node(node, "library")) { + return make_unique(node); + } else { + return make_unique(node); + } + } +} + //============================================================================== // IndependentSource implementation //============================================================================== @@ -81,32 +116,7 @@ IndependentSource::IndependentSource(pugi::xml_node node) // Spatial distribution for external source if (check_for_node(node, "space")) { - // Get pointer to spatial distribution - pugi::xml_node node_space = node.child("space"); - - // Check for type of spatial distribution and read - std::string type; - if (check_for_node(node_space, "type")) - type = get_node_value(node_space, "type", true, true); - if (type == "cartesian") { - space_ = UPtrSpace {new CartesianIndependent(node_space)}; - } else if (type == "cylindrical") { - 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") { - space_ = UPtrSpace {new SpatialBox(node_space, true)}; - } else if (type == "point") { - space_ = UPtrSpace {new SpatialPoint(node_space)}; - } else { - fatal_error(fmt::format( - "Invalid spatial distribution for external source: {}", type)); - } - + space_ = SpatialDistribution::create(node.child("space")); } else { // If no spatial distribution specified, make it a point source space_ = UPtrSpace {new SpatialPoint()}; @@ -114,24 +124,7 @@ IndependentSource::IndependentSource(pugi::xml_node node) // Determine external source angular distribution if (check_for_node(node, "angle")) { - // Get pointer to angular distribution - pugi::xml_node node_angle = node.child("angle"); - - // Check for type of angular distribution - std::string type; - if (check_for_node(node_angle, "type")) - type = get_node_value(node_angle, "type", true, true); - if (type == "isotropic") { - angle_ = UPtrAngle {new Isotropic()}; - } else if (type == "monodirectional") { - angle_ = UPtrAngle {new Monodirectional(node_angle)}; - } else if (type == "mu-phi") { - angle_ = UPtrAngle {new PolarAzimuthal(node_angle)}; - } else { - fatal_error(fmt::format( - "Invalid angular distribution for external source: {}", type)); - } - + angle_ = UnitSphereDistribution::create(node.child("angle")); } else { angle_ = UPtrAngle {new Isotropic()}; } @@ -290,8 +283,22 @@ SourceSite IndependentSource::sample(uint64_t* seed) const //============================================================================== // FileSource implementation //============================================================================== +FileSource::FileSource(pugi::xml_node node) +{ + auto path = get_node_value(node, "file", false, true); + if (ends_with(path, ".mcpl") || ends_with(path, ".mcpl.gz")) { + sites_ = mcpl_source_sites(path); + } else { + this->load_sites_from_file(path); + } +} -FileSource::FileSource(std::string path) +FileSource::FileSource(const std::string& path) +{ + load_sites_from_file(path); +} + +void FileSource::load_sites_from_file(const std::string& path) { // Check if source file exists if (!file_exists(path)) { @@ -328,9 +335,19 @@ SourceSite FileSource::sample(uint64_t* seed) const //============================================================================== // CompiledSourceWrapper implementation //============================================================================== +CompiledSourceWrapper::CompiledSourceWrapper(pugi::xml_node node) +{ + // Get shared library path and parameters + auto path = get_node_value(node, "library", false, true); + std::string parameters; + if (check_for_node(node, "parameters")) { + parameters = get_node_value(node, "parameters", false, true); + } + setup(path, parameters); +} -CompiledSourceWrapper::CompiledSourceWrapper( - std::string path, std::string parameters) +void CompiledSourceWrapper::setup( + const std::string& path, const std::string& parameters) { #ifdef HAS_DYNAMIC_LINKING // Open the library @@ -378,6 +395,50 @@ CompiledSourceWrapper::~CompiledSourceWrapper() #endif } +//============================================================================== +// MeshSource implementation +//============================================================================== + +MeshSource::MeshSource(pugi::xml_node node) +{ + int32_t mesh_id = stoi(get_node_value(node, "mesh")); + int32_t mesh_idx = model::mesh_map.at(mesh_id); + const auto& mesh = model::meshes[mesh_idx]; + + std::vector strengths; + // read all source distributions and populate strengths vector for MeshSpatial + // object + for (auto source_node : node.children("source")) { + sources_.emplace_back(Source::create(source_node)); + strengths.push_back(sources_.back()->strength()); + } + + // the number of source distributions should either be one or equal to the + // number of mesh elements + if (sources_.size() > 1 && sources_.size() != mesh->n_bins()) { + fatal_error(fmt::format("Incorrect number of source distributions ({}) for " + "mesh source with {} elements.", + sources_.size(), mesh->n_bins())); + } + + space_ = std::make_unique(mesh_idx, strengths); +} + +SourceSite MeshSource::sample(uint64_t* seed) const +{ + // sample location and element from mesh + auto mesh_location = space_->sample_mesh(seed); + + // Sample source for the chosen element + int32_t element = mesh_location.first; + SourceSite site = source(element)->sample(seed); + + // Replace spatial position with the one already sampled + site.r = mesh_location.second; + + return site; +} + //============================================================================== // Non-member functions //============================================================================== diff --git a/tests/unit_tests/test_source_mesh.py b/tests/unit_tests/test_source_mesh.py index fd8b04148e7..7bec2d0724b 100644 --- a/tests/unit_tests/test_source_mesh.py +++ b/tests/unit_tests/test_source_mesh.py @@ -1,6 +1,5 @@ from itertools import product from pathlib import Path -from subprocess import call import pytest import numpy as np @@ -8,9 +7,11 @@ import openmc.lib from tests import cdtemp -from tests.regression_tests import config +################### +# MeshSpatial Tests +################### TETS_PER_VOXEL = 12 # This test uses a geometry file with cells that match a regular mesh. Each cell @@ -48,9 +49,9 @@ def model(): settings.particles = 100 settings.batches = 2 - return openmc.model.Model(geometry=geometry, - materials=materials, - settings=settings) + return openmc.Model(geometry=geometry, + materials=materials, + settings=settings) ### Setup test cases ### param_values = (['libmesh', 'moab'], # mesh libraries @@ -174,6 +175,7 @@ def test_strengths_size_failure(request, model): 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.") @@ -201,3 +203,158 @@ def test_roundtrip(run_in_tmpdir, model, request): assert space_in.mesh.id == space_out.mesh.id assert space_in.volume_normalized == space_out.volume_normalized + + +################### +# MeshSource tests +################### +@pytest.mark.parametrize('mesh_type', ('rectangular', 'cylindrical')) +def test_mesh_source_independent(run_in_tmpdir, mesh_type): + """ + A void model containing a single box + """ + min, max = -10, 10 + box = openmc.model.RectangularParallelepiped( + min, max, min, max, min, max, boundary_type='vacuum') + + geometry = openmc.Geometry([openmc.Cell(region=-box)]) + + settings = openmc.Settings() + settings.particles = 100 + settings.batches = 10 + settings.run_mode = 'fixed source' + + model = openmc.Model(geometry=geometry, settings=settings) + + # define a 2 x 2 x 2 mesh + if mesh_type == 'rectangular': + mesh = openmc.RegularMesh.from_domain(model.geometry, (2, 2, 2)) + elif mesh_type == 'cylindrical': + mesh = openmc.CylindricalMesh.from_domain(model.geometry, (1, 4, 2)) + + energy = openmc.stats.Discrete([1.e6], [1.0]) + + # create sources with only one non-zero strength for the source in the mesh + # voxel occupying the lowest octant. Direct source particles straight out of + # the problem from there. This demonstrates that + # 1) particles are only being sourced within the intented mesh voxel based + # on source strength + # 2) particles are respecting the angle distributions assigned to each voxel + sources = np.empty(mesh.dimension, dtype=openmc.SourceBase) + centroids = mesh.centroids + x, y, z = np.swapaxes(mesh.centroids, -1, 0) + for i, j, k in mesh.indices: + # mesh.indices is currently one-indexed, adjust for Python arrays + ijk = (i-1, j-1, k-1) + + # get the centroid of the ijk mesh element and use it to set the + # direction of the source directly out of the problem + centroid = centroids[ijk] + vec = np.sign(centroid, dtype=float) + vec /= np.linalg.norm(vec) + angle = openmc.stats.Monodirectional(vec) + sources[ijk] = openmc.IndependentSource(energy=energy, angle=angle, strength=0.0) + + # create and apply the mesh source + mesh_source = openmc.MeshSource(mesh, sources) + model.settings.source = mesh_source + + # tally the flux on the mesh + mesh_filter = openmc.MeshFilter(mesh) + tally = openmc.Tally() + tally.filters = [mesh_filter] + tally.scores = ['flux'] + + model.tallies = openmc.Tallies([tally]) + + # for each element, set a single-non zero source with particles + # traveling out of the mesh (and geometry) w/o crossing any other + # mesh elements + for i, j, k in mesh.indices: + ijk = (i-1, j-1, k-1) + # zero-out all source strengths and set the strength + # on the element of interest + mesh_source.strength = 0.0 + mesh_source.sources[ijk].strength = 1.0 + + sp_file = model.run() + + with openmc.StatePoint(sp_file) as sp: + tally_out = sp.get_tally(id=tally.id) + mean = tally_out.get_reshaped_data(expand_dims=True) + + # remove nuclides and scores axes + mean = mean[..., 0, 0] + # the mesh elment with a non-zero source strength should have a value + assert mean[ijk] != 0 + # all other values should be zero + mean[ijk] = 0 + assert np.all(mean == 0), f'Failed on index {ijk} with centroid {mesh.centroids[ijk]}' + + # test roundtrip + xml_model = openmc.Model.from_model_xml() + xml_source = xml_model.settings.source[0] + assert isinstance(xml_source, openmc.MeshSource) + assert xml_source.strength == 1.0 + assert isinstance(xml_source.mesh, type(mesh_source.mesh)) + assert xml_source.mesh.dimension == mesh_source.mesh.dimension + assert xml_source.mesh.id == mesh_source.mesh.id + assert len(xml_source.sources) == len(mesh_source.sources) + + # check strength adjustment methods + assert mesh_source.strength == 1.0 + mesh_source.strength = 100.0 + assert mesh_source.strength == 100.0 + + mesh_source.normalize_source_strengths() + assert mesh_source.strength == 1.0 + + +def test_mesh_source_file(run_in_tmpdir): + # Creating a source file with a single particle + source_particle = openmc.SourceParticle(time=10.0) + openmc.write_source_file([source_particle], 'source.h5') + file_source = openmc.FileSource('source.h5') + + model = openmc.Model() + + rect_prism = openmc.model.RectangularParallelepiped( + -5.0, 5.0, -5.0, 5.0, -5.0, 5.0, boundary_type='vacuum') + + mat = openmc.Material() + mat.add_nuclide('H1', 1.0) + + model.geometry = openmc.Geometry([openmc.Cell(fill=mat, region=-rect_prism)]) + model.settings.particles = 1000 + model.settings.batches = 10 + model.settings.run_mode = 'fixed source' + + mesh = openmc.RegularMesh() + mesh.lower_left = (-1, -2, -3) + mesh.upper_right = (2, 3, 4) + mesh.dimension = (1, 1, 1) + + mesh_source_arr = np.asarray([file_source]).reshape(mesh.dimension) + source = openmc.MeshSource(mesh, mesh_source_arr) + + model.settings.source = source + + model.export_to_model_xml() + + openmc.lib.init() + openmc.lib.simulation_init() + sites = openmc.lib.sample_external_source(10) + openmc.lib.simulation_finalize() + openmc.lib.finalize() + + # The mesh bounds do not contain the point of the lone source site in the + # file source, so it should not appear in the set of source sites produced + # from the mesh source. Additionally, the source should be located within + # the mesh + bbox = mesh.bounding_box + for site in sites: + assert site.r != (0, 0, 0) + assert site.E == source_particle.E + assert site.u == source_particle.u + assert site.time == source_particle.time + assert site.r in bbox