22110756672

Committed 17 Feb 2026 06:31PM UTC coverage: 81.831% (+0.1%) from 81.721%

Build # 22110756672

Build Type

Pull #3813

github

Committed by

web-flow

Commit Message

Merge 47dc7194f into 977ade79a

Pull Request Pull Request #3813: Check for positive radii

Run Details

17268 of 24298 branches covered (71.07%)

Branch coverage included in aggregate %.

16 of 16 new or added lines in 1 file covered. (100.0%)

1065 existing lines in 28 files now uncovered.

57520 of 67095 relevant lines covered (85.73%)

45595325.95 hits per line

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

77.36

/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;
  }
}

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,233,788✔
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,233,788!
28	case 1:	10,455,423✔
29	case 11:
30	case 21:
31	return Interpolation::histogram;	10,455,423✔
32	case 2:	42,775,021✔
33	case 12:
34	case 22:
35	return Interpolation::lin_lin;	42,775,021✔
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)	2,147,483,647✔
54	{
55	return mt == N_FISSION \|\| mt == N_F \|\| mt == N_NF \|\| mt == N_2NF \|\|	2,147,483,647✔
56	mt == N_3NF;	2,147,483,647✔
57	}
58
59	bool is_disappearance(int mt)	1,123,522✔
60	{
61	if (mt >= N_DISAPPEAR && mt <= N_DA) {	1,123,522✔
62	return true;	156,990✔
63	} else if (mt >= N_P0 && mt <= N_AC) {	966,532!
64	return true;	328,038✔
65	} else if (mt == N_TA \|\| mt == N_DT \|\| mt == N_P3HE \|\| mt == N_D3HE \|\|	638,494!
66	mt == N_3HEA \|\| mt == N_3P) {	638,494!
UNCOV 67	return true;	×
68	} else {
69	return false;	638,494✔
70	}
71	}
72
73	bool is_inelastic_scatter(int mt)	1,192,031✔
74	{
75	if (mt < 100) {	1,192,031✔
76	if (is_fission(mt)) {	655,180✔
77	return false;	14,453✔
78	} else {
79	return mt >= MISC && mt != 27;	640,727!
80	}
81	} else if (mt <= 200) {	536,851✔
82	return !is_disappearance(mt);	91,130✔
83	} else if (mt >= N_2N0 && mt <= N_2NC) {	445,721!
UNCOV 84	return true;	×
85	} else {
86	return false;	445,721✔
87	}
88	}
89
90	unique_ptr<Function1D> read_function(hid_t group, const char* name)	3,118,308✔
91	{
92	hid_t obj_id = open_object(group, name);	3,118,308✔
93	std::string func_type;	3,118,308✔
94	read_attribute(obj_id, "type", func_type);	3,118,308✔
95	unique_ptr<Function1D> func;	3,118,308✔
96	if (func_type == "Tabulated1D") {	3,118,308✔
97	func = make_unique<Tabulated1D>(obj_id);	2,464,142✔
98	} else if (func_type == "Polynomial") {	654,166✔
99	func = make_unique<Polynomial>(obj_id);	654,002✔
100	} else if (func_type == "CoherentElastic") {	164✔
101	func = make_unique<CoherentElasticXS>(obj_id);	124✔
102	} else if (func_type == "IncoherentElastic") {	40!
UNCOV 103	func = make_unique<IncoherentElasticXS>(obj_id);	×
104	} else if (func_type == "Sum") {	40!
105	func = make_unique<Sum1D>(obj_id);	40✔
106	} else {
UNCOV 107	throw std::runtime_error {"Unknown function type " + func_type +	×
108	" for dataset " + object_name(obj_id)};	×
109	}
110	close_object(obj_id);	3,118,308✔
111	return func;	6,236,616✔
112	}	3,118,308✔
113
114	//==============================================================================
115	// Polynomial implementation
116	//==============================================================================
117
118	Polynomial::Polynomial(hid_t dset)	654,002✔
119	{
120	// Read coefficients into a vector
121	read_dataset(dset, coef_);	654,002✔
122	}	654,002✔
123
124	double Polynomial::operator()(double x) const	1,141,380,135✔
125	{
126	// Use Horner's rule to evaluate polynomial. Note that coefficients are
127	// ordered in increasing powers of x.
128	double y = 0.0;	1,141,380,135✔
129	for (auto c = coef_.crbegin(); c != coef_.crend(); ++c) {	2,147,483,647✔
130	y = y * x + *c;	2,147,483,647✔
131	}
132	return y;	1,141,380,135✔
133	}
134
135	//==============================================================================
136	// Tabulated1D implementation
137	//==============================================================================
138
139	Tabulated1D::Tabulated1D(hid_t dset)	2,473,962✔
140	{
141	read_attribute(dset, "breakpoints", nbt_);	2,473,962✔
142	n_regions_ = nbt_.size();	2,473,962✔
143
144	// Change 1-indexing to 0-indexing
145	for (auto& b : nbt_)	4,950,260✔
146	--b;	2,476,298✔
147
148	vector<int> int_temp;	2,473,962✔
149	read_attribute(dset, "interpolation", int_temp);	2,473,962✔
150
151	// Convert vector of ints into Interpolation
152	for (const auto i : int_temp)	4,950,260✔
153	int_.push_back(int2interp(i));	2,476,298✔
154
155	tensor::Tensor<double> arr;	2,473,962✔
156	read_dataset(dset, arr);	2,473,962✔
157
158	tensor::View<double> xs = arr.slice(0);	2,473,962✔
159	tensor::View<double> ys = arr.slice(1);	2,473,962✔
160
161	std::copy(xs.begin(), xs.end(), std::back_inserter(x_));	2,473,962✔
162	std::copy(ys.begin(), ys.end(), std::back_inserter(y_));	2,473,962✔
163	n_pairs_ = x_.size();	2,473,962✔
164	}	2,473,962✔
165
166	double Tabulated1D::operator()(double x) const	2,147,483,647✔
167	{
168	// find which bin the abscissa is in -- if the abscissa is outside the
169	// tabulated range, the first or last point is chosen, i.e. no interpolation
170	// is done outside the energy range
171	int i;
172	if (x < x_[0]) {	2,147,483,647✔
173	return y_[0];	1,621,757,048✔
174	} else if (x > x_[n_pairs_ - 1]) {	2,147,483,647✔
175	return y_[n_pairs_ - 1];	37,995,916✔
176	} else {
177	i = lower_bound_index(x_.begin(), x_.end(), x);	2,147,483,647✔
178	}
179
180	// determine interpolation scheme
181	Interpolation interp;
182	if (n_regions_ == 0) {	2,147,483,647!
UNCOV 183	interp = Interpolation::lin_lin;	×
184	} else {
185	interp = int_[0];	2,147,483,647✔
186	for (int j = 0; j < n_regions_; ++j) {	2,147,483,647!
187	if (i < nbt_[j]) {	2,147,483,647✔
188	interp = int_[j];	2,147,483,647✔
189	break;	2,147,483,647✔
190	}
191	}
192	}
193
194	// handle special case of histogram interpolation
195	if (interp == Interpolation::histogram)	2,147,483,647✔
196	return y_[i];	2,147,483,647✔
197
198	// determine bounding values
199	double x0 = x_[i];	2,147,483,647✔
200	double x1 = x_[i + 1];	2,147,483,647✔
201	double y0 = y_[i];	2,147,483,647✔
202	double y1 = y_[i + 1];	2,147,483,647✔
203
204	// determine interpolation factor and interpolated value
205	double r;
206	switch (interp) {	2,147,483,647!
207	case Interpolation::lin_lin:	2,147,483,647✔
208	r = (x - x0) / (x1 - x0);	2,147,483,647✔
209	return y0 + r * (y1 - y0);	2,147,483,647✔
UNCOV 210	case Interpolation::lin_log:	×
211	r = log(x / x0) / log(x1 / x0);	×
212	return y0 + r * (y1 - y0);	×
213	case Interpolation::log_lin:	×
214	r = (x - x0) / (x1 - x0);	×
215	return y0 * exp(r * log(y1 / y0));	×
216	case Interpolation::log_log:	21,280,626✔
217	r = log(x / x0) / log(x1 / x0);	21,280,626✔
218	return y0 * exp(r * log(y1 / y0));	21,280,626✔
UNCOV 219	default:	×
220	throw std::runtime_error {"Invalid interpolation scheme."};	×
221	}
222	}
223
224	//==============================================================================
225	// CoherentElasticXS implementation
226	//==============================================================================
227
228	CoherentElasticXS::CoherentElasticXS(hid_t dset)	124✔
229	{
230	// Read 2D array from dataset
231	tensor::Tensor<double> arr;	124✔
232	read_dataset(dset, arr);	124✔
233
234	// Get views for Bragg edges and structure factors
235	tensor::View<double> E = arr.slice(0);	124✔
236	tensor::View<double> s = arr.slice(1);	124✔
237
238	// Copy Bragg edges and partial sums of structure factors
239	std::copy(E.begin(), E.end(), std::back_inserter(bragg_edges_));	124✔
240	std::copy(s.begin(), s.end(), std::back_inserter(factors_));	124✔
241	}	124✔
242
243	double CoherentElasticXS::operator()(double E) const	4,490,870✔
244	{
245	if (E < bragg_edges_[0]) {	4,490,870✔
246	// If energy is below that of the lowest Bragg peak, the elastic cross
247	// section will be zero
248	return 0.0;	2,820✔
249	} else {
250	auto i_grid =
251	lower_bound_index(bragg_edges_.begin(), bragg_edges_.end(), E);	4,488,050✔
252	return factors_[i_grid] / E;	4,488,050✔
253	}
254	}
255
256	//==============================================================================
257	// IncoherentElasticXS implementation
258	//==============================================================================
259
UNCOV 260	IncoherentElasticXS::IncoherentElasticXS(hid_t dset)	×
261	{
262	array<double, 2> tmp;
UNCOV 263	read_dataset(dset, nullptr, tmp);	×
264	bound_xs_ = tmp[0];	×
265	debye_waller_ = tmp[1];	×
266	}	×
267
UNCOV 268	double IncoherentElasticXS::operator()(double E) const	×
269	{
270	// Determine cross section using ENDF-102, Eq. (7.5)
UNCOV 271	double W = debye_waller_;	×
272	return bound_xs_ / 2.0 * ((1 - std::exp(-4.0 * E * W)) / (2.0 * E * W));	×
273	}
274
275	//==============================================================================
276	// Sum1D implementation
277	//==============================================================================
278
279	Sum1D::Sum1D(hid_t group)	40✔
280	{
281	// Get number of functions
282	int n;
283	read_attribute(group, "n", n);	40✔
284
285	// Get each function
286	for (int i = 0; i < n; ++i) {	120✔
287	auto dset_name = fmt::format("func_{}", i + 1);	144✔
288	functions_.push_back(read_function(group, dset_name.c_str()));	80✔
289	}	80✔
290	}	40✔
291
292	double Sum1D::operator()(double x) const	3,681,200✔
293	{
294	double result = 0.0;	3,681,200✔
295	for (auto& func : functions_) {	11,043,600✔
296	result += (*func)(x);	7,362,400✔
297	}
298	return result;	3,681,200✔
299	}
300
301	} // namespace openmc

openmc-dev / openmc / 22110756672

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