Convert multivariate probability distributions

CMakeLists.txt

@@ -385,6 +385,8 @@ add_library(libopenmc SHARED
   src/tallies/trigger_header.F90
   src/cell.cpp
   src/distribution.cpp
+  src/distribution_multi.cpp
+  src/distribution_spatial.cpp
   src/initialize.cpp
   src/finalize.cpp
   src/geometry_aux.cpp

src/distribution.cpp

@@ -1,26 +1,40 @@
 #include "distribution.h"
 
-#include <algorithm>
-#include <cmath>
-#include <numeric>
+#include <algorithm> // for copy
+#include <cmath>     // for sqrt, floor, max
+#include <iterator>  // for back_inserter
+#include <numeric>   // for accumulate
+#include <string>    // for string, stod
 
 #include "error.h"
 #include "math_functions.h"
 #include "random_lcg.h"
+#include "xml_interface.h"
 
 namespace openmc {
 
+//==============================================================================
+// Discrete implementation
+//==============================================================================
+
+Discrete::Discrete(pugi::xml_node node)
+{
+  auto params = get_node_array<double>(node, "parameters");
+
+  std::size_t n = params.size();
+  std::copy(params.begin(), params.begin() + n/2, std::back_inserter(x_));
+  std::copy(params.begin() + n/2, params.end(), std::back_inserter(p_));
+
+  normalize();
+}
+
 Discrete::Discrete(const double* x, const double* p, int n)
   : x_{x, x+n}, p_{p, p+n}
 {
-  // Renormalize density function so that it sums to unity
-  double norm = std::accumulate(p_.begin(), p_.end(), 0.0);
-  for (auto& p_i : p_)
-    p_i /= norm;
+  normalize();
 }
 
-
-double Discrete::sample()
+double Discrete::sample() const
 {
   int n = x_.size();
   if (n > 1) {
@@ -36,28 +50,107 @@ double Discrete::sample()
   }
 }
 
+void Discrete::normalize()
+{
+  // Renormalize density function so that it sums to unity
+  double norm = std::accumulate(p_.begin(), p_.end(), 0.0);
+  for (auto& p_i : p_)
+    p_i /= norm;
+}
 
-double Uniform::sample()
+//==============================================================================
+// Uniform implementation
+//==============================================================================
+
+Uniform::Uniform(pugi::xml_node node)
+{
+  auto params = get_node_array<double>(node, "parameters");
+  if (params.size() != 2)
+    openmc::fatal_error("Uniform distribution must have two "
+                        "parameters specified.");
+
+  a_ = params.at(0);
+  b_ = params.at(1);
+}
+
+double Uniform::sample() const
 {
   return a_ + prn()*(b_ - a_);
 }
 
+//==============================================================================
+// Maxwell implementation
+//==============================================================================
 
-double Maxwell::sample()
+Maxwell::Maxwell(pugi::xml_node node)
+{
+  theta_ = std::stod(get_node_value(node, "parameters"));
+}
+
+double Maxwell::sample() const
 {
   return maxwell_spectrum_c(theta_);
 }
 
+//==============================================================================
+// Watt implementation
+//==============================================================================
 
-double Watt::sample()
+Watt::Watt(pugi::xml_node node)
+{
+  auto params = get_node_array<double>(node, "parameters");
+  if (params.size() != 2)
+    openmc::fatal_error("Watt energy distribution must have two "
+                        "parameters specified.");
+
+  a_ = params.at(0);
+  b_ = params.at(1);
+}
+
+double Watt::sample() const
 {
   return watt_spectrum_c(a_, b_);
 }
 
+//==============================================================================
+// Tabular implementation
+//==============================================================================
 
-Tabular::Tabular(const double* x, const double* p, int n, Interpolation interp)
-  : x_{x, x+n}, p_{p, p+n}, interp_{interp}, c_(n, 0.0)
+Tabular::Tabular(pugi::xml_node node)
 {
+  if (check_for_node(node, "interpolation")) {
+    std::string temp = get_node_value(node, "interpolation");
+    if (temp == "histogram") {
+      interp_ = Interpolation::histogram;
+    } else if (temp == "linear-linear") {
+      interp_ = Interpolation::lin_lin;
+    } else {
+      openmc::fatal_error("Unknown interpolation type for distribution: " + temp);
+    }
+  } else {
+    interp_ = Interpolation::histogram;
+  }
+
+  // Read and initialize tabular distribution
+  auto params = get_node_array<double>(node, "parameters");
+  std::size_t n = params.size() / 2;
+  const double* x = params.data();
+  const double* p = x + n;
+  init(x, p, n);
+}
+
+Tabular::Tabular(const double* x, const double* p, int n, Interpolation interp, const double* c)
+  : interp_{interp}
+{
+  init(x, p, n, c);
+}
+
+void Tabular::init(const double* x, const double* p, std::size_t n, const double* c)
+{
+  // Copy x/p arrays into vectors
+  std::copy(x, x + n, std::back_inserter(x_));
+  std::copy(p, p + n, std::back_inserter(p_));
+
   // Check interpolation parameter
   if (interp_ != Interpolation::histogram &&
       interp_ != Interpolation::lin_lin) {
@@ -66,11 +159,17 @@ Tabular::Tabular(const double* x, const double* p, int n, Interpolation interp)
   }
 
   // Calculate cumulative distribution function
-  for (int i = 1; i < n; ++i) {
-    if (interp_ == Interpolation::histogram) {
-      c_[i] = c_[i-1] + p_[i-1]*(x_[i] - x_[i-1]);
-    } else if (interp_ == Interpolation::lin_lin) {
-      c_[i] = c_[i-1] + 0.5*(p_[i-1] + p_[i]) * (x_[i] - x_[i-1]);
+  if (c) {
+    std::copy(c, c + n, std::back_inserter(c_));
+  } else {
+    c_.resize(n);
+    c_[0] = 0.0;
+    for (int i = 1; i < n; ++i) {
+      if (interp_ == Interpolation::histogram) {
+        c_[i] = c_[i-1] + p_[i-1]*(x_[i] - x_[i-1]);
+      } else if (interp_ == Interpolation::lin_lin) {
+        c_[i] = c_[i-1] + 0.5*(p_[i-1] + p_[i]) * (x_[i] - x_[i-1]);
+      }
     }
   }
 
@@ -81,8 +180,7 @@ Tabular::Tabular(const double* x, const double* p, int n, Interpolation interp)
   }
 }
 
-
-double Tabular::sample()
+double Tabular::sample() const
 {
   // Sample value of CDF
   double c = prn();
@@ -90,7 +188,7 @@ double Tabular::sample()
   // Find first CDF bin which is above the sampled value
   double c_i = c_[0];
   int i;
-  size_t n = c_.size();
+  std::size_t n = c_.size();
   for (i = 0; i < n - 1; ++i) {
     if (c <= c_[i+1]) break;
     c_i = c_[i+1];
@@ -121,10 +219,13 @@ double Tabular::sample()
   }
 }
 
+//==============================================================================
+// Equiprobable implementation
+//==============================================================================
 
-double Equiprobable::sample()
+double Equiprobable::sample() const
 {
-  size_t n = x_.size();
+  std::size_t n = x_.size();
 
   double r = prn();
   int i = std::floor((n - 1)*r);
@@ -134,4 +235,32 @@ double Equiprobable::sample()
   return xl + ((n - 1)*r - i) * (xr - xl);
 }
 
+//==============================================================================
+// Helper function
+//==============================================================================
+
+UPtrDist distribution_from_xml(pugi::xml_node node)
+{
+  if (!check_for_node(node, "type"))
+    openmc::fatal_error("Distribution type must be specified.");
+
+  // Determine type of distribution
+  std::string type = get_node_value(node, "type");
+
+  // Allocate extension of Distribution
+  if (type == "uniform") {
+    return UPtrDist{new Uniform(node)};
+  } else if (type == "maxwell") {
+    return UPtrDist{new Maxwell(node)};
+  } else if (type == "watt") {
+    return UPtrDist{new Watt(node)};
+  } else if (type == "discrete") {
+    return UPtrDist{new Discrete(node)};
+  } else if (type == "tabular") {
+    return UPtrDist{new Tabular(node)};
+  } else {
+    openmc::fatal_error("Invalid distribution type: " + type);
+  }
+}
+
 } // namespace openmc

src/distribution.h

@@ -1,78 +1,136 @@
 #ifndef DISTRIBUTION_H
 #define DISTRIBUTION_H
 
-#include <vector>
+#include <cstddef> // for size_t
+#include <memory> // for unique_ptr
+#include <vector> // for vector
 
 #include "constants.h"
+#include "pugixml.hpp"
 
 namespace openmc {
 
+//==============================================================================
+//! Abstract class representing a univariate probability distribution
+//==============================================================================
+
 class Distribution {
 public:
-  virtual double sample() = 0;
+  virtual double sample() const = 0;
   virtual ~Distribution() = default;
 };
 
+//==============================================================================
+//! A discrete distribution (probability mass function)
+//==============================================================================
 
 class Discrete : public Distribution {
 public:
+  explicit Discrete(pugi::xml_node node);
   Discrete(const double* x, const double* p, int n);
-  double sample();
+
+  //! Sample a value from the distribution
+  double sample() const;
 private:
   std::vector<double> x_;
   std::vector<double> p_;
+
+  void normalize();
 };
 
+//==============================================================================
+//! Uniform distribution over the interval [a,b]
+//==============================================================================
 
 class Uniform : public Distribution {
 public:
+  explicit Uniform(pugi::xml_node node);
   Uniform(double a, double b) : a_{a}, b_{b} {};
-  double sample();
+
+  //! Sample a value from the distribution
+  double sample() const;
 private:
   double a_;
   double b_;
 };
 
+//==============================================================================
+//! Maxwellian distribution of form c*E*exp(-E/a)
+//==============================================================================
 
 class Maxwell : public Distribution {
 public:
+  explicit Maxwell(pugi::xml_node node);
   Maxwell(double theta) : theta_{theta} { };
-  double sample();
+
+  //! Sample a value from the distribution
+  double sample() const;
 private:
   double theta_;
 };
 
+//==============================================================================
+//! Watt fission spectrum with form c*exp(-E/a)*sinh(sqrt(b*E))
+//==============================================================================
 
 class Watt : public Distribution {
 public:
+  explicit Watt(pugi::xml_node node);
   Watt(double a, double b) : a_{a}, b_{b} { };
-  double sample();
+
+  //! Sample a value from the distribution
+  double sample() const;
 private:
   double a_;
   double b_;
 };
 
+//==============================================================================
+//! Histogram or linear-linear interpolated tabular distribution
+//==============================================================================
 
 class Tabular : public Distribution {
 public:
-  Tabular(const double* x, const double* p, int n, Interpolation interp);
-  double sample();
+  explicit Tabular(pugi::xml_node node);
+  Tabular(const double* x, const double* p, int n, Interpolation interp,
+          const double* c=nullptr);
+
+  //! Sample a value from the distribution
+  double sample() const;
 private:
-  std::vector<double> x_;
-  std::vector<double> p_;
-  std::vector<double> c_;
-  Interpolation interp_;
+  std::vector<double> x_; //!< tabulated independent variable
+  std::vector<double> p_; //!< tabulated probability density
+  std::vector<double> c_; //!< cumulative distribution at tabulated values
+  Interpolation interp_;  //!< interpolation rule
+
+  //! Initialize tabulated probability density function
+  //! @param x Array of values for independent variable
+  //! @param p Array of tabulated probabilities
+  //! @param n Number of tabulated values
+  void init(const double* x, const double* p, std::size_t n,
+            const double* c=nullptr);
 };
 
+//==============================================================================
+//!
+//==============================================================================
 
 class Equiprobable : public Distribution {
 public:
+  explicit Equiprobable(pugi::xml_node node);
   Equiprobable(const double* x, int n) : x_{x, x+n} { };
-  double sample();
+
+  //! Sample a value from the distribution
+  double sample() const;
 private:
   std::vector<double> x_;
 };
 
+
+using UPtrDist = std::unique_ptr<Distribution>;
+
+UPtrDist distribution_from_xml(pugi::xml_node node);
+
 } // namespace openmc
 
 #endif // DISTRIBUTION_H

src/distribution_multi.cpp

@@ -0,0 +1,51 @@
+#include "distribution_multi.h"
+
+#include <algorithm> // for move
+#include <cmath>     // for sqrt, sin, cos, max
+
+#include "constants.h"
+#include "math_functions.h"
+#include "random_lcg.h"
+
+namespace openmc {
+
+//==============================================================================
+// PolarAzimuthal implementation
+//==============================================================================
+
+PolarAzimuthal::PolarAzimuthal(Direction u, UPtrDist mu, UPtrDist phi) :
+  UnitSphereDistribution{u}, mu_{std::move(mu)}, phi_{std::move(phi)} { }
+
+Direction PolarAzimuthal::sample() const
+{
+  // Sample cosine of polar angle
+  double mu = mu_->sample();
+  if (mu == 1.0) return u_ref;
+
+  // Sample azimuthal angle
+  double phi = phi_->sample();
+  return rotate_angle(u_ref, mu, &phi);
+}
+
+//==============================================================================
+// Isotropic implementation
+//==============================================================================
+
+Direction Isotropic::sample() const
+{
+  double phi = 2.0*PI*prn();
+  double mu = 2.0*prn() - 1.0;
+  return {mu, std::sqrt(1.0 - mu*mu) * std::cos(phi),
+      std::sqrt(1.0 - mu*mu) * std::sin(phi)};
+}
+
+//==============================================================================
+// Monodirectional implementation
+//==============================================================================
+
+Direction Monodirectional::sample() const
+{
+  return u_ref;
+}
+
+} // namespace openmc

src/distribution_multi.h

@@ -0,0 +1,73 @@
+#ifndef DISTRIBUTION_MULTI_H
+#define DISTRIBUTION_MULTI_H
+
+#include <memory>
+
+#include "distribution.h"
+#include "position.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Probability density function for points on the unit sphere. Extensions of
+//! this type are used to sample angular distributions for starting sources
+//==============================================================================
+
+class UnitSphereDistribution {
+public:
+  UnitSphereDistribution() { };
+  explicit UnitSphereDistribution(Direction u) : u_ref{u} { };
+  virtual ~UnitSphereDistribution() = default;
+
+  //! Sample a direction from the distribution
+  //! \return Direction sampled
+  virtual Direction sample() const = 0;
+
+  Direction u_ref {0.0, 0.0, 1.0};  //!< reference direction
+};
+
+//==============================================================================
+//! Explicit distribution of polar and azimuthal angles
+//==============================================================================
+
+class PolarAzimuthal : public UnitSphereDistribution {
+public:
+  PolarAzimuthal(Direction u, UPtrDist mu, UPtrDist phi);
+
+  //! Sample a direction from the distribution
+  //! \return Direction sampled
+  Direction sample() const;
+private:
+  UPtrDist mu_;  //!< Distribution of polar angle
+  UPtrDist phi_; //!< Distribution of azimuthal angle
+};
+
+//==============================================================================
+//! Uniform distribution on the unit sphere
+//==============================================================================
+
+class Isotropic : public UnitSphereDistribution {
+public:
+  Isotropic() { };
+
+  //! Sample a direction from the distribution
+  //! \return Sampled direction
+  Direction sample() const;
+};
+
+//==============================================================================
+//! Monodirectional distribution
+//==============================================================================
+
+class Monodirectional : public UnitSphereDistribution {
+public:
+  Monodirectional(Direction u) : UnitSphereDistribution{u} { };
+
+  //! Sample a direction from the distribution
+  //! \return Sampled direction
+  Direction sample() const;
+};
+
+} // namespace openmc
+
+#endif // DISTRIBUTION_MULTI_H

src/distribution_spatial.cpp

@@ -0,0 +1,97 @@
+#include "distribution_spatial.h"
+
+#include "error.h"
+#include "random_lcg.h"
+#include "xml_interface.h"
+
+namespace openmc {
+
+//==============================================================================
+// CartesianIndependent implementation
+//==============================================================================
+
+CartesianIndependent::CartesianIndependent(pugi::xml_node node)
+{
+  // Read distribution for x coordinate
+  if (check_for_node(node, "x")) {
+    pugi::xml_node node_dist = node.child("x");
+    x_ = distribution_from_xml(node_dist);
+  } else {
+    // If no distribution was specified, default to a single point at x=0
+    double x[] {0.0};
+    double p[] {1.0};
+    x_ = UPtrDist{new Discrete{x, p, 1}};
+  }
+
+  // Read distribution for y coordinate
+  if (check_for_node(node, "y")) {
+    pugi::xml_node node_dist = node.child("y");
+    y_ = distribution_from_xml(node_dist);
+  } else {
+    // If no distribution was specified, default to a single point at y=0
+    double x[] {0.0};
+    double p[] {1.0};
+    y_ = UPtrDist{new Discrete{x, p, 1}};
+  }
+
+  // Read distribution for z coordinate
+  if (check_for_node(node, "z")) {
+    pugi::xml_node node_dist = node.child("z");
+    z_ = distribution_from_xml(node_dist);
+  } else {
+    // If no distribution was specified, default to a single point at z=0
+    double x[] {0.0};
+    double p[] {1.0};
+    z_ = UPtrDist{new Discrete{x, p, 1}};
+  }
+}
+
+Position CartesianIndependent::sample() const
+{
+  return {x_->sample(), y_->sample(), z_->sample()};
+}
+
+//==============================================================================
+// SpatialBox implementation
+//==============================================================================
+
+SpatialBox::SpatialBox(pugi::xml_node node)
+{
+  // Read lower-right/upper-left coordinates
+  auto params = get_node_array<double>(node, "parameters");
+  if (params.size() != 6)
+    openmc::fatal_error("Box/fission spatial source must have six "
+                        "parameters specified.");
+
+  lower_left_ = Position{params[0], params[1], params[2]};
+  upper_right_ = Position{params[3], params[4], params[5]};
+}
+
+Position SpatialBox::sample() const
+{
+  Position xi {prn(), prn(), prn()};
+  return lower_left_ + xi*(upper_right_ - lower_left_);
+}
+
+//==============================================================================
+// SpatialPoint implementation
+//==============================================================================
+
+SpatialPoint::SpatialPoint(pugi::xml_node node)
+{
+  // Read location of point source
+  auto params = get_node_array<double>(node, "parameters");
+  if (params.size() != 3)
+    openmc::fatal_error("Point spatial source must have three "
+                        "parameters specified.");
+
+  // Set position
+  r_ = Position{params.data()};
+}
+
+Position SpatialPoint::sample() const
+{
+  return r_;
+}
+
+} // namespace openmc

src/distribution_spatial.h

@@ -0,0 +1,69 @@
+#ifndef OPENMC_DISTRIBTUION_SPATIAL_H
+#define OPENMC_DISTRIBUTION_SPATIAL_H
+
+#include "pugixml.hpp"
+
+#include "distribution.h"
+#include "position.h"
+
+namespace openmc {
+
+//==============================================================================
+//! Probability density function for points in Euclidean space
+//==============================================================================
+
+class SpatialDistribution {
+public:
+  virtual Position sample() const = 0;
+  virtual ~SpatialDistribution() = default;
+};
+
+//==============================================================================
+//! Distribution of points specified by independent distributions in x,y,z
+//==============================================================================
+
+class CartesianIndependent : public SpatialDistribution {
+public:
+  explicit CartesianIndependent(pugi::xml_node node);
+
+  //! Sample a position from the distribution
+  Position sample() const;
+private:
+  UPtrDist x_; //!< Distribution of x coordinates
+  UPtrDist y_; //!< Distribution of y coordinates
+  UPtrDist z_; //!< Distribution of z coordinates
+};
+
+//==============================================================================
+//! Uniform distribution of points over a box
+//==============================================================================
+
+class SpatialBox : public SpatialDistribution {
+public:
+  explicit SpatialBox(pugi::xml_node node);
+
+  //! Sample a position from the distribution
+  Position sample() const;
+private:
+  Position lower_left_;
+  Position upper_right_;
+  bool only_fissionable {false};
+};
+
+//==============================================================================
+//! Distribution at a single point
+//==============================================================================
+
+class SpatialPoint : public SpatialDistribution {
+public:
+  explicit SpatialPoint(pugi::xml_node node);
+
+  //! Sample a position from the distribution
+  Position sample() const;
+private:
+  Position r_;
+};
+
+} // namespace openmc
+
+#endif // OPENMC_DISTRIBUTION_SPATIAL_H

src/math_functions.cpp

@@ -623,6 +623,14 @@ void rotate_angle_c(double uvw[3], double mu, double* phi) {
 }
 
 
+Direction rotate_angle(Direction u, double mu, double* phi)
+{
+  double uvw[] {u.x, u.y, u.z};
+  rotate_angle_c(uvw, mu, phi);
+  return {uvw[0], uvw[1], uvw[2]};
+}
+
+
 double maxwell_spectrum_c(double T) {
   // Set the random numbers
   double r1 = prn();

src/math_functions.h

@@ -1,13 +1,14 @@
 //! \file math_functions.h
 //! A collection of elementary math functions.
 
-#ifndef MATH_FUNCTIONS_H
-#define MATH_FUNCTIONS_H
+#ifndef OPENMC_MATH_FUNCTIONS_H
+#define OPENMC_MATH_FUNCTIONS_H
 
 #include <cmath>
 #include <cstdlib>
 
 #include "constants.h"
+#include "position.h"
 #include "random_lcg.h"
 
 
@@ -106,6 +107,8 @@ extern "C" void calc_zn_c(int n, double rho, double phi, double zn[]);
 
 extern "C" void rotate_angle_c(double uvw[3], double mu, double* phi);
 
+Direction rotate_angle(Direction u, double mu, double* phi);
+
 //==============================================================================
 //! Samples an energy from the Maxwell fission distribution based on a direct
 //! sampling scheme.
@@ -202,4 +205,4 @@ extern "C" double spline_integrate_c(int n, const double x[], const double y[],
                                      const double z[], double xa, double xb);
 
 } // namespace openmc
-#endif // MATH_FUNCTIONS_H
+#endif // OPENMC_MATH_FUNCTIONS_H

src/xml_interface.cpp

@@ -9,7 +9,7 @@
 namespace openmc {
 
 std::string
-get_node_value(pugi::xml_node node, const char *name)
+get_node_value_str(pugi::xml_node node, const char* name)
 {
   // Search for either an attribute or child tag and get the data as a char*.
   const pugi::char_t *value_char;
@@ -23,9 +23,15 @@ get_node_value(pugi::xml_node node, const char *name)
             << node.name() << "\" XML node";
     fatal_error(err_msg);
   }
+  return value_char;
+}
 
-  // Convert to lowercase string.
-  std::string value(value_char);
+
+std::string
+get_node_value(pugi::xml_node node, const char *name)
+{
+  // Get char* and convert to lowercase string.
+  std::string value {get_node_value_str(node, name)};
   std::transform(value.begin(), value.end(), value.begin(), ::tolower);
 
   // Remove whitespace.

src/xml_interface.h

@@ -1,6 +1,7 @@
 #ifndef XML_INTERFACE_H
 #define XML_INTERFACE_H
 
+#include <sstream> // for stringstream
 #include <string>
 #include <vector>
 
@@ -15,7 +16,24 @@ check_for_node(pugi::xml_node node, const char *name)
   return node.attribute(name) || node.child(name);
 }
 
+std::string get_node_value_str(pugi::xml_node node, const char *name);
 std::string get_node_value(pugi::xml_node node, const char *name);
 
+template <typename T>
+std::vector<T> get_node_array(pugi::xml_node node, const char* name)
+{
+  // Get value of node attribute/child
+  std::string s {get_node_value_str(node, name)};
+
+  // Read values one by one into vector
+  std::stringstream iss {s};
+  T value;
+  std::vector<T> values;
+  while (iss >> value)
+    values.push_back(value);
+
+  return values;
+}
+
 } // namespace openmc
 #endif // XML_INTERFACE_H