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

openmc-dev / openmc / 16256055736

14 Jul 2025 01:49AM UTC coverage: 85.136% (-0.1%) from 85.251%
16256055736

Pull #3499

github

web-flow
Merge b0177d9f7 into d700d395d
Pull Request #3499: Add random ray info to statepoint file

84 of 122 new or added lines in 4 files covered. (68.85%)

85 existing lines in 15 files now uncovered.

52607 of 61792 relevant lines covered (85.14%)

36308708.86 hits per line

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

88.64
/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()
560✔
26
{
27
  //////////////////////////////////////////////////////////
28
  // Run forward simulation
29
  //////////////////////////////////////////////////////////
30

31
  // Check if adjoint calculation is needed. If it is, we will run the forward
32
  // calculation first and then the adjoint calculation later.
33
  bool adjoint_needed = FlatSourceDomain::adjoint_;
560✔
34

35
  // Configure the domain for forward simulation
36
  FlatSourceDomain::adjoint_ = false;
560✔
37

38
  // If we're going to do an adjoint simulation afterwards, report that this is
39
  // the initial forward flux solve.
40
  if (adjoint_needed && mpi::master)
560✔
41
    header("FORWARD FLUX SOLVE", 3);
55✔
42

43
  // Initialize OpenMC general data structures
44
  openmc_simulation_init();
560✔
45

46
  // Validate that inputs meet requirements for random ray mode
47
  if (mpi::master)
560✔
48
    validate_random_ray_inputs();
385✔
49

50
  // Declare forward flux so that it can be saved for later adjoint simulation
51
  vector<double> forward_flux;
560✔
52
  SourceRegionContainer forward_source_regions;
560✔
53
  SourceRegionContainer forward_base_source_regions;
560✔
54
  std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>
55
    forward_source_region_map;
560✔
56

57
  {
58
    // Initialize Random Ray Simulation Object
59
    RandomRaySimulation sim;
560✔
60

61
    // Initialize fixed sources, if present
62
    sim.apply_fixed_sources_and_mesh_domains();
560✔
63

64
    // Begin main simulation timer
65
    simulation::time_total.start();
560✔
66

67
    // Execute random ray simulation
68
    sim.simulate();
560✔
69

70
    // End main simulation timer
71
    simulation::time_total.stop();
560✔
72

73
    // Normalize and save the final forward flux
74
    sim.domain()->serialize_final_fluxes(forward_flux);
560✔
75

76
    double source_normalization_factor =
77
      sim.domain()->compute_fixed_source_normalization_factor() /
560✔
78
      (settings::n_batches - settings::n_inactive);
560✔
79

80
#pragma omp parallel for
315✔
81
    for (uint64_t i = 0; i < forward_flux.size(); i++) {
1,258,577✔
82
      forward_flux[i] *= source_normalization_factor;
1,258,332✔
83
    }
84

85
    forward_source_regions = sim.domain()->source_regions_;
560✔
86
    forward_source_region_map = sim.domain()->source_region_map_;
560✔
87
    forward_base_source_regions = sim.domain()->base_source_regions_;
560✔
88

89
    // Finalize OpenMC
90
    openmc_simulation_finalize();
560✔
91

92
    // Output all simulation results
93
    sim.output_simulation_results();
560✔
94
  }
560✔
95

96
  //////////////////////////////////////////////////////////
97
  // Run adjoint simulation (if enabled)
98
  //////////////////////////////////////////////////////////
99

100
  if (adjoint_needed) {
560✔
101
    reset_timers();
80✔
102

103
    // Configure the domain for adjoint simulation
104
    FlatSourceDomain::adjoint_ = true;
80✔
105

106
    if (mpi::master)
80✔
107
      header("ADJOINT FLUX SOLVE", 3);
55✔
108

109
    // Initialize OpenMC general data structures
110
    openmc_simulation_init();
80✔
111

112
    // Initialize Random Ray Simulation Object
113
    RandomRaySimulation adjoint_sim;
80✔
114

115
    // Initialize adjoint fixed sources, if present
116
    adjoint_sim.prepare_fixed_sources_adjoint(forward_flux,
80✔
117
      forward_source_regions, forward_base_source_regions,
118
      forward_source_region_map);
119

120
    // Transpose scattering matrix
121
    adjoint_sim.domain()->transpose_scattering_matrix();
80✔
122

123
    // Swap nu_sigma_f and chi
124
    adjoint_sim.domain()->nu_sigma_f_.swap(adjoint_sim.domain()->chi_);
80✔
125

126
    // Begin main simulation timer
127
    simulation::time_total.start();
80✔
128

129
    // Execute random ray simulation
130
    adjoint_sim.simulate();
80✔
131

132
    // End main simulation timer
133
    simulation::time_total.stop();
80✔
134

135
    // Finalize OpenMC
136
    openmc_simulation_finalize();
80✔
137

138
    // Output all simulation results
139
    adjoint_sim.output_simulation_results();
80✔
140
  }
80✔
141
}
560✔
142

143
// Enforces restrictions on inputs in random ray mode.  While there are
144
// many features that don't make sense in random ray mode, and are therefore
145
// unsupported, we limit our testing/enforcement operations only to inputs
146
// that may cause erroneous/misleading output or crashes from the solver.
147
void validate_random_ray_inputs()
385✔
148
{
149
  // Validate tallies
150
  ///////////////////////////////////////////////////////////////////
151
  for (auto& tally : model::tallies) {
1,232✔
152

153
    // Validate score types
154
    for (auto score_bin : tally->scores_) {
1,870✔
155
      switch (score_bin) {
1,023✔
156
      case SCORE_FLUX:
1,023✔
157
      case SCORE_TOTAL:
158
      case SCORE_FISSION:
159
      case SCORE_NU_FISSION:
160
      case SCORE_EVENTS:
161
        break;
1,023✔
162
      default:
×
163
        fatal_error(
×
164
          "Invalid score specified. Only flux, total, fission, nu-fission, and "
165
          "event scores are supported in random ray mode.");
166
      }
167
    }
168

169
    // Validate filter types
170
    for (auto f : tally->filters()) {
1,848✔
171
      auto& filter = *model::tally_filters[f];
1,001✔
172

173
      switch (filter.type()) {
1,001✔
174
      case FilterType::CELL:
1,001✔
175
      case FilterType::CELL_INSTANCE:
176
      case FilterType::DISTRIBCELL:
177
      case FilterType::ENERGY:
178
      case FilterType::MATERIAL:
179
      case FilterType::MESH:
180
      case FilterType::UNIVERSE:
181
      case FilterType::PARTICLE:
182
        break;
1,001✔
183
      default:
×
184
        fatal_error("Invalid filter specified. Only cell, cell_instance, "
×
185
                    "distribcell, energy, material, mesh, and universe filters "
186
                    "are supported in random ray mode.");
187
      }
188
    }
189
  }
190

191
  // Validate MGXS data
192
  ///////////////////////////////////////////////////////////////////
193
  for (auto& material : data::mg.macro_xs_) {
1,430✔
194
    if (!material.is_isotropic) {
1,045✔
195
      fatal_error("Anisotropic MGXS detected. Only isotropic XS data sets "
×
196
                  "supported in random ray mode.");
197
    }
198
    if (material.get_xsdata().size() > 1) {
1,045✔
199
      warning("Non-isothermal MGXS detected. Only isothermal XS data sets "
×
200
              "supported in random ray mode. Using lowest temperature.");
201
    }
202
    for (int g = 0; g < data::mg.num_energy_groups_; g++) {
5,522✔
203
      if (material.exists_in_model) {
4,477✔
204
        // Temperature and angle indices, if using multiple temperature
205
        // data sets and/or anisotropic data sets.
206
        // TODO: Currently assumes we are only using single temp/single angle
207
        // data.
208
        const int t = 0;
4,433✔
209
        const int a = 0;
4,433✔
210
        double sigma_t =
211
          material.get_xs(MgxsType::TOTAL, g, NULL, NULL, NULL, t, a);
4,433✔
212
        if (sigma_t <= 0.0) {
4,433✔
213
          fatal_error("No zero or negative total macroscopic cross sections "
×
214
                      "allowed in random ray mode. If the intention is to make "
215
                      "a void material, use a cell fill of 'None' instead.");
216
        }
217
      }
218
    }
219
  }
220

221
  // Validate ray source
222
  ///////////////////////////////////////////////////////////////////
223

224
  // Check for independent source
225
  IndependentSource* is =
226
    dynamic_cast<IndependentSource*>(RandomRay::ray_source_.get());
385✔
227
  if (!is) {
385✔
228
    fatal_error("Invalid ray source definition. Ray source must provided and "
×
229
                "be of type IndependentSource.");
230
  }
231

232
  // Check for box source
233
  SpatialDistribution* space_dist = is->space();
385✔
234
  SpatialBox* sb = dynamic_cast<SpatialBox*>(space_dist);
385✔
235
  if (!sb) {
385✔
236
    fatal_error(
×
237
      "Invalid ray source definition -- only box sources are allowed.");
238
  }
239

240
  // Check that box source is not restricted to fissionable areas
241
  if (sb->only_fissionable()) {
385✔
242
    fatal_error(
×
243
      "Invalid ray source definition -- fissionable spatial distribution "
244
      "not allowed.");
245
  }
246

247
  // Check for isotropic source
248
  UnitSphereDistribution* angle_dist = is->angle();
385✔
249
  Isotropic* id = dynamic_cast<Isotropic*>(angle_dist);
385✔
250
  if (!id) {
385✔
251
    fatal_error("Invalid ray source definition -- only isotropic sources are "
×
252
                "allowed.");
253
  }
254

255
  // Validate external sources
256
  ///////////////////////////////////////////////////////////////////
257
  if (settings::run_mode == RunMode::FIXED_SOURCE) {
385✔
258
    if (model::external_sources.size() < 1) {
264✔
259
      fatal_error("Must provide a particle source (in addition to ray source) "
×
260
                  "in fixed source random ray mode.");
261
    }
262

263
    for (int i = 0; i < model::external_sources.size(); i++) {
528✔
264
      Source* s = model::external_sources[i].get();
264✔
265

266
      // Check for independent source
267
      IndependentSource* is = dynamic_cast<IndependentSource*>(s);
264✔
268

269
      if (!is) {
264✔
270
        fatal_error(
×
271
          "Only IndependentSource external source types are allowed in "
272
          "random ray mode");
273
      }
274

275
      // Check for isotropic source
276
      UnitSphereDistribution* angle_dist = is->angle();
264✔
277
      Isotropic* id = dynamic_cast<Isotropic*>(angle_dist);
264✔
278
      if (!id) {
264✔
279
        fatal_error(
×
280
          "Invalid source definition -- only isotropic external sources are "
281
          "allowed in random ray mode.");
282
      }
283

284
      // Validate that a domain ID was specified OR that it is a point source
285
      auto sp = dynamic_cast<SpatialPoint*>(is->space());
264✔
286
      if (is->domain_ids().size() == 0 && !sp) {
264✔
287
        fatal_error("Fixed sources must be point source or spatially "
×
288
                    "constrained by domain id (cell, material, or universe) in "
289
                    "random ray mode.");
290
      } else if (is->domain_ids().size() > 0 && sp) {
264✔
291
        // If both a domain constraint and a non-default point source location
292
        // are specified, notify user that domain constraint takes precedence.
293
        if (sp->r().x == 0.0 && sp->r().y == 0.0 && sp->r().z == 0.0) {
231✔
294
          warning("Fixed source has both a domain constraint and a point "
231✔
295
                  "type spatial distribution. The domain constraint takes "
296
                  "precedence in random ray mode -- point source coordinate "
297
                  "will be ignored.");
298
        }
299
      }
300

301
      // Check that a discrete energy distribution was used
302
      Distribution* d = is->energy();
264✔
303
      Discrete* dd = dynamic_cast<Discrete*>(d);
264✔
304
      if (!dd) {
264✔
305
        fatal_error(
×
306
          "Only discrete (multigroup) energy distributions are allowed for "
307
          "external sources in random ray mode.");
308
      }
309
    }
310
  }
311

312
  // Validate plotting files
313
  ///////////////////////////////////////////////////////////////////
314
  for (int p = 0; p < model::plots.size(); p++) {
385✔
315

316
    // Get handle to OpenMC plot object
317
    const auto& openmc_plottable = model::plots[p];
×
318
    Plot* openmc_plot = dynamic_cast<Plot*>(openmc_plottable.get());
×
319

320
    // Random ray plots only support voxel plots
321
    if (!openmc_plot) {
×
322
      warning(fmt::format(
×
323
        "Plot {} will not be used for end of simulation data plotting -- only "
324
        "voxel plotting is allowed in random ray mode.",
325
        openmc_plottable->id()));
×
326
      continue;
×
327
    } else if (openmc_plot->type_ != Plot::PlotType::voxel) {
×
328
      warning(fmt::format(
×
329
        "Plot {} will not be used for end of simulation data plotting -- only "
330
        "voxel plotting is allowed in random ray mode.",
331
        openmc_plottable->id()));
×
332
      continue;
×
333
    }
334
  }
335

336
  // Warn about slow MPI domain replication, if detected
337
  ///////////////////////////////////////////////////////////////////
338
#ifdef OPENMC_MPI
339
  if (mpi::n_procs > 1) {
175✔
340
    warning(
175✔
341
      "MPI parallelism is not supported by the random ray solver. All work "
342
      "will be performed by rank 0. Domain decomposition may be implemented in "
343
      "the future to provide efficient MPI scaling.");
344
  }
345
#endif
346

347
  // Warn about instability resulting from linear sources in small regions
348
  // when generating weight windows with FW-CADIS and an overlaid mesh.
349
  ///////////////////////////////////////////////////////////////////
350
  if (RandomRay::source_shape_ == RandomRaySourceShape::LINEAR &&
132✔
351
      RandomRay::mesh_subdivision_enabled_ &&
517✔
352
      variance_reduction::weight_windows.size() > 0) {
66✔
353
    warning(
11✔
354
      "Linear sources may result in negative fluxes in small source regions "
355
      "generated by mesh subdivision. Negative sources may result in low "
356
      "quality FW-CADIS weight windows. We recommend you use flat source mode "
357
      "when generating weight windows with an overlaid mesh tally.");
358
  }
359
}
385✔
360

361
void openmc_reset_random_ray()
6,842✔
362
{
363
  FlatSourceDomain::volume_estimator_ = RandomRayVolumeEstimator::HYBRID;
6,842✔
364
  FlatSourceDomain::volume_normalized_flux_tallies_ = false;
6,842✔
365
  FlatSourceDomain::adjoint_ = false;
6,842✔
366
  FlatSourceDomain::mesh_domain_map_.clear();
6,842✔
367
  RandomRay::ray_source_.reset();
6,842✔
368
  RandomRay::source_shape_ = RandomRaySourceShape::FLAT;
6,842✔
369
  RandomRay::mesh_subdivision_enabled_ = false;
6,842✔
370
  RandomRay::sample_method_ = RandomRaySampleMethod::PRNG;
6,842✔
371
}
6,842✔
372

373
void write_random_ray_hdf5(hid_t group)
440✔
374
{
375
  switch (RandomRay::source_shape_) {
440✔
376
  case RandomRaySourceShape::FLAT:
264✔
377
    write_dataset(group, "source_shape", "flat");
264✔
378
    break;
264✔
379
  case RandomRaySourceShape::LINEAR:
143✔
380
    write_dataset(group, "source_shape", "linear");
143✔
381
    break;
143✔
382
  case RandomRaySourceShape::LINEAR_XY:
33✔
383
    write_dataset(group, "source_shape", "linear_xy");
33✔
384
    break;
33✔
NEW
385
  default:
×
NEW
386
    break;
×
387
  }
388
  std::string sample_method =
389
    (RandomRay::sample_method_ == RandomRaySampleMethod::PRNG) ? "prng"
440✔
390
                                                               : "halton";
880✔
391
  write_dataset(group, "sample_method", sample_method);
440✔
392
  switch (FlatSourceDomain::volume_estimator_) {
440✔
393
  case RandomRayVolumeEstimator::SIMULATION_AVERAGED:
22✔
394
    write_dataset(group, "volume_estimator", "simulation_averaged");
22✔
395
    break;
22✔
396
  case RandomRayVolumeEstimator::NAIVE:
66✔
397
    write_dataset(group, "volume_estimator", "naive");
66✔
398
    break;
66✔
399
  case RandomRayVolumeEstimator::HYBRID:
352✔
400
    write_dataset(group, "volume_estimator", "hybrid");
352✔
401
    break;
352✔
NEW
402
  default:
×
NEW
403
    break;
×
404
  }
405
  write_dataset(group, "distance_active", RandomRay::distance_active_);
440✔
406
  write_dataset(group, "distance_inactive", RandomRay::distance_inactive_);
440✔
407
  write_dataset(group, "adjoint_mode", FlatSourceDomain::adjoint_);
440✔
408
  write_dataset(group, "avg_miss_rate", RandomRay::avg_miss_rate);
440✔
409
  write_dataset(group, "n_source_regions", RandomRay::n_source_regions);
440✔
410
  write_dataset(
440✔
411
    group, "n_external_source_regions", RandomRay::n_external_source_regions);
412
  write_dataset(group, "n_geometric_intersections",
440✔
413
    RandomRay::total_geometric_intersections);
414
  int64_t n_integrations =
440✔
415
    RandomRay::total_geometric_intersections * RandomRay::negroups_;
440✔
416
  write_dataset(group, "n_integrations", n_integrations);
440✔
417
}
440✔
418

419
//==============================================================================
420
// RandomRaySimulation implementation
421
//==============================================================================
422

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

430
  // Random ray mode does not have an inner loop over generations within a
431
  // batch, so set the current gen to 1
432
  simulation::current_gen = 1;
640✔
433

434
  switch (RandomRay::source_shape_) {
640✔
435
  case RandomRaySourceShape::FLAT:
384✔
436
    domain_ = make_unique<FlatSourceDomain>();
384✔
437
    break;
384✔
438
  case RandomRaySourceShape::LINEAR:
256✔
439
  case RandomRaySourceShape::LINEAR_XY:
440
    domain_ = make_unique<LinearSourceDomain>();
256✔
441
    break;
256✔
442
  default:
×
443
    fatal_error("Unknown random ray source shape");
×
444
  }
445

446
  // Convert OpenMC native MGXS into a more efficient format
447
  // internal to the random ray solver
448
  domain_->flatten_xs();
640✔
449
}
640✔
450

451
void RandomRaySimulation::apply_fixed_sources_and_mesh_domains()
560✔
452
{
453
  domain_->apply_meshes();
560✔
454
  if (settings::run_mode == RunMode::FIXED_SOURCE) {
560✔
455
    // Transfer external source user inputs onto random ray source regions
456
    domain_->convert_external_sources();
384✔
457
    domain_->count_external_source_regions();
384✔
458
  }
459
}
560✔
460

461
void RandomRaySimulation::prepare_fixed_sources_adjoint(
80✔
462
  vector<double>& forward_flux, SourceRegionContainer& forward_source_regions,
463
  SourceRegionContainer& forward_base_source_regions,
464
  std::unordered_map<SourceRegionKey, int64_t, SourceRegionKey::HashFunctor>&
465
    forward_source_region_map)
466
{
467
  if (settings::run_mode == RunMode::FIXED_SOURCE) {
80✔
468
    if (RandomRay::mesh_subdivision_enabled_) {
64✔
469
      domain_->source_regions_ = forward_source_regions;
32✔
470
      domain_->source_region_map_ = forward_source_region_map;
32✔
471
      domain_->base_source_regions_ = forward_base_source_regions;
32✔
472
      domain_->source_regions_.adjoint_reset();
32✔
473
    }
474
    domain_->set_adjoint_sources(forward_flux);
64✔
475
  }
476
}
80✔
477

478
void RandomRaySimulation::simulate()
640✔
479
{
480
  // Random ray power iteration loop
481
  while (simulation::current_batch < settings::n_batches) {
17,440✔
482
    // Initialize the current batch
483
    initialize_batch();
16,800✔
484
    initialize_generation();
16,800✔
485

486
    // MPI not supported in random ray solver, so all work is done by rank 0
487
    // TODO: Implement domain decomposition for MPI parallelism
488
    if (mpi::master) {
16,800✔
489

490
      // Reset total starting particle weight used for normalizing tallies
491
      simulation::total_weight = 1.0;
11,550✔
492

493
      // Update source term (scattering + fission)
494
      domain_->update_neutron_source(k_eff_);
11,550✔
495

496
      // Reset scalar fluxes, iteration volume tallies, and region hit flags to
497
      // zero
498
      domain_->batch_reset();
11,550✔
499

500
      // At the beginning of the simulation, if mesh subvivision is in use, we
501
      // need to swap the main source region container into the base container,
502
      // as the main source region container will be used to hold the true
503
      // subdivided source regions. The base container will therefore only
504
      // contain the external source region information, the mesh indices,
505
      // material properties, and initial guess values for the flux/source.
506
      if (RandomRay::mesh_subdivision_enabled_ &&
11,550✔
507
          simulation::current_batch == 1 && !FlatSourceDomain::adjoint_) {
2,970✔
508
        domain_->prepare_base_source_regions();
110✔
509
      }
510

511
      // Start timer for transport
512
      simulation::time_transport.start();
11,550✔
513

514
// Transport sweep over all random rays for the iteration
515
#pragma omp parallel for schedule(dynamic)                                     \
6,300✔
516
  reduction(+ : total_geometric_intersections_)
6,300✔
517
      for (int i = 0; i < settings::n_particles; i++) {
849,250✔
518
        RandomRay ray(i, domain_.get());
844,000✔
519
        total_geometric_intersections_ +=
844,000✔
520
          ray.transport_history_based_single_ray();
844,000✔
521
      }
844,000✔
522

523
      simulation::time_transport.stop();
11,550✔
524

525
      // If using mesh subdivision, add any newly discovered source regions
526
      // to the main source region container.
527
      if (RandomRay::mesh_subdivision_enabled_) {
11,550✔
528
        domain_->finalize_discovered_source_regions();
2,970✔
529
      }
530

531
      // Normalize scalar flux and update volumes
532
      domain_->normalize_scalar_flux_and_volumes(
11,550✔
533
        settings::n_particles * RandomRay::distance_active_);
534

535
      // Add source to scalar flux, compute number of FSR hits
536
      int64_t n_hits = domain_->add_source_to_scalar_flux();
11,550✔
537

538
      // Apply transport stabilization factors
539
      domain_->apply_transport_stabilization();
11,550✔
540

541
      if (settings::run_mode == RunMode::EIGENVALUE) {
11,550✔
542
        // Compute random ray k-eff
543
        k_eff_ = domain_->compute_k_eff(k_eff_);
2,090✔
544

545
        // Store random ray k-eff into OpenMC's native k-eff variable
546
        global_tally_tracklength = k_eff_;
2,090✔
547
      }
548

549
      // Execute all tallying tasks, if this is an active batch
550
      if (simulation::current_batch > settings::n_inactive) {
11,550✔
551

552
        // Add this iteration's scalar flux estimate to final accumulated
553
        // estimate
554
        domain_->accumulate_iteration_flux();
4,675✔
555

556
        // Generate mapping between source regions and tallies
557
        if (!domain_->mapped_all_tallies_ &&
5,203✔
558
            !RandomRay::mesh_subdivision_enabled_) {
528✔
559
          domain_->convert_source_regions_to_tallies(0);
308✔
560
        }
561

562
        // Use above mapping to contribute FSR flux data to appropriate
563
        // tallies
564
        domain_->random_ray_tally();
4,675✔
565
      }
566

567
      // Set phi_old = phi_new
568
      domain_->flux_swap();
11,550✔
569

570
      // Check for any obvious insabilities/nans/infs
571
      instability_check(n_hits, k_eff_, avg_miss_rate_);
11,550✔
572
    } // End MPI master work
573

574
    // Set global variables for reporting
575
    RandomRay::avg_miss_rate = avg_miss_rate_ / settings::n_batches;
16,800✔
576
    RandomRay::total_geometric_intersections = total_geometric_intersections_;
16,800✔
577
    RandomRay::n_external_source_regions = domain_->n_external_source_regions_;
16,800✔
578
    RandomRay::n_source_regions = domain_->n_source_regions();
16,800✔
579

580
    // Finalize the current batch
581
    finalize_generation();
16,800✔
582
    finalize_batch();
16,800✔
583
  } // End random ray power iteration loop
584

585
  domain_->count_external_source_regions();
640✔
586
}
640✔
587

588
void RandomRaySimulation::output_simulation_results() const
640✔
589
{
590
  // Print random ray results
591
  if (mpi::master) {
640✔
592
    print_results_random_ray();
440✔
593
    if (model::plots.size() > 0) {
440✔
594
      domain_->output_to_vtk();
×
595
    }
596
  }
597
}
640✔
598

599
// Apply a few sanity checks to catch obvious cases of numerical instability.
600
// Instability typically only occurs if ray density is extremely low.
601
void RandomRaySimulation::instability_check(
11,550✔
602
  int64_t n_hits, double k_eff, double& avg_miss_rate) const
603
{
604
  double percent_missed = ((domain_->n_source_regions() - n_hits) /
11,550✔
605
                            static_cast<double>(domain_->n_source_regions())) *
11,550✔
606
                          100.0;
11,550✔
607
  avg_miss_rate += percent_missed;
11,550✔
608

609
  if (mpi::master) {
11,550✔
610
    if (percent_missed > 10.0) {
11,550✔
611
      warning(fmt::format(
957✔
612
        "Very high FSR miss rate detected ({:.3f}%). Instability may occur. "
613
        "Increase ray density by adding more rays and/or active distance.",
614
        percent_missed));
615
    } else if (percent_missed > 1.0) {
10,593✔
616
      warning(
×
617
        fmt::format("Elevated FSR miss rate detected ({:.3f}%). Increasing "
×
618
                    "ray density by adding more rays and/or active "
619
                    "distance may improve simulation efficiency.",
620
          percent_missed));
621
    }
622

623
    if (k_eff > 10.0 || k_eff < 0.01 || !(std::isfinite(k_eff))) {
11,550✔
624
      fatal_error("Instability detected");
×
625
    }
626
  }
627
}
11,550✔
628

629
// Print random ray simulation results
630
void RandomRaySimulation::print_results_random_ray() const
440✔
631
{
632
  using namespace simulation;
633

634
  if (settings::verbosity >= 6) {
440✔
635
    double total_integrations =
440✔
636
      RandomRay::total_geometric_intersections * RandomRay::negroups_;
440✔
637
    double time_per_integration =
638
      simulation::time_transport.elapsed() / total_integrations;
440✔
639
    double misc_time = time_total.elapsed() - time_update_src.elapsed() -
440✔
640
                       time_transport.elapsed() - time_tallies.elapsed() -
440✔
641
                       time_bank_sendrecv.elapsed();
440✔
642

643
    header("Simulation Statistics", 4);
440✔
644
    fmt::print(
440✔
645
      " Total Iterations                  = {}\n", settings::n_batches);
646
    fmt::print(
440✔
647
      " Number of Rays per Iteration      = {}\n", settings::n_particles);
648
    fmt::print(" Inactive Distance                 = {} cm\n",
440✔
649
      RandomRay::distance_inactive_);
650
    fmt::print(" Active Distance                   = {} cm\n",
440✔
651
      RandomRay::distance_active_);
652
    fmt::print(
440✔
653
      " Source Regions (SRs)              = {}\n", RandomRay::n_source_regions);
654
    fmt::print(" SRs Containing External Sources   = {}\n",
360✔
655
      RandomRay::n_external_source_regions);
656
    fmt::print(" Total Geometric Intersections     = {:.4e}\n",
360✔
657
      static_cast<double>(RandomRay::total_geometric_intersections));
440✔
658
    fmt::print("   Avg per Iteration               = {:.4e}\n",
360✔
659
      static_cast<double>(RandomRay::total_geometric_intersections) /
440✔
660
        settings::n_batches);
661
    fmt::print("   Avg per Iteration per SR        = {:.2f}\n",
360✔
662
      static_cast<double>(RandomRay::total_geometric_intersections) /
440✔
663
        static_cast<double>(settings::n_batches) / RandomRay::n_source_regions);
880✔
664
    fmt::print(" Avg SR Miss Rate per Iteration    = {:.4f}%\n",
440✔
665
      RandomRay::avg_miss_rate);
666
    fmt::print(
440✔
667
      " Energy Groups                     = {}\n", RandomRay::negroups_);
668
    fmt::print(
360✔
669
      " Total Integrations                = {:.4e}\n", total_integrations);
670
    fmt::print("   Avg per Iteration               = {:.4e}\n",
360✔
671
      total_integrations / settings::n_batches);
440✔
672

673
    std::string estimator;
440✔
674
    switch (domain_->volume_estimator_) {
440✔
675
    case RandomRayVolumeEstimator::SIMULATION_AVERAGED:
22✔
676
      estimator = "Simulation Averaged";
22✔
677
      break;
22✔
678
    case RandomRayVolumeEstimator::NAIVE:
66✔
679
      estimator = "Naive";
66✔
680
      break;
66✔
681
    case RandomRayVolumeEstimator::HYBRID:
352✔
682
      estimator = "Hybrid";
352✔
683
      break;
352✔
684
    default:
×
685
      fatal_error("Invalid volume estimator type");
×
686
    }
687
    fmt::print(" Volume Estimator Type             = {}\n", estimator);
360✔
688

689
    std::string adjoint_true = (FlatSourceDomain::adjoint_) ? "ON" : "OFF";
1,320✔
690
    fmt::print(" Adjoint Flux Mode                 = {}\n", adjoint_true);
360✔
691

692
    std::string shape;
880✔
693
    switch (RandomRay::source_shape_) {
440✔
694
    case RandomRaySourceShape::FLAT:
264✔
695
      shape = "Flat";
264✔
696
      break;
264✔
697
    case RandomRaySourceShape::LINEAR:
143✔
698
      shape = "Linear";
143✔
699
      break;
143✔
700
    case RandomRaySourceShape::LINEAR_XY:
33✔
701
      shape = "Linear XY";
33✔
702
      break;
33✔
703
    default:
×
704
      fatal_error("Invalid random ray source shape");
×
705
    }
706
    fmt::print(" Source Shape                      = {}\n", shape);
360✔
707
    std::string sample_method =
708
      (RandomRay::sample_method_ == RandomRaySampleMethod::PRNG) ? "PRNG"
440✔
709
                                                                 : "Halton";
1,320✔
710
    fmt::print(" Sample Method                     = {}\n", sample_method);
360✔
711

712
    if (domain_->is_transport_stabilization_needed_) {
440✔
713
      fmt::print(" Transport XS Stabilization Used   = YES (rho = {:.3f})\n",
11✔
714
        FlatSourceDomain::diagonal_stabilization_rho_);
715
    } else {
716
      fmt::print(" Transport XS Stabilization Used   = NO\n");
429✔
717
    }
718

719
    header("Timing Statistics", 4);
440✔
720
    show_time("Total time for initialization", time_initialize.elapsed());
440✔
721
    show_time("Reading cross sections", time_read_xs.elapsed(), 1);
440✔
722
    show_time("Total simulation time", time_total.elapsed());
440✔
723
    show_time("Transport sweep only", time_transport.elapsed(), 1);
440✔
724
    show_time("Source update only", time_update_src.elapsed(), 1);
440✔
725
    show_time("Tally conversion only", time_tallies.elapsed(), 1);
440✔
726
    show_time("MPI source reductions only", time_bank_sendrecv.elapsed(), 1);
440✔
727
    show_time("Other iteration routines", misc_time, 1);
440✔
728
    if (settings::run_mode == RunMode::EIGENVALUE) {
440✔
729
      show_time("Time in inactive batches", time_inactive.elapsed());
132✔
730
    }
731
    show_time("Time in active batches", time_active.elapsed());
440✔
732
    show_time("Time writing statepoints", time_statepoint.elapsed());
440✔
733
    show_time("Total time for finalization", time_finalize.elapsed());
440✔
734
    show_time("Time per integration", time_per_integration);
440✔
735
  }
440✔
736

737
  if (settings::verbosity >= 4 && settings::run_mode == RunMode::EIGENVALUE) {
440✔
738
    header("Results", 4);
132✔
739
    fmt::print(" k-effective                       = {:.5f} +/- {:.5f}\n",
132✔
740
      simulation::keff, simulation::keff_std);
741
  }
742
}
440✔
743

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