diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 6318a991601..ae3266dabed 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -534,8 +534,9 @@ attributes/sub-elements: *Default*: 1.0 :type: - 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. + 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``. @@ -716,6 +717,39 @@ attributes/sub-elements: mesh element and follows the format for :ref:`source_element`. The number of ```` sub-elements should correspond to the number of mesh elements. + :constraints: + This sub-element indicates the presence of constraints on sampled source + sites (see :ref:`usersguide_source_constraints` for details). It may have + the following sub-elements: + + :domain_ids: + The unique IDs of domains for which source sites must be within. + + *Default*: None + + :domain_type: + The type of each domain for source rejection ("cell", "material", or + "universe"). + + *Default*: None + + :fissionable: + A boolean indicating whether source sites must be sampled within a + material that is fissionable in order to be accepted. + + :time_bounds: + A pair of times in [s] indicating the lower and upper bound for a time + interval that source particles must be within. + + :energy_bounds: + A pair of energies in [eV] indicating the lower and upper bound for an + energy interval that source particles must be within. + + :rejection_strategy: + Either "resample", indicating that source sites should be resampled when + one is rejected, or "kill", indicating that a rejected source site is + assigned zero weight. + .. _univariate: Univariate Probability Distributions diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index e966423f0e3..eb02f654c30 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -420,6 +420,54 @@ string, which gets passed down to the ``openmc_create_source()`` function:: settings.source = openmc.CompiledSource('libsource.so', '3.5e6') +.. _usersguide_source_constraints: + +Source Constraints +------------------ + +All source classes in OpenMC have the ability to apply a set of "constraints" +that limit which sampled source sites are actually used for transport. The most +common use case is to sample source sites over some simple spatial distribution +(e.g., uniform over a box) and then only accept those that appear in a given +cell or material. This can be done with a domain constraint, which can be +specified as follows:: + + source_cell = openmc.Cell(...) + ... + + spatial_dist = openmc.stats.Box((-10., -10., -10.), (10., 10., 10.)) + source = openmc.IndependentSource( + space=spatial_dist, + constraints={'domains': [source_cell]} + ) + +For k-eigenvalue problems, a convenient constraint is available that limits +source sites to those sampled in a fissionable material:: + + source = openmc.IndependentSource( + space=spatial_dist, constraints={'fissionable': True} + ) + +Constraints can also be placed on a range of energies or times:: + + # Only use source sites between 500 keV and 1 MeV and with times under 1 sec + source = openmc.FileSource( + 'source.h5', + constraints={'energy_bounds': [500.0e3, 1.0e6], 'time_bounds': [0.0, 1.0]} + ) + +Normally, when a source site is rejected, a new one will be resampled until one +is found that meets the constraints. However, the rejection strategy can be +changed so that a rejected site will just not be simulated by specifying:: + + source = openmc.IndependentSource( + space=spatial_dist, + constraints={'domains': [cell], 'rejection_strategy': 'kill'} + ) + +In this case, the actual number of particles simulated may be less than what you +specified in :attr:`Settings.particles`. + .. _usersguide_entropy: --------------- diff --git a/examples/assembly/assembly.py b/examples/assembly/assembly.py index a370a361144..d94cb69bfc4 100644 --- a/examples/assembly/assembly.py +++ b/examples/assembly/assembly.py @@ -112,11 +112,10 @@ def assembly_model(): model.settings.batches = 150 model.settings.inactive = 50 model.settings.particles = 1000 - model.settings.source = openmc.IndependentSource(space=openmc.stats.Box( - (-pitch/2, -pitch/2, -1), - (pitch/2, pitch/2, 1), - only_fissionable=True - )) + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Box((-pitch/2, -pitch/2, -1), (pitch/2, pitch/2, 1)), + constraints={'fissionable': True} + ) # NOTE: We never actually created a Materials object. When you export/run # using the Model object, if no materials were assigned it will look through diff --git a/examples/pincell_depletion/restart_depletion.py b/examples/pincell_depletion/restart_depletion.py index ac065d3402f..de9fc16cb1f 100644 --- a/examples/pincell_depletion/restart_depletion.py +++ b/examples/pincell_depletion/restart_depletion.py @@ -26,8 +26,9 @@ # Create an initial uniform spatial source distribution over fissionable zones bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] -uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True) -settings.source = openmc.IndependentSource(space=uniform_dist) +uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:]) +settings.source = openmc.IndependentSource( + space=uniform_dist, constraints={'fissionable': True}) entropy_mesh = openmc.RegularMesh() entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50] diff --git a/examples/pincell_depletion/run_depletion.py b/examples/pincell_depletion/run_depletion.py index e1959938f05..013c83c86a5 100644 --- a/examples/pincell_depletion/run_depletion.py +++ b/examples/pincell_depletion/run_depletion.py @@ -71,8 +71,9 @@ # Create an initial uniform spatial source distribution over fissionable zones bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] -uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True) -settings.source = openmc.IndependentSource(space=uniform_dist) +uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:]) +settings.source = openmc.IndependentSource( + space=uniform_dist, constraints={'fissionable': True}) entropy_mesh = openmc.RegularMesh() entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50] diff --git a/include/openmc/distribution_spatial.h b/include/openmc/distribution_spatial.h index bd10a299605..28568c890d8 100644 --- a/include/openmc/distribution_spatial.h +++ b/include/openmc/distribution_spatial.h @@ -116,6 +116,11 @@ class MeshSpatial : public SpatialDistribution { //! \return Sampled element index and position within that element std::pair sample_mesh(uint64_t* seed) const; + //! Sample a mesh element + //! \param seed Pseudorandom number seed pointer + //! \return Sampled element index + int32_t sample_element_index(uint64_t* seed) const; + //! For unstructured meshes, ensure that elements are all linear tetrahedra void check_element_types() const; diff --git a/include/openmc/source.h b/include/openmc/source.h index d5ea9bd1d56..9ef8b925127 100644 --- a/include/openmc/source.h +++ b/include/openmc/source.h @@ -4,6 +4,7 @@ #ifndef OPENMC_SOURCE_H #define OPENMC_SOURCE_H +#include #include #include "pugixml.hpp" @@ -39,19 +40,72 @@ extern vector> external_sources; //============================================================================== //! Abstract source interface +// +//! The Source class provides the interface that must be implemented by derived +//! classes, namely the sample() method that returns a sampled source site. From +//! this base class, source rejection is handled within the +//! sample_with_constraints() method. However, note that some classes directly +//! check for constraints for efficiency reasons (like IndependentSource), in +//! which case the constraints_applied() method indicates that constraints +//! should not be checked a second time from the base class. //============================================================================== class Source { public: + // Constructors, destructors + Source() = default; + explicit Source(pugi::xml_node node); virtual ~Source() = default; - // Methods that must be implemented - virtual SourceSite sample(uint64_t* seed) const = 0; - // Methods that can be overridden - virtual double strength() const { return 1.0; } + virtual double strength() const { return strength_; } + + //! Sample a source site and apply constraints + // + //! \param[inout] seed Pseudorandom seed pointer + //! \return Sampled site + SourceSite sample_with_constraints(uint64_t* seed) const; + + //! Sample a source site (without applying constraints) + // + //! Sample from the external source distribution + //! \param[inout] seed Pseudorandom seed pointer + //! \return Sampled site + virtual SourceSite sample(uint64_t* seed) const = 0; static unique_ptr create(pugi::xml_node node); + +protected: + // Domain types + enum class DomainType { UNIVERSE, MATERIAL, CELL }; + + // Strategy used for rejecting sites when constraints are applied. KILL means + // that sites are always accepted but if they don't satisfy constraints, they + // are given weight 0. RESAMPLE means that a new source site will be sampled + // until constraints are met. + enum class RejectionStrategy { KILL, RESAMPLE }; + + // Indicates whether derived class already handles constraints + virtual bool constraints_applied() const { return false; } + + // Methods for constraints + void read_constraints(pugi::xml_node node); + bool satisfies_spatial_constraints(Position r) const; + bool satisfies_energy_constraints(double E) const; + bool satisfies_time_constraints(double time) const; + + // Data members + double strength_ {1.0}; //!< Source strength + std::unordered_set domain_ids_; //!< Domains to reject from + DomainType domain_type_; //!< Domain type for rejection + std::pair time_bounds_ {-std::numeric_limits::max(), + std::numeric_limits::max()}; //!< time limits + std::pair energy_bounds_ { + 0, std::numeric_limits::max()}; //!< energy limits + bool only_fissionable_ { + false}; //!< Whether site must be in fissionable material + RejectionStrategy rejection_strategy_ { + RejectionStrategy::RESAMPLE}; //!< Procedure for rejecting }; //============================================================================== @@ -73,7 +127,6 @@ class IndependentSource : public Source { // Properties ParticleType particle_type() const { return particle_; } - double strength() const override { return strength_; } // Make observing pointers available SpatialDistribution* space() const { return space_.get(); } @@ -81,19 +134,17 @@ class IndependentSource : public Source { Distribution* energy() const { return energy_.get(); } Distribution* time() const { return time_.get(); } -private: - // Domain types - enum class DomainType { UNIVERSE, MATERIAL, CELL }; +protected: + // Indicates whether derived class already handles constraints + bool constraints_applied() const override { return true; } +private: // Data members ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted - double strength_ {1.0}; //!< Source strength UPtrSpace space_; //!< Spatial distribution UPtrAngle angle_; //!< Angular distribution UPtrDist energy_; //!< Energy distribution UPtrDist time_; //!< Time distribution - DomainType domain_type_; //!< Domain type for rejection - std::unordered_set domain_ids_; //!< Domains to reject from }; //============================================================================== @@ -107,9 +158,12 @@ class FileSource : public Source { 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 + +protected: + SourceSite sample(uint64_t* seed) const override; + private: vector sites_; //!< Source sites from a file }; @@ -124,16 +178,17 @@ class CompiledSourceWrapper : public Source { CompiledSourceWrapper(pugi::xml_node node); ~CompiledSourceWrapper(); + double strength() const override { return compiled_source_->strength(); } + + void setup(const std::string& path, const std::string& parameters); + +protected: // Defer implementation to custom source library SourceSite sample(uint64_t* seed) const override { return compiled_source_->sample(seed); } - 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_; @@ -164,6 +219,9 @@ class MeshSource : public Source { return sources_.size() == 1 ? sources_[0] : sources_[i]; } +protected: + bool constraints_applied() const override { return true; } + private: // Data members unique_ptr space_; //!< Mesh spatial diff --git a/openmc/examples.py b/openmc/examples.py index cef86a47a87..fc94b8b528e 100644 --- a/openmc/examples.py +++ b/openmc/examples.py @@ -76,8 +76,10 @@ def pwr_pin_cell(): model.settings.batches = 10 model.settings.inactive = 5 model.settings.particles = 100 - model.settings.source = openmc.IndependentSource(space=openmc.stats.Box( - [-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1], only_fissionable=True)) + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Box([-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1]), + constraints={'fissionable': True} + ) plot = openmc.Plot.from_geometry(model.geometry) plot.pixels = (300, 300) @@ -527,8 +529,10 @@ def pwr_assembly(): model.settings.batches = 10 model.settings.inactive = 5 model.settings.particles = 100 - model.settings.source = openmc.IndependentSource(space=openmc.stats.Box( - [-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1], only_fissionable=True)) + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Box([-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1]), + constraints={'fissionable': True} + ) plot = openmc.Plot() plot.origin = (0.0, 0.0, 0) diff --git a/openmc/source.py b/openmc/source.py index 5deb6b05f50..91318f8e9ad 100644 --- a/openmc/source.py +++ b/openmc/source.py @@ -6,7 +6,7 @@ import warnings import typing # imported separately as py3.8 requires typing.Iterable # also required to prevent typing.Union namespace overwriting Union -from typing import Optional, Sequence +from typing import Optional, Sequence, Dict, Any import lxml.etree as ET import numpy as np @@ -28,6 +28,19 @@ class SourceBase(ABC): ---------- strength : float Strength of the source + constraints : dict + Constraints on sampled source particles. Valid keys include 'domains', + 'time_bounds', 'energy_bounds', 'fissionable', and 'rejection_strategy'. + For 'domains', the corresponding value is an iterable of + :class:`openmc.Cell`, :class:`openmc.Material`, or + :class:`openmc.Universe` for which sampled sites must be within. For + 'time_bounds' and 'energy_bounds', the corresponding value is a sequence + of floats giving the lower and upper bounds on time in [s] or energy in + [eV] that the sampled particle must be within. For 'fissionable', the + value is a bool indicating that only sites in fissionable material + should be accepted. The 'rejection_strategy' indicates what should + happen when a source particle is rejected: either 'resample' (pick a new + particle) or 'kill' (accept and terminate). Attributes ---------- @@ -35,11 +48,20 @@ class SourceBase(ABC): Indicator of source type. strength : float Strength of the source + constraints : dict + Constraints on sampled source particles. Valid keys include + 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', + 'fissionable', and 'rejection_strategy'. """ - def __init__(self, strength=1.0): + def __init__( + self, + strength: Optional[float] = 1.0, + constraints: Optional[Dict[str, Any]] = None + ): self.strength = strength + self.constraints = constraints @property def strength(self): @@ -47,10 +69,47 @@ def strength(self): @strength.setter def strength(self, strength): - cv.check_type('source strength', strength, Real) - cv.check_greater_than('source strength', strength, 0.0, True) + cv.check_type('source strength', strength, Real, none_ok=True) + if strength is not None: + cv.check_greater_than('source strength', strength, 0.0, True) self._strength = strength + @property + def constraints(self) -> Dict[str, Any]: + return self._constraints + + @constraints.setter + def constraints(self, constraints: Optional[Dict[str, Any]]): + self._constraints = {} + if constraints is None: + return + + for key, value in constraints.items(): + if key == 'domains': + cv.check_type('domains', value, Iterable, + (openmc.Cell, openmc.Material, openmc.Universe)) + if isinstance(value[0], openmc.Cell): + self._constraints['domain_type'] = 'cell' + elif isinstance(value[0], openmc.Material): + self._constraints['domain_type'] = 'material' + elif isinstance(value[0], openmc.Universe): + self._constraints['domain_type'] = 'universe' + self._constraints['domain_ids'] = [d.id for d in value] + elif key == 'time_bounds': + cv.check_type('time bounds', value, Iterable, Real) + self._constraints['time_bounds'] = tuple(value) + elif key == 'energy_bounds': + cv.check_type('energy bounds', value, Iterable, Real) + self._constraints['energy_bounds'] = tuple(value) + elif key == 'fissionable': + cv.check_type('fissionable', value, bool) + self._constraints['fissionable'] = value + elif key == 'rejection_strategy': + cv.check_value('rejection strategy', value, ('resample', 'kill')) + self._constraints['rejection_strategy'] = value + else: + raise ValueError('Unknown key in constraints dictionary: {key}') + @abstractmethod def populate_xml_element(self, element): """Add necessary source information to an XML element @@ -73,8 +132,30 @@ def to_xml_element(self) -> ET.Element: """ element = ET.Element("source") element.set("type", self.type) - element.set("strength", str(self.strength)) + if self.strength is not None: + element.set("strength", str(self.strength)) self.populate_xml_element(element) + constraints = self.constraints + if constraints: + constraints_elem = ET.SubElement(element, "constraints") + if "domain_ids" in constraints: + dt_elem = ET.SubElement(constraints_elem, "domain_type") + dt_elem.text = constraints["domain_type"] + id_elem = ET.SubElement(constraints_elem, "domain_ids") + id_elem.text = ' '.join(str(uid) for uid in constraints["domain_ids"]) + if "time_bounds" in constraints: + dt_elem = ET.SubElement(constraints_elem, "time_bounds") + dt_elem.text = ' '.join(str(t) for t in constraints["time_bounds"]) + if "energy_bounds" in constraints: + dt_elem = ET.SubElement(constraints_elem, "energy_bounds") + dt_elem.text = ' '.join(str(E) for E in constraints["energy_bounds"]) + if "fissionable" in constraints: + dt_elem = ET.SubElement(constraints_elem, "fissionable") + dt_elem.text = str(constraints["fissionable"]).lower() + if "rejection_strategy" in constraints: + dt_elem = ET.SubElement(constraints_elem, "rejection_strategy") + dt_elem.text = constraints["rejection_strategy"] + return element @classmethod @@ -118,6 +199,47 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase: else: raise ValueError(f'Source type {source_type} is not recognized') + @staticmethod + def _get_constraints(elem: ET.Element) -> Dict[str, Any]: + # Find element containing constraints + constraints_elem = elem.find("constraints") + elem = constraints_elem if constraints_elem is not None else elem + + constraints = {} + domain_type = get_text(elem, "domain_type") + if domain_type is not None: + domain_ids = [int(x) for x in get_text(elem, "domain_ids").split()] + + # Instantiate some throw-away domains that are used by the + # constructor to assign IDs + with warnings.catch_warnings(): + warnings.simplefilter('ignore', openmc.IDWarning) + if domain_type == 'cell': + domains = [openmc.Cell(uid) for uid in domain_ids] + elif domain_type == 'material': + domains = [openmc.Material(uid) for uid in domain_ids] + elif domain_type == 'universe': + domains = [openmc.Universe(uid) for uid in domain_ids] + constraints['domains'] = domains + + time_bounds = get_text(elem, "time_bounds") + if time_bounds is not None: + constraints['time_bounds'] = [float(x) for x in time_bounds.split()] + + energy_bounds = get_text(elem, "energy_bounds") + if energy_bounds is not None: + constraints['energy_bounds'] = [float(x) for x in energy_bounds.split()] + + fissionable = get_text(elem, "fissionable") + if fissionable is not None: + constraints['fissionable'] = fissionable in ('true', '1') + + rejection_strategy = get_text(elem, "rejection_strategy") + if rejection_strategy is not None: + constraints['rejection_strategy'] = rejection_strategy + + return constraints + class IndependentSource(SourceBase): """Distribution of phase space coordinates for source sites. @@ -142,6 +264,22 @@ class IndependentSource(SourceBase): Domains to reject based on, i.e., if a sampled spatial location is not within one of these domains, it will be rejected. + .. deprecated:: 0.14.1 + Use the `constraints` argument instead. + constraints : dict + Constraints on sampled source particles. Valid keys include 'domains', + 'time_bounds', 'energy_bounds', 'fissionable', and 'rejection_strategy'. + For 'domains', the corresponding value is an iterable of + :class:`openmc.Cell`, :class:`openmc.Material`, or + :class:`openmc.Universe` for which sampled sites must be within. For + 'time_bounds' and 'energy_bounds', the corresponding value is a sequence + of floats giving the lower and upper bounds on time in [s] or energy in + [eV] that the sampled particle must be within. For 'fissionable', the + value is a bool indicating that only sites in fissionable material + should be accepted. The 'rejection_strategy' indicates what should + happen when a source particle is rejected: either 'resample' (pick a new + particle) or 'kill' (accept and terminate). + Attributes ---------- space : openmc.stats.Spatial or None @@ -161,10 +299,10 @@ class IndependentSource(SourceBase): particle : {'neutron', 'photon'} Source particle type - ids : Iterable of int - IDs of domains to use for rejection - domain_type : {'cell', 'material', 'universe'} - Type of domain to use for rejection + constraints : dict + Constraints on sampled source particles. Valid keys include + 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', + 'fissionable', and 'rejection_strategy'. """ @@ -176,9 +314,15 @@ def __init__( time: Optional[openmc.stats.Univariate] = None, strength: float = 1.0, particle: str = 'neutron', - domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None + domains: Optional[Sequence[typing.Union[openmc.Cell, openmc.Material, openmc.Universe]]] = None, + constraints: Optional[Dict[str, Any]] = None ): - super().__init__(strength) + if domains is not None: + warnings.warn("The 'domains' arguments has been replaced by the " + "'constraints' argument.", FutureWarning) + constraints = {'domains': domains} + + super().__init__(strength=strength, constraints=constraints) self._space = None self._angle = None @@ -193,20 +337,8 @@ def __init__( self.energy = energy if time is not None: self.time = time - self.strength = strength self.particle = particle - self._domain_ids = [] - self._domain_type = None - if domains is not None: - if isinstance(domains[0], openmc.Cell): - self.domain_type = 'cell' - elif isinstance(domains[0], openmc.Material): - self.domain_type = 'material' - elif isinstance(domains[0], openmc.Universe): - self.domain_type = 'universe' - self.domain_ids = [d.id for d in domains] - @property def type(self) -> str: return 'independent' @@ -273,24 +405,6 @@ def particle(self, particle): cv.check_value('source particle', particle, ['neutron', 'photon']) self._particle = particle - @property - def domain_ids(self): - return self._domain_ids - - @domain_ids.setter - def domain_ids(self, ids): - cv.check_type('domain IDs', ids, Iterable, Real) - self._domain_ids = ids - - @property - def domain_type(self): - return self._domain_type - - @domain_type.setter - def domain_type(self, domain_type): - cv.check_value('domain type', domain_type, ('cell', 'material', 'universe')) - self._domain_type = domain_type - def populate_xml_element(self, element): """Add necessary source information to an XML element @@ -309,11 +423,6 @@ def populate_xml_element(self, element): element.append(self.energy.to_xml_element('energy')) if self.time is not None: element.append(self.time.to_xml_element('time')) - if self.domain_ids: - dt_elem = ET.SubElement(element, "domain_type") - dt_elem.text = self.domain_type - id_elem = ET.SubElement(element, "domain_ids") - id_elem.text = ' '.join(str(uid) for uid in self.domain_ids) @classmethod def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase: @@ -333,24 +442,8 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase: Source generated from XML element """ - domain_type = get_text(elem, "domain_type") - if domain_type is not None: - domain_ids = [int(x) for x in get_text(elem, "domain_ids").split()] - - # Instantiate some throw-away domains that are used by the - # constructor to assign IDs - with warnings.catch_warnings(): - warnings.simplefilter('ignore', openmc.IDWarning) - if domain_type == 'cell': - domains = [openmc.Cell(uid) for uid in domain_ids] - elif domain_type == 'material': - domains = [openmc.Material(uid) for uid in domain_ids] - elif domain_type == 'universe': - domains = [openmc.Universe(uid) for uid in domain_ids] - else: - domains = None - - source = cls(domains=domains) + constraints = cls._get_constraints(elem) + source = cls(constraints=constraints) strength = get_text(elem, 'strength') if strength is not None: @@ -360,10 +453,6 @@ def from_xml_element(cls, elem: ET.Element, meshes=None) -> SourceBase: if particle is not None: source.particle = particle - filename = get_text(elem, 'file') - if filename is not None: - source.file = filename - space = elem.find('space') if space is not None: source.space = Spatial.from_xml_element(space, meshes) @@ -405,6 +494,19 @@ class MeshSource(SourceBase): multidimensional array whose shape matches the mesh shape. If spatial distributions are set on any of the source objects, they will be ignored during source site sampling. + constraints : dict + Constraints on sampled source particles. Valid keys include 'domains', + 'time_bounds', 'energy_bounds', 'fissionable', and 'rejection_strategy'. + For 'domains', the corresponding value is an iterable of + :class:`openmc.Cell`, :class:`openmc.Material`, or + :class:`openmc.Universe` for which sampled sites must be within. For + 'time_bounds' and 'energy_bounds', the corresponding value is a sequence + of floats giving the lower and upper bounds on time in [s] or energy in + [eV] that the sampled particle must be within. For 'fissionable', the + value is a bool indicating that only sites in fissionable material + should be accepted. The 'rejection_strategy' indicates what should + happen when a source particle is rejected: either 'resample' (pick a new + particle) or 'kill' (accept and terminate). Attributes ---------- @@ -416,9 +518,19 @@ class MeshSource(SourceBase): Strength of the source type : str Indicator of source type: 'mesh' + constraints : dict + Constraints on sampled source particles. Valid keys include + 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', + 'fissionable', and 'rejection_strategy'. """ - def __init__(self, mesh: MeshBase, sources: Sequence[SourceBase]): + def __init__( + self, + mesh: MeshBase, + sources: Sequence[SourceBase], + constraints: Optional[Dict[str, Any]] = None, + ): + super().__init__(strength=None, constraints=constraints) self.mesh = mesh self.sources = sources @@ -473,8 +585,9 @@ def sources(self, s): @strength.setter def strength(self, val): - cv.check_type('mesh source strength', val, Real) - self.set_total_strength(val) + if val is not None: + 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. @@ -531,7 +644,8 @@ def from_xml_element(cls, elem: ET.Element, meshes) -> openmc.MeshSource: mesh = meshes[mesh_id] sources = [SourceBase.from_xml_element(e) for e in elem.iterchildren('source')] - return cls(mesh, sources) + constraints = cls._get_constraints(elem) + return cls(mesh, sources, constraints=constraints) def Source(*args, **kwargs): @@ -556,6 +670,19 @@ class CompiledSource(SourceBase): Parameters to be provided to the compiled shared library function strength : float Strength of the source + constraints : dict + Constraints on sampled source particles. Valid keys include 'domains', + 'time_bounds', 'energy_bounds', 'fissionable', and 'rejection_strategy'. + For 'domains', the corresponding value is an iterable of + :class:`openmc.Cell`, :class:`openmc.Material`, or + :class:`openmc.Universe` for which sampled sites must be within. For + 'time_bounds' and 'energy_bounds', the corresponding value is a sequence + of floats giving the lower and upper bounds on time in [s] or energy in + [eV] that the sampled particle must be within. For 'fissionable', the + value is a bool indicating that only sites in fissionable material + should be accepted. The 'rejection_strategy' indicates what should + happen when a source particle is rejected: either 'resample' (pick a new + particle) or 'kill' (accept and terminate). Attributes ---------- @@ -567,10 +694,20 @@ class CompiledSource(SourceBase): Strength of the source type : str Indicator of source type: 'compiled' + constraints : dict + Constraints on sampled source particles. Valid keys include + 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', + 'fissionable', and 'rejection_strategy'. """ - def __init__(self, library: Optional[str] = None, parameters: Optional[str] = None, strength=1.0) -> None: - super().__init__(strength=strength) + def __init__( + self, + library: Optional[str] = None, + parameters: Optional[str] = None, + strength: float = 1.0, + constraints: Optional[Dict[str, Any]] = None + ) -> None: + super().__init__(strength=strength, constraints=constraints) self._library = None if library is not None: @@ -634,9 +771,10 @@ def from_xml_element(cls, elem: ET.Element) -> openmc.CompiledSource: Source generated from XML element """ - library = get_text(elem, 'library') + kwargs = {'constraints': cls._get_constraints(elem)} + kwargs['library'] = get_text(elem, 'library') - source = cls(library) + source = cls(**kwargs) strength = get_text(elem, 'strength') if strength is not None: @@ -660,6 +798,19 @@ class FileSource(SourceBase): Path to the source file from which sites should be sampled strength : float Strength of the source (default is 1.0) + constraints : dict + Constraints on sampled source particles. Valid keys include 'domains', + 'time_bounds', 'energy_bounds', 'fissionable', and 'rejection_strategy'. + For 'domains', the corresponding value is an iterable of + :class:`openmc.Cell`, :class:`openmc.Material`, or + :class:`openmc.Universe` for which sampled sites must be within. For + 'time_bounds' and 'energy_bounds', the corresponding value is a sequence + of floats giving the lower and upper bounds on time in [s] or energy in + [eV] that the sampled particle must be within. For 'fissionable', the + value is a bool indicating that only sites in fissionable material + should be accepted. The 'rejection_strategy' indicates what should + happen when a source particle is rejected: either 'resample' (pick a new + particle) or 'kill' (accept and terminate). Attributes ---------- @@ -669,13 +820,21 @@ class FileSource(SourceBase): Strength of the source type : str Indicator of source type: 'file' + constraints : dict + Constraints on sampled source particles. Valid keys include + 'domain_type', 'domain_ids', 'time_bounds', 'energy_bounds', + 'fissionable', and 'rejection_strategy'. """ - def __init__(self, path: Optional[PathLike] = None, strength=1.0) -> None: - super().__init__(strength=strength) + def __init__( + self, + path: Optional[PathLike] = None, + strength: float = 1.0, + constraints: Optional[Dict[str, Any]] = None + ): + super().__init__(strength=strength, constraints=constraints) self._path = None - if path is not None: self.path = path @@ -722,16 +881,13 @@ def from_xml_element(cls, elem: ET.Element) -> openmc.FileSource: Source generated from XML element """ - - filename = get_text(elem, 'file') - - source = cls(filename) - + kwargs = {'constraints': cls._get_constraints(elem)} + kwargs['path'] = get_text(elem, 'file') strength = get_text(elem, 'strength') if strength is not None: - source.strength = float(strength) + kwargs['strength'] = float(strength) - return source + return cls(**kwargs) class ParticleType(IntEnum): diff --git a/openmc/stats/multivariate.py b/openmc/stats/multivariate.py index 3922d601aa8..212083c3cff 100644 --- a/openmc/stats/multivariate.py +++ b/openmc/stats/multivariate.py @@ -4,6 +4,7 @@ from collections.abc import Iterable from math import cos, pi from numbers import Real +from warnings import warn import lxml.etree as ET import numpy as np @@ -768,6 +769,9 @@ class Box(Spatial): Whether spatial sites should only be accepted if they occur in fissionable materials + .. deprecated:: 0.14.1 + Use the `constraints` argument when defining a source object instead. + Attributes ---------- lower_left : Iterable of float @@ -778,6 +782,9 @@ class Box(Spatial): Whether spatial sites should only be accepted if they occur in fissionable materials + .. deprecated:: 0.14.1 + Use the `constraints` argument when defining a source object instead. + """ def __init__( @@ -818,6 +825,10 @@ def only_fissionable(self): def only_fissionable(self, only_fissionable): cv.check_type('only fissionable', only_fissionable, bool) self._only_fissionable = only_fissionable + if only_fissionable: + warn("The 'only_fissionable' has been deprecated. Use the " + "'constraints' argument when defining a source instead.", + FutureWarning) def to_xml_element(self): """Return XML representation of the box distribution diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp index 8f34ea6b936..8dbd7d88706 100644 --- a/src/distribution_spatial.cpp +++ b/src/distribution_spatial.cpp @@ -281,10 +281,15 @@ void MeshSpatial::check_element_types() const } } +int32_t MeshSpatial::sample_element_index(uint64_t* seed) const +{ + return elem_idx_dist_.sample(seed); +} + std::pair MeshSpatial::sample_mesh(uint64_t* seed) const { // Sample the CDF defined in initialization above - int32_t elem_idx = elem_idx_dist_.sample(seed); + int32_t elem_idx = this->sample_element_index(seed); return {elem_idx, mesh()->sample_element(elem_idx, seed)}; } diff --git a/src/simulation.cpp b/src/simulation.cpp index 43d92060669..527d98db4f1 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -750,7 +750,7 @@ void free_memory_simulation() void transport_history_based_single_particle(Particle& p) { - while (true) { + while (p.alive()) { p.event_calculate_xs(); if (!p.alive()) break; @@ -761,8 +761,6 @@ void transport_history_based_single_particle(Particle& p) p.event_collide(); } p.event_revive_from_secondary(); - if (!p.alive()) - break; } p.event_death(); } diff --git a/src/source.cpp b/src/source.cpp index 5ddac7f4287..405de998c89 100644 --- a/src/source.cpp +++ b/src/source.cpp @@ -47,9 +47,23 @@ vector> external_sources; } //============================================================================== -// Source create implementation +// Source implementation //============================================================================== +Source::Source(pugi::xml_node node) +{ + // Check for source strength + if (check_for_node(node, "strength")) { + strength_ = std::stod(get_node_value(node, "strength")); + if (strength_ < 0.0) { + fatal_error("Source strength is negative."); + } + } + + // Check for additional defined constraints + read_constraints(node); +} + unique_ptr Source::create(pugi::xml_node node) { // if the source type is present, use it to determine the type @@ -79,6 +93,163 @@ unique_ptr Source::create(pugi::xml_node node) } } +void Source::read_constraints(pugi::xml_node node) +{ + // Check for constraints node. For backwards compatibility, if no constraints + // node is given, still try searching for domain constraints from top-level + // node. + pugi::xml_node constraints_node = node.child("constraints"); + if (constraints_node) { + node = constraints_node; + } + + // Check for domains to reject from + if (check_for_node(node, "domain_type")) { + std::string domain_type = get_node_value(node, "domain_type"); + if (domain_type == "cell") { + domain_type_ = DomainType::CELL; + } else if (domain_type == "material") { + domain_type_ = DomainType::MATERIAL; + } else if (domain_type == "universe") { + domain_type_ = DomainType::UNIVERSE; + } else { + fatal_error( + std::string("Unrecognized domain type for constraint: " + domain_type)); + } + + auto ids = get_node_array(node, "domain_ids"); + domain_ids_.insert(ids.begin(), ids.end()); + } + + if (check_for_node(node, "time_bounds")) { + auto ids = get_node_array(node, "time_bounds"); + if (ids.size() != 2) { + fatal_error("Time bounds must be represented by two numbers."); + } + time_bounds_ = std::make_pair(ids[0], ids[1]); + } + if (check_for_node(node, "energy_bounds")) { + auto ids = get_node_array(node, "energy_bounds"); + if (ids.size() != 2) { + fatal_error("Energy bounds must be represented by two numbers."); + } + energy_bounds_ = std::make_pair(ids[0], ids[1]); + } + + if (check_for_node(node, "fissionable")) { + only_fissionable_ = get_node_value_bool(node, "fissionable"); + } + + // Check for how to handle rejected particles + if (check_for_node(node, "rejection_strategy")) { + std::string rejection_strategy = get_node_value(node, "rejection_strategy"); + if (rejection_strategy == "kill") { + rejection_strategy_ = RejectionStrategy::KILL; + } else if (rejection_strategy == "resample") { + rejection_strategy_ = RejectionStrategy::RESAMPLE; + } else { + fatal_error(std::string( + "Unrecognized strategy source rejection: " + rejection_strategy)); + } + } +} + +SourceSite Source::sample_with_constraints(uint64_t* seed) const +{ + bool accepted = false; + static int n_reject = 0; + static int n_accept = 0; + SourceSite site; + + while (!accepted) { + // Sample a source site without considering constraints yet + site = this->sample(seed); + + if (constraints_applied()) { + accepted = true; + } else { + // Check whether sampled site satisfies constraints + accepted = satisfies_spatial_constraints(site.r) && + satisfies_energy_constraints(site.E) && + satisfies_time_constraints(site.time); + if (!accepted) { + ++n_reject; + if (n_reject >= EXTSRC_REJECT_THRESHOLD && + static_cast(n_accept) / n_reject <= + EXTSRC_REJECT_FRACTION) { + fatal_error("More than 95% of external source sites sampled were " + "rejected. Please check your source definition."); + } + + // For the "kill" strategy, accept particle but set weight to 0 so that + // it is terminated immediately + if (rejection_strategy_ == RejectionStrategy::KILL) { + accepted = true; + site.wgt = 0.0; + } + } + } + } + + // Increment number of accepted samples + ++n_accept; + + return site; +} + +bool Source::satisfies_energy_constraints(double E) const +{ + return E > energy_bounds_.first && E < energy_bounds_.second; +} + +bool Source::satisfies_time_constraints(double time) const +{ + return time > time_bounds_.first && time < time_bounds_.second; +} + +bool Source::satisfies_spatial_constraints(Position r) const +{ + GeometryState geom_state; + geom_state.r() = r; + + // Reject particle if it's not in the geometry at all + bool found = exhaustive_find_cell(geom_state); + if (!found) + return false; + + // Check the geometry state against specified domains + bool accepted = true; + if (!domain_ids_.empty()) { + if (domain_type_ == DomainType::MATERIAL) { + auto mat_index = geom_state.material(); + if (mat_index != MATERIAL_VOID) { + accepted = contains(domain_ids_, model::materials[mat_index]->id()); + } + } else { + for (int i = 0; i < geom_state.n_coord(); i++) { + auto id = (domain_type_ == DomainType::CELL) + ? model::cells[geom_state.coord(i).cell]->id_ + : model::universes[geom_state.coord(i).universe]->id_; + if ((accepted = contains(domain_ids_, id))) + break; + } + } + } + + // Check if spatial site is in fissionable material + if (accepted && only_fissionable_) { + // Determine material + auto mat_index = geom_state.material(); + if (mat_index == MATERIAL_VOID) { + accepted = false; + } else { + accepted = model::materials[mat_index]->fissionable(); + } + } + + return accepted; +} + //============================================================================== // IndependentSource implementation //============================================================================== @@ -89,7 +260,7 @@ IndependentSource::IndependentSource( energy_ {std::move(energy)}, time_ {std::move(time)} {} -IndependentSource::IndependentSource(pugi::xml_node node) +IndependentSource::IndependentSource(pugi::xml_node node) : Source(node) { // Check for particle type if (check_for_node(node, "particle")) { @@ -104,11 +275,6 @@ IndependentSource::IndependentSource(pugi::xml_node node) } } - // Check for source strength - if (check_for_node(node, "strength")) { - strength_ = std::stod(get_node_value(node, "strength")); - } - // Check for external source file if (check_for_node(node, "file")) { @@ -122,6 +288,15 @@ IndependentSource::IndependentSource(pugi::xml_node node) space_ = UPtrSpace {new SpatialPoint()}; } + // For backwards compatibility, check for only fissionable setting on box + // source + auto space_box = dynamic_cast(space_.get()); + if (space_box) { + if (!only_fissionable_) { + only_fissionable_ = space_box->only_fissionable(); + } + } + // Determine external source angular distribution if (check_for_node(node, "angle")) { angle_ = UnitSphereDistribution::create(node.child("angle")); @@ -148,24 +323,6 @@ IndependentSource::IndependentSource(pugi::xml_node node) double p[] {1.0}; time_ = UPtrDist {new Discrete {T, p, 1}}; } - - // Check for domains to reject from - if (check_for_node(node, "domain_type")) { - std::string domain_type = get_node_value(node, "domain_type"); - if (domain_type == "cell") { - domain_type_ = DomainType::CELL; - } else if (domain_type == "material") { - domain_type_ = DomainType::MATERIAL; - } else if (domain_type == "universe") { - domain_type_ = DomainType::UNIVERSE; - } else { - fatal_error(std::string( - "Unrecognized domain type for source rejection: " + domain_type)); - } - - auto ids = get_node_array(node, "domain_ids"); - domain_ids_.insert(ids.begin(), ids.end()); - } } } @@ -174,60 +331,21 @@ SourceSite IndependentSource::sample(uint64_t* seed) const SourceSite site; site.particle = particle_; - // Repeat sampling source location until a good site has been found - bool found = false; - int n_reject = 0; + // Repeat sampling source location until a good site has been accepted + bool accepted = false; + static int n_reject = 0; static int n_accept = 0; - while (!found) { - // Set particle type - Particle p; - p.type() = particle_; - p.u() = {0.0, 0.0, 1.0}; + while (!accepted) { // Sample spatial distribution - p.r() = space_->sample(seed); - - // Now search to see if location exists in geometry - found = exhaustive_find_cell(p); - - // Check if spatial site is in fissionable material - if (found) { - auto space_box = dynamic_cast(space_.get()); - if (space_box) { - if (space_box->only_fissionable()) { - // Determine material - auto mat_index = p.material(); - if (mat_index == MATERIAL_VOID) { - found = false; - } else { - found = model::materials[mat_index]->fissionable(); - } - } - } + site.r = space_->sample(seed); - // Rejection based on cells/materials/universes - if (!domain_ids_.empty()) { - found = false; - if (domain_type_ == DomainType::MATERIAL) { - auto mat_index = p.material(); - if (mat_index != MATERIAL_VOID) { - found = contains(domain_ids_, model::materials[mat_index]->id()); - } - } else { - for (int i = 0; i < p.n_coord(); i++) { - auto id = (domain_type_ == DomainType::CELL) - ? model::cells[p.coord(i).cell]->id_ - : model::universes[p.coord(i).universe]->id_; - if ((found = contains(domain_ids_, id))) - break; - } - } - } - } + // Check if sampled position satisfies spatial constraints + accepted = satisfies_spatial_constraints(site.r); // Check for rejection - if (!found) { + if (!accepted) { ++n_reject; if (n_reject >= EXTSRC_REJECT_THRESHOLD && static_cast(n_accept) / n_reject <= EXTSRC_REJECT_FRACTION) { @@ -236,8 +354,6 @@ SourceSite IndependentSource::sample(uint64_t* seed) const "definition."); } } - - site.r = p.r(); } // Sample angle @@ -261,7 +377,8 @@ SourceSite IndependentSource::sample(uint64_t* seed) const site.E = energy_->sample(seed); // Resample if energy falls above maximum particle energy - if (site.E < data::energy_max[p]) + if (site.E < data::energy_max[p] and + (satisfies_energy_constraints(site.E))) break; n_reject++; @@ -287,7 +404,8 @@ SourceSite IndependentSource::sample(uint64_t* seed) const //============================================================================== // FileSource implementation //============================================================================== -FileSource::FileSource(pugi::xml_node node) + +FileSource::FileSource(pugi::xml_node node) : Source(node) { auto path = get_node_value(node, "file", false, true); if (ends_with(path, ".mcpl") || ends_with(path, ".mcpl.gz")) { @@ -332,6 +450,7 @@ void FileSource::load_sites_from_file(const std::string& path) SourceSite FileSource::sample(uint64_t* seed) const { + // Sample a particle randomly from list size_t i_site = sites_.size() * prn(seed); return sites_[i_site]; } @@ -339,7 +458,8 @@ SourceSite FileSource::sample(uint64_t* seed) const //============================================================================== // CompiledSourceWrapper implementation //============================================================================== -CompiledSourceWrapper::CompiledSourceWrapper(pugi::xml_node node) + +CompiledSourceWrapper::CompiledSourceWrapper(pugi::xml_node node) : Source(node) { // Get shared library path and parameters auto path = get_node_value(node, "library", false, true); @@ -403,7 +523,7 @@ CompiledSourceWrapper::~CompiledSourceWrapper() // MeshSource implementation //============================================================================== -MeshSource::MeshSource(pugi::xml_node node) +MeshSource::MeshSource(pugi::xml_node node) : Source(node) { int32_t mesh_id = stoi(get_node_value(node, "mesh")); int32_t mesh_idx = model::mesh_map.at(mesh_id); @@ -430,15 +550,27 @@ MeshSource::MeshSource(pugi::xml_node node) SourceSite MeshSource::sample(uint64_t* seed) const { - // sample location and element from mesh - auto mesh_location = space_->sample_mesh(seed); + // Sample the CDF defined in initialization above + int32_t element = space_->sample_element_index(seed); - // Sample source for the chosen element - int32_t element = mesh_location.first; - SourceSite site = source(element)->sample(seed); + // Sample position and apply rejection on spatial domains + Position r; + do { + r = space_->mesh()->sample_element(element, seed); + } while (!this->satisfies_spatial_constraints(r)); - // Replace spatial position with the one already sampled - site.r = mesh_location.second; + SourceSite site; + while (true) { + // Sample source for the chosen element and replace the position + site = source(element)->sample_with_constraints(seed); + site.r = r; + + // Apply other rejections + if (satisfies_energy_constraints(site.E) && + satisfies_time_constraints(site.time)) { + break; + } + } return site; } @@ -493,7 +625,7 @@ SourceSite sample_external_source(uint64_t* seed) } // Sample source site from i-th source distribution - SourceSite site {model::external_sources[i]->sample(seed)}; + SourceSite site {model::external_sources[i]->sample_with_constraints(seed)}; // If running in MG, convert site.E to group if (!settings::run_CE) { diff --git a/tests/regression_tests/asymmetric_lattice/inputs_true.dat b/tests/regression_tests/asymmetric_lattice/inputs_true.dat index 3977220da61..302ef24c200 100644 --- a/tests/regression_tests/asymmetric_lattice/inputs_true.dat +++ b/tests/regression_tests/asymmetric_lattice/inputs_true.dat @@ -209,9 +209,12 @@ 10 5 - + -32 -32 0 32 32 32 + + true + diff --git a/tests/regression_tests/asymmetric_lattice/test.py b/tests/regression_tests/asymmetric_lattice/test.py index 0d17a6c98eb..fadf272ffd3 100644 --- a/tests/regression_tests/asymmetric_lattice/test.py +++ b/tests/regression_tests/asymmetric_lattice/test.py @@ -54,8 +54,10 @@ def __init__(self, *args, **kwargs): self._model.tallies.append(tally) # Specify summary output and correct source sampling box - self._model.settings.source = openmc.IndependentSource(space=openmc.stats.Box( - [-32, -32, 0], [32, 32, 32], only_fissionable = True)) + self._model.settings.source = openmc.IndependentSource( + space=openmc.stats.Box([-32, -32, 0], [32, 32, 32]), + constraints={'fissionable': True} + ) def _get_results(self, hash_output=True): """Digest info in statepoint and summary and return as a string.""" diff --git a/tests/regression_tests/mg_temperature_multi/inputs_true.dat b/tests/regression_tests/mg_temperature_multi/inputs_true.dat index b2ad5063c2b..0c452cd520e 100644 --- a/tests/regression_tests/mg_temperature_multi/inputs_true.dat +++ b/tests/regression_tests/mg_temperature_multi/inputs_true.dat @@ -28,9 +28,12 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + + true + multi-group diff --git a/tests/regression_tests/mg_temperature_multi/test.py b/tests/regression_tests/mg_temperature_multi/test.py index 404c8c0cd05..3117e29ba03 100755 --- a/tests/regression_tests/mg_temperature_multi/test.py +++ b/tests/regression_tests/mg_temperature_multi/test.py @@ -133,8 +133,9 @@ def test_mg_temperature_multi(): # Create an initial uniform spatial source distribution over fissionable zones lower_left = (-pitch/2, -pitch/2, -1) upper_right = (pitch/2, pitch/2, 1) - uniform_dist = openmc.stats.Box(lower_left, upper_right, only_fissionable=True) - settings.source = openmc.IndependentSource(space=uniform_dist) + uniform_dist = openmc.stats.Box(lower_left, upper_right) + settings.source = openmc.IndependentSource( + space=uniform_dist, constraints={'fissionable': True}) ############################################################################### # Define tallies diff --git a/tests/regression_tests/mgxs_library_ce_to_mg/inputs_true.dat b/tests/regression_tests/mgxs_library_ce_to_mg/inputs_true.dat index d9a3c52c840..6c071263f3d 100644 --- a/tests/regression_tests/mgxs_library_ce_to_mg/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_ce_to_mg/inputs_true.dat @@ -42,9 +42,12 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + + true + diff --git a/tests/regression_tests/mgxs_library_ce_to_mg_nuclides/inputs_true.dat b/tests/regression_tests/mgxs_library_ce_to_mg_nuclides/inputs_true.dat index f84adc32967..ad3d50974b0 100644 --- a/tests/regression_tests/mgxs_library_ce_to_mg_nuclides/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_ce_to_mg_nuclides/inputs_true.dat @@ -42,9 +42,12 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + + true + diff --git a/tests/regression_tests/mgxs_library_condense/inputs_true.dat b/tests/regression_tests/mgxs_library_condense/inputs_true.dat index 93519cd2a2a..8606e9fdda3 100644 --- a/tests/regression_tests/mgxs_library_condense/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_condense/inputs_true.dat @@ -42,9 +42,12 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + + true + diff --git a/tests/regression_tests/mgxs_library_correction/inputs_true.dat b/tests/regression_tests/mgxs_library_correction/inputs_true.dat index 80cbc4bb56a..ea2d1b714a2 100644 --- a/tests/regression_tests/mgxs_library_correction/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_correction/inputs_true.dat @@ -42,9 +42,12 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + + true + diff --git a/tests/regression_tests/mgxs_library_distribcell/inputs_true.dat b/tests/regression_tests/mgxs_library_distribcell/inputs_true.dat index fc2c06d1cfa..f1a60533491 100644 --- a/tests/regression_tests/mgxs_library_distribcell/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_distribcell/inputs_true.dat @@ -69,9 +69,12 @@ 10 5 - + -10.71 -10.71 -1 10.71 10.71 1 + + true + diff --git a/tests/regression_tests/mgxs_library_hdf5/inputs_true.dat b/tests/regression_tests/mgxs_library_hdf5/inputs_true.dat index 93519cd2a2a..8606e9fdda3 100644 --- a/tests/regression_tests/mgxs_library_hdf5/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_hdf5/inputs_true.dat @@ -42,9 +42,12 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + + true + diff --git a/tests/regression_tests/mgxs_library_histogram/inputs_true.dat b/tests/regression_tests/mgxs_library_histogram/inputs_true.dat index 178a44464e9..013fcd0fa5f 100644 --- a/tests/regression_tests/mgxs_library_histogram/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_histogram/inputs_true.dat @@ -42,9 +42,12 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + + true + diff --git a/tests/regression_tests/mgxs_library_no_nuclides/inputs_true.dat b/tests/regression_tests/mgxs_library_no_nuclides/inputs_true.dat index 0003413acb1..5cad36b7ea3 100644 --- a/tests/regression_tests/mgxs_library_no_nuclides/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_no_nuclides/inputs_true.dat @@ -42,9 +42,12 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + + true + diff --git a/tests/regression_tests/mgxs_library_nuclides/inputs_true.dat b/tests/regression_tests/mgxs_library_nuclides/inputs_true.dat index 97e227177e0..9b3a22510b1 100644 --- a/tests/regression_tests/mgxs_library_nuclides/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_nuclides/inputs_true.dat @@ -42,9 +42,12 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + + true + diff --git a/tests/regression_tests/mgxs_library_specific_nuclides/inputs_true.dat b/tests/regression_tests/mgxs_library_specific_nuclides/inputs_true.dat index b73ed179e5e..1372741c947 100644 --- a/tests/regression_tests/mgxs_library_specific_nuclides/inputs_true.dat +++ b/tests/regression_tests/mgxs_library_specific_nuclides/inputs_true.dat @@ -42,9 +42,12 @@ 10 5 - + -0.63 -0.63 -1 0.63 0.63 1 + + true + diff --git a/tests/regression_tests/surface_tally/inputs_true.dat b/tests/regression_tests/surface_tally/inputs_true.dat index abd01266d48..192c2e89f40 100644 --- a/tests/regression_tests/surface_tally/inputs_true.dat +++ b/tests/regression_tests/surface_tally/inputs_true.dat @@ -30,9 +30,12 @@ 10 0 - + -0.62992 -0.62992 -1 0.62992 0.62992 1 + + true + diff --git a/tests/regression_tests/surface_tally/test.py b/tests/regression_tests/surface_tally/test.py index 8968db9b6fb..e496ac0f65c 100644 --- a/tests/regression_tests/surface_tally/test.py +++ b/tests/regression_tests/surface_tally/test.py @@ -72,9 +72,9 @@ def __init__(self, *args, **kwargs): # Create an initial uniform spatial source distribution bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] - uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:],\ - only_fissionable=True) - settings_file.source = openmc.IndependentSource(space=uniform_dist) + uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:],) + settings_file.source = openmc.IndependentSource( + space=uniform_dist, constraints={'fissionable': True}) self._model.settings = settings_file # Tallies file diff --git a/tests/unit_tests/test_lib.py b/tests/unit_tests/test_lib.py index 1310a6a7fa0..009a9096831 100644 --- a/tests/unit_tests/test_lib.py +++ b/tests/unit_tests/test_lib.py @@ -67,8 +67,10 @@ def uo2_trigger_model(): model.settings.batches = 10 model.settings.inactive = 5 model.settings.particles = 100 - model.settings.source = openmc.IndependentSource(space=openmc.stats.Box( - [-0.5, -0.5, -1], [0.5, 0.5, 1], only_fissionable=True)) + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Box([-0.5, -0.5, -1], [0.5, 0.5, 1]), + constraints={'fissionable': True}, + ) model.settings.verbosity = 1 model.settings.keff_trigger = {'type': 'std_dev', 'threshold': 0.001} model.settings.trigger_active = True diff --git a/tests/unit_tests/test_model.py b/tests/unit_tests/test_model.py index 16fa18d453c..404235bef1f 100644 --- a/tests/unit_tests/test_model.py +++ b/tests/unit_tests/test_model.py @@ -57,9 +57,9 @@ def pin_model_attributes(): # Create a uniform spatial source distribution over fissionable zones bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1] - uniform_dist = openmc.stats.Box( - bounds[:3], bounds[3:], only_fissionable=True) - settings.source = openmc.IndependentSource(space=uniform_dist) + uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:]) + settings.source = openmc.IndependentSource( + space=uniform_dist, constraints={'fissionable': True}) entropy_mesh = openmc.RegularMesh() entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50] diff --git a/tests/unit_tests/test_source.py b/tests/unit_tests/test_source.py index 4783ff8f11d..9a19f6f24dd 100644 --- a/tests/unit_tests/test_source.py +++ b/tests/unit_tests/test_source.py @@ -7,6 +7,8 @@ import pytest from pytest import approx +from tests.regression_tests import config + def test_source(): space = openmc.stats.Point() @@ -100,7 +102,8 @@ def test_source_xml_roundtrip(): assert new_src.strength == approx(src.strength) -def test_rejection(run_in_tmpdir): +@pytest.fixture +def sphere_box_model(): # Model with two spheres inside a box mat = openmc.Material() mat.add_nuclide('H1', 1.0) @@ -119,24 +122,108 @@ def test_rejection(run_in_tmpdir): model.settings.batches = 10 model.settings.run_mode = 'fixed source' + return model, cell1, cell2, cell3 + + +def test_constraints_independent(sphere_box_model, run_in_tmpdir): + model, cell1, cell2, cell3 = sphere_box_model + # Set up a box source with rejection on the spherical cell - space = openmc.stats.Box(*cell3.bounding_box) - model.settings.source = openmc.IndependentSource(space=space, domains=[cell1, cell2]) + space = openmc.stats.Box((-4., -1., -1.), (4., 1., 1.)) + model.settings.source = openmc.IndependentSource( + space=space, constraints={'domains': [cell1, cell2]} + ) + + # Load up model via openmc.lib and sample source + model.export_to_model_xml() + openmc.lib.init() + particles = openmc.lib.sample_external_source(1000) + + # Make sure that all sampled sources are within one of the spheres + for p in particles: + assert p.r in (cell1.region | cell2.region) + assert p.r not in cell3.region + + openmc.lib.finalize() + + +def test_constraints_mesh(sphere_box_model, run_in_tmpdir): + model, cell1, cell2, cell3 = sphere_box_model + + bbox = cell3.bounding_box + mesh = openmc.RegularMesh() + mesh.lower_left = bbox.lower_left + mesh.upper_right = bbox.upper_right + mesh.dimension = (2, 1, 1) + + left_source = openmc.IndependentSource() + right_source = openmc.IndependentSource() + model.settings.source = openmc.MeshSource( + mesh, [left_source, right_source], constraints={'domains': [cell1, cell2]} + ) # Load up model via openmc.lib and sample source - model.export_to_xml() + model.export_to_model_xml() openmc.lib.init() particles = openmc.lib.sample_external_source(1000) # Make sure that all sampled sources are within one of the spheres - joint_region = cell1.region | cell2.region for p in particles: - assert p.r in joint_region - assert p.r not in non_source_region + assert p.r in (cell1.region | cell2.region) + assert p.r not in cell3.region openmc.lib.finalize() +def test_constraints_file(sphere_box_model, run_in_tmpdir): + model = sphere_box_model[0] + + # Create source file with randomly sampled source sites + rng = np.random.default_rng() + energy = rng.uniform(0., 1e6, 10_000) + time = rng.uniform(0., 1., 10_000) + particles = [openmc.SourceParticle(E=e, time=t) for e, t in zip(energy, time)] + openmc.write_source_file(particles, 'uniform_source.h5') + + # Use source file + model.settings.source = openmc.FileSource( + 'uniform_source.h5', + constraints={ + 'time_bounds': [0.25, 0.75], + 'energy_bounds': [500.e3, 1.0e6], + } + ) + + # Load up model via openmc.lib and sample source + model.export_to_model_xml() + openmc.lib.init() + particles = openmc.lib.sample_external_source(1000) + + # Make sure that all sampled sources are within energy/time bounds + for p in particles: + assert 0.25 <= p.time <= 0.75 + assert 500.e3 <= p.E <= 1.0e6 + + openmc.lib.finalize() + + +@pytest.mark.skipif(config['mpi'], reason='Not compatible with MPI') +def test_rejection_limit(sphere_box_model, run_in_tmpdir): + model, cell1 = sphere_box_model[:2] + + # Define a point source that will get rejected 100% of the time + model.settings.source = openmc.IndependentSource( + space=openmc.stats.Point((-3., 0., 0.)), + constraints={'domains': [cell1]} + ) + + # Confirm that OpenMC doesn't run in an infinite loop. Note that this may + # work when running with MPI since it won't necessarily capture the error + # message correctly + with pytest.raises(RuntimeError, match="rejected"): + model.run(openmc_exec=config['exe']) + + def test_exceptions(): with pytest.raises(AttributeError, match=r'Please use the FileSource class'): @@ -156,4 +243,4 @@ def test_exceptions(): with pytest.raises(AttributeError, match=r'has no attribute \'frisbee\''): s = openmc.IndependentSource() - s.frisbee \ No newline at end of file + s.frisbee