Restricted file source (openmc-dev#2916)

Co-authored-by: church89 <l.chierici89@gmail.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

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
     ``<source>`` 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

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:
 
 ---------------

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

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]

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]

include/openmc/distribution_spatial.h

@@ -116,6 +116,11 @@ class MeshSpatial : public SpatialDistribution {
   //! \return Sampled element index and position within that element
   std::pair<int32_t, Position> 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;
 

include/openmc/source.h

@@ -4,6 +4,7 @@
 #ifndef OPENMC_SOURCE_H
 #define OPENMC_SOURCE_H
 
+#include <limits>
 #include <unordered_set>
 
 #include "pugixml.hpp"
@@ -39,19 +40,72 @@ extern vector<unique_ptr<Source>> 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<Source> 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<int32_t> domain_ids_; //!< Domains to reject from
+  DomainType domain_type_;                 //!< Domain type for rejection
+  std::pair<double, double> time_bounds_ {-std::numeric_limits<double>::max(),
+    std::numeric_limits<double>::max()}; //!< time limits
+  std::pair<double, double> energy_bounds_ {
+    0, std::numeric_limits<double>::max()}; //!< energy limits
+  bool only_fissionable_ {
+    false}; //!< Whether site must be in fissionable material
+  RejectionStrategy rejection_strategy_ {
+    RejectionStrategy::RESAMPLE}; //!< Procedure for rejecting
 };
 
 //==============================================================================
@@ -73,27 +127,24 @@ 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(); }
   UnitSphereDistribution* angle() const { return angle_.get(); }
   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<int32_t> 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<SourceSite> 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<Source> 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<MeshSpatial> space_;           //!< Mesh spatial

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)