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

Merge pull request #2214 from paulromano/muir-function

Browse source 
Replace Muir classes with muir function
This commit is contained in:
Jonathan Shimwell 2022-09-22 16:46:26 +01:00 • committed by GitHub
commit cbd8376238
No known key found for this signature in database
GPG key ID: 4AEE18F83AFDEB23
7 changed files with 49 additions and 192 deletions
Download patch file Download diff file Expand all files Collapse all files

2
docs/source/pythonapi/stats.rst
View file

@ -22,7 +22,7 @@ Univariate Probability Distributions
openmc.stats.Legendre
openmc.stats.Mixture
openmc.stats.Normal
openmc.stats.Muir
openmc.stats.muir
Angular Distributions
---------------------

36
include/openmc/distribution.h
View file

@ -31,7 +31,6 @@ using UPtrDist = unique_ptr<Distribution>;
//! \return Unique pointer to distribution
UPtrDist distribution_from_xml(pugi::xml_node node);
//==============================================================================
//! A discrete distribution (probability mass function)
//==============================================================================
@ -99,11 +98,12 @@ public:
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)
double span_; //!< b^(n+1) - a^(n+1)
double ninv_; //!< 1/(n+1)
};
//==============================================================================
@ -172,35 +172,6 @@ private:
double std_dev_; //!< standard deviation [eV]
};
//==============================================================================
//! Muir (fusion) spectrum derived from Normal with extra params e0 is mean
//! std dev is sqrt(4*e0*kt/m)
//==============================================================================
class Muir : public Distribution {
public:
explicit Muir(pugi::xml_node node);
Muir(double e0, double m_rat, double kt)
: e0_ {e0}, m_rat_ {m_rat}, kt_ {kt} {};
//! Sample a value from the distribution
//! \param seed Pseudorandom number seed pointer
//! \return Sampled value
double sample(uint64_t* seed) const;
double e0() const { return e0_; }
double m_rat() const { return m_rat_; }
double kt() const { return kt_; }
private:
// example DT fusion m_rat = 5 (D = 2 + T = 3)
// ion temp = 20000 eV
// mean neutron energy 14.08e6 eV
double e0_; //!< mean neutron energy [eV]
double m_rat_; //!< ratio of reactant masses relative to atomic mass unit
double kt_; //!< ion temperature [eV]
};
//==============================================================================
//! Histogram or linear-linear interpolated tabular distribution
//==============================================================================
@ -277,7 +248,6 @@ private:
distribution_; //!< sub-distributions + cummulative probabilities
};
} // namespace openmc
#endif // OPENMC_DISTRIBUTION_H

17
include/openmc/random_dist.h
View file

@ -64,23 +64,6 @@ extern "C" double watt_spectrum(double a, double b, uint64_t* seed);
extern "C" double normal_variate(double mean, double std_dev, uint64_t* seed);
//==============================================================================
//! Samples an energy from the Muir (Gaussian) energy-dependent distribution.
//!
//! This is another form of the Gaussian distribution but with more easily
//! modifiable parameters
//! https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS
//!
//! \param e0 peak neutron energy [eV]
//! \param m_rat ratio of the fusion reactants to AMU
//! \param kt the ion temperature of the reactants [eV]
//! \param seed A pointer to the pseudorandom seed
//! \result The sampled outgoing energy
//==============================================================================
extern "C" double muir_spectrum(
double e0, double m_rat, double kt, uint64_t* seed);
} // namespace openmc
#endif // OPENMC_RANDOM_DIST_H

134
openmc/stats/univariate.py
View file

@ -2,7 +2,9 @@ from abc import ABC, abstractmethod
from collections import defaultdict
from collections.abc import Iterable
from copy import deepcopy
import math
from numbers import Real
from warnings import warn
from xml.etree import ElementTree as ET
import numpy as np
@ -53,7 +55,9 @@ class Univariate(EqualityMixin, ABC):
elif distribution == 'normal':
return Normal.from_xml_element(elem)
elif distribution == 'muir':
return Muir.from_xml_element(elem)
# Support older files where Muir had its own class
params = [float(x) for x in get_text(elem, 'parameters').split()]
return muir(*params)
elif distribution == 'tabular':
return Tabular.from_xml_element(elem)
elif distribution == 'legendre':
@ -717,119 +721,45 @@ class Normal(Univariate):
return cls(*map(float, params))
class Muir(Univariate):
"""Muir energy spectrum.
def muir(e0, m_rat, kt):
"""Generate a Muir energy spectrum
The Muir energy spectrum is a Gaussian spectrum, but for
convenience reasons allows the user 3 parameters to define
the distribution, e0 the mean energy of particles, the mass
of reactants m_rat, and the ion temperature kt.
The Muir energy spectrum is a normal distribution, but for convenience
reasons allows the user to specify three parameters to define the
distribution: the mean energy of particles ``e0``, the mass of reactants
``m_rat``, and the ion temperature ``kt``.
.. versionadded:: 0.14.0
Parameters
----------
e0 : float
Mean of the Muir distribution in units of eV
Mean of the Muir distribution in [eV]
m_rat : float
Ratio of the sum of the masses of the reaction inputs to an
AMU
Ratio of the sum of the masses of the reaction inputs to 1 amu
kt : float
Ion temperature for the Muir distribution in units of eV
Ion temperature for the Muir distribution in [eV]
Attributes
----------
e0 : float
Mean of the Muir distribution in units of eV
m_rat : float
Ratio of the sum of the masses of the reaction inputs to an
AMU
kt : float
Ion temperature for the Muir distribution in units of eV
Returns
-------
openmc.stats.Normal
Corresponding normal distribution
"""
# https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS
std_dev = math.sqrt(4 * e0 * kt / m_rat)
return Normal(e0, std_dev)
def __init__(self, e0=14.08e6, m_rat = 5., kt = 20000.):
self.e0 = e0
self.m_rat = m_rat
self.kt = kt
def __len__(self):
return 3
@property
def e0(self):
return self._e0
@property
def m_rat(self):
return self._m_rat
@property
def kt(self):
return self._kt
@e0.setter
def e0(self, e0):
cv.check_type('Muir e0', e0, Real)
cv.check_greater_than('Muir e0', e0, 0.0)
self._e0 = e0
@m_rat.setter
def m_rat(self, m_rat):
cv.check_type('Muir m_rat', m_rat, Real)
cv.check_greater_than('Muir m_rat', m_rat, 0.0)
self._m_rat = m_rat
@kt.setter
def kt(self, kt):
cv.check_type('Muir kt', kt, Real)
cv.check_greater_than('Muir kt', kt, 0.0)
self._kt = kt
@property
def std_dev(self):
return np.sqrt(4.*self.e0*self.kt/self.m_rat)
def sample(self, n_samples=1, seed=None):
# Based on LANL report LA-05411-MS
np.random.seed(seed)
return np.random.normal(self.e0, self.std_dev, n_samples)
def to_xml_element(self, element_name):
"""Return XML representation of the Watt distribution
Parameters
----------
element_name : str
XML element name
Returns
-------
element : xml.etree.ElementTree.Element
XML element containing Watt distribution data
"""
element = ET.Element(element_name)
element.set("type", "muir")
element.set("parameters", '{} {} {}'.format(self._e0, self._m_rat, self._kt))
return element
@classmethod
def from_xml_element(cls, elem):
"""Generate Muir distribution from an XML element
Parameters
----------
elem : xml.etree.ElementTree.Element
XML element
Returns
-------
openmc.stats.Muir
Muir distribution generated from XML element
"""
params = get_text(elem, 'parameters').split()
return cls(*map(float, params))
# Retain deprecated name for the time being
def Muir(*args, **kwargs):
# warn of name change
warn(
"The Muir(...) class has been replaced by the muir(...) function and "
"will be removed in a future version of OpenMC. Use muir(...) instead.",
FutureWarning
)
return muir(*args, **kwargs)
class Tabular(Univariate):

33
src/distribution.cpp
View file

@ -165,27 +165,6 @@ double Normal::sample(uint64_t* seed) const
return normal_variate(mean_value_, std_dev_, seed);
}
//==============================================================================
// Muir implementation
//==============================================================================
Muir::Muir(pugi::xml_node node)
{
auto params = get_node_array<double>(node, "parameters");
if (params.size() != 3) {
openmc::fatal_error("Muir energy distribution must have three "
"parameters specified.");
}
e0_ = params.at(0);
m_rat_ = params.at(1);
kt_ = params.at(2);
}
double Muir::sample(uint64_t* seed) const
{
return muir_spectrum(e0_, m_rat_, kt_, seed);
}
//==============================================================================
// Tabular implementation
//==============================================================================
@ -324,8 +303,10 @@ Mixture::Mixture(pugi::xml_node node)
double cumsum = 0.0;
for (pugi::xml_node pair : node.children("pair")) {
// Check that required data exists
if (!pair.attribute("probability")) fatal_error("Mixture pair element does not have probability.");
if (!pair.child("dist")) fatal_error("Mixture pair element does not have a distribution.");
if (!pair.attribute("probability"))
fatal_error("Mixture pair element does not have probability.");
if (!pair.child("dist"))
fatal_error("Mixture pair element does not have a distribution.");
// cummulative sum of probybilities
cumsum += std::stod(pair.attribute("probability").value());
@ -381,14 +362,16 @@ UPtrDist distribution_from_xml(pugi::xml_node node)
dist = UPtrDist {new Watt(node)};
} else if (type == "normal") {
dist = UPtrDist {new Normal(node)};
} else if (type == "muir") {
dist = UPtrDist {new Muir(node)};
} else if (type == "discrete") {
dist = UPtrDist {new Discrete(node)};
} else if (type == "tabular") {
dist = UPtrDist {new Tabular(node)};
} else if (type == "mixture") {
dist = UPtrDist {new Mixture(node)};
} else if (type == "muir") {
openmc::fatal_error(
"'muir' distributions are now specified using the openmc.stats.muir() "
"function in Python. Please regenerate your XML files.");
} else {
openmc::fatal_error("Invalid distribution type: " + type);
}

7
src/random_dist.cpp
View file

@ -46,11 +46,4 @@ double normal_variate(double mean, double standard_deviation, uint64_t* seed)
return mean + standard_deviation * z * x;
}
double muir_spectrum(double e0, double m_rat, double kt, uint64_t* seed)
{
// https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS
double sigma = std::sqrt(4. * e0 * kt / m_rat);
return normal_variate(e0, sigma, seed);
}
} // namespace openmc

12
tests/unit_tests/test_stats.py
View file

@ -376,16 +376,14 @@ def test_muir():
mean = 10.0
mass = 5.0
temp = 20000.
d = openmc.stats.Muir(mean,mass,temp)
d = openmc.stats.muir(mean, mass, temp)
assert isinstance(d, openmc.stats.Normal)
elem = d.to_xml_element('energy')
assert elem.attrib['type'] == 'muir'
assert elem.attrib['type'] == 'normal'
d = openmc.stats.Muir.from_xml_element(elem)
assert d.e0 == pytest.approx(mean)
assert d.m_rat == pytest.approx(mass)
assert d.kt == pytest.approx(temp)
assert len(d) == 3
d = openmc.stats.Univariate.from_xml_element(elem)
assert isinstance(d, openmc.stats.Normal)
# sample muir distribution
n_samples = 10000