Merge pull request #1131 from makeclean/muir_gauss

First commit of Muir/Gaussian sampling

docs/source/pythonapi/stats.rst

@@ -20,6 +20,8 @@ Univariate Probability Distributions
    openmc.stats.Tabular
    openmc.stats.Legendre
    openmc.stats.Mixture
+   openmc.stats.Normal
+   openmc.stats.Muir
 
 Angular Distributions
 ---------------------

include/openmc/distribution.h

@@ -98,6 +98,45 @@ private:
   double b_; //!< Factor in square root [1/eV]
 };
 
+//==============================================================================
+//! Normal distributions with form 1/2*std_dev*sqrt(pi) exp (-(e-E0)/2*std_dev)^2
+//==============================================================================
+
+class Normal : public Distribution {
+public:
+  explicit Normal(pugi::xml_node node);
+  Normal(double mean_value, double std_dev) : mean_value_{mean_value}, std_dev_{std_dev} { };
+
+  //! Sample a value from the distribution
+  //! \return Sampled value
+  double sample() const;
+private:
+  double mean_value_;    //!< middle of distribution [eV]
+  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
+  //! \return Sampled value
+  double sample() const;
+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
 //==============================================================================

include/openmc/math_functions.h

@@ -161,6 +161,39 @@ extern "C" double maxwell_spectrum(double T);
 
 extern "C" double watt_spectrum(double a, double b);
 
+//==============================================================================
+//! Samples an energy from the Gaussian energy-dependent fission distribution.
+//!
+//! Samples from a Normal distribution with a given mean and standard deviation
+//! The PDF is defined as s(x) = (1/2*sigma*sqrt(2) * e-((mu-x)/2*sigma)^2
+//! Its sampled according to 
+//! http://www-pdg.lbl.gov/2009/reviews/rpp2009-rev-monte-carlo-techniques.pdf
+//! section 33.4.4
+//!
+//! @param mean mean of the Gaussian distribution
+//! @param std_dev standard deviation of the Gaussian distribution
+//! @result The sampled outgoing energy
+//==============================================================================
+
+extern "C" double normal_variate(double mean, double std_dev);
+
+//==============================================================================
+//! 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]
+//! @result The sampled outgoing energy
+//==============================================================================
+
+extern "C" double muir_spectrum(double e0, double m_rat, double kt);
+
+
+
 //==============================================================================
 //! Doppler broadens the windowed multipole curvefit.
 //!

openmc/capi/math.py

@@ -37,6 +37,8 @@ _dll.broaden_wmp_polynomials_c.restype = None
 _dll.broaden_wmp_polynomials_c.argtypes = [c_double, c_double, c_int,
                                            ndpointer(c_double)]
 
+_dll.normal_variate.restype = c_double
+_dll.normal_variate.argtypes = [c_double, c_double]
 
 def t_percentile(p, df):
     """ Calculate the percentile of the Student's t distribution with a
@@ -253,6 +255,26 @@ def watt_spectrum(a, b):
     return _dll.watt_spectrum(a, b)
 
 
+def normal_variate(mean_value, std_dev):
+    """ Samples an energy from the Normal distribution.
+
+    Parameters
+    ----------
+    mean_value : float
+        Mean of the Normal distribution 
+    std_dev : float
+        Standard deviation of the normal distribution
+
+    Returns
+    -------
+    float
+        Sampled outgoing normally distributed value
+
+    """
+
+    return _dll.normal_variate(mean_value, std_dev)
+
+
 def broaden_wmp_polynomials(E, dopp, n):
     """ Doppler broadens the windowed multipole curvefit.  The curvefit is a
     polynomial of the form a/E + b/sqrt(E) + c + d sqrt(E) ...

openmc/stats/univariate.py

@@ -308,6 +308,163 @@ class Watt(Univariate):
         element.set("parameters", '{} {}'.format(self.a, self.b))
         return element
 
+class Normal(Univariate):
+    r"""Normally distributed sampling.
+
+    The Normal Distribution is characterized by two parameters
+    :math:`\mu` and :math:`\sigma` and has density function 
+    :math:`p(X) dX = 1/(\sqrt{2\pi}\sigma) e^{-(X-\mu)^2/(2\sigma^2)}`
+
+    Parameters
+    ----------
+    mean_value : float
+        Mean value of the  distribution
+    std_dev : float
+        Standard deviation of the Normal distribution
+
+    Attributes
+    ----------
+    mean_value : float
+        Mean of the Normal distribution 
+    std_dev : float
+        Standard deviation of the Normal distribution 
+    """
+
+    def __init__(self, mean_value, std_dev):
+        super().__init__()
+        self.mean_value = mean_value
+        self.std_dev = std_dev
+
+    def __len__(self):
+        return 2
+
+    @property
+    def mean_value(self):
+        return self._mean_value
+
+    @property
+    def std_dev(self):
+        return self._std_dev
+
+    @mean_value.setter
+    def mean_value(self, mean_value):
+        cv.check_type('Normal mean_value', mean_value, Real)
+        cv.check_greater_than('Normal mean_value', mean_value, 0.0)
+        self._mean_value = mean_value
+
+    @std_dev.setter
+    def std_dev(self, std_dev):
+        cv.check_type('Normal std_dev', std_dev, Real)
+        cv.check_greater_than('Normal std_dev', std_dev, 0.0)
+        self._std_dev = std_dev
+
+    def to_xml_element(self, element_name):
+        """Return XML representation of the Normal 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", "normal")
+        element.set("parameters", '{} {}'.format(self.mean_value, self.std_dev))
+        return element
+
+class Muir(Univariate):
+    """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. 
+
+    Parameters
+    ----------
+    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
+
+    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
+
+    """
+
+    def __init__(self, e0=14.08e6, m_rat = 5., kt = 20000.):
+        super().__init__()
+        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
+
+    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
+
 
 class Tabular(Univariate):
     """Piecewise continuous probability distribution.

src/distribution.cpp

@@ -113,6 +113,45 @@ double Watt::sample() const
   return watt_spectrum(a_, b_);
 }
 
+//==============================================================================
+// Normal implementation
+//==============================================================================
+Normal::Normal(pugi::xml_node node) 
+{
+  auto params = get_node_array<double>(node,"parameters");
+  if (params.size() != 2)
+    openmc::fatal_error("Normal energy distribution must have two "
+                        "parameters specified.");
+
+  mean_value_ = params.at(0);
+  std_dev_ = params.at(1);
+}
+
+double Normal::sample() const
+{
+  return normal_variate(mean_value_, std_dev_);
+}
+
+//==============================================================================
+// 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() const
+{
+  return muir_spectrum(e0_, m_rat_, kt_);
+}
+
 //==============================================================================
 // Tabular implementation
 //==============================================================================
@@ -256,6 +295,10 @@ UPtrDist distribution_from_xml(pugi::xml_node node)
     dist = UPtrDist{new Maxwell(node)};
   } else if (type == "watt") {
     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") {

src/math_functions.cpp

@@ -674,6 +674,30 @@ double maxwell_spectrum(double T) {
 }
 
 
+double normal_variate(double mean, double standard_deviation) {
+  // perhaps there should be a limit to the number of resamples
+  while ( true ) {
+    double v1 = 2 * prn() - 1.;
+    double v2 = 2 * prn() - 1.;
+    
+    double r = std::pow(v1, 2) + std::pow(v2, 2);
+    double r2 = std::pow(r, 2);
+    if (r2 < 1) {
+      double z = std::sqrt(-2.0 * std::log(r2)/r2);
+      z *= (prn() <= 0.5) ? v1 : v2;
+      return mean + standard_deviation*z;
+    }
+  }
+}
+
+double muir_spectrum(double e0, double m_rat, double kt) {
+  // note sigma here is a factor of 2 shy of equation
+  // 8 in https://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-05411-MS
+  double sigma = std::sqrt(2.*e0*kt/m_rat);
+  return normal_variate(e0, sigma);
+}  
+
+
 double watt_spectrum(double a, double b) {
   double w = maxwell_spectrum(a);
   double E_out = w + 0.25 * a * a * b + (2. * prn() - 1.) * std::sqrt(a * a * b * w);

tests/unit_tests/test_math.py

@@ -4,6 +4,7 @@ import scipy as sp
 import openmc
 import openmc.capi
 
+import pytest
 
 def test_t_percentile():
     # Permutations include 1 DoF, 2 DoF, and > 2 DoF
@@ -205,6 +206,25 @@ def test_watt_spectrum():
     assert ref_val == test_val
 
 
+def test_normal_dist():
+    settings = openmc.capi.settings
+    settings.seed = 1
+    a = 14.08
+    b = 0.0
+    ref_val = 14.08
+    test_val = openmc.capi.math.normal_variate(a, b)
+
+    assert ref_val == pytest.approx(test_val)
+
+    settings.seed = 1
+    a = 14.08
+    b = 1.0
+    ref_val = 16.436645416691427
+    test_val = openmc.capi.math.normal_variate(a, b)
+
+    assert ref_val == pytest.approx(test_val)
+
+
 def test_broaden_wmp_polynomials():
     # Two branches of the code to worry about, beta > 6 and otherwise
     # beta = sqrtE * dopp

tests/unit_tests/test_stats.py

@@ -177,3 +177,25 @@ def test_point():
     assert elem.tag == 'space'
     assert elem.attrib['type'] == 'point'
     assert elem.find('parameters') is not None
+
+def test_normal():
+    mean = 10.0
+    std_dev = 2.0
+    d = openmc.stats.Normal(mean,std_dev)
+    assert d.mean_value == pytest.approx(mean)
+    assert d.std_dev == pytest.approx(std_dev)
+    assert len(d) == 2
+    elem = d.to_xml_element('distribution')
+    assert elem.attrib['type'] == 'normal'
+
+def test_muir():
+    mean = 10.0
+    mass = 5.0
+    temp = 20000.
+    d = openmc.stats.Muir(mean,mass,temp)
+    assert d.e0 == pytest.approx(mean)
+    assert d.m_rat == pytest.approx(mass)
+    assert d.kt == pytest.approx(temp)
+    assert len(d) == 3
+    elem = d.to_xml_element('energy')
+    assert elem.attrib['type'] == 'muir'