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

openmc-dev / openmc / 22070080915

16 Feb 2026 04:15PM UTC coverage: 81.702% (-0.02%) from 81.72%
22070080915

Pull #3806

github

web-flow
Merge 4824efac8 into c6ef84d1d
Pull Request #3806: Introduce new C API function for slice plots

17579 of 24709 branches covered (71.14%)

Branch coverage included in aggregate %.

197 of 245 new or added lines in 4 files covered. (80.41%)

336 existing lines in 26 files now uncovered.

57190 of 66805 relevant lines covered (85.61%)

50553675.8 hits per line

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

79.2
/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,
5,592✔
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;
5,592✔
35
  awr = in_awr;
5,592✔
36
  // TODO: Remove adapt when in_KTs is an xtensor
37
  kTs = xt::adapt(in_kTs);
5,592✔
38
  fissionable = in_fissionable;
5,592✔
39
  scatter_format = in_scatter_format;
5,592✔
40
  xs.resize(in_kTs.size());
5,592✔
41
  is_isotropic = in_is_isotropic;
5,592✔
42
  n_pol = in_polar.size();
5,592✔
43
  n_azi = in_azimuthal.size();
5,592✔
44
  polar = in_polar;
5,592✔
45
  azimuthal = in_azimuthal;
5,592✔
46
}
5,592✔
47

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

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

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

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

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

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

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

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

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

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

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

157
          if (std::find(temps_to_read.begin(), temps_to_read.end(),
168✔
158
                temps_available[j + 1]) == temps_to_read.end()) {
336!
159
            temps_to_read.push_back(temps_available[j + 1]);
168✔
160
          }
161
          break;
168✔
162
        }
163
      }
164
    }
165
  }
166
  std::sort(temps_to_read.begin(), temps_to_read.end());
2,985✔
167

168
  // Get the library's temperatures
169
  int n_temperature = temps_to_read.size();
2,985✔
170
  vector<double> in_kTs(n_temperature);
5,970✔
171
  for (int i = 0; i < n_temperature; i++) {
6,110✔
172
    std::string temp_str(std::to_string(temps_to_read[i]) + "K");
3,125✔
173

174
    // read exact temperature value
175
    read_double(kT_group, temp_str.c_str(), &in_kTs[i], true);
3,125✔
176
  }
3,125✔
177
  close_group(kT_group);
2,985✔
178

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

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

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

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

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

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

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

263
  // Set the angular bins to use equally-spaced bins
264
  vector<double> in_polar(in_n_pol);
5,970✔
265
  double dangle = PI / in_n_pol;
2,985✔
266
  for (int p = 0; p < in_n_pol; p++) {
5,998✔
267
    in_polar[p] = (p + 0.5) * dangle;
3,013✔
268
  }
269
  vector<double> in_azimuthal(in_n_azi);
5,970✔
270
  dangle = 2. * PI / in_n_azi;
2,985✔
271
  for (int a = 0; a < in_n_azi; a++) {
5,998✔
272
    in_azimuthal[a] = (a + 0.5) * dangle - PI;
3,013✔
273
  }
274

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

280
//==============================================================================
281

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

291
  // Set number of energy and delayed groups
292
  AngleDistributionType final_scatter_format = scatter_format;
2,985✔
293
  if (settings::legendre_to_tabular) {
2,985✔
294
    if (scatter_format == AngleDistributionType::LEGENDRE)
2,761✔
295
      final_scatter_format = AngleDistributionType::TABULAR;
2,677✔
296
  }
297

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

306
    xs[t].from_hdf5(xsdata_grp, fissionable, scatter_format,
3,125✔
307
      final_scatter_format, order_data, is_isotropic, n_pol, n_azi);
3,125✔
308
    close_group(xsdata_grp);
3,125✔
309

310
  } // end temperature loop
3,125✔
311

312
  // Make sure the scattering format is updated to the final case
313
  scatter_format = final_scatter_format;
2,985✔
314
}
2,985✔
315

316
//==============================================================================
317

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

340
  init(in_name, in_awr, mat_kTs, in_fissionable, in_scatter_format,
2,607✔
341
    in_is_isotropic, in_polar, in_azimuthal);
342

343
  // Create the xs data for each temperature
344
  for (int t = 0; t < mat_kTs.size(); t++) {
5,242✔
345
    xs[t] = XsData(in_fissionable, in_scatter_format, in_polar.size(),
5,270✔
346
      in_azimuthal.size(), num_groups, num_delayed_groups);
5,270✔
347

348
    // Find the right temperature index to use
349
    double temp_desired = mat_kTs[t];
2,635✔
350

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

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

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

419
    // And finally, combine the data
420
    combine(mgxs_to_combine, interpolant, temp_indices, t);
2,635✔
421
  } // end temperature (t) loop
2,635✔
422
}
2,607✔
423

424
//==============================================================================
425

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

435
  xs[this_t].combine(those_xs, scalars);
2,635✔
436
}
2,635✔
437

438
//==============================================================================
439

440
double Mgxs::get_xs(MgxsType xstype, int gin, const int* gout, const double* mu,
2,090,094,052✔
441
  const int* dg, int t, int a)
442
{
443
  XsData* xs_t = &xs[t];
2,090,094,052✔
444
  double val;
445
  switch (xstype) {
2,090,094,052!
446
  case MgxsType::TOTAL:
18,757,418✔
447
    val = xs_t->total(a, gin);
18,757,418✔
448
    break;
18,757,418✔
449
  case MgxsType::NU_FISSION:
9,270,566✔
450
    val = fissionable ? xs_t->nu_fission(a, gin) : 0.;
9,270,566✔
451
    break;
9,270,566✔
452
  case MgxsType::ABSORPTION:
5,296,100✔
453
    val = xs_t->absorption(a, gin);
5,296,100✔
454
    ;
455
    break;
5,296,100✔
456
  case MgxsType::FISSION:
10,154,966✔
457
    val = fissionable ? xs_t->fission(a, gin) : 0.;
10,154,966✔
458
    break;
10,154,966✔
459
  case MgxsType::KAPPA_FISSION:
9,492,206✔
460
    val = fissionable ? xs_t->kappa_fission(a, gin) : 0.;
9,492,206✔
461
    break;
9,492,206✔
462
  case MgxsType::NU_SCATTER:
17,010,260✔
463
  case MgxsType::SCATTER:
464
  case MgxsType::NU_SCATTER_FMU:
465
  case MgxsType::SCATTER_FMU:
466
    val = xs_t->scatter[a]->get_xs(xstype, gin, gout, mu);
17,010,260✔
467
    break;
17,010,260✔
468
  case MgxsType::PROMPT_NU_FISSION:
9,262,360✔
469
    val = fissionable ? xs_t->prompt_nu_fission(a, gin) : 0.;
9,262,360!
470
    break;
9,262,360✔
471
  case MgxsType::DELAYED_NU_FISSION:
64,836,520✔
472
    if (fissionable) {
64,836,520!
473
      if (dg != nullptr) {
64,836,520✔
474
        val = xs_t->delayed_nu_fission(a, *dg, gin);
55,574,160✔
475
      } else {
476
        val = 0.;
9,262,360✔
477
        for (int d = 0; d < xs_t->delayed_nu_fission.shape()[1]; d++) {
64,836,520✔
478
          val += xs_t->delayed_nu_fission(a, d, gin);
55,574,160✔
479
        }
480
      }
481
    } else {
UNCOV
482
      val = 0.;
×
483
    }
484
    break;
64,836,520✔
485
  case MgxsType::CHI_PROMPT:
8,206✔
486
    if (fissionable) {
8,206✔
487
      if (gout != nullptr) {
2,842!
488
        val = xs_t->chi_prompt(a, gin, *gout);
2,842✔
489
      } else {
490
        // provide an outgoing group-wise sum
UNCOV
491
        val = 0.;
×
492
        for (int g = 0; g < xs_t->chi_prompt.shape()[2]; g++) {
×
493
          val += xs_t->chi_prompt(a, gin, g);
×
494
        }
495
      }
496
    } else {
497
      val = 0.;
5,364✔
498
    }
499
    break;
8,206✔
UNCOV
500
  case MgxsType::CHI_DELAYED:
×
501
    if (fissionable) {
×
502
      if (gout != nullptr) {
×
503
        if (dg != nullptr) {
×
504
          val = xs_t->chi_delayed(a, *dg, gin, *gout);
×
505
        } else {
UNCOV
506
          val = xs_t->chi_delayed(a, 0, gin, *gout);
×
507
        }
508
      } else {
509
        if (dg != nullptr) {
×
510
          val = 0.;
×
UNCOV
511
          for (int g = 0; g < xs_t->delayed_nu_fission.shape()[2]; g++) {
×
UNCOV
512
            val += xs_t->delayed_nu_fission(a, *dg, gin, g);
×
513
          }
514
        } else {
UNCOV
515
          val = 0.;
×
516
          for (int g = 0; g < xs_t->delayed_nu_fission.shape()[2]; g++) {
×
UNCOV
517
            for (int d = 0; d < xs_t->delayed_nu_fission.shape()[3]; d++) {
×
518
              val += xs_t->delayed_nu_fission(a, d, gin, g);
×
519
            }
520
          }
521
        }
522
      }
523
    } else {
UNCOV
524
      val = 0.;
×
525
    }
526
    break;
×
527
  case MgxsType::INVERSE_VELOCITY:
1,890,428,360✔
528
    val = xs_t->inverse_velocity(a, gin);
1,890,428,360✔
529
    break;
1,890,428,360✔
530
  case MgxsType::DECAY_RATE:
55,577,090✔
531
    if (dg != nullptr) {
55,577,090!
532
      val = xs_t->decay_rate(a, *dg);
55,577,090✔
533
    } else {
UNCOV
534
      val = xs_t->decay_rate(a, 0);
×
535
    }
536
    break;
55,577,090✔
UNCOV
537
  default:
×
UNCOV
538
    val = 0.;
×
539
  }
540
  return val;
2,090,094,052✔
541
}
542

543
//==============================================================================
544

545
void Mgxs::sample_fission_energy(
100,990,850✔
546
  int gin, int& dg, int& gout, uint64_t* seed, int t, int a)
547
{
548
  XsData* xs_t = &xs[t];
100,990,850✔
549
  double nu_fission = xs_t->nu_fission(a, gin);
100,990,850✔
550

551
  // Find the probability of having a prompt neutron
552
  double prob_prompt = xs_t->prompt_nu_fission(a, gin);
100,990,850✔
553

554
  // sample random numbers
555
  double xi_pd = prn(seed) * nu_fission;
100,990,850✔
556
  double xi_gout = prn(seed);
100,990,850✔
557

558
  // Select whether the neutron is prompt or delayed
559
  if (xi_pd <= prob_prompt) {
100,990,850✔
560
    // the neutron is prompt
561

562
    // set the delayed group for the particle to be -1, indicating prompt
563
    dg = -1;
100,989,480✔
564

565
    // sample the outgoing energy group
566
    double prob_gout = 0.;
100,989,480✔
567
    for (gout = 0; gout < num_groups; ++gout) {
101,036,500!
568
      prob_gout += xs_t->chi_prompt(a, gin, gout);
101,036,500✔
569
      if (xi_gout < prob_gout)
101,036,500✔
570
        break;
100,989,480✔
571
    }
572

573
  } else {
574
    // the neutron is delayed
575

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

583
    // adjust dg in case of round-off error
584
    dg = std::min(dg, num_delayed_groups - 1);
1,370✔
585

586
    // sample the outgoing energy group
587
    double prob_gout = 0.;
1,370✔
588
    for (gout = 0; gout < num_groups; ++gout) {
1,370!
589
      prob_gout += xs_t->chi_delayed(a, dg, gin, gout);
1,370✔
590
      if (xi_gout < prob_gout)
1,370!
591
        break;
1,370✔
592
    }
593
  }
594
}
100,990,850✔
595

596
//==============================================================================
597

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

605
//==============================================================================
606

607
void Mgxs::calculate_xs(Particle& p)
1,876,306,740✔
608
{
609
  // If the material is different, then we need to do a full lookup
610
  if (p.material() != p.mg_xs_cache().material) {
1,876,306,740✔
611
    set_temperature_index(p);
101,161,459✔
612
    set_angle_index(p);
101,161,459✔
613
    p.mg_xs_cache().material = p.material();
101,161,459✔
614
  } else {
615
    // If material is the same, but temperature is different, need to
616
    // find the new temperature index
617
    if (p.sqrtkT() != p.mg_xs_cache().sqrtkT) {
1,775,145,281✔
618
      set_temperature_index(p);
2,024,841✔
619
    }
620
    // If the material is the same, but angle is different, need to
621
    // find the new angle index
622
    if (p.u_local() != p.mg_xs_cache().u) {
1,775,145,281!
623
      set_angle_index(p);
1,775,145,281✔
624
    }
625
  }
626
  int temperature = p.mg_xs_cache().t;
1,876,306,740✔
627
  int angle = p.mg_xs_cache().a;
1,876,306,740✔
628
  p.macro_xs().total = xs[temperature].total(angle, p.g()) * p.density_mult();
1,876,306,740✔
629
  p.macro_xs().absorption =
2,522,744,995✔
630
    xs[temperature].absorption(angle, p.g()) * p.density_mult();
1,876,306,740✔
631
  p.macro_xs().nu_fission =
1,876,306,740✔
632
    fissionable ? xs[temperature].nu_fission(angle, p.g()) * p.density_mult()
1,876,306,740✔
633
                : 0.;
634
}
1,876,306,740✔
635

636
//==============================================================================
637

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

649
//==============================================================================
650

651
int Mgxs::get_temperature_index(double sqrtkT) const
114,575,513✔
652
{
653
  return xt::argmin(xt::abs(kTs - sqrtkT * sqrtkT))[0];
114,575,513✔
654
}
655

656
//==============================================================================
657

658
void Mgxs::set_temperature_index(Particle& p)
103,186,300✔
659
{
660
  p.mg_xs_cache().t = get_temperature_index(p.sqrtkT());
103,186,300✔
661
  p.mg_xs_cache().sqrtkT = p.sqrtkT();
103,186,300✔
662
}
103,186,300✔
663

664
//==============================================================================
665

666
int Mgxs::get_angle_index(const Direction& u) const
1,940,734,430✔
667
{
668
  if (is_isotropic) {
1,940,734,430✔
669
    return 0;
1,940,694,210✔
670
  } else {
671
    // convert direction to polar and azimuthal angles
672
    double my_pol = std::acos(u.z);
40,220✔
673
    double my_azi = std::atan2(u.y, u.x);
40,220✔
674

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

681
    return n_azi * p + a;
40,220✔
682
  }
683
}
684

685
//==============================================================================
686

687
void Mgxs::set_angle_index(Particle& p)
1,876,306,740✔
688
{
689
  // See if we need to find the new index
690
  if (!is_isotropic) {
1,876,306,740✔
691
    p.mg_xs_cache().a = get_angle_index(p.u_local());
20,110✔
692
    p.mg_xs_cache().u = p.u_local();
20,110✔
693
  }
694
}
1,876,306,740✔
695

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