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

openmc-dev / openmc / 9617074623

21 Jun 2024 04:53PM UTC coverage: 84.715% (+0.03%) from 84.688%
9617074623

Pull #3042

github

web-flow
Merge 47e50a3a9 into 4bd0b09e6
Pull Request #3042: Rely on std::filesystem for file_utils

26 of 31 new or added lines in 5 files covered. (83.87%)

921 existing lines in 34 files now uncovered.

48900 of 57723 relevant lines covered (84.71%)

31496763.14 hits per line

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

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

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

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

38
#include <fmt/format.h>
39

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

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

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

52
int openmc_run()
5,171✔
53
{
54
  openmc::simulation::time_total.start();
5,171✔
55
  openmc_simulation_init();
5,171✔
56

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

64
  int err = 0;
5,171✔
65
  while (status == 0 && err == 0) {
88,117✔
66
    err = openmc_next_batch(&status);
82,956✔
67
  }
68

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

74
int openmc_simulation_init()
5,655✔
75
{
76
  using namespace openmc;
77

78
  // Skip if simulation has already been initialized
79
  if (simulation::initialized)
5,655✔
80
    return 0;
12✔
81

82
  // Initialize nuclear data (energy limits, log grid)
83
  if (settings::run_CE) {
5,643✔
84
    initialize_data();
5,031✔
85
  }
86

87
  // Determine how much work each process should do
88
  calculate_work();
5,643✔
89

90
  // Allocate source, fission and surface source banks.
91
  allocate_banks();
5,643✔
92

93
  // Create track file if needed
94
  if (!settings::track_identifiers.empty() || settings::write_all_tracks) {
5,643✔
95
    open_track_file();
102✔
96
  }
97

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

106
  // Allocate tally results arrays if they're not allocated yet
107
  for (auto& t : model::tallies) {
31,626✔
108
    t->set_strides();
25,983✔
109
    t->init_results();
25,983✔
110
  }
111

112
  // Set up material nuclide index mapping
113
  for (auto& mat : model::materials) {
24,234✔
114
    mat->init_nuclide_index();
18,591✔
115
  }
116

117
  // Reset global variables -- this is done before loading state point (as that
118
  // will potentially populate k_generation and entropy)
119
  simulation::current_batch = 0;
5,643✔
120
  simulation::k_generation.clear();
5,643✔
121
  simulation::entropy.clear();
5,643✔
122
  openmc_reset();
5,643✔
123

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

137
  // Display header
138
  if (mpi::master) {
5,643✔
139
    if (settings::run_mode == RunMode::FIXED_SOURCE) {
4,860✔
140
      if (settings::solver_type == SolverType::MONTE_CARLO) {
1,670✔
141
        header("FIXED SOURCE TRANSPORT SIMULATION", 3);
1,634✔
142
      } else if (settings::solver_type == SolverType::RANDOM_RAY) {
36✔
143
        header("FIXED SOURCE TRANSPORT SIMULATION (RANDOM RAY SOLVER)", 3);
36✔
144
      }
145
    } else if (settings::run_mode == RunMode::EIGENVALUE) {
3,190✔
146
      if (settings::solver_type == SolverType::MONTE_CARLO) {
3,190✔
147
        header("K EIGENVALUE SIMULATION", 3);
3,166✔
148
      } else if (settings::solver_type == SolverType::RANDOM_RAY) {
24✔
149
        header("K EIGENVALUE SIMULATION (RANDOM RAY SOLVER)", 3);
24✔
150
      }
151
      if (settings::verbosity >= 7)
3,190✔
152
        print_columns();
2,935✔
153
    }
154
  }
155

156
  // load weight windows from file
157
  if (!settings::weight_windows_file.empty()) {
5,643✔
158
    openmc_weight_windows_import(settings::weight_windows_file.c_str());
×
159
  }
160

161
  // Set flag indicating initialization is done
162
  simulation::initialized = true;
5,643✔
163
  return 0;
5,643✔
164
}
165

166
int openmc_simulation_finalize()
5,633✔
167
{
168
  using namespace openmc;
169

170
  // Skip if simulation was never run
171
  if (!simulation::initialized)
5,633✔
172
    return 0;
×
173

174
  // Stop active batch timer and start finalization timer
175
  simulation::time_active.stop();
5,633✔
176
  simulation::time_finalize.start();
5,633✔
177

178
  // Clear material nuclide mapping
179
  for (auto& mat : model::materials) {
24,214✔
180
    mat->mat_nuclide_index_.clear();
18,581✔
181
  }
182

183
  // Close track file if open
184
  if (!settings::track_identifiers.empty() || settings::write_all_tracks) {
5,633✔
185
    close_track_file();
102✔
186
  }
187

188
  // Increment total number of generations
189
  simulation::total_gen += simulation::current_batch * settings::gen_per_batch;
5,633✔
190

191
#ifdef OPENMC_MPI
192
  broadcast_results();
2,934✔
193
#endif
194

195
  // Write tally results to tallies.out
196
  if (settings::output_tallies && mpi::master)
5,633✔
197
    write_tallies();
4,838✔
198

199
  // If weight window generators are present in this simulation,
200
  // write a weight windows file
201
  if (variance_reduction::weight_windows_generators.size() > 0) {
5,633✔
202
    openmc_weight_windows_export();
24✔
203
  }
204

205
  // Deactivate all tallies
206
  for (auto& t : model::tallies) {
31,616✔
207
    t->active_ = false;
25,983✔
208
  }
209

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

224
  // Reset flags
225
  simulation::initialized = false;
5,633✔
226
  return 0;
5,633✔
227
}
228

229
int openmc_next_batch(int* status)
89,736✔
230
{
231
  using namespace openmc;
232
  using openmc::simulation::current_gen;
233

234
  // Make sure simulation has been initialized
235
  if (!simulation::initialized) {
89,736✔
236
    set_errmsg("Simulation has not been initialized yet.");
12✔
237
    return OPENMC_E_ALLOCATE;
12✔
238
  }
239

240
  initialize_batch();
89,724✔
241

242
  // =======================================================================
243
  // LOOP OVER GENERATIONS
244
  for (current_gen = 1; current_gen <= settings::gen_per_batch; ++current_gen) {
179,676✔
245

246
    initialize_generation();
89,962✔
247

248
    // Start timer for transport
249
    simulation::time_transport.start();
89,962✔
250

251
    // Transport loop
252
    if (settings::event_based) {
89,962✔
253
      transport_event_based();
3,016✔
254
    } else {
255
      transport_history_based();
86,946✔
256
    }
257

258
    // Accumulate time for transport
259
    simulation::time_transport.stop();
89,952✔
260

261
    finalize_generation();
89,952✔
262
  }
263

264
  finalize_batch();
89,714✔
265

266
  // Check simulation ending criteria
267
  if (status) {
89,714✔
268
    if (simulation::current_batch >= settings::n_max_batches) {
89,714✔
269
      *status = STATUS_EXIT_MAX_BATCH;
5,523✔
270
    } else if (simulation::satisfy_triggers) {
84,191✔
271
      *status = STATUS_EXIT_ON_TRIGGER;
70✔
272
    } else {
273
      *status = STATUS_EXIT_NORMAL;
84,121✔
274
    }
275
  }
276
  return 0;
89,714✔
277
}
278

279
bool openmc_is_statepoint_batch()
5,700✔
280
{
281
  using namespace openmc;
282
  using openmc::simulation::current_gen;
283

284
  if (!simulation::initialized)
5,700✔
285
    return false;
×
286
  else
287
    return contains(settings::statepoint_batch, simulation::current_batch);
5,700✔
288
}
289

290
namespace openmc {
291

292
//==============================================================================
293
// Global variables
294
//==============================================================================
295

296
namespace simulation {
297

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

315
const RegularMesh* entropy_mesh {nullptr};
316
const RegularMesh* ufs_mesh {nullptr};
317

318
vector<double> k_generation;
319
vector<int64_t> work_index;
320

321
} // namespace simulation
322

323
//==============================================================================
324
// Non-member functions
325
//==============================================================================
326

327
void allocate_banks()
5,643✔
328
{
329
  if (settings::run_mode == RunMode::EIGENVALUE &&
5,643✔
330
      settings::solver_type == SolverType::MONTE_CARLO) {
3,834✔
331
    // Allocate source bank
332
    simulation::source_bank.resize(simulation::work_per_rank);
3,800✔
333

334
    // Allocate fission bank
335
    init_fission_bank(3 * simulation::work_per_rank);
3,800✔
336
  }
337

338
  if (settings::surf_source_write) {
5,643✔
339
    // Allocate surface source bank
340
    simulation::surf_source_bank.reserve(settings::max_surface_particles);
293✔
341
  }
342
}
5,643✔
343

344
void initialize_batch()
90,574✔
345
{
346
  // Increment current batch
347
  ++simulation::current_batch;
90,574✔
348

349
  if (settings::run_mode == RunMode::FIXED_SOURCE) {
90,574✔
350
    if (settings::solver_type == SolverType::RANDOM_RAY &&
14,619✔
351
        simulation::current_batch < settings::n_inactive + 1) {
510✔
352
      write_message(
255✔
353
        6, "Simulating batch {:<4} (inactive)", simulation::current_batch);
354
    } else {
355
      write_message(6, "Simulating batch {}", simulation::current_batch);
14,364✔
356
    }
357
  }
358

359
  // Reset total starting particle weight used for normalizing tallies
360
  simulation::total_weight = 0.0;
90,574✔
361

362
  // Determine if this batch is the first inactive or active batch.
363
  bool first_inactive = false;
90,574✔
364
  bool first_active = false;
90,574✔
365
  if (!settings::restart_run) {
90,574✔
366
    first_inactive = settings::n_inactive > 0 && simulation::current_batch == 1;
90,285✔
367
    first_active = simulation::current_batch == settings::n_inactive + 1;
90,285✔
368
  } else if (simulation::current_batch == simulation::restart_batch + 1) {
289✔
369
    first_inactive = simulation::restart_batch < settings::n_inactive;
66✔
370
    first_active = !first_inactive;
66✔
371
  }
372

373
  // Manage active/inactive timers and activate tallies if necessary.
374
  if (first_inactive) {
90,574✔
375
    simulation::time_inactive.start();
2,792✔
376
  } else if (first_active) {
87,782✔
377
    simulation::time_inactive.stop();
5,616✔
378
    simulation::time_active.start();
5,616✔
379
    for (auto& t : model::tallies) {
31,575✔
380
      t->active_ = true;
25,959✔
381
    }
382
  }
383

384
  // Add user tallies to active tallies list
385
  setup_active_tallies();
90,574✔
386
}
90,574✔
387

388
void finalize_batch()
90,564✔
389
{
390
  // Reduce tallies onto master process and accumulate
391
  simulation::time_tallies.start();
90,564✔
392
  accumulate_tallies();
90,564✔
393
  simulation::time_tallies.stop();
90,564✔
394

395
  // update weight windows if needed
396
  for (const auto& wwg : variance_reduction::weight_windows_generators) {
90,684✔
397
    wwg->update();
120✔
398
  }
399

400
  // Reset global tally results
401
  if (simulation::current_batch <= settings::n_inactive) {
90,564✔
402
    xt::view(simulation::global_tallies, xt::all()) = 0.0;
14,548✔
403
    simulation::n_realizations = 0;
14,548✔
404
  }
405

406
  // Check_triggers
407
  if (mpi::master)
90,564✔
408
    check_triggers();
74,204✔
409
#ifdef OPENMC_MPI
410
  MPI_Bcast(&simulation::satisfy_triggers, 1, MPI_C_BOOL, 0, mpi::intracomm);
49,113✔
411
#endif
412
  if (simulation::satisfy_triggers ||
90,564✔
413
      (settings::trigger_on &&
2,864✔
414
        simulation::current_batch == settings::n_max_batches)) {
2,864✔
415
    settings::statepoint_batch.insert(simulation::current_batch);
145✔
416
  }
417

418
  // Write out state point if it's been specified for this batch and is not
419
  // a CMFD run instance
420
  if (contains(settings::statepoint_batch, simulation::current_batch) &&
96,478✔
421
      !settings::cmfd_run) {
5,914✔
422
    if (contains(settings::sourcepoint_batch, simulation::current_batch) &&
10,963✔
423
        settings::source_write && !settings::source_separate) {
10,963✔
424
      bool b = (settings::run_mode == RunMode::EIGENVALUE);
5,216✔
425
      openmc_statepoint_write(nullptr, &b);
5,216✔
426
    } else {
427
      bool b = false;
378✔
428
      openmc_statepoint_write(nullptr, &b);
378✔
429
    }
430
  }
431

432
  if (settings::run_mode == RunMode::EIGENVALUE) {
90,564✔
433
    // Write out a separate source point if it's been specified for this batch
434
    if (contains(settings::sourcepoint_batch, simulation::current_batch) &&
79,860✔
435
        settings::source_write && settings::source_separate) {
79,860✔
436

437
      // Determine width for zero padding
438
      int w = std::to_string(settings::n_max_batches).size();
68✔
439
      std::string source_point_filename = fmt::format("{0}source.{1:0{2}}",
440
        settings::path_output, simulation::current_batch, w);
68✔
441
      gsl::span<SourceSite> bankspan(simulation::source_bank);
68✔
442
      if (settings::source_mcpl_write) {
68✔
443
        source_point_filename.append(".mcpl");
17✔
444
        write_mcpl_source_point(
17✔
445
          source_point_filename.c_str(), bankspan, simulation::work_index);
446
      } else {
447
        source_point_filename.append(".h5");
51✔
448
        write_source_point(
51✔
449
          source_point_filename.c_str(), bankspan, simulation::work_index);
450
      }
451
    }
68✔
452

453
    // Write a continously-overwritten source point if requested.
454
    if (settings::source_latest) {
75,955✔
455

456
      auto filename = settings::path_output + "source";
170✔
457
      gsl::span<SourceSite> bankspan(simulation::source_bank);
170✔
458
      if (settings::source_mcpl_write) {
170✔
NEW
459
        filename.append(".mcpl");
×
UNCOV
460
        write_mcpl_source_point(
×
461
          filename.c_str(), bankspan, simulation::work_index);
462
      } else {
463
        filename.append(".h5");
170✔
464
        write_source_point(filename.c_str(), bankspan, simulation::work_index);
170✔
465
      }
466
    }
170✔
467
  }
468

469
  // Write out surface source if requested.
470
  if (settings::surf_source_write &&
90,564✔
471
      simulation::current_batch == settings::n_batches) {
1,477✔
472
    auto filename = settings::path_output + "surface_source";
293✔
473
    auto surf_work_index =
474
      mpi::calculate_parallel_index_vector(simulation::surf_source_bank.size());
293✔
475
    gsl::span<SourceSite> surfbankspan(simulation::surf_source_bank.begin(),
476
      simulation::surf_source_bank.size());
293✔
477
    if (settings::surf_mcpl_write) {
293✔
NEW
478
      filename.append(".mcpl");
×
UNCOV
479
      write_mcpl_source_point(filename.c_str(), surfbankspan, surf_work_index);
×
480
    } else {
481
      filename.append(".h5");
293✔
482
      write_source_point(filename.c_str(), surfbankspan, surf_work_index);
293✔
483
    }
484
  }
293✔
485
}
90,564✔
486

487
void initialize_generation()
90,812✔
488
{
489
  if (settings::run_mode == RunMode::EIGENVALUE) {
90,812✔
490
    // Clear out the fission bank
491
    simulation::fission_bank.resize(0);
76,193✔
492

493
    // Count source sites if using uniform fission source weighting
494
    if (settings::ufs_on)
76,193✔
495
      ufs_count_sites();
170✔
496

497
    // Store current value of tracklength k
498
    simulation::keff_generation = simulation::global_tallies(
76,193✔
499
      GlobalTally::K_TRACKLENGTH, TallyResult::VALUE);
500
  }
501
}
90,812✔
502

503
void finalize_generation()
90,802✔
504
{
505
  auto& gt = simulation::global_tallies;
90,802✔
506

507
  // Update global tallies with the accumulation variables
508
  if (settings::run_mode == RunMode::EIGENVALUE) {
90,802✔
509
    gt(GlobalTally::K_COLLISION, TallyResult::VALUE) += global_tally_collision;
76,193✔
510
    gt(GlobalTally::K_ABSORPTION, TallyResult::VALUE) +=
76,193✔
511
      global_tally_absorption;
512
    gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) +=
76,193✔
513
      global_tally_tracklength;
514
  }
515
  gt(GlobalTally::LEAKAGE, TallyResult::VALUE) += global_tally_leakage;
90,802✔
516

517
  // reset tallies
518
  if (settings::run_mode == RunMode::EIGENVALUE) {
90,802✔
519
    global_tally_collision = 0.0;
76,193✔
520
    global_tally_absorption = 0.0;
76,193✔
521
    global_tally_tracklength = 0.0;
76,193✔
522
  }
523
  global_tally_leakage = 0.0;
90,802✔
524

525
  if (settings::run_mode == RunMode::EIGENVALUE &&
90,802✔
526
      settings::solver_type == SolverType::MONTE_CARLO) {
76,193✔
527
    // If using shared memory, stable sort the fission bank (by parent IDs)
528
    // so as to allow for reproducibility regardless of which order particles
529
    // are run in.
530
    sort_fission_bank();
75,853✔
531

532
    // Distribute fission bank across processors evenly
533
    synchronize_bank();
75,853✔
534
  }
535

536
  if (settings::run_mode == RunMode::EIGENVALUE) {
90,802✔
537

538
    // Calculate shannon entropy
539
    if (settings::entropy_on)
76,193✔
540
      shannon_entropy();
10,670✔
541

542
    // Collect results and statistics
543
    calculate_generation_keff();
76,193✔
544
    calculate_average_keff();
76,193✔
545

546
    // Write generation output
547
    if (mpi::master && settings::verbosity >= 7) {
76,193✔
548
      print_generation();
57,472✔
549
    }
550
  }
551
}
90,802✔
552

553
void initialize_history(Particle& p, int64_t index_source)
159,205,542✔
554
{
555
  // set defaults
556
  if (settings::run_mode == RunMode::EIGENVALUE) {
159,205,542✔
557
    // set defaults for eigenvalue simulations from primary bank
558
    p.from_source(&simulation::source_bank[index_source - 1]);
146,840,600✔
559
  } else if (settings::run_mode == RunMode::FIXED_SOURCE) {
12,364,942✔
560
    // initialize random number seed
561
    int64_t id = (simulation::total_gen + overall_generation() - 1) *
12,364,942✔
562
                   settings::n_particles +
12,364,942✔
563
                 simulation::work_index[mpi::rank] + index_source;
12,364,942✔
564
    uint64_t seed = init_seed(id, STREAM_SOURCE);
12,364,942✔
565
    // sample from external source distribution or custom library then set
566
    auto site = sample_external_source(&seed);
12,364,942✔
567
    p.from_source(&site);
12,364,938✔
568
  }
569
  p.current_work() = index_source;
159,205,538✔
570

571
  // set identifier for particle
572
  p.id() = simulation::work_index[mpi::rank] + index_source;
159,205,538✔
573

574
  // set progeny count to zero
575
  p.n_progeny() = 0;
159,205,538✔
576

577
  // Reset particle event counter
578
  p.n_event() = 0;
159,205,538✔
579

580
  // Reset split counter
581
  p.n_split() = 0;
159,205,538✔
582

583
  // Reset weight window ratio
584
  p.ww_factor() = 0.0;
159,205,538✔
585

586
  // Reset pulse_height_storage
587
  std::fill(p.pht_storage().begin(), p.pht_storage().end(), 0);
159,205,538✔
588

589
  // set random number seed
590
  int64_t particle_seed =
591
    (simulation::total_gen + overall_generation() - 1) * settings::n_particles +
159,205,538✔
592
    p.id();
159,205,538✔
593
  init_particle_seeds(particle_seed, p.seeds());
159,205,538✔
594

595
  // set particle trace
596
  p.trace() = false;
159,205,538✔
597
  if (simulation::current_batch == settings::trace_batch &&
318,423,076✔
598
      simulation::current_gen == settings::trace_gen &&
159,217,538✔
599
      p.id() == settings::trace_particle)
12,000✔
600
    p.trace() = true;
12✔
601

602
  // Set particle track.
603
  p.write_track() = check_track_criteria(p);
159,205,538✔
604

605
  // Display message if high verbosity or trace is on
606
  if (settings::verbosity >= 9 || p.trace()) {
159,205,538✔
607
    write_message("Simulating Particle {}", p.id());
12✔
608
  }
609

610
// Add paricle's starting weight to count for normalizing tallies later
611
#pragma omp atomic
79,865,650✔
612
  simulation::total_weight += p.wgt();
159,205,538✔
613

614
  // Force calculation of cross-sections by setting last energy to zero
615
  if (settings::run_CE) {
159,205,538✔
616
    p.invalidate_neutron_xs();
38,197,538✔
617
  }
618

619
  // Prepare to write out particle track.
620
  if (p.write_track())
159,205,538✔
621
    add_particle_track(p);
1,128✔
622
}
159,205,538✔
623

624
int overall_generation()
171,780,239✔
625
{
626
  using namespace simulation;
627
  return settings::gen_per_batch * (current_batch - 1) + current_gen;
171,780,239✔
628
}
629

630
void calculate_work()
5,643✔
631
{
632
  // Determine minimum amount of particles to simulate on each processor
633
  int64_t min_work = settings::n_particles / mpi::n_procs;
5,643✔
634

635
  // Determine number of processors that have one extra particle
636
  int64_t remainder = settings::n_particles % mpi::n_procs;
5,643✔
637

638
  int64_t i_bank = 0;
5,643✔
639
  simulation::work_index.resize(mpi::n_procs + 1);
5,643✔
640
  simulation::work_index[0] = 0;
5,643✔
641
  for (int i = 0; i < mpi::n_procs; ++i) {
12,851✔
642
    // Number of particles for rank i
643
    int64_t work_i = i < remainder ? min_work + 1 : min_work;
7,208✔
644

645
    // Set number of particles
646
    if (mpi::rank == i)
7,208✔
647
      simulation::work_per_rank = work_i;
5,643✔
648

649
    // Set index into source bank for rank i
650
    i_bank += work_i;
7,208✔
651
    simulation::work_index[i + 1] = i_bank;
7,208✔
652
  }
653
}
5,643✔
654

655
void initialize_data()
5,067✔
656
{
657
  // Determine minimum/maximum energy for incident neutron/photon data
658
  data::energy_max = {INFTY, INFTY};
5,067✔
659
  data::energy_min = {0.0, 0.0};
5,067✔
660
  for (const auto& nuc : data::nuclides) {
33,745✔
661
    if (nuc->grid_.size() >= 1) {
28,678✔
662
      int neutron = static_cast<int>(ParticleType::neutron);
28,678✔
663
      data::energy_min[neutron] =
57,356✔
664
        std::max(data::energy_min[neutron], nuc->grid_[0].energy.front());
28,678✔
665
      data::energy_max[neutron] =
28,678✔
666
        std::min(data::energy_max[neutron], nuc->grid_[0].energy.back());
28,678✔
667
    }
668
  }
669

670
  if (settings::photon_transport) {
5,067✔
671
    for (const auto& elem : data::elements) {
781✔
672
      if (elem->energy_.size() >= 1) {
539✔
673
        int photon = static_cast<int>(ParticleType::photon);
539✔
674
        int n = elem->energy_.size();
539✔
675
        data::energy_min[photon] =
1,078✔
676
          std::max(data::energy_min[photon], std::exp(elem->energy_(1)));
539✔
677
        data::energy_max[photon] =
1,078✔
678
          std::min(data::energy_max[photon], std::exp(elem->energy_(n - 1)));
539✔
679
      }
680
    }
681

682
    if (settings::electron_treatment == ElectronTreatment::TTB) {
242✔
683
      // Determine if minimum/maximum energy for bremsstrahlung is greater/less
684
      // than the current minimum/maximum
685
      if (data::ttb_e_grid.size() >= 1) {
242✔
686
        int photon = static_cast<int>(ParticleType::photon);
242✔
687
        int n_e = data::ttb_e_grid.size();
242✔
688
        data::energy_min[photon] =
484✔
689
          std::max(data::energy_min[photon], std::exp(data::ttb_e_grid(1)));
242✔
690
        data::energy_max[photon] = std::min(
242✔
691
          data::energy_max[photon], std::exp(data::ttb_e_grid(n_e - 1)));
484✔
692
      }
693
    }
694
  }
695

696
  // Show which nuclide results in lowest energy for neutron transport
697
  for (const auto& nuc : data::nuclides) {
6,297✔
698
    // If a nuclide is present in a material that's not used in the model, its
699
    // grid has not been allocated
700
    if (nuc->grid_.size() > 0) {
5,864✔
701
      double max_E = nuc->grid_[0].energy.back();
5,864✔
702
      int neutron = static_cast<int>(ParticleType::neutron);
5,864✔
703
      if (max_E == data::energy_max[neutron]) {
5,864✔
704
        write_message(7, "Maximum neutron transport energy: {} eV for {}",
4,634✔
705
          data::energy_max[neutron], nuc->name_);
4,634✔
706
        if (mpi::master && data::energy_max[neutron] < 20.0e6) {
4,634✔
707
          warning("Maximum neutron energy is below 20 MeV. This may bias "
×
708
                  "the results.");
709
        }
710
        break;
4,634✔
711
      }
712
    }
713
  }
714

715
  // Set up logarithmic grid for nuclides
716
  for (auto& nuc : data::nuclides) {
33,745✔
717
    nuc->init_grid();
28,678✔
718
  }
719
  int neutron = static_cast<int>(ParticleType::neutron);
5,067✔
720
  simulation::log_spacing =
5,067✔
721
    std::log(data::energy_max[neutron] / data::energy_min[neutron]) /
5,067✔
722
    settings::n_log_bins;
723
}
5,067✔
724

725
#ifdef OPENMC_MPI
726
void broadcast_results()
2,934✔
727
{
728
  // Broadcast tally results so that each process has access to results
729
  for (auto& t : model::tallies) {
17,867✔
730
    // Create a new datatype that consists of all values for a given filter
731
    // bin and then use that to broadcast. This is done to minimize the
732
    // chance of the 'count' argument of MPI_BCAST exceeding 2**31
733
    auto& results = t->results_;
14,933✔
734

735
    auto shape = results.shape();
14,933✔
736
    int count_per_filter = shape[1] * shape[2];
14,933✔
737
    MPI_Datatype result_block;
738
    MPI_Type_contiguous(count_per_filter, MPI_DOUBLE, &result_block);
14,933✔
739
    MPI_Type_commit(&result_block);
14,933✔
740
    MPI_Bcast(results.data(), shape[0], result_block, 0, mpi::intracomm);
14,933✔
741
    MPI_Type_free(&result_block);
14,933✔
742
  }
743

744
  // Also broadcast global tally results
745
  auto& gt = simulation::global_tallies;
2,934✔
746
  MPI_Bcast(gt.data(), gt.size(), MPI_DOUBLE, 0, mpi::intracomm);
2,934✔
747

748
  // These guys are needed so that non-master processes can calculate the
749
  // combined estimate of k-effective
750
  double temp[] {
751
    simulation::k_col_abs, simulation::k_col_tra, simulation::k_abs_tra};
2,934✔
752
  MPI_Bcast(temp, 3, MPI_DOUBLE, 0, mpi::intracomm);
2,934✔
753
  simulation::k_col_abs = temp[0];
2,934✔
754
  simulation::k_col_tra = temp[1];
2,934✔
755
  simulation::k_abs_tra = temp[2];
2,934✔
756
}
2,934✔
757

758
#endif
759

760
void free_memory_simulation()
6,249✔
761
{
762
  simulation::k_generation.clear();
6,249✔
763
  simulation::entropy.clear();
6,249✔
764
}
6,249✔
765

766
void transport_history_based_single_particle(Particle& p)
148,197,174✔
767
{
768
  while (p.alive()) {
2,147,483,647✔
769
    p.event_calculate_xs();
2,147,483,647✔
770
    if (p.alive()) {
2,147,483,647✔
771
      p.event_advance();
2,147,483,647✔
772
    }
773
    if (p.alive()) {
2,147,483,647✔
774
      if (p.collision_distance() > p.boundary().distance) {
2,147,483,647✔
775
        p.event_cross_surface();
1,317,901,458✔
776
      } else if (p.alive()) {
2,147,483,647✔
777
        p.event_collide();
2,147,483,647✔
778
      }
779
    }
780
    p.event_revive_from_secondary();
2,147,483,647✔
781
  }
782
  p.event_death();
148,197,168✔
783
}
148,197,168✔
784

785
void transport_history_based()
86,946✔
786
{
787
#pragma omp parallel for schedule(runtime)
788
  for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
79,482,928✔
789
    Particle p;
79,440,462✔
790
    initialize_history(p, i_work);
79,440,462✔
791
    transport_history_based_single_particle(p);
79,440,458✔
792
  }
79,440,452✔
793
}
86,936✔
794

795
void transport_event_based()
3,016✔
796
{
797
  int64_t remaining_work = simulation::work_per_rank;
3,016✔
798
  int64_t source_offset = 0;
3,016✔
799

800
  // To cap the total amount of memory used to store particle object data, the
801
  // number of particles in flight at any point in time can bet set. In the case
802
  // that the maximum in flight particle count is lower than the total number
803
  // of particles that need to be run this iteration, the event-based transport
804
  // loop is executed multiple times until all particles have been completed.
805
  while (remaining_work > 0) {
6,032✔
806
    // Figure out # of particles to run for this subiteration
807
    int64_t n_particles =
808
      std::min(remaining_work, settings::max_particles_in_flight);
3,016✔
809

810
    // Initialize all particle histories for this subiteration
811
    process_init_events(n_particles, source_offset);
3,016✔
812

813
    // Event-based transport loop
814
    while (true) {
815
      // Determine which event kernel has the longest queue
816
      int64_t max = std::max({simulation::calculate_fuel_xs_queue.size(),
4,732,848✔
817
        simulation::calculate_nonfuel_xs_queue.size(),
2,366,424✔
818
        simulation::advance_particle_queue.size(),
2,366,424✔
819
        simulation::surface_crossing_queue.size(),
2,366,424✔
820
        simulation::collision_queue.size()});
2,366,424✔
821

822
      // Execute event with the longest queue
823
      if (max == 0) {
2,366,424✔
824
        break;
3,016✔
825
      } else if (max == simulation::calculate_fuel_xs_queue.size()) {
2,363,408✔
826
        process_calculate_xs_events(simulation::calculate_fuel_xs_queue);
422,106✔
827
      } else if (max == simulation::calculate_nonfuel_xs_queue.size()) {
1,941,302✔
828
        process_calculate_xs_events(simulation::calculate_nonfuel_xs_queue);
362,839✔
829
      } else if (max == simulation::advance_particle_queue.size()) {
1,578,463✔
830
        process_advance_particle_events();
780,091✔
831
      } else if (max == simulation::surface_crossing_queue.size()) {
798,372✔
832
        process_surface_crossing_events();
255,910✔
833
      } else if (max == simulation::collision_queue.size()) {
542,462✔
834
        process_collision_events();
542,462✔
835
      }
836
    }
2,363,408✔
837

838
    // Execute death event for all particles
839
    process_death_events(n_particles);
3,016✔
840

841
    // Adjust remaining work and source offset variables
842
    remaining_work -= n_particles;
3,016✔
843
    source_offset += n_particles;
3,016✔
844
  }
845
}
3,016✔
846

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