• Home
  • Features
  • Pricing
  • Docs
  • Announcements
  • Sign In

openmc-dev / openmc / 22026428798

14 Feb 2026 11:58PM UTC coverage: 81.582% (-0.1%) from 81.72%
22026428798

Pull #3803

github

web-flow
Merge 9fe9422d1 into a35927aad
Pull Request #3803: Set upper and lower interpolation bounds for MGXS data.

17011 of 23656 branches covered (71.91%)

Branch coverage included in aggregate %.

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

332 existing lines in 25 files now uncovered.

56150 of 66022 relevant lines covered (85.05%)

42033857.66 hits per line

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

79.11
/src/mgxs.cpp
1
#include "openmc/mgxs.h"
2

3
#include <algorithm>
4
#include <cmath>
5
#include <cstdlib>
6
#include <sstream>
7

8
#include "xtensor/xadapt.hpp"
9
#include "xtensor/xmath.hpp"
10
#include "xtensor/xsort.hpp"
11
#include "xtensor/xview.hpp"
12
#include <fmt/core.h>
13

14
#include "openmc/error.h"
15
#include "openmc/math_functions.h"
16
#include "openmc/mgxs_interface.h"
17
#include "openmc/nuclide.h"
18
#include "openmc/random_lcg.h"
19
#include "openmc/settings.h"
20
#include "openmc/string_utils.h"
21

22
namespace openmc {
23

24
//==============================================================================
25
// Mgxs base-class methods
26
//==============================================================================
27

28
void Mgxs::init(const std::string& in_name, double in_awr,
4,800✔
29
  const vector<double>& in_kTs, bool in_fissionable,
30
  AngleDistributionType in_scatter_format, bool in_is_isotropic,
31
  const vector<double>& in_polar, const vector<double>& in_azimuthal)
32
{
33
  // Set the metadata
34
  name = in_name;
4,800✔
35
  awr = in_awr;
4,800✔
36
  // TODO: Remove adapt when in_KTs is an xtensor
37
  kTs = xt::adapt(in_kTs);
4,800✔
38
  fissionable = in_fissionable;
4,800✔
39
  scatter_format = in_scatter_format;
4,800✔
40
  xs.resize(in_kTs.size());
4,800✔
41
  is_isotropic = in_is_isotropic;
4,800✔
42
  n_pol = in_polar.size();
4,800✔
43
  n_azi = in_azimuthal.size();
4,800✔
44
  polar = in_polar;
4,800✔
45
  azimuthal = in_azimuthal;
4,800✔
46
}
4,800✔
47

48
//==============================================================================
49

50
void Mgxs::metadata_from_hdf5(hid_t xs_id, const vector<double>& temperature,
2,562✔
51
  vector<int>& temps_to_read, int& order_dim)
52
{
53
  // get name
54
  std::string in_name;
2,562✔
55
  get_name(xs_id, in_name);
2,562✔
56
  // remove the leading '/'
57
  in_name = in_name.substr(1);
2,562✔
58

59
  // Get the AWR
60
  double in_awr;
61
  if (attribute_exists(xs_id, "atomic_weight_ratio")) {
2,562✔
62
    read_attr_double(xs_id, "atomic_weight_ratio", &in_awr);
96✔
63
  } else {
64
    in_awr = MACROSCOPIC_AWR;
2,466✔
65
  }
66

67
  // Determine the available temperatures
68
  hid_t kT_group = open_group(xs_id, "kTs");
2,562✔
69
  size_t num_temps = get_num_datasets(kT_group);
2,562✔
70
  char** dset_names = new char*[num_temps];
2,562!
71
  for (int i = 0; i < num_temps; i++) {
5,988✔
72
    dset_names[i] = new char[151];
3,426✔
73
  }
74
  get_datasets(kT_group, dset_names);
2,562✔
75
  vector<size_t> shape = {num_temps};
2,562✔
76
  xt::xarray<double> temps_available(shape);
2,562✔
77
  for (int i = 0; i < num_temps; i++) {
5,988✔
78
    read_double(kT_group, dset_names[i], &temps_available[i], true);
3,426✔
79

80
    // convert eV to Kelvin
81
    temps_available[i] = std::round(temps_available[i] / K_BOLTZMANN);
3,426✔
82

83
    // Done with dset_names, so delete it
84
    delete[] dset_names[i];
3,426!
85
  }
86
  delete[] dset_names;
2,562!
87
  std::sort(temps_available.begin(), temps_available.end());
2,562✔
88

89
  // Set the global upper and lower interpolation bounds to avoid errors
90
  // involving C-API functions.
91
  data::temperature_min = std::min(data::temperature_min, *temps_available.begin());
2,562✔
92
  data::temperature_max = std::max(data::temperature_max, *temps_available.end());
2,562✔
93

94
  // If only one temperature is available, lets just use nearest temperature
95
  // interpolation
96
  if ((num_temps == 1) &&
2,562✔
97
      (settings::temperature_method == TemperatureMethod::INTERPOLATION)) {
2,058!
UNCOV
98
    warning("Cross sections for " + strtrim(name) + " are only available " +
×
UNCOV
99
            "at one temperature.  Reverting to the nearest temperature " +
×
100
            "method.");
101
    settings::temperature_method = TemperatureMethod::NEAREST;
×
102
  }
103

104
  switch (settings::temperature_method) {
2,562!
105
  case TemperatureMethod::NEAREST:
2,418✔
106
    // Determine actual temperatures to read
107
    for (const auto& T : temperature) {
4,812✔
108
      // Determine the closest temperature value
109
      // NOTE: the below block could be replaced with the following line,
110
      // though this gives a runtime error if using LLVM 20 or newer,
111
      // likely due to a bug in xtensor.
112
      // auto i_closest = xt::argmin(xt::abs(temps_available - T))[0];
113
      double closest = std::numeric_limits<double>::max();
2,394✔
114
      int i_closest = 0;
2,394✔
115
      for (int i = 0; i < temps_available.size(); i++) {
5,388✔
116
        double diff = std::abs(temps_available[i] - T);
2,994✔
117
        if (diff < closest) {
2,994✔
118
          closest = diff;
2,670✔
119
          i_closest = i;
2,670✔
120
        }
121
      }
122

123
      double temp_actual = temps_available[i_closest];
2,394✔
124
      if (std::fabs(temp_actual - T) < settings::temperature_tolerance) {
2,394!
125
        if (std::find(temps_to_read.begin(), temps_to_read.end(),
2,394✔
126
              std::round(temp_actual)) == temps_to_read.end()) {
4,788!
127
          temps_to_read.push_back(std::round(temp_actual));
2,394✔
128
        }
129
      } else {
UNCOV
130
        fatal_error(fmt::format(
×
131
          "MGXS library does not contain cross sections "
132
          "for {} at or near {} K. Available temperatures "
133
          "are {} K. Consider making use of openmc.Settings.temperature "
134
          "to specify how intermediate temperatures are treated.",
UNCOV
135
          in_name, std::round(T), concatenate(temps_available)));
×
136
      }
137
    }
138
    break;
2,418✔
139

140
  case TemperatureMethod::INTERPOLATION:
144✔
141
    for (int i = 0; i < temperature.size(); i++) {
288✔
142
      for (int j = 0; j < num_temps; j++) {
216!
143
        if (j == (num_temps - 1)) {
216!
UNCOV
144
          fatal_error("MGXS Library does not contain cross sections for " +
×
UNCOV
145
                      in_name + " at temperatures that bound " +
×
146
                      std::to_string(std::round(temperature[i])));
×
147
        }
148
        if ((temps_available[j] <= temperature[i]) &&
432!
149
            (temperature[i] < temps_available[j + 1])) {
216✔
150
          if (std::find(temps_to_read.begin(), temps_to_read.end(),
144✔
151
                temps_available[j]) == temps_to_read.end()) {
288!
152
            temps_to_read.push_back(temps_available[j]);
144✔
153
          }
154

155
          if (std::find(temps_to_read.begin(), temps_to_read.end(),
144✔
156
                temps_available[j + 1]) == temps_to_read.end()) {
288!
157
            temps_to_read.push_back(temps_available[j + 1]);
144✔
158
          }
159
          break;
144✔
160
        }
161
      }
162
    }
163
  }
164
  std::sort(temps_to_read.begin(), temps_to_read.end());
2,562✔
165

166
  // Get the library's temperatures
167
  int n_temperature = temps_to_read.size();
2,562✔
168
  vector<double> in_kTs(n_temperature);
5,124✔
169
  for (int i = 0; i < n_temperature; i++) {
5,244✔
170
    std::string temp_str(std::to_string(temps_to_read[i]) + "K");
2,682✔
171

172
    // read exact temperature value
173
    read_double(kT_group, temp_str.c_str(), &in_kTs[i], true);
2,682✔
174
  }
2,682✔
175
  close_group(kT_group);
2,562✔
176

177
  // Load the remaining metadata
178
  AngleDistributionType in_scatter_format;
179
  if (attribute_exists(xs_id, "scatter_format")) {
2,562!
180
    std::string temp_str;
2,562✔
181
    read_attribute(xs_id, "scatter_format", temp_str);
2,562✔
182
    to_lower(strtrim(temp_str));
2,562✔
183
    if (temp_str.compare(0, 8, "legendre") == 0) {
2,562✔
184
      in_scatter_format = AngleDistributionType::LEGENDRE;
2,466✔
185
    } else if (temp_str.compare(0, 9, "histogram") == 0) {
96✔
186
      in_scatter_format = AngleDistributionType::HISTOGRAM;
48✔
187
    } else if (temp_str.compare(0, 7, "tabular") == 0) {
48!
188
      in_scatter_format = AngleDistributionType::TABULAR;
48✔
189
    } else {
UNCOV
190
      fatal_error("Invalid scatter_format option!");
×
191
    }
192
  } else {
2,562✔
UNCOV
193
    in_scatter_format = AngleDistributionType::LEGENDRE;
×
194
  }
195

196
  if (attribute_exists(xs_id, "scatter_shape")) {
2,562!
197
    std::string temp_str;
2,562✔
198
    read_attribute(xs_id, "scatter_shape", temp_str);
2,562✔
199
    to_lower(strtrim(temp_str));
2,562✔
200
    if (temp_str.compare(0, 14, "[g][g\'][order]") != 0) {
2,562!
UNCOV
201
      fatal_error("Invalid scatter_shape option!");
×
202
    }
203
  }
2,562✔
204

205
  bool in_fissionable = false;
2,562✔
206
  if (attribute_exists(xs_id, "fissionable")) {
2,562!
207
    int int_fiss;
208
    read_attr_int(xs_id, "fissionable", &int_fiss);
2,562✔
209
    in_fissionable = int_fiss;
2,562✔
210
  } else {
UNCOV
211
    fatal_error("Fissionable element must be set!");
×
212
  }
213

214
  // Get the library's value for the order
215
  if (attribute_exists(xs_id, "order")) {
2,562!
216
    read_attr_int(xs_id, "order", &order_dim);
2,562✔
217
  } else {
UNCOV
218
    fatal_error("Order must be provided!");
×
219
  }
220

221
  // Store the dimensionality of the data in order_dim.
222
  // For Legendre data, we usually refer to it as Pn where n is the order.
223
  // However Pn has n+1 sets of points (since you need to count the P0
224
  // moment). Adjust for that. Histogram and Tabular formats dont need this
225
  // adjustment.
226
  if (in_scatter_format == AngleDistributionType::LEGENDRE) {
2,562✔
227
    ++order_dim;
2,466✔
228
  }
229

230
  // Get the angular information
231
  int in_n_pol;
232
  int in_n_azi;
233
  bool in_is_isotropic = true;
2,562✔
234
  if (attribute_exists(xs_id, "representation")) {
2,562!
235
    std::string temp_str;
2,562✔
236
    read_attribute(xs_id, "representation", temp_str);
2,562✔
237
    to_lower(strtrim(temp_str));
2,562✔
238
    if (temp_str.compare(0, 5, "angle") == 0) {
2,562✔
239
      in_is_isotropic = false;
24✔
240
    } else if (temp_str.compare(0, 9, "isotropic") != 0) {
2,538!
UNCOV
241
      fatal_error("Invalid Data Representation!");
×
242
    }
243
  }
2,562✔
244

245
  if (!in_is_isotropic) {
2,562✔
246
    if (attribute_exists(xs_id, "num_polar")) {
24!
247
      read_attr_int(xs_id, "num_polar", &in_n_pol);
24✔
248
    } else {
UNCOV
249
      fatal_error("num_polar must be provided!");
×
250
    }
251
    if (attribute_exists(xs_id, "num_azimuthal")) {
24!
252
      read_attr_int(xs_id, "num_azimuthal", &in_n_azi);
24✔
253
    } else {
UNCOV
254
      fatal_error("num_azimuthal must be provided!");
×
255
    }
256
  } else {
257
    in_n_pol = 1;
2,538✔
258
    in_n_azi = 1;
2,538✔
259
  }
260

261
  // Set the angular bins to use equally-spaced bins
262
  vector<double> in_polar(in_n_pol);
5,124✔
263
  double dangle = PI / in_n_pol;
2,562✔
264
  for (int p = 0; p < in_n_pol; p++) {
5,148✔
265
    in_polar[p] = (p + 0.5) * dangle;
2,586✔
266
  }
267
  vector<double> in_azimuthal(in_n_azi);
5,124✔
268
  dangle = 2. * PI / in_n_azi;
2,562✔
269
  for (int a = 0; a < in_n_azi; a++) {
5,148✔
270
    in_azimuthal[a] = (a + 0.5) * dangle - PI;
2,586✔
271
  }
272

273
  // Finally use this data to initialize the MGXS Object
274
  init(in_name, in_awr, in_kTs, in_fissionable, in_scatter_format,
2,562✔
275
    in_is_isotropic, in_polar, in_azimuthal);
276
}
2,562✔
277

278
//==============================================================================
279

280
Mgxs::Mgxs(
2,562✔
281
  hid_t xs_id, const vector<double>& temperature, int num_group, int num_delay)
2,562✔
282
  : num_groups(num_group), num_delayed_groups(num_delay)
2,562✔
283
{
284
  // Call generic data gathering routine (will populate the metadata)
285
  int order_data;
286
  vector<int> temps_to_read;
2,562✔
287
  metadata_from_hdf5(xs_id, temperature, temps_to_read, order_data);
2,562✔
288

289
  // Set number of energy and delayed groups
290
  AngleDistributionType final_scatter_format = scatter_format;
2,562✔
291
  if (settings::legendre_to_tabular) {
2,562✔
292
    if (scatter_format == AngleDistributionType::LEGENDRE)
2,370✔
293
      final_scatter_format = AngleDistributionType::TABULAR;
2,298✔
294
  }
295

296
  // Load the more specific XsData information
297
  for (int t = 0; t < temps_to_read.size(); t++) {
5,244✔
298
    xs[t] = XsData(fissionable, final_scatter_format, n_pol, n_azi, num_groups,
5,364✔
299
      num_delayed_groups);
5,364✔
300
    // Get the temperature as a string and then open the HDF5 group
301
    std::string temp_str = std::to_string(temps_to_read[t]) + "K";
2,682✔
302
    hid_t xsdata_grp = open_group(xs_id, temp_str.c_str());
2,682✔
303

304
    xs[t].from_hdf5(xsdata_grp, fissionable, scatter_format,
2,682✔
305
      final_scatter_format, order_data, is_isotropic, n_pol, n_azi);
2,682✔
306
    close_group(xsdata_grp);
2,682✔
307

308
  } // end temperature loop
2,682✔
309

310
  // Make sure the scattering format is updated to the final case
311
  scatter_format = final_scatter_format;
2,562✔
312
}
2,562✔
313

314
//==============================================================================
315

316
Mgxs::Mgxs(const std::string& in_name, const vector<double>& mat_kTs,
2,238✔
317
  const vector<Mgxs*>& micros, const vector<double>& atom_densities,
318
  int num_group, int num_delay)
2,238✔
319
  : num_groups(num_group), num_delayed_groups(num_delay)
2,238✔
320
{
321
  // Get the minimum data needed to initialize:
322
  // Dont need awr, but lets just initialize it anyways
323
  double in_awr = -1.;
2,238✔
324
  // start with the assumption it is not fissionable and set
325
  // the fissionable status if we learn differently
326
  bool in_fissionable = false;
2,238✔
327
  for (int m = 0; m < micros.size(); m++) {
4,764✔
328
    if (micros[m]->fissionable)
2,526✔
329
      in_fissionable = true;
897✔
330
  }
331
  // Force all of the following data to be the same; these will be verified
332
  // to be true later
333
  AngleDistributionType in_scatter_format = micros[0]->scatter_format;
2,238✔
334
  bool in_is_isotropic = micros[0]->is_isotropic;
2,238✔
335
  vector<double> in_polar = micros[0]->polar;
2,238✔
336
  vector<double> in_azimuthal = micros[0]->azimuthal;
2,238✔
337

338
  init(in_name, in_awr, mat_kTs, in_fissionable, in_scatter_format,
2,238✔
339
    in_is_isotropic, in_polar, in_azimuthal);
340

341
  // Create the xs data for each temperature
342
  for (int t = 0; t < mat_kTs.size(); t++) {
4,500✔
343
    xs[t] = XsData(in_fissionable, in_scatter_format, in_polar.size(),
4,524✔
344
      in_azimuthal.size(), num_groups, num_delayed_groups);
4,524✔
345

346
    // Find the right temperature index to use
347
    double temp_desired = mat_kTs[t];
2,262✔
348

349
    // Create the list of temperature indices and interpolation factors for
350
    // each microscopic data at the material temperature
351
    vector<int> micro_t(micros.size(), 0);
2,262✔
352
    vector<double> micro_t_interp(micros.size(), 0.);
2,262✔
353
    for (int m = 0; m < micros.size(); m++) {
4,812✔
354
      switch (settings::temperature_method) {
2,550!
355
      case TemperatureMethod::NEAREST: {
2,406✔
356
        micro_t[m] = xt::argmin(xt::abs(micros[m]->kTs - temp_desired))[0];
2,406✔
357
        auto temp_actual = micros[m]->kTs[micro_t[m]];
2,406✔
358

359
        if (std::abs(temp_actual - temp_desired) >=
2,406✔
360
            K_BOLTZMANN * settings::temperature_tolerance) {
2,406!
UNCOV
361
          fatal_error(fmt::format("MGXS Library does not contain cross section "
×
362
                                  "for {} at or near {} K.",
363
            name, std::round(temp_desired / K_BOLTZMANN)));
×
364
        }
365
      } break;
2,406✔
366
      case TemperatureMethod::INTERPOLATION:
144✔
367
        // Get a list of bounding temperatures for each actual temperature
368
        // present in the model
369
        for (int k = 0; k < micros[m]->kTs.shape()[0] - 1; k++) {
288✔
370
          if ((micros[m]->kTs[k] <= temp_desired) &&
288!
371
              (temp_desired < micros[m]->kTs[k + 1])) {
144!
372
            micro_t[m] = k;
144✔
373
            if (k == 0) {
144!
374
              micro_t_interp[m] = (temp_desired - micros[m]->kTs[k]) /
288✔
375
                                  (micros[m]->kTs[k + 1] - micros[m]->kTs[k]);
144✔
376
            } else {
UNCOV
377
              micro_t_interp[m] = 1.;
×
378
            }
379
          }
380
        }
381
      } // end switch
382
    }   // end microscopic temperature loop
383

384
    // Now combine the microscopic data at each relevant temperature
385
    // We will do this by treating the multiple temperatures of a nuclide as
386
    // a different nuclide. Mathematically this just means the temperature
387
    // interpolant is included in the number density.
388
    // These interpolants are contained within interpolant.
389
    vector<double> interpolant;    // the interpolant for the Mgxs
4,524✔
390
    vector<int> temp_indices;      // the temperature index for each Mgxs
4,524✔
391
    vector<Mgxs*> mgxs_to_combine; // The Mgxs to combine
4,524✔
392
    // Now go through and build the above vectors so that we can use them to
393
    // combine the data. We will step through each microscopic data and
394
    // add in its lower and upper temperature points
395
    for (int m = 0; m < micros.size(); m++) {
4,812✔
396
      if (settings::temperature_method == TemperatureMethod::NEAREST) {
2,550✔
397
        // Nearest interpolation only has one temperature point per isotope
398
        // and so we dont need to include a temperature interpolant in
399
        // the interpolant vector
400
        interpolant.push_back(atom_densities[m]);
2,406✔
401
        temp_indices.push_back(micro_t[m]);
2,406✔
402
        mgxs_to_combine.push_back(micros[m]);
2,406✔
403
      } else {
404
        // This will be an interpolation between two points so get both these
405
        // points
406
        // Start with the low point
407
        interpolant.push_back((1. - micro_t_interp[m]) * atom_densities[m]);
144✔
408
        temp_indices.push_back(micro_t[m]);
144✔
409
        mgxs_to_combine.push_back(micros[m]);
144✔
410
        // The higher point
411
        interpolant.push_back((micro_t_interp[m]) * atom_densities[m]);
144✔
412
        temp_indices.push_back(micro_t[m] + 1);
144✔
413
        mgxs_to_combine.push_back(micros[m]);
144✔
414
      }
415
    }
416

417
    // And finally, combine the data
418
    combine(mgxs_to_combine, interpolant, temp_indices, t);
2,262✔
419
  } // end temperature (t) loop
2,262✔
420
}
2,238✔
421

422
//==============================================================================
423

424
void Mgxs::combine(const vector<Mgxs*>& micros, const vector<double>& scalars,
2,262✔
425
  const vector<int>& micro_ts, int this_t)
426
{
427
  // Build the vector of pointers to the xs objects within micros
428
  vector<XsData*> those_xs(micros.size());
2,262✔
429
  for (int i = 0; i < micros.size(); i++) {
4,956✔
430
    those_xs[i] = &(micros[i]->xs[micro_ts[i]]);
2,694✔
431
  }
432

433
  xs[this_t].combine(those_xs, scalars);
2,262✔
434
}
2,262✔
435

436
//==============================================================================
437

438
double Mgxs::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu,
1,881,072,900✔
439
  const int* dg, int t, int a)
440
{
441
  XsData* xs_t = &xs[t];
1,881,072,900✔
442
  double val;
443
  switch (xstype) {
1,881,072,900!
444
  case MgxsType::TOTAL:
16,881,327✔
445
    val = xs_t->total(a, gin);
16,881,327✔
446
    break;
16,881,327✔
447
  case MgxsType::NU_FISSION:
8,343,162✔
448
    val = fissionable ? xs_t->nu_fission(a, gin) : 0.;
8,343,162✔
449
    break;
8,343,162✔
450
  case MgxsType::ABSORPTION:
4,766,490✔
451
    val = xs_t->absorption(a, gin);
4,766,490✔
452
    ;
453
    break;
4,766,490✔
454
  case MgxsType::FISSION:
9,139,122✔
455
    val = fissionable ? xs_t->fission(a, gin) : 0.;
9,139,122✔
456
    break;
9,139,122✔
457
  case MgxsType::KAPPA_FISSION:
8,542,638✔
458
    val = fissionable ? xs_t->kappa_fission(a, gin) : 0.;
8,542,638✔
459
    break;
8,542,638✔
460
  case MgxsType::NU_SCATTER:
15,299,226✔
461
  case MgxsType::SCATTER:
462
  case MgxsType::NU_SCATTER_FMU:
463
  case MgxsType::SCATTER_FMU:
464
    val = xs_t->scatter[a]->get_xs(xstype, gin, gout, mu);
15,299,226✔
465
    break;
15,299,226✔
466
  case MgxsType::PROMPT_NU_FISSION:
8,336,124✔
467
    val = fissionable ? xs_t->prompt_nu_fission(a, gin) : 0.;
8,336,124!
468
    break;
8,336,124✔
469
  case MgxsType::DELAYED_NU_FISSION:
58,352,868✔
470
    if (fissionable) {
58,352,868!
471
      if (dg != nullptr) {
58,352,868✔
472
        val = xs_t->delayed_nu_fission(a, *dg, gin);
50,016,744✔
473
      } else {
474
        val = 0.;
8,336,124✔
475
        for (int d = 0; d < xs_t->delayed_nu_fission.shape()[1]; d++) {
58,352,868✔
476
          val += xs_t->delayed_nu_fission(a, d, gin);
50,016,744✔
477
        }
478
      }
479
    } else {
UNCOV
480
      val = 0.;
×
481
    }
482
    break;
58,352,868✔
483
  case MgxsType::CHI_PROMPT:
7,038✔
484
    if (fissionable) {
7,038✔
485
      if (gout != nullptr) {
2,439!
486
        val = xs_t->chi_prompt(a, gin, *gout);
2,439✔
487
      } else {
488
        // provide an outgoing group-wise sum
UNCOV
489
        val = 0.;
×
UNCOV
490
        for (int g = 0; g < xs_t->chi_prompt.shape()[2]; g++) {
×
491
          val += xs_t->chi_prompt(a, gin, g);
×
492
        }
493
      }
494
    } else {
495
      val = 0.;
4,599✔
496
    }
497
    break;
7,038✔
UNCOV
498
  case MgxsType::CHI_DELAYED:
×
UNCOV
499
    if (fissionable) {
×
500
      if (gout != nullptr) {
×
501
        if (dg != nullptr) {
×
502
          val = xs_t->chi_delayed(a, *dg, gin, *gout);
×
503
        } else {
504
          val = xs_t->chi_delayed(a, 0, gin, *gout);
×
505
        }
506
      } else {
UNCOV
507
        if (dg != nullptr) {
×
UNCOV
508
          val = 0.;
×
509
          for (int g = 0; g < xs_t->delayed_nu_fission.shape()[2]; g++) {
×
510
            val += xs_t->delayed_nu_fission(a, *dg, gin, g);
×
511
          }
512
        } else {
UNCOV
513
          val = 0.;
×
UNCOV
514
          for (int g = 0; g < xs_t->delayed_nu_fission.shape()[2]; g++) {
×
515
            for (int d = 0; d < xs_t->delayed_nu_fission.shape()[3]; d++) {
×
516
              val += xs_t->delayed_nu_fission(a, d, gin, g);
×
517
            }
518
          }
519
        }
520
      }
521
    } else {
UNCOV
522
      val = 0.;
×
523
    }
524
    break;
×
525
  case MgxsType::INVERSE_VELOCITY:
1,701,385,524✔
526
    val = xs_t->inverse_velocity(a, gin);
1,701,385,524✔
527
    break;
1,701,385,524✔
528
  case MgxsType::DECAY_RATE:
50,019,381✔
529
    if (dg != nullptr) {
50,019,381!
530
      val = xs_t->decay_rate(a, *dg);
50,019,381✔
531
    } else {
UNCOV
532
      val = xs_t->decay_rate(a, 0);
×
533
    }
534
    break;
50,019,381✔
UNCOV
535
  default:
×
UNCOV
536
    val = 0.;
×
537
  }
538
  return val;
1,881,072,900✔
539
}
540

541
//==============================================================================
542

543
void Mgxs::sample_fission_energy(
90,891,765✔
544
  int gin, int& dg, int& gout, uint64_t* seed, int t, int a)
545
{
546
  XsData* xs_t = &xs[t];
90,891,765✔
547
  double nu_fission = xs_t->nu_fission(a, gin);
90,891,765✔
548

549
  // Find the probability of having a prompt neutron
550
  double prob_prompt = xs_t->prompt_nu_fission(a, gin);
90,891,765✔
551

552
  // sample random numbers
553
  double xi_pd = prn(seed) * nu_fission;
90,891,765✔
554
  double xi_gout = prn(seed);
90,891,765✔
555

556
  // Select whether the neutron is prompt or delayed
557
  if (xi_pd <= prob_prompt) {
90,891,765✔
558
    // the neutron is prompt
559

560
    // set the delayed group for the particle to be -1, indicating prompt
561
    dg = -1;
90,890,532✔
562

563
    // sample the outgoing energy group
564
    double prob_gout = 0.;
90,890,532✔
565
    for (gout = 0; gout < num_groups; ++gout) {
90,932,850!
566
      prob_gout += xs_t->chi_prompt(a, gin, gout);
90,932,850✔
567
      if (xi_gout < prob_gout)
90,932,850✔
568
        break;
90,890,532✔
569
    }
570

571
  } else {
572
    // the neutron is delayed
573

574
    // get the delayed group
575
    for (dg = 0; dg < num_delayed_groups; ++dg) {
3,105!
576
      prob_prompt += xs_t->delayed_nu_fission(a, dg, gin);
3,105✔
577
      if (xi_pd < prob_prompt)
3,105✔
578
        break;
1,233✔
579
    }
580

581
    // adjust dg in case of round-off error
582
    dg = std::min(dg, num_delayed_groups - 1);
1,233✔
583

584
    // sample the outgoing energy group
585
    double prob_gout = 0.;
1,233✔
586
    for (gout = 0; gout < num_groups; ++gout) {
1,233!
587
      prob_gout += xs_t->chi_delayed(a, dg, gin, gout);
1,233✔
588
      if (xi_gout < prob_gout)
1,233!
589
        break;
1,233✔
590
    }
591
  }
592
}
90,891,765✔
593

594
//==============================================================================
595

596
void Mgxs::sample_scatter(
1,365,885,576✔
597
  int gin, int& gout, double& mu, double& wgt, uint64_t* seed, int t, int a)
598
{
599
  // Sample the data
600
  xs[t].scatter[a]->sample(gin, gout, mu, wgt, seed);
1,365,885,576✔
601
}
1,365,885,576✔
602

603
//==============================================================================
604

605
void Mgxs::calculate_xs(Particle& p)
1,688,676,066✔
606
{
607
  // If the material is different, then we need to do a full lookup
608
  if (p.material() != p.mg_xs_cache().material) {
1,688,676,066✔
609
    set_temperature_index(p);
90,044,295✔
610
    set_angle_index(p);
90,044,295✔
611
    p.mg_xs_cache().material = p.material();
90,044,295✔
612
  } else {
613
    // If material is the same, but temperature is different, need to
614
    // find the new temperature index
615
    if (p.sqrtkT() != p.mg_xs_cache().sqrtkT) {
1,598,631,771✔
616
      set_temperature_index(p);
1,822,626✔
617
    }
618
    // If the material is the same, but angle is different, need to
619
    // find the new angle index
620
    if (p.u_local() != p.mg_xs_cache().u) {
1,598,631,771!
621
      set_angle_index(p);
1,598,631,771✔
622
    }
623
  }
624
  int temperature = p.mg_xs_cache().t;
1,688,676,066✔
625
  int angle = p.mg_xs_cache().a;
1,688,676,066✔
626
  p.macro_xs().total = xs[temperature].total(angle, p.g()) * p.density_mult();
1,688,676,066✔
627
  p.macro_xs().absorption =
2,147,483,647✔
628
    xs[temperature].absorption(angle, p.g()) * p.density_mult();
1,688,676,066✔
629
  p.macro_xs().nu_fission =
1,688,676,066✔
630
    fissionable ? xs[temperature].nu_fission(angle, p.g()) * p.density_mult()
1,688,676,066✔
631
                : 0.;
632
}
1,688,676,066✔
633

634
//==============================================================================
635

UNCOV
636
bool Mgxs::equiv(const Mgxs& that)
×
637
{
638
  return (
UNCOV
639
    (num_delayed_groups == that.num_delayed_groups) &&
×
UNCOV
640
    (num_groups == that.num_groups) && (n_pol == that.n_pol) &&
×
641
    (n_azi == that.n_azi) &&
×
642
    (std::equal(polar.begin(), polar.end(), that.polar.begin())) &&
×
643
    (std::equal(azimuthal.begin(), azimuthal.end(), that.azimuthal.begin())) &&
×
644
    (scatter_format == that.scatter_format));
×
645
}
646

647
//==============================================================================
648

649
int Mgxs::get_temperature_index(double sqrtkT) const
102,117,192✔
650
{
651
  return xt::argmin(xt::abs(kTs - sqrtkT * sqrtkT))[0];
102,117,192✔
652
}
653

654
//==============================================================================
655

656
void Mgxs::set_temperature_index(Particle& p)
91,866,921✔
657
{
658
  p.mg_xs_cache().t = get_temperature_index(p.sqrtkT());
91,866,921✔
659
  p.mg_xs_cache().sqrtkT = p.sqrtkT();
91,866,921✔
660
}
91,866,921✔
661

662
//==============================================================================
663

664
int Mgxs::get_angle_index(const Direction& u) const
1,746,660,987✔
665
{
666
  if (is_isotropic) {
1,746,660,987✔
667
    return 0;
1,746,624,789✔
668
  } else {
669
    // convert direction to polar and azimuthal angles
670
    double my_pol = std::acos(u.z);
36,198✔
671
    double my_azi = std::atan2(u.y, u.x);
36,198✔
672

673
    // Find the location, assuming equal-bin angles
674
    double delta_angle = PI / n_pol;
36,198✔
675
    int p = std::floor(my_pol / delta_angle);
36,198✔
676
    delta_angle = 2. * PI / n_azi;
36,198✔
677
    int a = std::floor((my_azi + PI) / delta_angle);
36,198✔
678

679
    return n_azi * p + a;
36,198✔
680
  }
681
}
682

683
//==============================================================================
684

685
void Mgxs::set_angle_index(Particle& p)
1,688,676,066✔
686
{
687
  // See if we need to find the new index
688
  if (!is_isotropic) {
1,688,676,066✔
689
    p.mg_xs_cache().a = get_angle_index(p.u_local());
18,099✔
690
    p.mg_xs_cache().u = p.u_local();
18,099✔
691
  }
692
}
1,688,676,066✔
693

694
} // namespace openmc
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2026 Coveralls, Inc