diff --git a/include/openmc/random_ray/random_ray_simulation.h b/include/openmc/random_ray/random_ray_simulation.h index 3dec48bf266..3db6510697c 100644 --- a/include/openmc/random_ray/random_ray_simulation.h +++ b/include/openmc/random_ray/random_ray_simulation.h @@ -19,9 +19,9 @@ class RandomRaySimulation { //---------------------------------------------------------------------------- // Methods - void compute_segment_correction_factors(); void apply_fixed_sources_and_mesh_domains(); void prepare_fixed_sources_adjoint(); + void prepare_adjoint_simulation(); void simulate(); void output_simulation_results() const; void instability_check( @@ -34,9 +34,15 @@ class RandomRaySimulation { // Accessors FlatSourceDomain* domain() const { return domain_.get(); } + //---------------------------------------------------------------------------- + // Public data members + + // Flag for adjoint simulation; + bool adjoint_needed_; + private: //---------------------------------------------------------------------------- - // Data members + // Private data members // Contains all flat source region data unique_ptr domain_; @@ -51,6 +57,9 @@ class RandomRaySimulation { // Number of energy groups int negroups_; + // Toggle for first simulation + bool is_first_simulation_; + }; // class RandomRaySimulation //============================================================================ @@ -60,6 +69,7 @@ class RandomRaySimulation { void openmc_run_random_ray(); void validate_random_ray_inputs(); void openmc_reset_random_ray(); +void print_adjoint_header(); } // namespace openmc diff --git a/openmc/model/model.py b/openmc/model/model.py index 6e4c1c5856d..ca78f282f6a 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -1179,8 +1179,8 @@ def plot( # Convert ID map to RGB image img = id_map_to_rgb( - id_map=id_map, - color_by=color_by, + id_map=id_map, + color_by=color_by, colors=colors, overlap_color=overlap_color ) @@ -1217,7 +1217,7 @@ def plot( extent=(x_min, x_max, y_min, y_max), **contour_kwargs ) - + # If only showing outline, set the axis limits and aspect explicitly if outline == 'only': axes.set_xlim(x_min, x_max) @@ -1685,6 +1685,85 @@ def differentiate_mats(self, diff_volume_method: str = None, depletable_only: bo self.geometry.get_all_materials().values() ) + def _auto_generate_mgxs_lib( + self, + model: openmc.model.model, + groups: openmc.mgxs.EnergyGroups, + correction: str | none, + directory: pathlike, + ) -> openmc.mgxs.Library: + """ + Automatically generate a multi-group cross section libray from a model + with the specified group structure. + + Parameters + ---------- + groups : openmc.mgxs.EnergyGroups + Energy group structure for the MGXS. + nparticles : int + Number of particles to simulate per batch when generating MGXS. + mgxs_path : str + Filename for the MGXS HDF5 file. + correction : str + Transport correction to apply to the MGXS. Options are None and + "P0". + directory : str + Directory to run the simulation in, so as to contain XML files. + + Returns + ------- + mgxs_lib : openmc.mgxs.Library + OpenMC MGXS Library object + """ + + # Initialize MGXS library with a finished OpenMC geometry object + mgxs_lib = openmc.mgxs.Library(model.geometry) + + # Pick energy group structure + mgxs_lib.energy_groups = groups + + # Disable transport correction + mgxs_lib.correction = correction + + # Specify needed cross sections for random ray + if correction == 'P0': + mgxs_lib.mgxs_types = [ + 'nu-transport', 'absorption', 'nu-fission', 'fission', + 'consistent nu-scatter matrix', 'multiplicity matrix', 'chi' + ] + elif correction is None: + mgxs_lib.mgxs_types = [ + 'total', 'absorption', 'nu-fission', 'fission', + 'consistent nu-scatter matrix', 'multiplicity matrix', 'chi' + ] + + # Specify a "material" domain type for the cross section tally filters + mgxs_lib.domain_type = "material" + + # Specify the domains over which to compute multi-group cross sections + mgxs_lib.domains = model.geometry.get_all_materials().values() + + # Do not compute cross sections on a nuclide-by-nuclide basis + mgxs_lib.by_nuclide = False + + # Check the library - if no errors are raised, then the library is satisfactory. + mgxs_lib.check_library_for_openmc_mgxs() + + # Construct all tallies needed for the multi-group cross section library + mgxs_lib.build_library() + + # Create a "tallies.xml" file for the MGXS Library + mgxs_lib.add_to_tallies(model.tallies, merge=True) + + # Run + statepoint_filename = model.run(cwd=directory) + + # Load MGXS + with openmc.StatePoint(statepoint_filename) as sp: + mgxs_lib.load_from_statepoint(sp) + + return mgxs_lib + def _create_mgxs_sources( self, groups: openmc.mgxs.EnergyGroups, @@ -1848,52 +1927,8 @@ def _generate_infinite_medium_mgxs( model.geometry.root_universe = infinite_universe # Add MGXS Tallies - - # Initialize MGXS library with a finished OpenMC geometry object - mgxs_lib = openmc.mgxs.Library(model.geometry) - - # Pick energy group structure - mgxs_lib.energy_groups = groups - - # Disable transport correction - mgxs_lib.correction = correction - - # Specify needed cross sections for random ray - if correction == 'P0': - mgxs_lib.mgxs_types = [ - 'nu-transport', 'absorption', 'nu-fission', 'fission', - 'consistent nu-scatter matrix', 'multiplicity matrix', 'chi' - ] - elif correction is None: - mgxs_lib.mgxs_types = [ - 'total', 'absorption', 'nu-fission', 'fission', - 'consistent nu-scatter matrix', 'multiplicity matrix', 'chi' - ] - - # Specify a "cell" domain type for the cross section tally filters - mgxs_lib.domain_type = "material" - - # Specify the cell domains over which to compute multi-group cross sections - mgxs_lib.domains = model.geometry.get_all_materials().values() - - # Do not compute cross sections on a nuclide-by-nuclide basis - mgxs_lib.by_nuclide = False - - # Check the library - if no errors are raised, then the library is satisfactory. - mgxs_lib.check_library_for_openmc_mgxs() - - # Construct all tallies needed for the multi-group cross section library - mgxs_lib.build_library() - - # Create a "tallies.xml" file for the MGXS Library - mgxs_lib.add_to_tallies(model.tallies, merge=True) - - # Run - statepoint_filename = model.run(cwd=directory) - - # Load MGXS - with openmc.StatePoint(statepoint_filename) as sp: - mgxs_lib.load_from_statepoint(sp) + mgxs_lib = self._auto_generate_mgxs_lib( + model, groups, correction, directory) # Create a MGXS File which can then be written to disk mgxs_set = mgxs_lib.get_xsdata(domain=material, xsdata_name=name) @@ -2055,48 +2090,8 @@ def _generate_stochastic_slab_mgxs( model.settings.output = {'summary': True, 'tallies': False} # Add MGXS Tallies - - # Initialize MGXS library with a finished OpenMC geometry object - mgxs_lib = openmc.mgxs.Library(model.geometry) - - # Pick energy group structure - mgxs_lib.energy_groups = groups - - # Disable transport correction - mgxs_lib.correction = correction - - # Specify needed cross sections for random ray - if correction == 'P0': - mgxs_lib.mgxs_types = ['nu-transport', 'absorption', 'nu-fission', 'fission', - 'consistent nu-scatter matrix', 'multiplicity matrix', 'chi'] - elif correction is None: - mgxs_lib.mgxs_types = ['total', 'absorption', 'nu-fission', 'fission', - 'consistent nu-scatter matrix', 'multiplicity matrix', 'chi'] - - # Specify a "cell" domain type for the cross section tally filters - mgxs_lib.domain_type = "material" - - # Specify the cell domains over which to compute multi-group cross sections - mgxs_lib.domains = model.geometry.get_all_materials().values() - - # Do not compute cross sections on a nuclide-by-nuclide basis - mgxs_lib.by_nuclide = False - - # Check the library - if no errors are raised, then the library is satisfactory. - mgxs_lib.check_library_for_openmc_mgxs() - - # Construct all tallies needed for the multi-group cross section library - mgxs_lib.build_library() - - # Create a "tallies.xml" file for the MGXS Library - mgxs_lib.add_to_tallies(model.tallies, merge=True) - - # Run - statepoint_filename = model.run(cwd=directory) - - # Load MGXS - with openmc.StatePoint(statepoint_filename) as sp: - mgxs_lib.load_from_statepoint(sp) + mgxs_lib = self._auto_generate_mgxs_lib( + model, groups, correction, directory) names = [mat.name for mat in mgxs_lib.domains] @@ -2146,52 +2141,8 @@ def _generate_material_wise_mgxs( model.settings.output = {'summary': True, 'tallies': False} # Add MGXS Tallies - - # Initialize MGXS library with a finished OpenMC geometry object - mgxs_lib = openmc.mgxs.Library(model.geometry) - - # Pick energy group structure - mgxs_lib.energy_groups = groups - - # Disable transport correction - mgxs_lib.correction = correction - - # Specify needed cross sections for random ray - if correction == 'P0': - mgxs_lib.mgxs_types = [ - 'nu-transport', 'absorption', 'nu-fission', 'fission', - 'consistent nu-scatter matrix', 'multiplicity matrix', 'chi' - ] - elif correction is None: - mgxs_lib.mgxs_types = [ - 'total', 'absorption', 'nu-fission', 'fission', - 'consistent nu-scatter matrix', 'multiplicity matrix', 'chi' - ] - - # Specify a "cell" domain type for the cross section tally filters - mgxs_lib.domain_type = "material" - - # Specify the cell domains over which to compute multi-group cross sections - mgxs_lib.domains = model.geometry.get_all_materials().values() - - # Do not compute cross sections on a nuclide-by-nuclide basis - mgxs_lib.by_nuclide = False - - # Check the library - if no errors are raised, then the library is satisfactory. - mgxs_lib.check_library_for_openmc_mgxs() - - # Construct all tallies needed for the multi-group cross section library - mgxs_lib.build_library() - - # Create a "tallies.xml" file for the MGXS Library - mgxs_lib.add_to_tallies(model.tallies, merge=True) - - # Run - statepoint_filename = model.run(cwd=directory) - - # Load MGXS - with openmc.StatePoint(statepoint_filename) as sp: - mgxs_lib.load_from_statepoint(sp) + mgxs_lib = self._auto_generate_mgxs_lib( + model, groups, correction, directory) names = [mat.name for mat in mgxs_lib.domains] @@ -2617,5 +2568,3 @@ def function_calls(self) -> int: def total_batches(self) -> int: """Total number of active batches used across all evaluations.""" return sum(self.batches) - - diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index d475b2593ea..2b54ae2eb79 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -28,17 +28,13 @@ void openmc_run_random_ray() // Run forward simulation ////////////////////////////////////////////////////////// - // Check if adjoint calculation is needed. If it is, we will run the forward - // calculation first and then the adjoint calculation later. - bool adjoint_needed = FlatSourceDomain::adjoint_; - - // Configure the domain for forward simulation - FlatSourceDomain::adjoint_ = false; - - // If we're going to do an adjoint simulation afterwards, report that this is - // the initial forward flux solve. - if (adjoint_needed && mpi::master) - header("FORWARD FLUX SOLVE", 3); + if (mpi::master) { + if (FlatSourceDomain::adjoint_) { + FlatSourceDomain::adjoint_ = false; + openmc::print_adjoint_header(); + FlatSourceDomain::adjoint_ = true; + } + } // Initialize OpenMC general data structures openmc_simulation_init(); @@ -53,76 +49,20 @@ void openmc_run_random_ray() // Initialize fixed sources, if present sim.apply_fixed_sources_and_mesh_domains(); - // Begin main simulation timer - simulation::time_total.start(); - - // Execute random ray simulation + // Run initial random ray simulation sim.simulate(); - // End main simulation timer - simulation::time_total.stop(); - - // 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 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(); - - // Output all simulation results - sim.output_simulation_results(); - ////////////////////////////////////////////////////////// // Run adjoint simulation (if enabled) ////////////////////////////////////////////////////////// - if (!adjoint_needed) { - return; - } - - reset_timers(); - - // Configure the domain for adjoint simulation - FlatSourceDomain::adjoint_ = true; - - if (mpi::master) - header("ADJOINT FLUX SOLVE", 3); - - // Initialize OpenMC general data structures - openmc_simulation_init(); - - sim.domain()->k_eff_ = 1.0; - - // Initialize adjoint fixed sources, if present - sim.prepare_fixed_sources_adjoint(); - - // Transpose scattering matrix - sim.domain()->transpose_scattering_matrix(); - - // Swap nu_sigma_f and chi - sim.domain()->nu_sigma_f_.swap(sim.domain()->chi_); - - // Begin main simulation timer - simulation::time_total.start(); - - // Execute random ray simulation - sim.simulate(); - - // End main simulation timer - simulation::time_total.stop(); - - // Finalize OpenMC - openmc_simulation_finalize(); + if (sim.adjoint_needed_) { + // Setup for adjoint simulation + sim.prepare_adjoint_simulation(); - // Output all simulation results - sim.output_simulation_results(); + // Run adjoint simulation + sim.simulate(); + } } // Enforces restrictions on inputs in random ray mode. While there are @@ -353,6 +293,17 @@ void openmc_reset_random_ray() RandomRay::sample_method_ = RandomRaySampleMethod::PRNG; } +void print_adjoint_header() +{ + if (!FlatSourceDomain::adjoint_) + // If we're going to do an adjoint simulation afterwards, report that this + // is the initial forward flux solve. + header("FORWARD FLUX SOLVE", 3); + else + // Otherwise report that we are doing the adjoint simulation + header("ADJOINT FLUX SOLVE", 3); +} + //============================================================================== // RandomRaySimulation implementation //============================================================================== @@ -383,6 +334,16 @@ RandomRaySimulation::RandomRaySimulation() // Convert OpenMC native MGXS into a more efficient format // internal to the random ray solver domain_->flatten_xs(); + + // Check if adjoint calculation is needed. If it is, we will run the forward + // calculation first and then the adjoint calculation later. + adjoint_needed_ = FlatSourceDomain::adjoint_; + + // Adjoint is always false for the forward calculation + FlatSourceDomain::adjoint_ = false; + + // The first simulation is run after initialization + is_first_simulation_ = true; } void RandomRaySimulation::apply_fixed_sources_and_mesh_domains() @@ -403,8 +364,41 @@ void RandomRaySimulation::prepare_fixed_sources_adjoint() } } +void RandomRaySimulation::prepare_adjoint_simulation() +{ + // Configure the domain for adjoint simulation + FlatSourceDomain::adjoint_ = true; + + // Reset k-eff + domain_->k_eff_ = 1.0; + + // Initialize adjoint fixed sources, if present + prepare_fixed_sources_adjoint(); + + // Transpose scattering matrix + domain_->transpose_scattering_matrix(); + + // Swap nu_sigma_f and chi + domain_->nu_sigma_f_.swap(domain_->chi_); +} + void RandomRaySimulation::simulate() { + if (!is_first_simulation_) { + if (mpi::master && adjoint_needed_) + openmc::print_adjoint_header(); + + // Reset the timers and reinitialize the general OpenMC datastructures if + // this is after the first simulation + reset_timers(); + + // Initialize OpenMC general data structures + openmc_simulation_init(); + } + + // Begin main simulation timer + simulation::time_total.start(); + // Random ray power iteration loop while (simulation::current_batch < settings::n_batches) { // Initialize the current batch @@ -493,6 +487,31 @@ void RandomRaySimulation::simulate() } // End random ray power iteration loop domain_->count_external_source_regions(); + + // End main simulation timer + simulation::time_total.stop(); + + // Normalize and save the final flux + double source_normalization_factor = + domain_->compute_fixed_source_normalization_factor() / + (settings::n_batches - settings::n_inactive); + +#pragma omp parallel for + for (uint64_t se = 0; se < domain_->n_source_elements(); se++) { + domain_->source_regions_.scalar_flux_final(se) *= + source_normalization_factor; + } + + // Finalize OpenMC + openmc_simulation_finalize(); + + // Output all simulation results + output_simulation_results(); + + // Toggle that the simulation object has been initialized after the first + // simulation + if (is_first_simulation_) + is_first_simulation_ = false; } void RandomRaySimulation::output_simulation_results() const