22263374650

Committed 21 Feb 2026 08:00PM UTC coverage: 81.66% (-0.2%) from 81.839%

Build # 22263374650

Build Type

push

github

Committed by

web-flow

Commit Message

Support arbitrary symmetry axis for CylindricalIndependent class (#3474)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

Run Details

16814 of 23454 branches covered (71.69%)

Branch coverage included in aggregate %.

50 of 64 new or added lines in 2 files covered. (78.13%)

696 existing lines in 31 files now uncovered.

56816 of 66713 relevant lines covered (85.16%)

43189734.94 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:
    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)	46,438,628✔
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) {	46,438,628!
28	case 1:	9,119,185✔
29	case 11:
30	case 21:
31	return Interpolation::histogram;	9,119,185✔
32	case 2:	37,316,484✔
33	case 12:
34	case 22:
35	return Interpolation::lin_lin;	37,316,484✔
36	case 3:	×
37	case 13:
38	case 23:
39	return Interpolation::lin_log;	×
40	case 4:	×
41	case 14:
42	case 24:
43	return Interpolation::log_lin;	×
44	case 5:	2,959✔
45	case 15:
46	case 25:
47	return Interpolation::log_log;	2,959✔
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,017,223✔
60	{
61	if (mt >= N_DISAPPEAR && mt <= N_DA) {	1,017,223✔
62	return true;	172,426✔
63	} else if (mt >= N_P0 && mt <= N_AC) {	844,797!
64	return true;	285,302✔
65	} else if (mt == N_TA \|\| mt == N_DT \|\| mt == N_P3HE \|\| mt == N_D3HE \|\|	559,495!
66	mt == N_3HEA \|\| mt == N_3P) {	559,495!
67	return true;	×
68	} else {
69	return false;	559,495✔
70	}
71	}
72
73	bool is_inelastic_scatter(int mt)	1,042,161✔
74	{
75	if (mt < 100) {	1,042,161✔
76	if (is_fission(mt)) {	574,087✔
77	return false;	12,725✔
78	} else {
79	return mt >= MISC && mt != 27;	561,362!
80	}
81	} else if (mt <= 200) {	468,074✔
82	return !is_disappearance(mt);	79,745✔
83	} else if (mt >= N_2N0 && mt <= N_2NC) {	388,329!
84	return true;	×
85	} else {
86	return false;	388,329✔
87	}
88	}
89
90	bool mt_matches(int event_mt, int target_mt)	1,901,655✔
91	{
92	// Direct match
93	if (event_mt == target_mt)	1,901,655✔
94	return true;	177,696✔
95
96	// Check if event_mt is a component of target_mt summation reaction
97	switch (target_mt) {	1,723,959!
98	case TOTAL_XS:	261,279✔
99	return event_mt == ELASTIC \|\| mt_matches(event_mt, N_NONELASTIC);	261,279!
100
101	case N_NONELASTIC: {	94,212✔
102	static constexpr int components[] = {4, 5, 11, 16, 17, 22, 23, 24, 25, 27,
103	28, 29, 30, 32, 33, 34, 35, 36, 37, 41, 42, 44, 45, 152, 153, 154, 156,
104	157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171,
105	172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 183, 184, 185, 186, 187,
106	188, 189, 190, 194, 195, 196, 198, 199, 200};
107	for (int mt : components) {	501,048!
108	if (mt_matches(event_mt, mt))	501,048✔
109	return true;	94,212✔
110	}
UNCOV 111	return false;	×
112	}
113
114	case N_LEVEL:	94,212✔
115	// Inelastic scattering levels
116	return event_mt >= 50 && event_mt <= N_NC;	94,212✔
117
118	case N_2N:	306,207✔
119	// (n,2n) to excited states
120	return event_mt >= N_2N0 && event_mt <= N_2NC;	306,207!
121
122	case N_FISSION:	251,064✔
123	return is_fission(event_mt);	251,064✔
124
125	case 27:	45,135✔
126	return is_fission(event_mt) \|\| is_disappearance(event_mt);	45,135!
127
UNCOV 128	case N_DISAPPEAR: {	×
UNCOV 129	return is_disappearance(event_mt);	×
130	}
131
UNCOV 132	case N_P:	×
133	// (n,p) to excited states
UNCOV 134	return event_mt >= N_P0 && event_mt <= N_PC;	×
135
UNCOV 136	case N_D:	×
137	// (n,d) to excited states
UNCOV 138	return event_mt >= N_D0 && event_mt <= N_DC;	×
139
UNCOV 140	case N_T:	×
141	// (n,t) to excited states
UNCOV 142	return event_mt >= N_T0 && event_mt <= N_TC;	×
143
UNCOV 144	case N_3HE:	×
145	// (n,3He) to excited states
UNCOV 146	return event_mt >= N_3HE0 && event_mt <= N_3HEC;	×
147
UNCOV 148	case N_A:	×
149	// (n,alpha) to excited states
UNCOV 150	return event_mt >= N_A0 && event_mt <= N_AC;	×
151
UNCOV 152	case 501:	×
UNCOV 153	return event_mt == 502 \|\| event_mt == 504 \|\| mt_matches(event_mt, 516) \|\|	×
UNCOV 154	mt_matches(event_mt, 522);	×
155
UNCOV 156	case PAIR_PROD:	×
UNCOV 157	return event_mt == PAIR_PROD_ELEC \|\| event_mt == PAIR_PROD_NUC;	×
158
UNCOV 159	case PHOTOELECTRIC:	×
UNCOV 160	return event_mt >= 534 && event_mt < 573;	×
161
162	default:	671,850✔
163	return false;	671,850✔
164	}
165	}
166
167	unique_ptr<Function1D> read_function(hid_t group, const char* name)	2,727,745✔
168	{
169	hid_t obj_id = open_object(group, name);	2,727,745✔
170	std::string func_type;	2,727,745✔
171	read_attribute(obj_id, "type", func_type);	2,727,745✔
172	unique_ptr<Function1D> func;	2,727,745✔
173	if (func_type == "Tabulated1D") {	2,727,745✔
174	func = make_unique<Tabulated1D>(obj_id);	2,154,597✔
175	} else if (func_type == "Polynomial") {	573,148✔
176	func = make_unique<Polynomial>(obj_id);	573,004✔
177	} else if (func_type == "CoherentElastic") {	144✔
178	func = make_unique<CoherentElasticXS>(obj_id);	108✔
179	} else if (func_type == "IncoherentElastic") {	36!
UNCOV 180	func = make_unique<IncoherentElasticXS>(obj_id);	×
181	} else if (func_type == "Sum") {	36!
182	func = make_unique<Sum1D>(obj_id);	36✔
183	} else {
UNCOV 184	throw std::runtime_error {"Unknown function type " + func_type +	×
UNCOV 185	" for dataset " + object_name(obj_id)};	×
186	}
187	close_object(obj_id);	2,727,745✔
188	return func;	5,455,490✔
189	}	2,727,745✔
190
191	//==============================================================================
192	// Polynomial implementation
193	//==============================================================================
194
195	Polynomial::Polynomial(hid_t dset)	573,004✔
196	{
197	// Read coefficients into a vector
198	read_dataset(dset, coef_);	573,004✔
199	}	573,004✔
200
201	double Polynomial::operator()(double x) const	998,377,528✔
202	{
203	// Use Horner's rule to evaluate polynomial. Note that coefficients are
204	// ordered in increasing powers of x.
205	double y = 0.0;	998,377,528✔
206	for (auto c = coef_.crbegin(); c != coef_.crend(); ++c) {	2,147,483,647✔
207	y = y * x + *c;	2,147,483,647✔
208	}
209	return y;	998,377,528✔
210	}
211
212	//==============================================================================
213	// Tabulated1D implementation
214	//==============================================================================
215
216	Tabulated1D::Tabulated1D(hid_t dset)	2,163,250✔
217	{
218	read_attribute(dset, "breakpoints", nbt_);	2,163,250✔
219	n_regions_ = nbt_.size();	2,163,250✔
220
221	// Change 1-indexing to 0-indexing
222	for (auto& b : nbt_)	4,328,558✔
223	--b;	2,165,308✔
224
225	vector<int> int_temp;	2,163,250✔
226	read_attribute(dset, "interpolation", int_temp);	2,163,250✔
227
228	// Convert vector of ints into Interpolation
229	for (const auto i : int_temp)	4,328,558✔
230	int_.push_back(int2interp(i));	2,165,308✔
231
232	tensor::Tensor<double> arr;	2,163,250✔
233	read_dataset(dset, arr);	2,163,250✔
234
235	tensor::View<double> xs = arr.slice(0);	2,163,250✔
236	tensor::View<double> ys = arr.slice(1);	2,163,250✔
237
238	std::copy(xs.begin(), xs.end(), std::back_inserter(x_));	2,163,250✔
239	std::copy(ys.begin(), ys.end(), std::back_inserter(y_));	2,163,250✔
240	n_pairs_ = x_.size();	2,163,250✔
241	}	2,163,250✔
242
243	double Tabulated1D::operator()(double x) const	2,147,483,647✔
244	{
245	// find which bin the abscissa is in -- if the abscissa is outside the
246	// tabulated range, the first or last point is chosen, i.e. no interpolation
247	// is done outside the energy range
248	int i;
249	if (x < x_[0]) {	2,147,483,647✔
250	return y_[0];	1,420,870,983✔
251	} else if (x > x_[n_pairs_ - 1]) {	2,147,483,647✔
252	return y_[n_pairs_ - 1];	33,149,481✔
253	} else {
254	i = lower_bound_index(x_.begin(), x_.end(), x);	2,147,483,647✔
255	}
256
257	// determine interpolation scheme
258	Interpolation interp;
259	if (n_regions_ == 0) {	2,147,483,647!
260	interp = Interpolation::lin_lin;	×
261	} else {
262	interp = int_[0];	2,147,483,647✔
263	for (int j = 0; j < n_regions_; ++j) {	2,147,483,647!
264	if (i < nbt_[j]) {	2,147,483,647✔
265	interp = int_[j];	2,147,483,647✔
266	break;	2,147,483,647✔
267	}
268	}
269	}
270
271	// handle special case of histogram interpolation
272	if (interp == Interpolation::histogram)	2,147,483,647✔
273	return y_[i];	2,147,483,647✔
274
275	// determine bounding values
276	double x0 = x_[i];	2,147,483,647✔
277	double x1 = x_[i + 1];	2,147,483,647✔
278	double y0 = y_[i];	2,147,483,647✔
279	double y1 = y_[i + 1];	2,147,483,647✔
280
281	// determine interpolation factor and interpolated value
282	double r;
283	switch (interp) {	2,147,483,647!
284	case Interpolation::lin_lin:	2,147,483,647✔
285	r = (x - x0) / (x1 - x0);	2,147,483,647✔
286	return y0 + r * (y1 - y0);	2,147,483,647✔
UNCOV 287	case Interpolation::lin_log:	×
UNCOV 288	r = log(x / x0) / log(x1 / x0);	×
UNCOV 289	return y0 + r * (y1 - y0);	×
UNCOV 290	case Interpolation::log_lin:	×
UNCOV 291	r = (x - x0) / (x1 - x0);	×
UNCOV 292	return y0 * exp(r * log(y1 / y0));	×
293	case Interpolation::log_log:	19,132,265✔
294	r = log(x / x0) / log(x1 / x0);	19,132,265✔
295	return y0 * exp(r * log(y1 / y0));	19,132,265✔
UNCOV 296	default:	×
UNCOV 297	throw std::runtime_error {"Invalid interpolation scheme."};	×
298	}
299	}
300
301	//==============================================================================
302	// CoherentElasticXS implementation
303	//==============================================================================
304
305	CoherentElasticXS::CoherentElasticXS(hid_t dset)	108✔
306	{
307	// Read 2D array from dataset
308	tensor::Tensor<double> arr;	108✔
309	read_dataset(dset, arr);	108✔
310
311	// Get views for Bragg edges and structure factors
312	tensor::View<double> E = arr.slice(0);	108✔
313	tensor::View<double> s = arr.slice(1);	108✔
314
315	// Copy Bragg edges and partial sums of structure factors
316	std::copy(E.begin(), E.end(), std::back_inserter(bragg_edges_));	108✔
317	std::copy(s.begin(), s.end(), std::back_inserter(factors_));	108✔
318	}	108✔
319
320	double CoherentElasticXS::operator()(double E) const	4,041,783✔
321	{
322	if (E < bragg_edges_[0]) {	4,041,783✔
323	// If energy is below that of the lowest Bragg peak, the elastic cross
324	// section will be zero
325	return 0.0;	2,538✔
326	} else {
327	auto i_grid =
328	lower_bound_index(bragg_edges_.begin(), bragg_edges_.end(), E);	4,039,245✔
329	return factors_[i_grid] / E;	4,039,245✔
330	}
331	}
332
333	//==============================================================================
334	// IncoherentElasticXS implementation
335	//==============================================================================
336
UNCOV 337	IncoherentElasticXS::IncoherentElasticXS(hid_t dset)	×
338	{
339	array<double, 2> tmp;
UNCOV 340	read_dataset(dset, nullptr, tmp);	×
UNCOV 341	bound_xs_ = tmp[0];	×
UNCOV 342	debye_waller_ = tmp[1];	×
UNCOV 343	}	×
344
UNCOV 345	double IncoherentElasticXS::operator()(double E) const	×
346	{
347	// Determine cross section using ENDF-102, Eq. (7.5)
UNCOV 348	double W = debye_waller_;	×
UNCOV 349	return bound_xs_ / 2.0 * ((1 - std::exp(-4.0 * E * W)) / (2.0 * E * W));	×
350	}
351
352	//==============================================================================
353	// Sum1D implementation
354	//==============================================================================
355
356	Sum1D::Sum1D(hid_t group)	36✔
357	{
358	// Get number of functions
359	int n;
360	read_attribute(group, "n", n);	36✔
361
362	// Get each function
363	for (int i = 0; i < n; ++i) {	108✔
364	auto dset_name = fmt::format("func_{}", i + 1);	128✔
365	functions_.push_back(read_function(group, dset_name.c_str()));	72✔
366	}	72✔
367	}	36✔
368
369	double Sum1D::operator()(double x) const	3,313,080✔
370	{
371	double result = 0.0;	3,313,080✔
372	for (auto& func : functions_) {	9,939,240✔
373	result += (*func)(x);	6,626,160✔
374	}
375	return result;	3,313,080✔
376	}
377
378	} // namespace openmc

openmc-dev / openmc / 22263374650

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