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/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index 78351fcc5f5..d4e80273460 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -27,8 +27,9 @@ class FlatSourceDomain { //---------------------------------------------------------------------------- // Methods - virtual void update_neutron_source(double k_eff); - double compute_k_eff(double k_eff_old) const; + 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); @@ -41,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; @@ -54,9 +55,8 @@ 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); + SourceRegionKey sr_key, Position r, 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 @@ -86,6 +90,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 @@ -110,14 +115,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. @@ -134,8 +131,17 @@ 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> 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 @@ -147,12 +153,11 @@ class FlatSourceDomain { //---------------------------------------------------------------------------- // Methods void apply_external_source_to_source_region( - Discrete* discrete, double strength_factor, SourceRegionHandle& srh); - void apply_external_source_to_cell_instances(int32_t i_cell, - Discrete* discrete, double strength_factor, 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, 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); 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..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(double k_eff) 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/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/include/openmc/random_ray/random_ray_simulation.h b/include/openmc/random_ray/random_ray_simulation.h index b94e7401b3e..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, - SourceRegionContainer& forward_base_source_regions, - std::unordered_map& - forward_source_region_map); + void prepare_fixed_sources_adjoint(); void simulate(); void output_simulation_results() const; void instability_check( @@ -45,9 +41,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/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/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 40923883088..1bf27e1eda1 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -53,24 +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 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_) { @@ -118,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(double k_eff) +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]; @@ -154,18 +126,30 @@ void FlatSourceDomain::update_neutron_source(double k_eff) 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(); } @@ -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 @@ -652,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 @@ -726,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; } @@ -794,23 +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(); - if (RandomRay::mesh_subdivision_enabled_) { - int mesh_idx = base_source_regions_.mesh(sr); - 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; - } + 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; } voxel_indices[z * Ny * Nx + y * Nx + x] = sr; @@ -967,13 +938,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 +956,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 +974,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 +988,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); } } } @@ -1070,36 +1041,17 @@ 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 (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; - } 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); - } + 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 + // 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. + external_point_source_map_[key].push_back(es); } else { // If not a point source, then use the volumetric domain constraints to @@ -1107,42 +1059,25 @@ 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; - } - } } void FlatSourceDomain::flux_swap() @@ -1159,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); @@ -1206,7 +1151,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, @@ -1220,7 +1165,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; } @@ -1230,7 +1175,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 { @@ -1239,6 +1184,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; } } @@ -1265,7 +1211,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 @@ -1326,13 +1271,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; } } } @@ -1402,18 +1348,9 @@ 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) + SourceRegionKey sr_key, Position r, 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 @@ -1475,9 +1412,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; @@ -1485,9 +1421,9 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( } // Sanity check on mesh bin - int mesh_idx = base_source_regions_.mesh(sr); + 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; @@ -1496,7 +1432,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; @@ -1508,26 +1444,60 @@ 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. - SourceRegion* sr_ptr = discovered_source_regions_.emplace( - sr_key, {base_source_regions_.get_source_region_handle(sr), sr}); - discovered_source_regions_.unlock(sr_key); + // Additionally, we need to determine the source region's material, initialize + // the starting scalar flux guess, and apply any known external sources. + + // 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}; - // 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) { + // 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()); + + // 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 + handle.mesh() = mesh_idx; + + 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); + } + } + + // 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); + } + } + + // 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; @@ -1535,6 +1505,21 @@ 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. + discovered_source_regions_.unlock(sr_key); + return handle; } @@ -1620,4 +1605,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/linear_source_domain.cpp b/src/random_ray/linear_source_domain.cpp index 81412164eca..e1ad68e3d87 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(double k_eff) +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(double k_eff) 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(double k_eff) } // 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 @@ -91,25 +84,21 @@ void LinearSourceDomain::update_neutron_source(double k_eff) // 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 && - source_regions_.source(sr, g_out) >= 0.0) { - source_regions_.source_gradients(sr, g_out) = + 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 { - 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.cpp b/src/random_ray/random_ray.cpp index 27f674c4278..89a91449df7 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() @@ -336,71 +335,60 @@ 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 - if (mesh_subdivision_enabled_) { - // Determine the mesh index for the base source region, if any - int mesh_idx = domain_->base_source_regions_.mesh(sr); - - 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(); - } - } + // Determine the mesh index for the base source region, if any + int mesh_idx = domain_->lookup_mesh_idx(sr); + + 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(); + } } } 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; - 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_key, r, u()); + if (srh.is_numerical_fp_artifact_) { + return; } 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); @@ -408,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); @@ -439,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++) { @@ -490,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 @@ -537,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); @@ -810,27 +798,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(), 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; - 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()); - } else { - srh = domain_->source_regions_.get_source_region_handle(sr); - } - if (!srh.is_numerical_fp_artifact_) { for (int g = 0; g < negroups_; g++) { angular_flux_[g] = srh.source(g); diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 388a778b840..d475b2593ea 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -47,97 +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; - SourceRegionContainer forward_base_source_regions; - std::unordered_map - forward_source_region_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; - } + for (uint64_t se = 0; se < sim.domain()->n_source_elements(); se++) { + sim.domain()->source_regions_.scalar_flux_final(se) *= + source_normalization_factor; + } - 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(); - // 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 Random Ray Simulation Object - RandomRaySimulation adjoint_sim; + // Initialize OpenMC general data structures + openmc_simulation_init(); - // Initialize adjoint fixed sources, if present - adjoint_sim.prepare_fixed_sources_adjoint(forward_flux, - forward_source_regions, forward_base_source_regions, - forward_source_region_map); + sim.domain()->k_eff_ = 1.0; - // Transpose scattering matrix - adjoint_sim.domain()->transpose_scattering_matrix(); + // Initialize adjoint fixed sources, if present + sim.prepare_fixed_sources_adjoint(); - // Swap nu_sigma_f and chi - adjoint_sim.domain()->nu_sigma_f_.swap(adjoint_sim.domain()->chi_); + // Transpose scattering matrix + sim.domain()->transpose_scattering_matrix(); - // Begin main simulation timer - simulation::time_total.start(); + // Swap nu_sigma_f and chi + sim.domain()->nu_sigma_f_.swap(sim.domain()->chi_); - // Execute random ray simulation - adjoint_sim.simulate(); + // Begin main simulation timer + simulation::time_total.start(); - // End main simulation timer - simulation::time_total.stop(); + // Execute random ray simulation + sim.simulate(); - // Finalize OpenMC - openmc_simulation_finalize(); + // End main simulation timer + simulation::time_total.stop(); - // Output all simulation results - adjoint_sim.output_simulation_results(); - } + // Finalize OpenMC + openmc_simulation_finalize(); + + // Output all simulation results + sim.output_simulation_results(); } // Enforces restrictions on inputs in random ray mode. While there are @@ -348,7 +333,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 +350,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; } @@ -412,20 +395,11 @@ 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) +void RandomRaySimulation::prepare_fixed_sources_adjoint() { + domain_->source_regions_.adjoint_reset(); 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_->set_adjoint_sources(forward_flux); + domain_->set_adjoint_sources(); } } @@ -445,22 +419,18 @@ void RandomRaySimulation::simulate() simulation::total_weight = 1.0; // Update source term (scattering + fission) - domain_->update_neutron_source(k_eff_); + domain_->update_all_neutron_sources(); - // Reset scalar fluxes, iteration volume tallies, and region hit flags to - // zero + // 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 + // 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 // 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(); @@ -476,11 +446,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_) { - domain_->finalize_discovered_source_regions(); - } + // 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( @@ -494,10 +462,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 @@ -507,12 +475,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(); @@ -522,7 +484,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 @@ -571,7 +533,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/random_ray/source_region.cpp b/src/random_ray/source_region.cpp index 1205b995a6f..3b06f0ed09f 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); @@ -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 //============================================================================== @@ -259,9 +240,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(), 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; } } } 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() 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