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

openmc-dev / openmc / 9272838221

28 May 2024 04:14PM UTC coverage: 84.646% (-0.004%) from 84.65%
9272838221

push

github

web-flow
Fix `CylinderSector` and `IsogonalOctagon` translations (#3018)

Co-authored-by: Paul Romano <paul.k.romano@gmail.com>

27 of 27 new or added lines in 1 file covered. (100.0%)

4 existing lines in 2 files now uncovered.

48537 of 57341 relevant lines covered (84.65%)

34990860.16 hits per line

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

97.89
/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,576✔
53
{
54
  openmc::simulation::time_total.start();
5,576✔
55
  openmc_simulation_init();
5,576✔
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,576✔
60
  if (openmc::simulation::current_batch >= openmc::settings::n_max_batches) {
5,576✔
61
    status = openmc::STATUS_EXIT_MAX_BATCH;
14✔
62
  }
63

64
  int err = 0;
5,576✔
65
  while (status == 0 && err == 0) {
97,932✔
66
    err = openmc_next_batch(&status);
92,370✔
67
  }
68

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

74
int openmc_simulation_init()
6,164✔
75
{
76
  using namespace openmc;
77

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

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

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

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

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

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

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

112
  // Set up material nuclide index mapping
113
  for (auto& mat : model::materials) {
28,086✔
114
    mat->init_nuclide_index();
21,936✔
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;
6,150✔
120
  simulation::k_generation.clear();
6,150✔
121
  simulation::entropy.clear();
6,150✔
122
  openmc_reset();
6,150✔
123

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

137
  // Display header
138
  if (mpi::master) {
6,150✔
139
    if (settings::run_mode == RunMode::FIXED_SOURCE) {
5,387✔
140
      header("FIXED SOURCE TRANSPORT SIMULATION", 3);
1,876✔
141
    } else if (settings::run_mode == RunMode::EIGENVALUE) {
3,511✔
142
      if (settings::solver_type == SolverType::MONTE_CARLO) {
3,511✔
143
        header("K EIGENVALUE SIMULATION", 3);
3,483✔
144
      } else if (settings::solver_type == SolverType::RANDOM_RAY) {
28✔
145
        header("K EIGENVALUE SIMULATION (RANDOM RAY SOLVER)", 3);
28✔
146
      }
147
      if (settings::verbosity >= 7)
3,511✔
148
        print_columns();
3,197✔
149
    }
150
  }
151

152
  // load weight windows from file
153
  if (!settings::weight_windows_file.empty()) {
6,150✔
154
    openmc_weight_windows_import(settings::weight_windows_file.c_str());
×
155
  }
156

157
  // Set flag indicating initialization is done
158
  simulation::initialized = true;
6,150✔
159
  return 0;
6,150✔
160
}
161

162
int openmc_simulation_finalize()
6,136✔
163
{
164
  using namespace openmc;
165

166
  // Skip if simulation was never run
167
  if (!simulation::initialized)
6,136✔
168
    return 0;
×
169

170
  // Stop active batch timer and start finalization timer
171
  simulation::time_active.stop();
6,136✔
172
  simulation::time_finalize.start();
6,136✔
173

174
  // Clear material nuclide mapping
175
  for (auto& mat : model::materials) {
28,058✔
176
    mat->mat_nuclide_index_.clear();
21,922✔
177
  }
178

179
  // Close track file if open
180
  if (!settings::track_identifiers.empty() || settings::write_all_tracks) {
6,136✔
181
    close_track_file();
114✔
182
  }
183

184
  // Increment total number of generations
185
  simulation::total_gen += simulation::current_batch * settings::gen_per_batch;
6,136✔
186

187
#ifdef OPENMC_MPI
188
  broadcast_results();
2,920✔
189
#endif
190

191
  // Write tally results to tallies.out
192
  if (settings::output_tallies && mpi::master)
6,136✔
193
    write_tallies();
5,359✔
194

195
  // If weight window generators are present in this simulation,
196
  // write a weight windows file
197
  if (variance_reduction::weight_windows_generators.size() > 0) {
6,136✔
198
    openmc_weight_windows_export();
28✔
199
  }
200

201
  // Deactivate all tallies
202
  for (auto& t : model::tallies) {
35,331✔
203
    t->active_ = false;
29,195✔
204
  }
205

206
  // Stop timers and show timing statistics
207
  simulation::time_finalize.stop();
6,136✔
208
  simulation::time_total.stop();
6,136✔
209
  if (mpi::master) {
6,136✔
210
    if (settings::solver_type != SolverType::RANDOM_RAY) {
5,373✔
211
      if (settings::verbosity >= 6)
5,345✔
212
        print_runtime();
5,031✔
213
      if (settings::verbosity >= 4)
5,345✔
214
        print_results();
5,031✔
215
    }
216
  }
217
  if (settings::check_overlaps)
6,136✔
218
    print_overlap_check();
×
219

220
  // Reset flags
221
  simulation::initialized = false;
6,136✔
222
  return 0;
6,136✔
223
}
224

225
int openmc_next_batch(int* status)
101,895✔
226
{
227
  using namespace openmc;
228
  using openmc::simulation::current_gen;
229

230
  // Make sure simulation has been initialized
231
  if (!simulation::initialized) {
101,895✔
232
    set_errmsg("Simulation has not been initialized yet.");
14✔
233
    return OPENMC_E_ALLOCATE;
14✔
234
  }
235

236
  initialize_batch();
101,881✔
237

238
  // =======================================================================
239
  // LOOP OVER GENERATIONS
240
  for (current_gen = 1; current_gen <= settings::gen_per_batch; ++current_gen) {
204,014✔
241

242
    initialize_generation();
102,147✔
243

244
    // Start timer for transport
245
    simulation::time_transport.start();
102,147✔
246

247
    // Transport loop
248
    if (settings::event_based) {
102,147✔
249
      transport_event_based();
2,986✔
250
    } else {
251
      transport_history_based();
99,161✔
252
    }
253

254
    // Accumulate time for transport
255
    simulation::time_transport.stop();
102,133✔
256

257
    finalize_generation();
102,133✔
258
  }
259

260
  finalize_batch();
101,867✔
261

262
  // Check simulation ending criteria
263
  if (status) {
101,867✔
264
    if (simulation::current_batch >= settings::n_max_batches) {
101,867✔
265
      *status = STATUS_EXIT_MAX_BATCH;
6,071✔
266
    } else if (simulation::satisfy_triggers) {
95,796✔
267
      *status = STATUS_EXIT_ON_TRIGGER;
80✔
268
    } else {
269
      *status = STATUS_EXIT_NORMAL;
95,716✔
270
    }
271
  }
272
  return 0;
101,867✔
273
}
274

275
bool openmc_is_statepoint_batch()
8,265✔
276
{
277
  using namespace openmc;
278
  using openmc::simulation::current_gen;
279

280
  if (!simulation::initialized)
8,265✔
281
    return false;
×
282
  else
283
    return contains(settings::statepoint_batch, simulation::current_batch);
8,265✔
284
}
285

286
namespace openmc {
287

288
//==============================================================================
289
// Global variables
290
//==============================================================================
291

292
namespace simulation {
293

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

311
const RegularMesh* entropy_mesh {nullptr};
312
const RegularMesh* ufs_mesh {nullptr};
313

314
vector<double> k_generation;
315
vector<int64_t> work_index;
316

317
} // namespace simulation
318

319
//==============================================================================
320
// Non-member functions
321
//==============================================================================
322

323
void allocate_banks()
6,150✔
324
{
325
  if (settings::run_mode == RunMode::EIGENVALUE &&
6,150✔
326
      settings::solver_type == SolverType::MONTE_CARLO) {
4,150✔
327
    // Allocate source bank
328
    simulation::source_bank.resize(simulation::work_per_rank);
4,112✔
329

330
    // Allocate fission bank
331
    init_fission_bank(3 * simulation::work_per_rank);
4,112✔
332
  }
333

334
  if (settings::surf_source_write) {
6,150✔
335
    // Allocate surface source bank
336
    simulation::surf_source_bank.reserve(settings::max_surface_particles);
16✔
337
  }
338
}
6,150✔
339

340
void initialize_batch()
102,261✔
341
{
342
  // Increment current batch
343
  ++simulation::current_batch;
102,261✔
344

345
  if (settings::run_mode == RunMode::FIXED_SOURCE) {
102,261✔
346
    write_message(6, "Simulating batch {}", simulation::current_batch);
16,211✔
347
  }
348

349
  // Reset total starting particle weight used for normalizing tallies
350
  simulation::total_weight = 0.0;
102,261✔
351

352
  // Determine if this batch is the first inactive or active batch.
353
  bool first_inactive = false;
102,261✔
354
  bool first_active = false;
102,261✔
355
  if (!settings::restart_run) {
102,261✔
356
    first_inactive = settings::n_inactive > 0 && simulation::current_batch == 1;
101,903✔
357
    first_active = simulation::current_batch == settings::n_inactive + 1;
101,903✔
358
  } else if (simulation::current_batch == simulation::restart_batch + 1) {
358✔
359
    first_inactive = simulation::restart_batch < settings::n_inactive;
81✔
360
    first_active = !first_inactive;
81✔
361
  }
362

363
  // Manage active/inactive timers and activate tallies if necessary.
364
  if (first_inactive) {
102,261✔
365
    simulation::time_inactive.start();
2,863✔
366
  } else if (first_active) {
99,398✔
367
    simulation::time_inactive.stop();
6,119✔
368
    simulation::time_active.start();
6,119✔
369
    for (auto& t : model::tallies) {
35,286✔
370
      t->active_ = true;
29,167✔
371
    }
372
  }
373

374
  // Add user tallies to active tallies list
375
  setup_active_tallies();
102,261✔
376
}
102,261✔
377

378
void finalize_batch()
102,247✔
379
{
380
  // Reduce tallies onto master process and accumulate
381
  simulation::time_tallies.start();
102,247✔
382
  accumulate_tallies();
102,247✔
383
  simulation::time_tallies.stop();
102,247✔
384

385
  // update weight windows if needed
386
  for (const auto& wwg : variance_reduction::weight_windows_generators) {
102,387✔
387
    wwg->update();
140✔
388
  }
389

390
  // Reset global tally results
391
  if (simulation::current_batch <= settings::n_inactive) {
102,247✔
392
    xt::view(simulation::global_tallies, xt::all()) = 0.0;
16,639✔
393
    simulation::n_realizations = 0;
16,639✔
394
  }
395

396
  // Check_triggers
397
  if (mpi::master)
102,247✔
398
    check_triggers();
86,112✔
399
#ifdef OPENMC_MPI
400
  MPI_Bcast(&simulation::satisfy_triggers, 1, MPI_C_BOOL, 0, mpi::intracomm);
50,538✔
401
#endif
402
  if (simulation::satisfy_triggers ||
102,247✔
403
      (settings::trigger_on &&
3,258✔
404
        simulation::current_batch == settings::n_max_batches)) {
3,258✔
405
    settings::statepoint_batch.insert(simulation::current_batch);
165✔
406
  }
407

408
  // Write out state point if it's been specified for this batch and is not
409
  // a CMFD run instance
410
  if (contains(settings::statepoint_batch, simulation::current_batch) &&
108,707✔
411
      !settings::cmfd_run) {
6,460✔
412
    if (contains(settings::sourcepoint_batch, simulation::current_batch) &&
11,737✔
413
        settings::source_write && !settings::source_separate) {
11,737✔
414
      bool b = (settings::run_mode == RunMode::EIGENVALUE);
5,627✔
415
      openmc_statepoint_write(nullptr, &b);
5,627✔
416
    } else {
417
      bool b = false;
369✔
418
      openmc_statepoint_write(nullptr, &b);
369✔
419
    }
420
  }
421

422
  if (settings::run_mode == RunMode::EIGENVALUE) {
102,247✔
423
    // Write out a separate source point if it's been specified for this batch
424
    if (contains(settings::sourcepoint_batch, simulation::current_batch) &&
90,286✔
425
        settings::source_write && settings::source_separate) {
90,286✔
426

427
      // Determine width for zero padding
428
      int w = std::to_string(settings::n_max_batches).size();
76✔
429
      std::string source_point_filename = fmt::format("{0}source.{1:0{2}}",
430
        settings::path_output, simulation::current_batch, w);
76✔
431
      gsl::span<SourceSite> bankspan(simulation::source_bank);
76✔
432
      if (settings::source_mcpl_write) {
76✔
433
        write_mcpl_source_point(
19✔
434
          source_point_filename.c_str(), bankspan, simulation::work_index);
435
      } else {
436
        write_source_point(
57✔
437
          source_point_filename.c_str(), bankspan, simulation::work_index);
438
      }
439
    }
76✔
440

441
    // Write a continously-overwritten source point if requested.
442
    if (settings::source_latest) {
86,050✔
443

444
      // note: correct file extension appended automatically
445
      auto filename = settings::path_output + "source";
190✔
446
      gsl::span<SourceSite> bankspan(simulation::source_bank);
190✔
447
      if (settings::source_mcpl_write) {
190✔
448
        write_mcpl_source_point(
×
449
          filename.c_str(), bankspan, simulation::work_index);
450
      } else {
451
        write_source_point(filename.c_str(), bankspan, simulation::work_index);
190✔
452
      }
453
    }
190✔
454
  }
455

456
  // Write out surface source if requested.
457
  if (settings::surf_source_write &&
102,247✔
458
      simulation::current_batch == settings::n_batches) {
150✔
459
    auto filename = settings::path_output + "surface_source";
16✔
460
    auto surf_work_index =
461
      mpi::calculate_parallel_index_vector(simulation::surf_source_bank.size());
16✔
462
    gsl::span<SourceSite> surfbankspan(simulation::surf_source_bank.begin(),
463
      simulation::surf_source_bank.size());
16✔
464
    if (settings::surf_mcpl_write) {
16✔
465
      write_mcpl_source_point(filename.c_str(), surfbankspan, surf_work_index);
×
466
    } else {
467
      write_source_point(filename.c_str(), surfbankspan, surf_work_index);
16✔
468
    }
469
  }
16✔
470
}
102,247✔
471

472
void initialize_generation()
102,527✔
473
{
474
  if (settings::run_mode == RunMode::EIGENVALUE) {
102,527✔
475
    // Clear out the fission bank
476
    simulation::fission_bank.resize(0);
86,316✔
477

478
    // Count source sites if using uniform fission source weighting
479
    if (settings::ufs_on)
86,316✔
480
      ufs_count_sites();
190✔
481

482
    // Store current value of tracklength k
483
    simulation::keff_generation = simulation::global_tallies(
86,316✔
484
      GlobalTally::K_TRACKLENGTH, TallyResult::VALUE);
485
  }
486
}
102,527✔
487

488
void finalize_generation()
102,513✔
489
{
490
  auto& gt = simulation::global_tallies;
102,513✔
491

492
  // Update global tallies with the accumulation variables
493
  if (settings::run_mode == RunMode::EIGENVALUE) {
102,513✔
494
    gt(GlobalTally::K_COLLISION, TallyResult::VALUE) += global_tally_collision;
86,316✔
495
    gt(GlobalTally::K_ABSORPTION, TallyResult::VALUE) +=
86,316✔
496
      global_tally_absorption;
497
    gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) +=
86,316✔
498
      global_tally_tracklength;
499
  }
500
  gt(GlobalTally::LEAKAGE, TallyResult::VALUE) += global_tally_leakage;
102,513✔
501

502
  // reset tallies
503
  if (settings::run_mode == RunMode::EIGENVALUE) {
102,513✔
504
    global_tally_collision = 0.0;
86,316✔
505
    global_tally_absorption = 0.0;
86,316✔
506
    global_tally_tracklength = 0.0;
86,316✔
507
  }
508
  global_tally_leakage = 0.0;
102,513✔
509

510
  if (settings::run_mode == RunMode::EIGENVALUE &&
102,513✔
511
      settings::solver_type == SolverType::MONTE_CARLO) {
86,316✔
512
    // If using shared memory, stable sort the fission bank (by parent IDs)
513
    // so as to allow for reproducibility regardless of which order particles
514
    // are run in.
515
    sort_fission_bank();
85,936✔
516

517
    // Distribute fission bank across processors evenly
518
    synchronize_bank();
85,936✔
519
  }
520

521
  if (settings::run_mode == RunMode::EIGENVALUE) {
102,513✔
522

523
    // Calculate shannon entropy
524
    if (settings::entropy_on)
86,316✔
525
      shannon_entropy();
14,055✔
526

527
    // Collect results and statistics
528
    calculate_generation_keff();
86,316✔
529
    calculate_average_keff();
86,316✔
530

531
    // Write generation output
532
    if (mpi::master && settings::verbosity >= 7) {
86,316✔
533
      print_generation();
66,964✔
534
    }
535
  }
536
}
102,513✔
537

538
void initialize_history(Particle& p, int64_t index_source)
186,549,130✔
539
{
540
  // set defaults
541
  if (settings::run_mode == RunMode::EIGENVALUE) {
186,549,130✔
542
    // set defaults for eigenvalue simulations from primary bank
543
    p.from_source(&simulation::source_bank[index_source - 1]);
172,328,700✔
544
  } else if (settings::run_mode == RunMode::FIXED_SOURCE) {
14,220,430✔
545
    // initialize random number seed
546
    int64_t id = (simulation::total_gen + overall_generation() - 1) *
14,220,430✔
547
                   settings::n_particles +
14,220,430✔
548
                 simulation::work_index[mpi::rank] + index_source;
14,220,430✔
549
    uint64_t seed = init_seed(id, STREAM_SOURCE);
14,220,430✔
550
    // sample from external source distribution or custom library then set
551
    auto site = sample_external_source(&seed);
14,220,430✔
552
    p.from_source(&site);
14,220,424✔
553
  }
554
  p.current_work() = index_source;
186,549,124✔
555

556
  // set identifier for particle
557
  p.id() = simulation::work_index[mpi::rank] + index_source;
186,549,124✔
558

559
  // set progeny count to zero
560
  p.n_progeny() = 0;
186,549,124✔
561

562
  // Reset particle event counter
563
  p.n_event() = 0;
186,549,124✔
564

565
  // Reset split counter
566
  p.n_split() = 0;
186,549,124✔
567

568
  // Reset weight window ratio
569
  p.ww_factor() = 0.0;
186,549,124✔
570

571
  // Reset pulse_height_storage
572
  std::fill(p.pht_storage().begin(), p.pht_storage().end(), 0);
186,549,124✔
573

574
  // set random number seed
575
  int64_t particle_seed =
576
    (simulation::total_gen + overall_generation() - 1) * settings::n_particles +
186,549,124✔
577
    p.id();
186,549,124✔
578
  init_particle_seeds(particle_seed, p.seeds());
186,549,124✔
579

580
  // set particle trace
581
  p.trace() = false;
186,549,124✔
582
  if (simulation::current_batch == settings::trace_batch &&
373,112,248✔
583
      simulation::current_gen == settings::trace_gen &&
186,563,124✔
584
      p.id() == settings::trace_particle)
14,000✔
585
    p.trace() = true;
14✔
586

587
  // Set particle track.
588
  p.write_track() = check_track_criteria(p);
186,549,124✔
589

590
  // Display message if high verbosity or trace is on
591
  if (settings::verbosity >= 9 || p.trace()) {
186,549,124✔
592
    write_message("Simulating Particle {}", p.id());
14✔
593
  }
594

595
// Add paricle's starting weight to count for normalizing tallies later
596
#pragma omp atomic
80,689,119✔
597
  simulation::total_weight += p.wgt();
186,549,124✔
598

599
  // Force calculation of cross-sections by setting last energy to zero
600
  if (settings::run_CE) {
186,549,124✔
601
    p.invalidate_neutron_xs();
45,373,124✔
602
  }
603

604
  // Prepare to write out particle track.
605
  if (p.write_track())
186,549,124✔
606
    add_particle_track(p);
1,266✔
607
}
186,549,124✔
608

609
int overall_generation()
201,009,057✔
610
{
611
  using namespace simulation;
612
  return settings::gen_per_batch * (current_batch - 1) + current_gen;
201,009,057✔
613
}
614

615
void calculate_work()
6,150✔
616
{
617
  // Determine minimum amount of particles to simulate on each processor
618
  int64_t min_work = settings::n_particles / mpi::n_procs;
6,150✔
619

620
  // Determine number of processors that have one extra particle
621
  int64_t remainder = settings::n_particles % mpi::n_procs;
6,150✔
622

623
  int64_t i_bank = 0;
6,150✔
624
  simulation::work_index.resize(mpi::n_procs + 1);
6,150✔
625
  simulation::work_index[0] = 0;
6,150✔
626
  for (int i = 0; i < mpi::n_procs; ++i) {
13,825✔
627
    // Number of particles for rank i
628
    int64_t work_i = i < remainder ? min_work + 1 : min_work;
7,675✔
629

630
    // Set number of particles
631
    if (mpi::rank == i)
7,675✔
632
      simulation::work_per_rank = work_i;
6,150✔
633

634
    // Set index into source bank for rank i
635
    i_bank += work_i;
7,675✔
636
    simulation::work_index[i + 1] = i_bank;
7,675✔
637
  }
638
}
6,150✔
639

640
void initialize_data()
5,565✔
641
{
642
  // Determine minimum/maximum energy for incident neutron/photon data
643
  data::energy_max = {INFTY, INFTY};
5,565✔
644
  data::energy_min = {0.0, 0.0};
5,565✔
645
  for (const auto& nuc : data::nuclides) {
38,143✔
646
    if (nuc->grid_.size() >= 1) {
32,578✔
647
      int neutron = static_cast<int>(ParticleType::neutron);
32,578✔
648
      data::energy_min[neutron] =
65,156✔
649
        std::max(data::energy_min[neutron], nuc->grid_[0].energy.front());
32,578✔
650
      data::energy_max[neutron] =
32,578✔
651
        std::min(data::energy_max[neutron], nuc->grid_[0].energy.back());
32,578✔
652
    }
653
  }
654

655
  if (settings::photon_transport) {
5,565✔
656
    for (const auto& elem : data::elements) {
887✔
657
      if (elem->energy_.size() >= 1) {
613✔
658
        int photon = static_cast<int>(ParticleType::photon);
613✔
659
        int n = elem->energy_.size();
613✔
660
        data::energy_min[photon] =
1,226✔
661
          std::max(data::energy_min[photon], std::exp(elem->energy_(1)));
613✔
662
        data::energy_max[photon] =
1,226✔
663
          std::min(data::energy_max[photon], std::exp(elem->energy_(n - 1)));
613✔
664
      }
665
    }
666

667
    if (settings::electron_treatment == ElectronTreatment::TTB) {
274✔
668
      // Determine if minimum/maximum energy for bremsstrahlung is greater/less
669
      // than the current minimum/maximum
670
      if (data::ttb_e_grid.size() >= 1) {
274✔
671
        int photon = static_cast<int>(ParticleType::photon);
274✔
672
        int n_e = data::ttb_e_grid.size();
274✔
673
        data::energy_min[photon] =
548✔
674
          std::max(data::energy_min[photon], std::exp(data::ttb_e_grid(1)));
274✔
675
        data::energy_max[photon] = std::min(
274✔
676
          data::energy_max[photon], std::exp(data::ttb_e_grid(n_e - 1)));
548✔
677
      }
678
    }
679
  }
680

681
  // Show which nuclide results in lowest energy for neutron transport
682
  for (const auto& nuc : data::nuclides) {
6,683✔
683
    // If a nuclide is present in a material that's not used in the model, its
684
    // grid has not been allocated
685
    if (nuc->grid_.size() > 0) {
6,180✔
686
      double max_E = nuc->grid_[0].energy.back();
6,180✔
687
      int neutron = static_cast<int>(ParticleType::neutron);
6,180✔
688
      if (max_E == data::energy_max[neutron]) {
6,180✔
689
        write_message(7, "Maximum neutron transport energy: {} eV for {}",
5,062✔
690
          data::energy_max[neutron], nuc->name_);
5,062✔
691
        if (mpi::master && data::energy_max[neutron] < 20.0e6) {
5,062✔
692
          warning("Maximum neutron energy is below 20 MeV. This may bias "
×
693
                  "the results.");
694
        }
695
        break;
5,062✔
696
      }
697
    }
698
  }
699

700
  // Set up logarithmic grid for nuclides
701
  for (auto& nuc : data::nuclides) {
38,143✔
702
    nuc->init_grid();
32,578✔
703
  }
704
  int neutron = static_cast<int>(ParticleType::neutron);
5,565✔
705
  simulation::log_spacing =
5,565✔
706
    std::log(data::energy_max[neutron] / data::energy_min[neutron]) /
5,565✔
707
    settings::n_log_bins;
708
}
5,565✔
709

710
#ifdef OPENMC_MPI
711
void broadcast_results()
2,920✔
712
{
713
  // Broadcast tally results so that each process has access to results
714
  for (auto& t : model::tallies) {
18,169✔
715
    // Create a new datatype that consists of all values for a given filter
716
    // bin and then use that to broadcast. This is done to minimize the
717
    // chance of the 'count' argument of MPI_BCAST exceeding 2**31
718
    auto& results = t->results_;
15,249✔
719

720
    auto shape = results.shape();
15,249✔
721
    int count_per_filter = shape[1] * shape[2];
15,249✔
722
    MPI_Datatype result_block;
723
    MPI_Type_contiguous(count_per_filter, MPI_DOUBLE, &result_block);
15,249✔
724
    MPI_Type_commit(&result_block);
15,249✔
725
    MPI_Bcast(results.data(), shape[0], result_block, 0, mpi::intracomm);
15,249✔
726
    MPI_Type_free(&result_block);
15,249✔
727
  }
728

729
  // Also broadcast global tally results
730
  auto& gt = simulation::global_tallies;
2,920✔
731
  MPI_Bcast(gt.data(), gt.size(), MPI_DOUBLE, 0, mpi::intracomm);
2,920✔
732

733
  // These guys are needed so that non-master processes can calculate the
734
  // combined estimate of k-effective
735
  double temp[] {
736
    simulation::k_col_abs, simulation::k_col_tra, simulation::k_abs_tra};
2,920✔
737
  MPI_Bcast(temp, 3, MPI_DOUBLE, 0, mpi::intracomm);
2,920✔
738
  simulation::k_col_abs = temp[0];
2,920✔
739
  simulation::k_col_tra = temp[1];
2,920✔
740
  simulation::k_abs_tra = temp[2];
2,920✔
741
}
2,920✔
742

743
#endif
744

745
void free_memory_simulation()
6,784✔
746
{
747
  simulation::k_generation.clear();
6,784✔
748
  simulation::entropy.clear();
6,784✔
749
}
6,784✔
750

751
void transport_history_based_single_particle(Particle& p)
175,572,266✔
752
{
753
  while (p.alive()) {
2,147,483,647✔
754
    p.event_calculate_xs();
2,147,483,647✔
755
    if (!p.alive())
2,147,483,647✔
UNCOV
756
      break;
×
757
    p.event_advance();
2,147,483,647✔
758
    if (p.collision_distance() > p.boundary().distance) {
2,147,483,647✔
759
      p.event_cross_surface();
1,596,175,622✔
760
    } else {
761
      p.event_collide();
2,147,483,647✔
762
    }
763
    p.event_revive_from_secondary();
2,147,483,647✔
764
  }
765
  p.event_death();
175,572,258✔
766
}
175,572,258✔
767

768
void transport_history_based()
99,161✔
769
{
770
#pragma omp parallel for schedule(runtime)
771
  for (int64_t i_work = 1; i_work <= simulation::work_per_rank; ++i_work) {
105,955,223✔
772
    Particle p;
105,901,150✔
773
    initialize_history(p, i_work);
105,901,150✔
774
    transport_history_based_single_particle(p);
105,901,144✔
775
  }
105,901,136✔
776
}
99,147✔
777

778
void transport_event_based()
2,986✔
779
{
780
  int64_t remaining_work = simulation::work_per_rank;
2,986✔
781
  int64_t source_offset = 0;
2,986✔
782

783
  // To cap the total amount of memory used to store particle object data, the
784
  // number of particles in flight at any point in time can bet set. In the case
785
  // that the maximum in flight particle count is lower than the total number
786
  // of particles that need to be run this iteration, the event-based transport
787
  // loop is executed multiple times until all particles have been completed.
788
  while (remaining_work > 0) {
5,972✔
789
    // Figure out # of particles to run for this subiteration
790
    int64_t n_particles =
791
      std::min(remaining_work, settings::max_particles_in_flight);
2,986✔
792

793
    // Initialize all particle histories for this subiteration
794
    process_init_events(n_particles, source_offset);
2,986✔
795

796
    // Event-based transport loop
797
    while (true) {
798
      // Determine which event kernel has the longest queue
799
      int64_t max = std::max({simulation::calculate_fuel_xs_queue.size(),
4,698,320✔
800
        simulation::calculate_nonfuel_xs_queue.size(),
2,349,160✔
801
        simulation::advance_particle_queue.size(),
2,349,160✔
802
        simulation::surface_crossing_queue.size(),
2,349,160✔
803
        simulation::collision_queue.size()});
2,349,160✔
804

805
      // Execute event with the longest queue
806
      if (max == 0) {
2,349,160✔
807
        break;
2,986✔
808
      } else if (max == simulation::calculate_fuel_xs_queue.size()) {
2,346,174✔
809
        process_calculate_xs_events(simulation::calculate_fuel_xs_queue);
421,740✔
810
      } else if (max == simulation::calculate_nonfuel_xs_queue.size()) {
1,924,434✔
811
        process_calculate_xs_events(simulation::calculate_nonfuel_xs_queue);
357,537✔
812
      } else if (max == simulation::advance_particle_queue.size()) {
1,566,897✔
813
        process_advance_particle_events();
774,513✔
814
      } else if (max == simulation::surface_crossing_queue.size()) {
792,384✔
815
        process_surface_crossing_events();
254,319✔
816
      } else if (max == simulation::collision_queue.size()) {
538,065✔
817
        process_collision_events();
538,065✔
818
      }
819
    }
2,346,174✔
820

821
    // Execute death event for all particles
822
    process_death_events(n_particles);
2,986✔
823

824
    // Adjust remaining work and source offset variables
825
    remaining_work -= n_particles;
2,986✔
826
    source_offset += n_particles;
2,986✔
827
  }
828
}
2,986✔
829

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