diff --git a/docs/source/io_formats/settings.rst b/docs/source/io_formats/settings.rst index c797212c1..38fc2d79f 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 "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` + 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..1d4ef0302 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.PowerLaw openmc.stats.Maxwell openmc.stats.Watt openmc.stats.Tabular diff --git a/include/openmc/distribution.h b/include/openmc/distribution.h index 8cdd78f78..8403f2c67 100644 --- a/include/openmc/distribution.h +++ b/include/openmc/distribution.h @@ -80,6 +80,32 @@ private: double b_; //!< Upper bound of distribution }; +//============================================================================== +//! PowerLaw distribution over the interval [a,b] with exponent n : p(x)=c x^n +//============================================================================== + +class PowerLaw : public Distribution { +public: + 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)} {}; + + //! 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 std::pow(offset_, ninv_); } + 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) +}; + //============================================================================== //! Maxwellian distribution of form c*sqrt(E)*exp(-E/theta) //============================================================================== diff --git a/openmc/stats/univariate.py b/openmc/stats/univariate.py index af6c9a6d1..8d10104fd 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 == 'powerlaw': + return PowerLaw.from_xml_element(elem) elif distribution == 'maxwell': return Maxwell.from_xml_element(elem) elif distribution == 'watt': @@ -242,6 +244,104 @@ class Uniform(Univariate): params = get_text(elem, 'parameters').split() return cls(*map(float, params)) +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`. + + 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, which is equivalent to a 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('interval lower bound', a, Real) + self._a = a + + @b.setter + def b(self, b): + cv.check_type('interval upper bound', b, Real) + self._b = b + + @n.setter + def n(self, n): + cv.check_type('power law exponent', n, Real) + self._n = n + + def to_xml_element(self, element_name): + """Return XML representation of the power law distribution + + Parameters + ---------- + element_name : str + XML element name + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing distribution data + + """ + element = ET.Element(element_name) + element.set("type", "powerlaw") + element.set("parameters", f'{self.a} {self.b} {self.n}') + return element + + @classmethod + def from_xml_element(cls, elem): + """Generate power law distribution from an XML element + + Parameters + ---------- + elem : xml.etree.ElementTree.Element + XML element + + Returns + ------- + openmc.stats.PowerLaw + 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..901852e0f 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); @@ -85,6 +85,32 @@ double Uniform::sample(uint64_t* seed) const return a_ + prn(seed) * (b_ - a_); } +//============================================================================== +// PowerLaw implementation +//============================================================================== + +PowerLaw::PowerLaw(pugi::xml_node node) +{ + auto params = get_node_array(node, "parameters"); + if (params.size() != 3) { + fatal_error("PowerLaw distribution must have three " + "parameters specified."); + } + + 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(b, n + 1) - offset_; + ninv_ = 1 / (n + 1); +} + +double PowerLaw::sample(uint64_t* seed) const +{ + return std::pow(offset_ + prn(seed) * span_, ninv_); +} + //============================================================================== // 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 == "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 eb60821c6..0568b6543 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..62eaf6eff 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 +2.865754E-01 6.762423E-03 diff --git a/tests/regression_tests/source/test.py b/tests/regression_tests/source/test.py index 006358ca7..41cd25856 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.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) @@ -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..bbcb12ff1 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_powerlaw(): + a, b, n = 10.0, 20.0, 2.0 + d = openmc.stats.PowerLaw(a, b, n) + elem = d.to_xml_element('distribution') + + d = openmc.stats.PowerLaw.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