diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 6ea8cfc56..b76b23e63 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -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: diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index dcd31f429..30640d07b 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -243,13 +243,17 @@ public: //! \return Sampled value double sample(uint64_t* seed) const; - // d and p property - const vector& d() const { return d_; } - const vector& p() const { return p_; } private: - vector d_; //!< Pointer to sub-distributions - vector p_; //!< tabulated probability density - vector 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 distribution_; //!< sub-distributions + cummulative probabilities }; diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index c7b6dadbe..b5997b1cd 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -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) diff --git a/src/distribution.cpp b/src/distribution.cpp index f473513ba..4b52a35e0 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -297,30 +297,43 @@ Mixture::Mixture(pugi::xml_node node) //auto params = get_node_array(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); } //============================================================================== diff --git a/tests/regression_tests/source/inputs_true.dat b/tests/regression_tests/source/inputs_true.dat index de81a4179..eb60821c6 100644 --- a/tests/regression_tests/source/inputs_true.dat +++ b/tests/regression_tests/source/inputs_true.dat @@ -76,4 +76,27 @@ 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 + + + + + + -2.0 0.0 2.0 0.2 0.3 0.2 + + + + + + + + + + + + + 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 + + + + diff --git a/tests/regression_tests/source/results_true.dat b/tests/regression_tests/source/results_true.dat index 66b610b3a..1815324f6 100644 --- a/tests/regression_tests/source/results_true.dat +++ b/tests/regression_tests/source/results_true.dat @@ -1,2 +1,2 @@ k-combined: -3.033600E-01 1.676108E-03 +2.980096E-01 9.632798E-04 diff --git a/tests/regression_tests/source/test.py b/tests/regression_tests/source/test.py index 91316acc6..229ac8383 100644 --- a/tests/regression_tests/source/test.py +++ b/tests/regression_tests/source/test.py @@ -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()