From eb0eac79c139049970c439ea79f31badbe548a34 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Tue, 21 Sep 2021 19:05:47 +0000 Subject: [PATCH 1/9] 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 8cdd78f78..6497078e4 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 af6c9a6d1..cdc2f4c30 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 62296cb25..8fedbf57a 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") { From 9eeff1ede98895c4b35a6adfd1f9791dc81dd733 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Tue, 5 Oct 2021 20:04:21 +0000 Subject: [PATCH 2/9] Update to rational distribution * store offset, span and ninv in class * accessor functions calculate a, b, n --- include/openmc/distribution.h | 16 +++++++--------- openmc/stats/univariate.py | 2 +- src/distribution.cpp | 23 +++++++++++------------ 3 files changed, 19 insertions(+), 22 deletions(-) diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index 6497078e4..48c167810 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -87,22 +87,20 @@ private: 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)} {}; + Rational(double a, double b, double n) : 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; - double a() const { return a_; } - double b() const { return b_; } - double n() const { return n_; } + double a() const { return std::pow(offset_, ninv_); } + double b() const { return std::pow(offset_+span_, ninv_); } + double n() const { return 1/ninv_-1; } 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 + double offset_; //!< Lower bound of distribution + double span_; //!< Upper bound of distribution + double ninv_; //!< inverse Exponent of distribution }; //============================================================================== diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index cdc2f4c30..cfbd0535b 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -254,7 +254,7 @@ class Rational(Univariate): b : float, optional Upper bound of the sampling interval. Defaults to unity. n : float, optional - power law exponent. Defaults to zero -> Uniform distribution + power law exponent. Defaults to zero, which is eqivalent to an uniform distribution Attributes ---------- diff --git a/src/distribution.cpp b/src/distribution.cpp index 8fedbf57a..362c237fc 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -72,8 +72,8 @@ Uniform::Uniform(pugi::xml_node node) { auto params = get_node_array(node, "parameters"); if (params.size() != 2) { - openmc::fatal_error("Uniform distribution must have two " - "parameters specified."); + fatal_error("Uniform distribution must have two " + "parameters specified."); } a_ = params.at(0); @@ -93,22 +93,21 @@ 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."); + 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); + const double a = params.at(0); + const double b = params.at(1); + const double n = params.at(2); + offset_ = std::pow(a, n+1); + span_ = std::pow(a, n+1) - offset_; + ninv_ = 1/(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; + return std::pow(offset_ + prn(seed) * span_, ninv_); } //============================================================================== From 557a5c66b2ab0afd601dc71a8e0a50c9032224da Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Tue, 5 Oct 2021 20:06:02 +0000 Subject: [PATCH 3/9] Add rational distribution to documentation --- docs/source/io_formats/settings.rst | 5 +++++ docs/source/pythonapi/stats.rst | 1 + 2 files changed, 6 insertions(+) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index c797212c1..313d064a7 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -638,6 +638,11 @@ variable and whose sub-elements/attributes are as follows: numbers :math:`a` and :math:`b` that define the interval :math:`[a,b]` over which random variates are sampled. + For a "rational" distribution, ``parameters`` should be given as three real + numbers :math:`a` and :math:`b` that define the interval :math:`[a,b]` over + which random variates are sampled and :math:`n` that defines the exponent of + the probability distribution :math:`p(x)=x^n` + For a "discrete" or "tabular" distribution, ``parameters`` provides the :math:`(x,p)` pairs defining the discrete/tabular distribution. All :math:`x` points are given first followed by corresponding :math:`p` points. diff --git a/docs/source/pythonapi/stats.rst b/docs/source/pythonapi/stats.rst index d14aff558..55f18f45d 100644 --- a/docs/source/pythonapi/stats.rst +++ b/docs/source/pythonapi/stats.rst @@ -15,6 +15,7 @@ Univariate Probability Distributions openmc.stats.Univariate openmc.stats.Discrete openmc.stats.Uniform + openmc.stats.Rational openmc.stats.Maxwell openmc.stats.Watt openmc.stats.Tabular From 29a515c7ba7866cf5a6357b8617a667cd8cb9383 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Tue, 5 Oct 2021 20:06:41 +0000 Subject: [PATCH 4/9] Add regression and unit test for rational distribution --- tests/regression_tests/source/inputs_true.dat | 50 ++++++++++++++++++- .../regression_tests/source/results_true.dat | 2 +- tests/regression_tests/source/test.py | 16 ++++-- tests/unit_tests/test_stats.py | 10 ++++ 4 files changed, 72 insertions(+), 6 deletions(-) diff --git a/tests/regression_tests/source/inputs_true.dat b/tests/regression_tests/source/inputs_true.dat index eb60821c6..7dac4a535 100644 --- a/tests/regression_tests/source/inputs_true.dat +++ b/tests/regression_tests/source/inputs_true.dat @@ -16,7 +16,7 @@ 1000 10 5 - + @@ -34,7 +34,7 @@ - + -4.0 -4.0 -4.0 4.0 4.0 4.0 @@ -99,4 +99,50 @@ + + + + + 0.7853981633974483 1.5707963267948966 2.356194490192345 0.3 0.4 0.3 + + + + + + + + + + + + + + 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 1815324f6..08383c09f 100644 --- a/tests/regression_tests/source/results_true.dat +++ b/tests/regression_tests/source/results_true.dat @@ -1,2 +1,2 @@ k-combined: -2.980096E-01 9.632798E-04 +3.007053E-01 3.738114E-03 diff --git a/tests/regression_tests/source/test.py b/tests/regression_tests/source/test.py index 006358ca7..3a194f577 100644 --- a/tests/regression_tests/source/test.py +++ b/tests/regression_tests/source/test.py @@ -30,6 +30,8 @@ class SourceTestHarness(PyAPITestHarness): y_dist = openmc.stats.Discrete([-4., -1., 3.], [0.2, 0.3, 0.5]) z_dist = openmc.stats.Tabular([-2., 0., 2.], [0.2, 0.3, 0.2]) r_dist = openmc.stats.Uniform(2., 3.) + r_dist1 = openmc.stats.Rational(2., 3., 1.) + r_dist2 = openmc.stats.Rational(2., 3., 2.) theta_dist = openmc.stats.Discrete([pi/4, pi/2, 3*pi/4], [0.3, 0.4, 0.3]) phi_dist = openmc.stats.Uniform(0.0, 2*pi) @@ -42,6 +44,12 @@ class SourceTestHarness(PyAPITestHarness): spatial5 = openmc.stats.CylindricalIndependent(r_dist, phi_dist, z_dist, origin=(1., 1., 0.)) + spatial6 = openmc.stats.SphericalIndependent(r_dist2, theta_dist, + phi_dist, + origin=(1., 1., 0.)) + spatial7 = openmc.stats.CylindricalIndependent(r_dist1, phi_dist, + z_dist, + origin=(1., 1., 0.)) mu_dist = openmc.stats.Discrete([-1., 0., 1.], [0.5, 0.25, 0.25]) phi_dist = openmc.stats.Uniform(0., 6.28318530718) @@ -57,18 +65,20 @@ class SourceTestHarness(PyAPITestHarness): 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) + source1 = openmc.Source(spatial1, angle1, energy1, strength=0.3) + source2 = openmc.Source(spatial2, angle2, energy2, strength=0.1) 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) + source7 = openmc.Source(spatial6, angle3, energy4, strength=0.1) + source8 = openmc.Source(spatial7, angle3, energy4, strength=0.1) settings = openmc.Settings() settings.batches = 10 settings.inactive = 5 settings.particles = 1000 - settings.source = [source1, source2, source3, source4, source5, source6] + settings.source = [source1, source2, source3, source4, source5, source6, source7, source8] settings.export_to_xml() diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py index 9e6237ba1..f741a9e5b 100644 --- a/tests/unit_tests/test_stats.py +++ b/tests/unit_tests/test_stats.py @@ -42,6 +42,16 @@ def test_uniform(): assert t.p == [1/(b-a), 1/(b-a)] assert t.interpolation == 'histogram' +def test_rational(): + a, b, n = 10.0, 20.0, 2.0 + d = openmc.stats.Rational(a, b, n) + elem = d.to_xml_element('distribution') + + d = openmc.stats.Rational.from_xml_element(elem) + assert d.a == a + assert d.b == b + assert d.n == n + assert len(d) == 3 def test_maxwell(): theta = 1.2895e6 From b776dfdc0bcff9db3108fd3009f24b580e2f4e9c Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Tue, 5 Oct 2021 20:20:53 +0000 Subject: [PATCH 5/9] Fixed typo in rational distribution constructor --- src/distribution.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/distribution.cpp b/src/distribution.cpp index 362c237fc..54bd67d80 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -100,8 +100,9 @@ Rational::Rational(pugi::xml_node node) const double a = params.at(0); const double b = params.at(1); const double n = params.at(2); + offset_ = std::pow(a, n+1); - span_ = std::pow(a, n+1) - offset_; + span_ = std::pow(b, n + 1) - offset_; ninv_ = 1/(n+1); } From 3028193f585fcc1ef236a7e7ccf9dd9445872c53 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Wed, 6 Oct 2021 07:19:31 +0000 Subject: [PATCH 6/9] Fix test after fixing bug --- tests/regression_tests/source/results_true.dat | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/regression_tests/source/results_true.dat b/tests/regression_tests/source/results_true.dat index 08383c09f..62eaf6eff 100644 --- a/tests/regression_tests/source/results_true.dat +++ b/tests/regression_tests/source/results_true.dat @@ -1,2 +1,2 @@ k-combined: -3.007053E-01 3.738114E-03 +2.865754E-01 6.762423E-03 From 358f64c53ac7dfd1ab11fb3a078bea2b2b2796a4 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Mon, 11 Oct 2021 21:20:30 +0200 Subject: [PATCH 7/9] Apply suggestions from code review Co-authored-by: Paul Romano --- docs/source/io_formats/settings.rst | 2 +- include/openmc/distribution.h | 14 ++++++++------ openmc/stats/univariate.py | 17 ++++++++++------- src/distribution.cpp | 4 ++-- 4 files changed, 21 insertions(+), 16 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 313d064a7..82c45f3ee 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -641,7 +641,7 @@ variable and whose sub-elements/attributes are as follows: For a "rational" distribution, ``parameters`` should be given as three real numbers :math:`a` and :math:`b` that define the interval :math:`[a,b]` over which random variates are sampled and :math:`n` that defines the exponent of - the probability distribution :math:`p(x)=x^n` + the probability distribution :math:`p(x)=c x^n` For a "discrete" or "tabular" distribution, ``parameters`` provides the :math:`(x,p)` pairs defining the discrete/tabular distribution. All :math:`x` diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index 48c167810..5d67f2a8a 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -87,7 +87,9 @@ private: class Rational : public Distribution { public: explicit Rational(pugi::xml_node node); - Rational(double a, double b, double n) : offset_{std::pow(a, n+1)}, span_{std::pow(b, n+1) - offset_}, ninv_{1/(n+1)} {}; + Rational(double a, double b, double n) + : 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 @@ -95,12 +97,12 @@ public: double sample(uint64_t* seed) const; double a() const { return std::pow(offset_, ninv_); } - double b() const { return std::pow(offset_+span_, ninv_); } - double n() const { return 1/ninv_-1; } + double b() const { return std::pow(offset_ + span_, ninv_); } + double n() const { return 1 / ninv_ - 1; } private: - double offset_; //!< Lower bound of distribution - double span_; //!< Upper bound of distribution - double ninv_; //!< inverse Exponent of distribution + double offset_; //!< a^(n+1) + double span_; //!< b^(n+1) - a^(n+1) + double ninv_; //!< 1/(n+1) }; //============================================================================== diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index cfbd0535b..64b94045c 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -246,6 +246,8 @@ class Uniform(Univariate): class Rational(Univariate): """Distribution with power law probability over a finite interval [a,b] + + The power law distribution has density function :math:`p(x) dx = c x^n dx`. Parameters ---------- @@ -254,7 +256,8 @@ class Rational(Univariate): b : float, optional Upper bound of the sampling interval. Defaults to unity. n : float, optional - power law exponent. Defaults to zero, which is eqivalent to an uniform distribution + Power law exponent. Defaults to zero, which is equivalent to a uniform + distribution. Attributes ---------- @@ -289,17 +292,17 @@ class Rational(Univariate): @a.setter def a(self, a): - cv.check_type('Uniform a', a, Real) + cv.check_type('interval lower bound', a, Real) self._a = a @b.setter def b(self, b): - cv.check_type('Uniform b', b, Real) + cv.check_type('interval upper bound', b, Real) self._b = b @n.setter def n(self, n): - cv.check_type('Uniform n', n, Real) + cv.check_type('power law exponent', n, Real) self._n = n def to_xml_element(self, element_name): @@ -313,12 +316,12 @@ class Rational(Univariate): Returns ------- element : xml.etree.ElementTree.Element - XML element containing uniform distribution data + XML element containing distribution data """ element = ET.Element(element_name) element.set("type", "rational") - element.set("parameters", '{} {} {}'.format(self.a, self.b, self.n)) + element.set("parameters", f'{self.a} {self.b} {self.n}') return element @classmethod @@ -333,7 +336,7 @@ class Rational(Univariate): Returns ------- openmc.stats.Rational - Uniform distribution generated from XML element + Distribution generated from XML element """ params = get_text(elem, 'parameters').split() diff --git a/src/distribution.cpp b/src/distribution.cpp index 54bd67d80..5b33d4a7e 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -101,9 +101,9 @@ Rational::Rational(pugi::xml_node node) const double b = params.at(1); const double n = params.at(2); - offset_ = std::pow(a, n+1); + offset_ = std::pow(a, n + 1); span_ = std::pow(b, n + 1) - offset_; - ninv_ = 1/(n+1); + ninv_ = 1 / (n + 1); } double Rational::sample(uint64_t* seed) const From c0bcc3f1cab69cfd088ae5b340566ab399589760 Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Mon, 11 Oct 2021 19:48:15 +0000 Subject: [PATCH 8/9] Renamed Rational to Power Law --- docs/source/io_formats/settings.rst | 2 +- docs/source/pythonapi/stats.rst | 2 +- include/openmc/distribution.h | 11 ++++++----- openmc/stats/univariate.py | 14 +++++++------- src/distribution.cpp | 12 ++++++------ tests/regression_tests/source/inputs_true.dat | 4 ++-- tests/regression_tests/source/test.py | 4 ++-- tests/unit_tests/test_stats.py | 6 +++--- 8 files changed, 28 insertions(+), 27 deletions(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index 82c45f3ee..ac2fda458 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -638,7 +638,7 @@ variable and whose sub-elements/attributes are as follows: numbers :math:`a` and :math:`b` that define the interval :math:`[a,b]` over which random variates are sampled. - For a "rational" distribution, ``parameters`` should be given as three real + For a "power law" distribution, ``parameters`` should be given as three real numbers :math:`a` and :math:`b` that define the interval :math:`[a,b]` over which random variates are sampled and :math:`n` that defines the exponent of the probability distribution :math:`p(x)=c x^n` diff --git a/docs/source/pythonapi/stats.rst b/docs/source/pythonapi/stats.rst index 55f18f45d..1d4ef0302 100644 --- a/docs/source/pythonapi/stats.rst +++ b/docs/source/pythonapi/stats.rst @@ -15,7 +15,7 @@ Univariate Probability Distributions openmc.stats.Univariate openmc.stats.Discrete openmc.stats.Uniform - openmc.stats.Rational + openmc.stats.PowerLaw openmc.stats.Maxwell openmc.stats.Watt openmc.stats.Tabular diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index 5d67f2a8a..8403f2c67 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -81,15 +81,15 @@ private: }; //============================================================================== -//! Rational distribution over the interval [a,b] with exponent n +//! PowerLaw distribution over the interval [a,b] with exponent n : p(x)=c x^n //============================================================================== -class Rational : public Distribution { +class PowerLaw : public Distribution { public: - explicit Rational(pugi::xml_node node); - Rational(double a, double b, double n) + explicit PowerLaw(pugi::xml_node node); + PowerLaw(double a, double b, double n) : offset_ {std::pow(a, n + 1)}, span_ {std::pow(b, n + 1) - offset_}, - ninv_ {1/(n+1)} {}; + ninv_ {1 / (n + 1)} {}; //! Sample a value from the distribution //! \param seed Pseudorandom number seed pointer @@ -100,6 +100,7 @@ public: double b() const { return std::pow(offset_ + span_, ninv_); } double n() const { return 1 / ninv_ - 1; } private: + //! Store processed values in object to allow for faster sampling double offset_; //!< a^(n+1) double span_; //!< b^(n+1) - a^(n+1) double ninv_; //!< 1/(n+1) diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index 64b94045c..8d10104fd 100644 --- a/openmc/stats/univariate.py +++ b/openmc/stats/univariate.py @@ -42,8 +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 == 'powerlaw': + return PowerLaw.from_xml_element(elem) elif distribution == 'maxwell': return Maxwell.from_xml_element(elem) elif distribution == 'watt': @@ -244,7 +244,7 @@ class Uniform(Univariate): params = get_text(elem, 'parameters').split() return cls(*map(float, params)) -class Rational(Univariate): +class PowerLaw(Univariate): """Distribution with power law probability over a finite interval [a,b] The power law distribution has density function :math:`p(x) dx = c x^n dx`. @@ -306,7 +306,7 @@ class Rational(Univariate): self._n = n def to_xml_element(self, element_name): - """Return XML representation of the rational distribution + """Return XML representation of the power law distribution Parameters ---------- @@ -320,13 +320,13 @@ class Rational(Univariate): """ element = ET.Element(element_name) - element.set("type", "rational") + element.set("type", "powerlaw") element.set("parameters", f'{self.a} {self.b} {self.n}') return element @classmethod def from_xml_element(cls, elem): - """Generate rational distribution from an XML element + """Generate power law distribution from an XML element Parameters ---------- @@ -335,7 +335,7 @@ class Rational(Univariate): Returns ------- - openmc.stats.Rational + openmc.stats.PowerLaw Distribution generated from XML element """ diff --git a/src/distribution.cpp b/src/distribution.cpp index 5b33d4a7e..901852e0f 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -86,14 +86,14 @@ double Uniform::sample(uint64_t* seed) const } //============================================================================== -// Rational implementation +// PowerLaw implementation //============================================================================== -Rational::Rational(pugi::xml_node node) +PowerLaw::PowerLaw(pugi::xml_node node) { auto params = get_node_array(node, "parameters"); if (params.size() != 3) { - fatal_error("Rational distribution must have three " + fatal_error("PowerLaw distribution must have three " "parameters specified."); } @@ -106,7 +106,7 @@ Rational::Rational(pugi::xml_node node) ninv_ = 1 / (n + 1); } -double Rational::sample(uint64_t* seed) const +double PowerLaw::sample(uint64_t* seed) const { return std::pow(offset_ + prn(seed) * span_, ninv_); } @@ -373,8 +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 == "powerlaw") { + dist = UPtrDist {new PowerLaw(node)}; } else if (type == "maxwell") { dist = UPtrDist {new Maxwell(node)}; } else if (type == "watt") { diff --git a/tests/regression_tests/source/inputs_true.dat b/tests/regression_tests/source/inputs_true.dat index 7dac4a535..0568b6543 100644 --- a/tests/regression_tests/source/inputs_true.dat +++ b/tests/regression_tests/source/inputs_true.dat @@ -101,7 +101,7 @@ - + 0.7853981633974483 1.5707963267948966 2.356194490192345 0.3 0.4 0.3 @@ -124,7 +124,7 @@ - + -2.0 0.0 2.0 0.2 0.3 0.2 diff --git a/tests/regression_tests/source/test.py b/tests/regression_tests/source/test.py index 3a194f577..41cd25856 100644 --- a/tests/regression_tests/source/test.py +++ b/tests/regression_tests/source/test.py @@ -30,8 +30,8 @@ class SourceTestHarness(PyAPITestHarness): y_dist = openmc.stats.Discrete([-4., -1., 3.], [0.2, 0.3, 0.5]) z_dist = openmc.stats.Tabular([-2., 0., 2.], [0.2, 0.3, 0.2]) r_dist = openmc.stats.Uniform(2., 3.) - r_dist1 = openmc.stats.Rational(2., 3., 1.) - r_dist2 = openmc.stats.Rational(2., 3., 2.) + r_dist1 = openmc.stats.PowerLaw(2., 3., 1.) + r_dist2 = openmc.stats.PowerLaw(2., 3., 2.) theta_dist = openmc.stats.Discrete([pi/4, pi/2, 3*pi/4], [0.3, 0.4, 0.3]) phi_dist = openmc.stats.Uniform(0.0, 2*pi) diff --git a/tests/unit_tests/test_stats.py b/tests/unit_tests/test_stats.py index f741a9e5b..bbcb12ff1 100644 --- a/tests/unit_tests/test_stats.py +++ b/tests/unit_tests/test_stats.py @@ -42,12 +42,12 @@ def test_uniform(): assert t.p == [1/(b-a), 1/(b-a)] assert t.interpolation == 'histogram' -def test_rational(): +def test_powerlaw(): a, b, n = 10.0, 20.0, 2.0 - d = openmc.stats.Rational(a, b, n) + d = openmc.stats.PowerLaw(a, b, n) elem = d.to_xml_element('distribution') - d = openmc.stats.Rational.from_xml_element(elem) + d = openmc.stats.PowerLaw.from_xml_element(elem) assert d.a == a assert d.b == b assert d.n == n From 2fe12062d29bcce7f09582e9131015dcba325a9b Mon Sep 17 00:00:00 2001 From: Olaf Schumann Date: Tue, 12 Oct 2021 08:55:37 +0200 Subject: [PATCH 9/9] Update docs/source/io_formats/settings.rst Co-authored-by: Paul Romano --- docs/source/io_formats/settings.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index ac2fda458..38fc2d79f 100644 --- a/docs/source/io_formats/settings.rst +++ b/docs/source/io_formats/settings.rst @@ -638,7 +638,7 @@ variable and whose sub-elements/attributes are as follows: numbers :math:`a` and :math:`b` that define the interval :math:`[a,b]` over which random variates are sampled. - For a "power law" distribution, ``parameters`` should be given as three real + For a "powerlaw" distribution, ``parameters`` should be given as three real numbers :math:`a` and :math:`b` that define the interval :math:`[a,b]` over which random variates are sampled and :math:`n` that defines the exponent of the probability distribution :math:`p(x)=c x^n`