diff --git a/docs/source/methods/index.rst b/docs/source/methods/index.rst index 121d04b1ded..c761e461308 100644 --- a/docs/source/methods/index.rst +++ b/docs/source/methods/index.rst @@ -17,6 +17,7 @@ Theory and Methodology charged_particles_physics tallies eigenvalue + subcritical_multiplication depletion energy_deposition parallelization diff --git a/docs/source/methods/subcritical_multiplication.rst b/docs/source/methods/subcritical_multiplication.rst new file mode 100644 index 00000000000..3ee4676f3bc --- /dev/null +++ b/docs/source/methods/subcritical_multiplication.rst @@ -0,0 +1,100 @@ +.. _methods_subcritical-multiplication: + +======================================= +Subcritical Multiplication Calculations +======================================= + +Subcritical multiplication problems are fixed source simulations with a large +amount of multiplication that are still subcritical with respect to +:math:`k_{eff}`. These problems are common in the design of accelerator driven +systems (ADS) which use an accelerator source to drive a subcritical +multiplication reaction in, e.g., spent fuel to transmute nuclear waste and +even generate power [Bowman]_. An ADS is inherently safe and allows for a much +more flexible fuel composition, hence their popularity for waste transmutation. + +For ADS's, the production of fission neutrons is central as these produced +neutrons amplify the external source. For the case of a proton spallation ADS, +source neutrons are produced from spallation reactions in a heavy metal target +bombarded by high-energy protons. These source neutrons then induce fission in +a subcritical core, producing additional neutrons. + + +.. _methods_subcritical-multiplication-factors: + +---------------------------------- +Subcritical Multiplication Factors +---------------------------------- +In a fixed source simulation, the total neutron production (per source particle) is used to define + +.. math:: + :label: integral_multiplicity + + \begin{align*} + M - 1 &= \frac{1}{N_{source}}\int d\boldsymbol{r} \int d\boldsymbol{\Omega} \int dE \nu \Sigma_f(\boldsymbol{r}, E) \psi(\boldsymbol{r}, \boldsymbol{\Omega}, E)\\ + &\coloneqq k + k^2 + k^3 + ... \\ + &= \frac{k}{1-k}\\ + \implies M &= \frac{1}{1-k} + \end{align*} + +Where :math:`M` is the subcritical multiplicity, and :math:`k` the subcritical +multiplication factor. The identification on the second line comes from the +picture of a single source neutron inducing several generations of fission +neutrons, producing on average :math:`k` neutrons in the first generation, +which in turn produce :math:`k^2` neutrons in the second generation, and so on. +However, the above picture cannot be taken literally, because the neutrons +born from the external source will have a different importance to neutron +production than will fission neutrons, and we have the following alternative +picture for :math:`M` [Kobayashi]_: + +.. math:: + :label: subcritical_k_factors + + \begin{align*} + M-1 &= k_q + k_q k_s + k_q k_s^2 + ... \\ + &= \frac{k_q}{1-k_s} \\ + \implies \frac{k}{1-k} &= \frac{k_q}{1-k_s}\\ + k &= \frac{k_q}{1 - k_s + k_q} + \end{align*} +Where :math:`k_q` is the multiplication factor of source neutrons, and :math:`k_s` +is the multiplication factor of fission neutrons, which together define an overall +subcritical multiplication factor :math:`k`. From the above it is clear that +:math:`k < 1 \iff k_s < 1`, and :math:`k_s <1` for :math:`k_{eff}<1`, so a +subcritical system will correctly have :math:`k < 1`. It is however not the case +that :math:`k_s = k_{eff}`, because the angular flux of fission neutrons is not +necessarily the fundamental mode calculated in eigenvalue mode, nor is it +true that :math:`k = k_{eff}`. In fact, for deeply subcritical systems, +:math:`k_{eff}` generally underestimates :math:`k` [Forget]_. In addition, we may +have :math:`k_q>1` despite :math:`k_s<1`, and in the limit, we may have :math:`k` +arbitrarily close to 1: an arbitrarily multiplying system that is still +subcritical with respect to fission neutrons. In fact, this is a primary design +consideration of ADS's, where :math:`k_s` is fixed :math:`<1` to ensure +subcritciality, while :math:`k_q` is maximized to achieve optimal multiplication. +It is therefore necessary to perform fixed source simulations to accurately determine +subcritical multiplication and the flux distribution in ADS's. + +.. _methods_subcritical-multiplication-methods: +---------------------------------- +Subcritical Multiplication Methods +---------------------------------- +The most straightforward method for performing a subcritical multiplication calculation is to perform a fixed source +simulation. For near critical systems, however, excessive secondary particle production can result in memory issues if +not properly managed. In very near critical systems, performing an eigenvalue simulation may be sufficient, +as the flux distribution will be close to the fundamental mode, and :math:`k` will be close to :math:`k_{eff}`, just note that +tally results must be multiplied by :math:`1/(1-k)`. For more deeply subcritical systems, there are several tailor-made methods +for subcritical multiplication calculations that limit excessive secondary particle production while equivalent results to fixed source simulations. +The method currently implemented in OpenMC is that of [Forget]_, which uses a modified power iteration to converge a source +distribution taking into account the external source and all generations of fission neutrons, which is then sampled from in generations +to accumulate statistics about the system. This method resembles the standard eigenvalue simulation (see :ref:`methods_eigenvalue`) in every way except +that tallies have a pre-defined normalization by :math:`1/(1-k)` (rather than normalized to an arbitrary power distribution), the :math:`k`-eigenvalue is the +subcritical multiplication factor (rather than the fundamental mode eigenvalue :math:`k_{eff}`), and the results are determined by the prescribed +source distribution. + + +.. [Bowman] Bowman, Charles D. "Accelerator-driven systems for nuclear waste + transmutation." *Annual Review of Nuclear and Particle Science* 48.1 (1998): + 505-556. +.. [Kobayashi] Kobayashi, Keisuke, and Kenji Nishihara. "Definition of + subcriticality using the importance function for the production of fission + neutrons." *Nuclear science and engineering* 136.2 (2000): 272-281. +.. [Forget] Forget, Benoit. "An Efficient Subcritical Multiplication Mode for + Monte Carlo Solvers." *Nuclear Science and Engineering* (2025): 1-11. \ No newline at end of file diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 5a04fedd70a..8469623d93e 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -33,6 +33,12 @@ be specified: Runs a fixed-source calculation with a specified external source, specified in the :attr:`Settings.source` attribute. +'subcritical multiplication' + Runs a subcritical multiplication calculation (as described in :ref:`methods_subcritical-multiplication`). + In this mode, the :attr:`Settings.source` specifies an external source distribution that is used for all + generations. Note that tallies must be manually normalized by :math:`1/(1-k)` for direct comparison with fixed + source results. + 'volume' Runs a stochastic volume calculation. diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 8185214fced..b3649d634ad 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -356,6 +356,7 @@ enum class RunMode { UNSET, // default value, OpenMC throws error if left to this FIXED_SOURCE, EIGENVALUE, + SUBCRITICAL_MULTIPLICATION, PLOTTING, PARTICLE, VOLUME diff --git a/include/openmc/settings.h b/include/openmc/settings.h index b369c99fef8..a1658c65d2b 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -198,6 +198,8 @@ extern "C" int verbosity; //!< How verbose to make output extern double weight_cutoff; //!< Weight cutoff for Russian roulette extern double weight_survive; //!< Survival weight after Russian roulette +bool eigenvalue_like(); + } // namespace settings //============================================================================== diff --git a/openmc/settings.py b/openmc/settings.py index 3cf662e311e..a7217eadefc 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -22,6 +22,7 @@ class RunMode(Enum): EIGENVALUE = 'eigenvalue' + SUBCRITICAL_MULTIPLICATION = 'subcritical multiplication' FIXED_SOURCE = 'fixed source' PLOT = 'plot' VOLUME = 'volume' @@ -238,7 +239,7 @@ class Settings: The 'nuclides' list indicates what nuclides the method should be applied to. In its absence, the method will be applied to all nuclides with 0 K elastic scattering data present. - run_mode : {'eigenvalue', 'fixed source', 'plot', 'volume', 'particle restart'} + run_mode : {'eigenvalue', 'fixed source', 'subcritical multiplication', 'plot', 'volume', 'particle restart'} The type of calculation to perform (default is 'eigenvalue') seed : int Seed for the linear congruential pseudorandom number generator diff --git a/openmc/statepoint.py b/openmc/statepoint.py index a10ec3a839d..70b6ce85181 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -212,7 +212,7 @@ def date_and_time(self): @property def entropy(self): - if self.run_mode == 'eigenvalue': + if 'entropy' in self._f: return self._f['entropy'][()] else: return None @@ -233,7 +233,7 @@ def filters(self): @property def generations_per_batch(self): - if self.run_mode == 'eigenvalue': + if 'generations_per_batch' in self._f: return self._f['generations_per_batch'][()] else: return None @@ -268,14 +268,14 @@ def k_cmfd(self): @property def k_generation(self): - if self.run_mode == 'eigenvalue': + if 'k_generation' in self._f: return self._f['k_generation'][()] else: return None @property def keff(self): - if self.run_mode == 'eigenvalue': + if 'k_combined' in self._f: return ufloat(*self._f['k_combined'][()]) else: return None @@ -290,21 +290,21 @@ def k_combined(self): @property def k_col_abs(self): - if self.run_mode == 'eigenvalue': + if 'k_col_abs' in self._f: return self._f['k_col_abs'][()] else: return None @property def k_col_tra(self): - if self.run_mode == 'eigenvalue': + if 'k_col_tra' in self._f: return self._f['k_col_tra'][()] else: return None @property def k_abs_tra(self): - if self.run_mode == 'eigenvalue': + if 'k_abs_tra' in self._f: return self._f['k_abs_tra'][()] else: return None @@ -329,7 +329,7 @@ def n_batches(self): @property def n_inactive(self): - if self.run_mode == 'eigenvalue': + if 'n_inactive' in self._f: return self._f['n_inactive'][()] else: return None diff --git a/src/main.cpp b/src/main.cpp index 88251ac7232..5f513b32b0d 100644 --- a/src/main.cpp +++ b/src/main.cpp @@ -42,6 +42,13 @@ int main(int argc, char* argv[]) break; } break; + case RunMode::SUBCRITICAL_MULTIPLICATION: + switch (settings::solver_type) { + case SolverType::MONTE_CARLO: + err = openmc_run(); + break; + } + break; case RunMode::PLOTTING: err = openmc_plot_geometry(); break; diff --git a/src/output.cpp b/src/output.cpp index ae2daaffc14..b33b008c366 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -444,11 +444,11 @@ void print_runtime() show_time("Collisions", time_event_collision.elapsed(), 2); show_time("Particle death", time_event_death.elapsed(), 2); } - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { show_time("Time in inactive batches", time_inactive.elapsed(), 1); } show_time("Time in active batches", time_active.elapsed(), 1); - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { show_time("Time synchronizing fission bank", time_bank.elapsed(), 1); show_time("Sampling source sites", time_bank_sample.elapsed(), 2); show_time("SEND/RECV source sites", time_bank_sendrecv.elapsed(), 2); @@ -535,7 +535,7 @@ void print_results() const auto& gt = simulation::global_tallies; double mean, stdev; if (n > 1) { - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { std::tie(mean, stdev) = mean_stdev(>(GlobalTally::K_COLLISION, 0), n); fmt::print(" k-effective (Collision) = {:.5f} +/- {:.5f}\n", mean, t_n1 * stdev); @@ -560,7 +560,7 @@ void print_results() warning("Could not compute uncertainties -- only one " "active batch simulated!"); - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { fmt::print(" k-effective (Collision) = {:.5f}\n", gt(GlobalTally::K_COLLISION, TallyResult::SUM) / n); fmt::print(" k-effective (Track-length) = {:.5f}\n", diff --git a/src/particle.cpp b/src/particle.cpp index a035e76b917..b2f55f4d1f6 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -279,7 +279,7 @@ void Particle::event_advance() } // Score track-length estimate of k-eff - if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) { + if (settings::eigenvalue_like() && type().is_neutron()) { keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission; } @@ -341,7 +341,7 @@ void Particle::event_collide() { // Score collision estimate of keff - if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) { + if (settings::eigenvalue_like() && type().is_neutron()) { keff_tally_collision() += wgt() * macro_xs().nu_fission / macro_xs().total; } @@ -518,7 +518,7 @@ void Particle::event_death() // Record the number of progeny created by this particle. // This data will be used to efficiently sort the fission bank. - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { int64_t offset = id() - 1 - simulation::work_index[mpi::rank]; simulation::progeny_per_particle[offset] = n_progeny(); } @@ -840,7 +840,7 @@ void Particle::write_restart() const write_dataset(file_id, "type", type().pdg_number()); int64_t i = current_work(); - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { // take source data from primary bank for eigenvalue simulation write_dataset(file_id, "weight", simulation::source_bank[i - 1].wgt); write_dataset(file_id, "energy", simulation::source_bank[i - 1].E); diff --git a/src/particle_restart.cpp b/src/particle_restart.cpp index f02fcb94b55..d5c39206594 100644 --- a/src/particle_restart.cpp +++ b/src/particle_restart.cpp @@ -51,6 +51,8 @@ void read_particle_restart(Particle& p, RunMode& previous_run_mode) previous_run_mode = RunMode::EIGENVALUE; } else if (mode == "fixed source") { previous_run_mode = RunMode::FIXED_SOURCE; + } else if (mode == "subcritical multiplication") { + previous_run_mode = RunMode::SUBCRITICAL_MULTIPLICATION; } read_dataset(file_id, "id", p.id()); int type; @@ -110,6 +112,7 @@ void run_particle_restart() int64_t particle_seed; switch (previous_run_mode) { case RunMode::EIGENVALUE: + case RunMode::SUBCRITICAL_MULTIPLICATION: case RunMode::FIXED_SOURCE: particle_seed = (simulation::total_gen + overall_generation() - 1) * settings::n_particles + diff --git a/src/physics.cpp b/src/physics.cpp index 1f72e4fac07..75a4943fc70 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -110,7 +110,7 @@ void sample_neutron_reaction(Particle& p) if (nuc->fissionable_ && p.neutron_xs(i_nuclide).fission > 0.0) { auto& rx = sample_fission(i_nuclide, p); - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { create_fission_sites(p, i_nuclide, rx); } else if (settings::run_mode == RunMode::FIXED_SOURCE && settings::create_fission_neutrons) { @@ -201,7 +201,7 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) // Determine whether to place fission sites into the shared fission bank // or the secondary particle bank. - bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE); + bool use_fission_bank = (settings::eigenvalue_like()); // Counter for the number of fission sites successfully stored to the shared // fission bank or the secondary particle bank @@ -649,7 +649,7 @@ void absorption(Particle& p, int i_nuclide) p.wgt() -= wgt_absorb; // Score implicit absorption estimate of keff - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { p.keff_tally_absorption() += wgt_absorb * p.neutron_xs(i_nuclide).nu_fission / p.neutron_xs(i_nuclide).absorption; @@ -659,7 +659,7 @@ void absorption(Particle& p, int i_nuclide) if (p.neutron_xs(i_nuclide).absorption > prn(p.current_seed()) * p.neutron_xs(i_nuclide).total) { // Score absorption estimate of keff - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { p.keff_tally_absorption() += p.wgt() * p.neutron_xs(i_nuclide).nu_fission / p.neutron_xs(i_nuclide).absorption; @@ -1212,7 +1212,7 @@ void sample_secondary_photons(Particle& p, int i_nuclide) // Stedry, "Self-consistent energy normalization for quasistatic reactor // calculations", Proc. PHYSOR, Cambridge, UK, Mar 29-Apr 2, 2020. double wgt = photon_wgt; - if (settings::run_mode == RunMode::EIGENVALUE && !is_fission(rx->mt_)) { + if (settings::eigenvalue_like() && !is_fission(rx->mt_)) { wgt *= simulation::keff; } diff --git a/src/physics_mg.cpp b/src/physics_mg.cpp index 58a21a67895..83d491e7108 100644 --- a/src/physics_mg.cpp +++ b/src/physics_mg.cpp @@ -49,7 +49,7 @@ void sample_reaction(Particle& p) // absorption (including fission) if (model::materials[p.material()]->fissionable()) { - if (settings::run_mode == RunMode::EIGENVALUE || + if (settings::eigenvalue_like() || (settings::run_mode == RunMode::FIXED_SOURCE && settings::create_fission_neutrons)) { create_fission_sites(p); @@ -127,7 +127,7 @@ void create_fission_sites(Particle& p) // Determine whether to place fission sites into the shared fission bank // or the secondary particle bank. - bool use_fission_bank = (settings::run_mode == RunMode::EIGENVALUE); + bool use_fission_bank = (settings::eigenvalue_like()); // Counter for the number of fission sites successfully stored to the shared // fission bank or the secondary particle bank diff --git a/src/settings.cpp b/src/settings.cpp index 5b472468fc4..b3d64bbd4f7 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -149,6 +149,12 @@ int verbosity {-1}; double weight_cutoff {0.25}; double weight_survive {1.0}; +bool eigenvalue_like() +{ + return run_mode == RunMode::EIGENVALUE || + run_mode == RunMode::SUBCRITICAL_MULTIPLICATION; +} + } // namespace settings //============================================================================== @@ -208,8 +214,7 @@ void get_run_parameters(pugi::xml_node node_base) } // Get number of inactive batches - if (run_mode == RunMode::EIGENVALUE || - solver_type == SolverType::RANDOM_RAY) { + if (settings::eigenvalue_like() || solver_type == SolverType::RANDOM_RAY) { if (check_for_node(node_base, "inactive")) { n_inactive = std::stoi(get_node_value(node_base, "inactive")); } @@ -490,6 +495,8 @@ void read_settings_xml(pugi::xml_node root) std::string temp_str = get_node_value(root, "run_mode", true, true); if (temp_str == "eigenvalue") { run_mode = RunMode::EIGENVALUE; + } else if (temp_str == "subcritical multiplication") { + run_mode = RunMode::SUBCRITICAL_MULTIPLICATION; } else if (temp_str == "fixed source") { run_mode = RunMode::FIXED_SOURCE; } else if (temp_str == "plot") { @@ -525,12 +532,16 @@ void read_settings_xml(pugi::xml_node root) // Check solver type if (check_for_node(root, "random_ray")) { solver_type = SolverType::RANDOM_RAY; + if (run_mode == RunMode::SUBCRITICAL_MULTIPLICATION) { + fatal_error("random ray solver not currently supported in subcritical " + "multiplication mode"); + } if (run_CE) fatal_error("multi-group energy mode must be specified in settings XML " "when using the random ray solver."); } - if (run_mode == RunMode::EIGENVALUE || run_mode == RunMode::FIXED_SOURCE) { + if (settings::eigenvalue_like() || run_mode == RunMode::FIXED_SOURCE) { // Read run parameters get_run_parameters(node_mode); diff --git a/src/simulation.cpp b/src/simulation.cpp index 18e40a8bc73..764b9e68880 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -130,8 +130,8 @@ int openmc_simulation_init() load_state_point(); write_message("Resuming simulation...", 6); } else { - // Only initialize primary source bank for eigenvalue simulations - if (settings::run_mode == RunMode::EIGENVALUE && + // Only initialize primary source bank for eigenvalue-like simulations + if (settings::eigenvalue_like() && settings::solver_type == SolverType::MONTE_CARLO) { initialize_source(); } @@ -150,6 +150,9 @@ int openmc_simulation_init() header("K EIGENVALUE SIMULATION", 3); } else if (settings::solver_type == SolverType::RANDOM_RAY) { header("K EIGENVALUE SIMULATION (RANDOM RAY SOLVER)", 3); + } else if (settings::run_mode == RunMode::SUBCRITICAL_MULTIPLICATION) { + header( + "FIXED SOURCE (SUBCRITICAL MULTIPLICATION) TRANSPORT SIMULATION", 3); } if (settings::verbosity >= 7) print_columns(); @@ -303,6 +306,7 @@ int current_batch; int current_gen; bool initialized {false}; double keff {1.0}; +double kold {1.0}; double keff_std; double k_col_abs {0.0}; double k_col_tra {0.0}; @@ -331,7 +335,7 @@ vector work_index; void allocate_banks() { - if (settings::run_mode == RunMode::EIGENVALUE && + if (settings::eigenvalue_like() && settings::solver_type == SolverType::MONTE_CARLO) { // Allocate source bank simulation::source_bank.resize(simulation::work_per_rank); @@ -435,7 +439,7 @@ void finalize_batch() !settings::cmfd_run) { if (contains(settings::sourcepoint_batch, simulation::current_batch) && settings::source_write && !settings::source_separate) { - bool b = (settings::run_mode == RunMode::EIGENVALUE); + bool b = (settings::eigenvalue_like()); openmc_statepoint_write(nullptr, &b); } else { bool b = false; @@ -443,7 +447,7 @@ void finalize_batch() } } - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { // Write out a separate source point if it's been specified for this batch if (contains(settings::sourcepoint_batch, simulation::current_batch) && settings::source_write && settings::source_separate) { @@ -505,7 +509,7 @@ void finalize_batch() void initialize_generation() { - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { // Clear out the fission bank simulation::fission_bank.resize(0); @@ -524,7 +528,7 @@ void finalize_generation() auto& gt = simulation::global_tallies; // Update global tallies with the accumulation variables - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { gt(GlobalTally::K_COLLISION, TallyResult::VALUE) += global_tally_collision; gt(GlobalTally::K_ABSORPTION, TallyResult::VALUE) += global_tally_absorption; @@ -534,14 +538,14 @@ void finalize_generation() gt(GlobalTally::LEAKAGE, TallyResult::VALUE) += global_tally_leakage; // reset tallies - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { global_tally_collision = 0.0; global_tally_absorption = 0.0; global_tally_tracklength = 0.0; } global_tally_leakage = 0.0; - if (settings::run_mode == RunMode::EIGENVALUE && + if (settings::eigenvalue_like() && settings::solver_type == SolverType::MONTE_CARLO) { // If using shared memory, stable sort the fission bank (by parent IDs) // so as to allow for reproducibility regardless of which order particles @@ -552,7 +556,7 @@ void finalize_generation() synchronize_bank(); } - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { // Calculate shannon entropy if (settings::entropy_on && @@ -563,6 +567,8 @@ void finalize_generation() calculate_generation_keff(); calculate_average_keff(); + simulation::kold = simulation::keff; + // Write generation output if (mpi::master && settings::verbosity >= 7) { print_generation(); @@ -576,15 +582,38 @@ void initialize_history(Particle& p, int64_t index_source) if (settings::run_mode == RunMode::EIGENVALUE) { // set defaults for eigenvalue simulations from primary bank p.from_source(&simulation::source_bank[index_source - 1]); - } else if (settings::run_mode == RunMode::FIXED_SOURCE) { + } else { // initialize random number seed int64_t id = (simulation::total_gen + overall_generation() - 1) * settings::n_particles + simulation::work_index[mpi::rank] + index_source; uint64_t seed = init_seed(id, STREAM_SOURCE); - // sample from external source distribution or custom library then set - auto site = sample_external_source(&seed); - p.from_source(&site); + if (settings::run_mode == RunMode::SUBCRITICAL_MULTIPLICATION) { + double rnd = prn(&seed); + double k_avg = 0.0; + int n = simulation::k_generation.size(); + if (n >= 2) { + // Average the last two values + double k_last = simulation::k_generation[n - 1]; + double k_prev = simulation::k_generation[n - 2]; + k_avg = (k_last + k_prev) / 2.0; + } else if (n == 1) { + // Only one generation exists, use it directly + k_avg = simulation::k_generation[0]; + } + if (rnd < k_avg) { + // sample from fission source bank + p.from_source(&simulation::source_bank[index_source - 1]); + } else { + // sample from external source + auto site = sample_external_source(&seed); + p.from_source(&site); + } + } else if (settings::run_mode == RunMode::FIXED_SOURCE) { + // sample from external source distribution or custom library then set + auto site = sample_external_source(&seed); + p.from_source(&site); + } } p.current_work() = index_source; diff --git a/src/state_point.cpp b/src/state_point.cpp index 455299529b3..af07807d70c 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -104,6 +104,9 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source) case RunMode::EIGENVALUE: write_dataset(file_id, "run_mode", "eigenvalue"); break; + case RunMode::SUBCRITICAL_MULTIPLICATION: + write_dataset(file_id, "run_mode", "subcritical multiplication"); + break; default: break; } @@ -118,7 +121,7 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source) write_attribute(file_id, "source_present", write_source_); // Write out information for eigenvalue run - if (settings::run_mode == RunMode::EIGENVALUE) + if (settings::eigenvalue_like()) write_eigenvalue_hdf5(file_id); hid_t tallies_group = create_group(file_id, "tallies"); @@ -312,11 +315,11 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source) write_dataset(runtime_group, "simulation", time_inactive.elapsed() + time_active.elapsed()); write_dataset(runtime_group, "transport", time_transport.elapsed()); - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { write_dataset(runtime_group, "inactive batches", time_inactive.elapsed()); } write_dataset(runtime_group, "active batches", time_active.elapsed()); - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { write_dataset( runtime_group, "synchronizing fission bank", time_bank.elapsed()); write_dataset( @@ -433,6 +436,8 @@ extern "C" int openmc_statepoint_load(const char* filename) settings::run_mode = RunMode::FIXED_SOURCE; } else if (word == "eigenvalue") { settings::run_mode = RunMode::EIGENVALUE; + } else if (word == "subcritical multiplication") { + settings::run_mode = RunMode::SUBCRITICAL_MULTIPLICATION; } read_attribute(file_id, "photon_transport", settings::photon_transport); read_dataset(file_id, "n_particles", settings::n_particles); @@ -459,7 +464,7 @@ extern "C" int openmc_statepoint_load(const char* filename) read_attribute(file_id, "source_present", source_present); // Read information specific to eigenvalue run - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { read_dataset(file_id, "n_inactive", temp); read_eigenvalue_hdf5(file_id); @@ -478,7 +483,7 @@ extern "C" int openmc_statepoint_load(const char* filename) // Set k_sum, keff, and current_batch based on whether restart file is part // of active cycle or inactive cycle - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { restart_set_keff(); } @@ -529,7 +534,7 @@ extern "C" int openmc_statepoint_load(const char* filename) } // Read source if in eigenvalue mode - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { // Check if source was written out separately if (!source_present) { diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 8cdf3a9050c..e6ea1939f6a 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -204,7 +204,7 @@ Tally::Tally(pugi::xml_node node) // Check for errors if (has_ifp_score) { - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { if (settings::ifp_n_generation < 0) { settings::ifp_n_generation = DEFAULT_IFP_N_GENERATION; warning(fmt::format( @@ -1074,7 +1074,7 @@ void accumulate_tallies() if (mpi::master || !settings::reduce_tallies) { auto& gt = simulation::global_tallies; - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { if (simulation::current_batch > settings::n_inactive) { // Accumulate products of different estimators of k double k_col = gt(GlobalTally::K_COLLISION, TallyResult::VALUE) / diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 51b5d9ffcc5..d29f366aea2 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -264,7 +264,7 @@ double get_nuclide_neutron_heating( if (kerma == 0.0) return 0.0; - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::eigenvalue_like()) { // Determine kerma for fission as (EFR + EB)*sigma_f double kerma_fission = nuc.fragments_ diff --git a/src/tallies/trigger.cpp b/src/tallies/trigger.cpp index f1f83e2982c..13ee434c75a 100644 --- a/src/tallies/trigger.cpp +++ b/src/tallies/trigger.cpp @@ -126,7 +126,7 @@ void check_tally_triggers(double& ratio, int& tally_id, int& score) double check_keff_trigger() { - if (settings::run_mode != RunMode::EIGENVALUE) + if (!settings::eigenvalue_like()) return 0.0; double k_combined[2]; diff --git a/tests/regression_tests/subcritical_multiplication/__init__.py b/tests/regression_tests/subcritical_multiplication/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/subcritical_multiplication/inputs_true.dat b/tests/regression_tests/subcritical_multiplication/inputs_true.dat new file mode 100644 index 00000000000..a86982b1e55 --- /dev/null +++ b/tests/regression_tests/subcritical_multiplication/inputs_true.dat @@ -0,0 +1,36 @@ + + + + mgxs.h5 + + + + + + + + + + + + subcritical multiplication + 100000 + 20 + 10 + + + 0 -1000 -1000 10 1000 1000 + + + 10000000.0 1.0 + + + + false + + multi-group + + false + + + diff --git a/tests/regression_tests/subcritical_multiplication/results_true.dat b/tests/regression_tests/subcritical_multiplication/results_true.dat new file mode 100644 index 00000000000..2a5e7474625 --- /dev/null +++ b/tests/regression_tests/subcritical_multiplication/results_true.dat @@ -0,0 +1,23 @@ +k +8.878177E-01 1.072332E-03 +k_generation: +1.380360E+00 +8.242283E-01 +8.322178E-01 +9.234869E-01 +8.951838E-01 +8.756053E-01 +8.894170E-01 +8.920616E-01 +8.891805E-01 +8.826478E-01 +8.909997E-01 +8.965785E-01 +8.917731E-01 +8.823861E-01 +8.912402E-01 +8.871419E-01 +8.850541E-01 +8.878719E-01 +8.930822E-01 +8.874183E-01 diff --git a/tests/regression_tests/subcritical_multiplication/test.py b/tests/regression_tests/subcritical_multiplication/test.py new file mode 100644 index 00000000000..294bad0a967 --- /dev/null +++ b/tests/regression_tests/subcritical_multiplication/test.py @@ -0,0 +1,95 @@ +"""This test is based on a simple 4-group slab model from +"MCNP Calculations of Subcritical Fixed and Fission Multiplication Factors", +LA-UR-10-00141 which can be found at https://mcnp.lanl.gov/pdf_files/TechReport_2010_LANL_LA-UR-10-00141_KiedrowskiBrown.pdf +""" +import openmc +from openmc.stats import delta_function +import numpy as np +import pytest +from openmc.examples import slab_mg +import os +import glob + +from tests.testing_harness import PyAPITestHarness + + +class MGXSTestHarness(PyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = 'mgxs.h5' + if os.path.exists(f): + os.remove(f) + def _get_results(self, hash_output=False): + outstr = super()._get_results(hash_output=hash_output) + # Read the statepoint file. + statepoint = glob.glob(self._sp_name)[0] + with openmc.StatePoint(statepoint) as sp: + # Write out k + outstr += 'k\n' + form = '{0:12.6E} {1:12.6E}\n' + k = sp.keff + outstr += form.format(k.n, k.s) + + # Write out k_generation + outstr += 'k_generation:\n' + form = '{0:12.6E}\n' + k_gen = sp.k_generation + for kg in k_gen: + outstr += form.format(kg) + + return outstr + + +@pytest.fixture() +def slab_model(): + model = slab_mg(mgxslib_name='mgxs.h5') + right_boundary = model.geometry.get_all_surfaces()[2] + right_boundary.coefficients['x0'] = 10.0 + + cell = model.geometry.get_all_cells()[1] + + mat = model.geometry.get_all_materials()[1] + mat.set_density('macro', 0.01) + ########################################################################### + # Create multigroup data + + # Instantiate the energy group data + ebins = np.geomspace(1e-5, 20.0e6, 5) + groups = openmc.mgxs.EnergyGroups(group_edges=ebins) + + nusigma_f = np.array([9.6,5.4,5.2,2.5]) + sigma_s = np.array([[0.5,0.5,0.5,0.5], + [0.0,1.0,0.5,0.5], + [0.0,0.0,1.5,0.5], + [0.0,0.0,0.0,2.0]]) + sigma_t = np.array([5.0,5.0,5.0,5.0]) + sigma_a = sigma_t - sigma_s.sum(axis=1) + chi = np.array([0.0,0.2,0.8,0.0]) + mat_data = openmc.XSdata('mat_1', groups) + mat_data.order = 0 + mat_data.set_total(sigma_t) + mat_data.set_nu_fission(nusigma_f) + mat_data.set_chi(chi) + mat_data.set_absorption(sigma_a) + mat_data.set_scatter_matrix(sigma_s[...,np.newaxis]) + + mg_cross_sections_file = openmc.MGXSLibrary(groups) + mg_cross_sections_file.add_xsdata(mat_data) + mg_cross_sections_file.export_to_hdf5() + + # Settings + model.settings.particles = 100000 + model.settings.inactive = 10 + model.settings.batches = 20 + model.settings.run_mode = 'subcritical multiplication' + + space = openmc.stats.Box([0,-1000,-1000],[10,1000,1000]) + model.settings.source = openmc.IndependentSource( + space=space, energy = delta_function(10e6)) + + return model + + +def test_multiplication(slab_model): + harness = MGXSTestHarness("statepoint.20.h5", model=slab_model) + harness.main()