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

87.83

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

#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
  tensor::Tensor<int> temp;
  read_attribute(dset, "interpolation", temp);

  tensor::View<int> temp_b = temp.slice(0); // breakpoints
  tensor::View<int> temp_i = temp.slice(1); // 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);

  tensor::Tensor<double> eout;
  read_dataset(dset, eout);
  close_dataset(dset);

  // Read angle distributions
  tensor::Tensor<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 = eout.slice(0, tensor::range(j, j + n));
    d.p = eout.slice(1, tensor::range(j, j + n));
    d.c = eout.slice(2, tensor::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);
      tensor::View<double> xs =
        mu.slice(0, tensor::range(offset_mu, offset_mu + m));
      tensor::View<double> ps =
        mu.slice(1, tensor::range(offset_mu, offset_mu + m));
      tensor::View<double> cs =
        mu.slice(2, tensor::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
}
Distribution& CorrelatedAngleEnergy::sample_dist(
  double E_in, double& E_out, 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) {
    return *distribution_[l].angle[k];
  } else {
    return *distribution_[l].angle[k + 1];
  }
}

void CorrelatedAngleEnergy::sample(
  double E_in, double& E_out, double& mu, uint64_t* seed) const
{
  mu = sample_dist(E_in, E_out, seed).sample(seed).first;
}

double CorrelatedAngleEnergy::sample_energy_and_pdf(double E_in, double mu,
  double& E_out, uint64_t* seed, bool is_com, double awr) const
{
  auto& dist = sample_dist(E_in, E_out, seed);
  double jac = 1.0;
  if (is_com)
    jac = get_jac_and_transform(E_in, mu, E_out, seed, awr);

  return jac * dist.evaluate(mu);
}

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