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

openmc-dev / openmc / 18554903056

16 Oct 2025 08:18AM UTC coverage: 81.836% (-0.08%) from 81.917%
18554903056

Pull #3547

github

web-flow
Merge 192e6b2c2 into b94b49611
Pull Request #3547: Tally spectrum of secondary particles

16598 of 23152 branches covered (71.69%)

Branch coverage included in aggregate %.

24 of 67 new or added lines in 10 files covered. (35.82%)

30 existing lines in 3 files now uncovered.

53729 of 62784 relevant lines covered (85.58%)

44550947.96 hits per line

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

84.44
/src/particle.cpp
1
#include "openmc/particle.h"
2

3
#include <algorithm> // copy, min
4
#include <cmath>     // log, abs
5

6
#include <fmt/core.h>
7

8
#include "openmc/bank.h"
9
#include "openmc/capi.h"
10
#include "openmc/cell.h"
11
#include "openmc/constants.h"
12
#include "openmc/dagmc.h"
13
#include "openmc/error.h"
14
#include "openmc/geometry.h"
15
#include "openmc/hdf5_interface.h"
16
#include "openmc/material.h"
17
#include "openmc/message_passing.h"
18
#include "openmc/mgxs_interface.h"
19
#include "openmc/nuclide.h"
20
#include "openmc/particle_data.h"
21
#include "openmc/photon.h"
22
#include "openmc/physics.h"
23
#include "openmc/physics_mg.h"
24
#include "openmc/random_lcg.h"
25
#include "openmc/settings.h"
26
#include "openmc/simulation.h"
27
#include "openmc/source.h"
28
#include "openmc/surface.h"
29
#include "openmc/tallies/derivative.h"
30
#include "openmc/tallies/tally.h"
31
#include "openmc/tallies/tally_scoring.h"
32
#include "openmc/track_output.h"
33
#include "openmc/weight_windows.h"
34

35
#ifdef OPENMC_DAGMC_ENABLED
36
#include "DagMC.hpp"
37
#endif
38

39
namespace openmc {
40

41
namespace simulation {
42
thread_local Particle tmp_particle;
43
}
44

45
//==============================================================================
46
// Particle implementation
47
//==============================================================================
48

49
double Particle::speed() const
2,147,483,647✔
50
{
51
  if (settings::run_CE) {
2,147,483,647✔
52
    // Determine mass in eV/c^2
53
    double mass;
54
    switch (this->type()) {
2,147,483,647!
55
    case ParticleType::neutron:
2,073,855,210✔
56
      mass = MASS_NEUTRON_EV;
2,073,855,210✔
57
      break;
2,073,855,210✔
58
    case ParticleType::photon:
21,525,750✔
59
      mass = 0.0;
21,525,750✔
60
      break;
21,525,750✔
61
    case ParticleType::electron:
53,146,649✔
62
    case ParticleType::positron:
63
      mass = MASS_ELECTRON_EV;
53,146,649✔
64
      break;
53,146,649✔
65
    }
66
    // Equivalent to C * sqrt(1-(m/(m+E))^2) without problem at E<<m:
67
    return C_LIGHT * std::sqrt(this->E() * (this->E() + 2 * mass)) /
2,147,483,647✔
68
           (this->E() + mass);
2,147,483,647✔
69
  } else {
70
    auto& macro_xs = data::mg.macro_xs_[this->material()];
2,063,937,414✔
71
    int macro_t = this->mg_xs_cache().t;
2,063,937,414✔
72
    int macro_a = macro_xs.get_angle_index(this->u());
2,063,937,414✔
73
    return 1.0 / macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, this->g(), nullptr,
2,063,937,414✔
74
                   nullptr, nullptr, macro_t, macro_a);
2,063,937,414✔
75
  }
76
}
77

78
bool Particle::create_secondary(
110,860,078✔
79
  double wgt, Direction u, double E, ParticleType type)
80
{
81
  // If energy is below cutoff for this particle, don't create secondary
82
  // particle
83
  if (E < settings::energy_cutoff[static_cast<int>(type)]) {
110,860,078✔
84
    return false;
52,995,193✔
85
  }
86

87
  auto& bank = secondary_bank().emplace_back();
57,864,885✔
88
  bank.particle = type;
57,864,885✔
89
  bank.wgt = wgt;
57,864,885✔
90
  bank.r = r();
57,864,885✔
91
  bank.u = u;
57,864,885✔
92
  bank.E = settings::run_CE ? E : g();
57,864,885!
93
  bank.time = time();
57,864,885✔
94
  bank_second_E() += bank.E;
57,864,885✔
95

96
  // Score tallies affected by secondary particles
97
  if (!model::active_particleout_analog_tallies.empty()) {
57,864,885!
98
    // Create secondary particle for tallying purposes only
NEW
99
    simulation::tmp_particle.from_source(&bank);
×
NEW
100
    simulation::tmp_particle.u_last() = this->u();
×
NEW
101
    simulation::tmp_particle.r_last() = this->r();
×
NEW
102
    simulation::tmp_particle.E_last() = this->E();
×
NEW
103
    simulation::tmp_particle.type_last() = this->type();
×
104

NEW
105
    if (settings::run_CE) {
×
NEW
106
      score_analog_tally_ce(
×
107
        simulation::tmp_particle, model::active_particleout_analog_tallies);
108
    } else {
NEW
109
      score_analog_tally_mg(
×
110
        simulation::tmp_particle, model::active_particleout_analog_tallies);
111
    }
112
  }
113
  return true;
57,864,885✔
114
}
115

116
void Particle::split(double wgt)
4,566,897✔
117
{
118
  auto& bank = secondary_bank().emplace_back();
4,566,897✔
119
  bank.particle = type();
4,566,897✔
120
  bank.wgt = wgt;
4,566,897✔
121
  bank.r = r();
4,566,897✔
122
  bank.u = u();
4,566,897✔
123
  bank.E = settings::run_CE ? E() : g();
4,566,897✔
124
  bank.time = time();
4,566,897✔
125

126
  // Convert signed index to a signed surface ID
127
  if (surface() == SURFACE_NONE) {
4,566,897✔
128
    bank.surf_id = SURFACE_NONE;
4,566,877✔
129
  } else {
130
    int surf_id = model::surfaces[surface_index()]->id_;
20✔
131
    bank.surf_id = (surface() > 0) ? surf_id : -surf_id;
20!
132
  }
133
}
4,566,897✔
134

135
void Particle::from_source(const SourceSite* src)
229,985,062✔
136
{
137
  // Reset some attributes
138
  clear();
229,985,062✔
139
  surface() = SURFACE_NONE;
229,985,062✔
140
  cell_born() = C_NONE;
229,985,062✔
141
  material() = C_NONE;
229,985,062✔
142
  n_collision() = 0;
229,985,062✔
143
  fission() = false;
229,985,062✔
144
  zero_flux_derivs();
229,985,062✔
145
  lifetime() = 0.0;
229,985,062✔
146

147
  // Copy attributes from source bank site
148
  type() = src->particle;
229,985,062✔
149
  type_last() = src->particle;
229,985,062✔
150
  wgt() = src->wgt;
229,985,062✔
151
  wgt_last() = src->wgt;
229,985,062✔
152
  r() = src->r;
229,985,062✔
153
  u() = src->u;
229,985,062✔
154
  r_born() = src->r;
229,985,062✔
155
  r_last_current() = src->r;
229,985,062✔
156
  r_last() = src->r;
229,985,062✔
157
  u_last() = src->u;
229,985,062✔
158
  if (settings::run_CE) {
229,985,062✔
159
    E() = src->E;
114,344,345✔
160
    g() = 0;
114,344,345✔
161
  } else {
162
    g() = static_cast<int>(src->E);
115,640,717✔
163
    g_last() = static_cast<int>(src->E);
115,640,717✔
164
    E() = data::mg.energy_bin_avg_[g()];
115,640,717✔
165
  }
166
  E_last() = E();
229,985,062✔
167
  time() = src->time;
229,985,062✔
168
  time_last() = src->time;
229,985,062✔
169
  parent_nuclide() = src->parent_nuclide;
229,985,062✔
170
  delayed_group() = src->delayed_group;
229,985,062✔
171

172
  // Convert signed surface ID to signed index
173
  if (src->surf_id != SURFACE_NONE) {
229,985,062✔
174
    int index_plus_one = model::surface_map[std::abs(src->surf_id)] + 1;
110,020✔
175
    surface() = (src->surf_id > 0) ? index_plus_one : -index_plus_one;
110,020!
176
  }
177
}
229,985,062✔
178

179
void Particle::event_calculate_xs()
2,147,483,647✔
180
{
181
  // Set the random number stream
182
  stream() = STREAM_TRACKING;
2,147,483,647✔
183

184
  // Store pre-collision particle properties
185
  wgt_last() = wgt();
2,147,483,647✔
186
  E_last() = E();
2,147,483,647✔
187
  type_last() = type();
2,147,483,647✔
188
  u_last() = u();
2,147,483,647✔
189
  r_last() = r();
2,147,483,647✔
190
  time_last() = time();
2,147,483,647✔
191

192
  // Reset event variables
193
  event() = TallyEvent::KILL;
2,147,483,647✔
194
  event_nuclide() = NUCLIDE_NONE;
2,147,483,647✔
195
  event_mt() = REACTION_NONE;
2,147,483,647✔
196

197
  // If the cell hasn't been determined based on the particle's location,
198
  // initiate a search for the current cell. This generally happens at the
199
  // beginning of the history and again for any secondary particles
200
  if (lowest_coord().cell() == C_NONE) {
2,147,483,647✔
201
    if (!exhaustive_find_cell(*this)) {
226,375,335!
202
      mark_as_lost(
×
203
        "Could not find the cell containing particle " + std::to_string(id()));
×
UNCOV
204
      return;
×
205
    }
206

207
    // Set birth cell attribute
208
    if (cell_born() == C_NONE)
226,375,335!
209
      cell_born() = lowest_coord().cell();
226,375,335✔
210

211
    // Initialize last cells from current cell
212
    for (int j = 0; j < n_coord(); ++j) {
470,933,240✔
213
      cell_last(j) = coord(j).cell();
244,557,905✔
214
    }
215
    n_coord_last() = n_coord();
226,375,335✔
216
  }
217

218
  // Write particle track.
219
  if (write_track())
2,147,483,647✔
220
    write_particle_track(*this);
10,814✔
221

222
  if (settings::check_overlaps)
2,147,483,647!
UNCOV
223
    check_cell_overlap(*this);
×
224

225
  // Calculate microscopic and macroscopic cross sections
226
  if (material() != MATERIAL_VOID) {
2,147,483,647✔
227
    if (settings::run_CE) {
2,147,483,647✔
228
      if (material() != material_last() || sqrtkT() != sqrtkT_last() ||
2,147,483,647✔
229
          density_mult() != density_mult_last()) {
433,638,426✔
230
        // If the material is the same as the last material and the
231
        // temperature hasn't changed, we don't need to lookup cross
232
        // sections again.
233
        model::materials[material()]->calculate_xs(*this);
1,619,917,797✔
234
      }
235
    } else {
236
      // Get the MG data; unlike the CE case above, we have to re-calculate
237
      // cross sections for every collision since the cross sections may
238
      // be angle-dependent
239
      data::mg.macro_xs_[material()].calculate_xs(*this);
2,063,937,414✔
240

241
      // Update the particle's group while we know we are multi-group
242
      g_last() = g();
2,063,937,414✔
243
    }
244
  } else {
245
    macro_xs().total = 0.0;
67,534,317✔
246
    macro_xs().absorption = 0.0;
67,534,317✔
247
    macro_xs().fission = 0.0;
67,534,317✔
248
    macro_xs().nu_fission = 0.0;
67,534,317✔
249
  }
250
}
251

252
void Particle::event_advance()
2,147,483,647✔
253
{
254
  // Find the distance to the nearest boundary
255
  boundary() = distance_to_boundary(*this);
2,147,483,647✔
256

257
  // Sample a distance to collision
258
  if (type() == ParticleType::electron || type() == ParticleType::positron) {
2,147,483,647✔
259
    collision_distance() = material() == MATERIAL_VOID ? INFINITY : 0.0;
53,146,649!
260
  } else if (macro_xs().total == 0.0) {
2,147,483,647✔
261
    collision_distance() = INFINITY;
67,534,317✔
262
  } else {
263
    collision_distance() = -std::log(prn(current_seed())) / macro_xs().total;
2,147,483,647✔
264
  }
265

266
  double speed = this->speed();
2,147,483,647✔
267
  double time_cutoff = settings::time_cutoff[static_cast<int>(type())];
2,147,483,647✔
268
  double distance_cutoff =
269
    (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY;
2,147,483,647✔
270

271
  // Select smaller of the three distances
272
  double distance =
273
    std::min({boundary().distance(), collision_distance(), distance_cutoff});
2,147,483,647✔
274

275
  // Advance particle in space and time
276
  this->move_distance(distance);
2,147,483,647✔
277
  double dt = distance / speed;
2,147,483,647✔
278
  this->time() += dt;
2,147,483,647✔
279
  this->lifetime() += dt;
2,147,483,647✔
280

281
  // Score timed track-length tallies
282
  if (!model::active_timed_tracklength_tallies.empty()) {
2,147,483,647✔
283
    score_timed_tracklength_tally(*this, distance);
3,628,317✔
284
  }
285

286
  // Score track-length tallies
287
  if (!model::active_tracklength_tallies.empty()) {
2,147,483,647✔
288
    score_tracklength_tally(*this, distance);
1,569,323,182✔
289
  }
290

291
  // Score track-length estimate of k-eff
292
  if (settings::run_mode == RunMode::EIGENVALUE &&
2,147,483,647✔
293
      type() == ParticleType::neutron) {
2,147,483,647✔
294
    keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission;
2,147,483,647✔
295
  }
296

297
  // Score flux derivative accumulators for differential tallies.
298
  if (!model::active_tallies.empty()) {
2,147,483,647✔
299
    score_track_derivative(*this, distance);
1,780,514,863✔
300
  }
301

302
  // Set particle weight to zero if it hit the time boundary
303
  if (distance == distance_cutoff) {
2,147,483,647✔
304
    wgt() = 0.0;
224,928✔
305
  }
306
}
2,147,483,647✔
307

308
void Particle::event_cross_surface()
2,147,483,647✔
309
{
310
  // Saving previous cell data
311
  for (int j = 0; j < n_coord(); ++j) {
2,147,483,647✔
312
    cell_last(j) = coord(j).cell();
2,147,483,647✔
313
  }
314
  n_coord_last() = n_coord();
2,147,483,647✔
315

316
  // Set surface that particle is on and adjust coordinate levels
317
  surface() = boundary().surface();
2,147,483,647✔
318
  n_coord() = boundary().coord_level();
2,147,483,647✔
319

320
  if (boundary().lattice_translation()[0] != 0 ||
2,147,483,647✔
321
      boundary().lattice_translation()[1] != 0 ||
2,147,483,647✔
322
      boundary().lattice_translation()[2] != 0) {
1,821,410,933✔
323
    // Particle crosses lattice boundary
324

325
    bool verbose = settings::verbosity >= 10 || trace();
687,031,349!
326
    cross_lattice(*this, boundary(), verbose);
687,031,349✔
327
    event() = TallyEvent::LATTICE;
687,031,349✔
328
  } else {
329
    // Particle crosses surface
330
    const auto& surf {model::surfaces[surface_index()].get()};
1,633,951,808✔
331
    // If BC, add particle to surface source before crossing surface
332
    if (surf->surf_source_ && surf->bc_) {
1,633,951,808✔
333
      add_surf_source_to_bank(*this, *surf);
736,393,988✔
334
    }
335
    this->cross_surface(*surf);
1,633,951,808✔
336
    // If no BC, add particle to surface source after crossing surface
337
    if (surf->surf_source_ && !surf->bc_) {
1,633,951,799✔
338
      add_surf_source_to_bank(*this, *surf);
896,319,984✔
339
    }
340
    if (settings::weight_window_checkpoint_surface) {
1,633,951,799✔
341
      apply_weight_windows(*this);
396!
342
    }
343
    event() = TallyEvent::SURFACE;
1,633,951,799✔
344
  }
345
  // Score cell to cell partial currents
346
  if (!model::active_surface_tallies.empty()) {
2,147,483,647✔
347
    score_surface_tally(*this, model::active_surface_tallies);
34,900,767✔
348
  }
349
}
2,147,483,647✔
350

351
void Particle::event_collide()
2,147,483,647✔
352
{
353
  // Score collision estimate of keff
354
  if (settings::run_mode == RunMode::EIGENVALUE &&
2,147,483,647✔
355
      type() == ParticleType::neutron) {
2,147,483,647✔
356
    keff_tally_collision() += wgt() * macro_xs().nu_fission / macro_xs().total;
2,147,483,647✔
357
  }
358

359
  // Score surface current tallies -- this has to be done before the collision
360
  // since the direction of the particle will change and we need to use the
361
  // pre-collision direction to figure out what mesh surfaces were crossed
362

363
  if (!model::active_meshsurf_tallies.empty())
2,147,483,647✔
364
    score_surface_tally(*this, model::active_meshsurf_tallies);
97,100,980✔
365

366
  // Clear surface component
367
  surface() = SURFACE_NONE;
2,147,483,647✔
368

369
  if (settings::run_CE) {
2,147,483,647✔
370
    collision(*this);
849,255,941✔
371
  } else {
372
    collision_mg(*this);
1,783,060,477✔
373
  }
374

375
  // Score collision estimator tallies -- this is done after a collision
376
  // has occurred rather than before because we need information on the
377
  // outgoing energy for any tallies with an outgoing energy filter
378
  if (!model::active_collision_tallies.empty())
2,147,483,647✔
379
    score_collision_tally(*this);
108,606,883✔
380
  if (!model::active_analog_tallies.empty()) {
2,147,483,647✔
381
    if (settings::run_CE) {
142,312,212✔
382
      score_analog_tally_ce(*this, model::active_analog_tallies);
141,103,950✔
383
    } else {
384
      score_analog_tally_mg(*this, model::active_analog_tallies);
1,208,262✔
385
    }
386
  }
387

388
  if (!model::active_pulse_height_tallies.empty() &&
2,147,483,647✔
389
      type() == ParticleType::photon) {
16,918✔
390
    pht_collision_energy();
2,024✔
391
  }
392

393
  // Reset banked weight during collision
394
  n_bank() = 0;
2,147,483,647✔
395
  bank_second_E() = 0.0;
2,147,483,647✔
396
  wgt_bank() = 0.0;
2,147,483,647✔
397
  zero_delayed_bank();
2,147,483,647✔
398

399
  // Reset fission logical
400
  fission() = false;
2,147,483,647✔
401

402
  // Save coordinates for tallying purposes
403
  r_last_current() = r();
2,147,483,647✔
404

405
  // Set last material to none since cross sections will need to be
406
  // re-evaluated
407
  material_last() = C_NONE;
2,147,483,647✔
408

409
  // Set all directions to base level -- right now, after a collision, only
410
  // the base level directions are changed
411
  for (int j = 0; j < n_coord() - 1; ++j) {
2,147,483,647✔
412
    if (coord(j + 1).rotated()) {
131,044,994✔
413
      // If next level is rotated, apply rotation matrix
414
      const auto& m {model::cells[coord(j).cell()]->rotation_};
10,426,614✔
415
      const auto& u {coord(j).u()};
10,426,614✔
416
      coord(j + 1).u() = u.rotate(m);
10,426,614✔
417
    } else {
418
      // Otherwise, copy this level's direction
419
      coord(j + 1).u() = coord(j).u();
120,618,380✔
420
    }
421
  }
422

423
  // Score flux derivative accumulators for differential tallies.
424
  if (!model::active_tallies.empty())
2,147,483,647✔
425
    score_collision_derivative(*this);
738,319,471✔
426

427
#ifdef OPENMC_DAGMC_ENABLED
428
  history().reset();
288,795,050✔
429
#endif
430
}
2,147,483,647✔
431

432
void Particle::event_revive_from_secondary()
2,147,483,647✔
433
{
434
  // If particle has too many events, display warning and kill it
435
  ++n_event();
2,147,483,647✔
436
  if (n_event() == settings::max_particle_events) {
2,147,483,647!
UNCOV
437
    warning("Particle " + std::to_string(id()) +
×
438
            " underwent maximum number of events.");
UNCOV
439
    wgt() = 0.0;
×
440
  }
441

442
  // Check for secondary particles if this particle is dead
443
  if (!alive()) {
2,147,483,647✔
444
    // Write final position for this particle
445
    if (write_track()) {
226,374,931✔
446
      write_particle_track(*this);
6,674✔
447
    }
448

449
    // If no secondary particles, break out of event loop
450
    if (secondary_bank().empty())
226,374,931✔
451
      return;
163,604,418✔
452

453
    from_source(&secondary_bank().back());
62,770,513✔
454
    secondary_bank().pop_back();
62,770,513✔
455
    n_event() = 0;
62,770,513✔
456
    bank_second_E() = 0.0;
62,770,513✔
457

458
    // Subtract secondary particle energy from interim pulse-height results
459
    if (!model::active_pulse_height_tallies.empty() &&
62,786,012✔
460
        this->type() == ParticleType::photon) {
15,499✔
461
      // Since the birth cell of the particle has not been set we
462
      // have to determine it before the energy of the secondary particle can be
463
      // removed from the pulse-height of this cell.
464
      if (lowest_coord().cell() == C_NONE) {
605!
465
        bool verbose = settings::verbosity >= 10 || trace();
605!
466
        if (!exhaustive_find_cell(*this, verbose)) {
605!
467
          mark_as_lost("Could not find the cell containing particle " +
×
468
                       std::to_string(id()));
×
UNCOV
469
          return;
×
470
        }
471
        // Set birth cell attribute
472
        if (cell_born() == C_NONE)
605!
473
          cell_born() = lowest_coord().cell();
605✔
474

475
        // Initialize last cells from current cell
476
        for (int j = 0; j < n_coord(); ++j) {
1,210✔
477
          cell_last(j) = coord(j).cell();
605✔
478
        }
479
        n_coord_last() = n_coord();
605✔
480
      }
481
      pht_secondary_particles();
605✔
482
    }
483

484
    // Enter new particle in particle track file
485
    if (write_track())
62,770,513✔
486
      add_particle_track(*this);
5,604✔
487
  }
488
}
489

490
void Particle::event_death()
163,605,418✔
491
{
492
#ifdef OPENMC_DAGMC_ENABLED
493
  history().reset();
16,997,402✔
494
#endif
495

496
  // Finish particle track output.
497
  if (write_track()) {
163,605,418✔
498
    finalize_particle_track(*this);
1,070✔
499
  }
500

501
// Contribute tally reduction variables to global accumulator
502
#pragma omp atomic
90,637,927✔
503
  global_tally_absorption += keff_tally_absorption();
163,605,418✔
504
#pragma omp atomic
91,057,760✔
505
  global_tally_collision += keff_tally_collision();
163,605,418✔
506
#pragma omp atomic
90,406,102✔
507
  global_tally_tracklength += keff_tally_tracklength();
163,605,418✔
508
#pragma omp atomic
90,108,009✔
509
  global_tally_leakage += keff_tally_leakage();
163,605,418✔
510

511
  // Reset particle tallies once accumulated
512
  keff_tally_absorption() = 0.0;
163,605,418✔
513
  keff_tally_collision() = 0.0;
163,605,418✔
514
  keff_tally_tracklength() = 0.0;
163,605,418✔
515
  keff_tally_leakage() = 0.0;
163,605,418✔
516

517
  if (!model::active_pulse_height_tallies.empty()) {
163,605,418✔
518
    score_pulse_height_tally(*this, model::active_pulse_height_tallies);
5,500✔
519
  }
520

521
  // Record the number of progeny created by this particle.
522
  // This data will be used to efficiently sort the fission bank.
523
  if (settings::run_mode == RunMode::EIGENVALUE) {
163,605,418✔
524
    int64_t offset = id() - 1 - simulation::work_index[mpi::rank];
138,212,000✔
525
    simulation::progeny_per_particle[offset] = n_progeny();
138,212,000✔
526
  }
527
}
163,605,418✔
528

529
void Particle::pht_collision_energy()
2,024✔
530
{
531
  // Adds the energy particles lose in a collision to the pulse-height
532

533
  // determine index of cell in pulse_height_cells
534
  auto it = std::find(model::pulse_height_cells.begin(),
2,024✔
535
    model::pulse_height_cells.end(), lowest_coord().cell());
2,024✔
536

537
  if (it != model::pulse_height_cells.end()) {
2,024!
538
    int index = std::distance(model::pulse_height_cells.begin(), it);
2,024✔
539
    pht_storage()[index] += E_last() - E();
2,024✔
540

541
    // If the energy of the particle is below the cutoff, it will not be sampled
542
    // so its energy is added to the pulse-height in the cell
543
    int photon = static_cast<int>(ParticleType::photon);
2,024✔
544
    if (E() < settings::energy_cutoff[photon]) {
2,024✔
545
      pht_storage()[index] += E();
825✔
546
    }
547
  }
548
}
2,024✔
549

550
void Particle::pht_secondary_particles()
605✔
551
{
552
  // Removes the energy of secondary produced particles from the pulse-height
553

554
  // determine index of cell in pulse_height_cells
555
  auto it = std::find(model::pulse_height_cells.begin(),
605✔
556
    model::pulse_height_cells.end(), cell_born());
605✔
557

558
  if (it != model::pulse_height_cells.end()) {
605!
559
    int index = std::distance(model::pulse_height_cells.begin(), it);
605✔
560
    pht_storage()[index] -= E();
605✔
561
  }
562
}
605✔
563

564
void Particle::cross_surface(const Surface& surf)
1,635,703,114✔
565
{
566

567
  if (settings::verbosity >= 10 || trace()) {
1,635,703,114✔
568
    write_message(1, "    Crossing surface {}", surf.id_);
33✔
569
  }
570

571
// if we're crossing a CSG surface, make sure the DAG history is reset
572
#ifdef OPENMC_DAGMC_ENABLED
573
  if (surf.geom_type() == GeometryType::CSG)
260,187,259✔
574
    history().reset();
260,142,482✔
575
#endif
576

577
  // Handle any applicable boundary conditions.
578
  if (surf.bc_ && settings::run_mode != RunMode::PLOTTING &&
2,147,483,647!
579
      settings::run_mode != RunMode::VOLUME) {
736,931,464✔
580
    surf.bc_->handle_particle(*this, surf);
736,746,096✔
581
    return;
736,746,096✔
582
  }
583

584
  // ==========================================================================
585
  // SEARCH NEIGHBOR LISTS FOR NEXT CELL
586

587
#ifdef OPENMC_DAGMC_ENABLED
588
  // in DAGMC, we know what the next cell should be
589
  if (surf.geom_type() == GeometryType::DAG) {
152,981,531✔
590
    int32_t i_cell = next_cell(surface_index(), cell_last(n_coord() - 1),
37,445✔
591
                       lowest_coord().universe()) -
37,445✔
592
                     1;
37,445✔
593
    // save material, temperature, and density multiplier
594
    material_last() = material();
37,445✔
595
    sqrtkT_last() = sqrtkT();
37,445✔
596
    density_mult_last() = density_mult();
37,445✔
597
    // set new cell value
598
    lowest_coord().cell() = i_cell;
37,445✔
599
    auto& cell = model::cells[i_cell];
37,445✔
600

601
    cell_instance() = 0;
37,445✔
602
    if (cell->distribcell_index_ >= 0)
37,445✔
603
      cell_instance() = cell_instance_at_level(*this, n_coord() - 1);
36,421✔
604

605
    material() = cell->material(cell_instance());
37,445✔
606
    sqrtkT() = cell->sqrtkT(cell_instance());
37,445✔
607
    density_mult() = cell->density_mult(cell_instance());
37,445✔
608
    return;
37,445✔
609
  }
610
#endif
611

612
  bool verbose = settings::verbosity >= 10 || trace();
898,919,573!
613
  if (neighbor_list_find_cell(*this, verbose)) {
898,919,573✔
614
    return;
898,889,662✔
615
  }
616

617
  // ==========================================================================
618
  // COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS
619

620
  // Remove lower coordinate levels
621
  n_coord() = 1;
29,911✔
622
  bool found = exhaustive_find_cell(*this, verbose);
29,911✔
623

624
  if (settings::run_mode != RunMode::PLOTTING && (!found)) {
29,911!
625
    // If a cell is still not found, there are two possible causes: 1) there is
626
    // a void in the model, and 2) the particle hit a surface at a tangent. If
627
    // the particle is really traveling tangent to a surface, if we move it
628
    // forward a tiny bit it should fix the problem.
629

630
    surface() = SURFACE_NONE;
5,799✔
631
    n_coord() = 1;
5,799✔
632
    r() += TINY_BIT * u();
5,799✔
633

634
    // Couldn't find next cell anywhere! This probably means there is an actual
635
    // undefined region in the geometry.
636

637
    if (!exhaustive_find_cell(*this, verbose)) {
5,799!
638
      mark_as_lost("After particle " + std::to_string(id()) +
17,388✔
639
                   " crossed surface " + std::to_string(surf.id_) +
23,178✔
640
                   " it could not be located in any cell and it did not leak.");
641
      return;
5,790✔
642
    }
643
  }
644
}
645

646
void Particle::cross_vacuum_bc(const Surface& surf)
33,218,983✔
647
{
648
  // Score any surface current tallies -- note that the particle is moved
649
  // forward slightly so that if the mesh boundary is on the surface, it is
650
  // still processed
651

652
  if (!model::active_meshsurf_tallies.empty()) {
33,218,983✔
653
    // TODO: Find a better solution to score surface currents than
654
    // physically moving the particle forward slightly
655

656
    r() += TINY_BIT * u();
1,387,174✔
657
    score_surface_tally(*this, model::active_meshsurf_tallies);
1,387,174✔
658
  }
659

660
  // Score to global leakage tally
661
  keff_tally_leakage() += wgt();
33,218,983✔
662

663
  // Kill the particle
664
  wgt() = 0.0;
33,218,983✔
665

666
  // Display message
667
  if (settings::verbosity >= 10 || trace()) {
33,218,983!
668
    write_message(1, "    Leaked out of surface {}", surf.id_);
11✔
669
  }
670
}
33,218,983✔
671

672
void Particle::cross_reflective_bc(const Surface& surf, Direction new_u)
703,870,956✔
673
{
674
  // Do not handle reflective boundary conditions on lower universes
675
  if (n_coord() != 1) {
703,870,956!
UNCOV
676
    mark_as_lost("Cannot reflect particle " + std::to_string(id()) +
×
677
                 " off surface in a lower universe.");
UNCOV
678
    return;
×
679
  }
680

681
  // Score surface currents since reflection causes the direction of the
682
  // particle to change. For surface filters, we need to score the tallies
683
  // twice, once before the particle's surface attribute has changed and
684
  // once after. For mesh surface filters, we need to artificially move
685
  // the particle slightly back in case the surface crossing is coincident
686
  // with a mesh boundary
687

688
  if (!model::active_surface_tallies.empty()) {
703,870,956✔
689
    score_surface_tally(*this, model::active_surface_tallies);
285,021✔
690
  }
691

692
  if (!model::active_meshsurf_tallies.empty()) {
703,870,956✔
693
    Position r {this->r()};
72,147,011✔
694
    this->r() -= TINY_BIT * u();
72,147,011✔
695
    score_surface_tally(*this, model::active_meshsurf_tallies);
72,147,011✔
696
    this->r() = r;
72,147,011✔
697
  }
698

699
  // Set the new particle direction
700
  u() = new_u;
703,870,956✔
701

702
  // Reassign particle's cell and surface
703
  coord(0).cell() = cell_last(0);
703,870,956✔
704
  surface() = -surface();
703,870,956✔
705

706
  // If a reflective surface is coincident with a lattice or universe
707
  // boundary, it is necessary to redetermine the particle's coordinates in
708
  // the lower universes.
709
  // (unless we're using a dagmc model, which has exactly one universe)
710
  n_coord() = 1;
703,870,956✔
711
  if (surf.geom_type() != GeometryType::DAG &&
1,407,739,154!
712
      !neighbor_list_find_cell(*this)) {
703,868,198!
713
    mark_as_lost("Couldn't find particle after reflecting from surface " +
×
714
                 std::to_string(surf.id_) + ".");
×
UNCOV
715
    return;
×
716
  }
717

718
  // Set previous coordinate going slightly past surface crossing
719
  r_last_current() = r() + TINY_BIT * u();
703,870,956✔
720

721
  // Diagnostic message
722
  if (settings::verbosity >= 10 || trace()) {
703,870,956!
UNCOV
723
    write_message(1, "    Reflected from surface {}", surf.id_);
×
724
  }
725
}
726

727
void Particle::cross_periodic_bc(
661,623✔
728
  const Surface& surf, Position new_r, Direction new_u, int new_surface)
729
{
730
  // Do not handle periodic boundary conditions on lower universes
731
  if (n_coord() != 1) {
661,623!
732
    mark_as_lost(
×
UNCOV
733
      "Cannot transfer particle " + std::to_string(id()) +
×
734
      " across surface in a lower universe. Boundary conditions must be "
735
      "applied to root universe.");
UNCOV
736
    return;
×
737
  }
738

739
  // Score surface currents since reflection causes the direction of the
740
  // particle to change -- artificially move the particle slightly back in
741
  // case the surface crossing is coincident with a mesh boundary
742
  if (!model::active_meshsurf_tallies.empty()) {
661,623!
743
    Position r {this->r()};
×
744
    this->r() -= TINY_BIT * u();
×
745
    score_surface_tally(*this, model::active_meshsurf_tallies);
×
UNCOV
746
    this->r() = r;
×
747
  }
748

749
  // Adjust the particle's location and direction.
750
  r() = new_r;
661,623✔
751
  u() = new_u;
661,623✔
752

753
  // Reassign particle's surface
754
  surface() = new_surface;
661,623✔
755

756
  // Figure out what cell particle is in now
757
  n_coord() = 1;
661,623✔
758

759
  if (!neighbor_list_find_cell(*this)) {
661,623!
760
    mark_as_lost("Couldn't find particle after hitting periodic "
×
761
                 "boundary on surface " +
×
UNCOV
762
                 std::to_string(surf.id_) +
×
763
                 ". The normal vector "
764
                 "of one periodic surface may need to be reversed.");
UNCOV
765
    return;
×
766
  }
767

768
  // Set previous coordinate going slightly past surface crossing
769
  r_last_current() = r() + TINY_BIT * u();
661,623✔
770

771
  // Diagnostic message
772
  if (settings::verbosity >= 10 || trace()) {
661,623!
UNCOV
773
    write_message(1, "    Hit periodic boundary on surface {}", surf.id_);
×
774
  }
775
}
776

777
void Particle::mark_as_lost(const char* message)
5,799✔
778
{
779
  // Print warning and write lost particle file
780
  warning(message);
5,799✔
781
  if (settings::max_write_lost_particles < 0 ||
5,799✔
782
      simulation::n_lost_particles < settings::max_write_lost_particles) {
5,500✔
783
    write_restart();
379✔
784
  }
785
  // Increment number of lost particles
786
  wgt() = 0.0;
5,799✔
787
#pragma omp atomic
3,154✔
788
  simulation::n_lost_particles += 1;
2,645✔
789

790
  // Count the total number of simulated particles (on this processor)
791
  auto n = simulation::current_batch * settings::gen_per_batch *
5,799✔
792
           simulation::work_per_rank;
793

794
  // Abort the simulation if the maximum number of lost particles has been
795
  // reached
796
  if (simulation::n_lost_particles >= settings::max_lost_particles &&
5,799✔
797
      simulation::n_lost_particles >= settings::rel_max_lost_particles * n) {
9!
798
    fatal_error("Maximum number of lost particles has been reached.");
9✔
799
  }
800
}
5,790✔
801

802
void Particle::write_restart() const
379✔
803
{
804
  // Dont write another restart file if in particle restart mode
805
  if (settings::run_mode == RunMode::PARTICLE)
379✔
806
    return;
22✔
807

808
  // Set up file name
809
  auto filename = fmt::format("{}particle_{}_{}.h5", settings::path_output,
810
    simulation::current_batch, id());
665✔
811

812
#pragma omp critical(WriteParticleRestart)
374✔
813
  {
814
    // Create file
815
    hid_t file_id = file_open(filename, 'w');
357✔
816

817
    // Write filetype and version info
818
    write_attribute(file_id, "filetype", "particle restart");
357✔
819
    write_attribute(file_id, "version", VERSION_PARTICLE_RESTART);
357✔
820
    write_attribute(file_id, "openmc_version", VERSION);
357✔
821
#ifdef GIT_SHA1
822
    write_attr_string(file_id, "git_sha1", GIT_SHA1);
823
#endif
824

825
    // Write data to file
826
    write_dataset(file_id, "current_batch", simulation::current_batch);
357✔
827
    write_dataset(file_id, "generations_per_batch", settings::gen_per_batch);
357✔
828
    write_dataset(file_id, "current_generation", simulation::current_gen);
357✔
829
    write_dataset(file_id, "n_particles", settings::n_particles);
357✔
830
    switch (settings::run_mode) {
357!
831
    case RunMode::FIXED_SOURCE:
225✔
832
      write_dataset(file_id, "run_mode", "fixed source");
225✔
833
      break;
225✔
834
    case RunMode::EIGENVALUE:
132✔
835
      write_dataset(file_id, "run_mode", "eigenvalue");
132✔
836
      break;
132✔
837
    case RunMode::PARTICLE:
×
838
      write_dataset(file_id, "run_mode", "particle restart");
×
839
      break;
×
840
    default:
×
UNCOV
841
      break;
×
842
    }
843
    write_dataset(file_id, "id", id());
357✔
844
    write_dataset(file_id, "type", static_cast<int>(type()));
357✔
845

846
    int64_t i = current_work();
357✔
847
    if (settings::run_mode == RunMode::EIGENVALUE) {
357✔
848
      // take source data from primary bank for eigenvalue simulation
849
      write_dataset(file_id, "weight", simulation::source_bank[i - 1].wgt);
132✔
850
      write_dataset(file_id, "energy", simulation::source_bank[i - 1].E);
132✔
851
      write_dataset(file_id, "xyz", simulation::source_bank[i - 1].r);
132✔
852
      write_dataset(file_id, "uvw", simulation::source_bank[i - 1].u);
132✔
853
      write_dataset(file_id, "time", simulation::source_bank[i - 1].time);
132✔
854
    } else if (settings::run_mode == RunMode::FIXED_SOURCE) {
225!
855
      // re-sample using rng random number seed used to generate source particle
856
      int64_t id = (simulation::total_gen + overall_generation() - 1) *
225✔
857
                     settings::n_particles +
225✔
858
                   simulation::work_index[mpi::rank] + i;
225✔
859
      uint64_t seed = init_seed(id, STREAM_SOURCE);
225✔
860
      // re-sample source site
861
      auto site = sample_external_source(&seed);
225✔
862
      write_dataset(file_id, "weight", site.wgt);
225✔
863
      write_dataset(file_id, "energy", site.E);
225✔
864
      write_dataset(file_id, "xyz", site.r);
225✔
865
      write_dataset(file_id, "uvw", site.u);
225✔
866
      write_dataset(file_id, "time", site.time);
225✔
867
    }
868

869
    // Close file
870
    file_close(file_id);
357✔
871
  } // #pragma omp critical
872
}
357✔
873

874
void Particle::update_neutron_xs(
2,147,483,647✔
875
  int i_nuclide, int i_grid, int i_sab, double sab_frac, double ncrystal_xs)
876
{
877
  // Get microscopic cross section cache
878
  auto& micro = this->neutron_xs(i_nuclide);
2,147,483,647✔
879

880
  // If the cache doesn't match, recalculate micro xs
881
  if (this->E() != micro.last_E || this->sqrtkT() != micro.last_sqrtkT ||
2,147,483,647✔
882
      i_sab != micro.index_sab || sab_frac != micro.sab_frac ||
2,147,483,647✔
883
      ncrystal_xs != micro.ncrystal_xs) {
2,147,483,647!
884
    data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, *this);
2,147,483,647✔
885

886
    // If NCrystal is being used, update micro cross section cache
887
    micro.ncrystal_xs = ncrystal_xs;
2,147,483,647✔
888
    if (ncrystal_xs >= 0.0) {
2,147,483,647✔
889
      data::nuclides[i_nuclide]->calculate_elastic_xs(*this);
11,018,953✔
890
      ncrystal_update_micro(ncrystal_xs, micro);
11,018,953✔
891
    }
892
  }
893
}
2,147,483,647✔
894

895
//==============================================================================
896
// Non-method functions
897
//==============================================================================
898

899
std::string particle_type_to_str(ParticleType type)
3,901,635✔
900
{
901
  switch (type) {
3,901,635!
902
  case ParticleType::neutron:
2,785,668✔
903
    return "neutron";
2,785,668✔
904
  case ParticleType::photon:
1,115,703✔
905
    return "photon";
1,115,703✔
906
  case ParticleType::electron:
132✔
907
    return "electron";
132✔
908
  case ParticleType::positron:
132✔
909
    return "positron";
132✔
910
  }
UNCOV
911
  UNREACHABLE();
×
912
}
913

914
ParticleType str_to_particle_type(std::string str)
3,596,352✔
915
{
916
  if (str == "neutron") {
3,596,352✔
917
    return ParticleType::neutron;
827,815✔
918
  } else if (str == "photon") {
2,768,537✔
919
    return ParticleType::photon;
2,768,451✔
920
  } else if (str == "electron") {
86✔
921
    return ParticleType::electron;
43✔
922
  } else if (str == "positron") {
43!
923
    return ParticleType::positron;
43✔
924
  } else {
UNCOV
925
    throw std::invalid_argument {fmt::format("Invalid particle name: {}", str)};
×
926
  }
927
}
928

929
void add_surf_source_to_bank(Particle& p, const Surface& surf)
1,632,713,972✔
930
{
931
  if (simulation::current_batch <= settings::n_inactive ||
2,147,483,647✔
932
      simulation::surf_source_bank.full()) {
1,294,295,598✔
933
    return;
1,632,584,319✔
934
  }
935

936
  // If a cell/cellfrom/cellto parameter is defined
937
  if (settings::ssw_cell_id != C_NONE) {
337,085✔
938

939
    // Retrieve cell index and storage type
940
    int cell_idx = model::cell_map[settings::ssw_cell_id];
254,440✔
941

942
    if (surf.bc_) {
254,440✔
943
      // Leave if cellto with vacuum boundary condition
944
      if (surf.bc_->type() == "vacuum" &&
182,560!
945
          settings::ssw_cell_type == SSWCellType::To) {
33,100✔
946
        return;
12,136✔
947
      }
948

949
      // Leave if other boundary condition than vacuum
950
      if (surf.bc_->type() != "vacuum") {
137,324✔
951
        return;
116,360✔
952
      }
953
    }
954

955
    // Check if the cell of interest has been exited
956
    bool exited = false;
125,944✔
957
    for (int i = 0; i < p.n_coord_last(); ++i) {
333,677✔
958
      if (p.cell_last(i) == cell_idx) {
207,733✔
959
        exited = true;
73,765✔
960
      }
961
    }
962

963
    // Check if the cell of interest has been entered
964
    bool entered = false;
125,944✔
965
    for (int i = 0; i < p.n_coord(); ++i) {
297,979✔
966
      if (p.coord(i).cell() == cell_idx) {
172,035✔
967
        entered = true;
57,517✔
968
      }
969
    }
970

971
    // Vacuum boundary conditions: return if cell is not exited
972
    if (surf.bc_) {
125,944✔
973
      if (surf.bc_->type() == "vacuum" && !exited) {
20,964!
974
        return;
14,664✔
975
      }
976
    } else {
977

978
      // If we both enter and exit the cell of interest
979
      if (entered && exited) {
104,980✔
980
        return;
27,203✔
981
      }
982

983
      // If we did not enter nor exit the cell of interest
984
      if (!entered && !exited) {
77,777✔
985
        return;
13,501✔
986
      }
987

988
      // If cellfrom and the cell before crossing is not the cell of
989
      // interest
990
      if (settings::ssw_cell_type == SSWCellType::From && !exited) {
64,276✔
991
        return;
11,543✔
992
      }
993

994
      // If cellto and the cell after crossing is not the cell of interest
995
      if (settings::ssw_cell_type == SSWCellType::To && !entered) {
52,733✔
996
        return;
12,025✔
997
      }
998
    }
999
  }
1000

1001
  SourceSite site;
129,653✔
1002
  site.r = p.r();
129,653✔
1003
  site.u = p.u();
129,653✔
1004
  site.E = p.E();
129,653✔
1005
  site.time = p.time();
129,653✔
1006
  site.wgt = p.wgt();
129,653✔
1007
  site.delayed_group = p.delayed_group();
129,653✔
1008
  site.surf_id = surf.id_;
129,653✔
1009
  site.particle = p.type();
129,653✔
1010
  site.parent_id = p.id();
129,653✔
1011
  site.progeny_id = p.n_progeny();
129,653✔
1012
  int64_t idx = simulation::surf_source_bank.thread_safe_append(site);
129,653✔
1013
}
1014

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