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

openmc-dev / openmc / 21034277178

15 Jan 2026 02:11PM UTC coverage: 81.871% (-0.2%) from 82.044%
21034277178

Pull #3702

github

web-flow
Merge d0e725296 into 179048b80
Pull Request #3702: Random Ray Kinetic Simulation Mode

17650 of 24558 branches covered (71.87%)

Branch coverage included in aggregate %.

1008 of 1149 new or added lines in 20 files covered. (87.73%)

166 existing lines in 15 files now uncovered.

56377 of 65861 relevant lines covered (85.6%)

47768466.83 hits per line

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

87.91
/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/material.h"
13
#include "openmc/message_passing.h"
14
#include "openmc/nuclide.h"
15
#include "openmc/output.h"
16
#include "openmc/particle.h"
17
#include "openmc/photon.h"
18
#include "openmc/random_lcg.h"
19
#include "openmc/settings.h"
20
#include "openmc/source.h"
21
#include "openmc/state_point.h"
22
#include "openmc/tallies/derivative.h"
23
#include "openmc/tallies/filter.h"
24
#include "openmc/tallies/tally.h"
25
#include "openmc/tallies/trigger.h"
26
#include "openmc/timer.h"
27
#include "openmc/track_output.h"
28
#include "openmc/weight_windows.h"
29

30
#ifdef _OPENMP
31
#include <omp.h>
32
#endif
33
#include "xtensor/xview.hpp"
34

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

39
#include <fmt/format.h>
40

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

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

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

53
int openmc_run()
4,702✔
54
{
55
  openmc::simulation::time_total.start();
4,702✔
56
  openmc_simulation_init();
4,702✔
57

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

65
  int err = 0;
4,702✔
66
  while (status == 0 && err == 0) {
104,270!
67
    err = openmc_next_batch(&status);
99,578✔
68
  }
69

70
  openmc_simulation_finalize();
4,692✔
71
  openmc::simulation::time_total.stop();
4,692✔
72
  return err;
4,692✔
73
}
74

75
int openmc_simulation_init()
6,609✔
76
{
77
  using namespace openmc;
78

79
  // Skip if simulation has already been initialized
80
  if (simulation::initialized)
6,609✔
81
    return 0;
8✔
82

83
  // Initialize nuclear data (energy limits, log grid)
84
  if (settings::run_CE) {
6,601✔
85
    initialize_data();
4,477✔
86
  }
87

88
  // Determine how much work each process should do
89
  calculate_work();
6,601✔
90

91
  // Allocate source, fission and surface source banks.
92
  allocate_banks();
6,601✔
93

94
  // Create track file if needed
95
  if (!settings::track_identifiers.empty() || settings::write_all_tracks) {
6,601✔
96
    open_track_file();
78✔
97
  }
98

99
  // If doing an event-based simulation, intialize the particle buffer
100
  // and event queues
101
  if (settings::event_based) {
6,601!
102
    int64_t event_buffer_length =
UNCOV
103
      std::min(simulation::work_per_rank, settings::max_particles_in_flight);
×
UNCOV
104
    init_event_queues(event_buffer_length);
×
105
  }
106

107
  // Allocate tally results arrays if they're not allocated yet
108
  for (auto& t : model::tallies) {
29,555✔
109
    t->set_strides();
22,954✔
110
    t->init_results();
22,954✔
111
  }
112

113
  // Set up material nuclide index mapping
114
  for (auto& mat : model::materials) {
23,349✔
115
    mat->init_nuclide_index();
16,748✔
116
  }
117

118
  // Reset global variables -- this is done before loading state point (as that
119
  // will potentially populate k_generation and entropy)
120
  simulation::current_batch = 0;
6,601✔
121
  simulation::ct_current_file = 1;
6,601✔
122
  simulation::ssw_current_file = 1;
6,601✔
123
  simulation::k_generation.clear();
6,601✔
124
  simulation::entropy.clear();
6,601✔
125
  openmc_reset();
6,601✔
126

127
  // If this is a restart run, load the state point data and binary source
128
  // file
129
  if (settings::restart_run) {
6,601✔
130
    load_state_point();
50✔
131
    write_message("Resuming simulation...", 6);
50✔
132
  } else {
133
    // Only initialize primary source bank for eigenvalue simulations
134
    if (settings::run_mode == RunMode::EIGENVALUE &&
6,551✔
135
        settings::solver_type == SolverType::MONTE_CARLO) {
4,414✔
136
      initialize_source();
3,101✔
137
    }
138
  }
139

140
  // Display header
141
  if (mpi::master) {
6,601✔
142
    if (settings::kinetic_simulation) {
5,061✔
143
      if (simulation::is_initial_condition) {
672✔
144
        if (simulation::k_eff_correction)
192✔
145
          header("KINETIC SIMULATION INITIAL CONDITION (K-EFF CORRECTION)", 3);
96✔
146
        else
147
          header("KINETIC SIMULATION INITIAL CONDITION", 3);
96✔
148
      } else {
149
        std::string message = fmt::format(
150
          "KINETIC SIMULATION TIME STEP {0}", simulation::current_timestep);
420✔
151
        const char* msg = message.c_str();
480✔
152
        header(msg, 3);
480✔
153
      }
480✔
154
    }
155
    if (settings::run_mode == RunMode::FIXED_SOURCE) {
5,061✔
156
      if (settings::solver_type == SolverType::MONTE_CARLO) {
1,833✔
157
        header("FIXED SOURCE TRANSPORT SIMULATION", 3);
1,591✔
158
      } else if (settings::solver_type == SolverType::RANDOM_RAY) {
242!
159
        header("FIXED SOURCE TRANSPORT SIMULATION (RANDOM RAY SOLVER)", 3);
242✔
160
      }
161
    } else if (settings::run_mode == RunMode::EIGENVALUE) {
3,228!
162
      if (settings::solver_type == SolverType::MONTE_CARLO) {
3,228✔
163
        header("K EIGENVALUE SIMULATION", 3);
2,420✔
164
      } else if (settings::solver_type == SolverType::RANDOM_RAY) {
808!
165
        header("K EIGENVALUE SIMULATION (RANDOM RAY SOLVER)", 3);
808✔
166
      }
167
      if (settings::verbosity >= 7)
3,228✔
168
        print_columns();
3,088✔
169
    }
170
  }
171

172
  // load weight windows from file
173
  if (!settings::weight_windows_file.empty()) {
6,601!
174
    openmc_weight_windows_import(settings::weight_windows_file.c_str());
×
175
  }
176

177
  // Set flag indicating initialization is done
178
  simulation::initialized = true;
6,601✔
179
  return 0;
6,601✔
180
}
181

182
int openmc_simulation_finalize()
6,591✔
183
{
184
  using namespace openmc;
185

186
  // Skip if simulation was never run
187
  if (!simulation::initialized)
6,591!
188
    return 0;
×
189

190
  // Stop active batch timer and start finalization timer
191
  simulation::time_active.stop();
6,591✔
192
  simulation::time_finalize.start();
6,591✔
193

194
  // Clear material nuclide mapping
195
  for (auto& mat : model::materials) {
23,329✔
196
    mat->mat_nuclide_index_.clear();
16,738✔
197
  }
198

199
  // Close track file if open
200
  if (!settings::track_identifiers.empty() || settings::write_all_tracks) {
6,591✔
201
    close_track_file();
78✔
202
  }
203

204
  // Increment total number of generations
205
  simulation::total_gen += simulation::current_batch * settings::gen_per_batch;
6,591✔
206

207
#ifdef OPENMC_MPI
208
  broadcast_results();
4,704✔
209
#endif
210

211
  // Write tally results to tallies.out
212
  if (settings::output_tallies && mpi::master)
6,591!
213
    write_tallies();
4,823✔
214

215
  // If weight window generators are present in this simulation,
216
  // write a weight windows file
217
  if (variance_reduction::weight_windows_generators.size() > 0) {
6,591✔
218
    openmc_weight_windows_export();
104✔
219
  }
220

221
  // Deactivate all tallies
222
  for (auto& t : model::tallies) {
29,545✔
223
    t->active_ = false;
22,954✔
224
  }
225

226
  // Stop timers and show timing statistics
227
  simulation::time_finalize.stop();
6,591✔
228
  simulation::time_total.stop();
6,591✔
229
  if (mpi::master) {
6,591✔
230
    if (settings::solver_type != SolverType::RANDOM_RAY) {
5,051✔
231
      if (settings::verbosity >= 6)
4,001✔
232
        print_runtime();
3,861✔
233
      if (settings::verbosity >= 4)
4,001✔
234
        print_results();
3,861✔
235
    }
236
  }
237
  if (settings::check_overlaps)
6,591!
238
    print_overlap_check();
×
239

240
  // Reset flags
241
  simulation::initialized = false;
6,591✔
242
  return 0;
6,591✔
243
}
244

245
int openmc_next_batch(int* status)
102,578✔
246
{
247
  using namespace openmc;
248
  using openmc::simulation::current_gen;
249

250
  // Make sure simulation has been initialized
251
  if (!simulation::initialized) {
102,578✔
252
    set_errmsg("Simulation has not been initialized yet.");
8✔
253
    return OPENMC_E_ALLOCATE;
8✔
254
  }
255

256
  initialize_batch();
102,570✔
257

258
  // =======================================================================
259
  // LOOP OVER GENERATIONS
260
  for (current_gen = 1; current_gen <= settings::gen_per_batch; ++current_gen) {
205,312✔
261

262
    initialize_generation();
102,752✔
263

264
    // Start timer for transport
265
    simulation::time_transport.start();
102,752✔
266

267
    // Transport loop
268
    if (settings::event_based) {
102,752!
UNCOV
269
      transport_event_based();
×
270
    } else {
271
      transport_history_based();
102,752✔
272
    }
273

274
    // Accumulate time for transport
275
    simulation::time_transport.stop();
102,742✔
276

277
    finalize_generation();
102,742✔
278
  }
279

280
  finalize_batch();
102,560✔
281

282
  // Check simulation ending criteria
283
  if (status) {
102,560!
284
    if (simulation::current_batch >= settings::n_max_batches) {
102,560✔
285
      *status = STATUS_EXIT_MAX_BATCH;
4,824✔
286
    } else if (simulation::satisfy_triggers) {
97,736✔
287
      *status = STATUS_EXIT_ON_TRIGGER;
76✔
288
    } else {
289
      *status = STATUS_EXIT_NORMAL;
97,660✔
290
    }
291
  }
292
  return 0;
102,560✔
293
}
294

295
bool openmc_is_statepoint_batch()
2,280✔
296
{
297
  using namespace openmc;
298
  using openmc::simulation::current_gen;
299

300
  if (!simulation::initialized)
2,280!
301
    return false;
×
302
  else
303
    return contains(settings::statepoint_batch, simulation::current_batch);
2,280✔
304
}
305

306
namespace openmc {
307

308
//==============================================================================
309
// Global variables
310
//==============================================================================
311

312
namespace simulation {
313

314
int ct_current_file;
315
int current_batch;
316
int current_gen;
317
bool initialized {false};
318
double keff {1.0};
319
double keff_std;
320
double k_col_abs {0.0};
321
double k_col_tra {0.0};
322
double k_abs_tra {0.0};
323
double log_spacing;
324
int n_lost_particles {0};
325
bool need_depletion_rx {false};
326
int restart_batch;
327
bool satisfy_triggers {false};
328
int ssw_current_file;
329
int total_gen {0};
330
double total_weight;
331
int64_t work_per_rank;
332

333
const RegularMesh* entropy_mesh {nullptr};
334
const RegularMesh* ufs_mesh {nullptr};
335

336
vector<double> k_generation;
337
vector<int64_t> work_index;
338

339
bool is_initial_condition {true};
340
int current_timestep;
341
double current_time;
342
bool k_eff_correction {false};
343

344
} // namespace simulation
345

346
//==============================================================================
347
// Non-member functions
348
//==============================================================================
349

350
void allocate_banks()
6,601✔
351
{
352
  if (settings::run_mode == RunMode::EIGENVALUE &&
6,601✔
353
      settings::solver_type == SolverType::MONTE_CARLO) {
4,464✔
354
    // Allocate source bank
355
    simulation::source_bank.resize(simulation::work_per_rank);
3,151✔
356

357
    // Allocate fission bank
358
    init_fission_bank(3 * simulation::work_per_rank);
3,151✔
359

360
    // Allocate IFP bank
361
    if (settings::ifp_on) {
3,151✔
362
      resize_simulation_ifp_banks();
58✔
363
    }
364
  }
365

366
  if (settings::surf_source_write) {
6,601✔
367
    // Allocate surface source bank
368
    simulation::surf_source_bank.reserve(settings::ssw_max_particles);
720✔
369
  }
370

371
  if (settings::collision_track) {
6,601✔
372
    // Allocate collision track bank
373
    collision_track_reserve_bank();
123✔
374
  }
375
}
6,601✔
376

377
void initialize_batch()
271,842✔
378
{
379
  // Increment current batch
380
  ++simulation::current_batch;
271,842✔
381
  if (settings::run_mode == RunMode::FIXED_SOURCE) {
271,842✔
382
    if (settings::solver_type == SolverType::RANDOM_RAY &&
43,735✔
383
        simulation::current_batch < settings::n_inactive + 1) {
11,452✔
384
      write_message(
6,961✔
385
        6, "Simulating batch {:<4} (inactive)", simulation::current_batch);
386
    } else {
387
      write_message(6, "Simulating batch {}", simulation::current_batch);
36,774✔
388
    }
389
  }
390

391
  // Reset total starting particle weight used for normalizing tallies
392
  simulation::total_weight = 0.0;
271,842✔
393

394
  // Determine if this batch is the first inactive or active batch.
395
  bool first_inactive = false;
271,842✔
396
  bool first_active = false;
271,842✔
397
  if (!settings::restart_run) {
271,842✔
398
    first_inactive = settings::n_inactive > 0 && simulation::current_batch == 1;
271,713✔
399
    first_active = simulation::current_batch == settings::n_inactive + 1;
271,713✔
400
  } else if (simulation::current_batch == simulation::restart_batch + 1) {
129✔
401
    first_inactive = simulation::restart_batch < settings::n_inactive;
42✔
402
    first_active = !first_inactive;
42✔
403
  }
404

405
  // Manage active/inactive timers and activate tallies if necessary.
406
  if (first_inactive) {
271,842✔
407
    simulation::time_inactive.start();
3,979✔
408
  } else if (first_active) {
267,863✔
409
    simulation::time_inactive.stop();
6,567✔
410
    simulation::time_active.start();
6,567✔
411
    for (auto& t : model::tallies) {
29,505✔
412
      t->active_ = true;
22,938✔
413
    }
414
  }
415

416
  // Add user tallies to active tallies list
417
  setup_active_tallies();
271,842✔
418
}
271,842✔
419

420
void finalize_batch()
271,832✔
421
{
422
  // Reduce tallies onto master process and accumulate
423
  simulation::time_tallies.start();
271,832✔
424
  accumulate_tallies();
271,832✔
425
  simulation::time_tallies.stop();
271,832✔
426

427
  // update weight windows if needed
428
  for (const auto& wwg : variance_reduction::weight_windows_generators) {
273,784✔
429
    wwg->update();
1,952✔
430
  }
431

432
  // Reset global tally results
433
  if (simulation::current_batch <= settings::n_inactive) {
271,832✔
434
    xt::view(simulation::global_tallies, xt::all()) = 0.0;
102,538✔
435
    simulation::n_realizations = 0;
102,538✔
436
  }
437

438
  // Check_triggers
439
  if (mpi::master)
271,832✔
440
    check_triggers();
189,782✔
441
#ifdef OPENMC_MPI
442
  MPI_Bcast(&simulation::satisfy_triggers, 1, MPI_C_BOOL, 0, mpi::intracomm);
200,756✔
443
#endif
444
  if (simulation::satisfy_triggers ||
271,832✔
445
      (settings::trigger_on &&
2,030✔
446
        simulation::current_batch == settings::n_max_batches)) {
2,030✔
447
    settings::statepoint_batch.insert(simulation::current_batch);
113✔
448
  }
449

450
  // Write out state point if it's been specified for this batch and is not
451
  // a CMFD run instance
452
  if (contains(settings::statepoint_batch, simulation::current_batch) &&
278,607✔
453
      !settings::cmfd_run) {
6,775✔
454
    if (contains(settings::sourcepoint_batch, simulation::current_batch) &&
13,121✔
455
        settings::source_write && !settings::source_separate) {
13,121✔
456
      bool b = (settings::run_mode == RunMode::EIGENVALUE);
4,709✔
457
      openmc_statepoint_write(nullptr, &b);
4,709✔
458
    } else {
459
      bool b = false;
1,938✔
460
      openmc_statepoint_write(nullptr, &b);
1,938✔
461
    }
462
  }
463

464
  if (settings::run_mode == RunMode::EIGENVALUE) {
271,832✔
465
    // Write out a separate source point if it's been specified for this batch
466
    if (contains(settings::sourcepoint_batch, simulation::current_batch) &&
232,608✔
467
        settings::source_write && settings::source_separate) {
232,608✔
468

469
      // Determine width for zero padding
470
      int w = std::to_string(settings::n_max_batches).size();
60✔
471
      std::string source_point_filename = fmt::format("{0}source.{1:0{2}}",
472
        settings::path_output, simulation::current_batch, w);
51✔
473
      span<SourceSite> bankspan(simulation::source_bank);
60✔
474
      write_source_point(source_point_filename, bankspan,
60✔
475
        simulation::work_index, settings::source_mcpl_write);
476
    }
60✔
477

478
    // Write a continously-overwritten source point if requested.
479
    if (settings::source_latest) {
228,107✔
480
      auto filename = settings::path_output + "source";
130✔
481
      span<SourceSite> bankspan(simulation::source_bank);
130✔
482
      write_source_point(filename, bankspan, simulation::work_index,
130✔
483
        settings::source_mcpl_write);
484
    }
130✔
485
  }
486

487
  // Write out surface source if requested.
488
  if (settings::surf_source_write &&
271,832✔
489
      simulation::ssw_current_file <= settings::ssw_max_files) {
6,772✔
490
    bool last_batch = (simulation::current_batch == settings::n_batches);
1,372✔
491
    if (simulation::surf_source_bank.full() || last_batch) {
1,372✔
492
      // Determine appropriate filename
493
      auto filename = fmt::format("{}surface_source.{}", settings::path_output,
494
        simulation::current_batch);
654✔
495
      if (settings::ssw_max_files == 1 ||
744✔
496
          (simulation::ssw_current_file == 1 && last_batch)) {
40!
497
        filename = settings::path_output + "surface_source";
704✔
498
      }
499

500
      // Get span of source bank and calculate parallel index vector
501
      auto surf_work_index = mpi::calculate_parallel_index_vector(
502
        simulation::surf_source_bank.size());
744✔
503
      span<SourceSite> surfbankspan(simulation::surf_source_bank.begin(),
504
        simulation::surf_source_bank.size());
744✔
505

506
      // Write surface source file
507
      write_source_point(
744✔
508
        filename, surfbankspan, surf_work_index, settings::surf_mcpl_write);
509

510
      // Reset surface source bank and increment counter
511
      simulation::surf_source_bank.clear();
744✔
512
      if (!last_batch && settings::ssw_max_files >= 1) {
744!
513
        simulation::surf_source_bank.reserve(settings::ssw_max_particles);
603✔
514
      }
515
      ++simulation::ssw_current_file;
744✔
516
    }
744✔
517
  }
518
  // Write collision track file if requested
519
  if (settings::collision_track) {
271,832✔
520
    collision_track_flush_bank();
519✔
521
  }
522
}
271,832✔
523

524
void initialize_generation()
272,024✔
525
{
526
  if (settings::run_mode == RunMode::EIGENVALUE) {
272,024✔
527
    // Clear out the fission bank
528
    simulation::fission_bank.resize(0);
228,289✔
529

530
    // Count source sites if using uniform fission source weighting
531
    if (settings::ufs_on)
228,289✔
532
      ufs_count_sites();
130✔
533

534
    // Store current value of tracklength k
535
    simulation::keff_generation = simulation::global_tallies(
228,289✔
536
      GlobalTally::K_TRACKLENGTH, TallyResult::VALUE);
537
  }
538
}
272,024✔
539

540
void finalize_generation()
272,014✔
541
{
542
  auto& gt = simulation::global_tallies;
272,014✔
543

544
  // Update global tallies with the accumulation variables
545
  if (settings::run_mode == RunMode::EIGENVALUE) {
272,014✔
546
    gt(GlobalTally::K_COLLISION, TallyResult::VALUE) += global_tally_collision;
228,289✔
547
    gt(GlobalTally::K_ABSORPTION, TallyResult::VALUE) +=
228,289✔
548
      global_tally_absorption;
549
    gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) +=
228,289✔
550
      global_tally_tracklength;
551
  }
552
  gt(GlobalTally::LEAKAGE, TallyResult::VALUE) += global_tally_leakage;
272,014✔
553

554
  // reset tallies
555
  if (settings::run_mode == RunMode::EIGENVALUE) {
272,014✔
556
    global_tally_collision = 0.0;
228,289✔
557
    global_tally_absorption = 0.0;
228,289✔
558
    global_tally_tracklength = 0.0;
228,289✔
559
  }
560
  global_tally_leakage = 0.0;
272,014✔
561

562
  if (settings::run_mode == RunMode::EIGENVALUE &&
272,014✔
563
      settings::solver_type == SolverType::MONTE_CARLO) {
228,289✔
564
    // If using shared memory, stable sort the fission bank (by parent IDs)
565
    // so as to allow for reproducibility regardless of which order particles
566
    // are run in.
567
    sort_fission_bank();
70,469✔
568

569
    // Distribute fission bank across processors evenly
570
    synchronize_bank();
70,469✔
571
  }
572

573
  if (settings::run_mode == RunMode::EIGENVALUE) {
272,014✔
574

575
    // Calculate shannon entropy
576
    if (settings::entropy_on &&
228,289✔
577
        settings::solver_type == SolverType::MONTE_CARLO)
163,430✔
578
      shannon_entropy();
5,610✔
579

580
    // Collect results and statistics
581
    calculate_generation_keff();
228,289✔
582
    calculate_average_keff();
228,289✔
583

584
    // Write generation output
585
    if (mpi::master && settings::verbosity >= 7) {
228,289✔
586
      print_generation();
149,829✔
587
    }
588
  }
589
}
272,014✔
590

591
void initialize_history(Particle& p, int64_t index_source)
121,386,074✔
592
{
593
  // set defaults
594
  if (settings::run_mode == RunMode::EIGENVALUE) {
121,386,074✔
595
    // set defaults for eigenvalue simulations from primary bank
596
    p.from_source(&simulation::source_bank[index_source - 1]);
102,528,100✔
597
  } else if (settings::run_mode == RunMode::FIXED_SOURCE) {
18,857,974!
598
    // initialize random number seed
599
    int64_t id = (simulation::total_gen + overall_generation() - 1) *
18,857,974✔
600
                   settings::n_particles +
18,857,974✔
601
                 simulation::work_index[mpi::rank] + index_source;
18,857,974✔
602
    uint64_t seed = init_seed(id, STREAM_SOURCE);
18,857,974✔
603
    // sample from external source distribution or custom library then set
604
    auto site = sample_external_source(&seed);
18,857,974✔
605
    p.from_source(&site);
18,857,971✔
606
  }
607
  p.current_work() = index_source;
121,386,071✔
608

609
  // set identifier for particle
610
  p.id() = simulation::work_index[mpi::rank] + index_source;
121,386,071✔
611

612
  // set progeny count to zero
613
  p.n_progeny() = 0;
121,386,071✔
614

615
  // Reset particle event counter
616
  p.n_event() = 0;
121,386,071✔
617

618
  // Reset split counter
619
  p.n_split() = 0;
121,386,071✔
620

621
  // Reset weight window ratio
622
  p.ww_factor() = 0.0;
121,386,071✔
623

624
  // set particle history start weight
625
  p.wgt_born() = p.wgt();
121,386,071✔
626

627
  // Reset pulse_height_storage
628
  std::fill(p.pht_storage().begin(), p.pht_storage().end(), 0);
121,386,071✔
629

630
  // set random number seed
631
  int64_t particle_seed =
632
    (simulation::total_gen + overall_generation() - 1) * settings::n_particles +
121,386,071✔
633
    p.id();
121,386,071✔
634
  init_particle_seeds(particle_seed, p.seeds());
121,386,071✔
635

636
  // set particle trace
637
  p.trace() = false;
121,386,071✔
638
  if (simulation::current_batch == settings::trace_batch &&
242,780,142✔
639
      simulation::current_gen == settings::trace_gen &&
121,394,071!
640
      p.id() == settings::trace_particle)
8,000✔
641
    p.trace() = true;
8✔
642

643
  // Set particle track.
644
  p.write_track() = check_track_criteria(p);
121,386,071✔
645

646
  // Set the particle's initial weight window value.
647
  p.wgt_ww_born() = -1.0;
121,386,071✔
648
  apply_weight_windows(p);
121,386,071✔
649

650
  // Display message if high verbosity or trace is on
651
  if (settings::verbosity >= 9 || p.trace()) {
121,386,071!
652
    write_message("Simulating Particle {}", p.id());
8✔
653
  }
654

655
// Add particle's starting weight to count for normalizing tallies later
656
#pragma omp atomic
45,654,001✔
657
  simulation::total_weight += p.wgt();
121,386,071✔
658

659
  // Force calculation of cross-sections by setting last energy to zero
660
  if (settings::run_CE) {
121,386,071✔
661
    p.invalidate_neutron_xs();
39,914,071✔
662
  }
663

664
  // Prepare to write out particle track.
665
  if (p.write_track())
121,386,071✔
666
    add_particle_track(p);
852✔
667
}
121,386,071✔
668

669
int overall_generation()
140,692,831✔
670
{
671
  using namespace simulation;
672
  return settings::gen_per_batch * (current_batch - 1) + current_gen;
140,692,831✔
673
}
674

675
void calculate_work()
6,601✔
676
{
677
  // Determine minimum amount of particles to simulate on each processor
678
  int64_t min_work = settings::n_particles / mpi::n_procs;
6,601✔
679

680
  // Determine number of processors that have one extra particle
681
  int64_t remainder = settings::n_particles % mpi::n_procs;
6,601✔
682

683
  int64_t i_bank = 0;
6,601✔
684
  simulation::work_index.resize(mpi::n_procs + 1);
6,601✔
685
  simulation::work_index[0] = 0;
6,601✔
686
  for (int i = 0; i < mpi::n_procs; ++i) {
16,281✔
687
    // Number of particles for rank i
688
    int64_t work_i = i < remainder ? min_work + 1 : min_work;
9,680!
689

690
    // Set number of particles
691
    if (mpi::rank == i)
9,680✔
692
      simulation::work_per_rank = work_i;
6,601✔
693

694
    // Set index into source bank for rank i
695
    i_bank += work_i;
9,680✔
696
    simulation::work_index[i + 1] = i_bank;
9,680✔
697
  }
698
}
6,601✔
699

700
void initialize_data()
4,501✔
701
{
702
  // Determine minimum/maximum energy for incident neutron/photon data
703
  data::energy_max = {INFTY, INFTY, INFTY, INFTY};
4,501✔
704
  data::energy_min = {0.0, 0.0, 0.0, 0.0};
4,501✔
705

706
  for (const auto& nuc : data::nuclides) {
26,597✔
707
    if (nuc->grid_.size() >= 1) {
22,096!
708
      int neutron = static_cast<int>(ParticleType::neutron);
22,096✔
709
      data::energy_min[neutron] =
44,192✔
710
        std::max(data::energy_min[neutron], nuc->grid_[0].energy.front());
22,096✔
711
      data::energy_max[neutron] =
22,096✔
712
        std::min(data::energy_max[neutron], nuc->grid_[0].energy.back());
22,096✔
713
    }
714
  }
715

716
  if (settings::photon_transport) {
4,501✔
717
    for (const auto& elem : data::elements) {
744✔
718
      if (elem->energy_.size() >= 1) {
505!
719
        int photon = static_cast<int>(ParticleType::photon);
505✔
720
        int n = elem->energy_.size();
505✔
721
        data::energy_min[photon] =
1,010✔
722
          std::max(data::energy_min[photon], std::exp(elem->energy_(1)));
505✔
723
        data::energy_max[photon] =
1,010✔
724
          std::min(data::energy_max[photon], std::exp(elem->energy_(n - 1)));
505✔
725
      }
726
    }
727

728
    if (settings::electron_treatment == ElectronTreatment::TTB) {
239✔
729
      // Determine if minimum/maximum energy for bremsstrahlung is greater/less
730
      // than the current minimum/maximum
731
      if (data::ttb_e_grid.size() >= 1) {
215!
732
        int photon = static_cast<int>(ParticleType::photon);
215✔
733
        int electron = static_cast<int>(ParticleType::electron);
215✔
734
        int positron = static_cast<int>(ParticleType::positron);
215✔
735
        int n_e = data::ttb_e_grid.size();
215✔
736

737
        const std::vector<int> charged = {electron, positron};
215✔
738
        for (auto t : charged) {
645✔
739
          data::energy_min[t] = std::exp(data::ttb_e_grid(1));
430✔
740
          data::energy_max[t] = std::exp(data::ttb_e_grid(n_e - 1));
430✔
741
        }
742

743
        data::energy_min[photon] =
430✔
744
          std::max(data::energy_min[photon], data::energy_min[electron]);
215✔
745

746
        data::energy_max[photon] =
430✔
747
          std::min(data::energy_max[photon], data::energy_max[electron]);
215✔
748
      }
215✔
749
    }
750
  }
751

752
  // Show which nuclide results in lowest energy for neutron transport
753
  for (const auto& nuc : data::nuclides) {
5,637✔
754
    // If a nuclide is present in a material that's not used in the model, its
755
    // grid has not been allocated
756
    if (nuc->grid_.size() > 0) {
5,281!
757
      double max_E = nuc->grid_[0].energy.back();
5,281✔
758
      int neutron = static_cast<int>(ParticleType::neutron);
5,281✔
759
      if (max_E == data::energy_max[neutron]) {
5,281✔
760
        write_message(7, "Maximum neutron transport energy: {} eV for {}",
4,145✔
761
          data::energy_max[neutron], nuc->name_);
4,145✔
762
        if (mpi::master && data::energy_max[neutron] < 20.0e6) {
4,145!
763
          warning("Maximum neutron energy is below 20 MeV. This may bias "
×
764
                  "the results.");
765
        }
766
        break;
4,145✔
767
      }
768
    }
769
  }
770

771
  // Set up logarithmic grid for nuclides
772
  for (auto& nuc : data::nuclides) {
26,597✔
773
    nuc->init_grid();
22,096✔
774
  }
775
  int neutron = static_cast<int>(ParticleType::neutron);
4,501✔
776
  simulation::log_spacing =
4,501✔
777
    std::log(data::energy_max[neutron] / data::energy_min[neutron]) /
4,501✔
778
    settings::n_log_bins;
779
}
4,501✔
780

781
#ifdef OPENMC_MPI
782
void broadcast_results()
4,704✔
783
{
784
  // Broadcast tally results so that each process has access to results
785
  for (auto& t : model::tallies) {
21,739✔
786
    // Create a new datatype that consists of all values for a given filter
787
    // bin and then use that to broadcast. This is done to minimize the
788
    // chance of the 'count' argument of MPI_BCAST exceeding 2**31
789
    auto& results = t->results_;
17,035✔
790

791
    auto shape = results.shape();
17,035✔
792
    int count_per_filter = shape[1] * shape[2];
17,035✔
793
    MPI_Datatype result_block;
794
    MPI_Type_contiguous(count_per_filter, MPI_DOUBLE, &result_block);
17,035✔
795
    MPI_Type_commit(&result_block);
17,035✔
796
    MPI_Bcast(results.data(), shape[0], result_block, 0, mpi::intracomm);
17,035✔
797
    MPI_Type_free(&result_block);
17,035✔
798
  }
799

800
  // Also broadcast global tally results
801
  auto& gt = simulation::global_tallies;
4,704✔
802
  MPI_Bcast(gt.data(), gt.size(), MPI_DOUBLE, 0, mpi::intracomm);
4,704✔
803

804
  // These guys are needed so that non-master processes can calculate the
805
  // combined estimate of k-effective
806
  double temp[] {
807
    simulation::k_col_abs, simulation::k_col_tra, simulation::k_abs_tra};
4,704✔
808
  MPI_Bcast(temp, 3, MPI_DOUBLE, 0, mpi::intracomm);
4,704✔
809
  simulation::k_col_abs = temp[0];
4,704✔
810
  simulation::k_col_tra = temp[1];
4,704✔
811
  simulation::k_abs_tra = temp[2];
4,704✔
812
}
4,704✔
813

814
#endif
815

816
void free_memory_simulation()
6,435✔
817
{
818
  simulation::k_generation.clear();
6,435✔
819
  simulation::entropy.clear();
6,435✔
820
}
6,435✔
821

822
void transport_history_based_single_particle(Particle& p)
121,386,095✔
823
{
824
  while (p.alive()) {
2,147,483,647✔
825
    p.event_calculate_xs();
2,147,483,647✔
826
    if (p.alive()) {
2,147,483,647!
827
      p.event_advance();
2,147,483,647✔
828
    }
829
    if (p.alive()) {
2,147,483,647✔
830
      if (p.collision_distance() > p.boundary().distance()) {
2,147,483,647✔
831
        p.event_cross_surface();
1,039,240,164✔
832
      } else if (p.alive()) {
1,967,986,469!
833
        p.event_collide();
1,967,986,469✔
834
      }
835
    }
836
    p.event_revive_from_secondary();
2,147,483,647✔
837
  }
838
  p.event_death();
121,386,088✔
839
}
121,386,088✔
840

841
void transport_history_based()
102,752✔
842
{
843
#pragma omp parallel for schedule(runtime)
844
  for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
75,841,174✔
845
    Particle p;
75,781,088✔
846
    initialize_history(p, i_work);
75,781,088✔
847
    transport_history_based_single_particle(p);
75,781,085✔
848
  }
75,781,080✔
849
}
102,744✔
850

UNCOV
851
void transport_event_based()
×
852
{
UNCOV
853
  int64_t remaining_work = simulation::work_per_rank;
×
UNCOV
854
  int64_t source_offset = 0;
×
855

856
  // To cap the total amount of memory used to store particle object data, the
857
  // number of particles in flight at any point in time can bet set. In the case
858
  // that the maximum in flight particle count is lower than the total number
859
  // of particles that need to be run this iteration, the event-based transport
860
  // loop is executed multiple times until all particles have been completed.
UNCOV
861
  while (remaining_work > 0) {
×
862
    // Figure out # of particles to run for this subiteration
863
    int64_t n_particles =
UNCOV
864
      std::min(remaining_work, settings::max_particles_in_flight);
×
865

866
    // Initialize all particle histories for this subiteration
UNCOV
867
    process_init_events(n_particles, source_offset);
×
868

869
    // Event-based transport loop
870
    while (true) {
871
      // Determine which event kernel has the longest queue
UNCOV
872
      int64_t max = std::max({simulation::calculate_fuel_xs_queue.size(),
×
UNCOV
873
        simulation::calculate_nonfuel_xs_queue.size(),
×
UNCOV
874
        simulation::advance_particle_queue.size(),
×
UNCOV
875
        simulation::surface_crossing_queue.size(),
×
UNCOV
876
        simulation::collision_queue.size()});
×
877

878
      // Execute event with the longest queue
UNCOV
879
      if (max == 0) {
×
UNCOV
880
        break;
×
UNCOV
881
      } else if (max == simulation::calculate_fuel_xs_queue.size()) {
×
UNCOV
882
        process_calculate_xs_events(simulation::calculate_fuel_xs_queue);
×
UNCOV
883
      } else if (max == simulation::calculate_nonfuel_xs_queue.size()) {
×
UNCOV
884
        process_calculate_xs_events(simulation::calculate_nonfuel_xs_queue);
×
UNCOV
885
      } else if (max == simulation::advance_particle_queue.size()) {
×
UNCOV
886
        process_advance_particle_events();
×
UNCOV
887
      } else if (max == simulation::surface_crossing_queue.size()) {
×
UNCOV
888
        process_surface_crossing_events();
×
UNCOV
889
      } else if (max == simulation::collision_queue.size()) {
×
UNCOV
890
        process_collision_events();
×
891
      }
UNCOV
892
    }
×
893

894
    // Execute death event for all particles
UNCOV
895
    process_death_events(n_particles);
×
896

897
    // Adjust remaining work and source offset variables
UNCOV
898
    remaining_work -= n_particles;
×
UNCOV
899
    source_offset += n_particles;
×
900
  }
UNCOV
901
}
×
902

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