Add random ray info to statepoint file

include/openmc/random_ray/random_ray.h

@@ -46,10 +46,20 @@ class RandomRay : public Particle {
   // Static data members
   static double distance_inactive_;          // Inactive (dead zone) ray length
   static double distance_active_;            // Active ray length
+  static double avg_miss_rate;               // Average ray miss rate per
+                                             // iteration for reporting
   static unique_ptr<Source> ray_source_;     // Starting source for ray sampling
   static RandomRaySourceShape source_shape_; // Flag for linear source
   static bool mesh_subdivision_enabled_;     // Flag for mesh subdivision
   static RandomRaySampleMethod sample_method_; // Flag for sampling method
+  static int64_t n_source_regions;             // Total number of source regions
+  static int64_t
+    n_external_source_regions; // Total number of source regions with
+                               // non-zero external source terms
+  static uint64_t total_geometric_intersections; // Tracks the total number of
+                                                 // geometric intersections by
+                                                 // all rays for reporting
+  static int negroups_;                          // Number of energy groups
 
   //----------------------------------------------------------------------------
   // Public data members
@@ -65,7 +75,6 @@ class RandomRay : public Particle {
   vector<int> mesh_bins_;
   vector<double> mesh_fractional_lengths_;
 
-  int negroups_;
   FlatSourceDomain* domain_ {nullptr}; // pointer to domain that has flat source
                                        // data needed for ray transport
   double distance_travelled_ {0};

include/openmc/random_ray/random_ray_simulation.h

@@ -30,9 +30,11 @@ class RandomRaySimulation {
   void output_simulation_results() const;
   void instability_check(
     int64_t n_hits, double k_eff, double& avg_miss_rate) const;
-  void 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;
+  void print_results_random_ray() const;
+  int64_t total_geometric_intersections() const
+  {
+    return total_geometric_intersections_;
+  }
 
   //----------------------------------------------------------------------------
   // Accessors
@@ -68,6 +70,10 @@ void openmc_run_random_ray();
 void validate_random_ray_inputs();
 void openmc_reset_random_ray();
 
+//! Write data related to randaom ray to statepoint
+//! \param[in] group HDF5 group
+void write_random_ray_hdf5(hid_t group);
+
 } // namespace openmc
 
 #endif // OPENMC_RANDOM_RAY_SIMULATION_H

openmc/settings.py

@@ -159,7 +159,7 @@ class Settings:
         Options for configuring the random ray solver. Acceptable keys are:
 
         :distance_inactive:
-            Indicates the total active distance in [cm] a ray should travel
+            Indicates the total inactive distance in [cm] a ray should travel
         :distance_active:
             Indicates the total active distance in [cm] a ray should travel
         :ray_source:

openmc/statepoint.py

@@ -92,13 +92,45 @@ class StatePoint:
         Working directory for simulation
     photon_transport : bool
         Indicate whether photon transport is active
+    random_ray : dict
+        Dictionary for random ray solver parameters and results.
+        Acceptable keys are:
+
+        :adjoint_mode:
+            Indicate whether random ray solve was in adjoint mode
+        :avg_miss_rate:
+            The random ray average source region miss rate per iteration
+            expressed as a percent
+        :distance_active:
+            Indicates the total active distance in [cm] for each ray
+        :distance_inactive:
+            Indicates the total inactive distance in [cm] for each ray
+        :sample_method:
+            Sampling method for the ray starting location and direction of
+            travel, e.g. `prng` or 'halton`
+        :source_shape:
+            Assumed shape of the source distribution within each source region,
+            e.g. 'flat' (default), 'linear', or 'linear_xy'
+        :n_external_source_regions:
+            Number of external source regions in random ray simulation
+        :n_geometric_intersections:
+            Total number of geometric intersections in random ray simulation
+        :n_integrations:
+            Total number of integrations in random ray simulation
+        :n_source_regions:
+            Number of source regions in random ray simulation
+        :volume_estimator:
+            Choice of volume estimator for the random ray solver, e.g.
+            'naive', 'simulation_averaged', or 'hybrid'
     run_mode : str
         Simulation run mode, e.g. 'eigenvalue'
     runtime : dict
         Dictionary whose keys are strings describing various runtime metrics
         and whose values are time values in seconds.
     seed : int
         Pseudorandom number generator seed
+    solver_type : str
+        Transport method, e.g. 'monte carlo' or 'random ray'
     stride : int
         Number of random numbers allocated for each particle history
     source : numpy.ndarray of compound datatype
@@ -132,6 +164,7 @@ def __init__(self, filepath, autolink=True):
         self._filters = {}
         self._tallies = {}
         self._derivs = {}
+        self._random_ray = {}
 
         # Check filetype and version
         cv.check_filetype_version(self._f, 'statepoint', _VERSION_STATEPOINT)
@@ -144,6 +177,7 @@ def __init__(self, filepath, autolink=True):
         self._global_tallies = None
         self._sparse = False
         self._derivs_read = False
+        self._random_ray_read = False
 
         # Automatically link in a summary file if one exists
         if autolink:
@@ -345,6 +379,29 @@ def path(self):
     def photon_transport(self):
         return self._f.attrs['photon_transport'] > 0
 
+    @property
+    def solver_type(self):
+        return self._f['solver_type'][()].decode()
+
+    @property
+    def random_ray(self):
+        if not self._random_ray_read:
+            self._random_ray['adjoint_mode'] = True if self._f['adjoint_mode'][()] else False
+            self._random_ray['avg_miss_rate'] = self._f['avg_miss_rate'][()]
+            self._random_ray['distance_active'] = self._f['distance_active'][()]
+            self._random_ray['distance_inactive'] = self._f['distance_inactive'][()]
+            self._random_ray['sample_method'] = self._f['sample_method'][()].decode()
+            self._random_ray['source_shape'] = self._f['source_shape'][()].decode()
+            self._random_ray['n_external_source_regions'] = self._f['n_external_source_regions'][()]
+            self._random_ray['n_geometric_intersections'] = self._f['n_geometric_intersections'][()]
+            self._random_ray['n_integrations'] = self._f['n_integrations'][()]
+            self._random_ray['n_source_regions'] = self._f['n_source_regions'][()]
+            self._random_ray['volume_estimator'] = self._f['volume_estimator'][()].decode()
+
+            self._random_ray_read = True
+
+        return self._random_ray
+
     @property
     def run_mode(self):
         return self._f['run_mode'][()].decode()

src/random_ray/random_ray.cpp

@@ -235,16 +235,21 @@ vector<double> rhalton(int dim, uint64_t* seed, int64_t skip = 0)
 // Static Variable Declarations
 double RandomRay::distance_inactive_;
 double RandomRay::distance_active_;
+double RandomRay::avg_miss_rate;
 unique_ptr<Source> RandomRay::ray_source_;
 RandomRaySourceShape RandomRay::source_shape_ {RandomRaySourceShape::FLAT};
 bool RandomRay::mesh_subdivision_enabled_ {false};
 RandomRaySampleMethod RandomRay::sample_method_ {RandomRaySampleMethod::PRNG};
+int64_t RandomRay::n_source_regions;
+int64_t RandomRay::n_external_source_regions;
+uint64_t RandomRay::total_geometric_intersections;
+int RandomRay::negroups_;
 
 RandomRay::RandomRay()
   : angular_flux_(data::mg.num_energy_groups_),
-    delta_psi_(data::mg.num_energy_groups_),
-    negroups_(data::mg.num_energy_groups_)
+    delta_psi_(data::mg.num_energy_groups_)
 {
+  negroups_ = data::mg.num_energy_groups_;
   if (source_shape_ == RandomRaySourceShape::LINEAR ||
       source_shape_ == RandomRaySourceShape::LINEAR_XY) {
     delta_moments_.resize(negroups_);

src/random_ray/random_ray_simulation.cpp

@@ -370,6 +370,52 @@ void openmc_reset_random_ray()
   RandomRay::sample_method_ = RandomRaySampleMethod::PRNG;
 }
 
+void write_random_ray_hdf5(hid_t group)
+{
+  switch (RandomRay::source_shape_) {
+  case RandomRaySourceShape::FLAT:
+    write_dataset(group, "source_shape", "flat");
+    break;
+  case RandomRaySourceShape::LINEAR:
+    write_dataset(group, "source_shape", "linear");
+    break;
+  case RandomRaySourceShape::LINEAR_XY:
+    write_dataset(group, "source_shape", "linear_xy");
+    break;
+  default:
+    break;
+  }
+  std::string sample_method =
+    (RandomRay::sample_method_ == RandomRaySampleMethod::PRNG) ? "prng"
+                                                               : "halton";
+  write_dataset(group, "sample_method", sample_method);
+  switch (FlatSourceDomain::volume_estimator_) {
+  case RandomRayVolumeEstimator::SIMULATION_AVERAGED:
+    write_dataset(group, "volume_estimator", "simulation_averaged");
+    break;
+  case RandomRayVolumeEstimator::NAIVE:
+    write_dataset(group, "volume_estimator", "naive");
+    break;
+  case RandomRayVolumeEstimator::HYBRID:
+    write_dataset(group, "volume_estimator", "hybrid");
+    break;
+  default:
+    break;
+  }
+  write_dataset(group, "distance_active", RandomRay::distance_active_);
+  write_dataset(group, "distance_inactive", RandomRay::distance_inactive_);
+  write_dataset(group, "adjoint_mode", FlatSourceDomain::adjoint_);
+  write_dataset(group, "avg_miss_rate", RandomRay::avg_miss_rate);
+  write_dataset(group, "n_source_regions", RandomRay::n_source_regions);
+  write_dataset(
+    group, "n_external_source_regions", RandomRay::n_external_source_regions);
+  write_dataset(group, "n_geometric_intersections",
+    RandomRay::total_geometric_intersections);
+  int64_t n_integrations =
+    RandomRay::total_geometric_intersections * RandomRay::negroups_;
+  write_dataset(group, "n_integrations", n_integrations);
+}
+
 //==============================================================================
 // RandomRaySimulation implementation
 //==============================================================================
@@ -525,6 +571,12 @@ void RandomRaySimulation::simulate()
       instability_check(n_hits, k_eff_, avg_miss_rate_);
     } // End MPI master work
 
+    // Set global variables for reporting
+    RandomRay::avg_miss_rate = avg_miss_rate_ / settings::n_batches;
+    RandomRay::total_geometric_intersections = total_geometric_intersections_;
+    RandomRay::n_external_source_regions = domain_->n_external_source_regions_;
+    RandomRay::n_source_regions = domain_->n_source_regions();
+
     // Finalize the current batch
     finalize_generation();
     finalize_batch();
@@ -537,9 +589,7 @@ 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_);
+    print_results_random_ray();
     if (model::plots.size() > 0) {
       domain_->output_to_vtk();
     }
@@ -577,14 +627,13 @@ void RandomRaySimulation::instability_check(
 }
 
 // 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
+void RandomRaySimulation::print_results_random_ray() const
 {
   using namespace simulation;
 
   if (settings::verbosity >= 6) {
-    double total_integrations = total_geometric_intersections * negroups;
+    double total_integrations =
+      RandomRay::total_geometric_intersections * RandomRay::negroups_;
     double time_per_integration =
       simulation::time_transport.elapsed() / total_integrations;
     double misc_time = time_total.elapsed() - time_update_src.elapsed() -
@@ -600,18 +649,22 @@ void RandomRaySimulation::print_results_random_ray(
       RandomRay::distance_inactive_);
     fmt::print(" Active Distance                   = {} cm\n",
       RandomRay::distance_active_);
-    fmt::print(" Source Regions (SRs)              = {}\n", n_source_regions);
     fmt::print(
-      " SRs Containing External Sources   = {}\n", n_external_source_regions);
+      " Source Regions (SRs)              = {}\n", RandomRay::n_source_regions);
+    fmt::print(" SRs Containing External Sources   = {}\n",
+      RandomRay::n_external_source_regions);
     fmt::print(" Total Geometric Intersections     = {:.4e}\n",
-      static_cast<double>(total_geometric_intersections));
+      static_cast<double>(RandomRay::total_geometric_intersections));
     fmt::print("   Avg per Iteration               = {:.4e}\n",
-      static_cast<double>(total_geometric_intersections) / settings::n_batches);
+      static_cast<double>(RandomRay::total_geometric_intersections) /
+        settings::n_batches);
     fmt::print("   Avg per Iteration per SR        = {:.2f}\n",
-      static_cast<double>(total_geometric_intersections) /
-        static_cast<double>(settings::n_batches) / n_source_regions);
-    fmt::print(" Avg SR Miss Rate per Iteration    = {:.4f}%\n", avg_miss_rate);
-    fmt::print(" Energy Groups                     = {}\n", negroups);
+      static_cast<double>(RandomRay::total_geometric_intersections) /
+        static_cast<double>(settings::n_batches) / RandomRay::n_source_regions);
+    fmt::print(" Avg SR Miss Rate per Iteration    = {:.4f}%\n",
+      RandomRay::avg_miss_rate);
+    fmt::print(
+      " Energy Groups                     = {}\n", RandomRay::negroups_);
     fmt::print(
       " Total Integrations                = {:.4e}\n", total_integrations);
     fmt::print("   Avg per Iteration               = {:.4e}\n",

src/state_point.cpp

@@ -21,6 +21,7 @@
 #include "openmc/mgxs_interface.h"
 #include "openmc/nuclide.h"
 #include "openmc/output.h"
+#include "openmc/random_ray/random_ray_simulation.h"
 #include "openmc/settings.h"
 #include "openmc/simulation.h"
 #include "openmc/tallies/derivative.h"
@@ -119,6 +120,16 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source)
     if (settings::run_mode == RunMode::EIGENVALUE)
       write_eigenvalue_hdf5(file_id);
 
+    switch (settings::solver_type) {
+    case SolverType::RANDOM_RAY:
+      write_dataset(file_id, "solver_type", "random ray");
+      write_random_ray_hdf5(file_id);
+      break;
+    default:
+      write_dataset(file_id, "solver_type", "monte carlo");
+      break;
+    }
+
     hid_t tallies_group = create_group(file_id, "tallies");
 
     // Write meshes