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

openmc-dev / openmc / 13461947318

21 Feb 2025 05:24PM UTC coverage: 85.019% (+0.05%) from 84.969%
13461947318

Pull #3140

github

web-flow
Merge ece247ada into 2b788ea6e
Pull Request #3140: Add Versioning Support from `version.txt`

8 of 8 new or added lines in 2 files covered. (100.0%)

1293 existing lines in 43 files now uncovered.

50627 of 59548 relevant lines covered (85.02%)

35584918.67 hits per line

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

84.33
/src/reaction.cpp
1
#include "openmc/reaction.h"
2

3
#include <string>
4
#include <unordered_map>
5
#include <utility> // for move
6

7
#include <fmt/core.h>
8

9
#include "openmc/constants.h"
10
#include "openmc/endf.h"
11
#include "openmc/hdf5_interface.h"
12
#include "openmc/random_lcg.h"
13
#include "openmc/search.h"
14
#include "openmc/secondary_uncorrelated.h"
15

16
namespace openmc {
17

18
//==============================================================================
19
// Reaction implementation
20
//==============================================================================
21

22
Reaction::Reaction(hid_t group, const vector<int>& temperatures)
1,326,956✔
23
{
24
  read_attribute(group, "Q_value", q_value_);
1,326,956✔
25
  read_attribute(group, "mt", mt_);
1,326,956✔
26
  int tmp;
27
  read_attribute(group, "center_of_mass", tmp);
1,326,956✔
28
  scatter_in_cm_ = (tmp == 1);
1,326,956✔
29

30
  // Checks if redudant attribute exists before loading
31
  // (for compatibiltiy with legacy .h5 libraries)
32
  if (attribute_exists(group, "redundant")) {
1,326,956✔
33
    read_attribute(group, "redundant", tmp);
1,326,956✔
34
    redundant_ = (tmp == 1);
1,326,956✔
35
  } else {
36
    redundant_ = false;
×
37
  }
38

39
  // Read cross section and threshold_idx data
40
  for (auto t : temperatures) {
2,654,164✔
41
    // Get group corresponding to temperature
42
    hid_t temp_group = open_group(group, fmt::format("{}K", t).c_str());
1,327,208✔
43
    hid_t dset = open_dataset(temp_group, "xs");
1,327,208✔
44

45
    // Get threshold index
46
    TemperatureXS xs;
1,327,208✔
47
    read_attribute(dset, "threshold_idx", xs.threshold);
1,327,208✔
48

49
    // Read cross section values
50
    read_dataset(dset, xs.value);
1,327,208✔
51
    close_dataset(dset);
1,327,208✔
52
    close_group(temp_group);
1,327,208✔
53

54
    // create new entry in xs vector
55
    xs_.push_back(std::move(xs));
1,327,208✔
56
  }
1,327,208✔
57

58
  // Read products
59
  for (const auto& name : group_names(group)) {
6,113,826✔
60
    if (name.rfind("product_", 0) == 0) {
4,786,870✔
61
      hid_t pgroup = open_group(group, name.c_str());
3,453,402✔
62
      products_.emplace_back(pgroup);
3,453,402✔
63
      close_group(pgroup);
3,453,402✔
64
    }
65
  }
1,326,956✔
66
}
1,326,956✔
67

68
double Reaction::xs(int64_t i_temp, int64_t i_grid, double interp_factor) const
815,690,508✔
69
{
70
  // If energy is below threshold, return 0. Otherwise interpolate between
71
  // nearest grid points
72
  const auto& x = xs_[i_temp];
815,690,508✔
73
  return (i_grid < x.threshold)
815,690,508✔
74
           ? 0.0
815,690,508✔
75
           : (1.0 - interp_factor) * x.value[i_grid - x.threshold] +
477,043,995✔
76
               interp_factor * x.value[i_grid - x.threshold + 1];
815,690,508✔
77
}
78

79
double Reaction::xs(const NuclideMicroXS& micro) const
815,690,508✔
80
{
81
  return this->xs(micro.index_temp, micro.index_grid, micro.interp_factor);
815,690,508✔
82
}
83

84
double Reaction::collapse_rate(int64_t i_temp, span<const double> energy,
816✔
85
  span<const double> flux, const vector<double>& grid) const
86
{
87
  // Find index corresponding to first energy
88
  const auto& xs = xs_[i_temp].value;
816✔
89
  int i_low = lower_bound_index(grid.cbegin(), grid.cend(), energy.front());
816✔
90

91
  // Check for threshold and adjust starting point if necessary
92
  int j_start = 0;
816✔
93
  int i_threshold = xs_[i_temp].threshold;
816✔
94
  if (i_low < i_threshold) {
816✔
95
    i_low = i_threshold;
792✔
96
    while (energy[j_start + 1] < grid[i_low]) {
1,116✔
97
      ++j_start;
324✔
98
      if (j_start + 1 == energy.size())
324✔
UNCOV
99
        return 0.0;
×
100
    }
101
  }
102

103
  double xs_flux_sum = 0.0;
816✔
104

105
  for (int j = j_start; j < flux.size(); ++j) {
5,604✔
106
    double E_group_low = energy[j];
4,788✔
107
    double E_group_high = energy[j + 1];
4,788✔
108
    double flux_per_eV = flux[j] / (E_group_high - E_group_low);
4,788✔
109

110
    // Determine energy grid index corresponding to group high
111
    int i_high = i_low;
4,788✔
112
    while (grid[i_high + 1] < E_group_high && i_high + 1 < grid.size() - 1)
30,743,544✔
113
      ++i_high;
30,738,756✔
114

115
    // Loop over energy grid points within [E_group_low, E_group_high]
116
    for (; i_low <= i_high; ++i_low) {
30,748,332✔
117
      // Determine bounding grid energies and cross sections
118
      double E_l = grid[i_low];
30,743,544✔
119
      double E_r = grid[i_low + 1];
30,743,544✔
120
      if (E_l == E_r)
30,743,544✔
UNCOV
121
        continue;
×
122

123
      double xs_l = xs[i_low - i_threshold];
30,743,544✔
124
      double xs_r = xs[i_low + 1 - i_threshold];
30,743,544✔
125

126
      // Determine actual energies
127
      double E_low = std::max(E_group_low, E_l);
30,743,544✔
128
      double E_high = std::min(E_group_high, E_r);
30,743,544✔
129

130
      // Determine average cross section across segment
131
      double m = (xs_r - xs_l) / (E_r - E_l);
30,743,544✔
132
      double xs_low = xs_l + m * (E_low - E_l);
30,743,544✔
133
      double xs_high = xs_l + m * (E_high - E_l);
30,743,544✔
134
      double xs_avg = 0.5 * (xs_low + xs_high);
30,743,544✔
135

136
      // Add contribution from segment
137
      double dE = (E_high - E_low);
30,743,544✔
138
      xs_flux_sum += flux_per_eV * xs_avg * dE;
30,743,544✔
139
    }
140

141
    i_low = i_high;
4,788✔
142

143
    // Check for end of energy grid
144
    if (i_low + 1 == grid.size())
4,788✔
UNCOV
145
      break;
×
146
  }
147

148
  return xs_flux_sum;
816✔
149
}
150

151
//==============================================================================
152
// Non-member functions
153
//==============================================================================
154

155
std::unordered_map<int, std::string> REACTION_NAME_MAP {
156
  {SCORE_FLUX, "flux"},
157
  {SCORE_TOTAL, "total"},
158
  {SCORE_SCATTER, "scatter"},
159
  {SCORE_NU_SCATTER, "nu-scatter"},
160
  {SCORE_ABSORPTION, "absorption"},
161
  {SCORE_FISSION, "fission"},
162
  {SCORE_NU_FISSION, "nu-fission"},
163
  {SCORE_DECAY_RATE, "decay-rate"},
164
  {SCORE_DELAYED_NU_FISSION, "delayed-nu-fission"},
165
  {SCORE_PROMPT_NU_FISSION, "prompt-nu-fission"},
166
  {SCORE_KAPPA_FISSION, "kappa-fission"},
167
  {SCORE_CURRENT, "current"},
168
  {SCORE_EVENTS, "events"},
169
  {SCORE_INVERSE_VELOCITY, "inverse-velocity"},
170
  {SCORE_FISS_Q_PROMPT, "fission-q-prompt"},
171
  {SCORE_FISS_Q_RECOV, "fission-q-recoverable"},
172
  {SCORE_PULSE_HEIGHT, "pulse-height"},
173
  // Normal ENDF-based reactions
174
  {TOTAL_XS, "(n,total)"},
175
  {ELASTIC, "(n,elastic)"},
176
  {N_NONELASTIC, "(n,nonelastic)"},
177
  {N_LEVEL, "(n,level)"},
178
  {N_2ND, "(n,2nd)"},
179
  {N_2N, "(n,2n)"},
180
  {N_3N, "(n,3n)"},
181
  {N_FISSION, "(n,fission)"},
182
  {N_F, "(n,f)"},
183
  {N_NF, "(n,nf)"},
184
  {N_2NF, "(n,2nf)"},
185
  {N_NA, "(n,na)"},
186
  {N_N3A, "(n,n3a)"},
187
  {N_2NA, "(n,2na)"},
188
  {N_3NA, "(n,3na)"},
189
  {N_NP, "(n,np)"},
190
  {N_N2A, "(n,n2a)"},
191
  {N_2N2A, "(n,2n2a)"},
192
  {N_ND, "(n,nd)"},
193
  {N_NT, "(n,nt)"},
194
  {N_N3HE, "(n,n3He)"},
195
  {N_ND2A, "(n,nd2a)"},
196
  {N_NT2A, "(n,nt2a)"},
197
  {N_4N, "(n,4n)"},
198
  {N_3NF, "(n,3nf)"},
199
  {N_2NP, "(n,2np)"},
200
  {N_3NP, "(n,3np)"},
201
  {N_N2P, "(n,n2p)"},
202
  {N_NPA, "(n,npa)"},
203
  {N_NC, "(n,nc)"},
204
  {N_DISAPPEAR, "(n,disappear)"},
205
  {N_GAMMA, "(n,gamma)"},
206
  {N_P, "(n,p)"},
207
  {N_D, "(n,d)"},
208
  {N_T, "(n,t)"},
209
  {N_3HE, "(n,3He)"},
210
  {N_A, "(n,a)"},
211
  {N_2A, "(n,2a)"},
212
  {N_3A, "(n,3a)"},
213
  {N_2P, "(n,2p)"},
214
  {N_PA, "(n,pa)"},
215
  {N_T2A, "(n,t2a)"},
216
  {N_D2A, "(n,d2a)"},
217
  {N_PD, "(n,pd)"},
218
  {N_PT, "(n,pt)"},
219
  {N_DA, "(n,da)"},
220
  {N_5N, "(n,5n)"},
221
  {N_6N, "(n,6n)"},
222
  {N_2NT, "(n,2nt)"},
223
  {N_TA, "(n,ta)"},
224
  {N_4NP, "(n,4np)"},
225
  {N_3ND, "(n,3nd)"},
226
  {N_NDA, "(n,nda)"},
227
  {N_2NPA, "(n,2npa)"},
228
  {N_7N, "(n,7n)"},
229
  {N_8N, "(n,8n)"},
230
  {N_5NP, "(n,5np)"},
231
  {N_6NP, "(n,6np)"},
232
  {N_7NP, "(n,7np)"},
233
  {N_4NA, "(n,4na)"},
234
  {N_5NA, "(n,5na)"},
235
  {N_6NA, "(n,6na)"},
236
  {N_7NA, "(n,7na)"},
237
  {N_4ND, "(n,4nd)"},
238
  {N_5ND, "(n,5nd)"},
239
  {N_6ND, "(n,6nd)"},
240
  {N_3NT, "(n,3nt)"},
241
  {N_4NT, "(n,4nt)"},
242
  {N_5NT, "(n,5nt)"},
243
  {N_6NT, "(n,6nt)"},
244
  {N_2N3HE, "(n,2n3He)"},
245
  {N_3N3HE, "(n,3n3He)"},
246
  {N_4N3HE, "(n,4n3He)"},
247
  {N_3N2P, "(n,3n2p)"},
248
  {N_3N2A, "(n,3n2a)"},
249
  {N_3NPA, "(n,3npa)"},
250
  {N_DT, "(n,dt)"},
251
  {N_NPD, "(n,npd)"},
252
  {N_NPT, "(n,npt)"},
253
  {N_NDT, "(n,ndt)"},
254
  {N_NP3HE, "(n,np3He)"},
255
  {N_ND3HE, "(n,nd3He)"},
256
  {N_NT3HE, "(n,nt3He)"},
257
  {N_NTA, "(n,nta)"},
258
  {N_2N2P, "(n,2n2p)"},
259
  {N_P3HE, "(n,p3He)"},
260
  {N_D3HE, "(n,d3He)"},
261
  {N_3HEA, "(n,3Hea)"},
262
  {N_4N2P, "(n,4n2p)"},
263
  {N_4N2A, "(n,4n2a)"},
264
  {N_4NPA, "(n,4npa)"},
265
  {N_3P, "(n,3p)"},
266
  {N_N3P, "(n,n3p)"},
267
  {N_3N2PA, "(n,3n2pa)"},
268
  {N_5N2P, "(n,5n2p)"},
269
  {201, "(n,Xn)"},
270
  {202, "(n,Xgamma)"},
271
  {N_XP, "(n,Xp)"},
272
  {N_XD, "(n,Xd)"},
273
  {N_XT, "(n,Xt)"},
274
  {N_X3HE, "(n,X3He)"},
275
  {N_XA, "(n,Xa)"},
276
  {HEATING, "heating"},
277
  {DAMAGE_ENERGY, "damage-energy"},
278
  {COHERENT, "coherent-scatter"},
279
  {INCOHERENT, "incoherent-scatter"},
280
  {PAIR_PROD_ELEC, "pair-production-electron"},
281
  {PAIR_PROD, "pair-production"},
282
  {PAIR_PROD_NUC, "pair-production-nuclear"},
283
  {PHOTOELECTRIC, "photoelectric"},
284
  {N_PC, "(n,pc)"},
285
  {N_DC, "(n,dc)"},
286
  {N_TC, "(n,tc)"},
287
  {N_3HEC, "(n,3Hec)"},
288
  {N_AC, "(n,ac)"},
289
  {N_2NC, "(n,2nc)"},
290
  {HEATING_LOCAL, "heating-local"},
291
};
292

293
std::unordered_map<std::string, int> REACTION_TYPE_MAP;
294

295
void initialize_maps()
3,110✔
296
{
297
  // Add level reactions to name map
298
  for (int level = 0; level <= 48; ++level) {
155,500✔
299
    if (level >= 1 && level <= 40) {
152,390✔
300
      REACTION_NAME_MAP[50 + level] = fmt::format("(n,n{})", level);
124,400✔
301
    }
302
    REACTION_NAME_MAP[600 + level] = fmt::format("(n,p{})", level);
152,390✔
303
    REACTION_NAME_MAP[650 + level] = fmt::format("(n,d{})", level);
152,390✔
304
    REACTION_NAME_MAP[700 + level] = fmt::format("(n,t{})", level);
152,390✔
305
    REACTION_NAME_MAP[750 + level] = fmt::format("(n,3He{})", level);
152,390✔
306
    REACTION_NAME_MAP[800 + level] = fmt::format("(n,a{})", level);
152,390✔
307
    if (level <= 15) {
152,390✔
308
      REACTION_NAME_MAP[875 + level] = fmt::format("(n,2n{})", level);
49,760✔
309
    }
310
  }
311

312
  // Create photoelectric subshells
313
  for (int mt = 534; mt <= 572; ++mt) {
124,400✔
314
    REACTION_NAME_MAP[mt] =
142,389✔
315
      fmt::format("photoelectric, {} subshell", SUBSHELLS[mt - 534]);
342,771✔
316
  }
317

318
  // Invert name map to create type map
319
  for (const auto& kv : REACTION_NAME_MAP) {
1,477,250✔
320
    REACTION_TYPE_MAP[kv.second] = kv.first;
1,474,140✔
321
  }
322
}
3,110✔
323

324
std::string reaction_name(int mt)
60,618✔
325
{
326
  // Initialize remainder of name map and all of type map
327
  if (REACTION_TYPE_MAP.empty())
60,618✔
328
    initialize_maps();
36✔
329

330
  // Get reaction name from map
331
  auto it = REACTION_NAME_MAP.find(mt);
60,618✔
332
  if (it != REACTION_NAME_MAP.end()) {
60,618✔
333
    return it->second;
60,618✔
334
  } else {
UNCOV
335
    return fmt::format("MT={}", mt);
×
336
  }
337
}
338

339
int reaction_type(std::string name)
37,736✔
340
{
341
  // Initialize remainder of name map and all of type map
342
  if (REACTION_TYPE_MAP.empty())
37,736✔
343
    initialize_maps();
3,074✔
344

345
  // (n,total) exists in REACTION_TYPE_MAP for MT=1, but we need this to return
346
  // the special SCORE_TOTAL score
347
  if (name == "(n,total)")
37,736✔
UNCOV
348
    return SCORE_TOTAL;
×
349

350
  // Check if type map has an entry for this reaction name
351
  auto it = REACTION_TYPE_MAP.find(name);
37,736✔
352
  if (it != REACTION_TYPE_MAP.end()) {
37,736✔
353
    return it->second;
37,448✔
354
  }
355

356
  // Alternate names for several reactions
357
  if (name == "elastic") {
288✔
358
    return ELASTIC;
179✔
359
  } else if (name == "n2n") {
109✔
UNCOV
360
    return N_2N;
×
361
  } else if (name == "n3n") {
109✔
362
    return N_3N;
×
363
  } else if (name == "n4n") {
109✔
364
    return N_4N;
×
365
  } else if (name == "H1-production") {
109✔
366
    return N_XP;
17✔
367
  } else if (name == "H2-production") {
92✔
368
    return N_XD;
17✔
369
  } else if (name == "H3-production") {
75✔
370
    return N_XT;
17✔
371
  } else if (name == "He3-production") {
58✔
372
    return N_X3HE;
41✔
373
  } else if (name == "He4-production") {
17✔
374
    return N_XA;
17✔
375
  }
376

377
  // Assume the given string is a reaction MT number.  Make sure it's a natural
378
  // number then return.
UNCOV
379
  int MT = 0;
×
380
  try {
381
    MT = std::stoi(name);
×
UNCOV
382
  } catch (const std::invalid_argument& ex) {
×
383
    throw std::invalid_argument(
×
384
      "Invalid tally score \"" + name +
×
385
      "\". See the docs "
386
      "for details: "
UNCOV
387
      "https://docs.openmc.org/en/stable/usersguide/tallies.html#scores");
×
UNCOV
388
  }
×
389
  if (MT < 1)
×
390
    throw std::invalid_argument(
×
391
      "Invalid tally score \"" + name +
×
392
      "\". See the docs "
393
      "for details: "
UNCOV
394
      "https://docs.openmc.org/en/stable/usersguide/tallies.html#scores");
×
UNCOV
395
  return MT;
×
396
}
397

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

© 2025 Coveralls, Inc