From 767db7e6a090b88a594be72b4971d71e3c2be243 Mon Sep 17 00:00:00 2001 From: Joffrey Dorville <54550047+JoffreyDorville@users.noreply.github.com> Date: Thu, 25 Sep 2025 14:58:29 -0500 Subject: [PATCH 1/6] Fix IFP implementation (#3580) --- include/openmc/ifp.h | 7 +++---- include/openmc/particle_data.h | 1 + src/ifp.cpp | 4 ++-- src/particle.cpp | 1 + src/physics.cpp | 9 +++------ tests/regression_tests/ifp/results_true.dat | 4 ++-- .../case-03/surface_source_true.h5 | Bin 33344 -> 33344 bytes .../case-04/surface_source_true.h5 | Bin 33344 -> 33344 bytes .../case-05/surface_source_true.h5 | Bin 33344 -> 33344 bytes .../case-07/surface_source_true.h5 | Bin 33344 -> 33344 bytes .../case-08/surface_source_true.h5 | Bin 33344 -> 33344 bytes .../case-09/surface_source_true.h5 | Bin 33344 -> 33344 bytes .../case-10/surface_source_true.h5 | Bin 33344 -> 33344 bytes .../case-12/surface_source_true.h5 | Bin 33344 -> 33344 bytes .../case-13/surface_source_true.h5 | Bin 33344 -> 33344 bytes .../case-14/surface_source_true.h5 | Bin 33344 -> 33344 bytes .../case-a01/surface_source_true.h5 | Bin 6720 -> 6720 bytes .../case-d07/surface_source_true.h5 | Bin 33344 -> 33344 bytes .../case-d08/surface_source_true.h5 | Bin 33344 -> 33344 bytes .../case-e01/surface_source_true.h5 | Bin 33344 -> 33344 bytes .../case-e02/surface_source_true.h5 | Bin 33344 -> 33344 bytes .../case-e03/surface_source_true.h5 | Bin 33344 -> 33344 bytes 22 files changed, 12 insertions(+), 14 deletions(-) diff --git a/include/openmc/ifp.h b/include/openmc/ifp.h index 633a262d5..01904d13c 100644 --- a/include/openmc/ifp.h +++ b/include/openmc/ifp.h @@ -68,15 +68,14 @@ vector _ifp(const T& value, const vector& data) //! //! Add the IFP information in the IFP banks using the same index //! as the one used to append the fission site to the fission bank. +//! The information stored are the delayed group number and lifetime +//! of the neutron that created the fission event. //! Multithreading protection is guaranteed by the index returned by the //! thread_safe_append call in physics.cpp. //! -//! Needs to be done after the delayed group is found. -//! //! \param[in] p Particle -//! \param[in] site Fission site //! \param[in] idx Bank index from the thread_safe_append call in physics.cpp -void ifp(const Particle& p, const SourceSite& site, int64_t idx); +void ifp(const Particle& p, int64_t idx); //! Resize the IFP banks used in the simulation void resize_simulation_ifp_banks(); diff --git a/include/openmc/particle_data.h b/include/openmc/particle_data.h index 1a22f5837..afcd56476 100644 --- a/include/openmc/particle_data.h +++ b/include/openmc/particle_data.h @@ -631,6 +631,7 @@ public: int& event_mt() { return event_mt_; } // MT number of collision const int& event_mt() const { return event_mt_; } int& delayed_group() { return delayed_group_; } // delayed group + const int& delayed_group() const { return delayed_group_; } const int& parent_nuclide() const { return parent_nuclide_; } int& parent_nuclide() { return parent_nuclide_; } // Parent nuclide diff --git a/src/ifp.cpp b/src/ifp.cpp index 1f81f26f6..cc4a76538 100644 --- a/src/ifp.cpp +++ b/src/ifp.cpp @@ -28,13 +28,13 @@ bool is_generation_time_or_both() return false; } -void ifp(const Particle& p, const SourceSite& site, int64_t idx) +void ifp(const Particle& p, int64_t idx) { if (is_beta_effective_or_both()) { const auto& delayed_groups = simulation::ifp_source_delayed_group_bank[p.current_work() - 1]; simulation::ifp_fission_delayed_group_bank[idx] = - _ifp(site.delayed_group, delayed_groups); + _ifp(p.delayed_group(), delayed_groups); } if (is_generation_time_or_both()) { const auto& lifetimes = diff --git a/src/particle.cpp b/src/particle.cpp index f5ad45d80..402af2498 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -144,6 +144,7 @@ void Particle::from_source(const SourceSite* src) time() = src->time; time_last() = src->time; parent_nuclide() = src->parent_nuclide; + delayed_group() = src->delayed_group; // Convert signed surface ID to signed index if (src->surf_id != SURFACE_NONE) { diff --git a/src/physics.cpp b/src/physics.cpp index f667fd586..e947fecbb 100644 --- a/src/physics.cpp +++ b/src/physics.cpp @@ -246,18 +246,15 @@ void create_fission_sites(Particle& p, int i_nuclide, const Reaction& rx) } // Iterated Fission Probability (IFP) method if (settings::ifp_on) { - ifp(p, site, idx); + ifp(p, idx); } } else { p.secondary_bank().push_back(site); } - // Set the delayed group on the particle as well - p.delayed_group() = site.delayed_group; - // Increment the number of neutrons born delayed - if (p.delayed_group() > 0) { - nu_d[p.delayed_group() - 1]++; + if (site.delayed_group > 0) { + nu_d[site.delayed_group - 1]++; } // Write fission particles to nuBank diff --git a/tests/regression_tests/ifp/results_true.dat b/tests/regression_tests/ifp/results_true.dat index 1d8f69e12..466ca1f01 100644 --- a/tests/regression_tests/ifp/results_true.dat +++ b/tests/regression_tests/ifp/results_true.dat @@ -3,7 +3,7 @@ k-combined: tally 1: 9.109384E-08 5.667165E-16 -6.500000E-02 -6.710000E-04 +5.200000E-02 +5.420000E-04 1.489000E+01 1.480036E+01 diff --git a/tests/regression_tests/surface_source_write/case-03/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-03/surface_source_true.h5 index 4dd5f821af92dd98395a8df31af2a81fab8a8ba4..228f5a7d0a462f38a1cb0b4274dd40d4307b5a12 100644 GIT binary patch delta 42 ycmX@m!gQd8X@dwmbDaH^&Eo8yDU2+W6SL1VvTZiZaTj7_nw%&gzj;R83T6N=*$u}4 delta 48 zcmX@m!gQd8X@dwm^EJ&=o5k5ZQzln%3QR7@J}|k1S75V4j)f42F}XlMVDpT+70dv* Csu83B diff --git a/tests/regression_tests/surface_source_write/case-04/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-04/surface_source_true.h5 index ee00af3895c349d60d786af1aedef39e05e24d20..c276af40f4302f4e0feff883bc3c472bc4df7ca1 100644 GIT binary patch delta 38 ucmX@m!gQd8X@dwmbDaH^&Eo8xzKkrB7lxc?WZPUAI-ie`X|rI8z8nArx(svx delta 44 xcmX@m!gQd8X@dwm^L5Qro5k5XeJ59N3QXP*a$vFtm%!$R&|0h5VN3AkPef$3?Q>mU+^N7 CA`&M6 diff --git a/tests/regression_tests/surface_source_write/case-12/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-12/surface_source_true.h5 index ca5eda8e2f347a1fc88b05db6690ad5b6d3a2a51..6c6925daba145c5db52a281023dd017f066e0732 100644 GIT binary patch delta 117 zcmX@m!gQd8X@dwmbAtVq&Eo8lQj<>z2yAXpOweRx+5FH;oYOG=2mfb3*B|@WHpD(y z^yT0FTL*$W*`i%Qt;HW$zbo>Udt8lS9X5&5&?bMwZEqk?BA) UUPkiE0XaX*)4?VeR9xT$0J0-9HUIzs delta 124 zcmV-?0E7R)gaW{X0a}p%0E$1%ExYoLO_SJDO6Vllc)#~laN_h elc)#~qxe}KlPD=5lYkHilOT9N0RppOc;5**VKsyR diff --git a/tests/regression_tests/surface_source_write/case-13/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-13/surface_source_true.h5 index ca5eda8e2f347a1fc88b05db6690ad5b6d3a2a51..6c6925daba145c5db52a281023dd017f066e0732 100644 GIT binary patch delta 117 zcmX@m!gQd8X@dwmbAtVq&Eo8lQj<>z2yAXpOweRx+5FH;oYOG=2mfb3*B|@WHpD(y z^yT0FTL*$W*`i%Qt;HW$zbo>Udt8lS9X5&5&?bMwZEqk?BA) UUPkiE0XaX*)4?VeR9xT$0J0-9HUIzs delta 124 zcmV-?0E7R)gaW{X0a}p%0E$1%ExYoLO_SJDO6Vllc)#~laN_h elc)#~qxe}KlPD=5lYkHilOT9N0RppOc;5**VKsyR diff --git a/tests/regression_tests/surface_source_write/case-14/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-14/surface_source_true.h5 index c9900bd9c0469d3b171f863f0518d6c94e966027..6c6925daba145c5db52a281023dd017f066e0732 100644 GIT binary patch delta 117 zcmX@m!gQd8X@dwmbAtVq&Eo8lQj<>z2yAXpOweRx+5FH;oYOG=2mfb3*B|@WHpD(y z^yT0FTL*$W*`i%Qt;HW$zbo>Udt8lS9X5&5&?bMwZEqk?BA) UUPkiE0XaX*)4?VeR9xT$0J0-9HUIzs delta 124 zcmV-?0E7R)gaW{X0a}p%0E$1%ExYoLO_SJDO6Vllc)#~laN_h elc)#~qxe}KlPD=5lYkHilOT9N0RppOc;5**(=~+v diff --git a/tests/regression_tests/surface_source_write/case-a01/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-a01/surface_source_true.h5 index 8dc148ee07eb33d6374387e557b35914f5fa2003..da36fc505eb67f3da0ce48511464b699a4fb16e3 100644 GIT binary patch delta 70 zcmX?La=>JR2s?9<{gut)?19XaFK}f{&fv<}Y{8wtJxM%c@(i(z$rph10tuivkhTEQ W1(F$)IV3$8StcLklb?J;G6MjFG#L5- delta 68 zcmX?La=>JR2s`s*%~PAj*#ns;U*O7^oWYf`*@8QPdy;s@X>K%LUYL4oBy*u_xg7`DD%%1Q6vVXpy?7Ij1?@Ya&Zg=K`&foCGUf}-ZSJHay^ESNO9MxAl*L;*XO>X>K%KZYL4oBy*u{Hg7`DD%%1Q6ynnu+?7Ij1?@Ya&ZgX>K%LUYL4oBy*u_xg7`DD%%1Q6vVXpy?7Ij1?@Ya&Zg=K`&foCGUf}-ZSJHay^ESNO9MxAl*L;*XO>X>K%KZYL4oBy*u{Hg7`DD%%1Q6ynnu+?7Ij1?@Ya&ZgLSiQX diff --git a/tests/regression_tests/surface_source_write/case-e01/surface_source_true.h5 b/tests/regression_tests/surface_source_write/case-e01/surface_source_true.h5 index 0dda15ed2192739435e8f50e126831021d620f53..bbfbd152bc6543871750cb9bd4cc55fe0899754a 100644 GIT binary patch delta 39 vcmX@m!gQd8X@dwm^L6Jdo5k5Vtr*!hD>|xkF|us-^iEf2WZJy3sGlDI{ml#e delta 45 zcmV+|0Mh@!gaW{X0|xkF|us-^iEf2WZJy3sGlDI{ml#e delta 45 zcmV+|0Mh@!gaW{X0|fgu z`(V+RfBSD82<~KyeqpaT&*);s!pHXSSC}u~^zD?rU%;v3d2LP(9h)bH)H5-%Y!(bR wXPmsFL|}77#0D`CW3ohk#N->+5}P~ne=vh(CKPKhf^3?6!Bk-Ljgkyr09VyO4FCWD delta 148 zcmX@m!gQd8X@dwm^Bj*eo5k5BRTx<|8*2A+P4 Date: Fri, 26 Sep 2025 12:27:05 -0400 Subject: [PATCH 2/6] Multi-group capability for kinetics parameter calculations with Iterated Fission Probability (#3425) Co-authored-by: GuySten Co-authored-by: Paul Romano --- docs/source/usersguide/kinetics.rst | 25 +++++++- openmc/model/model.py | 41 +++++++++++-- openmc/statepoint.py | 58 +++++++++++++++++++ src/tallies/tally.cpp | 3 +- src/tallies/tally_scoring.cpp | 9 +++ .../ifp/groupwise/__init__.py | 0 .../ifp/groupwise/inputs_true.dat | 43 ++++++++++++++ .../ifp/groupwise/results_true.dat | 21 +++++++ tests/regression_tests/ifp/groupwise/test.py | 40 +++++++++++++ tests/regression_tests/ifp/total/__init__.py | 0 .../ifp/{ => total}/inputs_true.dat | 0 .../ifp/{ => total}/results_true.dat | 0 .../regression_tests/ifp/{ => total}/test.py | 1 - tests/unit_tests/test_ifp.py | 39 +++++++++++++ 14 files changed, 272 insertions(+), 8 deletions(-) create mode 100644 tests/regression_tests/ifp/groupwise/__init__.py create mode 100644 tests/regression_tests/ifp/groupwise/inputs_true.dat create mode 100644 tests/regression_tests/ifp/groupwise/results_true.dat create mode 100644 tests/regression_tests/ifp/groupwise/test.py create mode 100644 tests/regression_tests/ifp/total/__init__.py rename tests/regression_tests/ifp/{ => total}/inputs_true.dat (100%) rename tests/regression_tests/ifp/{ => total}/results_true.dat (100%) rename tests/regression_tests/ifp/{ => total}/test.py (99%) diff --git a/docs/source/usersguide/kinetics.rst b/docs/source/usersguide/kinetics.rst index bdf26d341..9024ff822 100644 --- a/docs/source/usersguide/kinetics.rst +++ b/docs/source/usersguide/kinetics.rst @@ -67,6 +67,23 @@ are needed to compute kinetics parameters in OpenMC: Obtaining kinetics parameters ----------------------------- +The ``Model`` class can be used to automatically generate all IFP tallies using +the Python API with :attr:`openmc.Settings.ifp_n_generation` greater than 0 and +the :meth:`openmc.Model.add_ifp_kinetics_tallies` method:: + + model = openmc.Model(geometry, settings=settings) + model.add_kinetics_parameters_tallies(num_groups=6) # Add 6 precursor groups + +Alternatively, each of the tallies can be manually defined using group-wise or +total :math:`\beta_{\text{eff}}` specified by providing a 6-group +:class:`openmc.DelayedGroupFilter`:: + + beta_tally = openmc.Tally(name="group-beta-score") + beta_tally.scores = ["ifp-beta-numerator"] + + # Add DelayedGroupFilter to enable group-wise tallies + beta_tally.filters = [openmc.DelayedGroupFilter(list(range(1, 7)))] + Here is an example showing how to declare the three available IFP scores in a single tally:: @@ -95,6 +112,12 @@ for ``ifp-denominator``: \beta_{\text{eff}} = \frac{S_{\text{ifp-beta-numerator}}}{S_{\text{ifp-denominator}}} +The kinetics parameters can be retrieved directly from a statepoint file using +the :meth:`openmc.StatePoint.ifp_results` method:: + + with openmc.StatePoint(output_path) as sp: + generation_time, beta_eff = sp.get_kinetics_parameters() + .. only:: html .. rubric:: References @@ -107,4 +130,4 @@ for ``ifp-denominator``: of the Iterated Fission Probability Method in OpenMC to Compute Adjoint-Weighted Kinetics Parameters", International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering (M&C 2025), Denver, April 27-30, - 2025 (to be presented). + 2025. diff --git a/openmc/model/model.py b/openmc/model/model.py index 59e6fa511..c1ffafafd 100644 --- a/openmc/model/model.py +++ b/openmc/model/model.py @@ -1,7 +1,7 @@ from __future__ import annotations from collections.abc import Iterable, Sequence import copy -from functools import lru_cache +from functools import cache from pathlib import Path import math from numbers import Integral, Real @@ -160,7 +160,7 @@ class Model: return False @property - @lru_cache(maxsize=None) + @cache def _materials_by_id(self) -> dict: """Dictionary mapping material ID --> material""" if self.materials: @@ -170,14 +170,14 @@ class Model: return {mat.id: mat for mat in mats} @property - @lru_cache(maxsize=None) + @cache def _cells_by_id(self) -> dict: """Dictionary mapping cell ID --> cell""" cells = self.geometry.get_all_cells() return {cell.id: cell for cell in cells.values()} @property - @lru_cache(maxsize=None) + @cache def _cells_by_name(self) -> dict[int, openmc.Cell]: # Get the names maps, but since names are not unique, store a set for # each name key. In this way when the user requests a change by a name, @@ -190,7 +190,7 @@ class Model: return result @property - @lru_cache(maxsize=None) + @cache def _materials_by_name(self) -> dict[int, openmc.Material]: if self.materials is None: mats = self.geometry.get_all_materials().values() @@ -203,6 +203,37 @@ class Model: result[mat.name].add(mat) return result + def add_kinetics_parameters_tallies(self, num_groups: int | None = None): + """Add tallies for calculating kinetics parameters using the IFP method. + + This method adds tallies to the model for calculating two kinetics + parameters, the generation time and the effective delayed neutron + fraction (beta effective). After a model is run, these parameters can be + determined through the :meth:`openmc.StatePoint.ifp_results` method. + + Parameters + ---------- + num_groups : int, optional + Number of precursor groups to filter the delayed neutron fraction. + If None, only the total effective delayed neutron fraction is + tallied. + + """ + if not any('ifp-time-numerator' in t.scores for t in self.tallies): + gen_time_tally = openmc.Tally(name='IFP time numerator') + gen_time_tally.scores = ['ifp-time-numerator'] + self.tallies.append(gen_time_tally) + if not any('ifp-beta-numerator' in t.scores for t in self.tallies): + beta_tally = openmc.Tally(name='IFP beta numerator') + beta_tally.scores = ['ifp-beta-numerator'] + if num_groups is not None: + beta_tally.filters = [openmc.DelayedGroupFilter(list(range(1, num_groups + 1)))] + self.tallies.append(beta_tally) + if not any('ifp-denominator' in t.scores for t in self.tallies): + denom_tally = openmc.Tally(name='IFP denominator') + denom_tally.scores = ['ifp-denominator'] + self.tallies.append(denom_tally) + @classmethod def from_xml( cls, diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 715becf48..29c11921c 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -1,4 +1,5 @@ from datetime import datetime +from collections import namedtuple import glob import re import os @@ -8,6 +9,7 @@ import h5py import numpy as np from pathlib import Path from uncertainties import ufloat +from uncertainties.unumpy import uarray import openmc import openmc.checkvalue as cv @@ -15,6 +17,9 @@ import openmc.checkvalue as cv _VERSION_STATEPOINT = 18 +KineticsParameters = namedtuple("KineticsParameters", ["generation_time", "beta_effective"]) + + class StatePoint: """State information on a simulation at a certain point in time (at the end of a given batch). Statepoints can be used to analyze tally results as well @@ -710,3 +715,56 @@ class StatePoint: tally_filter.paths = cell.paths self._summary = summary + + def get_kinetics_parameters(self) -> KineticsParameters: + """Get kinetics parameters from IFP tallies. + + This method searches the tallies in the statepoint for the tallies + required to compute kinetics parameters using the Iterated Fission + Probability (IFP) method. + + Returns + ------- + KineticsParameters + A named tuple containing the generation time and effective delayed + neutron fraction. If the necessary tallies for one or both + parameters are not found, that parameter is returned as None. + + """ + + denom_tally = None + gen_time_tally = None + beta_tally = None + for tally in self.tallies.values(): + if 'ifp-denominator' in tally.scores: + denom_tally = self.get_tally(scores=['ifp-denominator']) + if 'ifp-time-numerator' in tally.scores: + gen_time_tally = self.get_tally(scores=['ifp-time-numerator']) + if 'ifp-beta-numerator' in tally.scores: + beta_tally = self.get_tally(scores=['ifp-beta-numerator']) + + if denom_tally is None: + return KineticsParameters(None, None) + + def get_ufloat(tally, score): + return uarray(tally.get_values(scores=[score]), + tally.get_values(scores=[score], value='std_dev')) + + denom_values = get_ufloat(denom_tally, 'ifp-denominator') + if gen_time_tally is None: + generation_time = None + else: + gen_time_values = get_ufloat(gen_time_tally, 'ifp-time-numerator') + gen_time_values /= denom_values*self.keff + generation_time = gen_time_values.flatten()[0] + + if beta_tally is None: + beta_effective = None + else: + beta_values = get_ufloat(beta_tally, 'ifp-beta-numerator') + beta_values /= denom_values + beta_effective = beta_values.flatten() + if beta_effective.size == 1: + beta_effective = beta_effective[0] + + return KineticsParameters(generation_time, beta_effective) diff --git a/src/tallies/tally.cpp b/src/tallies/tally.cpp index ae0bffe6e..b9c615ecb 100644 --- a/src/tallies/tally.cpp +++ b/src/tallies/tally.cpp @@ -560,7 +560,8 @@ void Tally::set_scores(const vector& scores) // Make sure a delayed group filter wasn't used with an incompatible // score. if (delayedgroup_filter_ != C_NONE) { - if (score_str != "delayed-nu-fission" && score_str != "decay-rate") + if (score_str != "delayed-nu-fission" && score_str != "decay-rate" && + score_str != "ifp-beta-numerator") fatal_error("Cannot tally " + score_str + "with a delayedgroup filter"); } diff --git a/src/tallies/tally_scoring.cpp b/src/tallies/tally_scoring.cpp index 0df80a239..67e851644 100644 --- a/src/tallies/tally_scoring.cpp +++ b/src/tallies/tally_scoring.cpp @@ -964,6 +964,15 @@ void score_general_ce_nonanalog(Particle& p, int i_tally, int start_index, if (delayed_groups.size() == settings::ifp_n_generation) { if (delayed_groups[0] > 0) { score = p.wgt_last(); + if (tally.delayedgroup_filter_ != C_NONE) { + auto i_dg_filt = tally.filters()[tally.delayedgroup_filter_]; + const DelayedGroupFilter& filt { + *dynamic_cast( + model::tally_filters[i_dg_filt].get())}; + score_fission_delayed_dg(i_tally, delayed_groups[0] - 1, + score, score_index, p.filter_matches()); + continue; + } } } } diff --git a/tests/regression_tests/ifp/groupwise/__init__.py b/tests/regression_tests/ifp/groupwise/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/ifp/groupwise/inputs_true.dat b/tests/regression_tests/ifp/groupwise/inputs_true.dat new file mode 100644 index 000000000..6d7e20717 --- /dev/null +++ b/tests/regression_tests/ifp/groupwise/inputs_true.dat @@ -0,0 +1,43 @@ + + + + + + + + + + + + + + eigenvalue + 1000 + 20 + 5 + + + -10.0 -10.0 -10.0 10.0 10.0 10.0 + + + true + + + 5 + + + + 1 2 3 4 5 6 + + + ifp-time-numerator + + + 1 + ifp-beta-numerator + + + ifp-denominator + + + diff --git a/tests/regression_tests/ifp/groupwise/results_true.dat b/tests/regression_tests/ifp/groupwise/results_true.dat new file mode 100644 index 000000000..ea66a8de3 --- /dev/null +++ b/tests/regression_tests/ifp/groupwise/results_true.dat @@ -0,0 +1,21 @@ +k-combined: +1.006559E+00 5.389391E-03 +tally 1: +9.109384E-08 +5.667165E-16 +tally 2: +3.000000E-03 +9.000000E-06 +0.000000E+00 +0.000000E+00 +2.100000E-02 +1.370000E-04 +2.800000E-02 +2.220000E-04 +0.000000E+00 +0.000000E+00 +0.000000E+00 +0.000000E+00 +tally 3: +1.489000E+01 +1.480036E+01 diff --git a/tests/regression_tests/ifp/groupwise/test.py b/tests/regression_tests/ifp/groupwise/test.py new file mode 100644 index 000000000..a1a0ebefb --- /dev/null +++ b/tests/regression_tests/ifp/groupwise/test.py @@ -0,0 +1,40 @@ +"""Test the Iterated Fission Probability (IFP) method to compute adjoint-weighted +kinetics parameters using dedicated tallies.""" + +import openmc +import pytest + +from tests.testing_harness import PyAPITestHarness + +@pytest.fixture() +def ifp_model(): + # Material + material = openmc.Material(name="core") + material.add_nuclide("U235", 1.0) + material.set_density('g/cm3', 16.0) + + # Geometry + radius = 10.0 + sphere = openmc.Sphere(r=radius, boundary_type="vacuum") + cell = openmc.Cell(region=-sphere, fill=material) + geometry = openmc.Geometry([cell]) + + # Settings + settings = openmc.Settings() + settings.particles = 1000 + settings.batches = 20 + settings.inactive = 5 + settings.ifp_n_generation = 5 + + model = openmc.Model(settings=settings, geometry=geometry) + + space = openmc.stats.Box(*cell.bounding_box) + model.settings.source = openmc.IndependentSource( + space=space, constraints={'fissionable': True}) + model.add_kinetics_parameters_tallies(num_groups=6) + return model + + +def test_iterated_fission_probability(ifp_model): + harness = PyAPITestHarness("statepoint.20.h5", model=ifp_model) + harness.main() diff --git a/tests/regression_tests/ifp/total/__init__.py b/tests/regression_tests/ifp/total/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/ifp/inputs_true.dat b/tests/regression_tests/ifp/total/inputs_true.dat similarity index 100% rename from tests/regression_tests/ifp/inputs_true.dat rename to tests/regression_tests/ifp/total/inputs_true.dat diff --git a/tests/regression_tests/ifp/results_true.dat b/tests/regression_tests/ifp/total/results_true.dat similarity index 100% rename from tests/regression_tests/ifp/results_true.dat rename to tests/regression_tests/ifp/total/results_true.dat diff --git a/tests/regression_tests/ifp/test.py b/tests/regression_tests/ifp/total/test.py similarity index 99% rename from tests/regression_tests/ifp/test.py rename to tests/regression_tests/ifp/total/test.py index 6969a54c4..18b89cfc0 100644 --- a/tests/regression_tests/ifp/test.py +++ b/tests/regression_tests/ifp/total/test.py @@ -6,7 +6,6 @@ import pytest from tests.testing_harness import PyAPITestHarness - @pytest.fixture() def ifp_model(): model = openmc.Model() diff --git a/tests/unit_tests/test_ifp.py b/tests/unit_tests/test_ifp.py index 8d0fd9801..e527f1624 100644 --- a/tests/unit_tests/test_ifp.py +++ b/tests/unit_tests/test_ifp.py @@ -47,3 +47,42 @@ def test_exceptions(options, error, run_in_tmpdir, geometry): tallies = openmc.Tallies([tally]) model = openmc.Model(geometry=geometry, settings=settings, tallies=tallies) model.run() + + +@pytest.mark.parametrize( + "num_groups, use_auto_tallies", + [ + (None, True), + (None, False), + (6, True), + (6, False), + ], +) +def test_get_kinetics_parameters(run_in_tmpdir, geometry, num_groups, use_auto_tallies): + # Create basic model + model = openmc.Model(geometry=geometry) + model.settings.particles = 1000 + model.settings.batches = 20 + model.settings.inactive = 5 + model.settings.ifp_n_generation = 5 + + # Add IFP tallies either via the convenience method or manually + if use_auto_tallies: + model.add_kinetics_parameters_tallies(num_groups=num_groups) + else: + for score in ["ifp-time-numerator", "ifp-beta-numerator", "ifp-denominator"]: + tally = openmc.Tally() + tally.scores = [score] + if score == "ifp-beta-numerator" and num_groups is not None: + tally.filters = [openmc.DelayedGroupFilter(list(range(1, num_groups + 1)))] + model.tallies.append(tally) + + # Run and get kinetics parameters + sp_file = model.run() + with openmc.StatePoint(sp_file) as sp: + params = sp.get_kinetics_parameters() + assert isinstance(params, openmc.KineticsParameters) + assert params.generation_time is not None + assert params.beta_effective is not None + if num_groups is not None: + assert len(params.beta_effective) == num_groups From 4011b7a5515e4098f179d825313f574680d3a68f Mon Sep 17 00:00:00 2001 From: Perry <100789850+yrrepy@users.noreply.github.com> Date: Fri, 26 Sep 2025 09:46:23 -0700 Subject: [PATCH 3/6] Optional separation of mesh-material-volume calc from get_homogenized_materials (#3581) Co-authored-by: Paul Romano --- openmc/mesh.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/openmc/mesh.py b/openmc/mesh.py index 2e9abd1b6..9601207e9 100644 --- a/openmc/mesh.py +++ b/openmc/mesh.py @@ -287,6 +287,7 @@ class MeshBase(IDManagerMixin, ABC): model: openmc.Model, n_samples: int | tuple[int, int, int] = 10_000, include_void: bool = True, + material_volumes: MeshMaterialVolumes | None = None, **kwargs ) -> list[openmc.Material]: """Generate homogenized materials over each element in a mesh. @@ -305,8 +306,12 @@ class MeshBase(IDManagerMixin, ABC): the x, y, and z dimensions. include_void : bool, optional Whether homogenization should include voids. + material_volumes : MeshMaterialVolumes, optional + Previously computed mesh material volumes to use for homogenization. + If not provided, they will be computed by calling + :meth:`material_volumes`. **kwargs - Keyword-arguments passed to :meth:`MeshBase.material_volumes`. + Keyword-arguments passed to :meth:`material_volumes`. Returns ------- @@ -314,7 +319,10 @@ class MeshBase(IDManagerMixin, ABC): Homogenized material in each mesh element """ - vols = self.material_volumes(model, n_samples, **kwargs) + if material_volumes is None: + vols = self.material_volumes(model, n_samples, **kwargs) + else: + vols = material_volumes mat_volume_by_element = [vols.by_element(i) for i in range(vols.num_elements)] # Create homogenized material for each element @@ -424,7 +432,6 @@ class MeshBase(IDManagerMixin, ABC): # Restore original tallies model.tallies = original_tallies - return volumes From feefcc671307fc84e7357430f5f07343503919e5 Mon Sep 17 00:00:00 2001 From: Jonathan Shimwell Date: Wed, 1 Oct 2025 00:09:15 +0200 Subject: [PATCH 4/6] Adding tally filter type option to statepoint get_tally (#3584) Co-authored-by: Jon Shimwell Co-authored-by: GuySten <62616591+GuySten@users.noreply.github.com> --- openmc/statepoint.py | 9 +++- tests/unit_tests/test_statepoint.py | 65 +++++++++++++++++++++++++++++ 2 files changed, 73 insertions(+), 1 deletion(-) create mode 100644 tests/unit_tests/test_statepoint.py diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 29c11921c..a763db397 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -536,7 +536,7 @@ class StatePoint: def get_tally(self, scores=[], filters=[], nuclides=[], name=None, id=None, estimator=None, exact_filters=False, exact_nuclides=False, exact_scores=False, - multiply_density=None, derivative=None): + multiply_density=None, derivative=None, filter_type=None): """Finds and returns a Tally object with certain properties. This routine searches the list of Tallies and returns the first Tally @@ -580,6 +580,9 @@ class StatePoint: to the same value as this parameter. derivative : openmc.TallyDerivative, optional TallyDerivative object to match. + filter_type : type, optional + If not None, the Tally must have at least one Filter that is an + instance of this type. For example `openmc.MeshFilter`. Returns ------- @@ -653,6 +656,10 @@ class StatePoint: if not contains_filters: continue + if filter_type is not None: + if not any(isinstance(f, filter_type) for f in test_tally.filters): + continue + # Determine if Tally has the queried Nuclide(s) if nuclides: if not all(nuclide in test_tally.nuclides for nuclide in nuclides): diff --git a/tests/unit_tests/test_statepoint.py b/tests/unit_tests/test_statepoint.py new file mode 100644 index 000000000..7ffaf7ec2 --- /dev/null +++ b/tests/unit_tests/test_statepoint.py @@ -0,0 +1,65 @@ +import openmc + + +def test_get_tally_filter_type(run_in_tmpdir): + """Test various ways of retrieving tallies from a StatePoint object.""" + + mat = openmc.Material() + mat.add_nuclide("H1", 1.0) + mat.set_density("g/cm3", 10.0) + + sphere = openmc.Sphere(r=10.0, boundary_type="vacuum") + cell = openmc.Cell(fill=mat, region=-sphere) + geometry = openmc.Geometry([cell]) + + settings = openmc.Settings() + settings.particles = 10 + settings.batches = 2 + settings.run_mode = "fixed source" + + reg_mesh = openmc.RegularMesh().from_domain(cell) + tally1 = openmc.Tally(tally_id=1) + mesh_filter = openmc.MeshFilter(reg_mesh) + tally1.filters = [mesh_filter] + tally1.scores = ["flux"] + + tally2 = openmc.Tally(tally_id=2, name="heating tally") + cell_filter = openmc.CellFilter(cell) + tally2.filters = [cell_filter] + tally2.scores = ["heating"] + + tallies = openmc.Tallies([tally1, tally2]) + model = openmc.Model( + geometry=geometry, materials=[mat], settings=settings, tallies=tallies + ) + + sp_filename = model.run() + + sp = openmc.StatePoint(sp_filename) + + tally_found = sp.get_tally(filter_type=openmc.MeshFilter) + assert tally_found.id == 1 + + tally_found = sp.get_tally(filter_type=openmc.CellFilter) + assert tally_found.id == 2 + + tally_found = sp.get_tally(filters=[mesh_filter]) + assert tally_found.id == 1 + + tally_found = sp.get_tally(filters=[cell_filter]) + assert tally_found.id == 2 + + tally_found = sp.get_tally(scores=["heating"]) + assert tally_found.id == 2 + + tally_found = sp.get_tally(name="heating tally") + assert tally_found.id == 2 + + tally_found = sp.get_tally(name=None) + assert tally_found.id == 1 + + tally_found = sp.get_tally(id=1) + assert tally_found.id == 1 + + tally_found = sp.get_tally(id=2) + assert tally_found.id == 2 From 3ac64d9a01d8d10ff126d62d8a9f618fbd803f4b Mon Sep 17 00:00:00 2001 From: John Tramm Date: Thu, 2 Oct 2025 11:02:43 -0500 Subject: [PATCH 5/6] Random Ray Base Source Region Refactor (#3576) --- include/openmc/constants.h | 5 + .../openmc/random_ray/flat_source_domain.h | 47 ++- .../openmc/random_ray/linear_source_domain.h | 2 +- include/openmc/random_ray/random_ray.h | 1 - .../openmc/random_ray/random_ray_simulation.h | 9 +- include/openmc/random_ray/source_region.h | 1 - src/random_ray/flat_source_domain.cpp | 373 ++++++++++-------- src/random_ray/linear_source_domain.cpp | 47 +-- src/random_ray/random_ray.cpp | 121 +++--- src/random_ray/random_ray_simulation.cpp | 190 ++++----- src/random_ray/source_region.cpp | 28 +- src/settings.cpp | 1 - .../random_ray_low_density/__init__.py | 0 .../random_ray_low_density/inputs_true.dat | 244 ++++++++++++ .../random_ray_low_density/results_true.dat | 9 + .../random_ray_low_density/test.py | 60 +++ .../results_true.dat | 8 +- 17 files changed, 700 insertions(+), 446 deletions(-) create mode 100644 tests/regression_tests/random_ray_low_density/__init__.py create mode 100644 tests/regression_tests/random_ray_low_density/inputs_true.dat create mode 100644 tests/regression_tests/random_ray_low_density/results_true.dat create mode 100644 tests/regression_tests/random_ray_low_density/test.py diff --git a/include/openmc/constants.h b/include/openmc/constants.h index df13da370..a0d164613 100644 --- a/include/openmc/constants.h +++ b/include/openmc/constants.h @@ -68,6 +68,11 @@ constexpr double MIN_HITS_PER_BATCH {1.5}; // prevent extremely large adjoint source terms from being generated. constexpr double ZERO_FLUX_CUTOFF {1e-22}; +// The minimum macroscopic cross section value considered non-void for the +// random ray solver. Materials with any group with a cross section below this +// value will be converted to pure void. +constexpr double MINIMUM_MACRO_XS {1e-6}; + // ============================================================================ // MATH AND PHYSICAL CONSTANTS diff --git a/include/openmc/random_ray/flat_source_domain.h b/include/openmc/random_ray/flat_source_domain.h index 78351fcc5..d4e802734 100644 --- a/include/openmc/random_ray/flat_source_domain.h +++ b/include/openmc/random_ray/flat_source_domain.h @@ -27,8 +27,9 @@ public: //---------------------------------------------------------------------------- // Methods - virtual void update_neutron_source(double k_eff); - double compute_k_eff(double k_eff_old) const; + virtual void update_single_neutron_source(SourceRegionHandle& srh); + virtual void update_all_neutron_sources(); + void compute_k_eff(); virtual void normalize_scalar_flux_and_volumes( double total_active_distance_per_iteration); @@ -41,7 +42,7 @@ public: void output_to_vtk() const; void convert_external_sources(); void count_external_source_regions(); - void set_adjoint_sources(const vector& forward_flux); + void set_adjoint_sources(); void flux_swap(); virtual double evaluate_flux_at_point(Position r, int64_t sr, int g) const; double compute_fixed_source_normalization_factor() const; @@ -54,9 +55,8 @@ public: bool is_target_void); void apply_mesh_to_cell_and_children(int32_t i_cell, int32_t mesh_idx, int32_t target_material_id, bool is_target_void); - void prepare_base_source_regions(); SourceRegionHandle get_subdivided_source_region_handle( - int64_t sr, int mesh_bin, Position r, double dist, Direction u); + SourceRegionKey sr_key, Position r, Direction u); void finalize_discovered_source_regions(); void apply_transport_stabilization(); int64_t n_source_regions() const @@ -67,6 +67,10 @@ public: { return source_regions_.n_source_regions() * negroups_; } + int64_t lookup_base_source_region_idx(const GeometryState& p) const; + SourceRegionKey lookup_source_region_key(const GeometryState& p) const; + int64_t lookup_mesh_bin(int64_t sr, Position r) const; + int lookup_mesh_idx(int64_t sr) const; //---------------------------------------------------------------------------- // Static Data members @@ -86,6 +90,7 @@ public: //---------------------------------------------------------------------------- // Public Data members + double k_eff_ {1.0}; // Eigenvalue bool mapped_all_tallies_ {false}; // If all source regions have been visited int64_t n_external_source_regions_ {0}; // Total number of source regions with @@ -110,14 +115,6 @@ public: // The abstract container holding all source region-specific data SourceRegionContainer source_regions_; - // Base source region container. When source region subdivision via mesh - // is in use, this container holds the original (non-subdivided) material - // filled cell instance source regions. These are useful as they can be - // initialized with external source and mesh domain information ahead of time. - // Then, dynamically discovered source regions can be initialized by cloning - // their base region. - SourceRegionContainer base_source_regions_; - // Parallel hash map holding all source regions discovered during // a single iteration. This is a threadsafe data structure that is cleaned // out after each iteration and stored in the "source_regions_" container. @@ -134,8 +131,17 @@ public: // Map that relates a SourceRegionKey to the external source index. This map // is used to check if there are any point sources within a subdivided source // region at the time it is discovered. - std::unordered_map - point_source_map_; + std::unordered_map, SourceRegionKey::HashFunctor> + external_point_source_map_; + + // Map that relates a base source region index to the external source index. + // This map is used to check if there are any volumetric sources within a + // subdivided source region at the time it is discovered. + std::unordered_map> external_volumetric_source_map_; + + // Map that relates a base source region index to a mesh index. This map + // is used to check which subdivision mesh is present in a source region. + std::unordered_map mesh_map_; // If transport corrected MGXS data is being used, there may be negative // in-group scattering cross sections that can result in instability in MOC @@ -147,12 +153,11 @@ protected: //---------------------------------------------------------------------------- // Methods void apply_external_source_to_source_region( - Discrete* discrete, double strength_factor, SourceRegionHandle& srh); - void apply_external_source_to_cell_instances(int32_t i_cell, - Discrete* discrete, double strength_factor, int target_material_id, - const vector& instances); - void apply_external_source_to_cell_and_children(int32_t i_cell, - Discrete* discrete, double strength_factor, int32_t target_material_id); + int src_idx, SourceRegionHandle& srh); + void apply_external_source_to_cell_instances(int32_t i_cell, int src_idx, + int target_material_id, const vector& instances); + void apply_external_source_to_cell_and_children( + int32_t i_cell, int src_idx, int32_t target_material_id); virtual void set_flux_to_flux_plus_source(int64_t sr, double volume, int g); void set_flux_to_source(int64_t sr, int g); virtual void set_flux_to_old_flux(int64_t sr, int g); diff --git a/include/openmc/random_ray/linear_source_domain.h b/include/openmc/random_ray/linear_source_domain.h index 67fdd99f8..0098c7820 100644 --- a/include/openmc/random_ray/linear_source_domain.h +++ b/include/openmc/random_ray/linear_source_domain.h @@ -20,7 +20,7 @@ class LinearSourceDomain : public FlatSourceDomain { public: //---------------------------------------------------------------------------- // Methods - void update_neutron_source(double k_eff) override; + void update_single_neutron_source(SourceRegionHandle& srh) override; void normalize_scalar_flux_and_volumes( double total_active_distance_per_iteration) override; diff --git a/include/openmc/random_ray/random_ray.h b/include/openmc/random_ray/random_ray.h index abf2a2688..40c67ef95 100644 --- a/include/openmc/random_ray/random_ray.h +++ b/include/openmc/random_ray/random_ray.h @@ -48,7 +48,6 @@ public: static double distance_active_; // Active ray length static unique_ptr ray_source_; // Starting source for ray sampling static RandomRaySourceShape source_shape_; // Flag for linear source - static bool mesh_subdivision_enabled_; // Flag for mesh subdivision static RandomRaySampleMethod sample_method_; // Flag for sampling method //---------------------------------------------------------------------------- diff --git a/include/openmc/random_ray/random_ray_simulation.h b/include/openmc/random_ray/random_ray_simulation.h index b94e7401b..3dec48bf2 100644 --- a/include/openmc/random_ray/random_ray_simulation.h +++ b/include/openmc/random_ray/random_ray_simulation.h @@ -21,11 +21,7 @@ public: // Methods void compute_segment_correction_factors(); void apply_fixed_sources_and_mesh_domains(); - void prepare_fixed_sources_adjoint(vector& forward_flux, - SourceRegionContainer& forward_source_regions, - SourceRegionContainer& forward_base_source_regions, - std::unordered_map& - forward_source_region_map); + void prepare_fixed_sources_adjoint(); void simulate(); void output_simulation_results() const; void instability_check( @@ -45,9 +41,6 @@ private: // Contains all flat source region data unique_ptr domain_; - // Random ray eigenvalue - double k_eff_ {1.0}; - // Tracks the average FSR miss rate for analysis and reporting double avg_miss_rate_ {0.0}; diff --git a/include/openmc/random_ray/source_region.h b/include/openmc/random_ray/source_region.h index 5c5b31f39..0f5a747ff 100644 --- a/include/openmc/random_ray/source_region.h +++ b/include/openmc/random_ray/source_region.h @@ -308,7 +308,6 @@ public: //---------------------------------------------------------------------------- // Constructors SourceRegion(int negroups, bool is_linear); - SourceRegion(const SourceRegionHandle& handle, int64_t parent_sr); SourceRegion() = default; //---------------------------------------------------------------------------- diff --git a/src/random_ray/flat_source_domain.cpp b/src/random_ray/flat_source_domain.cpp index 409238830..1bf27e1ed 100644 --- a/src/random_ray/flat_source_domain.cpp +++ b/src/random_ray/flat_source_domain.cpp @@ -53,24 +53,6 @@ FlatSourceDomain::FlatSourceDomain() : negroups_(data::mg.num_energy_groups_) // Initialize source regions. bool is_linear = RandomRay::source_shape_ != RandomRaySourceShape::FLAT; source_regions_ = SourceRegionContainer(negroups_, is_linear); - source_regions_.assign( - base_source_regions, SourceRegion(negroups_, is_linear)); - - // Initialize materials - int64_t source_region_id = 0; - for (int i = 0; i < model::cells.size(); i++) { - Cell& cell = *model::cells[i]; - if (cell.type_ == Fill::MATERIAL) { - for (int j = 0; j < cell.n_instances(); j++) { - source_regions_.material(source_region_id++) = cell.material(j); - } - } - } - - // Sanity check - if (source_region_id != base_source_regions) { - fatal_error("Unexpected number of source regions"); - } // Initialize tally volumes if (volume_normalized_flux_tallies_) { @@ -118,34 +100,24 @@ void FlatSourceDomain::accumulate_iteration_flux() } } -// Compute new estimate of scattering + fission sources in each source region -// based on the flux estimate from the previous iteration. -void FlatSourceDomain::update_neutron_source(double k_eff) +void FlatSourceDomain::update_single_neutron_source(SourceRegionHandle& srh) { - simulation::time_update_src.start(); - - double inverse_k_eff = 1.0 / k_eff; - -// Reset all source regions to zero (important for void regions) -#pragma omp parallel for - for (int64_t se = 0; se < n_source_elements(); se++) { - source_regions_.source(se) = 0.0; + // Reset all source regions to zero (important for void regions) + for (int g = 0; g < negroups_; g++) { + srh.source(g) = 0.0; } // Add scattering + fission source -#pragma omp parallel for - for (int64_t sr = 0; sr < n_source_regions(); sr++) { - int material = source_regions_.material(sr); - if (material == MATERIAL_VOID) { - continue; - } + int material = srh.material(); + if (material != MATERIAL_VOID) { + double inverse_k_eff = 1.0 / k_eff_; for (int g_out = 0; g_out < negroups_; g_out++) { double sigma_t = sigma_t_[material * negroups_ + g_out]; double scatter_source = 0.0; double fission_source = 0.0; for (int g_in = 0; g_in < negroups_; g_in++) { - double scalar_flux = source_regions_.scalar_flux_old(sr, g_in); + double scalar_flux = srh.scalar_flux_old(g_in); double sigma_s = sigma_s_[material * negroups_ * negroups_ + g_out * negroups_ + g_in]; double nu_sigma_f = nu_sigma_f_[material * negroups_ + g_in]; @@ -154,18 +126,30 @@ void FlatSourceDomain::update_neutron_source(double k_eff) scatter_source += sigma_s * scalar_flux; fission_source += nu_sigma_f * scalar_flux * chi; } - source_regions_.source(sr, g_out) = + srh.source(g_out) = (scatter_source + fission_source * inverse_k_eff) / sigma_t; } } // Add external source if in fixed source mode if (settings::run_mode == RunMode::FIXED_SOURCE) { -#pragma omp parallel for - for (int64_t se = 0; se < n_source_elements(); se++) { - source_regions_.source(se) += source_regions_.external_source(se); + for (int g = 0; g < negroups_; g++) { + srh.source(g) += srh.external_source(g); } } +} + +// Compute new estimate of scattering + fission sources in each source region +// based on the flux estimate from the previous iteration. +void FlatSourceDomain::update_all_neutron_sources() +{ + simulation::time_update_src.start(); + +#pragma omp parallel for + for (int64_t sr = 0; sr < n_source_regions(); sr++) { + SourceRegionHandle srh = source_regions_.get_source_region_handle(sr); + update_single_neutron_source(srh); + } simulation::time_update_src.stop(); } @@ -320,7 +304,7 @@ int64_t FlatSourceDomain::add_source_to_scalar_flux() // Generates new estimate of k_eff based on the differences between this // iteration's estimate of the scalar flux and the last iteration's estimate. -double FlatSourceDomain::compute_k_eff(double k_eff_old) const +void FlatSourceDomain::compute_k_eff() { double fission_rate_old = 0; double fission_rate_new = 0; @@ -365,7 +349,7 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const p[sr] = sr_fission_source_new; } - double k_eff_new = k_eff_old * (fission_rate_new / fission_rate_old); + double k_eff_new = k_eff_ * (fission_rate_new / fission_rate_old); double H = 0.0; // defining an inverse sum for better performance @@ -385,7 +369,7 @@ double FlatSourceDomain::compute_k_eff(double k_eff_old) const // Adds entropy value to shared entropy vector in openmc namespace. simulation::entropy.push_back(H); - return k_eff_new; + k_eff_ = k_eff_new; } // This function is responsible for generating a mapping between random @@ -652,7 +636,6 @@ void FlatSourceDomain::random_ray_tally() "random ray mode."); break; } - // Apply score to the appropriate tally bin Tally& tally {*model::tallies[task.tally_idx]}; #pragma omp atomic @@ -726,21 +709,21 @@ void FlatSourceDomain::output_to_vtk() const print_plot(); // Outer loop over plots - for (int p = 0; p < model::plots.size(); p++) { + for (int plt = 0; plt < model::plots.size(); plt++) { // Get handle to OpenMC plot object and extract params - Plot* openmc_plot = dynamic_cast(model::plots[p].get()); + Plot* openmc_plot = dynamic_cast(model::plots[plt].get()); // Random ray plots only support voxel plots if (!openmc_plot) { warning(fmt::format("Plot {} is invalid plot type -- only voxel plotting " "is allowed in random ray mode.", - p)); + plt)); continue; } else if (openmc_plot->type_ != Plot::PlotType::voxel) { warning(fmt::format("Plot {} is invalid plot type -- only voxel plotting " "is allowed in random ray mode.", - p)); + plt)); continue; } @@ -794,23 +777,11 @@ void FlatSourceDomain::output_to_vtk() const continue; } - int i_cell = p.lowest_coord().cell(); - int64_t sr = source_region_offsets_[i_cell] + p.cell_instance(); - if (RandomRay::mesh_subdivision_enabled_) { - int mesh_idx = base_source_regions_.mesh(sr); - int mesh_bin; - if (mesh_idx == C_NONE) { - mesh_bin = 0; - } else { - mesh_bin = model::meshes[mesh_idx]->get_bin(p.r()); - } - SourceRegionKey sr_key {sr, mesh_bin}; - auto it = source_region_map_.find(sr_key); - if (it != source_region_map_.end()) { - sr = it->second; - } else { - sr = -1; - } + SourceRegionKey sr_key = lookup_source_region_key(p); + int64_t sr = -1; + auto it = source_region_map_.find(sr_key); + if (it != source_region_map_.end()) { + sr = it->second; } voxel_indices[z * Ny * Nx + y * Nx + x] = sr; @@ -967,13 +938,17 @@ void FlatSourceDomain::output_to_vtk() const } void FlatSourceDomain::apply_external_source_to_source_region( - Discrete* discrete, double strength_factor, SourceRegionHandle& srh) + int src_idx, SourceRegionHandle& srh) { - srh.external_source_present() = 1; - + auto s = model::external_sources[src_idx].get(); + auto is = dynamic_cast(s); + auto discrete = dynamic_cast(is->energy()); + double strength_factor = is->strength(); const auto& discrete_energies = discrete->x(); const auto& discrete_probs = discrete->prob(); + srh.external_source_present() = 1; + for (int i = 0; i < discrete_energies.size(); i++) { int g = data::mg.get_group_index(discrete_energies[i]); srh.external_source(g) += discrete_probs[i] * strength_factor; @@ -981,8 +956,7 @@ void FlatSourceDomain::apply_external_source_to_source_region( } void FlatSourceDomain::apply_external_source_to_cell_instances(int32_t i_cell, - Discrete* discrete, double strength_factor, int target_material_id, - const vector& instances) + int src_idx, int target_material_id, const vector& instances) { Cell& cell = *model::cells[i_cell]; @@ -1000,16 +974,13 @@ void FlatSourceDomain::apply_external_source_to_cell_instances(int32_t i_cell, if (target_material_id == C_NONE || cell_material_id == target_material_id) { int64_t source_region = source_region_offsets_[i_cell] + j; - SourceRegionHandle srh = - source_regions_.get_source_region_handle(source_region); - apply_external_source_to_source_region(discrete, strength_factor, srh); + external_volumetric_source_map_[source_region].push_back(src_idx); } } } void FlatSourceDomain::apply_external_source_to_cell_and_children( - int32_t i_cell, Discrete* discrete, double strength_factor, - int32_t target_material_id) + int32_t i_cell, int src_idx, int32_t target_material_id) { Cell& cell = *model::cells[i_cell]; @@ -1017,14 +988,14 @@ void FlatSourceDomain::apply_external_source_to_cell_and_children( vector instances(cell.n_instances()); std::iota(instances.begin(), instances.end(), 0); apply_external_source_to_cell_instances( - i_cell, discrete, strength_factor, target_material_id, instances); + i_cell, src_idx, target_material_id, instances); } else if (target_material_id == C_NONE) { std::unordered_map> cell_instance_list = cell.get_contained_cells(0, nullptr); for (const auto& pair : cell_instance_list) { int32_t i_child_cell = pair.first; - apply_external_source_to_cell_instances(i_child_cell, discrete, - strength_factor, target_material_id, pair.second); + apply_external_source_to_cell_instances( + i_child_cell, src_idx, target_material_id, pair.second); } } } @@ -1070,36 +1041,17 @@ void FlatSourceDomain::convert_external_sources() "point source at {}", sp->r())); } - int i_cell = gs.lowest_coord().cell(); - int64_t sr = source_region_offsets_[i_cell] + gs.cell_instance(); + SourceRegionKey key = lookup_source_region_key(gs); - if (RandomRay::mesh_subdivision_enabled_) { - // If mesh subdivision is enabled, we need to determine which subdivided - // mesh bin the point source coordinate is in as well - int mesh_idx = source_regions_.mesh(sr); - int mesh_bin; - if (mesh_idx == C_NONE) { - mesh_bin = 0; - } else { - mesh_bin = model::meshes[mesh_idx]->get_bin(gs.r()); - } - // With the source region and mesh bin known, we can use the - // accompanying SourceRegionKey as a key into a map that stores the - // corresponding external source index for the point source. Notably, we - // do not actually apply the external source to any source regions here, - // as if mesh subdivision is enabled, they haven't actually been - // discovered & initilized yet. When discovered, they will read from the - // point_source_map to determine if there are any point source terms - // that should be applied. - SourceRegionKey key {sr, mesh_bin}; - point_source_map_[key] = es; - } else { - // If we are not using mesh subdivision, we can apply the external - // source directly to the source region as we do for volumetric domain - // constraint sources. - SourceRegionHandle srh = source_regions_.get_source_region_handle(sr); - apply_external_source_to_source_region(energy, strength_factor, srh); - } + // With the source region and mesh bin known, we can use the + // accompanying SourceRegionKey as a key into a map that stores the + // corresponding external source index for the point source. Notably, we + // do not actually apply the external source to any source regions here, + // as if mesh subdivision is enabled, they haven't actually been + // discovered & initilized yet. When discovered, they will read from the + // external_source_map to determine if there are any external source + // terms that should be applied. + external_point_source_map_[key].push_back(es); } else { // If not a point source, then use the volumetric domain constraints to @@ -1107,42 +1059,25 @@ void FlatSourceDomain::convert_external_sources() if (is->domain_type() == Source::DomainType::MATERIAL) { for (int32_t material_id : domain_ids) { for (int i_cell = 0; i_cell < model::cells.size(); i_cell++) { - apply_external_source_to_cell_and_children( - i_cell, energy, strength_factor, material_id); + apply_external_source_to_cell_and_children(i_cell, es, material_id); } } } else if (is->domain_type() == Source::DomainType::CELL) { for (int32_t cell_id : domain_ids) { int32_t i_cell = model::cell_map[cell_id]; - apply_external_source_to_cell_and_children( - i_cell, energy, strength_factor, C_NONE); + apply_external_source_to_cell_and_children(i_cell, es, C_NONE); } } else if (is->domain_type() == Source::DomainType::UNIVERSE) { for (int32_t universe_id : domain_ids) { int32_t i_universe = model::universe_map[universe_id]; Universe& universe = *model::universes[i_universe]; for (int32_t i_cell : universe.cells_) { - apply_external_source_to_cell_and_children( - i_cell, energy, strength_factor, C_NONE); + apply_external_source_to_cell_and_children(i_cell, es, C_NONE); } } } } } // End loop over external sources - -// Divide the fixed source term by sigma t (to save time when applying each -// iteration) -#pragma omp parallel for - for (int64_t sr = 0; sr < n_source_regions(); sr++) { - int material = source_regions_.material(sr); - if (material == MATERIAL_VOID) { - continue; - } - for (int g = 0; g < negroups_; g++) { - double sigma_t = sigma_t_[material * negroups_ + g]; - source_regions_.external_source(sr, g) /= sigma_t; - } - } } void FlatSourceDomain::flux_swap() @@ -1159,13 +1094,23 @@ void FlatSourceDomain::flatten_xs() const int a = 0; n_materials_ = data::mg.macro_xs_.size(); - for (auto& m : data::mg.macro_xs_) { + for (int i = 0; i < n_materials_; i++) { + auto& m = data::mg.macro_xs_[i]; for (int g_out = 0; g_out < negroups_; g_out++) { if (m.exists_in_model) { double sigma_t = m.get_xs(MgxsType::TOTAL, g_out, NULL, NULL, NULL, t, a); sigma_t_.push_back(sigma_t); + if (sigma_t < MINIMUM_MACRO_XS) { + Material* mat = model::materials[i].get(); + warning(fmt::format( + "Material \"{}\" (id: {}) has a group {} total cross section " + "({:.3e}) below the minimum threshold " + "({:.3e}). Material will be treated as pure void.", + mat->name(), mat->id(), g_out, sigma_t, MINIMUM_MACRO_XS)); + } + double nu_sigma_f = m.get_xs(MgxsType::NU_FISSION, g_out, NULL, NULL, NULL, t, a); nu_sigma_f_.push_back(nu_sigma_f); @@ -1206,7 +1151,7 @@ void FlatSourceDomain::flatten_xs() } } -void FlatSourceDomain::set_adjoint_sources(const vector& forward_flux) +void FlatSourceDomain::set_adjoint_sources() { // Set the adjoint external source to 1/forward_flux. If the forward flux is // negative, zero, or extremely close to zero, set the adjoint source to zero, @@ -1220,7 +1165,7 @@ void FlatSourceDomain::set_adjoint_sources(const vector& forward_flux) double max_flux = 0.0; #pragma omp parallel for reduction(max : max_flux) for (int64_t se = 0; se < n_source_elements(); se++) { - double flux = forward_flux[se]; + double flux = source_regions_.scalar_flux_final(se); if (flux > max_flux) { max_flux = flux; } @@ -1230,7 +1175,7 @@ void FlatSourceDomain::set_adjoint_sources(const vector& forward_flux) #pragma omp parallel for for (int64_t sr = 0; sr < n_source_regions(); sr++) { for (int g = 0; g < negroups_; g++) { - double flux = forward_flux[sr * negroups_ + g]; + double flux = source_regions_.scalar_flux_final(sr, g); if (flux <= ZERO_FLUX_CUTOFF * max_flux) { source_regions_.external_source(sr, g) = 0.0; } else { @@ -1239,6 +1184,7 @@ void FlatSourceDomain::set_adjoint_sources(const vector& forward_flux) if (flux > 0.0) { source_regions_.external_source_present(sr) = 1; } + source_regions_.scalar_flux_final(sr, g) = 0.0; } } @@ -1265,7 +1211,6 @@ void FlatSourceDomain::set_adjoint_sources(const vector& forward_flux) source_regions_.external_source_present(sr) = 0; } } - // Divide the fixed source term by sigma t (to save time when applying each // iteration) #pragma omp parallel for @@ -1326,13 +1271,14 @@ void FlatSourceDomain::apply_mesh_to_cell_instances(int32_t i_cell, if ((target_material_id == C_NONE && !is_target_void) || cell_material_id == target_material_id) { int64_t sr = source_region_offsets_[i_cell] + j; - if (source_regions_.mesh(sr) != C_NONE) { - // print out the source region that is broken: + // Check if the key is already present in the mesh_map_ + if (mesh_map_.find(sr) != mesh_map_.end()) { fatal_error(fmt::format("Source region {} already has mesh idx {} " "applied, but trying to apply mesh idx {}", - sr, source_regions_.mesh(sr), mesh_idx)); + sr, mesh_map_[sr], mesh_idx)); } - source_regions_.mesh(sr) = mesh_idx; + // If the SR has not already been assigned, then we can write to it + mesh_map_[sr] = mesh_idx; } } } @@ -1402,18 +1348,9 @@ void FlatSourceDomain::apply_meshes() } } -void FlatSourceDomain::prepare_base_source_regions() -{ - std::swap(source_regions_, base_source_regions_); - source_regions_.negroups() = base_source_regions_.negroups(); - source_regions_.is_linear() = base_source_regions_.is_linear(); -} - SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( - int64_t sr, int mesh_bin, Position r, double dist, Direction u) + SourceRegionKey sr_key, Position r, Direction u) { - SourceRegionKey sr_key {sr, mesh_bin}; - // Case 1: Check if the source region key is already present in the permanent // map. This is the most common condition, as any source region visited in a // previous power iteration will already be present in the permanent map. If @@ -1475,9 +1412,8 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( gs.r() = r + TINY_BIT * u; gs.u() = {1.0, 0.0, 0.0}; exhaustive_find_cell(gs); - int gs_i_cell = gs.lowest_coord().cell(); - int64_t sr_found = source_region_offsets_[gs_i_cell] + gs.cell_instance(); - if (sr_found != sr) { + int64_t sr_found = lookup_base_source_region_idx(gs); + if (sr_found != sr_key.base_source_region_id) { discovered_source_regions_.unlock(sr_key); SourceRegionHandle handle; handle.is_numerical_fp_artifact_ = true; @@ -1485,9 +1421,9 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( } // Sanity check on mesh bin - int mesh_idx = base_source_regions_.mesh(sr); + int mesh_idx = lookup_mesh_idx(sr_key.base_source_region_id); if (mesh_idx == C_NONE) { - if (mesh_bin != 0) { + if (sr_key.mesh_bin != 0) { discovered_source_regions_.unlock(sr_key); SourceRegionHandle handle; handle.is_numerical_fp_artifact_ = true; @@ -1496,7 +1432,7 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( } else { Mesh* mesh = model::meshes[mesh_idx].get(); int bin_found = mesh->get_bin(r + TINY_BIT * u); - if (bin_found != mesh_bin) { + if (bin_found != sr_key.mesh_bin) { discovered_source_regions_.unlock(sr_key); SourceRegionHandle handle; handle.is_numerical_fp_artifact_ = true; @@ -1508,26 +1444,60 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( // condition only occurs the first time the source region is discovered // (typically in the first power iteration). In this case, we need to handle // creation of the new source region and its storage into the parallel map. - // The new source region is created by copying the base source region, so as - // to inherit material, external source, and some flux properties etc. We - // also pass the base source region id to allow the new source region to - // know which base source region it is derived from. - SourceRegion* sr_ptr = discovered_source_regions_.emplace( - sr_key, {base_source_regions_.get_source_region_handle(sr), sr}); - discovered_source_regions_.unlock(sr_key); + // Additionally, we need to determine the source region's material, initialize + // the starting scalar flux guess, and apply any known external sources. + + // Call the basic constructor for the source region and store in the parallel + // map. + bool is_linear = RandomRay::source_shape_ != RandomRaySourceShape::FLAT; + SourceRegion* sr_ptr = + discovered_source_regions_.emplace(sr_key, {negroups_, is_linear}); SourceRegionHandle handle {*sr_ptr}; - // Check if the new source region contains a point source and apply it if so - auto it2 = point_source_map_.find(sr_key); - if (it2 != point_source_map_.end()) { - int es = it2->second; - auto s = model::external_sources[es].get(); - auto is = dynamic_cast(s); - auto energy = dynamic_cast(is->energy()); - double strength_factor = is->strength(); - apply_external_source_to_source_region(energy, strength_factor, handle); - int material = handle.material(); - if (material != MATERIAL_VOID) { + // Determine the material + int gs_i_cell = gs.lowest_coord().cell(); + Cell& cell = *model::cells[gs_i_cell]; + int material = cell.material(gs.cell_instance()); + + // If material total XS is extremely low, just set it to void to avoid + // problems with 1/Sigma_t + for (int g = 0; g < negroups_; g++) { + double sigma_t = sigma_t_[material * negroups_ + g]; + if (sigma_t < MINIMUM_MACRO_XS) { + material = MATERIAL_VOID; + break; + } + } + + handle.material() = material; + + // Store the mesh index (if any) assigned to this source region + handle.mesh() = mesh_idx; + + if (settings::run_mode == RunMode::FIXED_SOURCE) { + // Determine if there are any volumetric sources, and apply them. + // Volumetric sources are specifc only to the base SR idx. + auto it_vol = + external_volumetric_source_map_.find(sr_key.base_source_region_id); + if (it_vol != external_volumetric_source_map_.end()) { + const vector& vol_sources = it_vol->second; + for (int src_idx : vol_sources) { + apply_external_source_to_source_region(src_idx, handle); + } + } + + // Determine if there are any point sources, and apply them. + // Point sources are specific to the source region key. + auto it_point = external_point_source_map_.find(sr_key); + if (it_point != external_point_source_map_.end()) { + const vector& point_sources = it_point->second; + for (int src_idx : point_sources) { + apply_external_source_to_source_region(src_idx, handle); + } + } + + // Divide external source term by sigma_t + if (material != C_NONE) { for (int g = 0; g < negroups_; g++) { double sigma_t = sigma_t_[material * negroups_ + g]; handle.external_source(g) /= sigma_t; @@ -1535,6 +1505,21 @@ SourceRegionHandle FlatSourceDomain::get_subdivided_source_region_handle( } } + // Compute the combined source term + update_single_neutron_source(handle); + + // Unlock the parallel map. Note: we may be tempted to release + // this lock earlier, and then just use the source region's lock to protect + // the flux/source initialization stages above. However, the rest of the code + // only protects updates to the new flux and volume fields, and assumes that + // the source is constant for the duration of transport. Thus, using just the + // source region's lock by itself would result in other threads potentially + // reading from the source before it is computed, as they won't use the lock + // when only reading from the SR's source. It would be expensive to protect + // those operations, whereas generating the SR is only done once, so we just + // hold the map's bucket lock until the source region is fully initialized. + discovered_source_regions_.unlock(sr_key); + return handle; } @@ -1620,4 +1605,52 @@ void FlatSourceDomain::apply_transport_stabilization() } } +// Determines the base source region index (i.e., a material filled cell +// instance) that corresponds to a particular location in the geometry. Requires +// that the "gs" object passed in has already been initialized and has called +// find_cell etc. +int64_t FlatSourceDomain::lookup_base_source_region_idx( + const GeometryState& gs) const +{ + int i_cell = gs.lowest_coord().cell(); + int64_t sr = source_region_offsets_[i_cell] + gs.cell_instance(); + return sr; +} + +// Determines the index of the mesh (if any) that has been applied +// to a particular base source region index. +int FlatSourceDomain::lookup_mesh_idx(int64_t sr) const +{ + int mesh_idx = C_NONE; + auto mesh_it = mesh_map_.find(sr); + if (mesh_it != mesh_map_.end()) { + mesh_idx = mesh_it->second; + } + return mesh_idx; +} + +// Determines the source region key that corresponds to a particular location in +// the geometry. This takes into account both the base source region index as +// well as the mesh bin if a mesh is applied to this source region for +// subdivision. +SourceRegionKey FlatSourceDomain::lookup_source_region_key( + const GeometryState& gs) const +{ + int64_t sr = lookup_base_source_region_idx(gs); + int64_t mesh_bin = lookup_mesh_bin(sr, gs.r()); + return SourceRegionKey {sr, mesh_bin}; +} + +// Determines the mesh bin that corresponds to a particular base source region +// index and position. +int64_t FlatSourceDomain::lookup_mesh_bin(int64_t sr, Position r) const +{ + int mesh_idx = lookup_mesh_idx(sr); + int mesh_bin = 0; + if (mesh_idx != C_NONE) { + mesh_bin = model::meshes[mesh_idx]->get_bin(r); + } + return mesh_bin; +} + } // namespace openmc diff --git a/src/random_ray/linear_source_domain.cpp b/src/random_ray/linear_source_domain.cpp index 81412164e..e1ad68e3d 100644 --- a/src/random_ray/linear_source_domain.cpp +++ b/src/random_ray/linear_source_domain.cpp @@ -34,25 +34,18 @@ void LinearSourceDomain::batch_reset() } } -void LinearSourceDomain::update_neutron_source(double k_eff) +void LinearSourceDomain::update_single_neutron_source(SourceRegionHandle& srh) { - simulation::time_update_src.start(); - - double inverse_k_eff = 1.0 / k_eff; - -// Reset all source regions to zero (important for void regions) -#pragma omp parallel for - for (int64_t se = 0; se < n_source_elements(); se++) { - source_regions_.source(se) = 0.0; + // Reset all source regions to zero (important for void regions) + for (int g = 0; g < negroups_; g++) { + srh.source(g) = 0.0; } -#pragma omp parallel for - for (int64_t sr = 0; sr < n_source_regions(); sr++) { - int material = source_regions_.material(sr); - if (material == MATERIAL_VOID) { - continue; - } - MomentMatrix invM = source_regions_.mom_matrix(sr).inverse(); + // Add scattering + fission source + int material = srh.material(); + if (material != MATERIAL_VOID) { + double inverse_k_eff = 1.0 / k_eff_; + MomentMatrix invM = srh.mom_matrix().inverse(); for (int g_out = 0; g_out < negroups_; g_out++) { double sigma_t = sigma_t_[material * negroups_ + g_out]; @@ -64,8 +57,8 @@ void LinearSourceDomain::update_neutron_source(double k_eff) for (int g_in = 0; g_in < negroups_; g_in++) { // Handles for the flat and linear components of the flux - double flux_flat = source_regions_.scalar_flux_old(sr, g_in); - MomentArray flux_linear = source_regions_.flux_moments_old(sr, g_in); + double flux_flat = srh.scalar_flux_old(g_in); + MomentArray flux_linear = srh.flux_moments_old(g_in); // Handles for cross sections double sigma_s = @@ -81,7 +74,7 @@ void LinearSourceDomain::update_neutron_source(double k_eff) } // Compute the flat source term - source_regions_.source(sr, g_out) = + srh.source(g_out) = (scatter_flat + fission_flat * inverse_k_eff) / sigma_t; // Compute the linear source terms. In the first 10 iterations when the @@ -91,25 +84,21 @@ void LinearSourceDomain::update_neutron_source(double k_eff) // very small/noisy or have poorly developed spatial moments, so we zero // the source gradients (effectively making this a flat source region // temporarily), so as to improve stability. - if (simulation::current_batch > 10 && - source_regions_.source(sr, g_out) >= 0.0) { - source_regions_.source_gradients(sr, g_out) = + if (simulation::current_batch > 10 && srh.source(g_out) >= 0.0) { + srh.source_gradients(g_out) = invM * ((scatter_linear + fission_linear * inverse_k_eff) / sigma_t); } else { - source_regions_.source_gradients(sr, g_out) = {0.0, 0.0, 0.0}; + srh.source_gradients(g_out) = {0.0, 0.0, 0.0}; } } } + // Add external source if in fixed source mode if (settings::run_mode == RunMode::FIXED_SOURCE) { -// Add external source to flat source term if in fixed source mode -#pragma omp parallel for - for (int64_t se = 0; se < n_source_elements(); se++) { - source_regions_.source(se) += source_regions_.external_source(se); + for (int g = 0; g < negroups_; g++) { + srh.source(g) += srh.external_source(g); } } - - simulation::time_update_src.stop(); } void LinearSourceDomain::normalize_scalar_flux_and_volumes( diff --git a/src/random_ray/random_ray.cpp b/src/random_ray/random_ray.cpp index 27f674c42..89a91449d 100644 --- a/src/random_ray/random_ray.cpp +++ b/src/random_ray/random_ray.cpp @@ -237,7 +237,6 @@ double RandomRay::distance_inactive_; double RandomRay::distance_active_; unique_ptr RandomRay::ray_source_; RandomRaySourceShape RandomRay::source_shape_ {RandomRaySourceShape::FLAT}; -bool RandomRay::mesh_subdivision_enabled_ {false}; RandomRaySampleMethod RandomRay::sample_method_ {RandomRaySampleMethod::PRNG}; RandomRay::RandomRay() @@ -336,71 +335,60 @@ void RandomRay::event_advance_ray() void RandomRay::attenuate_flux(double distance, bool is_active, double offset) { - // Determine source region index etc. - int i_cell = lowest_coord().cell(); - - // The base source region is the spatial region index - int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance(); + // Lookup base source region index + int64_t sr = domain_->lookup_base_source_region_idx(*this); // Perform ray tracing across mesh - if (mesh_subdivision_enabled_) { - // Determine the mesh index for the base source region, if any - int mesh_idx = domain_->base_source_regions_.mesh(sr); + // Determine the mesh index for the base source region, if any + int mesh_idx = domain_->lookup_mesh_idx(sr); - if (mesh_idx == C_NONE) { - // If there's no mesh being applied to this cell, then - // we just attenuate the flux as normal, and set - // the mesh bin to 0 - attenuate_flux_inner(distance, is_active, sr, 0, r()); - } else { - // If there is a mesh being applied to this cell, then - // we loop over all the bin crossings and attenuate - // separately. - Mesh* mesh = model::meshes[mesh_idx].get(); - - // We adjust the start and end positions of the ray slightly - // to accomodate for floating point precision issues that tend - // to occur at mesh boundaries that overlap with geometry lattice - // boundaries. - Position start = r() + (offset + TINY_BIT) * u(); - Position end = start + (distance - 2.0 * TINY_BIT) * u(); - double reduced_distance = (end - start).norm(); - - // Ray trace through the mesh and record bins and lengths - mesh_bins_.resize(0); - mesh_fractional_lengths_.resize(0); - mesh->bins_crossed(start, end, u(), mesh_bins_, mesh_fractional_lengths_); - - // Loop over all mesh bins and attenuate flux - for (int b = 0; b < mesh_bins_.size(); b++) { - double physical_length = reduced_distance * mesh_fractional_lengths_[b]; - attenuate_flux_inner( - physical_length, is_active, sr, mesh_bins_[b], start); - start += physical_length * u(); - } - } + if (mesh_idx == C_NONE) { + // If there's no mesh being applied to this cell, then + // we just attenuate the flux as normal, and set + // the mesh bin to 0 + attenuate_flux_inner(distance, is_active, sr, 0, r()); } else { - attenuate_flux_inner(distance, is_active, sr, C_NONE, r()); + // If there is a mesh being applied to this cell, then + // we loop over all the bin crossings and attenuate + // separately. + Mesh* mesh = model::meshes[mesh_idx].get(); + + // We adjust the start and end positions of the ray slightly + // to accomodate for floating point precision issues that tend + // to occur at mesh boundaries that overlap with geometry lattice + // boundaries. + Position start = r() + (offset + TINY_BIT) * u(); + Position end = start + (distance - 2.0 * TINY_BIT) * u(); + double reduced_distance = (end - start).norm(); + + // Ray trace through the mesh and record bins and lengths + mesh_bins_.resize(0); + mesh_fractional_lengths_.resize(0); + mesh->bins_crossed(start, end, u(), mesh_bins_, mesh_fractional_lengths_); + + // Loop over all mesh bins and attenuate flux + for (int b = 0; b < mesh_bins_.size(); b++) { + double physical_length = reduced_distance * mesh_fractional_lengths_[b]; + attenuate_flux_inner( + physical_length, is_active, sr, mesh_bins_[b], start); + start += physical_length * u(); + } } } void RandomRay::attenuate_flux_inner( double distance, bool is_active, int64_t sr, int mesh_bin, Position r) { + SourceRegionKey sr_key {sr, mesh_bin}; SourceRegionHandle srh; - if (mesh_subdivision_enabled_) { - srh = domain_->get_subdivided_source_region_handle( - sr, mesh_bin, r, distance, u()); - if (srh.is_numerical_fp_artifact_) { - return; - } - } else { - srh = domain_->source_regions_.get_source_region_handle(sr); + srh = domain_->get_subdivided_source_region_handle(sr_key, r, u()); + if (srh.is_numerical_fp_artifact_) { + return; } switch (source_shape_) { case RandomRaySourceShape::FLAT: - if (this->material() == MATERIAL_VOID) { + if (srh.material() == MATERIAL_VOID) { attenuate_flux_flat_source_void(srh, distance, is_active, r); } else { attenuate_flux_flat_source(srh, distance, is_active, r); @@ -408,7 +396,7 @@ void RandomRay::attenuate_flux_inner( break; case RandomRaySourceShape::LINEAR: case RandomRaySourceShape::LINEAR_XY: - if (this->material() == MATERIAL_VOID) { + if (srh.material() == MATERIAL_VOID) { attenuate_flux_linear_source_void(srh, distance, is_active, r); } else { attenuate_flux_linear_source(srh, distance, is_active, r); @@ -439,7 +427,7 @@ void RandomRay::attenuate_flux_flat_source( n_event()++; // Get material - int material = this->material(); + int material = srh.material(); // MOC incoming flux attenuation + source contribution/attenuation equation for (int g = 0; g < negroups_; g++) { @@ -490,7 +478,7 @@ void RandomRay::attenuate_flux_flat_source_void( // The number of geometric intersections is counted for reporting purposes n_event()++; - int material = this->material(); + int material = srh.material(); // If ray is in the active phase (not in dead zone), make contributions to // source region bookkeeping @@ -537,7 +525,7 @@ void RandomRay::attenuate_flux_linear_source( // The number of geometric intersections is counted for reporting purposes n_event()++; - int material = this->material(); + int material = srh.material(); Position& centroid = srh.centroid(); Position midpoint = r + u() * (distance / 2.0); @@ -810,27 +798,12 @@ void RandomRay::initialize_ray(uint64_t ray_id, FlatSourceDomain* domain) cell_born() = lowest_coord().cell(); } + SourceRegionKey sr_key = domain_->lookup_source_region_key(*this); + SourceRegionHandle srh = + domain_->get_subdivided_source_region_handle(sr_key, r(), u()); + // Initialize ray's starting angular flux to starting location's isotropic // source - int i_cell = lowest_coord().cell(); - int64_t sr = domain_->source_region_offsets_[i_cell] + cell_instance(); - - SourceRegionHandle srh; - if (mesh_subdivision_enabled_) { - int mesh_idx = domain_->base_source_regions_.mesh(sr); - int mesh_bin; - if (mesh_idx == C_NONE) { - mesh_bin = 0; - } else { - Mesh* mesh = model::meshes[mesh_idx].get(); - mesh_bin = mesh->get_bin(r()); - } - srh = - domain_->get_subdivided_source_region_handle(sr, mesh_bin, r(), 0.0, u()); - } else { - srh = domain_->source_regions_.get_source_region_handle(sr); - } - if (!srh.is_numerical_fp_artifact_) { for (int g = 0; g < negroups_; g++) { angular_flux_[g] = srh.source(g); diff --git a/src/random_ray/random_ray_simulation.cpp b/src/random_ray/random_ray_simulation.cpp index 388a778b8..d475b2593 100644 --- a/src/random_ray/random_ray_simulation.cpp +++ b/src/random_ray/random_ray_simulation.cpp @@ -47,97 +47,82 @@ void openmc_run_random_ray() if (mpi::master) validate_random_ray_inputs(); - // Declare forward flux so that it can be saved for later adjoint simulation - vector forward_flux; - SourceRegionContainer forward_source_regions; - SourceRegionContainer forward_base_source_regions; - std::unordered_map - forward_source_region_map; + // Initialize Random Ray Simulation Object + RandomRaySimulation sim; - { - // Initialize Random Ray Simulation Object - RandomRaySimulation sim; + // Initialize fixed sources, if present + sim.apply_fixed_sources_and_mesh_domains(); - // Initialize fixed sources, if present - sim.apply_fixed_sources_and_mesh_domains(); + // Begin main simulation timer + simulation::time_total.start(); - // Begin main simulation timer - simulation::time_total.start(); + // Execute random ray simulation + sim.simulate(); - // Execute random ray simulation - sim.simulate(); + // End main simulation timer + simulation::time_total.stop(); - // End main simulation timer - simulation::time_total.stop(); - - // Normalize and save the final forward flux - sim.domain()->serialize_final_fluxes(forward_flux); - - double source_normalization_factor = - sim.domain()->compute_fixed_source_normalization_factor() / - (settings::n_batches - settings::n_inactive); + // Normalize and save the final forward flux + double source_normalization_factor = + sim.domain()->compute_fixed_source_normalization_factor() / + (settings::n_batches - settings::n_inactive); #pragma omp parallel for - for (uint64_t i = 0; i < forward_flux.size(); i++) { - forward_flux[i] *= source_normalization_factor; - } - - forward_source_regions = sim.domain()->source_regions_; - forward_source_region_map = sim.domain()->source_region_map_; - forward_base_source_regions = sim.domain()->base_source_regions_; - - // Finalize OpenMC - openmc_simulation_finalize(); - - // Output all simulation results - sim.output_simulation_results(); + for (uint64_t se = 0; se < sim.domain()->n_source_elements(); se++) { + sim.domain()->source_regions_.scalar_flux_final(se) *= + source_normalization_factor; } + // Finalize OpenMC + openmc_simulation_finalize(); + + // Output all simulation results + sim.output_simulation_results(); + ////////////////////////////////////////////////////////// // Run adjoint simulation (if enabled) ////////////////////////////////////////////////////////// - if (adjoint_needed) { - reset_timers(); - - // Configure the domain for adjoint simulation - FlatSourceDomain::adjoint_ = true; - - if (mpi::master) - header("ADJOINT FLUX SOLVE", 3); - - // Initialize OpenMC general data structures - openmc_simulation_init(); - - // Initialize Random Ray Simulation Object - RandomRaySimulation adjoint_sim; - - // Initialize adjoint fixed sources, if present - adjoint_sim.prepare_fixed_sources_adjoint(forward_flux, - forward_source_regions, forward_base_source_regions, - forward_source_region_map); - - // Transpose scattering matrix - adjoint_sim.domain()->transpose_scattering_matrix(); - - // Swap nu_sigma_f and chi - adjoint_sim.domain()->nu_sigma_f_.swap(adjoint_sim.domain()->chi_); - - // Begin main simulation timer - simulation::time_total.start(); - - // Execute random ray simulation - adjoint_sim.simulate(); - - // End main simulation timer - simulation::time_total.stop(); - - // Finalize OpenMC - openmc_simulation_finalize(); - - // Output all simulation results - adjoint_sim.output_simulation_results(); + if (!adjoint_needed) { + return; } + + reset_timers(); + + // Configure the domain for adjoint simulation + FlatSourceDomain::adjoint_ = true; + + if (mpi::master) + header("ADJOINT FLUX SOLVE", 3); + + // Initialize OpenMC general data structures + openmc_simulation_init(); + + sim.domain()->k_eff_ = 1.0; + + // Initialize adjoint fixed sources, if present + sim.prepare_fixed_sources_adjoint(); + + // Transpose scattering matrix + sim.domain()->transpose_scattering_matrix(); + + // Swap nu_sigma_f and chi + sim.domain()->nu_sigma_f_.swap(sim.domain()->chi_); + + // Begin main simulation timer + simulation::time_total.start(); + + // Execute random ray simulation + sim.simulate(); + + // End main simulation timer + simulation::time_total.stop(); + + // Finalize OpenMC + openmc_simulation_finalize(); + + // Output all simulation results + sim.output_simulation_results(); } // Enforces restrictions on inputs in random ray mode. While there are @@ -348,7 +333,6 @@ void validate_random_ray_inputs() // when generating weight windows with FW-CADIS and an overlaid mesh. /////////////////////////////////////////////////////////////////// if (RandomRay::source_shape_ == RandomRaySourceShape::LINEAR && - RandomRay::mesh_subdivision_enabled_ && variance_reduction::weight_windows.size() > 0) { warning( "Linear sources may result in negative fluxes in small source regions " @@ -366,7 +350,6 @@ void openmc_reset_random_ray() FlatSourceDomain::mesh_domain_map_.clear(); RandomRay::ray_source_.reset(); RandomRay::source_shape_ = RandomRaySourceShape::FLAT; - RandomRay::mesh_subdivision_enabled_ = false; RandomRay::sample_method_ = RandomRaySampleMethod::PRNG; } @@ -412,20 +395,11 @@ void RandomRaySimulation::apply_fixed_sources_and_mesh_domains() } } -void RandomRaySimulation::prepare_fixed_sources_adjoint( - vector& forward_flux, SourceRegionContainer& forward_source_regions, - SourceRegionContainer& forward_base_source_regions, - std::unordered_map& - forward_source_region_map) +void RandomRaySimulation::prepare_fixed_sources_adjoint() { + domain_->source_regions_.adjoint_reset(); if (settings::run_mode == RunMode::FIXED_SOURCE) { - if (RandomRay::mesh_subdivision_enabled_) { - domain_->source_regions_ = forward_source_regions; - domain_->source_region_map_ = forward_source_region_map; - domain_->base_source_regions_ = forward_base_source_regions; - domain_->source_regions_.adjoint_reset(); - } - domain_->set_adjoint_sources(forward_flux); + domain_->set_adjoint_sources(); } } @@ -445,22 +419,18 @@ void RandomRaySimulation::simulate() simulation::total_weight = 1.0; // Update source term (scattering + fission) - domain_->update_neutron_source(k_eff_); + domain_->update_all_neutron_sources(); - // Reset scalar fluxes, iteration volume tallies, and region hit flags to - // zero + // Reset scalar fluxes, iteration volume tallies, and region hit flags + // to zero domain_->batch_reset(); - // At the beginning of the simulation, if mesh subvivision is in use, we + // At the beginning of the simulation, if mesh subdivision is in use, we // need to swap the main source region container into the base container, // as the main source region container will be used to hold the true // subdivided source regions. The base container will therefore only // contain the external source region information, the mesh indices, // material properties, and initial guess values for the flux/source. - if (RandomRay::mesh_subdivision_enabled_ && - simulation::current_batch == 1 && !FlatSourceDomain::adjoint_) { - domain_->prepare_base_source_regions(); - } // Start timer for transport simulation::time_transport.start(); @@ -476,11 +446,9 @@ void RandomRaySimulation::simulate() simulation::time_transport.stop(); - // If using mesh subdivision, add any newly discovered source regions - // to the main source region container. - if (RandomRay::mesh_subdivision_enabled_) { - domain_->finalize_discovered_source_regions(); - } + // Add any newly discovered source regions to the main source region + // container. + domain_->finalize_discovered_source_regions(); // Normalize scalar flux and update volumes domain_->normalize_scalar_flux_and_volumes( @@ -494,10 +462,10 @@ void RandomRaySimulation::simulate() if (settings::run_mode == RunMode::EIGENVALUE) { // Compute random ray k-eff - k_eff_ = domain_->compute_k_eff(k_eff_); + domain_->compute_k_eff(); // Store random ray k-eff into OpenMC's native k-eff variable - global_tally_tracklength = k_eff_; + global_tally_tracklength = domain_->k_eff_; } // Execute all tallying tasks, if this is an active batch @@ -507,12 +475,6 @@ void RandomRaySimulation::simulate() // estimate domain_->accumulate_iteration_flux(); - // Generate mapping between source regions and tallies - if (!domain_->mapped_all_tallies_ && - !RandomRay::mesh_subdivision_enabled_) { - domain_->convert_source_regions_to_tallies(0); - } - // Use above mapping to contribute FSR flux data to appropriate // tallies domain_->random_ray_tally(); @@ -522,7 +484,7 @@ void RandomRaySimulation::simulate() domain_->flux_swap(); // Check for any obvious insabilities/nans/infs - instability_check(n_hits, k_eff_, avg_miss_rate_); + instability_check(n_hits, domain_->k_eff_, avg_miss_rate_); } // End MPI master work // Finalize the current batch @@ -571,7 +533,7 @@ void RandomRaySimulation::instability_check( } if (k_eff > 10.0 || k_eff < 0.01 || !(std::isfinite(k_eff))) { - fatal_error("Instability detected"); + fatal_error(fmt::format("Instability detected: k-eff = {:.5f}", k_eff)); } } } diff --git a/src/random_ray/source_region.cpp b/src/random_ray/source_region.cpp index 1205b995a..3b06f0ed0 100644 --- a/src/random_ray/source_region.cpp +++ b/src/random_ray/source_region.cpp @@ -48,7 +48,7 @@ SourceRegion::SourceRegion(int negroups, bool is_linear) } scalar_flux_new_.assign(negroups, 0.0); - source_.resize(negroups); + source_.assign(negroups, 0.0); scalar_flux_final_.assign(negroups, 0.0); tally_task_.resize(negroups); @@ -60,25 +60,6 @@ SourceRegion::SourceRegion(int negroups, bool is_linear) } } -SourceRegion::SourceRegion(const SourceRegionHandle& handle, int64_t parent_sr) - : SourceRegion(handle.negroups_, handle.is_linear_) -{ - material_ = handle.material(); - mesh_ = handle.mesh(); - parent_sr_ = parent_sr; - for (int g = 0; g < scalar_flux_new_.size(); g++) { - scalar_flux_old_[g] = handle.scalar_flux_old(g); - source_[g] = handle.source(g); - } - - if (settings::run_mode == RunMode::FIXED_SOURCE) { - external_source_present_ = handle.external_source_present(); - for (int g = 0; g < scalar_flux_new_.size(); g++) { - external_source_[g] = handle.external_source(g); - } - } -} - //============================================================================== // SourceRegionContainer implementation //============================================================================== @@ -259,9 +240,12 @@ void SourceRegionContainer::adjoint_reset() MomentMatrix {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}); std::fill(mom_matrix_t_.begin(), mom_matrix_t_.end(), MomentMatrix {0.0, 0.0, 0.0, 0.0, 0.0, 0.0}); - std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 0.0); + if (settings::run_mode == RunMode::FIXED_SOURCE) { + std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 0.0); + } else { + std::fill(scalar_flux_old_.begin(), scalar_flux_old_.end(), 1.0); + } std::fill(scalar_flux_new_.begin(), scalar_flux_new_.end(), 0.0); - std::fill(scalar_flux_final_.begin(), scalar_flux_final_.end(), 0.0); std::fill(source_.begin(), source_.end(), 0.0f); std::fill(external_source_.begin(), external_source_.end(), 0.0f); std::fill(source_gradients_.begin(), source_gradients_.end(), diff --git a/src/settings.cpp b/src/settings.cpp index 03d42bb9d..f9a2469db 100644 --- a/src/settings.cpp +++ b/src/settings.cpp @@ -346,7 +346,6 @@ void get_run_parameters(pugi::xml_node node_base) } FlatSourceDomain::mesh_domain_map_[mesh_id].emplace_back( type, domain_id); - RandomRay::mesh_subdivision_enabled_ = true; } } } diff --git a/tests/regression_tests/random_ray_low_density/__init__.py b/tests/regression_tests/random_ray_low_density/__init__.py new file mode 100644 index 000000000..e69de29bb diff --git a/tests/regression_tests/random_ray_low_density/inputs_true.dat b/tests/regression_tests/random_ray_low_density/inputs_true.dat new file mode 100644 index 000000000..c4fd06f42 --- /dev/null +++ b/tests/regression_tests/random_ray_low_density/inputs_true.dat @@ -0,0 +1,244 @@ + + + + mgxs.h5 + + + + + + + + + + + + + + + + + + + + + 2.5 2.5 2.5 + 12 12 12 + 0.0 0.0 0.0 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 +1 1 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 +2 2 2 2 2 2 2 2 2 2 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 + +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 +3 3 3 3 3 3 3 3 3 3 3 3 + + + + + + + + + + fixed source + 90 + 10 + 5 + + + 100.0 1.0 + + + universe + 1 + + + multi-group + + 500.0 + 100.0 + + + 0.0 0.0 0.0 30.0 30.0 30.0 + + + True + + + + + 1 + + + 2 + + + 3 + + + 3 + flux + tracklength + + + 2 + flux + tracklength + + + 1 + flux + tracklength + + + diff --git a/tests/regression_tests/random_ray_low_density/results_true.dat b/tests/regression_tests/random_ray_low_density/results_true.dat new file mode 100644 index 000000000..a4b3ee1bc --- /dev/null +++ b/tests/regression_tests/random_ray_low_density/results_true.dat @@ -0,0 +1,9 @@ +tally 1: +5.973607E-01 +7.155477E-02 +tally 2: +3.206216E-02 +2.063375E-04 +tally 3: +2.096415E-03 +8.804963E-07 diff --git a/tests/regression_tests/random_ray_low_density/test.py b/tests/regression_tests/random_ray_low_density/test.py new file mode 100644 index 000000000..1b4ffb781 --- /dev/null +++ b/tests/regression_tests/random_ray_low_density/test.py @@ -0,0 +1,60 @@ +import os + +import numpy as np +import openmc +from openmc.examples import random_ray_three_region_cube + +from tests.testing_harness import TolerantPyAPITestHarness + + +class MGXSTestHarness(TolerantPyAPITestHarness): + def _cleanup(self): + super()._cleanup() + f = 'mgxs.h5' + if os.path.exists(f): + os.remove(f) + + +def test_random_ray_low_density(): + model = random_ray_three_region_cube() + + # Rebuild the MGXS library to have a material with very + # low macroscopic cross sections + ebins = [1e-5, 20.0e6] + groups = openmc.mgxs.EnergyGroups(group_edges=ebins) + + void_sigma_a = 4.0e-6 + void_sigma_s = 3.0e-4 + void_mat_data = openmc.XSdata('void', groups) + void_mat_data.order = 0 + void_mat_data.set_total([void_sigma_a + void_sigma_s]) + void_mat_data.set_absorption([void_sigma_a]) + void_mat_data.set_scatter_matrix( + np.rollaxis(np.array([[[void_sigma_s]]]), 0, 3)) + + absorber_sigma_a = 0.75 + absorber_sigma_s = 0.25 + absorber_mat_data = openmc.XSdata('absorber', groups) + absorber_mat_data.order = 0 + absorber_mat_data.set_total([absorber_sigma_a + absorber_sigma_s]) + absorber_mat_data.set_absorption([absorber_sigma_a]) + absorber_mat_data.set_scatter_matrix( + np.rollaxis(np.array([[[absorber_sigma_s]]]), 0, 3)) + + multiplier = 0.0000001 + source_sigma_a = void_sigma_a * multiplier + source_sigma_s = void_sigma_s * multiplier + source_mat_data = openmc.XSdata('source', groups) + source_mat_data.order = 0 + source_mat_data.set_total([source_sigma_a + source_sigma_s]) + source_mat_data.set_absorption([source_sigma_a]) + source_mat_data.set_scatter_matrix( + np.rollaxis(np.array([[[source_sigma_s]]]), 0, 3)) + + mg_cross_sections_file = openmc.MGXSLibrary(groups) + mg_cross_sections_file.add_xsdatas( + [source_mat_data, void_mat_data, absorber_mat_data]) + mg_cross_sections_file.export_to_hdf5() + + harness = MGXSTestHarness('statepoint.10.h5', model) + harness.main() diff --git a/tests/regression_tests/random_ray_point_source_locator/results_true.dat b/tests/regression_tests/random_ray_point_source_locator/results_true.dat index 1785dda57..8c6f358dd 100644 --- a/tests/regression_tests/random_ray_point_source_locator/results_true.dat +++ b/tests/regression_tests/random_ray_point_source_locator/results_true.dat @@ -1,9 +1,9 @@ tally 1: -2.633900E+00 -2.948207E+00 +2.633923E+00 +2.948228E+00 tally 2: -1.440463E-01 -3.294032E-03 +1.440456E-01 +3.293984E-03 tally 3: 9.425207E-03 1.089748E-05 From 1dacf4fd2bba3582b575b6abf92b5d201eabab37 Mon Sep 17 00:00:00 2001 From: Patrick Shriwise Date: Thu, 2 Oct 2025 11:48:02 -0500 Subject: [PATCH 6/6] Add missing documentation on in depletion chain file format (#3590) Co-authored-by: Paul Romano Co-authored-by: April Novak --- docs/source/io_formats/depletion_chain.rst | 21 +++++++++++++++++++++ openmc/data/decay.py | 2 +- 2 files changed, 22 insertions(+), 1 deletion(-) diff --git a/docs/source/io_formats/depletion_chain.rst b/docs/source/io_formats/depletion_chain.rst index 89c76525f..74413e7b6 100644 --- a/docs/source/io_formats/depletion_chain.rst +++ b/docs/source/io_formats/depletion_chain.rst @@ -56,6 +56,27 @@ attributes: .. _io_chain_reaction: +-------------------- +```` Element +-------------------- + +The ```` element represents photon and electron sources associated with +the decay of a nuclide and contains information to construct an +:class:`openmc.stats.Univariate` object that represents this emission as an +energy distribution. This element has the following attributes: + + :type: + The type of :class:`openmc.stats.Univariate` source term. + + :particle: + The type of particle emitted, e.g., 'photon' or 'electron' + + :parameters: + The parameters of the source term, e.g., for a + :class:`openmc.stats.Discrete` source, the energies (in [eV]) at which the + particles are emitted and their relative intensities in [Bq/atom] (in other + words, decay constants). + ---------------------- ```` Element ---------------------- diff --git a/openmc/data/decay.py b/openmc/data/decay.py index 1a11d3614..c8a0bb5e7 100644 --- a/openmc/data/decay.py +++ b/openmc/data/decay.py @@ -591,7 +591,7 @@ def decay_photon_energy(nuclide: str) -> Univariate | None: openmc.stats.Univariate or None Distribution of energies in [eV] of photons emitted from decay, or None if no photon source exists. Note that the probabilities represent - intensities, given as [Bq]. + intensities, given as [Bq/atom] (in other words, decay constants). """ if not _DECAY_PHOTON_ENERGY: chain_file = openmc.config.get('chain_file')