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..5cca628615b --- /dev/null +++ b/docs/source/methods/subcritical_multiplication.rst @@ -0,0 +1,121 @@ +.. _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. + +Source Effectiveness +-------------------- +In addition, one may quantify the relative effectiveness of a source by the amount of +multiplication it induces relative to the fundamental mode of the system. This is given +by the source effectiveness factor :math:`g^*` (or :math:`\phi^*`), defined as [Ye]_: + +.. math:: + M = \frac{g^*}{1 - k_{eff}} + +which corrects for the error in the multiplicity one would get from using the +fundamental mode multiplication factor. As :math:`k_{eff} \to 1`, :math:`g^* \to 1`, +since as the system approaches criticality, the importance of the source diminishes +relative to subsequent fisison neutrons. + +.. _methods_subcritical-multiplication-estimating: + +----------------------------------------------------------------------- +Estimating :math:`k`, :math:`k_q`, and :math:`k_s` in Fixed Source Mode +----------------------------------------------------------------------- +The total multiplication factor :math:`k` can be estimated through :eq:`integral_multiplicity`. +The total fission production can be tallied and estiamted using standard collision, absorption, +and track-length estimators over a neutron history, giving :math:`M-1`, which can be used to +compute :math:`k`. + +To estimate :math:`k_q`, we may use its interpretation as the multiplication +factor of source neutrons. For a given source history, we may tally the neutron production +estimators, and simply stop before simulating any of the secondary fission neutrons. This gives +an estimate of the neutron production due to source neutrons alone, which can be used to compute +:math:`k_q`. :math:`k_s` can then be computed from :math:`k` and :math:`k_q` using +:eq:`subcritical_k_factors`. + +.. [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. + +.. [Ye] Ye, Bin, Chao-Wen Yang, and Chun Zheng. "Measurement of k eff by delayed + neutron multiplication in subcritical systems." Nuclear Science and + Techniques 29.2 (2018): 29. \ No newline at end of file diff --git a/docs/source/usersguide/settings.rst b/docs/source/usersguide/settings.rst index 5a04fedd70a..d7aa7a96721 100644 --- a/docs/source/usersguide/settings.rst +++ b/docs/source/usersguide/settings.rst @@ -851,3 +851,39 @@ the inputs were not modified before the restart run), no particles will be transported and OpenMC will exit immediately. .. note:: A statepoint file must match the input model to be successfully used in a restart simulation. + +---------------------------------------------- +Calculating Subcritical Multiplication Factors +---------------------------------------------- +In addition to the standard effective multiplication factor, :math:`k_{eff}`, calculated using eigenvalue +mode, OpenMC can also calculate subcritical multiplication factors: :math:`k`, :math:`k_q`, and :math:`k_s` +in fixed source mode or eigenvalue mode, but note these values are only strictly well-defined for fixed source simulations. +Please see :ref:`methods_subcritical-multiplication` for theoretical details on +subcritical multiplication factors. + +Fixed Source Mode +----------------- +To enable the calculation of subcritical multiplication factors in +fixed source mode, set + +.. code-block:: python + + settings.calculate_subcritical_k = True + +this will enable printing of :math:`k`, :math:`k_q`, and :math:`k_s` both generation-wise and cumulative +averages throughout the simulation analogous to eigenvalue mode. The generation statistics as well as +combined estimates for :math:`k`, :math:`k_q`, and :math:`k_s` will also be stored in the statepoint file. +To print out all of the subcritical multiplication factors during the transport calculation, set + +.. code-block:: python + + settings.print_all_k_factors = True + +Eigenvalue Mode +--------------- +These quantities may also be computed using the eigensolver if `settings.source` is set to the source +distribution for your corresponding fixed source problem. These quantities are computed by default, and there's not need +to set `settings.calculate_subcritical_k` to True. The generation statistics as well as combined estimates for :math:`k`, +:math:`k_q`, and :math:`k_s` will be stored in the statepoint file. In addition, the source effectiveness factor will also +be computed, since :math:`k_{eff}` is available. Please interpret these values with caution, as they are only strictly +well-defined when :math:`k_{eff} < 1`. \ No newline at end of file diff --git a/include/openmc/constants.h b/include/openmc/constants.h index 8185214fced..dc482412a2b 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -330,8 +330,15 @@ enum TallyScore { }; // Global tally parameters -constexpr int N_GLOBAL_TALLIES {4}; -enum class GlobalTally { K_COLLISION, K_ABSORPTION, K_TRACKLENGTH, LEAKAGE }; +enum class GlobalTally { + K_COLLISION, + K_ABSORPTION, + K_TRACKLENGTH, + LEAKAGE, + K_TRACKLENGTH_SQ, // for calculation of stddev for generational k + SIZE +}; +constexpr int N_GLOBAL_TALLIES {static_cast(GlobalTally::SIZE)}; // Miscellaneous constexpr int C_NONE {-1}; @@ -361,6 +368,9 @@ enum class RunMode { VOLUME }; +// Eigenvalue calculation parameters +enum class KeffType { k, kq, ks }; + enum class SolverType { MONTE_CARLO, RANDOM_RAY }; enum class RandomRayVolumeEstimator { NAIVE, SIMULATION_AVERAGED, HYBRID }; diff --git a/include/openmc/eigenvalue.h b/include/openmc/eigenvalue.h index b456fee21eb..2cb52aeeab1 100644 --- a/include/openmc/eigenvalue.h +++ b/include/openmc/eigenvalue.h @@ -21,8 +21,12 @@ namespace openmc { namespace simulation { -extern double keff_generation; //!< Single-generation k on each processor +extern array + keff_generation; //!< Single-generation k on each processor +extern array + kq_generation_val; //!< Single-generation kq on each processor extern array k_sum; //!< Used to reduce sum and sum_sq +extern array kq_sum; extern vector entropy; //!< Shannon entropy at each generation extern xt::xtensor source_frac; //!< Source fraction for UFS @@ -32,15 +36,25 @@ extern xt::xtensor source_frac; //!< Source fraction for UFS // Non-member functions //============================================================================== -//! Collect/normalize the tracklength keff from each process +//! Collect/normalize the tracklength keff from each process (keff) void calculate_generation_keff(); +//! Collect/normalize the tracklength keff from each process +//! +//! Depending on KeffType, this function will calculate either the kq or ks +//! generation keff +void calculate_generation_keff(KeffType type); + +std::pair convert_m_to_k(double m, double m_std); +std::pair convert_k_to_m(double k, double k_std); + //! Calculate mean/standard deviation of keff during active generations //! //! This function sets the global variables keff and keff_std which represent //! the mean and standard deviation of the mean of k-effective over active //! generations. It also broadcasts the value from the master process. void calculate_average_keff(); +void calculate_average_keff(KeffType type); //! Calculates a minimum variance estimate of k-effective //! @@ -55,6 +69,11 @@ void calculate_average_keff(); //! \param[out] k_combined Estimate of k-effective and its standard deviation //! \return Error status extern "C" int openmc_get_keff(double* k_combined); +extern "C" int openmc_get_kq(double* kq_combined); +int get_combined_k_from_tallies(double* k_combined, + std::array& k_combined_weights, double keff, double keff_std, + xt::xtensor_fixed> global_tallies, + double k_col_abs, double k_col_tra, double k_abs_tra); //! Sample/redistribute source sites from accumulated fission sites void synchronize_bank(); diff --git a/include/openmc/output.h b/include/openmc/output.h index 940ea78ceb9..8aa9b5ed0f6 100644 --- a/include/openmc/output.h +++ b/include/openmc/output.h @@ -5,6 +5,7 @@ #define OPENMC_OUTPUT_H #include +#include #include "openmc/particle.h" @@ -48,6 +49,8 @@ void print_columns(); //! Display information about a generation of neutrons void print_generation(); +void print_generation_values(int idx, int n, std::string batch_and_gen, + array k_gen, double k, double k_std); //! Display time elapsed for various stages of a run void print_runtime(); diff --git a/include/openmc/settings.h b/include/openmc/settings.h index b369c99fef8..4987f666183 100644 --- a/include/openmc/settings.h +++ b/include/openmc/settings.h @@ -169,7 +169,10 @@ extern double source_rejection_fraction; //!< Minimum fraction of source sites //!< that must be accepted extern double free_gas_threshold; //!< Threshold multiplier for free gas //!< scattering treatment - +extern bool calculate_subcritical_k; //!< Calculate subcritical k in fixed + //!< source mode +extern bool + print_all_k_factors; //!< Print k-effective for all batches and generations extern int max_history_splits; //!< maximum number of particle splits for weight windows extern int max_secondaries; //!< maximum number of secondaries in the bank diff --git a/include/openmc/simulation.h b/include/openmc/simulation.h index 9a6cf1b2131..9ec7ff3262f 100644 --- a/include/openmc/simulation.h +++ b/include/openmc/simulation.h @@ -4,11 +4,13 @@ #ifndef OPENMC_SIMULATION_H #define OPENMC_SIMULATION_H +#include "openmc/array.h" #include "openmc/mesh.h" #include "openmc/particle.h" #include "openmc/vector.h" #include +#include namespace openmc { @@ -28,13 +30,30 @@ extern "C" int current_gen; //!< current fission generation extern "C" bool initialized; //!< has simulation been initialized? extern "C" double keff; //!< average k over batches extern "C" double keff_std; //!< standard deviation of average k +extern "C" double + k; //!< average k over batches, used for subcritical multiplication problems +extern "C" double k_std; //!< standard deviation of average k, used for + //!< subcritical multiplication problems +extern "C" double + kq; //!< average kq over batches, used for subcritical multiplication problems +extern "C" double kq_std; //!< standard deviation of average kq, used for + //!< subcritical multiplication problems +extern "C" double + ks; //!< average ks over batches, used for subcritical multiplication problems +extern "C" double ks_std; //!< standard deviation of average ks, used for + //!< subcritical multiplication problems extern "C" double k_col_abs; //!< sum over batches of k_collision * k_absorption extern "C" double k_col_tra; //!< sum over batches of k_collision * k_tracklength extern "C" double - k_abs_tra; //!< sum over batches of k_absorption * k_tracklength -extern double log_spacing; //!< lethargy spacing for energy grid searches -extern "C" int n_lost_particles; //!< cumulative number of lost particles + k_abs_tra; //!< sum over batches of k_absorption * k_tracklength +extern "C" double + kq_col_abs; //!< sum over batches of kq_collision * kq_absorption +extern "C" double + kq_col_tra; //!< sum over batches of kq_collision * kq_tracklength +extern "C" double kq_abs_tra; +extern double log_spacing; //!< lethargy spacing for energy grid searches +extern "C" int n_lost_particles; //!< cumulative number of lost particles extern "C" bool need_depletion_rx; //!< need to calculate depletion rx? extern "C" int restart_batch; //!< batch at which a restart job resumed extern "C" bool satisfy_triggers; //!< have tally triggers been satisfied? @@ -46,9 +65,25 @@ extern int64_t work_per_rank; //!< number of particles per MPI rank extern const RegularMesh* entropy_mesh; extern const RegularMesh* ufs_mesh; -extern vector k_generation; +extern vector> k_generation; +extern vector> kq_generation; +extern vector> ks_generation; extern vector work_index; +enum class KEstimator : int { + K_COLLISION = 0, + K_ABSORPTION = 1, + K_TRACKLENGTH = 2, + N_ESTIMATORS +}; + +constexpr int N_K_EST = static_cast(KEstimator::N_ESTIMATORS); +extern std::array, N_K_EST> k_kq_products; +extern std::array, N_K_EST> k_kq_product; + +extern std::array k_combined_weights; +extern std::array kq_combined_weights; + } // namespace simulation //============================================================================== diff --git a/include/openmc/state_point.h b/include/openmc/state_point.h index fb1aaf7b985..00df7d4419e 100644 --- a/include/openmc/state_point.h +++ b/include/openmc/state_point.h @@ -49,6 +49,8 @@ void write_source_bank(hid_t group_id, span source_bank, void read_source_bank( hid_t group_id, vector& sites, bool distribute); +void write_global_tallies(hid_t file_id); +void read_global_tallies(hid_t file_id); void write_tally_results_nr(hid_t file_id); void restart_set_keff(); void write_unstructured_mesh_results(); diff --git a/include/openmc/tallies/tally.h b/include/openmc/tallies/tally.h index 374daff92a0..3437cf8356b 100644 --- a/include/openmc/tallies/tally.h +++ b/include/openmc/tallies/tally.h @@ -222,6 +222,9 @@ namespace simulation { //! Global tallies (such as k-effective estimators) extern xt::xtensor_fixed> global_tallies; +//! Global tallies for first generation (for subcritical multiplication) +extern xt::xtensor_fixed> + global_tallies_first_gen; //! Number of realizations for global tallies extern "C" int32_t n_realizations; @@ -230,8 +233,14 @@ extern "C" int32_t n_realizations; extern double global_tally_absorption; extern double global_tally_collision; extern double global_tally_tracklength; +extern double global_tally_tracklength_sq; extern double global_tally_leakage; +extern double global_tally_absorption_first_gen; +extern double global_tally_collision_first_gen; +extern double global_tally_tracklength_first_gen; +extern double global_tally_tracklength_sq_first_gen; + //============================================================================== // Non-member functions //============================================================================== diff --git a/openmc/lib/tally.py b/openmc/lib/tally.py index c17b16597f9..fe8e83e3e2e 100644 --- a/openmc/lib/tally.py +++ b/openmc/lib/tally.py @@ -127,11 +127,11 @@ def global_tallies(): """ ptr = POINTER(c_double)() _dll.openmc_global_tallies(ptr) - array = as_array(ptr, (4, 3)) + array = as_array(ptr, (5, 3)) # Get sum, sum-of-squares, and number of realizations - sum_ = array[:, 1] - sum_sq = array[:, 2] + sum_ = array[:4, 1] + sum_sq = array[:4, 2] n = num_realizations() # Determine mean diff --git a/openmc/settings.py b/openmc/settings.py index 3cf662e311e..82c7036f5c9 100644 --- a/openmc/settings.py +++ b/openmc/settings.py @@ -366,6 +366,13 @@ class Settings: .. versionadded::0.14.0 write_initial_source : bool Indicate whether to write the initial source distribution to file + + .. versionadded::0.15.4 + calculate_subcritical_k : bool + Indicate whether to calculate and report the subcritical multiplication + factor k in fixed source simulations + print_all_k_factors : bool + In subcritical multipliction factor calculations, indicate whether to print all k factors (k, ks, kq) during the transport calculation """ def __init__(self, **kwargs): @@ -378,6 +385,8 @@ def __init__(self, **kwargs): self._max_write_lost_particles = None self._particles = None self._keff_trigger = None + self._calculate_subcritical_k = False + self._print_all_k_factors = False # Energy mode subelement self._energy_mode = None @@ -1396,6 +1405,30 @@ def free_gas_threshold(self, free_gas_threshold: float | None): cv.check_greater_than('free gas threshold', free_gas_threshold, 0.0) self._free_gas_threshold = free_gas_threshold + @property + def calculate_subcritical_k(self) -> bool: + return self._calculate_subcritical_k + + @calculate_subcritical_k.setter + def calculate_subcritical_k(self, calculate_subcritical_k: bool): + cv.check_type('calculate subcritical k', calculate_subcritical_k, bool) + if not self._run_mode == RunMode.FIXED_SOURCE: + raise ValueError("calculate_subcritical_k can only be set when " + "run_mode is 'fixed source'") + self._calculate_subcritical_k = calculate_subcritical_k + + @property + def print_all_k_factors(self) -> bool: + return self._print_all_k_factors + + @print_all_k_factors.setter + def print_all_k_factors(self, print_all_k_factors: bool): + cv.check_type('print all k factors', print_all_k_factors, bool) + if not self._run_mode == RunMode.FIXED_SOURCE or not self._calculate_subcritical_k: + raise ValueError("print_all_k_factors can only be set when " + "run_mode is 'fixed source' and calculate_subcritical_k is True") + self._print_all_k_factors = print_all_k_factors + def _create_run_mode_subelement(self, root): elem = ET.SubElement(root, "run_mode") elem.text = self._run_mode.value @@ -1928,6 +1961,16 @@ def _create_free_gas_threshold_subelement(self, root): element = ET.SubElement(root, "free_gas_threshold") element.text = str(self._free_gas_threshold) + def _create_calculate_subcritical_k_subelement(self, root): + if self._calculate_subcritical_k: + elem = ET.SubElement(root, "calculate_subcritical_k") + elem.text = str(self._calculate_subcritical_k).lower() + + def _create_print_all_k_factors_subelement(self, root): + if self._print_all_k_factors: + elem = ET.SubElement(root, "print_all_k_factors") + elem.text = str(self._print_all_k_factors).lower() + def _eigenvalue_from_xml_element(self, root): elem = root.find('eigenvalue') if elem is not None: @@ -2394,6 +2437,16 @@ def _free_gas_threshold_from_xml_element(self, root): if text is not None: self.free_gas_threshold = float(text) + def _calculate_subcritical_k_from_xml_element(self, root): + text = get_text(root, 'calculate_subcritical_k') + if text is not None: + self.calculate_subcritical_k = text in ('true', '1') + + def _print_all_k_factors_from_xml_element(self, root): + text = get_text(root, 'print_all_k_factors') + if text is not None: + self.print_all_k_factors = text in ('true', '1') + def to_xml_element(self, mesh_memo=None): """Create a 'settings' element to be written to an XML file. @@ -2466,6 +2519,8 @@ def to_xml_element(self, mesh_memo=None): self._create_use_decay_photons_subelement(element) self._create_source_rejection_fraction_subelement(element) self._create_free_gas_threshold_subelement(element) + self._create_calculate_subcritical_k_subelement(element) + self._create_print_all_k_factors_subelement(element) # Clean the indentation in the file to be user-readable clean_indentation(element) diff --git a/openmc/statepoint.py b/openmc/statepoint.py index a10ec3a839d..5bd21339d13 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -36,6 +36,8 @@ class StatePoint: Attributes ---------- + calculate_subcritical_k : bool + Whether subcritical multiplication factors were calculated in fixed source mode cmfd_on : bool Indicate whether CMFD is active cmfd_balance : numpy.ndarray @@ -79,12 +81,21 @@ class StatePoint: Cross-product of absorption and tracklength estimates of k-effective k_generation : numpy.ndarray Estimate of k-effective for each batch/generation + kq_generation : numpy.ndarray + Estimate of kq for each batch/generation + ks_generation : numpy.ndarray + Estimate of ks for each batch/generation keff : uncertainties.UFloat Combined estimator for k-effective - .. versionadded:: 0.13.1 + kq : uncertainties.UFloat + Combined estimator for kq + ks : uncertainaties.UFloat + Combined estimator for ks meshes : dict Dictionary whose keys are mesh IDs and whose values are MeshBase objects + multiplication : uncertainties.UFloat + Estimator for multiplication of initial source neutrons n_batches : int Number of batches n_inactive : int @@ -111,6 +122,8 @@ class StatePoint: 'E', 'wgt', 'delayed_group', 'surf_id', and 'particle', corresponding to the position, direction, energy, weight, delayed group, surface ID and particle type of the source site, respectively. + source_effectiveness : uncertainties.UFloat + Source effectiveness factor. source_present : bool Indicate whether source sites are present sparse : bool @@ -169,6 +182,15 @@ def __enter__(self): def __exit__(self, *exc): self.close() + @property + def calculate_subcritical_k(self): + """Indicate whether subcritical multiplication was calculated.""" + # If we are in fixed source mode and k_generation exists, + # then calculate_subcritical_k was True during the simulation. + if self.run_mode == 'fixed source' and 'k_generation' in self._f: + return True + return False + @property def cmfd_on(self): return self._f.attrs['cmfd_on'] > 0 @@ -247,8 +269,8 @@ def global_tallies(self): ('mean', 'f8'), ('std_dev', 'f8')]) gt['name'] = ['k-collision', 'k-absorption', 'k-tracklength', 'leakage'] - gt['sum'] = data[:,1] - gt['sum_sq'] = data[:,2] + gt['sum'] = data[:,0] + gt['sum_sq'] = data[:,1] # Calculate mean and sample standard deviation of mean n = self.n_realizations @@ -268,17 +290,63 @@ def k_cmfd(self): @property def k_generation(self): - if self.run_mode == 'eigenvalue': - return self._f['k_generation'][()] + if 'k_generation' in self._f: + arr = self._f['k_generation'][()] + if arr.ndim==1: + return arr + elif arr.ndim==2: + return uarray(arr[:,0],arr[:,1]) + else: + raise ValueError(f'k_generation shape ({arr.shape}) must be either 1d or 2d') + else: + return None + + @property + def kq_generation(self): + if 'kq_generation' in self._f: + arr = self._f['kq_generation'][()] + if arr.ndim==1: + return arr + elif arr.ndim==2: + return uarray(arr[:,0],arr[:,1]) + else: + raise ValueError(f'kq_generation shape ({arr.shape}) must be either 1d or 2d') else: return None + @property + def ks_generation(self): + if 'ks_generation' in self._f: + arr = self._f['ks_generation'][()] + if arr.ndim==1: + return arr + elif arr.ndim==2: + return uarray(arr[:,0],arr[:,1]) + else: + raise ValueError(f'ks_generation shape ({arr.shape}) must be either 1d or 2d') + 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 + + @property + def kq(self): + if 'kq_combined' in self._f: + return ufloat(*self._f['kq_combined'][()]) + else: + return None + + @property + def ks(self): + if 'ks_combined' in self._f: + return ufloat(*self._f['ks_combined'][()]) + else: + return None @property def k_combined(self): @@ -322,6 +390,19 @@ def meshes(self): self._meshes_read = True return self._meshes + + @property + def multiplication(self): + if self.run_mode == 'eigenvalue': + k_gen = self.k_generation + k_inactive = k_gen[:self.n_inactive] + cumprod = np.cumprod(k_inactive) + cumprod[-1] *= 1/(1-self.keff) + return 1.0+np.sum(cumprod) + elif self.calculate_subcritical_k: + return 1/(1 - self.keff) + else: + return None @property def n_batches(self): @@ -370,6 +451,15 @@ def stride(self): @property def source(self): return self._f['source_bank'][()] if self.source_present else None + + @property + def source_effectiveness(self): + if self.run_mode == 'eigenvalue': + k_gen = self.k_generation + k_inactive = k_gen[:self.n_inactive] + return np.prod(k_inactive/self.keff) + else: + return None @property def source_present(self): diff --git a/src/eigenvalue.cpp b/src/eigenvalue.cpp index dc4fbb420b8..593e5046433 100644 --- a/src/eigenvalue.cpp +++ b/src/eigenvalue.cpp @@ -24,8 +24,9 @@ #include // for min #include // for sqrt, abs, pow -#include // for back_inserter -#include //for infinity +#include +#include // for back_inserter +#include //for infinity #include namespace openmc { @@ -36,8 +37,12 @@ namespace openmc { namespace simulation { -double keff_generation; +array keff_generation; +array kq_generation_val; +array ks_generation_val; array k_sum; +array kq_sum; +array ks_sum; vector entropy; xt::xtensor source_frac; @@ -47,20 +52,76 @@ xt::xtensor source_frac; // Non-member functions //============================================================================== +double compute_cov_M_kq(int i, int j) +{ + // Number of active generations + int N = settings::n_particles; + int idx = static_cast(GlobalTally::K_TRACKLENGTH); + + // Note before conversion, k_generation is M + auto [mean_m, m_std] = simulation::k_generation.back(); + double mean_kq = simulation::kq_generation.back()[0]; + + double sum_k_kq = + simulation::k_kq_products[idx][idx] - simulation::k_kq_product[idx][idx]; + return (sum_k_kq - mean_m * mean_kq) / (N - 1); +} + +void calculate_generation_ks() +{ + auto [m, m_std] = simulation::k_generation.back(); + auto [kq, kq_std] = simulation::kq_generation.back(); + double ks_mean = 1 - kq / (m - 1); + double cov_M_kq = compute_cov_M_kq(0, 0); + double rho = cov_M_kq / (m_std * kq_std); + double ks_std = + std::sqrt(std::pow(kq_std, 2) / std::pow(m - 1, 2) + + std::pow(kq, 2) * std::pow(m_std, 2) / std::pow(m - 1, 4) - + 2 * kq / std::pow(m - 1, 3) * rho * m_std * kq_std); + simulation::ks_generation.push_back({ks_mean, ks_std}); +} + void calculate_generation_keff() { - const auto& gt = simulation::global_tallies; + calculate_generation_keff(KeffType::k); +} + +void calculate_generation_keff(KeffType type) +{ + // Initialize variables + xt::xtensor_fixed> gt; + array* keff_generation_ptr; + vector>* k_generation_ptr; + switch (type) { + case KeffType::k: + gt = simulation::global_tallies; + keff_generation_ptr = &simulation::keff_generation; + k_generation_ptr = &simulation::k_generation; + break; + case KeffType::kq: + gt = simulation::global_tallies_first_gen; + keff_generation_ptr = &simulation::kq_generation_val; + k_generation_ptr = &simulation::kq_generation; + break; + case KeffType::ks: + calculate_generation_ks(); + return; + } - // Get keff for this generation by subtracting off the starting value - simulation::keff_generation = + (*keff_generation_ptr)[0] = gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) - - simulation::keff_generation; + (*keff_generation_ptr)[0]; + (*keff_generation_ptr)[1] = + gt(GlobalTally::K_TRACKLENGTH_SQ, TallyResult::VALUE) - + (*keff_generation_ptr)[1]; - double keff_reduced; + array keff_reduced; #ifdef OPENMC_MPI if (settings::solver_type != SolverType::RANDOM_RAY) { // Combine values across all processors - MPI_Allreduce(&simulation::keff_generation, &keff_reduced, 1, MPI_DOUBLE, + MPI_Allreduce(&(*keff_generation_ptr)[0], &keff_reduced[0], 1, MPI_DOUBLE, + MPI_SUM, mpi::intracomm); + MPI_Allreduce(&(*keff_generation_ptr)[1], &keff_reduced[1], 1, MPI_DOUBLE, MPI_SUM, mpi::intracomm); } else { // If using random ray, MPI parallelism is provided by domain replication. @@ -68,19 +129,36 @@ void calculate_generation_keff() // such that all ranks have identical scalar flux vectors, and will all // independently compute the same value of k. Thus, there is no need to // perform any additional MPI reduction here. - keff_reduced = simulation::keff_generation; + keff_reduced = *keff_generation_ptr; } #else - keff_reduced = simulation::keff_generation; + keff_reduced = *keff_generation_ptr; #endif // Normalize single batch estimate of k // TODO: This should be normalized by total_weight, not by n_particles if (settings::solver_type != SolverType::RANDOM_RAY) { - keff_reduced /= settings::n_particles; + keff_reduced[0] /= settings::n_particles; + keff_reduced[1] /= settings::n_particles; } + double k_mean = keff_reduced[0]; + double k_std = std::sqrt( + (keff_reduced[1] - std::pow(k_mean, 2)) / (settings::n_particles - 1)); + k_generation_ptr->push_back({k_mean, k_std}); +} + +std::pair convert_m_to_k(double m, double m_std) +{ + double k = 1.0 - 1.0 / m; + double k_std = m_std / (m * m); + return {k, k_std}; +} - simulation::k_generation.push_back(keff_reduced); +std::pair convert_k_to_m(double k, double k_std) +{ + double m = 1.0 / (1.0 - k); + double m_std = k_std / std::pow(1.0 - k, 2); + return {m, m_std}; } void synchronize_bank() @@ -376,6 +454,11 @@ void synchronize_bank() } void calculate_average_keff() +{ + calculate_average_keff(KeffType::k); +} + +void calculate_average_keff(KeffType type) { // Determine overall generation and number of active generations int i = overall_generation() - 1; @@ -386,19 +469,50 @@ void calculate_average_keff() } else { n = 0; } + // Initialize variables + double keff; + double keff_std; + + array* keff_generation_ptr; + vector>* k_generation_ptr; + array* k_sum_ptr; + double* k_ptr; + double* k_std_ptr; + switch (type) { + case KeffType::k: + keff_generation_ptr = &simulation::keff_generation; + k_generation_ptr = &simulation::k_generation; + k_sum_ptr = &simulation::k_sum; + k_ptr = &simulation::k; + k_std_ptr = &simulation::k_std; + break; + case KeffType::kq: + keff_generation_ptr = &simulation::kq_generation_val; + k_generation_ptr = &simulation::kq_generation; + k_sum_ptr = &simulation::kq_sum; + k_ptr = &simulation::kq; + k_std_ptr = &simulation::kq_std; + break; + case KeffType::ks: + keff_generation_ptr = &simulation::ks_generation_val; + k_generation_ptr = &simulation::ks_generation; + k_sum_ptr = &simulation::ks_sum; + k_ptr = &simulation::ks; + k_std_ptr = &simulation::ks_std; + break; + } if (n <= 0) { // For inactive generations, use current generation k as estimate for next // generation - simulation::keff = simulation::k_generation[i]; + keff = (*k_generation_ptr)[i][0]; } else { // Sample mean of keff - simulation::k_sum[0] += simulation::k_generation[i]; - simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2); + (*k_sum_ptr)[0] += (*k_generation_ptr)[i][0]; + (*k_sum_ptr)[1] += std::pow((*k_generation_ptr)[i][0], 2); // Determine mean - simulation::keff = simulation::k_sum[0] / n; - + keff = (*k_sum_ptr)[0] / n; if (n > 1) { double t_value; if (settings::confidence_intervals) { @@ -410,10 +524,8 @@ void calculate_average_keff() } // Standard deviation of the sample mean of k - simulation::keff_std = - t_value * - std::sqrt( - (simulation::k_sum[1] / n - std::pow(simulation::keff, 2)) / (n - 1)); + keff_std = t_value * + std::sqrt(((*k_sum_ptr)[1] / n - std::pow(keff, 2)) / (n - 1)); // In some cases (such as an infinite medium problem), random ray // may estimate k exactly and in an unvarying manner between iterations. @@ -421,14 +533,25 @@ void calculate_average_keff() // power operations may cause an extremely small negative value to occur // inside the sqrt operation, leading to NaN. If this occurs, we check for // it and set the std dev to zero. - if (!std::isfinite(simulation::keff_std)) { - simulation::keff_std = 0.0; + if (!std::isfinite(keff_std)) { + keff_std = 0.0; } } } + (*k_ptr) = keff; + (*k_std_ptr) = keff_std; + if (settings::run_mode == RunMode::EIGENVALUE and type == KeffType::k) { + // Only set simulation::keff for eigenvalue mode, since it's used to bias + // physics + simulation::keff = keff; + simulation::keff_std = keff_std; + } } -int openmc_get_keff(double* k_combined) +int get_combined_k_from_tallies(double* k_combined, + std::array& k_combined_weights, double keff, double keff_std, + xt::xtensor_fixed> gt, + double k_col_abs, double k_col_tra, double k_abs_tra) { k_combined[0] = 0.0; k_combined[1] = 0.0; @@ -437,8 +560,8 @@ int openmc_get_keff(double* k_combined) // there is a N-3 term in a denominator. if (simulation::n_realizations <= 3 || settings::solver_type == SolverType::RANDOM_RAY) { - k_combined[0] = simulation::keff; - k_combined[1] = simulation::keff_std; + k_combined[0] = keff; + k_combined[1] = keff_std; if (simulation::n_realizations <= 1) { k_combined[1] = std::numeric_limits::infinity(); } @@ -449,8 +572,6 @@ int openmc_get_keff(double* k_combined) int64_t n = simulation::n_realizations; // Copy estimates of k-effective and its variance (not variance of the mean) - const auto& gt = simulation::global_tallies; - array kv {}; xt::xtensor cov = xt::zeros({3, 3}); kv[0] = gt(GlobalTally::K_COLLISION, TallyResult::SUM) / n; @@ -467,9 +588,9 @@ int openmc_get_keff(double* k_combined) (n - 1); // Calculate covariances based on sums with Bessel's correction - cov(0, 1) = (simulation::k_col_abs - n * kv[0] * kv[1]) / (n - 1); - cov(0, 2) = (simulation::k_col_tra - n * kv[0] * kv[2]) / (n - 1); - cov(1, 2) = (simulation::k_abs_tra - n * kv[1] * kv[2]) / (n - 1); + cov(0, 1) = (k_col_abs - n * kv[0] * kv[1]) / (n - 1); + cov(0, 2) = (k_col_tra - n * kv[0] * kv[2]) / (n - 1); + cov(1, 2) = (k_abs_tra - n * kv[1] * kv[2]) / (n - 1); cov(1, 0) = cov(0, 1); cov(2, 0) = cov(0, 2); cov(2, 1) = cov(1, 2); @@ -548,9 +669,13 @@ int openmc_get_keff(double* k_combined) S[2] += (cov(k, k) + cov(i, j) - cov(j, k) - cov(i, k)) * kv[l] * kv[j]; // Add to sum for combined k-effective + k_combined_weights[l] += f; k_combined[0] += f * kv[l]; g += f; } + for (auto& w : k_combined_weights) { + w /= g; + } // Complete calculations of S sums for (auto& S_i : S) { @@ -587,6 +712,48 @@ int openmc_get_keff(double* k_combined) return 0; } +int openmc_get_keff(double* k_combined) +{ + return get_combined_k_from_tallies(k_combined, simulation::k_combined_weights, + simulation::k, simulation::k_std, simulation::global_tallies, + simulation::k_col_abs, simulation::k_col_tra, simulation::k_abs_tra); +} + +int openmc_get_kq(double* kq_combined) +{ + return get_combined_k_from_tallies(kq_combined, + simulation::kq_combined_weights, simulation::kq, simulation::kq_std, + simulation::global_tallies_first_gen, simulation::kq_col_abs, + simulation::kq_col_tra, simulation::kq_abs_tra); +} + +int openmc_get_ks(double* ks_combined, double* k_combined, double* kq_combined) +{ + ks_combined[0] = 1.0 - kq_combined[0] / (k_combined[0] - 1.0); + double total_cov {0.0}; + array M; + array kq; + int n = simulation::n_realizations; + for (int i = 0; i < simulation::N_K_EST; ++i) { + for (int j = 0; j < simulation::N_K_EST; ++j) { + double cov = + (simulation::k_kq_products[i][j] - + n * simulation::global_tallies(i, TallyResult::SUM) * + simulation::global_tallies_first_gen(j, TallyResult::SUM) / + std::pow(n, 2)) / + (n * (n - 1)); // Note extra division by n to get standard error + total_cov += simulation::k_combined_weights[i] * + simulation::kq_combined_weights[j] * cov; + } + } + ks_combined[1] = + std::sqrt(std::pow(kq_combined[1], 2) / (std::pow(k_combined[0] - 1, 2)) + + std::pow(kq_combined[0] * k_combined[1], 2) / + std::pow(k_combined[0] - 1, 4) - + 2 * kq_combined[0] / std::pow(k_combined[0] - 1, 3) * total_cov); + return 0; +} + void shannon_entropy() { // Get source weight in each mesh bin @@ -679,26 +846,83 @@ double ufs_get_weight(const Particle& p) void write_eigenvalue_hdf5(hid_t group) { + auto n = simulation::k_generation.size(); + xt::xtensor k_generation({n, 2}); + xt::xtensor kq_generation({n, 2}); + xt::xtensor ks_generation({n, 2}); + for (int i = 0; i < n; ++i) { + double k, k_std; + if (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k) { + // Temporary fix until formula for combined estimator can be implemented + // correctly + auto [k0, k1] = convert_m_to_k( + simulation::k_generation[i][0], simulation::k_generation[i][1]); + k = k0; + k_std = k1; + + kq_generation(i, 0) = simulation::kq_generation[i][0]; + kq_generation(i, 1) = simulation::kq_generation[i][1]; + + ks_generation(i, 0) = simulation::ks_generation[i][0]; + ks_generation(i, 1) = simulation::ks_generation[i][1]; + } else { + k = simulation::k_generation[i][0]; + k_std = simulation::k_generation[i][1]; + } + k_generation(i, 0) = k; + k_generation(i, 1) = k_std; + } write_dataset(group, "n_inactive", settings::n_inactive); write_dataset(group, "generations_per_batch", settings::gen_per_batch); - write_dataset(group, "k_generation", simulation::k_generation); + write_dataset(group, "k_generation", k_generation); if (settings::entropy_on) { write_dataset(group, "entropy", simulation::entropy); } - write_dataset(group, "k_col_abs", simulation::k_col_abs); - write_dataset(group, "k_col_tra", simulation::k_col_tra); - write_dataset(group, "k_abs_tra", simulation::k_abs_tra); + if (settings::run_mode == RunMode::EIGENVALUE) { + write_dataset(group, "k_col_abs", simulation::k_col_abs); + write_dataset(group, "k_col_tra", simulation::k_col_tra); + write_dataset(group, "k_abs_tra", simulation::k_abs_tra); + } + array k_combined; openmc_get_keff(k_combined.data()); + if (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k) { + auto [k0, k1] = convert_m_to_k(k_combined[0], k_combined[1]); + k_combined[0] = k0; + k_combined[1] = k1; + } write_dataset(group, "k_combined", k_combined); + if (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k) { + write_dataset(group, "kq_generation", kq_generation); + array kq_combined; + openmc_get_kq(kq_combined.data()); + fmt::print("kq_combined: {} +/- {}\n", kq_combined[0], kq_combined[1]); + write_dataset(group, "kq_combined", kq_combined); + + // Convert back to m for calculation of ks + array ks_combined; + std::tie(k_combined[0], k_combined[1]) = + convert_k_to_m(k_combined[0], k_combined[1]); + openmc_get_ks(ks_combined.data(), k_combined.data(), kq_combined.data()); + fmt::print("ks_combined: {} +/- {}\n", ks_combined[0], ks_combined[1]); + write_dataset(group, "ks_generation", ks_generation); + write_dataset(group, "ks_combined", ks_combined); + } } void read_eigenvalue_hdf5(hid_t group) { read_dataset(group, "generations_per_batch", settings::gen_per_batch); - int n = simulation::restart_batch * settings::gen_per_batch; + size_t n = simulation::restart_batch * settings::gen_per_batch; + xt::xtensor k_generation({n, 2}); + read_dataset(group, "k_generation", k_generation); simulation::k_generation.resize(n); - read_dataset(group, "k_generation", simulation::k_generation); + for (int i = 0; i < n; ++i) { + simulation::k_generation[i] = {k_generation(i, 0), k_generation(i, 1)}; + } if (settings::entropy_on) { read_dataset(group, "entropy", simulation::entropy); } diff --git a/src/finalize.cpp b/src/finalize.cpp index 5ca9d574627..824656241e0 100644 --- a/src/finalize.cpp +++ b/src/finalize.cpp @@ -89,6 +89,7 @@ int openmc_finalize() settings::entropy_on = false; settings::event_based = false; settings::free_gas_threshold = 400.0; + settings::calculate_subcritical_k = false; settings::gen_per_batch = 1; settings::legendre_to_tabular = true; settings::legendre_to_tabular_points = -1; @@ -208,6 +209,10 @@ int openmc_reset() simulation::k_col_abs = 0.0; simulation::k_col_tra = 0.0; simulation::k_abs_tra = 0.0; + simulation::kq_col_abs = 0.0; + simulation::kq_col_tra = 0.0; + simulation::kq_abs_tra = 0.0; + simulation::k_kq_products.fill({}); simulation::k_sum = {0.0, 0.0}; simulation::satisfy_triggers = false; diff --git a/src/output.cpp b/src/output.cpp index ae2daaffc14..8471417d9b1 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -19,6 +19,7 @@ #endif #include "xtensor/xview.hpp" +#include "openmc/array.h" #include "openmc/capi.h" #include "openmc/cell.h" #include "openmc/constants.h" @@ -371,11 +372,12 @@ void print_build_info() void print_columns() { if (settings::entropy_on) { - fmt::print(" Bat./Gen. k Entropy Average k \n" - " ========= ======== ======== ====================\n"); + fmt::print( + " Bat./Gen. k Entropy Average k \n" + " ========= ==================== ======== ====================\n"); } else { - fmt::print(" Bat./Gen. k Average k\n" - " ========= ======== ====================\n"); + fmt::print(" Bat./Gen. k Average k\n" + " ========= ==================== ====================\n"); } } @@ -393,20 +395,60 @@ void print_generation() // write out batch/generation and generation k-effective auto batch_and_gen = std::to_string(simulation::current_batch) + "/" + std::to_string(simulation::current_gen); - fmt::print(" {:>9} {:8.5f}", batch_and_gen, simulation::k_generation[idx]); + array k_gen; + k_gen[0] = simulation::k_generation[idx][0]; + k_gen[1] = simulation::k_generation[idx][1]; + if (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k) { + // convert from multiplication factor to subcritical k + auto [k0, k1] = convert_m_to_k(k_gen[0], k_gen[1]); + k_gen[0] = k0; + k_gen[1] = k1; + } + + double k, k_std; + if (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k) { + // convert from multiplication factor to subcritical k + auto [k0, k1] = convert_m_to_k(simulation::k, simulation::k_std); + k = k0; + k_std = k1; + } else { + k = simulation::k; + k_std = simulation::k_std; + } + + if (settings::print_all_k_factors) { + print_generation_values(idx, n, batch_and_gen + " k:", k_gen, k, k_std); + array kq_gen; + array ks_gen; + kq_gen = simulation::kq_generation[idx]; + ks_gen = simulation::ks_generation[idx]; + std::string spaces(batch_and_gen.size(), ' '); + print_generation_values( + idx, n, spaces + "kq:", kq_gen, simulation::kq, simulation::kq_std); + print_generation_values( + idx, n, spaces + "ks:", ks_gen, simulation::ks, simulation::ks_std); + } else { + print_generation_values(idx, n, batch_and_gen, k_gen, k, k_std); + } +} + +void print_generation_values(int idx, int n, std::string batch_and_gen, + array k_gen, double k, double k_std) +{ + fmt::print(" {:>9} {:8.5f} +/-{:8.5f}", batch_and_gen, k_gen[0], k_gen[1]); // write out entropy info if (settings::entropy_on) { fmt::print(" {:8.5f}", simulation::entropy[idx]); } - if (n > 1) { - fmt::print(" {:8.5f} +/-{:8.5f}", simulation::keff, simulation::keff_std); + fmt::print(" {:8.5f} +/-{:8.5f}", k, k_std); } fmt::print("\n"); std::fflush(stdout); } - //============================================================================== void show_time(const char* label, double secs, int indent_level) diff --git a/src/particle.cpp b/src/particle.cpp index a035e76b917..a83d364e877 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -279,7 +279,10 @@ void Particle::event_advance() } // Score track-length estimate of k-eff - if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) { + if ((settings::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) && + type().is_neutron()) { keff_tally_tracklength() += wgt() * distance * macro_xs().nu_fission; } @@ -341,7 +344,10 @@ void Particle::event_collide() { // Score collision estimate of keff - if (settings::run_mode == RunMode::EIGENVALUE && type().is_neutron()) { + if ((settings::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) && + type().is_neutron()) { keff_tally_collision() += wgt() * macro_xs().nu_fission / macro_xs().total; } @@ -501,8 +507,11 @@ void Particle::event_death() global_tally_absorption += keff_tally_absorption(); #pragma omp atomic global_tally_collision += keff_tally_collision(); + double val = keff_tally_tracklength(); #pragma omp atomic - global_tally_tracklength += keff_tally_tracklength(); + global_tally_tracklength += val; +#pragma omp atomic + global_tally_tracklength_sq += val * val; #pragma omp atomic global_tally_leakage += keff_tally_leakage(); diff --git a/src/physics.cpp b/src/physics.cpp index 1f72e4fac07..3c47f303e72 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -177,9 +177,14 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) double weight = settings::ufs_on ? ufs_get_weight(p) : 1.0; // Determine the expected number of neutrons produced - double nu_t = p.wgt() / simulation::keff * weight * - p.neutron_xs(i_nuclide).nu_fission / - p.neutron_xs(i_nuclide).total; + double nu_t {0}; + if (settings::run_mode == RunMode::EIGENVALUE) { + nu_t = p.wgt() / simulation::keff * weight * + p.neutron_xs(i_nuclide).nu_fission / p.neutron_xs(i_nuclide).total; + } else { + nu_t = p.wgt() * weight * p.neutron_xs(i_nuclide).nu_fission / + p.neutron_xs(i_nuclide).total; + } // Sample the number of neutrons produced int nu = static_cast(nu_t); @@ -649,7 +654,9 @@ 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::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) { p.keff_tally_absorption() += wgt_absorb * p.neutron_xs(i_nuclide).nu_fission / p.neutron_xs(i_nuclide).absorption; @@ -659,7 +666,9 @@ 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::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) { p.keff_tally_absorption() += p.wgt() * p.neutron_xs(i_nuclide).nu_fission / p.neutron_xs(i_nuclide).absorption; diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 5970dc65a2b..ff803f64343 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -404,6 +404,7 @@ void RandomRaySimulation::simulate() // Store random ray k-eff into OpenMC's native k-eff variable global_tally_tracklength = domain_->k_eff_; + global_tally_tracklength_sq = std::pow(domain_->k_eff_, 2); } // Execute all tallying tasks, if this is an active batch diff --git a/src/settings.cpp b/src/settings.cpp index 5b472468fc4..d7ff0cf85bc 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -130,6 +130,8 @@ std::unordered_set sourcepoint_batch; std::unordered_set statepoint_batch; double source_rejection_fraction {0.05}; double free_gas_threshold {400.0}; +bool calculate_subcritical_k {false}; +bool print_all_k_factors {false}; std::unordered_set source_write_surf_id; CollisionTrackConfig collision_track_config {}; int64_t ssw_max_particles; @@ -676,6 +678,35 @@ void read_settings_xml(pugi::xml_node root) free_gas_threshold = std::stod(get_node_value(root, "free_gas_threshold")); } + if (check_for_node(root, "calculate_subcritical_k")) { + if (run_mode == RunMode::FIXED_SOURCE) { + if (solver_type != SolverType::MONTE_CARLO) { + fatal_error("The 'calculate_subcritical_k' setting is only valid in " + "fixed source mode with the Monte Carlo solver."); + } + calculate_subcritical_k = + get_node_value_bool(root, "calculate_subcritical_k"); + } else { + fatal_error("The 'calculate_subcritical_k' setting is only valid in " + "fixed source mode."); + } + } + + if (check_for_node(root, "print_all_k_factors")) { + if (run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k) { + if (solver_type != SolverType::MONTE_CARLO) { + fatal_error("The 'print_all_k_factors' setting is only valid in " + "fixed source mode with the Monte Carlo solver."); + } + print_all_k_factors = get_node_value_bool(root, "print_all_k_factors"); + } else { + fatal_error( + "The 'print_all_k_factors' setting is only valid in " + "fixed source mode with calculate_subcritical_k set to true."); + } + } + // Survival biasing if (check_for_node(root, "survival_biasing")) { survival_biasing = get_node_value_bool(root, "survival_biasing"); diff --git a/src/simulation.cpp b/src/simulation.cpp index 18e40a8bc73..7fd52c3e22b 100644 --- a/src/simulation.cpp +++ b/src/simulation.cpp @@ -1,8 +1,10 @@ #include "openmc/simulation.h" +#include "openmc/array.h" #include "openmc/bank.h" #include "openmc/capi.h" #include "openmc/collision_track.h" +#include "openmc/constants.h" #include "openmc/container_util.h" #include "openmc/eigenvalue.h" #include "openmc/error.h" @@ -121,6 +123,7 @@ int openmc_simulation_init() simulation::ct_current_file = 1; simulation::ssw_current_file = 1; simulation::k_generation.clear(); + simulation::kq_generation.clear(); simulation::entropy.clear(); openmc_reset(); @@ -142,6 +145,9 @@ int openmc_simulation_init() if (settings::run_mode == RunMode::FIXED_SOURCE) { if (settings::solver_type == SolverType::MONTE_CARLO) { header("FIXED SOURCE TRANSPORT SIMULATION", 3); + if (settings::verbosity >= 7 && settings::calculate_subcritical_k) { + print_columns(); + } } else if (settings::solver_type == SolverType::RANDOM_RAY) { header("FIXED SOURCE TRANSPORT SIMULATION (RANDOM RAY SOLVER)", 3); } @@ -304,9 +310,18 @@ int current_gen; bool initialized {false}; double keff {1.0}; double keff_std; +double k; +double k_std; +double kq; +double kq_std; +double ks; +double ks_std; double k_col_abs {0.0}; double k_col_tra {0.0}; double k_abs_tra {0.0}; +double kq_col_abs {0.0}; +double kq_col_tra {0.0}; +double kq_abs_tra {0.0}; double log_spacing; int n_lost_particles {0}; bool need_depletion_rx {false}; @@ -320,9 +335,19 @@ int64_t work_per_rank; const RegularMesh* entropy_mesh {nullptr}; const RegularMesh* ufs_mesh {nullptr}; -vector k_generation; +vector> k_generation; +vector> kq_generation; +vector> ks_generation; vector work_index; +// k estimator × kq estimator products +std::array, N_K_EST> k_kq_products; +std::array, N_K_EST> k_kq_product; + +// Combined k estimator weights +std::array k_combined_weights; +std::array kq_combined_weights; + } // namespace simulation //============================================================================== @@ -366,7 +391,9 @@ void initialize_batch() write_message( 6, "Simulating batch {:<4} (inactive)", simulation::current_batch); } else { - write_message(6, "Simulating batch {}", simulation::current_batch); + if (!settings::calculate_subcritical_k) { + write_message(6, "Simulating batch {}", simulation::current_batch); + } } } @@ -403,6 +430,7 @@ void finalize_batch() { // Reduce tallies onto master process and accumulate simulation::time_tallies.start(); + simulation::k_kq_product = simulation::k_kq_products; accumulate_tallies(); simulation::time_tallies.stop(); @@ -414,6 +442,7 @@ void finalize_batch() // Reset global tally results if (simulation::current_batch <= settings::n_inactive) { xt::view(simulation::global_tallies, xt::all()) = 0.0; + xt::view(simulation::global_tallies_first_gen, xt::all()) = 0.0; simulation::n_realizations = 0; } @@ -505,7 +534,9 @@ void finalize_batch() void initialize_generation() { - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) { // Clear out the fission bank simulation::fission_bank.resize(0); @@ -514,14 +545,25 @@ void initialize_generation() ufs_count_sites(); // Store current value of tracklength k - simulation::keff_generation = simulation::global_tallies( - GlobalTally::K_TRACKLENGTH, TallyResult::VALUE); + auto& gt = simulation::global_tallies; + simulation::keff_generation = { + gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE), + gt(GlobalTally::K_TRACKLENGTH_SQ, TallyResult::VALUE)}; + } + if (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k) { + // Store current value of tracklength kq + auto& gt_first_gen = simulation::global_tallies_first_gen; + simulation::kq_generation_val = { + gt_first_gen(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE), + gt_first_gen(GlobalTally::K_TRACKLENGTH_SQ, TallyResult::VALUE)}; } } void finalize_generation() { auto& gt = simulation::global_tallies; + auto& gt_first_gen = simulation::global_tallies_first_gen; // Update global tallies with the accumulation variables if (settings::run_mode == RunMode::EIGENVALUE) { @@ -530,30 +572,63 @@ void finalize_generation() global_tally_absorption; gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) += global_tally_tracklength; + gt(GlobalTally::K_TRACKLENGTH_SQ, TallyResult::VALUE) += + global_tally_tracklength_sq; + } else if ((settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) { + gt(GlobalTally::K_COLLISION, TallyResult::VALUE) += + settings::n_particles + global_tally_collision; + gt(GlobalTally::K_ABSORPTION, TallyResult::VALUE) += + settings::n_particles + global_tally_absorption; + gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) += + settings::n_particles + global_tally_tracklength; + gt(GlobalTally::K_TRACKLENGTH_SQ, TallyResult::VALUE) += + settings::n_particles + global_tally_tracklength_sq; + + // Update first generation tallies + gt_first_gen(GlobalTally::K_ABSORPTION, TallyResult::VALUE) += + global_tally_absorption_first_gen; + gt_first_gen(GlobalTally::K_COLLISION, TallyResult::VALUE) += + global_tally_collision_first_gen; + gt_first_gen(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) += + global_tally_tracklength_first_gen; + gt_first_gen(GlobalTally::K_TRACKLENGTH_SQ, TallyResult::VALUE) += + global_tally_tracklength_sq_first_gen; } gt(GlobalTally::LEAKAGE, TallyResult::VALUE) += global_tally_leakage; // reset tallies - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) { global_tally_collision = 0.0; global_tally_absorption = 0.0; global_tally_tracklength = 0.0; + global_tally_tracklength_sq = 0.0; + } + if (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k) { + global_tally_absorption_first_gen = 0.0; + global_tally_collision_first_gen = 0.0; + global_tally_tracklength_first_gen = 0.0; + global_tally_tracklength_sq_first_gen = 0.0; } global_tally_leakage = 0.0; - if (settings::run_mode == RunMode::EIGENVALUE && - 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 - // are run in. - sort_fission_bank(); + // For fixed source mode, we need different handling + if (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k) { - // Distribute fission bank across processors evenly + } else if (settings::run_mode == RunMode::EIGENVALUE && + settings::solver_type == SolverType::MONTE_CARLO) { + + sort_fission_bank(); synchronize_bank(); } - if (settings::run_mode == RunMode::EIGENVALUE) { - + if (settings::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) { // Calculate shannon entropy if (settings::entropy_on && settings::solver_type == SolverType::MONTE_CARLO) @@ -562,6 +637,22 @@ void finalize_generation() // Collect results and statistics calculate_generation_keff(); calculate_average_keff(); + if (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k) { + calculate_generation_keff(KeffType::kq); + calculate_average_keff(KeffType::kq); + calculate_generation_keff(KeffType::ks); + calculate_average_keff(KeffType::ks); + } + + // fmt::print("Kq {} Kq std {}\n", simulation::kq, simulation::kq_std); + // fmt::print("Kq_generation {} Kq std {}\n", + // simulation::kq_generation_val[0], + // simulation::kq_generation_val[1]); + // fmt::print("Ks {} Ks std {}\n", simulation::ks, simulation::ks_std); + // fmt::print("Ks_generation {} Ks std {}\n", + // simulation::ks_generation.back()[0], + // simulation::ks_generation.back()[1]); // Write generation output if (mpi::master && settings::verbosity >= 7) { @@ -782,15 +873,50 @@ void broadcast_results() // Also broadcast global tally results auto& gt = simulation::global_tallies; MPI_Bcast(gt.data(), gt.size(), MPI_DOUBLE, 0, mpi::intracomm); + auto& gt_first_gen = simulation::global_tallies_first_gen; + MPI_Bcast( + gt_first_gen.data(), gt_first_gen.size(), MPI_DOUBLE, 0, mpi::intracomm); // These guys are needed so that non-master processes can calculate the // combined estimate of k-effective - double temp[] { + double temp_k[] { simulation::k_col_abs, simulation::k_col_tra, simulation::k_abs_tra}; - MPI_Bcast(temp, 3, MPI_DOUBLE, 0, mpi::intracomm); - simulation::k_col_abs = temp[0]; - simulation::k_col_tra = temp[1]; - simulation::k_abs_tra = temp[2]; + MPI_Bcast(temp_k, 3, MPI_DOUBLE, 0, mpi::intracomm); + simulation::k_col_abs = temp_k[0]; + simulation::k_col_tra = temp_k[1]; + simulation::k_abs_tra = temp_k[2]; + + // These guys are needed so that non-master processes can calculate the + // combined estimate of kq + double temp_kq[] { + simulation::kq_col_abs, simulation::kq_col_tra, simulation::kq_abs_tra}; + MPI_Bcast(temp_kq, 3, MPI_DOUBLE, 0, mpi::intracomm); + simulation::kq_col_abs = temp_kq[0]; + simulation::kq_col_tra = temp_kq[1]; + simulation::kq_abs_tra = temp_kq[2]; + + double temp_products[simulation::N_K_EST * simulation::N_K_EST]; + + if (mpi::master) { + int idx = 0; + for (int i = 0; i < simulation::N_K_EST; ++i) { + for (int j = 0; j < simulation::N_K_EST; ++j) { + temp_products[idx++] = simulation::k_kq_products[i][j]; + } + } + } + + MPI_Bcast(temp_products, simulation::N_K_EST * simulation::N_K_EST, + MPI_DOUBLE, 0, mpi::intracomm); + + if (!mpi::master) { + int idx = 0; + for (int i = 0; i < simulation::N_K_EST; ++i) { + for (int j = 0; j < simulation::N_K_EST; ++j) { + simulation::k_kq_products[i][j] = temp_products[idx++]; + } + } + } } #endif @@ -798,11 +924,16 @@ void broadcast_results() void free_memory_simulation() { simulation::k_generation.clear(); + simulation::kq_generation.clear(); simulation::entropy.clear(); } void transport_history_based_single_particle(Particle& p) { + bool tally_first_generation = (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k) + ? true + : false; while (p.alive()) { p.event_calculate_xs(); if (p.alive()) { @@ -815,7 +946,25 @@ void transport_history_based_single_particle(Particle& p) p.event_collide(); } } - p.event_revive_from_secondary(); + // Check for first generation completion + if (!p.alive() && tally_first_generation) { + if (settings::calculate_subcritical_k) { + // Protect global updates with atomic to prevent data races +#pragma omp atomic + global_tally_absorption_first_gen += p.keff_tally_absorption(); +#pragma omp atomic + global_tally_collision_first_gen += p.keff_tally_collision(); +#pragma omp atomic + global_tally_tracklength_first_gen += p.keff_tally_tracklength(); +#pragma omp atomic + global_tally_tracklength_sq_first_gen += + std::pow(p.keff_tally_tracklength(), 2); + } + tally_first_generation = false; + } + if (!tally_first_generation) { + p.event_revive_from_secondary(); + } } p.event_death(); } diff --git a/src/state_point.cpp b/src/state_point.cpp index 455299529b3..ebc4b8655ab 100644 --- a/src/state_point.cpp +++ b/src/state_point.cpp @@ -118,9 +118,11 @@ 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::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) { write_eigenvalue_hdf5(file_id); - + } hid_t tallies_group = create_group(file_id, "tallies"); // Write meshes @@ -261,7 +263,7 @@ extern "C" int openmc_statepoint_write(const char* filename, bool* write_source) if (settings::reduce_tallies) { // Write global tallies - write_dataset(file_id, "global_tallies", simulation::global_tallies); + write_global_tallies(file_id); // Write tallies if (model::active_tallies.size() > 0) { @@ -363,13 +365,13 @@ void restart_set_keff() { if (simulation::restart_batch > settings::n_inactive) { for (int i = settings::n_inactive; i < simulation::restart_batch; ++i) { - simulation::k_sum[0] += simulation::k_generation[i]; - simulation::k_sum[1] += std::pow(simulation::k_generation[i], 2); + simulation::k_sum[0] += simulation::k_generation[i][0]; + simulation::k_sum[1] += std::pow(simulation::k_generation[i][0], 2); } int n = settings::gen_per_batch * simulation::n_realizations; simulation::keff = simulation::k_sum[0] / n; } else { - simulation::keff = simulation::k_generation.back(); + simulation::keff = simulation::k_generation.back()[0]; } } @@ -493,8 +495,7 @@ extern "C" int openmc_statepoint_load(const char* filename) if (mpi::master) { #endif // Read global tally data - read_dataset_lowlevel(file_id, "global_tallies", H5T_NATIVE_DOUBLE, H5S_ALL, - false, simulation::global_tallies.data()); + read_global_tallies(file_id); // Check if tally results are present bool present; @@ -870,6 +871,34 @@ void write_unstructured_mesh_results() } } } +void write_global_tallies(hid_t file_id) +{ + // Get global tallies + auto& gt = simulation::global_tallies; + auto write_view = xt::view(gt, + xt::range(static_cast(GlobalTally::K_COLLISION), + static_cast(GlobalTally::LEAKAGE) + 1), + xt::range(static_cast(TallyResult::SUM), + static_cast(TallyResult::SUM_SQ) + 1)); + auto gt_reduced = xt::empty_like(write_view); + gt_reduced = write_view; + write_dataset(file_id, "global_tallies", gt_reduced); +} + +void read_global_tallies(hid_t file_id) +{ + // Get global tallies + auto& gt = simulation::global_tallies; + auto gt_view = xt::view(gt, + xt::range(static_cast(GlobalTally::K_COLLISION), + static_cast(GlobalTally::LEAKAGE) + 1), + xt::range(static_cast(TallyResult::SUM), + static_cast(TallyResult::SUM_SQ) + 1)); + auto gt_reduced = xt::empty_like(gt_view); + read_dataset_lowlevel(file_id, "global_tallies", H5T_NATIVE_DOUBLE, H5S_ALL, + false, gt_reduced.data()); + gt_view = gt_reduced; +} void write_tally_results_nr(hid_t file_id) { @@ -904,7 +933,7 @@ void write_tally_results_nr(hid_t file_id) // Write out global tallies sum and sum_sq if (mpi::master) { - write_dataset(file_id, "global_tallies", gt); + write_global_tallies(file_id); } for (const auto& t : model::tallies) { @@ -977,4 +1006,4 @@ void write_tally_results_nr(hid_t file_id) } } -} // namespace openmc +} // namespace openmc \ No newline at end of file diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index 8cdf3a9050c..b8e5ce6924b 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -70,14 +70,23 @@ vector time_grid; namespace simulation { xt::xtensor_fixed> global_tallies; +xt::xtensor_fixed> + global_tallies_first_gen; int32_t n_realizations {0}; } // namespace simulation double global_tally_absorption; double global_tally_collision; double global_tally_tracklength; +double global_tally_tracklength_sq; double global_tally_leakage; +// Tallies for first generation quantities +double global_tally_absorption_first_gen; +double global_tally_collision_first_gen; +double global_tally_tracklength_first_gen; +double global_tally_tracklength_sq_first_gen; + //============================================================================== // Tally object implementation //============================================================================== @@ -1073,16 +1082,25 @@ void accumulate_tallies() // Accumulate on master only unless run is not reduced then do it on all if (mpi::master || !settings::reduce_tallies) { auto& gt = simulation::global_tallies; + double k_col = 0.0; + double k_abs = 0.0; + double k_tra = 0.0; + + double kq_col = 0.0; + double kq_abs = 0.0; + double kq_tra = 0.0; - if (settings::run_mode == RunMode::EIGENVALUE) { + if (settings::run_mode == RunMode::EIGENVALUE || + (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k)) { if (simulation::current_batch > settings::n_inactive) { // Accumulate products of different estimators of k - double k_col = gt(GlobalTally::K_COLLISION, TallyResult::VALUE) / - simulation::total_weight; - double k_abs = gt(GlobalTally::K_ABSORPTION, TallyResult::VALUE) / - simulation::total_weight; - double k_tra = gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) / - simulation::total_weight; + k_col = gt(GlobalTally::K_COLLISION, TallyResult::VALUE) / + simulation::total_weight; + k_abs = gt(GlobalTally::K_ABSORPTION, TallyResult::VALUE) / + simulation::total_weight; + k_tra = gt(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) / + simulation::total_weight; simulation::k_col_abs += k_col * k_abs; simulation::k_col_tra += k_col * k_tra; simulation::k_abs_tra += k_abs * k_tra; @@ -1090,12 +1108,54 @@ void accumulate_tallies() } // Accumulate results for global tallies + double k_vals[3]; + double kq_vals[3]; for (int i = 0; i < N_GLOBAL_TALLIES; ++i) { double val = gt(i, TallyResult::VALUE) / simulation::total_weight; + if (i < 3) { + k_vals[i] = val; + } gt(i, TallyResult::VALUE) = 0.0; gt(i, TallyResult::SUM) += val; gt(i, TallyResult::SUM_SQ) += val * val; } + + // Accumulate subcritical multiplication tallies + if (settings::run_mode == RunMode::FIXED_SOURCE && + settings::calculate_subcritical_k) { + auto& gt_first_gen = simulation::global_tallies_first_gen; + if (mpi::master || !settings::reduce_tallies) { + if (mpi::master || !settings::reduce_tallies) { + // Accumulate products of different estimators of k + kq_col = gt_first_gen(GlobalTally::K_COLLISION, TallyResult::VALUE) / + simulation::total_weight; + kq_abs = gt_first_gen(GlobalTally::K_ABSORPTION, TallyResult::VALUE) / + simulation::total_weight; + kq_tra = + gt_first_gen(GlobalTally::K_TRACKLENGTH, TallyResult::VALUE) / + simulation::total_weight; + simulation::kq_col_abs += kq_col * kq_abs; + simulation::kq_col_tra += kq_col * kq_tra; + simulation::kq_abs_tra += kq_abs * kq_tra; + } + } + for (int i = 0; i < N_GLOBAL_TALLIES; ++i) { + double val_first_gen = + gt_first_gen(i, TallyResult::VALUE) / simulation::total_weight; + if (i < 3) { + kq_vals[i] = val_first_gen; + } + gt_first_gen(i, TallyResult::VALUE) = 0.0; + gt_first_gen(i, TallyResult::SUM) += val_first_gen; + gt_first_gen(i, TallyResult::SUM_SQ) += val_first_gen * val_first_gen; + } + } + + for (int i = 0; i < simulation::N_K_EST; ++i) { + for (int j = 0; j < simulation::N_K_EST; ++j) { + simulation::k_kq_products[i][j] += k_vals[i] * kq_vals[j]; + } + } } // Accumulate results for each tally @@ -1551,8 +1611,8 @@ extern "C" int openmc_tally_get_n_realizations(int32_t index, int32_t* n) return 0; } -//! \brief Returns a pointer to a tally results array along with its shape. -//! This allows a user to obtain in-memory tally results from Python directly. +//! \brief Returns a pointer to a tally results array along with its shape. This +//! allows a user to obtain in-memory tally results from Python directly. extern "C" int openmc_tally_results( int32_t index, double** results, size_t* shape) { diff --git a/tests/geometry.xml b/tests/geometry.xml new file mode 100644 index 00000000000..b5e29c3bbaa --- /dev/null +++ b/tests/geometry.xml @@ -0,0 +1,65 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 1.0 1.0 +
0.0 0.0 0.0
+ + 9 + 20 10 +19 3 11 + 8 4 +18 2 12 + 7 5 +17 6 13 + 16 14 + 15 + 29 + 40 30 +39 23 31 + 28 24 +38 22 32 + 27 25 +37 26 33 + 36 34 + 35 + + diff --git a/tests/regression_tests/fixed_source_k/__init__.py b/tests/regression_tests/fixed_source_k/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/fixed_source_k/inputs_true.dat b/tests/regression_tests/fixed_source_k/inputs_true.dat new file mode 100644 index 00000000000..faa3b96d9a2 --- /dev/null +++ b/tests/regression_tests/fixed_source_k/inputs_true.dat @@ -0,0 +1,38 @@ + + + + mgxs.h5 + + + + + + + + + + + + fixed source + 100000 + 20 + 5 + + + 0 -1000 -1000 10 1000 1000 + + + 10000000.0 1.0 + + + + false + + multi-group + + false + + true + true + + diff --git a/tests/regression_tests/fixed_source_k/results_true.dat b/tests/regression_tests/fixed_source_k/results_true.dat new file mode 100644 index 00000000000..2ad72e7c6d8 --- /dev/null +++ b/tests/regression_tests/fixed_source_k/results_true.dat @@ -0,0 +1,71 @@ +multiplication: +8.979362E+00 6.965776E-03 +k +8.886335E-01 8.639299E-05 +k_generation: +8.891825E-01 5.167920E-04 +8.891873E-01 5.144271E-04 +8.880890E-01 5.217593E-04 +8.892142E-01 5.171018E-04 +8.877472E-01 5.215507E-04 +8.886940E-01 5.236110E-04 +8.886075E-01 5.127534E-04 +8.886461E-01 5.158635E-04 +8.886604E-01 5.211064E-04 +8.887879E-01 5.161212E-04 +8.895620E-01 5.260243E-04 +8.883862E-01 5.198265E-04 +8.873593E-01 5.108713E-04 +8.886636E-01 5.178553E-04 +8.882488E-01 5.197064E-04 +8.894376E-01 5.177350E-04 +8.891000E-01 5.186993E-04 +8.877115E-01 5.158857E-04 +8.889858E-01 5.174576E-04 +8.884413E-01 5.214052E-04 +ks: +8.267411E-01 1.412098E-04 +ks_generation: +8.282394E-01 1.317049E-03 +8.276995E-01 1.037028E-03 +8.253728E-01 1.051372E-03 +8.270288E-01 1.047493E-03 +8.255261E-01 1.044401E-03 +8.272328E-01 1.051199E-03 +8.267898E-01 1.035525E-03 +8.270472E-01 1.035300E-03 +8.263453E-01 1.051913E-03 +8.274425E-01 1.037003E-03 +8.286009E-01 1.052082E-03 +8.257329E-01 1.048266E-03 +8.232758E-01 1.042139E-03 +8.270810E-01 1.042270E-03 +8.257970E-01 1.049061E-03 +8.276237E-01 1.046188E-03 +8.275681E-01 1.041544E-03 +8.251495E-01 1.040380E-03 +8.264155E-01 1.053247E-03 +8.254596E-01 1.051565E-03 +kq: +1.382496E+00 6.246398E-04 +kq_generation: +1.378182E+00 4.087893E-03 +1.382580E+00 4.123522E-03 +1.385784E+00 4.124623E-03 +1.388341E+00 4.130962E-03 +1.379820E+00 4.111702E-03 +1.379414E+00 4.096006E-03 +1.381743E+00 4.119172E-03 +1.380228E+00 4.064636E-03 +1.386030E+00 4.125097E-03 +1.379049E+00 4.122997E-03 +1.380596E+00 4.101632E-03 +1.387073E+00 4.122476E-03 +1.392196E+00 4.123907E-03 +1.380203E+00 4.110304E-03 +1.384644E+00 4.122543E-03 +1.386710E+00 4.124031E-03 +1.382409E+00 4.135859E-03 +1.382304E+00 4.124394E-03 +1.390039E+00 4.157437E-03 +1.390021E+00 4.107879E-03 diff --git a/tests/regression_tests/fixed_source_k/test.py b/tests/regression_tests/fixed_source_k/test.py new file mode 100644 index 00000000000..b3b019e1a2d --- /dev/null +++ b/tests/regression_tests/fixed_source_k/test.py @@ -0,0 +1,127 @@ +"""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 multiplication. + outstr += 'multiplication:\n' + form = '{0:12.6E} {1:12.6E}\n' + M = sp.multiplication + outstr += form.format(M.n, M.s) + + # 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} {1:12.6E}\n' + k_gen = sp.k_generation + for kg in k_gen: + outstr += form.format(kg.n, kg.s) + + # Write out ks + outstr += 'ks:\n' + form = '{0:12.6E} {1:12.6E}\n' + ks = sp.ks + outstr += form.format(ks.n, ks.s) + + # Write out ks_generation + outstr += 'ks_generation:\n' + form = '{0:12.6E} {1:12.6E}\n' + ks_gen = sp.ks_generation + for ksg in ks_gen: + outstr += form.format(ksg.n, ksg.s) + + # Write out kq. + outstr += 'kq:\n' + form = '{0:12.6E} {1:12.6E}\n' + kq = sp.kq + outstr += form.format(kq.n, kq.s) + + # Write out kq_generation + outstr += 'kq_generation:\n' + form = '{0:12.6E} {1:12.6E}\n' + kq_gen = sp.kq_generation + for kqg in kq_gen: + outstr += form.format(kqg.n, kqg.s) + 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.batches = 20 + model.settings.run_mode = 'fixed source' + model.settings.calculate_subcritical_k = True + model.settings.print_all_k_factors = True + + 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() \ No newline at end of file diff --git a/tests/regression_tests/multiplication/__init__.py b/tests/regression_tests/multiplication/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/multiplication/inputs_true.dat b/tests/regression_tests/multiplication/inputs_true.dat new file mode 100644 index 00000000000..b98a95ebb32 --- /dev/null +++ b/tests/regression_tests/multiplication/inputs_true.dat @@ -0,0 +1,36 @@ + + + + mgxs.h5 + + + + + + + + + + + + eigenvalue + 100000 + 20 + 10 + + + 0 -1000 -1000 10 1000 1000 + + + 10000000.0 1.0 + + + + false + + multi-group + + false + + + diff --git a/tests/regression_tests/multiplication/results_true.dat b/tests/regression_tests/multiplication/results_true.dat new file mode 100644 index 00000000000..b8a2310dfd5 --- /dev/null +++ b/tests/regression_tests/multiplication/results_true.dat @@ -0,0 +1,10 @@ +k-combined: +8.271126E-01 4.868738E-04 +multiplication: +8.963514E+00 4.210376E-02 +k1: +1.378182E+00 4.087893E-03 +kf: +8.269380E-01 1.049151E-03 +g*: +1.663722E+00 1.856058E-02 diff --git a/tests/regression_tests/multiplication/test.py b/tests/regression_tests/multiplication/test.py new file mode 100644 index 00000000000..a39ab9839cd --- /dev/null +++ b/tests/regression_tests/multiplication/test.py @@ -0,0 +1,104 @@ +"""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 multiplication. + outstr += 'multiplication:\n' + form = '{0:12.6E} {1:12.6E}\n' + M = sp.multiplication + outstr += form.format(M.n, M.s) + + # Write out k1. + outstr += 'k1:\n' + form = '{0:12.6E} {1:12.6E}\n' + k1 = sp.k_generation[0] + outstr += form.format(k1.n, k1.s) + + # Write out kf. + outstr += 'kf:\n' + form = '{0:12.6E} {1:12.6E}\n' + kf = 1-k1/(M-1) + outstr += form.format(kf.n, kf.s) + + # Write out g*. + outstr += 'g*:\n' + form = '{0:12.6E} {1:12.6E}\n' + g_star = sp.source_effectiveness + outstr += form.format(g_star.n, g_star.s) + 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.batches = 20 + model.settings.inactive = 10 + + 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() diff --git a/tests/regression_tests/source_effectiveness/__init__.py b/tests/regression_tests/source_effectiveness/__init__.py new file mode 100644 index 00000000000..e69de29bb2d diff --git a/tests/regression_tests/source_effectiveness/inputs_true.dat b/tests/regression_tests/source_effectiveness/inputs_true.dat new file mode 100644 index 00000000000..94e7fdcd963 --- /dev/null +++ b/tests/regression_tests/source_effectiveness/inputs_true.dat @@ -0,0 +1,44 @@ + + + + + + + + + + + + + + + + + eigenvalue + 100000 + 20 + 10 + + + 0.0 0.0 0.0 + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/regression_tests/source_effectiveness/results_true.dat b/tests/regression_tests/source_effectiveness/results_true.dat new file mode 100644 index 00000000000..b4980e9290d --- /dev/null +++ b/tests/regression_tests/source_effectiveness/results_true.dat @@ -0,0 +1,10 @@ +k-combined: +9.960995E-01 5.611810E-04 +multiplication: +3.646446E+02 5.066747E+01 +k1: +1.185417E+00 2.963267E-03 +kf: +9.967402E-01 4.542722E-04 +g*: +1.425744E+00 1.461587E-02 diff --git a/tests/regression_tests/source_effectiveness/test.py b/tests/regression_tests/source_effectiveness/test.py new file mode 100644 index 00000000000..48c70024214 --- /dev/null +++ b/tests/regression_tests/source_effectiveness/test.py @@ -0,0 +1,95 @@ +"""This model is based on a numerical example at doi:10.1016/s0306-4549(98)00048-6 +""" +import openmc +from openmc.stats import spherical_uniform, Point, Mixture, Watt +import numpy as np +import pytest +import os +import glob + +from tests.testing_harness import PyAPITestHarness + + +class SourceTestHarness(PyAPITestHarness): + 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 multiplication. + outstr += 'multiplication:\n' + form = '{0:12.6E} {1:12.6E}\n' + M = sp.multiplication + outstr += form.format(M.n, M.s) + + # Write out k1. + outstr += 'k1:\n' + form = '{0:12.6E} {1:12.6E}\n' + k1 = sp.k_generation[0] + outstr += form.format(k1.n, k1.s) + + # Write out kf. + outstr += 'kf:\n' + form = '{0:12.6E} {1:12.6E}\n' + kf = 1-k1/(M-1) + outstr += form.format(kf.n, kf.s) + + # Write out g*. + outstr += 'g*:\n' + form = '{0:12.6E} {1:12.6E}\n' + g_star = sp.source_effectiveness + outstr += form.format(g_star.n, g_star.s) + return outstr + +@pytest.fixture() +def sources(): + point = Point((0.0,0.0,0.0)) + uniform = spherical_uniform(r_outer=8.85) + cf252_spec = Watt(1.0250e6,2.926e-6) + u235_spec = Watt(0.7747e6,4.852e-6) + u238_spec = Watt(0.6483e6,6.811e-6) + cf252 = openmc.IndependentSource(space=point, + energy=cf252_spec, + strength=100) + u235 = openmc.IndependentSource(space=uniform, + energy=u235_spec, + strength=0.5) + u238 = openmc.IndependentSource(space=uniform, + energy=u238_spec, + strength=57.3) + return [cf252,u235,u238] + + +@pytest.fixture() +def sphere_model(sources): + + mat = openmc.Material(name='heu') + mat.add_element('U',1,percent_type="ao",enrichment=92.2) + mat.set_density('g/cm3',18.6) + + materials = openmc.Materials([mat]) + + + + surface = openmc.Sphere(r=8.85,boundary_type='vacuum') + cell = openmc.Cell(fill=mat,region=-surface) + + root = openmc.Universe() + root.add_cell(cell) + + geometry = openmc.Geometry(root) + + model = openmc.Model(geometry=geometry,materials=materials) + + # Settings + model.settings.particles = 100000 + model.settings.batches = 20 + model.settings.inactive = 10 + model.settings.source = sources + + return model + + +def test_source_effectiveness(sphere_model): + harness = SourceTestHarness("statepoint.20.h5", model=sphere_model) + harness.main()