22111819275

Committed 17 Feb 2026 07:03PM UTC coverage: 81.802% (+0.08%) from 81.721%

Build # 22111819275

Build Type

Pull #3809

github

Committed by

web-flow

Commit Message

Merge f0a3f78bd into 977ade79a

Pull Request Pull Request #3809: Implement tally filter for filtering by reaction

Run Details

17321 of 24402 branches covered (70.98%)

Branch coverage included in aggregate %.

97 of 119 new or added lines in 8 files covered. (81.51%)

1361 existing lines in 48 files now uncovered.

57642 of 67238 relevant lines covered (85.73%)

52441039.54 hits per line

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

66.34

/src/endf.cpp

#include "openmc/endf.h"

#include <algorithm> // for copy
#include <cmath>     // for log, exp
#include <iterator>  // for back_inserter
#include <stdexcept> // for runtime_error

#include "openmc/tensor.h"

#include "openmc/array.h"
#include "openmc/constants.h"
#include "openmc/hdf5_interface.h"
#include "openmc/search.h"

namespace openmc {

//==============================================================================
// Functions
//==============================================================================

Interpolation int2interp(int i)
{
  // TODO: We are ignoring specification of two-dimensional interpolation
  // schemes (method of corresponding points and unit base interpolation). Those
  // should be accounted for in the distribution classes somehow.

  switch (i) {
  case 1:
  case 11:
  case 21:
    return Interpolation::histogram;
  case 2:
  case 12:
  case 22:
    return Interpolation::lin_lin;
  case 3:
  case 13:
  case 23:
    return Interpolation::lin_log;
  case 4:
  case 14:
  case 24:
    return Interpolation::log_lin;
  case 5:
  case 15:
  case 25:
    return Interpolation::log_log;
  default:
    throw std::runtime_error {"Invalid interpolation code."};
  }
}

bool is_fission(int mt)
{
  return mt == N_FISSION || mt == N_F || mt == N_NF || mt == N_2NF ||
         mt == N_3NF;
}

bool is_disappearance(int mt)
{
  if (mt >= N_DISAPPEAR && mt <= N_DA) {
    return true;
  } else if (mt >= N_P0 && mt <= N_AC) {
    return true;
  } else if (mt == N_TA || mt == N_DT || mt == N_P3HE || mt == N_D3HE ||
             mt == N_3HEA || mt == N_3P) {
    return true;
  } else {
    return false;
  }
}

bool is_inelastic_scatter(int mt)
{
  if (mt < 100) {
    if (is_fission(mt)) {
      return false;
    } else {
      return mt >= MISC && mt != 27;
    }
  } else if (mt <= 200) {
    return !is_disappearance(mt);
  } else if (mt >= N_2N0 && mt <= N_2NC) {
    return true;
  } else {
    return false;
  }
}

bool mt_matches(int event_mt, int target_mt)
{
  // Direct match
  if (event_mt == target_mt)
    return true;

  // Check if event_mt is a component of target_mt summation reaction
  switch (target_mt) {
  case TOTAL_XS:
  case SCORE_TOTAL:
    return event_mt == ELASTIC || mt_matches(event_mt, N_NONELASTIC);

  case N_NONELASTIC: {
    static constexpr int components[] = {4, 5, 11, 16, 17, 22, 23, 24, 25, 27,
      28, 29, 30, 32, 33, 34, 35, 36, 37, 41, 42, 44, 45, 152, 153, 154, 156,
      157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171,
      172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 183, 184, 185, 186, 187,
      188, 189, 190, 194, 195, 196, 198, 199, 200};
    for (int mt : components) {
      if (mt_matches(event_mt, mt))
        return true;
    }
    return false;
  }

  case N_LEVEL:
    // Inelastic scattering levels
    return event_mt >= 50 && event_mt <= N_NC;

  case N_2N:
    // (n,2n) to excited states
    return event_mt >= N_2N0 && event_mt <= N_2NC;

  case N_FISSION:
    return is_fission(event_mt);

  case 27:
    return is_fission(event_mt) || is_disappearance(event_mt);

  case N_DISAPPEAR: {
    return is_disappearance(event_mt);
  }

  case N_P:
    // (n,p) to excited states
    return event_mt >= N_P0 && event_mt <= N_PC;

  case N_D:
    // (n,d) to excited states
    return event_mt >= N_D0 && event_mt <= N_DC;

  case N_T:
    // (n,t) to excited states
    return event_mt >= N_T0 && event_mt <= N_TC;

  case N_3HE:
    // (n,3He) to excited states
    return event_mt >= N_3HE0 && event_mt <= N_3HEC;

  case N_A:
    // (n,alpha) to excited states
    return event_mt >= N_A0 && event_mt <= N_AC;

  case 501:
    return event_mt == 502 || event_mt == 504 || mt_matches(event_mt, 516) ||
           mt_matches(event_mt, 522);

  case PAIR_PROD:
    return event_mt == PAIR_PROD_ELEC || event_mt == PAIR_PROD_NUC;

  case PHOTOELECTRIC:
    return event_mt >= 534 && event_mt < 573;

  default:
    return false;
  }
}

unique_ptr<Function1D> read_function(hid_t group, const char* name)
{
  hid_t obj_id = open_object(group, name);
  std::string func_type;
  read_attribute(obj_id, "type", func_type);
  unique_ptr<Function1D> func;
  if (func_type == "Tabulated1D") {
    func = make_unique<Tabulated1D>(obj_id);
  } else if (func_type == "Polynomial") {
    func = make_unique<Polynomial>(obj_id);
  } else if (func_type == "CoherentElastic") {
    func = make_unique<CoherentElasticXS>(obj_id);
  } else if (func_type == "IncoherentElastic") {
    func = make_unique<IncoherentElasticXS>(obj_id);
  } else if (func_type == "Sum") {
    func = make_unique<Sum1D>(obj_id);
  } else {
    throw std::runtime_error {"Unknown function type " + func_type +
                              " for dataset " + object_name(obj_id)};
  }
  close_object(obj_id);
  return func;
}

//==============================================================================
// Polynomial implementation
//==============================================================================

Polynomial::Polynomial(hid_t dset)
{
  // Read coefficients into a vector
  read_dataset(dset, coef_);
}

double Polynomial::operator()(double x) const
{
  // Use Horner's rule to evaluate polynomial. Note that coefficients are
  // ordered in increasing powers of x.
  double y = 0.0;
  for (auto c = coef_.crbegin(); c != coef_.crend(); ++c) {
    y = y * x + *c;
  }
  return y;
}

//==============================================================================
// Tabulated1D implementation
//==============================================================================

Tabulated1D::Tabulated1D(hid_t dset)
{
  read_attribute(dset, "breakpoints", nbt_);
  n_regions_ = nbt_.size();

  // Change 1-indexing to 0-indexing
  for (auto& b : nbt_)
    --b;

  vector<int> int_temp;
  read_attribute(dset, "interpolation", int_temp);

  // Convert vector of ints into Interpolation
  for (const auto i : int_temp)
    int_.push_back(int2interp(i));

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

  tensor::View<double> xs = arr.slice(0);
  tensor::View<double> ys = arr.slice(1);

  std::copy(xs.begin(), xs.end(), std::back_inserter(x_));
  std::copy(ys.begin(), ys.end(), std::back_inserter(y_));
  n_pairs_ = x_.size();
}

double Tabulated1D::operator()(double x) const
{
  // find which bin the abscissa is in -- if the abscissa is outside the
  // tabulated range, the first or last point is chosen, i.e. no interpolation
  // is done outside the energy range
  int i;
  if (x < x_[0]) {
    return y_[0];
  } else if (x > x_[n_pairs_ - 1]) {
    return y_[n_pairs_ - 1];
  } else {
    i = lower_bound_index(x_.begin(), x_.end(), x);
  }

  // determine interpolation scheme
  Interpolation interp;
  if (n_regions_ == 0) {
    interp = Interpolation::lin_lin;
  } else {
    interp = int_[0];
    for (int j = 0; j < n_regions_; ++j) {
      if (i < nbt_[j]) {
        interp = int_[j];
        break;
      }
    }
  }

  // handle special case of histogram interpolation
  if (interp == Interpolation::histogram)
    return y_[i];

  // determine bounding values
  double x0 = x_[i];
  double x1 = x_[i + 1];
  double y0 = y_[i];
  double y1 = y_[i + 1];

  // determine interpolation factor and interpolated value
  double r;
  switch (interp) {
  case Interpolation::lin_lin:
    r = (x - x0) / (x1 - x0);
    return y0 + r * (y1 - y0);
  case Interpolation::lin_log:
    r = log(x / x0) / log(x1 / x0);
    return y0 + r * (y1 - y0);
  case Interpolation::log_lin:
    r = (x - x0) / (x1 - x0);
    return y0 * exp(r * log(y1 / y0));
  case Interpolation::log_log:
    r = log(x / x0) / log(x1 / x0);
    return y0 * exp(r * log(y1 / y0));
  default:
    throw std::runtime_error {"Invalid interpolation scheme."};
  }
}

//==============================================================================
// CoherentElasticXS implementation
//==============================================================================

CoherentElasticXS::CoherentElasticXS(hid_t dset)
{
  // Read 2D array from dataset
  tensor::Tensor<double> arr;
  read_dataset(dset, arr);

  // Get views for Bragg edges and structure factors
  tensor::View<double> E = arr.slice(0);
  tensor::View<double> s = arr.slice(1);

  // Copy Bragg edges and partial sums of structure factors
  std::copy(E.begin(), E.end(), std::back_inserter(bragg_edges_));
  std::copy(s.begin(), s.end(), std::back_inserter(factors_));
}

double CoherentElasticXS::operator()(double E) const
{
  if (E < bragg_edges_[0]) {
    // If energy is below that of the lowest Bragg peak, the elastic cross
    // section will be zero
    return 0.0;
  } else {
    auto i_grid =
      lower_bound_index(bragg_edges_.begin(), bragg_edges_.end(), E);
    return factors_[i_grid] / E;
  }
}

//==============================================================================
// IncoherentElasticXS implementation
//==============================================================================

IncoherentElasticXS::IncoherentElasticXS(hid_t dset)
{
  array<double, 2> tmp;
  read_dataset(dset, nullptr, tmp);
  bound_xs_ = tmp[0];
  debye_waller_ = tmp[1];
}

double IncoherentElasticXS::operator()(double E) const
{
  // Determine cross section using ENDF-102, Eq. (7.5)
  double W = debye_waller_;
  return bound_xs_ / 2.0 * ((1 - std::exp(-4.0 * E * W)) / (2.0 * E * W));
}

//==============================================================================
// Sum1D implementation
//==============================================================================

Sum1D::Sum1D(hid_t group)
{
  // Get number of functions
  int n;
  read_attribute(group, "n", n);

  // Get each function
  for (int i = 0; i < n; ++i) {
    auto dset_name = fmt::format("func_{}", i + 1);
    functions_.push_back(read_function(group, dset_name.c_str()));
  }
}

double Sum1D::operator()(double x) const
{
  double result = 0.0;
  for (auto& func : functions_) {
    result += (*func)(x);
  }
  return result;
}

} // namespace openmc

1	#include "openmc/endf.h"
2
3	#include <algorithm> // for copy
4	#include <cmath> // for log, exp
5	#include <iterator> // for back_inserter
6	#include <stdexcept> // for runtime_error
7
8	#include "openmc/tensor.h"
9
10	#include "openmc/array.h"
11	#include "openmc/constants.h"
12	#include "openmc/hdf5_interface.h"
13	#include "openmc/search.h"
14
15	namespace openmc {
16
17	//==============================================================================
18	// Functions
19	//==============================================================================
20
21	Interpolation int2interp(int i)	53,248,474✔
22	{
23	// TODO: We are ignoring specification of two-dimensional interpolation
24	// schemes (method of corresponding points and unit base interpolation). Those
25	// should be accounted for in the distribution classes somehow.
26
27	switch (i) {	53,248,474!
28	case 1:	10,457,873✔
29	case 11:
30	case 21:
31	return Interpolation::histogram;	10,457,873✔
32	case 2:	42,787,257✔
33	case 12:
34	case 22:
35	return Interpolation::lin_lin;	42,787,257✔
UNCOV 36	case 3:	×
37	case 13:
38	case 23:
UNCOV 39	return Interpolation::lin_log;	×
40	case 4:	×
41	case 14:
42	case 24:
UNCOV 43	return Interpolation::log_lin;	×
44	case 5:	3,344✔
45	case 15:
46	case 25:
47	return Interpolation::log_log;	3,344✔
UNCOV 48	default:	×
49	throw std::runtime_error {"Invalid interpolation code."};	×
50	}
51	}
52
53	bool is_fission(int mt)	3,569,304,828✔
54	{
55	return mt == N_FISSION \|\| mt == N_F \|\| mt == N_NF \|\| mt == N_2NF \|\|	3,569,304,828✔
56	mt == N_3NF;	3,569,304,828✔
57	}
58
59	bool is_disappearance(int mt)	1,162,994✔
60	{
61	if (mt >= N_DISAPPEAR && mt <= N_DA) {	1,162,994✔
62	return true;	195,818✔
63	} else if (mt >= N_P0 && mt <= N_AC) {	967,176!
64	return true;	328,038✔
65	} else if (mt == N_TA \|\| mt == N_DT \|\| mt == N_P3HE \|\| mt == N_D3HE \|\|	639,138!
66	mt == N_3HEA \|\| mt == N_3P) {	639,138!
UNCOV 67	return true;	×
68	} else {
69	return false;	639,138✔
70	}
71	}
72
73	bool is_inelastic_scatter(int mt)	1,192,745✔
74	{
75	if (mt < 100) {	1,192,745✔
76	if (is_fission(mt)) {	655,852✔
77	return false;	14,467✔
78	} else {
79	return mt >= MISC && mt != 27;	641,385!
80	}
81	} else if (mt <= 200) {	536,893✔
82	return !is_disappearance(mt);	91,144✔
83	} else if (mt >= N_2N0 && mt <= N_2NC) {	445,749!
UNCOV 84	return true;	×
85	} else {
86	return false;	445,749✔
87	}
88	}
89
90	bool mt_matches(int event_mt, int target_mt)	2,112,950✔
91	{
92	// Direct match
93	if (event_mt == target_mt)	2,112,950✔
94	return true;	197,440✔
95
96	// Check if event_mt is a component of target_mt summation reaction
97	switch (target_mt) {	1,915,510!
98	case TOTAL_XS:	290,310✔
99	case SCORE_TOTAL:
100	return event_mt == ELASTIC \|\| mt_matches(event_mt, N_NONELASTIC);	290,310!
101
102	case N_NONELASTIC: {	104,680✔
103	static constexpr int components[] = {4, 5, 11, 16, 17, 22, 23, 24, 25, 27,
104	28, 29, 30, 32, 33, 34, 35, 36, 37, 41, 42, 44, 45, 152, 153, 154, 156,
105	157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171,
106	172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 183, 184, 185, 186, 187,
107	188, 189, 190, 194, 195, 196, 198, 199, 200};
108	for (int mt : components) {	556,720!
109	if (mt_matches(event_mt, mt))	556,720✔
110	return true;	104,680✔
111	}
NEW 112	return false;	×
113	}
114
115	case N_LEVEL:	104,680✔
116	// Inelastic scattering levels
117	return event_mt >= 50 && event_mt <= N_NC;	104,680✔
118
119	case N_2N:	340,230✔
120	// (n,2n) to excited states
121	return event_mt >= N_2N0 && event_mt <= N_2NC;	340,230!
122
123	case N_FISSION:	278,960✔
124	return is_fission(event_mt);	278,960✔
125
126	case 27:	50,150✔
127	return is_fission(event_mt) \|\| is_disappearance(event_mt);	50,150!
128
NEW 129	case N_DISAPPEAR: {	×
NEW 130	return is_disappearance(event_mt);	×
131	}
132
NEW 133	case N_P:	×
134	// (n,p) to excited states
NEW 135	return event_mt >= N_P0 && event_mt <= N_PC;	×
136
NEW 137	case N_D:	×
138	// (n,d) to excited states
NEW 139	return event_mt >= N_D0 && event_mt <= N_DC;	×
140
NEW 141	case N_T:	×
142	// (n,t) to excited states
NEW 143	return event_mt >= N_T0 && event_mt <= N_TC;	×
144
NEW 145	case N_3HE:	×
146	// (n,3He) to excited states
NEW 147	return event_mt >= N_3HE0 && event_mt <= N_3HEC;	×
148
NEW 149	case N_A:	×
150	// (n,alpha) to excited states
NEW 151	return event_mt >= N_A0 && event_mt <= N_AC;	×
152
NEW 153	case 501:	×
NEW 154	return event_mt == 502 \|\| event_mt == 504 \|\| mt_matches(event_mt, 516) \|\|	×
NEW 155	mt_matches(event_mt, 522);	×
156
NEW 157	case PAIR_PROD:	×
NEW 158	return event_mt == PAIR_PROD_ELEC \|\| event_mt == PAIR_PROD_NUC;	×
159
NEW 160	case PHOTOELECTRIC:	×
NEW 161	return event_mt >= 534 && event_mt < 573;	×
162
163	default:	746,500✔
164	return false;	746,500✔
165	}
166	}
167
168	unique_ptr<Function1D> read_function(hid_t group, const char* name)	3,119,596✔
169	{
170	hid_t obj_id = open_object(group, name);	3,119,596✔
171	std::string func_type;	3,119,596✔
172	read_attribute(obj_id, "type", func_type);	3,119,596✔
173	unique_ptr<Function1D> func;	3,119,596✔
174	if (func_type == "Tabulated1D") {	3,119,596✔
175	func = make_unique<Tabulated1D>(obj_id);	2,464,716✔
176	} else if (func_type == "Polynomial") {	654,880✔
177	func = make_unique<Polynomial>(obj_id);	654,716✔
178	} else if (func_type == "CoherentElastic") {	164✔
179	func = make_unique<CoherentElasticXS>(obj_id);	124✔
180	} else if (func_type == "IncoherentElastic") {	40!
UNCOV 181	func = make_unique<IncoherentElasticXS>(obj_id);	×
182	} else if (func_type == "Sum") {	40!
183	func = make_unique<Sum1D>(obj_id);	40✔
184	} else {
UNCOV 185	throw std::runtime_error {"Unknown function type " + func_type +	×
186	" for dataset " + object_name(obj_id)};	×
187	}
188	close_object(obj_id);	3,119,596✔
189	return func;	6,239,192✔
190	}	3,119,596✔
191
192	//==============================================================================
193	// Polynomial implementation
194	//==============================================================================
195
196	Polynomial::Polynomial(hid_t dset)	654,716✔
197	{
198	// Read coefficients into a vector
199	read_dataset(dset, coef_);	654,716✔
200	}	654,716✔
201
202	double Polynomial::operator()(double x) const	1,143,705,328✔
203	{
204	// Use Horner's rule to evaluate polynomial. Note that coefficients are
205	// ordered in increasing powers of x.
206	double y = 0.0;	1,143,705,328✔
207	for (auto c = coef_.crbegin(); c != coef_.crend(); ++c) {	2,707,889,817✔
208	y = y * x + *c;	2,562,522,875✔
209	}
210	return y;	1,143,705,328✔
211	}
212
213	//==============================================================================
214	// Tabulated1D implementation
215	//==============================================================================
216
217	Tabulated1D::Tabulated1D(hid_t dset)	2,474,536✔
218	{
219	read_attribute(dset, "breakpoints", nbt_);	2,474,536✔
220	n_regions_ = nbt_.size();	2,474,536✔
221
222	// Change 1-indexing to 0-indexing
223	for (auto& b : nbt_)	4,951,408✔
224	--b;	2,476,872✔
225
226	vector<int> int_temp;	2,474,536✔
227	read_attribute(dset, "interpolation", int_temp);	2,474,536✔
228
229	// Convert vector of ints into Interpolation
230	for (const auto i : int_temp)	4,951,408✔
231	int_.push_back(int2interp(i));	2,476,872✔
232
233	tensor::Tensor<double> arr;	2,474,536✔
234	read_dataset(dset, arr);	2,474,536✔
235
236	tensor::View<double> xs = arr.slice(0);	2,474,536✔
237	tensor::View<double> ys = arr.slice(1);	2,474,536✔
238
239	std::copy(xs.begin(), xs.end(), std::back_inserter(x_));	2,474,536✔
240	std::copy(ys.begin(), ys.end(), std::back_inserter(y_));	2,474,536✔
241	n_pairs_ = x_.size();	2,474,536✔
242	}	2,474,536✔
243
244	double Tabulated1D::operator()(double x) const	3,702,632,883✔
245	{
246	// find which bin the abscissa is in -- if the abscissa is outside the
247	// tabulated range, the first or last point is chosen, i.e. no interpolation
248	// is done outside the energy range
249	int i;
250	if (x < x_[0]) {	3,702,632,883✔
251	return y_[0];	1,622,922,744✔
252	} else if (x > x_[n_pairs_ - 1]) {	3,500,581,122✔
253	return y_[n_pairs_ - 1];	37,995,916✔
254	} else {
255	i = lower_bound_index(x_.begin(), x_.end(), x);	3,495,734,687✔
256	}
257
258	// determine interpolation scheme
259	Interpolation interp;
260	if (n_regions_ == 0) {	3,495,734,687!
UNCOV 261	interp = Interpolation::lin_lin;	×
262	} else {
263	interp = int_[0];	3,495,734,687✔
264	for (int j = 0; j < n_regions_; ++j) {	3,501,056,548!
265	if (i < nbt_[j]) {	3,501,056,548✔
266	interp = int_[j];	3,495,734,687✔
267	break;	3,495,734,687✔
268	}
269	}
270	}
271
272	// handle special case of histogram interpolation
273	if (interp == Interpolation::histogram)	3,495,734,687✔
274	return y_[i];	3,121,346,226✔
275
276	// determine bounding values
277	double x0 = x_[i];	2,521,872,108✔
278	double x1 = x_[i + 1];	2,521,872,108✔
279	double y0 = y_[i];	2,521,872,108✔
280	double y1 = y_[i + 1];	2,521,872,108✔
281
282	// determine interpolation factor and interpolated value
283	double r;
284	switch (interp) {	2,521,872,108!
285	case Interpolation::lin_lin:	2,519,725,024✔
286	r = (x - x0) / (x1 - x0);	2,519,725,024✔
287	return y0 + r * (y1 - y0);	2,519,725,024✔
UNCOV 288	case Interpolation::lin_log:	×
289	r = log(x / x0) / log(x1 / x0);	×
290	return y0 + r * (y1 - y0);	×
291	case Interpolation::log_lin:	×
292	r = (x - x0) / (x1 - x0);	×
293	return y0 * exp(r * log(y1 / y0));	×
294	case Interpolation::log_log:	21,280,379✔
295	r = log(x / x0) / log(x1 / x0);	21,280,379✔
296	return y0 * exp(r * log(y1 / y0));	21,280,379✔
UNCOV 297	default:	×
298	throw std::runtime_error {"Invalid interpolation scheme."};	×
299	}
300	}
301
302	//==============================================================================
303	// CoherentElasticXS implementation
304	//==============================================================================
305
306	CoherentElasticXS::CoherentElasticXS(hid_t dset)	124✔
307	{
308	// Read 2D array from dataset
309	tensor::Tensor<double> arr;	124✔
310	read_dataset(dset, arr);	124✔
311
312	// Get views for Bragg edges and structure factors
313	tensor::View<double> E = arr.slice(0);	124✔
314	tensor::View<double> s = arr.slice(1);	124✔
315
316	// Copy Bragg edges and partial sums of structure factors
317	std::copy(E.begin(), E.end(), std::back_inserter(bragg_edges_));	124✔
318	std::copy(s.begin(), s.end(), std::back_inserter(factors_));	124✔
319	}	124✔
320
321	double CoherentElasticXS::operator()(double E) const	4,490,870✔
322	{
323	if (E < bragg_edges_[0]) {	4,490,870✔
324	// If energy is below that of the lowest Bragg peak, the elastic cross
325	// section will be zero
326	return 0.0;	2,820✔
327	} else {
328	auto i_grid =
329	lower_bound_index(bragg_edges_.begin(), bragg_edges_.end(), E);	4,488,050✔
330	return factors_[i_grid] / E;	4,488,050✔
331	}
332	}
333
334	//==============================================================================
335	// IncoherentElasticXS implementation
336	//==============================================================================
337
UNCOV 338	IncoherentElasticXS::IncoherentElasticXS(hid_t dset)	×
339	{
340	array<double, 2> tmp;
UNCOV 341	read_dataset(dset, nullptr, tmp);	×
342	bound_xs_ = tmp[0];	×
343	debye_waller_ = tmp[1];	×
344	}	×
345
UNCOV 346	double IncoherentElasticXS::operator()(double E) const	×
347	{
348	// Determine cross section using ENDF-102, Eq. (7.5)
UNCOV 349	double W = debye_waller_;	×
350	return bound_xs_ / 2.0 * ((1 - std::exp(-4.0 * E * W)) / (2.0 * E * W));	×
351	}
352
353	//==============================================================================
354	// Sum1D implementation
355	//==============================================================================
356
357	Sum1D::Sum1D(hid_t group)	40✔
358	{
359	// Get number of functions
360	int n;
361	read_attribute(group, "n", n);	40✔
362
363	// Get each function
364	for (int i = 0; i < n; ++i) {	120✔
365	auto dset_name = fmt::format("func_{}", i + 1);	144✔
366	functions_.push_back(read_function(group, dset_name.c_str()));	80✔
367	}	80✔
368	}	40✔
369
370	double Sum1D::operator()(double x) const	3,681,200✔
371	{
372	double result = 0.0;	3,681,200✔
373	for (auto& func : functions_) {	11,043,600✔
374	result += (*func)(x);	7,362,400✔
375	}
376	return result;	3,681,200✔
377	}
378
379	} // namespace openmc

openmc-dev / openmc / 22111819275

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