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

openmc-dev / openmc / 23919792556

02 Apr 2026 08:10PM UTC coverage: 81.336% (+0.01%) from 81.324%
23919792556

Pull #3903

github

web-flow
Merge 05ee6d675 into d9b30bbbd
Pull Request #3903: Add from_bounding_box classmethod to structured mesh classes (Issue #3871)

17622 of 25463 branches covered (69.21%)

Branch coverage included in aggregate %.

33 of 40 new or added lines in 1 file covered. (82.5%)

69 existing lines in 3 files now uncovered.

58119 of 67658 relevant lines covered (85.9%)

44576082.99 hits per line

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

86.17
/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/collision_track.h"
12
#include "openmc/constants.h"
13
#include "openmc/dagmc.h"
14
#include "openmc/error.h"
15
#include "openmc/geometry.h"
16
#include "openmc/hdf5_interface.h"
17
#include "openmc/lattice.h"
18
#include "openmc/material.h"
19
#include "openmc/message_passing.h"
20
#include "openmc/mgxs_interface.h"
21
#include "openmc/nuclide.h"
22
#include "openmc/particle_data.h"
23
#include "openmc/photon.h"
24
#include "openmc/physics.h"
25
#include "openmc/physics_mg.h"
26
#include "openmc/random_lcg.h"
27
#include "openmc/settings.h"
28
#include "openmc/simulation.h"
29
#include "openmc/source.h"
30
#include "openmc/surface.h"
31
#include "openmc/tallies/derivative.h"
32
#include "openmc/tallies/tally.h"
33
#include "openmc/tallies/tally_scoring.h"
34
#include "openmc/track_output.h"
35
#include "openmc/weight_windows.h"
36

37
#ifdef OPENMC_DAGMC_ENABLED
38
#include "DagMC.hpp"
39
#endif
40

41
namespace openmc {
42

43
//==============================================================================
44
// Particle implementation
45
//==============================================================================
46

47
double Particle::speed() const
2,147,483,647✔
48
{
49
  if (settings::run_CE) {
2,147,483,647✔
50
    // Determine mass in eV/c^2
51
    double mass = this->mass();
2,147,483,647✔
52

53
    // Equivalent to C * sqrt(1-(m/(m+E))^2) without problem at E<<m:
54
    return C_LIGHT * std::sqrt(this->E() * (this->E() + 2 * mass)) /
2,147,483,647✔
55
           (this->E() + mass);
2,147,483,647✔
56
  } else {
57
    auto mat = this->material();
2,063,937,744✔
58
    if (mat == MATERIAL_VOID)
2,063,937,744!
UNCOV
59
      return 1.0 / data::mg.default_inverse_velocity_[this->g()];
×
60
    auto& macro_xs = data::mg.macro_xs_[mat];
2,063,937,744✔
61
    int macro_t = this->mg_xs_cache().t;
2,063,937,744✔
62
    int macro_a = macro_xs.get_angle_index(this->u());
2,063,937,744✔
63
    return 1.0 / macro_xs.get_xs(
2,147,483,647✔
64
                   MgxsType::INVERSE_VELOCITY, this->g(), macro_t, macro_a);
2,063,937,744✔
65
  }
66
}
67

68
double Particle::mass() const
2,147,483,647✔
69
{
70
  switch (type().pdg_number()) {
2,147,483,647✔
71
  case PDG_NEUTRON:
72
    return MASS_NEUTRON_EV;
73
  case PDG_ELECTRON:
54,141,467✔
74
  case PDG_POSITRON:
54,141,467✔
75
    return MASS_ELECTRON_EV;
54,141,467✔
76
  default:
22,300,322✔
77
    return this->type().mass() * AMU_EV;
22,300,322✔
78
  }
79
}
80

81
bool Particle::create_secondary(
112,642,192✔
82
  double wgt, Direction u, double E, ParticleType type)
83
{
84
  // If energy is below cutoff for this particle, don't create secondary
85
  // particle
86
  int idx = type.transport_index();
112,642,192✔
87
  if (idx == C_NONE) {
112,642,192!
88
    return false;
89
  }
90
  if (E < settings::energy_cutoff[idx]) {
112,642,192✔
91
    return false;
92
  }
93

94
  // Increment number of secondaries created (for ParticleProductionFilter)
95
  n_secondaries()++;
58,894,055✔
96

97
  auto& bank = secondary_bank().emplace_back();
58,894,055✔
98
  bank.particle = type;
58,894,055✔
99
  bank.wgt = wgt;
58,894,055✔
100
  bank.r = r();
58,894,055!
101
  bank.u = u;
58,894,055✔
102
  bank.E = settings::run_CE ? E : g();
58,894,055!
103
  bank.time = time();
58,894,055✔
104
  bank_second_E() += bank.E;
58,894,055✔
105
  return true;
58,894,055✔
106
}
107

108
void Particle::split(double wgt)
4,250,141✔
109
{
110
  auto& bank = secondary_bank().emplace_back();
4,250,141✔
111
  bank.particle = type();
4,250,141✔
112
  bank.wgt = wgt;
4,250,141✔
113
  bank.r = r();
4,250,141✔
114
  bank.u = u();
4,250,141✔
115
  bank.E = settings::run_CE ? E() : g();
4,250,141✔
116
  bank.time = time();
4,250,141✔
117

118
  // Convert signed index to a signed surface ID
119
  if (surface() == SURFACE_NONE) {
4,250,141✔
120
    bank.surf_id = SURFACE_NONE;
4,249,596✔
121
  } else {
122
    int surf_id = model::surfaces[surface_index()]->id_;
545✔
123
    bank.surf_id = (surface() > 0) ? surf_id : -surf_id;
545✔
124
  }
125
}
4,250,141✔
126

127
void Particle::from_source(const SourceSite* src)
240,931,329✔
128
{
129
  // Reset some attributes
130
  clear();
240,931,329✔
131
  surface() = SURFACE_NONE;
240,931,329✔
132
  cell_born() = C_NONE;
240,931,329✔
133
  material() = C_NONE;
240,931,329✔
134
  n_collision() = 0;
240,931,329✔
135
  fission() = false;
240,931,329✔
136
  zero_flux_derivs();
240,931,329✔
137
  lifetime() = 0.0;
240,931,329✔
138
#ifdef OPENMC_DAGMC_ENABLED
139
  history().reset();
22,065,494✔
140
#endif
141

142
  // Copy attributes from source bank site
143
  type() = src->particle;
240,931,329✔
144
  wgt() = src->wgt;
240,931,329✔
145
  wgt_last() = src->wgt;
240,931,329✔
146
  r() = src->r;
240,931,329✔
147
  u() = src->u;
240,931,329✔
148
  r_born() = src->r;
240,931,329✔
149
  r_last_current() = src->r;
240,931,329✔
150
  r_last() = src->r;
240,931,329✔
151
  u_last() = src->u;
240,931,329✔
152
  if (settings::run_CE) {
240,931,329✔
153
    E() = src->E;
125,005,382✔
154
    g() = 0;
125,005,382✔
155
  } else {
156
    g() = static_cast<int>(src->E);
115,925,947✔
157
    g_last() = static_cast<int>(src->E);
115,925,947✔
158
    E() = data::mg.energy_bin_avg_[g()];
115,925,947✔
159
  }
160
  E_last() = E();
240,931,329✔
161
  time() = src->time;
240,931,329✔
162
  time_last() = src->time;
240,931,329✔
163
  parent_nuclide() = src->parent_nuclide;
240,931,329✔
164
  delayed_group() = src->delayed_group;
240,931,329✔
165

166
  // Convert signed surface ID to signed index
167
  if (src->surf_id != SURFACE_NONE) {
240,931,329✔
168
    int index_plus_one = model::surface_map[std::abs(src->surf_id)] + 1;
110,545✔
169
    surface() = (src->surf_id > 0) ? index_plus_one : -index_plus_one;
110,545✔
170
  }
171
}
240,931,329✔
172

173
void Particle::event_calculate_xs()
2,147,483,647✔
174
{
175
  // Set the random number stream
176
  stream() = STREAM_TRACKING;
2,147,483,647✔
177

178
  // Store pre-collision particle properties
179
  wgt_last() = wgt();
2,147,483,647✔
180
  E_last() = E();
2,147,483,647✔
181
  u_last() = u();
2,147,483,647✔
182
  r_last() = r();
2,147,483,647✔
183
  time_last() = time();
2,147,483,647✔
184

185
  // Reset event variables
186
  event() = TallyEvent::KILL;
2,147,483,647✔
187
  event_nuclide() = NUCLIDE_NONE;
2,147,483,647✔
188
  event_mt() = REACTION_NONE;
2,147,483,647✔
189

190
  // If the cell hasn't been determined based on the particle's location,
191
  // initiate a search for the current cell. This generally happens at the
192
  // beginning of the history and again for any secondary particles
193
  if (lowest_coord().cell() == C_NONE) {
2,147,483,647✔
194
    if (!exhaustive_find_cell(*this)) {
232,260,420!
UNCOV
195
      mark_as_lost(
×
UNCOV
196
        "Could not find the cell containing particle " + std::to_string(id()));
×
UNCOV
197
      return;
×
198
    }
199

200
    // Set birth cell attribute
201
    if (cell_born() == C_NONE)
232,260,420!
202
      cell_born() = lowest_coord().cell();
232,260,420✔
203

204
    // Initialize last cells from current cell
205
    for (int j = 0; j < n_coord(); ++j) {
481,992,752✔
206
      cell_last(j) = coord(j).cell();
249,732,332✔
207
    }
208
    n_coord_last() = n_coord();
232,260,420✔
209
  }
210

211
  // Write particle track.
212
  if (write_track())
2,147,483,647✔
213
    write_particle_track(*this);
10,313✔
214

215
  if (settings::check_overlaps)
2,147,483,647!
UNCOV
216
    check_cell_overlap(*this);
×
217

218
  // Calculate microscopic and macroscopic cross sections
219
  if (material() != MATERIAL_VOID) {
2,147,483,647✔
220
    if (settings::run_CE) {
2,147,483,647✔
221
      if (material() != material_last() || sqrtkT() != sqrtkT_last() ||
2,147,483,647✔
222
          density_mult() != density_mult_last()) {
385,350,345✔
223
        // If the material is the same as the last material and the
224
        // temperature hasn't changed, we don't need to lookup cross
225
        // sections again.
226
        model::materials[material()]->calculate_xs(*this);
1,781,909,653✔
227
      }
228
    } else {
229
      // Get the MG data; unlike the CE case above, we have to re-calculate
230
      // cross sections for every collision since the cross sections may
231
      // be angle-dependent
232
      data::mg.macro_xs_[material()].calculate_xs(*this);
2,063,937,744✔
233

234
      // Update the particle's group while we know we are multi-group
235
      g_last() = g();
2,063,937,744✔
236
    }
237
  } else {
238
    macro_xs().total = 0.0;
111,838,608✔
239
    macro_xs().absorption = 0.0;
111,838,608✔
240
    macro_xs().fission = 0.0;
111,838,608✔
241
    macro_xs().nu_fission = 0.0;
111,838,608✔
242
  }
243
}
244

245
void Particle::event_advance()
2,147,483,647✔
246
{
247
  // Find the distance to the nearest boundary
248
  boundary() = distance_to_boundary(*this);
2,147,483,647✔
249

250
  // Sample a distance to collision
251
  if (type() == ParticleType::electron() ||
2,147,483,647✔
252
      type() == ParticleType::positron()) {
2,147,483,647✔
253
    collision_distance() = material() == MATERIAL_VOID ? INFINITY : 0.0;
108,282,934!
254
  } else if (macro_xs().total == 0.0) {
2,147,483,647✔
255
    collision_distance() = INFINITY;
111,838,608✔
256
  } else {
257
    collision_distance() = -std::log(prn(current_seed())) / macro_xs().total;
2,147,483,647✔
258
  }
259

260
  double speed = this->speed();
2,147,483,647✔
261
  double time_cutoff = settings::time_cutoff[type().transport_index()];
2,147,483,647✔
262
  double distance_cutoff =
2,147,483,647✔
263
    (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY;
2,147,483,647✔
264

265
  // Select smaller of the three distances
266
  double distance =
2,147,483,647✔
267
    std::min({boundary().distance(), collision_distance(), distance_cutoff});
2,147,483,647✔
268

269
  // Advance particle in space and time
270
  this->move_distance(distance);
2,147,483,647✔
271
  double dt = distance / speed;
2,147,483,647✔
272
  this->time() += dt;
2,147,483,647✔
273
  this->lifetime() += dt;
2,147,483,647✔
274

275
  // Score timed track-length tallies
276
  if (!model::active_timed_tracklength_tallies.empty()) {
2,147,483,647✔
277
    score_timed_tracklength_tally(*this, distance);
3,628,317✔
278
  }
279

280
  // Score track-length tallies
281
  if (!model::active_tracklength_tallies.empty()) {
2,147,483,647✔
282
    score_tracklength_tally(*this, distance);
1,723,200,332✔
283
  }
284

285
  // Score track-length estimate of k-eff
286
  if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {
2,147,483,647✔
287
    keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission;
2,147,483,647✔
288
  }
289

290
  // Score flux derivative accumulators for differential tallies.
291
  if (!model::active_tallies.empty()) {
2,147,483,647✔
292
    score_track_derivative(*this, distance);
1,893,466,173✔
293
  }
294

295
  // Set particle weight to zero if it hit the time boundary
296
  if (distance == distance_cutoff) {
2,147,483,647✔
297
    wgt() = 0.0;
224,928✔
298
  }
299
}
2,147,483,647✔
300

301
void Particle::event_cross_surface()
2,147,483,647✔
302
{
303
  // Saving previous cell data
304
  for (int j = 0; j < n_coord(); ++j) {
2,147,483,647✔
305
    cell_last(j) = coord(j).cell();
2,147,483,647✔
306
  }
307
  n_coord_last() = n_coord();
2,147,483,647✔
308

309
  // Set surface that particle is on and adjust coordinate levels
310
  surface() = boundary().surface();
2,147,483,647✔
311
  n_coord() = boundary().coord_level();
2,147,483,647✔
312

313
  if (boundary().lattice_translation()[0] != 0 ||
2,147,483,647✔
314
      boundary().lattice_translation()[1] != 0 ||
2,147,483,647✔
315
      boundary().lattice_translation()[2] != 0) {
1,860,089,195✔
316
    // Particle crosses lattice boundary
317

318
    bool verbose = settings::verbosity >= 10 || trace();
750,250,914!
319
    cross_lattice(*this, boundary(), verbose);
750,250,914✔
320
    event() = TallyEvent::LATTICE;
750,250,914✔
321

322
    // Score cell to cell partial currents
323
    if (!model::active_surface_tallies.empty()) {
750,250,914!
UNCOV
324
      auto& lat {*model::lattices[lowest_coord().lattice()]};
×
UNCOV
325
      bool is_valid;
×
UNCOV
326
      Direction normal =
×
UNCOV
327
        lat.get_normal(boundary().lattice_translation(), is_valid);
×
UNCOV
328
      if (is_valid) {
×
UNCOV
329
        normal /= normal.norm();
×
UNCOV
330
        score_surface_tally(*this, model::active_surface_tallies, normal);
×
331
      }
332
    }
333

334
  } else {
335

336
    const auto& surf {*model::surfaces[surface_index()].get()};
1,671,547,301✔
337

338
    // Particle crosses surface
339
    // If BC, add particle to surface source before crossing surface
340
    if (surf.surf_source_ && surf.bc_) {
1,671,547,301✔
341
      add_surf_source_to_bank(*this, surf);
720,304,561✔
342
    }
343
    this->cross_surface(surf);
1,671,547,301✔
344
    // If no BC, add particle to surface source after crossing surface
345
    if (surf.surf_source_ && !surf.bc_) {
1,671,547,292✔
346
      add_surf_source_to_bank(*this, surf);
950,004,904✔
347
    }
348
    if (settings::weight_window_checkpoint_surface) {
1,671,547,292✔
349
      apply_weight_windows(*this);
74,912✔
350
    }
351
    event() = TallyEvent::SURFACE;
1,671,547,292✔
352

353
    // Score cell to cell partial currents
354
    if (!model::active_surface_tallies.empty()) {
1,671,547,292✔
355
      Direction normal = surf.normal(r());
34,931,567✔
356
      normal /= normal.norm();
34,931,567✔
357
      score_surface_tally(*this, model::active_surface_tallies, normal);
34,931,567✔
358
    }
359
  }
360
}
2,147,483,647✔
361

362
void Particle::event_collide()
2,147,483,647✔
363
{
364

365
  // Score collision estimate of keff
366
  if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {
2,147,483,647✔
367
    keff_tally_collision() += wgt() * macro_xs().nu_fission / macro_xs().total;
2,139,390,524✔
368
  }
369

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

374
  if (!model::active_meshsurf_tallies.empty())
2,147,483,647✔
375
    score_meshsurface_tally(*this, model::active_meshsurf_tallies);
63,098,926✔
376

377
  // Clear surface component
378
  surface() = SURFACE_NONE;
2,147,483,647✔
379

380
  if (settings::run_CE) {
2,147,483,647✔
381
    collision(*this);
1,072,606,270✔
382
  } else {
383
    collision_mg(*this);
1,783,060,477✔
384
  }
385

386
  // Collision track feature to recording particle interaction
387
  if (settings::collision_track) {
2,147,483,647✔
388
    collision_track_record(*this);
150,087✔
389
  }
390

391
  // Score collision estimator tallies -- this is done after a collision
392
  // has occurred rather than before because we need information on the
393
  // outgoing energy for any tallies with an outgoing energy filter
394
  if (!model::active_collision_tallies.empty())
2,147,483,647✔
395
    score_collision_tally(*this);
107,843,712✔
396
  if (!model::active_analog_tallies.empty()) {
2,147,483,647✔
397
    if (settings::run_CE) {
406,034,945✔
398
      score_analog_tally_ce(*this);
404,826,683✔
399
    } else {
400
      score_analog_tally_mg(*this);
1,208,262✔
401
    }
402
  }
403

404
  if (!model::active_pulse_height_tallies.empty() && type().is_photon()) {
2,147,483,647✔
405
    pht_collision_energy();
2,024✔
406
  }
407

408
  // Reset banked weight during collision
409
  n_bank() = 0;
2,147,483,647✔
410
  bank_second_E() = 0.0;
2,147,483,647✔
411
  wgt_bank() = 0.0;
2,147,483,647✔
412

413
  // Clear number of secondaries in this collision. This is
414
  // distinct from the number of created neutrons n_bank() above!
415
  n_secondaries() = 0;
2,147,483,647✔
416

417
  zero_delayed_bank();
2,147,483,647✔
418

419
  // Reset fission logical
420
  fission() = false;
2,147,483,647✔
421

422
  // Save coordinates for tallying purposes
423
  r_last_current() = r();
2,147,483,647✔
424

425
  // Set last material to none since cross sections will need to be
426
  // re-evaluated
427
  material_last() = C_NONE;
2,147,483,647✔
428

429
  // Set all directions to base level -- right now, after a collision, only
430
  // the base level directions are changed
431
  for (int j = 0; j < n_coord() - 1; ++j) {
2,147,483,647✔
432
    if (coord(j + 1).rotated()) {
161,247,365✔
433
      // If next level is rotated, apply rotation matrix
434
      const auto& m {model::cells[coord(j).cell()]->rotation_};
10,426,614✔
435
      const auto& u {coord(j).u()};
10,426,614✔
436
      coord(j + 1).u() = u.rotate(m);
10,426,614✔
437
    } else {
438
      // Otherwise, copy this level's direction
439
      coord(j + 1).u() = coord(j).u();
150,820,751✔
440
    }
441
  }
442

443
  // Score flux derivative accumulators for differential tallies.
444
  if (!model::active_tallies.empty())
2,147,483,647✔
445
    score_collision_derivative(*this);
937,993,381✔
446

447
#ifdef OPENMC_DAGMC_ENABLED
448
  history().reset();
261,721,037✔
449
#endif
450
}
2,147,483,647✔
451

452
void Particle::event_revive_from_secondary()
2,147,483,647✔
453
{
454
  // If particle has too many events, display warning and kill it
455
  ++n_event();
2,147,483,647✔
456
  if (n_event() == settings::max_particle_events) {
2,147,483,647!
UNCOV
457
    warning("Particle " + std::to_string(id()) +
×
458
            " underwent maximum number of events.");
UNCOV
459
    wgt() = 0.0;
×
460
  }
461

462
  // Check for secondary particles if this particle is dead
463
  if (!alive()) {
2,147,483,647✔
464
    // Write final position for this particle
465
    if (write_track()) {
232,260,016✔
466
      write_particle_track(*this);
6,244✔
467
    }
468

469
    // If no secondary particles, break out of event loop
470
    if (secondary_bank().empty())
232,260,016✔
471
      return;
472

473
    from_source(&secondary_bank().back());
63,347,715✔
474
    secondary_bank().pop_back();
63,347,715✔
475
    n_event() = 0;
63,347,715✔
476
    bank_second_E() = 0.0;
63,347,715✔
477

478
    // Subtract secondary particle energy from interim pulse-height results
479
    if (!model::active_pulse_height_tallies.empty() &&
63,347,715✔
480
        this->type().is_photon()) {
15,499✔
481
      // Since the birth cell of the particle has not been set we
482
      // have to determine it before the energy of the secondary particle can be
483
      // removed from the pulse-height of this cell.
484
      if (lowest_coord().cell() == C_NONE) {
605!
485
        bool verbose = settings::verbosity >= 10 || trace();
605!
486
        if (!exhaustive_find_cell(*this, verbose)) {
605!
UNCOV
487
          mark_as_lost("Could not find the cell containing particle " +
×
UNCOV
488
                       std::to_string(id()));
×
UNCOV
489
          return;
×
490
        }
491
        // Set birth cell attribute
492
        if (cell_born() == C_NONE)
605!
493
          cell_born() = lowest_coord().cell();
605✔
494

495
        // Initialize last cells from current cell
496
        for (int j = 0; j < n_coord(); ++j) {
1,210✔
497
          cell_last(j) = coord(j).cell();
605✔
498
        }
499
        n_coord_last() = n_coord();
605✔
500
      }
501
      pht_secondary_particles();
605✔
502
    }
503

504
    // Enter new particle in particle track file
505
    if (write_track())
63,347,715✔
506
      add_particle_track(*this);
5,234✔
507
  }
508
}
509

510
void Particle::event_death()
168,913,301✔
511
{
512
#ifdef OPENMC_DAGMC_ENABLED
513
  history().reset();
15,432,415✔
514
#endif
515

516
  // Finish particle track output.
517
  if (write_track()) {
168,913,301✔
518
    finalize_particle_track(*this);
1,010✔
519
  }
520

521
// Contribute tally reduction variables to global accumulator
522
#pragma omp atomic
92,994,579✔
523
  global_tally_absorption += keff_tally_absorption();
168,913,301✔
524
#pragma omp atomic
93,305,087✔
525
  global_tally_collision += keff_tally_collision();
168,913,301✔
526
#pragma omp atomic
92,869,844✔
527
  global_tally_tracklength += keff_tally_tracklength();
168,913,301✔
528
#pragma omp atomic
92,427,053✔
529
  global_tally_leakage += keff_tally_leakage();
168,913,301✔
530

531
  // Reset particle tallies once accumulated
532
  keff_tally_absorption() = 0.0;
168,913,301✔
533
  keff_tally_collision() = 0.0;
168,913,301✔
534
  keff_tally_tracklength() = 0.0;
168,913,301✔
535
  keff_tally_leakage() = 0.0;
168,913,301✔
536

537
  if (!model::active_pulse_height_tallies.empty()) {
168,913,301✔
538
    score_pulse_height_tally(*this, model::active_pulse_height_tallies);
5,500✔
539
  }
540

541
  // Record the number of progeny created by this particle.
542
  // This data will be used to efficiently sort the fission bank.
543
  if (settings::run_mode == RunMode::EIGENVALUE) {
168,913,301✔
544
    int64_t offset = id() - 1 - simulation::work_index[mpi::rank];
141,534,700✔
545
    simulation::progeny_per_particle[offset] = n_progeny();
141,534,700✔
546
  }
547
}
168,913,301✔
548

549
void Particle::pht_collision_energy()
2,024✔
550
{
551
  // Adds the energy particles lose in a collision to the pulse-height
552

553
  // determine index of cell in pulse_height_cells
554
  auto it = std::find(model::pulse_height_cells.begin(),
2,024✔
555
    model::pulse_height_cells.end(), lowest_coord().cell());
2,024!
556

557
  if (it != model::pulse_height_cells.end()) {
2,024!
558
    int index = std::distance(model::pulse_height_cells.begin(), it);
2,024✔
559
    pht_storage()[index] += E_last() - E();
2,024✔
560

561
    // If the energy of the particle is below the cutoff, it will not be sampled
562
    // so its energy is added to the pulse-height in the cell
563
    int photon = ParticleType::photon().transport_index();
2,024✔
564
    if (E() < settings::energy_cutoff[photon]) {
2,024✔
565
      pht_storage()[index] += E();
825✔
566
    }
567
  }
568
}
2,024✔
569

570
void Particle::pht_secondary_particles()
605✔
571
{
572
  // Removes the energy of secondary produced particles from the pulse-height
573

574
  // determine index of cell in pulse_height_cells
575
  auto it = std::find(model::pulse_height_cells.begin(),
605✔
576
    model::pulse_height_cells.end(), cell_born());
605!
577

578
  if (it != model::pulse_height_cells.end()) {
605!
579
    int index = std::distance(model::pulse_height_cells.begin(), it);
605✔
580
    pht_storage()[index] -= E();
605✔
581
  }
582
}
605✔
583

584
void Particle::cross_surface(const Surface& surf)
1,673,434,599✔
585
{
586

587
  if (settings::verbosity >= 10 || trace()) {
1,673,434,599✔
588
    write_message(1, "    Crossing surface {}", surf.id_);
66✔
589
  }
590

591
// if we're crossing a CSG surface, make sure the DAG history is reset
592
#ifdef OPENMC_DAGMC_ENABLED
593
  if (surf.geom_type() == GeometryType::CSG)
152,810,101✔
594
    history().reset();
152,752,942✔
595
#endif
596

597
  // Handle any applicable boundary conditions.
598
  if (surf.bc_ && settings::run_mode != RunMode::PLOTTING &&
1,673,434,599!
599
      settings::run_mode != RunMode::VOLUME) {
600
    surf.bc_->handle_particle(*this, surf);
720,656,669✔
601
    return;
720,656,669✔
602
  }
603

604
  // ==========================================================================
605
  // SEARCH NEIGHBOR LISTS FOR NEXT CELL
606

607
#ifdef OPENMC_DAGMC_ENABLED
608
  // in DAGMC, we know what the next cell should be
609
  if (surf.geom_type() == GeometryType::DAG) {
86,853,401✔
610
    int32_t i_cell = next_cell(surface_index(), cell_last(n_coord() - 1),
46,350✔
611
                       lowest_coord().universe()) -
46,350✔
612
                     1;
46,350✔
613
    // save material, temperature, and density multiplier
614
    material_last() = material();
46,350✔
615
    sqrtkT_last() = sqrtkT();
46,350✔
616
    density_mult_last() = density_mult();
46,350✔
617
    // set new cell value
618
    lowest_coord().cell() = i_cell;
46,350✔
619
    auto& cell = model::cells[i_cell];
46,350✔
620

621
    cell_instance() = 0;
46,350✔
622
    if (cell->distribcell_index_ >= 0)
46,350✔
623
      cell_instance() = cell_instance_at_level(*this, n_coord() - 1);
45,326✔
624

625
    material() = cell->material(cell_instance());
46,350!
626
    sqrtkT() = cell->sqrtkT(cell_instance());
46,350!
627
    density_mult() = cell->density_mult(cell_instance());
46,350✔
628
    return;
46,350✔
629
  }
630
#endif
631

632
  bool verbose = settings::verbosity >= 10 || trace();
952,731,580!
633
  if (neighbor_list_find_cell(*this, verbose)) {
952,731,580✔
634
    return;
635
  }
636

637
  // ==========================================================================
638
  // COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS
639

640
  // Remove lower coordinate levels
641
  n_coord() = 1;
29,911✔
642
  bool found = exhaustive_find_cell(*this, verbose);
29,911✔
643

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

650
    surface() = SURFACE_NONE;
5,799✔
651
    n_coord() = 1;
5,799✔
652
    r() += TINY_BIT * u();
5,799✔
653

654
    // Couldn't find next cell anywhere! This probably means there is an actual
655
    // undefined region in the geometry.
656

657
    if (!exhaustive_find_cell(*this, verbose)) {
5,799!
658
      mark_as_lost("After particle " + std::to_string(id()) +
17,388✔
659
                   " crossed surface " + std::to_string(surf.id_) +
17,388✔
660
                   " it could not be located in any cell and it did not leak.");
661
      return;
5,790✔
662
    }
663
  }
664
}
665

666
void Particle::cross_vacuum_bc(const Surface& surf)
35,093,373✔
667
{
668
  // Score any surface current tallies -- note that the particle is moved
669
  // forward slightly so that if the mesh boundary is on the surface, it is
670
  // still processed
671

672
  if (!model::active_meshsurf_tallies.empty()) {
35,093,373✔
673
    // TODO: Find a better solution to score surface currents than
674
    // physically moving the particle forward slightly
675

676
    r() += TINY_BIT * u();
937,222✔
677
    score_meshsurface_tally(*this, model::active_meshsurf_tallies);
937,222✔
678
  }
679

680
  // Score to global leakage tally
681
  keff_tally_leakage() += wgt();
35,093,373✔
682

683
  // Kill the particle
684
  wgt() = 0.0;
35,093,373✔
685

686
  // Display message
687
  if (settings::verbosity >= 10 || trace()) {
35,093,373!
688
    write_message(1, "    Leaked out of surface {}", surf.id_);
22✔
689
  }
690
}
35,093,373✔
691

692
void Particle::cross_reflective_bc(const Surface& surf, Direction new_u)
684,318,903✔
693
{
694
  // Do not handle reflective boundary conditions on lower universes
695
  if (n_coord() != 1) {
684,318,903!
UNCOV
696
    mark_as_lost("Cannot reflect particle " + std::to_string(id()) +
×
697
                 " off surface in a lower universe.");
UNCOV
698
    return;
×
699
  }
700

701
  // Score surface currents since reflection causes the direction of the
702
  // particle to change. For surface filters, we need to score the tallies
703
  // twice, once before the particle's surface attribute has changed and
704
  // once after. For mesh surface filters, we need to artificially move
705
  // the particle slightly back in case the surface crossing is coincident
706
  // with a mesh boundary
707

708
  if (!model::active_surface_tallies.empty()) {
684,318,903✔
709
    Direction normal = surf.normal(r());
285,021✔
710
    normal /= normal.norm();
285,021✔
711
    score_surface_tally(*this, model::active_surface_tallies, normal);
285,021✔
712
  }
713

714
  if (!model::active_meshsurf_tallies.empty()) {
684,318,903✔
715
    Position r {this->r()};
46,885,487✔
716
    this->r() -= TINY_BIT * u();
46,885,487✔
717
    score_meshsurface_tally(*this, model::active_meshsurf_tallies);
46,885,487✔
718
    this->r() = r;
46,885,487✔
719
  }
720

721
  // Set the new particle direction
722
  u() = new_u;
684,318,903✔
723

724
  // Reassign particle's cell and surface
725
  coord(0).cell() = cell_last(0);
684,318,903✔
726
  surface() = -surface();
684,318,903✔
727

728
  // If a reflective surface is coincident with a lattice or universe
729
  // boundary, it is necessary to redetermine the particle's coordinates in
730
  // the lower universes.
731
  // (unless we're using a dagmc model, which has exactly one universe)
732
  n_coord() = 1;
684,318,903✔
733
  if (surf.geom_type() != GeometryType::DAG &&
1,368,635,048!
734
      !neighbor_list_find_cell(*this)) {
684,316,145✔
UNCOV
735
    mark_as_lost("Couldn't find particle after reflecting from surface " +
×
UNCOV
736
                 std::to_string(surf.id_) + ".");
×
UNCOV
737
    return;
×
738
  }
739

740
  // Set previous coordinate going slightly past surface crossing
741
  r_last_current() = r() + TINY_BIT * u();
684,318,903✔
742

743
  // Diagnostic message
744
  if (settings::verbosity >= 10 || trace()) {
684,318,903!
UNCOV
745
    write_message(1, "    Reflected from surface {}", surf.id_);
×
746
  }
747
}
748

749
void Particle::cross_periodic_bc(
2,249,859✔
750
  const Surface& surf, Position new_r, Direction new_u, int new_surface)
751
{
752
  // Do not handle periodic boundary conditions on lower universes
753
  if (n_coord() != 1) {
2,249,859!
UNCOV
754
    mark_as_lost(
×
UNCOV
755
      "Cannot transfer particle " + std::to_string(id()) +
×
756
      " across surface in a lower universe. Boundary conditions must be "
757
      "applied to root universe.");
758
    return;
×
759
  }
760

761
  // Score surface currents since reflection causes the direction of the
762
  // particle to change -- artificially move the particle slightly back in
763
  // case the surface crossing is coincident with a mesh boundary
764
  if (!model::active_meshsurf_tallies.empty()) {
2,249,859!
UNCOV
765
    Position r {this->r()};
×
UNCOV
766
    this->r() -= TINY_BIT * u();
×
UNCOV
767
    score_meshsurface_tally(*this, model::active_meshsurf_tallies);
×
UNCOV
768
    this->r() = r;
×
769
  }
770

771
  // Adjust the particle's location and direction.
772
  r() = new_r;
2,249,859✔
773
  u() = new_u;
2,249,859✔
774

775
  // Reassign particle's surface
776
  surface() = new_surface;
2,249,859✔
777

778
  // Figure out what cell particle is in now
779
  n_coord() = 1;
2,249,859✔
780

781
  if (!neighbor_list_find_cell(*this)) {
2,249,859!
UNCOV
782
    mark_as_lost("Couldn't find particle after hitting periodic "
×
UNCOV
783
                 "boundary on surface " +
×
UNCOV
784
                 std::to_string(surf.id_) + ".");
×
UNCOV
785
    return;
×
786
  }
787

788
  // Set previous coordinate going slightly past surface crossing
789
  r_last_current() = r() + TINY_BIT * u();
2,249,859✔
790

791
  // Diagnostic message
792
  if (settings::verbosity >= 10 || trace()) {
2,249,859!
UNCOV
793
    write_message(1, "    Hit periodic boundary on surface {}", surf.id_);
×
794
  }
795
}
796

797
void Particle::mark_as_lost(const char* message)
5,799✔
798
{
799
  // Print warning and write lost particle file
800
  warning(message);
5,799✔
801
  if (settings::max_write_lost_particles < 0 ||
5,799✔
802
      simulation::n_lost_particles < settings::max_write_lost_particles) {
5,500✔
803
    write_restart();
374✔
804
  }
805
  // Increment number of lost particles
806
  wgt() = 0.0;
5,799✔
807
#pragma omp atomic
3,154✔
808
  simulation::n_lost_particles += 1;
2,645✔
809

810
  // Count the total number of simulated particles (on this processor)
811
  auto n = simulation::current_batch * settings::gen_per_batch *
5,799✔
812
           simulation::work_per_rank;
813

814
  // Abort the simulation if the maximum number of lost particles has been
815
  // reached
816
  if (simulation::n_lost_particles >= settings::max_lost_particles &&
5,799✔
817
      simulation::n_lost_particles >= settings::rel_max_lost_particles * n) {
9!
818
    fatal_error("Maximum number of lost particles has been reached.");
9✔
819
  }
820
}
5,790✔
821

822
void Particle::write_restart() const
374✔
823
{
824
  // Dont write another restart file if in particle restart mode
825
  if (settings::run_mode == RunMode::PARTICLE)
374✔
826
    return;
22✔
827

828
  // Set up file name
829
  auto filename = fmt::format("{}particle_{}_{}.h5", settings::path_output,
352✔
830
    simulation::current_batch, id());
352✔
831

832
#pragma omp critical(WriteParticleRestart)
187✔
833
  {
352✔
834
    // Create file
835
    hid_t file_id = file_open(filename, 'w');
352✔
836

837
    // Write filetype and version info
838
    write_attribute(file_id, "filetype", "particle restart");
352✔
839
    write_attribute(file_id, "version", VERSION_PARTICLE_RESTART);
352✔
840
    write_attribute(file_id, "openmc_version", VERSION);
352✔
841
#ifdef GIT_SHA1
842
    write_attr_string(file_id, "git_sha1", GIT_SHA1);
843
#endif
844

845
    // Write data to file
846
    write_dataset(file_id, "current_batch", simulation::current_batch);
352✔
847
    write_dataset(file_id, "generations_per_batch", settings::gen_per_batch);
352✔
848
    write_dataset(file_id, "current_generation", simulation::current_gen);
352✔
849
    write_dataset(file_id, "n_particles", settings::n_particles);
352✔
850
    switch (settings::run_mode) {
352!
851
    case RunMode::FIXED_SOURCE:
220✔
852
      write_dataset(file_id, "run_mode", "fixed source");
220✔
853
      break;
115✔
854
    case RunMode::EIGENVALUE:
132✔
855
      write_dataset(file_id, "run_mode", "eigenvalue");
132✔
856
      break;
72✔
UNCOV
857
    case RunMode::PARTICLE:
×
UNCOV
858
      write_dataset(file_id, "run_mode", "particle restart");
×
859
      break;
860
    default:
861
      break;
862
    }
863
    write_dataset(file_id, "id", id());
352✔
864
    write_dataset(file_id, "type", type().pdg_number());
352✔
865

866
    int64_t i = current_work();
352✔
867
    if (settings::run_mode == RunMode::EIGENVALUE) {
352✔
868
      // take source data from primary bank for eigenvalue simulation
869
      write_dataset(file_id, "weight", simulation::source_bank[i - 1].wgt);
132✔
870
      write_dataset(file_id, "energy", simulation::source_bank[i - 1].E);
132✔
871
      write_dataset(file_id, "xyz", simulation::source_bank[i - 1].r);
132✔
872
      write_dataset(file_id, "uvw", simulation::source_bank[i - 1].u);
132✔
873
      write_dataset(file_id, "time", simulation::source_bank[i - 1].time);
132✔
874
    } else if (settings::run_mode == RunMode::FIXED_SOURCE) {
220!
875
      // re-sample using rng random number seed used to generate source particle
876
      int64_t id = (simulation::total_gen + overall_generation() - 1) *
220✔
877
                     settings::n_particles +
220✔
878
                   simulation::work_index[mpi::rank] + i;
220✔
879
      uint64_t seed = init_seed(id, STREAM_SOURCE);
220✔
880
      // re-sample source site
881
      auto site = sample_external_source(&seed);
220✔
882
      write_dataset(file_id, "weight", site.wgt);
220✔
883
      write_dataset(file_id, "energy", site.E);
220✔
884
      write_dataset(file_id, "xyz", site.r);
220✔
885
      write_dataset(file_id, "uvw", site.u);
220✔
886
      write_dataset(file_id, "time", site.time);
220✔
887
    }
888

889
    // Close file
890
    file_close(file_id);
352✔
891
  } // #pragma omp critical
892
}
352✔
893

894
void Particle::update_neutron_xs(
2,147,483,647✔
895
  int i_nuclide, int i_grid, int i_sab, double sab_frac, double ncrystal_xs)
896
{
897
  // Get microscopic cross section cache
898
  auto& micro = this->neutron_xs(i_nuclide);
2,147,483,647✔
899

900
  // If the cache doesn't match, recalculate micro xs
901
  if (this->E() != micro.last_E || this->sqrtkT() != micro.last_sqrtkT ||
2,147,483,647✔
902
      i_sab != micro.index_sab || sab_frac != micro.sab_frac ||
2,147,483,647✔
903
      ncrystal_xs != micro.ncrystal_xs) {
2,147,483,647!
904
    data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, *this);
2,147,483,647✔
905

906
    // If NCrystal is being used, update micro cross section cache
907
    micro.ncrystal_xs = ncrystal_xs;
2,147,483,647✔
908
    if (ncrystal_xs >= 0.0) {
2,147,483,647✔
909
      data::nuclides[i_nuclide]->calculate_elastic_xs(*this);
11,018,953✔
910
      ncrystal_update_micro(ncrystal_xs, micro);
11,018,953✔
911
    }
912
  }
913
}
2,147,483,647✔
914

915
//==============================================================================
916
// Non-method functions
917
//==============================================================================
918
void add_surf_source_to_bank(Particle& p, const Surface& surf)
1,670,309,465✔
919
{
920
  if (simulation::current_batch <= settings::n_inactive ||
1,670,309,465✔
921
      simulation::surf_source_bank.full()) {
1,284,837,127✔
922
    return;
1,670,179,812✔
923
  }
924

925
  // If a cell/cellfrom/cellto parameter is defined
926
  if (settings::ssw_cell_id != C_NONE) {
337,080✔
927

928
    // Retrieve cell index and storage type
929
    int cell_idx = model::cell_map[settings::ssw_cell_id];
254,435✔
930

931
    if (surf.bc_) {
254,435✔
932
      // Leave if cellto with vacuum boundary condition
933
      if (surf.bc_->type() == "vacuum" &&
298,918✔
934
          settings::ssw_cell_type == SSWCellType::To) {
33,099✔
935
        return;
936
      }
937

938
      // Leave if other boundary condition than vacuum
939
      if (surf.bc_->type() != "vacuum") {
274,648✔
940
        return;
941
      }
942
    }
943

944
    // Check if the cell of interest has been exited
945
    bool exited = false;
946
    for (int i = 0; i < p.n_coord_last(); ++i) {
333,669✔
947
      if (p.cell_last(i) == cell_idx) {
207,729✔
948
        exited = true;
73,763✔
949
      }
950
    }
951

952
    // Check if the cell of interest has been entered
953
    bool entered = false;
954
    for (int i = 0; i < p.n_coord(); ++i) {
297,971✔
955
      if (p.coord(i).cell() == cell_idx) {
172,031✔
956
        entered = true;
57,517✔
957
      }
958
    }
959

960
    // Vacuum boundary conditions: return if cell is not exited
961
    if (surf.bc_) {
125,940✔
962
      if (surf.bc_->type() == "vacuum" && !exited) {
41,928!
963
        return;
964
      }
965
    } else {
966

967
      // If we both enter and exit the cell of interest
968
      if (entered && exited) {
104,976✔
969
        return;
970
      }
971

972
      // If we did not enter nor exit the cell of interest
973
      if (!entered && !exited) {
77,773✔
974
        return;
975
      }
976

977
      // If cellfrom and the cell before crossing is not the cell of
978
      // interest
979
      if (settings::ssw_cell_type == SSWCellType::From && !exited) {
64,274✔
980
        return;
981
      }
982

983
      // If cellto and the cell after crossing is not the cell of interest
984
      if (settings::ssw_cell_type == SSWCellType::To && !entered) {
52,731✔
985
        return;
986
      }
987
    }
988
  }
989

990
  SourceSite site;
129,653✔
991
  site.r = p.r();
129,653✔
992
  site.u = p.u();
129,653✔
993
  site.E = p.E();
129,653✔
994
  site.time = p.time();
129,653✔
995
  site.wgt = p.wgt();
129,653✔
996
  site.delayed_group = p.delayed_group();
129,653✔
997
  site.surf_id = surf.id_;
129,653✔
998
  site.particle = p.type();
129,653✔
999
  site.parent_id = p.id();
129,653✔
1000
  site.progeny_id = p.n_progeny();
129,653✔
1001
  int64_t idx = simulation::surf_source_bank.thread_safe_append(site);
129,653✔
1002
}
1003

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