13461947318

Committed 21 Feb 2025 05:24PM UTC coverage: 85.019% (+0.05%) from 84.969%

Build # 13461947318

Build Type

Pull #3140

github

Committed by

web-flow

Commit Message

Merge ece247ada into 2b788ea6e

Pull Request Pull Request #3140: Add Versioning Support from `version.txt`

Run Details

8 of 8 new or added lines in 2 files covered. (100.0%)

1293 existing lines in 43 files now uncovered.

50627 of 59548 relevant lines covered (85.02%)

35584918.67 hits per line

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

92.9

/src/secondary_thermal.cpp

#include "openmc/secondary_thermal.h"

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

#include "xtensor/xview.hpp"

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

namespace openmc {

// Helper function to get index on incident energy grid
void get_energy_index(
  const vector<double>& energies, double E, int& i, double& f)
{
  // Get index and interpolation factor for elastic grid
  i = 0;
  f = 0.0;
  if (E >= energies.front()) {
    i = lower_bound_index(energies.begin(), energies.end(), E);
    if (i + 1 < energies.size())
      f = (E - energies[i]) / (energies[i + 1] - energies[i]);
  }
}

//==============================================================================
// 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;
}

//==============================================================================
// 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
{
  // Sample angle by inverting the distribution in ENDF-102, Eq. 7.4
  double c = 2 * E_in * debye_waller_;
  mu = std::log(1.0 + prn(seed) * (std::exp(2.0 * c) - 1)) / c - 1.0;

  // Energy doesn't change in elastic scattering (ENDF-102, Eq. 7.4)
  E_out = E_in;
}

//==============================================================================
// 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;
}

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

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

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

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

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

  // Determine outgoing energy bin
  // (First reset n_energy_out to the right value)
  auto n = distribution_[l].n_e_out;
  double r1 = prn(seed);
  double c_j = distribution_[l].e_out_cdf[0];
  double c_j1;
  std::size_t j;
  for (j = 0; j < n - 1; ++j) {
    c_j1 = distribution_[l].e_out_cdf[j + 1];
    if (r1 < c_j1)
      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;
  }

  // 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;
  f = (r1 - c_j) / (c_j1 - c_j);

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

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

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

void MixedElasticAE::sample(
  double E_in, double& E_out, double& mu, uint64_t* seed) const
{
  // Evaluate coherent and incoherent elastic cross sections
  double xs_coh = coherent_xs_(E_in);
  double xs_incoh = incoherent_xs_(E_in);

  if (prn(seed) * (xs_coh + xs_incoh) < xs_coh) {
    coherent_dist_.sample(E_in, E_out, mu, seed);
  } else {
    incoherent_dist_->sample(E_in, E_out, mu, seed);
  }
}

} // namespace openmc

1	#include "openmc/secondary_thermal.h"
2
3	#include "openmc/hdf5_interface.h"
4	#include "openmc/random_lcg.h"
5	#include "openmc/search.h"
6
7	#include "xtensor/xview.hpp"
8
9	#include <cassert>
10	#include <cmath> // for log, exp
11
12	namespace openmc {
13
14	// Helper function to get index on incident energy grid
15	void get_energy_index(	78,377,469✔
16	const vector<double>& energies, double E, int& i, double& f)
17	{
18	// Get index and interpolation factor for elastic grid
19	i = 0;	78,377,469✔
20	f = 0.0;	78,377,469✔
21	if (E >= energies.front()) {	78,377,469✔
22	i = lower_bound_index(energies.begin(), energies.end(), E);	78,377,457✔
23	if (i + 1 < energies.size())	78,377,457✔
24	f = (E - energies[i]) / (energies[i + 1] - energies[i]);	78,375,993✔
25	}
26	}	78,377,469✔
27
28	//==============================================================================
29	// CoherentElasticAE implementation
30	//==============================================================================
31
32	CoherentElasticAE::CoherentElasticAE(const CoherentElasticXS& xs) : xs_ {xs} {}	150✔
33
34	void CoherentElasticAE::sample(	198,492✔
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;	198,492✔
39
40	const auto& energies {xs_.bragg_edges()};	198,492✔
41
42	assert(E_in >= energies.front());	165,410✔
43
44	const int i = lower_bound_index(energies.begin(), energies.end(), E_in);	198,492✔
45
46	// Sample a Bragg edge between 1 and i
47	// E[0] < E_in < E[i+1] -> can scatter in bragg edges 0..i
48	const auto& factors = xs_.factors();	198,492✔
49	const double prob = prn(seed) * factors[i];	198,492✔
50
51	const int k = std::lower_bound(factors.begin(), factors.begin() + i, prob) -	198,492✔
52	factors.begin();	198,492✔
53
54	// Characteristic scattering cosine for this Bragg edge (ENDF-102, Eq. 7-2)
55	mu = 1.0 - 2.0 * energies[k] / E_in;	198,492✔
56	}	198,492✔
57
58	//==============================================================================
59	// IncoherentElasticAE implementation
60	//==============================================================================
61
UNCOV 62	IncoherentElasticAE::IncoherentElasticAE(hid_t group)	×
63	{
UNCOV 64	read_dataset(group, "debye_waller", debye_waller_);	×
65	}
66
UNCOV 67	void IncoherentElasticAE::sample(	×
68	double E_in, double& E_out, double& mu, uint64_t* seed) const
69	{
70	// Sample angle by inverting the distribution in ENDF-102, Eq. 7.4
UNCOV 71	double c = 2 * E_in * debye_waller_;	×
72	mu = std::log(1.0 + prn(seed) * (std::exp(2.0 * c) - 1)) / c - 1.0;	×
73
74	// Energy doesn't change in elastic scattering (ENDF-102, Eq. 7.4)
UNCOV 75	E_out = E_in;	×
76	}
77
78	//==============================================================================
79	// IncoherentElasticAEDiscrete implementation
80	//==============================================================================
81
82	IncoherentElasticAEDiscrete::IncoherentElasticAEDiscrete(	48✔
83	hid_t group, const vector<double>& energy)	48✔
84	: energy_ {energy}	48✔
85	{
86	read_dataset(group, "mu_out", mu_out_);	48✔
87	}	48✔
88
89	void IncoherentElasticAEDiscrete::sample(	1,464✔
90	double E_in, double& E_out, double& mu, uint64_t* seed) const
91	{
92	// Get index and interpolation factor for elastic grid
93	int i;
94	double f;
95	get_energy_index(energy_, E_in, i, f);	1,464✔
96
97	// Interpolate between two discrete cosines corresponding to neighboring
98	// incoming energies.
99
100	// Sample outgoing cosine bin
101	int n_mu = mu_out_.shape()[1];	1,464✔
102	int k = prn(seed) * n_mu;	1,464✔
103
104	// Rather than use the sampled discrete mu directly, it is smeared over
105	// a bin of width 0.5*min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
106	// discrete mu value itself.
107
108	// Interpolate kth mu value between distributions at energies i and i+1
109	mu = mu_out_(i, k) + f * (mu_out_(i + 1, k) - mu_out_(i, k));	1,464✔
110
111	// Inteprolate (k-1)th mu value between distributions at energies i and i+1.
112	// When k==0, pick a value that will smear the cosine out to a minimum of -1.
113	double mu_left = (k == 0) ? -1.0 - (mu + 1.0)	1,464✔
114	: mu_out_(i, k - 1) +	1,080✔
115	f * (mu_out_(i + 1, k - 1) - mu_out_(i, k - 1));	1,080✔
116
117	// Inteprolate (k+1)th mu value between distributions at energies i and i+1.
118	// When k is the last discrete value, pick a value that will smear the cosine
119	// out to a maximum of 1.
120	double mu_right =
121	(k == n_mu - 1)	1,464✔
122	? 1.0 + (1.0 - mu)	1,464✔
123	: mu_out_(i, k + 1) + f * (mu_out_(i + 1, k + 1) - mu_out_(i, k + 1));	1,080✔
124
125	// Smear cosine
126	mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5);	1,464✔
127
128	// Energy doesn't change in elastic scattering
129	E_out = E_in;	1,464✔
130	}	1,464✔
131
132	//==============================================================================
133	// IncoherentInelasticAEDiscrete implementation
134	//==============================================================================
135
136	IncoherentInelasticAEDiscrete::IncoherentInelasticAEDiscrete(	1,086✔
137	hid_t group, const vector<double>& energy)	1,086✔
138	: energy_ {energy}	1,086✔
139	{
140	read_dataset(group, "energy_out", energy_out_);	1,086✔
141	read_dataset(group, "mu_out", mu_out_);	1,086✔
142	read_dataset(group, "skewed", skewed_);	1,086✔
143	}	1,086✔
144
145	void IncoherentInelasticAEDiscrete::sample(	74,432,481✔
146	double E_in, double& E_out, double& mu, uint64_t* seed) const
147	{
148	// Get index and interpolation factor for inelastic grid
149	int i;
150	double f;
151	get_energy_index(energy_, E_in, i, f);	74,432,481✔
152
153	// Now that we have an incoming energy bin, we need to determine the outgoing
154	// energy bin. This will depend on whether the outgoing energy distribution is
155	// skewed. If it is skewed, then the first two and last two bins have lower
156	// probabilities than the other bins (0.1 for the first and last bins and 0.4
157	// for the second and second to last bins, relative to a normal bin
158	// probability of 1). Otherwise, each bin is equally probable.
159
160	int j;
161	int n = energy_out_.shape()[1];	74,432,481✔
162	if (!skewed_) {	74,432,481✔
163	// All bins equally likely
UNCOV 164	j = prn(seed) * n;	×
165	} else {
166	// Distribution skewed away from edge points
167	double r = prn(seed) * (n - 3);	74,432,481✔
168	if (r > 1.0) {	74,432,481✔
169	// equally likely N-4 middle bins
170	j = r + 1;	73,213,584✔
171	} else if (r > 0.6) {	1,218,897✔
172	// second to last bin has relative probability of 0.4
173	j = n - 2;	488,495✔
174	} else if (r > 0.5) {	730,402✔
175	// last bin has relative probability of 0.1
176	j = n - 1;	121,516✔
177	} else if (r > 0.1) {	608,886✔
178	// second bin has relative probability of 0.4
179	j = 1;	488,352✔
180	} else {
181	// first bin has relative probability of 0.1
182	j = 0;	120,534✔
183	}
184	}
185
186	// Determine outgoing energy corresponding to E_in[i] and E_in[i+1]
187	double E_ij = energy_out_(i, j);	74,432,481✔
188	double E_i1j = energy_out_(i + 1, j);	74,432,481✔
189
190	// Outgoing energy
191	E_out = (1 - f) * E_ij + f * E_i1j;	74,432,481✔
192
193	// Sample outgoing cosine bin
194	int m = mu_out_.shape()[2];	74,432,481✔
195	int k = prn(seed) * m;	74,432,481✔
196
197	// Determine outgoing cosine corresponding to E_in[i] and E_in[i+1]
198	double mu_ijk = mu_out_(i, j, k);	74,432,481✔
199	double mu_i1jk = mu_out_(i + 1, j, k);	74,432,481✔
200
201	// Cosine of angle between incoming and outgoing neutron
202	mu = (1 - f) * mu_ijk + f * mu_i1jk;	74,432,481✔
203	}	74,432,481✔
204
205	//==============================================================================
206	// IncoherentInelasticAE implementation
207	//==============================================================================
208
209	IncoherentInelasticAE::IncoherentInelasticAE(hid_t group)	48✔
210	{
211	// Read correlated angle-energy distribution
212	CorrelatedAngleEnergy dist {group};	48✔
213
214	// Copy incident energies
215	energy_ = dist.energy();	48✔
216
217	// Convert to S(a,b) native format
218	for (const auto& edist : dist.distribution()) {	192✔
219	// Create temporary distribution
220	DistEnergySab d;	144✔
221
222	// Copy outgoing energy distribution
223	d.n_e_out = edist.e_out.size();	144✔
224	d.e_out = edist.e_out;	144✔
225	d.e_out_pdf = edist.p;	144✔
226	d.e_out_cdf = edist.c;	144✔
227
228	for (int j = 0; j < d.n_e_out; ++j) {	720✔
229	auto adist = dynamic_cast<Tabular*>(edist.angle[j].get());	576✔
230	if (adist) {	576✔
231	// On first pass, allocate space for angles
232	if (j == 0) {	576✔
233	auto n_mu = adist->x().size();	144✔
234	d.mu = xt::empty<double>({d.n_e_out, n_mu});	144✔
235	}
236
237	// Copy outgoing angles
238	auto mu_j = xt::view(d.mu, j);	576✔
239	std::copy(adist->x().begin(), adist->x().end(), mu_j.begin());	576✔
240	}	576✔
241	}
242
243	distribution_.emplace_back(std::move(d));	144✔
244	}	144✔
245	}	48✔
246
247	void IncoherentInelasticAE::sample(	3,943,524✔
248	double E_in, double& E_out, double& mu, uint64_t* seed) const
249	{
250	// Get index and interpolation factor for inelastic grid
251	int i;
252	double f;
253	get_energy_index(energy_, E_in, i, f);	3,943,524✔
254
255	// Pick closer energy based on interpolation factor
256	int l = f > 0.5 ? i + 1 : i;	3,943,524✔
257
258	// Determine outgoing energy bin
259	// (First reset n_energy_out to the right value)
260	auto n = distribution_[l].n_e_out;	3,943,524✔
261	double r1 = prn(seed);	3,943,524✔
262	double c_j = distribution_[l].e_out_cdf[0];	3,943,524✔
263	double c_j1;
264	std::size_t j;
265	for (j = 0; j < n - 1; ++j) {	9,514,404✔
266	c_j1 = distribution_[l].e_out_cdf[j + 1];	9,514,404✔
267	if (r1 < c_j1)	9,514,404✔
268	break;	3,943,524✔
269	c_j = c_j1;	5,570,880✔
270	}
271
272	// check to make sure j is <= n_energy_out - 2
273	j = std::min(j, n - 2);	3,943,524✔
274
275	// Get the data to interpolate between
276	double E_l_j = distribution_[l].e_out[j];	3,943,524✔
277	double p_l_j = distribution_[l].e_out_pdf[j];	3,943,524✔
278
279	// Next part assumes linear-linear interpolation in standard
280	double E_l_j1 = distribution_[l].e_out[j + 1];	3,943,524✔
281	double p_l_j1 = distribution_[l].e_out_pdf[j + 1];	3,943,524✔
282
283	// Find secondary energy (variable E)
284	double frac = (p_l_j1 - p_l_j) / (E_l_j1 - E_l_j);	3,943,524✔
285	if (frac == 0.0) {	3,943,524✔
286	E_out = E_l_j + (r1 - c_j) / p_l_j;	3,943,524✔
287	} else {
UNCOV 288	E_out = E_l_j +	×
289	(std::sqrt(std::max(0.0, p_l_j * p_l_j + 2.0 * frac * (r1 - c_j))) -	×
290	p_l_j) /	×
291	frac;
292	}
293
294	// Adjustment of outgoing energy
295	double E_l = energy_[l];	3,943,524✔
296	if (E_out < 0.5 * E_l) {	3,943,524✔
297	E_out = 2.0 E_in / E_l - 1.0;	375,792✔
298	} else {
299	E_out += E_in - E_l;	3,567,732✔
300	}
301
302	// Sample outgoing cosine bin
303	int n_mu = distribution_[l].mu.shape()[1];	3,943,524✔
304	std::size_t k = prn(seed) * n_mu;	3,943,524✔
305
306	// Rather than use the sampled discrete mu directly, it is smeared over
307	// a bin of width 0.5*min(mu[k] - mu[k-1], mu[k+1] - mu[k]) centered on the
308	// discrete mu value itself.
309	const auto& mu_l = distribution_[l].mu;	3,943,524✔
310	f = (r1 - c_j) / (c_j1 - c_j);	3,943,524✔
311
312	// Interpolate kth mu value between distributions at energies j and j+1
313	mu = mu_l(j, k) + f * (mu_l(j + 1, k) - mu_l(j, k));	3,943,524✔
314
315	// Inteprolate (k-1)th mu value between distributions at energies j and j+1.
316	// When k==0, pick a value that will smear the cosine out to a minimum of -1.
317	double mu_left =
318	(k == 0)
319	? mu_left = -1.0 - (mu + 1.0)	3,943,524✔
320	: mu_left = mu_l(j, k - 1) + f * (mu_l(j + 1, k - 1) - mu_l(j, k - 1));	3,446,688✔
321
322	// Inteprolate (k+1)th mu value between distributions at energies j and j+1.
323	// When k is the last discrete value, pick a value that will smear the cosine
324	// out to a maximum of 1.
325	double mu_right =
326	(k == n_mu - 1)	3,943,524✔
327	? mu_right = 1.0 + (1.0 - mu)	3,943,524✔
328	: mu_right = mu_l(j, k + 1) + f * (mu_l(j + 1, k + 1) - mu_l(j, k + 1));	3,448,104✔
329
330	// Smear cosine
331	mu += std::min(mu - mu_left, mu_right - mu) * (prn(seed) - 0.5);	3,943,524✔
332	}	3,943,524✔
333
334	//==============================================================================
335	// MixedElasticAE implementation
336	//==============================================================================
337
338	MixedElasticAE::MixedElasticAE(	48✔
339	hid_t group, const CoherentElasticXS& coh_xs, const Function1D& incoh_xs)	48✔
340	: coherent_dist_(coh_xs), coherent_xs_(coh_xs), incoherent_xs_(incoh_xs)	48✔
341	{
342	// Read incoherent elastic distribution
343	hid_t incoherent_group = open_group(group, "incoherent");	48✔
344	std::string temp;	48✔
345	read_attribute(incoherent_group, "type", temp);	48✔
346	if (temp == "incoherent_elastic") {	48✔
UNCOV 347	incoherent_dist_ = make_unique<IncoherentElasticAE>(incoherent_group);	×
348	} else if (temp == "incoherent_elastic_discrete") {	48✔
349	auto xs = dynamic_cast<const Tabulated1D*>(&incoh_xs);	48✔
350	incoherent_dist_ =
351	make_unique<IncoherentElasticAEDiscrete>(incoherent_group, xs->x());	48✔
352	}
353	close_group(incoherent_group);	48✔
354	}	48✔
355
356	void MixedElasticAE::sample(	50,712✔
357	double E_in, double& E_out, double& mu, uint64_t* seed) const
358	{
359	// Evaluate coherent and incoherent elastic cross sections
360	double xs_coh = coherent_xs_(E_in);	50,712✔
361	double xs_incoh = incoherent_xs_(E_in);	50,712✔
362
363	if (prn(seed) * (xs_coh + xs_incoh) < xs_coh) {	50,712✔
364	coherent_dist_.sample(E_in, E_out, mu, seed);	49,248✔
365	} else {
366	incoherent_dist_->sample(E_in, E_out, mu, seed);	1,464✔
367	}
368	}	50,712✔
369
370	} // namespace openmc

openmc-dev / openmc / 13461947318

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