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

openmc-dev / openmc / 21691294911

04 Feb 2026 10:44PM UTC coverage: 81.754% (-0.009%) from 81.763%
21691294911

Pull #3765

github

web-flow
Merge a50f9a6c7 into a22238e06
Pull Request #3765: Store atomic mass in ParticleType.

17301 of 24254 branches covered (71.33%)

Branch coverage included in aggregate %.

5 of 6 new or added lines in 2 files covered. (83.33%)

28 existing lines in 1 file now uncovered.

55897 of 65280 relevant lines covered (85.63%)

44329169.07 hits per line

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

85.34
/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 = this->type().mass() * AMU_EV;
2,114,290,008✔
51
    // Equivalent to C * sqrt(1-(m/(m+E))^2) without problem at E<<m:
52
    return C_LIGHT * std::sqrt(this->E() * (this->E() + 2 * mass)) /
2,114,290,008✔
53
           (this->E() + mass);
2,114,290,008✔
54
  } else {
55
    auto& macro_xs = data::mg.macro_xs_[this->material()];
2,063,937,414✔
56
    int macro_t = this->mg_xs_cache().t;
2,063,937,414✔
57
    int macro_a = macro_xs.get_angle_index(this->u());
2,063,937,414✔
58
    return 1.0 / macro_xs.get_xs(MgxsType::INVERSE_VELOCITY, this->g(), nullptr,
2,063,937,414✔
59
                   nullptr, nullptr, macro_t, macro_a);
2,063,937,414✔
60
  }
61
}
62

63
bool Particle::create_secondary(
111,458,605✔
64
  double wgt, Direction u, double E, ParticleType type)
65
{
66
  // If energy is below cutoff for this particle, don't create secondary
67
  // particle
68
  int idx = type.transport_index();
111,458,605✔
69
  if (idx == C_NONE) {
111,458,605!
UNCOV
70
    return false;
×
71
  }
72
  if (E < settings::energy_cutoff[idx]) {
111,458,605✔
73
    return false;
53,276,669✔
74
  }
75

76
  auto& bank = secondary_bank().emplace_back();
58,181,936✔
77
  bank.particle = type;
58,181,936✔
78
  bank.wgt = wgt;
58,181,936✔
79
  bank.r = r();
58,181,936✔
80
  bank.u = u;
58,181,936✔
81
  bank.E = settings::run_CE ? E : g();
58,181,936!
82
  bank.time = time();
58,181,936✔
83
  bank_second_E() += bank.E;
58,181,936✔
84
  return true;
58,181,936✔
85
}
86

87
void Particle::split(double wgt)
4,084,524✔
88
{
89
  auto& bank = secondary_bank().emplace_back();
4,084,524✔
90
  bank.particle = type();
4,084,524✔
91
  bank.wgt = wgt;
4,084,524✔
92
  bank.r = r();
4,084,524✔
93
  bank.u = u();
4,084,524✔
94
  bank.E = settings::run_CE ? E() : g();
4,084,524✔
95
  bank.time = time();
4,084,524✔
96

97
  // Convert signed index to a signed surface ID
98
  if (surface() == SURFACE_NONE) {
4,084,524✔
99
    bank.surf_id = SURFACE_NONE;
4,084,188✔
100
  } else {
101
    int surf_id = model::surfaces[surface_index()]->id_;
336✔
102
    bank.surf_id = (surface() > 0) ? surf_id : -surf_id;
336!
103
  }
104
}
4,084,524✔
105

106
void Particle::from_source(const SourceSite* src)
238,030,373✔
107
{
108
  // Reset some attributes
109
  clear();
238,030,373✔
110
  surface() = SURFACE_NONE;
238,030,373✔
111
  cell_born() = C_NONE;
238,030,373✔
112
  material() = C_NONE;
238,030,373✔
113
  n_collision() = 0;
238,030,373✔
114
  fission() = false;
238,030,373✔
115
  zero_flux_derivs();
238,030,373✔
116
  lifetime() = 0.0;
238,030,373✔
117
#ifdef OPENMC_DAGMC_ENABLED
118
  history().reset();
21,803,685✔
119
#endif
120

121
  // Copy attributes from source bank site
122
  type() = src->particle;
238,030,373✔
123
  wgt() = src->wgt;
238,030,373✔
124
  wgt_last() = src->wgt;
238,030,373✔
125
  r() = src->r;
238,030,373✔
126
  u() = src->u;
238,030,373✔
127
  r_born() = src->r;
238,030,373✔
128
  r_last_current() = src->r;
238,030,373✔
129
  r_last() = src->r;
238,030,373✔
130
  u_last() = src->u;
238,030,373✔
131
  if (settings::run_CE) {
238,030,373✔
132
    E() = src->E;
122,291,756✔
133
    g() = 0;
122,291,756✔
134
  } else {
135
    g() = static_cast<int>(src->E);
115,738,617✔
136
    g_last() = static_cast<int>(src->E);
115,738,617✔
137
    E() = data::mg.energy_bin_avg_[g()];
115,738,617✔
138
  }
139
  E_last() = E();
238,030,373✔
140
  time() = src->time;
238,030,373✔
141
  time_last() = src->time;
238,030,373✔
142
  parent_nuclide() = src->parent_nuclide;
238,030,373✔
143
  delayed_group() = src->delayed_group;
238,030,373✔
144

145
  // Convert signed surface ID to signed index
146
  if (src->surf_id != SURFACE_NONE) {
238,030,373✔
147
    int index_plus_one = model::surface_map[std::abs(src->surf_id)] + 1;
110,336✔
148
    surface() = (src->surf_id > 0) ? index_plus_one : -index_plus_one;
110,336!
149
  }
150
}
238,030,373✔
151

152
void Particle::event_calculate_xs()
2,147,483,647✔
153
{
154
  // Set the random number stream
155
  stream() = STREAM_TRACKING;
2,147,483,647✔
156

157
  // Store pre-collision particle properties
158
  wgt_last() = wgt();
2,147,483,647✔
159
  E_last() = E();
2,147,483,647✔
160
  u_last() = u();
2,147,483,647✔
161
  r_last() = r();
2,147,483,647✔
162
  time_last() = time();
2,147,483,647✔
163

164
  // Reset event variables
165
  event() = TallyEvent::KILL;
2,147,483,647✔
166
  event_nuclide() = NUCLIDE_NONE;
2,147,483,647✔
167
  event_mt() = REACTION_NONE;
2,147,483,647✔
168

169
  // If the cell hasn't been determined based on the particle's location,
170
  // initiate a search for the current cell. This generally happens at the
171
  // beginning of the history and again for any secondary particles
172
  if (lowest_coord().cell() == C_NONE) {
2,147,483,647✔
173
    if (!exhaustive_find_cell(*this)) {
229,546,464!
UNCOV
174
      mark_as_lost(
×
UNCOV
175
        "Could not find the cell containing particle " + std::to_string(id()));
×
UNCOV
176
      return;
×
177
    }
178

179
    // Set birth cell attribute
180
    if (cell_born() == C_NONE)
229,546,464!
181
      cell_born() = lowest_coord().cell();
229,546,464✔
182

183
    // Initialize last cells from current cell
184
    for (int j = 0; j < n_coord(); ++j) {
476,117,608✔
185
      cell_last(j) = coord(j).cell();
246,571,144✔
186
    }
187
    n_coord_last() = n_coord();
229,546,464✔
188
  }
189

190
  // Write particle track.
191
  if (write_track())
2,147,483,647✔
192
    write_particle_track(*this);
10,818✔
193

194
  if (settings::check_overlaps)
2,147,483,647!
UNCOV
195
    check_cell_overlap(*this);
×
196

197
  // Calculate microscopic and macroscopic cross sections
198
  if (material() != MATERIAL_VOID) {
2,147,483,647✔
199
    if (settings::run_CE) {
2,147,483,647✔
200
      if (material() != material_last() || sqrtkT() != sqrtkT_last() ||
2,147,483,647✔
201
          density_mult() != density_mult_last()) {
365,033,992✔
202
        // If the material is the same as the last material and the
203
        // temperature hasn't changed, we don't need to lookup cross
204
        // sections again.
205
        model::materials[material()]->calculate_xs(*this);
1,609,993,853✔
206
      }
207
    } else {
208
      // Get the MG data; unlike the CE case above, we have to re-calculate
209
      // cross sections for every collision since the cross sections may
210
      // be angle-dependent
211
      data::mg.macro_xs_[material()].calculate_xs(*this);
2,063,937,414✔
212

213
      // Update the particle's group while we know we are multi-group
214
      g_last() = g();
2,063,937,414✔
215
    }
216
  } else {
217
    macro_xs().total = 0.0;
111,825,094✔
218
    macro_xs().absorption = 0.0;
111,825,094✔
219
    macro_xs().fission = 0.0;
111,825,094✔
220
    macro_xs().nu_fission = 0.0;
111,825,094✔
221
  }
222
}
223

224
void Particle::event_advance()
2,147,483,647✔
225
{
226
  // Find the distance to the nearest boundary
227
  boundary() = distance_to_boundary(*this);
2,147,483,647✔
228

229
  // Sample a distance to collision
230
  if (type() == ParticleType::electron() ||
2,147,483,647✔
231
      type() == ParticleType::positron()) {
2,147,483,647✔
232
    collision_distance() = material() == MATERIAL_VOID ? INFINITY : 0.0;
53,430,397!
233
  } else if (macro_xs().total == 0.0) {
2,147,483,647✔
234
    collision_distance() = INFINITY;
111,825,094✔
235
  } else {
236
    collision_distance() = -std::log(prn(current_seed())) / macro_xs().total;
2,147,483,647✔
237
  }
238

239
  double speed = this->speed();
2,147,483,647✔
240
  double time_cutoff = settings::time_cutoff[type().transport_index()];
2,147,483,647✔
241
  double distance_cutoff =
242
    (time_cutoff < INFTY) ? (time_cutoff - time()) * speed : INFTY;
2,147,483,647✔
243

244
  // Select smaller of the three distances
245
  double distance =
246
    std::min({boundary().distance(), collision_distance(), distance_cutoff});
2,147,483,647✔
247

248
  // Advance particle in space and time
249
  this->move_distance(distance);
2,147,483,647✔
250
  double dt = distance / speed;
2,147,483,647✔
251
  this->time() += dt;
2,147,483,647✔
252
  this->lifetime() += dt;
2,147,483,647✔
253

254
  // Score timed track-length tallies
255
  if (!model::active_timed_tracklength_tallies.empty()) {
2,147,483,647✔
256
    score_timed_tracklength_tally(*this, distance);
3,628,317✔
257
  }
258

259
  // Score track-length tallies
260
  if (!model::active_tracklength_tallies.empty()) {
2,147,483,647✔
261
    score_tracklength_tally(*this, distance);
1,562,801,097✔
262
  }
263

264
  // Score track-length estimate of k-eff
265
  if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {
2,147,483,647✔
266
    keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission;
2,147,483,647✔
267
  }
268

269
  // Score flux derivative accumulators for differential tallies.
270
  if (!model::active_tallies.empty()) {
2,147,483,647✔
271
    score_track_derivative(*this, distance);
1,732,243,720✔
272
  }
273

274
  // Set particle weight to zero if it hit the time boundary
275
  if (distance == distance_cutoff) {
2,147,483,647✔
276
    wgt() = 0.0;
224,928✔
277
  }
278
}
2,147,483,647✔
279

280
void Particle::event_cross_surface()
2,147,483,647✔
281
{
282
  // Saving previous cell data
283
  for (int j = 0; j < n_coord(); ++j) {
2,147,483,647✔
284
    cell_last(j) = coord(j).cell();
2,147,483,647✔
285
  }
286
  n_coord_last() = n_coord();
2,147,483,647✔
287

288
  // Set surface that particle is on and adjust coordinate levels
289
  surface() = boundary().surface();
2,147,483,647✔
290
  n_coord() = boundary().coord_level();
2,147,483,647✔
291

292
  if (boundary().lattice_translation()[0] != 0 ||
2,147,483,647✔
293
      boundary().lattice_translation()[1] != 0 ||
2,147,483,647✔
294
      boundary().lattice_translation()[2] != 0) {
1,695,547,466✔
295
    // Particle crosses lattice boundary
296

297
    bool verbose = settings::verbosity >= 10 || trace();
732,435,563!
298
    cross_lattice(*this, boundary(), verbose);
732,435,563✔
299
    event() = TallyEvent::LATTICE;
732,435,563✔
300
  } else {
301
    // Particle crosses surface
302
    const auto& surf {model::surfaces[surface_index()].get()};
1,507,005,571✔
303
    // If BC, add particle to surface source before crossing surface
304
    if (surf->surf_source_ && surf->bc_) {
1,507,005,571✔
305
      add_surf_source_to_bank(*this, *surf);
688,546,037✔
306
    }
307
    this->cross_surface(*surf);
1,507,005,571✔
308
    // If no BC, add particle to surface source after crossing surface
309
    if (surf->surf_source_ && !surf->bc_) {
1,507,005,562✔
310
      add_surf_source_to_bank(*this, *surf);
817,221,698✔
311
    }
312
    if (settings::weight_window_checkpoint_surface) {
1,507,005,562✔
313
      apply_weight_windows(*this);
10,738!
314
    }
315
    event() = TallyEvent::SURFACE;
1,507,005,562✔
316
  }
317
  // Score cell to cell partial currents
318
  if (!model::active_surface_tallies.empty()) {
2,147,483,647✔
319
    score_surface_tally(*this, model::active_surface_tallies);
34,922,767✔
320
  }
321
}
2,147,483,647✔
322

323
void Particle::event_collide()
2,147,483,647✔
324
{
325
  // Score collision estimate of keff
326
  if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) {
2,147,483,647✔
327
    keff_tally_collision() += wgt() * macro_xs().nu_fission / macro_xs().total;
2,125,051,490✔
328
  }
329

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

334
  if (!model::active_meshsurf_tallies.empty())
2,147,483,647✔
335
    score_surface_tally(*this, model::active_meshsurf_tallies);
63,098,926✔
336

337
  // Clear surface component
338
  surface() = SURFACE_NONE;
2,147,483,647✔
339

340
  if (settings::run_CE) {
2,147,483,647✔
341
    collision(*this);
944,050,036✔
342
  } else {
343
    collision_mg(*this);
1,783,060,477✔
344
  }
345

346
  // Collision track feature to recording particle interaction
347
  if (settings::collision_track) {
2,147,483,647✔
348
    collision_track_record(*this);
150,087✔
349
  }
350

351
  // Score collision estimator tallies -- this is done after a collision
352
  // has occurred rather than before because we need information on the
353
  // outgoing energy for any tallies with an outgoing energy filter
354
  if (!model::active_collision_tallies.empty())
2,147,483,647✔
355
    score_collision_tally(*this);
107,002,479✔
356
  if (!model::active_analog_tallies.empty()) {
2,147,483,647✔
357
    if (settings::run_CE) {
291,821,708✔
358
      score_analog_tally_ce(*this);
290,613,446✔
359
    } else {
360
      score_analog_tally_mg(*this);
1,208,262✔
361
    }
362
  }
363

364
  if (!model::active_pulse_height_tallies.empty() && type().is_photon()) {
2,147,483,647✔
365
    pht_collision_energy();
2,024✔
366
  }
367

368
  // Reset banked weight during collision
369
  n_bank() = 0;
2,147,483,647✔
370
  bank_second_E() = 0.0;
2,147,483,647✔
371
  wgt_bank() = 0.0;
2,147,483,647✔
372
  zero_delayed_bank();
2,147,483,647✔
373

374
  // Reset fission logical
375
  fission() = false;
2,147,483,647✔
376

377
  // Save coordinates for tallying purposes
378
  r_last_current() = r();
2,147,483,647✔
379

380
  // Set last material to none since cross sections will need to be
381
  // re-evaluated
382
  material_last() = C_NONE;
2,147,483,647✔
383

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

398
  // Score flux derivative accumulators for differential tallies.
399
  if (!model::active_tallies.empty())
2,147,483,647✔
400
    score_collision_derivative(*this);
816,414,278✔
401

402
#ifdef OPENMC_DAGMC_ENABLED
403
  history().reset();
250,231,364✔
404
#endif
405
}
2,147,483,647✔
406

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

417
  // Check for secondary particles if this particle is dead
418
  if (!alive()) {
2,147,483,647✔
419
    // Write final position for this particle
420
    if (write_track()) {
229,546,060✔
421
      write_particle_track(*this);
6,676✔
422
    }
423

424
    // If no secondary particles, break out of event loop
425
    if (secondary_bank().empty())
229,546,060✔
426
      return;
167,076,081✔
427

428
    from_source(&secondary_bank().back());
62,469,979✔
429
    secondary_bank().pop_back();
62,469,979✔
430
    n_event() = 0;
62,469,979✔
431
    bank_second_E() = 0.0;
62,469,979✔
432

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

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

459
    // Enter new particle in particle track file
460
    if (write_track())
62,469,979✔
461
      add_particle_track(*this);
5,606✔
462
  }
463
}
464

465
void Particle::event_death()
167,077,081✔
466
{
467
#ifdef OPENMC_DAGMC_ENABLED
468
  history().reset();
15,264,135✔
469
#endif
470

471
  // Finish particle track output.
472
  if (write_track()) {
167,077,081✔
473
    finalize_particle_track(*this);
1,070✔
474
  }
475

476
// Contribute tally reduction variables to global accumulator
477
#pragma omp atomic
91,683,199✔
478
  global_tally_absorption += keff_tally_absorption();
167,077,081✔
479
#pragma omp atomic
92,125,625✔
480
  global_tally_collision += keff_tally_collision();
167,077,081✔
481
#pragma omp atomic
91,450,975✔
482
  global_tally_tracklength += keff_tally_tracklength();
167,077,081✔
483
#pragma omp atomic
91,345,669✔
484
  global_tally_leakage += keff_tally_leakage();
167,077,081✔
485

486
  // Reset particle tallies once accumulated
487
  keff_tally_absorption() = 0.0;
167,077,081✔
488
  keff_tally_collision() = 0.0;
167,077,081✔
489
  keff_tally_tracklength() = 0.0;
167,077,081✔
490
  keff_tally_leakage() = 0.0;
167,077,081✔
491

492
  if (!model::active_pulse_height_tallies.empty()) {
167,077,081✔
493
    score_pulse_height_tally(*this, model::active_pulse_height_tallies);
5,500✔
494
  }
495

496
  // Record the number of progeny created by this particle.
497
  // This data will be used to efficiently sort the fission bank.
498
  if (settings::run_mode == RunMode::EIGENVALUE) {
167,077,081✔
499
    int64_t offset = id() - 1 - simulation::work_index[mpi::rank];
140,925,700✔
500
    simulation::progeny_per_particle[offset] = n_progeny();
140,925,700✔
501
  }
502
}
167,077,081✔
503

504
void Particle::pht_collision_energy()
2,024✔
505
{
506
  // Adds the energy particles lose in a collision to the pulse-height
507

508
  // determine index of cell in pulse_height_cells
509
  auto it = std::find(model::pulse_height_cells.begin(),
2,024✔
510
    model::pulse_height_cells.end(), lowest_coord().cell());
2,024✔
511

512
  if (it != model::pulse_height_cells.end()) {
2,024!
513
    int index = std::distance(model::pulse_height_cells.begin(), it);
2,024✔
514
    pht_storage()[index] += E_last() - E();
2,024✔
515

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

525
void Particle::pht_secondary_particles()
605✔
526
{
527
  // Removes the energy of secondary produced particles from the pulse-height
528

529
  // determine index of cell in pulse_height_cells
530
  auto it = std::find(model::pulse_height_cells.begin(),
605✔
531
    model::pulse_height_cells.end(), cell_born());
605✔
532

533
  if (it != model::pulse_height_cells.end()) {
605!
534
    int index = std::distance(model::pulse_height_cells.begin(), it);
605✔
535
    pht_storage()[index] -= E();
605✔
536
  }
537
}
605✔
538

539
void Particle::cross_surface(const Surface& surf)
1,508,874,989✔
540
{
541

542
  if (settings::verbosity >= 10 || trace()) {
1,508,874,989✔
543
    write_message(1, "    Crossing surface {}", surf.id_);
33✔
544
  }
545

546
// if we're crossing a CSG surface, make sure the DAG history is reset
547
#ifdef OPENMC_DAGMC_ENABLED
548
  if (surf.geom_type() == GeometryType::CSG)
137,898,723✔
549
    history().reset();
137,843,604✔
550
#endif
551

552
  // Handle any applicable boundary conditions.
553
  if (surf.bc_ && settings::run_mode != RunMode::PLOTTING &&
2,147,483,647!
554
      settings::run_mode != RunMode::VOLUME) {
689,018,089✔
555
    surf.bc_->handle_particle(*this, surf);
688,898,145✔
556
    return;
688,898,145✔
557
  }
558

559
  // ==========================================================================
560
  // SEARCH NEIGHBOR LISTS FOR NEXT CELL
561

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

576
    cell_instance() = 0;
44,310✔
577
    if (cell->distribcell_index_ >= 0)
44,310✔
578
      cell_instance() = cell_instance_at_level(*this, n_coord() - 1);
43,286✔
579

580
    material() = cell->material(cell_instance());
44,310✔
581
    sqrtkT() = cell->sqrtkT(cell_instance());
44,310✔
582
    density_mult() = cell->density_mult(cell_instance());
44,310✔
583
    return;
44,310✔
584
  }
585
#endif
586

587
  bool verbose = settings::verbosity >= 10 || trace();
819,932,534!
588
  if (neighbor_list_find_cell(*this, verbose)) {
819,932,534✔
589
    return;
819,902,623✔
590
  }
591

592
  // ==========================================================================
593
  // COULDN'T FIND PARTICLE IN NEIGHBORING CELLS, SEARCH ALL CELLS
594

595
  // Remove lower coordinate levels
596
  n_coord() = 1;
29,911✔
597
  bool found = exhaustive_find_cell(*this, verbose);
29,911✔
598

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

605
    surface() = SURFACE_NONE;
5,799✔
606
    n_coord() = 1;
5,799✔
607
    r() += TINY_BIT * u();
5,799✔
608

609
    // Couldn't find next cell anywhere! This probably means there is an actual
610
    // undefined region in the geometry.
611

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

621
void Particle::cross_vacuum_bc(const Surface& surf)
35,018,314✔
622
{
623
  // Score any surface current tallies -- note that the particle is moved
624
  // forward slightly so that if the mesh boundary is on the surface, it is
625
  // still processed
626

627
  if (!model::active_meshsurf_tallies.empty()) {
35,018,314✔
628
    // TODO: Find a better solution to score surface currents than
629
    // physically moving the particle forward slightly
630

631
    r() += TINY_BIT * u();
937,222✔
632
    score_surface_tally(*this, model::active_meshsurf_tallies);
937,222✔
633
  }
634

635
  // Score to global leakage tally
636
  keff_tally_leakage() += wgt();
35,018,314✔
637

638
  // Kill the particle
639
  wgt() = 0.0;
35,018,314✔
640

641
  // Display message
642
  if (settings::verbosity >= 10 || trace()) {
35,018,314!
643
    write_message(1, "    Leaked out of surface {}", surf.id_);
11✔
644
  }
645
}
35,018,314✔
646

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

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

663
  if (!model::active_surface_tallies.empty()) {
652,640,003✔
664
    score_surface_tally(*this, model::active_surface_tallies);
285,021✔
665
  }
666

667
  if (!model::active_meshsurf_tallies.empty()) {
652,640,003✔
668
    Position r {this->r()};
46,885,487✔
669
    this->r() -= TINY_BIT * u();
46,885,487✔
670
    score_surface_tally(*this, model::active_meshsurf_tallies);
46,885,487✔
671
    this->r() = r;
46,885,487✔
672
  }
673

674
  // Set the new particle direction
675
  u() = new_u;
652,640,003✔
676

677
  // Reassign particle's cell and surface
678
  coord(0).cell() = cell_last(0);
652,640,003✔
679
  surface() = -surface();
652,640,003✔
680

681
  // If a reflective surface is coincident with a lattice or universe
682
  // boundary, it is necessary to redetermine the particle's coordinates in
683
  // the lower universes.
684
  // (unless we're using a dagmc model, which has exactly one universe)
685
  n_coord() = 1;
652,640,003✔
686
  if (surf.geom_type() != GeometryType::DAG &&
1,305,277,248!
687
      !neighbor_list_find_cell(*this)) {
652,637,245!
UNCOV
688
    mark_as_lost("Couldn't find particle after reflecting from surface " +
×
UNCOV
689
                 std::to_string(surf.id_) + ".");
×
UNCOV
690
    return;
×
691
  }
692

693
  // Set previous coordinate going slightly past surface crossing
694
  r_last_current() = r() + TINY_BIT * u();
652,640,003✔
695

696
  // Diagnostic message
697
  if (settings::verbosity >= 10 || trace()) {
652,640,003!
698
    write_message(1, "    Reflected from surface {}", surf.id_);
×
699
  }
700
}
701

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

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

724
  // Adjust the particle's location and direction.
725
  r() = new_r;
2,245,294✔
726
  u() = new_u;
2,245,294✔
727

728
  // Reassign particle's surface
729
  surface() = new_surface;
2,245,294✔
730

731
  // Figure out what cell particle is in now
732
  n_coord() = 1;
2,245,294✔
733

734
  if (!neighbor_list_find_cell(*this)) {
2,245,294!
UNCOV
735
    mark_as_lost("Couldn't find particle after hitting periodic "
×
UNCOV
736
                 "boundary on surface " +
×
UNCOV
737
                 std::to_string(surf.id_) + ".");
×
UNCOV
738
    return;
×
739
  }
740

741
  // Set previous coordinate going slightly past surface crossing
742
  r_last_current() = r() + TINY_BIT * u();
2,245,294✔
743

744
  // Diagnostic message
745
  if (settings::verbosity >= 10 || trace()) {
2,245,294!
746
    write_message(1, "    Hit periodic boundary on surface {}", surf.id_);
×
747
  }
748
}
749

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

763
  // Count the total number of simulated particles (on this processor)
764
  auto n = simulation::current_batch * settings::gen_per_batch *
5,799✔
765
           simulation::work_per_rank;
766

767
  // Abort the simulation if the maximum number of lost particles has been
768
  // reached
769
  if (simulation::n_lost_particles >= settings::max_lost_particles &&
5,799✔
770
      simulation::n_lost_particles >= settings::rel_max_lost_particles * n) {
9!
771
    fatal_error("Maximum number of lost particles has been reached.");
9✔
772
  }
773
}
5,790✔
774

775
void Particle::write_restart() const
379✔
776
{
777
  // Dont write another restart file if in particle restart mode
778
  if (settings::run_mode == RunMode::PARTICLE)
379✔
779
    return;
22✔
780

781
  // Set up file name
782
  auto filename = fmt::format("{}particle_{}_{}.h5", settings::path_output,
783
    simulation::current_batch, id());
665✔
784

785
#pragma omp critical(WriteParticleRestart)
374✔
786
  {
787
    // Create file
788
    hid_t file_id = file_open(filename, 'w');
357✔
789

790
    // Write filetype and version info
791
    write_attribute(file_id, "filetype", "particle restart");
357✔
792
    write_attribute(file_id, "version", VERSION_PARTICLE_RESTART);
357✔
793
    write_attribute(file_id, "openmc_version", VERSION);
357✔
794
#ifdef GIT_SHA1
795
    write_attr_string(file_id, "git_sha1", GIT_SHA1);
796
#endif
797

798
    // Write data to file
799
    write_dataset(file_id, "current_batch", simulation::current_batch);
357✔
800
    write_dataset(file_id, "generations_per_batch", settings::gen_per_batch);
357✔
801
    write_dataset(file_id, "current_generation", simulation::current_gen);
357✔
802
    write_dataset(file_id, "n_particles", settings::n_particles);
357✔
803
    switch (settings::run_mode) {
357!
804
    case RunMode::FIXED_SOURCE:
225✔
805
      write_dataset(file_id, "run_mode", "fixed source");
225✔
806
      break;
225✔
807
    case RunMode::EIGENVALUE:
132✔
808
      write_dataset(file_id, "run_mode", "eigenvalue");
132✔
809
      break;
132✔
UNCOV
810
    case RunMode::PARTICLE:
×
UNCOV
811
      write_dataset(file_id, "run_mode", "particle restart");
×
UNCOV
812
      break;
×
UNCOV
813
    default:
×
UNCOV
814
      break;
×
815
    }
816
    write_dataset(file_id, "id", id());
357✔
817
    write_dataset(file_id, "type", type().pdg_number());
357✔
818

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

842
    // Close file
843
    file_close(file_id);
357✔
844
  } // #pragma omp critical
845
}
357✔
846

847
void Particle::update_neutron_xs(
2,147,483,647✔
848
  int i_nuclide, int i_grid, int i_sab, double sab_frac, double ncrystal_xs)
849
{
850
  // Get microscopic cross section cache
851
  auto& micro = this->neutron_xs(i_nuclide);
2,147,483,647✔
852

853
  // If the cache doesn't match, recalculate micro xs
854
  if (this->E() != micro.last_E || this->sqrtkT() != micro.last_sqrtkT ||
2,147,483,647✔
855
      i_sab != micro.index_sab || sab_frac != micro.sab_frac ||
2,147,483,647✔
856
      ncrystal_xs != micro.ncrystal_xs) {
2,147,483,647!
857
    data::nuclides[i_nuclide]->calculate_xs(i_sab, i_grid, sab_frac, *this);
2,147,483,647✔
858

859
    // If NCrystal is being used, update micro cross section cache
860
    micro.ncrystal_xs = ncrystal_xs;
2,147,483,647✔
861
    if (ncrystal_xs >= 0.0) {
2,147,483,647✔
862
      data::nuclides[i_nuclide]->calculate_elastic_xs(*this);
11,018,953✔
863
      ncrystal_update_micro(ncrystal_xs, micro);
11,018,953✔
864
    }
865
  }
866
}
2,147,483,647✔
867

868
//==============================================================================
869
// Non-method functions
870
//==============================================================================
871
void add_surf_source_to_bank(Particle& p, const Surface& surf)
1,505,767,735✔
872
{
873
  if (simulation::current_batch <= settings::n_inactive ||
2,147,483,647✔
874
      simulation::surf_source_bank.full()) {
1,178,696,153✔
875
    return;
1,505,638,082✔
876
  }
877

878
  // If a cell/cellfrom/cellto parameter is defined
879
  if (settings::ssw_cell_id != C_NONE) {
337,083✔
880

881
    // Retrieve cell index and storage type
882
    int cell_idx = model::cell_map[settings::ssw_cell_id];
254,438✔
883

884
    if (surf.bc_) {
254,438✔
885
      // Leave if cellto with vacuum boundary condition
886
      if (surf.bc_->type() == "vacuum" &&
182,558!
887
          settings::ssw_cell_type == SSWCellType::To) {
33,099✔
888
        return;
12,136✔
889
      }
890

891
      // Leave if other boundary condition than vacuum
892
      if (surf.bc_->type() != "vacuum") {
137,323✔
893
        return;
116,360✔
894
      }
895
    }
896

897
    // Check if the cell of interest has been exited
898
    bool exited = false;
125,942✔
899
    for (int i = 0; i < p.n_coord_last(); ++i) {
333,673✔
900
      if (p.cell_last(i) == cell_idx) {
207,731✔
901
        exited = true;
73,764✔
902
      }
903
    }
904

905
    // Check if the cell of interest has been entered
906
    bool entered = false;
125,942✔
907
    for (int i = 0; i < p.n_coord(); ++i) {
297,975✔
908
      if (p.coord(i).cell() == cell_idx) {
172,033✔
909
        entered = true;
57,517✔
910
      }
911
    }
912

913
    // Vacuum boundary conditions: return if cell is not exited
914
    if (surf.bc_) {
125,942✔
915
      if (surf.bc_->type() == "vacuum" && !exited) {
20,963!
916
        return;
14,663✔
917
      }
918
    } else {
919

920
      // If we both enter and exit the cell of interest
921
      if (entered && exited) {
104,979✔
922
        return;
27,203✔
923
      }
924

925
      // If we did not enter nor exit the cell of interest
926
      if (!entered && !exited) {
77,776✔
927
        return;
13,501✔
928
      }
929

930
      // If cellfrom and the cell before crossing is not the cell of
931
      // interest
932
      if (settings::ssw_cell_type == SSWCellType::From && !exited) {
64,275✔
933
        return;
11,543✔
934
      }
935

936
      // If cellto and the cell after crossing is not the cell of interest
937
      if (settings::ssw_cell_type == SSWCellType::To && !entered) {
52,732✔
938
        return;
12,024✔
939
      }
940
    }
941
  }
942

943
  SourceSite site;
129,653✔
944
  site.r = p.r();
129,653✔
945
  site.u = p.u();
129,653✔
946
  site.E = p.E();
129,653✔
947
  site.time = p.time();
129,653✔
948
  site.wgt = p.wgt();
129,653✔
949
  site.delayed_group = p.delayed_group();
129,653✔
950
  site.surf_id = surf.id_;
129,653✔
951
  site.particle = p.type();
129,653✔
952
  site.parent_id = p.id();
129,653✔
953
  site.progeny_id = p.n_progeny();
129,653✔
954
  int64_t idx = simulation::surf_source_bank.thread_safe_append(site);
129,653✔
955
}
956

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