From 5b8db460cbb4afee28f79d4c1f782d5d3ec3b4f1 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Thu, 18 Sep 2025 15:21:57 -0500 Subject: [PATCH 01/19] initial implementation complete. Now for debugging. --- .../openmc/random_ray/flat_source_domain.h | 28 +- .../openmc/random_ray/linear_source_domain.h | 2 +- .../openmc/random_ray/random_ray_simulation.h | 3 - src/random_ray/flat_source_domain.cpp | 249 +++++++++++------- src/random_ray/linear_source_domain.cpp | 4 +- src/random_ray/random_ray.cpp | 32 ++- src/random_ray/random_ray_simulation.cpp | 21 +- src/random_ray/source_region.cpp | 2 +- 8 files changed, 203 insertions(+), 138 deletions(-) diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index 78351fcc5f5..e61f37abe22 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -27,8 +27,8 @@ class FlatSourceDomain { //---------------------------------------------------------------------------- // Methods - virtual void update_neutron_source(double k_eff); - double compute_k_eff(double k_eff_old) const; + virtual void update_neutron_source(); + void compute_k_eff(); virtual void normalize_scalar_flux_and_volumes( double total_active_distance_per_iteration); @@ -86,6 +86,7 @@ class FlatSourceDomain { //---------------------------------------------------------------------------- // Public Data members + double k_eff_ {1.0}; // Eigenvalue bool mapped_all_tallies_ {false}; // If all source regions have been visited int64_t n_external_source_regions_ {0}; // Total number of source regions with @@ -134,8 +135,19 @@ class FlatSourceDomain { // Map that relates a SourceRegionKey to the external source index. This map // is used to check if there are any point sources within a subdivided source // region at the time it is discovered. - std::unordered_map - point_source_map_; + std::unordered_map, SourceRegionKey::HashFunctor> + external_point_source_map_; + + // Map that relates a base source region index to the external source index. + // This map is used to check if there are any volumetric sources within a + // subdivided source region at the time it is discovered. + std::unordered_map, SourceRegionKey::HashFunctor> + external_volumetric_source_map_; + + // Map that relates a base source region index to a mesh index. This map + // is used to check which subdivision mesh is present in a source region. + std::unordered_map + mesh_map_; // If transport corrected MGXS data is being used, there may be negative // in-group scattering cross sections that can result in instability in MOC @@ -146,13 +158,13 @@ class FlatSourceDomain { protected: //---------------------------------------------------------------------------- // Methods - void apply_external_source_to_source_region( - Discrete* discrete, double strength_factor, SourceRegionHandle& srh); + void FlatSourceDomain::apply_external_source_to_source_region( + int src_idx, SourceRegionHandle& srh); void apply_external_source_to_cell_instances(int32_t i_cell, - Discrete* discrete, double strength_factor, int target_material_id, + int src_idx, int target_material_id, const vector& instances); void apply_external_source_to_cell_and_children(int32_t i_cell, - Discrete* discrete, double strength_factor, int32_t target_material_id); + int src_idx, int32_t target_material_id); virtual void set_flux_to_flux_plus_source(int64_t sr, double volume, int g); void set_flux_to_source(int64_t sr, int g); virtual void set_flux_to_old_flux(int64_t sr, int g); diff --git a/include/openmc/random_ray/linear_source_domain.h b/include/openmc/random_ray/linear_source_domain.h index 67fdd99f880..7b5fd86dc52 100644 --- a/include/openmc/random_ray/linear_source_domain.h +++ b/include/openmc/random_ray/linear_source_domain.h @@ -20,7 +20,7 @@ class LinearSourceDomain : public FlatSourceDomain { public: //---------------------------------------------------------------------------- // Methods - void update_neutron_source(double k_eff) override; + void update_neutron_source() override; void normalize_scalar_flux_and_volumes( double total_active_distance_per_iteration) override; diff --git a/include/openmc/random_ray/random_ray_simulation.h b/include/openmc/random_ray/random_ray_simulation.h index b94e7401b3e..13bc68834ec 100644 --- a/include/openmc/random_ray/random_ray_simulation.h +++ b/include/openmc/random_ray/random_ray_simulation.h @@ -45,9 +45,6 @@ class RandomRaySimulation { // Contains all flat source region data unique_ptr domain_; - // Random ray eigenvalue - double k_eff_ {1.0}; - // Tracks the average FSR miss rate for analysis and reporting double avg_miss_rate_ {0.0}; diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 40923883088..69db68aef11 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -56,22 +56,6 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_) source_regions_.assign( base_source_regions, SourceRegion(negroups_, is_linear)); - // Initialize materials - int64_t source_region_id = 0; - for (int i = 0; i < model::cells.size(); i++) { - Cell& cell = *model::cells[i]; - if (cell.type_ == Fill::MATERIAL) { - for (int j = 0; j < cell.n_instances(); j++) { - source_regions_.material(source_region_id++) = cell.material(j); - } - } - } - - // Sanity check - if (source_region_id != base_source_regions) { - fatal_error("Unexpected number of source regions"); - } - // Initialize tally volumes if (volume_normalized_flux_tallies_) { tally_volumes_.resize(model::tallies.size()); @@ -120,11 +104,11 @@ void FlatSourceDomain::accumulate_iteration_flux() // Compute new estimate of scattering + fission sources in each source region // based on the flux estimate from the previous iteration. -void FlatSourceDomain::update_neutron_source(double k_eff) +void FlatSourceDomain::update_neutron_source() { simulation::time_update_src.start(); - double inverse_k_eff = 1.0 / k_eff; + double inverse_k_eff = 1.0 / k_eff_; // Reset all source regions to zero (important for void regions) #pragma omp parallel for @@ -320,7 +304,7 @@ int64_t FlatSourceDomain::add_source_to_scalar_flux() // Generates new estimate of k_eff based on the differences between this // iteration's estimate of the scalar flux and the last iteration's estimate. -double FlatSourceDomain::compute_k_eff(double k_eff_old) const +void FlatSourceDomain::compute_k_eff() { double fission_rate_old = 0; double fission_rate_new = 0; @@ -365,7 +349,7 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const p[sr] = sr_fission_source_new; } - double k_eff_new = k_eff_old * (fission_rate_new / fission_rate_old); + double k_eff_new = k_eff_ * (fission_rate_new / fission_rate_old); double H = 0.0; // defining an inverse sum for better performance @@ -385,7 +369,7 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const // Adds entropy value to shared entropy vector in openmc namespace. simulation::entropy.push_back(H); - return k_eff_new; + k_eff_ = k_eff_new; } // This function is responsible for generating a mapping between random @@ -797,7 +781,11 @@ void FlatSourceDomain::output_to_vtk() const int i_cell = p.lowest_coord().cell(); int64_t sr = source_region_offsets_[i_cell] + p.cell_instance(); if (RandomRay::mesh_subdivision_enabled_) { - int mesh_idx = base_source_regions_.mesh(sr); + int mesh_idx = C_NONE; + auto mesh_it = mesh_map_.find(sr); + if (mesh_it != mesh_map_.end()) { + mesh_idx = mesh_it->second; + } int mesh_bin; if (mesh_idx == C_NONE) { mesh_bin = 0; @@ -967,13 +955,17 @@ void FlatSourceDomain::output_to_vtk() const } void FlatSourceDomain::apply_external_source_to_source_region( - Discrete* discrete, double strength_factor, SourceRegionHandle& srh) + int src_idx, SourceRegionHandle& srh) { - srh.external_source_present() = 1; - + auto s = model::external_sources[src_idx].get(); + auto is = dynamic_cast(s); + auto discrete = dynamic_cast(is->energy()); + double strength_factor = is->strength(); const auto& discrete_energies = discrete->x(); const auto& discrete_probs = discrete->prob(); + srh.external_source_present() = 1; + for (int i = 0; i < discrete_energies.size(); i++) { int g = data::mg.get_group_index(discrete_energies[i]); srh.external_source(g) += discrete_probs[i] * strength_factor; @@ -981,8 +973,7 @@ void FlatSourceDomain::apply_external_source_to_source_region( } void FlatSourceDomain::apply_external_source_to_cell_instances(int32_t i_cell, - Discrete* discrete, double strength_factor, int target_material_id, - const vector& instances) + int src_idx, int target_material_id, const vector& instances) { Cell& cell = *model::cells[i_cell]; @@ -1000,16 +991,13 @@ void FlatSourceDomain::apply_external_source_to_cell_instances(int32_t i_cell, if (target_material_id == C_NONE || cell_material_id == target_material_id) { int64_t source_region = source_region_offsets_[i_cell] + j; - SourceRegionHandle srh = - source_regions_.get_source_region_handle(source_region); - apply_external_source_to_source_region(discrete, strength_factor, srh); + external_volumetric_source_map_[source_region].push_back(src_idx); } } } void FlatSourceDomain::apply_external_source_to_cell_and_children( - int32_t i_cell, Discrete* discrete, double strength_factor, - int32_t target_material_id) + int32_t i_cell, int src_idx, int32_t target_material_id) { Cell& cell = *model::cells[i_cell]; @@ -1017,14 +1005,14 @@ void FlatSourceDomain::apply_external_source_to_cell_and_children( vector instances(cell.n_instances()); std::iota(instances.begin(), instances.end(), 0); apply_external_source_to_cell_instances( - i_cell, discrete, strength_factor, target_material_id, instances); + i_cell, src_idx, target_material_id, instances); } else if (target_material_id == C_NONE) { std::unordered_map> cell_instance_list = cell.get_contained_cells(0, nullptr); for (const auto& pair : cell_instance_list) { int32_t i_child_cell = pair.first; - apply_external_source_to_cell_instances(i_child_cell, discrete, - strength_factor, target_material_id, pair.second); + apply_external_source_to_cell_instances( + i_child_cell, src_idx, target_material_id, pair.second); } } } @@ -1073,33 +1061,25 @@ void FlatSourceDomain::convert_external_sources() int i_cell = gs.lowest_coord().cell(); int64_t sr = source_region_offsets_[i_cell] + gs.cell_instance(); - if (RandomRay::mesh_subdivision_enabled_) { - // If mesh subdivision is enabled, we need to determine which subdivided - // mesh bin the point source coordinate is in as well - int mesh_idx = source_regions_.mesh(sr); - int mesh_bin; - if (mesh_idx == C_NONE) { - mesh_bin = 0; - } else { - mesh_bin = model::meshes[mesh_idx]->get_bin(gs.r()); - } - // With the source region and mesh bin known, we can use the - // accompanying SourceRegionKey as a key into a map that stores the - // corresponding external source index for the point source. Notably, we - // do not actually apply the external source to any source regions here, - // as if mesh subdivision is enabled, they haven't actually been - // discovered & initilized yet. When discovered, they will read from the - // point_source_map to determine if there are any point source terms - // that should be applied. - SourceRegionKey key {sr, mesh_bin}; - point_source_map_[key] = es; + // If mesh subdivision is enabled, we need to determine which subdivided + // mesh bin the point source coordinate is in as well + int mesh_idx = source_regions_.mesh(sr); + int mesh_bin; + if (mesh_idx == C_NONE) { + mesh_bin = 0; } else { - // If we are not using mesh subdivision, we can apply the external - // source directly to the source region as we do for volumetric domain - // constraint sources. - SourceRegionHandle srh = source_regions_.get_source_region_handle(sr); - apply_external_source_to_source_region(energy, strength_factor, srh); + mesh_bin = model::meshes[mesh_idx]->get_bin(gs.r()); } + // With the source region and mesh bin known, we can use the + // accompanying SourceRegionKey as a key into a map that stores the + // corresponding external source index for the point source. Notably, we + // do not actually apply the external source to any source regions here, + // as if mesh subdivision is enabled, they haven't actually been + // discovered & initilized yet. When discovered, they will read from the + // external_source_map to determine if there are any external source + // terms that should be applied. + SourceRegionKey key {sr, mesh_bin}; + external_point_source_map_[key].push_back(es); } else { // If not a point source, then use the volumetric domain constraints to @@ -1107,42 +1087,41 @@ void FlatSourceDomain::convert_external_sources() if (is->domain_type() == Source::DomainType::MATERIAL) { for (int32_t material_id : domain_ids) { for (int i_cell = 0; i_cell < model::cells.size(); i_cell++) { - apply_external_source_to_cell_and_children( - i_cell, energy, strength_factor, material_id); + apply_external_source_to_cell_and_children(i_cell, es, material_id); } } } else if (is->domain_type() == Source::DomainType::CELL) { for (int32_t cell_id : domain_ids) { int32_t i_cell = model::cell_map[cell_id]; - apply_external_source_to_cell_and_children( - i_cell, energy, strength_factor, C_NONE); + apply_external_source_to_cell_and_children(i_cell, es, C_NONE); } } else if (is->domain_type() == Source::DomainType::UNIVERSE) { for (int32_t universe_id : domain_ids) { int32_t i_universe = model::universe_map[universe_id]; Universe& universe = *model::universes[i_universe]; for (int32_t i_cell : universe.cells_) { - apply_external_source_to_cell_and_children( - i_cell, energy, strength_factor, C_NONE); + apply_external_source_to_cell_and_children(i_cell, es, C_NONE); } } } } } // End loop over external sources -// Divide the fixed source term by sigma t (to save time when applying each -// iteration) -#pragma omp parallel for - for (int64_t sr = 0; sr < n_source_regions(); sr++) { - int material = source_regions_.material(sr); - if (material == MATERIAL_VOID) { - continue; - } - for (int g = 0; g < negroups_; g++) { - double sigma_t = sigma_t_[material * negroups_ + g]; - source_regions_.external_source(sr, g) /= sigma_t; + // Divide the fixed source term by sigma t (to save time when applying each + // iteration) + /* + #pragma omp parallel for + for (int64_t sr = 0; sr < n_source_regions(); sr++) { + int material = source_regions_.material(sr); + if (material == MATERIAL_VOID) { + continue; + } + for (int g = 0; g < negroups_; g++) { + double sigma_t = sigma_t_[material * negroups_ + g]; + source_regions_.external_source(sr, g) /= sigma_t; + } } - } + */ } void FlatSourceDomain::flux_swap() @@ -1326,13 +1305,14 @@ void FlatSourceDomain::apply_mesh_to_cell_instances(int32_t i_cell, if ((target_material_id == C_NONE && !is_target_void) || cell_material_id == target_material_id) { int64_t sr = source_region_offsets_[i_cell] + j; - if (source_regions_.mesh(sr) != C_NONE) { - // print out the source region that is broken: + // Check if the key is already present in the mesh_map_ + if (mesh_map_.find(sr) != mesh_map_.end()) { fatal_error(fmt::format("Source region {} already has mesh idx {} " "applied, but trying to apply mesh idx {}", - sr, source_regions_.mesh(sr), mesh_idx)); + sr, mesh_map_[sr], mesh_idx)); } - source_regions_.mesh(sr) = mesh_idx; + // If the SR has not already been assigned, then we can write to it + mesh_map_[sr] = mesh_idx; } } } @@ -1485,7 +1465,11 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( } // Sanity check on mesh bin - int mesh_idx = base_source_regions_.mesh(sr); + int mesh_idx = C_NONE; + auto mesh_it = mesh_map_.find(sr); + if (mesh_it != mesh_map_.end()) { + mesh_idx = mesh_it->second; + } if (mesh_idx == C_NONE) { if (mesh_bin != 0) { discovered_source_regions_.unlock(sr_key); @@ -1512,29 +1496,94 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( // to inherit material, external source, and some flux properties etc. We // also pass the base source region id to allow the new source region to // know which base source region it is derived from. - SourceRegion* sr_ptr = discovered_source_regions_.emplace( - sr_key, {base_source_regions_.get_source_region_handle(sr), sr}); - discovered_source_regions_.unlock(sr_key); + + // 0. Call the basic constructor for the source region + bool is_linear = RandomRay::source_shape_ != RandomRaySourceShape::FLAT; + SourceRegion* sr_ptr = + discovered_source_regions_.emplace(sr_key, {negroups_, is_linear}); SourceRegionHandle handle {*sr_ptr}; - // Check if the new source region contains a point source and apply it if so - auto it2 = point_source_map_.find(sr_key); - if (it2 != point_source_map_.end()) { - int es = it2->second; - auto s = model::external_sources[es].get(); - auto is = dynamic_cast(s); - auto energy = dynamic_cast(is->energy()); - double strength_factor = is->strength(); - apply_external_source_to_source_region(energy, strength_factor, handle); - int material = handle.material(); - if (material != MATERIAL_VOID) { - for (int g = 0; g < negroups_; g++) { - double sigma_t = sigma_t_[material * negroups_ + g]; - handle.external_source(g) /= sigma_t; + // 1. Lock the SR (before writing anything to it) + handle.lock(); + discovered_source_regions_.unlock(sr_key); + + // 2. Determine the material + Cell& cell = *model::cells[gs_i_cell]; + int material = cell.material(gs.cell_instance()); + handle.material() = material; + + // 3. Determine if there are any meshes assigned to this + handle.mesh() = mesh_idx; + + // 4. Determine if there are any point sources, and apply them. + // Point sources are specific to the SRK. + auto it4 = external_point_source_map_.find(sr_key); + if (it4 != external_point_source_map_.end()) { + handle.external_source_present() = 1; + const vector& point_sources = it4->second; + for (int src_idx : point_sources) { + apply_external_source_to_source_region(src_idx, handle); + } + } + + // 5. Determine if there are any volumetric sources, and apply them. + // Volumetric sources are specifc only to the base SR idx. + auto it5 = external_volumetric_source_map_.find(sr); + if (it5 != external_volumetric_source_map_.end()) { + handle.external_source_present() = 1; + const vector& vol_sources = it5->second; + for (int src_idx : vol_sources) { + apply_external_source_to_source_region(src_idx, handle); + } + } + + // 6. Initialize scalar flux based off of 1.0 or 0.0 depending on k vs. fixed, + // and if external source is present or not. + for (int g = 0; g < negroups_; g++) { + if (settings::run_mode == RunMode::FIXED_SOURCE && + handle.external_source_present()) { + handle.scalar_flux_old(g) = 1.0; + } + } + + // 7. Divide external source by sigma_t + if (material != C_NONE) { + for (int g = 0; g < negroups_; g++) { + double sigma_t = sigma_t_[material * negroups_ + g]; + handle.external_source(g) /= sigma_t; + } + } + + // 8. Compute the scattering + fission and divide by sigma_t + if (material != MATERIAL_VOID) { + for (int g_out = 0; g_out < negroups_; g_out++) { + double sigma_t = sigma_t_[material * negroups_ + g_out]; + double scatter_source = 0.0; + double fission_source = 0.0; + + for (int g_in = 0; g_in < negroups_; g_in++) { + double scalar_flux = handle.scalar_flux_old(g_in); + double sigma_s = + sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in]; + double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in]; + double chi = chi_[material * negroups_ + g_out]; + + scatter_source += sigma_s * scalar_flux; + fission_source += nu_sigma_f * scalar_flux * chi; } + handle.source(g_out) = + (scatter_source + fission_source / k_eff_) / sigma_t; } } + // 9. Add in external source + for (int g = 0; g < negroups_; g++) { + handle.source(g) += handle.external_source(g); + } + + // 10. Unlock the SR + handle.unlock(); + return handle; } diff --git a/src/random_ray/linear_source_domain.cpp b/src/random_ray/linear_source_domain.cpp index 81412164eca..3a34478aab2 100644 --- a/src/random_ray/linear_source_domain.cpp +++ b/src/random_ray/linear_source_domain.cpp @@ -34,11 +34,11 @@ void LinearSourceDomain::batch_reset() } } -void LinearSourceDomain::update_neutron_source(double k_eff) +void LinearSourceDomain::update_neutron_source() { simulation::time_update_src.start(); - double inverse_k_eff = 1.0 / k_eff; + double inverse_k_eff = 1.0 / k_eff_; // Reset all source regions to zero (important for void regions) #pragma omp parallel for diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index 27f674c4278..8c72f06f995 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -345,7 +345,11 @@ void RandomRay::attenuate_flux(double distance, bool is_active, double offset) // Perform ray tracing across mesh if (mesh_subdivision_enabled_) { // Determine the mesh index for the base source region, if any - int mesh_idx = domain_->base_source_regions_.mesh(sr); + int mesh_idx = C_NONE; + auto it = domain_->mesh_map_.find(sr); + if (it != domain_->mesh_map_.end()) { + mesh_idx = it->second; + } if (mesh_idx == C_NONE) { // If there's no mesh being applied to this cell, then @@ -816,20 +820,22 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain) int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance(); SourceRegionHandle srh; - if (mesh_subdivision_enabled_) { - int mesh_idx = domain_->base_source_regions_.mesh(sr); - int mesh_bin; - if (mesh_idx == C_NONE) { - mesh_bin = 0; - } else { - Mesh* mesh = model::meshes[mesh_idx].get(); - mesh_bin = mesh->get_bin(r()); - } - srh = - domain_->get_subdivided_source_region_handle(sr, mesh_bin, r(), 0.0, u()); + + // Determine if there are any meshes assigned to this + int mesh_idx = C_NONE; + auto it = domain_->mesh_map_.find(sr); + if (it != domain_->mesh_map_.end()) { + mesh_idx = it->second; + } + int mesh_bin; + if (mesh_idx == C_NONE) { + mesh_bin = 0; } else { - srh = domain_->source_regions_.get_source_region_handle(sr); + Mesh* mesh = model::meshes[mesh_idx].get(); + mesh_bin = mesh->get_bin(r()); } + srh = + domain_->get_subdivided_source_region_handle(sr, mesh_bin, r(), 0.0, u()); if (!srh.is_numerical_fp_artifact_) { for (int g = 0; g < negroups_; g++) { diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 388a778b840..cf9706675f6 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -444,12 +444,13 @@ void RandomRaySimulation::simulate() // 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(); + if (simulation::current_batch > 1) { + // Update source term (scattering + fission) + domain_->update_neutron_source(); + // Reset scalar fluxes, iteration volume tallies, and region hit flags + // to zero + domain_->batch_reset(); + } // At the beginning of the simulation, if mesh subvivision is in use, we // need to swap the main source region container into the base container, @@ -459,7 +460,7 @@ void RandomRaySimulation::simulate() // material properties, and initial guess values for the flux/source. if (RandomRay::mesh_subdivision_enabled_ && simulation::current_batch == 1 && !FlatSourceDomain::adjoint_) { - domain_->prepare_base_source_regions(); + // domain_->prepare_base_source_regions(); } // Start timer for transport @@ -494,10 +495,10 @@ void RandomRaySimulation::simulate() if (settings::run_mode == RunMode::EIGENVALUE) { // Compute random ray k-eff - k_eff_ = domain_->compute_k_eff(k_eff_); + domain_->compute_k_eff(); // Store random ray k-eff into OpenMC's native k-eff variable - global_tally_tracklength = k_eff_; + global_tally_tracklength = domain_->k_eff_; } // Execute all tallying tasks, if this is an active batch @@ -522,7 +523,7 @@ void RandomRaySimulation::simulate() domain_->flux_swap(); // Check for any obvious insabilities/nans/infs - instability_check(n_hits, k_eff_, avg_miss_rate_); + instability_check(n_hits, domain_->k_eff_, avg_miss_rate_); } // End MPI master work // Finalize the current batch diff --git a/src/random_ray/source_region.cpp b/src/random_ray/source_region.cpp index 1205b995a6f..e0eb090900b 100644 --- a/src/random_ray/source_region.cpp +++ b/src/random_ray/source_region.cpp @@ -48,7 +48,7 @@ SourceRegion::SourceRegion(int negroups, bool is_linear) } scalar_flux_new_.assign(negroups, 0.0); - source_.resize(negroups); + source_.assign(negroups, 0.0); scalar_flux_final_.assign(negroups, 0.0); tally_task_.resize(negroups); From 8719e23f2c1a3a8f89d0b09368e646cde1520bbc Mon Sep 17 00:00:00 2001 From: John Tramm Date: Thu, 18 Sep 2025 15:46:08 -0500 Subject: [PATCH 02/19] seems to be working --- .../openmc/random_ray/flat_source_domain.h | 18 ++-- include/openmc/random_ray/random_ray.h | 1 - src/random_ray/flat_source_domain.cpp | 48 +++++----- src/random_ray/random_ray.cpp | 89 +++++++++---------- src/random_ray/random_ray_simulation.cpp | 30 ++----- src/settings.cpp | 1 - 6 files changed, 80 insertions(+), 107 deletions(-) diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index e61f37abe22..20bd0f1307f 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -86,7 +86,7 @@ class FlatSourceDomain { //---------------------------------------------------------------------------- // Public Data members - double k_eff_ {1.0}; // Eigenvalue + double k_eff_ {1.0}; // Eigenvalue bool mapped_all_tallies_ {false}; // If all source regions have been visited int64_t n_external_source_regions_ {0}; // Total number of source regions with @@ -141,13 +141,12 @@ class FlatSourceDomain { // Map that relates a base source region index to the external source index. // This map is used to check if there are any volumetric sources within a // subdivided source region at the time it is discovered. - std::unordered_map, SourceRegionKey::HashFunctor> + std::unordered_map> external_volumetric_source_map_; // Map that relates a base source region index to a mesh index. This map // is used to check which subdivision mesh is present in a source region. - std::unordered_map - mesh_map_; + std::unordered_map mesh_map_; // If transport corrected MGXS data is being used, there may be negative // in-group scattering cross sections that can result in instability in MOC @@ -158,13 +157,12 @@ class FlatSourceDomain { protected: //---------------------------------------------------------------------------- // Methods - void FlatSourceDomain::apply_external_source_to_source_region( + void apply_external_source_to_source_region( int src_idx, SourceRegionHandle& srh); - void apply_external_source_to_cell_instances(int32_t i_cell, - int src_idx, int target_material_id, - const vector& instances); - void apply_external_source_to_cell_and_children(int32_t i_cell, - int src_idx, int32_t target_material_id); + void apply_external_source_to_cell_instances(int32_t i_cell, int src_idx, + int target_material_id, const vector& instances); + void apply_external_source_to_cell_and_children( + int32_t i_cell, int src_idx, int32_t target_material_id); virtual void set_flux_to_flux_plus_source(int64_t sr, double volume, int g); void set_flux_to_source(int64_t sr, int g); virtual void set_flux_to_old_flux(int64_t sr, int g); diff --git a/include/openmc/random_ray/random_ray.h b/include/openmc/random_ray/random_ray.h index abf2a26881f..40c67ef9549 100644 --- a/include/openmc/random_ray/random_ray.h +++ b/include/openmc/random_ray/random_ray.h @@ -48,7 +48,6 @@ class RandomRay : public Particle { static double distance_active_; // Active ray length 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 //---------------------------------------------------------------------------- diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 69db68aef11..f9b63e586a5 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -53,8 +53,8 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_) // Initialize source regions. bool is_linear = RandomRay::source_shape_ != RandomRaySourceShape::FLAT; source_regions_ = SourceRegionContainer(negroups_, is_linear); - source_regions_.assign( - base_source_regions, SourceRegion(negroups_, is_linear)); + //source_regions_.assign( + // base_source_regions, SourceRegion(negroups_, is_linear)); // Initialize tally volumes if (volume_normalized_flux_tallies_) { @@ -780,25 +780,23 @@ void FlatSourceDomain::output_to_vtk() const int i_cell = p.lowest_coord().cell(); int64_t sr = source_region_offsets_[i_cell] + p.cell_instance(); - if (RandomRay::mesh_subdivision_enabled_) { - int mesh_idx = C_NONE; - auto mesh_it = mesh_map_.find(sr); - if (mesh_it != mesh_map_.end()) { - mesh_idx = mesh_it->second; - } - int mesh_bin; - if (mesh_idx == C_NONE) { - mesh_bin = 0; - } else { - mesh_bin = model::meshes[mesh_idx]->get_bin(p.r()); - } - SourceRegionKey sr_key {sr, mesh_bin}; - auto it = source_region_map_.find(sr_key); - if (it != source_region_map_.end()) { - sr = it->second; - } else { - sr = -1; - } + int mesh_idx = C_NONE; + auto mesh_it = mesh_map_.find(sr); + if (mesh_it != mesh_map_.end()) { + mesh_idx = mesh_it->second; + } + int mesh_bin; + if (mesh_idx == C_NONE) { + mesh_bin = 0; + } else { + mesh_bin = model::meshes[mesh_idx]->get_bin(p.r()); + } + SourceRegionKey sr_key {sr, mesh_bin}; + auto it = source_region_map_.find(sr_key); + if (it != source_region_map_.end()) { + sr = it->second; + } else { + sr = -1; } voxel_indices[z * Ny * Nx + y * Nx + x] = sr; @@ -1547,7 +1545,7 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( } // 7. Divide external source by sigma_t - if (material != C_NONE) { + if (material != C_NONE && settings::run_mode == RunMode::FIXED_SOURCE) { for (int g = 0; g < negroups_; g++) { double sigma_t = sigma_t_[material * negroups_ + g]; handle.external_source(g) /= sigma_t; @@ -1577,8 +1575,10 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( } // 9. Add in external source - for (int g = 0; g < negroups_; g++) { - handle.source(g) += handle.external_source(g); + if (settings::run_mode == RunMode::FIXED_SOURCE) { + for (int g = 0; g < negroups_; g++) { + handle.source(g) += handle.external_source(g); + } } // 10. Unlock the SR diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index 8c72f06f995..c0e42a88547 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -237,7 +237,6 @@ double RandomRay::distance_inactive_; double RandomRay::distance_active_; unique_ptr RandomRay::ray_source_; RandomRaySourceShape RandomRay::source_shape_ {RandomRaySourceShape::FLAT}; -bool RandomRay::mesh_subdivision_enabled_ {false}; RandomRaySampleMethod RandomRay::sample_method_ {RandomRaySampleMethod::PRNG}; RandomRay::RandomRay() @@ -343,48 +342,44 @@ void RandomRay::attenuate_flux(double distance, bool is_active, double offset) int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance(); // Perform ray tracing across mesh - if (mesh_subdivision_enabled_) { - // Determine the mesh index for the base source region, if any - int mesh_idx = C_NONE; - auto it = domain_->mesh_map_.find(sr); - if (it != domain_->mesh_map_.end()) { - mesh_idx = it->second; - } + // Determine the mesh index for the base source region, if any + int mesh_idx = C_NONE; + auto it = domain_->mesh_map_.find(sr); + if (it != domain_->mesh_map_.end()) { + mesh_idx = it->second; + } - if (mesh_idx == C_NONE) { - // If there's no mesh being applied to this cell, then - // we just attenuate the flux as normal, and set - // the mesh bin to 0 - attenuate_flux_inner(distance, is_active, sr, 0, r()); - } else { - // If there is a mesh being applied to this cell, then - // we loop over all the bin crossings and attenuate - // separately. - Mesh* mesh = model::meshes[mesh_idx].get(); - - // We adjust the start and end positions of the ray slightly - // to accomodate for floating point precision issues that tend - // to occur at mesh boundaries that overlap with geometry lattice - // boundaries. - Position start = r() + (offset + TINY_BIT) * u(); - Position end = start + (distance - 2.0 * TINY_BIT) * u(); - double reduced_distance = (end - start).norm(); - - // Ray trace through the mesh and record bins and lengths - mesh_bins_.resize(0); - mesh_fractional_lengths_.resize(0); - mesh->bins_crossed(start, end, u(), mesh_bins_, mesh_fractional_lengths_); - - // Loop over all mesh bins and attenuate flux - for (int b = 0; b < mesh_bins_.size(); b++) { - double physical_length = reduced_distance * mesh_fractional_lengths_[b]; - attenuate_flux_inner( - physical_length, is_active, sr, mesh_bins_[b], start); - start += physical_length * u(); - } - } + if (mesh_idx == C_NONE) { + // If there's no mesh being applied to this cell, then + // we just attenuate the flux as normal, and set + // the mesh bin to 0 + attenuate_flux_inner(distance, is_active, sr, 0, r()); } else { - attenuate_flux_inner(distance, is_active, sr, C_NONE, r()); + // If there is a mesh being applied to this cell, then + // we loop over all the bin crossings and attenuate + // separately. + Mesh* mesh = model::meshes[mesh_idx].get(); + + // We adjust the start and end positions of the ray slightly + // to accomodate for floating point precision issues that tend + // to occur at mesh boundaries that overlap with geometry lattice + // boundaries. + Position start = r() + (offset + TINY_BIT) * u(); + Position end = start + (distance - 2.0 * TINY_BIT) * u(); + double reduced_distance = (end - start).norm(); + + // Ray trace through the mesh and record bins and lengths + mesh_bins_.resize(0); + mesh_fractional_lengths_.resize(0); + mesh->bins_crossed(start, end, u(), mesh_bins_, mesh_fractional_lengths_); + + // Loop over all mesh bins and attenuate flux + for (int b = 0; b < mesh_bins_.size(); b++) { + double physical_length = reduced_distance * mesh_fractional_lengths_[b]; + attenuate_flux_inner( + physical_length, is_active, sr, mesh_bins_[b], start); + start += physical_length * u(); + } } } @@ -392,14 +387,10 @@ void RandomRay::attenuate_flux_inner( double distance, bool is_active, int64_t sr, int mesh_bin, Position r) { SourceRegionHandle srh; - if (mesh_subdivision_enabled_) { - srh = domain_->get_subdivided_source_region_handle( - sr, mesh_bin, r, distance, u()); - if (srh.is_numerical_fp_artifact_) { - return; - } - } else { - srh = domain_->source_regions_.get_source_region_handle(sr); + srh = domain_->get_subdivided_source_region_handle( + sr, mesh_bin, r, distance, u()); + if (srh.is_numerical_fp_artifact_) { + return; } switch (source_shape_) { diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index cf9706675f6..d9a60a41885 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -348,7 +348,6 @@ void validate_random_ray_inputs() // when generating weight windows with FW-CADIS and an overlaid mesh. /////////////////////////////////////////////////////////////////// if (RandomRay::source_shape_ == RandomRaySourceShape::LINEAR && - RandomRay::mesh_subdivision_enabled_ && variance_reduction::weight_windows.size() > 0) { warning( "Linear sources may result in negative fluxes in small source regions " @@ -366,7 +365,6 @@ void openmc_reset_random_ray() FlatSourceDomain::mesh_domain_map_.clear(); RandomRay::ray_source_.reset(); RandomRay::source_shape_ = RandomRaySourceShape::FLAT; - RandomRay::mesh_subdivision_enabled_ = false; RandomRay::sample_method_ = RandomRaySampleMethod::PRNG; } @@ -419,12 +417,10 @@ void RandomRaySimulation::prepare_fixed_sources_adjoint( forward_source_region_map) { if (settings::run_mode == RunMode::FIXED_SOURCE) { - if (RandomRay::mesh_subdivision_enabled_) { - domain_->source_regions_ = forward_source_regions; - domain_->source_region_map_ = forward_source_region_map; - domain_->base_source_regions_ = forward_base_source_regions; - domain_->source_regions_.adjoint_reset(); - } + domain_->source_regions_ = forward_source_regions; + domain_->source_region_map_ = forward_source_region_map; + domain_->base_source_regions_ = forward_base_source_regions; + domain_->source_regions_.adjoint_reset(); domain_->set_adjoint_sources(forward_flux); } } @@ -458,10 +454,6 @@ void RandomRaySimulation::simulate() // subdivided source regions. The base container will therefore only // contain the external source region information, the mesh indices, // material properties, and initial guess values for the flux/source. - if (RandomRay::mesh_subdivision_enabled_ && - simulation::current_batch == 1 && !FlatSourceDomain::adjoint_) { - // domain_->prepare_base_source_regions(); - } // Start timer for transport simulation::time_transport.start(); @@ -479,9 +471,9 @@ void RandomRaySimulation::simulate() // If using mesh subdivision, add any newly discovered source regions // to the main source region container. - if (RandomRay::mesh_subdivision_enabled_) { - domain_->finalize_discovered_source_regions(); - } + // if (RandomRay::mesh_subdivision_enabled_) { + domain_->finalize_discovered_source_regions(); + //} // Normalize scalar flux and update volumes domain_->normalize_scalar_flux_and_volumes( @@ -508,12 +500,6 @@ void RandomRaySimulation::simulate() // estimate domain_->accumulate_iteration_flux(); - // Generate mapping between source regions and tallies - if (!domain_->mapped_all_tallies_ && - !RandomRay::mesh_subdivision_enabled_) { - domain_->convert_source_regions_to_tallies(0); - } - // Use above mapping to contribute FSR flux data to appropriate // tallies domain_->random_ray_tally(); @@ -572,7 +558,7 @@ void RandomRaySimulation::instability_check( } if (k_eff > 10.0 || k_eff < 0.01 || !(std::isfinite(k_eff))) { - fatal_error("Instability detected"); + fatal_error(fmt::format("Instability detected: k-eff = {:.5f}", k_eff)); } } } diff --git a/src/settings.cpp b/src/settings.cpp index afa42a3c5e6..2b703c684c6 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -345,7 +345,6 @@ void get_run_parameters(pugi::xml_node node_base) } FlatSourceDomain::mesh_domain_map_[mesh_id].emplace_back( type, domain_id); - RandomRay::mesh_subdivision_enabled_ = true; } } } From e702bbec85002fa491f85dc9b0e778f85a4351d1 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Fri, 19 Sep 2025 09:04:50 -0500 Subject: [PATCH 03/19] finished deleting all references to base source regions --- include/openmc/random_ray/flat_source_domain.h | 9 --------- include/openmc/random_ray/random_ray_simulation.h | 1 - src/random_ray/flat_source_domain.cpp | 9 --------- src/random_ray/random_ray_simulation.cpp | 6 +----- 4 files changed, 1 insertion(+), 24 deletions(-) diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index 20bd0f1307f..2e40ecaeb7c 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -54,7 +54,6 @@ class FlatSourceDomain { bool is_target_void); void apply_mesh_to_cell_and_children(int32_t i_cell, int32_t mesh_idx, int32_t target_material_id, bool is_target_void); - void prepare_base_source_regions(); SourceRegionHandle get_subdivided_source_region_handle( int64_t sr, int mesh_bin, Position r, double dist, Direction u); void finalize_discovered_source_regions(); @@ -111,14 +110,6 @@ class FlatSourceDomain { // The abstract container holding all source region-specific data SourceRegionContainer source_regions_; - // Base source region container. When source region subdivision via mesh - // is in use, this container holds the original (non-subdivided) material - // filled cell instance source regions. These are useful as they can be - // initialized with external source and mesh domain information ahead of time. - // Then, dynamically discovered source regions can be initialized by cloning - // their base region. - SourceRegionContainer base_source_regions_; - // Parallel hash map holding all source regions discovered during // a single iteration. This is a threadsafe data structure that is cleaned // out after each iteration and stored in the "source_regions_" container. diff --git a/include/openmc/random_ray/random_ray_simulation.h b/include/openmc/random_ray/random_ray_simulation.h index 13bc68834ec..3d1e6158a0d 100644 --- a/include/openmc/random_ray/random_ray_simulation.h +++ b/include/openmc/random_ray/random_ray_simulation.h @@ -23,7 +23,6 @@ class RandomRaySimulation { void apply_fixed_sources_and_mesh_domains(); void prepare_fixed_sources_adjoint(vector& forward_flux, SourceRegionContainer& forward_source_regions, - SourceRegionContainer& forward_base_source_regions, std::unordered_map& forward_source_region_map); void simulate(); diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index f9b63e586a5..8db512a8120 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -53,8 +53,6 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_) // Initialize source regions. bool is_linear = RandomRay::source_shape_ != RandomRaySourceShape::FLAT; source_regions_ = SourceRegionContainer(negroups_, is_linear); - //source_regions_.assign( - // base_source_regions, SourceRegion(negroups_, is_linear)); // Initialize tally volumes if (volume_normalized_flux_tallies_) { @@ -1380,13 +1378,6 @@ void FlatSourceDomain::apply_meshes() } } -void FlatSourceDomain::prepare_base_source_regions() -{ - std::swap(source_regions_, base_source_regions_); - source_regions_.negroups() = base_source_regions_.negroups(); - source_regions_.is_linear() = base_source_regions_.is_linear(); -} - SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( int64_t sr, int mesh_bin, Position r, double dist, Direction u) { diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index d9a60a41885..8877c7e111c 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -50,7 +50,6 @@ void openmc_run_random_ray() // Declare forward flux so that it can be saved for later adjoint simulation vector forward_flux; SourceRegionContainer forward_source_regions; - SourceRegionContainer forward_base_source_regions; std::unordered_map forward_source_region_map; @@ -84,7 +83,6 @@ void openmc_run_random_ray() forward_source_regions = sim.domain()->source_regions_; forward_source_region_map = sim.domain()->source_region_map_; - forward_base_source_regions = sim.domain()->base_source_regions_; // Finalize OpenMC openmc_simulation_finalize(); @@ -114,7 +112,7 @@ void openmc_run_random_ray() // Initialize adjoint fixed sources, if present adjoint_sim.prepare_fixed_sources_adjoint(forward_flux, - forward_source_regions, forward_base_source_regions, + forward_source_regions, forward_source_region_map); // Transpose scattering matrix @@ -412,14 +410,12 @@ void RandomRaySimulation::apply_fixed_sources_and_mesh_domains() void RandomRaySimulation::prepare_fixed_sources_adjoint( vector& forward_flux, SourceRegionContainer& forward_source_regions, - SourceRegionContainer& forward_base_source_regions, std::unordered_map& forward_source_region_map) { if (settings::run_mode == RunMode::FIXED_SOURCE) { domain_->source_regions_ = forward_source_regions; domain_->source_region_map_ = forward_source_region_map; - domain_->base_source_regions_ = forward_base_source_regions; domain_->source_regions_.adjoint_reset(); domain_->set_adjoint_sources(forward_flux); } From 2990038f19c76fb487b4fc8fa2e098a3ba66c720 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Fri, 19 Sep 2025 09:33:43 -0500 Subject: [PATCH 04/19] fixed issue with mesh map for point source locator --- src/random_ray/flat_source_domain.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 8db512a8120..35467aaf3fa 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -1059,7 +1059,11 @@ void FlatSourceDomain::convert_external_sources() // If mesh subdivision is enabled, we need to determine which subdivided // mesh bin the point source coordinate is in as well - int mesh_idx = source_regions_.mesh(sr); + int mesh_idx = C_NONE; + auto mesh_it = mesh_map_.find(sr); + if (mesh_it != mesh_map_.end()) { + mesh_idx = mesh_it->second; + } int mesh_bin; if (mesh_idx == C_NONE) { mesh_bin = 0; From c7df7ebf79d9bab4fb72f9618a327a485d1db13d Mon Sep 17 00:00:00 2001 From: John Tramm Date: Fri, 19 Sep 2025 10:18:45 -0500 Subject: [PATCH 05/19] moved source updates into single SR function, to allow for calling at either batch or init on single source regions --- .../openmc/random_ray/flat_source_domain.h | 3 +- .../openmc/random_ray/linear_source_domain.h | 2 +- include/openmc/random_ray/moment_matrix.h | 12 +-- src/random_ray/flat_source_domain.cpp | 76 +++++++------------ src/random_ray/linear_source_domain.cpp | 46 +++++------ src/random_ray/random_ray_simulation.cpp | 3 +- 6 files changed, 55 insertions(+), 87 deletions(-) diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index 2e40ecaeb7c..1a167bb64e6 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -27,7 +27,8 @@ class FlatSourceDomain { //---------------------------------------------------------------------------- // Methods - virtual void update_neutron_source(); + virtual void update_single_neutron_source(SourceRegionHandle& srh); + virtual void update_all_neutron_sources(); void compute_k_eff(); virtual void normalize_scalar_flux_and_volumes( double total_active_distance_per_iteration); diff --git a/include/openmc/random_ray/linear_source_domain.h b/include/openmc/random_ray/linear_source_domain.h index 7b5fd86dc52..0098c782001 100644 --- a/include/openmc/random_ray/linear_source_domain.h +++ b/include/openmc/random_ray/linear_source_domain.h @@ -20,7 +20,7 @@ class LinearSourceDomain : public FlatSourceDomain { public: //---------------------------------------------------------------------------- // Methods - void update_neutron_source() override; + void update_single_neutron_source(SourceRegionHandle& srh) override; void normalize_scalar_flux_and_volumes( double total_active_distance_per_iteration) override; diff --git a/include/openmc/random_ray/moment_matrix.h b/include/openmc/random_ray/moment_matrix.h index c95bb2c1286..508b562a478 100644 --- a/include/openmc/random_ray/moment_matrix.h +++ b/include/openmc/random_ray/moment_matrix.h @@ -26,12 +26,12 @@ class MomentMatrix { public: //---------------------------------------------------------------------------- // Public data members - double a; - double b; - double c; - double d; - double e; - double f; + double a {0}; + double b {0}; + double c {0}; + double d {0}; + double e {0}; + double f {0}; //---------------------------------------------------------------------------- // Constructors diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 35467aaf3fa..58663552094 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -100,34 +100,24 @@ void FlatSourceDomain::accumulate_iteration_flux() } } -// Compute new estimate of scattering + fission sources in each source region -// based on the flux estimate from the previous iteration. -void FlatSourceDomain::update_neutron_source() +void FlatSourceDomain::update_single_neutron_source(SourceRegionHandle& srh) { - simulation::time_update_src.start(); - - double inverse_k_eff = 1.0 / k_eff_; - -// Reset all source regions to zero (important for void regions) -#pragma omp parallel for - for (int64_t se = 0; se < n_source_elements(); se++) { - source_regions_.source(se) = 0.0; + // Reset all source regions to zero (important for void regions) + for (int g = 0; g < negroups_; g++) { + srh.source(g) = 0.0; } // Add scattering + fission source -#pragma omp parallel for - for (int64_t sr = 0; sr < n_source_regions(); sr++) { - int material = source_regions_.material(sr); - if (material == MATERIAL_VOID) { - continue; - } + int material = srh.material(); + if (material != MATERIAL_VOID) { + double inverse_k_eff = 1.0 / k_eff_; for (int g_out = 0; g_out < negroups_; g_out++) { double sigma_t = sigma_t_[material * negroups_ + g_out]; double scatter_source = 0.0; double fission_source = 0.0; for (int g_in = 0; g_in < negroups_; g_in++) { - double scalar_flux = source_regions_.scalar_flux_old(sr, g_in); + double scalar_flux = srh.scalar_flux_old(g_in); double sigma_s = sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in]; double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in]; @@ -136,18 +126,30 @@ void FlatSourceDomain::update_neutron_source() scatter_source += sigma_s * scalar_flux; fission_source += nu_sigma_f * scalar_flux * chi; } - source_regions_.source(sr, g_out) = + srh.source(g_out) = (scatter_source + fission_source * inverse_k_eff) / sigma_t; } } // Add external source if in fixed source mode if (settings::run_mode == RunMode::FIXED_SOURCE) { -#pragma omp parallel for - for (int64_t se = 0; se < n_source_elements(); se++) { - source_regions_.source(se) += source_regions_.external_source(se); + for (int g = 0; g < negroups_; g++) { + srh.source(g) += srh.external_source(g); } } +} + +// Compute new estimate of scattering + fission sources in each source region +// based on the flux estimate from the previous iteration. +void FlatSourceDomain::update_all_neutron_sources() +{ + simulation::time_update_src.start(); + + #pragma omp parallel for + for (int64_t sr = 0; sr < n_source_regions(); sr++) { + SourceRegionHandle srh = source_regions_.get_source_region_handle(sr); + update_single_neutron_source(srh); + } simulation::time_update_src.stop(); } @@ -1547,34 +1549,8 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( } } - // 8. Compute the scattering + fission and divide by sigma_t - if (material != MATERIAL_VOID) { - for (int g_out = 0; g_out < negroups_; g_out++) { - double sigma_t = sigma_t_[material * negroups_ + g_out]; - double scatter_source = 0.0; - double fission_source = 0.0; - - for (int g_in = 0; g_in < negroups_; g_in++) { - double scalar_flux = handle.scalar_flux_old(g_in); - double sigma_s = - sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in]; - double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in]; - double chi = chi_[material * negroups_ + g_out]; - - scatter_source += sigma_s * scalar_flux; - fission_source += nu_sigma_f * scalar_flux * chi; - } - handle.source(g_out) = - (scatter_source + fission_source / k_eff_) / sigma_t; - } - } - - // 9. Add in external source - if (settings::run_mode == RunMode::FIXED_SOURCE) { - for (int g = 0; g < negroups_; g++) { - handle.source(g) += handle.external_source(g); - } - } + // Compute the combined source term + update_single_neutron_source(handle); // 10. Unlock the SR handle.unlock(); diff --git a/src/random_ray/linear_source_domain.cpp b/src/random_ray/linear_source_domain.cpp index 3a34478aab2..ea97e61be59 100644 --- a/src/random_ray/linear_source_domain.cpp +++ b/src/random_ray/linear_source_domain.cpp @@ -34,25 +34,18 @@ void LinearSourceDomain::batch_reset() } } -void LinearSourceDomain::update_neutron_source() +void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh) { - simulation::time_update_src.start(); - - double inverse_k_eff = 1.0 / k_eff_; - -// Reset all source regions to zero (important for void regions) -#pragma omp parallel for - for (int64_t se = 0; se < n_source_elements(); se++) { - source_regions_.source(se) = 0.0; + // Reset all source regions to zero (important for void regions) + for (int g = 0; g < negroups_; g++) { + srh.source(g) = 0.0; } -#pragma omp parallel for - for (int64_t sr = 0; sr < n_source_regions(); sr++) { - int material = source_regions_.material(sr); - if (material == MATERIAL_VOID) { - continue; - } - MomentMatrix invM = source_regions_.mom_matrix(sr).inverse(); + // Add scattering + fission source + int material = srh.material(); + if (material != MATERIAL_VOID) { + double inverse_k_eff = 1.0 / k_eff_; + MomentMatrix invM = srh.mom_matrix().inverse(); for (int g_out = 0; g_out < negroups_; g_out++) { double sigma_t = sigma_t_[material * negroups_ + g_out]; @@ -64,8 +57,8 @@ void LinearSourceDomain::update_neutron_source() for (int g_in = 0; g_in < negroups_; g_in++) { // Handles for the flat and linear components of the flux - double flux_flat = source_regions_.scalar_flux_old(sr, g_in); - MomentArray flux_linear = source_regions_.flux_moments_old(sr, g_in); + double flux_flat = srh.scalar_flux_old(g_in); + MomentArray flux_linear = srh.flux_moments_old(g_in); // Handles for cross sections double sigma_s = @@ -81,7 +74,7 @@ void LinearSourceDomain::update_neutron_source() } // Compute the flat source term - source_regions_.source(sr, g_out) = + srh.source(g_out) = (scatter_flat + fission_flat * inverse_k_eff) / sigma_t; // Compute the linear source terms. In the first 10 iterations when the @@ -92,24 +85,21 @@ void LinearSourceDomain::update_neutron_source() // the source gradients (effectively making this a flat source region // temporarily), so as to improve stability. if (simulation::current_batch > 10 && - source_regions_.source(sr, g_out) >= 0.0) { - source_regions_.source_gradients(sr, g_out) = + srh.source(g_out) >= 0.0) { + srh.source_gradients(g_out) = invM * ((scatter_linear + fission_linear * inverse_k_eff) / sigma_t); } else { - source_regions_.source_gradients(sr, g_out) = {0.0, 0.0, 0.0}; + srh.source_gradients(g_out) = {0.0, 0.0, 0.0}; } } } + // Add external source if in fixed source mode if (settings::run_mode == RunMode::FIXED_SOURCE) { -// Add external source to flat source term if in fixed source mode -#pragma omp parallel for - for (int64_t se = 0; se < n_source_elements(); se++) { - source_regions_.source(se) += source_regions_.external_source(se); + for (int g = 0; g < negroups_; g++) { + srh.source(g) += srh.external_source(g); } } - - simulation::time_update_src.stop(); } void LinearSourceDomain::normalize_scalar_flux_and_volumes( diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 8877c7e111c..3068be2cb1c 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -413,6 +413,7 @@ void RandomRaySimulation::prepare_fixed_sources_adjoint( std::unordered_map& forward_source_region_map) { + domain_->k_eff_ = 1.0; if (settings::run_mode == RunMode::FIXED_SOURCE) { domain_->source_regions_ = forward_source_regions; domain_->source_region_map_ = forward_source_region_map; @@ -438,7 +439,7 @@ void RandomRaySimulation::simulate() if (simulation::current_batch > 1) { // Update source term (scattering + fission) - domain_->update_neutron_source(); + domain_->update_all_neutron_sources(); // Reset scalar fluxes, iteration volume tallies, and region hit flags // to zero domain_->batch_reset(); From d3fd941781ce68d9b15f52287ef287271f76e7e7 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Fri, 19 Sep 2025 10:54:43 -0500 Subject: [PATCH 06/19] moved volume src before pt source to be consistent with original ordering --- src/random_ray/flat_source_domain.cpp | 58 ++++++++++-------------- src/random_ray/random_ray_simulation.cpp | 11 ++--- 2 files changed, 27 insertions(+), 42 deletions(-) diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 58663552094..ee566bb676b 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -1108,22 +1108,6 @@ void FlatSourceDomain::convert_external_sources() } } } // End loop over external sources - - // Divide the fixed source term by sigma t (to save time when applying each - // iteration) - /* - #pragma omp parallel for - for (int64_t sr = 0; sr < n_source_regions(); sr++) { - int material = source_regions_.material(sr); - if (material == MATERIAL_VOID) { - continue; - } - for (int g = 0; g < negroups_; g++) { - double sigma_t = sigma_t_[material * negroups_ + g]; - source_regions_.external_source(sr, g) /= sigma_t; - } - } - */ } void FlatSourceDomain::flux_swap() @@ -1510,37 +1494,41 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( // 3. Determine if there are any meshes assigned to this handle.mesh() = mesh_idx; - // 4. Determine if there are any point sources, and apply them. - // Point sources are specific to the SRK. - auto it4 = external_point_source_map_.find(sr_key); - if (it4 != external_point_source_map_.end()) { - handle.external_source_present() = 1; - const vector& point_sources = it4->second; - for (int src_idx : point_sources) { - apply_external_source_to_source_region(src_idx, handle); - } - } - - // 5. Determine if there are any volumetric sources, and apply them. + // 4. Determine if there are any volumetric sources, and apply them. // Volumetric sources are specifc only to the base SR idx. auto it5 = external_volumetric_source_map_.find(sr); if (it5 != external_volumetric_source_map_.end()) { - handle.external_source_present() = 1; const vector& vol_sources = it5->second; for (int src_idx : vol_sources) { apply_external_source_to_source_region(src_idx, handle); } } - // 6. Initialize scalar flux based off of 1.0 or 0.0 depending on k vs. fixed, - // and if external source is present or not. - for (int g = 0; g < negroups_; g++) { - if (settings::run_mode == RunMode::FIXED_SOURCE && - handle.external_source_present()) { - handle.scalar_flux_old(g) = 1.0; + // 5. Determine if there are any point sources, and apply them. + // Point sources are specific to the SRK. + auto it4 = external_point_source_map_.find(sr_key); + if (it4 != external_point_source_map_.end()) { + const vector& point_sources = it4->second; + for (int src_idx : point_sources) { + apply_external_source_to_source_region(src_idx, handle); } } + // 6. Initialize scalar flux based off of 1.0 or 0.0 depending on k vs. fixed, + // and if external source is present or not. + // TODO: This changes the solutions slightly, despite probably helping + // things to converge more rapidly? + // That said, perhaps this is not the best idea. In void regions, the + // assumption that there is a ton of flux present is maybe not the best + // idea? Maybe better to just assume everything starts at zero. + + // for (int g = 0; g < negroups_; g++) { + // if (settings::run_mode == RunMode::FIXED_SOURCE && + // handle.external_source_present()) { + // handle.scalar_flux_old(g) = 1.0; + // } + // } + // 7. Divide external source by sigma_t if (material != C_NONE && settings::run_mode == RunMode::FIXED_SOURCE) { for (int g = 0; g < negroups_; g++) { diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 3068be2cb1c..cefa42a4998 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -111,9 +111,8 @@ void openmc_run_random_ray() RandomRaySimulation adjoint_sim; // Initialize adjoint fixed sources, if present - adjoint_sim.prepare_fixed_sources_adjoint(forward_flux, - forward_source_regions, - forward_source_region_map); + adjoint_sim.prepare_fixed_sources_adjoint( + forward_flux, forward_source_regions, forward_source_region_map); // Transpose scattering matrix adjoint_sim.domain()->transpose_scattering_matrix(); @@ -466,11 +465,9 @@ void RandomRaySimulation::simulate() simulation::time_transport.stop(); - // If using mesh subdivision, add any newly discovered source regions - // to the main source region container. - // if (RandomRay::mesh_subdivision_enabled_) { + // Add any newly discovered source regions to the main source region + // container. domain_->finalize_discovered_source_regions(); - //} // Normalize scalar flux and update volumes domain_->normalize_scalar_flux_and_volumes( From 81d3516a5b466a9638c6ca020efbdde28e541e90 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Fri, 19 Sep 2025 11:36:08 -0500 Subject: [PATCH 07/19] fixed bug with adjoints, caused by update_neutron_source not being called in first batch of adjoint solve --- src/random_ray/flat_source_domain.cpp | 2 +- src/random_ray/random_ray_simulation.cpp | 17 ++++++++--------- 2 files changed, 9 insertions(+), 10 deletions(-) diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index ee566bb676b..c84051dfc13 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -636,7 +636,7 @@ void FlatSourceDomain::random_ray_tally() "random ray mode."); break; } - + fmt::print("Score: {} Flux: {} Volume: {}\n", score, flux, volume); // Apply score to the appropriate tally bin Tally& tally {*model::tallies[task.tally_idx]}; #pragma omp atomic diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index cefa42a4998..46065ba52eb 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -436,15 +436,14 @@ void RandomRaySimulation::simulate() // Reset total starting particle weight used for normalizing tallies simulation::total_weight = 1.0; - if (simulation::current_batch > 1) { - // Update source term (scattering + fission) - domain_->update_all_neutron_sources(); - // Reset scalar fluxes, iteration volume tallies, and region hit flags - // to zero - domain_->batch_reset(); - } - - // At the beginning of the simulation, if mesh subvivision is in use, we + // Update source term (scattering + fission) + domain_->update_all_neutron_sources(); + + // Reset scalar fluxes, iteration volume tallies, and region hit flags + // to zero + domain_->batch_reset(); + + // At the beginning of the simulation, if mesh subdivision is in use, we // need to swap the main source region container into the base container, // as the main source region container will be used to hold the true // subdivided source regions. The base container will therefore only From 615e319f79395ebb86e0069330acbc5e687b17a8 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Fri, 19 Sep 2025 12:00:20 -0500 Subject: [PATCH 08/19] removed debug prints --- src/random_ray/flat_source_domain.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index c84051dfc13..23e0f18d040 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -636,7 +636,7 @@ void FlatSourceDomain::random_ray_tally() "random ray mode."); break; } - fmt::print("Score: {} Flux: {} Volume: {}\n", score, flux, volume); + // Apply score to the appropriate tally bin Tally& tally {*model::tallies[task.tally_idx]}; #pragma omp atomic From e3aa71f64f4dcf7636d931772292a5d7556e81f3 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Fri, 19 Sep 2025 13:46:59 -0500 Subject: [PATCH 09/19] fixed final bug. Issue was with adjoint and meshes, the mesh_map was not getting copied to the adjoint Simulation domain. --- include/openmc/random_ray/random_ray_simulation.h | 3 ++- src/random_ray/flat_source_domain.cpp | 2 +- src/random_ray/random_ray_simulation.cpp | 10 +++++++--- .../random_ray_point_source_locator/results_true.dat | 8 ++++---- 4 files changed, 14 insertions(+), 9 deletions(-) diff --git a/include/openmc/random_ray/random_ray_simulation.h b/include/openmc/random_ray/random_ray_simulation.h index 3d1e6158a0d..e18d34342c6 100644 --- a/include/openmc/random_ray/random_ray_simulation.h +++ b/include/openmc/random_ray/random_ray_simulation.h @@ -24,7 +24,8 @@ class RandomRaySimulation { void prepare_fixed_sources_adjoint(vector& forward_flux, SourceRegionContainer& forward_source_regions, std::unordered_map& - forward_source_region_map); + forward_source_region_map, + std::unordered_map& mesh_map); void simulate(); void output_simulation_results() const; void instability_check( diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 23e0f18d040..54c9f3d2b7d 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -1230,7 +1230,6 @@ void FlatSourceDomain::set_adjoint_sources(const vector& forward_flux) source_regions_.external_source_present(sr) = 0; } } - // Divide the fixed source term by sigma t (to save time when applying each // iteration) #pragma omp parallel for @@ -1385,6 +1384,7 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( return source_regions_.get_source_region_handle(sr); } + // Case 2: Check if the source region key is present in the temporary (thread // safe) map. This is a common occurrence in the first power iteration when // the source region has already been visited already by some other ray. We diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 46065ba52eb..35624590272 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -52,6 +52,7 @@ void openmc_run_random_ray() SourceRegionContainer forward_source_regions; std::unordered_map forward_source_region_map; + std::unordered_map mesh_map; { // Initialize Random Ray Simulation Object @@ -83,6 +84,7 @@ void openmc_run_random_ray() forward_source_regions = sim.domain()->source_regions_; forward_source_region_map = sim.domain()->source_region_map_; + mesh_map = sim.domain()->mesh_map_; // Finalize OpenMC openmc_simulation_finalize(); @@ -112,7 +114,7 @@ void openmc_run_random_ray() // Initialize adjoint fixed sources, if present adjoint_sim.prepare_fixed_sources_adjoint( - forward_flux, forward_source_regions, forward_source_region_map); + forward_flux, forward_source_regions, forward_source_region_map, mesh_map); // Transpose scattering matrix adjoint_sim.domain()->transpose_scattering_matrix(); @@ -410,12 +412,13 @@ void RandomRaySimulation::apply_fixed_sources_and_mesh_domains() void RandomRaySimulation::prepare_fixed_sources_adjoint( vector& forward_flux, SourceRegionContainer& forward_source_regions, std::unordered_map& - forward_source_region_map) + forward_source_region_map, std::unordered_map& mesh_map) { domain_->k_eff_ = 1.0; if (settings::run_mode == RunMode::FIXED_SOURCE) { domain_->source_regions_ = forward_source_regions; domain_->source_region_map_ = forward_source_region_map; + domain_->mesh_map_ = mesh_map; domain_->source_regions_.adjoint_reset(); domain_->set_adjoint_sources(forward_flux); } @@ -433,12 +436,13 @@ void RandomRaySimulation::simulate() // TODO: Implement domain decomposition for MPI parallelism if (mpi::master) { + // Reset total starting particle weight used for normalizing tallies simulation::total_weight = 1.0; // Update source term (scattering + fission) domain_->update_all_neutron_sources(); - + // Reset scalar fluxes, iteration volume tallies, and region hit flags // to zero domain_->batch_reset(); diff --git a/tests/regression_tests/random_ray_point_source_locator/results_true.dat b/tests/regression_tests/random_ray_point_source_locator/results_true.dat index 1785dda574a..8c6f358dd31 100644 --- a/tests/regression_tests/random_ray_point_source_locator/results_true.dat +++ b/tests/regression_tests/random_ray_point_source_locator/results_true.dat @@ -1,9 +1,9 @@ tally 1: -2.633900E+00 -2.948207E+00 +2.633923E+00 +2.948228E+00 tally 2: -1.440463E-01 -3.294032E-03 +1.440456E-01 +3.293984E-03 tally 3: 9.425207E-03 1.089748E-05 From 1125eb850a596fbd92afeff6675cd73334275c36 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Fri, 19 Sep 2025 14:17:04 -0500 Subject: [PATCH 10/19] removal of forward vs. adjoint top level simulation objects. Now only one is reused --- .../openmc/random_ray/flat_source_domain.h | 2 +- .../openmc/random_ray/random_ray_simulation.h | 6 +- src/random_ray/flat_source_domain.cpp | 7 +- src/random_ray/random_ray_simulation.cpp | 128 ++++++++---------- src/random_ray/source_region.cpp | 7 +- 5 files changed, 64 insertions(+), 86 deletions(-) diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index 1a167bb64e6..e04568e42be 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -42,7 +42,7 @@ class FlatSourceDomain { void output_to_vtk() const; void convert_external_sources(); void count_external_source_regions(); - void set_adjoint_sources(const vector& forward_flux); + void set_adjoint_sources(); void flux_swap(); virtual double evaluate_flux_at_point(Position r, int64_t sr, int g) const; double compute_fixed_source_normalization_factor() const; diff --git a/include/openmc/random_ray/random_ray_simulation.h b/include/openmc/random_ray/random_ray_simulation.h index e18d34342c6..3dec48bf266 100644 --- a/include/openmc/random_ray/random_ray_simulation.h +++ b/include/openmc/random_ray/random_ray_simulation.h @@ -21,11 +21,7 @@ class RandomRaySimulation { // Methods void compute_segment_correction_factors(); void apply_fixed_sources_and_mesh_domains(); - void prepare_fixed_sources_adjoint(vector& forward_flux, - SourceRegionContainer& forward_source_regions, - std::unordered_map& - forward_source_region_map, - std::unordered_map& mesh_map); + void prepare_fixed_sources_adjoint(); void simulate(); void output_simulation_results() const; void instability_check( diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 54c9f3d2b7d..59c6bddb4da 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -1171,7 +1171,7 @@ void FlatSourceDomain::flatten_xs() } } -void FlatSourceDomain::set_adjoint_sources(const vector& forward_flux) +void FlatSourceDomain::set_adjoint_sources() { // Set the adjoint external source to 1/forward_flux. If the forward flux is // negative, zero, or extremely close to zero, set the adjoint source to zero, @@ -1185,7 +1185,7 @@ void FlatSourceDomain::set_adjoint_sources(const vector& forward_flux) double max_flux = 0.0; #pragma omp parallel for reduction(max : max_flux) for (int64_t se = 0; se < n_source_elements(); se++) { - double flux = forward_flux[se]; + double flux = source_regions_.scalar_flux_final(se); if (flux > max_flux) { max_flux = flux; } @@ -1195,7 +1195,7 @@ void FlatSourceDomain::set_adjoint_sources(const vector& forward_flux) #pragma omp parallel for for (int64_t sr = 0; sr < n_source_regions(); sr++) { for (int g = 0; g < negroups_; g++) { - double flux = forward_flux[sr * negroups_ + g]; + double flux = source_regions_.scalar_flux_final(sr,g); if (flux <= ZERO_FLUX_CUTOFF * max_flux) { source_regions_.external_source(sr, g) = 0.0; } else { @@ -1204,6 +1204,7 @@ void FlatSourceDomain::set_adjoint_sources(const vector& forward_flux) if (flux > 0.0) { source_regions_.external_source_present(sr) = 1; } + source_regions_.scalar_flux_final(sr,g) = 0.0; } } diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 35624590272..d475b2593ea 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -47,96 +47,82 @@ void openmc_run_random_ray() if (mpi::master) validate_random_ray_inputs(); - // Declare forward flux so that it can be saved for later adjoint simulation - vector forward_flux; - SourceRegionContainer forward_source_regions; - std::unordered_map - forward_source_region_map; - std::unordered_map mesh_map; + // Initialize Random Ray Simulation Object + RandomRaySimulation sim; - { - // Initialize Random Ray Simulation Object - RandomRaySimulation sim; + // Initialize fixed sources, if present + sim.apply_fixed_sources_and_mesh_domains(); - // Initialize fixed sources, if present - sim.apply_fixed_sources_and_mesh_domains(); + // Begin main simulation timer + simulation::time_total.start(); - // Begin main simulation timer - simulation::time_total.start(); + // Execute random ray simulation + sim.simulate(); - // Execute random ray simulation - sim.simulate(); + // End main simulation timer + simulation::time_total.stop(); - // End main simulation timer - simulation::time_total.stop(); - - // Normalize and save the final forward flux - sim.domain()->serialize_final_fluxes(forward_flux); - - double source_normalization_factor = - sim.domain()->compute_fixed_source_normalization_factor() / - (settings::n_batches - settings::n_inactive); + // Normalize and save the final forward flux + double source_normalization_factor = + sim.domain()->compute_fixed_source_normalization_factor() / + (settings::n_batches - settings::n_inactive); #pragma omp parallel for - for (uint64_t i = 0; i < forward_flux.size(); i++) { - forward_flux[i] *= source_normalization_factor; - } - - forward_source_regions = sim.domain()->source_regions_; - forward_source_region_map = sim.domain()->source_region_map_; - mesh_map = sim.domain()->mesh_map_; + for (uint64_t se = 0; se < sim.domain()->n_source_elements(); se++) { + sim.domain()->source_regions_.scalar_flux_final(se) *= + source_normalization_factor; + } - // Finalize OpenMC - openmc_simulation_finalize(); + // Finalize OpenMC + openmc_simulation_finalize(); - // Output all simulation results - sim.output_simulation_results(); - } + // Output all simulation results + sim.output_simulation_results(); ////////////////////////////////////////////////////////// // Run adjoint simulation (if enabled) ////////////////////////////////////////////////////////// - if (adjoint_needed) { - reset_timers(); + if (!adjoint_needed) { + return; + } - // Configure the domain for adjoint simulation - FlatSourceDomain::adjoint_ = true; + reset_timers(); - if (mpi::master) - header("ADJOINT FLUX SOLVE", 3); + // Configure the domain for adjoint simulation + FlatSourceDomain::adjoint_ = true; - // Initialize OpenMC general data structures - openmc_simulation_init(); + if (mpi::master) + header("ADJOINT FLUX SOLVE", 3); + + // Initialize OpenMC general data structures + openmc_simulation_init(); - // Initialize Random Ray Simulation Object - RandomRaySimulation adjoint_sim; + sim.domain()->k_eff_ = 1.0; - // Initialize adjoint fixed sources, if present - adjoint_sim.prepare_fixed_sources_adjoint( - forward_flux, forward_source_regions, forward_source_region_map, mesh_map); + // Initialize adjoint fixed sources, if present + sim.prepare_fixed_sources_adjoint(); - // Transpose scattering matrix - adjoint_sim.domain()->transpose_scattering_matrix(); + // Transpose scattering matrix + sim.domain()->transpose_scattering_matrix(); - // Swap nu_sigma_f and chi - adjoint_sim.domain()->nu_sigma_f_.swap(adjoint_sim.domain()->chi_); + // Swap nu_sigma_f and chi + sim.domain()->nu_sigma_f_.swap(sim.domain()->chi_); - // Begin main simulation timer - simulation::time_total.start(); + // Begin main simulation timer + simulation::time_total.start(); - // Execute random ray simulation - adjoint_sim.simulate(); + // Execute random ray simulation + sim.simulate(); - // End main simulation timer - simulation::time_total.stop(); + // End main simulation timer + simulation::time_total.stop(); - // Finalize OpenMC - openmc_simulation_finalize(); + // Finalize OpenMC + openmc_simulation_finalize(); - // Output all simulation results - adjoint_sim.output_simulation_results(); - } + // Output all simulation results + sim.output_simulation_results(); } // Enforces restrictions on inputs in random ray mode. While there are @@ -409,18 +395,11 @@ void RandomRaySimulation::apply_fixed_sources_and_mesh_domains() } } -void RandomRaySimulation::prepare_fixed_sources_adjoint( - vector& forward_flux, SourceRegionContainer& forward_source_regions, - std::unordered_map& - forward_source_region_map, std::unordered_map& mesh_map) +void RandomRaySimulation::prepare_fixed_sources_adjoint() { - domain_->k_eff_ = 1.0; + domain_->source_regions_.adjoint_reset(); if (settings::run_mode == RunMode::FIXED_SOURCE) { - domain_->source_regions_ = forward_source_regions; - domain_->source_region_map_ = forward_source_region_map; - domain_->mesh_map_ = mesh_map; - domain_->source_regions_.adjoint_reset(); - domain_->set_adjoint_sources(forward_flux); + domain_->set_adjoint_sources(); } } @@ -436,7 +415,6 @@ void RandomRaySimulation::simulate() // TODO: Implement domain decomposition for MPI parallelism if (mpi::master) { - // Reset total starting particle weight used for normalizing tallies simulation::total_weight = 1.0; diff --git a/src/random_ray/source_region.cpp b/src/random_ray/source_region.cpp index e0eb090900b..f2a36f1ec26 100644 --- a/src/random_ray/source_region.cpp +++ b/src/random_ray/source_region.cpp @@ -259,9 +259,12 @@ void SourceRegionContainer::adjoint_reset() MomentMatrix {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}); std::fill(mom_matrix_t_.begin(), mom_matrix_t_.end(), MomentMatrix {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}); - std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 0.0); + if (settings::run_mode == RunMode::FIXED_SOURCE) { + std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 0.0); + } else { + std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 1.0); + } std::fill(scalar_flux_new_.begin(), scalar_flux_new_.end(), 0.0); - std::fill(scalar_flux_final_.begin(), scalar_flux_final_.end(), 0.0); std::fill(source_.begin(), source_.end(), 0.0f); std::fill(external_source_.begin(), external_source_.end(), 0.0f); std::fill(source_gradients_.begin(), source_gradients_.end(), From 9862252ccb8bd48b1a2cac7d023742dcd6e3f10f Mon Sep 17 00:00:00 2001 From: John Tramm Date: Fri, 19 Sep 2025 15:04:31 -0500 Subject: [PATCH 11/19] Moving all source region indexing logic into functions to avoid reproducing same logic everywhere --- .../openmc/random_ray/flat_source_domain.h | 9 +- src/random_ray/flat_source_domain.cpp | 125 ++++++++++-------- src/random_ray/random_ray.cpp | 41 ++---- 3 files changed, 83 insertions(+), 92 deletions(-) diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index e04568e42be..7fdfe4cf78b 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -56,7 +56,7 @@ class FlatSourceDomain { void apply_mesh_to_cell_and_children(int32_t i_cell, int32_t mesh_idx, int32_t target_material_id, bool is_target_void); SourceRegionHandle get_subdivided_source_region_handle( - int64_t sr, int mesh_bin, Position r, double dist, Direction u); + SourceRegionKey sr_key, Position r, double dist, Direction u); void finalize_discovered_source_regions(); void apply_transport_stabilization(); int64_t n_source_regions() const @@ -67,6 +67,10 @@ class FlatSourceDomain { { return source_regions_.n_source_regions() * negroups_; } + int64_t lookup_base_source_region_idx(const GeometryState& p) const; + SourceRegionKey lookup_source_region_key(const GeometryState& p) const; + int64_t lookup_mesh_bin(int64_t sr, Position r) const; + int lookup_mesh_idx(int64_t sr) const; //---------------------------------------------------------------------------- // Static Data members @@ -133,8 +137,7 @@ class FlatSourceDomain { // Map that relates a base source region index to the external source index. // This map is used to check if there are any volumetric sources within a // subdivided source region at the time it is discovered. - std::unordered_map> - external_volumetric_source_map_; + std::unordered_map> external_volumetric_source_map_; // Map that relates a base source region index to a mesh index. This map // is used to check which subdivision mesh is present in a source region. diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 59c6bddb4da..ade9070e5dd 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -145,7 +145,7 @@ void FlatSourceDomain::update_all_neutron_sources() { simulation::time_update_src.start(); - #pragma omp parallel for +#pragma omp parallel for for (int64_t sr = 0; sr < n_source_regions(); sr++) { SourceRegionHandle srh = source_regions_.get_source_region_handle(sr); update_single_neutron_source(srh); @@ -636,7 +636,6 @@ void FlatSourceDomain::random_ray_tally() "random ray mode."); break; } - // Apply score to the appropriate tally bin Tally& tally {*model::tallies[task.tally_idx]}; #pragma omp atomic @@ -710,21 +709,21 @@ void FlatSourceDomain::output_to_vtk() const print_plot(); // Outer loop over plots - for (int p = 0; p < model::plots.size(); p++) { + for (int plt = 0; plt < model::plots.size(); plt++) { // Get handle to OpenMC plot object and extract params - Plot* openmc_plot = dynamic_cast(model::plots[p].get()); + Plot* openmc_plot = dynamic_cast(model::plots[plt].get()); // Random ray plots only support voxel plots if (!openmc_plot) { warning(fmt::format("Plot {} is invalid plot type -- only voxel plotting " "is allowed in random ray mode.", - p)); + plt)); continue; } else if (openmc_plot->type_ != Plot::PlotType::voxel) { warning(fmt::format("Plot {} is invalid plot type -- only voxel plotting " "is allowed in random ray mode.", - p)); + plt)); continue; } @@ -778,25 +777,11 @@ void FlatSourceDomain::output_to_vtk() const continue; } - int i_cell = p.lowest_coord().cell(); - int64_t sr = source_region_offsets_[i_cell] + p.cell_instance(); - int mesh_idx = C_NONE; - auto mesh_it = mesh_map_.find(sr); - if (mesh_it != mesh_map_.end()) { - mesh_idx = mesh_it->second; - } - int mesh_bin; - if (mesh_idx == C_NONE) { - mesh_bin = 0; - } else { - mesh_bin = model::meshes[mesh_idx]->get_bin(p.r()); - } - SourceRegionKey sr_key {sr, mesh_bin}; + SourceRegionKey sr_key = lookup_source_region_key(p); + int64_t sr = -1; auto it = source_region_map_.find(sr_key); if (it != source_region_map_.end()) { sr = it->second; - } else { - sr = -1; } voxel_indices[z * Ny * Nx + y * Nx + x] = sr; @@ -1056,22 +1041,8 @@ void FlatSourceDomain::convert_external_sources() "point source at {}", sp->r())); } - int i_cell = gs.lowest_coord().cell(); - int64_t sr = source_region_offsets_[i_cell] + gs.cell_instance(); - - // If mesh subdivision is enabled, we need to determine which subdivided - // mesh bin the point source coordinate is in as well - int mesh_idx = C_NONE; - auto mesh_it = mesh_map_.find(sr); - if (mesh_it != mesh_map_.end()) { - mesh_idx = mesh_it->second; - } - int mesh_bin; - if (mesh_idx == C_NONE) { - mesh_bin = 0; - } else { - mesh_bin = model::meshes[mesh_idx]->get_bin(gs.r()); - } + SourceRegionKey key = lookup_source_region_key(gs); + // With the source region and mesh bin known, we can use the // accompanying SourceRegionKey as a key into a map that stores the // corresponding external source index for the point source. Notably, we @@ -1080,7 +1051,6 @@ void FlatSourceDomain::convert_external_sources() // discovered & initilized yet. When discovered, they will read from the // external_source_map to determine if there are any external source // terms that should be applied. - SourceRegionKey key {sr, mesh_bin}; external_point_source_map_[key].push_back(es); } else { @@ -1195,7 +1165,7 @@ void FlatSourceDomain::set_adjoint_sources() #pragma omp parallel for for (int64_t sr = 0; sr < n_source_regions(); sr++) { for (int g = 0; g < negroups_; g++) { - double flux = source_regions_.scalar_flux_final(sr,g); + double flux = source_regions_.scalar_flux_final(sr, g); if (flux <= ZERO_FLUX_CUTOFF * max_flux) { source_regions_.external_source(sr, g) = 0.0; } else { @@ -1204,7 +1174,7 @@ void FlatSourceDomain::set_adjoint_sources() if (flux > 0.0) { source_regions_.external_source_present(sr) = 1; } - source_regions_.scalar_flux_final(sr,g) = 0.0; + source_regions_.scalar_flux_final(sr, g) = 0.0; } } @@ -1369,10 +1339,8 @@ void FlatSourceDomain::apply_meshes() } SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( - int64_t sr, int mesh_bin, Position r, double dist, Direction u) + SourceRegionKey sr_key, Position r, double dist, Direction u) { - SourceRegionKey sr_key {sr, mesh_bin}; - // Case 1: Check if the source region key is already present in the permanent // map. This is the most common condition, as any source region visited in a // previous power iteration will already be present in the permanent map. If @@ -1385,7 +1353,6 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( return source_regions_.get_source_region_handle(sr); } - // Case 2: Check if the source region key is present in the temporary (thread // safe) map. This is a common occurrence in the first power iteration when // the source region has already been visited already by some other ray. We @@ -1435,9 +1402,8 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( gs.r() = r + TINY_BIT * u; gs.u() = {1.0, 0.0, 0.0}; exhaustive_find_cell(gs); - int gs_i_cell = gs.lowest_coord().cell(); - int64_t sr_found = source_region_offsets_[gs_i_cell] + gs.cell_instance(); - if (sr_found != sr) { + int64_t sr_found = lookup_base_source_region_idx(gs); + if (sr_found != sr_key.base_source_region_id) { discovered_source_regions_.unlock(sr_key); SourceRegionHandle handle; handle.is_numerical_fp_artifact_ = true; @@ -1445,13 +1411,9 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( } // Sanity check on mesh bin - int mesh_idx = C_NONE; - auto mesh_it = mesh_map_.find(sr); - if (mesh_it != mesh_map_.end()) { - mesh_idx = mesh_it->second; - } + int mesh_idx = lookup_mesh_idx(sr_key.base_source_region_id); if (mesh_idx == C_NONE) { - if (mesh_bin != 0) { + if (sr_key.mesh_bin != 0) { discovered_source_regions_.unlock(sr_key); SourceRegionHandle handle; handle.is_numerical_fp_artifact_ = true; @@ -1460,7 +1422,7 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( } else { Mesh* mesh = model::meshes[mesh_idx].get(); int bin_found = mesh->get_bin(r + TINY_BIT * u); - if (bin_found != mesh_bin) { + if (bin_found != sr_key.mesh_bin) { discovered_source_regions_.unlock(sr_key); SourceRegionHandle handle; handle.is_numerical_fp_artifact_ = true; @@ -1488,6 +1450,7 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( discovered_source_regions_.unlock(sr_key); // 2. Determine the material + int gs_i_cell = gs.lowest_coord().cell(); Cell& cell = *model::cells[gs_i_cell]; int material = cell.material(gs.cell_instance()); handle.material() = material; @@ -1497,7 +1460,7 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( // 4. Determine if there are any volumetric sources, and apply them. // Volumetric sources are specifc only to the base SR idx. - auto it5 = external_volumetric_source_map_.find(sr); + auto it5 = external_volumetric_source_map_.find(sr_key.base_source_region_id); if (it5 != external_volumetric_source_map_.end()) { const vector& vol_sources = it5->second; for (int src_idx : vol_sources) { @@ -1517,7 +1480,7 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( // 6. Initialize scalar flux based off of 1.0 or 0.0 depending on k vs. fixed, // and if external source is present or not. - // TODO: This changes the solutions slightly, despite probably helping + // TODO: This changes the solutions slightly, despite probably helping // things to converge more rapidly? // That said, perhaps this is not the best idea. In void regions, the // assumption that there is a ton of flux present is maybe not the best @@ -1629,4 +1592,52 @@ void FlatSourceDomain::apply_transport_stabilization() } } +// Determines the base source region index (i.e., a material filled cell +// instance) that corresponds to a particular location in the geometry. Requires +// that the "gs" object passed in has already been initialized and has called +// find_cell etc. +int64_t FlatSourceDomain::lookup_base_source_region_idx( + const GeometryState& gs) const +{ + int i_cell = gs.lowest_coord().cell(); + int64_t sr = source_region_offsets_[i_cell] + gs.cell_instance(); + return sr; +} + +// Determines the index of the mesh (if any) that has been applied +// to a particular base source region index. +int FlatSourceDomain::lookup_mesh_idx(int64_t sr) const +{ + int mesh_idx = C_NONE; + auto mesh_it = mesh_map_.find(sr); + if (mesh_it != mesh_map_.end()) { + mesh_idx = mesh_it->second; + } + return mesh_idx; +} + +// Determines the source region key that corresponds to a particular location in +// the geometry. This takes into account both the base source region index as +// well as the mesh bin if a mesh is applied to this source region for +// subdivision. +SourceRegionKey FlatSourceDomain::lookup_source_region_key( + const GeometryState& gs) const +{ + int64_t sr = lookup_base_source_region_idx(gs); + int64_t mesh_bin = lookup_mesh_bin(sr, gs.r()); + return SourceRegionKey {sr, mesh_bin}; +} + +// Determines the mesh bin that corresponds to a particular base source region +// index and position. +int64_t FlatSourceDomain::lookup_mesh_bin(int64_t sr, Position r) const +{ + int mesh_idx = lookup_mesh_idx(sr); + int mesh_bin = 0; + if (mesh_idx != C_NONE) { + mesh_bin = model::meshes[mesh_idx]->get_bin(r); + } + return mesh_bin; +} + } // namespace openmc diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index c0e42a88547..3db148da057 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -335,19 +335,12 @@ void RandomRay::event_advance_ray() void RandomRay::attenuate_flux(double distance, bool is_active, double offset) { - // Determine source region index etc. - int i_cell = lowest_coord().cell(); - - // The base source region is the spatial region index - int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance(); + // Lookup base source region index + int64_t sr = domain_->lookup_base_source_region_idx(*this); // Perform ray tracing across mesh // Determine the mesh index for the base source region, if any - int mesh_idx = C_NONE; - auto it = domain_->mesh_map_.find(sr); - if (it != domain_->mesh_map_.end()) { - mesh_idx = it->second; - } + int mesh_idx = domain_->lookup_mesh_idx(sr); if (mesh_idx == C_NONE) { // If there's no mesh being applied to this cell, then @@ -386,9 +379,10 @@ void RandomRay::attenuate_flux(double distance, bool is_active, double offset) void RandomRay::attenuate_flux_inner( double distance, bool is_active, int64_t sr, int mesh_bin, Position r) { + SourceRegionKey sr_key {sr, mesh_bin}; SourceRegionHandle srh; srh = domain_->get_subdivided_source_region_handle( - sr, mesh_bin, r, distance, u()); + sr_key, r, distance, u()); if (srh.is_numerical_fp_artifact_) { return; } @@ -805,29 +799,12 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain) cell_born() = lowest_coord().cell(); } + SourceRegionKey sr_key = domain_->lookup_source_region_key(*this); + SourceRegionHandle srh = + domain_->get_subdivided_source_region_handle(sr_key, r(), 0.0, u()); + // Initialize ray's starting angular flux to starting location's isotropic // source - int i_cell = lowest_coord().cell(); - int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance(); - - SourceRegionHandle srh; - - // Determine if there are any meshes assigned to this - int mesh_idx = C_NONE; - auto it = domain_->mesh_map_.find(sr); - if (it != domain_->mesh_map_.end()) { - mesh_idx = it->second; - } - int mesh_bin; - if (mesh_idx == C_NONE) { - mesh_bin = 0; - } else { - Mesh* mesh = model::meshes[mesh_idx].get(); - mesh_bin = mesh->get_bin(r()); - } - srh = - domain_->get_subdivided_source_region_handle(sr, mesh_bin, r(), 0.0, u()); - if (!srh.is_numerical_fp_artifact_) { for (int g = 0; g < negroups_; g++) { angular_flux_[g] = srh.source(g); From fb037a210acb950bde3811f8136403e3642fff18 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Fri, 19 Sep 2025 15:20:22 -0500 Subject: [PATCH 12/19] cleanup of comments on source region handle finder --- include/openmc/random_ray/moment_matrix.h | 12 ++-- src/random_ray/flat_source_domain.cpp | 80 ++++++++++------------- 2 files changed, 41 insertions(+), 51 deletions(-) diff --git a/include/openmc/random_ray/moment_matrix.h b/include/openmc/random_ray/moment_matrix.h index 508b562a478..c95bb2c1286 100644 --- a/include/openmc/random_ray/moment_matrix.h +++ b/include/openmc/random_ray/moment_matrix.h @@ -26,12 +26,12 @@ class MomentMatrix { public: //---------------------------------------------------------------------------- // Public data members - double a {0}; - double b {0}; - double c {0}; - double d {0}; - double e {0}; - double f {0}; + double a; + double b; + double c; + double d; + double e; + double f; //---------------------------------------------------------------------------- // Constructors diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index ade9070e5dd..8c1de5dce40 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -1434,77 +1434,67 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( // condition only occurs the first time the source region is discovered // (typically in the first power iteration). In this case, we need to handle // creation of the new source region and its storage into the parallel map. - // The new source region is created by copying the base source region, so as - // to inherit material, external source, and some flux properties etc. We - // also pass the base source region id to allow the new source region to - // know which base source region it is derived from. + // Additionally, we need to determine the source region's material, initialize + // the starting scalar flux guess, and apply any known external sources. - // 0. Call the basic constructor for the source region + // Call the basic constructor for the source region and store in the parallel + // map. bool is_linear = RandomRay::source_shape_ != RandomRaySourceShape::FLAT; SourceRegion* sr_ptr = discovered_source_regions_.emplace(sr_key, {negroups_, is_linear}); SourceRegionHandle handle {*sr_ptr}; - // 1. Lock the SR (before writing anything to it) + // Lock the SR (before writing anything to it), and release the lock + // on the parallel map object. Locking the handle must happen first, + // so as to ensure another thread does not visit this source region + // before it is fully initialized. handle.lock(); discovered_source_regions_.unlock(sr_key); - // 2. Determine the material + // Determine the material int gs_i_cell = gs.lowest_coord().cell(); Cell& cell = *model::cells[gs_i_cell]; int material = cell.material(gs.cell_instance()); handle.material() = material; - // 3. Determine if there are any meshes assigned to this + // Store the mesh index (if any) assigned to this source region handle.mesh() = mesh_idx; - // 4. Determine if there are any volumetric sources, and apply them. - // Volumetric sources are specifc only to the base SR idx. - auto it5 = external_volumetric_source_map_.find(sr_key.base_source_region_id); - if (it5 != external_volumetric_source_map_.end()) { - const vector& vol_sources = it5->second; - for (int src_idx : vol_sources) { - apply_external_source_to_source_region(src_idx, handle); + if (settings::run_mode == RunMode::FIXED_SOURCE) { + // Determine if there are any volumetric sources, and apply them. + // Volumetric sources are specifc only to the base SR idx. + auto it_vol = + external_volumetric_source_map_.find(sr_key.base_source_region_id); + if (it_vol != external_volumetric_source_map_.end()) { + const vector& vol_sources = it_vol->second; + for (int src_idx : vol_sources) { + apply_external_source_to_source_region(src_idx, handle); + } } - } - // 5. Determine if there are any point sources, and apply them. - // Point sources are specific to the SRK. - auto it4 = external_point_source_map_.find(sr_key); - if (it4 != external_point_source_map_.end()) { - const vector& point_sources = it4->second; - for (int src_idx : point_sources) { - apply_external_source_to_source_region(src_idx, handle); + // Determine if there are any point sources, and apply them. + // Point sources are specific to the source region key. + auto it_point = external_point_source_map_.find(sr_key); + if (it_point != external_point_source_map_.end()) { + const vector& point_sources = it_point->second; + for (int src_idx : point_sources) { + apply_external_source_to_source_region(src_idx, handle); + } } - } - // 6. Initialize scalar flux based off of 1.0 or 0.0 depending on k vs. fixed, - // and if external source is present or not. - // TODO: This changes the solutions slightly, despite probably helping - // things to converge more rapidly? - // That said, perhaps this is not the best idea. In void regions, the - // assumption that there is a ton of flux present is maybe not the best - // idea? Maybe better to just assume everything starts at zero. - - // for (int g = 0; g < negroups_; g++) { - // if (settings::run_mode == RunMode::FIXED_SOURCE && - // handle.external_source_present()) { - // handle.scalar_flux_old(g) = 1.0; - // } - // } - - // 7. Divide external source by sigma_t - if (material != C_NONE && settings::run_mode == RunMode::FIXED_SOURCE) { - for (int g = 0; g < negroups_; g++) { - double sigma_t = sigma_t_[material * negroups_ + g]; - handle.external_source(g) /= sigma_t; + // Divide external source term by sigma_t + if (material != C_NONE) { + for (int g = 0; g < negroups_; g++) { + double sigma_t = sigma_t_[material * negroups_ + g]; + handle.external_source(g) /= sigma_t; + } } } // Compute the combined source term update_single_neutron_source(handle); - // 10. Unlock the SR + // Unlock the source region handle.unlock(); return handle; From 5a6f2f2832cd3a5c1810cd381e0a91a34bcd9bf2 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Fri, 19 Sep 2025 16:20:15 -0500 Subject: [PATCH 13/19] ran clang format --- src/random_ray/linear_source_domain.cpp | 3 +-- src/random_ray/random_ray.cpp | 3 +-- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/src/random_ray/linear_source_domain.cpp b/src/random_ray/linear_source_domain.cpp index ea97e61be59..e1ad68e3d87 100644 --- a/src/random_ray/linear_source_domain.cpp +++ b/src/random_ray/linear_source_domain.cpp @@ -84,8 +84,7 @@ void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh) // very small/noisy or have poorly developed spatial moments, so we zero // the source gradients (effectively making this a flat source region // temporarily), so as to improve stability. - if (simulation::current_batch > 10 && - srh.source(g_out) >= 0.0) { + if (simulation::current_batch > 10 && srh.source(g_out) >= 0.0) { srh.source_gradients(g_out) = invM * ((scatter_linear + fission_linear * inverse_k_eff) / sigma_t); } else { diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index 3db148da057..25fe6a7062a 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -381,8 +381,7 @@ void RandomRay::attenuate_flux_inner( { SourceRegionKey sr_key {sr, mesh_bin}; SourceRegionHandle srh; - srh = domain_->get_subdivided_source_region_handle( - sr_key, r, distance, u()); + srh = domain_->get_subdivided_source_region_handle(sr_key, r, distance, u()); if (srh.is_numerical_fp_artifact_) { return; } From f6644512e65d28b4304918edd4feab5d036a7476 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Mon, 22 Sep 2025 14:55:04 -0500 Subject: [PATCH 14/19] removed unneeded distance argument, and corrected bug in locking of parallel map --- .../openmc/random_ray/flat_source_domain.h | 2 +- src/random_ray/flat_source_domain.cpp | 24 ++++++++++--------- src/random_ray/random_ray.cpp | 4 ++-- 3 files changed, 16 insertions(+), 14 deletions(-) diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index 7fdfe4cf78b..d4e80273460 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -56,7 +56,7 @@ class FlatSourceDomain { void apply_mesh_to_cell_and_children(int32_t i_cell, int32_t mesh_idx, int32_t target_material_id, bool is_target_void); SourceRegionHandle get_subdivided_source_region_handle( - SourceRegionKey sr_key, Position r, double dist, Direction u); + SourceRegionKey sr_key, Position r, Direction u); void finalize_discovered_source_regions(); void apply_transport_stabilization(); int64_t n_source_regions() const diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 8c1de5dce40..31f4b02eca6 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -1339,7 +1339,7 @@ void FlatSourceDomain::apply_meshes() } SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( - SourceRegionKey sr_key, Position r, double dist, Direction u) + SourceRegionKey sr_key, Position r, Direction u) { // Case 1: Check if the source region key is already present in the permanent // map. This is the most common condition, as any source region visited in a @@ -1444,13 +1444,6 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( discovered_source_regions_.emplace(sr_key, {negroups_, is_linear}); SourceRegionHandle handle {*sr_ptr}; - // Lock the SR (before writing anything to it), and release the lock - // on the parallel map object. Locking the handle must happen first, - // so as to ensure another thread does not visit this source region - // before it is fully initialized. - handle.lock(); - discovered_source_regions_.unlock(sr_key); - // Determine the material int gs_i_cell = gs.lowest_coord().cell(); Cell& cell = *model::cells[gs_i_cell]; @@ -1493,9 +1486,18 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( // Compute the combined source term update_single_neutron_source(handle); - - // Unlock the source region - handle.unlock(); + + // Unlock the parallel map. Note: we may be tempted to release + // this lock earlier, and then just use the source region's lock to protect + // the flux/source initialization stages above. However, the rest of the code + // only protects updates to the new flux and volume fields, and assumes that the + // source is constant for the duration of transport. Thus, using just the source + // region's lock by itself would result in other threads potentially reading from + // the source before it is computed, as they won't use the lock when only reading + // from the SR's source. It would be expensive to protect those operations, + // whereas generating the SR is only done once, so we just hold the map's bucket + // lock until the source region is fully initialized. + discovered_source_regions_.unlock(sr_key); return handle; } diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index 25fe6a7062a..6b38618b471 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -381,7 +381,7 @@ void RandomRay::attenuate_flux_inner( { SourceRegionKey sr_key {sr, mesh_bin}; SourceRegionHandle srh; - srh = domain_->get_subdivided_source_region_handle(sr_key, r, distance, u()); + srh = domain_->get_subdivided_source_region_handle(sr_key, r, u()); if (srh.is_numerical_fp_artifact_) { return; } @@ -800,7 +800,7 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain) SourceRegionKey sr_key = domain_->lookup_source_region_key(*this); SourceRegionHandle srh = - domain_->get_subdivided_source_region_handle(sr_key, r(), 0.0, u()); + domain_->get_subdivided_source_region_handle(sr_key, r(), u()); // Initialize ray's starting angular flux to starting location's isotropic // source From 14adca47ff02906676fb9121a7c7886dd98a9dba Mon Sep 17 00:00:00 2001 From: John Tramm Date: Mon, 22 Sep 2025 14:56:20 -0500 Subject: [PATCH 15/19] ran clang format --- src/random_ray/flat_source_domain.cpp | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 31f4b02eca6..b8ce1bc9b76 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -1486,17 +1486,17 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( // Compute the combined source term update_single_neutron_source(handle); - + // Unlock the parallel map. Note: we may be tempted to release // this lock earlier, and then just use the source region's lock to protect // the flux/source initialization stages above. However, the rest of the code - // only protects updates to the new flux and volume fields, and assumes that the - // source is constant for the duration of transport. Thus, using just the source - // region's lock by itself would result in other threads potentially reading from - // the source before it is computed, as they won't use the lock when only reading - // from the SR's source. It would be expensive to protect those operations, - // whereas generating the SR is only done once, so we just hold the map's bucket - // lock until the source region is fully initialized. + // only protects updates to the new flux and volume fields, and assumes that + // the source is constant for the duration of transport. Thus, using just the + // source region's lock by itself would result in other threads potentially + // reading from the source before it is computed, as they won't use the lock + // when only reading from the SR's source. It would be expensive to protect + // those operations, whereas generating the SR is only done once, so we just + // hold the map's bucket lock until the source region is fully initialized. discovered_source_regions_.unlock(sr_key); return handle; From a2e2adf20a9a24c816e6a154f28c88d815ea54b5 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Mon, 22 Sep 2025 16:36:03 -0500 Subject: [PATCH 16/19] removed unneded SourceRegion constructor --- include/openmc/random_ray/source_region.h | 1 - src/random_ray/source_region.cpp | 19 ------------------- 2 files changed, 20 deletions(-) diff --git a/include/openmc/random_ray/source_region.h b/include/openmc/random_ray/source_region.h index 5c5b31f392e..0f5a747fff8 100644 --- a/include/openmc/random_ray/source_region.h +++ b/include/openmc/random_ray/source_region.h @@ -308,7 +308,6 @@ class SourceRegion { //---------------------------------------------------------------------------- // Constructors SourceRegion(int negroups, bool is_linear); - SourceRegion(const SourceRegionHandle& handle, int64_t parent_sr); SourceRegion() = default; //---------------------------------------------------------------------------- diff --git a/src/random_ray/source_region.cpp b/src/random_ray/source_region.cpp index f2a36f1ec26..3b06f0ed09f 100644 --- a/src/random_ray/source_region.cpp +++ b/src/random_ray/source_region.cpp @@ -60,25 +60,6 @@ SourceRegion::SourceRegion(int negroups, bool is_linear) } } -SourceRegion::SourceRegion(const SourceRegionHandle& handle, int64_t parent_sr) - : SourceRegion(handle.negroups_, handle.is_linear_) -{ - material_ = handle.material(); - mesh_ = handle.mesh(); - parent_sr_ = parent_sr; - for (int g = 0; g < scalar_flux_new_.size(); g++) { - scalar_flux_old_[g] = handle.scalar_flux_old(g); - source_[g] = handle.source(g); - } - - if (settings::run_mode == RunMode::FIXED_SOURCE) { - external_source_present_ = handle.external_source_present(); - for (int g = 0; g < scalar_flux_new_.size(); g++) { - external_source_[g] = handle.external_source(g); - } - } -} - //============================================================================== // SourceRegionContainer implementation //============================================================================== From 4a02d1657e7503c0872e49b3b36f7df1afde074a Mon Sep 17 00:00:00 2001 From: John Tramm Date: Thu, 25 Sep 2025 09:45:00 -0500 Subject: [PATCH 17/19] added check to convert low XS materials to void materials --- include/openmc/constants.h | 5 +++++ src/random_ray/flat_source_domain.cpp | 23 ++++++++++++++++++++++- 2 files changed, 27 insertions(+), 1 deletion(-) diff --git a/include/openmc/constants.h b/include/openmc/constants.h index ae705607958..597816066f3 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -68,6 +68,11 @@ constexpr double MIN_HITS_PER_BATCH {1.5}; // prevent extremely large adjoint source terms from being generated. constexpr double ZERO_FLUX_CUTOFF {1e-22}; +// The minimum macroscopic cross section value considered non-void for the +// random ray solver. Materials with any group with a cross section below this +// value will be converted to pure void. +constexpr double MINIMUM_MACRO_XS {1e-6}; + // ============================================================================ // MATH AND PHYSICAL CONSTANTS diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index b8ce1bc9b76..1bf27e1eda1 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -1094,13 +1094,23 @@ void FlatSourceDomain::flatten_xs() const int a = 0; n_materials_ = data::mg.macro_xs_.size(); - for (auto& m : data::mg.macro_xs_) { + for (int i = 0; i < n_materials_; i++) { + auto& m = data::mg.macro_xs_[i]; for (int g_out = 0; g_out < negroups_; g_out++) { if (m.exists_in_model) { double sigma_t = m.get_xs(MgxsType::TOTAL, g_out, NULL, NULL, NULL, t, a); sigma_t_.push_back(sigma_t); + if (sigma_t < MINIMUM_MACRO_XS) { + Material* mat = model::materials[i].get(); + warning(fmt::format( + "Material \"{}\" (id: {}) has a group {} total cross section " + "({:.3e}) below the minimum threshold " + "({:.3e}). Material will be treated as pure void.", + mat->name(), mat->id(), g_out, sigma_t, MINIMUM_MACRO_XS)); + } + double nu_sigma_f = m.get_xs(MgxsType::NU_FISSION, g_out, NULL, NULL, NULL, t, a); nu_sigma_f_.push_back(nu_sigma_f); @@ -1448,6 +1458,17 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( int gs_i_cell = gs.lowest_coord().cell(); Cell& cell = *model::cells[gs_i_cell]; int material = cell.material(gs.cell_instance()); + + // If material total XS is extremely low, just set it to void to avoid + // problems with 1/Sigma_t + for (int g = 0; g < negroups_; g++) { + double sigma_t = sigma_t_[material * negroups_ + g]; + if (sigma_t < MINIMUM_MACRO_XS) { + material = MATERIAL_VOID; + break; + } + } + handle.material() = material; // Store the mesh index (if any) assigned to this source region From 423e106fddc302962095f467a47f93841856b3c0 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Fri, 26 Sep 2025 10:42:39 -0500 Subject: [PATCH 18/19] fixed issue with ray using particle material vs. source region material --- src/random_ray/random_ray.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index 6b38618b471..89a91449df7 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -388,7 +388,7 @@ void RandomRay::attenuate_flux_inner( switch (source_shape_) { case RandomRaySourceShape::FLAT: - if (this->material() == MATERIAL_VOID) { + if (srh.material() == MATERIAL_VOID) { attenuate_flux_flat_source_void(srh, distance, is_active, r); } else { attenuate_flux_flat_source(srh, distance, is_active, r); @@ -396,7 +396,7 @@ void RandomRay::attenuate_flux_inner( break; case RandomRaySourceShape::LINEAR: case RandomRaySourceShape::LINEAR_XY: - if (this->material() == MATERIAL_VOID) { + if (srh.material() == MATERIAL_VOID) { attenuate_flux_linear_source_void(srh, distance, is_active, r); } else { attenuate_flux_linear_source(srh, distance, is_active, r); @@ -427,7 +427,7 @@ void RandomRay::attenuate_flux_flat_source( n_event()++; // Get material - int material = this->material(); + int material = srh.material(); // MOC incoming flux attenuation + source contribution/attenuation equation for (int g = 0; g < negroups_; g++) { @@ -478,7 +478,7 @@ void RandomRay::attenuate_flux_flat_source_void( // The number of geometric intersections is counted for reporting purposes n_event()++; - int material = this->material(); + int material = srh.material(); // If ray is in the active phase (not in dead zone), make contributions to // source region bookkeeping @@ -525,7 +525,7 @@ void RandomRay::attenuate_flux_linear_source( // The number of geometric intersections is counted for reporting purposes n_event()++; - int material = this->material(); + int material = srh.material(); Position& centroid = srh.centroid(); Position midpoint = r + u() * (distance / 2.0); From 19d7046afd6c353b63d692fa55cb9034434dfc99 Mon Sep 17 00:00:00 2001 From: John Tramm Date: Fri, 26 Sep 2025 11:11:08 -0500 Subject: [PATCH 19/19] added test for random ray low density flux conversion --- .../random_ray_low_density/__init__.py | 0 .../random_ray_low_density/inputs_true.dat | 244 ++++++++++++++++++ .../random_ray_low_density/results_true.dat | 9 + .../random_ray_low_density/test.py | 60 +++++ 4 files changed, 313 insertions(+) create mode 100644 tests/regression_tests/random_ray_low_density/__init__.py create mode 100644 tests/regression_tests/random_ray_low_density/inputs_true.dat create mode 100644 tests/regression_tests/random_ray_low_density/results_true.dat create mode 100644 tests/regression_tests/random_ray_low_density/test.py diff --git a/tests/regression_tests/random_ray_low_density/__init__.py b/tests/regression_tests/random_ray_low_density/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/random_ray_low_density/inputs_true.dat b/tests/regression_tests/random_ray_low_density/inputs_true.dat new file mode 100644 index 00000000000..c4fd06f421e --- /dev/null +++ b/tests/regression_tests/random_ray_low_density/inputs_true.dat @@ -0,0 +1,244 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + + + + + + 2.5 2.5 2.5 + 12 12 12 + 0.0 0.0 0.0 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 + + + + + + + + + + fixed source + 90 + 10 + 5 + + + 100.0 1.0 + + + universe + 1 + + + multi-group + + 500.0 + 100.0 + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + True + + + + + 1 + + + 2 + + + 3 + + + 3 + flux + tracklength + + + 2 + flux + tracklength + + + 1 + flux + tracklength + + + diff --git a/tests/regression_tests/random_ray_low_density/results_true.dat b/tests/regression_tests/random_ray_low_density/results_true.dat new file mode 100644 index 00000000000..a4b3ee1bcde --- /dev/null +++ b/tests/regression_tests/random_ray_low_density/results_true.dat @@ -0,0 +1,9 @@ +tally 1: +5.973607E-01 +7.155477E-02 +tally 2: +3.206216E-02 +2.063375E-04 +tally 3: +2.096415E-03 +8.804963E-07 diff --git a/tests/regression_tests/random_ray_low_density/test.py b/tests/regression_tests/random_ray_low_density/test.py new file mode 100644 index 00000000000..1b4ffb78183 --- /dev/null +++ b/tests/regression_tests/random_ray_low_density/test.py @@ -0,0 +1,60 @@ +import os + +import numpy as np +import openmc +from openmc.examples import random_ray_three_region_cube + +from tests.testing_harness import TolerantPyAPITestHarness + + +class MGXSTestHarness(TolerantPyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = 'mgxs.h5' + if os.path.exists(f): + os.remove(f) + + +def test_random_ray_low_density(): + model = random_ray_three_region_cube() + + # Rebuild the MGXS library to have a material with very + # low macroscopic cross sections + ebins = [1e-5, 20.0e6] + groups = openmc.mgxs.EnergyGroups(group_edges=ebins) + + void_sigma_a = 4.0e-6 + void_sigma_s = 3.0e-4 + void_mat_data = openmc.XSdata('void', groups) + void_mat_data.order = 0 + void_mat_data.set_total([void_sigma_a + void_sigma_s]) + void_mat_data.set_absorption([void_sigma_a]) + void_mat_data.set_scatter_matrix( + np.rollaxis(np.array([[[void_sigma_s]]]), 0, 3)) + + absorber_sigma_a = 0.75 + absorber_sigma_s = 0.25 + absorber_mat_data = openmc.XSdata('absorber', groups) + absorber_mat_data.order = 0 + absorber_mat_data.set_total([absorber_sigma_a + absorber_sigma_s]) + absorber_mat_data.set_absorption([absorber_sigma_a]) + absorber_mat_data.set_scatter_matrix( + np.rollaxis(np.array([[[absorber_sigma_s]]]), 0, 3)) + + multiplier = 0.0000001 + source_sigma_a = void_sigma_a * multiplier + source_sigma_s = void_sigma_s * multiplier + source_mat_data = openmc.XSdata('source', groups) + source_mat_data.order = 0 + source_mat_data.set_total([source_sigma_a + source_sigma_s]) + source_mat_data.set_absorption([source_sigma_a]) + source_mat_data.set_scatter_matrix( + np.rollaxis(np.array([[[source_sigma_s]]]), 0, 3)) + + mg_cross_sections_file = openmc.MGXSLibrary(groups) + mg_cross_sections_file.add_xsdatas( + [source_mat_data, void_mat_data, absorber_mat_data]) + mg_cross_sections_file.export_to_hdf5() + + harness = MGXSTestHarness('statepoint.10.h5', model) + harness.main()