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

openmc-dev / openmc / 22669057039

04 Mar 2026 12:19PM UTC coverage: 81.555% (-0.003%) from 81.558%
22669057039

Pull #3849

github

web-flow
Merge 6ca13f81e into 70be65000
Pull Request #3849: Fixed a bug in cell last update

17537 of 25247 branches covered (69.46%)

Branch coverage included in aggregate %.

8 of 8 new or added lines in 1 file covered. (100.0%)

26 existing lines in 2 files now uncovered.

57926 of 67283 relevant lines covered (86.09%)

44918489.62 hits per line

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

87.1
/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/material.h"
18
#include "openmc/message_passing.h"
19
#include "openmc/mgxs_interface.h"
20
#include "openmc/nuclide.h"
21
#include "openmc/particle_data.h"
22
#include "openmc/photon.h"
23
#include "openmc/physics.h"
24
#include "openmc/physics_mg.h"
25
#include "openmc/random_lcg.h"
26
#include "openmc/settings.h"
27
#include "openmc/simulation.h"
28
#include "openmc/source.h"
29
#include "openmc/surface.h"
30
#include "openmc/tallies/derivative.h"
31
#include "openmc/tallies/tally.h"
32
#include "openmc/tallies/tally_scoring.h"
33
#include "openmc/track_output.h"
34
#include "openmc/weight_windows.h"
35

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

40
namespace openmc {
41

42
//==============================================================================
43
// Particle implementation
44
//==============================================================================
45

46
double Particle::speed() const
2,147,483,647✔
47
{
48
  if (settings::run_CE) {
2,147,483,647✔
49
    // Determine mass in eV/c^2
50
    double mass;
2,147,483,647✔
51
    switch (type().pdg_number()) {
2,147,483,647✔
52
    case PDG_NEUTRON:
53
      mass = MASS_NEUTRON_EV;
54
    case PDG_ELECTRON:
55
    case PDG_POSITRON:
56
      mass = MASS_ELECTRON_EV;
2,147,483,647✔
57
    default:
2,147,483,647✔
58
      mass = this->type().mass() * AMU_EV;
2,147,483,647✔
59
    }
60

61
    // Equivalent to C * sqrt(1-(m/(m+E))^2) without problem at E<<m:
62
    return C_LIGHT * std::sqrt(this->E() * (this->E() + 2 * mass)) /
2,147,483,647✔
63
           (this->E() + mass);
2,147,483,647✔
64
  } else {
65
    auto& macro_xs = data::mg.macro_xs_[this->material()];
2,063,937,414✔
66
    int macro_t = this->mg_xs_cache().t;
2,063,937,414✔
67
    int macro_a = macro_xs.get_angle_index(this->u());
2,063,937,414✔
68
    return 1.0 / macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, this->g(), nullptr,
2,063,937,414✔
69
                   nullptr, nullptr, macro_t, macro_a);
2,063,937,414✔
70
  }
71
}
72

73
bool Particle::create_secondary(
112,449,300✔
74
  double wgt, Direction u, double E, ParticleType type)
75
{
76
  // If energy is below cutoff for this particle, don't create secondary
77
  // particle
78
  int idx = type.transport_index();
112,449,300✔
79
  if (idx == C_NONE) {
112,449,300!
80
    return false;
81
  }
82
  if (E < settings::energy_cutoff[idx]) {
112,449,300✔
83
    return false;
84
  }
85

86
  // Increment number of secondaries created (for ParticleProductionFilter)
87
  n_secondaries()++;
58,678,232✔
88

89
  auto& bank = secondary_bank().emplace_back();
58,678,232✔
90
  bank.particle = type;
58,678,232✔
91
  bank.wgt = wgt;
58,678,232✔
92
  bank.r = r();
58,678,232!
93
  bank.u = u;
58,678,232✔
94
  bank.E = settings::run_CE ? E : g();
58,678,232!
95
  bank.time = time();
58,678,232✔
96
  bank_second_E() += bank.E;
58,678,232✔
97
  return true;
58,678,232✔
98
}
99

100
void Particle::split(double wgt)
4,250,186✔
101
{
102
  auto& bank = secondary_bank().emplace_back();
4,250,186✔
103
  bank.particle = type();
4,250,186✔
104
  bank.wgt = wgt;
4,250,186✔
105
  bank.r = r();
4,250,186✔
106
  bank.u = u();
4,250,186✔
107
  bank.E = settings::run_CE ? E() : g();
4,250,186✔
108
  bank.time = time();
4,250,186✔
109

110
  // Convert signed index to a signed surface ID
111
  if (surface() == SURFACE_NONE) {
4,250,186✔
112
    bank.surf_id = SURFACE_NONE;
4,249,641✔
113
  } else {
114
    int surf_id = model::surfaces[surface_index()]->id_;
545✔
115
    bank.surf_id = (surface() > 0) ? surf_id : -surf_id;
545✔
116
  }
117
}
4,250,186✔
118

119
void Particle::from_source(const SourceSite* src)
240,605,221✔
120
{
121
  // Reset some attributes
122
  clear();
240,605,221✔
123
  surface() = SURFACE_NONE;
240,605,221✔
124
  cell_born() = C_NONE;
240,605,221✔
125
  material() = C_NONE;
240,605,221✔
126
  n_collision() = 0;
240,605,221✔
127
  fission() = false;
240,605,221✔
128
  zero_flux_derivs();
240,605,221✔
129
  lifetime() = 0.0;
240,605,221✔
130
#ifdef OPENMC_DAGMC_ENABLED
131
  history().reset();
22,035,388✔
132
#endif
133

134
  // Copy attributes from source bank site
135
  type() = src->particle;
240,605,221✔
136
  wgt() = src->wgt;
240,605,221✔
137
  wgt_last() = src->wgt;
240,605,221✔
138
  r() = src->r;
240,605,221✔
139
  u() = src->u;
240,605,221✔
140
  r_born() = src->r;
240,605,221✔
141
  r_last_current() = src->r;
240,605,221✔
142
  r_last() = src->r;
240,605,221✔
143
  u_last() = src->u;
240,605,221✔
144
  if (settings::run_CE) {
240,605,221✔
145
    E() = src->E;
124,679,604✔
146
    g() = 0;
124,679,604✔
147
  } else {
148
    g() = static_cast<int>(src->E);
115,925,617✔
149
    g_last() = static_cast<int>(src->E);
115,925,617✔
150
    E() = data::mg.energy_bin_avg_[g()];
115,925,617✔
151
  }
152
  E_last() = E();
240,605,221✔
153
  time() = src->time;
240,605,221✔
154
  time_last() = src->time;
240,605,221✔
155
  parent_nuclide() = src->parent_nuclide;
240,605,221✔
156
  delayed_group() = src->delayed_group;
240,605,221✔
157

158
  // Convert signed surface ID to signed index
159
  if (src->surf_id != SURFACE_NONE) {
240,605,221✔
160
    int index_plus_one = model::surface_map[std::abs(src->surf_id)] + 1;
110,545✔
161
    surface() = (src->surf_id > 0) ? index_plus_one : -index_plus_one;
110,545✔
162
  }
163
}
240,605,221✔
164

165
void Particle::event_calculate_xs()
2,147,483,647✔
166
{
167
  // Set the random number stream
168
  stream() = STREAM_TRACKING;
2,147,483,647✔
169

170
  // Store pre-collision particle properties
171
  wgt_last() = wgt();
2,147,483,647✔
172
  E_last() = E();
2,147,483,647✔
173
  u_last() = u();
2,147,483,647✔
174
  r_last() = r();
2,147,483,647✔
175
  time_last() = time();
2,147,483,647✔
176

177
  // If the cell hasn't been determined based on the particle's location,
178
  // initiate a search for the current cell. This generally happens at the
179
  // beginning of the history and again for any secondary particles
180
  if (lowest_coord().cell() == C_NONE) {
2,147,483,647✔
181
    if (!exhaustive_find_cell(*this)) {
231,934,312!
182
      mark_as_lost(
×
183
        "Could not find the cell containing particle " + std::to_string(id()));
×
184
      return;
×
185
    }
186

187
    // Set birth cell attribute
188
    if (cell_born() == C_NONE)
231,934,312!
189
      cell_born() = lowest_coord().cell();
231,934,312✔
190
  }
191

192
  // Initialize last cells from current cell
193
  for (int j = 0; j < n_coord(); ++j) {
2,147,483,647✔
194
    cell_last(j) = coord(j).cell();
2,147,483,647✔
195
  }
196
  n_coord_last() = n_coord();
2,147,483,647✔
197

198
  // Reset event variables
199
  event() = TallyEvent::KILL;
2,147,483,647✔
200
  event_nuclide() = NUCLIDE_NONE;
2,147,483,647✔
201
  event_mt() = REACTION_NONE;
2,147,483,647✔
202

203
  // Write particle track.
204
  if (write_track())
2,147,483,647✔
205
    write_particle_track(*this);
10,309✔
206

207
  if (settings::check_overlaps)
2,147,483,647!
208
    check_cell_overlap(*this);
×
209

210
  // Calculate microscopic and macroscopic cross sections
211
  if (material() != MATERIAL_VOID) {
2,147,483,647✔
212
    if (settings::run_CE) {
2,147,483,647✔
213
      if (material() != material_last() || sqrtkT() != sqrtkT_last() ||
2,147,483,647✔
214
          density_mult() != density_mult_last()) {
385,346,733✔
215
        // If the material is the same as the last material and the
216
        // temperature hasn't changed, we don't need to lookup cross
217
        // sections again.
218
        model::materials[material()]->calculate_xs(*this);
1,781,410,774✔
219
      }
220
    } else {
221
      // Get the MG data; unlike the CE case above, we have to re-calculate
222
      // cross sections for every collision since the cross sections may
223
      // be angle-dependent
224
      data::mg.macro_xs_[material()].calculate_xs(*this);
2,063,937,414✔
225

226
      // Update the particle's group while we know we are multi-group
227
      g_last() = g();
2,063,937,414✔
228
    }
229
  } else {
230
    macro_xs().total = 0.0;
111,838,554✔
231
    macro_xs().absorption = 0.0;
111,838,554✔
232
    macro_xs().fission = 0.0;
111,838,554✔
233
    macro_xs().nu_fission = 0.0;
111,838,554✔
234
  }
235
}
236

237
void Particle::event_advance()
2,147,483,647✔
238
{
239
  // Find the distance to the nearest boundary
240
  boundary() = distance_to_boundary(*this);
2,147,483,647✔
241

242
  // Sample a distance to collision
243
  if (type() == ParticleType::electron() ||
2,147,483,647✔
244
      type() == ParticleType::positron()) {
2,147,483,647✔
245
    collision_distance() = material() == MATERIAL_VOID ? INFINITY : 0.0;
107,847,232!
246
  } else if (macro_xs().total == 0.0) {
2,147,483,647✔
247
    collision_distance() = INFINITY;
111,838,554✔
248
  } else {
249
    collision_distance() = -std::log(prn(current_seed())) / macro_xs().total;
2,147,483,647✔
250
  }
251

252
  double speed = this->speed();
2,147,483,647✔
253
  double time_cutoff = settings::time_cutoff[type().transport_index()];
2,147,483,647✔
254
  double distance_cutoff =
2,147,483,647✔
255
    (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY;
2,147,483,647✔
256

257
  // Select smaller of the three distances
258
  double distance =
2,147,483,647✔
259
    std::min({boundary().distance(), collision_distance(), distance_cutoff});
2,147,483,647✔
260

261
  // Advance particle in space and time
262
  this->move_distance(distance);
2,147,483,647✔
263
  double dt = distance / speed;
2,147,483,647✔
264
  this->time() += dt;
2,147,483,647✔
265
  this->lifetime() += dt;
2,147,483,647✔
266

267
  // Score timed track-length tallies
268
  if (!model::active_timed_tracklength_tallies.empty()) {
2,147,483,647✔
269
    score_timed_tracklength_tally(*this, distance);
3,628,317✔
270
  }
271

272
  // Score track-length tallies
273
  if (!model::active_tracklength_tallies.empty()) {
2,147,483,647✔
274
    score_tracklength_tally(*this, distance);
1,723,195,749✔
275
  }
276

277
  // Score track-length estimate of k-eff
278
  if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {
2,147,483,647✔
279
    keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission;
2,147,483,647✔
280
  }
281

282
  // Score flux derivative accumulators for differential tallies.
283
  if (!model::active_tallies.empty()) {
2,147,483,647✔
284
    score_track_derivative(*this, distance);
1,892,966,601✔
285
  }
286

287
  // Set particle weight to zero if it hit the time boundary
288
  if (distance == distance_cutoff) {
2,147,483,647✔
289
    wgt() = 0.0;
224,928✔
290
  }
291
}
2,147,483,647✔
292

293
void Particle::event_cross_surface()
2,147,483,647✔
294
{
295
  // Saving previous cell data
296
  for (int j = 0; j < n_coord(); ++j) {
2,147,483,647✔
297
    cell_last(j) = coord(j).cell();
2,147,483,647✔
298
  }
299
  n_coord_last() = n_coord();
2,147,483,647✔
300

301
  // Set surface that particle is on and adjust coordinate levels
302
  surface() = boundary().surface();
2,147,483,647✔
303
  n_coord() = boundary().coord_level();
2,147,483,647✔
304

305
  const auto& surf {*model::surfaces[surface_index()].get()};
2,147,483,647✔
306

307
  if (boundary().lattice_translation()[0] != 0 ||
2,147,483,647✔
308
      boundary().lattice_translation()[1] != 0 ||
2,147,483,647✔
309
      boundary().lattice_translation()[2] != 0) {
1,860,085,623✔
310
    // Particle crosses lattice boundary
311

312
    bool verbose = settings::verbosity >= 10 || trace();
750,250,860!
313
    cross_lattice(*this, boundary(), verbose);
750,250,860✔
314
    event() = TallyEvent::LATTICE;
750,250,860✔
315
  } else {
316
    // Particle crosses surface
317
    // If BC, add particle to surface source before crossing surface
318
    if (surf.surf_source_ && surf.bc_) {
1,671,543,739✔
319
      add_surf_source_to_bank(*this, surf);
720,301,239✔
320
    }
321
    this->cross_surface(surf);
1,671,543,739✔
322
    // If no BC, add particle to surface source after crossing surface
323
    if (surf.surf_source_ && !surf.bc_) {
1,671,543,730✔
324
      add_surf_source_to_bank(*this, surf);
950,004,664✔
325
    }
326
    if (settings::weight_window_checkpoint_surface) {
1,671,543,730✔
327
      apply_weight_windows(*this);
74,912✔
328
    }
329
    event() = TallyEvent::SURFACE;
1,671,543,730✔
330
  }
331
  // Score cell to cell partial currents
332
  if (!model::active_surface_tallies.empty()) {
2,147,483,647✔
333
    score_surface_tally(*this, model::active_surface_tallies, surf);
34,931,567✔
334
  }
335
}
2,147,483,647✔
336

337
void Particle::event_collide()
2,147,483,647✔
338
{
339

340
  // Score collision estimate of keff
341
  if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {
2,147,483,647✔
342
    keff_tally_collision() += wgt() * macro_xs().nu_fission / macro_xs().total;
2,139,390,648✔
343
  }
344

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

349
  if (!model::active_meshsurf_tallies.empty())
2,147,483,647✔
350
    score_meshsurface_tally(*this, model::active_meshsurf_tallies);
63,098,926✔
351

352
  // Clear surface component
353
  surface() = SURFACE_NONE;
2,147,483,647✔
354

355
  if (settings::run_CE) {
2,147,483,647✔
356
    collision(*this);
1,072,107,011✔
357
  } else {
358
    collision_mg(*this);
1,783,060,477✔
359
  }
360

361
  // Collision track feature to recording particle interaction
362
  if (settings::collision_track) {
2,147,483,647✔
363
    collision_track_record(*this);
150,087✔
364
  }
365

366
  // Score collision estimator tallies -- this is done after a collision
367
  // has occurred rather than before because we need information on the
368
  // outgoing energy for any tallies with an outgoing energy filter
369
  if (!model::active_collision_tallies.empty())
2,147,483,647✔
370
    score_collision_tally(*this);
107,372,900✔
371
  if (!model::active_analog_tallies.empty()) {
2,147,483,647✔
372
    if (settings::run_CE) {
406,034,945✔
373
      score_analog_tally_ce(*this);
404,826,683✔
374
    } else {
375
      score_analog_tally_mg(*this);
1,208,262✔
376
    }
377
  }
378

379
  if (!model::active_pulse_height_tallies.empty() && type().is_photon()) {
2,147,483,647✔
380
    pht_collision_energy();
2,024✔
381
  }
382

383
  // Reset banked weight during collision
384
  n_bank() = 0;
2,147,483,647✔
385
  bank_second_E() = 0.0;
2,147,483,647✔
386
  wgt_bank() = 0.0;
2,147,483,647✔
387

388
  // Clear number of secondaries in this collision. This is
389
  // distinct from the number of created neutrons n_bank() above!
390
  n_secondaries() = 0;
2,147,483,647✔
391

392
  zero_delayed_bank();
2,147,483,647✔
393

394
  // Reset fission logical
395
  fission() = false;
2,147,483,647✔
396

397
  // Save coordinates for tallying purposes
398
  r_last_current() = r();
2,147,483,647✔
399

400
  // Set last material to none since cross sections will need to be
401
  // re-evaluated
402
  material_last() = C_NONE;
2,147,483,647✔
403

404
  // Set all directions to base level -- right now, after a collision, only
405
  // the base level directions are changed
406
  for (int j = 0; j < n_coord() - 1; ++j) {
2,147,483,647✔
407
    if (coord(j + 1).rotated()) {
161,247,365✔
408
      // If next level is rotated, apply rotation matrix
409
      const auto& m {model::cells[coord(j).cell()]->rotation_};
10,426,614✔
410
      const auto& u {coord(j).u()};
10,426,614✔
411
      coord(j + 1).u() = u.rotate(m);
10,426,614✔
412
    } else {
413
      // Otherwise, copy this level's direction
414
      coord(j + 1).u() = coord(j).u();
150,820,751✔
415
    }
416
  }
417

418
  // Score flux derivative accumulators for differential tallies.
419
  if (!model::active_tallies.empty())
2,147,483,647✔
420
    score_collision_derivative(*this);
937,493,998✔
421

422
#ifdef OPENMC_DAGMC_ENABLED
423
  history().reset();
261,407,959✔
424
#endif
425
}
2,147,483,647✔
426

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

437
  // Check for secondary particles if this particle is dead
438
  if (!alive()) {
2,147,483,647✔
439
    // Write final position for this particle
440
    if (write_track()) {
231,933,908✔
441
      write_particle_track(*this);
6,244✔
442
    }
443

444
    // If no secondary particles, break out of event loop
445
    if (secondary_bank().empty())
231,933,908✔
446
      return;
447

448
    from_source(&secondary_bank().back());
63,131,937✔
449
    secondary_bank().pop_back();
63,131,937✔
450
    n_event() = 0;
63,131,937✔
451
    bank_second_E() = 0.0;
63,131,937✔
452

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

470
        // Initialize last cells from current cell
471
        for (int j = 0; j < n_coord(); ++j) {
1,210✔
472
          cell_last(j) = coord(j).cell();
605✔
473
        }
474
        n_coord_last() = n_coord();
605✔
475
      }
476
      pht_secondary_particles();
605✔
477
    }
478

479
    // Enter new particle in particle track file
480
    if (write_track())
63,131,937✔
481
      add_particle_track(*this);
5,234✔
482
  }
483
}
484

485
void Particle::event_death()
168,802,971✔
486
{
487
#ifdef OPENMC_DAGMC_ENABLED
488
  history().reset();
15,422,385✔
489
#endif
490

491
  // Finish particle track output.
492
  if (write_track()) {
168,802,971✔
493
    finalize_particle_track(*this);
1,010✔
494
  }
495

496
// Contribute tally reduction variables to global accumulator
497
#pragma omp atomic
93,021,287✔
498
  global_tally_absorption += keff_tally_absorption();
168,802,971✔
499
#pragma omp atomic
93,170,415✔
500
  global_tally_collision += keff_tally_collision();
168,802,971✔
501
#pragma omp atomic
93,019,002✔
502
  global_tally_tracklength += keff_tally_tracklength();
168,802,971✔
503
#pragma omp atomic
92,365,059✔
504
  global_tally_leakage += keff_tally_leakage();
168,802,971✔
505

506
  // Reset particle tallies once accumulated
507
  keff_tally_absorption() = 0.0;
168,802,971✔
508
  keff_tally_collision() = 0.0;
168,802,971✔
509
  keff_tally_tracklength() = 0.0;
168,802,971✔
510
  keff_tally_leakage() = 0.0;
168,802,971✔
511

512
  if (!model::active_pulse_height_tallies.empty()) {
168,802,971✔
513
    score_pulse_height_tally(*this, model::active_pulse_height_tallies);
5,500✔
514
  }
515

516
  // Record the number of progeny created by this particle.
517
  // This data will be used to efficiently sort the fission bank.
518
  if (settings::run_mode == RunMode::EIGENVALUE) {
168,802,971✔
519
    int64_t offset = id() - 1 - simulation::work_index[mpi::rank];
141,534,700✔
520
    simulation::progeny_per_particle[offset] = n_progeny();
141,534,700✔
521
  }
522
}
168,802,971✔
523

524
void Particle::pht_collision_energy()
2,024✔
525
{
526
  // Adds the energy particles lose in a collision to the pulse-height
527

528
  // determine index of cell in pulse_height_cells
529
  auto it = std::find(model::pulse_height_cells.begin(),
2,024✔
530
    model::pulse_height_cells.end(), lowest_coord().cell());
2,024!
531

532
  if (it != model::pulse_height_cells.end()) {
2,024!
533
    int index = std::distance(model::pulse_height_cells.begin(), it);
2,024✔
534
    pht_storage()[index] += E_last() - E();
2,024✔
535

536
    // If the energy of the particle is below the cutoff, it will not be sampled
537
    // so its energy is added to the pulse-height in the cell
538
    int photon = ParticleType::photon().transport_index();
2,024✔
539
    if (E() < settings::energy_cutoff[photon]) {
2,024✔
540
      pht_storage()[index] += E();
825✔
541
    }
542
  }
543
}
2,024✔
544

545
void Particle::pht_secondary_particles()
605✔
546
{
547
  // Removes the energy of secondary produced particles from the pulse-height
548

549
  // determine index of cell in pulse_height_cells
550
  auto it = std::find(model::pulse_height_cells.begin(),
605✔
551
    model::pulse_height_cells.end(), cell_born());
605!
552

553
  if (it != model::pulse_height_cells.end()) {
605!
554
    int index = std::distance(model::pulse_height_cells.begin(), it);
605✔
555
    pht_storage()[index] -= E();
605✔
556
  }
557
}
605✔
558

559
void Particle::cross_surface(const Surface& surf)
1,673,451,005✔
560
{
561

562
  if (settings::verbosity >= 10 || trace()) {
1,673,451,005✔
563
    write_message(1, "    Crossing surface {}", surf.id_);
66✔
564
  }
565

566
// if we're crossing a CSG surface, make sure the DAG history is reset
567
#ifdef OPENMC_DAGMC_ENABLED
568
  if (surf.geom_type() == GeometryType::CSG)
152,817,327✔
569
    history().reset();
152,760,168✔
570
#endif
571

572
  // Handle any applicable boundary conditions.
573
  if (surf.bc_ && settings::run_mode != RunMode::PLOTTING &&
1,673,451,005!
574
      settings::run_mode != RunMode::VOLUME) {
575
    surf.bc_->handle_particle(*this, surf);
720,653,347✔
576
    return;
720,653,347✔
577
  }
578

579
  // ==========================================================================
580
  // SEARCH NEIGHBOR LISTS FOR NEXT CELL
581

582
#ifdef OPENMC_DAGMC_ENABLED
583
  // in DAGMC, we know what the next cell should be
584
  if (surf.geom_type() == GeometryType::DAG) {
86,858,510✔
585
    int32_t i_cell = next_cell(surface_index(), cell_last(n_coord() - 1),
46,350✔
586
                       lowest_coord().universe()) -
46,350✔
587
                     1;
46,350✔
588
    // save material, temperature, and density multiplier
589
    material_last() = material();
46,350✔
590
    sqrtkT_last() = sqrtkT();
46,350✔
591
    density_mult_last() = density_mult();
46,350✔
592
    // set new cell value
593
    lowest_coord().cell() = i_cell;
46,350✔
594
    auto& cell = model::cells[i_cell];
46,350✔
595

596
    cell_instance() = 0;
46,350✔
597
    if (cell->distribcell_index_ >= 0)
46,350✔
598
      cell_instance() = cell_instance_at_level(*this, n_coord() - 1);
45,326✔
599

600
    material() = cell->material(cell_instance());
46,350!
601
    sqrtkT() = cell->sqrtkT(cell_instance());
46,350!
602
    density_mult() = cell->density_mult(cell_instance());
46,350✔
603
    return;
46,350✔
604
  }
605
#endif
606

607
  bool verbose = settings::verbosity >= 10 || trace();
952,751,308!
608
  if (neighbor_list_find_cell(*this, verbose)) {
952,751,308✔
609
    return;
610
  }
611

612
  // ==========================================================================
613
  // COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS
614

615
  // Remove lower coordinate levels
616
  n_coord() = 1;
29,911✔
617
  bool found = exhaustive_find_cell(*this, verbose);
29,911✔
618

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

625
    surface() = SURFACE_NONE;
5,799✔
626
    n_coord() = 1;
5,799✔
627
    r() += TINY_BIT * u();
5,799✔
628

629
    // Couldn't find next cell anywhere! This probably means there is an actual
630
    // undefined region in the geometry.
631

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

641
void Particle::cross_vacuum_bc(const Surface& surf)
35,093,663✔
642
{
643
  // Score any surface current tallies -- note that the particle is moved
644
  // forward slightly so that if the mesh boundary is on the surface, it is
645
  // still processed
646

647
  if (!model::active_meshsurf_tallies.empty()) {
35,093,663✔
648
    // TODO: Find a better solution to score surface currents than
649
    // physically moving the particle forward slightly
650

651
    r() += TINY_BIT * u();
937,222✔
652
    score_meshsurface_tally(*this, model::active_meshsurf_tallies);
937,222✔
653
  }
654

655
  // Score to global leakage tally
656
  keff_tally_leakage() += wgt();
35,093,663✔
657

658
  // Kill the particle
659
  wgt() = 0.0;
35,093,663✔
660

661
  // Display message
662
  if (settings::verbosity >= 10 || trace()) {
35,093,663!
663
    write_message(1, "    Leaked out of surface {}", surf.id_);
22✔
664
  }
665
}
35,093,663✔
666

667
void Particle::cross_reflective_bc(const Surface& surf, Direction new_u)
684,318,903✔
668
{
669
  // Do not handle reflective boundary conditions on lower universes
670
  if (n_coord() != 1) {
684,318,903!
UNCOV
671
    mark_as_lost("Cannot reflect particle " + std::to_string(id()) +
×
672
                 " off surface in a lower universe.");
UNCOV
673
    return;
×
674
  }
675

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

683
  if (!model::active_surface_tallies.empty()) {
684,318,903✔
684
    score_surface_tally(*this, model::active_surface_tallies, surf);
285,021✔
685
  }
686

687
  if (!model::active_meshsurf_tallies.empty()) {
684,318,903✔
688
    Position r {this->r()};
46,885,487✔
689
    this->r() -= TINY_BIT * u();
46,885,487✔
690
    score_meshsurface_tally(*this, model::active_meshsurf_tallies);
46,885,487✔
691
    this->r() = r;
46,885,487✔
692
  }
693

694
  // Set the new particle direction
695
  u() = new_u;
684,318,903✔
696

697
  // Reassign particle's cell and surface
698
  coord(0).cell() = cell_last(0);
684,318,903✔
699
  surface() = -surface();
684,318,903✔
700

701
  // If a reflective surface is coincident with a lattice or universe
702
  // boundary, it is necessary to redetermine the particle's coordinates in
703
  // the lower universes.
704
  // (unless we're using a dagmc model, which has exactly one universe)
705
  n_coord() = 1;
684,318,903✔
706
  if (surf.geom_type() != GeometryType::DAG &&
1,368,635,048!
707
      !neighbor_list_find_cell(*this)) {
684,316,145✔
UNCOV
708
    mark_as_lost("Couldn't find particle after reflecting from surface " +
×
UNCOV
709
                 std::to_string(surf.id_) + ".");
×
UNCOV
710
    return;
×
711
  }
712

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

716
  // Diagnostic message
717
  if (settings::verbosity >= 10 || trace()) {
684,318,903!
UNCOV
718
    write_message(1, "    Reflected from surface {}", surf.id_);
×
719
  }
720
}
721

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

734
  // Score surface currents since reflection causes the direction of the
735
  // particle to change -- artificially move the particle slightly back in
736
  // case the surface crossing is coincident with a mesh boundary
737
  if (!model::active_meshsurf_tallies.empty()) {
2,246,247!
UNCOV
738
    Position r {this->r()};
×
UNCOV
739
    this->r() -= TINY_BIT * u();
×
UNCOV
740
    score_meshsurface_tally(*this, model::active_meshsurf_tallies);
×
UNCOV
741
    this->r() = r;
×
742
  }
743

744
  // Adjust the particle's location and direction.
745
  r() = new_r;
2,246,247✔
746
  u() = new_u;
2,246,247✔
747

748
  // Reassign particle's surface
749
  surface() = new_surface;
2,246,247✔
750

751
  // Figure out what cell particle is in now
752
  n_coord() = 1;
2,246,247✔
753

754
  if (!neighbor_list_find_cell(*this)) {
2,246,247!
UNCOV
755
    mark_as_lost("Couldn't find particle after hitting periodic "
×
UNCOV
756
                 "boundary on surface " +
×
UNCOV
757
                 std::to_string(surf.id_) + ".");
×
UNCOV
758
    return;
×
759
  }
760

761
  // Set previous coordinate going slightly past surface crossing
762
  r_last_current() = r() + TINY_BIT * u();
2,246,247✔
763

764
  // Diagnostic message
765
  if (settings::verbosity >= 10 || trace()) {
2,246,247!
UNCOV
766
    write_message(1, "    Hit periodic boundary on surface {}", surf.id_);
×
767
  }
768
}
769

770
void Particle::mark_as_lost(const char* message)
5,799✔
771
{
772
  // Print warning and write lost particle file
773
  warning(message);
5,799✔
774
  if (settings::max_write_lost_particles < 0 ||
5,799✔
775
      simulation::n_lost_particles < settings::max_write_lost_particles) {
5,500✔
776
    write_restart();
374✔
777
  }
778
  // Increment number of lost particles
779
  wgt() = 0.0;
5,799✔
780
#pragma omp atomic
3,154✔
781
  simulation::n_lost_particles += 1;
2,645✔
782

783
  // Count the total number of simulated particles (on this processor)
784
  auto n = simulation::current_batch * settings::gen_per_batch *
5,799✔
785
           simulation::work_per_rank;
786

787
  // Abort the simulation if the maximum number of lost particles has been
788
  // reached
789
  if (simulation::n_lost_particles >= settings::max_lost_particles &&
5,799✔
790
      simulation::n_lost_particles >= settings::rel_max_lost_particles * n) {
9!
791
    fatal_error("Maximum number of lost particles has been reached.");
9✔
792
  }
793
}
5,790✔
794

795
void Particle::write_restart() const
374✔
796
{
797
  // Dont write another restart file if in particle restart mode
798
  if (settings::run_mode == RunMode::PARTICLE)
374✔
799
    return;
22✔
800

801
  // Set up file name
802
  auto filename = fmt::format("{}particle_{}_{}.h5", settings::path_output,
352✔
803
    simulation::current_batch, id());
352✔
804

805
#pragma omp critical(WriteParticleRestart)
187✔
806
  {
352✔
807
    // Create file
808
    hid_t file_id = file_open(filename, 'w');
352✔
809

810
    // Write filetype and version info
811
    write_attribute(file_id, "filetype", "particle restart");
352✔
812
    write_attribute(file_id, "version", VERSION_PARTICLE_RESTART);
352✔
813
    write_attribute(file_id, "openmc_version", VERSION);
352✔
814
#ifdef GIT_SHA1
815
    write_attr_string(file_id, "git_sha1", GIT_SHA1);
816
#endif
817

818
    // Write data to file
819
    write_dataset(file_id, "current_batch", simulation::current_batch);
352✔
820
    write_dataset(file_id, "generations_per_batch", settings::gen_per_batch);
352✔
821
    write_dataset(file_id, "current_generation", simulation::current_gen);
352✔
822
    write_dataset(file_id, "n_particles", settings::n_particles);
352✔
823
    switch (settings::run_mode) {
352!
824
    case RunMode::FIXED_SOURCE:
220✔
825
      write_dataset(file_id, "run_mode", "fixed source");
220✔
826
      break;
115✔
827
    case RunMode::EIGENVALUE:
132✔
828
      write_dataset(file_id, "run_mode", "eigenvalue");
132✔
829
      break;
72✔
UNCOV
830
    case RunMode::PARTICLE:
×
UNCOV
831
      write_dataset(file_id, "run_mode", "particle restart");
×
832
      break;
833
    default:
834
      break;
835
    }
836
    write_dataset(file_id, "id", id());
352✔
837
    write_dataset(file_id, "type", type().pdg_number());
352✔
838

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

862
    // Close file
863
    file_close(file_id);
352✔
864
  } // #pragma omp critical
865
}
352✔
866

867
void Particle::update_neutron_xs(
2,147,483,647✔
868
  int i_nuclide, int i_grid, int i_sab, double sab_frac, double ncrystal_xs)
869
{
870
  // Get microscopic cross section cache
871
  auto& micro = this->neutron_xs(i_nuclide);
2,147,483,647✔
872

873
  // If the cache doesn't match, recalculate micro xs
874
  if (this->E() != micro.last_E || this->sqrtkT() != micro.last_sqrtkT ||
2,147,483,647✔
875
      i_sab != micro.index_sab || sab_frac != micro.sab_frac ||
2,147,483,647✔
876
      ncrystal_xs != micro.ncrystal_xs) {
2,147,483,647!
877
    data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, *this);
2,147,483,647✔
878

879
    // If NCrystal is being used, update micro cross section cache
880
    micro.ncrystal_xs = ncrystal_xs;
2,147,483,647✔
881
    if (ncrystal_xs >= 0.0) {
2,147,483,647✔
882
      data::nuclides[i_nuclide]->calculate_elastic_xs(*this);
11,018,953✔
883
      ncrystal_update_micro(ncrystal_xs, micro);
11,018,953✔
884
    }
885
  }
886
}
2,147,483,647✔
887

888
//==============================================================================
889
// Non-method functions
890
//==============================================================================
891
void add_surf_source_to_bank(Particle& p, const Surface& surf)
1,670,305,903✔
892
{
893
  if (simulation::current_batch <= settings::n_inactive ||
1,670,305,903✔
894
      simulation::surf_source_bank.full()) {
1,284,835,522✔
895
    return;
1,670,176,250✔
896
  }
897

898
  // If a cell/cellfrom/cellto parameter is defined
899
  if (settings::ssw_cell_id != C_NONE) {
337,080✔
900

901
    // Retrieve cell index and storage type
902
    int cell_idx = model::cell_map[settings::ssw_cell_id];
254,435✔
903

904
    if (surf.bc_) {
254,435✔
905
      // Leave if cellto with vacuum boundary condition
906
      if (surf.bc_->type() == "vacuum" &&
298,918✔
907
          settings::ssw_cell_type == SSWCellType::To) {
33,099✔
908
        return;
909
      }
910

911
      // Leave if other boundary condition than vacuum
912
      if (surf.bc_->type() != "vacuum") {
274,648✔
913
        return;
914
      }
915
    }
916

917
    // Check if the cell of interest has been exited
918
    bool exited = false;
919
    for (int i = 0; i < p.n_coord_last(); ++i) {
333,669✔
920
      if (p.cell_last(i) == cell_idx) {
207,729✔
921
        exited = true;
73,765✔
922
      }
923
    }
924

925
    // Check if the cell of interest has been entered
926
    bool entered = false;
927
    for (int i = 0; i < p.n_coord(); ++i) {
297,971✔
928
      if (p.coord(i).cell() == cell_idx) {
172,031✔
929
        entered = true;
57,515✔
930
      }
931
    }
932

933
    // Vacuum boundary conditions: return if cell is not exited
934
    if (surf.bc_) {
125,940✔
935
      if (surf.bc_->type() == "vacuum" && !exited) {
41,928!
936
        return;
937
      }
938
    } else {
939

940
      // If we both enter and exit the cell of interest
941
      if (entered && exited) {
104,976✔
942
        return;
943
      }
944

945
      // If we did not enter nor exit the cell of interest
946
      if (!entered && !exited) {
77,773✔
947
        return;
948
      }
949

950
      // If cellfrom and the cell before crossing is not the cell of
951
      // interest
952
      if (settings::ssw_cell_type == SSWCellType::From && !exited) {
64,274✔
953
        return;
954
      }
955

956
      // If cellto and the cell after crossing is not the cell of interest
957
      if (settings::ssw_cell_type == SSWCellType::To && !entered) {
52,733✔
958
        return;
959
      }
960
    }
961
  }
962

963
  SourceSite site;
129,653✔
964
  site.r = p.r();
129,653✔
965
  site.u = p.u();
129,653✔
966
  site.E = p.E();
129,653✔
967
  site.time = p.time();
129,653✔
968
  site.wgt = p.wgt();
129,653✔
969
  site.delayed_group = p.delayed_group();
129,653✔
970
  site.surf_id = surf.id_;
129,653✔
971
  site.particle = p.type();
129,653✔
972
  site.parent_id = p.id();
129,653✔
973
  site.progeny_id = p.n_progeny();
129,653✔
974
  int64_t idx = simulation::surf_source_bank.thread_safe_append(site);
129,653✔
975
}
976

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