30880334341

Committed 04 Aug 2026 05:19AM UTC coverage: 81.393% (-0.03%) from 81.425%

Build # 30880334341

Build Type

Pull #3675

github

Committed by

web-flow

Commit Message

Merge 553ee9e20 into 8202ef6fb

Pull Request Pull Request #3675: Extend level scattering to support incident photons

Coverage Stats

18534 of 26869 branches covered (68.98%)

Branch coverage included in aggregate %.

55 of 77 new or added lines in 4 files covered. (71.43%)

60286 of 69970 relevant lines covered (86.16%)

49998368.1 hits per line

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

82.7

/src/distribution_energy.cpp

#include "openmc/distribution_energy.h"

#include <algorithm> // for max, min, copy, move
#include <cmath>     // for sqrt, abs
#include <cstddef>   // for size_t
#include <iterator>  // for back_inserter

#include "openmc/tensor.h"

#include "openmc/constants.h"
#include "openmc/endf.h"
#include "openmc/hdf5_interface.h"
#include "openmc/math_functions.h"
#include "openmc/particle.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)
{
  // for backwards compatibility:
  if (attribute_exists(group, "mass_ratio")) {
    read_attribute(group, "threshold", b_);
    read_attribute(group, "mass_ratio", a_);
    c_ = 0.0;
  } else {
    double A, Q;
    std::string temp;
    read_attribute(group, "mass", A);
    read_attribute(group, "q_value", Q);
    read_attribute(group, "particle", temp);
    auto type = ParticleType(temp);
    if (type.is_neutron()) {
      a_ = (A / (A + 1.0)) * (A / (A + 1.0));
      b_ = (A + 1.0) / A * std::abs(Q);
      c_ = 0.0;
    } else if (type.is_photon()) {
      a_ = (A - 1.0) / A;
      b_ = std::abs(Q);
      c_ = 1.0 / (2.0 * MASS_NEUTRON_EV * (A - 1.0));
    } else {
      fatal_error("Unrecognized particle: " + type.str());
    }
  }
}

double LevelInelastic::sample(double E, uint64_t* seed) const
{
  return a_ * (E - b_ - c_ * (E * E));
}

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

openmc-dev / openmc / 30880334341

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