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

openmc-dev / openmc / 18875412062

28 Oct 2025 12:53PM UTC coverage: 81.805% (+0.001%) from 81.804%
18875412062

Pull #3618

github

web-flow
Merge d145e08f9 into f10d7d9f6
Pull Request #3618: Updates all `check_type('filename', filename, str)` calls to accept both `str` and `os.PathLike` objects.

16651 of 23259 branches covered (71.59%)

Branch coverage included in aggregate %.

8 of 11 new or added lines in 5 files covered. (72.73%)

89 existing lines in 2 files now uncovered.

53799 of 62860 relevant lines covered (85.59%)

43335951.24 hits per line

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

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

3
#include "openmc/bank.h"
4
#include "openmc/capi.h"
5
#include "openmc/container_util.h"
6
#include "openmc/eigenvalue.h"
7
#include "openmc/error.h"
8
#include "openmc/event.h"
9
#include "openmc/geometry_aux.h"
10
#include "openmc/ifp.h"
11
#include "openmc/material.h"
12
#include "openmc/mcpl_interface.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()
5,759✔
54
{
55
  openmc::simulation::time_total.start();
5,759✔
56
  openmc_simulation_init();
5,759✔
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;
5,759✔
61
  if (openmc::simulation::current_batch >= openmc::settings::n_max_batches) {
5,759✔
62
    status = openmc::STATUS_EXIT_MAX_BATCH;
11✔
63
  }
64

65
  int err = 0;
5,759✔
66
  while (status == 0 && err == 0) {
101,776!
67
    err = openmc_next_batch(&status);
96,029✔
68
  }
69

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

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

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

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

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

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

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

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

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

113
  // Set up material nuclide index mapping
114
  for (auto& mat : model::materials) {
25,914✔
115
    mat->init_nuclide_index();
19,176✔
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,738✔
121
  simulation::ssw_current_file = 1;
6,738✔
122
  simulation::k_generation.clear();
6,738✔
123
  simulation::entropy.clear();
6,738✔
124
  openmc_reset();
6,738✔
125

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

139
  // Display header
140
  if (mpi::master) {
6,738✔
141
    if (settings::run_mode == RunMode::FIXED_SOURCE) {
5,728✔
142
      if (settings::solver_type == SolverType::MONTE_CARLO) {
2,523✔
143
        header("FIXED SOURCE TRANSPORT SIMULATION", 3);
2,202✔
144
      } else if (settings::solver_type == SolverType::RANDOM_RAY) {
321!
145
        header("FIXED SOURCE TRANSPORT SIMULATION (RANDOM RAY SOLVER)", 3);
321✔
146
      }
147
    } else if (settings::run_mode == RunMode::EIGENVALUE) {
3,205!
148
      if (settings::solver_type == SolverType::MONTE_CARLO) {
3,205✔
149
        header("K EIGENVALUE SIMULATION", 3);
3,073✔
150
      } else if (settings::solver_type == SolverType::RANDOM_RAY) {
132!
151
        header("K EIGENVALUE SIMULATION (RANDOM RAY SOLVER)", 3);
132✔
152
      }
153
      if (settings::verbosity >= 7)
3,205✔
154
        print_columns();
2,957✔
155
    }
156
  }
157

158
  // load weight windows from file
159
  if (!settings::weight_windows_file.empty()) {
6,738!
160
    openmc_weight_windows_import(settings::weight_windows_file.c_str());
×
161
  }
162

163
  // Set flag indicating initialization is done
164
  simulation::initialized = true;
6,738✔
165
  return 0;
6,738✔
166
}
167

168
int openmc_simulation_finalize()
6,726✔
169
{
170
  using namespace openmc;
171

172
  // Skip if simulation was never run
173
  if (!simulation::initialized)
6,726!
174
    return 0;
×
175

176
  // Stop active batch timer and start finalization timer
177
  simulation::time_active.stop();
6,726✔
178
  simulation::time_finalize.start();
6,726✔
179

180
  // Clear material nuclide mapping
181
  for (auto& mat : model::materials) {
25,890✔
182
    mat->mat_nuclide_index_.clear();
19,164✔
183
  }
184

185
  // Close track file if open
186
  if (!settings::track_identifiers.empty() || settings::write_all_tracks) {
6,726✔
187
    close_track_file();
96✔
188
  }
189

190
  // Increment total number of generations
191
  simulation::total_gen += simulation::current_batch * settings::gen_per_batch;
6,726✔
192

193
#ifdef OPENMC_MPI
194
  broadcast_results();
3,880✔
195
#endif
196

197
  // Write tally results to tallies.out
198
  if (settings::output_tallies && mpi::master)
6,726!
199
    write_tallies();
5,625✔
200

201
  // If weight window generators are present in this simulation,
202
  // write a weight windows file
203
  if (variance_reduction::weight_windows_generators.size() > 0) {
6,726✔
204
    openmc_weight_windows_export();
131✔
205
  }
206

207
  // Deactivate all tallies
208
  for (auto& t : model::tallies) {
33,807✔
209
    t->active_ = false;
27,081✔
210
  }
211

212
  // Stop timers and show timing statistics
213
  simulation::time_finalize.stop();
6,726✔
214
  simulation::time_total.stop();
6,726✔
215
  if (mpi::master) {
6,726✔
216
    if (settings::solver_type != SolverType::RANDOM_RAY) {
5,716✔
217
      if (settings::verbosity >= 6)
5,263✔
218
        print_runtime();
5,015✔
219
      if (settings::verbosity >= 4)
5,263✔
220
        print_results();
5,015✔
221
    }
222
  }
223
  if (settings::check_overlaps)
6,726!
224
    print_overlap_check();
×
225

226
  // Reset flags
227
  simulation::initialized = false;
6,726✔
228
  return 0;
6,726✔
229
}
230

231
int openmc_next_batch(int* status)
100,994✔
232
{
233
  using namespace openmc;
234
  using openmc::simulation::current_gen;
235

236
  // Make sure simulation has been initialized
237
  if (!simulation::initialized) {
100,994✔
238
    set_errmsg("Simulation has not been initialized yet.");
14✔
239
    return OPENMC_E_ALLOCATE;
14✔
240
  }
241

242
  initialize_batch();
100,980✔
243

244
  // =======================================================================
245
  // LOOP OVER GENERATIONS
246
  for (current_gen = 1; current_gen <= settings::gen_per_batch; ++current_gen) {
202,172✔
247

248
    initialize_generation();
101,204✔
249

250
    // Start timer for transport
251
    simulation::time_transport.start();
101,204✔
252

253
    // Transport loop
254
    if (settings::event_based) {
101,204✔
255
      transport_event_based();
3,102✔
256
    } else {
257
      transport_history_based();
98,102✔
258
    }
259

260
    // Accumulate time for transport
261
    simulation::time_transport.stop();
101,192✔
262

263
    finalize_generation();
101,192✔
264
  }
265

266
  finalize_batch();
100,968✔
267

268
  // Check simulation ending criteria
269
  if (status) {
100,968!
270
    if (simulation::current_batch >= settings::n_max_batches) {
100,968✔
271
      *status = STATUS_EXIT_MAX_BATCH;
6,002✔
272
    } else if (simulation::satisfy_triggers) {
94,966✔
273
      *status = STATUS_EXIT_ON_TRIGGER;
97✔
274
    } else {
275
      *status = STATUS_EXIT_NORMAL;
94,869✔
276
    }
277
  }
278
  return 0;
100,968✔
279
}
280

281
bool openmc_is_statepoint_batch()
3,705✔
282
{
283
  using namespace openmc;
284
  using openmc::simulation::current_gen;
285

286
  if (!simulation::initialized)
3,705!
287
    return false;
×
288
  else
289
    return contains(settings::statepoint_batch, simulation::current_batch);
3,705✔
290
}
291

292
namespace openmc {
293

294
//==============================================================================
295
// Global variables
296
//==============================================================================
297

298
namespace simulation {
299

300
int current_batch;
301
int current_gen;
302
bool initialized {false};
303
double keff {1.0};
304
double keff_std;
305
double k_col_abs {0.0};
306
double k_col_tra {0.0};
307
double k_abs_tra {0.0};
308
double log_spacing;
309
int n_lost_particles {0};
310
bool need_depletion_rx {false};
311
int restart_batch;
312
bool satisfy_triggers {false};
313
int ssw_current_file;
314
int total_gen {0};
315
double total_weight;
316
int64_t work_per_rank;
317

318
const RegularMesh* entropy_mesh {nullptr};
319
const RegularMesh* ufs_mesh {nullptr};
320

321
vector<double> k_generation;
322
vector<int64_t> work_index;
323

324
} // namespace simulation
325

326
//==============================================================================
327
// Non-member functions
328
//==============================================================================
329

330
void allocate_banks()
6,738✔
331
{
332
  if (settings::run_mode == RunMode::EIGENVALUE &&
6,738✔
333
      settings::solver_type == SolverType::MONTE_CARLO) {
3,931✔
334
    // Allocate source bank
335
    simulation::source_bank.resize(simulation::work_per_rank);
3,739✔
336

337
    // Allocate fission bank
338
    init_fission_bank(3 * simulation::work_per_rank);
3,739✔
339

340
    // Allocate IFP bank
341
    if (settings::ifp_on) {
3,739✔
342
      resize_simulation_ifp_banks();
76✔
343
    }
344
  }
345

346
  if (settings::surf_source_write) {
6,738✔
347
    // Allocate surface source bank
348
    simulation::surf_source_bank.reserve(settings::ssw_max_particles);
1,009✔
349
  }
350
}
6,738✔
351

352
void initialize_batch()
117,952✔
353
{
354
  // Increment current batch
355
  ++simulation::current_batch;
117,952✔
356
  if (settings::run_mode == RunMode::FIXED_SOURCE) {
117,952✔
357
    if (settings::solver_type == SolverType::RANDOM_RAY &&
35,770✔
358
        simulation::current_batch < settings::n_inactive + 1) {
13,932✔
359
      write_message(
8,486✔
360
        6, "Simulating batch {:<4} (inactive)", simulation::current_batch);
361
    } else {
362
      write_message(6, "Simulating batch {}", simulation::current_batch);
27,284✔
363
    }
364
  }
365

366
  // Reset total starting particle weight used for normalizing tallies
367
  simulation::total_weight = 0.0;
117,952✔
368

369
  // Determine if this batch is the first inactive or active batch.
370
  bool first_inactive = false;
117,952✔
371
  bool first_active = false;
117,952✔
372
  if (!settings::restart_run) {
117,952✔
373
    first_inactive = settings::n_inactive > 0 && simulation::current_batch == 1;
117,765✔
374
    first_active = simulation::current_batch == settings::n_inactive + 1;
117,765✔
375
  } else if (simulation::current_batch == simulation::restart_batch + 1) {
187✔
376
    first_inactive = simulation::restart_batch < settings::n_inactive;
59✔
377
    first_active = !first_inactive;
59✔
378
  }
379

380
  // Manage active/inactive timers and activate tallies if necessary.
381
  if (first_inactive) {
117,952✔
382
    simulation::time_inactive.start();
3,364✔
383
  } else if (first_active) {
114,588✔
384
    simulation::time_inactive.stop();
6,685✔
385
    simulation::time_active.start();
6,685✔
386
    for (auto& t : model::tallies) {
33,744✔
387
      t->active_ = true;
27,059✔
388
    }
389
  }
390

391
  // Add user tallies to active tallies list
392
  setup_active_tallies();
117,952✔
393
}
117,952✔
394

395
void finalize_batch()
117,940✔
396
{
397
  // Reduce tallies onto master process and accumulate
398
  simulation::time_tallies.start();
117,940✔
399
  accumulate_tallies();
117,940✔
400
  simulation::time_tallies.stop();
117,940✔
401

402
  // update weight windows if needed
403
  for (const auto& wwg : variance_reduction::weight_windows_generators) {
120,357✔
404
    wwg->update();
2,417✔
405
  }
406

407
  // Reset global tally results
408
  if (simulation::current_batch <= settings::n_inactive) {
117,940✔
409
    xt::view(simulation::global_tallies, xt::all()) = 0.0;
26,834✔
410
    simulation::n_realizations = 0;
26,834✔
411
  }
412

413
  // Check_triggers
414
  if (mpi::master)
117,940✔
415
    check_triggers();
96,140✔
416
#ifdef OPENMC_MPI
417
  MPI_Bcast(&simulation::satisfy_triggers, 1, MPI_C_BOOL, 0, mpi::intracomm);
68,719✔
418
#endif
419
  if (simulation::satisfy_triggers ||
117,940✔
420
      (settings::trigger_on &&
2,705✔
421
        simulation::current_batch == settings::n_max_batches)) {
2,705✔
422
    settings::statepoint_batch.insert(simulation::current_batch);
149✔
423
  }
424

425
  // Write out state point if it's been specified for this batch and is not
426
  // a CMFD run instance
427
  if (contains(settings::statepoint_batch, simulation::current_batch) &&
124,911✔
428
      !settings::cmfd_run) {
6,971✔
429
    if (contains(settings::sourcepoint_batch, simulation::current_batch) &&
13,290✔
430
        settings::source_write && !settings::source_separate) {
13,290✔
431
      bool b = (settings::run_mode == RunMode::EIGENVALUE);
5,794✔
432
      openmc_statepoint_write(nullptr, &b);
5,794✔
433
    } else {
434
      bool b = false;
969✔
435
      openmc_statepoint_write(nullptr, &b);
969✔
436
    }
437
  }
438

439
  if (settings::run_mode == RunMode::EIGENVALUE) {
117,940✔
440
    // Write out a separate source point if it's been specified for this batch
441
    if (contains(settings::sourcepoint_batch, simulation::current_batch) &&
86,164✔
442
        settings::source_write && settings::source_separate) {
86,164✔
443

444
      // Determine width for zero padding
445
      int w = std::to_string(settings::n_max_batches).size();
75✔
446
      std::string source_point_filename = fmt::format("{0}source.{1:0{2}}",
447
        settings::path_output, simulation::current_batch, w);
61✔
448
      span<SourceSite> bankspan(simulation::source_bank);
75✔
449
      write_source_point(source_point_filename, bankspan,
75✔
450
        simulation::work_index, settings::source_mcpl_write);
451
    }
75✔
452

453
    // Write a continously-overwritten source point if requested.
454
    if (settings::source_latest) {
82,182✔
455
      auto filename = settings::path_output + "source";
160✔
456
      span<SourceSite> bankspan(simulation::source_bank);
160✔
457
      write_source_point(filename, bankspan, simulation::work_index,
160✔
458
        settings::source_mcpl_write);
459
    }
160✔
460
  }
461

462
  // Write out surface source if requested.
463
  if (settings::surf_source_write &&
117,940✔
464
      simulation::ssw_current_file <= settings::ssw_max_files) {
9,889✔
465
    bool last_batch = (simulation::current_batch == settings::n_batches);
1,831✔
466
    if (simulation::surf_source_bank.full() || last_batch) {
1,831✔
467
      // Determine appropriate filename
468
      auto filename = fmt::format("{}surface_source.{}", settings::path_output,
469
        simulation::current_batch);
866✔
470
      if (settings::ssw_max_files == 1 ||
1,042✔
471
          (simulation::ssw_current_file == 1 && last_batch)) {
55!
472
        filename = settings::path_output + "surface_source";
987✔
473
      }
474

475
      // Get span of source bank and calculate parallel index vector
476
      auto surf_work_index = mpi::calculate_parallel_index_vector(
477
        simulation::surf_source_bank.size());
1,042✔
478
      span<SourceSite> surfbankspan(simulation::surf_source_bank.begin(),
479
        simulation::surf_source_bank.size());
1,042✔
480

481
      // Write surface source file
482
      write_source_point(
1,042✔
483
        filename, surfbankspan, surf_work_index, settings::surf_mcpl_write);
484

485
      // Reset surface source bank and increment counter
486
      simulation::surf_source_bank.clear();
1,042✔
487
      if (!last_batch && settings::ssw_max_files >= 1) {
1,042!
488
        simulation::surf_source_bank.reserve(settings::ssw_max_particles);
860✔
489
      }
490
      ++simulation::ssw_current_file;
1,042✔
491
    }
1,042✔
492
  }
493
}
117,940✔
494

495
void initialize_generation()
118,176✔
496
{
497
  if (settings::run_mode == RunMode::EIGENVALUE) {
118,176✔
498
    // Clear out the fission bank
499
    simulation::fission_bank.resize(0);
82,406✔
500

501
    // Count source sites if using uniform fission source weighting
502
    if (settings::ufs_on)
82,406✔
503
      ufs_count_sites();
160✔
504

505
    // Store current value of tracklength k
506
    simulation::keff_generation = simulation::global_tallies(
82,406✔
507
      GlobalTally::K_TRACKLENGTH, TallyResult::VALUE);
508
  }
509
}
118,176✔
510

511
void finalize_generation()
118,164✔
512
{
513
  auto& gt = simulation::global_tallies;
118,164✔
514

515
  // Update global tallies with the accumulation variables
516
  if (settings::run_mode == RunMode::EIGENVALUE) {
118,164✔
517
    gt(GlobalTally::K_COLLISION, TallyResult::VALUE) += global_tally_collision;
82,406✔
518
    gt(GlobalTally::K_ABSORPTION, TallyResult::VALUE) +=
82,406✔
519
      global_tally_absorption;
520
    gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) +=
82,406✔
521
      global_tally_tracklength;
522
  }
523
  gt(GlobalTally::LEAKAGE, TallyResult::VALUE) += global_tally_leakage;
118,164✔
524

525
  // reset tallies
526
  if (settings::run_mode == RunMode::EIGENVALUE) {
118,164✔
527
    global_tally_collision = 0.0;
82,406✔
528
    global_tally_absorption = 0.0;
82,406✔
529
    global_tally_tracklength = 0.0;
82,406✔
530
  }
531
  global_tally_leakage = 0.0;
118,164✔
532

533
  if (settings::run_mode == RunMode::EIGENVALUE &&
118,164✔
534
      settings::solver_type == SolverType::MONTE_CARLO) {
82,406✔
535
    // If using shared memory, stable sort the fission bank (by parent IDs)
536
    // so as to allow for reproducibility regardless of which order particles
537
    // are run in.
538
    sort_fission_bank();
79,366✔
539

540
    // Distribute fission bank across processors evenly
541
    synchronize_bank();
79,366✔
542
  }
543

544
  if (settings::run_mode == RunMode::EIGENVALUE) {
118,164✔
545

546
    // Calculate shannon entropy
547
    if (settings::entropy_on &&
82,406✔
548
        settings::solver_type == SolverType::MONTE_CARLO)
11,905✔
549
      shannon_entropy();
8,865✔
550

551
    // Collect results and statistics
552
    calculate_generation_keff();
82,406✔
553
    calculate_average_keff();
82,406✔
554

555
    // Write generation output
556
    if (mpi::master && settings::verbosity >= 7) {
82,406✔
557
      print_generation();
62,619✔
558
    }
559
  }
560
}
118,164✔
561

562
void initialize_history(Particle& p, int64_t index_source)
162,146,497✔
563
{
564
  // set defaults
565
  if (settings::run_mode == RunMode::EIGENVALUE) {
162,146,497✔
566
    // set defaults for eigenvalue simulations from primary bank
567
    p.from_source(&simulation::source_bank[index_source - 1]);
135,868,300✔
568
  } else if (settings::run_mode == RunMode::FIXED_SOURCE) {
26,278,197!
569
    // initialize random number seed
570
    int64_t id = (simulation::total_gen + overall_generation() - 1) *
26,278,197✔
571
                   settings::n_particles +
26,278,197✔
572
                 simulation::work_index[mpi::rank] + index_source;
26,278,197✔
573
    uint64_t seed = init_seed(id, STREAM_SOURCE);
26,278,197✔
574
    // sample from external source distribution or custom library then set
575
    auto site = sample_external_source(&seed);
26,278,197✔
576
    p.from_source(&site);
26,278,194✔
577
  }
578
  p.current_work() = index_source;
162,146,494✔
579

580
  // set identifier for particle
581
  p.id() = simulation::work_index[mpi::rank] + index_source;
162,146,494✔
582

583
  // set progeny count to zero
584
  p.n_progeny() = 0;
162,146,494✔
585

586
  // Reset particle event counter
587
  p.n_event() = 0;
162,146,494✔
588

589
  // Reset split counter
590
  p.n_split() = 0;
162,146,494✔
591

592
  // Reset weight window ratio
593
  p.ww_factor() = 0.0;
162,146,494✔
594

595
  // set particle history start weight
596
  p.wgt_born() = p.wgt();
162,146,494✔
597

598
  // Reset pulse_height_storage
599
  std::fill(p.pht_storage().begin(), p.pht_storage().end(), 0);
162,146,494✔
600

601
  // set random number seed
602
  int64_t particle_seed =
603
    (simulation::total_gen + overall_generation() - 1) * settings::n_particles +
162,146,494✔
604
    p.id();
162,146,494✔
605
  init_particle_seeds(particle_seed, p.seeds());
162,146,494✔
606

607
  // set particle trace
608
  p.trace() = false;
162,146,494✔
609
  if (simulation::current_batch == settings::trace_batch &&
324,303,988✔
610
      simulation::current_gen == settings::trace_gen &&
162,157,494!
611
      p.id() == settings::trace_particle)
11,000✔
612
    p.trace() = true;
11✔
613

614
  // Set particle track.
615
  p.write_track() = check_track_criteria(p);
162,146,494✔
616

617
  // Set the particle's initial weight window value.
618
  p.wgt_ww_born() = -1.0;
162,146,494✔
619
  apply_weight_windows(p);
162,146,494✔
620

621
  // Display message if high verbosity or trace is on
622
  if (settings::verbosity >= 9 || p.trace()) {
162,146,494!
623
    write_message("Simulating Particle {}", p.id());
11✔
624
  }
625

626
// Add particle's starting weight to count for normalizing tallies later
627
#pragma omp atomic
88,170,928✔
628
  simulation::total_weight += p.wgt();
162,146,494✔
629

630
  // Force calculation of cross-sections by setting last energy to zero
631
  if (settings::run_CE) {
162,146,494✔
632
    p.invalidate_neutron_xs();
50,122,494✔
633
  }
634

635
  // Prepare to write out particle track.
636
  if (p.write_track())
162,146,494✔
637
    add_particle_track(p);
1,059✔
638
}
162,146,494✔
639

640
int overall_generation()
188,649,340✔
641
{
642
  using namespace simulation;
643
  return settings::gen_per_batch * (current_batch - 1) + current_gen;
188,649,340✔
644
}
645

646
void calculate_work()
6,738✔
647
{
648
  // Determine minimum amount of particles to simulate on each processor
649
  int64_t min_work = settings::n_particles / mpi::n_procs;
6,738✔
650

651
  // Determine number of processors that have one extra particle
652
  int64_t remainder = settings::n_particles % mpi::n_procs;
6,738✔
653

654
  int64_t i_bank = 0;
6,738✔
655
  simulation::work_index.resize(mpi::n_procs + 1);
6,738✔
656
  simulation::work_index[0] = 0;
6,738✔
657
  for (int i = 0; i < mpi::n_procs; ++i) {
15,495✔
658
    // Number of particles for rank i
659
    int64_t work_i = i < remainder ? min_work + 1 : min_work;
8,757!
660

661
    // Set number of particles
662
    if (mpi::rank == i)
8,757✔
663
      simulation::work_per_rank = work_i;
6,738✔
664

665
    // Set index into source bank for rank i
666
    i_bank += work_i;
8,757✔
667
    simulation::work_index[i + 1] = i_bank;
8,757✔
668
  }
669
}
6,738✔
670

671
void initialize_data()
5,595✔
672
{
673
  // Determine minimum/maximum energy for incident neutron/photon data
674
  data::energy_max = {INFTY, INFTY, INFTY, INFTY};
5,595✔
675
  data::energy_min = {0.0, 0.0, 0.0, 0.0};
5,595✔
676

677
  for (const auto& nuc : data::nuclides) {
33,120✔
678
    if (nuc->grid_.size() >= 1) {
27,525!
679
      int neutron = static_cast<int>(ParticleType::neutron);
27,525✔
680
      data::energy_min[neutron] =
55,050✔
681
        std::max(data::energy_min[neutron], nuc->grid_[0].energy.front());
27,525✔
682
      data::energy_max[neutron] =
27,525✔
683
        std::min(data::energy_max[neutron], nuc->grid_[0].energy.back());
27,525✔
684
    }
685
  }
686

687
  if (settings::photon_transport) {
5,595✔
688
    for (const auto& elem : data::elements) {
875✔
689
      if (elem->energy_.size() >= 1) {
589!
690
        int photon = static_cast<int>(ParticleType::photon);
589✔
691
        int n = elem->energy_.size();
589✔
692
        data::energy_min[photon] =
1,178✔
693
          std::max(data::energy_min[photon], std::exp(elem->energy_(1)));
589✔
694
        data::energy_max[photon] =
1,178✔
695
          std::min(data::energy_max[photon], std::exp(elem->energy_(n - 1)));
589✔
696
      }
697
    }
698

699
    if (settings::electron_treatment == ElectronTreatment::TTB) {
286✔
700
      // Determine if minimum/maximum energy for bremsstrahlung is greater/less
701
      // than the current minimum/maximum
702
      if (data::ttb_e_grid.size() >= 1) {
275!
703
        int photon = static_cast<int>(ParticleType::photon);
275✔
704
        int electron = static_cast<int>(ParticleType::electron);
275✔
705
        int positron = static_cast<int>(ParticleType::positron);
275✔
706
        int n_e = data::ttb_e_grid.size();
275✔
707

708
        const std::vector<int> charged = {electron, positron};
275✔
709
        for (auto t : charged) {
825✔
710
          data::energy_min[t] = std::exp(data::ttb_e_grid(1));
550✔
711
          data::energy_max[t] = std::exp(data::ttb_e_grid(n_e - 1));
550✔
712
        }
713

714
        data::energy_min[photon] =
550✔
715
          std::max(data::energy_min[photon], data::energy_min[electron]);
275✔
716

717
        data::energy_max[photon] =
550✔
718
          std::min(data::energy_max[photon], data::energy_max[electron]);
275✔
719
      }
275✔
720
    }
721
  }
722

723
  // Show which nuclide results in lowest energy for neutron transport
724
  for (const auto& nuc : data::nuclides) {
6,865✔
725
    // If a nuclide is present in a material that's not used in the model, its
726
    // grid has not been allocated
727
    if (nuc->grid_.size() > 0) {
6,418!
728
      double max_E = nuc->grid_[0].energy.back();
6,418✔
729
      int neutron = static_cast<int>(ParticleType::neutron);
6,418✔
730
      if (max_E == data::energy_max[neutron]) {
6,418✔
731
        write_message(7, "Maximum neutron transport energy: {} eV for {}",
5,148✔
732
          data::energy_max[neutron], nuc->name_);
5,148✔
733
        if (mpi::master && data::energy_max[neutron] < 20.0e6) {
5,148!
UNCOV
734
          warning("Maximum neutron energy is below 20 MeV. This may bias "
×
735
                  "the results.");
736
        }
737
        break;
5,148✔
738
      }
739
    }
740
  }
741

742
  // Set up logarithmic grid for nuclides
743
  for (auto& nuc : data::nuclides) {
33,120✔
744
    nuc->init_grid();
27,525✔
745
  }
746
  int neutron = static_cast<int>(ParticleType::neutron);
5,595✔
747
  simulation::log_spacing =
5,595✔
748
    std::log(data::energy_max[neutron] / data::energy_min[neutron]) /
5,595✔
749
    settings::n_log_bins;
750
}
5,595✔
751

752
#ifdef OPENMC_MPI
753
void broadcast_results()
3,880✔
754
{
755
  // Broadcast tally results so that each process has access to results
756
  for (auto& t : model::tallies) {
20,487✔
757
    // Create a new datatype that consists of all values for a given filter
758
    // bin and then use that to broadcast. This is done to minimize the
759
    // chance of the 'count' argument of MPI_BCAST exceeding 2**31
760
    auto& results = t->results_;
16,607✔
761

762
    auto shape = results.shape();
16,607✔
763
    int count_per_filter = shape[1] * shape[2];
16,607✔
764
    MPI_Datatype result_block;
765
    MPI_Type_contiguous(count_per_filter, MPI_DOUBLE, &result_block);
16,607✔
766
    MPI_Type_commit(&result_block);
16,607✔
767
    MPI_Bcast(results.data(), shape[0], result_block, 0, mpi::intracomm);
16,607✔
768
    MPI_Type_free(&result_block);
16,607✔
769
  }
770

771
  // Also broadcast global tally results
772
  auto& gt = simulation::global_tallies;
3,880✔
773
  MPI_Bcast(gt.data(), gt.size(), MPI_DOUBLE, 0, mpi::intracomm);
3,880✔
774

775
  // These guys are needed so that non-master processes can calculate the
776
  // combined estimate of k-effective
777
  double temp[] {
778
    simulation::k_col_abs, simulation::k_col_tra, simulation::k_abs_tra};
3,880✔
779
  MPI_Bcast(temp, 3, MPI_DOUBLE, 0, mpi::intracomm);
3,880✔
780
  simulation::k_col_abs = temp[0];
3,880✔
781
  simulation::k_col_tra = temp[1];
3,880✔
782
  simulation::k_abs_tra = temp[2];
3,880✔
783
}
3,880✔
784

785
#endif
786

787
void free_memory_simulation()
8,063✔
788
{
789
  simulation::k_generation.clear();
8,063✔
790
  simulation::entropy.clear();
8,063✔
791
}
8,063✔
792

793
void transport_history_based_single_particle(Particle& p)
150,048,327✔
794
{
795
  while (p.alive()) {
2,147,483,647✔
796
    p.event_calculate_xs();
2,147,483,647✔
797
    if (p.alive()) {
2,147,483,647!
798
      p.event_advance();
2,147,483,647✔
799
    }
800
    if (p.alive()) {
2,147,483,647✔
801
      if (p.collision_distance() > p.boundary().distance()) {
2,147,483,647✔
802
        p.event_cross_surface();
1,333,101,640✔
803
      } else if (p.alive()) {
2,147,483,647!
804
        p.event_collide();
2,147,483,647✔
805
      }
806
    }
807
    p.event_revive_from_secondary();
2,147,483,647✔
808
  }
809
  p.event_death();
150,048,318✔
810
}
150,048,318✔
811

812
void transport_history_based()
98,102✔
813
{
814
#pragma omp parallel for schedule(runtime)
815
  for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
74,170,477✔
816
    Particle p;
74,124,823✔
817
    initialize_history(p, i_work);
74,124,823✔
818
    transport_history_based_single_particle(p);
74,124,820✔
819
  }
74,124,815✔
820
}
98,094✔
821

822
void transport_event_based()
3,102✔
823
{
824
  int64_t remaining_work = simulation::work_per_rank;
3,102✔
825
  int64_t source_offset = 0;
3,102✔
826

827
  // To cap the total amount of memory used to store particle object data, the
828
  // number of particles in flight at any point in time can bet set. In the case
829
  // that the maximum in flight particle count is lower than the total number
830
  // of particles that need to be run this iteration, the event-based transport
831
  // loop is executed multiple times until all particles have been completed.
832
  while (remaining_work > 0) {
6,204✔
833
    // Figure out # of particles to run for this subiteration
834
    int64_t n_particles =
835
      std::min(remaining_work, settings::max_particles_in_flight);
3,102✔
836

837
    // Initialize all particle histories for this subiteration
838
    process_init_events(n_particles, source_offset);
3,102!
839

840
    // Event-based transport loop
841
    while (true) {
842
      // Determine which event kernel has the longest queue
843
      int64_t max = std::max({simulation::calculate_fuel_xs_queue.size(),
4,698,088!
844
        simulation::calculate_nonfuel_xs_queue.size(),
2,349,044✔
845
        simulation::advance_particle_queue.size(),
2,349,044✔
846
        simulation::surface_crossing_queue.size(),
2,349,044✔
847
        simulation::collision_queue.size()});
2,349,044✔
848

849
      // Execute event with the longest queue
850
      if (max == 0) {
2,349,044✔
851
        break;
3,102✔
852
      } else if (max == simulation::calculate_fuel_xs_queue.size()) {
2,345,942✔
853
        process_calculate_xs_events(simulation::calculate_fuel_xs_queue);
423,192!
854
      } else if (max == simulation::calculate_nonfuel_xs_queue.size()) {
1,922,750✔
855
        process_calculate_xs_events(simulation::calculate_nonfuel_xs_queue);
356,217!
856
      } else if (max == simulation::advance_particle_queue.size()) {
1,566,533✔
857
        process_advance_particle_events();
774,269!
858
      } else if (max == simulation::surface_crossing_queue.size()) {
792,264✔
859
        process_surface_crossing_events();
256,949!
860
      } else if (max == simulation::collision_queue.size()) {
535,315!
861
        process_collision_events();
535,315!
862
      }
863
    }
2,345,942✔
864

865
    // Execute death event for all particles
866
    process_death_events(n_particles);
3,102!
867

868
    // Adjust remaining work and source offset variables
869
    remaining_work -= n_particles;
3,102✔
870
    source_offset += n_particles;
3,102✔
871
  }
872
}
3,102✔
873

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

© 2025 Coveralls, Inc