21489819490

Committed 29 Jan 2026 06:21PM UTC coverage: 80.077% (-1.9%) from 81.953%

Build # 21489819490

Build Type

Pull #3757

github

Committed by

web-flow

Commit Message

Merge d08626053 into f7a734189

Pull Request Pull Request #3757: Testing point detectors

Run Details

16004 of 22621 branches covered (70.75%)

Branch coverage included in aggregate %.

94 of 518 new or added lines in 26 files covered. (18.15%)

1021 existing lines in 52 files now uncovered.

53779 of 64524 relevant lines covered (83.35%)

8016833.26 hits per line

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

60.14

/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, double mu_0)
{
  DoubleVector w {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 n = upper_bound_index(energies.begin(), energies.end(), E_in);

  vector<double> mu_vector;
  mu_vector.reserve(n);

  std::transform(energies.rbegin() + n - 1, energies.rend(),
    std::back_inserter(mu_vector),
    [E_in](double Ei) { return 1 - 2 * Ei / E_in; });

  vector<double> weights;
  weights.reserve(n);

  weights.emplace_back(factors[0] / factors[n]);
  for (int i = 1; i <= n; ++i) {
    weights.emplace_back((factors[i] - factors[i - 1]) / factors[n]);
  }
  return get_pdf_discrete(mu_vector, weights, 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;
  mu_vector.reserve(n_mu);

  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;
  mu_vector.reserve(m);

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

openmc-dev / openmc / 21489819490

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