18656608910

Committed 20 Oct 2025 03:16PM UTC coverage: 81.844% (-3.4%) from 85.218%

Build # 18656608910

Build Type

Pull #3454

github

Committed by

web-flow

Commit Message

Merge 5eee478a5 into 055ea15a2

Pull Request Pull Request #3454: Adding variance of variance and normality tests for tally statistics

Run Details

16610 of 23118 branches covered (71.85%)

Branch coverage included in aggregate %.

188 of 312 new or added lines in 7 files covered. (60.26%)

2157 existing lines in 77 files now uncovered.

53896 of 63029 relevant lines covered (85.51%)

42554448.35 hits per line

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

76.9

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

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

  xt::xarray<double> arr;
  read_dataset(dset, arr);

  auto xs = xt::view(arr, 0);
  auto ys = xt::view(arr, 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
  xt::xarray<double> arr;
  read_dataset(dset, arr);

  // Get views for Bragg edges and structure factors
  auto E = xt::view(arr, 0);
  auto s = xt::view(arr, 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 "xtensor/xarray.hpp"
9	#include "xtensor/xview.hpp"
10
11	#include "openmc/array.h"
12	#include "openmc/constants.h"
13	#include "openmc/hdf5_interface.h"
14	#include "openmc/search.h"
15
16	namespace openmc {
17
18	//==============================================================================
19	// Functions
20	//==============================================================================
21
22	Interpolation int2interp(int i)	56,851,546✔
23	{
24	// TODO: We are ignoring specification of two-dimensional interpolation
25	// schemes (method of corresponding points and unit base interpolation). Those
26	// should be accounted for in the distribution classes somehow.
27
28	switch (i) {	56,851,546!
29	case 1:	10,679,473✔
30	case 11:
31	case 21:
32	return Interpolation::histogram;	10,679,473✔
33	case 2:	46,169,001✔
34	case 12:
35	case 22:
36	return Interpolation::lin_lin;	46,169,001✔
37	case 3:	×
38	case 13:
39	case 23:
40	return Interpolation::lin_log;	×
41	case 4:	×
42	case 14:
43	case 24:
44	return Interpolation::log_lin;	×
45	case 5:	3,072✔
46	case 15:
47	case 25:
48	return Interpolation::log_log;	3,072✔
49	default:	×
50	throw std::runtime_error {"Invalid interpolation code."};	×
51	}
52	}
53
54	bool is_fission(int mt)	2,147,483,647✔
55	{
56	return mt == N_FISSION \|\| mt == N_F \|\| mt == N_NF \|\| mt == N_2NF \|\|	2,147,483,647✔
57	mt == N_3NF;	2,147,483,647✔
58	}
59
60	bool is_disappearance(int mt)	1,178,139✔
61	{
62	if (mt >= N_DISAPPEAR && mt <= N_DA) {	1,178,139✔
63	return true;	165,281✔
64	} else if (mt >= N_P0 && mt <= N_AC) {	1,012,858!
65	return true;	335,778✔
66	} else if (mt == N_TA \|\| mt == N_DT \|\| mt == N_P3HE \|\| mt == N_D3HE \|\|	677,080!
67	mt == N_3HEA \|\| mt == N_3P) {	677,080!
68	return true;	×
69	} else {
70	return false;	677,080✔
71	}
72	}
73
74	bool is_inelastic_scatter(int mt)	1,249,513✔
75	{
76	if (mt < 100) {	1,249,513✔
77	if (is_fission(mt)) {	694,117✔
78	return false;	14,304✔
79	} else {
80	return mt >= MISC && mt != 27;	679,813!
81	}
82	} else if (mt <= 200) {	555,396✔
83	return !is_disappearance(mt);	95,760✔
84	} else if (mt >= N_2N0 && mt <= N_2NC) {	459,636!
85	return true;	×
86	} else {
87	return false;	459,636✔
88	}
89	}
90
91	unique_ptr<Function1D> read_function(hid_t group, const char* name)	3,286,402✔
92	{
93	hid_t obj_id = open_object(group, name);	3,286,402✔
94	std::string func_type;	3,286,402✔
95	read_attribute(obj_id, "type", func_type);	3,286,402✔
96	unique_ptr<Function1D> func;	3,286,402✔
97	if (func_type == "Tabulated1D") {	3,286,402✔
98	func = make_unique<Tabulated1D>(obj_id);	2,591,937✔
99	} else if (func_type == "Polynomial") {	694,465✔
100	func = make_unique<Polynomial>(obj_id);	694,281✔
101	} else if (func_type == "CoherentElastic") {	184✔
102	func = make_unique<CoherentElasticXS>(obj_id);	140✔
103	} else if (func_type == "IncoherentElastic") {	44!
104	func = make_unique<IncoherentElasticXS>(obj_id);	×
105	} else if (func_type == "Sum") {	44!
106	func = make_unique<Sum1D>(obj_id);	44✔
107	} else {
108	throw std::runtime_error {"Unknown function type " + func_type +	×
109	" for dataset " + object_name(obj_id)};	×
110	}
111	close_object(obj_id);	3,286,402✔
112	return func;	6,572,804✔
113	}	3,286,402✔
114
115	//==============================================================================
116	// Polynomial implementation
117	//==============================================================================
118
119	Polynomial::Polynomial(hid_t dset)	694,281✔
120	{
121	// Read coefficients into a vector
122	read_dataset(dset, coef_);	694,281✔
123	}	694,281✔
124
125	double Polynomial::operator()(double x) const	1,218,609,493✔
126	{
127	// Use Horner's rule to evaluate polynomial. Note that coefficients are
128	// ordered in increasing powers of x.
129	double y = 0.0;	1,218,609,493✔
130	for (auto c = coef_.crbegin(); c != coef_.crend(); ++c) {	2,147,483,647✔
131	y = y * x + *c;	2,147,483,647✔
132	}
133	return y;	1,218,609,493✔
134	}
135
136	//==============================================================================
137	// Tabulated1D implementation
138	//==============================================================================
139
140	Tabulated1D::Tabulated1D(hid_t dset)	2,602,415✔
141	{
142	read_attribute(dset, "breakpoints", nbt_);	2,602,415✔
143	n_regions_ = nbt_.size();	2,602,415✔
144
145	// Change 1-indexing to 0-indexing
146	for (auto& b : nbt_)	5,207,345✔
147	--b;	2,604,930✔
148
149	vector<int> int_temp;	2,602,415✔
150	read_attribute(dset, "interpolation", int_temp);	2,602,415✔
151
152	// Convert vector of ints into Interpolation
153	for (const auto i : int_temp)	5,207,345✔
154	int_.push_back(int2interp(i));	2,604,930✔
155
156	xt::xarray<double> arr;	2,602,415✔
157	read_dataset(dset, arr);	2,602,415✔
158
159	auto xs = xt::view(arr, 0);	2,602,415✔
160	auto ys = xt::view(arr, 1);	2,602,415✔
161
162	std::copy(xs.begin(), xs.end(), std::back_inserter(x_));	2,602,415✔
163	std::copy(ys.begin(), ys.end(), std::back_inserter(y_));	2,602,415✔
164	n_pairs_ = x_.size();	2,602,415✔
165	}	2,602,415✔
166
167	double Tabulated1D::operator()(double x) const	2,147,483,647✔
168	{
169	// find which bin the abscissa is in -- if the abscissa is outside the
170	// tabulated range, the first or last point is chosen, i.e. no interpolation
171	// is done outside the energy range
172	int i;
173	if (x < x_[0]) {	2,147,483,647✔
174	return y_[0];	1,732,351,644✔
175	} else if (x > x_[n_pairs_ - 1]) {	2,147,483,647✔
176	return y_[n_pairs_ - 1];	38,425,163✔
177	} else {
178	i = lower_bound_index(x_.begin(), x_.end(), x);	2,147,483,647✔
179	}
180
181	// determine interpolation scheme
182	Interpolation interp;
183	if (n_regions_ == 0) {	2,147,483,647!
184	interp = Interpolation::lin_lin;	×
185	} else {
186	interp = int_[0];	2,147,483,647✔
187	for (int j = 0; j < n_regions_; ++j) {	2,147,483,647!
188	if (i < nbt_[j]) {	2,147,483,647✔
189	interp = int_[j];	2,147,483,647✔
190	break;	2,147,483,647✔
191	}
192	}
193	}
194
195	// handle special case of histogram interpolation
196	if (interp == Interpolation::histogram)	2,147,483,647✔
197	return y_[i];	2,147,483,647✔
198
199	// determine bounding values
200	double x0 = x_[i];	2,147,483,647✔
201	double x1 = x_[i + 1];	2,147,483,647✔
202	double y0 = y_[i];	2,147,483,647✔
203	double y1 = y_[i + 1];	2,147,483,647✔
204
205	// determine interpolation factor and interpolated value
206	double r;
207	switch (interp) {	2,147,483,647!
208	case Interpolation::lin_lin:	2,147,483,647✔
209	r = (x - x0) / (x1 - x0);	2,147,483,647✔
210	return y0 + r * (y1 - y0);	2,147,483,647✔
211	case Interpolation::lin_log:	×
212	r = log(x / x0) / log(x1 / x0);	×
213	return y0 + r * (y1 - y0);	×
214	case Interpolation::log_lin:	×
215	r = (x - x0) / (x1 - x0);	×
216	return y0 * exp(r * log(y1 / y0));	×
217	case Interpolation::log_log:	22,808,627✔
218	r = log(x / x0) / log(x1 / x0);	22,808,627✔
219	return y0 * exp(r * log(y1 / y0));	22,808,627✔
220	default:	×
221	throw std::runtime_error {"Invalid interpolation scheme."};	×
222	}
223	}
224
225	//==============================================================================
226	// CoherentElasticXS implementation
227	//==============================================================================
228
229	CoherentElasticXS::CoherentElasticXS(hid_t dset)	140✔
230	{
231	// Read 2D array from dataset
232	xt::xarray<double> arr;	140✔
233	read_dataset(dset, arr);	140✔
234
235	// Get views for Bragg edges and structure factors
236	auto E = xt::view(arr, 0);	140✔
237	auto s = xt::view(arr, 1);	140✔
238
239	// Copy Bragg edges and partial sums of structure factors
240	std::copy(E.begin(), E.end(), std::back_inserter(bragg_edges_));	140✔
241	std::copy(s.begin(), s.end(), std::back_inserter(factors_));	140✔
242	}	140✔
243
244	double CoherentElasticXS::operator()(double E) const	4,939,957✔
245	{
246	if (E < bragg_edges_[0]) {	4,939,957✔
247	// If energy is below that of the lowest Bragg peak, the elastic cross
248	// section will be zero
249	return 0.0;	3,102✔
250	} else {
251	auto i_grid =
252	lower_bound_index(bragg_edges_.begin(), bragg_edges_.end(), E);	4,936,855✔
253	return factors_[i_grid] / E;	4,936,855✔
254	}
255	}
256
257	//==============================================================================
258	// IncoherentElasticXS implementation
259	//==============================================================================
260
261	IncoherentElasticXS::IncoherentElasticXS(hid_t dset)	×
262	{
263	array<double, 2> tmp;
264	read_dataset(dset, nullptr, tmp);	×
265	bound_xs_ = tmp[0];	×
266	debye_waller_ = tmp[1];	×
UNCOV 267	}	×
268
269	double IncoherentElasticXS::operator()(double E) const	×
270	{
271	// Determine cross section using ENDF-102, Eq. (7.5)
272	double W = debye_waller_;	×
273	return bound_xs_ / 2.0 * ((1 - std::exp(-4.0 * E * W)) / (2.0 * E * W));	×
274	}
275
276	//==============================================================================
277	// Sum1D implementation
278	//==============================================================================
279
280	Sum1D::Sum1D(hid_t group)	44✔
281	{
282	// Get number of functions
283	int n;
284	read_attribute(group, "n", n);	44✔
285
286	// Get each function
287	for (int i = 0; i < n; ++i) {	132✔
288	auto dset_name = fmt::format("func_{}", i + 1);	160✔
289	functions_.push_back(read_function(group, dset_name.c_str()));	88✔
290	}	88✔
291	}	44✔
292
293	double Sum1D::operator()(double x) const	4,049,320✔
294	{
295	double result = 0.0;	4,049,320✔
296	for (auto& func : functions_) {	12,147,960✔
297	result += (*func)(x);	8,098,640✔
298	}
299	return result;	4,049,320✔
300	}
301
302	} // namespace openmc

openmc-dev / openmc / 18656608910

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