25472760865

Committed 07 May 2026 02:32AM UTC coverage: 80.952% (-0.4%) from 81.374%

Build # 25472760865

Build Type

Pull #3757

github

Committed by

web-flow

Commit Message

Merge 4928ac803 into e542b2f03

Pull Request Pull Request #3757: Testing point detectors

Coverage Stats

17726 of 25786 branches covered (68.74%)

Branch coverage included in aggregate %.

44 of 386 new or added lines in 25 files covered. (11.4%)

81 existing lines in 3 files now uncovered.

58569 of 68461 relevant lines covered (85.55%)

47138839.61 hits per line

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

69.73

/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 "openmc/tensor.h"

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

namespace openmc {

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

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

  // 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, bool is_com, double awr) const
{
  // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7-1)
  E_out = E_in;

  double jac = 1.0;
  if (is_com)
    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);

  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 jac * 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
//==============================================================================

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, bool is_com, double awr) const
{
  E_out = E_in;

  double jac = 1.0;
  if (is_com)
    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);

  // 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 jac * 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, bool is_com, double awr) 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;

  double jac = 1.0;
  if (is_com)
    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);

  return jac * get_pdf_discrete_interpolated(mu_out_.slice(i, tensor::all),
                 mu_out_.slice(i + 1, tensor::all), f, 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, bool is_com, double awr) 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);

  double jac = 1.0;
  if (is_com)
    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);

  return jac * get_pdf_discrete_interpolated(mu_out_.slice(i, j, tensor::all),
                 mu_out_.slice(i + 1, j, tensor::all), f, 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 = tensor::Tensor<double>({d.n_e_out, n_mu});
        }

        // Copy outgoing angles
        tensor::View<double> mu_j = d.mu.slice(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, bool is_com, double awr) const
{
  double f;
  int l, j;
  sample_params(E_in, E_out, f, l, j, seed);

  const auto& mu_l = distribution_[l].mu;

  double jac = 1.0;
  if (is_com)
    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);

  return jac * get_pdf_discrete_interpolated(mu_l.slice(j, tensor::all),
                 mu_l.slice(j + 1, tensor::all), f, 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, bool is_com, double awr) const
{
  return sample_dist(E_in, seed)
    .sample_energy_and_pdf(E_in, mu, E_out, seed, is_com, awr);
}

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

openmc-dev / openmc / 25472760865

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