From eb0eac79c139049970c439ea79f31badbe548a34 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Tue, 21 Sep 2021 19:05:47 +0000 Subject: [PATCH] Rational univariate distribution, that sampled a random variable with a pdf of p = x^n Usefull for samplig homogenous cylindrical or spherical sources Simmilar to MCNPs "sp -21 n" functionality --- include/openmc/distribution.h | 25 +++++++++ openmc/stats/univariate.py | 97 +++++++++++++++++++++++++++++++++++ src/distribution.cpp | 28 ++++++++++ 3 files changed, 150 insertions(+) diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index 8cdd78f788..6497078e47 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -80,6 +80,31 @@ private: double b_; //!< Upper bound of distribution }; +//============================================================================== +//! Rational distribution over the interval [a,b] with exponent n +//============================================================================== + +class Rational : public Distribution { +public: + explicit Rational(pugi::xml_node node); + Rational(double a, double b, double n) : a_{a}, b_{b}, n_{n}, an_{std::pow(a, n+1)}, bn_{std::pow(b, n+1)} {}; + + //! Sample a value from the distribution + //! \param seed Pseudorandom number seed pointer + //! \return Sampled value + double sample(uint64_t* seed) const; + + double a() const { return a_; } + double b() const { return b_; } + double n() const { return n_; } +private: + double a_; //!< Lower bound of distribution + double b_; //!< Upper bound of distribution + double n_; //!< Exponent of distribution + double an_; //!< Lower bound of distribution + double bn_; //!< Upper bound of distribution +}; + //============================================================================== //! Maxwellian distribution of form c*sqrt(E)*exp(-E/theta) //============================================================================== diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index af6c9a6d12..cdc2f4c300 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -42,6 +42,8 @@ class Univariate(EqualityMixin, ABC): return Discrete.from_xml_element(elem) elif distribution == 'uniform': return Uniform.from_xml_element(elem) + elif distribution == 'rational': + return Rational.from_xml_element(elem) elif distribution == 'maxwell': return Maxwell.from_xml_element(elem) elif distribution == 'watt': @@ -242,6 +244,101 @@ class Uniform(Univariate): params = get_text(elem, 'parameters').split() return cls(*map(float, params)) +class Rational(Univariate): + """Distribution with power law probability over a finite interval [a,b] + + Parameters + ---------- + a : float, optional + Lower bound of the sampling interval. Defaults to zero. + b : float, optional + Upper bound of the sampling interval. Defaults to unity. + n : float, optional + power law exponent. Defaults to zero -> Uniform distribution + + Attributes + ---------- + a : float + Lower bound of the sampling interval + b : float + Upper bound of the sampling interval + n : float + Power law exponent + + """ + + def __init__(self, a=0.0, b=1.0, n=0): + self.a = a + self.b = b + self.n = n + + def __len__(self): + return 3 + + @property + def a(self): + return self._a + + @property + def b(self): + return self._b + + @property + def n(self): + return self._n + + @a.setter + def a(self, a): + cv.check_type('Uniform a', a, Real) + self._a = a + + @b.setter + def b(self, b): + cv.check_type('Uniform b', b, Real) + self._b = b + + @n.setter + def n(self, n): + cv.check_type('Uniform n', n, Real) + self._n = n + + def to_xml_element(self, element_name): + """Return XML representation of the rational distribution + + Parameters + ---------- + element_name : str + XML element name + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing uniform distribution data + + """ + element = ET.Element(element_name) + element.set("type", "rational") + element.set("parameters", '{} {} {}'.format(self.a, self.b, self.n)) + return element + + @classmethod + def from_xml_element(cls, elem): + """Generate rational distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.stats.Rational + Uniform distribution generated from XML element + + """ + params = get_text(elem, 'parameters').split() + return cls(*map(float, params)) + class Maxwell(Univariate): r"""Maxwellian distribution in energy. diff --git a/src/distribution.cpp b/src/distribution.cpp index 62296cb25d..8fedbf57a6 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -85,6 +85,32 @@ double Uniform::sample(uint64_t* seed) const return a_ + prn(seed) * (b_ - a_); } +//============================================================================== +// Rational implementation +//============================================================================== + +Rational::Rational(pugi::xml_node node) +{ + auto params = get_node_array(node, "parameters"); + if (params.size() != 3) { + openmc::fatal_error("Rational distribution must have three " + "parameters specified."); + } + + a_ = params.at(0); + b_ = params.at(1); + n_ = params.at(2); + an_ = std::pow(a_, n_+1); + bn_ = std::pow(b_, n_+1); +} + +double Rational::sample(uint64_t* seed) const +{ + const double u = prn(seed); + const double r = std::pow(an_+u*(bn_-an_), 1./(n_+1)); + return r; +} + //============================================================================== // Maxwell implementation //============================================================================== @@ -347,6 +373,8 @@ UPtrDist distribution_from_xml(pugi::xml_node node) UPtrDist dist; if (type == "uniform") { dist = UPtrDist {new Uniform(node)}; + } else if (type == "rational") { + dist = UPtrDist{new Rational(node)}; } else if (type == "maxwell") { dist = UPtrDist {new Maxwell(node)}; } else if (type == "watt") {