diff --git a/docs/source/usersguide/random_ray.rst b/docs/source/usersguide/random_ray.rst index d35aff83b..06d4e50ed 100644 --- a/docs/source/usersguide/random_ray.rst +++ b/docs/source/usersguide/random_ray.rst @@ -1105,10 +1105,15 @@ external source is present in the problem. Simulation settings (e.g., number of rays, batches, etc.) will be identical for both calculations. At the conclusion of the run, all results (e.g., tallies, plots, etc.) will be derived from the adjoint flux rather than the forward flux but are not labeled -any differently. The initial forward flux solution will not be stored or -available in the final statepoint file. Those wishing to do analysis requiring -both the forward and adjoint solutions will need to run two separate -simulations and load both statepoint files. +any differently. When an initial forward solve is performed (i.e., when no +user-specified adjoint source is present), its output files are also written to +disk with a ``forward`` infix, so they are not overwritten by the subsequent +adjoint solve. This applies to the statepoint, ``tallies.out``, and any voxel +plots, e.g., ``statepoint.forward.N.h5`` and ``tallies.forward.out``; the +adjoint solve keeps the usual file names. This allows analyses requiring both +the forward and adjoint solutions to be performed from a single run. When +generating FW-CADIS weight windows, no weight window file is written for the +forward solve, as only the final adjoint-derived weight windows are meaningful. .. note:: Use of the automated diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 7baed25c5..26b4a224c 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -368,6 +368,7 @@ enum class SolverType { MONTE_CARLO, RANDOM_RAY }; enum class RandomRayVolumeEstimator { NAIVE, SIMULATION_AVERAGED, HYBRID }; enum class RandomRaySourceShape { FLAT, LINEAR, LINEAR_XY }; enum class RandomRaySampleMethod { PRNG, HALTON, S2 }; +enum class RandomRaySolve { FORWARD, FORWARD_FOR_ADJOINT, ADJOINT }; //============================================================================== // Geometry Constants diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index 6f51af34d..09414fd44 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -76,7 +76,10 @@ public: //---------------------------------------------------------------------------- // Static Data members static bool volume_normalized_flux_tallies_; - static bool adjoint_; // If the user wants outputs based on the adjoint flux + // If the user wants outputs based on the adjoint flux + static bool adjoint_requested_; + // The solve currently being executed + static RandomRaySolve solve_; static bool fw_cadis_local_; static double diagonal_stabilization_rho_; // Adjusts strength of diagonal stabilization diff --git a/include/openmc/random_ray/random_ray_simulation.h b/include/openmc/random_ray/random_ray_simulation.h index ccd2cbe47..e186c549f 100644 --- a/include/openmc/random_ray/random_ray_simulation.h +++ b/include/openmc/random_ray/random_ray_simulation.h @@ -22,7 +22,7 @@ public: void apply_fixed_sources_and_mesh_domains(); void prepare_fw_fixed_sources_adjoint(); void prepare_local_fixed_sources_adjoint(); - void prepare_adjoint_simulation(bool fw_adjoint); + void prepare_adjoint_simulation(bool from_forward); void simulate(); void output_simulation_results() const; void instability_check( diff --git a/src/output.cpp b/src/output.cpp index 725e8d069..04a8caa60 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -621,8 +621,15 @@ void write_tallies() if (model::tallies.empty()) return; + // Tag tallies.out written during the forward solve of an adjoint run + const char* forward = + (FlatSourceDomain::solve_ == RandomRaySolve::FORWARD_FOR_ADJOINT) + ? "forward." + : ""; + // Set filename for tallies_out - std::string filename = fmt::format("{}tallies.out", settings::path_output); + std::string filename = + fmt::format("{}tallies.{}out", settings::path_output, forward); // Open the tallies.out file. std::ofstream tallies_out; diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 06c6ef14d..895e0663e 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -30,7 +30,8 @@ namespace openmc { RandomRayVolumeEstimator FlatSourceDomain::volume_estimator_ { RandomRayVolumeEstimator::HYBRID}; bool FlatSourceDomain::volume_normalized_flux_tallies_ {false}; -bool FlatSourceDomain::adjoint_ {false}; +bool FlatSourceDomain::adjoint_requested_ {false}; +RandomRaySolve FlatSourceDomain::solve_ {RandomRaySolve::FORWARD}; bool FlatSourceDomain::fw_cadis_local_ {false}; double FlatSourceDomain::diagonal_stabilization_rho_ {1.0}; std::unordered_map>> @@ -556,7 +557,7 @@ double FlatSourceDomain::compute_fixed_source_normalization_factor() const // If we are in adjoint mode of a fixed source problem, the external // source is already normalized, such that all resulting fluxes are // also normalized. - if (adjoint_) { + if (solve_ == RandomRaySolve::ADJOINT) { return 1.0; } @@ -795,6 +796,12 @@ void FlatSourceDomain::output_to_vtk() const double z_delta = width.z / Nz; std::string filename = openmc_plot->path_plot(); + // Tag plots written during the forward solve of an adjoint run + if (solve_ == RandomRaySolve::FORWARD_FOR_ADJOINT) { + auto dot = filename.find_last_of('.'); + filename = filename.substr(0, dot) + ".forward" + filename.substr(dot); + } + // Perform sanity checks on file size uint64_t bytes = Nx * Ny * Nz * (negroups_ + 1 + 1 + 1) * sizeof(float); write_message(5, "Processing plot {}: {}... (Estimated size is {} MB)", @@ -1002,9 +1009,10 @@ void FlatSourceDomain::output_to_vtk() const void FlatSourceDomain::apply_external_source_to_source_region( int src_idx, SourceRegionHandle& srh) { - auto s = (adjoint_ && !model::adjoint_sources.empty()) - ? model::adjoint_sources[src_idx].get() - : model::external_sources[src_idx].get(); + auto s = + (solve_ == RandomRaySolve::ADJOINT && !model::adjoint_sources.empty()) + ? model::adjoint_sources[src_idx].get() + : model::external_sources[src_idx].get(); auto is = dynamic_cast(s); auto discrete = dynamic_cast(is->energy()); double strength_factor = is->strength(); diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 24650afb8..00fff99a7 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -189,7 +189,7 @@ void validate_random_ray_inputs() // Validate adjoint sources /////////////////////////////////////////////////////////////////// - if (FlatSourceDomain::adjoint_ && !model::adjoint_sources.empty()) { + if (FlatSourceDomain::adjoint_requested_ && !model::adjoint_sources.empty()) { for (int i = 0; i < model::adjoint_sources.size(); i++) { Source* s = model::adjoint_sources[i].get(); @@ -289,7 +289,8 @@ void openmc_finalize_random_ray() { FlatSourceDomain::volume_estimator_ = RandomRayVolumeEstimator::HYBRID; FlatSourceDomain::volume_normalized_flux_tallies_ = false; - FlatSourceDomain::adjoint_ = false; + FlatSourceDomain::adjoint_requested_ = false; + FlatSourceDomain::solve_ = RandomRaySolve::FORWARD; FlatSourceDomain::fw_cadis_local_ = false; FlatSourceDomain::fw_cadis_local_targets_.clear(); FlatSourceDomain::mesh_domain_map_.clear(); @@ -356,19 +357,17 @@ void RandomRaySimulation::prepare_local_fixed_sources_adjoint() } } -void RandomRaySimulation::prepare_adjoint_simulation(bool fw_adjoint) +void RandomRaySimulation::prepare_adjoint_simulation(bool from_forward) { reset_timers(); if (mpi::master) header("ADJOINT FLUX SOLVE", 3); - if (fw_adjoint) { - // Forward simulation has already been run; - // Configure the domain for adjoint simulation and - // re-initialize OpenMC general data structures - FlatSourceDomain::adjoint_ = true; - + if (from_forward) { + // The forward solve has already run. Re-initialize OpenMC's general data + // structures for the adjoint solve and derive the adjoint source from the + // forward flux. openmc_simulation_init(); prepare_fw_fixed_sources_adjoint(); @@ -603,7 +602,8 @@ void RandomRaySimulation::print_results_random_ray( } fmt::print(" Volume Estimator Type = {}\n", estimator); - std::string adjoint_true = (FlatSourceDomain::adjoint_) ? "ON" : "OFF"; + std::string adjoint_true = + (FlatSourceDomain::solve_ == RandomRaySolve::ADJOINT) ? "ON" : "OFF"; fmt::print(" Adjoint Flux Mode = {}\n", adjoint_true); std::string shape; @@ -675,60 +675,49 @@ void RandomRaySimulation::print_results_random_ray( void openmc_run_random_ray() { - ////////////////////////////////////////////////////////// - // Run forward simulation - ////////////////////////////////////////////////////////// + using namespace openmc; - // Check if adjoint calculation is needed, and if local adjoint source(s) - // are present. If an adjoint calculation is needed and no sources are - // specified, we will run a forward calculation first to calculate adjoint - // sources for global variance reduction, then perform an adjoint - // calculation later. - bool adjoint_needed = openmc::FlatSourceDomain::adjoint_; - bool fw_adjoint = openmc::model::adjoint_sources.empty() && adjoint_needed; + // Determine which solves to run. If adjoint results are requested and no + // user-defined adjoint source is present, an initial forward solve is needed + // to construct the adjoint source from the forward flux (FW-CADIS). If the + // user has defined an adjoint source, the forward solve is skipped and only + // the adjoint solve is run. + const bool run_adjoint = FlatSourceDomain::adjoint_requested_; + const bool have_adjoint_source = !model::adjoint_sources.empty(); + const bool run_forward = !(run_adjoint && have_adjoint_source); - // If we're going to do an adjoint simulation with forward-weighted adjoint - // sources afterwards, report that this is the initial forward flux solve. - if (!adjoint_needed || fw_adjoint) { - // Configure the domain for forward simulation - openmc::FlatSourceDomain::adjoint_ = false; - - if (adjoint_needed && openmc::mpi::master) - openmc::header("FORWARD FLUX SOLVE", 3); + // Set the initial solve type + if (!run_forward) { + FlatSourceDomain::solve_ = RandomRaySolve::ADJOINT; + } else if (run_adjoint) { + FlatSourceDomain::solve_ = RandomRaySolve::FORWARD_FOR_ADJOINT; } else { - // Configure domain for adjoint simulation (later) - openmc::FlatSourceDomain::adjoint_ = true; + FlatSourceDomain::solve_ = RandomRaySolve::FORWARD; } // Initialize OpenMC general data structures openmc_simulation_init(); // Validate that inputs meet requirements for random ray mode - if (openmc::mpi::master) - openmc::validate_random_ray_inputs(); + if (mpi::master) + validate_random_ray_inputs(); // Initialize Random Ray Simulation Object - openmc::RandomRaySimulation sim; + RandomRaySimulation sim; - if (!adjoint_needed || fw_adjoint) { - // Initialize fixed sources, if present + // Run the forward solve + if (run_forward) { + // When an adjoint solve follows, report this as the initial forward solve + if (run_adjoint && mpi::master) + header("FORWARD FLUX SOLVE", 3); sim.apply_fixed_sources_and_mesh_domains(); - - // Execute random ray simulation sim.simulate(); } - ////////////////////////////////////////////////////////// - // Run adjoint simulation (if enabled) - ////////////////////////////////////////////////////////// - - if (!adjoint_needed) { - return; + // Run the adjoint solve + if (run_adjoint) { + FlatSourceDomain::solve_ = RandomRaySolve::ADJOINT; + sim.prepare_adjoint_simulation(run_forward); + sim.simulate(); } - - // Setup for adjoint simulation - sim.prepare_adjoint_simulation(fw_adjoint); - - // Execute random ray simulation - sim.simulate(); } diff --git a/src/settings.cpp b/src/settings.cpp index 8ae252ae1..58bced940 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -322,7 +322,7 @@ void get_run_parameters(pugi::xml_node node_base) get_node_value_bool(random_ray_node, "volume_normalized_flux_tallies"); } if (check_for_node(random_ray_node, "adjoint")) { - FlatSourceDomain::adjoint_ = + FlatSourceDomain::adjoint_requested_ = get_node_value_bool(random_ray_node, "adjoint"); } if (check_for_node(random_ray_node, "sample_method")) { diff --git a/src/simulation.cpp b/src/simulation.cpp index 89aa9ca0e..900a383ab 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -16,6 +16,7 @@ #include "openmc/particle.h" #include "openmc/photon.h" #include "openmc/random_lcg.h" +#include "openmc/random_ray/flat_source_domain.h" #include "openmc/settings.h" #include "openmc/source.h" #include "openmc/state_point.h" @@ -200,9 +201,12 @@ int openmc_simulation_finalize() if (settings::output_tallies && mpi::master) write_tallies(); - // If weight window generators are present in this simulation, - // write a weight windows file - if (variance_reduction::weight_windows_generators.size() > 0) { + // If weight window generators are present in this simulation, write a + // weight windows file. This is skipped during the forward solve of an + // adjoint (FW-CADIS) run, where only the adjoint-derived weight windows + // are meaningful. + if (variance_reduction::weight_windows_generators.size() > 0 && + FlatSourceDomain::solve_ != RandomRaySolve::FORWARD_FOR_ADJOINT) { openmc_weight_windows_export(); } diff --git a/src/state_point.cpp b/src/state_point.cpp index da1c141a2..02c68bf60 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -22,6 +22,7 @@ #include "openmc/nuclide.h" #include "openmc/output.h" #include "openmc/particle_type.h" +#include "openmc/random_ray/flat_source_domain.h" #include "openmc/settings.h" #include "openmc/simulation.h" #include "openmc/tallies/derivative.h" @@ -46,9 +47,15 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source) // Determine width for zero padding int w = std::to_string(settings::n_max_batches).size(); + // Tag statepoints written during the forward solve of an adjoint run + const char* forward = + (FlatSourceDomain::solve_ == RandomRaySolve::FORWARD_FOR_ADJOINT) + ? "forward." + : ""; + // Set filename for state point - filename_ = fmt::format("{0}statepoint.{1:0{2}}.h5", settings::path_output, - simulation::current_batch, w); + filename_ = fmt::format("{0}statepoint.{3}{1:0{2}}.h5", + settings::path_output, simulation::current_batch, w, forward); } // If a file name was specified, ensure it has .h5 file extension diff --git a/src/weight_windows.cpp b/src/weight_windows.cpp index 0614110cd..39e46026c 100644 --- a/src/weight_windows.cpp +++ b/src/weight_windows.cpp @@ -791,7 +791,7 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node) if (method_string == "magic") { method_ = WeightWindowUpdateMethod::MAGIC; if (settings::solver_type == SolverType::RANDOM_RAY && - FlatSourceDomain::adjoint_) { + FlatSourceDomain::adjoint_requested_) { fatal_error("Random ray weight window generation with MAGIC cannot be " "done in adjoint mode."); } @@ -800,7 +800,7 @@ WeightWindowsGenerator::WeightWindowsGenerator(pugi::xml_node node) if (settings::solver_type != SolverType::RANDOM_RAY) { fatal_error("FW-CADIS can only be run in random ray solver mode."); } - FlatSourceDomain::adjoint_ = true; + FlatSourceDomain::adjoint_requested_ = true; if (check_for_node(node, "targets")) { FlatSourceDomain::fw_cadis_local_ = true; targets_ = get_node_array(node, "targets"); diff --git a/tests/testing_harness.py b/tests/testing_harness.py index 8d156bd64..c1977556d 100644 --- a/tests/testing_harness.py +++ b/tests/testing_harness.py @@ -143,7 +143,8 @@ class TestHarness: def _cleanup(self): """Delete statepoints, tally, and test files.""" output = glob.glob('statepoint.*.h5') - output += ['tallies.out', 'results_test.dat', 'summary.h5'] + output += ['tallies.out', 'tallies.forward.out'] + output += ['results_test.dat', 'summary.h5'] output += glob.glob('volume_*.h5') for f in output: if os.path.exists(f):