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

openmc-dev / openmc / 27444237628

12 Jun 2026 09:31PM UTC coverage: 81.304% (+0.02%) from 81.281%
27444237628

Pull #3971

github

web-flow
Merge 9487507b3 into 02eb999af
Pull Request #3971: Delta tracking

18480 of 26797 branches covered (68.96%)

Branch coverage included in aggregate %.

592 of 644 new or added lines in 20 files covered. (91.93%)

31 existing lines in 3 files now uncovered.

59819 of 69507 relevant lines covered (86.06%)

49514105.74 hits per line

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

93.74
/src/simulation.cpp
1
#include "openmc/simulation.h"
2

3
#include "openmc/bank.h"
4
#include "openmc/capi.h"
5
#include "openmc/collision_track.h"
6
#include "openmc/container_util.h"
7
#include "openmc/eigenvalue.h"
8
#include "openmc/error.h"
9
#include "openmc/event.h"
10
#include "openmc/geometry_aux.h"
11
#include "openmc/ifp.h"
12
#include "openmc/majorant.h"
13
#include "openmc/material.h"
14
#include "openmc/message_passing.h"
15
#include "openmc/nuclide.h"
16
#include "openmc/output.h"
17
#include "openmc/particle.h"
18
#include "openmc/photon.h"
19
#include "openmc/random_lcg.h"
20
#include "openmc/settings.h"
21
#include "openmc/source.h"
22
#include "openmc/state_point.h"
23
#include "openmc/tallies/derivative.h"
24
#include "openmc/tallies/filter.h"
25
#include "openmc/tallies/tally.h"
26
#include "openmc/tallies/trigger.h"
27
#include "openmc/timer.h"
28
#include "openmc/track_output.h"
29
#include "openmc/weight_windows.h"
30

31
#ifdef _OPENMP
32
#include <omp.h>
33
#endif
34
#include "openmc/tensor.h"
35

36
#ifdef OPENMC_MPI
37
#include <mpi.h>
38
#endif
39

40
#include <fmt/format.h>
41

42
#include <algorithm>
43
#include <cmath>
44
#include <string>
45

46
//==============================================================================
47
// C API functions
48
//==============================================================================
49

50
// OPENMC_RUN encompasses all the main logic where iterations are performed
51
// over the batches, generations, and histories in a fixed source or
52
// k-eigenvalue calculation.
53

54
int openmc_run()
6,726✔
55
{
56
  openmc::simulation::time_total.start();
6,726✔
57
  openmc_simulation_init();
6,726✔
58

59
  // Ensure that a batch isn't executed in the case that the maximum number of
60
  // batches has already been run in a restart statepoint file
61
  int status = 0;
6,726✔
62
  if (openmc::simulation::current_batch >= openmc::settings::n_max_batches) {
6,726✔
63
    status = openmc::STATUS_EXIT_MAX_BATCH;
11✔
64
  }
65

66
  int err = 0;
67
  while (status == 0 && err == 0) {
149,019✔
68
    err = openmc_next_batch(&status);
142,306✔
69
  }
70

71
  openmc_simulation_finalize();
6,713✔
72
  openmc::simulation::time_total.stop();
6,713✔
73
  return err;
6,713✔
74
}
75

76
int openmc_simulation_init()
7,898✔
77
{
78
  using namespace openmc;
7,898✔
79

80
  // Skip if simulation has already been initialized
81
  if (simulation::initialized)
7,898✔
82
    return 0;
83

84
  // Initialize nuclear data (energy limits, log grid)
85
  if (settings::run_CE) {
7,876✔
86
    initialize_data();
6,465✔
87
  }
88

89
  // Create the majorant cross sections for delta tracking.
90
  if (settings::delta_tracking) {
7,876✔
91
    create_majorants();
150✔
92
  }
93

94
  // Determine how much work each process should do
95
  calculate_work(settings::n_particles);
7,876✔
96

97
  // Allocate source, fission and surface source banks.
98
  allocate_banks();
7,876✔
99

100
  // Create track file if needed
101
  if (!settings::track_identifiers.empty() || settings::write_all_tracks) {
7,876✔
102
    open_track_file();
90✔
103
  }
104

105
  // If doing an event-based simulation, intialize the particle buffer
106
  // and event queues
107
  if (settings::event_based) {
7,876✔
108
    int64_t event_buffer_length =
259!
109
      std::min(simulation::work_per_rank, settings::max_particles_in_flight);
259✔
110
    init_event_queues(event_buffer_length);
259✔
111
  }
112

113
  // Allocate tally results arrays if they're not allocated yet
114
  for (auto& t : model::tallies) {
35,527✔
115
    t->set_strides();
27,651✔
116
    t->init_results();
27,651✔
117
  }
118

119
  // Set up material nuclide index mapping
120
  for (auto& mat : model::materials) {
27,823✔
121
    mat->init_nuclide_index();
19,947✔
122
  }
123

124
  // Reset global variables -- this is done before loading state point (as that
125
  // will potentially populate k_generation and entropy)
126
  simulation::current_batch = 0;
7,876✔
127
  simulation::ct_current_file = 1;
7,876✔
128
  simulation::ssw_current_file = 1;
7,876✔
129
  simulation::k_generation.clear();
7,876✔
130
  simulation::entropy.clear();
7,876✔
131
  reset_source_rejection_counters();
7,876✔
132
  openmc_reset();
7,876✔
133

134
  // If this is a restart run, load the state point data and binary source
135
  // file
136
  if (settings::restart_run) {
7,876✔
137
    load_state_point();
63✔
138
    write_message("Resuming simulation...", 6);
126✔
139
  } else {
140
    // Only initialize primary source bank for eigenvalue simulations
141
    if (settings::run_mode == RunMode::EIGENVALUE &&
7,813✔
142
        settings::solver_type == SolverType::MONTE_CARLO) {
4,568✔
143
      initialize_source();
4,197✔
144
    }
145
  }
146

147
  // Display header
148
  if (mpi::master) {
7,876✔
149
    if (settings::run_mode == RunMode::FIXED_SOURCE) {
6,828✔
150
      if (settings::solver_type == SolverType::MONTE_CARLO) {
2,931✔
151
        header("FIXED SOURCE TRANSPORT SIMULATION", 3);
2,555✔
152
      } else if (settings::solver_type == SolverType::RANDOM_RAY) {
376!
153
        header("FIXED SOURCE TRANSPORT SIMULATION (RANDOM RAY SOLVER)", 3);
376✔
154
      }
155
    } else if (settings::run_mode == RunMode::EIGENVALUE) {
3,897!
156
      if (settings::solver_type == SolverType::MONTE_CARLO) {
3,897✔
157
        header("K EIGENVALUE SIMULATION", 3);
3,622✔
158
      } else if (settings::solver_type == SolverType::RANDOM_RAY) {
275!
159
        header("K EIGENVALUE SIMULATION (RANDOM RAY SOLVER)", 3);
275✔
160
      }
161
      if (settings::verbosity >= 7)
3,897✔
162
        print_columns();
3,517✔
163
    }
164
  }
165

166
  // load weight windows from file
167
  if (!settings::weight_windows_file.empty()) {
7,876!
168
    openmc_weight_windows_import(settings::weight_windows_file.c_str());
×
169
  }
170

171
  // Set flag indicating initialization is done
172
  simulation::initialized = true;
7,876✔
173
  return 0;
7,876✔
174
}
175

176
int openmc_simulation_finalize()
7,863✔
177
{
178
  using namespace openmc;
7,863✔
179

180
  // Skip if simulation was never run
181
  if (!simulation::initialized)
7,863!
182
    return 0;
183

184
  // Stop active batch timer and start finalization timer
185
  simulation::time_active.stop();
7,863✔
186
  simulation::time_finalize.start();
7,863✔
187

188
  // Clear material nuclide mapping
189
  for (auto& mat : model::materials) {
27,797✔
190
    mat->mat_nuclide_index_.clear();
39,868!
191
  }
192

193
  // Close track file if open
194
  if (!settings::track_identifiers.empty() || settings::write_all_tracks) {
7,863✔
195
    close_track_file();
90✔
196
  }
197

198
  // Increment total number of generations
199
  simulation::total_gen += simulation::current_batch * settings::gen_per_batch;
7,863✔
200

201
#ifdef OPENMC_MPI
202
  broadcast_results();
3,546✔
203
#endif
204

205
  // Write tally results to tallies.out
206
  if (settings::output_tallies && mpi::master)
7,863!
207
    write_tallies();
6,469✔
208

209
  // If weight window generators are present in this simulation,
210
  // write a weight windows file
211
  if (variance_reduction::weight_windows_generators.size() > 0) {
7,863✔
212
    openmc_weight_windows_export();
166✔
213
  }
214

215
  // Deactivate all tallies
216
  for (auto& t : model::tallies) {
35,514✔
217
    t->active_ = false;
27,651✔
218
  }
219

220
  // Stop timers and show timing statistics
221
  simulation::time_finalize.stop();
7,863✔
222
  simulation::time_total.stop();
7,863✔
223

224
#ifdef OPENMC_MPI
225
  // Reduce track count across ranks for correct reporting. In shared secondary
226
  // bank mode, all ranks already have the global count; in non-shared mode,
227
  // each rank only has its own count.
228
  if (settings::weight_windows_on && !settings::use_shared_secondary_bank) {
3,546✔
229
    int64_t total_tracks;
80✔
230
    MPI_Reduce(&simulation::simulation_tracks_completed, &total_tracks, 1,
80✔
231
      MPI_INT64_T, MPI_SUM, 0, mpi::intracomm);
232
    if (mpi::master)
80✔
233
      simulation::simulation_tracks_completed = total_tracks;
64✔
234
  }
235
#endif
236

237
  if (mpi::master) {
7,863✔
238
    if (settings::solver_type != SolverType::RANDOM_RAY) {
6,815✔
239
      if (settings::verbosity >= 6)
6,164✔
240
        print_runtime();
5,784✔
241
      if (settings::verbosity >= 4)
6,164✔
242
        print_results();
5,784✔
243
    }
244
  }
245
  if (settings::check_overlaps)
7,863!
246
    print_overlap_check();
×
247

248
  // Clear majorants as they could change if OpenMC is run again.
249
  reset_majorants();
7,863✔
250

251
  // Reset flags
252
  simulation::initialized = false;
7,863✔
253
  return 0;
7,863✔
254
}
255

256
int openmc_next_batch(int* status)
146,431✔
257
{
258
  using namespace openmc;
146,431✔
259
  using openmc::simulation::current_gen;
146,431✔
260

261
  // Make sure simulation has been initialized
262
  if (!simulation::initialized) {
146,431✔
263
    set_errmsg("Simulation has not been initialized yet.");
11✔
264
    return OPENMC_E_ALLOCATE;
11✔
265
  }
266

267
  initialize_batch();
146,420✔
268

269
  // =======================================================================
270
  // LOOP OVER GENERATIONS
271
  for (current_gen = 1; current_gen <= settings::gen_per_batch; ++current_gen) {
293,037✔
272

273
    initialize_generation();
146,630✔
274

275
    // Start timer for transport
276
    simulation::time_transport.start();
146,630✔
277

278
    // Transport loop
279
    if (settings::event_based) {
146,630✔
280
      if (settings::use_shared_secondary_bank) {
3,601✔
281
        transport_event_based_shared_secondary();
11✔
282
      } else {
283
        transport_event_based();
3,590✔
284
      }
285
    } else {
286
      if (settings::use_shared_secondary_bank) {
143,029✔
287
        transport_history_based_shared_secondary();
667✔
288
      } else {
289
        transport_history_based();
142,362✔
290
      }
291
    }
292

293
    // Accumulate time for transport
294
    simulation::time_transport.stop();
146,617✔
295

296
    finalize_generation();
146,617✔
297
  }
298

299
  finalize_batch();
146,407✔
300

301
  // Check simulation ending criteria
302
  if (status) {
146,407!
303
    if (simulation::current_batch >= settings::n_max_batches) {
146,407✔
304
      *status = STATUS_EXIT_MAX_BATCH;
6,906✔
305
    } else if (simulation::satisfy_triggers) {
139,501✔
306
      *status = STATUS_EXIT_ON_TRIGGER;
93✔
307
    } else {
308
      *status = STATUS_EXIT_NORMAL;
139,408✔
309
    }
310
  }
311
  return 0;
312
}
313

314
bool openmc_is_statepoint_batch()
3,135✔
315
{
316
  using namespace openmc;
3,135✔
317
  using openmc::simulation::current_gen;
3,135✔
318

319
  if (!simulation::initialized)
3,135!
320
    return false;
321
  else
322
    return contains(settings::statepoint_batch, simulation::current_batch);
6,270✔
323
}
324

325
namespace openmc {
326

327
//==============================================================================
328
// Global variables
329
//==============================================================================
330

331
namespace simulation {
332

333
int ct_current_file;
334
int current_batch;
335
int current_gen;
336
bool initialized {false};
337
double keff {1.0};
338
double keff_std;
339
double k_col_abs {0.0};
340
double k_col_tra {0.0};
341
double k_abs_tra {0.0};
342
double log_spacing;
343
int n_lost_particles {0};
344
bool need_depletion_rx {false};
345
int restart_batch;
346
bool satisfy_triggers {false};
347
int ssw_current_file;
348
int total_gen {0};
349
double total_weight;
350
int64_t work_per_rank;
351

352
const RegularMesh* entropy_mesh {nullptr};
353
const RegularMesh* ufs_mesh {nullptr};
354

355
vector<double> k_generation;
356
vector<int64_t> work_index;
357

358
int64_t simulation_tracks_completed {0};
359

360
} // namespace simulation
361

362
//==============================================================================
363
// Non-member functions
364
//==============================================================================
365

366
void allocate_banks()
7,876✔
367
{
368
  if (settings::run_mode == RunMode::EIGENVALUE &&
7,876✔
369
      settings::solver_type == SolverType::MONTE_CARLO) {
4,631✔
370
    // Allocate source bank
371
    simulation::source_bank.resize(simulation::work_per_rank);
4,260✔
372

373
    // Allocate fission bank
374
    init_fission_bank(3 * simulation::work_per_rank);
4,260✔
375

376
    // Allocate IFP bank
377
    if (settings::ifp_on) {
4,260✔
378
      resize_simulation_ifp_banks();
74✔
379
    }
380
  }
381

382
  if (settings::surf_source_write) {
7,876✔
383
    // Allocate surface source bank
384
    simulation::surf_source_bank.reserve(settings::ssw_max_particles);
1,154✔
385
  }
386

387
  if (settings::collision_track) {
7,876✔
388
    // Allocate collision track bank
389
    collision_track_reserve_bank();
160✔
390
  }
391
}
7,876✔
392

393
void initialize_batch()
167,382✔
394
{
395
  // Increment current batch
396
  ++simulation::current_batch;
167,382✔
397
  if (settings::run_mode == RunMode::FIXED_SOURCE) {
167,382✔
398
    if (settings::solver_type == SolverType::RANDOM_RAY &&
64,939✔
399
        simulation::current_batch < settings::n_inactive + 1) {
14,112✔
400
      write_message(
16,812✔
401
        6, "Simulating batch {:<4} (inactive)", simulation::current_batch);
402
    } else {
403
      write_message(6, "Simulating batch {}", simulation::current_batch);
113,066✔
404
    }
405
  }
406

407
  // Reset total starting particle weight used for normalizing tallies
408
  simulation::total_weight = 0.0;
167,382✔
409

410
  // Determine if this batch is the first inactive or active batch.
411
  bool first_inactive = false;
167,382✔
412
  bool first_active = false;
167,382✔
413
  if (!settings::restart_run) {
167,382✔
414
    first_inactive = settings::n_inactive > 0 && simulation::current_batch == 1;
167,219✔
415
    first_active = simulation::current_batch == settings::n_inactive + 1;
167,219✔
416
  } else if (simulation::current_batch == simulation::restart_batch + 1) {
163✔
417
    first_inactive = simulation::restart_batch < settings::n_inactive;
52✔
418
    first_active = !first_inactive;
52✔
419
  }
420

421
  // Manage active/inactive timers and activate tallies if necessary.
422
  if (first_inactive) {
167,271✔
423
    simulation::time_inactive.start();
4,015✔
424
  } else if (first_active) {
163,367✔
425
    simulation::time_inactive.stop();
7,829✔
426
    simulation::time_active.start();
7,829✔
427
    for (auto& t : model::tallies) {
35,458✔
428
      t->active_ = true;
27,629✔
429
    }
430
  }
431

432
  // Add user tallies to active tallies list
433
  setup_active_tallies();
167,382✔
434
}
167,382✔
435

436
void finalize_batch()
167,369✔
437
{
438
  // Reduce tallies onto master process and accumulate
439
  simulation::time_tallies.start();
167,369✔
440
  accumulate_tallies();
167,369✔
441
  simulation::time_tallies.stop();
167,369✔
442

443
  // update weight windows if needed
444
  for (const auto& wwg : variance_reduction::weight_windows_generators) {
170,001✔
445
    wwg->update();
2,632✔
446
  }
447

448
  // Reset global tally results
449
  if (simulation::current_batch <= settings::n_inactive) {
167,369✔
450
    simulation::global_tallies.fill(0.0);
33,353✔
451
    simulation::n_realizations = 0;
33,353✔
452
  }
453

454
  // Check_triggers
455
  if (mpi::master)
167,369✔
456
    check_triggers();
148,050✔
457
#ifdef OPENMC_MPI
458
  MPI_Bcast(&simulation::satisfy_triggers, 1, MPI_C_BOOL, 0, mpi::intracomm);
73,480✔
459
#endif
460
  if (simulation::satisfy_triggers ||
167,369✔
461
      (settings::trigger_on &&
2,567✔
462
        simulation::current_batch == settings::n_max_batches)) {
2,567✔
463
    settings::statepoint_batch.insert(simulation::current_batch);
141✔
464
  }
465

466
  // Write out state point if it's been specified for this batch and is not
467
  // a CMFD run instance
468
  if (contains(settings::statepoint_batch, simulation::current_batch) &&
334,738✔
469
      !settings::cmfd_run) {
8,131✔
470
    if (contains(settings::sourcepoint_batch, simulation::current_batch) &&
15,632✔
471
        settings::source_write && !settings::source_separate) {
14,749✔
472
      bool b = (settings::run_mode == RunMode::EIGENVALUE);
6,723✔
473
      openmc_statepoint_write(nullptr, &b);
6,723✔
474
    } else {
475
      bool b = false;
1,232✔
476
      openmc_statepoint_write(nullptr, &b);
1,232✔
477
    }
478
  }
479

480
  if (settings::run_mode == RunMode::EIGENVALUE) {
167,369✔
481
    // Write out a separate source point if it's been specified for this batch
482
    if (contains(settings::sourcepoint_batch, simulation::current_batch) &&
107,100✔
483
        settings::source_write && settings::source_separate) {
106,729✔
484

485
      // Determine width for zero padding
486
      int w = std::to_string(settings::n_max_batches).size();
71✔
487
      std::string source_point_filename = fmt::format("{0}source.{1:0{2}}",
71✔
488
        settings::path_output, simulation::current_batch, w);
71✔
489
      span<SourceSite> bankspan(simulation::source_bank);
71✔
490
      write_source_point(source_point_filename, bankspan,
142✔
491
        simulation::work_index, settings::source_mcpl_write);
492
    }
71✔
493

494
    // Write a continously-overwritten source point if requested.
495
    if (settings::source_latest) {
102,443✔
496
      auto filename = settings::path_output + "source";
150✔
497
      span<SourceSite> bankspan(simulation::source_bank);
150✔
498
      write_source_point(filename, bankspan, simulation::work_index,
300✔
499
        settings::source_mcpl_write);
500
    }
150✔
501
  }
502

503
  // Write out surface source if requested.
504
  if (settings::surf_source_write &&
167,369✔
505
      simulation::ssw_current_file <= settings::ssw_max_files) {
17,669✔
506
    bool last_batch = (simulation::current_batch == settings::n_batches);
1,976✔
507
    if (simulation::surf_source_bank.full() || last_batch) {
1,976✔
508
      // Determine appropriate filename
509
      auto filename = fmt::format("{}surface_source.{}", settings::path_output,
1,187✔
510
        simulation::current_batch);
1,187✔
511
      if (settings::ssw_max_files == 1 ||
1,187✔
512
          (simulation::ssw_current_file == 1 && last_batch)) {
55!
513
        filename = settings::path_output + "surface_source";
1,132✔
514
      }
515

516
      // Get span of source bank and calculate parallel index vector
517
      auto surf_work_index = mpi::calculate_parallel_index_vector(
1,187✔
518
        simulation::surf_source_bank.size());
1,187✔
519
      span<SourceSite> surfbankspan(simulation::surf_source_bank.begin(),
1,187✔
520
        simulation::surf_source_bank.size());
1,187✔
521

522
      // Write surface source file
523
      write_source_point(
1,187✔
524
        filename, surfbankspan, surf_work_index, settings::surf_mcpl_write);
525

526
      // Reset surface source bank and increment counter
527
      simulation::surf_source_bank.clear();
1,187✔
528
      if (!last_batch && settings::ssw_max_files >= 1) {
1,187!
529
        simulation::surf_source_bank.reserve(settings::ssw_max_particles);
1,005✔
530
      }
531
      ++simulation::ssw_current_file;
1,187✔
532
    }
1,187✔
533
  }
534
  // Write collision track file if requested
535
  if (settings::collision_track) {
167,369✔
536
    collision_track_flush_bank();
580✔
537
  }
538
}
167,369✔
539

540
void initialize_generation()
167,592✔
541
{
542
  if (settings::run_mode == RunMode::EIGENVALUE) {
167,592✔
543
    // Clear out the fission bank
544
    simulation::fission_bank.resize(0);
102,653✔
545

546
    // Count source sites if using uniform fission source weighting
547
    if (settings::ufs_on)
102,653✔
548
      ufs_count_sites();
150✔
549

550
    // Store current value of tracklength k
551
    if (settings::delta_tracking) {
102,653✔
552
      simulation::keff_generation = simulation::global_tallies(
1,500✔
553
        GlobalTally::K_COLLISION, TallyResult::VALUE);
554
    } else {
555
      simulation::keff_generation = simulation::global_tallies(
101,153✔
556
        GlobalTally::K_TRACKLENGTH, TallyResult::VALUE);
557
    }
558
  }
559
}
167,592✔
560

561
void finalize_generation()
167,579✔
562
{
563
  auto& gt = simulation::global_tallies;
167,579✔
564

565
  // Update global tallies with the accumulation variables
566
  if (settings::run_mode == RunMode::EIGENVALUE) {
167,579✔
567
    gt(GlobalTally::K_COLLISION, TallyResult::VALUE) += global_tally_collision;
102,653✔
568
    gt(GlobalTally::K_ABSORPTION, TallyResult::VALUE) +=
102,653✔
569
      global_tally_absorption;
570
    gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) +=
102,653✔
571
      global_tally_tracklength;
572
  }
573
  gt(GlobalTally::LEAKAGE, TallyResult::VALUE) += global_tally_leakage;
167,579✔
574

575
  // reset tallies
576
  if (settings::run_mode == RunMode::EIGENVALUE) {
167,579✔
577
    global_tally_collision = 0.0;
102,653✔
578
    global_tally_absorption = 0.0;
102,653✔
579
    global_tally_tracklength = 0.0;
102,653✔
580
  }
581
  global_tally_leakage = 0.0;
167,579✔
582

583
  if (settings::run_mode == RunMode::EIGENVALUE &&
167,579✔
584
      settings::solver_type == SolverType::MONTE_CARLO) {
102,653✔
585
    // If using shared memory, stable sort the fission bank (by parent IDs)
586
    // so as to allow for reproducibility regardless of which order particles
587
    // are run in.
588
    sort_bank(simulation::fission_bank, true);
95,803✔
589

590
    // Distribute fission bank across processors evenly
591
    synchronize_bank();
95,803✔
592
  }
593

594
  if (settings::run_mode == RunMode::EIGENVALUE) {
167,579✔
595

596
    // Calculate shannon entropy
597
    if (settings::entropy_on &&
102,653✔
598
        settings::solver_type == SolverType::MONTE_CARLO)
14,535✔
599
      shannon_entropy();
7,685✔
600

601
    // Collect results and statistics
602
    calculate_generation_keff();
102,653✔
603
    calculate_average_keff();
102,653✔
604

605
    // Write generation output
606
    if (mpi::master && settings::verbosity >= 7) {
102,653✔
607
      print_generation();
77,188✔
608
    }
609
  }
610
}
167,579✔
611

612
void sample_source_particle(Particle& p, int64_t index_source)
178,937,692✔
613
{
614
  // Sample a particle from the source bank
615
  if (settings::run_mode == RunMode::EIGENVALUE) {
178,937,692✔
616
    p.from_source(&simulation::source_bank[index_source - 1]);
151,004,000✔
617
  } else if (settings::run_mode == RunMode::FIXED_SOURCE) {
27,933,692!
618
    // initialize random number seed
619
    int64_t id = compute_transport_seed(compute_particle_id(index_source));
27,933,692✔
620
    uint64_t seed = init_seed(id, STREAM_SOURCE);
27,933,692✔
621
    // sample from external source distribution or custom library then set
622
    auto site = sample_external_source(&seed);
27,933,692✔
623
    p.from_source(&site);
27,933,688✔
624
  }
625
}
178,937,688✔
626

627
void initialize_particle_track(
200,270,276✔
628
  Particle& p, int64_t index_source, bool is_secondary)
629
{
630
  // Note: index_source is 1-based (first particle = 1), but current_work() is
631
  // stored as 0-based for direct use as an array index into
632
  // progeny_per_particle, source_bank, ifp banks, etc.
633
  if (!is_secondary) {
200,270,276✔
634
    sample_source_particle(p, index_source);
178,937,692✔
635
  }
636

637
  p.current_work() = index_source - 1;
200,270,272✔
638

639
  // set identifier for particle
640
  p.id() = compute_particle_id(index_source);
200,270,272✔
641

642
  // set progeny count to zero
643
  p.n_progeny() = 0;
200,270,272✔
644

645
  // Reset particle event counter
646
  p.n_event() = 0;
200,270,272✔
647

648
  // Initialize track counter (1 for this primary/secondary track)
649
  p.n_tracks() = 1;
200,270,272✔
650

651
  // Reset split counter
652
  p.n_split() = 0;
200,270,272✔
653

654
  // Reset weight window ratio
655
  p.ww_factor() = 0.0;
200,270,272✔
656

657
  // set particle history start weight
658
  p.wgt_born() = p.wgt();
200,270,272✔
659

660
  // Reset pulse_height_storage
661
  std::fill(p.pht_storage().begin(), p.pht_storage().end(), 0);
200,270,272✔
662

663
  // set random number seed
664
  int64_t particle_seed = compute_transport_seed(p.id());
200,270,272✔
665
  init_particle_seeds(particle_seed, p.seeds());
200,270,272✔
666

667
  // set particle trace
668
  p.trace() = false;
200,270,272✔
669
  if (simulation::current_batch == settings::trace_batch &&
200,281,272✔
670
      simulation::current_gen == settings::trace_gen &&
200,270,272!
671
      p.id() == settings::trace_particle)
11,000✔
672
    p.trace() = true;
11✔
673

674
  // Set particle track.
675
  p.write_track() = check_track_criteria(p);
200,270,272✔
676

677
  // Set the particle's initial weight window value.
678
  if (!is_secondary) {
200,270,272✔
679
    p.wgt_ww_born() = -1.0;
178,937,688✔
680
    apply_weight_windows(p);
178,937,688✔
681
  }
682

683
  // Display message if high verbosity or trace is on
684
  if (settings::verbosity >= 9 || p.trace()) {
200,270,272!
685
    write_message("Simulating Particle {}", p.id());
22✔
686
  }
687

688
  // Compute the majorant.
689
  if (settings::delta_tracking) {
200,270,272✔
690
    p.update_majorant();
1,100,000✔
691
  }
692

693
  // Add particle's starting weight to count for normalizing tallies later
694
  if (!is_secondary) {
200,270,272✔
695
#pragma omp atomic
98,098,109✔
696
    simulation::total_weight += p.wgt();
178,937,688✔
697
  }
698

699
  // Force calculation of cross-sections by setting last energy to zero
700
  if (settings::run_CE) {
200,270,272✔
701
    p.invalidate_neutron_xs();
85,721,728✔
702
  }
703

704
  // Prepare to write out particle track.
705
  if (p.write_track())
200,270,272✔
706
    add_particle_track(p);
999✔
707
}
200,270,272✔
708

709
int overall_generation()
205,578,307✔
710
{
711
  using namespace simulation;
205,578,307✔
712
  return settings::gen_per_batch * (current_batch - 1) + current_gen;
205,578,307✔
713
}
714

715
int64_t compute_particle_id(int64_t index_source)
228,204,239✔
716
{
717
  if (settings::use_shared_secondary_bank) {
228,204,239✔
718
    return simulation::work_index[mpi::rank] + index_source +
22,901,609✔
719
           simulation::simulation_tracks_completed;
22,901,609✔
720
  } else {
721
    return simulation::work_index[mpi::rank] + index_source;
205,302,630✔
722
  }
723
}
724

725
int64_t compute_transport_seed(int64_t particle_id)
228,204,283✔
726
{
727
  if (settings::use_shared_secondary_bank) {
228,204,283✔
728
    return particle_id;
729
  } else {
730
    return (simulation::total_gen + overall_generation() - 1) *
205,302,663✔
731
             settings::n_particles +
732
           particle_id;
205,302,663✔
733
  }
734
}
735

736
void calculate_work(int64_t n_particles)
17,202✔
737
{
738
  // Determine minimum amount of particles to simulate on each processor
739
  int64_t min_work = n_particles / mpi::n_procs;
17,202✔
740

741
  // Determine number of processors that have one extra particle
742
  int64_t remainder = n_particles % mpi::n_procs;
17,202✔
743

744
  int64_t i_bank = 0;
17,202✔
745
  simulation::work_index.resize(mpi::n_procs + 1);
17,202✔
746
  simulation::work_index[0] = 0;
17,202✔
747
  for (int i = 0; i < mpi::n_procs; ++i) {
40,027✔
748
    // Number of particles for rank i
749
    int64_t work_i = i < remainder ? min_work + 1 : min_work;
22,825✔
750

751
    // Set number of particles
752
    if (mpi::rank == i)
22,825✔
753
      simulation::work_per_rank = work_i;
17,202✔
754

755
    // Set index into source bank for rank i
756
    i_bank += work_i;
22,825✔
757
    simulation::work_index[i + 1] = i_bank;
22,825✔
758
  }
759
}
17,202✔
760

761
void initialize_data()
6,509✔
762
{
763
  // Determine minimum/maximum energy for incident neutron/photon data
764
  data::energy_max = {INFTY, INFTY, INFTY, INFTY};
6,509✔
765
  data::energy_min = {0.0, 0.0, 0.0, 0.0};
6,509✔
766
  int neutron = ParticleType::neutron().transport_index();
6,509✔
767
  int photon = ParticleType::photon().transport_index();
6,509✔
768
  int electron = ParticleType::electron().transport_index();
6,509✔
769
  int positron = ParticleType::positron().transport_index();
6,509✔
770

771
  for (const auto& nuc : data::nuclides) {
39,751✔
772
    if (nuc->grid_.size() >= 1) {
33,242!
773
      data::energy_min[neutron] =
33,242✔
774
        std::max(data::energy_min[neutron], nuc->grid_[0].energy.front());
39,230✔
775
      data::energy_max[neutron] =
33,242✔
776
        std::min(data::energy_max[neutron], nuc->grid_[0].energy.back());
40,794✔
777
    }
778
  }
779

780
  if (settings::photon_transport) {
6,509✔
781
    for (const auto& elem : data::elements) {
1,805✔
782
      if (elem->energy_.size() >= 1) {
1,281!
783
        int n = elem->energy_.size();
1,281✔
784
        data::energy_min[photon] =
2,562✔
785
          std::max(data::energy_min[photon], std::exp(elem->energy_(1)));
2,120✔
786
        data::energy_max[photon] =
1,281✔
787
          std::min(data::energy_max[photon], std::exp(elem->energy_(n - 1)));
1,805✔
788
      }
789
    }
790

791
    if (settings::electron_treatment == ElectronTreatment::TTB) {
524✔
792
      // Determine if minimum/maximum energy for bremsstrahlung is greater/less
793
      // than the current minimum/maximum
794
      if (data::ttb_e_grid.size() >= 1) {
465!
795
        int n_e = data::ttb_e_grid.size();
465✔
796

797
        const std::vector<int> charged = {electron, positron};
465✔
798
        for (auto t : charged) {
1,395✔
799
          data::energy_min[t] = std::exp(data::ttb_e_grid(1));
930✔
800
          data::energy_max[t] = std::exp(data::ttb_e_grid(n_e - 1));
930✔
801
        }
802

803
        data::energy_min[photon] =
930✔
804
          std::max(data::energy_min[photon], data::energy_min[electron]);
930!
805

806
        data::energy_max[photon] =
930✔
807
          std::min(data::energy_max[photon], data::energy_max[electron]);
930!
808
      }
465✔
809
    }
810
  }
811

812
  // Show which nuclide results in lowest energy for neutron transport
813
  for (const auto& nuc : data::nuclides) {
8,133✔
814
    // If a nuclide is present in a material that's not used in the model, its
815
    // grid has not been allocated
816
    if (nuc->grid_.size() > 0) {
7,612!
817
      double max_E = nuc->grid_[0].energy.back();
7,612✔
818
      if (max_E == data::energy_max[neutron]) {
7,612✔
819
        write_message(7, "Maximum neutron transport energy: {} eV for {}",
5,988✔
820
          data::energy_max[neutron], nuc->name_);
5,988✔
821
        if (mpi::master && data::energy_max[neutron] < 20.0e6) {
5,988!
822
          warning("Maximum neutron energy is below 20 MeV. This may bias "
×
823
                  "the results.");
824
        }
825
        break;
826
      }
827
    }
828
  }
829

830
  // Set up logarithmic grid for nuclides
831
  for (auto& nuc : data::nuclides) {
39,751✔
832
    nuc->init_grid();
33,242✔
833
  }
834
  simulation::log_spacing =
13,018✔
835
    std::log(data::energy_max[neutron] / data::energy_min[neutron]) /
6,509✔
836
    settings::n_log_bins;
837
}
6,509✔
838

839
#ifdef OPENMC_MPI
840
void broadcast_results()
3,546✔
841
{
842
  // Broadcast tally results so that each process has access to results
843
  for (auto& t : model::tallies) {
17,291✔
844
    // Create a new datatype that consists of all values for a given filter
845
    // bin and then use that to broadcast. This is done to minimize the
846
    // chance of the 'count' argument of MPI_BCAST exceeding 2**31
847
    auto& results = t->results_;
13,745✔
848

849
    auto shape = results.shape();
13,745✔
850
    int count_per_filter = shape[1] * shape[2];
13,745✔
851
    MPI_Datatype result_block;
13,745✔
852
    MPI_Type_contiguous(count_per_filter, MPI_DOUBLE, &result_block);
13,745✔
853
    MPI_Type_commit(&result_block);
13,745✔
854
    MPI_Bcast(results.data(), shape[0], result_block, 0, mpi::intracomm);
13,745✔
855
    MPI_Type_free(&result_block);
13,745✔
856
  }
13,745✔
857

858
  // Also broadcast global tally results
859
  auto& gt = simulation::global_tallies;
3,546✔
860
  MPI_Bcast(gt.data(), gt.size(), MPI_DOUBLE, 0, mpi::intracomm);
3,546✔
861

862
  // These guys are needed so that non-master processes can calculate the
863
  // combined estimate of k-effective
864
  double temp[] {
3,546✔
865
    simulation::k_col_abs, simulation::k_col_tra, simulation::k_abs_tra};
3,546✔
866
  MPI_Bcast(temp, 3, MPI_DOUBLE, 0, mpi::intracomm);
3,546✔
867
  simulation::k_col_abs = temp[0];
3,546✔
868
  simulation::k_col_tra = temp[1];
3,546✔
869
  simulation::k_abs_tra = temp[2];
3,546✔
870
}
3,546✔
871

872
#endif
873

874
void free_memory_simulation()
9,004✔
875
{
876
  simulation::k_generation.clear();
9,004✔
877
  simulation::entropy.clear();
9,004✔
878
}
9,004✔
879

880
void transport_history_based_single_particle(Particle& p)
186,834,506✔
881
{
882
  while (p.alive()) {
2,147,483,647✔
883
    p.event_calculate_xs();
2,147,483,647✔
884
    if (p.alive()) {
2,147,483,647!
885
      p.event_advance();
2,147,483,647✔
886
    }
887
    if (p.alive()) {
2,147,483,647✔
888
      if (p.collision_distance() > p.boundary().distance()) {
2,147,483,647✔
889
        p.event_cross_surface();
2,147,483,647✔
890
      } else if (p.alive()) {
2,147,483,647✔
891
        p.event_collide();
2,147,483,647✔
892
      }
893
    }
894
    p.event_check_limit_and_revive();
2,147,483,647✔
895
  }
896
  p.event_death();
186,834,497✔
897
}
186,834,497✔
898

899
void transport_delta_history_based_single_particle(Particle& p)
800,000✔
900
{
901
  p.delta_tracking() = true;
800,000✔
902

903
  while (p.alive()) {
166,263,420✔
904
    p.event_delta_advance();
165,463,420✔
905

906
    if (p.alive()) {
165,463,420!
907
      // Electrons and positrons collide in-place, no need to rejection sample.
908
      if (p.type() == ParticleType::electron() ||
165,463,420✔
909
          p.type() == ParticleType::positron()) {
72,980,960✔
910
        p.event_collide();
92,565,270✔
911
      }
912

913
      if (p.collision_distance() < p.boundary().distance()) {
165,463,420✔
914
        // Collided before hitting an external boundary. Rejection sample the
915
        // majorant.
916
        p.event_calculate_xs();
156,078,460✔
917
        if (p.alive() && (p.macro_xs().total / p.majorant() > 1.0)) {
156,078,460!
NEW
918
          p.mark_as_lost(fmt::format(
×
919
            "Ratio of the total cross section ({}) to the majorant "
920
            "cross section ({}) for particle {} ({}) with energy {} is "
921
            "greater than unity!",
NEW
922
            p.macro_xs().total, p.majorant(), p.id(), p.type().str(), p.E()));
×
NEW
923
          break;
×
924
        }
925
        if (p.alive() &&
156,078,460✔
926
            (prn(p.current_seed()) < (p.macro_xs().total / p.majorant()))) {
63,513,190✔
927
          p.event_collide();
22,562,390✔
928
        }
929
      } else {
930
        // Crossed an external boundary before colliding.
931
        p.event_cross_surface();
9,384,960✔
932
      }
933
    }
934

935
    p.event_check_limit_and_revive();
165,463,420✔
936
  }
937
  p.event_death();
800,000✔
938
}
800,000✔
939

940
void transport_history_based()
142,362✔
941
{
942
#pragma omp parallel for schedule(runtime)
79,337✔
943
  for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
80,861,175✔
944
    Particle p;
80,798,159✔
945
    initialize_particle_track(p, i_work, false);
80,798,159✔
946
    if (settings::delta_tracking) {
80,798,155✔
947
      transport_delta_history_based_single_particle(p);
400,000✔
948
    } else {
949
      transport_history_based_single_particle(p);
80,398,155✔
950
    }
951
  }
80,798,150✔
952
}
142,353✔
953

954
// The shared secondary bank transport algorithm works in two phases. In the
955
// first phase, all primary particles are sampled then transported, and their
956
// secondary particles are deposited into a shared secondary bank. The second
957
// phase occurs in a loop, where all secondary tracks in the shared secondary
958
// bank are transported. Any secondary particles generated during this phase are
959
// deposited back into the shared secondary bank. The shared secondary bank is
960
// sorted for consistent ordering and load balanced across MPI ranks. This loop
961
// continues until there are no more secondary tracks left to transport.
962
void transport_history_based_shared_secondary()
667✔
963
{
964
  // Clear shared secondary banks from any prior use
965
  simulation::shared_secondary_bank_read.clear();
667✔
966
  simulation::shared_secondary_bank_write.clear();
667✔
967

968
  if (mpi::master) {
667✔
969
    write_message(fmt::format(" Primary source          particles: {}",
1,150✔
970
                    settings::n_particles),
971
      6);
972
  }
973

974
  simulation::progeny_per_particle.resize(simulation::work_per_rank);
667✔
975
  std::fill(simulation::progeny_per_particle.begin(),
1,334✔
976
    simulation::progeny_per_particle.end(), 0);
667✔
977

978
  // Phase 1: Transport primary particles and deposit first generation of
979
  // secondaries in the shared secondary bank
980
#pragma omp parallel
384✔
981
  {
283✔
982
    vector<SourceSite> thread_bank;
283✔
983

984
#pragma omp for schedule(runtime)
985
    for (int64_t i = 1; i <= simulation::work_per_rank; i++) {
356,058✔
986
      Particle p;
355,775✔
987
      initialize_particle_track(p, i, false);
355,775✔
988
      if (settings::delta_tracking) {
355,775!
989
        transport_delta_history_based_single_particle(p);
×
990
      } else {
991
        transport_history_based_single_particle(p);
355,775✔
992
      }
993
      for (auto& site : p.local_secondary_bank()) {
1,087,225✔
994
        thread_bank.push_back(site);
731,450✔
995
      }
996
    }
355,775✔
997

998
    // Drain thread-local bank into the shared secondary bank (once per thread)
999
#pragma omp critical(SharedSecondaryBank)
1000
    {
283✔
1001
      for (auto& site : thread_bank) {
731,733✔
1002
        simulation::shared_secondary_bank_write.thread_unsafe_append(site);
731,450✔
1003
      }
1004
    }
1005
  }
1006

1007
  simulation::simulation_tracks_completed += settings::n_particles;
667✔
1008

1009
  // Phase 2: Now that the secondary bank has been populated, enter loop over
1010
  // all secondary generations
1011
  int n_generation_depth = 1;
667✔
1012
  int64_t alive_secondary = 1;
667✔
1013
  while (alive_secondary) {
8,909✔
1014

1015
    // Sort the shared secondary bank by parent ID then progeny ID to
1016
    // ensure reproducibility.
1017
    sort_bank(simulation::shared_secondary_bank_write, false);
8,242✔
1018

1019
    // Synchronize the shared secondary bank amongst all MPI ranks, such
1020
    // that each MPI rank has an approximately equal number of secondary
1021
    // tracks. Also reports the total number of secondaries alive across
1022
    // all MPI ranks.
1023
    alive_secondary = synchronize_global_secondary_bank(
8,242✔
1024
      simulation::shared_secondary_bank_write);
1025

1026
    // Recalculate work for each MPI rank based on number of alive secondary
1027
    // tracks
1028
    calculate_work(alive_secondary);
8,242✔
1029

1030
    // Display the number of secondary tracks in this generation. This
1031
    // is useful for user monitoring so as to see if the secondary population is
1032
    // exploding and to determine how many generations of secondaries are being
1033
    // transported.
1034
    if (mpi::master) {
8,242✔
1035
      write_message(fmt::format(" Secondary generation {:<2}    tracks: {}",
13,140✔
1036
                      n_generation_depth, alive_secondary),
1037
        6);
1038
    }
1039

1040
    simulation::shared_secondary_bank_read =
8,242✔
1041
      std::move(simulation::shared_secondary_bank_write);
8,242✔
1042
    simulation::shared_secondary_bank_write = SharedArray<SourceSite>();
8,242!
1043
    simulation::progeny_per_particle.resize(
8,242✔
1044
      simulation::shared_secondary_bank_read.size());
8,242✔
1045
    std::fill(simulation::progeny_per_particle.begin(),
16,484✔
1046
      simulation::progeny_per_particle.end(), 0);
8,242✔
1047

1048
    // Transport all secondary tracks from the shared secondary bank
1049
#pragma omp parallel
4,649✔
1050
    {
3,593✔
1051
      vector<SourceSite> thread_bank;
3,593✔
1052

1053
#pragma omp for schedule(runtime)
1054
      for (int64_t i = 1; i <= simulation::shared_secondary_bank_read.size();
9,727,253✔
1055
           i++) {
1056
        Particle p;
9,723,660✔
1057
        initialize_particle_track(p, i, true);
9,723,660✔
1058
        SourceSite& site = simulation::shared_secondary_bank_read[i - 1];
9,723,660✔
1059
        p.event_revive_from_secondary(site);
9,723,660✔
1060
        if (settings::delta_tracking) {
9,723,660!
1061
          transport_delta_history_based_single_particle(p);
×
1062
        } else {
1063
          transport_history_based_single_particle(p);
9,723,660✔
1064
        }
1065
        for (auto& secondary_site : p.local_secondary_bank()) {
18,715,870✔
1066
          thread_bank.push_back(secondary_site);
8,992,210✔
1067
        }
1068
      }
9,723,660✔
1069

1070
      // Drain thread-local bank into the shared secondary bank (once per
1071
      // thread)
1072
#pragma omp critical(SharedSecondaryBank)
1073
      {
3,593✔
1074
        for (auto& secondary_site : thread_bank) {
8,995,803✔
1075
          simulation::shared_secondary_bank_write.thread_unsafe_append(
8,992,210✔
1076
            secondary_site);
1077
        }
1078
      }
1079
    } // End of transport loop over tracks in shared secondary bank
3,593✔
1080
    n_generation_depth++;
8,242✔
1081
    simulation::simulation_tracks_completed += alive_secondary;
8,242✔
1082
  } // End of loop over secondary generations
1083

1084
  // Reset work so that fission bank etc works correctly
1085
  calculate_work(settings::n_particles);
667✔
1086
}
667✔
1087

1088
void transport_event_based()
3,590✔
1089
{
1090
  int64_t remaining_work = simulation::work_per_rank;
3,590✔
1091
  int64_t source_offset = 0;
3,590✔
1092

1093
  // To cap the total amount of memory used to store particle object data, the
1094
  // number of particles in flight at any point in time can bet set. In the case
1095
  // that the maximum in flight particle count is lower than the total number
1096
  // of particles that need to be run this iteration, the event-based transport
1097
  // loop is executed multiple times until all particles have been completed.
1098
  while (remaining_work > 0) {
7,180✔
1099
    // Figure out # of particles to run for this subiteration
1100
    int64_t n_particles =
3,590!
1101
      std::min(remaining_work, settings::max_particles_in_flight);
3,590✔
1102

1103
    // Initialize all particle histories for this subiteration
1104
    if (settings::delta_tracking) {
3,590✔
1105
      process_delta_init_events(n_particles, source_offset);
380✔
1106
      process_delta_transport_events();
380✔
1107
    } else {
1108
      process_init_events(n_particles, source_offset);
3,210✔
1109
      process_transport_events();
3,210✔
1110
    }
1111
    process_death_events(n_particles);
3,590✔
1112

1113
    // Adjust remaining work and source offset variables
1114
    remaining_work -= n_particles;
3,590✔
1115
    source_offset += n_particles;
3,590✔
1116
  }
1117
}
3,590✔
1118

1119
void transport_event_based_shared_secondary()
11✔
1120
{
1121
  // Clear shared secondary banks from any prior use
1122
  simulation::shared_secondary_bank_read.clear();
11✔
1123
  simulation::shared_secondary_bank_write.clear();
11✔
1124

1125
  if (mpi::master) {
11!
1126
    write_message(fmt::format(" Primary source          particles: {}",
22!
1127
                    settings::n_particles),
1128
      6);
1129
  }
1130

1131
  simulation::progeny_per_particle.resize(simulation::work_per_rank);
11✔
1132
  std::fill(simulation::progeny_per_particle.begin(),
22✔
1133
    simulation::progeny_per_particle.end(), 0);
11✔
1134

1135
  // Phase 1: Transport primary particles using event-based processing and
1136
  // deposit first generation of secondaries in the shared secondary bank
1137
  int64_t remaining_work = simulation::work_per_rank;
11✔
1138
  int64_t source_offset = 0;
11✔
1139

1140
  while (remaining_work > 0) {
22✔
1141
    int64_t n_particles =
11!
1142
      std::min(remaining_work, settings::max_particles_in_flight);
11✔
1143

1144
    if (settings::delta_tracking) {
11!
NEW
1145
      process_delta_init_events(n_particles, source_offset);
×
NEW
UNCOV
1146
      process_delta_transport_events();
×
1147
    } else {
1148
      process_init_events(n_particles, source_offset);
11✔
1149
      process_transport_events();
11✔
1150
    }
1151
    process_death_events(n_particles);
11✔
1152

1153
    // Collect secondaries from all particle buffers into shared bank
1154
    for (int64_t i = 0; i < n_particles; i++) {
1,661✔
1155
      for (auto& site : simulation::particles[i].local_secondary_bank()) {
6,132✔
1156
        simulation::shared_secondary_bank_write.thread_unsafe_append(site);
4,482✔
1157
      }
1158
      simulation::particles[i].local_secondary_bank().clear();
3,017✔
1159
    }
1160

1161
    remaining_work -= n_particles;
11✔
1162
    source_offset += n_particles;
11✔
1163
  }
1164

1165
  simulation::simulation_tracks_completed += settings::n_particles;
11✔
1166

1167
  // Phase 2: Now that the secondary bank has been populated, enter loop over
1168
  // all secondary generations
1169
  int n_generation_depth = 1;
11✔
1170
  int64_t alive_secondary = 1;
11✔
1171
  while (alive_secondary) {
417✔
1172

1173
    // Sort the shared secondary bank by parent ID then progeny ID to
1174
    // ensure reproducibility.
1175
    sort_bank(simulation::shared_secondary_bank_write, false);
406✔
1176

1177
    // Synchronize the shared secondary bank amongst all MPI ranks, such
1178
    // that each MPI rank has an approximately equal number of secondary
1179
    // tracks.
1180
    alive_secondary = synchronize_global_secondary_bank(
406✔
1181
      simulation::shared_secondary_bank_write);
1182

1183
    // Recalculate work for each MPI rank based on number of alive secondary
1184
    // tracks
1185
    calculate_work(alive_secondary);
406✔
1186

1187
    if (mpi::master) {
406!
1188
      write_message(fmt::format(" Secondary generation {:<2}    tracks: {}",
812!
1189
                      n_generation_depth, alive_secondary),
1190
        6);
1191
    }
1192

1193
    simulation::shared_secondary_bank_read =
406✔
1194
      std::move(simulation::shared_secondary_bank_write);
406✔
1195
    simulation::shared_secondary_bank_write = SharedArray<SourceSite>();
406!
1196
    simulation::progeny_per_particle.resize(
406✔
1197
      simulation::shared_secondary_bank_read.size());
406✔
1198
    std::fill(simulation::progeny_per_particle.begin(),
812✔
1199
      simulation::progeny_per_particle.end(), 0);
406✔
1200

1201
    // Ensure particle buffer is large enough for this secondary generation
1202
    int64_t sec_buffer_length = std::min(
406!
1203
      static_cast<int64_t>(simulation::shared_secondary_bank_read.size()),
406!
1204
      settings::max_particles_in_flight);
406✔
1205
    if (sec_buffer_length >
406✔
1206
        static_cast<int64_t>(simulation::particles.size())) {
406✔
1207
      init_event_queues(sec_buffer_length);
34✔
1208
    }
1209

1210
    // Transport secondary tracks using event-based processing
1211
    int64_t sec_remaining = simulation::shared_secondary_bank_read.size();
406✔
1212
    int64_t sec_offset = 0;
406✔
1213

1214
    while (sec_remaining > 0) {
801✔
1215
      int64_t n_particles =
395!
1216
        std::min(sec_remaining, settings::max_particles_in_flight);
395✔
1217

1218
      if (settings::delta_tracking) {
395!
NEW
1219
        process_delta_init_secondary_events(
×
1220
          n_particles, sec_offset, simulation::shared_secondary_bank_read);
NEW
1221
        process_delta_transport_events();
×
1222
      } else {
1223
        process_init_secondary_events(
395✔
1224
          n_particles, sec_offset, simulation::shared_secondary_bank_read);
1225
        process_transport_events();
395✔
1226
      }
1227
      process_death_events(n_particles);
395✔
1228

1229
      // Collect secondaries from all particle buffers into shared bank
1230
      for (int64_t i = 0; i < n_particles; i++) {
180,005✔
1231
        for (auto& site : simulation::particles[i].local_secondary_bank()) {
354,738✔
1232
          simulation::shared_secondary_bank_write.thread_unsafe_append(site);
175,128✔
1233
        }
1234
        simulation::particles[i].local_secondary_bank().clear();
226,133✔
1235
      }
1236

1237
      sec_remaining -= n_particles;
395✔
1238
      sec_offset += n_particles;
395✔
1239
    } // End of subiteration loop over secondary tracks
1240
    n_generation_depth++;
406✔
1241
    simulation::simulation_tracks_completed += alive_secondary;
406✔
1242
  } // End of loop over secondary generations
1243

1244
  // Reset work so that fission bank etc works correctly
1245
  calculate_work(settings::n_particles);
11✔
1246
}
11✔
1247

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