crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-21 14:35:27 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Source biasing capabilities (#3460)

Browse source 
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
Jack Fletcher 2026-01-12 10:51:12 -05:00 • committed by GitHub
parent 37e2feb34b
commit 0486e433d2
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
24 changed files with 2003 additions and 252 deletions
Download patch file Download diff file Expand all files Collapse all files

34
docs/source/io_formats/settings.rst
View file

@ -844,13 +844,18 @@ attributes/sub-elements:
relative source strength of each mesh element or each point in the cloud.
:volume_normalized:
For "mesh" spatial distrubtions, this optional boolean element specifies
For "mesh" spatial distributions, this optional boolean element specifies
whether the vector of relative strengths should be multiplied by the mesh
element volume. This is most common if the strengths represent a source
per unit volume.
*Default*: false
:bias:
For "mesh" and "cloud" spatial distributions, this optional element
specifies floating point values corresponding to alternative probabilities
for each value/component to use for biased sampling.
:angle:
An element specifying the angular distribution of source sites. This element
has the following attributes:
@ -883,6 +888,10 @@ attributes/sub-elements:
are those of a univariate probability distribution (see the description in
:ref:`univariate`).
:bias:
For "isotropic" angular distributions, this optional element specifies a
"mu-phi" angular distribution used for biased sampling.
:energy:
An element specifying the energy distribution of source sites. The necessary
sub-elements/attributes are those of a univariate probability distribution
@ -906,6 +915,10 @@ 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.
.. note:: Biased sampling can be applied to the spatial and energy distributions
of a source by using the ``<bias>`` sub-element (see
:ref:`univariate` for details on how to specify bias distributions).
:constraints:
This sub-element indicates the presence of constraints on sampled source
sites (see :ref:`usersguide_source_constraints` for details). It may have
@ -998,13 +1011,26 @@ variable and whose sub-elements/attributes are as follows:
*Default*: histogram
:pair:
For a "mixture" distribution, this element provides a distribution and its corresponding probability.
For a "mixture" distribution, this element provides a distribution and its
corresponding probability.
:probability:
An attribute or ``pair`` that provides the probability of a univariate distribution within a "mixture" distribution.
An attribute or ``pair`` that provides the probability of a univariate
distribution within a "mixture" distribution.
:dist:
This sub-element of a ``pair`` element provides information on the corresponding univariate distribution.
This sub-element of a ``pair`` element provides information on the
corresponding univariate distribution.
:bias:
This optional element specifies a biased distribution for importance sampling.
For continuous distributions, the ``bias`` element should contain another
univariate distribution with the same support (interval) as the parent
distribution. For discrete distributions, the ``bias`` element should contain
floating point values corresponding to alternative probabilities for each
value/component to be used for biased sampling.
*Default*: None
---------------------------------------
``<source_rejection_fraction>`` Element

82
docs/source/methods/variance_reduction.rst
View file

@ -22,12 +22,14 @@ not experience a single scoring event, even after billions of analog histories.
Variance reduction techniques aim to either flatten the global uncertainty
distribution, such that all regions of phase space have a fairly similar
uncertainty, or to reduce the uncertainty in specific locations (such as a
detector). There are two strategies available in OpenMC for variance reduction:
the Monte Carlo MAGIC method and the FW-CADIS method. Both strategies work by
developing a weight window mesh that can be utilized by subsequent Monte Carlo
solves to split particles heading towards areas of lower flux densities while
terminating particles in higher flux regions---all while maintaining a fair
game.
detector). There are three strategies available in OpenMC for variance
reduction: weight windows generated via the MAGIC method or the FW-CADIS method,
and source biasing. Both weight windowing strategies work by developing a mesh
that can be utilized by subsequent Monte Carlo solves to split particles heading
towards areas of lower flux densities while terminating particles in higher flux
regions. In contrast, source biasing modifies source site sampling behavior to
preferentially track particles more likely to reach phase space regions of
interest.
------------
MAGIC Method
@ -132,3 +134,71 @@ aware of this.
:label: variance_fom
\text{FOM} = \frac{1}{\text{Time} \times \sigma^2}
.. _methods_source_biasing:
--------------
Source Biasing
--------------
In contrast to the previous two methods that introduce population controls
during transport, source biasing modifies the sampling of the external source
distribution. The basic premise of the technique is that for each spatial,
angular, energy, or time distribution of a source, an additional distribution
can be specified provided that the two share a common support (set of points
where the distribution is nonzero). Samples are then drawn from this "bias"
distribution, which can be chosen to preferentially direct particles towards
phase space regions of interest. In order to avoid biasing the tally results,
however, a weight adjustment is applied to each sampled site as described below.
Assume that the unbiased probability density function of a random variable
:math:`X:x \rightarrow \mathbb{R}` is given by :math:`f(x)`, but that using the
biased distribution :math:`g(x)` will result in a greater number of particle
trajectories reaching some phase space region of interest. Then a sample
:math:`x_0` may be drawn from :math:`g(x)` while maintaining a fair game,
provided that its weight is adjusted as:
.. math::
:label: source_bias
w = w_0 \times \frac{f(x_0)}{g(x_0)}
where :math:`w_0` is the weight of an unbiased sample from :math:`f(x)`,
typically unity.
Returning now to Equation :eq:`source_bias`, the requirement for common support
becomes evident. If :math:`\mathrm{supp} (g)` fully contains but is not
identical to :math:`\mathrm{supp} (f)`, then some samples from :math:`g(x)` will
correspond to points where :math:`f(x) = 0`. Thus these source sites would be
assigned a starting weight of 0, meaning the particles would be killed
immediately upon transport, effectively wasting computation time. Conversely, if
:math:`\mathrm{supp} (g)` is fully contained by but not identical to
:math:`\mathrm{supp} (f)`, the contributions of some regions outside
:math:`\mathrm{supp} (g)` will not be counted towards the integral, potentially
biasing the tally. The weight assigned to such points would be undefined since
:math:`g(x) = \mathbf{0}` at these points.
When an independent source is sampled in OpenMC, the particle's coordinate in
each variable of phase space :math:`(\mathbf{r},\mathbf{\Omega},E,t)` is
successively drawn from an independent probability distribution. Multiple
variables can be biased, in which case the resultant weight :math:`w` applied to
the particle is the product of the weights assigned from all sampled
distributions: space, angle, energy, and time, as shown in Equation
:eq:`tot_wgt`.
.. math::
:label: tot_wgt
w = w_r \times w_{\Omega} \times w_E \times w_t
Finally, source biasing and weight windows serve different purposes. Source
biasing changes how particles are born, allowing the initial source sites to be
sampled preferentially from important regions of phase space (space, angle,
energy, and time) with an accompanying weight adjustment. Weight windows, by
contrast, apply population control during transport (splitting and Russian
roulette) to help particles reach and contribute in important regions as they
move through the system. Because particle transport proceeds as usual after a
biased source is sampled, particle attenuation in optically thick regions
outside the source volume will not be affected by source biasing; in such
scenarios, transport biasing techniques such as weight windows are often more
effective.

6
docs/source/usersguide/settings.rst
View file

@ -272,6 +272,12 @@ option::
settings.source = [src1, src2]
settings.uniform_source_sampling = True
Additionally, sampling from an :class:`openmc.IndependentSource` may be biased
for local or global variance reduction by modifying the
:attr:`~openmc.IndependentSource.bias` attribute of each of its four main
distributions. Further discussion of source biasing can be found in
:ref:`source_biasing`.
Finally, the :attr:`IndependentSource.particle` attribute can be used to
indicate the source should be composed of particles other than neutrons. For
example, the following would generate a photon source::

155
docs/source/usersguide/variance_reduction.rst
View file

@ -5,10 +5,12 @@ Variance Reduction
==================
Global variance reduction in OpenMC is accomplished by weight windowing
techniques. OpenMC is capable of generating weight windows using either the
MAGIC or FW-CADIS methods. Both techniques will produce a ``weight_windows.h5``
file that can be loaded and used later on. In this section, we break down the
steps required to both generate and then apply weight windows.
or source biasing techniques, the latter of which additionally provides a
local variance reduction capability. OpenMC is capable of generating weight
windows using either the MAGIC or FW-CADIS methods. Both techniques will
produce a ``weight_windows.h5`` file that can be loaded and used later on. In
this section, we first break down the steps required to generate and apply
weight windows, then describe how source biasing may be applied.
.. _ww_generator:
@ -172,3 +174,148 @@ Weight window mesh information is embedded into the weight window file, so the
mesh does not need to be redefined. Monte Carlo solves that load a weight window
file as above will utilize weight windows to reduce the variance of the
simulation.
.. _source_biasing:
--------------
Source Biasing
--------------
In fixed source problems, source biasing provides a means to reduce the variance
on global or localized responses, depending on the biasing scheme. In either
case, the premise of the method is to sample source sites from a biased
distribution that directs a larger fraction of the simulated histories towards
phase space regions of interest than would be found there under analog sampling.
In order to preserve an unbiased estimate of the tally mean, the weight of these
with analog sampling, divided by the probability assigned by the biased
distribution. While the assignment of statistical weights is outlined in the
:ref:`methods section <methods_source_biasing>`, this section demonstrates the
implementation of source biasing to problems in OpenMC.
Source biasing in OpenMC is accomplished by applying a distribution to the
:attr:`bias` attribute of one or more of the univariate or independent
multivariate distributions which make up an :class:`~openmc.IndependentSource`
instance as follows::
# First create the biased distribution
biased_dist = openmc.stats.PowerLaw(a=0, b=3, n=3)
# Construct a new distribution with the bias applied
dist = openmc.stats.PowerLaw(a=0, b=3, n=2, bias=biased_dist)
# The bias attribute can also be set on an existing "analog" distribution:
sphere_dist = openmc.stats.spherical_uniform(r_outer=3)
sphere_dist.r.bias = biased_dist
Univariate distributions may be sampled via the Python API, returning the
sample(s) along with the associated weight(s)::
sample_vec, wgt_vec = dist.sample(n_samples=100)
Here, if the distribution is unbiased, the weight of each sample will be unity.
Finally, :class:`~openmc.IndependentSource` instances can be constructed with
biased distributions::
# Create a source with a biased spatial distribution
source = openmc.IndependentSource(space=sphere_dist)
During the simulation, source sites are then sampled using the biased
distributions where available and given starting statistical weights
corresponding to the cumulative product of the weights assigned by each
distribution in the source object. Hence multiple source variables (e.g.,
direction and energy) may be biased and the resulting source sites will have
their weights adjusted accordingly.
.. note::
Combining source biasing with weight windows can be a powerful variance
reduction technique if each is constructed appropriately for the response
of interest. For example, if a source biasing scheme is devised for
variance reduction of a specific localized response, the user may be able
to specify their own weight window structure that results in more efficient
transport than if weight windows were generated by either of OpenMC's
automatic weight window generators, which are intended for global variance
reduction.
Biased distributions that could result in degenerate weight mappings are not
recommended; this is most commonly seen when biasing the :math:`\phi`-coordinate
of spherical or cylindrical independent multivariate distributions. In such
cases degenerate behavior will be observed at the pole about which :math:`\phi`
is measured, with all values of :math:`\phi` (hence many possible statistical
weights) mapping to the same point for :math:`r=0` or :math:`\mu=0`, and large
weight gradients in the vicinity. In most cases requiring a spherical
independent source, it would be preferable to reorient the reference vector of
the distribution such that biasing could be applied to the
:math:`\mu`-coordinate instead.
When biasing a distribution, care should also be taken to ensure that both the
unbiased and biased distribution share a common support---that is, every region
of phase space mapped to a nonzero probability density by the unbiased
distribution should likewise map to nonzero probability under the biased
distribution, and vice versa. In OpenMC, this places restrictions on the set of
compatible distributions that may be used to bias sampling of each distribution
type. The following table summarizes the method for each distribution in OpenMC
that permits biased sampling.
.. list-table:: **Distributions that support biased sampling**
:header-rows: 1
:widths: 35 65
* - Discrete Univariate PDFs
- Biasing Method
* - :class:`openmc.stats.Discrete`
- Apply a vector of alternative probabilities to the :attr:`bias`
attribute
.. list-table::
:header-rows: 1
:widths: 35 65
* - Continuous Univariate PDFs
- Biasing Method
* - :class:`openmc.stats.Uniform`,
:class:`openmc.stats.PowerLaw`,
:class:`openmc.stats.Maxwell`,
:class:`openmc.stats.Watt`,
:class:`openmc.stats.Normal`,
:class:`openmc.stats.Tabular`
- Apply a second, unbiased continous univariate PDF to the :attr:`bias`
attribute, ensuring that the :attr:`support` attribute of each
distribution is the same
.. list-table::
:header-rows: 1
:widths: 35 65
* - Mixed Univariate PDFs
- Biasing Method
* - :class:`openmc.stats.Mixture`
- May be constructed from multiple biased univariate distributions, or a
second, unbiased continous univariate PDF may be applied to the
:attr:`bias` attribute
.. list-table::
:header-rows: 1
:widths: 35 65
* - Discrete Multivariate PDFs
- Biasing Method
* - :class:`openmc.stats.PointCloud`,
:class:`openmc.stats.MeshSpatial`
- Apply a vector of the new relative probabilities of each point or mesh
element under biased sampling to the :attr:`bias` attribute
.. list-table::
:header-rows: 1
:widths: 35 65
* - Continuous Multivariate PDFs
- Biasing Method
* - :class:`openmc.stats.CartesianIndependent`,
:class:`openmc.stats.CylindricalIndependent`,
:class:`openmc.stats.SphericalIndependent`,
:class:`openmc.stats.PolarAzimuthal`
- Construct from biased univariate distributions for :attr:`x`, :attr:`y`,
:attr:`z`, etc.
* - :class:`openmc.stats.Isotropic`
- Apply an unbiased :class:`openmc.stats.PolarAzimuthal` to the
:attr:`bias` attribute

188
include/openmc/distribution.h
View file

@ -15,6 +15,17 @@
namespace openmc {
//==============================================================================
// Helper function for computing importance weights from biased sampling
//==============================================================================
//! Compute importance weights for biased sampling
//! \param p Unnormalized original probability vector
//! \param b Unnormalized bias probability vector
//! \return Vector of importance weights (p_norm[i] / b_norm[i])
vector<double> compute_importance_weights(
const vector<double>& p, const vector<double>& b);
//==============================================================================
//! Abstract class representing a univariate probability distribution
//==============================================================================
@ -22,11 +33,41 @@ namespace openmc {
class Distribution {
public:
virtual ~Distribution() = default;
virtual double sample(uint64_t* seed) const = 0;
//! Sample a value from the distribution, handling biasing automatically
//! \param seed Pseudorandom number seed pointer
//! \return (sampled value, importance weight)
virtual std::pair<double, double> sample(uint64_t* seed) const;
//! Evaluate probability density, f(x), at a point
//! \param x Point to evaluate f(x)
//! \return f(x)
virtual double evaluate(double x) const;
//! Return integral of distribution
//! \return Integral of distribution
virtual double integral() const { return 1.0; };
//! Set bias distribution
virtual void set_bias(std::unique_ptr<Distribution> bias)
{
bias_ = std::move(bias);
}
const Distribution* bias() const { return bias_.get(); }
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
virtual double sample_unbiased(uint64_t* seed) const = 0;
//! Read bias distribution from XML
//! \param node XML node that may contain a bias child element
void read_bias_from_xml(pugi::xml_node node);
// Biasing distribution
unique_ptr<Distribution> bias_;
};
using UPtrDist = unique_ptr<Distribution>;
@ -50,7 +91,7 @@ public:
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
//! \return Sampled index
size_t sample(uint64_t* seed) const;
// Properties
@ -67,7 +108,7 @@ private:
//! Normalize distribution so that probabilities sum to unity
void normalize();
//! Initialize alias tables for distribution
//! Initialize alias table for sampling
void init_alias();
};
@ -82,20 +123,30 @@ public:
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t* seed) const override;
//! \return (sampled value, sample weight)
std::pair<double, double> sample(uint64_t* seed) const override;
double integral() const override { return di_.integral(); };
//! Override set_bias as no-op (bias handled in constructor)
void set_bias(std::unique_ptr<Distribution> bias) override {}
// Properties
const vector<double>& x() const { return x_; }
const vector<double>& prob() const { return di_.prob(); }
const vector<size_t>& alias() const { return di_.alias(); }
const vector<double>& weight() const { return weight_; }
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
vector<double> x_; //!< Possible outcomes
DiscreteIndex di_; //!< discrete probability distribution of
//!< outcome indices
vector<double> x_; //!< Possible outcomes
vector<double> weight_; //!< Importance weights (empty if unbiased)
DiscreteIndex di_; //!< Discrete probability distribution of outcome indices
};
//==============================================================================
@ -107,14 +158,20 @@ public:
explicit Uniform(pugi::xml_node node);
Uniform(double a, double b) : a_ {a}, b_ {b} {};
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t* seed) const override;
//! Evaluate probability density, f(x), at a point
//! \param x Point to evaluate f(x)
//! \return f(x)
double evaluate(double x) const override;
double a() const { return a_; }
double b() const { return b_; }
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
double a_; //!< Lower bound of distribution
double b_; //!< Upper bound of distribution
@ -131,15 +188,21 @@ public:
: offset_ {std::pow(a, n + 1)}, span_ {std::pow(b, n + 1) - offset_},
ninv_ {1 / (n + 1)} {};
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t* seed) const override;
//! Evaluate probability density, f(x), at a point
//! \param x Point to evaluate f(x)
//! \return f(x)
double evaluate(double x) const override;
double a() const { return std::pow(offset_, ninv_); }
double b() const { return std::pow(offset_ + span_, ninv_); }
double n() const { return 1 / ninv_ - 1; }
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
//! Store processed values in object to allow for faster sampling
double offset_; //!< a^(n+1)
@ -156,13 +219,19 @@ public:
explicit Maxwell(pugi::xml_node node);
Maxwell(double theta) : theta_ {theta} {};
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t* seed) const override;
//! Evaluate probability density, f(x), at a point
//! \param x Point to evaluate f(x)
//! \return f(x)
double evaluate(double x) const override;
double theta() const { return theta_; }
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
double theta_; //!< Factor in exponential [eV]
};
@ -176,14 +245,20 @@ public:
explicit Watt(pugi::xml_node node);
Watt(double a, double b) : a_ {a}, b_ {b} {};
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t* seed) const override;
//! Evaluate probability density, f(x), at a point
//! \param x Point to evaluate f(x)
//! \return f(x)
double evaluate(double x) const override;
double a() const { return a_; }
double b() const { return b_; }
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
double a_; //!< Factor in exponential [eV]
double b_; //!< Factor in square root [1/eV]
@ -200,14 +275,20 @@ public:
Normal(double mean_value, double std_dev)
: mean_value_ {mean_value}, std_dev_ {std_dev} {};
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t* seed) const override;
//! Evaluate probability density, f(x), at a point
//! \param x Point to evaluate f(x)
//! \return f(x)
double evaluate(double x) const override;
double mean_value() const { return mean_value_; }
double std_dev() const { return std_dev_; }
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
double mean_value_; //!< middle of distribution [eV]
double std_dev_; //!< standard deviation [eV]
@ -223,10 +304,10 @@ public:
Tabular(const double* x, const double* p, int n, Interpolation interp,
const double* c = nullptr);
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t* seed) const override;
//! Evaluate probability density, f(x), at a point
//! \param x Point to evaluate f(x)
//! \return f(x)
double evaluate(double x) const override;
// properties
vector<double>& x() { return x_; }
@ -235,6 +316,12 @@ public:
Interpolation interp() const { return interp_; }
double integral() const override { return integral_; };
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
vector<double> x_; //!< tabulated independent variable
vector<double> p_; //!< tabulated probability density
@ -259,13 +346,19 @@ public:
explicit Equiprobable(pugi::xml_node node);
Equiprobable(const double* x, int n) : x_ {x, x + n} {};
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t* seed) const override;
//! Evaluate probability density, f(x), at a point
//! \param x Point to evaluate f(x)
//! \return f(x)
double evaluate(double x) const override;
const vector<double>& x() const { return x_; }
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
vector<double> x_; //! Possible outcomes
};
@ -280,18 +373,25 @@ public:
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t* seed) const override;
//! \return (sampled value, sample weight)
std::pair<double, double> sample(uint64_t* seed) const override;
double integral() const override { return integral_; }
private:
// Storrage for probability + distribution
using DistPair = std::pair<double, UPtrDist>;
//! Override set_bias as no-op (bias handled in constructor)
void set_bias(std::unique_ptr<Distribution> bias) override {}
vector<DistPair>
distribution_; //!< sub-distributions + cummulative probabilities
double integral_; //!< integral of distribution
protected:
//! Sample a value (unbiased) from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample_unbiased(uint64_t* seed) const override;
private:
vector<UPtrDist> distribution_; //!< Sub-distributions
vector<double> weight_; //!< Importance weights for component selection
DiscreteIndex di_; //!< Discrete probability distribution of indices
double integral_; //!< Integral of distribution
};
} // namespace openmc

43
include/openmc/distribution_multi.h
View file

@ -26,8 +26,8 @@ public:
//! Sample a direction from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Direction sampled
virtual Direction sample(uint64_t* seed) const = 0;
//! \return (sampled Direction, sample weight)
virtual std::pair<Direction, double> sample(uint64_t* seed) const = 0;
Direction u_ref_ {0.0, 0.0, 1.0}; //!< reference direction
};
@ -43,14 +43,28 @@ public:
//! Sample a direction from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Direction sampled
Direction sample(uint64_t* seed) const override;
//! \return (sampled Direction, sample weight)
std::pair<Direction, double> sample(uint64_t* seed) const override;
//! Sample a direction and return evaluation of the PDF for biased sampling.
//! Note that bias distributions are intended to return unit-weight samples.
//! \param seed Pseudorandom number seed points
//! \return (sampled Direction, value of the PDF at this Direction)
std::pair<Direction, double> sample_as_bias(uint64_t* seed) const;
// Observing pointers
Distribution* mu() const { return mu_.get(); }
Distribution* phi() const { return phi_.get(); }
private:
//! Common sampling implementation
//! \param seed Pseudorandom number seed pointer
//! \param return_pdf If true, return PDF evaluation; if false, return
//! importance weight
//! \return (sampled Direction, weight or PDF value)
std::pair<Direction, double> sample_impl(
uint64_t* seed, bool return_pdf) const;
Direction v_ref_ {1.0, 0.0, 0.0}; //!< reference direction
Direction w_ref_;
UPtrDist mu_; //!< Distribution of polar angle
@ -66,11 +80,24 @@ Direction isotropic_direction(uint64_t* seed);
class Isotropic : public UnitSphereDistribution {
public:
Isotropic() {};
explicit Isotropic(pugi::xml_node node);
//! Sample a direction from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled direction
Direction sample(uint64_t* seed) const override;
//! \return (sampled direction, sample weight)
std::pair<Direction, double> sample(uint64_t* seed) const override;
// Set or get bias distribution
void set_bias(std::unique_ptr<PolarAzimuthal> bias)
{
bias_ = std::move(bias);
}
const PolarAzimuthal* bias() const { return bias_.get(); }
protected:
// Biasing distribution
unique_ptr<PolarAzimuthal> bias_;
};
//==============================================================================
@ -85,8 +112,8 @@ public:
//! Sample a direction from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled direction
Direction sample(uint64_t* seed) const override;
//! \return (sampled direction, sample weight)
std::pair<Direction, double> sample(uint64_t* seed) const override;
};
using UPtrAngle = unique_ptr<UnitSphereDistribution>;

37
include/openmc/distribution_spatial.h
View file

@ -19,7 +19,9 @@ public:
virtual ~SpatialDistribution() = default;
//! Sample a position from the distribution
virtual Position sample(uint64_t* seed) const = 0;
//! \param seed Pseudorandom number seed pointer
//! \return Sampled (position, importance weight)
virtual std::pair<Position, double> sample(uint64_t* seed) const = 0;
static unique_ptr<SpatialDistribution> create(pugi::xml_node node);
};
@ -34,8 +36,8 @@ public:
//! Sample a position from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const override;
//! \return Sampled (position, importance weight)
std::pair<Position, double> sample(uint64_t* seed) const override;
// Observer pointers
Distribution* x() const { return x_.get(); }
@ -58,8 +60,8 @@ public:
//! Sample a position from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const override;
//! \return Sampled (position, importance weight)
std::pair<Position, double> sample(uint64_t* seed) const override;
Distribution* r() const { return r_.get(); }
Distribution* phi() const { return phi_.get(); }
@ -83,8 +85,8 @@ public:
//! Sample a position from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const override;
//! \return Sampled (position, importance weight)
std::pair<Position, double> sample(uint64_t* seed) const override;
Distribution* r() const { return r_.get(); }
Distribution* cos_theta() const { return cos_theta_.get(); }
@ -109,8 +111,8 @@ public:
//! Sample a position from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const override;
//! \return Sampled (position, importance weight)
std::pair<Position, double> sample(uint64_t* seed) const override;
//! Sample the mesh for an element and position within that element
//! \param seed Pseudorandom number seed pointer
@ -133,8 +135,8 @@ public:
private:
int32_t mesh_idx_ {C_NONE};
DiscreteIndex elem_idx_dist_; //!< Distribution of
//!< mesh element indices
DiscreteIndex elem_idx_dist_; //!< Distribution of mesh element indices
vector<double> weight_; //!< Importance weights (empty if unbiased)
};
//==============================================================================
@ -149,12 +151,13 @@ public:
//! Sample a position from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const override;
//! \return Sampled (position, importance weight)
std::pair<Position, double> sample(uint64_t* seed) const override;
private:
std::vector<Position> point_cloud_;
DiscreteIndex point_idx_dist_; //!< Distribution of Position indices
vector<double> weight_; //!< Importance weights (empty if unbiased)
};
//==============================================================================
@ -167,8 +170,8 @@ public:
//! Sample a position from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const override;
//! \return Sampled (position, importance weight)
std::pair<Position, double> sample(uint64_t* seed) const override;
// Properties
bool only_fissionable() const { return only_fissionable_; }
@ -193,8 +196,8 @@ public:
//! Sample a position from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const override;
//! \return Sampled (position, importance weight)
std::pair<Position, double> sample(uint64_t* seed) const override;
Position r() const { return r_; }

4
include/openmc/source.h
View file

@ -217,8 +217,8 @@ public:
//! Sample a position from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled position
Position sample(uint64_t* seed) const override;
//! \return (sampled position, importance weight)
std::pair<Position, double> sample(uint64_t* seed) const override;
private:
int32_t mesh_index_ {C_NONE}; //!< Index in global meshes array

2
openmc/filter.py
View file

@ -815,7 +815,7 @@ class ParentNuclideFilter(ParticleFilter):
class MeshFilter(Filter):
"""Bins tally event locations by mesh elements.
r"""Bins tally event locations by mesh elements.
Parameters
----------

9
openmc/settings.py
View file

@ -10,7 +10,7 @@ import warnings
import openmc
import openmc.checkvalue as cv
from openmc.checkvalue import PathLike
from openmc.stats.multivariate import MeshSpatial
from openmc.stats.multivariate import MeshSpatial, Box, PolarAzimuthal, Isotropic
from ._xml import clean_indentation, get_elem_list, get_text
from .mesh import _read_meshes, RegularMesh, MeshBase
from .source import SourceBase, MeshSource, IndependentSource
@ -1877,7 +1877,11 @@ class Settings:
for key, value in self._random_ray.items():
if key == 'ray_source' and isinstance(value, SourceBase):
source_element = value.to_xml_element()
if source_element.find('bias') is not None:
raise RuntimeError(
"Ray source distributions should not be biased.")
element.append(source_element)
elif key == 'source_region_meshes':
subelement = ET.SubElement(element, 'source_region_meshes')
for mesh, domains in value:
@ -2324,6 +2328,9 @@ class Settings:
self.random_ray[child.tag] = float(child.text)
elif child.tag == 'source':
source = SourceBase.from_xml_element(child)
if child.find('bias') is not None:
raise RuntimeError(
"Ray source distributions should not be biased.")
self.random_ray['ray_source'] = source
elif child.tag == 'volume_estimator':
self.random_ray['volume_estimator'] = child.text

142
openmc/stats/multivariate.py
View file

@ -81,7 +81,7 @@ class PolarAzimuthal(UnitSphere):
z-direction.
reference_vwu : Iterable of float
Direction from which azimuthal angle is measured. Defaults to the positive
x-direction.
x-direction.
Attributes
----------
@ -104,7 +104,7 @@ class PolarAzimuthal(UnitSphere):
self.phi = phi
else:
self.phi = Uniform(0., 2*pi)
@property
def reference_vwu(self):
return self._reference_vwu
@ -114,8 +114,8 @@ class PolarAzimuthal(UnitSphere):
cv.check_type('reference v direction', vwu, Iterable, Real)
vwu = np.asarray(vwu)
uvw = self.reference_uvw
cv.check_greater_than('reference v direction must not be parallel to reference u direction', np.linalg.norm(np.cross(vwu,uvw)), 1e-6*np.linalg.norm(vwu))
vwu -= vwu.dot(uvw)*uvw
cv.check_greater_than('reference v direction must not be parallel to reference u direction', np.linalg.norm(np.cross(vwu,uvw)), 1e-6*np.linalg.norm(vwu))
vwu -= vwu.dot(uvw)*uvw
cv.check_less_than('reference v direction must be orthogonal to reference u direction', np.abs(vwu.dot(uvw)), 1e-6)
self._reference_vwu = vwu/np.linalg.norm(vwu)
@ -137,21 +137,30 @@ class PolarAzimuthal(UnitSphere):
cv.check_type('azimuthal angle', phi, Univariate)
self._phi = phi
def to_xml_element(self):
def to_xml_element(self, element_name: str = None):
"""Return XML representation of the angular distribution
Parameters
----------
element_name : str, optional
XML element name
Returns
-------
element : lxml.etree._Element
XML element containing angular distribution data
"""
element = ET.Element('angle')
if element_name is not None:
element = ET.Element(element_name)
else:
element = ET.Element('angle')
element.set("type", "mu-phi")
if self.reference_uvw is not None:
element.set("reference_uvw", ' '.join(map(str, self.reference_uvw)))
if self.reference_vwu is not None:
element.set("reference_vwu", ' '.join(map(str, self.reference_vwu)))
element.set("reference_vwu", ' '.join(map(str, self.reference_vwu)))
element.append(self.mu.to_xml_element('mu'))
element.append(self.phi.to_xml_element('phi'))
return element
@ -177,17 +186,48 @@ class PolarAzimuthal(UnitSphere):
mu_phi.reference_uvw = uvw
vwu = get_elem_list(elem, "reference_vwu", float)
if vwu is not None:
mu_phi.reference_vwu = vwu
mu_phi.reference_vwu = vwu
mu_phi.mu = Univariate.from_xml_element(elem.find('mu'))
mu_phi.phi = Univariate.from_xml_element(elem.find('phi'))
return mu_phi
class Isotropic(UnitSphere):
"""Isotropic angular distribution."""
"""Isotropic angular distribution.
def __init__(self):
Parameters
----------
bias : openmc.stats.PolarAzimuthal, optional
Distribution for biased sampling.
Attributes
----------
bias : openmc.stats.PolarAzimuthal or None
Distribution for biased sampling
"""
def __init__(self, bias: PolarAzimuthal | None = None):
super().__init__()
self.bias = bias
@property
def bias(self):
return self._bias
@bias.setter
def bias(self, bias):
cv.check_type('Biasing distribution', bias, PolarAzimuthal, none_ok=True)
if bias is not None:
if (bias.mu.bias is not None) or (bias.phi.bias is not None):
raise RuntimeError('Biasing distributions should not have their own bias.')
elif (bias.mu.support != (-1., 1.)
or not np.all(np.isclose(bias.phi.support, (0., 2*np.pi)))):
raise ValueError("Biasing distribution for an isotropic "
"distribution should be supported on "
"mu=(-1.0,1.0) and phi=(0.0,2*pi).")
self._bias = bias
def to_xml_element(self):
"""Return XML representation of the isotropic distribution
@ -200,6 +240,15 @@ class Isotropic(UnitSphere):
"""
element = ET.Element('angle')
element.set("type", "isotropic")
if self.bias is not None:
bias_dist = self.bias
if (bias_dist.mu.bias is not None) or (bias_dist.phi.bias is not None):
raise RuntimeError('Biasing distributions should not have their own bias!')
else:
bias_elem = self.bias.to_xml_element("bias")
element.append(bias_elem)
return element
@classmethod
@ -217,7 +266,13 @@ class Isotropic(UnitSphere):
Isotropic distribution generated from XML element
"""
return cls()
bias_elem = elem.find('bias')
if bias_elem is not None:
bias_dist = PolarAzimuthal.from_xml_element(bias_elem)
return cls(bias=bias_dist)
else:
return cls()
class Monodirectional(UnitSphere):
@ -672,6 +727,9 @@ class MeshSpatial(Spatial):
volume_normalized : bool, optional
Whether or not the strengths will be multiplied by element volumes at
runtime. Default is True.
bias : iterable of float, optional
An iterable of values giving the selection weights assigned to each
element during biased sampling.
Attributes
----------
@ -682,12 +740,16 @@ class MeshSpatial(Spatial):
volume_normalized : bool
Whether or not the strengths will be multiplied by element volumes at
runtime.
bias : numpy.ndarray or None
Distribution for biased sampling
"""
def __init__(self, mesh, strengths=None, volume_normalized=True):
def __init__(self, mesh, strengths=None, volume_normalized=True,
bias: Sequence[float] | None = None):
self.mesh = mesh
self.strengths = strengths
self.volume_normalized = volume_normalized
self.bias = bias
@property
def mesh(self):
@ -720,6 +782,23 @@ class MeshSpatial(Spatial):
else:
self._strengths = None
@property
def bias(self):
return self._bias
@bias.setter
def bias(self, given_bias):
if given_bias is not None:
cv.check_type('Biasing strengths array', given_bias, Iterable, Real)
bias_array = np.asarray(given_bias, dtype=float).flatten()
if bias_array.size != self.strengths.size:
raise ValueError(
'Bias strengths array must have same size as strengths array.')
else:
self._bias = bias_array
else:
self._bias = None
@property
def num_strength_bins(self):
if self.strengths is None:
@ -745,6 +824,9 @@ class MeshSpatial(Spatial):
subelement = ET.SubElement(element, 'strengths')
subelement.text = ' '.join(str(e) for e in self.strengths)
if self.bias is not None:
Univariate._append_array_bias_to_xml(self, element)
return element
@classmethod
@ -774,7 +856,8 @@ class MeshSpatial(Spatial):
volume_normalized = get_text(elem, 'volume_normalized').lower() == 'true'
strengths = get_elem_list(elem, 'strengths', float)
return cls(meshes[mesh_id], strengths, volume_normalized)
bias_strengths = Univariate._read_array_bias_from_xml(elem)
return cls(meshes[mesh_id], strengths, volume_normalized, bias=bias_strengths)
class PointCloud(Spatial):
@ -792,6 +875,9 @@ class PointCloud(Spatial):
strengths : iterable of float, optional
An iterable of values that represents the relative probabilty of each
point.
bias : iterable of float, optional
An iterable of values representing the relative probability of each
point under biased sampling.
Attributes
----------
@ -799,15 +885,19 @@ class PointCloud(Spatial):
The points in space to be sampled with shape (N, 3)
strengths : numpy.ndarray or None
An array of relative probabilities for each mesh point
bias : numpy.ndarray or None
An array of relative probabilities for biased sampling of mesh points
"""
def __init__(
self,
positions: Sequence[Sequence[float]],
strengths: Sequence[float] | None = None
strengths: Sequence[float] | None = None,
bias: Sequence[float] | None = None
):
self.positions = positions
self.strengths = strengths
self.bias = bias
@property
def positions(self) -> np.ndarray:
@ -836,6 +926,23 @@ class PointCloud(Spatial):
raise ValueError('strengths must have the same length as positions')
self._strengths = strengths
@property
def bias(self):
return self._bias
@bias.setter
def bias(self, given_bias):
if given_bias is not None:
cv.check_type('Biasing strengths array', given_bias, Iterable, Real)
bias_array = np.asarray(given_bias, dtype=float).flatten()
if bias_array.size != self.strengths.size:
raise ValueError(
'Bias strengths array must have same size as strengths array.')
else:
self._bias = bias_array
else:
self._bias = None
@property
def num_strength_bins(self) -> int:
if self.strengths is None:
@ -861,6 +968,9 @@ class PointCloud(Spatial):
subelement = ET.SubElement(element, 'strengths')
subelement.text = ' '.join(str(e) for e in self.strengths)
if self.bias is not None:
Univariate._append_array_bias_to_xml(self, element)
return element
@classmethod
@ -883,8 +993,8 @@ class PointCloud(Spatial):
positions = np.array(coord_data).reshape((-1, 3))
strengths = get_elem_list(elem, 'strengths', float)
return cls(positions, strengths)
bias_strengths = Univariate._read_array_bias_from_xml(elem)
return cls(positions, strengths, bias=bias_strengths)
class Box(Spatial):

574
openmc/stats/univariate.py
View file

File diff suppressed because it is too large Load diff

4
src/chain.cpp
View file

@ -70,8 +70,8 @@ ChainNuclide::~ChainNuclide()
void DecayPhotonAngleEnergy::sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const
{
E_out = photon_energy_->sample(seed);
mu = Uniform(-1., 1.).sample(seed);
E_out = photon_energy_->sample(seed).first;
mu = Uniform(-1., 1.).sample(seed).first;
}
//==============================================================================

288
src/distribution.cpp
View file

@ -8,6 +8,7 @@
#include <stdexcept> // for runtime_error
#include <string> // for string, stod
#include "openmc/constants.h"
#include "openmc/error.h"
#include "openmc/math_functions.h"
#include "openmc/random_dist.h"
@ -16,6 +17,74 @@
namespace openmc {
//==============================================================================
// Helper function for computing importance weights from biased sampling
//==============================================================================
vector<double> compute_importance_weights(
const vector<double>& p, const vector<double>& b)
{
std::size_t n = p.size();
// Normalize original probabilities
double sum_p = std::accumulate(p.begin(), p.end(), 0.0);
vector<double> p_norm(n);
for (std::size_t i = 0; i < n; ++i) {
p_norm[i] = p[i] / sum_p;
}
// Normalize bias probabilities
double sum_b = std::accumulate(b.begin(), b.end(), 0.0);
vector<double> b_norm(n);
for (std::size_t i = 0; i < n; ++i) {
b_norm[i] = b[i] / sum_b;
}
// Compute importance weights
vector<double> weights(n);
for (std::size_t i = 0; i < n; ++i) {
weights[i] = (b_norm[i] == 0.0) ? INFTY : p_norm[i] / b_norm[i];
}
return weights;
}
std::pair<double, double> Distribution::sample(uint64_t* seed) const
{
if (bias_) {
// Sample from the bias distribution and compute importance weight
double val = bias_->sample_unbiased(seed);
double wgt = this->evaluate(val) / bias_->evaluate(val);
return {val, wgt};
} else {
// Unbiased sampling: return sampled value with weight 1.0
double val = sample_unbiased(seed);
return {val, 1.0};
}
}
// PDF evaluation not supported for all distribution types
double Distribution::evaluate(double x) const
{
throw std::runtime_error(
"PDF evaluation not implemented for this distribution type.");
}
void Distribution::read_bias_from_xml(pugi::xml_node node)
{
if (check_for_node(node, "bias")) {
pugi::xml_node bias_node = node.child("bias");
if (check_for_node(bias_node, "bias")) {
openmc::fatal_error(
"Distribution has a bias distribution with its own bias distribution. "
"Please ensure bias distributions do not have their own bias.");
}
UPtrDist bias = distribution_from_xml(bias_node);
this->set_bias(std::move(bias));
}
}
//==============================================================================
// DiscreteIndex implementation
//==============================================================================
@ -36,7 +105,6 @@ DiscreteIndex::DiscreteIndex(span<const double> p)
void DiscreteIndex::assign(span<const double> p)
{
prob_.assign(p.begin(), p.end());
this->init_alias();
}
@ -115,24 +183,55 @@ void DiscreteIndex::normalize()
// Discrete implementation
//==============================================================================
Discrete::Discrete(pugi::xml_node node) : di_(node)
Discrete::Discrete(pugi::xml_node node)
{
auto params = get_node_array<double>(node, "parameters");
std::size_t n = params.size() / 2;
// First half is x values, second half is probabilities
x_.assign(params.begin(), params.begin() + n);
const double* p = params.data() + n;
// Check for bias
if (check_for_node(node, "bias")) {
// Get bias probabilities
auto bias_params = get_node_array<double>(node, "bias");
if (bias_params.size() != n) {
openmc::fatal_error(
"Size mismatch: Attempted to bias Discrete distribution with " +
std::to_string(n) + " probability entries using a bias with " +
std::to_string(bias_params.size()) +
" entries. Please ensure distributions have the same size.");
}
// Compute importance weights
vector<double> p_vec(p, p + n);
weight_ = compute_importance_weights(p_vec, bias_params);
// Initialize DiscreteIndex with bias probabilities for sampling
di_.assign(bias_params);
} else {
// Unbiased case: weight_ stays empty
di_.assign({p, n});
}
}
Discrete::Discrete(const double* x, const double* p, size_t n) : di_({p, n})
{
x_.assign(x, x + n);
}
double Discrete::sample(uint64_t* seed) const
std::pair<double, double> Discrete::sample(uint64_t* seed) const
{
return x_[di_.sample(seed)];
size_t idx = di_.sample(seed);
double wgt = weight_.empty() ? 1.0 : weight_[idx];
return {x_[idx], wgt};
}
double Discrete::sample_unbiased(uint64_t* seed) const
{
size_t idx = di_.sample(seed);
return x_[idx];
}
//==============================================================================
@ -149,13 +248,26 @@ Uniform::Uniform(pugi::xml_node node)
a_ = params.at(0);
b_ = params.at(1);
read_bias_from_xml(node);
}
double Uniform::sample(uint64_t* seed) const
double Uniform::sample_unbiased(uint64_t* seed) const
{
return a_ + prn(seed) * (b_ - a_);
}
double Uniform::evaluate(double x) const
{
if (x <= a()) {
return 0.0;
} else if (x >= b()) {
return 0.0;
} else {
return 1 / (b() - a());
}
}
//==============================================================================
// PowerLaw implementation
//==============================================================================
@ -175,9 +287,24 @@ PowerLaw::PowerLaw(pugi::xml_node node)
offset_ = std::pow(a, n + 1);
span_ = std::pow(b, n + 1) - offset_;
ninv_ = 1 / (n + 1);
read_bias_from_xml(node);
}
double PowerLaw::sample(uint64_t* seed) const
double PowerLaw::evaluate(double x) const
{
if (x <= a()) {
return 0.0;
} else if (x >= b()) {
return 0.0;
} else {
int pwr = n() + 1;
double norm = pwr / span_;
return norm * std::pow(std::fabs(x), n());
}
}
double PowerLaw::sample_unbiased(uint64_t* seed) const
{
return std::pow(offset_ + prn(seed) * span_, ninv_);
}
@ -189,13 +316,21 @@ double PowerLaw::sample(uint64_t* seed) const
Maxwell::Maxwell(pugi::xml_node node)
{
theta_ = std::stod(get_node_value(node, "parameters"));
read_bias_from_xml(node);
}
double Maxwell::sample(uint64_t* seed) const
double Maxwell::sample_unbiased(uint64_t* seed) const
{
return maxwell_spectrum(theta_, seed);
}
double Maxwell::evaluate(double x) const
{
double c = (2.0 / SQRT_PI) * std::pow(theta_, -1.5);
return c * std::sqrt(x) * std::exp(-x / theta_);
}
//==============================================================================
// Watt implementation
//==============================================================================
@ -209,13 +344,22 @@ Watt::Watt(pugi::xml_node node)
a_ = params.at(0);
b_ = params.at(1);
read_bias_from_xml(node);
}
double Watt::sample(uint64_t* seed) const
double Watt::sample_unbiased(uint64_t* seed) const
{
return watt_spectrum(a_, b_, seed);
}
double Watt::evaluate(double x) const
{
double c =
2.0 / (std::sqrt(PI * b_) * std::pow(a_, 1.5) * std::exp(a_ * b_ / 4.0));
return c * std::exp(-x / a_) * std::sinh(std::sqrt(b_ * x));
}
//==============================================================================
// Normal implementation
//==============================================================================
@ -229,13 +373,22 @@ Normal::Normal(pugi::xml_node node)
mean_value_ = params.at(0);
std_dev_ = params.at(1);
read_bias_from_xml(node);
}
double Normal::sample(uint64_t* seed) const
double Normal::sample_unbiased(uint64_t* seed) const
{
return normal_variate(mean_value_, std_dev_, seed);
}
double Normal::evaluate(double x) const
{
return (1.0 / (std::sqrt(2.0 / PI) * std_dev_)) *
std::exp(-(std::pow((x - mean_value_), 2.0)) /
(2.0 * std::pow(std_dev_, 2.0)));
}
//==============================================================================
// Tabular implementation
//==============================================================================
@ -266,6 +419,8 @@ Tabular::Tabular(pugi::xml_node node)
const double* x = params.data();
const double* p = x + n;
init(x, p, n);
read_bias_from_xml(node);
}
Tabular::Tabular(const double* x, const double* p, int n, Interpolation interp,
@ -314,7 +469,7 @@ void Tabular::init(
}
}
double Tabular::sample(uint64_t* seed) const
double Tabular::sample_unbiased(uint64_t* seed) const
{
// Sample value of CDF
double c = prn(seed);
@ -356,11 +511,39 @@ double Tabular::sample(uint64_t* seed) const
}
}
double Tabular::evaluate(double x) const
{
int i;
if (interp_ == Interpolation::histogram) {
i = std::upper_bound(x_.begin(), x_.end(), x) - x_.begin() - 1;
if (i < 0 || i >= static_cast<int>(p_.size())) {
return 0.0;
} else {
return p_[i];
}
} else {
i = std::lower_bound(x_.begin(), x_.end(), x) - x_.begin() - 1;
if (i < 0 || i >= static_cast<int>(p_.size()) - 1) {
return 0.0;
} else {
double x0 = x_[i];
double x1 = x_[i + 1];
double p0 = p_[i];
double p1 = p_[i + 1];
double t = (x - x0) / (x1 - x0);
return (1 - t) * p0 + t * p1;
}
}
}
//==============================================================================
// Equiprobable implementation
//==============================================================================
double Equiprobable::sample(uint64_t* seed) const
double Equiprobable::sample_unbiased(uint64_t* seed) const
{
std::size_t n = x_.size();
@ -372,13 +555,27 @@ double Equiprobable::sample(uint64_t* seed) const
return xl + ((n - 1) * r - i) * (xr - xl);
}
double Equiprobable::evaluate(double x) const
{
double x_min = *std::min_element(x_.begin(), x_.end());
double x_max = *std::max_element(x_.begin(), x_.end());
if (x < x_min || x > x_max) {
return 0.0;
} else {
return 1.0 / (x_max - x_min);
}
}
//==============================================================================
// Mixture implementation
//==============================================================================
Mixture::Mixture(pugi::xml_node node)
{
double cumsum = 0.0;
vector<double> probabilities;
// First pass: collect distributions and their probabilities
for (pugi::xml_node pair : node.children("pair")) {
// Check that required data exists
if (!pair.attribute("probability"))
@ -386,39 +583,60 @@ Mixture::Mixture(pugi::xml_node node)
if (!pair.child("dist"))
fatal_error("Mixture pair element does not have a distribution.");
// cummulative sum of probabilities
// Get probability and distribution
double p = std::stod(pair.attribute("probability").value());
// Save cummulative probability and distribution
auto dist = distribution_from_xml(pair.child("dist"));
cumsum += p * dist->integral();
distribution_.push_back(std::make_pair(cumsum, std::move(dist)));
// Weight probability by the distribution's integral
double weighted_prob = p * dist->integral();
probabilities.push_back(weighted_prob);
distribution_.push_back(std::move(dist));
}
// Save integral of distribution
integral_ = cumsum;
// Save sum of weighted probabilities
integral_ = std::accumulate(probabilities.begin(), probabilities.end(), 0.0);
// Normalize cummulative probabilities to 1
for (auto& pair : distribution_) {
pair.first /= cumsum;
std::size_t n = probabilities.size();
// Check for bias
if (check_for_node(node, "bias")) {
// Get bias probabilities
auto bias_params = get_node_array<double>(node, "bias");
if (bias_params.size() != n) {
openmc::fatal_error(
"Size mismatch: Attempted to bias Mixture distribution with " +
std::to_string(n) + " components using a bias with " +
std::to_string(bias_params.size()) +
" entries. Please ensure distributions have the same size.");
}
// Compute importance weights
weight_ = compute_importance_weights(probabilities, bias_params);
// Initialize DiscreteIndex with bias probabilities for sampling
di_.assign(bias_params);
} else {
// Unbiased case: weight_ stays empty
di_.assign(probabilities);
}
}
double Mixture::sample(uint64_t* seed) const
std::pair<double, double> Mixture::sample(uint64_t* seed) const
{
// Sample value of CDF
const double p = prn(seed);
// find matching distribution
const auto it = std::lower_bound(distribution_.cbegin(), distribution_.cend(),
p, [](const DistPair& pair, double p) { return pair.first < p; });
// This should not happen. Catch it
assert(it != distribution_.cend());
size_t idx = di_.sample(seed);
// Sample the chosen distribution
return it->second->sample(seed);
auto [val, sub_wgt] = distribution_[idx]->sample(seed);
// Multiply by component selection weight
double mix_wgt = weight_.empty() ? 1.0 : weight_[idx];
return {val, mix_wgt * sub_wgt};
}
double Mixture::sample_unbiased(uint64_t* seed) const
{
size_t idx = di_.sample(seed);
return distribution_[idx]->sample(seed).first;
}
//==============================================================================

2
src/distribution_angle.cpp
View file

@ -75,7 +75,7 @@ double AngleDistribution::sample(double E, uint64_t* seed) const
++i;
// Sample i-th distribution
double mu = distribution_[i]->sample(seed);
double mu = distribution_[i]->sample(seed).first;
// Make sure mu is in range [-1,1] and return
if (std::abs(mu) > 1.0)

69
src/distribution_multi.cpp
View file

@ -20,7 +20,7 @@ unique_ptr<UnitSphereDistribution> UnitSphereDistribution::create(
if (check_for_node(node, "type"))
type = get_node_value(node, "type", true, true);
if (type == "isotropic") {
return UPtrAngle {new Isotropic()};
return UPtrAngle {new Isotropic(node)};
} else if (type == "monodirectional") {
return UPtrAngle {new Monodirectional(node)};
} else if (type == "mu-phi") {
@ -82,27 +82,63 @@ PolarAzimuthal::PolarAzimuthal(pugi::xml_node node)
}
}
Direction PolarAzimuthal::sample(uint64_t* seed) const
std::pair<Direction, double> PolarAzimuthal::sample(uint64_t* seed) const
{
return sample_impl(seed, false);
}
std::pair<Direction, double> PolarAzimuthal::sample_as_bias(
uint64_t* seed) const
{
return sample_impl(seed, true);
}
std::pair<Direction, double> PolarAzimuthal::sample_impl(
uint64_t* seed, bool return_pdf) const
{
// Sample cosine of polar angle
double mu = mu_->sample(seed);
if (mu == 1.0)
return u_ref_;
if (mu == -1.0)
return -u_ref_;
auto [mu, mu_wgt] = mu_->sample(seed);
// Sample azimuthal angle
double phi = phi_->sample(seed);
auto [phi, phi_wgt] = phi_->sample(seed);
// Compute either the PDF value or the importance weight
double weight =
return_pdf ? (mu_->evaluate(mu) * phi_->evaluate(phi)) : (mu_wgt * phi_wgt);
if (mu == 1.0)
return {u_ref_, weight};
if (mu == -1.0)
return {-u_ref_, weight};
double f = std::sqrt(1 - mu * mu);
return mu * u_ref_ + f * std::cos(phi) * v_ref_ + f * std::sin(phi) * w_ref_;
return {mu * u_ref_ + f * std::cos(phi) * v_ref_ + f * std::sin(phi) * w_ref_,
weight};
}
//==============================================================================
// Isotropic implementation
//==============================================================================
Isotropic::Isotropic(pugi::xml_node node) : UnitSphereDistribution {node}
{
if (check_for_node(node, "bias")) {
pugi::xml_node bias_node = node.child("bias");
std::string bias_type = get_node_value(bias_node, "type", true, true);
if (bias_type != "mu-phi") {
openmc::fatal_error(
"Isotropic distributions may only be biased by a PolarAzimuthal.");
}
auto bias = std::make_unique<PolarAzimuthal>(bias_node);
if (bias->mu()->bias() || bias->phi()->bias()) {
openmc::fatal_error(
"Attempted to bias Isotropic distribution with a biased PolarAzimuthal "
"distribution. Please ensure bias distributions are unbiased.");
}
this->set_bias(std::move(bias));
}
}
Direction isotropic_direction(uint64_t* seed)
{
double phi = uniform_distribution(0., 2.0 * PI, seed);
@ -111,18 +147,23 @@ Direction isotropic_direction(uint64_t* seed)
std::sqrt(1.0 - mu * mu) * std::sin(phi)};
}
Direction Isotropic::sample(uint64_t* seed) const
std::pair<Direction, double> Isotropic::sample(uint64_t* seed) const
{
return isotropic_direction(seed);
if (bias()) {
auto [val, eval] = bias()->sample_as_bias(seed);
return {val, 1.0 / (4.0 * PI * eval)};
} else {
return {isotropic_direction(seed), 1.0};
}
}
//==============================================================================
// Monodirectional implementation
//==============================================================================
Direction Monodirectional::sample(uint64_t* seed) const
std::pair<Direction, double> Monodirectional::sample(uint64_t* seed) const
{
return u_ref_;
return {u_ref_, 1.0};
}
} // namespace openmc

106
src/distribution_spatial.cpp
View file

@ -80,9 +80,13 @@ CartesianIndependent::CartesianIndependent(pugi::xml_node node)
}
}
Position CartesianIndependent::sample(uint64_t* seed) const
std::pair<Position, double> CartesianIndependent::sample(uint64_t* seed) const
{
return {x_->sample(seed), y_->sample(seed), z_->sample(seed)};
auto [x_val, x_wgt] = x_->sample(seed);
auto [y_val, y_wgt] = y_->sample(seed);
auto [z_val, z_wgt] = z_->sample(seed);
Position xi {x_val, y_val, z_val};
return {xi, x_wgt * y_wgt * z_wgt};
}
//==============================================================================
@ -139,14 +143,16 @@ CylindricalIndependent::CylindricalIndependent(pugi::xml_node node)
}
}
Position CylindricalIndependent::sample(uint64_t* seed) const
std::pair<Position, double> CylindricalIndependent::sample(uint64_t* seed) const
{
double r = r_->sample(seed);
double phi = phi_->sample(seed);
auto [r, r_wgt] = r_->sample(seed);
auto [phi, phi_wgt] = phi_->sample(seed);
auto [z, z_wgt] = z_->sample(seed);
double x = r * cos(phi) + origin_.x;
double y = r * sin(phi) + origin_.y;
double z = z_->sample(seed) + origin_.z;
return {x, y, z};
z += origin_.z;
Position xi {x, y, z};
return {xi, r_wgt * phi_wgt * z_wgt};
}
//==============================================================================
@ -203,16 +209,17 @@ SphericalIndependent::SphericalIndependent(pugi::xml_node node)
}
}
Position SphericalIndependent::sample(uint64_t* seed) const
std::pair<Position, double> SphericalIndependent::sample(uint64_t* seed) const
{
double r = r_->sample(seed);
double cos_theta = cos_theta_->sample(seed);
double phi = phi_->sample(seed);
auto [r, r_wgt] = r_->sample(seed);
auto [cos_theta, cos_theta_wgt] = cos_theta_->sample(seed);
auto [phi, phi_wgt] = phi_->sample(seed);
// sin(theta) by sin**2 + cos**2 = 1
double x = r * std::sqrt(1 - cos_theta * cos_theta) * cos(phi) + origin_.x;
double y = r * std::sqrt(1 - cos_theta * cos_theta) * sin(phi) + origin_.y;
double z = r * cos_theta + origin_.z;
return {x, y, z};
Position xi {x, y, z};
return {xi, r_wgt * cos_theta_wgt * phi_wgt};
}
//==============================================================================
@ -260,6 +267,37 @@ MeshSpatial::MeshSpatial(pugi::xml_node node)
}
elem_idx_dist_.assign(strengths);
if (check_for_node(node, "bias")) {
pugi::xml_node bias_node = node.child("bias");
if (check_for_node(bias_node, "strengths")) {
std::vector<double> bias_strengths(n_bins, 1.0);
bias_strengths = get_node_array<double>(node, "strengths");
if (bias_strengths.size() != n_bins) {
fatal_error(
fmt::format("Number of entries in the bias strengths array {} does "
"not match the number of entities in mesh {} ({}).",
bias_strengths.size(), mesh_id, n_bins));
}
if (get_node_value_bool(node, "volume_normalized")) {
for (int i = 0; i < n_bins; i++) {
bias_strengths[i] *= this->mesh()->volume(i);
}
}
// Compute importance weights
weight_ = compute_importance_weights(strengths, bias_strengths);
// Re-initialize DiscreteIndex with bias strengths for sampling
elem_idx_dist_.assign(bias_strengths);
} else {
fatal_error(fmt::format(
"Bias node for mesh {} found without strengths array.", mesh_id));
}
}
}
MeshSpatial::MeshSpatial(int32_t mesh_idx, span<const double> strengths)
@ -295,9 +333,11 @@ std::pair<int32_t, Position> MeshSpatial::sample_mesh(uint64_t* seed) const
return {elem_idx, mesh()->sample_element(elem_idx, seed)};
}
Position MeshSpatial::sample(uint64_t* seed) const
std::pair<Position, double> MeshSpatial::sample(uint64_t* seed) const
{
return this->sample_mesh(seed).second;
auto [elem_idx, u] = this->sample_mesh(seed);
double wgt = weight_.empty() ? 1.0 : weight_[elem_idx];
return {u, wgt};
}
//==============================================================================
@ -328,6 +368,31 @@ PointCloud::PointCloud(pugi::xml_node node)
}
point_idx_dist_.assign(strengths);
if (check_for_node(node, "bias")) {
pugi::xml_node bias_node = node.child("bias");
if (check_for_node(bias_node, "strengths")) {
std::vector<double> bias_strengths(point_cloud_.size(), 1.0);
bias_strengths = get_node_array<double>(node, "strengths");
if (bias_strengths.size() != point_cloud_.size()) {
fatal_error(
fmt::format("Number of entries in the bias strengths array {} does "
"not match the number of spatial points provided {}.",
bias_strengths.size(), point_cloud_.size()));
}
// Compute importance weights
weight_ = compute_importance_weights(strengths, bias_strengths);
// Re-initialize DiscreteIndex with bias strengths for sampling
point_idx_dist_.assign(bias_strengths);
} else {
fatal_error(
fmt::format("Bias node for PointCloud found without strengths array."));
}
}
}
PointCloud::PointCloud(
@ -337,10 +402,11 @@ PointCloud::PointCloud(
point_idx_dist_.assign(strengths);
}
Position PointCloud::sample(uint64_t* seed) const
std::pair<Position, double> PointCloud::sample(uint64_t* seed) const
{
int32_t index = point_idx_dist_.sample(seed);
return point_cloud_[index];
double wgt = weight_.empty() ? 1.0 : weight_[index];
return {point_cloud_[index], wgt};
}
//==============================================================================
@ -360,10 +426,10 @@ SpatialBox::SpatialBox(pugi::xml_node node, bool fission)
upper_right_ = Position {params[3], params[4], params[5]};
}
Position SpatialBox::sample(uint64_t* seed) const
std::pair<Position, double> SpatialBox::sample(uint64_t* seed) const
{
Position xi {prn(seed), prn(seed), prn(seed)};
return lower_left_ + xi * (upper_right_ - lower_left_);
return {lower_left_ + xi * (upper_right_ - lower_left_), 1.0};
}
//==============================================================================
@ -382,9 +448,9 @@ SpatialPoint::SpatialPoint(pugi::xml_node node)
r_ = Position {params.data()};
}
Position SpatialPoint::sample(uint64_t* seed) const
std::pair<Position, double> SpatialPoint::sample(uint64_t* seed) const
{
return r_;
return {r_, 1.0};
}
} // namespace openmc

4
src/secondary_correlated.cpp
View file

@ -247,9 +247,9 @@ void CorrelatedAngleEnergy::sample(
// Find correlated angular distribution for closest outgoing energy bin
if (r1 - c_k < c_k1 - r1 ||
distribution_[l].interpolation == Interpolation::histogram) {
mu = distribution_[l].angle[k]->sample(seed);
mu = distribution_[l].angle[k]->sample(seed).first;
} else {
mu = distribution_[l].angle[k + 1]->sample(seed);
mu = distribution_[l].angle[k + 1]->sample(seed).first;
}
}

24
src/source.cpp
View file

@ -356,6 +356,8 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
{
SourceSite site;
site.particle = particle_;
double r_wgt = 1.0;
double E_wgt = 1.0;
// Repeat sampling source location until a good site has been accepted
bool accepted = false;
@ -365,7 +367,9 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
while (!accepted) {
// Sample spatial distribution
site.r = space_->sample(seed);
auto [r, r_wgt_temp] = space_->sample(seed);
site.r = r;
r_wgt = r_wgt_temp;
// Check if sampled position satisfies spatial constraints
accepted = satisfies_spatial_constraints(site.r);
@ -378,7 +382,10 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
}
// Sample angle
site.u = angle_->sample(seed);
auto [u, u_wgt] = angle_->sample(seed);
site.u = u;
site.wgt = r_wgt * u_wgt;
// Sample energy and time for neutron and photon sources
if (settings::solver_type != SolverType::RANDOM_RAY) {
@ -395,7 +402,9 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
while (true) {
// Sample energy spectrum
site.E = energy_->sample(seed);
auto [E, E_wgt_temp] = energy_->sample(seed);
site.E = E;
E_wgt = E_wgt_temp;
// Resample if energy falls above maximum particle energy
if (site.E < data::energy_max[p] &&
@ -407,7 +416,10 @@ SourceSite IndependentSource::sample(uint64_t* seed) const
}
// Sample particle creation time
site.time = time_->sample(seed);
auto [time, time_wgt] = time_->sample(seed);
site.time = time;
site.wgt *= (E_wgt * time_wgt);
}
// Increment number of accepted samples
@ -537,9 +549,9 @@ CompiledSourceWrapper::~CompiledSourceWrapper()
// MeshElementSpatial implementation
//==============================================================================
Position MeshElementSpatial::sample(uint64_t* seed) const
std::pair<Position, double> MeshElementSpatial::sample(uint64_t* seed) const
{
return model::meshes[mesh_index_]->sample_element(elem_index_, seed);
return {model::meshes[mesh_index_]->sample_element(elem_index_, seed), 1.0};
}
//==============================================================================

4
tests/cpp_unit_tests/test_distribution.cpp
View file

@ -28,7 +28,7 @@ TEST_CASE("Test alias method sampling of a discrete distribution")
int counter = 0;
for (size_t i = 0; i < n_samples; i++) {
auto sample = dist.sample(&seed);
auto sample = dist.sample(&seed).first;
std += sample * sample / n_samples;
dist_mean += sample;
@ -61,7 +61,7 @@ TEST_CASE("Test alias sampling method for pugixml constructor")
// Initialize discrete distribution and seed
openmc::Discrete dist(energy);
uint64_t seed = openmc::init_seed(0, 0);
auto sample = dist.sample(&seed);
auto sample = dist.sample(&seed).first;
// Assertions
REQUIRE(dist.x().size() == 3);

2
tests/regression_tests/source/results_true.dat
View file

@ -1,2 +1,2 @@
k-combined:
3.080655E-01 4.837707E-03
2.942254E-01 2.571435E-03

9
tests/unit_tests/__init__.py
View file

@ -14,3 +14,12 @@ def assert_unbounded(obj):
ll, ur = obj.bounding_box
assert ll == pytest.approx((-np.inf, -np.inf, -np.inf))
assert ur == pytest.approx((np.inf, np.inf, np.inf))
def assert_sample_mean(samples, expected_mean):
# Calculate sample standard deviation
std_dev = samples.std() / np.sqrt(samples.size - 1)
# Means should agree within 4 sigma 99.993% of the time. Note that this is
# expected to fail about 1 out of 16,000 times
assert np.abs(expected_mean - samples.mean()) < 4*std_dev

276
tests/unit_tests/test_source_biasing.py Normal file
View file

@ -0,0 +1,276 @@
"""Tests for source biasing using C++ sampling routines via openmc.lib
This test module validates that the C++ distribution sampling implementations
correctly handle both unbiased and biased sampling when used in source
definitions. Each test:
1. Creates a minimal model with a source using a specific energy distribution
2. Uses model.sample_external_source() to generate samples via openmc.lib
3. Extracts energies from the returned particle list
4. Validates that:
- Unbiased sampling produces the expected mean
- Biased sampling with importance weighting produces the expected mean
- Weights are correctly applied (non-unity for biased case)
These tests complement the Python-level tests in test_stats.py by exercising
the full C++ sampling codepath that is used during actual simulations.
"""
import numpy as np
import pytest
import openmc
from tests.unit_tests import assert_sample_mean
@pytest.fixture
def model():
"""Create a minimal model for source sampling tests."""
sphere = openmc.Sphere(r=100.0, boundary_type='vacuum')
cell = openmc.Cell(region=-sphere)
geometry = openmc.Geometry([cell])
settings = openmc.Settings(particles=100, batches=1)
space = openmc.stats.Point()
angle = openmc.stats.Monodirectional((1.0, 0.0, 0.0))
settings.source = openmc.IndependentSource(space=space, angle=angle)
return openmc.Model(geometry=geometry, settings=settings)
@pytest.mark.flaky(reruns=1)
def test_discrete(run_in_tmpdir, model):
"""Test Discrete distribution sampling via C++ routines."""
vals = np.array([1.0, 2.0, 3.0])
probs = np.array([0.1, 0.7, 0.2])
exp_mean = (vals * probs).sum()
# Create source with discrete energy distribution
model.settings.source[0].energy = energy_dist = openmc.stats.Discrete(vals, probs)
# Sample using C++ routines and extract energies
n_samples = 10_000
particles = model.sample_external_source(n_samples)
energies = np.array([p.E for p in particles])
# Check unbiased mean
assert_sample_mean(energies, exp_mean)
# Sample from biased distribution
energy_dist.bias = np.array([0.2, 0.1, 0.7])
particles = model.sample_external_source(n_samples)
energies = np.array([p.E for p in particles])
weights = np.array([p.wgt for p in particles])
# Check biased weighted mean
weighted_energies = energies * weights
assert_sample_mean(weighted_energies, exp_mean)
assert np.any(weights != 1.0)
@pytest.mark.flaky(reruns=1)
def test_uniform(run_in_tmpdir, model):
"""Test Uniform distribution sampling via C++ routines."""
a, b = 5.0, 10.0
exp_mean = 0.5 * (a + b)
# Create source with uniform energy distribution
model.settings.source[0].energy = energy_dist = openmc.stats.Uniform(a, b)
# Sample using C++ routines and extract energies
n_samples = 10_000
particles = model.sample_external_source(n_samples)
energies = np.array([p.E for p in particles])
# Check unbiased mean
assert_sample_mean(energies, exp_mean)
# Sample from biased distribution
energy_dist.bias = openmc.stats.PowerLaw(a, b, 2)
particles = model.sample_external_source(n_samples)
energies = np.array([p.E for p in particles])
weights = np.array([p.wgt for p in particles])
# Check biased weighted mean
weighted_energies = energies * weights
assert_sample_mean(weighted_energies, exp_mean)
assert np.any(weights != 1.0)
@pytest.mark.flaky(reruns=1)
def test_powerlaw(run_in_tmpdir, model):
"""Test PowerLaw distribution sampling via C++ routines."""
a, b, n = 1.0, 20.0, 2.0
# Determine mean of distribution
exp_mean = (n+1)*(b**(n+2) - a**(n+2))/((n+2)*(b**(n+1) - a**(n+1)))
# Create source with powerlaw energy distribution
model.settings.source[0].energy = energy_dist = openmc.stats.PowerLaw(a, b, n)
# Sample using C++ routines and extract energies
n_samples = 10_000
particles = model.sample_external_source(n_samples)
energies = np.array([p.E for p in particles])
# Check unbiased mean
assert_sample_mean(energies, exp_mean)
# Sample from biased distribution
energy_dist.bias = openmc.stats.Uniform(a, b)
particles = model.sample_external_source(n_samples)
energies = np.array([p.E for p in particles])
weights = np.array([p.wgt for p in particles])
# Check biased weighted mean
weighted_energies = energies * weights
assert_sample_mean(weighted_energies, exp_mean)
assert np.any(weights != 1.0)
@pytest.mark.flaky(reruns=1)
def test_maxwell(run_in_tmpdir, model):
"""Test Maxwell distribution sampling via C++ routines."""
theta = 1.2895e6
exp_mean = 3/2 * theta
# Create source with Maxwell energy distribution
model.settings.source[0].energy = energy_dist = openmc.stats.Maxwell(theta)
# Sample using C++ routines and extract energies
n_samples = 10_000
particles = model.sample_external_source(n_samples)
energies = np.array([p.E for p in particles])
# Check unbiased mean
assert_sample_mean(energies, exp_mean)
# Sample from biased distribution
energy_dist.bias = openmc.stats.Maxwell(theta * 1.1)
particles = model.sample_external_source(n_samples)
energies = np.array([p.E for p in particles])
weights = np.array([p.wgt for p in particles])
# Check biased weighted mean
weighted_energies = energies * weights
assert_sample_mean(weighted_energies, exp_mean)
assert np.any(weights != 1.0)
@pytest.mark.flaky(reruns=1)
def test_watt(run_in_tmpdir, model):
"""Test Watt distribution sampling via C++ routines."""
a, b = 0.965e6, 2.29e-6
exp_mean = 3/2 * a + a**2 * b / 4
# Create source with Watt energy distribution
model.settings.source[0].energy = energy_dist = openmc.stats.Watt(a, b)
# Sample using C++ routines and extract energies
n_samples = 10_000
particles = model.sample_external_source(n_samples)
energies = np.array([p.E for p in particles])
# Check unbiased mean
assert_sample_mean(energies, exp_mean)
# Sample from biased distribution
energy_dist.bias = openmc.stats.Watt(a*1.05, b)
particles = model.sample_external_source(n_samples)
energies = np.array([p.E for p in particles])
weights = np.array([p.wgt for p in particles])
# Check biased weighted mean
weighted_energies = energies * weights
assert_sample_mean(weighted_energies, exp_mean)
assert np.any(weights != 1.0)
@pytest.mark.flaky(reruns=1)
def test_tabular(run_in_tmpdir, model):
"""Test Tabular distribution sampling via C++ routines."""
# Test linear-linear sampling
x = np.array([0.0, 5.0, 7.0, 10.0])
p = np.array([10.0, 20.0, 5.0, 6.0])
# Create tabular distribution and normalize to get expected mean
model.settings.source[0].energy = energy_dist = openmc.stats.Tabular(x, p, 'linear-linear')
energy_dist.normalize()
exp_mean = energy_dist.mean()
# Sample using C++ routines and extract energies
n_samples = 10_000
particles = model.sample_external_source(n_samples)
energies = np.array([p.E for p in particles])
# Check unbiased mean
assert_sample_mean(energies, exp_mean)
# Sample from biased distribution
energy_dist.bias = openmc.stats.Uniform(x[0], x[-1])
particles = model.sample_external_source(n_samples)
energies = np.array([p.E for p in particles])
weights = np.array([p.wgt for p in particles])
# Check biased weighted mean
weighted_energies = energies * weights
assert_sample_mean(weighted_energies, exp_mean)
assert np.any(weights != 1.0)
@pytest.mark.flaky(reruns=1)
def test_mixture(run_in_tmpdir, model):
"""Test Mixture distribution sampling via C++ routines."""
d1 = openmc.stats.Uniform(0, 5)
d2 = openmc.stats.Uniform(3, 7)
p = [0.5, 0.5]
# Create mixture energy distribution
model.settings.source[0].energy = energy_dist = openmc.stats.Mixture(p, [d1, d2])
exp_mean = (2.5 + 5.0) / 2
# Sample using C++ routines and extract energies
n_samples = 10_000
particles = model.sample_external_source(n_samples)
energies = np.array([p.E for p in particles])
# Check unbiased mean
assert_sample_mean(energies, exp_mean)
# Sample using biased sub-distribution
energy_dist.distribution[0].bias = openmc.stats.PowerLaw(0, 5, 2)
particles = model.sample_external_source(n_samples)
energies = np.array([p.E for p in particles])
weights = np.array([p.wgt for p in particles])
# Check biased weighted mean
weighted_energies = energies * weights
assert_sample_mean(weighted_energies, exp_mean)
assert np.any(weights != 1.0)
@pytest.mark.flaky(reruns=1)
def test_normal(run_in_tmpdir, model):
"""Test Normal distribution sampling via C++ routines."""
mean_val = 25.0
std_dev = 2.0
# Create source with normal energy distribution
model.settings.source[0].energy = energy_dist = openmc.stats.Normal(mean_val, std_dev)
# Sample using C++ routines and extract energies
n_samples = 10_000
particles = model.sample_external_source(n_samples)
energies = np.array([p.E for p in particles])
# Check unbiased mean
assert_sample_mean(energies, mean_val)
# Sample from biased distribution
energy_dist.bias = openmc.stats.Normal(mean_val * 1.1, std_dev)
particles = model.sample_external_source(n_samples)
energies = np.array([p.E for p in particles])
weights = np.array([p.wgt for p in particles])
# Check biased weighted mean
weighted_energies = energies * weights
assert_sample_mean(weighted_energies, mean_val)
assert np.any(weights != 1.0)

191
tests/unit_tests/test_stats.py
View file

@ -6,14 +6,7 @@ import openmc
import openmc.stats
from scipy.integrate import trapezoid
def assert_sample_mean(samples, expected_mean):
# Calculate sample standard deviation
std_dev = samples.std() / np.sqrt(samples.size - 1)
# Means should agree within 4 sigma 99.993% of the time. Note that this is
# expected to fail about 1 out of 16,000 times
assert np.abs(expected_mean - samples.mean()) < 4*std_dev
from tests.unit_tests import assert_sample_mean
@pytest.mark.flaky(reruns=1)
@ -47,8 +40,19 @@ def test_discrete():
# sample discrete distribution and check that the mean of the samples is
# within 4 std. dev. of the expected mean
n_samples = 1_000_000
samples = d3.sample(n_samples)
samples, weights = d3.sample(n_samples)
assert_sample_mean(samples, exp_mean)
assert np.all(weights == 1.0)
# Test biased distribution
d3.bias = np.array([0.2, 0.1, 0.7])
bias_elem = d3.to_xml_element('distribution')
d4 = openmc.stats.Univariate.from_xml_element(bias_elem)
np.testing.assert_array_equal(d4.bias, [0.2, 0.1, 0.7])
samples, weights = d4.sample(n_samples)
weighted_sample = samples * weights
assert_sample_mean(weighted_sample, exp_mean)
assert np.all(weights != 1.0)
def test_delta_function():
@ -83,6 +87,28 @@ def test_merge_discrete():
assert triple.integral() == pytest.approx(6.0)
def test_merge_discrete_with_bias():
# Two discrete distributions with different biases
d1 = openmc.stats.Discrete([1.0, 2.0], [0.5, 0.5])
d2 = openmc.stats.Discrete([2.0, 3.0], [0.3, 0.7], bias=[0.1, 0.9])
merged = openmc.stats.Discrete.merge([d1, d2], [0.6, 0.4])
exp_mean = 0.6 * d1.mean() + 0.4 * d2.mean()
# Verify merged distribution has correct x values
assert set(merged.x) == {1.0, 2.0, 3.0}
# Bias should not be changed in original distributions
assert d1.bias is None
assert np.all(d2.bias == [0.1, 0.9])
# Sample and verify bias is applied correctly
samples, weights = merged.sample(10_000)
# Verify weighted mean matches expected unbiased mean
assert_sample_mean(samples*weights, exp_mean)
def test_clip_discrete():
# Create discrete distribution with two points that are not important, one
# because the x value is very small, and one because the p value is very
@ -125,8 +151,19 @@ def test_uniform():
# std. dev. of the expected mean
exp_mean = 0.5 * (a + b)
n_samples = 1_000_000
samples = d.sample(n_samples)
samples, weights = d.sample(n_samples)
assert_sample_mean(samples, exp_mean)
assert np.all(weights == 1.0)
# Test biased distribution
d.bias = openmc.stats.PowerLaw(a, b, 2)
bias_elem = d.to_xml_element('distribution')
d2 = openmc.stats.Univariate.from_xml_element(bias_elem)
assert isinstance (d2.bias, openmc.stats.PowerLaw)
samples, weights = d2.sample(n_samples)
weighted_sample = samples * weights
assert_sample_mean(weighted_sample, exp_mean)
assert np.all(weights != 1.0)
@pytest.mark.flaky(reruns=1)
@ -147,8 +184,19 @@ def test_powerlaw():
# sample power law distribution and check that the mean of the samples is
# within 4 std. dev. of the expected mean
n_samples = 1_000_000
samples = d.sample(n_samples)
samples, weights = d.sample(n_samples)
assert_sample_mean(samples, exp_mean)
assert np.all(weights == 1.0)
# Test biased distribution
d.bias = openmc.stats.Uniform(a, b)
bias_elem = d.to_xml_element('distribution')
d2 = openmc.stats.Univariate.from_xml_element(bias_elem)
assert isinstance (d2.bias, openmc.stats.Uniform)
samples, weights = d2.sample(n_samples)
weighted_sample = samples * weights
assert_sample_mean(weighted_sample, exp_mean)
assert np.all(weights != 1.0)
@pytest.mark.flaky(reruns=1)
@ -166,13 +214,25 @@ def test_maxwell():
# sample maxwell distribution and check that the mean of the samples is
# within 4 std. dev. of the expected mean
n_samples = 1_000_000
samples = d.sample(n_samples)
samples, weights = d.sample(n_samples)
assert_sample_mean(samples, exp_mean)
assert np.all(weights == 1.0)
# A second sample starting from a different seed
samples_2 = d.sample(n_samples)
samples_2, weights_2 = d.sample(n_samples)
assert_sample_mean(samples_2, exp_mean)
assert samples_2.mean() != samples.mean()
assert np.all(weights_2 == 1)
# Test biased distribution
d.bias = openmc.stats.Maxwell((theta * 1.1))
bias_elem = d.to_xml_element('distribution')
d2 = openmc.stats.Univariate.from_xml_element(bias_elem)
assert isinstance (d2.bias, openmc.stats.Maxwell)
samples, weights = d2.sample(n_samples)
weighted_sample = samples * weights
assert_sample_mean(weighted_sample, exp_mean)
assert np.all(weights != 1.0)
@pytest.mark.flaky(reruns=1)
@ -195,8 +255,19 @@ def test_watt():
# sample Watt distribution and check that the mean of the samples is within
# 4 std. dev. of the expected mean
n_samples = 1_000_000
samples = d.sample(n_samples)
samples, weights = d.sample(n_samples)
assert_sample_mean(samples, exp_mean)
assert np.all(weights == 1.0)
# Test biased distribution with 5 percent higher T_e
d.bias = openmc.stats.Watt(a*1.05, b)
bias_elem = d.to_xml_element('distribution')
d2 = openmc.stats.Univariate.from_xml_element(bias_elem)
assert isinstance (d2.bias, openmc.stats.Watt)
samples, weights = d2.sample(n_samples)
weighted_sample = samples * weights
assert_sample_mean(weighted_sample, exp_mean)
assert np.all(weights != 1.0)
@pytest.mark.flaky(reruns=1)
@ -206,8 +277,9 @@ def test_tabular():
p = np.array([10.0, 20.0, 5.0, 6.0])
d = openmc.stats.Tabular(x, p, 'linear-linear')
n_samples = 100_000
samples = d.sample(n_samples)
samples, weights = d.sample(n_samples)
assert_sample_mean(samples, d.mean())
assert np.all(weights == 1.0)
# test linear-linear normalization
d.normalize()
@ -215,8 +287,9 @@ def test_tabular():
# test histogram sampling
d = openmc.stats.Tabular(x, p, interpolation='histogram')
samples = d.sample(n_samples)
samples, weights = d.sample(n_samples)
assert_sample_mean(samples, d.mean())
assert np.all(weights == 1.0)
d.normalize()
assert d.integral() == pytest.approx(1.0)
@ -226,7 +299,7 @@ def test_tabular():
d = openmc.stats.Tabular(x, p[:-1], interpolation='histogram')
d.cdf()
d.mean()
assert_sample_mean(d.sample(n_samples), d.mean())
assert_sample_mean(d.sample(n_samples)[0], d.mean())
# passing a shorter probability set should raise an error for linear-linear
with pytest.raises(ValueError):
@ -238,6 +311,16 @@ def test_tabular():
d = openmc.stats.Tabular(x, p, interpolation='linear-linear')
d.cdf()
# Test biased distribution
d.bias = openmc.stats.Uniform(x[0], x[-1])
bias_elem = d.to_xml_element('distribution')
d2 = openmc.stats.Univariate.from_xml_element(bias_elem)
assert isinstance (d2.bias, openmc.stats.Uniform)
samples, weights = d2.sample(n_samples)
weighted_sample = samples * weights
assert_sample_mean(weighted_sample, d2.mean())
assert np.all(weights != 1.0)
def test_tabular_from_xml():
x = np.array([0.0, 5.0, 7.0, 10.0])
@ -288,8 +371,9 @@ def test_mixture():
# Sample and make sure sample mean is close to expected mean
n_samples = 1_000_000
samples = mix.sample(n_samples)
samples, weights = mix.sample(n_samples)
assert_sample_mean(samples, (2.5 + 5.0)/2)
assert np.all(weights == 1.0)
elem = mix.to_xml_element('distribution')
@ -298,6 +382,26 @@ def test_mixture():
assert d.distribution == [d1, d2]
assert len(d) == 4
# Test biased sub-distribution
d.distribution[0].bias = openmc.stats.PowerLaw(0, 5, 2)
bias_elem = d.to_xml_element('distribution')
d3 = openmc.stats.Univariate.from_xml_element(bias_elem)
assert isinstance (d3.distribution[0].bias, openmc.stats.PowerLaw)
samples, weights = d3.sample(n_samples)
weighted_sample = samples * weights
assert_sample_mean(weighted_sample, (2.5 + 5.0)/2)
# Test biased meta-probability
d.distribution[0].bias = None
d.bias = [0.25, 0.75]
bias_elem_2 = d.to_xml_element('distribution')
d4 = openmc.stats.Univariate.from_xml_element(bias_elem_2)
assert isinstance (d4.bias, np.ndarray)
samples, weights = d4.sample(n_samples)
weighted_sample = samples * weights
assert_sample_mean(weighted_sample, (2.5 + 5.0)/2)
assert np.all(weights != 1.0)
def test_mixture_clip():
# Create mixture distribution containing a discrete distribution with two
@ -330,6 +434,13 @@ def test_mixture_clip():
with pytest.warns(UserWarning):
mix_clip = mix.clip(1e-6)
# Make sure warning is raised if a biased Discrete is clipped
d3 = openmc.stats.Discrete([1.0, 1.001], [1.0, 0.7e-8])
d3.bias = [0.9, 0.1]
mix = openmc.stats.Mixture([1.0, 1.0], [d3, d2])
with pytest.raises(RuntimeError):
mix_clip = mix.clip(1e-6)
def test_polar_azimuthal():
# default polar-azimuthal should be uniform in mu and phi
@ -365,12 +476,17 @@ def test_polar_azimuthal():
def test_isotropic():
d = openmc.stats.Isotropic()
mu = openmc.stats.Uniform(-1.0, 1.0)
phi = openmc.stats.PowerLaw(0., 2*np.pi, 2)
d2 = openmc.stats.PolarAzimuthal(mu, phi)
d.bias = d2
elem = d.to_xml_element()
assert elem.tag == 'angle'
assert elem.attrib['type'] == 'isotropic'
d = openmc.stats.Isotropic.from_xml_element(elem)
assert isinstance(d, openmc.stats.Isotropic)
assert isinstance(d.bias, openmc.stats.PolarAzimuthal)
def test_monodirectional():
@ -387,6 +503,7 @@ def test_cartesian():
x = openmc.stats.Uniform(-10., 10.)
y = openmc.stats.Uniform(-10., 10.)
z = openmc.stats.Uniform(0., 20.)
z.bias = openmc.stats.PowerLaw(0., 20., 3)
d = openmc.stats.CartesianIndependent(x, y, z)
elem = d.to_xml_element()
@ -402,6 +519,7 @@ def test_cartesian():
d = openmc.stats.Spatial.from_xml_element(elem)
assert isinstance(d, openmc.stats.CartesianIndependent)
assert isinstance (d.z.bias, openmc.stats.PowerLaw)
def test_box():
@ -448,8 +566,19 @@ def test_normal():
# sample normal distribution
n_samples = 100_000
samples = d.sample(n_samples)
samples, weights = d.sample(n_samples)
assert_sample_mean(samples, mean)
assert np.all(weights == 1.0)
# Test biased distribution
d.bias = openmc.stats.Normal(10.0, 4.0)
bias_elem = d.to_xml_element('distribution')
d2 = openmc.stats.Univariate.from_xml_element(bias_elem)
assert isinstance (d2.bias, openmc.stats.Normal)
samples, weights = d2.sample(n_samples)
weighted_sample = samples * weights
assert_sample_mean(weighted_sample, mean)
assert np.all(weights != 1.0)
@pytest.mark.flaky(reruns=1)
@ -468,8 +597,9 @@ def test_muir():
# sample muir distribution
n_samples = 100_000
samples = d.sample(n_samples)
samples, weights = d.sample(n_samples)
assert_sample_mean(samples, mean)
assert np.all(weights == 1.0)
@pytest.mark.flaky(reruns=1)
@ -514,8 +644,27 @@ def test_combine_distributions():
# Sample the combined distribution and make sure the sample mean is within
# uncertainty of the expected value
samples = combined.sample(10_000)
samples, weights = combined.sample(10_000)
assert_sample_mean(samples, 0.25)
assert np.all(weights == 1.0)
# If biased/unbiased Discrete distributions are combined, unbiased probability
# should be conserved and points from both original distributions should be
# assigned bias probabilities.
x1 = [0.0, 1.0, 10.0]
p1 = [0.3, 0.2, 0.5]
b1 = [0.2, 0.5, 0.3]
d1 = openmc.stats.Discrete(x1, p1, b1)
x2 = [0.5, 1.0, 5.0]
p2 = [0.4, 0.5, 0.1]
d2 = openmc.stats.Discrete(x2, p2)
combined = openmc.stats.combine_distributions([d1, d2, t1], [0.25, 0.25, 0.5])
p3 = [0.075, 0.1, 0.175, 0.025, 0.125]
b3 = [0.05, 0.1, 0.25, 0.025, 0.075]
assert all(combined.distribution[-1].p == p3)
assert all(combined.distribution[-1].bias == b3)
def test_reference_vwu_projection():
"""When a non-orthogonal vector is provided, the setter should project out