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

openmc-dev / openmc / 21077426970

16 Jan 2026 06:56PM UTC coverage: 82.029% (-0.04%) from 82.069%
21077426970

Pull #3740

github

web-flow
Merge ea498539e into 51ea89ccc
Pull Request #3740: Add `chi` cross sections to MGXS machinery

17267 of 23987 branches covered (71.98%)

Branch coverage included in aggregate %.

34 of 44 new or added lines in 3 files covered. (77.27%)

7 existing lines in 2 files now uncovered.

55664 of 64922 relevant lines covered (85.74%)

43515974.52 hits per line

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

93.68
/src/xsdata.cpp
1
#include "openmc/xsdata.h"
2

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

8
#include "xtensor/xbuilder.hpp"
9
#include "xtensor/xindex_view.hpp"
10
#include "xtensor/xmath.hpp"
11
#include "xtensor/xview.hpp"
12

13
#include "openmc/constants.h"
14
#include "openmc/error.h"
15
#include "openmc/math_functions.h"
16
#include "openmc/mgxs_interface.h"
17
#include "openmc/random_lcg.h"
18
#include "openmc/settings.h"
19

20
namespace openmc {
21

22
//==============================================================================
23
// XsData class methods
24
//==============================================================================
25

26
XsData::XsData(bool fissionable, AngleDistributionType scatter_format,
5,814✔
27
  int n_pol, int n_azi, size_t n_groups, size_t n_d_groups)
5,814✔
28
  : n_g_(n_groups), n_dg_(n_d_groups)
5,814✔
29
{
30
  size_t n_ang = n_pol * n_azi;
5,814✔
31

32
  // check to make sure scatter format is OK before we allocate
33
  if (scatter_format != AngleDistributionType::HISTOGRAM &&
5,814✔
34
      scatter_format != AngleDistributionType::TABULAR &&
448!
35
      scatter_format != AngleDistributionType::LEGENDRE) {
36
    fatal_error("Invalid scatter_format!");
×
37
  }
38
  // allocate all [temperature][angle][in group] quantities
39
  vector<size_t> shape {n_ang, n_g_};
5,814✔
40
  total = xt::zeros<double>(shape);
5,814✔
41
  absorption = xt::zeros<double>(shape);
5,814✔
42
  inverse_velocity = xt::zeros<double>(shape);
5,814✔
43
  if (fissionable) {
5,814✔
44
    fission = xt::zeros<double>(shape);
2,160✔
45
    nu_fission = xt::zeros<double>(shape);
2,160✔
46
    prompt_nu_fission = xt::zeros<double>(shape);
2,160✔
47
    kappa_fission = xt::zeros<double>(shape);
2,160✔
48
  }
49

50
  // allocate decay_rate; [temperature][angle][delayed group]
51
  shape[1] = n_dg_;
5,814✔
52
  decay_rate = xt::zeros<double>(shape);
5,814✔
53

54
  if (fissionable) {
5,814✔
55
    shape = {n_ang, n_dg_, n_g_};
2,160✔
56
    // allocate delayed_nu_fission; [temperature][angle][delay group][in group]
57
    delayed_nu_fission = xt::zeros<double>(shape);
2,160✔
58

59
    // chi and chi_prompt; [temperature][angle][in group][out group]
60
    shape = {n_ang, n_g_, n_g_};
2,160✔
61
    chi = xt::zeros<double>(shape);
2,160✔
62
    chi_prompt = xt::zeros<double>(shape);
2,160✔
63

64
    // chi_delayed; [temperature][angle][delay group][in group][out group]
65
    shape = {n_ang, n_dg_, n_g_, n_g_};
2,160✔
66
    chi_delayed = xt::zeros<double>(shape);
2,160✔
67
  }
68

69
  for (int a = 0; a < n_ang; a++) {
11,820✔
70
    if (scatter_format == AngleDistributionType::HISTOGRAM) {
6,006✔
71
      scatter.emplace_back(new ScattDataHistogram);
128✔
72
    } else if (scatter_format == AngleDistributionType::TABULAR) {
5,878✔
73
      scatter.emplace_back(new ScattDataTabular);
5,334✔
74
    } else if (scatter_format == AngleDistributionType::LEGENDRE) {
544!
75
      scatter.emplace_back(new ScattDataLegendre);
544✔
76
    }
77
  }
78
}
5,814✔
79

80
//==============================================================================
81

82
void XsData::from_hdf5(hid_t xsdata_grp, bool fissionable,
3,187✔
83
  AngleDistributionType scatter_format,
84
  AngleDistributionType final_scatter_format, int order_data, bool is_isotropic,
85
  int n_pol, int n_azi)
86
{
87
  // Reconstruct the dimension information so it doesn't need to be passed
88
  size_t n_ang = n_pol * n_azi;
3,187✔
89
  size_t energy_groups = total.shape()[1];
3,187✔
90

91
  // Set the fissionable-specific data
92
  if (fissionable) {
3,187✔
93
    fission_from_hdf5(xsdata_grp, n_ang, is_isotropic);
1,168✔
94
  }
95
  // Get the non-fission-specific data
96
  read_nd_vector(xsdata_grp, "decay-rate", decay_rate);
3,187✔
97
  read_nd_vector(xsdata_grp, "absorption", absorption, true);
3,187✔
98
  read_nd_vector(xsdata_grp, "inverse-velocity", inverse_velocity);
3,187✔
99

100
  // Get scattering data
101
  scatter_from_hdf5(
3,187✔
102
    xsdata_grp, n_ang, scatter_format, final_scatter_format, order_data);
103

104
  // Check absorption to ensure it is not 0 since it is often the
105
  // denominator in tally methods
106
  xt::filtration(absorption, xt::equal(absorption, 0.)) = 1.e-10;
3,187✔
107

108
  // Get or calculate the total x/s
109
  if (object_exists(xsdata_grp, "total")) {
3,187!
110
    read_nd_vector(xsdata_grp, "total", total);
3,187✔
111
  } else {
112
    for (size_t a = 0; a < n_ang; a++) {
×
113
      for (size_t gin = 0; gin < energy_groups; gin++) {
×
114
        total(a, gin) = absorption(a, gin) + scatter[a]->scattxs[gin];
×
115
      }
116
    }
117
  }
118

119
  // Fix if total is 0, since it is in the denominator when tallying
120
  xt::filtration(total, xt::equal(total, 0.)) = 1.e-10;
3,187✔
121
}
3,187✔
122

123
//==============================================================================
124

125
void XsData::fission_vector_beta_from_hdf5(
48✔
126
  hid_t xsdata_grp, size_t n_ang, bool is_isotropic)
127
{
128
  // Data is provided as nu-fission and chi with a beta for delayed info
129

130
  // Get chi
131
  xt::xtensor<double, 2> temp_chi({n_ang, n_g_}, 0.);
48✔
132
  read_nd_vector(xsdata_grp, "chi", temp_chi, true);
48✔
133

134
  // Normalize chi by summing over the outgoing groups for each incoming angle
135
  temp_chi /= xt::view(xt::sum(temp_chi, {1}), xt::all(), xt::newaxis());
48✔
136

137
  // Store chi in chi
138
  chi = xt::view(temp_chi, xt::all(), xt::newaxis(), xt::all());
48✔
139

140
  // Now every incoming group in prompt_chi and delayed_chi is the normalized
141
  // chi we just made
142
  // TODO: This is incorrect. This makes chi_prompt and chi_delayed identical
143
  // to chi.
144
  chi_prompt = xt::view(temp_chi, xt::all(), xt::newaxis(), xt::all());
48✔
145
  chi_delayed =
146
    xt::view(temp_chi, xt::all(), xt::newaxis(), xt::newaxis(), xt::all());
48✔
147

148
  // Get nu-fission
149
  xt::xtensor<double, 2> temp_nufiss({n_ang, n_g_}, 0.);
48✔
150
  read_nd_vector(xsdata_grp, "nu-fission", temp_nufiss, true);
48✔
151

152
  // Get beta (strategy will depend upon the number of dimensions in beta)
153
  hid_t beta_dset = open_dataset(xsdata_grp, "beta");
48✔
154
  int beta_ndims = dataset_ndims(beta_dset);
48✔
155
  close_dataset(beta_dset);
48✔
156
  int ndim_target = 1;
48✔
157
  if (!is_isotropic)
48!
158
    ndim_target += 2;
×
159
  if (beta_ndims == ndim_target) {
48✔
160
    xt::xtensor<double, 2> temp_beta({n_ang, n_dg_}, 0.);
32✔
161
    read_nd_vector(xsdata_grp, "beta", temp_beta, true);
32✔
162

163
    // Set prompt_nu_fission = (1. - beta_total)*nu_fission
164
    prompt_nu_fission = temp_nufiss * (1. - xt::sum(temp_beta, {1}));
32✔
165

166
    // Set delayed_nu_fission as beta * nu_fission
167
    delayed_nu_fission =
168
      xt::view(temp_beta, xt::all(), xt::all(), xt::newaxis()) *
64✔
169
      xt::view(temp_nufiss, xt::all(), xt::newaxis(), xt::all());
96✔
170
  } else if (beta_ndims == ndim_target + 1) {
48!
171
    xt::xtensor<double, 3> temp_beta({n_ang, n_dg_, n_g_}, 0.);
16✔
172
    read_nd_vector(xsdata_grp, "beta", temp_beta, true);
16✔
173

174
    // Set prompt_nu_fission = (1. - beta_total)*nu_fission
175
    prompt_nu_fission = temp_nufiss * (1. - xt::sum(temp_beta, {1}));
16✔
176

177
    // Set delayed_nu_fission as beta * nu_fission
178
    delayed_nu_fission =
179
      temp_beta * xt::view(temp_nufiss, xt::all(), xt::newaxis(), xt::all());
16✔
180
  }
16✔
181
}
48✔
182

183
void XsData::fission_vector_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang)
16✔
184
{
185
  // If chi is included in the dataset, we should store it!
186
  if (object_exists(xsdata_grp, "chi")) {
16!
NEW
187
    xt::xtensor<double, 2> temp_chi({n_ang, n_g_}, 0.);
×
NEW
188
    read_nd_vector(xsdata_grp, "chi", temp_chi, true);
×
189
    // Normalize chi by summing over the outgoing groups for each incoming angle
NEW
190
    temp_chi /= xt::view(xt::sum(temp_chi, {1}), xt::all(), xt::newaxis());
×
191
    // Now every incoming group in self.chi is the normalized chi we just made
NEW
192
    chi = xt::view(temp_chi, xt::all(), xt::newaxis(), xt::all());
×
NEW
193
  }
×
194

195
  // Data is provided separately as prompt + delayed nu-fission and chi
196
  // Get chi-prompt
197
  xt::xtensor<double, 2> temp_chi_p({n_ang, n_g_}, 0.);
16✔
198
  read_nd_vector(xsdata_grp, "chi-prompt", temp_chi_p, true);
16✔
199

200
  // Normalize chi-prompt by summing over the outgoing groups for each incoming
201
  // angle
202
  temp_chi_p /= xt::view(xt::sum(temp_chi_p, {1}), xt::all(), xt::newaxis());
16✔
203

204
  // Get chi-delayed
205
  xt::xtensor<double, 3> temp_chi_d({n_ang, n_dg_, n_g_}, 0.);
16✔
206
  read_nd_vector(xsdata_grp, "chi-delayed", temp_chi_d, true);
16✔
207

208
  // Normalize chi-delayed by summing over the outgoing groups for each incoming
209
  // angle
210
  temp_chi_d /=
211
    xt::view(xt::sum(temp_chi_d, {2}), xt::all(), xt::all(), xt::newaxis());
16✔
212

213
  // Now assign the prompt and delayed chis by replicating for each incoming
214
  // group
215
  chi_prompt = xt::view(temp_chi_p, xt::all(), xt::newaxis(), xt::all());
16✔
216
  chi_delayed =
217
    xt::view(temp_chi_d, xt::all(), xt::all(), xt::newaxis(), xt::all());
16✔
218

219
  // Get prompt and delayed nu-fission directly
220
  read_nd_vector(xsdata_grp, "prompt-nu-fission", prompt_nu_fission, true);
16✔
221
  read_nd_vector(xsdata_grp, "delayed-nu-fission", delayed_nu_fission, true);
16✔
222
}
16✔
223

224
void XsData::fission_vector_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang)
992✔
225
{
226
  // No beta is provided and there is no prompt/delay distinction.
227
  // Therefore, the code only considers the data as prompt.
228

229
  // Get chi
230
  xt::xtensor<double, 2> temp_chi({n_ang, n_g_}, 0.);
992✔
231
  read_nd_vector(xsdata_grp, "chi", temp_chi, true);
992✔
232

233
  // Normalize chi by summing over the outgoing groups for each incoming angle
234
  temp_chi /= xt::view(xt::sum(temp_chi, {1}), xt::all(), xt::newaxis());
992✔
235

236
  // Now every incoming group in self.chi is the normalized chi we just made
237
  chi = xt::view(temp_chi, xt::all(), xt::newaxis(), xt::all());
992✔
238
  chi_prompt = xt::view(temp_chi, xt::all(), xt::newaxis(), xt::all());
992✔
239

240
  // Get nu-fission directly
241
  read_nd_vector(xsdata_grp, "nu-fission", prompt_nu_fission, true);
992✔
242
}
992✔
243

244
//==============================================================================
245

246
void XsData::fission_matrix_beta_from_hdf5(
32✔
247
  hid_t xsdata_grp, size_t n_ang, bool is_isotropic)
248
{
249
  // Data is provided as nu-fission and chi with a beta for delayed info
250

251
  // Get nu-fission matrix
252
  xt::xtensor<double, 3> temp_matrix({n_ang, n_g_, n_g_}, 0.);
32✔
253
  read_nd_vector(xsdata_grp, "nu-fission", temp_matrix, true);
32✔
254

255
  // Get beta (strategy will depend upon the number of dimensions in beta)
256
  hid_t beta_dset = open_dataset(xsdata_grp, "beta");
32✔
257
  int beta_ndims = dataset_ndims(beta_dset);
32✔
258
  close_dataset(beta_dset);
32✔
259
  int ndim_target = 1;
32✔
260
  if (!is_isotropic)
32!
261
    ndim_target += 2;
×
262
  if (beta_ndims == ndim_target) {
32✔
263
    xt::xtensor<double, 2> temp_beta({n_ang, n_dg_}, 0.);
16✔
264
    read_nd_vector(xsdata_grp, "beta", temp_beta, true);
16✔
265

266
    xt::xtensor<double, 1> temp_beta_sum({n_ang}, 0.);
16✔
267
    temp_beta_sum = xt::sum(temp_beta, {1});
16✔
268

269
    // prompt_nu_fission is the sum of this matrix over outgoing groups and
270
    // multiplied by (1 - beta_sum)
271
    prompt_nu_fission = xt::sum(temp_matrix, {2}) * (1. - temp_beta_sum);
16✔
272

273
    // Store chi-prompt
274
    chi_prompt =
275
      xt::view(1.0 - temp_beta_sum, xt::all(), xt::newaxis(), xt::newaxis()) *
32✔
276
      temp_matrix;
16✔
277

278
    // delayed_nu_fission is the sum of this matrix over outgoing groups and
279
    // multiplied by beta
280
    delayed_nu_fission =
281
      xt::view(temp_beta, xt::all(), xt::all(), xt::newaxis()) *
32✔
282
      xt::view(xt::sum(temp_matrix, {2}), xt::all(), xt::newaxis(), xt::all());
48✔
283

284
    // Store chi-delayed
285
    chi_delayed =
286
      xt::view(temp_beta, xt::all(), xt::all(), xt::newaxis(), xt::newaxis()) *
32✔
287
      xt::view(temp_matrix, xt::all(), xt::newaxis(), xt::all(), xt::all());
48✔
288

289
    // Store chi
290
    chi =
291
      chi_prompt * (1. - temp_beta_sum) + xt::sum(temp_beta * chi_delayed, {1});
16✔
292

293
  } else if (beta_ndims == ndim_target + 1) {
32!
294
    xt::xtensor<double, 3> temp_beta({n_ang, n_dg_, n_g_}, 0.);
16✔
295
    read_nd_vector(xsdata_grp, "beta", temp_beta, true);
16✔
296

297
    xt::xtensor<double, 2> temp_beta_sum({n_ang, n_g_}, 0.);
16✔
298
    temp_beta_sum = xt::sum(temp_beta, {1});
16✔
299

300
    // prompt_nu_fission is the sum of this matrix over outgoing groups and
301
    // multiplied by (1 - beta_sum)
302
    prompt_nu_fission = xt::sum(temp_matrix, {2}) * (1. - temp_beta_sum);
16✔
303

304
    // Store chi-prompt
305
    chi_prompt =
306
      xt::view(1.0 - temp_beta_sum, xt::all(), xt::all(), xt::newaxis()) *
32✔
307
      temp_matrix;
16✔
308

309
    // delayed_nu_fission is the sum of this matrix over outgoing groups and
310
    // multiplied by beta
311
    delayed_nu_fission = temp_beta * xt::view(xt::sum(temp_matrix, {2}),
32✔
312
                                       xt::all(), xt::newaxis(), xt::all());
32✔
313

314
    // Store chi-delayed
315
    chi_delayed =
316
      xt::view(temp_beta, xt::all(), xt::all(), xt::all(), xt::newaxis()) *
32✔
317
      xt::view(temp_matrix, xt::all(), xt::newaxis(), xt::all(), xt::all());
48✔
318

319
    // Store chi
320
    chi =
321
      chi_prompt * (1. - temp_beta_sum) + xt::sum(temp_beta * chi_delayed, {1});
16✔
322
  }
16✔
323

324
  // Normalize chis
325
  chi_prompt /=
326
    xt::view(xt::sum(chi_prompt, {2}), xt::all(), xt::all(), xt::newaxis());
32✔
327

328
  chi_delayed /= xt::view(
64✔
329
    xt::sum(chi_delayed, {3}), xt::all(), xt::all(), xt::all(), xt::newaxis());
96✔
330

331
  chi /= xt::view(xt::sum(chi, {2}), xt::all(), xt::all(), xt::newaxis());
32✔
332
}
32✔
333

334
void XsData::fission_matrix_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang)
16✔
335
{
336
  // Data is provided separately as prompt + delayed nu-fission and chi
337

338
  // Get the prompt nu-fission matrix
339
  xt::xtensor<double, 3> temp_matrix_p({n_ang, n_g_, n_g_}, 0.);
16✔
340
  read_nd_vector(xsdata_grp, "prompt-nu-fission", temp_matrix_p, true);
16✔
341

342
  // prompt_nu_fission is the sum over outgoing groups
343
  prompt_nu_fission = xt::sum(temp_matrix_p, {2});
16✔
344

345
  // chi_prompt is this matrix but normalized over outgoing groups, which we
346
  // have already stored in prompt_nu_fission
347
  chi_prompt = temp_matrix_p /
32✔
348
               xt::view(prompt_nu_fission, xt::all(), xt::all(), xt::newaxis());
48✔
349

350
  // Get the delayed nu-fission matrix
351
  xt::xtensor<double, 4> temp_matrix_d({n_ang, n_dg_, n_g_, n_g_}, 0.);
16✔
352
  read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_matrix_d, true);
16✔
353

354
  // delayed_nu_fission is the sum over outgoing groups
355
  delayed_nu_fission = xt::sum(temp_matrix_d, {3});
16✔
356

357
  // chi_delayee is this matrix but normalized over outgoing groups, which we
358
  // have already stored in delayed_nu_fission
359
  chi_delayed = temp_matrix_d / xt::view(delayed_nu_fission, xt::all(),
32✔
360
                                  xt::all(), xt::all(), xt::newaxis());
32✔
361
}
16✔
362

363
void XsData::fission_matrix_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang)
64✔
364
{
365
  // No beta is provided and there is no prompt/delay distinction.
366
  // Therefore, the code only considers the data as prompt.
367

368
  // Get nu-fission matrix
369
  xt::xtensor<double, 3> temp_matrix({n_ang, n_g_, n_g_}, 0.);
64✔
370
  read_nd_vector(xsdata_grp, "nu-fission", temp_matrix, true);
64✔
371

372
  // prompt_nu_fission is the sum over outgoing groups
373
  prompt_nu_fission = xt::sum(temp_matrix, {2});
64✔
374

375
  // chi_prompt is this matrix but normalized over outgoing groups, which we
376
  // have already stored in prompt_nu_fission
377
  chi = temp_matrix /
128✔
378
        xt::view(prompt_nu_fission, xt::all(), xt::all(), xt::newaxis());
192✔
379
  chi_prompt = temp_matrix /
128✔
380
               xt::view(prompt_nu_fission, xt::all(), xt::all(), xt::newaxis());
192✔
381
}
64✔
382

383
//==============================================================================
384

385
void XsData::fission_from_hdf5(
1,168✔
386
  hid_t xsdata_grp, size_t n_ang, bool is_isotropic)
387
{
388
  // Get the fission and kappa_fission data xs; these are optional
389
  read_nd_vector(xsdata_grp, "fission", fission);
1,168✔
390
  read_nd_vector(xsdata_grp, "kappa-fission", kappa_fission);
1,168✔
391

392
  // Get the data; the strategy for doing so depends on if the data is provided
393
  // as a nu-fission matrix or a set of chi and nu-fission vectors
394
  if (object_exists(xsdata_grp, "chi") ||
1,296✔
395
      object_exists(xsdata_grp, "chi-prompt")) {
128✔
396
    if (n_dg_ == 0) {
1,056✔
397
      fission_vector_no_delayed_from_hdf5(xsdata_grp, n_ang);
992✔
398
    } else {
399
      if (object_exists(xsdata_grp, "beta")) {
64✔
400
        fission_vector_beta_from_hdf5(xsdata_grp, n_ang, is_isotropic);
48✔
401
      } else {
402
        fission_vector_no_beta_from_hdf5(xsdata_grp, n_ang);
16✔
403
      }
404
    }
405
  } else {
406
    if (n_dg_ == 0) {
112✔
407
      fission_matrix_no_delayed_from_hdf5(xsdata_grp, n_ang);
64✔
408
    } else {
409
      if (object_exists(xsdata_grp, "beta")) {
48✔
410
        fission_matrix_beta_from_hdf5(xsdata_grp, n_ang, is_isotropic);
32✔
411
      } else {
412
        fission_matrix_no_beta_from_hdf5(xsdata_grp, n_ang);
16✔
413
      }
414
    }
415
  }
416

417
  // Combine prompt_nu_fission and delayed_nu_fission into nu_fission
418
  if (n_dg_ == 0) {
1,168✔
419
    nu_fission = prompt_nu_fission;
1,056✔
420
  } else {
421
    nu_fission = prompt_nu_fission + xt::sum(delayed_nu_fission, {1});
112✔
422
  }
423
}
1,168✔
424

425
//==============================================================================
426

427
void XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang,
3,187✔
428
  AngleDistributionType scatter_format,
429
  AngleDistributionType final_scatter_format, int order_data)
430
{
431
  if (!object_exists(xsdata_grp, "scatter_data")) {
3,187!
432
    fatal_error("Must provide scatter_data group!");
×
433
  }
434
  hid_t scatt_grp = open_group(xsdata_grp, "scatter_data");
3,187✔
435

436
  // Get the outgoing group boundary indices
437
  xt::xtensor<int, 2> gmin({n_ang, n_g_}, 0.);
3,187✔
438
  read_nd_vector(scatt_grp, "g_min", gmin, true);
3,187✔
439
  xt::xtensor<int, 2> gmax({n_ang, n_g_}, 0.);
3,187✔
440
  read_nd_vector(scatt_grp, "g_max", gmax, true);
3,187✔
441

442
  // Make gmin and gmax start from 0 vice 1 as they do in the library
443
  gmin -= 1;
3,187✔
444
  gmax -= 1;
3,187✔
445

446
  // Now use this info to find the length of a vector to hold the flattened
447
  // data.
448
  size_t length = order_data * xt::sum(gmax - gmin + 1)();
3,187✔
449

450
  double_4dvec input_scatt(n_ang, double_3dvec(n_g_));
6,374✔
451
  xt::xtensor<double, 1> temp_arr({length}, 0.);
3,187✔
452
  read_nd_vector(scatt_grp, "scatter_matrix", temp_arr, true);
3,187✔
453

454
  // Compare the number of orders given with the max order of the problem;
455
  // strip off the superfluous orders if needed
456
  int order_dim;
457
  if (scatter_format == AngleDistributionType::LEGENDRE) {
3,187✔
458
    order_dim = std::min(order_data - 1, settings::max_order) + 1;
3,059✔
459
  } else {
460
    order_dim = order_data;
128✔
461
  }
462

463
  // convert the flattened temp_arr to a jagged array for passing to
464
  // scatt data
465
  size_t temp_idx = 0;
3,187✔
466
  for (size_t a = 0; a < n_ang; a++) {
6,470✔
467
    for (size_t gin = 0; gin < n_g_; gin++) {
14,695✔
468
      input_scatt[a][gin].resize(gmax(a, gin) - gmin(a, gin) + 1);
11,412✔
469
      for (size_t i_gout = 0; i_gout < input_scatt[a][gin].size(); i_gout++) {
61,946✔
470
        input_scatt[a][gin][i_gout].resize(order_dim);
50,534✔
471
        for (size_t l = 0; l < order_dim; l++) {
111,772✔
472
          input_scatt[a][gin][i_gout][l] = temp_arr[temp_idx++];
61,238✔
473
        }
474
        // Adjust index for the orders we didnt take
475
        temp_idx += (order_data - order_dim);
50,534✔
476
      }
477
    }
478
  }
479

480
  // Get multiplication matrix
481
  double_3dvec temp_mult(n_ang, double_2dvec(n_g_));
6,374✔
482
  if (object_exists(scatt_grp, "multiplicity_matrix")) {
3,187✔
483
    temp_arr.resize({length / order_data});
481✔
484
    read_nd_vector(scatt_grp, "multiplicity_matrix", temp_arr);
481✔
485

486
    // convert the flat temp_arr to a jagged array for passing to scatt data
487
    size_t temp_idx = 0;
481✔
488
    for (size_t a = 0; a < n_ang; a++) {
962✔
489
      for (size_t gin = 0; gin < n_g_; gin++) {
5,859✔
490
        temp_mult[a][gin].resize(gmax(a, gin) - gmin(a, gin) + 1);
5,378✔
491
        for (size_t i_gout = 0; i_gout < temp_mult[a][gin].size(); i_gout++) {
43,222✔
492
          temp_mult[a][gin][i_gout] = temp_arr[temp_idx++];
37,844✔
493
        }
494
      }
495
    }
496
  } else {
497
    // Use a default: multiplicities are 1.0.
498
    for (size_t a = 0; a < n_ang; a++) {
5,508✔
499
      for (size_t gin = 0; gin < n_g_; gin++) {
8,836✔
500
        temp_mult[a][gin].resize(gmax(a, gin) - gmin(a, gin) + 1);
6,034✔
501
        for (size_t i_gout = 0; i_gout < temp_mult[a][gin].size(); i_gout++) {
18,724✔
502
          temp_mult[a][gin][i_gout] = 1.;
12,690✔
503
        }
504
      }
505
    }
506
  }
507
  close_group(scatt_grp);
3,187✔
508

509
  // Finally, convert the Legendre data to tabular, if needed
510
  if (scatter_format == AngleDistributionType::LEGENDRE &&
3,187✔
511
      final_scatter_format == AngleDistributionType::TABULAR) {
512
    for (size_t a = 0; a < n_ang; a++) {
5,718✔
513
      ScattDataLegendre legendre_scatt;
2,883✔
514
      xt::xtensor<int, 1> in_gmin = xt::view(gmin, a, xt::all());
2,883✔
515
      xt::xtensor<int, 1> in_gmax = xt::view(gmax, a, xt::all());
2,883✔
516

517
      legendre_scatt.init(in_gmin, in_gmax, temp_mult[a], input_scatt[a]);
2,883✔
518

519
      // Now create a tabular version of legendre_scatt
520
      convert_legendre_to_tabular(
2,883✔
521
        legendre_scatt, *static_cast<ScattDataTabular*>(scatter[a].get()));
2,883✔
522

523
      scatter_format = final_scatter_format;
2,883✔
524
    }
2,883✔
525
  } else {
2,835✔
526
    // We are sticking with the current representation
527
    // Initialize the ScattData object with this data
528
    for (size_t a = 0; a < n_ang; a++) {
752✔
529
      xt::xtensor<int, 1> in_gmin = xt::view(gmin, a, xt::all());
400✔
530
      xt::xtensor<int, 1> in_gmax = xt::view(gmax, a, xt::all());
400✔
531
      scatter[a]->init(in_gmin, in_gmax, temp_mult[a], input_scatt[a]);
400✔
532
    }
400✔
533
  }
534
}
3,187✔
535

536
//==============================================================================
537

538
void XsData::combine(
2,627✔
539
  const vector<XsData*>& those_xs, const vector<double>& scalars)
540
{
541
  // Combine the non-scattering data
542
  for (size_t i = 0; i < those_xs.size(); i++) {
5,830✔
543
    XsData* that = those_xs[i];
3,203✔
544
    if (!equiv(*that))
3,203!
545
      fatal_error("Cannot combine the XsData objects!");
×
546
    double scalar = scalars[i];
3,203✔
547
    total += scalar * that->total;
3,203✔
548
    absorption += scalar * that->absorption;
3,203✔
549
    if (i == 0) {
3,203✔
550
      inverse_velocity = that->inverse_velocity;
2,627✔
551
    }
552
    if (that->prompt_nu_fission.shape()[0] > 0) {
3,203✔
553
      nu_fission += scalar * that->nu_fission;
1,184✔
554
      prompt_nu_fission += scalar * that->prompt_nu_fission;
1,184✔
555
      kappa_fission += scalar * that->kappa_fission;
1,184✔
556
      fission += scalar * that->fission;
1,184✔
557
      delayed_nu_fission += scalar * that->delayed_nu_fission;
1,184✔
558
      chi += scalar * that->chi;
1,184✔
559
      chi_prompt += scalar *
1,184✔
560
                    xt::view(xt::sum(that->prompt_nu_fission, {1}), xt::all(),
2,368✔
561
                      xt::newaxis(), xt::newaxis()) *
3,552✔
562
                    that->chi_prompt;
2,368✔
563
      chi_delayed += scalar *
1,184✔
564
                     xt::view(xt::sum(that->delayed_nu_fission, {2}), xt::all(),
2,368✔
565
                       xt::all(), xt::newaxis(), xt::newaxis()) *
4,736✔
566
                     that->chi_delayed;
2,368✔
567
    }
568
    decay_rate += scalar * that->decay_rate;
3,203✔
569
  }
570

571
  // Ensure chi, chi_prompt and chi_delayed are normalized to 1 for each
572
  // azimuthal angle and delayed group (for chi_delayed)
573
  chi /= xt::view(xt::sum(chi, {2}), xt::all(), xt::all(), xt::newaxis());
2,627✔
574
  chi_prompt /=
575
    xt::view(xt::sum(chi_prompt, {2}), xt::all(), xt::all(), xt::newaxis());
2,627✔
576
  chi_delayed /= xt::view(
5,254✔
577
    xt::sum(chi_delayed, {3}), xt::all(), xt::all(), xt::all(), xt::newaxis());
7,881✔
578

579
  // Allow the ScattData object to combine itself
580
  for (size_t a = 0; a < total.shape()[0]; a++) {
5,350✔
581
    // Build vector of the scattering objects to incorporate
582
    vector<ScattData*> those_scatts(those_xs.size());
2,723✔
583
    for (size_t i = 0; i < those_xs.size(); i++) {
6,022✔
584
      those_scatts[i] = those_xs[i]->scatter[a].get();
3,299✔
585
    }
586

587
    // Now combine these guys
588
    scatter[a]->combine(those_scatts, scalars);
2,723✔
589
  }
2,723✔
590
}
2,627✔
591

592
//==============================================================================
593

594
bool XsData::equiv(const XsData& that)
3,203✔
595
{
596
  return (absorption.shape() == that.absorption.shape());
3,203✔
597
}
598

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