crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 13:15:39 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Update mixture distribution

Browse source 
update storeage of c++ class
  update sampling -> uses std::lower_bound for O(log(N)) sampling of the
distribution
  implemented from_xml_element
  updated settings.rst
  include a mixture in the regression tests
This commit is contained in:
Olaf Schumann 2021-09-23 09:34:11 +00:00
parent 1c219600a1
commit e7fa663b3d
7 changed files with 98 additions and 30 deletions
Download patch file Download diff file Expand all files Collapse all files

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

@ -619,16 +619,17 @@ variable and whose sub-elements/attributes are as follows:
:type:
The type of the distribution. Valid options are "uniform", "discrete",
"tabular", "maxwell", and "watt". The "uniform" option produces variates
sampled from a uniform distribution over a finite interval. The "discrete"
option produces random variates that can assume a finite number of values
(i.e., a distribution characterized by a probability mass function). The
"tabular" option produces random variates sampled from a tabulated
"tabular", "maxwell", "watt", and "mixture". The "uniform" option producess
variates sampled from a uniform distribution over a finite interval. The
"discrete" option produces random variates that can assume a finite number
of values (i.e., a distribution characterized by a probability mass function).
The "tabular" option produces random variates sampled from a tabulated
distribution where the density function is either a histogram or
linearly-interpolated between tabulated points. The "watt" option produces
random variates is sampled from a Watt fission spectrum (only used for
energies). The "maxwell" option produce variates sampled from a Maxwell
fission spectrum (only used for energies).
fission spectrum (only used for energies). The "mixture" option produces samples
from univariate sub-distributions with given probabilities.
*Default*: None
@ -649,6 +650,11 @@ variable and whose sub-elements/attributes are as follows:
number :math:`a` that parameterizes the distribution :math:`p(x) dx = c x
e^{-x/a} dx`.
For a "mixture" distribution, ``parameters`` provide the :math:'(p,d)' pairs
connecting the probabilites :math:'p' with the different sub-distributions
:math:'d'. All probabilities :math:'p' are given first followed by the corresponding
distributions :math:'d'.
.. note:: The above format should be used even when using the multi-group
:ref:`energy_mode`.
:interpolation:

16
include/openmc/distribution.h
View file

@ -243,13 +243,17 @@ public:
//! \return Sampled value
double sample(uint64_t* seed) const;
// d and p property
const vector<UPtrDist>& d() const { return d_; }
const vector<double>& p() const { return p_; }
private:
vector<UPtrDist> d_; //!< Pointer to sub-distributions
vector<double> p_; //!< tabulated probability density
vector<double> c_; //!< cumulative distribution
struct Pair {
double cummulative_probability_ {0.0};
UPtrDist distribution_;
// bool operator<(const Pair& o) const {
// return cummulative_probability_ < o.cummulative_probability_;
//}
bool operator<(double p) const { return cummulative_probability_ < p; }
};
vector<Pair> distribution_; //!< sub-distributions + cummulative probabilities
};

22
openmc/stats/univariate.py
View file

@ -837,4 +837,24 @@ class Mixture(Univariate):
@classmethod
def from_xml_element(cls, elem):
raise NotImplementedError
"""Generate mixture distribution from an XML element
Parameters
----------
elem : xml.etree.ElementTree.Element
XML element
Returns
-------
openmc.stats.Mixture
Mixture distribution generated from XML element
"""
P = []
D = []
for pair in elem:
if pair.tag == "pair":
P.append( float(get_text(pair, 'probability')) )
D.append( Univariate.from_xml_element(pair.find("dist")) )
return cls(P,D)

43
src/distribution.cpp
View file

@ -297,30 +297,43 @@ Mixture::Mixture(pugi::xml_node node)
//auto params = get_node_array<double>(node, "parameters");
double cumsum = 0.0;
for (pugi::xml_node pair = node.child("pair"); pair; pair = pair.next_sibling("pair")) {
for (pugi::xml_node pair = node.child("pair"); pair;
pair = pair.next_sibling("pair")) {
// Check that required data exists
if (!pair.attribute("probability")) openmc::fatal_error("Mixture pair element does not have probability.");
if (!pair.child("dist")) openmc::fatal_error("Mixture pair element does not have a distribution.");
double p = std::stod(pair.attribute("probability").value());
p_.push_back(p);
cumsum += p;
c_.push_back(cumsum);
d_.push_back(distribution_from_xml(pair.child("dist")));
// cummulative sum of probybilities
cumsum += std::stod(pair.attribute("probability").value());
// Save cummulative probybility and distrubution
distribution_.push_back(
{cumsum, distribution_from_xml(pair.child("dist"))});
}
// Normalize cummulative probabilities to 1
for (auto& pair : distribution_) {
pair.cummulative_probability_ /= cumsum;
}
std::transform(c_.begin(), c_.end(), c_.begin(), [cumsum](double c){ return c/cumsum;});
}
double Mixture::sample(uint64_t* seed) const {
//
double Mixture::sample(uint64_t* seed) const
{
// Sample value of CDF
double c = prn(seed);
const double p = prn(seed);
// Find first CDF bin which is above the sampled value and sample corresponding sub-distribution
for (size_t i = 0; i < c_.size(); i++) {
if (c<=c_[i]) return d_[i]->sample(seed);
}
// find matching distribution
const auto it =
std::lower_bound(distribution_.cbegin(), distribution_.cend(), p);
// This should not happen. Catch it
openmc::fatal_error("Bad Sampling in Mixture Distribution.");
// TODO: Remove this check, when the code is trusted to always operate
if (it == distribution_.cend()) {
fatal_error("Bad Sampling in Mixture Distribution.");
}
// Sample the choosen distribution
return it->distribution_->sample(seed);
}
//==============================================================================

23
tests/regression_tests/source/inputs_true.dat
View file

@ -76,4 +76,27 @@
<parameters>1.0 1.3894954943731377 1.93069772888325 2.6826957952797255 3.72759372031494 5.17947467923121 7.196856730011519 10.0 13.894954943731374 19.306977288832496 26.826957952797247 37.2759372031494 51.7947467923121 71.96856730011518 100.0 138.94954943731375 193.06977288832496 268.26957952797244 372.7593720314938 517.9474679231207 719.6856730011514 1000.0 1389.4954943731375 1930.6977288832495 2682.6957952797247 3727.593720314938 5179.474679231207 7196.856730011514 10000.0 13894.95494373136 19306.977288832495 26826.95795279722 37275.93720314938 51794.74679231213 71968.56730011514 100000.0 138949.5494373136 193069.77288832495 268269.5795279722 372759.3720314938 517947.4679231202 719685.6730011514 1000000.0 1389495.494373136 1930697.7288832497 2682695.7952797217 3727593.720314938 5179474.679231202 7196856.730011513 10000000.0 0.0 2.9086439299358713e-08 5.80533561806147e-08 8.67817193689187e-08 1.1515347785771536e-07 1.4305204600565115e-07 1.7036278261198208e-07 1.9697346200185813e-07 2.227747351856934e-07 2.4766057919761985e-07 2.715287327665956e-07 2.9428111652990295e-07 3.1582423606228735e-07 3.360695660646056e-07 3.549339141332686e-07 3.723397626156721e-07 3.882155871468592e-07 4.024961505584776e-07 4.151227709522976e-07 4.260435628367196e-07 4.3521365033538783e-07 4.4259535159179273e-07 4.4815833361210174e-07 4.5187973690993757e-07 4.5374426944091084e-07 4.5374426944091084e-07 4.5187973690993757e-07 4.4815833361210174e-07 4.4259535159179273e-07 4.352136503353879e-07 4.2604356283671966e-07 4.1512277095229767e-07 4.0249615055847764e-07 3.8821558714685926e-07 3.723397626156722e-07 3.5493391413326864e-07 3.360695660646057e-07 3.158242360622874e-07 2.942811165299031e-07 2.715287327665957e-07 2.4766057919762e-07 2.2277473518569352e-07 1.9697346200185819e-07 1.7036278261198226e-07 1.4305204600565126e-07 1.1515347785771556e-07 8.678171936891881e-08 5.805335618061493e-08 2.9086439299358858e-08 5.559621115282002e-23</parameters>
</energy>
</source>
<source strength="0.1">
<space origin="1.0 1.0 0.0" type="cylindrical">
<r parameters="2.0 3.0" type="uniform" />
<phi parameters="0.0 6.283185307179586" type="uniform" />
<z interpolation="linear-linear" type="tabular">
<parameters>-2.0 0.0 2.0 0.2 0.3 0.2</parameters>
</z>
</space>
<angle type="isotropic" />
<energy type="mixture">
<pair probability="1">
<dist parameters="1289500.0" type="maxwell" />
</pair>
<pair probability="2">
<dist parameters="988000.0 2.249e-06" type="watt" />
</pair>
<pair probability="3">
<dist interpolation="histogram" type="tabular">
<parameters>1.0 1.3894954943731377 1.93069772888325 2.6826957952797255 3.72759372031494 5.17947467923121 7.196856730011519 10.0 13.894954943731374 19.306977288832496 26.826957952797247 37.2759372031494 51.7947467923121 71.96856730011518 100.0 138.94954943731375 193.06977288832496 268.26957952797244 372.7593720314938 517.9474679231207 719.6856730011514 1000.0 1389.4954943731375 1930.6977288832495 2682.6957952797247 3727.593720314938 5179.474679231207 7196.856730011514 10000.0 13894.95494373136 19306.977288832495 26826.95795279722 37275.93720314938 51794.74679231213 71968.56730011514 100000.0 138949.5494373136 193069.77288832495 268269.5795279722 372759.3720314938 517947.4679231202 719685.6730011514 1000000.0 1389495.494373136 1930697.7288832497 2682695.7952797217 3727593.720314938 5179474.679231202 7196856.730011513 10000000.0 0.0 2.9086439299358713e-08 5.80533561806147e-08 8.67817193689187e-08 1.1515347785771536e-07 1.4305204600565115e-07 1.7036278261198208e-07 1.9697346200185813e-07 2.227747351856934e-07 2.4766057919761985e-07 2.715287327665956e-07 2.9428111652990295e-07 3.1582423606228735e-07 3.360695660646056e-07 3.549339141332686e-07 3.723397626156721e-07 3.882155871468592e-07 4.024961505584776e-07 4.151227709522976e-07 4.260435628367196e-07 4.3521365033538783e-07 4.4259535159179273e-07 4.4815833361210174e-07 4.5187973690993757e-07 4.5374426944091084e-07 4.5374426944091084e-07 4.5187973690993757e-07 4.4815833361210174e-07 4.4259535159179273e-07 4.352136503353879e-07 4.2604356283671966e-07 4.1512277095229767e-07 4.0249615055847764e-07 3.8821558714685926e-07 3.723397626156722e-07 3.5493391413326864e-07 3.360695660646057e-07 3.158242360622874e-07 2.942811165299031e-07 2.715287327665957e-07 2.4766057919762e-07 2.2277473518569352e-07 1.9697346200185819e-07 1.7036278261198226e-07 1.4305204600565126e-07 1.1515347785771556e-07 8.678171936891881e-08 5.805335618061493e-08 2.9086439299358858e-08 5.559621115282002e-23</parameters>
</dist>
</pair>
</energy>
</source>
</settings>

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

@ -1,2 +1,2 @@
k-combined:
3.033600E-01 1.676108E-03
2.980096E-01 9.632798E-04

4
tests/regression_tests/source/test.py
View file

@ -55,18 +55,20 @@ class SourceTestHarness(PyAPITestHarness):
energy1 = openmc.stats.Maxwell(1.2895e6)
energy2 = openmc.stats.Watt(0.988e6, 2.249e-6)
energy3 = openmc.stats.Tabular(E, p, interpolation='histogram')
energy4 = openmc.stats.Mixture([1,2,3], [energy1, energy2, energy3])
source1 = openmc.Source(spatial1, angle1, energy1, strength=0.5)
source2 = openmc.Source(spatial2, angle2, energy2, strength=0.3)
source3 = openmc.Source(spatial3, angle3, energy3, strength=0.1)
source4 = openmc.Source(spatial4, angle3, energy3, strength=0.1)
source5 = openmc.Source(spatial5, angle3, energy3, strength=0.1)
source6 = openmc.Source(spatial5, angle3, energy4, strength=0.1)
settings = openmc.Settings()
settings.batches = 10
settings.inactive = 5
settings.particles = 1000
settings.source = [source1, source2, source3, source4, source5]
settings.source = [source1, source2, source3, source4, source5, source6]
settings.export_to_xml()