diff --git a/CMakeLists.txt b/CMakeLists.txt index 10484f161e..24b0bef569 100644 --- a/CMakeLists.txt +++ b/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 diff --git a/src/distribution.cpp b/src/distribution.cpp index 323a03de50..01255f2d75 100644 --- a/src/distribution.cpp +++ b/src/distribution.cpp @@ -1,26 +1,40 @@ #include "distribution.h" -#include -#include -#include +#include // for copy +#include // for sqrt, floor, max +#include // for back_inserter +#include // for accumulate +#include // 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(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(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(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(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 diff --git a/src/distribution.h b/src/distribution.h index 675acf9761..705db53dab 100644 --- a/src/distribution.h +++ b/src/distribution.h @@ -1,78 +1,136 @@ #ifndef DISTRIBUTION_H #define DISTRIBUTION_H -#include +#include // for size_t +#include // for unique_ptr +#include // 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 x_; std::vector 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 x_; - std::vector p_; - std::vector c_; - Interpolation interp_; + std::vector x_; //!< tabulated independent variable + std::vector p_; //!< tabulated probability density + std::vector 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 x_; }; + +using UPtrDist = std::unique_ptr; + +UPtrDist distribution_from_xml(pugi::xml_node node); + } // namespace openmc #endif // DISTRIBUTION_H diff --git a/src/distribution_multi.cpp b/src/distribution_multi.cpp new file mode 100644 index 0000000000..36184f7527 --- /dev/null +++ b/src/distribution_multi.cpp @@ -0,0 +1,51 @@ +#include "distribution_multi.h" + +#include // for move +#include // 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 diff --git a/src/distribution_multi.h b/src/distribution_multi.h new file mode 100644 index 0000000000..2d4dab42c1 --- /dev/null +++ b/src/distribution_multi.h @@ -0,0 +1,73 @@ +#ifndef DISTRIBUTION_MULTI_H +#define DISTRIBUTION_MULTI_H + +#include + +#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 diff --git a/src/distribution_spatial.cpp b/src/distribution_spatial.cpp new file mode 100644 index 0000000000..d870a5f93a --- /dev/null +++ b/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(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(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 diff --git a/src/distribution_spatial.h b/src/distribution_spatial.h new file mode 100644 index 0000000000..9b32a6344f --- /dev/null +++ b/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 diff --git a/src/math_functions.cpp b/src/math_functions.cpp index 3342e3a912..ddb7506af7 100644 --- a/src/math_functions.cpp +++ b/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(); diff --git a/src/math_functions.h b/src/math_functions.h index 566e68ea60..78aa166b32 100644 --- a/src/math_functions.h +++ b/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 #include #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 diff --git a/src/xml_interface.cpp b/src/xml_interface.cpp index b4fcc524c6..7f3222a524 100644 --- a/src/xml_interface.cpp +++ b/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. diff --git a/src/xml_interface.h b/src/xml_interface.h index de89018efe..b12e6185b6 100644 --- a/src/xml_interface.h +++ b/src/xml_interface.h @@ -1,6 +1,7 @@ #ifndef XML_INTERFACE_H #define XML_INTERFACE_H +#include // for stringstream #include #include @@ -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 +std::vector 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 values; + while (iss >> value) + values.push_back(value); + + return values; +} + } // namespace openmc #endif // XML_INTERFACE_H