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

openmc-dev / openmc / 20932149956

12 Jan 2026 07:24PM UTC coverage: 81.996% (-0.2%) from 82.195%
20932149956

Pull #3702

github

web-flow
Merge 80e4cff3d into 0486e433d
Pull Request #3702: Random Ray Kinetic Simulation Mode

17686 of 24568 branches covered (71.99%)

Branch coverage included in aggregate %.

1002 of 1145 new or added lines in 20 files covered. (87.51%)

484 existing lines in 23 files now uncovered.

56448 of 65844 relevant lines covered (85.73%)

52849177.45 hits per line

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

85.15
/src/random_ray/random_ray_simulation.cpp
1
#include "openmc/random_ray/random_ray_simulation.h"
2

3
#include "openmc/eigenvalue.h"
4
#include "openmc/geometry.h"
5
#include "openmc/message_passing.h"
6
#include "openmc/mgxs_interface.h"
7
#include "openmc/output.h"
8
#include "openmc/plot.h"
9
#include "openmc/random_ray/flat_source_domain.h"
10
#include "openmc/random_ray/random_ray.h"
11
#include "openmc/simulation.h"
12
#include "openmc/source.h"
13
#include "openmc/tallies/filter.h"
14
#include "openmc/tallies/tally.h"
15
#include "openmc/tallies/tally_scoring.h"
16
#include "openmc/timer.h"
17
#include "openmc/weight_windows.h"
18

19
namespace openmc {
20

21
//==============================================================================
22
// Non-member functions
23
//==============================================================================
24

25
void openmc_run_random_ray()
833✔
26
{
27
  //////////////////////////////////////////////////////////
28
  // Run forward simulation
29
  //////////////////////////////////////////////////////////
30

31
  if (mpi::master) {
833✔
32
    bool adjoint_needed = FlatSourceDomain::adjoint_;
573✔
33
    print_random_ray_headers(adjoint_needed);
573✔
34
  }
35

36
  // Initialize OpenMC general data structures
37
  openmc_simulation_init();
833✔
38

39
  // Validate that inputs meet requirements for random ray mode
40
  if (mpi::master)
833✔
41
    validate_random_ray_inputs();
573✔
42

43
  // Initialize Random Ray Simulation Object
44
  RandomRaySimulation sim;
833✔
45

46
  // Initialize fixed sources, if present
47
  sim.apply_fixed_sources_and_mesh_domains();
833✔
48

49
  // Run initial random ray simulation
50
  sim.run_single_simulation();
833✔
51

52
  if (settings::kinetic_simulation) {
833✔
53
    // Second steady state simulation to correct the batchwise k-eff
54
    sim.kinetic_initial_condition();
192✔
55

56
    // Timestepping loop
57
    for (int i = 0; i < settings::n_timesteps; i++)
1,152✔
58
      sim.kinetic_single_time_step(i);
960✔
59
  }
60

61
  //////////////////////////////////////////////////////////
62
  // Run adjoint simulation (if enabled)
63
  //////////////////////////////////////////////////////////
64

65
  sim.random_ray_adjoint();
833✔
66
}
833✔
67

68
// Enforces restrictions on inputs in random ray mode.  While there are
69
// many features that don't make sense in random ray mode, and are therefore
70
// unsupported, we limit our testing/enforcement operations only to inputs
71
// that may cause erroneous/misleading output or crashes from the solver.
72
void validate_random_ray_inputs()
573✔
73
{
74
  // Validate tallies
75
  ///////////////////////////////////////////////////////////////////
76
  for (auto& tally : model::tallies) {
1,542✔
77

78
    // Validate score types
79
    for (auto score_bin : tally->scores_) {
2,202✔
80
      switch (score_bin) {
1,233!
81
      case SCORE_FLUX:
1,189✔
82
      case SCORE_TOTAL:
83
      case SCORE_FISSION:
84
      case SCORE_NU_FISSION:
85
      case SCORE_EVENTS:
86
        break;
1,189✔
87

88
      // TODO: add support for prompt and delayed fission
89
      case SCORE_PRECURSORS: {
44✔
90
        if (settings::kinetic_simulation && settings::create_delayed_neutrons) {
44!
91
          break;
92
        } else {
NEW
93
          fatal_error("Invalid score specified in tallies.xml. Precursors can "
×
94
                      "only be scored for kinetic simulations.");
95
        }
96
      }
97
      default:
98
        fatal_error(
×
99
          "Invalid score specified. Only flux, total, fission, nu-fission, and "
100
          "event scores are supported in random ray mode. (precursors are "
101
          "supported for "
102
          "kinetic simulations when delayed neutrons are turned on).");
103
      }
104
    }
105

106
    // Validate filter types
107
    for (auto f : tally->filters()) {
2,116✔
108
      auto& filter = *model::tally_filters[f];
1,147✔
109

110
      switch (filter.type()) {
1,147!
111
      case FilterType::CELL:
1,103✔
112
      case FilterType::CELL_INSTANCE:
113
      case FilterType::DISTRIBCELL:
114
      case FilterType::ENERGY:
115
      case FilterType::MATERIAL:
116
      case FilterType::MESH:
117
      case FilterType::UNIVERSE:
118
      case FilterType::PARTICLE:
119
        break;
1,103✔
120
      case FilterType::DELAYED_GROUP:
44✔
121
        if (settings::kinetic_simulation) {
44!
122
          break;
44✔
123
        } else {
NEW
124
          fatal_error("Invalid filter specified in tallies.xml. Kinetic "
×
125
                      "simulations is required "
126
                      "to tally with a delayed_group filter.");
127
        }
128
      default:
×
129
        fatal_error("Invalid filter specified. Only cell, cell_instance, "
×
130
                    "distribcell, energy, material, mesh, and universe filters "
131
                    "are supported in random ray mode (delayed_group is "
132
                    "supported for kinetic simulations).");
133
      }
134
    }
135
  }
136

137
  // TODO: validate kinetic data is present
138
  //  Validate MGXS data
139
  ///////////////////////////////////////////////////////////////////
140
  for (auto& material : data::mg.macro_xs_) {
2,158✔
141
    if (!material.is_isotropic) {
1,585!
142
      fatal_error("Anisotropic MGXS detected. Only isotropic XS data sets "
×
143
                  "supported in random ray mode.");
144
    }
145
    if (material.get_xsdata().size() > 1) {
1,585!
146
      warning("Non-isothermal MGXS detected. Only isothermal XS data sets "
×
147
              "supported in random ray mode. Using lowest temperature.");
148
    }
149
    for (int g = 0; g < data::mg.num_energy_groups_; g++) {
12,939✔
150
      if (material.exists_in_model) {
11,354✔
151
        // Temperature and angle indices, if using multiple temperature
152
        // data sets and/or anisotropic data sets.
153
        // TODO: Currently assumes we are only using single temp/single angle
154
        // data.
155
        const int t = 0;
11,310✔
156
        const int a = 0;
11,310✔
157
        double sigma_t =
158
          material.get_xs(MgxsType::TOTAL, g, NULL, NULL, NULL, t, a);
11,310✔
159
        if (sigma_t <= 0.0) {
11,310!
160
          fatal_error("No zero or negative total macroscopic cross sections "
×
161
                      "allowed in random ray mode. If the intention is to make "
162
                      "a void material, use a cell fill of 'None' instead.");
163
        }
164
      }
165
    }
166
  }
167

168
  // Validate ray source
169
  ///////////////////////////////////////////////////////////////////
170

171
  // Check for independent source
172
  IndependentSource* is =
173
    dynamic_cast<IndependentSource*>(RandomRay::ray_source_.get());
573!
174
  if (!is) {
573!
175
    fatal_error("Invalid ray source definition. Ray source must provided and "
×
176
                "be of type IndependentSource.");
177
  }
178

179
  // Check for box source
180
  SpatialDistribution* space_dist = is->space();
573✔
181
  SpatialBox* sb = dynamic_cast<SpatialBox*>(space_dist);
573!
182
  if (!sb) {
573!
183
    fatal_error(
×
184
      "Invalid ray source definition -- only box sources are allowed.");
185
  }
186

187
  // Check that box source is not restricted to fissionable areas
188
  if (sb->only_fissionable()) {
573!
189
    fatal_error(
×
190
      "Invalid ray source definition -- fissionable spatial distribution "
191
      "not allowed.");
192
  }
193

194
  // Check for isotropic source
195
  UnitSphereDistribution* angle_dist = is->angle();
573✔
196
  Isotropic* id = dynamic_cast<Isotropic*>(angle_dist);
573!
197
  if (!id) {
573!
198
    fatal_error("Invalid ray source definition -- only isotropic sources are "
×
199
                "allowed.");
200
  }
201

202
  // Validate external sources
203
  ///////////////////////////////////////////////////////////////////
204
  if (settings::run_mode == RunMode::FIXED_SOURCE) {
573✔
205
    if (model::external_sources.size() < 1) {
276!
206
      fatal_error("Must provide a particle source (in addition to ray source) "
×
207
                  "in fixed source random ray mode.");
208
    }
209

210
    for (int i = 0; i < model::external_sources.size(); i++) {
552✔
211
      Source* s = model::external_sources[i].get();
276✔
212

213
      // Check for independent source
214
      IndependentSource* is = dynamic_cast<IndependentSource*>(s);
276!
215

216
      if (!is) {
276!
217
        fatal_error(
×
218
          "Only IndependentSource external source types are allowed in "
219
          "random ray mode");
220
      }
221

222
      // Check for isotropic source
223
      UnitSphereDistribution* angle_dist = is->angle();
276✔
224
      Isotropic* id = dynamic_cast<Isotropic*>(angle_dist);
276!
225
      if (!id) {
276!
226
        fatal_error(
×
227
          "Invalid source definition -- only isotropic external sources are "
228
          "allowed in random ray mode.");
229
      }
230

231
      // Validate that a domain ID was specified OR that it is a point source
232
      auto sp = dynamic_cast<SpatialPoint*>(is->space());
276!
233
      if (is->domain_ids().size() == 0 && !sp) {
276!
234
        fatal_error("Fixed sources must be point source or spatially "
×
235
                    "constrained by domain id (cell, material, or universe) in "
236
                    "random ray mode.");
237
      } else if (is->domain_ids().size() > 0 && sp) {
276✔
238
        // If both a domain constraint and a non-default point source location
239
        // are specified, notify user that domain constraint takes precedence.
240
        if (sp->r().x == 0.0 && sp->r().y == 0.0 && sp->r().z == 0.0) {
242!
241
          warning("Fixed source has both a domain constraint and a point "
242✔
242
                  "type spatial distribution. The domain constraint takes "
243
                  "precedence in random ray mode -- point source coordinate "
244
                  "will be ignored.");
245
        }
246
      }
247

248
      // Check that a discrete energy distribution was used
249
      Distribution* d = is->energy();
276✔
250
      Discrete* dd = dynamic_cast<Discrete*>(d);
276!
251
      if (!dd) {
276!
252
        fatal_error(
×
253
          "Only discrete (multigroup) energy distributions are allowed for "
254
          "external sources in random ray mode.");
255
      }
256
    }
257
  }
258

259
  // Validate plotting files
260
  ///////////////////////////////////////////////////////////////////
261
  for (int p = 0; p < model::plots.size(); p++) {
573!
262

263
    // Get handle to OpenMC plot object
264
    const auto& openmc_plottable = model::plots[p];
×
265
    Plot* openmc_plot = dynamic_cast<Plot*>(openmc_plottable.get());
×
266

267
    // Random ray plots only support voxel plots
268
    if (!openmc_plot) {
×
269
      warning(fmt::format(
×
270
        "Plot {} will not be used for end of simulation data plotting -- only "
271
        "voxel plotting is allowed in random ray mode.",
272
        openmc_plottable->id()));
×
273
      continue;
×
274
    } else if (openmc_plot->type_ != Plot::PlotType::voxel) {
×
275
      warning(fmt::format(
×
276
        "Plot {} will not be used for end of simulation data plotting -- only "
277
        "voxel plotting is allowed in random ray mode.",
278
        openmc_plottable->id()));
×
279
      continue;
×
280
    }
281
  }
282

283
  // Warn about slow MPI domain replication, if detected
284
  ///////////////////////////////////////////////////////////////////
285
#ifdef OPENMC_MPI
286
  if (mpi::n_procs > 1) {
261✔
287
    warning(
260✔
288
      "MPI parallelism is not supported by the random ray solver. All work "
289
      "will be performed by rank 0. Domain decomposition may be implemented in "
290
      "the future to provide efficient MPI scaling.");
291
  }
292
#endif
293

294
  // Warn about instability resulting from linear sources in small regions
295
  // when generating weight windows with FW-CADIS and an overlaid mesh.
296
  ///////////////////////////////////////////////////////////////////
297
  if (RandomRay::source_shape_ == RandomRaySourceShape::LINEAR &&
749✔
298
      variance_reduction::weight_windows.size() > 0) {
176✔
299
    warning(
11✔
300
      "Linear sources may result in negative fluxes in small source regions "
301
      "generated by mesh subdivision. Negative sources may result in low "
302
      "quality FW-CADIS weight windows. We recommend you use flat source mode "
303
      "when generating weight windows with an overlaid mesh tally.");
304
  }
305
}
573✔
306

307
void openmc_reset_random_ray()
8,342✔
308
{
309
  FlatSourceDomain::volume_estimator_ = RandomRayVolumeEstimator::HYBRID;
8,342✔
310
  FlatSourceDomain::volume_normalized_flux_tallies_ = false;
8,342✔
311
  FlatSourceDomain::adjoint_ = false;
8,342✔
312
  FlatSourceDomain::mesh_domain_map_.clear();
8,342✔
313
  RandomRay::ray_source_.reset();
8,342✔
314
  RandomRay::source_shape_ = RandomRaySourceShape::FLAT;
8,342✔
315
  RandomRay::sample_method_ = RandomRaySampleMethod::PRNG;
8,342✔
316
  RandomRay::bd_order_ = 3;
8,342✔
317
  RandomRay::time_method_ = RandomRayTimeMethod::ISOTROPIC;
8,342✔
318
}
8,342✔
319

320
void write_random_ray_hdf5(hid_t group)
1,421✔
321
{
322
  hid_t random_ray_group = create_group(group, "random_ray");
1,421✔
323
  switch (RandomRay::source_shape_) {
1,421!
324
  case RandomRaySourceShape::FLAT:
1,201✔
325
    write_dataset(random_ray_group, "source_shape", "flat");
1,201✔
326
    break;
1,201✔
327
  case RandomRaySourceShape::LINEAR:
187✔
328
    write_dataset(random_ray_group, "source_shape", "linear");
187✔
329
    break;
187✔
330
  case RandomRaySourceShape::LINEAR_XY:
33✔
331
    write_dataset(random_ray_group, "source_shape", "linear xy");
33✔
332
    break;
33✔
NEW
333
  default:
×
NEW
334
    break;
×
335
  }
336

337
  switch (FlatSourceDomain::volume_estimator_) {
1,421!
338
  case RandomRayVolumeEstimator::SIMULATION_AVERAGED:
22✔
339
    write_dataset(random_ray_group, "volume_estimator", "simulation averaged");
22✔
340
    break;
22✔
341
  case RandomRayVolumeEstimator::NAIVE:
68✔
342
    write_dataset(random_ray_group, "volume_estimator", "naive");
68✔
343
    break;
68✔
344
  case RandomRayVolumeEstimator::HYBRID:
1,331✔
345
    write_dataset(random_ray_group, "volume_estimator", "hybrid");
1,331✔
346
    break;
1,331✔
NEW
347
  default:
×
NEW
348
    break;
×
349
  }
350

351
  write_dataset(
1,421✔
352
    random_ray_group, "distance_active", RandomRay::distance_active_);
353
  write_dataset(
1,421✔
354
    random_ray_group, "distance_inactive", RandomRay::distance_inactive_);
355
  write_dataset(random_ray_group, "volume_normalized_flux_tallies",
1,421✔
356
    FlatSourceDomain::volume_normalized_flux_tallies_);
357
  write_dataset(random_ray_group, "adjoint_mode", FlatSourceDomain::adjoint_);
1,421✔
358

359
  write_dataset(random_ray_group, "avg_miss_rate", RandomRay::avg_miss_rate_);
1,421✔
360
  write_dataset(
1,421✔
361
    random_ray_group, "n_source_regions", RandomRay::n_source_regions_);
362
  write_dataset(random_ray_group, "n_external_source_regions",
1,421✔
363
    RandomRay::n_external_source_regions_);
364
  write_dataset(random_ray_group, "n_geometric_intersections",
1,421✔
365
    RandomRay::total_geometric_intersections_);
366
  int64_t n_integrations =
1,421✔
367
    RandomRay::total_geometric_intersections_ * data::mg.num_energy_groups_;
1,421✔
368
  write_dataset(random_ray_group, "n_integrations", n_integrations);
1,421✔
369

370
  if (settings::kinetic_simulation && !simulation::is_initial_condition) {
1,421✔
371
    write_dataset(random_ray_group, "bd_order", RandomRay::bd_order_);
660✔
372
    switch (RandomRay::time_method_) {
660!
373
    case RandomRayTimeMethod::ISOTROPIC:
440✔
374
      write_dataset(random_ray_group, "time_method", "isotropic");
440✔
375
      break;
440✔
376
    case RandomRayTimeMethod::PROPAGATION:
220✔
377
      write_dataset(random_ray_group, "time_method", "propogation");
220✔
378
      break;
220✔
NEW
379
    default:
×
NEW
380
      break;
×
381
    }
382
  }
383
  close_group(random_ray_group);
1,421✔
384
}
1,421✔
385

386
void print_random_ray_headers(bool& adjoint_needed)
1,421✔
387
{
388
  // If we're going to do an adjoint simulation afterwards, report that this is
389
  // the initial forward flux solve.
390
  if (adjoint_needed && !FlatSourceDomain::adjoint_)
1,421!
NEW
391
    header("FORWARD FLUX SOLVE", 3);
×
392

393
  // Otherwise report that we are doing the adjoint simulation
394
  if (adjoint_needed && FlatSourceDomain::adjoint_)
1,421!
395
    header("ADJOINT FLUX SOLVE", 3);
112✔
396
}
1,421✔
397

398
//-----------------------------------------------------------------------------
399
// Non-member functions for kinetic simulations
400

401
void rename_time_step_file(
2,304✔
402
  std::string base_filename, std::string extension, int i)
403
{
404
  // Rename file
405
  std::string old_filename_ =
406
    fmt::format("{0}{1}{2}", settings::path_output, base_filename, extension);
2,304✔
407
  std::string new_filename_ = fmt::format(
408
    "{0}{1}_{2}{3}", settings::path_output, base_filename, i, extension);
1,872✔
409

410
  const char* old_fname = old_filename_.c_str();
2,304✔
411
  const char* new_fname = new_filename_.c_str();
2,304✔
412
  std::rename(old_fname, new_fname);
2,304✔
413
}
2,304✔
414

415
//==============================================================================
416
// RandomRaySimulation implementation
417
//==============================================================================
418

419
RandomRaySimulation::RandomRaySimulation()
833✔
420
  : negroups_(data::mg.num_energy_groups_),
833✔
421
    ndgroups_(data::mg.num_delayed_groups_)
833✔
422
{
423
  // There are no source sites in random ray mode, so be sure to disable to
424
  // ensure we don't attempt to write source sites to statepoint
425
  settings::source_write = false;
833✔
426

427
  // Random ray mode does not have an inner loop over generations within a
428
  // batch, so set the current gen to 1
429
  simulation::current_gen = 1;
833✔
430

431
  switch (RandomRay::source_shape_) {
833!
432
  case RandomRaySourceShape::FLAT:
529✔
433
    domain_ = make_unique<FlatSourceDomain>();
529✔
434
    break;
529✔
435
  case RandomRaySourceShape::LINEAR:
304✔
436
  case RandomRaySourceShape::LINEAR_XY:
437
    domain_ = make_unique<LinearSourceDomain>();
304✔
438
    break;
304✔
439
  default:
×
440
    fatal_error("Unknown random ray source shape");
×
441
  }
442

443
  // Convert OpenMC native MGXS into a more efficient format
444
  // internal to the random ray solver
445
  domain_->flatten_xs();
833✔
446

447
  // Check if adjoint calculation is needed. If it is, we will run the forward
448
  // calculation first and then the adjoint calculation later.
449
  adjoint_needed_ = FlatSourceDomain::adjoint_;
833✔
450

451
  // Adjoint is always false for the forward calculation
452
  FlatSourceDomain::adjoint_ = false;
833✔
453

454
  // The first simulation is run after initialization
455
  is_first_simulation_ = true;
833✔
456

457
  // Initialize vectors used for baking in the initial condition during time
458
  // stepping
459
  if (settings::kinetic_simulation) {
833✔
460
    // Initialize vars used for time-consistent seed approach
461
    static_avg_k_eff_;
462
    static_k_eff_;
463
    static_fission_rate_;
464
  }
465
}
833✔
466

467
void RandomRaySimulation::apply_fixed_sources_and_mesh_domains()
833✔
468
{
469
  domain_->apply_meshes();
833✔
470
  if (settings::run_mode == RunMode::FIXED_SOURCE) {
833✔
471
    // Transfer external source user inputs onto random ray source regions
472
    domain_->convert_external_sources();
401✔
473
    domain_->count_external_source_regions();
401✔
474
  }
475
}
833✔
476

477
void RandomRaySimulation::prepare_fixed_sources_adjoint()
81✔
478
{
479
  domain_->source_regions_.adjoint_reset();
81✔
480
  if (settings::run_mode == RunMode::FIXED_SOURCE) {
81✔
481
    domain_->set_adjoint_sources();
65✔
482
  }
483
}
81✔
484

485
void RandomRaySimulation::print_random_ray_headers()
848✔
486
{
487
  openmc::print_random_ray_headers(adjoint_needed_);
848✔
488
}
848✔
489

490
void RandomRaySimulation::run_single_simulation()
2,066✔
491
{
492
  if (!is_first_simulation_) {
2,066✔
493
    if (mpi::master)
1,233✔
494
      print_random_ray_headers();
848✔
495

496
    // Reset the timers and reinitialize the general OpenMC datastructures if
497
    // this is after the first simulation
498
    reset_timers();
1,233✔
499

500
    // Initialize OpenMC general data structures
501
    openmc_simulation_init();
1,233✔
502
  }
503

504
  // Begin main simulation timer
505
  simulation::time_total.start();
2,066✔
506

507
  // Execute random ray simulation
508
  simulate();
2,066✔
509

510
  // End main simulation timer
511
  simulation::time_total.stop();
2,066✔
512

513
  // Normalize and save the final forward quantities
514
  domain_->normalize_final_quantities();
2,066✔
515

516
  // Finalize OpenMC
517
  openmc_simulation_finalize();
2,066✔
518

519
  // Output all simulation results
520
  output_simulation_results();
2,066✔
521

522
  // Toggle that the simulation object has been initialized after the first
523
  // simulation
524
  if (is_first_simulation_)
2,066✔
525
    is_first_simulation_ = false;
833✔
526
}
2,066✔
527

528
void RandomRaySimulation::random_ray_adjoint()
833✔
529
{
530
  if (!adjoint_needed_) {
833✔
531
    return;
752✔
532
  }
533

534
  // Configure the domain for adjoint simulation
535
  FlatSourceDomain::adjoint_ = true;
81✔
536

537
  // Reset k-eff
538
  domain_->k_eff_ = 1.0;
81✔
539

540
  // Initialize adjoint fixed sources, if present
541
  prepare_fixed_sources_adjoint();
81✔
542

543
  // Transpose scattering matrix
544
  domain_->transpose_scattering_matrix();
81✔
545

546
  // Swap nu_sigma_f and chi
547
  domain_->nu_sigma_f_.swap(domain_->chi_);
81✔
548

549
  // Run a single simulation
550
  run_single_simulation();
81✔
551
}
552

553
void RandomRaySimulation::kinetic_initial_condition()
192✔
554
{
555
  // Set flag for k_eff correction
556
  simulation::k_eff_correction = true;
192✔
557

558
  static_avg_k_eff_ = simulation::keff;
192✔
559
  domain_->k_eff_ = static_avg_k_eff_;
192✔
560
  domain_->source_regions_.adjoint_reset();
192✔
561
  domain_->propagate_final_quantities();
192✔
562
  domain_->source_regions_.time_step_reset();
192✔
563

564
  // Run the initial condition
565
  run_single_simulation();
192✔
566

567
  // Initialize the BD arrays
568
  domain_->store_time_step_quantities(false);
192✔
569

570
  // Store k-eff corrected initial condition statepoints
571
  rename_time_step_file(
540✔
572
    fmt::format("statepoint.{0}", settings::n_batches), ".h5", 0);
348✔
573
  if (settings::output_tallies) {
192!
574
    rename_time_step_file("tallies", ".out", 0);
192✔
575
  }
576

577
  // Set flags for kinetic simulation
578
  simulation::is_initial_condition = false;
192✔
579
  simulation::k_eff_correction = false;
192✔
580

581
  // Set starting time as zero
582
  simulation::current_time = 0;
192✔
583
}
192✔
584

585
// TODO: Add support for time-dependent restart
586
void RandomRaySimulation::kinetic_single_time_step(int i)
960✔
587
{
588
  // Increment time step
589
  simulation::current_timestep = i + 1;
960✔
590
  simulation::current_time += settings::dt;
960✔
591

592
  // Propogate results of previous simulation
593
  domain_->k_eff_ = static_avg_k_eff_;
960✔
594
  domain_->source_regions_.adjoint_reset();
960✔
595
  domain_->propagate_final_quantities();
960✔
596
  domain_->source_regions_.time_step_reset();
960✔
597

598
  // Compute RHS backward differences
599
  domain_->compute_rhs_bd_quantities();
960✔
600

601
  // Update time dependent cross section based on the density
602
  domain_->update_material_density(i);
960✔
603

604
  // Run the simulation for the current time step
605
  run_single_simulation();
960✔
606

607
  // Rename statepoint and tallies file for the current time step
608
  rename_time_step_file(fmt::format("statepoint.{0}", settings::n_batches),
1,920✔
609
    ".h5", simulation::current_timestep);
610
  if (settings::output_tallies) {
960!
611
    rename_time_step_file("tallies", ".out", simulation::current_timestep);
960✔
612
  }
613

614
  // Store final quantities for the current time step
615
  domain_->store_time_step_quantities();
960✔
616
}
960✔
617

618
void RandomRaySimulation::simulate()
2,066✔
619
{
620
  // Random ray power iteration loop
621
  while (simulation::current_batch < settings::n_batches) {
210,078✔
622
    // Initialize the current batch
623
    initialize_batch();
208,012✔
624
    initialize_generation();
208,012✔
625

626
    // MPI not supported in random ray solver, so all work is done by rank 0
627
    // TODO: Implement domain decomposition for MPI parallelism
628
    if (mpi::master) {
208,012✔
629

630
      // Reset total starting particle weight used for normalizing tallies
631
      simulation::total_weight = 1.0;
143,012✔
632

633
      // Update source term (scattering + fission (+ delayed if kinetic))
634
      domain_->update_all_neutron_sources();
143,012✔
635

636
      // Reset scalar fluxes, iteration volume tallies, and region hit flags
637
      // to zero
638
      domain_->batch_reset();
143,012✔
639

640
      // At the beginning of the simulation, if mesh subdivision is in use, we
641
      // need to swap the main source region container into the base container,
642
      // as the main source region container will be used to hold the true
643
      // subdivided source regions. The base container will therefore only
644
      // contain the external source region information, the mesh indices,
645
      // material properties, and initial guess values for the flux/source.
646

647
      // Start timer for transport
648
      simulation::time_transport.start();
143,012✔
649

650
// Transport sweep over all random rays for the iteration
651
#pragma omp parallel for schedule(dynamic)                                     \
78,012✔
652
  reduction(+ : total_geometric_intersections_)
78,012✔
653
      for (int i = 0; i < settings::n_particles; i++) {
6,883,500✔
654
        RandomRay ray(i, domain_.get());
6,818,500✔
655
        total_geometric_intersections_ +=
6,818,500✔
656
          ray.transport_history_based_single_ray();
6,818,500✔
657
      }
6,818,500✔
658

659
      simulation::time_transport.stop();
143,012✔
660

661
      // Add any newly discovered source regions to the main source region
662
      // container.
663
      domain_->finalize_discovered_source_regions();
143,012✔
664

665
      // Normalize scalar flux and update volumes
666
      domain_->normalize_scalar_flux_and_volumes(
143,012✔
667
        settings::n_particles * RandomRay::distance_active_);
668

669
      // Add source to scalar flux, compute number of FSR hits
670
      int64_t n_hits = domain_->add_source_to_scalar_flux();
143,012✔
671

672
      // Apply transport stabilization factors
673
      domain_->apply_transport_stabilization();
143,012✔
674

675
      if (settings::run_mode == RunMode::EIGENVALUE) {
143,012✔
676
        // Compute random ray k-eff
677
        if (!settings::kinetic_simulation ||
133,430!
678
            settings::kinetic_simulation && simulation::is_initial_condition) {
130,900✔
679
          domain_->compute_k_eff();
39,930✔
680
          if (simulation::k_eff_correction) {
39,930✔
681
            static_fission_rate_.push_back(domain_->fission_rate_);
18,700✔
682
            static_k_eff_.push_back(domain_->k_eff_);
18,700✔
683
          }
684
        } else {
685
          domain_->k_eff_ = static_k_eff_[simulation::current_batch - 1];
93,500✔
686
          domain_->fission_rate_ =
93,500✔
687
            static_fission_rate_[simulation::current_batch - 1];
93,500✔
688
        }
689

690
        // Store random ray k-eff into OpenMC's native k-eff variable
691
        global_tally_tracklength = domain_->k_eff_;
133,430✔
692
      }
693

694
      // Compute precursors if delayed neutrons are turned on
695
      if (settings::kinetic_simulation && settings::create_delayed_neutrons)
143,012!
696
        domain_->compute_all_precursors();
130,900✔
697

698
      // Execute all tallying tasks, if this is an active batch
699
      if (simulation::current_batch > settings::n_inactive) {
143,012✔
700

701
        // Add this iteration's scalar flux estimate to final accumulated
702
        // estimate
703
        domain_->accumulate_iteration_quantities();
69,636✔
704

705
        // Use above mapping to contribute FSR flux data to appropriate
706
        // tallies
707
        domain_->random_ray_tally();
69,636✔
708
      }
709

710
      // Set phi_old = phi_new
711
      domain_->flux_swap();
143,012✔
712
      if (settings::kinetic_simulation && settings::create_delayed_neutrons) {
143,012!
713
        domain_->precursors_swap();
130,900✔
714
      }
715

716
      // Check for any obvious insabilities/nans/infs
717
      instability_check(n_hits, domain_->k_eff_, avg_miss_rate_);
143,012✔
718
    } // End MPI master work
719

720
    // Store simulation metrics
721
    RandomRay::avg_miss_rate_ = avg_miss_rate_ / settings::n_batches;
208,012✔
722
    RandomRay::total_geometric_intersections_ = total_geometric_intersections_;
208,012✔
723
    RandomRay::n_external_source_regions_ = domain_->n_external_source_regions_;
208,012✔
724
    RandomRay::n_source_regions_ = domain_->n_source_regions();
208,012✔
725

726
    // Finalize the current batch
727
    finalize_generation();
208,012✔
728
    finalize_batch();
208,012✔
729
  } // End random ray power iteration loop
730

731
  domain_->count_external_source_regions();
2,066✔
732
}
2,066✔
733

734
void RandomRaySimulation::output_simulation_results() const
2,066✔
735
{
736
  // Print random ray results
737
  if (mpi::master) {
2,066✔
738
    print_results_random_ray();
1,421✔
739
    if (model::plots.size() > 0) {
1,421!
740
      domain_->output_to_vtk();
×
741
    }
742
  }
743
}
2,066✔
744

745
// Apply a few sanity checks to catch obvious cases of numerical instability.
746
// Instability typically only occurs if ray density is extremely low.
747
void RandomRaySimulation::instability_check(
143,012✔
748
  int64_t n_hits, double k_eff, double& avg_miss_rate) const
749
{
750
  double percent_missed = ((domain_->n_source_regions() - n_hits) /
143,012✔
751
                            static_cast<double>(domain_->n_source_regions())) *
143,012✔
752
                          100.0;
143,012✔
753
  avg_miss_rate += percent_missed;
143,012✔
754

755
  if (mpi::master) {
143,012!
756
    if (percent_missed > 10.0) {
143,012✔
757
      warning(fmt::format(
957✔
758
        "Very high FSR miss rate detected ({:.3f}%). Instability may occur. "
759
        "Increase ray density by adding more rays and/or active distance.",
760
        percent_missed));
761
    } else if (percent_missed > 1.0) {
142,055✔
762
      warning(
2!
763
        fmt::format("Elevated FSR miss rate detected ({:.3f}%). Increasing "
4!
764
                    "ray density by adding more rays and/or active "
765
                    "distance may improve simulation efficiency.",
766
          percent_missed));
767
    }
768

769
    if (k_eff > 10.0 || k_eff < 0.01 || !(std::isfinite(k_eff))) {
143,012!
770
      fatal_error(fmt::format("Instability detected: k-eff = {:.5f}", k_eff));
×
771
    }
772
  }
773
}
143,012✔
774

775
// Print random ray simulation results
776
void RandomRaySimulation::print_results_random_ray() const
1,421✔
777
{
778
  using namespace simulation;
779

780
  if (settings::verbosity >= 6) {
1,421!
781
    double total_integrations =
1,421✔
782
      RandomRay::total_geometric_intersections_ * negroups_;
1,421✔
783
    double time_per_integration =
784
      simulation::time_transport.elapsed() / total_integrations;
1,421✔
785
    double misc_time = time_total.elapsed() - time_update_src.elapsed() -
1,421✔
786
                       time_transport.elapsed() - time_tallies.elapsed() -
1,421✔
787
                       time_bank_sendrecv.elapsed();
1,421✔
788

789
    if (settings::kinetic_simulation && !simulation::is_initial_condition) {
1,421✔
790
      misc_time -= time_update_bd_vectors.elapsed();
660✔
791
    }
792
    header("Simulation Statistics", 4);
1,421✔
793
    fmt::print(
1,421✔
794
      " Total Iterations                  = {}\n", settings::n_batches);
795
    fmt::print(
1,421✔
796
      " Number of Rays per Iteration      = {}\n", settings::n_particles);
797
    fmt::print(" Inactive Distance                 = {} cm\n",
1,421✔
798
      RandomRay::distance_inactive_);
799
    fmt::print(" Active Distance                   = {} cm\n",
1,421✔
800
      RandomRay::distance_active_);
801
    fmt::print(" Source Regions (SRs)              = {}\n",
1,421✔
802
      RandomRay::n_source_regions_);
803
    fmt::print(" SRs Containing External Sources   = {}\n",
1,163✔
804
      RandomRay::n_external_source_regions_);
805
    fmt::print(" Total Geometric Intersections     = {:.4e}\n",
1,163✔
806
      static_cast<double>(RandomRay::total_geometric_intersections_));
1,421✔
807
    fmt::print("   Avg per Iteration               = {:.4e}\n",
1,163✔
808
      static_cast<double>(RandomRay::total_geometric_intersections_) /
1,421✔
809
        settings::n_batches);
810
    fmt::print("   Avg per Iteration per SR        = {:.2f}\n",
1,163✔
811
      static_cast<double>(RandomRay::total_geometric_intersections_) /
1,421✔
812
        static_cast<double>(settings::n_batches) /
2,842✔
813
        RandomRay::n_source_regions_);
814
    fmt::print(" Avg SR Miss Rate per Iteration    = {:.4f}%\n",
1,163✔
815
      RandomRay::avg_miss_rate_);
816
    fmt::print(" Energy Groups                     = {}\n", negroups_);
2,584✔
817
    if (settings::kinetic_simulation)
1,421✔
818
      fmt::print(" Delay Groups                      = {}\n", ndgroups_);
1,848✔
819
    fmt::print(
1,163✔
820
      " Total Integrations                = {:.4e}\n", total_integrations);
821
    fmt::print("   Avg per Iteration               = {:.4e}\n",
1,163✔
822
      total_integrations / settings::n_batches);
1,421✔
823

824
    std::string estimator;
1,421✔
825
    switch (domain_->volume_estimator_) {
1,421!
826
    case RandomRayVolumeEstimator::SIMULATION_AVERAGED:
22✔
827
      estimator = "Simulation Averaged";
22✔
828
      break;
22✔
829
    case RandomRayVolumeEstimator::NAIVE:
68✔
830
      estimator = "Naive";
68✔
831
      break;
68✔
832
    case RandomRayVolumeEstimator::HYBRID:
1,331✔
833
      estimator = "Hybrid";
1,331✔
834
      break;
1,331✔
835
    default:
×
836
      fatal_error("Invalid volume estimator type");
×
837
    }
838
    fmt::print(" Volume Estimator Type             = {}\n", estimator);
1,163✔
839

840
    std::string adjoint_true = (FlatSourceDomain::adjoint_) ? "ON" : "OFF";
4,263✔
841
    fmt::print(" Adjoint Flux Mode                 = {}\n", adjoint_true);
1,163✔
842

843
    std::string shape;
2,842✔
844
    switch (RandomRay::source_shape_) {
1,421!
845
    case RandomRaySourceShape::FLAT:
1,201✔
846
      shape = "Flat";
1,201✔
847
      break;
1,201✔
848
    case RandomRaySourceShape::LINEAR:
187✔
849
      shape = "Linear";
187✔
850
      break;
187✔
851
    case RandomRaySourceShape::LINEAR_XY:
33✔
852
      shape = "Linear XY";
33✔
853
      break;
33✔
854
    default:
×
855
      fatal_error("Invalid random ray source shape");
×
856
    }
857
    fmt::print(" Source Shape                      = {}\n", shape);
1,163✔
858
    std::string sample_method =
859
      (RandomRay::sample_method_ == RandomRaySampleMethod::PRNG) ? "PRNG"
1,421✔
860
                                                                 : "Halton";
4,263✔
861
    fmt::print(" Sample Method                     = {}\n", sample_method);
1,163✔
862

863
    if (domain_->is_transport_stabilization_needed_) {
1,421✔
864
      fmt::print(" Transport XS Stabilization Used   = YES (rho = {:.3f})\n",
165✔
865
        FlatSourceDomain::diagonal_stabilization_rho_);
866
    } else {
867
      fmt::print(" Transport XS Stabilization Used   = NO\n");
1,256✔
868
    }
869
    if (settings::kinetic_simulation && !simulation::is_initial_condition) {
1,421✔
870
      std::string time_method =
871
        (RandomRay::time_method_ == RandomRayTimeMethod::ISOTROPIC)
660✔
872
          ? "ISOTROPIC"
873
          : "PROPAGATION";
1,320✔
874
      fmt::print(" Time Method                       = {}\n", time_method);
660✔
875
      fmt::print(
540✔
876
        " Backwards Difference Order        = {}\n", RandomRay::bd_order_);
877
    }
660✔
878

879
    header("Timing Statistics", 4);
1,421✔
880
    show_time("Total time for initialization", time_initialize.elapsed());
1,421✔
881
    show_time("Reading cross sections", time_read_xs.elapsed(), 1);
1,421✔
882
    show_time("Total simulation time", time_total.elapsed());
1,421✔
883
    show_time("Transport sweep only", time_transport.elapsed(), 1);
1,421✔
884
    show_time("Source update only", time_update_src.elapsed(), 1);
1,421✔
885
    if (settings::kinetic_simulation && settings::create_delayed_neutrons) {
1,421!
886
      show_time(
924✔
887
        "Precursor computation only", time_compute_precursors.elapsed(), 1);
888
      misc_time -= time_compute_precursors.elapsed();
924✔
889
    }
890
    show_time("Tally conversion only", time_tallies.elapsed(), 1);
1,421✔
891
    show_time("MPI source reductions only", time_bank_sendrecv.elapsed(), 1);
1,421✔
892
    show_time("Other iteration routines", misc_time, 1);
1,421✔
893
    if (settings::run_mode == RunMode::EIGENVALUE) {
1,421✔
894
      show_time("Time in inactive batches", time_inactive.elapsed());
1,100✔
895
    }
896
    show_time("Time in active batches", time_active.elapsed());
1,421✔
897
    show_time("Time writing statepoints", time_statepoint.elapsed());
1,421✔
898
    show_time("Total time for finalization", time_finalize.elapsed());
1,421✔
899
    show_time("Time per integration", time_per_integration);
1,421✔
900
  }
1,421✔
901

902
  if (settings::verbosity >= 4 && settings::run_mode == RunMode::EIGENVALUE) {
1,421!
903
    header("Results", 4);
1,100✔
904
    fmt::print(" k-effective                       = {:.5f} +/- {:.5f}\n",
1,100✔
905
      simulation::keff, simulation::keff_std);
906
  }
907
}
1,421✔
908

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