crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-21 06:25:30 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Implement angular PDF evaluation for angle-energy distributions (#3550)

Browse source 
Co-authored-by: Your Name <you@example.com>
Co-authored-by: GuySten <62616591+GuySten@users.noreply.github.com>
Co-authored-by: GuySten <guyste@post.bgu.ac.il>
Co-authored-by: Eliezer214 <110336440+Eliezer214@users.noreply.github.com>
Co-authored-by: Paul Romano <paul.k.romano@gmail.com>
This commit is contained in:
itay-space 2026-03-14 23:58:17 +02:00 • committed by GitHub
parent bc9c31e0f9
commit 1578698129
No known key found for this signature in database
GPG key ID: B5690EEEBB952194
21 changed files with 632 additions and 51 deletions
Download patch file Download diff file Expand all files Collapse all files

14
include/openmc/angle_energy.h
View file

@ -14,8 +14,22 @@ namespace openmc {
class AngleEnergy {
public:
//! Sample an outgoing energy and scattering cosine
//! \param[in] E_in Incoming energy in [eV]
//! \param[out] E_out Outgoing energy in [eV]
//! \param[out] mu Outgoing cosine with respect to current direction
//! \param[inout] seed Pseudorandom seed pointer
virtual void sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const = 0;
//! Sample an outgoing energy and evaluate the angular PDF
//! \param[in] E_in Incoming energy in [eV]
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
virtual double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const = 0;
virtual ~AngleEnergy() = default;
};

9
include/openmc/chain.h
View file

@ -71,6 +71,15 @@ public:
void sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const override;
//! Sample an outgoing energy and evaluate the angular PDF
//! \param[in] E_in Incoming energy in [eV]
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
private:
const Distribution* photon_energy_;
};

6
include/openmc/distribution_angle.h
View file

@ -26,6 +26,12 @@ public:
//! \return Cosine of the angle in the range [-1,1]
double sample(double E, uint64_t* seed) const;
//! Evaluate the angular PDF at a given energy and cosine
//! \param[in] E Particle energy in [eV]
//! \param[in] mu Cosine of the scattering angle
//! \return Probability density for the scattering cosine
double evaluate(double E, double mu) const;
//! Determine whether angle distribution is empty
//! \return Whether distribution is empty
bool empty() const { return energy_.empty(); }

19
include/openmc/distribution_multi.h
View file

@ -6,6 +6,7 @@
#include "pugixml.hpp"
#include "openmc/distribution.h"
#include "openmc/error.h"
#include "openmc/position.h"
namespace openmc {
@ -29,6 +30,14 @@ public:
//! \return (sampled Direction, sample weight)
virtual std::pair<Direction, double> sample(uint64_t* seed) const = 0;
//! Evaluate the probability density for a given direction
//! \param[in] u Direction on the unit sphere
//! \return Probability density at the given direction
virtual double evaluate(Direction u) const
{
fatal_error("evaluate not available for this UnitSphereDistribution type");
}
Direction u_ref_ {0.0, 0.0, 1.0}; //!< reference direction
};
@ -52,6 +61,11 @@ public:
//! \return (sampled Direction, value of the PDF at this Direction)
std::pair<Direction, double> sample_as_bias(uint64_t* seed) const;
//! Evaluate the probability density for a given direction
//! \param[in] u Direction on the unit sphere
//! \return Probability density at the given direction
double evaluate(Direction u) const override;
// Observing pointers
Distribution* mu() const { return mu_.get(); }
Distribution* phi() const { return phi_.get(); }
@ -87,6 +101,11 @@ public:
//! \return (sampled direction, sample weight)
std::pair<Direction, double> sample(uint64_t* seed) const override;
//! Evaluate the probability density for a given direction
//! \param[in] u Direction on the unit sphere
//! \return Probability density at the given direction
double evaluate(Direction u) const override;
// Set or get bias distribution
void set_bias(std::unique_ptr<PolarAzimuthal> bias)
{

15
include/openmc/reaction_product.h
View file

@ -49,6 +49,21 @@ public:
//! \param[inout] seed Pseudorandom seed pointer
void sample(double E_in, double& E_out, double& mu, uint64_t* seed) const;
//! Select which angle-energy distribution to sample
//! \param[in] E_in Incoming energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Reference to the selected angle-energy distribution
AngleEnergy& sample_dist(double E_in, uint64_t* seed) const;
//! Sample an outgoing energy and evaluate the angular PDF
//! \param[in] E_in Incoming energy in [eV]
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const;
ParticleType particle_; //!< Particle type
EmissionMode emission_mode_; //!< Emission mode
double decay_rate_; //!< Decay rate (for delayed neutron precursors) in [1/s]

16
include/openmc/secondary_correlated.h
View file

@ -41,6 +41,22 @@ public:
void sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const override;
//! Sample the outgoing energy and return the angular distribution
//! \param[in] E_in Incoming energy in [eV]
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Reference to the angular distribution at the sampled energy bin
Distribution& sample_dist(double E_in, double& E_out, uint64_t* seed) const;
//! Sample an outgoing energy and evaluate the angular PDF
//! \param[in] E_in Incoming energy in [eV]
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
// energy property
vector<double>& energy() { return energy_; }
const vector<double>& energy() const { return energy_; }

18
include/openmc/secondary_kalbach.h
View file

@ -32,6 +32,24 @@ public:
void sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const override;
//! Sample outgoing energy and Kalbach-Mann parameters
//! \param[in] E_in Incoming energy in [eV]
//! \param[out] E_out Outgoing energy in [eV]
//! \param[out] km_a Kalbach-Mann 'a' parameter
//! \param[out] km_r Kalbach-Mann pre-compound fraction 'r'
//! \param[inout] seed Pseudorandom seed pointer
void sample_params(double E_in, double& E_out, double& km_a, double& km_r,
uint64_t* seed) const;
//! Sample an outgoing energy and evaluate the angular PDF
//! \param[in] E_in Incoming energy in [eV]
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
private:
//! Outgoing energy/angle at a single incoming energy
struct KMTable {

15
include/openmc/secondary_nbody.h
View file

@ -28,6 +28,21 @@ public:
void sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const override;
//! Sample an outgoing energy from the N-body phase space distribution
//! \param[in] E_in Incoming energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Sampled outgoing energy in [eV]
double sample_energy(double E_in, uint64_t* seed) const;
//! Sample an outgoing energy and evaluate the angular PDF
//! \param[in] E_in Incoming energy in [eV]
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
private:
int n_bodies_; //!< Number of particles distributed
double mass_ratio_; //!< Total mass of particles [neutron mass]

207
include/openmc/secondary_thermal.h
View file

@ -6,6 +6,7 @@
#include "openmc/angle_energy.h"
#include "openmc/endf.h"
#include "openmc/search.h"
#include "openmc/secondary_correlated.h"
#include "openmc/vector.h"
@ -33,8 +34,20 @@ public:
void sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const override;
//! Sample an outgoing energy and evaluate the angular PDF
//! \param[in] E_in Incoming energy in [eV]
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
private:
const CoherentElasticXS& xs_; //!< Coherent elastic scattering cross section
tensor::Tensor<double> bragg_edges_; //!< Copy of Bragg edges for slicing
tensor::Tensor<double>
factors_diff_; //!< Differences over elastic scattering factors
};
//==============================================================================
@ -56,6 +69,15 @@ public:
void sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const override;
//! Sample an outgoing energy and evaluate the angular PDF
//! \param[in] E_in Incoming energy in [eV]
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
private:
double debye_waller_;
};
@ -81,6 +103,15 @@ public:
void sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const override;
//! Sample an outgoing energy and evaluate the angular PDF
//! \param[in] E_in Incoming energy in [eV]
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
private:
const vector<double>& energy_; //!< Energies at which cosines are tabulated
tensor::Tensor<double> mu_out_; //!< Cosines for each incident energy
@ -106,6 +137,21 @@ public:
//! \param[inout] seed Pseudorandom number seed pointer
void sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const override;
//! Sample outgoing energy bin parameters
//! \param[in] E_in Incoming energy in [eV]
//! \param[out] E_out Outgoing energy in [eV]
//! \param[out] j Sampled outgoing energy bin index
//! \param[inout] seed Pseudorandom seed pointer
void sample_params(double E_in, double& E_out, int& j, uint64_t* seed) const;
//! Sample an outgoing energy and evaluate the angular PDF
//! \param[in] E_in Incoming energy in [eV]
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
private:
const vector<double>& energy_; //!< Incident energies
@ -135,6 +181,25 @@ public:
void sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const override;
//! Sample outgoing energy bin parameters
//! \param[in] E_in Incoming energy in [eV]
//! \param[out] E_out Outgoing energy in [eV]
//! \param[out] f Interpolation factor within sampled energy bin
//! \param[out] l Index of the closer incident energy
//! \param[out] j Sampled outgoing energy bin index
//! \param[inout] seed Pseudorandom seed pointer
void sample_params(double E_in, double& E_out, double& f, int& l, int& j,
uint64_t* seed) const;
//! Sample an outgoing energy and evaluate the angular PDF
//! \param[in] E_in Incoming energy in [eV]
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
private:
//! Secondary energy/angle distribution
struct DistEnergySab {
@ -170,6 +235,21 @@ public:
void sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const override;
//! Select the coherent or incoherent elastic distribution to sample
//! \param[in] E_in Incoming energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Reference to the selected angle-energy distribution
const AngleEnergy& sample_dist(double E_in, uint64_t* seed) const;
//! Sample an outgoing energy and evaluate the angular PDF
//! \param[in] E_in Incoming energy in [eV]
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
private:
CoherentElasticAE coherent_dist_; //!< Coherent distribution
unique_ptr<AngleEnergy> incoherent_dist_; //!< Incoherent distribution
@ -178,6 +258,133 @@ private:
const Function1D& incoherent_xs_; //!< Polymorphic ref. to incoherent XS
};
//! Internal helper for evaluating a piecewise-constant PDF on discrete points.
//!
//! The underlying discrete points are represented implicitly through a
//! monotonically increasing `center(i)` function and corresponding per-point
//! `weight(i)` values. Each point contributes a rectangular bin whose
//! half-width is half the distance to its nearest neighboring center.
//!
//! \tparam CenterFn Callable returning the location of the i-th discrete value
//! \tparam WeightFn Callable returning the weight of the i-th discrete value
//! \param[in] n Number of discrete values
//! \param[in] mu_0 Point at which to evaluate the PDF
//! \param[in] a Lower bound of the domain (default: -1)
//! \param[in] b Upper bound of the domain (default: 1)
//! \return Probability density at mu_0
template<typename CenterFn, typename WeightFn>
double get_pdf_discrete_impl(std::size_t n, double mu_0, double a, double b,
CenterFn center, WeightFn weight)
{
if (n == 0 || mu_0 < a || mu_0 > b)
return 0.0;
auto evaluate_bin = [&](std::size_t i) {
double x = center(i);
double left_span = (i == 0) ? 2.0 * (x - a) : x - center(i - 1);
double right_span = (i + 1 == n) ? 2.0 * (b - x) : center(i + 1) - x;
double delta = 0.5 * std::min(left_span, right_span);
if (delta <= 0.0)
return 0.0;
double left = x - delta;
double right = x + delta;
bool in_bin =
(mu_0 >= left) && ((i + 1 == n) ? (mu_0 <= right) : (mu_0 < right));
return in_bin ? weight(i) / (2.0 * delta) : 0.0;
};
// This is effectively a lower_bound over the sequence center(i), but the
// sequence is implicit rather than stored in a container, so the STL
// algorithms can not be used.
std::size_t low = 0;
std::size_t high = n;
while (low < high) {
std::size_t mid = low + (high - low) / 2;
if (center(mid) < mu_0) {
low = mid + 1;
} else {
high = mid;
}
}
if (low < n) {
double pdf = evaluate_bin(low);
if (pdf > 0.0)
return pdf;
}
if (low > 0)
return evaluate_bin(low - 1);
return 0.0;
}
//! Evaluate the PDF of a weighted discrete distribution at a given point.
//!
//! Given a set of discrete values mu[i] with weights w[i], this function
//! computes the probability density at mu_0 by treating each discrete value
//! as a rectangular bin. The bin half-width around each discrete value is
//! half the distance to its nearest neighbor.
//!
//! \tparam T1 Container type for discrete cosine values (must support
//! operator[], size())
//! \tparam T2 Container type for weights (must support operator[])
//! \param[in] mu Sorted array of discrete cosine values
//! \param[in] w Weights for each discrete value (need not be normalized)
//! \param[in] mu_0 Point at which to evaluate the PDF
//! \param[in] a Lower bound of the domain (default: -1)
//! \param[in] b Upper bound of the domain (default: 1)
//! \return Probability density at mu_0
template<typename T1, typename T2>
double get_pdf_discrete(
const T1 mu, const T2& w, double mu_0, double a = -1.0, double b = 1.0)
{
// Returns the location of the discrete value for a given index
auto center = [&](std::size_t i) { return mu[i]; };
auto weight = [&](std::size_t i) { return w[i]; };
return get_pdf_discrete_impl(mu.size(), mu_0, a, b, center, weight);
}
//! Evaluate the PDF of a discrete distribution with uniform weights
//!
//! \tparam T1 Container type for discrete cosine values
//! \param[in] mu Sorted array of discrete cosine values
//! \param[in] mu_0 Point at which to evaluate the PDF
//! \param[in] a Lower bound of the domain (default: -1)
//! \param[in] b Upper bound of the domain (default: 1)
//! \return Probability density at mu_0
template<typename T1>
double get_pdf_discrete(
const T1 mu, double mu_0, double a = -1.0, double b = 1.0)
{
auto center = [&](std::size_t i) { return mu[i]; };
auto weight = [&](std::size_t i) { return 1.0 / mu.size(); };
return get_pdf_discrete_impl(mu.size(), mu_0, a, b, center, weight);
}
//! Evaluate the PDF of a uniformly weighted distribution on interpolated points
//!
//! \tparam T1 Container type for the lower tabulated cosine values
//! \tparam T2 Container type for the upper tabulated cosine values
//! \param[in] mu0 Sorted array of discrete cosine values at the lower grid
//! \param[in] mu1 Sorted array of discrete cosine values at the upper grid
//! \param[in] f Interpolation factor between mu0 and mu1
//! \param[in] mu_0 Point at which to evaluate the PDF
//! \param[in] a Lower bound of the domain (default: -1)
//! \param[in] b Upper bound of the domain (default: 1)
//! \return Probability density at mu_0
template<typename T1, typename T2>
double get_pdf_discrete_interpolated(const T1 mu0, const T2 mu1, double f,
double mu_0, double a = -1.0, double b = 1.0)
{
if (mu0.size() != mu1.size())
return 0.0;
// Returns interpolated discrete value for a given index
auto center = [&](std::size_t i) { return mu0[i] + f * (mu1[i] - mu0[i]); };
auto weight = [&](std::size_t i) { return 1.0 / mu0.size(); };
return get_pdf_discrete_impl(mu0.size(), mu_0, a, b, center, weight);
}
} // namespace openmc
#endif // OPENMC_SECONDARY_THERMAL_H

9
include/openmc/secondary_uncorrelated.h
View file

@ -32,6 +32,15 @@ public:
void sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const override;
//! Sample an outgoing energy and evaluate the angular PDF
//! \param[in] E_in Incoming energy in [eV]
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const override;
// Accessors
AngleDistribution& angle() { return angle_; }

20
include/openmc/thermal.h
View file

@ -51,7 +51,25 @@ public:
//! \param[out] mu Outgoing scattering angle cosine
//! \param[inout] seed Pseudorandom seed pointer
void sample(const NuclideMicroXS& micro_xs, double E_in, double* E_out,
double* mu, uint64_t* seed);
double* mu, uint64_t* seed) const;
//! Select the elastic or inelastic distribution to sample
//! \param[in] micro_xs Microscopic cross sections
//! \param[in] E Incident neutron energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Reference to the selected angle-energy distribution
AngleEnergy& sample_dist(
const NuclideMicroXS& micro_xs, double E, uint64_t* seed) const;
//! Sample an outgoing energy and evaluate the angular PDF
//! \param[in] micro_xs Microscopic cross sections
//! \param[in] E_in Incoming energy in [eV]
//! \param[in] mu Scattering cosine with respect to current direction
//! \param[out] E_out Outgoing energy in [eV]
//! \param[inout] seed Pseudorandom seed pointer
//! \return Probability density for the scattering cosine
double sample_energy_and_pdf(const NuclideMicroXS& micro_xs, double E_in,
double mu, double& E_out, uint64_t* seed) const;
private:
struct Reaction {

7
src/chain.cpp
View file

@ -74,6 +74,13 @@ void DecayPhotonAngleEnergy::sample(
mu = Uniform(-1., 1.).sample(seed).first;
}
double DecayPhotonAngleEnergy::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
{
E_out = photon_energy_->sample(seed).first;
return 0.5;
}
//==============================================================================
// Global variables
//==============================================================================

15
src/distribution_angle.cpp
View file

@ -82,4 +82,19 @@ double AngleDistribution::sample(double E, uint64_t* seed) const
return mu;
}
double AngleDistribution::evaluate(double E, double mu) const
{
// Find energy bin and calculate interpolation factor
int i;
double r;
get_energy_index(energy_, E, i, r);
double pdf = 0.0;
if (r > 0.0)
pdf += r * distribution_[i + 1]->evaluate(mu);
if (r < 1.0)
pdf += (1.0 - r) * distribution_[i]->evaluate(mu);
return pdf;
}
} // namespace openmc

25
src/distribution_multi.cpp
View file

@ -1,6 +1,6 @@
#include "openmc/distribution_multi.h"
#include <algorithm> // for move
#include <algorithm> // for move, clamp
#include <cmath> // for sqrt, sin, cos, max
#include "openmc/constants.h"
@ -44,6 +44,7 @@ UnitSphereDistribution::UnitSphereDistribution(pugi::xml_node node)
fatal_error("Angular distribution reference direction must have "
"three parameters specified.");
u_ref_ = Direction(u_ref.data());
u_ref_ /= u_ref_.norm();
}
}
@ -65,6 +66,7 @@ PolarAzimuthal::PolarAzimuthal(pugi::xml_node node)
fatal_error("Angular distribution reference v direction must have "
"three parameters specified.");
v_ref_ = Direction(v_ref.data());
v_ref_ /= v_ref_.norm();
}
w_ref_ = u_ref_.cross(v_ref_);
if (check_for_node(node, "mu")) {
@ -116,6 +118,22 @@ std::pair<Direction, double> PolarAzimuthal::sample_impl(
weight};
}
double PolarAzimuthal::evaluate(Direction u) const
{
double mu = std::clamp(u.dot(u_ref_), -1.0, 1.0);
double phi = 0.0;
double sin_theta_sq = std::max(0.0, 1.0 - mu * mu);
if (sin_theta_sq > 0.0) {
double sin_theta = std::sqrt(sin_theta_sq);
double cos_phi = u.dot(v_ref_) / sin_theta;
double sin_phi = u.dot(w_ref_) / sin_theta;
phi = std::atan2(sin_phi, cos_phi);
if (phi < 0.0)
phi += 2.0 * PI;
}
return mu_->evaluate(mu) * phi_->evaluate(phi);
}
//==============================================================================
// Isotropic implementation
//==============================================================================
@ -157,6 +175,11 @@ std::pair<Direction, double> Isotropic::sample(uint64_t* seed) const
}
}
double Isotropic::evaluate(Direction u) const
{
return 1.0 / (4.0 * PI);
}
//==============================================================================
// Monodirectional implementation
//==============================================================================

29
src/reaction_product.cpp
View file

@ -1,5 +1,6 @@
#include "openmc/reaction_product.h"
#include <cassert>
#include <string> // for string
#include <fmt/core.h>
@ -106,9 +107,10 @@ ReactionProduct::ReactionProduct(const ChainNuclide::Product& product)
make_unique<DecayPhotonAngleEnergy>(chain_nuc->photon_energy()));
}
void ReactionProduct::sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const
AngleEnergy& ReactionProduct::sample_dist(double E_in, uint64_t* seed) const
{
assert(!distribution_.empty());
auto n = applicability_.size();
if (n > 1) {
double prob = 0.0;
@ -118,15 +120,24 @@ void ReactionProduct::sample(
prob += applicability_[i](E_in);
// If i-th distribution is sampled, sample energy from the distribution
if (c <= prob) {
distribution_[i]->sample(E_in, E_out, mu, seed);
break;
}
if (c <= prob)
return *distribution_[i];
}
} else {
// If only one distribution is present, go ahead and sample it
distribution_[0]->sample(E_in, E_out, mu, seed);
}
return *distribution_.back();
}
void ReactionProduct::sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const
{
sample_dist(E_in, seed).sample(E_in, E_out, mu, seed);
}
double ReactionProduct::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
{
return sample_dist(E_in, seed).sample_energy_and_pdf(E_in, mu, E_out, seed);
}
} // namespace openmc

21
src/secondary_correlated.cpp
View file

@ -155,9 +155,8 @@ CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group)
distribution_.push_back(std::move(d));
} // incoming energies
}
void CorrelatedAngleEnergy::sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const
Distribution& CorrelatedAngleEnergy::sample_dist(
double E_in, double& E_out, uint64_t* seed) const
{
// Find energy bin and calculate interpolation factor
int i;
@ -249,10 +248,22 @@ void CorrelatedAngleEnergy::sample(
// Find correlated angular distribution for closest outgoing energy bin
if (r1 - c_k < c_k1 - r1 ||
distribution_[l].interpolation == Interpolation::histogram) {
mu = distribution_[l].angle[k]->sample(seed).first;
return *distribution_[l].angle[k];
} else {
mu = distribution_[l].angle[k + 1]->sample(seed).first;
return *distribution_[l].angle[k + 1];
}
}
void CorrelatedAngleEnergy::sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const
{
mu = sample_dist(E_in, E_out, seed).sample(seed).first;
}
double CorrelatedAngleEnergy::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
{
return sample_dist(E_in, E_out, seed).evaluate(mu);
}
} // namespace openmc

22
src/secondary_kalbach.cpp
View file

@ -114,8 +114,8 @@ KalbachMann::KalbachMann(hid_t group)
} // incoming energies
}
void KalbachMann::sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const
void KalbachMann::sample_params(
double E_in, double& E_out, double& km_a, double& km_r, uint64_t* seed) const
{
// Find energy bin and calculate interpolation factor
int i;
@ -170,7 +170,6 @@ void KalbachMann::sample(
double E_l_k = distribution_[l].e_out[k];
double p_l_k = distribution_[l].p[k];
double km_r, km_a;
if (distribution_[l].interpolation == Interpolation::histogram) {
// Histogram interpolation
if (p_l_k > 0.0 && k >= n_discrete) {
@ -216,6 +215,13 @@ void KalbachMann::sample(
E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);
}
}
}
void KalbachMann::sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const
{
double km_r, km_a;
sample_params(E_in, E_out, km_a, km_r, seed);
// Sampled correlated angle from Kalbach-Mann parameters
if (prn(seed) > km_r) {
@ -226,5 +232,15 @@ void KalbachMann::sample(
mu = std::log(r1 * std::exp(km_a) + (1.0 - r1) * std::exp(-km_a)) / km_a;
}
}
double KalbachMann::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
{
double km_r, km_a;
sample_params(E_in, E_out, km_a, km_r, seed);
// https://docs.openmc.org/en/latest/methods/neutron_physics.html#equation-KM-pdf-angle
return km_a / (2 * std::sinh(km_a)) *
(std::cosh(km_a * mu) + km_r * std::sinh(km_a * mu));
}
} // namespace openmc

28
src/secondary_nbody.cpp
View file

@ -22,13 +22,8 @@ NBodyPhaseSpace::NBodyPhaseSpace(hid_t group)
read_attribute(group, "q_value", Q_);
}
void NBodyPhaseSpace::sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const
double NBodyPhaseSpace::sample_energy(double E_in, uint64_t* seed) const
{
// By definition, the distribution of the angle is isotropic for an N-body
// phase space distribution
mu = uniform_distribution(-1., 1., seed);
// Determine E_max parameter
double Ap = mass_ratio_;
double E_max = (Ap - 1.0) / Ap * (A_ / (A_ + 1.0) * E_in + Q_);
@ -59,12 +54,29 @@ void NBodyPhaseSpace::sample(
std::log(r5) * std::pow(std::cos(PI / 2.0 * r6), 2);
break;
default:
throw std::runtime_error {"N-body phase space with >5 bodies."};
fatal_error("N-body phase space with >5 bodies.");
}
// Now determine v and E_out
double v = x / (x + y);
E_out = E_max * v;
return E_max * v;
}
void NBodyPhaseSpace::sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const
{
// By definition, the distribution of the angle is isotropic for an N-body
// phase space distribution
mu = uniform_distribution(-1., 1., seed);
E_out = sample_energy(E_in, seed);
}
double NBodyPhaseSpace::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
{
E_out = sample_energy(E_in, seed);
return 0.5;
}
} // namespace openmc

150
src/secondary_thermal.cpp
View file

@ -4,6 +4,7 @@
#include "openmc/math_functions.h"
#include "openmc/random_lcg.h"
#include "openmc/search.h"
#include "openmc/vector.h"
#include "openmc/tensor.h"
@ -16,16 +17,26 @@ namespace openmc {
// CoherentElasticAE implementation
//==============================================================================
CoherentElasticAE::CoherentElasticAE(const CoherentElasticXS& xs) : xs_ {xs} {}
CoherentElasticAE::CoherentElasticAE(const CoherentElasticXS& xs) : xs_ {xs}
{
const auto& bragg = xs_.bragg_edges();
auto n = bragg.size();
bragg_edges_ = tensor::Tensor<double>(bragg.data(), n);
const auto& factors = xs_.factors();
factors_diff_ = tensor::zeros<double>({n});
factors_diff_.slice(0) = factors[0];
for (int i = 1; i < n; ++i) {
factors_diff_.slice(i) = factors[i] - factors[i - 1];
}
}
void CoherentElasticAE::sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const
{
// Energy doesn't change in elastic scattering (ENDF-102, Eq. 7-1)
E_out = E_in;
const auto& energies {xs_.bragg_edges()};
assert(E_in >= energies.front());
const int i = lower_bound_index(energies.begin(), energies.end(), E_in);
@ -42,6 +53,25 @@ void CoherentElasticAE::sample(
mu = 1.0 - 2.0 * energies[k] / E_in;
}
double CoherentElasticAE::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
{
// Energy doesn't change in elastic scattering (ENDF-102, Eq. 7-1)
E_out = E_in;
const auto& factors = xs_.factors();
if (E_in < bragg_edges_.front())
return 0.0;
const int i =
lower_bound_index(bragg_edges_.begin(), bragg_edges_.end(), E_in);
double E = 0.5 * (1 - mu) * E_in;
double C = 0.5 * E_in / factors[i];
return C * get_pdf_discrete(bragg_edges_.slice(tensor::range(i + 1)),
factors_diff_.slice(tensor::range(i + 1)), E, 0.0, E_in);
}
//==============================================================================
// IncoherentElasticAE implementation
//==============================================================================
@ -54,12 +84,21 @@ IncoherentElasticAE::IncoherentElasticAE(hid_t group)
void IncoherentElasticAE::sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const
{
E_out = E_in;
// Sample angle by inverting the distribution in ENDF-102, Eq. 7.4
double c = 2 * E_in * debye_waller_;
mu = std::log(1.0 + prn(seed) * (std::exp(2.0 * c) - 1)) / c - 1.0;
// Energy doesn't change in elastic scattering (ENDF-102, Eq. 7.4)
}
double IncoherentElasticAE::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
{
E_out = E_in;
// Sample angle by inverting the distribution in ENDF-102, Eq. 7.4
double c = 2 * E_in * debye_waller_;
double A = c / (1 - std::exp(-2.0 * c)); // normalization factor
return A * std::exp(-c * (1 - mu));
}
//==============================================================================
@ -116,6 +155,20 @@ void IncoherentElasticAEDiscrete::sample(
E_out = E_in;
}
double IncoherentElasticAEDiscrete::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
{
// Get index and interpolation factor for elastic grid
int i;
double f;
get_energy_index(energy_, E_in, i, f);
// Energy doesn't change in elastic scattering
E_out = E_in;
return get_pdf_discrete_interpolated(
mu_out_.slice(i, tensor::all), mu_out_.slice(i + 1, tensor::all), f, mu);
}
//==============================================================================
// IncoherentInelasticAEDiscrete implementation
//==============================================================================
@ -129,8 +182,8 @@ IncoherentInelasticAEDiscrete::IncoherentInelasticAEDiscrete(
read_dataset(group, "skewed", skewed_);
}
void IncoherentInelasticAEDiscrete::sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const
void IncoherentInelasticAEDiscrete::sample_params(
double E_in, double& E_out, int& j, uint64_t* seed) const
{
// Get index and interpolation factor for inelastic grid
int i;
@ -144,7 +197,6 @@ void IncoherentInelasticAEDiscrete::sample(
// for the second and second to last bins, relative to a normal bin
// probability of 1). Otherwise, each bin is equally probable.
int j;
int n = energy_out_.shape(1);
if (!skewed_) {
// All bins equally likely
@ -176,6 +228,18 @@ void IncoherentInelasticAEDiscrete::sample(
// Outgoing energy
E_out = (1 - f) * E_ij + f * E_i1j;
}
void IncoherentInelasticAEDiscrete::sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const
{
// Get index and interpolation factor for inelastic grid
int i;
double f;
get_energy_index(energy_, E_in, i, f);
int j;
sample_params(E_in, E_out, j, seed);
// Sample outgoing cosine bin
int m = mu_out_.shape(2);
@ -189,6 +253,20 @@ void IncoherentInelasticAEDiscrete::sample(
mu = (1 - f) * mu_ijk + f * mu_i1jk;
}
double IncoherentInelasticAEDiscrete::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
{
// Get index and interpolation factor for inelastic grid
int i;
double f;
get_energy_index(energy_, E_in, i, f);
int j;
sample_params(E_in, E_out, j, seed);
return get_pdf_discrete_interpolated(mu_out_.slice(i, j, tensor::all),
mu_out_.slice(i + 1, j, tensor::all), f, mu);
}
//==============================================================================
// IncoherentInelasticAE implementation
//==============================================================================
@ -231,24 +309,23 @@ IncoherentInelasticAE::IncoherentInelasticAE(hid_t group)
}
}
void IncoherentInelasticAE::sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const
void IncoherentInelasticAE::sample_params(
double E_in, double& E_out, double& f, int& l, int& j, uint64_t* seed) const
{
// Get index and interpolation factor for inelastic grid
int i;
double f;
get_energy_index(energy_, E_in, i, f);
double f0;
get_energy_index(energy_, E_in, i, f0);
// Pick closer energy based on interpolation factor
int l = f > 0.5 ? i + 1 : i;
l = f0 > 0.5 ? i + 1 : i;
// Determine outgoing energy bin
// (First reset n_energy_out to the right value)
auto n = distribution_[l].n_e_out;
int n = distribution_[l].n_e_out;
double r1 = prn(seed);
double c_j = distribution_[l].e_out_cdf[0];
double c_j1;
std::size_t j;
for (j = 0; j < n - 1; ++j) {
c_j1 = distribution_[l].e_out_cdf[j + 1];
if (r1 < c_j1)
@ -286,6 +363,15 @@ void IncoherentInelasticAE::sample(
E_out += E_in - E_l;
}
f = (r1 - c_j) / (c_j1 - c_j);
}
void IncoherentInelasticAE::sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const
{
double f;
int l, j;
sample_params(E_in, E_out, f, l, j, seed);
// Sample outgoing cosine bin
int n_mu = distribution_[l].mu.shape(1);
std::size_t k = prn(seed) * n_mu;
@ -294,7 +380,6 @@ void IncoherentInelasticAE::sample(
// a bin of width 0.5*min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
// discrete mu value itself.
const auto& mu_l = distribution_[l].mu;
f = (r1 - c_j) / (c_j1 - c_j);
// Interpolate kth mu value between distributions at energies j and j+1
mu = mu_l(j, k) + f * (mu_l(j + 1, k) - mu_l(j, k));
@ -318,6 +403,19 @@ void IncoherentInelasticAE::sample(
mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5);
}
double IncoherentInelasticAE::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
{
double f;
int l, j;
sample_params(E_in, E_out, f, l, j, seed);
const auto& mu_l = distribution_[l].mu;
return get_pdf_discrete_interpolated(
mu_l.slice(j, tensor::all), mu_l.slice(j + 1, tensor::all), f, mu);
}
//==============================================================================
// MixedElasticAE implementation
//==============================================================================
@ -340,18 +438,30 @@ MixedElasticAE::MixedElasticAE(
close_group(incoherent_group);
}
void MixedElasticAE::sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const
const AngleEnergy& MixedElasticAE::sample_dist(
double E_in, uint64_t* seed) const
{
// Evaluate coherent and incoherent elastic cross sections
double xs_coh = coherent_xs_(E_in);
double xs_incoh = incoherent_xs_(E_in);
if (prn(seed) * (xs_coh + xs_incoh) < xs_coh) {
coherent_dist_.sample(E_in, E_out, mu, seed);
return coherent_dist_;
} else {
incoherent_dist_->sample(E_in, E_out, mu, seed);
return *incoherent_dist_;
}
}
void MixedElasticAE::sample(
double E_in, double& E_out, double& mu, uint64_t* seed) const
{
sample_dist(E_in, seed).sample(E_in, E_out, mu, seed);
}
double MixedElasticAE::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
{
return sample_dist(E_in, seed).sample_energy_and_pdf(E_in, mu, E_out, seed);
}
} // namespace openmc

18
src/secondary_uncorrelated.cpp
View file

@ -65,4 +65,22 @@ void UncorrelatedAngleEnergy::sample(
E_out = energy_->sample(E_in, seed);
}
double UncorrelatedAngleEnergy::sample_energy_and_pdf(
double E_in, double mu, double& E_out, uint64_t* seed) const
{
// Sample outgoing energy
if (energy_ != nullptr) {
E_out = energy_->sample(E_in, seed);
} else {
E_out = E_in;
}
if (!angle_.empty()) {
return angle_.evaluate(E_in, mu);
} else {
// no angle distribution given => assume isotropic for all energies
return 0.5;
}
}
} // namespace openmc

20
src/thermal.cpp
View file

@ -291,16 +291,21 @@ void ThermalData::calculate_xs(
*inelastic = (*inelastic_.xs)(E);
}
void ThermalData::sample(const NuclideMicroXS& micro_xs, double E,
double* E_out, double* mu, uint64_t* seed)
AngleEnergy& ThermalData::sample_dist(
const NuclideMicroXS& micro_xs, double E, uint64_t* seed) const
{
// Determine whether inelastic or elastic scattering will occur
if (prn(seed) < micro_xs.thermal_elastic / micro_xs.thermal) {
elastic_.distribution->sample(E, *E_out, *mu, seed);
return *elastic_.distribution;
} else {
inelastic_.distribution->sample(E, *E_out, *mu, seed);
return *inelastic_.distribution;
}
}
void ThermalData::sample(const NuclideMicroXS& micro_xs, double E,
double* E_out, double* mu, uint64_t* seed) const
{
sample_dist(micro_xs, E, seed).sample(E, *E_out, *mu, seed);
// Because of floating-point roundoff, it may be possible for mu to be
// outside of the range [-1,1). In these cases, we just set mu to exactly
// -1 or 1
@ -308,6 +313,13 @@ void ThermalData::sample(const NuclideMicroXS& micro_xs, double E,
*mu = std::copysign(1.0, *mu);
}
double ThermalData::sample_energy_and_pdf(const NuclideMicroXS& micro_xs,
double E_in, double mu, double& E_out, uint64_t* seed) const
{
return sample_dist(micro_xs, E_in, seed)
.sample_energy_and_pdf(E_in, mu, E_out, seed);
}
void free_memory_thermal()
{
data::thermal_scatt.clear();