20470828925

Committed 23 Dec 2025 08:24PM UTC coverage: 81.898% (-0.04%) from 81.94%

Build # 20470828925

Build Type

Pull #3550

github

Committed by

web-flow

Commit Message

Merge b067e7f24 into 3f06a42ab

Pull Request Pull Request #3550: [Point Detector] Add distribution get_pdf functionality

Run Details

17078 of 23776 branches covered (71.83%)

Branch coverage included in aggregate %.

48 of 212 new or added lines in 13 files covered. (22.64%)

4 existing lines in 2 files now uncovered.

55208 of 64487 relevant lines covered (85.61%)

43370386.44 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

54.97

/src/secondary_thermal.cpp

#include "openmc/secondary_thermal.h"

#include "openmc/hdf5_interface.h"
#include "openmc/math_functions.h"
#include "openmc/random_lcg.h"
#include "openmc/search.h"
#include "openmc/vector.h"

#include "xtensor/xview.hpp"

#include <cassert>
#include <cmath> // for log, exp

namespace openmc {

double get_pdf_discrete(
  const vector<double>& mu, const vector<double>& w, double mu_0)
{
  // Make sure mu is in range [-1,1]
  if (std::abs(mu_0) > 1.0)
    mu_0 = std::copysign(1.0, mu_0);
  double a0;
  double a1;
  double b0;
  double b1;
  int32_t ai = -1;
  int32_t bi = -1;
  if (mu_0 > mu[0]) {
    ai = lower_bound_index(mu.begin(), mu.end(), mu_0);
    a0 = mu[ai];
    a1 = (ai > 1) ? mu[ai - 1] : -1.0;
  } else {
    a0 = -1.0;
    a1 = -1.0;
  }
  if (mu_0 < mu[mu.size() - 1]) {
    bi = upper_bound_index(mu.begin(), mu.end(), mu_0);
    b0 = mu[bi];
    b1 = (bi < mu.size() - 1) ? mu[bi + 1] : 1.0;
  } else {
    b0 = 1.0;
    b1 = 1.0;
  }

  //  Calculate Delta_a and Delta_b
  double delta_a = 0.5 * std::min(b0 - a0, a0 - a1);
  double delta_b = 0.5 * std::min(b1 - b0, b0 - a0);

  if (mu_0 < a0 + delta_a)
    return w[ai] / (2.0 * delta_a);
  else if (mu_0 + delta_b < b0)
    return w[bi] / (2.0 * delta_b);
  else
    return 0.0;
}

double get_pdf_discrete(const vector<double>& mu, double mu_0)
{
  vector<double> w(mu.size(), 1.0 / mu.size());
  return get_pdf_discrete(mu, w, mu_0);
}

//==============================================================================
// CoherentElasticAE implementation
//==============================================================================

CoherentElasticAE::CoherentElasticAE(const CoherentElasticXS& xs) : xs_ {xs} {}

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);

  // Sample a Bragg edge between 1 and i
  // E[0] < E_in < E[i+1] -> can scatter in bragg edges 0..i
  const auto& factors = xs_.factors();
  const double prob = prn(seed) * factors[i];

  const int k = std::lower_bound(factors.begin(), factors.begin() + i, prob) -
                factors.begin();

  // Characteristic scattering cosine for this Bragg edge (ENDF-102, Eq. 7-2)
  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)

  double pdf;
  E_out = E_in;
  const auto& energies {xs_.bragg_edges()};
  const auto& factors = xs_.factors();

  if (E_in < energies.front() || E_in > energies.back()) {
    return 0;
  }

  const int i = lower_bound_index(energies.begin(), energies.end(), E_in);
  vector<double> energies_cut(energies.begin(), energies.begin() + i + 1);
  vector<double> factors_cut(factors.begin(), factors.begin() + i + 1);

  vector<double> mu_vector_rev;
  std::transform(energies_cut.begin(), energies_cut.end(),
    std::back_inserter(mu_vector_rev),
    [E_in](double Ei) { return 1 - 2 * Ei / E_in; });
  vector<double> mu_vector(mu_vector_rev.rbegin(), mu_vector_rev.rend());

  auto f = xt::adapt(factors_cut, {
                                    factors_cut.size(),
                                  });
  auto weights = xt::diff(f);
  weights /= xt::sum(weights);
  vector<double> w(weights.begin(), weights.end());
  return get_pdf_discrete(mu_vector, w, mu);
}

//==============================================================================
// IncoherentElasticAE implementation
//==============================================================================

IncoherentElasticAE::IncoherentElasticAE(hid_t group)
{
  read_dataset(group, "debye_waller", debye_waller_);
}

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;
}
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));
}

//==============================================================================
// IncoherentElasticAEDiscrete implementation
//==============================================================================

IncoherentElasticAEDiscrete::IncoherentElasticAEDiscrete(
  hid_t group, const vector<double>& energy)
  : energy_ {energy}
{
  read_dataset(group, "mu_out", mu_out_);
}

void IncoherentElasticAEDiscrete::sample(
  double E_in, double& E_out, double& mu, uint64_t* seed) const
{
  // Get index and interpolation factor for elastic grid
  int i;
  double f;
  get_energy_index(energy_, E_in, i, f);

  // Interpolate between two discrete cosines corresponding to neighboring
  // incoming energies.

  // Sample outgoing cosine bin
  int n_mu = mu_out_.shape()[1];
  int k = prn(seed) * n_mu;

  // Rather than use the sampled discrete mu directly, it is smeared over
  // a bin of width 0.5*min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
  // discrete mu value itself.

  // Interpolate kth mu value between distributions at energies i and i+1
  mu = mu_out_(i, k) + f * (mu_out_(i + 1, k) - mu_out_(i, k));

  // Inteprolate (k-1)th mu value between distributions at energies i and i+1.
  // When k==0, pick a value that will smear the cosine out to a minimum of -1.
  double mu_left = (k == 0) ? -1.0 - (mu + 1.0)
                            : mu_out_(i, k - 1) +
                                f * (mu_out_(i + 1, k - 1) - mu_out_(i, k - 1));

  // Inteprolate (k+1)th mu value between distributions at energies i and i+1.
  // When k is the last discrete value, pick a value that will smear the cosine
  // out to a maximum of 1.
  double mu_right =
    (k == n_mu - 1)
      ? 1.0 + (1.0 - mu)
      : mu_out_(i, k + 1) + f * (mu_out_(i + 1, k + 1) - mu_out_(i, k + 1));

  // Smear cosine
  mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5);

  // Energy doesn't change in elastic scattering
  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;
  int n_mu = mu_out_.shape()[1];

  std::vector<double> mu_vector;

  for (int k = 0; k < n_mu; ++k) {
    double mu_k = mu_out_(i, k) + f * (mu_out_(i + 1, k) - mu_out_(i, k));
    mu_vector.push_back(mu_k);
  }

  return get_pdf_discrete(mu_vector, mu);
}

//==============================================================================
// IncoherentInelasticAEDiscrete implementation
//==============================================================================

IncoherentInelasticAEDiscrete::IncoherentInelasticAEDiscrete(
  hid_t group, const vector<double>& energy)
  : energy_ {energy}
{
  read_dataset(group, "energy_out", energy_out_);
  read_dataset(group, "mu_out", mu_out_);
  read_dataset(group, "skewed", skewed_);
}

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;
  double f;
  get_energy_index(energy_, E_in, i, f);

  // Now that we have an incoming energy bin, we need to determine the outgoing
  // energy bin. This will depend on whether the outgoing energy distribution is
  // skewed. If it is skewed, then the first two and last two bins have lower
  // probabilities than the other bins (0.1 for the first and last bins and 0.4
  // for the second and second to last bins, relative to a normal bin
  // probability of 1). Otherwise, each bin is equally probable.

  int n = energy_out_.shape()[1];
  if (!skewed_) {
    // All bins equally likely
    j = prn(seed) * n;
  } else {
    // Distribution skewed away from edge points
    double r = prn(seed) * (n - 3);
    if (r > 1.0) {
      // equally likely N-4 middle bins
      j = r + 1;
    } else if (r > 0.6) {
      // second to last bin has relative probability of 0.4
      j = n - 2;
    } else if (r > 0.5) {
      // last bin has relative probability of 0.1
      j = n - 1;
    } else if (r > 0.1) {
      // second bin has relative probability of 0.4
      j = 1;
    } else {
      // first bin has relative probability of 0.1
      j = 0;
    }
  }

  // Determine outgoing energy corresponding to E_in[i] and E_in[i+1]
  double E_ij = energy_out_(i, j);
  double E_i1j = energy_out_(i + 1, j);

  // 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];
  int k = prn(seed) * m;

  // Determine outgoing cosine corresponding to E_in[i] and E_in[i+1]
  double mu_ijk = mu_out_(i, j, k);
  double mu_i1jk = mu_out_(i + 1, j, k);

  // Cosine of angle between incoming and outgoing neutron
  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);

  int m = mu_out_.shape()[2];
  std::vector<double> mu_vector;

  for (int k = 0; k < m; ++k) {
    double mu_ijk = mu_out_(i, j, k);
    double mu_i1jk = mu_out_(i + 1, j, k);
    double mu_k = (1 - f) * mu_ijk + f * mu_i1jk;
    mu_vector.push_back(mu_k);
  }

  return get_pdf_discrete(mu_vector, mu);
}

//==============================================================================
// IncoherentInelasticAE implementation
//==============================================================================

IncoherentInelasticAE::IncoherentInelasticAE(hid_t group)
{
  // Read correlated angle-energy distribution
  CorrelatedAngleEnergy dist {group};

  // Copy incident energies
  energy_ = dist.energy();

  // Convert to S(a,b) native format
  for (const auto& edist : dist.distribution()) {
    // Create temporary distribution
    DistEnergySab d;

    // Copy outgoing energy distribution
    d.n_e_out = edist.e_out.size();
    d.e_out = edist.e_out;
    d.e_out_pdf = edist.p;
    d.e_out_cdf = edist.c;

    for (int j = 0; j < d.n_e_out; ++j) {
      auto adist = dynamic_cast<Tabular*>(edist.angle[j].get());
      if (adist) {
        // On first pass, allocate space for angles
        if (j == 0) {
          auto n_mu = adist->x().size();
          d.mu = xt::empty<double>({d.n_e_out, n_mu});
        }

        // Copy outgoing angles
        auto mu_j = xt::view(d.mu, j);
        std::copy(adist->x().begin(), adist->x().end(), mu_j.begin());
      }
    }

    distribution_.emplace_back(std::move(d));
  }
}

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 f0;
  get_energy_index(energy_, E_in, i, f0);

  // Pick closer energy based on interpolation factor
  l = f0 > 0.5 ? i + 1 : i;

  // Determine outgoing energy bin
  // (First reset n_energy_out to the right value)
  int n = distribution_[l].n_e_out;
  double r1 = prn(seed);
  double c_j = distribution_[l].e_out_cdf[0];
  double c_j1;
  for (j = 0; j < n - 1; ++j) {
    c_j1 = distribution_[l].e_out_cdf[j + 1];
    if (r1 < c_j1)
      break;
    c_j = c_j1;
  }

  // check to make sure j is <= n_energy_out - 2
  j = std::min(j, n - 2);

  // Get the data to interpolate between
  double E_l_j = distribution_[l].e_out[j];
  double p_l_j = distribution_[l].e_out_pdf[j];

  // Next part assumes linear-linear interpolation in standard
  double E_l_j1 = distribution_[l].e_out[j + 1];
  double p_l_j1 = distribution_[l].e_out_pdf[j + 1];

  // Find secondary energy (variable E)
  double frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j);
  if (frac == 0.0) {
    E_out = E_l_j + (r1 - c_j) / p_l_j;
  } else {
    E_out = E_l_j +
            (std::sqrt(std::max(0.0, p_l_j * p_l_j + 2.0 * frac * (r1 - c_j))) -
              p_l_j) /
              frac;
  }

  // Adjustment of outgoing energy
  double E_l = energy_[l];
  if (E_out < 0.5 * E_l) {
    E_out *= 2.0 * E_in / E_l - 1.0;
  } else {
    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;

  // Rather than use the sampled discrete mu directly, it is smeared over
  // 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;

  // 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));

  // Inteprolate (k-1)th mu value between distributions at energies j and j+1.
  // When k==0, pick a value that will smear the cosine out to a minimum of -1.
  double mu_left =
    (k == 0)
      ? mu_left = -1.0 - (mu + 1.0)
      : mu_left = mu_l(j, k - 1) + f * (mu_l(j + 1, k - 1) - mu_l(j, k - 1));

  // Inteprolate (k+1)th mu value between distributions at energies j and j+1.
  // When k is the last discrete value, pick a value that will smear the cosine
  // out to a maximum of 1.
  double mu_right =
    (k == n_mu - 1)
      ? mu_right = 1.0 + (1.0 - mu)
      : mu_right = mu_l(j, k + 1) + f * (mu_l(j + 1, k + 1) - mu_l(j, k + 1));

  // Smear cosine
  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);

  int n_mu = distribution_[l].mu.shape()[1];
  const auto& mu_l = distribution_[l].mu;
  std::vector<double> mu_vector;

  for (int k = 0; k < n_mu; ++k) {
    double mu_k = mu_l(j, k) + f * (mu_l(j + 1, k) - mu_l(j, k));
    mu_vector.push_back(mu_k);
  }

  return get_pdf_discrete(mu_vector, mu);
}

//==============================================================================
// MixedElasticAE implementation
//==============================================================================

MixedElasticAE::MixedElasticAE(
  hid_t group, const CoherentElasticXS& coh_xs, const Function1D& incoh_xs)
  : coherent_dist_(coh_xs), coherent_xs_(coh_xs), incoherent_xs_(incoh_xs)
{
  // Read incoherent elastic distribution
  hid_t incoherent_group = open_group(group, "incoherent");
  std::string temp;
  read_attribute(incoherent_group, "type", temp);
  if (temp == "incoherent_elastic") {
    incoherent_dist_ = make_unique<IncoherentElasticAE>(incoherent_group);
  } else if (temp == "incoherent_elastic_discrete") {
    auto xs = dynamic_cast<const Tabulated1D*>(&incoh_xs);
    incoherent_dist_ =
      make_unique<IncoherentElasticAEDiscrete>(incoherent_group, xs->x());
  }
  close_group(incoherent_group);
}

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) {
    return coherent_dist_;
  } else {
    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

1	#include "openmc/secondary_thermal.h"
2
3	#include "openmc/hdf5_interface.h"
4	#include "openmc/math_functions.h"
5	#include "openmc/random_lcg.h"
6	#include "openmc/search.h"
7	#include "openmc/vector.h"
8
9	#include "xtensor/xview.hpp"
10
11	#include <cassert>
12	#include <cmath> // for log, exp
13
14	namespace openmc {
15
NEW 16	double get_pdf_discrete(	×
17	const vector<double>& mu, const vector<double>& w, double mu_0)
18	{
19	// Make sure mu is in range [-1,1]
NEW 20	if (std::abs(mu_0) > 1.0)	×
NEW 21	mu_0 = std::copysign(1.0, mu_0);	×
22	double a0;
23	double a1;
24	double b0;
25	double b1;
NEW 26	int32_t ai = -1;	×
NEW 27	int32_t bi = -1;	×
NEW 28	if (mu_0 > mu[0]) {	×
NEW 29	ai = lower_bound_index(mu.begin(), mu.end(), mu_0);	×
NEW 30	a0 = mu[ai];	×
NEW 31	a1 = (ai > 1) ? mu[ai - 1] : -1.0;	×
32	} else {
NEW 33	a0 = -1.0;	×
NEW 34	a1 = -1.0;	×
35	}
NEW 36	if (mu_0 < mu[mu.size() - 1]) {	×
NEW 37	bi = upper_bound_index(mu.begin(), mu.end(), mu_0);	×
NEW 38	b0 = mu[bi];	×
NEW 39	b1 = (bi < mu.size() - 1) ? mu[bi + 1] : 1.0;	×
40	} else {
NEW 41	b0 = 1.0;	×
NEW 42	b1 = 1.0;	×
43	}
44
45	// Calculate Delta_a and Delta_b
NEW 46	double delta_a = 0.5 * std::min(b0 - a0, a0 - a1);	×
NEW 47	double delta_b = 0.5 * std::min(b1 - b0, b0 - a0);	×
48
NEW 49	if (mu_0 < a0 + delta_a)	×
NEW 50	return w[ai] / (2.0 * delta_a);	×
NEW 51	else if (mu_0 + delta_b < b0)	×
NEW 52	return w[bi] / (2.0 * delta_b);	×
53	else
NEW 54	return 0.0;	×
55	}
56
NEW 57	double get_pdf_discrete(const vector<double>& mu, double mu_0)	×
58	{
NEW 59	vector<double> w(mu.size(), 1.0 / mu.size());	×
NEW 60	return get_pdf_discrete(mu, w, mu_0);	×
NEW 61	}	×
62
63	//==============================================================================
64	// CoherentElasticAE implementation
65	//==============================================================================
66
67	CoherentElasticAE::CoherentElasticAE(const CoherentElasticXS& xs) : xs_ {xs} {}	140✔
68
69	void CoherentElasticAE::sample(	178,013✔
70	double E_in, double& E_out, double& mu, uint64_t* seed) const
71	{
72	// Energy doesn't change in elastic scattering (ENDF-102, Eq. 7-1)
73	E_out = E_in;	178,013✔
74
75	const auto& energies {xs_.bragg_edges()};	178,013✔
76
77	assert(E_in >= energies.front());	145,647!
78
79	const int i = lower_bound_index(energies.begin(), energies.end(), E_in);	178,013✔
80
81	// Sample a Bragg edge between 1 and i
82	// E[0] < E_in < E[i+1] -> can scatter in bragg edges 0..i
83	const auto& factors = xs_.factors();	178,013✔
84	const double prob = prn(seed) * factors[i];	178,013✔
85
86	const int k = std::lower_bound(factors.begin(), factors.begin() + i, prob) -	178,013✔
87	factors.begin();	178,013✔
88
89	// Characteristic scattering cosine for this Bragg edge (ENDF-102, Eq. 7-2)
90	mu = 1.0 - 2.0 * energies[k] / E_in;	178,013✔
91	}	178,013✔
92
NEW 93	double CoherentElasticAE::sample_energy_and_pdf(	×
94	double E_in, double mu, double& E_out, uint64_t* seed) const
95	{
96	// Energy doesn't change in elastic scattering (ENDF-102, Eq. 7-1)
97
98	double pdf;
NEW 99	E_out = E_in;	×
NEW 100	const auto& energies {xs_.bragg_edges()};	×
NEW 101	const auto& factors = xs_.factors();	×
102
NEW 103	if (E_in < energies.front() \|\| E_in > energies.back()) {	×
NEW 104	return 0;	×
105	}
106
NEW 107	const int i = lower_bound_index(energies.begin(), energies.end(), E_in);	×
NEW 108	vector<double> energies_cut(energies.begin(), energies.begin() + i + 1);	×
NEW 109	vector<double> factors_cut(factors.begin(), factors.begin() + i + 1);	×
110
NEW 111	vector<double> mu_vector_rev;	×
NEW 112	std::transform(energies_cut.begin(), energies_cut.end(),	×
113	std::back_inserter(mu_vector_rev),
NEW 114	[E_in](double Ei) { return 1 - 2 * Ei / E_in; });	×
NEW 115	vector<double> mu_vector(mu_vector_rev.rbegin(), mu_vector_rev.rend());	×
116
117	auto f = xt::adapt(factors_cut, {
NEW 118	factors_cut.size(),	×
NEW 119	});	×
NEW 120	auto weights = xt::diff(f);	×
NEW 121	weights /= xt::sum(weights);	×
NEW 122	vector<double> w(weights.begin(), weights.end());	×
NEW 123	return get_pdf_discrete(mu_vector, w, mu);	×
NEW 124	}	×
125
126	//==============================================================================
127	// IncoherentElasticAE implementation
128	//==============================================================================
129
130	IncoherentElasticAE::IncoherentElasticAE(hid_t group)	×
131	{
132	read_dataset(group, "debye_waller", debye_waller_);	×
133	}	×
134
135	void IncoherentElasticAE::sample(	×
136	double E_in, double& E_out, double& mu, uint64_t* seed) const
137	{
NEW 138	E_out = E_in;	×
139
140	// Sample angle by inverting the distribution in ENDF-102, Eq. 7.4
141	double c = 2 * E_in * debye_waller_;	×
142	mu = std::log(1.0 + prn(seed) * (std::exp(2.0 * c) - 1)) / c - 1.0;	×
NEW 143	}	×
NEW 144	double IncoherentElasticAE::sample_energy_and_pdf(	×
145	double E_in, double mu, double& E_out, uint64_t* seed) const
146	{
147	E_out = E_in;	×
148
149	// Sample angle by inverting the distribution in ENDF-102, Eq. 7.4
NEW 150	double c = 2 * E_in * debye_waller_;	×
NEW 151	double A = c / (1 - std::exp(-2.0 * c)); // normalization factor	×
NEW 152	return A * std::exp(-c * (1 - mu));	×
153	}
154
155	//==============================================================================
156	// IncoherentElasticAEDiscrete implementation
157	//==============================================================================
158
159	IncoherentElasticAEDiscrete::IncoherentElasticAEDiscrete(	44✔
160	hid_t group, const vector<double>& energy)	44✔
161	: energy_ {energy}	44✔
162	{
163	read_dataset(group, "mu_out", mu_out_);	44✔
164	}	44✔
165
166	void IncoherentElasticAEDiscrete::sample(	1,342✔
167	double E_in, double& E_out, double& mu, uint64_t* seed) const
168	{
169	// Get index and interpolation factor for elastic grid
170	int i;
171	double f;
172	get_energy_index(energy_, E_in, i, f);	1,342✔
173
174	// Interpolate between two discrete cosines corresponding to neighboring
175	// incoming energies.
176
177	// Sample outgoing cosine bin
178	int n_mu = mu_out_.shape()[1];	1,342✔
179	int k = prn(seed) * n_mu;	1,342✔
180
181	// Rather than use the sampled discrete mu directly, it is smeared over
182	// a bin of width 0.5*min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
183	// discrete mu value itself.
184
185	// Interpolate kth mu value between distributions at energies i and i+1
186	mu = mu_out_(i, k) + f * (mu_out_(i + 1, k) - mu_out_(i, k));	1,342✔
187
188	// Inteprolate (k-1)th mu value between distributions at energies i and i+1.
189	// When k==0, pick a value that will smear the cosine out to a minimum of -1.
190	double mu_left = (k == 0) ? -1.0 - (mu + 1.0)	1,342✔
191	: mu_out_(i, k - 1) +	990✔
192	f * (mu_out_(i + 1, k - 1) - mu_out_(i, k - 1));	990✔
193
194	// Inteprolate (k+1)th mu value between distributions at energies i and i+1.
195	// When k is the last discrete value, pick a value that will smear the cosine
196	// out to a maximum of 1.
197	double mu_right =
198	(k == n_mu - 1)	1,342✔
199	? 1.0 + (1.0 - mu)	1,342✔
200	: mu_out_(i, k + 1) + f * (mu_out_(i + 1, k + 1) - mu_out_(i, k + 1));	990✔
201
202	// Smear cosine
203	mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5);	1,342✔
204
205	// Energy doesn't change in elastic scattering
206	E_out = E_in;	1,342✔
207	}	1,342✔
208
NEW 209	double IncoherentElasticAEDiscrete::sample_energy_and_pdf(	×
210	double E_in, double mu, double& E_out, uint64_t* seed) const
211	{
212	// Get index and interpolation factor for elastic grid
213	int i;
214	double f;
NEW 215	get_energy_index(energy_, E_in, i, f);	×
216	// Energy doesn't change in elastic scattering
NEW 217	E_out = E_in;	×
NEW 218	int n_mu = mu_out_.shape()[1];	×
219
NEW 220	std::vector<double> mu_vector;	×
221
NEW 222	for (int k = 0; k < n_mu; ++k) {	×
NEW 223	double mu_k = mu_out_(i, k) + f * (mu_out_(i + 1, k) - mu_out_(i, k));	×
NEW 224	mu_vector.push_back(mu_k);	×
225	}
226
NEW 227	return get_pdf_discrete(mu_vector, mu);	×
NEW 228	}	×
229
230	//==============================================================================
231	// IncoherentInelasticAEDiscrete implementation
232	//==============================================================================
233
234	IncoherentInelasticAEDiscrete::IncoherentInelasticAEDiscrete(	1,133✔
235	hid_t group, const vector<double>& energy)	1,133✔
236	: energy_ {energy}	1,133✔
237	{
238	read_dataset(group, "energy_out", energy_out_);	1,133✔
239	read_dataset(group, "mu_out", mu_out_);	1,133✔
240	read_dataset(group, "skewed", skewed_);	1,133✔
241	}	1,133✔
242
243	void IncoherentInelasticAEDiscrete::sample_params(	89,943,685✔
244	double E_in, double& E_out, int& j, uint64_t* seed) const
245	{
246	// Get index and interpolation factor for inelastic grid
247	int i;
248	double f;
249	get_energy_index(energy_, E_in, i, f);	89,943,685✔
250
251	// Now that we have an incoming energy bin, we need to determine the outgoing
252	// energy bin. This will depend on whether the outgoing energy distribution is
253	// skewed. If it is skewed, then the first two and last two bins have lower
254	// probabilities than the other bins (0.1 for the first and last bins and 0.4
255	// for the second and second to last bins, relative to a normal bin
256	// probability of 1). Otherwise, each bin is equally probable.
257
258	int n = energy_out_.shape()[1];	89,943,685✔
259	if (!skewed_) {	89,943,685!
260	// All bins equally likely
261	j = prn(seed) * n;	×
262	} else {
263	// Distribution skewed away from edge points
264	double r = prn(seed) * (n - 3);	89,943,685✔
265	if (r > 1.0) {	89,943,685✔
266	// equally likely N-4 middle bins
267	j = r + 1;	88,464,103✔
268	} else if (r > 0.6) {	1,479,582✔
269	// second to last bin has relative probability of 0.4
270	j = n - 2;	587,603✔
271	} else if (r > 0.5) {	891,979✔
272	// last bin has relative probability of 0.1
273	j = n - 1;	148,507✔
274	} else if (r > 0.1) {	743,472✔
275	// second bin has relative probability of 0.4
276	j = 1;	595,673✔
277	} else {
278	// first bin has relative probability of 0.1
279	j = 0;	147,799✔
280	}
281	}
282
283	// Determine outgoing energy corresponding to E_in[i] and E_in[i+1]
284	double E_ij = energy_out_(i, j);	89,943,685✔
285	double E_i1j = energy_out_(i + 1, j);	89,943,685✔
286
287	// Outgoing energy
288	E_out = (1 - f) * E_ij + f * E_i1j;	89,943,685✔
289	}	89,943,685✔
290
291	void IncoherentInelasticAEDiscrete::sample(	89,943,685✔
292	double E_in, double& E_out, double& mu, uint64_t* seed) const
293	{
294	// Get index and interpolation factor for inelastic grid
295	int i;
296	double f;
297	get_energy_index(energy_, E_in, i, f);	89,943,685✔
298
299	int j;
300	sample_params(E_in, E_out, j, seed);	89,943,685✔
301
302	// Sample outgoing cosine bin
303	int m = mu_out_.shape()[2];	89,943,685✔
304	int k = prn(seed) * m;	89,943,685✔
305
306	// Determine outgoing cosine corresponding to E_in[i] and E_in[i+1]
307	double mu_ijk = mu_out_(i, j, k);	89,943,685✔
308	double mu_i1jk = mu_out_(i + 1, j, k);	89,943,685✔
309
310	// Cosine of angle between incoming and outgoing neutron
311	mu = (1 - f) * mu_ijk + f * mu_i1jk;	89,943,685✔
312	}	89,943,685✔
313
NEW 314	double IncoherentInelasticAEDiscrete::sample_energy_and_pdf(	×
315	double E_in, double mu, double& E_out, uint64_t* seed) const
316	{
317	// Get index and interpolation factor for inelastic grid
318	int i;
319	double f;
NEW 320	get_energy_index(energy_, E_in, i, f);	×
321	int j;
NEW 322	sample_params(E_in, E_out, j, seed);	×
323
NEW 324	int m = mu_out_.shape()[2];	×
NEW 325	std::vector<double> mu_vector;	×
326
NEW 327	for (int k = 0; k < m; ++k) {	×
NEW 328	double mu_ijk = mu_out_(i, j, k);	×
NEW 329	double mu_i1jk = mu_out_(i + 1, j, k);	×
NEW 330	double mu_k = (1 - f) * mu_ijk + f * mu_i1jk;	×
NEW 331	mu_vector.push_back(mu_k);	×
332	}
333
NEW 334	return get_pdf_discrete(mu_vector, mu);	×
NEW 335	}	×
336
337	//==============================================================================
338	// IncoherentInelasticAE implementation
339	//==============================================================================
340
341	IncoherentInelasticAE::IncoherentInelasticAE(hid_t group)	44✔
342	{
343	// Read correlated angle-energy distribution
344	CorrelatedAngleEnergy dist {group};	44✔
345
346	// Copy incident energies
347	energy_ = dist.energy();	44✔
348
349	// Convert to S(a,b) native format
350	for (const auto& edist : dist.distribution()) {	176✔
351	// Create temporary distribution
352	DistEnergySab d;	132✔
353
354	// Copy outgoing energy distribution
355	d.n_e_out = edist.e_out.size();	132✔
356	d.e_out = edist.e_out;	132✔
357	d.e_out_pdf = edist.p;	132✔
358	d.e_out_cdf = edist.c;	132✔
359
360	for (int j = 0; j < d.n_e_out; ++j) {	660✔
361	auto adist = dynamic_cast<Tabular*>(edist.angle[j].get());	528✔
362	if (adist) {	528!
363	// On first pass, allocate space for angles
364	if (j == 0) {	528✔
365	auto n_mu = adist->x().size();	132✔
366	d.mu = xt::empty<double>({d.n_e_out, n_mu});	132✔
367	}
368
369	// Copy outgoing angles
370	auto mu_j = xt::view(d.mu, j);	528✔
371	std::copy(adist->x().begin(), adist->x().end(), mu_j.begin());	528✔
372	}	528✔
373	}
374
375	distribution_.emplace_back(std::move(d));	132✔
376	}	132✔
377	}	44✔
378
379	void IncoherentInelasticAE::sample_params(	3,614,897✔
380	double E_in, double& E_out, double& f, int& l, int& j, uint64_t* seed) const
381	{
382	// Get index and interpolation factor for inelastic grid
383	int i;
384	double f0;
385	get_energy_index(energy_, E_in, i, f0);	3,614,897✔
386
387	// Pick closer energy based on interpolation factor
388	l = f0 > 0.5 ? i + 1 : i;	3,614,897✔
389
390	// Determine outgoing energy bin
391	// (First reset n_energy_out to the right value)
392	int n = distribution_[l].n_e_out;	3,614,897✔
393	double r1 = prn(seed);	3,614,897✔
394	double c_j = distribution_[l].e_out_cdf[0];	3,614,897✔
395	double c_j1;
396	for (j = 0; j < n - 1; ++j) {	8,721,537!
397	c_j1 = distribution_[l].e_out_cdf[j + 1];	8,721,537✔
398	if (r1 < c_j1)	8,721,537✔
399	break;	3,614,897✔
400	c_j = c_j1;	5,106,640✔
401	}
402
403	// check to make sure j is <= n_energy_out - 2
404	j = std::min(j, n - 2);	3,614,897✔
405
406	// Get the data to interpolate between
407	double E_l_j = distribution_[l].e_out[j];	3,614,897✔
408	double p_l_j = distribution_[l].e_out_pdf[j];	3,614,897✔
409
410	// Next part assumes linear-linear interpolation in standard
411	double E_l_j1 = distribution_[l].e_out[j + 1];	3,614,897✔
412	double p_l_j1 = distribution_[l].e_out_pdf[j + 1];	3,614,897✔
413
414	// Find secondary energy (variable E)
415	double frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j);	3,614,897✔
416	if (frac == 0.0) {	3,614,897!
417	E_out = E_l_j + (r1 - c_j) / p_l_j;	3,614,897✔
418	} else {
419	E_out = E_l_j +	×
420	(std::sqrt(std::max(0.0, p_l_j * p_l_j + 2.0 * frac * (r1 - c_j))) -	×
421	p_l_j) /	×
422	frac;
423	}
424
425	// Adjustment of outgoing energy
426	double E_l = energy_[l];	3,614,897✔
427	if (E_out < 0.5 * E_l) {	3,614,897✔
428	E_out = 2.0 E_in / E_l - 1.0;	344,476✔
429	} else {
430	E_out += E_in - E_l;	3,270,421✔
431	}
432
433	f = (r1 - c_j) / (c_j1 - c_j);	3,614,897✔
434	}	3,614,897✔
435	void IncoherentInelasticAE::sample(	3,614,897✔
436	double E_in, double& E_out, double& mu, uint64_t* seed) const
437	{
438	double f;
439	int l, j;
440	sample_params(E_in, E_out, f, l, j, seed);	3,614,897✔
441
442	// Sample outgoing cosine bin
443	int n_mu = distribution_[l].mu.shape()[1];	3,614,897✔
444	std::size_t k = prn(seed) * n_mu;	3,614,897✔
445
446	// Rather than use the sampled discrete mu directly, it is smeared over
447	// a bin of width 0.5*min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
448	// discrete mu value itself.
449	const auto& mu_l = distribution_[l].mu;	3,614,897✔
450
451	// Interpolate kth mu value between distributions at energies j and j+1
452	mu = mu_l(j, k) + f * (mu_l(j + 1, k) - mu_l(j, k));	3,614,897✔
453
454	// Inteprolate (k-1)th mu value between distributions at energies j and j+1.
455	// When k==0, pick a value that will smear the cosine out to a minimum of -1.
456	double mu_left =
457	(k == 0)
458	? mu_left = -1.0 - (mu + 1.0)	3,614,897✔
459	: mu_left = mu_l(j, k - 1) + f * (mu_l(j + 1, k - 1) - mu_l(j, k - 1));	3,159,464✔
460
461	// Inteprolate (k+1)th mu value between distributions at energies j and j+1.
462	// When k is the last discrete value, pick a value that will smear the cosine
463	// out to a maximum of 1.
464	double mu_right =
465	(k == n_mu - 1)	3,614,897✔
466	? mu_right = 1.0 + (1.0 - mu)	3,614,897✔
467	: mu_right = mu_l(j, k + 1) + f * (mu_l(j + 1, k + 1) - mu_l(j, k + 1));	3,160,762✔
468
469	// Smear cosine
470	mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5);	3,614,897✔
471	}	3,614,897✔
472
NEW 473	double IncoherentInelasticAE::sample_energy_and_pdf(	×
474	double E_in, double mu, double& E_out, uint64_t* seed) const
475	{
476	double f;
477	int l, j;
NEW 478	sample_params(E_in, E_out, f, l, j, seed);	×
479
NEW 480	int n_mu = distribution_[l].mu.shape()[1];	×
NEW 481	const auto& mu_l = distribution_[l].mu;	×
NEW 482	std::vector<double> mu_vector;	×
483
NEW 484	for (int k = 0; k < n_mu; ++k) {	×
NEW 485	double mu_k = mu_l(j, k) + f * (mu_l(j + 1, k) - mu_l(j, k));	×
NEW 486	mu_vector.push_back(mu_k);	×
487	}
488
NEW 489	return get_pdf_discrete(mu_vector, mu);	×
NEW 490	}	×
491
492	//==============================================================================
493	// MixedElasticAE implementation
494	//==============================================================================
495
496	MixedElasticAE::MixedElasticAE(	44✔
497	hid_t group, const CoherentElasticXS& coh_xs, const Function1D& incoh_xs)	44✔
498	: coherent_dist_(coh_xs), coherent_xs_(coh_xs), incoherent_xs_(incoh_xs)	44✔
499	{
500	// Read incoherent elastic distribution
501	hid_t incoherent_group = open_group(group, "incoherent");	44✔
502	std::string temp;	44✔
503	read_attribute(incoherent_group, "type", temp);	44✔
504	if (temp == "incoherent_elastic") {	44!
505	incoherent_dist_ = make_unique<IncoherentElasticAE>(incoherent_group);	×
506	} else if (temp == "incoherent_elastic_discrete") {	44!
507	auto xs = dynamic_cast<const Tabulated1D*>(&incoh_xs);	44!
508	incoherent_dist_ =
509	make_unique<IncoherentElasticAEDiscrete>(incoherent_group, xs->x());	44✔
510	}
511	close_group(incoherent_group);	44✔
512	}	44✔
513
514	const AngleEnergy& MixedElasticAE::sample_dist(	46,486✔
515	double E_in, uint64_t* seed) const
516	{
517	// Evaluate coherent and incoherent elastic cross sections
518	double xs_coh = coherent_xs_(E_in);	46,486✔
519	double xs_incoh = incoherent_xs_(E_in);	46,486✔
520
521	if (prn(seed) * (xs_coh + xs_incoh) < xs_coh) {	46,486✔
522	return coherent_dist_;	45,144✔
523	} else {
524	return *incoherent_dist_;	1,342✔
525	}
526	}
527
528	void MixedElasticAE::sample(	46,486✔
529	double E_in, double& E_out, double& mu, uint64_t* seed) const
530	{
531	sample_dist(E_in, seed).sample(E_in, E_out, mu, seed);	46,486✔
532	}	46,486✔
533
NEW 534	double MixedElasticAE::sample_energy_and_pdf(	×
535	double E_in, double mu, double& E_out, uint64_t* seed) const
536	{
NEW 537	return sample_dist(E_in, seed).sample_energy_and_pdf(E_in, mu, E_out, seed);	×
538	}
539
540	} // namespace openmc

openmc-dev / openmc / 20470828925

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous