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

openmc-dev / openmc / 19288072129

12 Nov 2025 06:05AM UTC coverage: 82.024% (+0.02%) from 82.006%
19288072129

Pull #3624

github

web-flow
Merge 2c40f0c4b into 8cd3911cb
Pull Request #3624: Avoid divide-by-zero in `from_multigroup_flux` when flux is zero

16723 of 23233 branches covered (71.98%)

Branch coverage included in aggregate %.

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

29 existing lines in 1 file now uncovered.

54251 of 63295 relevant lines covered (85.71%)

43808177.67 hits per line

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

85.35
/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
//==============================================================================
42
// Particle implementation
43
//==============================================================================
44

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

74
bool Particle::create_secondary(
113,024,402✔
75
  double wgt, Direction u, double E, ParticleType type)
76
{
77
  // If energy is below cutoff for this particle, don't create secondary
78
  // particle
79
  if (E < settings::energy_cutoff[static_cast<int>(type)]) {
113,024,402✔
80
    return false;
54,031,267✔
81
  }
82

83
  auto& bank = secondary_bank().emplace_back();
58,993,135✔
84
  bank.particle = type;
58,993,135✔
85
  bank.wgt = wgt;
58,993,135✔
86
  bank.r = r();
58,993,135✔
87
  bank.u = u;
58,993,135✔
88
  bank.E = settings::run_CE ? E : g();
58,993,135!
89
  bank.time = time();
58,993,135✔
90
  bank_second_E() += bank.E;
58,993,135✔
91
  return true;
58,993,135✔
92
}
93

94
void Particle::split(double wgt)
4,410,895✔
95
{
96
  auto& bank = secondary_bank().emplace_back();
4,410,895✔
97
  bank.particle = type();
4,410,895✔
98
  bank.wgt = wgt;
4,410,895✔
99
  bank.r = r();
4,410,895✔
100
  bank.u = u();
4,410,895✔
101
  bank.E = settings::run_CE ? E() : g();
4,410,895✔
102
  bank.time = time();
4,410,895✔
103

104
  // Convert signed index to a signed surface ID
105
  if (surface() == SURFACE_NONE) {
4,410,895✔
106
    bank.surf_id = SURFACE_NONE;
4,410,559✔
107
  } else {
108
    int surf_id = model::surfaces[surface_index()]->id_;
336✔
109
    bank.surf_id = (surface() > 0) ? surf_id : -surf_id;
336!
110
  }
111
}
4,410,895✔
112

113
void Particle::from_source(const SourceSite* src)
235,982,200✔
114
{
115
  // Reset some attributes
116
  clear();
235,982,200✔
117
  surface() = SURFACE_NONE;
235,982,200✔
118
  cell_born() = C_NONE;
235,982,200✔
119
  material() = C_NONE;
235,982,200✔
120
  n_collision() = 0;
235,982,200✔
121
  fission() = false;
235,982,200✔
122
  zero_flux_derivs();
235,982,200✔
123
  lifetime() = 0.0;
235,982,200✔
124
#ifdef OPENMC_DAGMC_ENABLED
125
  history().reset();
21,476,341✔
126
#endif
127

128
  // Copy attributes from source bank site
129
  type() = src->particle;
235,982,200✔
130
  wgt() = src->wgt;
235,982,200✔
131
  wgt_last() = src->wgt;
235,982,200✔
132
  r() = src->r;
235,982,200✔
133
  u() = src->u;
235,982,200✔
134
  r_born() = src->r;
235,982,200✔
135
  r_last_current() = src->r;
235,982,200✔
136
  r_last() = src->r;
235,982,200✔
137
  u_last() = src->u;
235,982,200✔
138
  if (settings::run_CE) {
235,982,200✔
139
    E() = src->E;
120,341,483✔
140
    g() = 0;
120,341,483✔
141
  } else {
142
    g() = static_cast<int>(src->E);
115,640,717✔
143
    g_last() = static_cast<int>(src->E);
115,640,717✔
144
    E() = data::mg.energy_bin_avg_[g()];
115,640,717✔
145
  }
146
  E_last() = E();
235,982,200✔
147
  time() = src->time;
235,982,200✔
148
  time_last() = src->time;
235,982,200✔
149
  parent_nuclide() = src->parent_nuclide;
235,982,200✔
150
  delayed_group() = src->delayed_group;
235,982,200✔
151

152
  // Convert signed surface ID to signed index
153
  if (src->surf_id != SURFACE_NONE) {
235,982,200✔
154
    int index_plus_one = model::surface_map[std::abs(src->surf_id)] + 1;
110,336✔
155
    surface() = (src->surf_id > 0) ? index_plus_one : -index_plus_one;
110,336!
156
  }
157
}
235,982,200✔
158

159
void Particle::event_calculate_xs()
2,147,483,647✔
160
{
161
  // Set the random number stream
162
  stream() = STREAM_TRACKING;
2,147,483,647✔
163

164
  // Store pre-collision particle properties
165
  wgt_last() = wgt();
2,147,483,647✔
166
  E_last() = E();
2,147,483,647✔
167
  u_last() = u();
2,147,483,647✔
168
  r_last() = r();
2,147,483,647✔
169
  time_last() = time();
2,147,483,647✔
170

171
  // Reset event variables
172
  event() = TallyEvent::KILL;
2,147,483,647✔
173
  event_nuclide() = NUCLIDE_NONE;
2,147,483,647✔
174
  event_mt() = REACTION_NONE;
2,147,483,647✔
175

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

186
    // Set birth cell attribute
187
    if (cell_born() == C_NONE)
232,426,003!
188
      cell_born() = lowest_coord().cell();
232,426,003✔
189

190
    // Initialize last cells from current cell
191
    for (int j = 0; j < n_coord(); ++j) {
482,320,818✔
192
      cell_last(j) = coord(j).cell();
249,894,815✔
193
    }
194
    n_coord_last() = n_coord();
232,426,003✔
195
  }
196

197
  // Write particle track.
198
  if (write_track())
2,147,483,647✔
199
    write_particle_track(*this);
10,814✔
200

201
  if (settings::check_overlaps)
2,147,483,647!
UNCOV
202
    check_cell_overlap(*this);
×
203

204
  // Calculate microscopic and macroscopic cross sections
205
  if (material() != MATERIAL_VOID) {
2,147,483,647✔
206
    if (settings::run_CE) {
2,147,483,647✔
207
      if (material() != material_last() || sqrtkT() != sqrtkT_last() ||
2,147,483,647✔
208
          density_mult() != density_mult_last()) {
400,750,860✔
209
        // If the material is the same as the last material and the
210
        // temperature hasn't changed, we don't need to lookup cross
211
        // sections again.
212
        model::materials[material()]->calculate_xs(*this);
1,549,830,756✔
213
      }
214
    } else {
215
      // Get the MG data; unlike the CE case above, we have to re-calculate
216
      // cross sections for every collision since the cross sections may
217
      // be angle-dependent
218
      data::mg.macro_xs_[material()].calculate_xs(*this);
2,063,937,414✔
219

220
      // Update the particle's group while we know we are multi-group
221
      g_last() = g();
2,063,937,414✔
222
    }
223
  } else {
224
    macro_xs().total = 0.0;
67,553,055✔
225
    macro_xs().absorption = 0.0;
67,553,055✔
226
    macro_xs().fission = 0.0;
67,553,055✔
227
    macro_xs().nu_fission = 0.0;
67,553,055✔
228
  }
229
}
230

231
void Particle::event_advance()
2,147,483,647✔
232
{
233
  // Find the distance to the nearest boundary
234
  boundary() = distance_to_boundary(*this);
2,147,483,647✔
235

236
  // Sample a distance to collision
237
  if (type() == ParticleType::electron || type() == ParticleType::positron) {
2,147,483,647✔
238
    collision_distance() = material() == MATERIAL_VOID ? INFINITY : 0.0;
54,187,846!
239
  } else if (macro_xs().total == 0.0) {
2,147,483,647✔
240
    collision_distance() = INFINITY;
67,553,055✔
241
  } else {
242
    collision_distance() = -std::log(prn(current_seed())) / macro_xs().total;
2,147,483,647✔
243
  }
244

245
  double speed = this->speed();
2,147,483,647✔
246
  double time_cutoff = settings::time_cutoff[static_cast<int>(type())];
2,147,483,647✔
247
  double distance_cutoff =
248
    (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY;
2,147,483,647✔
249

250
  // Select smaller of the three distances
251
  double distance =
252
    std::min({boundary().distance(), collision_distance(), distance_cutoff});
2,147,483,647✔
253

254
  // Advance particle in space and time
255
  this->move_distance(distance);
2,147,483,647✔
256
  double dt = distance / speed;
2,147,483,647✔
257
  this->time() += dt;
2,147,483,647✔
258
  this->lifetime() += dt;
2,147,483,647✔
259

260
  // Score timed track-length tallies
261
  if (!model::active_timed_tracklength_tallies.empty()) {
2,147,483,647✔
262
    score_timed_tracklength_tally(*this, distance);
3,628,317✔
263
  }
264

265
  // Score track-length tallies
266
  if (!model::active_tracklength_tallies.empty()) {
2,147,483,647✔
267
    score_tracklength_tally(*this, distance);
1,473,589,016✔
268
  }
269

270
  // Score track-length estimate of k-eff
271
  if (settings::run_mode == RunMode::EIGENVALUE &&
2,147,483,647✔
272
      type() == ParticleType::neutron) {
2,147,483,647✔
273
    keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission;
2,147,483,647✔
274
  }
275

276
  // Score flux derivative accumulators for differential tallies.
277
  if (!model::active_tallies.empty()) {
2,147,483,647✔
278
    score_track_derivative(*this, distance);
1,670,879,011✔
279
  }
280

281
  // Set particle weight to zero if it hit the time boundary
282
  if (distance == distance_cutoff) {
2,147,483,647✔
283
    wgt() = 0.0;
224,928✔
284
  }
285
}
2,147,483,647✔
286

287
void Particle::event_cross_surface()
2,147,483,647✔
288
{
289
  // Saving previous cell data
290
  for (int j = 0; j < n_coord(); ++j) {
2,147,483,647✔
291
    cell_last(j) = coord(j).cell();
2,147,483,647✔
292
  }
293
  n_coord_last() = n_coord();
2,147,483,647✔
294

295
  // Set surface that particle is on and adjust coordinate levels
296
  surface() = boundary().surface();
2,147,483,647✔
297
  n_coord() = boundary().coord_level();
2,147,483,647✔
298

299
  if (boundary().lattice_translation()[0] != 0 ||
2,147,483,647✔
300
      boundary().lattice_translation()[1] != 0 ||
2,147,483,647✔
301
      boundary().lattice_translation()[2] != 0) {
1,736,747,164✔
302
    // Particle crosses lattice boundary
303

304
    bool verbose = settings::verbosity >= 10 || trace();
686,024,086!
305
    cross_lattice(*this, boundary(), verbose);
686,024,086✔
306
    event() = TallyEvent::LATTICE;
686,024,086✔
307
  } else {
308
    // Particle crosses surface
309
    const auto& surf {model::surfaces[surface_index()].get()};
1,549,289,525✔
310
    // If BC, add particle to surface source before crossing surface
311
    if (surf->surf_source_ && surf->bc_) {
1,549,289,525✔
312
      add_surf_source_to_bank(*this, *surf);
708,444,281✔
313
    }
314
    this->cross_surface(*surf);
1,549,289,525✔
315
    // If no BC, add particle to surface source after crossing surface
316
    if (surf->surf_source_ && !surf->bc_) {
1,549,289,516✔
317
      add_surf_source_to_bank(*this, *surf);
839,607,408✔
318
    }
319
    if (settings::weight_window_checkpoint_surface) {
1,549,289,516✔
320
      apply_weight_windows(*this);
10,738!
321
    }
322
    event() = TallyEvent::SURFACE;
1,549,289,516✔
323
  }
324
  // Score cell to cell partial currents
325
  if (!model::active_surface_tallies.empty()) {
2,147,483,647✔
326
    score_surface_tally(*this, model::active_surface_tallies);
34,922,767✔
327
  }
328
}
2,147,483,647✔
329

330
void Particle::event_collide()
2,147,483,647✔
331
{
332
  // Score collision estimate of keff
333
  if (settings::run_mode == RunMode::EIGENVALUE &&
2,147,483,647✔
334
      type() == ParticleType::neutron) {
2,147,483,647✔
335
    keff_tally_collision() += wgt() * macro_xs().nu_fission / macro_xs().total;
2,147,483,647✔
336
  }
337

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

342
  if (!model::active_meshsurf_tallies.empty())
2,147,483,647✔
343
    score_surface_tally(*this, model::active_meshsurf_tallies);
85,766,962✔
344

345
  // Clear surface component
346
  surface() = SURFACE_NONE;
2,147,483,647✔
347

348
  if (settings::run_CE) {
2,147,483,647✔
349
    collision(*this);
831,969,618✔
350
  } else {
351
    collision_mg(*this);
1,783,060,477✔
352
  }
353

354
  // Score collision estimator tallies -- this is done after a collision
355
  // has occurred rather than before because we need information on the
356
  // outgoing energy for any tallies with an outgoing energy filter
357
  if (!model::active_collision_tallies.empty())
2,147,483,647✔
358
    score_collision_tally(*this);
108,083,562✔
359
  if (!model::active_analog_tallies.empty()) {
2,147,483,647✔
360
    if (settings::run_CE) {
130,978,194✔
361
      score_analog_tally_ce(*this);
129,769,932✔
362
    } else {
363
      score_analog_tally_mg(*this);
1,208,262✔
364
    }
365
  }
366

367
  if (!model::active_pulse_height_tallies.empty() &&
2,147,483,647✔
368
      type() == ParticleType::photon) {
16,918✔
369
    pht_collision_energy();
2,024✔
370
  }
371

372
  // Reset banked weight during collision
373
  n_bank() = 0;
2,147,483,647✔
374
  bank_second_E() = 0.0;
2,147,483,647✔
375
  wgt_bank() = 0.0;
2,147,483,647✔
376
  zero_delayed_bank();
2,147,483,647✔
377

378
  // Reset fission logical
379
  fission() = false;
2,147,483,647✔
380

381
  // Save coordinates for tallying purposes
382
  r_last_current() = r();
2,147,483,647✔
383

384
  // Set last material to none since cross sections will need to be
385
  // re-evaluated
386
  material_last() = C_NONE;
2,147,483,647✔
387

388
  // Set all directions to base level -- right now, after a collision, only
389
  // the base level directions are changed
390
  for (int j = 0; j < n_coord() - 1; ++j) {
2,147,483,647✔
391
    if (coord(j + 1).rotated()) {
125,005,150✔
392
      // If next level is rotated, apply rotation matrix
393
      const auto& m {model::cells[coord(j).cell()]->rotation_};
10,426,614✔
394
      const auto& u {coord(j).u()};
10,426,614✔
395
      coord(j + 1).u() = u.rotate(m);
10,426,614✔
396
    } else {
397
      // Otherwise, copy this level's direction
398
      coord(j + 1).u() = coord(j).u();
114,578,536✔
399
    }
400
  }
401

402
  // Score flux derivative accumulators for differential tallies.
403
  if (!model::active_tallies.empty())
2,147,483,647✔
404
    score_collision_derivative(*this);
708,906,088✔
405

406
#ifdef OPENMC_DAGMC_ENABLED
407
  history().reset();
234,145,623✔
408
#endif
409
}
2,147,483,647✔
410

411
void Particle::event_revive_from_secondary()
2,147,483,647✔
412
{
413
  // If particle has too many events, display warning and kill it
414
  ++n_event();
2,147,483,647✔
415
  if (n_event() == settings::max_particle_events) {
2,147,483,647!
UNCOV
416
    warning("Particle " + std::to_string(id()) +
×
417
            " underwent maximum number of events.");
UNCOV
418
    wgt() = 0.0;
×
419
  }
420

421
  // Check for secondary particles if this particle is dead
422
  if (!alive()) {
2,147,483,647✔
423
    // Write final position for this particle
424
    if (write_track()) {
232,425,599✔
425
      write_particle_track(*this);
6,674✔
426
    }
427

428
    // If no secondary particles, break out of event loop
429
    if (secondary_bank().empty())
232,425,599✔
430
      return;
168,682,838✔
431

432
    from_source(&secondary_bank().back());
63,742,761✔
433
    secondary_bank().pop_back();
63,742,761✔
434
    n_event() = 0;
63,742,761✔
435
    bank_second_E() = 0.0;
63,742,761✔
436

437
    // Subtract secondary particle energy from interim pulse-height results
438
    if (!model::active_pulse_height_tallies.empty() &&
63,758,260✔
439
        this->type() == ParticleType::photon) {
15,499✔
440
      // Since the birth cell of the particle has not been set we
441
      // have to determine it before the energy of the secondary particle can be
442
      // removed from the pulse-height of this cell.
443
      if (lowest_coord().cell() == C_NONE) {
605!
444
        bool verbose = settings::verbosity >= 10 || trace();
605!
445
        if (!exhaustive_find_cell(*this, verbose)) {
605!
UNCOV
446
          mark_as_lost("Could not find the cell containing particle " +
×
UNCOV
447
                       std::to_string(id()));
×
UNCOV
448
          return;
×
449
        }
450
        // Set birth cell attribute
451
        if (cell_born() == C_NONE)
605!
452
          cell_born() = lowest_coord().cell();
605✔
453

454
        // Initialize last cells from current cell
455
        for (int j = 0; j < n_coord(); ++j) {
1,210✔
456
          cell_last(j) = coord(j).cell();
605✔
457
        }
458
        n_coord_last() = n_coord();
605✔
459
      }
460
      pht_secondary_particles();
605✔
461
    }
462

463
    // Enter new particle in particle track file
464
    if (write_track())
63,742,761✔
465
      add_particle_track(*this);
5,604✔
466
  }
467
}
468

469
void Particle::event_death()
168,683,838✔
470
{
471
#ifdef OPENMC_DAGMC_ENABLED
472
  history().reset();
15,351,522✔
473
#endif
474

475
  // Finish particle track output.
476
  if (write_track()) {
168,683,838✔
477
    finalize_particle_track(*this);
1,070✔
478
  }
479

480
// Contribute tally reduction variables to global accumulator
481
#pragma omp atomic
92,810,087✔
482
  global_tally_absorption += keff_tally_absorption();
168,683,838✔
483
#pragma omp atomic
92,623,729✔
484
  global_tally_collision += keff_tally_collision();
168,683,838✔
485
#pragma omp atomic
92,585,773✔
486
  global_tally_tracklength += keff_tally_tracklength();
168,683,838✔
487
#pragma omp atomic
92,302,582✔
488
  global_tally_leakage += keff_tally_leakage();
168,683,838✔
489

490
  // Reset particle tallies once accumulated
491
  keff_tally_absorption() = 0.0;
168,683,838✔
492
  keff_tally_collision() = 0.0;
168,683,838✔
493
  keff_tally_tracklength() = 0.0;
168,683,838✔
494
  keff_tally_leakage() = 0.0;
168,683,838✔
495

496
  if (!model::active_pulse_height_tallies.empty()) {
168,683,838✔
497
    score_pulse_height_tally(*this, model::active_pulse_height_tallies);
5,500✔
498
  }
499

500
  // Record the number of progeny created by this particle.
501
  // This data will be used to efficiently sort the fission bank.
502
  if (settings::run_mode == RunMode::EIGENVALUE) {
168,683,838✔
503
    int64_t offset = id() - 1 - simulation::work_index[mpi::rank];
142,040,200✔
504
    simulation::progeny_per_particle[offset] = n_progeny();
142,040,200✔
505
  }
506
}
168,683,838✔
507

508
void Particle::pht_collision_energy()
2,024✔
509
{
510
  // Adds the energy particles lose in a collision to the pulse-height
511

512
  // determine index of cell in pulse_height_cells
513
  auto it = std::find(model::pulse_height_cells.begin(),
2,024✔
514
    model::pulse_height_cells.end(), lowest_coord().cell());
2,024✔
515

516
  if (it != model::pulse_height_cells.end()) {
2,024!
517
    int index = std::distance(model::pulse_height_cells.begin(), it);
2,024✔
518
    pht_storage()[index] += E_last() - E();
2,024✔
519

520
    // If the energy of the particle is below the cutoff, it will not be sampled
521
    // so its energy is added to the pulse-height in the cell
522
    int photon = static_cast<int>(ParticleType::photon);
2,024✔
523
    if (E() < settings::energy_cutoff[photon]) {
2,024✔
524
      pht_storage()[index] += E();
825✔
525
    }
526
  }
527
}
2,024✔
528

529
void Particle::pht_secondary_particles()
605✔
530
{
531
  // Removes the energy of secondary produced particles from the pulse-height
532

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

537
  if (it != model::pulse_height_cells.end()) {
605!
538
    int index = std::distance(model::pulse_height_cells.begin(), it);
605✔
539
    pht_storage()[index] -= E();
605✔
540
  }
541
}
605✔
542

543
void Particle::cross_surface(const Surface& surf)
1,550,885,255✔
544
{
545

546
  if (settings::verbosity >= 10 || trace()) {
1,550,885,255✔
547
    write_message(1, "    Crossing surface {}", surf.id_);
33✔
548
  }
549

550
// if we're crossing a CSG surface, make sure the DAG history is reset
551
#ifdef OPENMC_DAGMC_ENABLED
552
  if (surf.geom_type() == GeometryType::CSG)
125,365,187✔
553
    history().reset();
125,310,068✔
554
#endif
555

556
  // Handle any applicable boundary conditions.
557
  if (surf.bc_ && settings::run_mode != RunMode::PLOTTING &&
2,147,483,647!
558
      settings::run_mode != RunMode::VOLUME) {
708,959,949✔
559
    surf.bc_->handle_particle(*this, surf);
708,796,389✔
560
    return;
708,796,389✔
561
  }
562

563
  // ==========================================================================
564
  // SEARCH NEIGHBOR LISTS FOR NEXT CELL
565

566
#ifdef OPENMC_DAGMC_ENABLED
567
  // in DAGMC, we know what the next cell should be
568
  if (surf.geom_type() == GeometryType::DAG) {
66,401,671✔
569
    int32_t i_cell = next_cell(surface_index(), cell_last(n_coord() - 1),
44,310✔
570
                       lowest_coord().universe()) -
44,310✔
571
                     1;
44,310✔
572
    // save material, temperature, and density multiplier
573
    material_last() = material();
44,310✔
574
    sqrtkT_last() = sqrtkT();
44,310✔
575
    density_mult_last() = density_mult();
44,310✔
576
    // set new cell value
577
    lowest_coord().cell() = i_cell;
44,310✔
578
    auto& cell = model::cells[i_cell];
44,310✔
579

580
    cell_instance() = 0;
44,310✔
581
    if (cell->distribcell_index_ >= 0)
44,310✔
582
      cell_instance() = cell_instance_at_level(*this, n_coord() - 1);
43,286✔
583

584
    material() = cell->material(cell_instance());
44,310✔
585
    sqrtkT() = cell->sqrtkT(cell_instance());
44,310✔
586
    density_mult() = cell->density_mult(cell_instance());
44,310✔
587
    return;
44,310✔
588
  }
589
#endif
590

591
  bool verbose = settings::verbosity >= 10 || trace();
842,044,556!
592
  if (neighbor_list_find_cell(*this, verbose)) {
842,044,556✔
593
    return;
842,014,645✔
594
  }
595

596
  // ==========================================================================
597
  // COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS
598

599
  // Remove lower coordinate levels
600
  n_coord() = 1;
29,911✔
601
  bool found = exhaustive_find_cell(*this, verbose);
29,911✔
602

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

609
    surface() = SURFACE_NONE;
5,799✔
610
    n_coord() = 1;
5,799✔
611
    r() += TINY_BIT * u();
5,799✔
612

613
    // Couldn't find next cell anywhere! This probably means there is an actual
614
    // undefined region in the geometry.
615

616
    if (!exhaustive_find_cell(*this, verbose)) {
5,799!
617
      mark_as_lost("After particle " + std::to_string(id()) +
17,388✔
618
                   " crossed surface " + std::to_string(surf.id_) +
23,178✔
619
                   " it could not be located in any cell and it did not leak.");
620
      return;
5,790✔
621
    }
622
  }
623
}
624

625
void Particle::cross_vacuum_bc(const Surface& surf)
37,168,105✔
626
{
627
  // Score any surface current tallies -- note that the particle is moved
628
  // forward slightly so that if the mesh boundary is on the surface, it is
629
  // still processed
630

631
  if (!model::active_meshsurf_tallies.empty()) {
37,168,105✔
632
    // TODO: Find a better solution to score surface currents than
633
    // physically moving the particle forward slightly
634

635
    r() += TINY_BIT * u();
1,237,190✔
636
    score_surface_tally(*this, model::active_meshsurf_tallies);
1,237,190✔
637
  }
638

639
  // Score to global leakage tally
640
  keff_tally_leakage() += wgt();
37,168,105✔
641

642
  // Kill the particle
643
  wgt() = 0.0;
37,168,105✔
644

645
  // Display message
646
  if (settings::verbosity >= 10 || trace()) {
37,168,105!
647
    write_message(1, "    Leaked out of surface {}", surf.id_);
11✔
648
  }
649
}
37,168,105✔
650

651
void Particle::cross_reflective_bc(const Surface& surf, Direction new_u)
671,972,127✔
652
{
653
  // Do not handle reflective boundary conditions on lower universes
654
  if (n_coord() != 1) {
671,972,127!
UNCOV
655
    mark_as_lost("Cannot reflect particle " + std::to_string(id()) +
×
656
                 " off surface in a lower universe.");
UNCOV
657
    return;
×
658
  }
659

660
  // Score surface currents since reflection causes the direction of the
661
  // particle to change. For surface filters, we need to score the tallies
662
  // twice, once before the particle's surface attribute has changed and
663
  // once after. For mesh surface filters, we need to artificially move
664
  // the particle slightly back in case the surface crossing is coincident
665
  // with a mesh boundary
666

667
  if (!model::active_surface_tallies.empty()) {
671,972,127✔
668
    score_surface_tally(*this, model::active_surface_tallies);
285,021✔
669
  }
670

671
  if (!model::active_meshsurf_tallies.empty()) {
671,972,127✔
672
    Position r {this->r()};
63,726,503✔
673
    this->r() -= TINY_BIT * u();
63,726,503✔
674
    score_surface_tally(*this, model::active_meshsurf_tallies);
63,726,503✔
675
    this->r() = r;
63,726,503✔
676
  }
677

678
  // Set the new particle direction
679
  u() = new_u;
671,972,127✔
680

681
  // Reassign particle's cell and surface
682
  coord(0).cell() = cell_last(0);
671,972,127✔
683
  surface() = -surface();
671,972,127✔
684

685
  // If a reflective surface is coincident with a lattice or universe
686
  // boundary, it is necessary to redetermine the particle's coordinates in
687
  // the lower universes.
688
  // (unless we're using a dagmc model, which has exactly one universe)
689
  n_coord() = 1;
671,972,127✔
690
  if (surf.geom_type() != GeometryType::DAG &&
1,343,941,496!
691
      !neighbor_list_find_cell(*this)) {
671,969,369!
UNCOV
692
    mark_as_lost("Couldn't find particle after reflecting from surface " +
×
UNCOV
693
                 std::to_string(surf.id_) + ".");
×
UNCOV
694
    return;
×
695
  }
696

697
  // Set previous coordinate going slightly past surface crossing
698
  r_last_current() = r() + TINY_BIT * u();
671,972,127✔
699

700
  // Diagnostic message
701
  if (settings::verbosity >= 10 || trace()) {
671,972,127!
UNCOV
702
    write_message(1, "    Reflected from surface {}", surf.id_);
×
703
  }
704
}
705

706
void Particle::cross_periodic_bc(
661,623✔
707
  const Surface& surf, Position new_r, Direction new_u, int new_surface)
708
{
709
  // Do not handle periodic boundary conditions on lower universes
710
  if (n_coord() != 1) {
661,623!
UNCOV
711
    mark_as_lost(
×
712
      "Cannot transfer particle " + std::to_string(id()) +
×
713
      " across surface in a lower universe. Boundary conditions must be "
714
      "applied to root universe.");
UNCOV
715
    return;
×
716
  }
717

718
  // Score surface currents since reflection causes the direction of the
719
  // particle to change -- artificially move the particle slightly back in
720
  // case the surface crossing is coincident with a mesh boundary
721
  if (!model::active_meshsurf_tallies.empty()) {
661,623!
722
    Position r {this->r()};
×
UNCOV
723
    this->r() -= TINY_BIT * u();
×
UNCOV
724
    score_surface_tally(*this, model::active_meshsurf_tallies);
×
UNCOV
725
    this->r() = r;
×
726
  }
727

728
  // Adjust the particle's location and direction.
729
  r() = new_r;
661,623✔
730
  u() = new_u;
661,623✔
731

732
  // Reassign particle's surface
733
  surface() = new_surface;
661,623✔
734

735
  // Figure out what cell particle is in now
736
  n_coord() = 1;
661,623✔
737

738
  if (!neighbor_list_find_cell(*this)) {
661,623!
UNCOV
739
    mark_as_lost("Couldn't find particle after hitting periodic "
×
UNCOV
740
                 "boundary on surface " +
×
741
                 std::to_string(surf.id_) +
×
742
                 ". The normal vector "
743
                 "of one periodic surface may need to be reversed.");
UNCOV
744
    return;
×
745
  }
746

747
  // Set previous coordinate going slightly past surface crossing
748
  r_last_current() = r() + TINY_BIT * u();
661,623✔
749

750
  // Diagnostic message
751
  if (settings::verbosity >= 10 || trace()) {
661,623!
UNCOV
752
    write_message(1, "    Hit periodic boundary on surface {}", surf.id_);
×
753
  }
754
}
755

756
void Particle::mark_as_lost(const char* message)
5,799✔
757
{
758
  // Print warning and write lost particle file
759
  warning(message);
5,799✔
760
  if (settings::max_write_lost_particles < 0 ||
5,799✔
761
      simulation::n_lost_particles < settings::max_write_lost_particles) {
5,500✔
762
    write_restart();
379✔
763
  }
764
  // Increment number of lost particles
765
  wgt() = 0.0;
5,799✔
766
#pragma omp atomic
3,154✔
767
  simulation::n_lost_particles += 1;
2,645✔
768

769
  // Count the total number of simulated particles (on this processor)
770
  auto n = simulation::current_batch * settings::gen_per_batch *
5,799✔
771
           simulation::work_per_rank;
772

773
  // Abort the simulation if the maximum number of lost particles has been
774
  // reached
775
  if (simulation::n_lost_particles >= settings::max_lost_particles &&
5,799✔
776
      simulation::n_lost_particles >= settings::rel_max_lost_particles * n) {
9!
777
    fatal_error("Maximum number of lost particles has been reached.");
9✔
778
  }
779
}
5,790✔
780

781
void Particle::write_restart() const
379✔
782
{
783
  // Dont write another restart file if in particle restart mode
784
  if (settings::run_mode == RunMode::PARTICLE)
379✔
785
    return;
22✔
786

787
  // Set up file name
788
  auto filename = fmt::format("{}particle_{}_{}.h5", settings::path_output,
789
    simulation::current_batch, id());
665✔
790

791
#pragma omp critical(WriteParticleRestart)
374✔
792
  {
793
    // Create file
794
    hid_t file_id = file_open(filename, 'w');
357✔
795

796
    // Write filetype and version info
797
    write_attribute(file_id, "filetype", "particle restart");
357✔
798
    write_attribute(file_id, "version", VERSION_PARTICLE_RESTART);
357✔
799
    write_attribute(file_id, "openmc_version", VERSION);
357✔
800
#ifdef GIT_SHA1
801
    write_attr_string(file_id, "git_sha1", GIT_SHA1);
802
#endif
803

804
    // Write data to file
805
    write_dataset(file_id, "current_batch", simulation::current_batch);
357✔
806
    write_dataset(file_id, "generations_per_batch", settings::gen_per_batch);
357✔
807
    write_dataset(file_id, "current_generation", simulation::current_gen);
357✔
808
    write_dataset(file_id, "n_particles", settings::n_particles);
357✔
809
    switch (settings::run_mode) {
357!
810
    case RunMode::FIXED_SOURCE:
225✔
811
      write_dataset(file_id, "run_mode", "fixed source");
225✔
812
      break;
225✔
813
    case RunMode::EIGENVALUE:
132✔
814
      write_dataset(file_id, "run_mode", "eigenvalue");
132✔
815
      break;
132✔
816
    case RunMode::PARTICLE:
×
817
      write_dataset(file_id, "run_mode", "particle restart");
×
UNCOV
818
      break;
×
UNCOV
819
    default:
×
UNCOV
820
      break;
×
821
    }
822
    write_dataset(file_id, "id", id());
357✔
823
    write_dataset(file_id, "type", static_cast<int>(type()));
357✔
824

825
    int64_t i = current_work();
357✔
826
    if (settings::run_mode == RunMode::EIGENVALUE) {
357✔
827
      // take source data from primary bank for eigenvalue simulation
828
      write_dataset(file_id, "weight", simulation::source_bank[i - 1].wgt);
132✔
829
      write_dataset(file_id, "energy", simulation::source_bank[i - 1].E);
132✔
830
      write_dataset(file_id, "xyz", simulation::source_bank[i - 1].r);
132✔
831
      write_dataset(file_id, "uvw", simulation::source_bank[i - 1].u);
132✔
832
      write_dataset(file_id, "time", simulation::source_bank[i - 1].time);
132✔
833
    } else if (settings::run_mode == RunMode::FIXED_SOURCE) {
225!
834
      // re-sample using rng random number seed used to generate source particle
835
      int64_t id = (simulation::total_gen + overall_generation() - 1) *
225✔
836
                     settings::n_particles +
225✔
837
                   simulation::work_index[mpi::rank] + i;
225✔
838
      uint64_t seed = init_seed(id, STREAM_SOURCE);
225✔
839
      // re-sample source site
840
      auto site = sample_external_source(&seed);
225✔
841
      write_dataset(file_id, "weight", site.wgt);
225✔
842
      write_dataset(file_id, "energy", site.E);
225✔
843
      write_dataset(file_id, "xyz", site.r);
225✔
844
      write_dataset(file_id, "uvw", site.u);
225✔
845
      write_dataset(file_id, "time", site.time);
225✔
846
    }
847

848
    // Close file
849
    file_close(file_id);
357✔
850
  } // #pragma omp critical
851
}
357✔
852

853
void Particle::update_neutron_xs(
2,147,483,647✔
854
  int i_nuclide, int i_grid, int i_sab, double sab_frac, double ncrystal_xs)
855
{
856
  // Get microscopic cross section cache
857
  auto& micro = this->neutron_xs(i_nuclide);
2,147,483,647✔
858

859
  // If the cache doesn't match, recalculate micro xs
860
  if (this->E() != micro.last_E || this->sqrtkT() != micro.last_sqrtkT ||
2,147,483,647✔
861
      i_sab != micro.index_sab || sab_frac != micro.sab_frac ||
2,147,483,647✔
862
      ncrystal_xs != micro.ncrystal_xs) {
2,147,483,647!
863
    data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, *this);
2,147,483,647✔
864

865
    // If NCrystal is being used, update micro cross section cache
866
    micro.ncrystal_xs = ncrystal_xs;
2,147,483,647✔
867
    if (ncrystal_xs >= 0.0) {
2,147,483,647✔
868
      data::nuclides[i_nuclide]->calculate_elastic_xs(*this);
11,018,953✔
869
      ncrystal_update_micro(ncrystal_xs, micro);
11,018,953✔
870
    }
871
  }
872
}
2,147,483,647✔
873

874
//==============================================================================
875
// Non-method functions
876
//==============================================================================
877

878
std::string particle_type_to_str(ParticleType type)
3,644,489✔
879
{
880
  switch (type) {
3,644,489!
881
  case ParticleType::neutron:
2,657,092✔
882
    return "neutron";
2,657,092✔
883
  case ParticleType::photon:
987,133✔
884
    return "photon";
987,133✔
885
  case ParticleType::electron:
132✔
886
    return "electron";
132✔
887
  case ParticleType::positron:
132✔
888
    return "positron";
132✔
889
  }
UNCOV
890
  UNREACHABLE();
×
891
}
892

893
ParticleType str_to_particle_type(std::string str)
3,617,092✔
894
{
895
  if (str == "neutron") {
3,617,092✔
896
    return ParticleType::neutron;
826,657✔
897
  } else if (str == "photon") {
2,790,435✔
898
    return ParticleType::photon;
2,790,349✔
899
  } else if (str == "electron") {
86✔
900
    return ParticleType::electron;
43✔
901
  } else if (str == "positron") {
43!
902
    return ParticleType::positron;
43✔
903
  } else {
UNCOV
904
    throw std::invalid_argument {fmt::format("Invalid particle name: {}", str)};
×
905
  }
906
}
907

908
void add_surf_source_to_bank(Particle& p, const Surface& surf)
1,548,051,689✔
909
{
910
  if (simulation::current_batch <= settings::n_inactive ||
2,147,483,647✔
911
      simulation::surf_source_bank.full()) {
1,226,115,598✔
912
    return;
1,547,922,035✔
913
  }
914

915
  // If a cell/cellfrom/cellto parameter is defined
916
  if (settings::ssw_cell_id != C_NONE) {
337,086✔
917

918
    // Retrieve cell index and storage type
919
    int cell_idx = model::cell_map[settings::ssw_cell_id];
254,440✔
920

921
    if (surf.bc_) {
254,440✔
922
      // Leave if cellto with vacuum boundary condition
923
      if (surf.bc_->type() == "vacuum" &&
182,558!
924
          settings::ssw_cell_type == SSWCellType::To) {
33,099✔
925
        return;
12,136✔
926
      }
927

928
      // Leave if other boundary condition than vacuum
929
      if (surf.bc_->type() != "vacuum") {
137,323✔
930
        return;
116,360✔
931
      }
932
    }
933

934
    // Check if the cell of interest has been exited
935
    bool exited = false;
125,944✔
936
    for (int i = 0; i < p.n_coord_last(); ++i) {
333,677✔
937
      if (p.cell_last(i) == cell_idx) {
207,733✔
938
        exited = true;
73,765✔
939
      }
940
    }
941

942
    // Check if the cell of interest has been entered
943
    bool entered = false;
125,944✔
944
    for (int i = 0; i < p.n_coord(); ++i) {
297,979✔
945
      if (p.coord(i).cell() == cell_idx) {
172,035✔
946
        entered = true;
57,517✔
947
      }
948
    }
949

950
    // Vacuum boundary conditions: return if cell is not exited
951
    if (surf.bc_) {
125,944✔
952
      if (surf.bc_->type() == "vacuum" && !exited) {
20,963!
953
        return;
14,663✔
954
      }
955
    } else {
956

957
      // If we both enter and exit the cell of interest
958
      if (entered && exited) {
104,981✔
959
        return;
27,203✔
960
      }
961

962
      // If we did not enter nor exit the cell of interest
963
      if (!entered && !exited) {
77,778✔
964
        return;
13,502✔
965
      }
966

967
      // If cellfrom and the cell before crossing is not the cell of
968
      // interest
969
      if (settings::ssw_cell_type == SSWCellType::From && !exited) {
64,276✔
970
        return;
11,543✔
971
      }
972

973
      // If cellto and the cell after crossing is not the cell of interest
974
      if (settings::ssw_cell_type == SSWCellType::To && !entered) {
52,733✔
975
        return;
12,025✔
976
      }
977
    }
978
  }
979

980
  SourceSite site;
129,654✔
981
  site.r = p.r();
129,654✔
982
  site.u = p.u();
129,654✔
983
  site.E = p.E();
129,654✔
984
  site.time = p.time();
129,654✔
985
  site.wgt = p.wgt();
129,654✔
986
  site.delayed_group = p.delayed_group();
129,654✔
987
  site.surf_id = surf.id_;
129,654✔
988
  site.particle = p.type();
129,654✔
989
  site.parent_id = p.id();
129,654✔
990
  site.progeny_id = p.n_progeny();
129,654✔
991
  int64_t idx = simulation::surf_source_bank.thread_safe_append(site);
129,654✔
992
}
993

994
} // namespace openmc
STATUS · Troubleshooting · Open an Issue · Sales · Support · CAREERS · ENTERPRISE · START FREE · SCHEDULE DEMO
ANNOUNCEMENTS · TWITTER · TOS & SLA · Supported CI Services · What's a CI service? · Automated Testing

© 2025 Coveralls, Inc