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

openmc-dev / openmc / 23020451435

12 Mar 2026 07:39PM UTC coverage: 81.638% (+0.07%) from 81.566%
23020451435

Pull #3863

github

web-flow
Merge 5541c0bf0 into 27522fe85
Pull Request #3863: Shared Secondary Particle Bank

17761 of 25507 branches covered (69.63%)

Branch coverage included in aggregate %.

418 of 435 new or added lines in 17 files covered. (96.09%)

10 existing lines in 4 files now uncovered.

58266 of 67620 relevant lines covered (86.17%)

46325054.24 hits per line

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

84.13
/src/physics.cpp
1
#include "openmc/physics.h"
2

3
#include "openmc/bank.h"
4
#include "openmc/bremsstrahlung.h"
5
#include "openmc/chain.h"
6
#include "openmc/constants.h"
7
#include "openmc/distribution_multi.h"
8
#include "openmc/eigenvalue.h"
9
#include "openmc/endf.h"
10
#include "openmc/error.h"
11
#include "openmc/ifp.h"
12
#include "openmc/material.h"
13
#include "openmc/math_functions.h"
14
#include "openmc/message_passing.h"
15
#include "openmc/ncrystal_interface.h"
16
#include "openmc/nuclide.h"
17
#include "openmc/photon.h"
18
#include "openmc/physics_common.h"
19
#include "openmc/random_dist.h"
20
#include "openmc/random_lcg.h"
21
#include "openmc/reaction.h"
22
#include "openmc/search.h"
23
#include "openmc/secondary_uncorrelated.h"
24
#include "openmc/settings.h"
25
#include "openmc/simulation.h"
26
#include "openmc/string_utils.h"
27
#include "openmc/tallies/tally.h"
28
#include "openmc/thermal.h"
29
#include "openmc/weight_windows.h"
30

31
#include <fmt/core.h>
32

33
#include "openmc/tensor.h"
34
#include <algorithm> // for max, min, max_element
35
#include <cmath>     // for sqrt, exp, log, abs, copysign
36

37
namespace openmc {
38

39
//==============================================================================
40
// Non-member functions
41
//==============================================================================
42

43
void collision(Particle& p)
1,261,107,495✔
44
{
45
  // Add to collision counter for particle
46
  ++(p.n_collision());
1,261,107,495✔
47
  p.secondary_bank_index() = p.local_secondary_bank().size();
1,261,107,495!
48

49
  // Sample reaction for the material the particle is in
50
  switch (p.type().pdg_number()) {
1,261,107,495!
51
  case PDG_NEUTRON:
1,139,872,269✔
52
    sample_neutron_reaction(p);
1,139,872,269✔
53
    break;
1,139,872,269✔
54
  case PDG_PHOTON:
35,048,652✔
55
    sample_photon_reaction(p);
35,048,652✔
56
    break;
35,048,652✔
57
  case PDG_ELECTRON:
86,006,682✔
58
    sample_electron_reaction(p);
86,006,682✔
59
    break;
86,006,682✔
60
  case PDG_POSITRON:
179,892✔
61
    sample_positron_reaction(p);
179,892✔
62
    break;
179,892✔
63
  default:
×
64
    fatal_error("Unsupported particle PDG for collision sampling.");
×
65
  }
66

67
  if (settings::weight_windows_on) {
1,261,107,495✔
68
    auto [ww_found, ww] = search_weight_window(p);
269,941,573✔
69
    if (!ww_found && p.type() == ParticleType::neutron()) {
269,941,573✔
70
      // if the weight window is not valid, apply russian roulette for neutrons
71
      // (regardless of weight window collision checkpoint setting)
72
      apply_russian_roulette(p);
66,697✔
73
    } else if (settings::weight_window_checkpoint_collision) {
269,874,876!
74
      // if collision checkpointing is on, apply weight window
75
      apply_weight_window(p, ww);
269,874,876✔
76
    }
77
  }
78

79
  // Kill particle if energy falls below cutoff
80
  int type = p.type().transport_index();
1,261,107,495!
81
  if (type != C_NONE && p.E() < settings::energy_cutoff[type]) {
1,261,107,495!
82
    p.wgt() = 0.0;
7,849,553✔
83
  }
84

85
  // Display information about collision
86
  if (settings::verbosity >= 10 || p.trace()) {
1,261,107,495!
87
    std::string msg;
121!
88
    if (p.event() == TallyEvent::KILL) {
121!
89
      msg = fmt::format("    Killed. Energy = {} eV.", p.E());
×
90
    } else if (p.type().is_neutron()) {
121!
91
      msg = fmt::format("    {} with {}. Energy = {} eV.",
242✔
92
        reaction_name(p.event_mt()), data::nuclides[p.event_nuclide()]->name_,
242✔
93
        p.E());
121✔
94
    } else if (p.type().is_photon()) {
×
95
      msg = fmt::format("    {} with {}. Energy = {} eV.",
×
96
        reaction_name(p.event_mt()),
×
97
        to_element(data::nuclides[p.event_nuclide()]->name_), p.E());
×
98
    } else {
99
      msg = fmt::format("    Disappeared. Energy = {} eV.", p.E());
×
100
    }
101
    write_message(msg, 1);
121✔
102
  }
121✔
103
}
1,261,107,495✔
104

105
void sample_neutron_reaction(Particle& p)
1,139,872,269✔
106
{
107
  // Sample a nuclide within the material
108
  int i_nuclide = sample_nuclide(p);
1,139,872,269✔
109

110
  // Save which nuclide particle had collision with
111
  p.event_nuclide() = i_nuclide;
1,139,872,269✔
112

113
  // Create fission bank sites. Note that while a fission reaction is sampled,
114
  // it never actually "happens", i.e. the weight of the particle does not
115
  // change when sampling fission sites. The following block handles all
116
  // absorption (including fission)
117

118
  const auto& nuc {data::nuclides[i_nuclide]};
1,139,872,269✔
119

120
  if (nuc->fissionable_ && p.neutron_xs(i_nuclide).fission > 0.0) {
1,139,872,269✔
121
    auto& rx = sample_fission(i_nuclide, p);
125,773,799✔
122
    if (settings::run_mode == RunMode::EIGENVALUE) {
125,773,799✔
123
      create_fission_sites(p, i_nuclide, rx);
100,745,246✔
124
    } else if (settings::run_mode == RunMode::FIXED_SOURCE &&
25,028,553✔
125
               settings::create_fission_neutrons) {
126
      create_fission_sites(p, i_nuclide, rx);
582,505✔
127

128
      // Make sure particle population doesn't grow out of control for
129
      // subcritical multiplication problems.
130
      if (p.local_secondary_bank().size() >= settings::max_secondaries &&
582,505!
NEW
131
          !settings::use_shared_secondary_bank) {
×
UNCOV
132
        fatal_error(
×
133
          "The secondary particle bank appears to be growing without "
134
          "bound. You are likely running a subcritical multiplication problem "
135
          "with k-effective close to or greater than one.");
136
      }
137
    }
138
    p.event_mt() = rx.mt_;
125,773,799✔
139
  }
140

141
  // Create secondary photons
142
  if (settings::photon_transport) {
1,139,872,269✔
143
    sample_secondary_photons(p, i_nuclide);
66,629,112✔
144
  }
145

146
  // If survival biasing is being used, the following subroutine adjusts the
147
  // weight of the particle. Otherwise, it checks to see if absorption occurs
148

149
  if (p.neutron_xs(i_nuclide).absorption > 0.0) {
1,139,872,269✔
150
    absorption(p, i_nuclide);
1,139,868,749✔
151
  }
152
  if (!p.alive())
1,139,872,269✔
153
    return;
154

155
  // Sample a scattering reaction and determine the secondary energy of the
156
  // exiting neutron
157
  const auto& ncrystal_mat = model::materials[p.material()]->ncrystal_mat();
1,114,262,908✔
158
  if (ncrystal_mat && p.E() < NCRYSTAL_MAX_ENERGY) {
1,114,262,908!
159
    ncrystal_mat.scatter(p);
158,829✔
160
  } else {
161
    scatter(p, i_nuclide);
1,114,104,079✔
162
  }
163

164
  // Advance URR seed stream 'N' times after energy changes
165
  if (p.E() != p.E_last()) {
1,114,262,908✔
166
    advance_prn_seed(data::nuclides.size(), &p.seeds(STREAM_URR_PTABLE));
1,113,943,677✔
167
  }
168

169
  // Play russian roulette if there are no weight windows
170
  if (!settings::weight_windows_on)
1,114,262,908✔
171
    apply_russian_roulette(p);
914,907,189✔
172
}
173

174
void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx)
101,327,751✔
175
{
176
  // If uniform fission source weighting is turned on, we increase or decrease
177
  // the expected number of fission sites produced
178
  double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0;
101,327,751✔
179

180
  // Determine the expected number of neutrons produced
181
  double nu_t = p.wgt() / simulation::keff * weight *
101,327,751✔
182
                p.neutron_xs(i_nuclide).nu_fission /
101,327,751✔
183
                p.neutron_xs(i_nuclide).total;
101,327,751✔
184

185
  // Sample the number of neutrons produced
186
  int nu = static_cast<int>(nu_t);
101,327,751✔
187
  if (prn(p.current_seed()) <= (nu_t - nu))
101,327,751✔
188
    ++nu;
20,801,279✔
189

190
  // If no neutrons were produced then don't continue
191
  if (nu == 0)
101,327,751✔
192
    return;
76,510,364✔
193

194
  // Initialize the counter of delayed neutrons encountered for each delayed
195
  // group.
196
  double nu_d[MAX_DELAYED_GROUPS] = {0.};
24,817,387✔
197

198
  // Clear out particle's nu fission bank
199
  p.nu_bank().clear();
24,817,387✔
200

201
  p.fission() = true;
24,817,387✔
202

203
  // Determine whether to place fission sites into the shared fission bank
204
  // or the secondary particle bank.
205
  bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE);
24,817,387✔
206

207
  // Counter for the number of fission sites successfully stored to the shared
208
  // fission bank or the secondary particle bank
209
  int n_sites_stored;
24,817,387✔
210

211
  for (n_sites_stored = 0; n_sites_stored < nu; n_sites_stored++) {
55,451,397✔
212
    // Initialize fission site object with particle data
213
    SourceSite site;
30,634,010✔
214
    site.r = p.r();
30,634,010✔
215
    site.particle = ParticleType::neutron();
30,634,010✔
216
    site.time = p.time();
30,634,010✔
217
    site.wgt = 1. / weight;
30,634,010✔
218
    site.surf_id = 0;
30,634,010✔
219

220
    // Sample delayed group and angle/energy for fission reaction
221
    sample_fission_neutron(i_nuclide, rx, &site, p);
30,634,010✔
222

223
    // Reject site if it exceeds time cutoff
224
    if (site.delayed_group > 0) {
30,634,010✔
225
      double t_cutoff = settings::time_cutoff[site.particle.transport_index()];
192,175!
226
      if (site.time > t_cutoff) {
192,175!
227
        continue;
×
228
      }
229
    }
230

231
    // Set parent and progeny IDs
232
    site.parent_id = p.current_work();
30,634,010✔
233
    site.progeny_id = p.n_progeny()++;
30,634,010✔
234

235
    // Store fission site in bank
236
    if (use_fission_bank) {
30,634,010✔
237
      int64_t idx = simulation::fission_bank.thread_safe_append(site);
30,419,777✔
238
      if (idx == -1) {
30,419,777!
239
        warning(
×
240
          "The shared fission bank is full. Additional fission sites created "
241
          "in this generation will not be banked. Results may be "
242
          "non-deterministic.");
243

244
        // Decrement number of particle progeny as storage was unsuccessful.
245
        // This step is needed so that the sum of all progeny is equal to the
246
        // size of the shared fission bank.
247
        p.n_progeny()--;
×
248

249
        // Break out of loop as no more sites can be added to fission bank
250
        break;
×
251
      }
252
      // Iterated Fission Probability (IFP) method
253
      if (settings::ifp_on) {
30,419,777✔
254
        ifp(p, idx);
1,352,626✔
255
      }
256
    } else {
257
      site.wgt_born = p.wgt_born();
214,233✔
258
      site.wgt_ww_born = p.wgt_ww_born();
214,233✔
259
      site.n_split = p.n_split();
214,233✔
260
      p.local_secondary_bank().push_back(site);
214,233✔
261
      p.n_secondaries()++;
214,233✔
262
    }
263

264
    // Increment the number of neutrons born delayed
265
    if (site.delayed_group > 0) {
30,634,010✔
266
      nu_d[site.delayed_group - 1]++;
192,175✔
267
    }
268

269
    // Write fission particles to nuBank
270
    NuBank& nu_bank_entry = p.nu_bank().emplace_back();
30,634,010✔
271
    nu_bank_entry.wgt = site.wgt;
30,634,010✔
272
    nu_bank_entry.E = site.E;
30,634,010✔
273
    nu_bank_entry.delayed_group = site.delayed_group;
30,634,010✔
274
  }
275

276
  // If shared fission bank was full, and no fissions could be added,
277
  // set the particle fission flag to false.
278
  if (n_sites_stored == 0) {
24,817,387!
279
    p.fission() = false;
×
280
    return;
×
281
  }
282

283
  // Set nu to the number of fission sites successfully stored. If the fission
284
  // bank was not found to be full then these values are already equivalent.
285
  nu = n_sites_stored;
24,817,387✔
286

287
  // Store the total weight banked for analog fission tallies
288
  p.n_bank() = nu;
24,817,387✔
289
  p.wgt_bank() = nu / weight;
24,817,387✔
290
  for (size_t d = 0; d < MAX_DELAYED_GROUPS; d++) {
223,356,483✔
291
    p.n_delayed_bank(d) = nu_d[d];
198,539,096✔
292
  }
293
}
294

295
void sample_photon_reaction(Particle& p)
35,048,652✔
296
{
297
  // Kill photon if below energy cutoff -- an extra check is made here because
298
  // photons with energy below the cutoff may have been produced by neutrons
299
  // reactions or atomic relaxation
300
  int photon = ParticleType::photon().transport_index();
35,048,652✔
301
  if (p.E() < settings::energy_cutoff[photon]) {
35,048,652!
302
    p.E() = 0.0;
×
303
    p.wgt() = 0.0;
×
304
    return;
×
305
  }
306

307
  // Sample element within material
308
  int i_element = sample_element(p);
35,048,652✔
309
  const auto& micro {p.photon_xs(i_element)};
35,048,652✔
310
  const auto& element {*data::elements[i_element]};
35,048,652✔
311

312
  // Calculate photon energy over electron rest mass equivalent
313
  double alpha = p.E() / MASS_ELECTRON_EV;
35,048,652✔
314

315
  // For tallying purposes, this routine might be called directly. In that
316
  // case, we need to sample a reaction via the cutoff variable
317
  double prob = 0.0;
35,048,652✔
318
  double cutoff = prn(p.current_seed()) * micro.total;
35,048,652✔
319

320
  // Coherent (Rayleigh) scattering
321
  prob += micro.coherent;
35,048,652✔
322
  if (prob > cutoff) {
35,048,652✔
323
    p.mu() = element.rayleigh_scatter(alpha, p.current_seed());
1,738,644✔
324
    p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed());
1,738,644✔
325
    p.event() = TallyEvent::SCATTER;
1,738,644✔
326
    p.event_mt() = COHERENT;
1,738,644✔
327
    return;
1,738,644✔
328
  }
329

330
  // Incoherent (Compton) scattering
331
  prob += micro.incoherent;
33,310,008✔
332
  if (prob > cutoff) {
33,310,008✔
333
    double alpha_out;
25,471,378✔
334
    int i_shell;
25,471,378✔
335
    element.compton_scatter(
25,471,378✔
336
      alpha, true, &alpha_out, &p.mu(), &i_shell, p.current_seed());
25,471,378✔
337

338
    // Determine binding energy of shell. The binding energy is 0.0 if
339
    // doppler broadening is not used.
340
    double e_b;
25,471,378✔
341
    if (i_shell == -1) {
25,471,378!
342
      e_b = 0.0;
343
    } else {
344
      e_b = element.binding_energy_[i_shell];
25,471,378✔
345
    }
346

347
    // Create Compton electron
348
    double phi = uniform_distribution(0., 2.0 * PI, p.current_seed());
25,471,378✔
349
    double E_electron = (alpha - alpha_out) * MASS_ELECTRON_EV - e_b;
25,471,378✔
350
    int electron = ParticleType::electron().transport_index();
25,471,378✔
351
    if (E_electron >= settings::energy_cutoff[electron]) {
25,471,378✔
352
      double mu_electron = (alpha - alpha_out * p.mu()) /
25,274,368✔
353
                           std::sqrt(alpha * alpha + alpha_out * alpha_out -
25,274,368✔
354
                                     2.0 * alpha * alpha_out * p.mu());
25,274,368✔
355
      Direction u = rotate_angle(p.u(), mu_electron, &phi, p.current_seed());
25,274,368✔
356
      p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron());
25,274,368✔
357
    }
358

359
    // Allow electrons to fill orbital and produce Auger electrons and
360
    // fluorescent photons. Since Compton subshell data does not match atomic
361
    // relaxation data, use the mapping between the data to find the subshell
362
    if (settings::atomic_relaxation && i_shell >= 0 &&
25,471,378!
363
        element.subshell_map_[i_shell] >= 0) {
25,339,675!
364
      element.atomic_relaxation(element.subshell_map_[i_shell], p);
25,339,675✔
365
    }
366

367
    phi += PI;
25,471,378✔
368
    p.E() = alpha_out * MASS_ELECTRON_EV;
25,471,378✔
369
    p.u() = rotate_angle(p.u(), p.mu(), &phi, p.current_seed());
25,471,378✔
370
    p.event() = TallyEvent::SCATTER;
25,471,378✔
371
    p.event_mt() = INCOHERENT;
25,471,378✔
372
    return;
25,471,378✔
373
  }
374

375
  // Photoelectric effect
376
  double prob_after = prob + micro.photoelectric;
7,838,630✔
377

378
  if (prob_after > cutoff) {
7,838,630✔
379
    // Get grid index, interpolation factor, and bounding subshell
380
    // cross sections
381
    int i_grid = micro.index_grid;
7,658,738✔
382
    double f = micro.interp_factor;
7,658,738✔
383
    tensor::View<const double> xs_lower = element.cross_sections_.slice(i_grid);
7,658,738✔
384
    tensor::View<const double> xs_upper =
7,658,738✔
385
      element.cross_sections_.slice(i_grid + 1);
7,658,738✔
386

387
    for (int i_shell = 0; i_shell < element.shells_.size(); ++i_shell) {
27,275,524!
388
      const auto& shell {element.shells_[i_shell]};
27,275,524✔
389

390
      // Check threshold of reaction
391
      if (xs_lower(i_shell) == 0)
27,275,524✔
392
        continue;
10,317,592✔
393

394
      //  Evaluation subshell photoionization cross section
395
      prob += std::exp(
16,957,932✔
396
        xs_lower(i_shell) + f * (xs_upper(i_shell) - xs_lower(i_shell)));
16,957,932✔
397

398
      if (prob > cutoff) {
16,957,932✔
399
        // Determine binding energy based on whether atomic relaxation data is
400
        // present (if not, use value from Compton profile data)
401
        double binding_energy = element.has_atomic_relaxation_
7,658,738✔
402
                                  ? shell.binding_energy
7,658,738!
403
                                  : element.binding_energy_[i_shell];
×
404

405
        // Determine energy of secondary electron
406
        double E_electron = p.E() - binding_energy;
7,658,738✔
407

408
        // Sample mu using non-relativistic Sauter distribution.
409
        // See Eqns 3.19 and 3.20 in "Implementing a photon physics
410
        // model in Serpent 2" by Toni Kaltiaisenaho
411
        double mu;
11,476,635✔
412
        while (true) {
11,476,635✔
413
          double r = prn(p.current_seed());
11,476,635✔
414
          if (4.0 * (1.0 - r) * r >= prn(p.current_seed())) {
11,476,635✔
415
            double rel_vel =
7,658,738✔
416
              std::sqrt(E_electron * (E_electron + 2.0 * MASS_ELECTRON_EV)) /
7,658,738✔
417
              (E_electron + MASS_ELECTRON_EV);
7,658,738✔
418
            mu =
7,658,738✔
419
              (2.0 * r + rel_vel - 1.0) / (2.0 * rel_vel * r - rel_vel + 1.0);
7,658,738✔
420
            break;
7,658,738✔
421
          }
422
        }
423

424
        double phi = uniform_distribution(0., 2.0 * PI, p.current_seed());
7,658,738✔
425
        Direction u;
7,658,738✔
426
        u.x = mu;
7,658,738✔
427
        u.y = std::sqrt(1.0 - mu * mu) * std::cos(phi);
7,658,738✔
428
        u.z = std::sqrt(1.0 - mu * mu) * std::sin(phi);
7,658,738✔
429

430
        // Create secondary electron
431
        p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron());
7,658,738✔
432

433
        // Allow electrons to fill orbital and produce auger electrons
434
        // and fluorescent photons
435
        if (settings::atomic_relaxation) {
7,658,738✔
436
          element.atomic_relaxation(i_shell, p);
7,548,738✔
437
        }
438
        p.event() = TallyEvent::ABSORB;
7,658,738✔
439
        p.event_mt() = 533 + shell.index_subshell;
7,658,738✔
440
        p.wgt() = 0.0;
7,658,738✔
441
        p.E() = 0.0;
7,658,738✔
442
        return;
7,658,738✔
443
      }
444
    }
445
  }
15,317,476✔
446
  prob = prob_after;
179,892✔
447

448
  // Pair production
449
  prob += micro.pair_production;
179,892✔
450
  if (prob > cutoff) {
179,892!
451
    double E_electron, E_positron;
179,892✔
452
    double mu_electron, mu_positron;
179,892✔
453
    element.pair_production(alpha, &E_electron, &E_positron, &mu_electron,
179,892✔
454
      &mu_positron, p.current_seed());
455

456
    // Create secondary electron
457
    Direction u = rotate_angle(p.u(), mu_electron, nullptr, p.current_seed());
179,892✔
458
    p.create_secondary(p.wgt(), u, E_electron, ParticleType::electron());
179,892✔
459

460
    // Create secondary positron
461
    u = rotate_angle(p.u(), mu_positron, nullptr, p.current_seed());
179,892✔
462
    p.create_secondary(p.wgt(), u, E_positron, ParticleType::positron());
179,892✔
463
    p.event() = TallyEvent::ABSORB;
179,892✔
464
    p.event_mt() = PAIR_PROD;
179,892✔
465
    p.wgt() = 0.0;
179,892✔
466
    p.E() = 0.0;
179,892✔
467
  }
468
}
469

470
void sample_electron_reaction(Particle& p)
86,006,682✔
471
{
472
  // TODO: create reaction types
473

474
  if (settings::electron_treatment == ElectronTreatment::TTB) {
86,006,682✔
475
    double E_lost;
85,462,292✔
476
    thick_target_bremsstrahlung(p, &E_lost);
85,462,292✔
477
  }
478

479
  p.E() = 0.0;
86,006,682✔
480
  p.wgt() = 0.0;
86,006,682✔
481
  p.event() = TallyEvent::ABSORB;
86,006,682✔
482
}
86,006,682✔
483

484
void sample_positron_reaction(Particle& p)
179,892✔
485
{
486
  // TODO: create reaction types
487

488
  if (settings::electron_treatment == ElectronTreatment::TTB) {
179,892✔
489
    double E_lost;
178,627✔
490
    thick_target_bremsstrahlung(p, &E_lost);
178,627✔
491
  }
492

493
  // Sample angle isotropically
494
  Direction u = isotropic_direction(p.current_seed());
179,892✔
495

496
  // Create annihilation photon pair traveling in opposite directions
497
  p.create_secondary(p.wgt(), u, MASS_ELECTRON_EV, ParticleType::photon());
179,892✔
498
  p.create_secondary(p.wgt(), -u, MASS_ELECTRON_EV, ParticleType::photon());
179,892✔
499

500
  p.E() = 0.0;
179,892✔
501
  p.wgt() = 0.0;
179,892✔
502
  p.event() = TallyEvent::ABSORB;
179,892✔
503
}
179,892✔
504

505
int sample_nuclide(Particle& p)
1,139,872,269✔
506
{
507
  // Sample cumulative distribution function
508
  double cutoff = prn(p.current_seed()) * p.macro_xs().total;
1,139,872,269✔
509

510
  // Get pointers to nuclide/density arrays
511
  const auto& mat {model::materials[p.material()]};
1,139,872,269✔
512
  int n = mat->nuclide_.size();
1,139,872,269✔
513

514
  double prob = 0.0;
1,139,872,269✔
515
  for (int i = 0; i < n; ++i) {
2,147,483,647!
516
    // Get atom density
517
    int i_nuclide = mat->nuclide_[i];
2,147,483,647✔
518
    double atom_density = mat->atom_density(i, p.density_mult());
2,147,483,647✔
519

520
    // Increment probability to compare to cutoff
521
    prob += atom_density * p.neutron_xs(i_nuclide).total;
2,147,483,647✔
522
    if (prob >= cutoff)
2,147,483,647✔
523
      return i_nuclide;
1,139,872,269✔
524
  }
525

526
  // If we reach here, no nuclide was sampled
527
  p.write_restart();
×
528
  throw std::runtime_error {"Did not sample any nuclide during collision."};
×
529
}
530

531
int sample_element(Particle& p)
35,048,652✔
532
{
533
  // Sample cumulative distribution function
534
  double cutoff = prn(p.current_seed()) * p.macro_xs().total;
35,048,652✔
535

536
  // Get pointers to elements, densities
537
  const auto& mat {model::materials[p.material()]};
35,048,652✔
538

539
  double prob = 0.0;
35,048,652✔
540
  for (int i = 0; i < mat->element_.size(); ++i) {
133,185,082!
541
    // Find atom density
542
    int i_element = mat->element_[i];
133,185,082✔
543
    double atom_density = mat->atom_density(i, p.density_mult());
133,185,082✔
544

545
    // Determine microscopic cross section
546
    double sigma = atom_density * p.photon_xs(i_element).total;
133,185,082✔
547

548
    // Increment probability to compare to cutoff
549
    prob += sigma;
133,185,082✔
550
    if (prob > cutoff) {
133,185,082✔
551
      // Save which nuclide particle had collision with for tally purpose
552
      p.event_nuclide() = mat->nuclide_[i];
35,048,652✔
553

554
      return i_element;
35,048,652✔
555
    }
556
  }
557

558
  // If we made it here, no element was sampled
559
  p.write_restart();
×
560
  fatal_error("Did not sample any element during collision.");
×
561
}
562

563
Reaction& sample_fission(int i_nuclide, Particle& p)
125,773,799✔
564
{
565
  // Get pointer to nuclide
566
  const auto& nuc {data::nuclides[i_nuclide]};
125,773,799✔
567

568
  // If we're in the URR, by default use the first fission reaction. We also
569
  // default to the first reaction if we know that there are no partial fission
570
  // reactions
571
  if (p.neutron_xs(i_nuclide).use_ptable || !nuc->has_partial_fission_) {
125,773,799✔
572
    return *nuc->fission_rx_[0];
125,743,166✔
573
  }
574

575
  // Check to see if we are in a windowed multipole range.  WMP only supports
576
  // the first fission reaction.
577
  if (nuc->multipole_) {
30,633✔
578
    if (p.E() >= nuc->multipole_->E_min_ && p.E() <= nuc->multipole_->E_max_) {
2,849!
579
      return *nuc->fission_rx_[0];
1,991✔
580
    }
581
  }
582

583
  // Get grid index and interpolation factor and sample fission cdf
584
  const auto& micro = p.neutron_xs(i_nuclide);
28,642✔
585
  double cutoff = prn(p.current_seed()) * p.neutron_xs(i_nuclide).fission;
28,642✔
586
  double prob = 0.0;
28,642✔
587

588
  // Loop through each partial fission reaction type
589
  for (auto& rx : nuc->fission_rx_) {
28,702!
590
    // add to cumulative probability
591
    prob += rx->xs(micro);
28,702✔
592

593
    // Create fission bank sites if fission occurs
594
    if (prob > cutoff)
28,702✔
595
      return *rx;
28,642✔
596
  }
597

598
  // If we reached here, no reaction was sampled
599
  throw std::runtime_error {
×
600
    "No fission reaction was sampled for " + nuc->name_};
×
601
}
602

603
void sample_photon_product(
2,667,896✔
604
  int i_nuclide, Particle& p, int* i_rx, int* i_product)
605
{
606
  // Get grid index and interpolation factor and sample photon production cdf
607
  const auto& micro = p.neutron_xs(i_nuclide);
2,667,896✔
608
  double cutoff = prn(p.current_seed()) * micro.photon_prod;
2,667,896✔
609
  double prob = 0.0;
2,667,896✔
610

611
  // Loop through each reaction type
612
  const auto& nuc {data::nuclides[i_nuclide]};
2,667,896✔
613
  for (int i = 0; i < nuc->reactions_.size(); ++i) {
48,809,772!
614
    // Evaluate neutron cross section
615
    const auto& rx = nuc->reactions_[i];
48,809,772✔
616
    double xs = rx->xs(micro);
48,809,772✔
617

618
    // if cross section is zero for this reaction, skip it
619
    if (xs == 0.0)
48,809,772✔
620
      continue;
30,015,007✔
621

622
    for (int j = 0; j < rx->products_.size(); ++j) {
134,132,570✔
623
      if (rx->products_[j].particle_.is_photon()) {
118,005,701✔
624
        // For fission, artificially increase the photon yield to account
625
        // for delayed photons
626
        double f = 1.0;
103,484,601✔
627
        if (settings::delayed_photon_scaling) {
103,484,601!
628
          if (is_fission(rx->mt_)) {
103,484,601✔
629
            if (nuc->prompt_photons_ && nuc->delayed_photons_) {
540,221!
630
              double energy_prompt = (*nuc->prompt_photons_)(p.E());
540,221✔
631
              double energy_delayed = (*nuc->delayed_photons_)(p.E());
540,221✔
632
              f = (energy_prompt + energy_delayed) / (energy_prompt);
540,221✔
633
            }
634
          }
635
        }
636

637
        // add to cumulative probability
638
        prob += f * (*rx->products_[j].yield_)(p.E()) * xs;
103,484,601✔
639

640
        *i_rx = i;
103,484,601✔
641
        *i_product = j;
103,484,601✔
642
        if (prob > cutoff)
103,484,601✔
643
          return;
644
      }
645
    }
646
  }
647
}
648

649
void absorption(Particle& p, int i_nuclide)
1,139,868,749✔
650
{
651
  if (settings::survival_biasing) {
1,139,868,749✔
652
    // Determine weight absorbed in survival biasing
653
    const double wgt_absorb = p.wgt() * p.neutron_xs(i_nuclide).absorption /
9,679,604✔
654
                              p.neutron_xs(i_nuclide).total;
9,679,604✔
655

656
    // Adjust weight of particle by probability of absorption
657
    p.wgt() -= wgt_absorb;
9,679,604✔
658

659
    // Score implicit absorption estimate of keff
660
    if (settings::run_mode == RunMode::EIGENVALUE) {
9,679,604✔
661
      p.keff_tally_absorption() += wgt_absorb *
499,950✔
662
                                   p.neutron_xs(i_nuclide).nu_fission /
499,950✔
663
                                   p.neutron_xs(i_nuclide).absorption;
499,950✔
664
    }
665
  } else {
666
    // See if disappearance reaction happens
667
    if (p.neutron_xs(i_nuclide).absorption >
1,130,189,145✔
668
        prn(p.current_seed()) * p.neutron_xs(i_nuclide).total) {
1,130,189,145✔
669
      // Score absorption estimate of keff
670
      if (settings::run_mode == RunMode::EIGENVALUE) {
25,607,931✔
671
        p.keff_tally_absorption() += p.wgt() *
17,779,621✔
672
                                     p.neutron_xs(i_nuclide).nu_fission /
17,779,621✔
673
                                     p.neutron_xs(i_nuclide).absorption;
17,779,621✔
674
      }
675

676
      p.wgt() = 0.0;
25,607,931✔
677
      p.event() = TallyEvent::ABSORB;
25,607,931✔
678
      if (!p.fission()) {
25,607,931✔
679
        p.event_mt() = N_DISAPPEAR;
16,805,243✔
680
      }
681
    }
682
  }
683
}
1,139,868,749✔
684

685
void scatter(Particle& p, int i_nuclide)
1,114,104,079✔
686
{
687
  // copy incoming direction
688
  Direction u_old {p.u()};
1,114,104,079✔
689

690
  // Get pointer to nuclide and grid index/interpolation factor
691
  const auto& nuc {data::nuclides[i_nuclide]};
1,114,104,079✔
692
  const auto& micro {p.neutron_xs(i_nuclide)};
1,114,104,079✔
693
  int i_temp = micro.index_temp;
1,114,104,079✔
694

695
  // For tallying purposes, this routine might be called directly. In that
696
  // case, we need to sample a reaction via the cutoff variable
697
  double cutoff = prn(p.current_seed()) * (micro.total - micro.absorption);
1,114,104,079✔
698
  bool sampled = false;
1,114,104,079✔
699

700
  // Calculate elastic cross section if it wasn't precalculated
701
  if (micro.elastic == CACHE_INVALID) {
1,114,104,079✔
702
    nuc->calculate_elastic_xs(p);
864,641,634✔
703
  }
704

705
  double prob = micro.elastic - micro.thermal;
1,114,104,079✔
706
  if (prob > cutoff) {
1,114,104,079✔
707
    // =======================================================================
708
    // NON-S(A,B) ELASTIC SCATTERING
709

710
    // Determine temperature
711
    double kT = nuc->multipole_ ? p.sqrtkT() * p.sqrtkT() : nuc->kTs_[i_temp];
967,617,453✔
712

713
    // Perform collision physics for elastic scattering
714
    elastic_scatter(i_nuclide, *nuc->reactions_[0], kT, p);
967,617,453✔
715

716
    p.event_mt() = ELASTIC;
967,617,453✔
717
    sampled = true;
967,617,453✔
718
  }
719

720
  prob = micro.elastic;
1,114,104,079✔
721
  if (prob > cutoff && !sampled) {
1,114,104,079✔
722
    // =======================================================================
723
    // S(A,B) SCATTERING
724

725
    sab_scatter(i_nuclide, micro.index_sab, p);
127,478,075✔
726

727
    p.event_mt() = ELASTIC;
127,478,075✔
728
    sampled = true;
127,478,075✔
729
  }
730

731
  if (!sampled) {
1,114,104,079✔
732
    // =======================================================================
733
    // INELASTIC SCATTERING
734

735
    int n = nuc->index_inelastic_scatter_.size();
19,008,551✔
736
    int i = 0;
19,008,551✔
737
    for (int j = 0; j < n && prob < cutoff; ++j) {
355,831,861✔
738
      i = nuc->index_inelastic_scatter_[j];
336,823,310✔
739

740
      // add to cumulative probability
741
      prob += nuc->reactions_[i]->xs(micro);
336,823,310✔
742
    }
743

744
    // Perform collision physics for inelastic scattering
745
    const auto& rx {nuc->reactions_[i]};
19,008,551✔
746
    inelastic_scatter(*nuc, *rx, p);
19,008,551✔
747
    p.event_mt() = rx->mt_;
19,008,551✔
748
  }
749

750
  // Set event component
751
  p.event() = TallyEvent::SCATTER;
1,114,104,079✔
752

753
  // Sample new outgoing angle for isotropic-in-lab scattering
754
  const auto& mat {model::materials[p.material()]};
1,114,104,079!
755
  if (!mat->p0_.empty()) {
1,114,104,079!
756
    int i_nuc_mat = mat->mat_nuclide_index_[i_nuclide];
326,370✔
757
    if (mat->p0_[i_nuc_mat]) {
326,370!
758
      // Sample isotropic-in-lab outgoing direction
759
      p.u() = isotropic_direction(p.current_seed());
326,370✔
760
      p.mu() = u_old.dot(p.u());
326,370✔
761
    }
762
  }
763
}
1,114,104,079✔
764

765
void elastic_scatter(int i_nuclide, const Reaction& rx, double kT, Particle& p)
967,617,453✔
766
{
767
  // get pointer to nuclide
768
  const auto& nuc {data::nuclides[i_nuclide]};
967,617,453✔
769

770
  double vel = std::sqrt(p.E());
967,617,453✔
771
  double awr = nuc->awr_;
967,617,453✔
772

773
  // Neutron velocity in LAB
774
  Direction v_n = vel * p.u();
967,617,453✔
775

776
  // Sample velocity of target nucleus
777
  Direction v_t {};
967,617,453✔
778
  if (!p.neutron_xs(i_nuclide).use_ptable) {
967,617,453✔
779
    v_t = sample_target_velocity(*nuc, p.E(), p.u(), v_n,
926,052,798✔
780
      p.neutron_xs(i_nuclide).elastic, kT, p.current_seed());
926,052,798✔
781
  }
782

783
  // Velocity of center-of-mass
784
  Direction v_cm = (v_n + awr * v_t) / (awr + 1.0);
967,617,453✔
785

786
  // Transform to CM frame
787
  v_n -= v_cm;
967,617,453✔
788

789
  // Find speed of neutron in CM
790
  vel = v_n.norm();
967,617,453✔
791

792
  // Sample scattering angle, checking if angle distribution is present (assume
793
  // isotropic otherwise)
794
  double mu_cm;
967,617,453✔
795
  auto& d = rx.products_[0].distribution_[0];
967,617,453!
796
  auto d_ = dynamic_cast<UncorrelatedAngleEnergy*>(d.get());
967,617,453!
797
  if (!d_->angle().empty()) {
967,617,453!
798
    mu_cm = d_->angle().sample(p.E(), p.current_seed());
967,617,453✔
799
  } else {
800
    mu_cm = uniform_distribution(-1., 1., p.current_seed());
×
801
  }
802

803
  // Determine direction cosines in CM
804
  Direction u_cm = v_n / vel;
967,617,453✔
805

806
  // Rotate neutron velocity vector to new angle -- note that the speed of the
807
  // neutron in CM does not change in elastic scattering. However, the speed
808
  // will change when we convert back to LAB
809
  v_n = vel * rotate_angle(u_cm, mu_cm, nullptr, p.current_seed());
967,617,453✔
810

811
  // Transform back to LAB frame
812
  v_n += v_cm;
967,617,453✔
813

814
  p.E() = v_n.dot(v_n);
967,617,453✔
815
  vel = std::sqrt(p.E());
967,617,453✔
816

817
  // compute cosine of scattering angle in LAB frame by taking dot product of
818
  // neutron's pre- and post-collision angle
819
  p.mu() = p.u().dot(v_n) / vel;
967,617,453✔
820

821
  // Set energy and direction of particle in LAB frame
822
  p.u() = v_n / vel;
967,617,453!
823

824
  // Because of floating-point roundoff, it may be possible for mu_lab to be
825
  // outside of the range [-1,1). In these cases, we just set mu_lab to exactly
826
  // -1 or 1
827
  if (std::abs(p.mu()) > 1.0)
967,617,453!
828
    p.mu() = std::copysign(1.0, p.mu());
×
829
}
967,617,453✔
830

831
void sab_scatter(int i_nuclide, int i_sab, Particle& p)
127,478,075✔
832
{
833
  // Determine temperature index
834
  const auto& micro {p.neutron_xs(i_nuclide)};
127,478,075✔
835
  int i_temp = micro.index_temp_sab;
127,478,075✔
836

837
  // Sample energy and angle
838
  double E_out;
127,478,075✔
839
  data::thermal_scatt[i_sab]->data_[i_temp].sample(
127,478,075✔
840
    micro, p.E(), &E_out, &p.mu(), p.current_seed());
127,478,075✔
841

842
  // Set energy to outgoing, change direction of particle
843
  p.E() = E_out;
127,478,075✔
844
  p.u() = rotate_angle(p.u(), p.mu(), nullptr, p.current_seed());
127,478,075✔
845
}
127,478,075✔
846

847
Direction sample_target_velocity(const Nuclide& nuc, double E, Direction u,
926,052,798✔
848
  Direction v_neut, double xs_eff, double kT, uint64_t* seed)
849
{
850
  // check if nuclide is a resonant scatterer
851
  ResScatMethod sampling_method;
926,052,798✔
852
  if (nuc.resonant_) {
926,052,798✔
853

854
    // sampling method to use
855
    sampling_method = settings::res_scat_method;
84,557✔
856

857
    // upper resonance scattering energy bound (target is at rest above this E)
858
    if (E > settings::res_scat_energy_max) {
84,557✔
859
      return {};
40,755✔
860

861
      // lower resonance scattering energy bound (should be no resonances below)
862
    } else if (E < settings::res_scat_energy_min) {
43,802✔
863
      sampling_method = ResScatMethod::cxs;
864
    }
865

866
    // otherwise, use free gas model
867
  } else {
868
    if (E >= settings::free_gas_threshold * kT && nuc.awr_ > 1.0) {
925,968,241✔
869
      return {};
420,531,468✔
870
    } else {
871
      sampling_method = ResScatMethod::cxs;
872
    }
873
  }
874

875
  // use appropriate target velocity sampling method
876
  switch (sampling_method) {
18,810!
877
  case ResScatMethod::cxs:
505,461,765✔
878

879
    // sample target velocity with the constant cross section (cxs) approx.
880
    return sample_cxs_target_velocity(nuc.awr_, E, u, kT, seed);
505,461,765✔
881

882
  case ResScatMethod::dbrc:
18,810✔
883
  case ResScatMethod::rvs: {
18,810✔
884
    double E_red = std::sqrt(nuc.awr_ * E / kT);
18,810✔
885
    double E_low = std::pow(std::max(0.0, E_red - 4.0), 2) * kT / nuc.awr_;
37,620!
886
    double E_up = (E_red + 4.0) * (E_red + 4.0) * kT / nuc.awr_;
18,810✔
887

888
    // find lower and upper energy bound indices
889
    // lower index
890
    int i_E_low;
18,810✔
891
    if (E_low < nuc.energy_0K_.front()) {
18,810!
892
      i_E_low = 0;
893
    } else if (E_low > nuc.energy_0K_.back()) {
18,810!
894
      i_E_low = nuc.energy_0K_.size() - 2;
×
895
    } else {
896
      i_E_low =
18,810✔
897
        lower_bound_index(nuc.energy_0K_.begin(), nuc.energy_0K_.end(), E_low);
18,810✔
898
    }
899

900
    // upper index
901
    int i_E_up;
18,810✔
902
    if (E_up < nuc.energy_0K_.front()) {
18,810!
903
      i_E_up = 0;
904
    } else if (E_up > nuc.energy_0K_.back()) {
18,810!
905
      i_E_up = nuc.energy_0K_.size() - 2;
×
906
    } else {
907
      i_E_up =
18,810✔
908
        lower_bound_index(nuc.energy_0K_.begin(), nuc.energy_0K_.end(), E_up);
18,810✔
909
    }
910

911
    if (i_E_up == i_E_low) {
18,810✔
912
      // Handle degenerate case -- if the upper/lower bounds occur for the same
913
      // index, then using cxs is probably a good approximation
914
      return sample_cxs_target_velocity(nuc.awr_, E, u, kT, seed);
18,810✔
915
    }
916

917
    if (sampling_method == ResScatMethod::dbrc) {
15,532!
918
      // interpolate xs since we're not exactly at the energy indices
919
      double xs_low = nuc.elastic_0K_[i_E_low];
×
920
      double m = (nuc.elastic_0K_[i_E_low + 1] - xs_low) /
×
921
                 (nuc.energy_0K_[i_E_low + 1] - nuc.energy_0K_[i_E_low]);
×
922
      xs_low += m * (E_low - nuc.energy_0K_[i_E_low]);
×
923
      double xs_up = nuc.elastic_0K_[i_E_up];
×
924
      m = (nuc.elastic_0K_[i_E_up + 1] - xs_up) /
×
925
          (nuc.energy_0K_[i_E_up + 1] - nuc.energy_0K_[i_E_up]);
×
926
      xs_up += m * (E_up - nuc.energy_0K_[i_E_up]);
×
927

928
      // get max 0K xs value over range of practical relative energies
929
      double xs_max = *std::max_element(
×
930
        &nuc.elastic_0K_[i_E_low + 1], &nuc.elastic_0K_[i_E_up + 1]);
×
931
      xs_max = std::max({xs_low, xs_max, xs_up});
×
932

933
      while (true) {
×
934
        double E_rel;
×
935
        Direction v_target;
×
936
        while (true) {
×
937
          // sample target velocity with the constant cross section (cxs)
938
          // approx.
939
          v_target = sample_cxs_target_velocity(nuc.awr_, E, u, kT, seed);
×
940
          Direction v_rel = v_neut - v_target;
×
941
          E_rel = v_rel.dot(v_rel);
×
942
          if (E_rel < E_up)
×
943
            break;
944
        }
945

946
        // perform Doppler broadening rejection correction (dbrc)
947
        double xs_0K = nuc.elastic_xs_0K(E_rel);
×
948
        double R = xs_0K / xs_max;
×
949
        if (prn(seed) < R)
×
950
          return v_target;
×
951
      }
952

953
    } else if (sampling_method == ResScatMethod::rvs) {
15,532✔
954
      // interpolate xs CDF since we're not exactly at the energy indices
955
      // cdf value at lower bound attainable energy
956
      double cdf_low = 0.0;
15,532✔
957
      if (E_low > nuc.energy_0K_.front()) {
15,532!
958
        double m = (nuc.xs_cdf_[i_E_low + 1] - nuc.xs_cdf_[i_E_low]) /
15,532✔
959
                   (nuc.energy_0K_[i_E_low + 1] - nuc.energy_0K_[i_E_low]);
15,532✔
960
        cdf_low = nuc.xs_cdf_[i_E_low] + m * (E_low - nuc.energy_0K_[i_E_low]);
15,532✔
961
      }
962

963
      // cdf value at upper bound attainable energy
964
      double m = (nuc.xs_cdf_[i_E_up + 1] - nuc.xs_cdf_[i_E_up]) /
15,532✔
965
                 (nuc.energy_0K_[i_E_up + 1] - nuc.energy_0K_[i_E_up]);
15,532✔
966
      double cdf_up = nuc.xs_cdf_[i_E_up] + m * (E_up - nuc.energy_0K_[i_E_up]);
15,532✔
967

968
      while (true) {
325,908✔
969
        // directly sample Maxwellian
970
        double E_t = -kT * std::log(prn(seed));
170,720✔
971

972
        // sample a relative energy using the xs cdf
973
        double cdf_rel = cdf_low + prn(seed) * (cdf_up - cdf_low);
170,720✔
974
        int i_E_rel = lower_bound_index(nuc.xs_cdf_.begin() + i_E_low,
170,720✔
975
          nuc.xs_cdf_.begin() + i_E_up + 2, cdf_rel);
170,720✔
976
        double E_rel = nuc.energy_0K_[i_E_low + i_E_rel];
170,720✔
977
        double m = (nuc.xs_cdf_[i_E_low + i_E_rel + 1] -
170,720✔
978
                     nuc.xs_cdf_[i_E_low + i_E_rel]) /
170,720✔
979
                   (nuc.energy_0K_[i_E_low + i_E_rel + 1] -
170,720✔
980
                     nuc.energy_0K_[i_E_low + i_E_rel]);
170,720✔
981
        E_rel += (cdf_rel - nuc.xs_cdf_[i_E_low + i_E_rel]) / m;
170,720✔
982

983
        // perform rejection sampling on cosine between
984
        // neutron and target velocities
985
        double mu = (E_t + nuc.awr_ * (E - E_rel)) /
170,720✔
986
                    (2.0 * std::sqrt(nuc.awr_ * E * E_t));
170,720✔
987

988
        if (std::abs(mu) < 1.0) {
170,720✔
989
          // set and accept target velocity
990
          E_t /= nuc.awr_;
15,532✔
991
          return std::sqrt(E_t) * rotate_angle(u, mu, nullptr, seed);
15,532✔
992
        }
993
      }
155,188✔
994
    }
995
  } // case RVS, DBRC
996
  } // switch (sampling_method)
997

998
  UNREACHABLE();
×
999
}
1000

1001
Direction sample_cxs_target_velocity(
505,465,043✔
1002
  double awr, double E, Direction u, double kT, uint64_t* seed)
1003
{
1004
  double beta_vn = std::sqrt(awr * E / kT);
505,465,043✔
1005
  double alpha = 1.0 / (1.0 + std::sqrt(PI) * beta_vn / 2.0);
505,465,043✔
1006

1007
  double beta_vt_sq;
608,812,749✔
1008
  double mu;
608,812,749✔
1009
  while (true) {
608,812,749✔
1010
    // Sample two random numbers
1011
    double r1 = prn(seed);
608,812,749✔
1012
    double r2 = prn(seed);
608,812,749✔
1013

1014
    if (prn(seed) < alpha) {
608,812,749✔
1015
      // With probability alpha, we sample the distribution p(y) =
1016
      // y*e^(-y). This can be done with sampling scheme C45 from the Monte
1017
      // Carlo sampler
1018

1019
      beta_vt_sq = -std::log(r1 * r2);
167,048,526✔
1020

1021
    } else {
1022
      // With probability 1-alpha, we sample the distribution p(y) = y^2 *
1023
      // e^(-y^2). This can be done with sampling scheme C61 from the Monte
1024
      // Carlo sampler
1025

1026
      double c = std::cos(PI / 2.0 * prn(seed));
441,764,223✔
1027
      beta_vt_sq = -std::log(r1) - std::log(r2) * c * c;
441,764,223✔
1028
    }
1029

1030
    // Determine beta * vt
1031
    double beta_vt = std::sqrt(beta_vt_sq);
608,812,749✔
1032

1033
    // Sample cosine of angle between neutron and target velocity
1034
    mu = uniform_distribution(-1., 1., seed);
608,812,749✔
1035

1036
    // Determine rejection probability
1037
    double accept_prob =
608,812,749✔
1038
      std::sqrt(beta_vn * beta_vn + beta_vt_sq - 2 * beta_vn * beta_vt * mu) /
608,812,749✔
1039
      (beta_vn + beta_vt);
608,812,749✔
1040

1041
    // Perform rejection sampling on vt and mu
1042
    if (prn(seed) < accept_prob)
608,812,749✔
1043
      break;
1044
  }
1045

1046
  // Determine speed of target nucleus
1047
  double vt = std::sqrt(beta_vt_sq * kT / awr);
505,465,043✔
1048

1049
  // Determine velocity vector of target nucleus based on neutron's velocity
1050
  // and the sampled angle between them
1051
  return vt * rotate_angle(u, mu, nullptr, seed);
505,465,043✔
1052
}
1053

1054
void sample_fission_neutron(
30,634,010✔
1055
  int i_nuclide, const Reaction& rx, SourceSite* site, Particle& p)
1056
{
1057
  // Get attributes of particle
1058
  double E_in = p.E();
30,634,010✔
1059
  uint64_t* seed = p.current_seed();
30,634,010✔
1060

1061
  // Determine total nu, delayed nu, and delayed neutron fraction
1062
  const auto& nuc {data::nuclides[i_nuclide]};
30,634,010✔
1063
  double nu_t = nuc->nu(E_in, Nuclide::EmissionMode::total);
30,634,010✔
1064
  double nu_d = nuc->nu(E_in, Nuclide::EmissionMode::delayed);
30,634,010✔
1065
  double beta = nu_d / nu_t;
30,634,010✔
1066

1067
  if (prn(seed) < beta) {
30,634,010✔
1068
    // ====================================================================
1069
    // DELAYED NEUTRON SAMPLED
1070

1071
    // sampled delayed precursor group
1072
    double xi = prn(seed) * nu_d;
192,175✔
1073
    double prob = 0.0;
192,175✔
1074
    int group;
192,175✔
1075
    for (group = 1; group < nuc->n_precursor_; ++group) {
715,699✔
1076
      // determine delayed neutron precursor yield for group j
1077
      double yield = (*rx.products_[group].yield_)(E_in);
702,134✔
1078

1079
      // Check if this group is sampled
1080
      prob += yield;
702,134✔
1081
      if (xi < prob)
702,134✔
1082
        break;
1083
    }
1084

1085
    // if the sum of the probabilities is slightly less than one and the
1086
    // random number is greater, j will be greater than nuc %
1087
    // n_precursor -- check for this condition
1088
    group = std::min(group, nuc->n_precursor_);
192,175!
1089

1090
    // set the delayed group for the particle born from fission
1091
    site->delayed_group = group;
192,175✔
1092

1093
    // Sample time of emission based on decay constant of precursor
1094
    double decay_rate = rx.products_[site->delayed_group].decay_rate_;
192,175✔
1095
    site->time -= std::log(prn(p.current_seed())) / decay_rate;
192,175✔
1096

1097
  } else {
1098
    // ====================================================================
1099
    // PROMPT NEUTRON SAMPLED
1100

1101
    // set the delayed group for the particle born from fission to 0
1102
    site->delayed_group = 0;
30,441,835✔
1103
  }
1104

1105
  // sample from prompt neutron energy distribution
1106
  int n_sample = 0;
1107
  double mu;
30,634,010✔
1108
  while (true) {
30,634,010✔
1109
    rx.products_[site->delayed_group].sample(E_in, site->E, mu, seed);
30,634,010✔
1110

1111
    // resample if energy is greater than maximum neutron energy
1112
    int neutron = ParticleType::neutron().transport_index();
30,634,010✔
1113
    if (site->E < data::energy_max[neutron])
30,634,010!
1114
      break;
1115

1116
    // check for large number of resamples
1117
    ++n_sample;
×
1118
    if (n_sample == MAX_SAMPLE) {
×
1119
      // particle_write_restart(p)
1120
      fatal_error("Resampled energy distribution maximum number of times "
×
1121
                  "for nuclide " +
×
1122
                  nuc->name_);
×
1123
    }
1124
  }
1125

1126
  // Sample azimuthal angle uniformly in [0, 2*pi) and assign angle
1127
  site->u = rotate_angle(p.u(), mu, nullptr, seed);
30,634,010✔
1128
}
30,634,010✔
1129

1130
void inelastic_scatter(const Nuclide& nuc, const Reaction& rx, Particle& p)
19,008,551✔
1131
{
1132
  // copy energy of neutron
1133
  double E_in = p.E();
19,008,551✔
1134

1135
  // sample outgoing energy and scattering cosine
1136
  double E;
19,008,551✔
1137
  double mu;
19,008,551✔
1138
  rx.products_[0].sample(E_in, E, mu, p.current_seed());
19,008,551✔
1139

1140
  // if scattering system is in center-of-mass, transfer cosine of scattering
1141
  // angle and outgoing energy from CM to LAB
1142
  if (rx.scatter_in_cm_) {
19,008,551✔
1143
    double E_cm = E;
18,948,660✔
1144

1145
    // determine outgoing energy in lab
1146
    double A = nuc.awr_;
18,948,660✔
1147
    E = E_cm + (E_in + 2.0 * mu * (A + 1.0) * std::sqrt(E_in * E_cm)) /
18,948,660✔
1148
                 ((A + 1.0) * (A + 1.0));
18,948,660✔
1149

1150
    // determine outgoing angle in lab
1151
    mu = mu * std::sqrt(E_cm / E) + 1.0 / (A + 1.0) * std::sqrt(E_in / E);
18,948,660✔
1152
  }
1153

1154
  // Because of floating-point roundoff, it may be possible for mu to be
1155
  // outside of the range [-1,1). In these cases, we just set mu to exactly -1
1156
  // or 1
1157
  if (std::abs(mu) > 1.0)
19,008,551!
1158
    mu = std::copysign(1.0, mu);
×
1159

1160
  // Set outgoing energy and scattering angle
1161
  p.E() = E;
19,008,551✔
1162
  p.mu() = mu;
19,008,551✔
1163

1164
  // change direction of particle
1165
  p.u() = rotate_angle(p.u(), mu, nullptr, p.current_seed());
19,008,551✔
1166

1167
  // evaluate yield
1168
  double yield = (*rx.products_[0].yield_)(E_in);
19,008,551✔
1169
  if (std::floor(yield) == yield && yield > 0) {
19,008,551!
1170
    // If yield is integral, create exactly that many secondary particles
1171
    for (int i = 0; i < static_cast<int>(std::round(yield)) - 1; ++i) {
19,121,066✔
1172
      p.create_secondary(p.wgt(), p.u(), p.E(), ParticleType::neutron());
112,569✔
1173
    }
1174
  } else {
1175
    // Otherwise, change weight of particle based on yield
1176
    p.wgt() *= yield;
54✔
1177
  }
1178
}
19,008,551✔
1179

1180
void sample_secondary_photons(Particle& p, int i_nuclide)
66,629,112✔
1181
{
1182
  // Sample the number of photons produced
1183
  double y_t =
66,629,112✔
1184
    p.neutron_xs(i_nuclide).photon_prod / p.neutron_xs(i_nuclide).total;
66,629,112✔
1185
  double photon_wgt = p.wgt();
66,629,112✔
1186
  int y = 1;
66,629,112✔
1187

1188
  if (settings::use_decay_photons) {
66,629,112✔
1189
    // For decay photons, sample a single photon and modify the weight
1190
    if (y_t <= 0.0)
72,006✔
1191
      return;
1192
    photon_wgt *= y_t;
54,725✔
1193
  } else {
1194
    // For prompt photons, sample an integral number of photons with weight
1195
    // equal to the neutron's weight
1196
    y = static_cast<int>(y_t);
66,557,106✔
1197
    if (prn(p.current_seed()) <= y_t - y)
66,557,106✔
1198
      ++y;
1,869,263✔
1199
  }
1200

1201
  // Sample each secondary photon
1202
  for (int i = 0; i < y; ++i) {
69,279,727✔
1203
    // Sample the reaction and product
1204
    int i_rx;
2,667,896✔
1205
    int i_product;
2,667,896✔
1206
    sample_photon_product(i_nuclide, p, &i_rx, &i_product);
2,667,896✔
1207

1208
    // Sample the outgoing energy and angle
1209
    auto& rx = data::nuclides[i_nuclide]->reactions_[i_rx];
2,667,896✔
1210
    double E;
2,667,896✔
1211
    double mu;
2,667,896✔
1212
    rx->products_[i_product].sample(p.E(), E, mu, p.current_seed());
2,667,896✔
1213

1214
    // Sample the new direction
1215
    Direction u = rotate_angle(p.u(), mu, nullptr, p.current_seed());
2,667,896✔
1216

1217
    // In a k-eigenvalue simulation, it's necessary to provide higher weight to
1218
    // secondary photons from non-fission reactions to properly balance energy
1219
    // release and deposition. See D. P. Griesheimer, S. J. Douglass, and M. H.
1220
    // Stedry, "Self-consistent energy normalization for quasistatic reactor
1221
    // calculations", Proc. PHYSOR, Cambridge, UK, Mar 29-Apr 2, 2020.
1222
    double wgt = photon_wgt;
2,667,896✔
1223
    if (settings::run_mode == RunMode::EIGENVALUE && !is_fission(rx->mt_)) {
2,667,896✔
1224
      wgt *= simulation::keff;
348,128✔
1225
    }
1226

1227
    // Create the secondary photon
1228
    bool created_photon = p.create_secondary(wgt, u, E, ParticleType::photon());
2,667,896✔
1229

1230
    // Tag secondary particle with parent nuclide
1231
    if (created_photon && settings::use_decay_photons) {
2,667,896✔
1232
      p.local_secondary_bank().back().parent_nuclide =
52,844✔
1233
        rx->products_[i_product].parent_nuclide_;
52,844✔
1234
    }
1235
  }
1236
}
1237

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