34416149688

Committed 09 Sep 2026 11:15PM UTC coverage: 81.48% (+0.1%) from 81.347%

Build # 34416149688

Build Type

Pull #4087

github

Committed by

web-flow

Commit Message

Merge 93138e589 into 5260b9a0f

Pull Request Pull Request #4087: Compute bounding boxes for general planes and tori in C++

Coverage Stats

18756 of 27197 branches covered (68.96%)

Branch coverage included in aggregate %.

43 of 46 new or added lines in 2 files covered. (93.48%)

959 existing lines in 29 files now uncovered.

60863 of 70519 relevant lines covered (86.31%)

49704801.91 hits per line

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

94.44

/src/distribution_energy.cpp

#include "openmc/distribution_energy.h"

#include <algorithm> // for max, min, copy, move
#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_dist.h"
#include "openmc/random_lcg.h"

namespace openmc {

//==============================================================================
// DiscretePhoton implementation
//==============================================================================

DiscretePhoton::DiscretePhoton(hid_t group)
{
  read_attribute(group, "primary_flag", primary_flag_);
  read_attribute(group, "energy", energy_);
  read_attribute(group, "atomic_weight_ratio", A_);
}

double DiscretePhoton::sample(double E, uint64_t* seed) const
{
  if (primary_flag_ == 2) {
    return energy_ + A_ / (A_ + 1) * E;
  } else {
    return energy_;
  }
}

//==============================================================================
// LevelInelastic implementation
//==============================================================================

LevelInelastic::LevelInelastic(hid_t group)
{
  read_attribute(group, "threshold", threshold_);
  read_attribute(group, "mass_ratio", mass_ratio_);
}

double LevelInelastic::sample(double E, uint64_t* seed) const
{
  return mass_ratio_ * (E - threshold_);
}

//==============================================================================
// ContinuousTabular implementation
//==============================================================================

ContinuousTabular::ContinuousTabular(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, "distribution");
  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);

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

    // 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 (true) {
      d.c = eout.slice(2, tensor::range(j, j + n));
    } else {
      // 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];
    }

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

double ContinuousTabular::sample(double E, uint64_t* seed) const
{
  // Read number of interpolation regions and incoming energies
  bool histogram_interp;
  if (n_region_ == 1) {
    histogram_interp = (interpolation_[0] == Interpolation::histogram);
  } else {
    histogram_interp = false;
  }

  // Find energy bin and calculate interpolation factor -- if the energy is
  // outside the range of the tabulated energies, choose the first or last bins
  int i;
  double r;
  get_energy_index(energy_, E, i, r);

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

  // Determine outgoing energy bin
  int n_energy_out = distribution_[l].e_out.size();
  int 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];

  if (k < n_discrete) {
    // Discrete case
    return E_l_k;
  } else {
    // Continuous case
    double p_l_k = distribution_[l].p[k];
    double E_out;
    if (distribution_[l].interpolation == Interpolation::histogram) {
      // Histogram interpolation
      if (p_l_k > 0.0) {
        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];

      if (E_l_k != E_l_k1) {
        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;
        }
      } else {
        E_out = E_l_k;
      }
    } else {
      throw std::runtime_error {
        "Unexpected interpolation for continuous energy "
        "distribution."};
    }

    // Now interpolate between incident energy bins i and i + 1
    if (!histogram_interp && n_energy_out > 1) {
      // Interpolation for energy E1 and EK
      n_energy_out = distribution_[i].e_out.size();
      n_discrete = distribution_[i].n_discrete;
      const double E_i_1 = distribution_[i].e_out[n_discrete];
      const 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;
      const double E_i1_1 = distribution_[i + 1].e_out[n_discrete];
      const double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];

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

      if (l == i) {
        return E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);
      } else {
        return E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);
      }
    } else {
      return E_out;
    }
  }
}

//==============================================================================
// MaxwellEnergy implementation
//==============================================================================

MaxwellEnergy::MaxwellEnergy(hid_t group)
{
  read_attribute(group, "u", u_);
  hid_t dset = open_dataset(group, "theta");
  theta_ = Tabulated1D {dset};
  close_dataset(dset);
}

double MaxwellEnergy::sample(double E, uint64_t* seed) const
{
  // Get temperature corresponding to incoming energy
  double theta = theta_(E);

  while (true) {
    // Sample maxwell fission spectrum
    double E_out = maxwell_spectrum(theta, seed);

    // Accept energy based on restriction energy
    if (E_out <= E - u_)
      return E_out;
  }
}

//==============================================================================
// Evaporation implementation
//==============================================================================

Evaporation::Evaporation(hid_t group)
{
  read_attribute(group, "u", u_);
  hid_t dset = open_dataset(group, "theta");
  theta_ = Tabulated1D {dset};
  close_dataset(dset);
}

double Evaporation::sample(double E, uint64_t* seed) const
{
  // Get temperature corresponding to incoming energy
  double theta = theta_(E);

  double y = (E - u_) / theta;
  double v = 1.0 - std::exp(-y);

  // Sample outgoing energy based on evaporation spectrum probability
  // density function
  double x;
  while (true) {
    x = -std::log((1.0 - v * prn(seed)) * (1.0 - v * prn(seed)));
    if (x <= y)
      break;
  }

  return x * theta;
}

//==============================================================================
// WattEnergy implementation
//==============================================================================

WattEnergy::WattEnergy(hid_t group)
{
  // Read restriction energy
  read_attribute(group, "u", u_);

  // Read tabulated functions
  hid_t dset = open_dataset(group, "a");
  a_ = Tabulated1D {dset};
  close_dataset(dset);
  dset = open_dataset(group, "b");
  b_ = Tabulated1D {dset};
  close_dataset(dset);
}

double WattEnergy::sample(double E, uint64_t* seed) const
{
  // Determine Watt parameters at incident energy
  double a = a_(E);
  double b = b_(E);

  while (true) {
    // Sample energy-dependent Watt fission spectrum
    double E_out = watt_spectrum(a, b, seed);

    // Accept energy based on restriction energy
    if (E_out <= E - u_)
      return E_out;
  }
}

} // namespace openmc

1	#include "openmc/distribution_energy.h"
2
3	#include <algorithm> // for max, min, copy, move
4	#include <cstddef> // for size_t
5	#include <iterator> // for back_inserter
6
7	#include "openmc/tensor.h"
8
9	#include "openmc/endf.h"
10	#include "openmc/hdf5_interface.h"
11	#include "openmc/math_functions.h"
12	#include "openmc/random_dist.h"
13	#include "openmc/random_lcg.h"
14
15	namespace openmc {
16
17	//==============================================================================
18	// DiscretePhoton implementation
19	//==============================================================================
20
21	DiscretePhoton::DiscretePhoton(hid_t group)	2,734,640 ✔
22	{
23	read_attribute(group, "primary_flag", primary_flag_);	2,734,640 ✔
24	read_attribute(group, "energy", energy_);	2,734,640 ✔
25	read_attribute(group, "atomic_weight_ratio", A_);	2,734,640 ✔
26	}	2,734,640 ✔
27
28	double DiscretePhoton::sample(double E, uint64_t* seed) const	1,235,663 ✔
29	{
30	if (primary_flag_ == 2) {	1,235,663 ✔
31	return energy_ + A_ / (A_ + 1) * E;	282,073 ✔
32	} else {
33	return energy_;	953,590 ✔
34	}
35	}
36
37	//==============================================================================
38	// LevelInelastic implementation
39	//==============================================================================
40
41	LevelInelastic::LevelInelastic(hid_t group)	579,419 ✔
42	{
43	read_attribute(group, "threshold", threshold_);	579,419 ✔
44	read_attribute(group, "mass_ratio", mass_ratio_);	579,419 ✔
45	}	579,419 ✔
46
47	double LevelInelastic::sample(double E, uint64_t* seed) const	17,954,629 ✔
48	{
49	return mass_ratio_ * (E - threshold_);	17,954,629 ✔
50	}
51
52	//==============================================================================
53	// ContinuousTabular implementation
54	//==============================================================================
55
56	ContinuousTabular::ContinuousTabular(hid_t group)	376,026 ✔
57	{
58	// Open incoming energy dataset
59	hid_t dset = open_dataset(group, "energy");	376,026 ✔
60
61	// Get interpolation parameters
62	tensor::Tensor<int> temp;	376,026 ✔
63	read_attribute(dset, "interpolation", temp);	376,026 ✔
64
65	tensor::View<int> temp_b = temp.slice(0); // breakpoints	376,026 ✔
66	tensor::View<int> temp_i = temp.slice(1); // interpolation parameters	376,026 ✔
67
68	std::copy(temp_b.begin(), temp_b.end(), std::back_inserter(breakpoints_));	376,026 ✔
69	for (const auto i : temp_i)	752,284 ✔
70	interpolation_.push_back(int2interp(i));	376,258 ✔
71	n_region_ = breakpoints_.size();	376,026 ✔
72
73	// Get incoming energies
74	read_dataset(dset, energy_);	376,026 ✔
75	std::size_t n_energy = energy_.size();	376,026 ✔
76	close_dataset(dset);	376,026 ✔
77
78	// Get outgoing energy distribution data
79	dset = open_dataset(group, "distribution");	376,026 ✔
80	vector<int> offsets;	376,026 ✔
81	vector<int> interp;	376,026 ✔
82	vector<int> n_discrete;	376,026 ✔
83	read_attribute(dset, "offsets", offsets);	376,026 ✔
84	read_attribute(dset, "interpolation", interp);	376,026 ✔
85	read_attribute(dset, "n_discrete_lines", n_discrete);	376,026 ✔
86
87	tensor::Tensor<double> eout;	376,026 ✔
88	read_dataset(dset, eout);	376,026 ✔
89	close_dataset(dset);	376,026 ✔
90
91	for (int i = 0; i < n_energy; ++i) {	5,356,230 ✔
92	// Determine number of outgoing energies
93	int j = offsets[i];	4,980,204 ✔
94	int n;	4,980,204 ✔
95	if (i < n_energy - 1) {	4,980,204 ✔
96	n = offsets[i + 1] - j;	4,604,178 ✔
97	} else {
98	n = eout.shape(1) - j;	752,052 !
99	}
100
101	// Assign interpolation scheme and number of discrete lines
102	CTTable d;	4,980,204 ✔
103	d.interpolation = int2interp(interp[i]);	4,980,204 ✔
104	d.n_discrete = n_discrete[i];	4,980,204 ✔
105
106	// Copy data
107	d.e_out = eout.slice(0, tensor::range(j, j + n));	4,980,204 ✔
108	d.p = eout.slice(1, tensor::range(j, j + n));	4,980,204 ✔
109
110	// To get answers that match ACE data, for now we still use the tabulated
111	// CDF values that were passed through to the HDF5 library. At a later
112	// time, we can remove the CDF values from the HDF5 library and
113	// reconstruct them using the PDF
114	if (true) {	4,980,204 ✔
115	d.c = eout.slice(2, tensor::range(j, j + n));	4,980,204 ✔
116	} else {
117	// Calculate cumulative distribution function -- discrete portion
118	for (int k = 0; k < d.n_discrete; ++k) {
119	if (k == 0) {
120	d.c[k] = d.p[k];
121	} else {
122	d.c[k] = d.c[k - 1] + d.p[k];
123	}
124	}
125
126	// Continuous portion
127	for (int k = d.n_discrete; k < n; ++k) {
128	if (k == d.n_discrete) {
129	d.c[k] = d.c[k - 1] + d.p[k];
130	} else {
131	if (d.interpolation == Interpolation::histogram) {
132	d.c[k] = d.c[k - 1] + d.p[k - 1] * (d.e_out[k] - d.e_out[k - 1]);
133	} else if (d.interpolation == Interpolation::lin_lin) {
134	d.c[k] = d.c[k - 1] + 0.5 * (d.p[k - 1] + d.p[k]) *
135	(d.e_out[k] - d.e_out[k - 1]);
136	}
137	}
138	}
139
140	// Normalize density and distribution functions
141	d.p /= d.c[n - 1];
142	d.c /= d.c[n - 1];
143	}
144
145	distribution_.push_back(std::move(d));	4,980,204 ✔
146	} // incoming energies	4,980,204 ✔
147	}	1,504,104 ✔
148
149	double ContinuousTabular::sample(double E, uint64_t* seed) const	41,113,499 ✔
150	{
151	// Read number of interpolation regions and incoming energies
152	bool histogram_interp;	41,113,499 ✔
153	if (n_region_ == 1) {	41,113,499 !
154	histogram_interp = (interpolation_[0] == Interpolation::histogram);	41,113,499 ✔
155	} else {
156	histogram_interp = false;
157	}
158
159	// Find energy bin and calculate interpolation factor -- if the energy is
160	// outside the range of the tabulated energies, choose the first or last bins
161	int i;	41,113,499 ✔
162	double r;	41,113,499 ✔
163	get_energy_index(energy_, E, i, r);	41,113,499 ✔
164
165	// Sample between the ith and [i+1]th bin
166	int l;	41,113,499 ✔
167	if (histogram_interp) {	41,113,499 ✔
168	l = i;	12,309 ✔
169	} else {
170	l = r > prn(seed) ? i + 1 : i;	41,101,190 ✔
171	}
172
173	// Determine outgoing energy bin
174	int n_energy_out = distribution_[l].e_out.size();	41,113,499 ✔
175	int n_discrete = distribution_[l].n_discrete;	41,113,499 ✔
176	double r1 = prn(seed);	41,113,499 ✔
177	double c_k = distribution_[l].c[0];	41,113,499 ✔
178	int k = 0;	41,113,499 ✔
179	int end = n_energy_out - 2;	41,113,499 ✔
180
181	// Discrete portion
182	for (int j = 0; j < n_discrete; ++j) {	43,092,685 ✔
183	k = j;	2,589,422 ✔
184	c_k = distribution_[l].c[k];	2,589,422 ✔
185	if (r1 < c_k) {	2,589,422 ✔
186	end = j;
187	break;
188	}
189	}
190
191	// Continuous portion
192	double c_k1;	41,113,499 ✔
193	for (int j = n_discrete; j < end; ++j) {	2,147,483,647 ✔
194	k = j;	2,147,483,647 ✔
195	c_k1 = distribution_[l].c[k + 1];	2,147,483,647 ✔
196	if (r1 < c_k1)	2,147,483,647 ✔
197	break;
198	k = j + 1;
199	c_k = c_k1;
200	}
201
202	double E_l_k = distribution_[l].e_out[k];	41,113,499 ✔
203
204	if (k < n_discrete) {	41,113,499 ✔
205	// Discrete case
206	return E_l_k;
207	} else {
208	// Continuous case
209	double p_l_k = distribution_[l].p[k];	40,503,263 ✔
210	double E_out;	40,503,263 ✔
211	if (distribution_[l].interpolation == Interpolation::histogram) {	40,503,263 ✔
212	// Histogram interpolation
213	if (p_l_k > 0.0) {	767,890 !
214	E_out = E_l_k + (r1 - c_k) / p_l_k;	767,890 ✔
215	} else {
216	E_out = E_l_k;
217	}
218	} else if (distribution_[l].interpolation == Interpolation::lin_lin) {	39,735,373 !
219	// Linear-linear interpolation
220	double E_l_k1 = distribution_[l].e_out[k + 1];	39,735,373 !
221	double p_l_k1 = distribution_[l].p[k + 1];	39,735,373 ✔
222
223	if (E_l_k != E_l_k1) {	39,735,373 !
224	double frac = (p_l_k1 - p_l_k) / (E_l_k1 - E_l_k);	39,735,373 ✔
225	if (frac == 0.0) {	39,735,373 !
UNCOV 226	E_out = E_l_k + (r1 - c_k) / p_l_k;	×
227	} else {
228	E_out =	79,470,746 ✔
229	E_l_k +
230	(std::sqrt(std::max(0.0, p_l_k * p_l_k + 2.0 * frac * (r1 - c_k))) -	79,470,746 !
231	p_l_k) /	39,735,373 ✔
232	frac;
233	}
234	} else {
235	E_out = E_l_k;
236	}
237	} else {
UNCOV 238	throw std::runtime_error {	×
239	"Unexpected interpolation for continuous energy "
UNCOV 240	"distribution."};	×
241	}
242
243	// Now interpolate between incident energy bins i and i + 1
244	if (!histogram_interp && n_energy_out > 1) {	40,503,263 ✔
245	// Interpolation for energy E1 and EK
246	n_energy_out = distribution_[i].e_out.size();	40,490,954 ✔
247	n_discrete = distribution_[i].n_discrete;	40,490,954 ✔
248	const double E_i_1 = distribution_[i].e_out[n_discrete];	40,490,954 ✔
249	const double E_i_K = distribution_[i].e_out[n_energy_out - 1];	40,490,954 ✔
250
251	n_energy_out = distribution_[i + 1].e_out.size();	40,490,954 ✔
252	n_discrete = distribution_[i + 1].n_discrete;	40,490,954 ✔
253	const double E_i1_1 = distribution_[i + 1].e_out[n_discrete];	40,490,954 ✔
254	const double E_i1_K = distribution_[i + 1].e_out[n_energy_out - 1];	40,490,954 ✔
255
256	const double E_1 = E_i_1 + r * (E_i1_1 - E_i_1);	40,490,954 ✔
257	const double E_K = E_i_K + r * (E_i1_K - E_i_K);	40,490,954 ✔
258
259	if (l == i) {	40,490,954 ✔
260	return E_1 + (E_out - E_i_1) * (E_K - E_1) / (E_i_K - E_i_1);	29,075,914 ✔
261	} else {
262	return E_1 + (E_out - E_i1_1) * (E_K - E_1) / (E_i1_K - E_i1_1);	11,415,040 ✔
263	}
264	} else {
265	return E_out;
266	}
267	}
268	}
269
270	//==============================================================================
271	// MaxwellEnergy implementation
272	//==============================================================================
273
274	MaxwellEnergy::MaxwellEnergy(hid_t group)	2,709 ✔
275	{
276	read_attribute(group, "u", u_);	2,709 ✔
277	hid_t dset = open_dataset(group, "theta");	2,709 ✔
278	theta_ = Tabulated1D {dset};	2,709 ✔
279	close_dataset(dset);	2,709 ✔
280	}	2,709 ✔
281
282	double MaxwellEnergy::sample(double E, uint64_t* seed) const	143,803 ✔
283	{
284	// Get temperature corresponding to incoming energy
285	double theta = theta_(E);	143,803 ✔
286
287	while (true) {	143,803 ✔
288	// Sample maxwell fission spectrum
289	double E_out = maxwell_spectrum(theta, seed);	143,803 ✔
290
291	// Accept energy based on restriction energy
292	if (E_out <= E - u_)	143,803 !
293	return E_out;	143,803 ✔
294	}
295	}
296
297	//==============================================================================
298	// Evaporation implementation
299	//==============================================================================
300
301	Evaporation::Evaporation(hid_t group)	5,653 ✔
302	{
303	read_attribute(group, "u", u_);	5,653 ✔
304	hid_t dset = open_dataset(group, "theta");	5,653 ✔
305	theta_ = Tabulated1D {dset};	5,653 ✔
306	close_dataset(dset);	5,653 ✔
307	}	5,653 ✔
308
309	double Evaporation::sample(double E, uint64_t* seed) const	29,139 ✔
310	{
311	// Get temperature corresponding to incoming energy
312	double theta = theta_(E);	29,139 ✔
313
314	double y = (E - u_) / theta;	29,139 ✔
315	double v = 1.0 - std::exp(-y);	29,139 ✔
316
317	// Sample outgoing energy based on evaporation spectrum probability
318	// density function
319	double x;	33,990 ✔
320	while (true) {	33,990 ✔
321	x = -std::log((1.0 - v * prn(seed)) * (1.0 - v * prn(seed)));	33,990 ✔
322	if (x <= y)	33,990 ✔
323	break;
324	}
325
326	return x * theta;	29,139 ✔
327	}
328
329	//==============================================================================
330	// WattEnergy implementation
331	//==============================================================================
332
333	WattEnergy::WattEnergy(hid_t group)	195 ✔
334	{
335	// Read restriction energy
336	read_attribute(group, "u", u_);	195 ✔
337
338	// Read tabulated functions
339	hid_t dset = open_dataset(group, "a");	195 ✔
340	a_ = Tabulated1D {dset};	195 ✔
341	close_dataset(dset);	195 ✔
342	dset = open_dataset(group, "b");	195 ✔
343	b_ = Tabulated1D {dset};	195 ✔
344	close_dataset(dset);	195 ✔
345	}	195 ✔
346
347	double WattEnergy::sample(double E, uint64_t* seed) const	1,831,599 ✔
348	{
349	// Determine Watt parameters at incident energy
350	double a = a_(E);	1,831,599 ✔
351	double b = b_(E);	1,831,599 ✔
352
353	while (true) {	1,831,599 ✔
354	// Sample energy-dependent Watt fission spectrum
355	double E_out = watt_spectrum(a, b, seed);	1,831,599 ✔
356
357	// Accept energy based on restriction energy
358	if (E_out <= E - u_)	1,831,599 !
359	return E_out;	1,831,599 ✔
360	}
361	}
362
363	} // namespace openmc

openmc-dev / openmc / 34416149688

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