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

92.31

/src/secondary_correlated.cpp

#include "openmc/secondary_correlated.h"

#include <algorithm> // for copy
#include <cmath>
#include <cstddef>  // for size_t
#include <iterator> // for back_inserter

#include "xtensor/xarray.hpp"
#include "xtensor/xview.hpp"

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

namespace openmc {

//==============================================================================
//! CorrelatedAngleEnergy implementation
//==============================================================================

CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group)
{
  // Open incoming energy dataset
  hid_t dset = open_dataset(group, "energy");

  // Get interpolation parameters
  xt::xarray<int> temp;
  read_attribute(dset, "interpolation", temp);

  auto temp_b = xt::view(temp, 0); // view of breakpoints
  auto temp_i = xt::view(temp, 1); // view of interpolation parameters

  std::copy(temp_b.begin(), temp_b.end(), std::back_inserter(breakpoints_));
  for (const auto i : temp_i)
    interpolation_.push_back(int2interp(i));
  n_region_ = breakpoints_.size();

  // Get incoming energies
  read_dataset(dset, energy_);
  std::size_t n_energy = energy_.size();
  close_dataset(dset);

  // Get outgoing energy distribution data
  dset = open_dataset(group, "energy_out");
  vector<int> offsets;
  vector<int> interp;
  vector<int> n_discrete;
  read_attribute(dset, "offsets", offsets);
  read_attribute(dset, "interpolation", interp);
  read_attribute(dset, "n_discrete_lines", n_discrete);

  xt::xarray<double> eout;
  read_dataset(dset, eout);
  close_dataset(dset);

  // Read angle distributions
  xt::xarray<double> mu;
  read_dataset(group, "mu", mu);

  for (int i = 0; i < n_energy; ++i) {
    // Determine number of outgoing energies
    int j = offsets[i];
    int n;
    if (i < n_energy - 1) {
      n = offsets[i + 1] - j;
    } else {
      n = eout.shape()[1] - j;
    }

    // Assign interpolation scheme and number of discrete lines
    CorrTable d;
    d.interpolation = int2interp(interp[i]);
    d.n_discrete = n_discrete[i];

    // Copy data
    d.e_out = xt::view(eout, 0, xt::range(j, j + n));
    d.p = xt::view(eout, 1, xt::range(j, j + n));
    d.c = xt::view(eout, 2, xt::range(j, j + n));

    // To get answers that match ACE data, for now we still use the tabulated
    // CDF values that were passed through to the HDF5 library. At a later
    // time, we can remove the CDF values from the HDF5 library and
    // reconstruct them using the PDF
    if (false) {
      // Calculate cumulative distribution function -- discrete portion
      for (int k = 0; k < d.n_discrete; ++k) {
        if (k == 0) {
          d.c[k] = d.p[k];
        } else {
          d.c[k] = d.c[k - 1] + d.p[k];
        }
      }

      // Continuous portion
      for (int k = d.n_discrete; k < n; ++k) {
        if (k == d.n_discrete) {
          d.c[k] = d.c[k - 1] + d.p[k];
        } else {
          if (d.interpolation == Interpolation::histogram) {
            d.c[k] = d.c[k - 1] + d.p[k - 1] * (d.e_out[k] - d.e_out[k - 1]);
          } else if (d.interpolation == Interpolation::lin_lin) {
            d.c[k] = d.c[k - 1] + 0.5 * (d.p[k - 1] + d.p[k]) *
                                    (d.e_out[k] - d.e_out[k - 1]);
          }
        }
      }

      // Normalize density and distribution functions
      d.p /= d.c[n - 1];
      d.c /= d.c[n - 1];
    }

    for (j = 0; j < n; ++j) {
      // Get interpolation scheme
      int interp_mu = std::lround(eout(3, offsets[i] + j));

      // Determine offset and size of distribution
      int offset_mu = std::lround(eout(4, offsets[i] + j));
      int m;
      if (offsets[i] + j + 1 < eout.shape()[1]) {
        m = std::lround(eout(4, offsets[i] + j + 1)) - offset_mu;
      } else {
        m = mu.shape()[1] - offset_mu;
      }

      // For incoherent inelastic thermal scattering, the angle distributions
      // may be given as discrete mu values. In this case, interpolation values
      // of zero appear in the HDF5 file. Here we change it to a 1 so that
      // int2interp doesn't fail.
      if (interp_mu == 0)
        interp_mu = 1;

      auto interp = int2interp(interp_mu);
      auto xs = xt::view(mu, 0, xt::range(offset_mu, offset_mu + m));
      auto ps = xt::view(mu, 1, xt::range(offset_mu, offset_mu + m));
      auto cs = xt::view(mu, 2, xt::range(offset_mu, offset_mu + m));

      vector<double> x {xs.begin(), xs.end()};
      vector<double> p {ps.begin(), ps.end()};
      vector<double> c {cs.begin(), cs.end()};

      // To get answers that match ACE data, for now we still use the tabulated
      // CDF values that were passed through to the HDF5 library. At a later
      // time, we can remove the CDF values from the HDF5 library and
      // reconstruct them using the PDF
      Tabular* mudist = new Tabular {x.data(), p.data(), m, interp, c.data()};

      d.angle.emplace_back(mudist);
    } // outgoing energies

    distribution_.push_back(std::move(d));
  } // incoming energies
}

void CorrelatedAngleEnergy::sample(
  double E_in, double& E_out, double& mu, uint64_t* seed) const
{
  // Find energy bin and calculate interpolation factor
  int i;
  double r;
  get_energy_index(energy_, E_in, i, r);

  // Sample between the ith and [i+1]th bin
  int l = r > prn(seed) ? i + 1 : i;

  // Interpolation for energy E1 and EK
  int n_energy_out = distribution_[i].e_out.size();
  int n_discrete = distribution_[i].n_discrete;
  double E_i_1 = distribution_[i].e_out[n_discrete];
  double E_i_K = distribution_[i].e_out[n_energy_out - 1];

  n_energy_out = distribution_[i + 1].e_out.size();
  n_discrete = distribution_[i + 1].n_discrete;
  double E_i1_1 = distribution_[i + 1].e_out[n_discrete];
  double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];

  double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);
  double E_K = E_i_K + r * (E_i1_K - E_i_K);

  // Determine outgoing energy bin
  n_energy_out = distribution_[l].e_out.size();
  n_discrete = distribution_[l].n_discrete;
  double r1 = prn(seed);
  double c_k = distribution_[l].c[0];
  int k = 0;
  int end = n_energy_out - 2;

  // Discrete portion
  for (int j = 0; j < n_discrete; ++j) {
    k = j;
    c_k = distribution_[l].c[k];
    if (r1 < c_k) {
      end = j;
      break;
    }
  }

  // Continuous portion
  double c_k1;
  for (int j = n_discrete; j < end; ++j) {
    k = j;
    c_k1 = distribution_[l].c[k + 1];
    if (r1 < c_k1)
      break;
    k = j + 1;
    c_k = c_k1;
  }

  double E_l_k = distribution_[l].e_out[k];
  double p_l_k = distribution_[l].p[k];
  if (distribution_[l].interpolation == Interpolation::histogram) {
    // Histogram interpolation
    if (p_l_k > 0.0 && k >= n_discrete) {
      E_out = E_l_k + (r1 - c_k) / p_l_k;
    } else {
      E_out = E_l_k;
    }

  } else if (distribution_[l].interpolation == Interpolation::lin_lin) {
    // Linear-linear interpolation
    double E_l_k1 = distribution_[l].e_out[k + 1];
    double p_l_k1 = distribution_[l].p[k + 1];

    double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k);
    if (frac == 0.0) {
      E_out = E_l_k + (r1 - c_k) / p_l_k;
    } else {
      E_out =
        E_l_k +
        (std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) -
          p_l_k) /
          frac;
    }
  }

  // Now interpolate between incident energy bins i and i + 1
  if (k >= n_discrete) {
    if (l == i) {
      E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);
    } else {
      E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);
    }
  }

  // Find correlated angular distribution for closest outgoing energy bin
  if (r1 - c_k < c_k1 - r1 ||
      distribution_[l].interpolation == Interpolation::histogram) {
    mu = distribution_[l].angle[k]->sample(seed);
  } else {
    mu = distribution_[l].angle[k + 1]->sample(seed);
  }
}

} // namespace openmc

1	#include "openmc/secondary_correlated.h"
2
3	#include <algorithm> // for copy
4	#include <cmath>
5	#include <cstddef> // for size_t
6	#include <iterator> // for back_inserter
7
8	#include "xtensor/xarray.hpp"
9	#include "xtensor/xview.hpp"
10
11	#include "openmc/endf.h"
12	#include "openmc/hdf5_interface.h"
13	#include "openmc/math_functions.h"
14	#include "openmc/random_lcg.h"
15	#include "openmc/search.h"
16
17	namespace openmc {
18
19	//==============================================================================
20	//! CorrelatedAngleEnergy implementation
21	//==============================================================================
22
23	CorrelatedAngleEnergy::CorrelatedAngleEnergy(hid_t group)	46,579✔
24	{
25	// Open incoming energy dataset
26	hid_t dset = open_dataset(group, "energy");	46,579✔
27
28	// Get interpolation parameters
29	xt::xarray<int> temp;	46,579✔
30	read_attribute(dset, "interpolation", temp);	46,579✔
31
32	auto temp_b = xt::view(temp, 0); // view of breakpoints	46,579✔
33	auto temp_i = xt::view(temp, 1); // view of interpolation parameters	46,579✔
34
35	std::copy(temp_b.begin(), temp_b.end(), std::back_inserter(breakpoints_));	46,579✔
36	for (const auto i : temp_i)	93,158✔
37	interpolation_.push_back(int2interp(i));	46,579✔
38	n_region_ = breakpoints_.size();	46,579✔
39
40	// Get incoming energies
41	read_dataset(dset, energy_);	46,579✔
42	std::size_t n_energy = energy_.size();	46,579✔
43	close_dataset(dset);	46,579✔
44
45	// Get outgoing energy distribution data
46	dset = open_dataset(group, "energy_out");	46,579✔
47	vector<int> offsets;	46,579✔
48	vector<int> interp;	46,579✔
49	vector<int> n_discrete;	46,579✔
50	read_attribute(dset, "offsets", offsets);	46,579✔
51	read_attribute(dset, "interpolation", interp);	46,579✔
52	read_attribute(dset, "n_discrete_lines", n_discrete);	46,579✔
53
54	xt::xarray<double> eout;	46,579✔
55	read_dataset(dset, eout);	46,579✔
56	close_dataset(dset);	46,579✔
57
58	// Read angle distributions
59	xt::xarray<double> mu;	46,579✔
60	read_dataset(group, "mu", mu);	46,579✔
61
62	for (int i = 0; i < n_energy; ++i) {	957,418✔
63	// Determine number of outgoing energies
64	int j = offsets[i];	910,839✔
65	int n;
66	if (i < n_energy - 1) {	910,839✔
67	n = offsets[i + 1] - j;	864,260✔
68	} else {
69	n = eout.shape()[1] - j;	46,579✔
70	}
71
72	// Assign interpolation scheme and number of discrete lines
73	CorrTable d;	910,839✔
74	d.interpolation = int2interp(interp[i]);	910,839✔
75	d.n_discrete = n_discrete[i];	910,839✔
76
77	// Copy data
78	d.e_out = xt::view(eout, 0, xt::range(j, j + n));	910,839✔
79	d.p = xt::view(eout, 1, xt::range(j, j + n));	910,839✔
80	d.c = xt::view(eout, 2, xt::range(j, j + n));	910,839✔
81
82	// To get answers that match ACE data, for now we still use the tabulated
83	// CDF values that were passed through to the HDF5 library. At a later
84	// time, we can remove the CDF values from the HDF5 library and
85	// reconstruct them using the PDF
86	if (false) {
87	// Calculate cumulative distribution function -- discrete portion
88	for (int k = 0; k < d.n_discrete; ++k) {
89	if (k == 0) {
90	d.c[k] = d.p[k];
91	} else {
92	d.c[k] = d.c[k - 1] + d.p[k];
93	}
94	}
95
96	// Continuous portion
97	for (int k = d.n_discrete; k < n; ++k) {
98	if (k == d.n_discrete) {
99	d.c[k] = d.c[k - 1] + d.p[k];
100	} else {
101	if (d.interpolation == Interpolation::histogram) {
102	d.c[k] = d.c[k - 1] + d.p[k - 1] * (d.e_out[k] - d.e_out[k - 1]);
103	} else if (d.interpolation == Interpolation::lin_lin) {
104	d.c[k] = d.c[k - 1] + 0.5 * (d.p[k - 1] + d.p[k]) *
105	(d.e_out[k] - d.e_out[k - 1]);
106	}
107	}
108	}
109
110	// Normalize density and distribution functions
111	d.p /= d.c[n - 1];
112	d.c /= d.c[n - 1];
113	}
114
115	for (j = 0; j < n; ++j) {	21,672,375✔
116	// Get interpolation scheme
117	int interp_mu = std::lround(eout(3, offsets[i] + j));	20,761,536✔
118
119	// Determine offset and size of distribution
120	int offset_mu = std::lround(eout(4, offsets[i] + j));	20,761,536✔
121	int m;
122	if (offsets[i] + j + 1 < eout.shape()[1]) {	20,761,536✔
123	m = std::lround(eout(4, offsets[i] + j + 1)) - offset_mu;	20,714,957✔
124	} else {
125	m = mu.shape()[1] - offset_mu;	46,579✔
126	}
127
128	// For incoherent inelastic thermal scattering, the angle distributions
129	// may be given as discrete mu values. In this case, interpolation values
130	// of zero appear in the HDF5 file. Here we change it to a 1 so that
131	// int2interp doesn't fail.
132	if (interp_mu == 0)	20,761,536✔
133	interp_mu = 1;	528✔
134
135	auto interp = int2interp(interp_mu);	20,761,536✔
136	auto xs = xt::view(mu, 0, xt::range(offset_mu, offset_mu + m));	20,761,536✔
137	auto ps = xt::view(mu, 1, xt::range(offset_mu, offset_mu + m));	20,761,536✔
138	auto cs = xt::view(mu, 2, xt::range(offset_mu, offset_mu + m));	20,761,536✔
139
140	vector<double> x {xs.begin(), xs.end()};	20,761,536✔
141	vector<double> p {ps.begin(), ps.end()};	20,761,536✔
142	vector<double> c {cs.begin(), cs.end()};	20,761,536✔
143
144	// To get answers that match ACE data, for now we still use the tabulated
145	// CDF values that were passed through to the HDF5 library. At a later
146	// time, we can remove the CDF values from the HDF5 library and
147	// reconstruct them using the PDF
148	Tabular* mudist = new Tabular {x.data(), p.data(), m, interp, c.data()};	20,761,536✔
149
150	d.angle.emplace_back(mudist);	20,761,536✔
151	} // outgoing energies	20,761,536✔
152
153	distribution_.push_back(std::move(d));	910,839✔
154	} // incoming energies	910,839✔
155	}	46,579✔
156
157	void CorrelatedAngleEnergy::sample(	303,642✔
158	double E_in, double& E_out, double& mu, uint64_t* seed) const
159	{
160	// Find energy bin and calculate interpolation factor
161	int i;
162	double r;
163	get_energy_index(energy_, E_in, i, r);	303,642✔
164
165	// Sample between the ith and [i+1]th bin
166	int l = r > prn(seed) ? i + 1 : i;	303,642✔
167
168	// Interpolation for energy E1 and EK
169	int n_energy_out = distribution_[i].e_out.size();	303,642✔
170	int n_discrete = distribution_[i].n_discrete;	303,642✔
171	double E_i_1 = distribution_[i].e_out[n_discrete];	303,642✔
172	double E_i_K = distribution_[i].e_out[n_energy_out - 1];	303,642✔
173
174	n_energy_out = distribution_[i + 1].e_out.size();	303,642✔
175	n_discrete = distribution_[i + 1].n_discrete;	303,642✔
176	double E_i1_1 = distribution_[i + 1].e_out[n_discrete];	303,642✔
177	double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];	303,642✔
178
179	double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);	303,642✔
180	double E_K = E_i_K + r * (E_i1_K - E_i_K);	303,642✔
181
182	// Determine outgoing energy bin
183	n_energy_out = distribution_[l].e_out.size();	303,642✔
184	n_discrete = distribution_[l].n_discrete;	303,642✔
185	double r1 = prn(seed);	303,642✔
186	double c_k = distribution_[l].c[0];	303,642✔
187	int k = 0;	303,642✔
188	int end = n_energy_out - 2;	303,642✔
189
190	// Discrete portion
191	for (int j = 0; j < n_discrete; ++j) {	303,642!
UNCOV 192	k = j;	×
UNCOV 193	c_k = distribution_[l].c[k];	×
UNCOV 194	if (r1 < c_k) {	×
UNCOV 195	end = j;	×
UNCOV 196	break;	×
197	}
198	}
199
200	// Continuous portion
201	double c_k1;
202	for (int j = n_discrete; j < end; ++j) {	2,105,189✔
203	k = j;	2,104,167✔
204	c_k1 = distribution_[l].c[k + 1];	2,104,167✔
205	if (r1 < c_k1)	2,104,167✔
206	break;	302,620✔
207	k = j + 1;	1,801,547✔
208	c_k = c_k1;	1,801,547✔
209	}
210
211	double E_l_k = distribution_[l].e_out[k];	303,642✔
212	double p_l_k = distribution_[l].p[k];	303,642✔
213	if (distribution_[l].interpolation == Interpolation::histogram) {	303,642✔
214	// Histogram interpolation
215	if (p_l_k > 0.0 && k >= n_discrete) {	21,923!
216	E_out = E_l_k + (r1 - c_k) / p_l_k;	21,923✔
217	} else {
UNCOV 218	E_out = E_l_k;	×
219	}
220
221	} else if (distribution_[l].interpolation == Interpolation::lin_lin) {	281,719!
222	// Linear-linear interpolation
223	double E_l_k1 = distribution_[l].e_out[k + 1];	281,719✔
224	double p_l_k1 = distribution_[l].p[k + 1];	281,719✔
225
226	double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k);	281,719✔
227	if (frac == 0.0) {	281,719!
228	E_out = E_l_k + (r1 - c_k) / p_l_k;	×
229	} else {
230	E_out =	281,719✔
231	E_l_k +	281,719✔
232	(std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) -	281,719✔
233	p_l_k) /	281,719✔
234	frac;
235	}
236	}
237
238	// Now interpolate between incident energy bins i and i + 1
239	if (k >= n_discrete) {	303,642!
240	if (l == i) {	303,642✔
241	E_out = E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);	84,856✔
242	} else {
243	E_out = E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);	218,786✔
244	}
245	}
246
247	// Find correlated angular distribution for closest outgoing energy bin
248	if (r1 - c_k < c_k1 - r1 \|\|	455,040✔
249	distribution_[l].interpolation == Interpolation::histogram) {	151,398✔
250	mu = distribution_[l].angle[k]->sample(seed);	163,684✔
251	} else {
252	mu = distribution_[l].angle[k + 1]->sample(seed);	139,958✔
253	}
254	}	303,642✔
255
256	} // 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