Skip to content

Commit 1702b45

Browse files
Erik B Knudsenchurch89paulromano
authored
Restricted file source (#2916)
Co-authored-by: church89 <l.chierici89@gmail.com> Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
1 parent a7d6939 commit 1702b45

32 files changed

Lines changed: 823 additions & 240 deletions

File tree

‎docs/source/io_formats/settings.rst‎

Lines changed: 36 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -534,8 +534,9 @@ attributes/sub-elements:
534534
*Default*: 1.0
535535

536536
:type:
537-
Indicator of source type. One of ``independent``, ``file``, ``compiled``, or ``mesh``.
538-
The type of the source will be determined by this attribute if it is present.
537+
Indicator of source type. One of ``independent``, ``file``, ``compiled``, or
538+
``mesh``. The type of the source will be determined by this attribute if it
539+
is present.
539540

540541
:particle:
541542
The source particle type, either ``neutron`` or ``photon``.
@@ -716,6 +717,39 @@ attributes/sub-elements:
716717
mesh element and follows the format for :ref:`source_element`. The number of
717718
``<source>`` sub-elements should correspond to the number of mesh elements.
718719

720+
:constraints:
721+
This sub-element indicates the presence of constraints on sampled source
722+
sites (see :ref:`usersguide_source_constraints` for details). It may have
723+
the following sub-elements:
724+
725+
:domain_ids:
726+
The unique IDs of domains for which source sites must be within.
727+
728+
*Default*: None
729+
730+
:domain_type:
731+
The type of each domain for source rejection ("cell", "material", or
732+
"universe").
733+
734+
*Default*: None
735+
736+
:fissionable:
737+
A boolean indicating whether source sites must be sampled within a
738+
material that is fissionable in order to be accepted.
739+
740+
:time_bounds:
741+
A pair of times in [s] indicating the lower and upper bound for a time
742+
interval that source particles must be within.
743+
744+
:energy_bounds:
745+
A pair of energies in [eV] indicating the lower and upper bound for an
746+
energy interval that source particles must be within.
747+
748+
:rejection_strategy:
749+
Either "resample", indicating that source sites should be resampled when
750+
one is rejected, or "kill", indicating that a rejected source site is
751+
assigned zero weight.
752+
719753
.. _univariate:
720754

721755
Univariate Probability Distributions

‎docs/source/usersguide/settings.rst‎

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -420,6 +420,54 @@ string, which gets passed down to the ``openmc_create_source()`` function::
420420

421421
settings.source = openmc.CompiledSource('libsource.so', '3.5e6')
422422

423+
.. _usersguide_source_constraints:
424+
425+
Source Constraints
426+
------------------
427+
428+
All source classes in OpenMC have the ability to apply a set of "constraints"
429+
that limit which sampled source sites are actually used for transport. The most
430+
common use case is to sample source sites over some simple spatial distribution
431+
(e.g., uniform over a box) and then only accept those that appear in a given
432+
cell or material. This can be done with a domain constraint, which can be
433+
specified as follows::
434+
435+
source_cell = openmc.Cell(...)
436+
...
437+
438+
spatial_dist = openmc.stats.Box((-10., -10., -10.), (10., 10., 10.))
439+
source = openmc.IndependentSource(
440+
space=spatial_dist,
441+
constraints={'domains': [source_cell]}
442+
)
443+
444+
For k-eigenvalue problems, a convenient constraint is available that limits
445+
source sites to those sampled in a fissionable material::
446+
447+
source = openmc.IndependentSource(
448+
space=spatial_dist, constraints={'fissionable': True}
449+
)
450+
451+
Constraints can also be placed on a range of energies or times::
452+
453+
# Only use source sites between 500 keV and 1 MeV and with times under 1 sec
454+
source = openmc.FileSource(
455+
'source.h5',
456+
constraints={'energy_bounds': [500.0e3, 1.0e6], 'time_bounds': [0.0, 1.0]}
457+
)
458+
459+
Normally, when a source site is rejected, a new one will be resampled until one
460+
is found that meets the constraints. However, the rejection strategy can be
461+
changed so that a rejected site will just not be simulated by specifying::
462+
463+
source = openmc.IndependentSource(
464+
space=spatial_dist,
465+
constraints={'domains': [cell], 'rejection_strategy': 'kill'}
466+
)
467+
468+
In this case, the actual number of particles simulated may be less than what you
469+
specified in :attr:`Settings.particles`.
470+
423471
.. _usersguide_entropy:
424472

425473
---------------

‎examples/assembly/assembly.py‎

Lines changed: 4 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -112,11 +112,10 @@ def assembly_model():
112112
model.settings.batches = 150
113113
model.settings.inactive = 50
114114
model.settings.particles = 1000
115-
model.settings.source = openmc.IndependentSource(space=openmc.stats.Box(
116-
(-pitch/2, -pitch/2, -1),
117-
(pitch/2, pitch/2, 1),
118-
only_fissionable=True
119-
))
115+
model.settings.source = openmc.IndependentSource(
116+
space=openmc.stats.Box((-pitch/2, -pitch/2, -1), (pitch/2, pitch/2, 1)),
117+
constraints={'fissionable': True}
118+
)
120119

121120
# NOTE: We never actually created a Materials object. When you export/run
122121
# using the Model object, if no materials were assigned it will look through

‎examples/pincell_depletion/restart_depletion.py‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -26,8 +26,9 @@
2626

2727
# Create an initial uniform spatial source distribution over fissionable zones
2828
bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1]
29-
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
30-
settings.source = openmc.IndependentSource(space=uniform_dist)
29+
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:])
30+
settings.source = openmc.IndependentSource(
31+
space=uniform_dist, constraints={'fissionable': True})
3132

3233
entropy_mesh = openmc.RegularMesh()
3334
entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50]

‎examples/pincell_depletion/run_depletion.py‎

Lines changed: 3 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -71,8 +71,9 @@
7171

7272
# Create an initial uniform spatial source distribution over fissionable zones
7373
bounds = [-0.62992, -0.62992, -1, 0.62992, 0.62992, 1]
74-
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
75-
settings.source = openmc.IndependentSource(space=uniform_dist)
74+
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:])
75+
settings.source = openmc.IndependentSource(
76+
space=uniform_dist, constraints={'fissionable': True})
7677

7778
entropy_mesh = openmc.RegularMesh()
7879
entropy_mesh.lower_left = [-0.39218, -0.39218, -1.e50]

‎include/openmc/distribution_spatial.h‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -116,6 +116,11 @@ class MeshSpatial : public SpatialDistribution {
116116
//! \return Sampled element index and position within that element
117117
std::pair<int32_t, Position> sample_mesh(uint64_t* seed) const;
118118

119+
//! Sample a mesh element
120+
//! \param seed Pseudorandom number seed pointer
121+
//! \return Sampled element index
122+
int32_t sample_element_index(uint64_t* seed) const;
123+
119124
//! For unstructured meshes, ensure that elements are all linear tetrahedra
120125
void check_element_types() const;
121126

‎include/openmc/source.h‎

Lines changed: 74 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,7 @@
44
#ifndef OPENMC_SOURCE_H
55
#define OPENMC_SOURCE_H
66

7+
#include <limits>
78
#include <unordered_set>
89

910
#include "pugixml.hpp"
@@ -39,19 +40,72 @@ extern vector<unique_ptr<Source>> external_sources;
3940

4041
//==============================================================================
4142
//! Abstract source interface
43+
//
44+
//! The Source class provides the interface that must be implemented by derived
45+
//! classes, namely the sample() method that returns a sampled source site. From
46+
//! this base class, source rejection is handled within the
47+
//! sample_with_constraints() method. However, note that some classes directly
48+
//! check for constraints for efficiency reasons (like IndependentSource), in
49+
//! which case the constraints_applied() method indicates that constraints
50+
//! should not be checked a second time from the base class.
4251
//==============================================================================
4352

4453
class Source {
4554
public:
55+
// Constructors, destructors
56+
Source() = default;
57+
explicit Source(pugi::xml_node node);
4658
virtual ~Source() = default;
4759

48-
// Methods that must be implemented
49-
virtual SourceSite sample(uint64_t* seed) const = 0;
50-
5160
// Methods that can be overridden
52-
virtual double strength() const { return 1.0; }
61+
virtual double strength() const { return strength_; }
62+
63+
//! Sample a source site and apply constraints
64+
//
65+
//! \param[inout] seed Pseudorandom seed pointer
66+
//! \return Sampled site
67+
SourceSite sample_with_constraints(uint64_t* seed) const;
68+
69+
//! Sample a source site (without applying constraints)
70+
//
71+
//! Sample from the external source distribution
72+
//! \param[inout] seed Pseudorandom seed pointer
73+
//! \return Sampled site
74+
virtual SourceSite sample(uint64_t* seed) const = 0;
5375

5476
static unique_ptr<Source> create(pugi::xml_node node);
77+
78+
protected:
79+
// Domain types
80+
enum class DomainType { UNIVERSE, MATERIAL, CELL };
81+
82+
// Strategy used for rejecting sites when constraints are applied. KILL means
83+
// that sites are always accepted but if they don't satisfy constraints, they
84+
// are given weight 0. RESAMPLE means that a new source site will be sampled
85+
// until constraints are met.
86+
enum class RejectionStrategy { KILL, RESAMPLE };
87+
88+
// Indicates whether derived class already handles constraints
89+
virtual bool constraints_applied() const { return false; }
90+
91+
// Methods for constraints
92+
void read_constraints(pugi::xml_node node);
93+
bool satisfies_spatial_constraints(Position r) const;
94+
bool satisfies_energy_constraints(double E) const;
95+
bool satisfies_time_constraints(double time) const;
96+
97+
// Data members
98+
double strength_ {1.0}; //!< Source strength
99+
std::unordered_set<int32_t> domain_ids_; //!< Domains to reject from
100+
DomainType domain_type_; //!< Domain type for rejection
101+
std::pair<double, double> time_bounds_ {-std::numeric_limits<double>::max(),
102+
std::numeric_limits<double>::max()}; //!< time limits
103+
std::pair<double, double> energy_bounds_ {
104+
0, std::numeric_limits<double>::max()}; //!< energy limits
105+
bool only_fissionable_ {
106+
false}; //!< Whether site must be in fissionable material
107+
RejectionStrategy rejection_strategy_ {
108+
RejectionStrategy::RESAMPLE}; //!< Procedure for rejecting
55109
};
56110

57111
//==============================================================================
@@ -73,27 +127,24 @@ class IndependentSource : public Source {
73127

74128
// Properties
75129
ParticleType particle_type() const { return particle_; }
76-
double strength() const override { return strength_; }
77130

78131
// Make observing pointers available
79132
SpatialDistribution* space() const { return space_.get(); }
80133
UnitSphereDistribution* angle() const { return angle_.get(); }
81134
Distribution* energy() const { return energy_.get(); }
82135
Distribution* time() const { return time_.get(); }
83136

84-
private:
85-
// Domain types
86-
enum class DomainType { UNIVERSE, MATERIAL, CELL };
137+
protected:
138+
// Indicates whether derived class already handles constraints
139+
bool constraints_applied() const override { return true; }
87140

141+
private:
88142
// Data members
89143
ParticleType particle_ {ParticleType::neutron}; //!< Type of particle emitted
90-
double strength_ {1.0}; //!< Source strength
91144
UPtrSpace space_; //!< Spatial distribution
92145
UPtrAngle angle_; //!< Angular distribution
93146
UPtrDist energy_; //!< Energy distribution
94147
UPtrDist time_; //!< Time distribution
95-
DomainType domain_type_; //!< Domain type for rejection
96-
std::unordered_set<int32_t> domain_ids_; //!< Domains to reject from
97148
};
98149

99150
//==============================================================================
@@ -107,9 +158,12 @@ class FileSource : public Source {
107158
explicit FileSource(const std::string& path);
108159

109160
// Methods
110-
SourceSite sample(uint64_t* seed) const override;
111161
void load_sites_from_file(
112162
const std::string& path); //!< Load source sites from file
163+
164+
protected:
165+
SourceSite sample(uint64_t* seed) const override;
166+
113167
private:
114168
vector<SourceSite> sites_; //!< Source sites from a file
115169
};
@@ -124,16 +178,17 @@ class CompiledSourceWrapper : public Source {
124178
CompiledSourceWrapper(pugi::xml_node node);
125179
~CompiledSourceWrapper();
126180

181+
double strength() const override { return compiled_source_->strength(); }
182+
183+
void setup(const std::string& path, const std::string& parameters);
184+
185+
protected:
127186
// Defer implementation to custom source library
128187
SourceSite sample(uint64_t* seed) const override
129188
{
130189
return compiled_source_->sample(seed);
131190
}
132191

133-
double strength() const override { return compiled_source_->strength(); }
134-
135-
void setup(const std::string& path, const std::string& parameters);
136-
137192
private:
138193
void* shared_library_; //!< library from dlopen
139194
unique_ptr<Source> compiled_source_;
@@ -164,6 +219,9 @@ class MeshSource : public Source {
164219
return sources_.size() == 1 ? sources_[0] : sources_[i];
165220
}
166221

222+
protected:
223+
bool constraints_applied() const override { return true; }
224+
167225
private:
168226
// Data members
169227
unique_ptr<MeshSpatial> space_; //!< Mesh spatial

‎openmc/examples.py‎

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -76,8 +76,10 @@ def pwr_pin_cell():
7676
model.settings.batches = 10
7777
model.settings.inactive = 5
7878
model.settings.particles = 100
79-
model.settings.source = openmc.IndependentSource(space=openmc.stats.Box(
80-
[-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1], only_fissionable=True))
79+
model.settings.source = openmc.IndependentSource(
80+
space=openmc.stats.Box([-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1]),
81+
constraints={'fissionable': True}
82+
)
8183

8284
plot = openmc.Plot.from_geometry(model.geometry)
8385
plot.pixels = (300, 300)
@@ -527,8 +529,10 @@ def pwr_assembly():
527529
model.settings.batches = 10
528530
model.settings.inactive = 5
529531
model.settings.particles = 100
530-
model.settings.source = openmc.IndependentSource(space=openmc.stats.Box(
531-
[-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1], only_fissionable=True))
532+
model.settings.source = openmc.IndependentSource(
533+
space=openmc.stats.Box([-pitch/2, -pitch/2, -1], [pitch/2, pitch/2, 1]),
534+
constraints={'fissionable': True}
535+
)
532536

533537
plot = openmc.Plot()
534538
plot.origin = (0.0, 0.0, 0)

0 commit comments

Comments
 (0)