20501343020

Committed 25 Dec 2025 07:38AM UTC coverage: 82.185% (+0.05%) from 82.139%

Build # 20501343020

Build Type

Pull #3692

github

Committed by

web-flow

Commit Message

Merge fc35e8018 into 3f06a42ab

Pull Request Pull Request #3692: Fix a bug in rotational periodic boundary conditions

Run Details

17096 of 23665 branches covered (72.24%)

Branch coverage included in aggregate %.

36 of 37 new or added lines in 3 files covered. (97.3%)

198 existing lines in 11 files now uncovered.

55209 of 64313 relevant lines covered (85.84%)

43488756.38 hits per line

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

89.77

/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 "xtensor/xview.hpp"

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

namespace openmc {

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

openmc-dev / openmc / 20501343020

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