21139931957

Committed 19 Jan 2026 01:50PM UTC coverage: 81.835% (-0.2%) from 82.058%

Build # 21139931957

Build Type

Pull #2693

github

Committed by

web-flow

Commit Message

Merge 1258e263b into 5847b0de2

Pull Request Pull Request #2693: Add reactivity control to coupled transport-depletion analyses

Run Details

17180 of 23968 branches covered (71.68%)

Branch coverage included in aggregate %.

78 of 84 new or added lines in 4 files covered. (92.86%)

345 existing lines in 27 files now uncovered.

55616 of 64987 relevant lines covered (85.58%)

50201535.92 hits per line

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

90.98

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

#include "openmc/endf.h"
#include "openmc/hdf5_interface.h"
#include "openmc/math_functions.h"
#include "openmc/random_dist.h"
#include "openmc/random_lcg.h"
#include "openmc/search.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
  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, "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);

  xt::xarray<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 = xt::view(eout, 0, xt::range(j, j + n));
    d.p = xt::view(eout, 1, 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 (true) {
      d.c = xt::view(eout, 2, xt::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
  auto n_energy_in = energy_.size();
  int i;
  double r;
  if (E < energy_[0]) {
    i = 0;
    r = 0.0;
  } else if (E > energy_[n_energy_in - 1]) {
    i = n_energy_in - 2;
    r = 1.0;
  } else {
    i = lower_bound_index(energy_.begin(), energy_.end(), E);
    r = (E - energy_[i]) / (energy_[i + 1] - energy_[i]);
  }

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

openmc-dev / openmc / 21139931957

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