From b0177d9f7dcd243a4f871b318dadbebeed707e4f Mon Sep 17 00:00:00 2001 From: spasmann Date: Sun, 13 Jul 2025 18:14:25 -0700 Subject: [PATCH 1/2] Add random ray info to statepoint file --- include/openmc/random_ray/random_ray.h | 11 +- .../openmc/random_ray/random_ray_simulation.h | 12 +- openmc/statepoint.py | 105 ++++++++++++++++++ src/random_ray/random_ray.cpp | 9 +- src/random_ray/random_ray_simulation.cpp | 83 +++++++++++--- src/state_point.cpp | 11 ++ 6 files changed, 210 insertions(+), 21 deletions(-) diff --git a/include/openmc/random_ray/random_ray.h b/include/openmc/random_ray/random_ray.h index abf2a26881f..5e1545cbdf5 100644 --- a/include/openmc/random_ray/random_ray.h +++ b/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 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 mesh_bins_; vector mesh_fractional_lengths_; - int negroups_; FlatSourceDomain* domain_ {nullptr}; // pointer to domain that has flat source // data needed for ray transport double distance_travelled_ {0}; diff --git a/include/openmc/random_ray/random_ray_simulation.h b/include/openmc/random_ray/random_ray_simulation.h index b94e7401b3e..63f79e85fbb 100644 --- a/include/openmc/random_ray/random_ray_simulation.h +++ b/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 diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 715becf4882..6f3beacbec2 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -31,6 +31,10 @@ class StatePoint: Attributes ---------- + adjoint_mode : bool + Indicate whether random ray solve was in adjoint mode + avg_miss_rate : float + The random ray average source region miss rate per iteration cmfd_on : bool Indicate whether CMFD is active cmfd_balance : numpy.ndarray @@ -52,6 +56,10 @@ class StatePoint: Number of batches simulated date_and_time : datetime.datetime Date and time at which statepoint was written + distance_active : float + Active distance for each ray in Random ray + distance_inactive : float + Inactive distance for each ray in Random ray entropy : numpy.ndarray Shannon entropy of fission source at each batch filters : dict @@ -82,12 +90,20 @@ class StatePoint: Dictionary whose keys are mesh IDs and whose values are MeshBase objects n_batches : int Number of batches + n_external_source_regions : int + Number of external source regions in random ray simulation + n_geometric_intersections : int + Total number of geometric intersections in random ray simulation n_inactive : int Number of inactive batches + n_integrations : int + Total number of integrations in random ray simulation n_particles : int Number of particles per generation n_realizations : int Number of tally realizations + n_source_regions : int + Number of source regions in random ray simulation path : str Working directory for simulation photon_transport : bool @@ -97,8 +113,14 @@ class StatePoint: runtime : dict Dictionary whose keys are strings describing various runtime metrics and whose values are time values in seconds. + sample_method : str + Random ray sample method, e.g. 'prng' or 'halton' seed : int Pseudorandom number generator seed + solver_type : str + Transport method, e.g. 'monte carlo' or 'random ray' + source_shape : str + Random ray source region approximation, e.g. 'flat' or 'linear' stride : int Number of random numbers allocated for each particle history source : numpy.ndarray of compound datatype @@ -120,6 +142,8 @@ class StatePoint: TallyDerivative objects version: tuple of Integral Version of OpenMC + volume_estimator : str + Random ray source region volume estimatation method, e.g. 'naive' summary : None or openmc.Summary A summary object if the statepoint has been linked with a summary file @@ -304,6 +328,83 @@ def k_abs_tra(self): else: return None + @property + def source_shape(self): + if self.solver_type == 'random ray': + return self._f['source_shape'][()].decode() + else: + return None + + @property + def sample_method(self): + if self.solver_type == 'random ray': + return self._f['sample_method'][()].decode() + else: + return None + + @property + def volume_estimator(self): + if self.solver_type == 'random ray': + return self._f['volume_estimator'][()].decode() + else: + return None + + @property + def distance_active(self): + if self.solver_type == 'random ray': + return self._f['distance_active'][()] + else: + return None + + @property + def distance_inactive(self): + if self.solver_type == 'random ray': + return self._f['distance_inactive'][()] + else: + return None + + @property + def adjoint_mode(self): + if self.solver_type == 'random ray': + return True if self._f['adjoint_mode'][()] else False + else: + return None + + @property + def avg_miss_rate(self): + if self.solver_type == 'random ray': + return self._f['avg_miss_rate'][()] + else: + return None + + @property + def n_source_regions(self): + if self.solver_type == 'random ray': + return self._f['n_source_regions'][()] + else: + return None + + @property + def n_external_source_regions(self): + if self.solver_type == 'random ray': + return self._f['n_external_source_regions'][()] + else: + return None + + @property + def n_geometric_intersections(self): + if self.solver_type == 'random ray': + return self._f['n_geometric_intersections'][()] + else: + return None + + @property + def n_integrations(self): + if self.solver_type == 'random ray': + return self._f['n_integrations'][()] + else: + return None + @property def meshes(self): if not self._meshes_read: @@ -345,6 +446,10 @@ 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 run_mode(self): return self._f['run_mode'][()].decode() diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index 36eab062846..64e68c13008 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -235,16 +235,21 @@ vector 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 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_); diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 388a778b840..1787bad4340 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/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(total_geometric_intersections)); + static_cast(RandomRay::total_geometric_intersections)); fmt::print(" Avg per Iteration = {:.4e}\n", - static_cast(total_geometric_intersections) / settings::n_batches); + static_cast(RandomRay::total_geometric_intersections) / + settings::n_batches); fmt::print(" Avg per Iteration per SR = {:.2f}\n", - static_cast(total_geometric_intersections) / - static_cast(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(RandomRay::total_geometric_intersections) / + static_cast(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", diff --git a/src/state_point.cpp b/src/state_point.cpp index 12dfeb82ac6..f6dc6dea426 100644 --- a/src/state_point.cpp +++ b/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 From 9fec33dd7a2ccf3eb9cfa09e7e0d99da59391f85 Mon Sep 17 00:00:00 2001 From: spasmann Date: Mon, 14 Jul 2025 09:01:26 -0700 Subject: [PATCH 2/2] Return random ray variables in dictionary --- openmc/settings.py | 2 +- openmc/statepoint.py | 150 +++++++++++++++---------------------------- 2 files changed, 52 insertions(+), 100 deletions(-) diff --git a/openmc/settings.py b/openmc/settings.py index 36eaec6cb45..2f75906f729 100644 --- a/openmc/settings.py +++ b/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: diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 6f3beacbec2..408a1b956de 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -31,10 +31,6 @@ class StatePoint: Attributes ---------- - adjoint_mode : bool - Indicate whether random ray solve was in adjoint mode - avg_miss_rate : float - The random ray average source region miss rate per iteration cmfd_on : bool Indicate whether CMFD is active cmfd_balance : numpy.ndarray @@ -56,10 +52,6 @@ class StatePoint: Number of batches simulated date_and_time : datetime.datetime Date and time at which statepoint was written - distance_active : float - Active distance for each ray in Random ray - distance_inactive : float - Inactive distance for each ray in Random ray entropy : numpy.ndarray Shannon entropy of fission source at each batch filters : dict @@ -90,37 +82,55 @@ class StatePoint: Dictionary whose keys are mesh IDs and whose values are MeshBase objects n_batches : int Number of batches - n_external_source_regions : int - Number of external source regions in random ray simulation - n_geometric_intersections : int - Total number of geometric intersections in random ray simulation n_inactive : int Number of inactive batches - n_integrations : int - Total number of integrations in random ray simulation n_particles : int Number of particles per generation n_realizations : int Number of tally realizations - n_source_regions : int - Number of source regions in random ray simulation path : str 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. - sample_method : str - Random ray sample method, e.g. 'prng' or 'halton' seed : int Pseudorandom number generator seed solver_type : str Transport method, e.g. 'monte carlo' or 'random ray' - source_shape : str - Random ray source region approximation, e.g. 'flat' or 'linear' stride : int Number of random numbers allocated for each particle history source : numpy.ndarray of compound datatype @@ -142,8 +152,6 @@ class StatePoint: TallyDerivative objects version: tuple of Integral Version of OpenMC - volume_estimator : str - Random ray source region volume estimatation method, e.g. 'naive' summary : None or openmc.Summary A summary object if the statepoint has been linked with a summary file @@ -156,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) @@ -168,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: @@ -328,83 +338,6 @@ def k_abs_tra(self): else: return None - @property - def source_shape(self): - if self.solver_type == 'random ray': - return self._f['source_shape'][()].decode() - else: - return None - - @property - def sample_method(self): - if self.solver_type == 'random ray': - return self._f['sample_method'][()].decode() - else: - return None - - @property - def volume_estimator(self): - if self.solver_type == 'random ray': - return self._f['volume_estimator'][()].decode() - else: - return None - - @property - def distance_active(self): - if self.solver_type == 'random ray': - return self._f['distance_active'][()] - else: - return None - - @property - def distance_inactive(self): - if self.solver_type == 'random ray': - return self._f['distance_inactive'][()] - else: - return None - - @property - def adjoint_mode(self): - if self.solver_type == 'random ray': - return True if self._f['adjoint_mode'][()] else False - else: - return None - - @property - def avg_miss_rate(self): - if self.solver_type == 'random ray': - return self._f['avg_miss_rate'][()] - else: - return None - - @property - def n_source_regions(self): - if self.solver_type == 'random ray': - return self._f['n_source_regions'][()] - else: - return None - - @property - def n_external_source_regions(self): - if self.solver_type == 'random ray': - return self._f['n_external_source_regions'][()] - else: - return None - - @property - def n_geometric_intersections(self): - if self.solver_type == 'random ray': - return self._f['n_geometric_intersections'][()] - else: - return None - - @property - def n_integrations(self): - if self.solver_type == 'random ray': - return self._f['n_integrations'][()] - else: - return None - @property def meshes(self): if not self._meshes_read: @@ -450,6 +383,25 @@ def photon_transport(self): 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()