10586562087

Committed 27 Aug 2024 10:05PM UTC coverage: 84.707% (-0.2%) from 84.9%

Build # 10586562087

Build Type

Pull #3112

github

Committed by

web-flow

Commit Message

Merge f7f32bf18 into 5bc04b5d7

Pull Request Pull Request #3112: Revamp CI with dependency and Python caching for efficient installs

Run Details

49553 of 58499 relevant lines covered (84.71%)

34324762.08 hits per line

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

83.17

/src/random_ray/random_ray_simulation.cpp

#include "openmc/random_ray/random_ray_simulation.h"

#include "openmc/eigenvalue.h"
#include "openmc/geometry.h"
#include "openmc/message_passing.h"
#include "openmc/mgxs_interface.h"
#include "openmc/output.h"
#include "openmc/plot.h"
#include "openmc/random_ray/flat_source_domain.h"
#include "openmc/random_ray/random_ray.h"
#include "openmc/simulation.h"
#include "openmc/source.h"
#include "openmc/tallies/filter.h"
#include "openmc/tallies/tally.h"
#include "openmc/tallies/tally_scoring.h"
#include "openmc/timer.h"

namespace openmc {

//==============================================================================
// Non-member functions
//==============================================================================

void openmc_run_random_ray()
{
  // Initialize OpenMC general data structures
  openmc_simulation_init();

  // Validate that inputs meet requirements for random ray mode
  if (mpi::master)
    validate_random_ray_inputs();

  // Initialize Random Ray Simulation Object
  RandomRaySimulation sim;

  // Begin main simulation timer
  simulation::time_total.start();

  // Execute random ray simulation
  sim.simulate();

  // End main simulation timer
  openmc::simulation::time_total.stop();

  // Finalize OpenMC
  openmc_simulation_finalize();

  // Reduce variables across MPI ranks
  sim.reduce_simulation_statistics();

  // Output all simulation results
  sim.output_simulation_results();
}

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

    // Validate score types
    for (auto score_bin : tally->scores_) {
      switch (score_bin) {
      case SCORE_FLUX:
      case SCORE_TOTAL:
      case SCORE_FISSION:
      case SCORE_NU_FISSION:
      case SCORE_EVENTS:
        break;
      default:
        fatal_error(
          "Invalid score specified. Only flux, total, fission, nu-fission, and "
          "event scores are supported in random ray mode.");
      }
    }

    // Validate filter types
    for (auto f : tally->filters()) {
      auto& filter = *model::tally_filters[f];

      switch (filter.type()) {
      case FilterType::CELL:
      case FilterType::CELL_INSTANCE:
      case FilterType::DISTRIBCELL:
      case FilterType::ENERGY:
      case FilterType::MATERIAL:
      case FilterType::MESH:
      case FilterType::UNIVERSE:
        break;
      default:
        fatal_error("Invalid filter specified. Only cell, cell_instance, "
                    "distribcell, energy, material, mesh, and universe filters "
                    "are supported in random ray mode.");
      }
    }
  }

  // Validate MGXS data
  ///////////////////////////////////////////////////////////////////
  for (auto& material : data::mg.macro_xs_) {
    if (!material.is_isotropic) {
      fatal_error("Anisotropic MGXS detected. Only isotropic XS data sets "
                  "supported in random ray mode.");
    }
    if (material.get_xsdata().size() > 1) {
      fatal_error("Non-isothermal MGXS detected. Only isothermal XS data sets "
                  "supported in random ray mode.");
    }
  }

  // Validate ray source
  ///////////////////////////////////////////////////////////////////

  // Check for independent source
  IndependentSource* is =
    dynamic_cast<IndependentSource*>(RandomRay::ray_source_.get());
  if (!is) {
    fatal_error("Invalid ray source definition. Ray source must provided and "
                "be of type IndependentSource.");
  }

  // Check for box source
  SpatialDistribution* space_dist = is->space();
  SpatialBox* sb = dynamic_cast<SpatialBox*>(space_dist);
  if (!sb) {
    fatal_error(
      "Invalid ray source definition -- only box sources are allowed.");
  }

  // Check that box source is not restricted to fissionable areas
  if (sb->only_fissionable()) {
    fatal_error(
      "Invalid ray source definition -- fissionable spatial distribution "
      "not allowed.");
  }

  // Check for isotropic source
  UnitSphereDistribution* angle_dist = is->angle();
  Isotropic* id = dynamic_cast<Isotropic*>(angle_dist);
  if (!id) {
    fatal_error("Invalid ray source definition -- only isotropic sources are "
                "allowed.");
  }

  // Validate external sources
  ///////////////////////////////////////////////////////////////////
  if (settings::run_mode == RunMode::FIXED_SOURCE) {
    if (model::external_sources.size() < 1) {
      fatal_error("Must provide a particle source (in addition to ray source) "
                  "in fixed source random ray mode.");
    }

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

      // Check for independent source
      IndependentSource* is = dynamic_cast<IndependentSource*>(s);

      if (!is) {
        fatal_error(
          "Only IndependentSource external source types are allowed in "
          "random ray mode");
      }

      // Check for isotropic source
      UnitSphereDistribution* angle_dist = is->angle();
      Isotropic* id = dynamic_cast<Isotropic*>(angle_dist);
      if (!id) {
        fatal_error(
          "Invalid source definition -- only isotropic external sources are "
          "allowed in random ray mode.");
      }

      // Validate that a domain ID was specified
      if (is->domain_ids().size() == 0) {
        fatal_error("Fixed sources must be specified by domain "
                    "id (cell, material, or universe) in random ray mode.");
      }

      // Check that a discrete energy distribution was used
      Distribution* d = is->energy();
      Discrete* dd = dynamic_cast<Discrete*>(d);
      if (!dd) {
        fatal_error(
          "Only discrete (multigroup) energy distributions are allowed for "
          "external sources in random ray mode.");
      }
    }
  }

  // Validate plotting files
  ///////////////////////////////////////////////////////////////////
  for (int p = 0; p < model::plots.size(); p++) {

    // Get handle to OpenMC plot object
    Plot* openmc_plot = dynamic_cast<Plot*>(model::plots[p].get());

    // Random ray plots only support voxel plots
    if (!openmc_plot) {
      warning(fmt::format(
        "Plot {} will not be used for end of simulation data plotting -- only "
        "voxel plotting is allowed in random ray mode.",
        p));
      continue;
    } else if (openmc_plot->type_ != Plot::PlotType::voxel) {
      warning(fmt::format(
        "Plot {} will not be used for end of simulation data plotting -- only "
        "voxel plotting is allowed in random ray mode.",
        p));
      continue;
    }
  }

  // Warn about slow MPI domain replication, if detected
  ///////////////////////////////////////////////////////////////////
#ifdef OPENMC_MPI
  if (mpi::n_procs > 1) {
    warning(
      "Domain replication in random ray is supported, but suffers from poor "
      "scaling of source all-reduce operations. Performance may severely "
      "degrade beyond just a few MPI ranks. Domain decomposition may be "
      "implemented in the future to provide efficient scaling.");
  }
#endif
}

//==============================================================================
// RandomRaySimulation implementation
//==============================================================================

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

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

  switch (RandomRay::source_shape_) {
  case RandomRaySourceShape::FLAT:
    domain_ = make_unique<FlatSourceDomain>();
    break;
  case RandomRaySourceShape::LINEAR:
  case RandomRaySourceShape::LINEAR_XY:
    domain_ = make_unique<LinearSourceDomain>();
    break;
  default:
    fatal_error("Unknown random ray source shape");
  }
}

void RandomRaySimulation::simulate()
{
  if (settings::run_mode == RunMode::FIXED_SOURCE) {
    // Transfer external source user inputs onto random ray source regions
    domain_->convert_external_sources();
    domain_->count_external_source_regions();
  }

  // Random ray power iteration loop
  while (simulation::current_batch < settings::n_batches) {

    // Initialize the current batch
    initialize_batch();
    initialize_generation();

    // Reset total starting particle weight used for normalizing tallies
    simulation::total_weight = 1.0;

    // Update source term (scattering + fission)
    domain_->update_neutron_source(k_eff_);

    // Reset scalar fluxes, iteration volume tallies, and region hit flags to
    // zero
    domain_->batch_reset();

    // Start timer for transport
    simulation::time_transport.start();

// Transport sweep over all random rays for the iteration
#pragma omp parallel for schedule(dynamic)                                     \
  reduction(+ : total_geometric_intersections_)
    for (int i = 0; i < simulation::work_per_rank; i++) {
      RandomRay ray(i, domain_.get());
      total_geometric_intersections_ +=
        ray.transport_history_based_single_ray();
    }

    simulation::time_transport.stop();

    // If using multiple MPI ranks, perform all reduce on all transport results
    domain_->all_reduce_replicated_source_regions();

    // Normalize scalar flux and update volumes
    domain_->normalize_scalar_flux_and_volumes(
      settings::n_particles * RandomRay::distance_active_);

    // Add source to scalar flux, compute number of FSR hits
    int64_t n_hits = domain_->add_source_to_scalar_flux();

    if (settings::run_mode == RunMode::EIGENVALUE) {
      // Compute random ray k-eff
      k_eff_ = domain_->compute_k_eff(k_eff_);

      // Store random ray k-eff into OpenMC's native k-eff variable
      global_tally_tracklength = k_eff_;
    }

    // Execute all tallying tasks, if this is an active batch
    if (simulation::current_batch > settings::n_inactive && mpi::master) {

      // Generate mapping between source regions and tallies
      if (!domain_->mapped_all_tallies_) {
        domain_->convert_source_regions_to_tallies();
      }

      // Use above mapping to contribute FSR flux data to appropriate tallies
      domain_->random_ray_tally();

      // Add this iteration's scalar flux estimate to final accumulated estimate
      domain_->accumulate_iteration_flux();
    }

    // Set phi_old = phi_new
    domain_->flux_swap();

    // Check for any obvious insabilities/nans/infs
    instability_check(n_hits, k_eff_, avg_miss_rate_);

    // Finalize the current batch
    finalize_generation();
    finalize_batch();
  } // End random ray power iteration loop
}

void RandomRaySimulation::reduce_simulation_statistics()
{
  // Reduce number of intersections
#ifdef OPENMC_MPI
  if (mpi::n_procs > 1) {
    uint64_t total_geometric_intersections_reduced = 0;
    MPI_Reduce(&total_geometric_intersections_,
      &total_geometric_intersections_reduced, 1, MPI_UNSIGNED_LONG, MPI_SUM, 0,
      mpi::intracomm);
    total_geometric_intersections_ = total_geometric_intersections_reduced;
  }
#endif
}

void RandomRaySimulation::output_simulation_results() const
{
  // Print random ray results
  if (mpi::master) {
    print_results_random_ray(total_geometric_intersections_,
      avg_miss_rate_ / settings::n_batches, negroups_,
      domain_->n_source_regions_, domain_->n_external_source_regions_);
    if (model::plots.size() > 0) {
      domain_->output_to_vtk();
    }
  }
}

// Apply a few sanity checks to catch obvious cases of numerical instability.
// Instability typically only occurs if ray density is extremely low.
void RandomRaySimulation::instability_check(
  int64_t n_hits, double k_eff, double& avg_miss_rate) const
{
  double percent_missed = ((domain_->n_source_regions_ - n_hits) /
                            static_cast<double>(domain_->n_source_regions_)) *
                          100.0;
  avg_miss_rate += percent_missed;

  if (mpi::master) {
    if (percent_missed > 10.0) {
      warning(fmt::format(
        "Very high FSR miss rate detected ({:.3f}%). Instability may occur. "
        "Increase ray density by adding more rays and/or active distance.",
        percent_missed));
    } else if (percent_missed > 1.0) {
      warning(
        fmt::format("Elevated FSR miss rate detected ({:.3f}%). Increasing "
                    "ray density by adding more rays and/or active "
                    "distance may improve simulation efficiency.",
          percent_missed));
    }

    if (k_eff > 10.0 || k_eff < 0.01 || !(std::isfinite(k_eff))) {
      fatal_error("Instability detected");
    }
  }
}

// Print random ray simulation results
void RandomRaySimulation::print_results_random_ray(
  uint64_t total_geometric_intersections, double avg_miss_rate, int negroups,
  int64_t n_source_regions, int64_t n_external_source_regions) const
{
  using namespace simulation;

  if (settings::verbosity >= 6) {
    double total_integrations = total_geometric_intersections * negroups;
    double time_per_integration =
      simulation::time_transport.elapsed() / total_integrations;
    double misc_time = time_total.elapsed() - time_update_src.elapsed() -
                       time_transport.elapsed() - time_tallies.elapsed() -
                       time_bank_sendrecv.elapsed();

    header("Simulation Statistics", 4);
    fmt::print(
      " Total Iterations                  = {}\n", settings::n_batches);
    fmt::print(" Flat Source Regions (FSRs)        = {}\n", n_source_regions);
    fmt::print(
      " FSRs Containing External Sources  = {}\n", n_external_source_regions);
    fmt::print(" Total Geometric Intersections     = {:.4e}\n",
      static_cast<double>(total_geometric_intersections));
    fmt::print("   Avg per Iteration               = {:.4e}\n",
      static_cast<double>(total_geometric_intersections) / settings::n_batches);
    fmt::print("   Avg per Iteration per FSR       = {:.2f}\n",
      static_cast<double>(total_geometric_intersections) /
        static_cast<double>(settings::n_batches) / n_source_regions);
    fmt::print(" Avg FSR Miss Rate per Iteration   = {:.4f}%\n", avg_miss_rate);
    fmt::print(" Energy Groups                     = {}\n", negroups);
    fmt::print(
      " Total Integrations                = {:.4e}\n", total_integrations);
    fmt::print("   Avg per Iteration               = {:.4e}\n",
      total_integrations / settings::n_batches);

    std::string estimator;
    switch (domain_->volume_estimator_) {
    case RandomRayVolumeEstimator::SIMULATION_AVERAGED:
      estimator = "Simulation Averaged";
      break;
    case RandomRayVolumeEstimator::NAIVE:
      estimator = "Naive";
      break;
    case RandomRayVolumeEstimator::HYBRID:
      estimator = "Hybrid";
      break;
    default:
      fatal_error("Invalid volume estimator type");
    }
    fmt::print(" Volume Estimator Type             = {}\n", estimator);

    header("Timing Statistics", 4);
    show_time("Total time for initialization", time_initialize.elapsed());
    show_time("Reading cross sections", time_read_xs.elapsed(), 1);
    show_time("Total simulation time", time_total.elapsed());
    show_time("Transport sweep only", time_transport.elapsed(), 1);
    show_time("Source update only", time_update_src.elapsed(), 1);
    show_time("Tally conversion only", time_tallies.elapsed(), 1);
    show_time("MPI source reductions only", time_bank_sendrecv.elapsed(), 1);
    show_time("Other iteration routines", misc_time, 1);
    if (settings::run_mode == RunMode::EIGENVALUE) {
      show_time("Time in inactive batches", time_inactive.elapsed());
    }
    show_time("Time in active batches", time_active.elapsed());
    show_time("Time writing statepoints", time_statepoint.elapsed());
    show_time("Total time for finalization", time_finalize.elapsed());
    show_time("Time per integration", time_per_integration);
  }

  if (settings::verbosity >= 4 && settings::run_mode == RunMode::EIGENVALUE) {
    header("Results", 4);
    fmt::print(" k-effective                       = {:.5f} +/- {:.5f}\n",
      simulation::keff, simulation::keff_std);
  }
}

} // namespace openmc

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
18	namespace openmc {
19
20	//==============================================================================
21	// Non-member functions
22	//==============================================================================
23
24	void openmc_run_random_ray()	323✔
25	{
26	// Initialize OpenMC general data structures
27	openmc_simulation_init();	323✔
28
29	// Validate that inputs meet requirements for random ray mode
30	if (mpi::master)	323✔
31	validate_random_ray_inputs();	228✔
32
33	// Initialize Random Ray Simulation Object
34	RandomRaySimulation sim;	323✔
35
36	// Begin main simulation timer
37	simulation::time_total.start();	323✔
38
39	// Execute random ray simulation
40	sim.simulate();	323✔
41
42	// End main simulation timer
43	openmc::simulation::time_total.stop();	323✔
44
45	// Finalize OpenMC
46	openmc_simulation_finalize();	323✔
47
48	// Reduce variables across MPI ranks
49	sim.reduce_simulation_statistics();	323✔
50
51	// Output all simulation results
52	sim.output_simulation_results();	323✔
53	}	323✔
54
55	// Enforces restrictions on inputs in random ray mode. While there are
56	// many features that don't make sense in random ray mode, and are therefore
57	// unsupported, we limit our testing/enforcement operations only to inputs
58	// that may cause erroneous/misleading output or crashes from the solver.
59	void validate_random_ray_inputs()	228✔
60	{
61	// Validate tallies
62	///////////////////////////////////////////////////////////////////
63	for (auto& tally : model::tallies) {	756✔
64
65	// Validate score types
66	for (auto score_bin : tally->scores_) {	1,176✔
67	switch (score_bin) {	648✔
68	case SCORE_FLUX:	648✔
69	case SCORE_TOTAL:
70	case SCORE_FISSION:
71	case SCORE_NU_FISSION:
72	case SCORE_EVENTS:
73	break;	648✔
74	default:	×
75	fatal_error(	×
76	"Invalid score specified. Only flux, total, fission, nu-fission, and "
77	"event scores are supported in random ray mode.");
78	}
79	}
80
81	// Validate filter types
82	for (auto f : tally->filters()) {	1,116✔
83	auto& filter = *model::tally_filters[f];	588✔
84
85	switch (filter.type()) {	588✔
86	case FilterType::CELL:	588✔
87	case FilterType::CELL_INSTANCE:
88	case FilterType::DISTRIBCELL:
89	case FilterType::ENERGY:
90	case FilterType::MATERIAL:
91	case FilterType::MESH:
92	case FilterType::UNIVERSE:
93	break;	588✔
94	default:	×
95	fatal_error("Invalid filter specified. Only cell, cell_instance, "	×
96	"distribcell, energy, material, mesh, and universe filters "
97	"are supported in random ray mode.");
98	}
99	}
100	}
101
102	// Validate MGXS data
103	///////////////////////////////////////////////////////////////////
104	for (auto& material : data::mg.macro_xs_) {	828✔
105	if (!material.is_isotropic) {	600✔
106	fatal_error("Anisotropic MGXS detected. Only isotropic XS data sets "	×
107	"supported in random ray mode.");
108	}
109	if (material.get_xsdata().size() > 1) {	600✔
110	fatal_error("Non-isothermal MGXS detected. Only isothermal XS data sets "	×
111	"supported in random ray mode.");
112	}
113	}
114
115	// Validate ray source
116	///////////////////////////////////////////////////////////////////
117
118	// Check for independent source
119	IndependentSource* is =
120	dynamic_cast<IndependentSource*>(RandomRay::ray_source_.get());	228✔
121	if (!is) {	228✔
122	fatal_error("Invalid ray source definition. Ray source must provided and "	×
123	"be of type IndependentSource.");
124	}
125
126	// Check for box source
127	SpatialDistribution* space_dist = is->space();	228✔
128	SpatialBox* sb = dynamic_cast<SpatialBox*>(space_dist);	228✔
129	if (!sb) {	228✔
130	fatal_error(	×
131	"Invalid ray source definition -- only box sources are allowed.");
132	}
133
134	// Check that box source is not restricted to fissionable areas
135	if (sb->only_fissionable()) {	228✔
136	fatal_error(	×
137	"Invalid ray source definition -- fissionable spatial distribution "
138	"not allowed.");
139	}
140
141	// Check for isotropic source
142	UnitSphereDistribution* angle_dist = is->angle();	228✔
143	Isotropic* id = dynamic_cast<Isotropic*>(angle_dist);	228✔
144	if (!id) {	228✔
145	fatal_error("Invalid ray source definition -- only isotropic sources are "	×
146	"allowed.");
147	}
148
149	// Validate external sources
150	///////////////////////////////////////////////////////////////////
151	if (settings::run_mode == RunMode::FIXED_SOURCE) {	228✔
152	if (model::external_sources.size() < 1) {	180✔
153	fatal_error("Must provide a particle source (in addition to ray source) "	×
154	"in fixed source random ray mode.");
155	}
156
157	for (int i = 0; i < model::external_sources.size(); i++) {	360✔
158	Source* s = model::external_sources[i].get();	180✔
159
160	// Check for independent source
161	IndependentSource* is = dynamic_cast<IndependentSource*>(s);	180✔
162
163	if (!is) {	180✔
164	fatal_error(	×
165	"Only IndependentSource external source types are allowed in "
166	"random ray mode");
167	}
168
169	// Check for isotropic source
170	UnitSphereDistribution* angle_dist = is->angle();	180✔
171	Isotropic* id = dynamic_cast<Isotropic*>(angle_dist);	180✔
172	if (!id) {	180✔
173	fatal_error(	×
174	"Invalid source definition -- only isotropic external sources are "
175	"allowed in random ray mode.");
176	}
177
178	// Validate that a domain ID was specified
179	if (is->domain_ids().size() == 0) {	180✔
180	fatal_error("Fixed sources must be specified by domain "	×
181	"id (cell, material, or universe) in random ray mode.");
182	}
183
184	// Check that a discrete energy distribution was used
185	Distribution* d = is->energy();	180✔
186	Discrete* dd = dynamic_cast<Discrete*>(d);	180✔
187	if (!dd) {	180✔
188	fatal_error(	×
189	"Only discrete (multigroup) energy distributions are allowed for "
190	"external sources in random ray mode.");
191	}
192	}
193	}
194
195	// Validate plotting files
196	///////////////////////////////////////////////////////////////////
197	for (int p = 0; p < model::plots.size(); p++) {	228✔
198
199	// Get handle to OpenMC plot object
200	Plot* openmc_plot = dynamic_cast<Plot*>(model::plots[p].get());	×
201
202	// Random ray plots only support voxel plots
203	if (!openmc_plot) {	×
204	warning(fmt::format(	×
205	"Plot {} will not be used for end of simulation data plotting -- only "
206	"voxel plotting is allowed in random ray mode.",
207	p));
208	continue;	×
209	} else if (openmc_plot->type_ != Plot::PlotType::voxel) {	×
210	warning(fmt::format(	×
211	"Plot {} will not be used for end of simulation data plotting -- only "
212	"voxel plotting is allowed in random ray mode.",
213	p));
214	continue;	×
215	}
216	}
217
218	// Warn about slow MPI domain replication, if detected
219	///////////////////////////////////////////////////////////////////
220	#ifdef OPENMC_MPI
221	if (mpi::n_procs > 1) {	95✔
222	warning(	95✔
223	"Domain replication in random ray is supported, but suffers from poor "
224	"scaling of source all-reduce operations. Performance may severely "
225	"degrade beyond just a few MPI ranks. Domain decomposition may be "
226	"implemented in the future to provide efficient scaling.");
227	}
228	#endif
229	}	228✔
230
231	//==============================================================================
232	// RandomRaySimulation implementation
233	//==============================================================================
234
235	RandomRaySimulation::RandomRaySimulation()	323✔
236	: negroups_(data::mg.num_energy_groups_)	323✔
237	{
238	// There are no source sites in random ray mode, so be sure to disable to
239	// ensure we don't attempt to write source sites to statepoint
240	settings::source_write = false;	323✔
241
242	// Random ray mode does not have an inner loop over generations within a
243	// batch, so set the current gen to 1
244	simulation::current_gen = 1;	323✔
245
246	switch (RandomRay::source_shape_) {	323✔
247	case RandomRaySourceShape::FLAT:	187✔
248	domain_ = make_unique<FlatSourceDomain>();	187✔
249	break;	187✔
250	case RandomRaySourceShape::LINEAR:	136✔
251	case RandomRaySourceShape::LINEAR_XY:
252	domain_ = make_unique<LinearSourceDomain>();	136✔
253	break;	136✔
254	default:	×
255	fatal_error("Unknown random ray source shape");	×
256	}
257	}	323✔
258
259	void RandomRaySimulation::simulate()	323✔
260	{
261	if (settings::run_mode == RunMode::FIXED_SOURCE) {	323✔
262	// Transfer external source user inputs onto random ray source regions
263	domain_->convert_external_sources();	255✔
264	domain_->count_external_source_regions();	255✔
265	}
266
267	// Random ray power iteration loop
268	while (simulation::current_batch < settings::n_batches) {	11,033✔
269
270	// Initialize the current batch
271	initialize_batch();	10,710✔
272	initialize_generation();	10,710✔
273
274	// Reset total starting particle weight used for normalizing tallies
275	simulation::total_weight = 1.0;	10,710✔
276
277	// Update source term (scattering + fission)
278	domain_->update_neutron_source(k_eff_);	10,710✔
279
280	// Reset scalar fluxes, iteration volume tallies, and region hit flags to
281	// zero
282	domain_->batch_reset();	10,710✔
283
284	// Start timer for transport
285	simulation::time_transport.start();	10,710✔
286
287	// Transport sweep over all random rays for the iteration
288	#pragma omp parallel for schedule(dynamic) \	5,670✔
289	reduction(+ : total_geometric_intersections_)	5,670✔
290	for (int i = 0; i < simulation::work_per_rank; i++) {	261,240✔
291	RandomRay ray(i, domain_.get());	256,200✔
292	total_geometric_intersections_ +=	256,200✔
293	ray.transport_history_based_single_ray();	256,200✔
294	}	256,200✔
295
296	simulation::time_transport.stop();	10,710✔
297
298	// If using multiple MPI ranks, perform all reduce on all transport results
299	domain_->all_reduce_replicated_source_regions();	10,710✔
300
301	// Normalize scalar flux and update volumes
302	domain_->normalize_scalar_flux_and_volumes(	10,710✔
303	settings::n_particles * RandomRay::distance_active_);
304
305	// Add source to scalar flux, compute number of FSR hits
306	int64_t n_hits = domain_->add_source_to_scalar_flux();	10,710✔
307
308	if (settings::run_mode == RunMode::EIGENVALUE) {	10,710✔
309	// Compute random ray k-eff
310	k_eff_ = domain_->compute_k_eff(k_eff_);	1,700✔
311
312	// Store random ray k-eff into OpenMC's native k-eff variable
313	global_tally_tracklength = k_eff_;	1,700✔
314	}
315
316	// Execute all tallying tasks, if this is an active batch
317	if (simulation::current_batch > settings::n_inactive && mpi::master) {	10,710✔
318
319	// Generate mapping between source regions and tallies
320	if (!domain_->mapped_all_tallies_) {	2,880✔
321	domain_->convert_source_regions_to_tallies();	228✔
322	}
323
324	// Use above mapping to contribute FSR flux data to appropriate tallies
325	domain_->random_ray_tally();	2,880✔
326
327	// Add this iteration's scalar flux estimate to final accumulated estimate
328	domain_->accumulate_iteration_flux();	2,880✔
329	}
330
331	// Set phi_old = phi_new
332	domain_->flux_swap();	10,710✔
333
334	// Check for any obvious insabilities/nans/infs
335	instability_check(n_hits, k_eff_, avg_miss_rate_);	10,710✔
336
337	// Finalize the current batch
338	finalize_generation();	10,710✔
339	finalize_batch();	10,710✔
340	} // End random ray power iteration loop
341	}	323✔
342
343	void RandomRaySimulation::reduce_simulation_statistics()	323✔
344	{
345	// Reduce number of intersections
346	#ifdef OPENMC_MPI
347	if (mpi::n_procs > 1) {	190✔
348	uint64_t total_geometric_intersections_reduced = 0;	190✔
349	MPI_Reduce(&total_geometric_intersections_,	190✔
350	&total_geometric_intersections_reduced, 1, MPI_UNSIGNED_LONG, MPI_SUM, 0,
351	mpi::intracomm);
352	total_geometric_intersections_ = total_geometric_intersections_reduced;	190✔
353	}
354	#endif
355	}	323✔
356
357	void RandomRaySimulation::output_simulation_results() const	323✔
358	{
359	// Print random ray results
360	if (mpi::master) {	323✔
361	print_results_random_ray(total_geometric_intersections_,	456✔
362	avg_miss_rate_ / settings::n_batches, negroups_,	228✔
363	domain_->n_source_regions_, domain_->n_external_source_regions_);	228✔
364	if (model::plots.size() > 0) {	228✔
365	domain_->output_to_vtk();	×
366	}
367	}
368	}	323✔
369
370	// Apply a few sanity checks to catch obvious cases of numerical instability.
371	// Instability typically only occurs if ray density is extremely low.
372	void RandomRaySimulation::instability_check(	10,710✔
373	int64_t n_hits, double k_eff, double& avg_miss_rate) const
374	{
375	double percent_missed = ((domain_->n_source_regions_ - n_hits) /	10,710✔
376	static_cast<double>(domain_->n_source_regions_)) *	10,710✔
377	100.0;	10,710✔
378	avg_miss_rate += percent_missed;	10,710✔
379
380	if (mpi::master) {	10,710✔
381	if (percent_missed > 10.0) {	7,560✔
382	warning(fmt::format(	×
383	"Very high FSR miss rate detected ({:.3f}%). Instability may occur. "
384	"Increase ray density by adding more rays and/or active distance.",
385	percent_missed));
386	} else if (percent_missed > 1.0) {	7,560✔
387	warning(	×
388	fmt::format("Elevated FSR miss rate detected ({:.3f}%). Increasing "	×
389	"ray density by adding more rays and/or active "
390	"distance may improve simulation efficiency.",
391	percent_missed));
392	}
393
394	if (k_eff > 10.0 \|\| k_eff < 0.01 \|\| !(std::isfinite(k_eff))) {	7,560✔
395	fatal_error("Instability detected");	×
396	}
397	}
398	}	10,710✔
399
400	// Print random ray simulation results
401	void RandomRaySimulation::print_results_random_ray(	228✔
402	uint64_t total_geometric_intersections, double avg_miss_rate, int negroups,
403	int64_t n_source_regions, int64_t n_external_source_regions) const
404	{
405	using namespace simulation;
406
407	if (settings::verbosity >= 6) {	228✔
408	double total_integrations = total_geometric_intersections * negroups;	228✔
409	double time_per_integration =
410	simulation::time_transport.elapsed() / total_integrations;	228✔
411	double misc_time = time_total.elapsed() - time_update_src.elapsed() -	228✔
412	time_transport.elapsed() - time_tallies.elapsed() -	228✔
413	time_bank_sendrecv.elapsed();	228✔
414
415	header("Simulation Statistics", 4);	228✔
416	fmt::print(	228✔
417	" Total Iterations = {}\n", settings::n_batches);
418	fmt::print(" Flat Source Regions (FSRs) = {}\n", n_source_regions);	228✔
419	fmt::print(	228✔
420	" FSRs Containing External Sources = {}\n", n_external_source_regions);
421	fmt::print(" Total Geometric Intersections = {:.4e}\n",	×
422	static_cast<double>(total_geometric_intersections));	228✔
423	fmt::print(" Avg per Iteration = {:.4e}\n",	×
424	static_cast<double>(total_geometric_intersections) / settings::n_batches);	228✔
425	fmt::print(" Avg per Iteration per FSR = {:.2f}\n",	×
426	static_cast<double>(total_geometric_intersections) /	228✔
427	static_cast<double>(settings::n_batches) / n_source_regions);	228✔
428	fmt::print(" Avg FSR Miss Rate per Iteration = {:.4f}%\n", avg_miss_rate);	228✔
429	fmt::print(" Energy Groups = {}\n", negroups);	228✔
430	fmt::print(	228✔
431	" Total Integrations = {:.4e}\n", total_integrations);
432	fmt::print(" Avg per Iteration = {:.4e}\n",	×
433	total_integrations / settings::n_batches);	228✔
434
435	std::string estimator;	228✔
436	switch (domain_->volume_estimator_) {	228✔
437	case RandomRayVolumeEstimator::SIMULATION_AVERAGED:	24✔
438	estimator = "Simulation Averaged";	24✔
439	break;	24✔
440	case RandomRayVolumeEstimator::NAIVE:	24✔
441	estimator = "Naive";	24✔
442	break;	24✔
443	case RandomRayVolumeEstimator::HYBRID:	180✔
444	estimator = "Hybrid";	180✔
445	break;	180✔
446	default:	×
447	fatal_error("Invalid volume estimator type");	×
448	}
449	fmt::print(" Volume Estimator Type = {}\n", estimator);	228✔
450
451	header("Timing Statistics", 4);	228✔
452	show_time("Total time for initialization", time_initialize.elapsed());	228✔
453	show_time("Reading cross sections", time_read_xs.elapsed(), 1);	228✔
454	show_time("Total simulation time", time_total.elapsed());	228✔
455	show_time("Transport sweep only", time_transport.elapsed(), 1);	228✔
456	show_time("Source update only", time_update_src.elapsed(), 1);	228✔
457	show_time("Tally conversion only", time_tallies.elapsed(), 1);	228✔
458	show_time("MPI source reductions only", time_bank_sendrecv.elapsed(), 1);	228✔
459	show_time("Other iteration routines", misc_time, 1);	228✔
460	if (settings::run_mode == RunMode::EIGENVALUE) {	228✔
461	show_time("Time in inactive batches", time_inactive.elapsed());	48✔
462	}
463	show_time("Time in active batches", time_active.elapsed());	228✔
464	show_time("Time writing statepoints", time_statepoint.elapsed());	228✔
465	show_time("Total time for finalization", time_finalize.elapsed());	228✔
466	show_time("Time per integration", time_per_integration);	228✔
467	}	228✔
468
469	if (settings::verbosity >= 4 && settings::run_mode == RunMode::EIGENVALUE) {	228✔
470	header("Results", 4);	48✔
471	fmt::print(" k-effective = {:.5f} +/- {:.5f}\n",	96✔
472	simulation::keff, simulation::keff_std);
473	}
474	}	228✔
475
476	} // namespace openmc

openmc-dev / openmc / 10586562087

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

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