Issue 
EPJ Nuclear Sci. Technol.
Volume 9, 2023
Euratom Research and Training in 2022: the Awards collection



Article Number  11  
Number of page(s)  9  
Section  Part 1: Safety research and training of reactor systems  
DOI  https://doi.org/10.1051/epjn/2023002  
Published online  10 February 2023 
https://doi.org/10.1051/epjn/2023002
Regular Article
Coupling the GUARDYAN code to subchanflow
Budapest University of Technology and Economics, Institute of Nuclear Techniques, Department of Nuclear Techniques, Műegyetem rkp. 3–9, H1111 Budapest, Hungary
^{*} email: legrady@reak.bme.hu
Received:
2
September
2022
Received in final form:
19
December
2022
Accepted:
12
January
2023
Published online: 10 February 2023
GUARDYAN is a GPUbased dynamic Monte Carlo code developed at the Budapest University of Technology and Economics, Hungary. Dynamic Monte Carlo computes the neutron population evolution by calculating the direct time dependence of the neutron histories in multiplying systems. Some wellestablished Monte Carlo codes have DMC versions with coupling to ThermalHydraulic solvers. GUARDYAN has the computational advantage of applying GPUs, thus calculation burden can be carried by commonly available hardware, and is capable of handling power plant size systems for kinetics problems. GUARDYAN has been recently coupled to the subchannel thermalhydraulics code SUBCHANFLOW in order to carry out dynamic calculations with TH feedback. This paper describes some convergence studies regarding reaching the initial equilibrium state. A literaturesuggested set of TH input settings and high sample numbers resulted in very low statistical errors of the power estimates and stable global measures (L_{2}) of power release, fuel, and coolant temperatures for both static and dynamic convergence. Dynamic mode lowsample simulations provided surprisingly stable global L_{2} measures, correct fuel temperatures, and power release, while coolant temperatures were off, without any indication of the incorrectness of the result. Static convergence showed an alternating, fluctuating L_{2} behavior that did not affect the final stable state.
© D. Legrady et al., Published by EDP Sciences, 2023
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1. Introduction
GUARDYAN (GpU Assisted Reactor Dynamic ANalysis) is a GPUbased dynamic Monte Carlo code developed at the Budapest University of Technology and Economics, Hungary [1–3]. As an emerging calculation technology, Dynamic Monte Carlo (DMC) [4] computes the neutron population evolution by calculating the direct time dependence of the neutron histories in multiplying systems. Some wellestablished Monte Carlo codes have DMC versions [5, 6] with coupling to ThermalHydraulic (TH) solvers. GUARDYAN has the computational advantage of applying GPUs, thus calculation burden can be carried by commonly available hardware and is capable of handling power plant size systems for kinetics problems. GUARDYAN has been recently coupled to the subchannel thermalhydraulics code SUBCHANFLOW in order to carry out dynamic calculations with TH feedback. This paper is the first to report on this coupling, further, we describe some convergence studies regarding reaching the initial equilibrium state. Comparison is given of reaching the steady state using dynamic or static (kmode) calculations.
2. The main features of GUARDYAN
GUARDYAN is specifically written to perform dynamic Monte Carlo calculations using GPUs. Geometry setup uses constructive solid geometry, where cells with homogenous material distributions are bounded by secondorder surfaces. Within a cell, temperatures may continuously change thanks to the Woodcock (delta) tracking algorithm for free path sampling. The code is written for the Linux environment in C++ and CUDA. Code verification involved rigorous comparison to wellestablished Monte Carlo and deterministic codes. Dynamic validation involved the reproduction of rapid reactivity insertions and withdrawals at the BME Training Reactor, while the current highlight being the reproduction of excore detector signals of a roddrop experiment of a VVER440/V213 nuclear power plant of unit 4 at the Paks NPP at the beginning of campaign number 32 with realistic burnup data considered [2], see Figure 1
Fig. 1. Power distribution in the core of Paks NPP Unit 4 campaign 32rod drop experiment beginning state (left), halfway (middle), and end state (right). 
The dynamic versions of the CPUbased codes TRIPOLI4 [5] and Serpent2 [6] are also coupled to SUBCHANFLOW and transient simulations for realistic small reactors with comparison to measurements or benchmark scenarios for larger, NPPsized cores have been reported. Performance comparison is challenging due to different computing architectures.
3. Coupling to thermalhydraulics
GUARDYAN has been coupled to the thermalhydraulic subchannel code SUBCHANFLOW (SCF) [7]. Based on the COBRATF set of codes [8, 9] and COBRAEN, the code SCF was developed at the Karlsruhe Institute of Technology. The development included a modernized software implementation in Fortran 95, substantial improvement in userfriendliness, and more robust error handling but foremost a refined twophase flow handling. Empirical correlations of pressure drop, void generation, heat transfer, etc. are incorporated into the basic conservation equations of energy, momentum, and mass and implemented in SCF making it capable of accurate modeling of light water reactors. Water properties are updated to comply with IAPWS97 (International Association for the Properties of Water and Steam). Other liquid and gas coolants can also be specified such as liquid sodium or lead, or gases like Helium or air. Hexagonal and rectangular lattices can be specified for regular plant geometries.
In the course of the McSafe [10] project SCF was coupled to the DMC version of TRIPOLI and Serpent2. The SCF features developed to ease coupling with these external codes provided an almost readily applicable way for coupling to GUARDYAN. For a flexible solution only subroutinelevel modifications (like ISO_C_BINDING) were used and we did not utilize the SALOME^{1} platform. The two codes are compiled together, thus data is directly exchanged between SCF and GUARDYAN without the need for generating data files for data swapping.
3.1. Coupling scheme
The thermalhydraulics solver expects to receive input data on the power release per pin axial node resulting from the Monte Carlo simulation. GUARDYAN considers interactions for fission power release other sources of energy releases were not considered. The probability of fission interaction probability was favorably biased towards more frequent sampling with proper weight compensation to keep the expected value intact.
The coupling of thermalhydraulics and neutronics follows the Operator Splitting SemiImplicit (OSSI) method. Essentially, fuel temperature (T_{f}), coolant temperature (T_{c}), and coolant density (ρ) data of the previous time step (i − 1) are used in the MC simulation for the calculation (F_{N}) of the pinwise power release:
$$\begin{array}{c}\hfill {P}^{k}={F}_{N}\left({T}_{f}^{i1},{T}_{c}^{i1},{\rho}^{i1}\right).\end{array}$$(1)
The thermalhydraulics module returns (F _{TH}) the new temperature and density data based on the power release:
$$\begin{array}{c}\hfill \left[{T}_{f}^{i},{T}_{c}^{i},{\rho}^{i}\right]={F}_{\mathrm{TH}}\left({P}^{i}\right).\end{array}$$(2)
The relaxation parameter was not introduced to average previous time step results to obtain the new one, though this may have improved the convergence rate, especially with the GK mode. The thermalhydraulic calculations are executed separately for each subchannel. Fuel and coolant temperatures are used for calculating temperaturedependent crosssections.
GUARDYAN coupled dynamic calculations are executed in three distinct phases, see Figure 2, GD stands for GUARDYAN Dynamic mode. In the first, sourceconvergence phase the equilibrium dynamic flux shape is calculated. Starting from uniform distribution dynamic transport is pursued with instant decay of precursors for a millisecond (in thermal systems). At the end of the time step the nodewise power distribution is normalized to the nominal target value and the thermalhydraulics module is called. A fully converging static TH calculation is carried out until the TH equilibrium is reached. After these temperature and density values are returned, the loop of these two steps repeats until convergence is reached. Criticality is set by the negative feedback of the system controlled by the moderator Boron concentration.
Fig. 2. Coupling scheme of GUARDYAN Dynamic (GD) mode to thermalhydraulics (TH) solvers. 
In the second step, the delayed neutron precursors are generated in a pure neutronics run, during which neutron histories are followed and at fission sites precursor samples are banked until the designated memory space is filled. Temperatures and densities are kept constant, nodewise power release is not estimated and the thermalhydraulic module is not called.
In the third step, the transient calculation starts with precursor decaying with properly sampled delayed neutron time dependence, and at each coupling time step the TH module is called and temperature and density parameters are exchanged. Nodewise power release is estimated but it is not normalized anymore to the nominal power. The TH module now only calculates the length of the coupling time step in contrast to the first step where static TH results are obtained.
General purpose MC codes usually operate in keigenvalue mode. For aiding comparison with these codes GUARDYAN can also calculate static k _{eff}. In this GUARDYAN kmode (GK) calculation neutrons are followed from their birth at a fission site until the next induced fission or a termination (escape, failed Russian Roulette) event. For the coupling scheme, see Figure 3 inactive cycles are run first (GK I) to reach neutron source convergence, now the power release is not estimated. Thereafter, during the active cycles (GK A) nodewise power release is estimated. After normalization to nominal power, data is given to the TH solver. With the new temperature and density data the next inactive cycle starts. This loop repeats until stopped, the last step offers the estimation for k_{eff}.
Fig. 3. Coupling scheme of static calculation (GUARDYAN, GK) and thermalhydraulics (TH). 
4. Convergence studies on the TMI Minicore benchmark
4.1. The TMI Minicore model
Based on the Three Mile Island TMI1 reactor a 3 × 3assembly simplified core model is presented in the NEA benchmark set [11] intended for benchmarking kinetic calculations. The 9 fuel assemblies contain 15 × 15 fuel pins made of 4.12% enriched UOX fuel. In each assembly, four Gdcontaining burnable poison pins are located. The fuel length is 353 cm, and the assembly is 21.64 cm wide. At the center of the assemblies, void tubes can be found. The core is surrounded by a 21.64 cm thick layer of a mixture of steel and borated water. The central assembly has the specialty of housing 16 control rods made of Ag–In–Cd. At Hot Zero Power control rods are fully inserted, see Figure 4.
Fig. 4. TMI Minicore crosssection material models. 
The SCF model of the TMI Minicore places a subchannel between each set of four fuel pins, filled with water. The inlet temperature is kept constant at 565 K, outlet pressure is 15.5 MPa, and mass flow rate is set to 773.64 kg/s. The power is defined as 140.94 MW. The core axial length was divided into 30 axial nodes, within a pin further radial divisions were not created.
4.2. Convergence to steady state in dynamic mode
To scout, the convergence properties of the dynamic coupling scheme, calculations for the TMI Minicore benchmark were carried out with two different parameter settings. Run #1 used more samples (2^{20} neutrons, about 1 000 000) but with less frequent coupling with an interval size of 5 ms. Run #1 was based on the settings found in référence [12]. Run #2 employed a significantly smaller sample size (2^{17} neutrons, equaling about 100 000) with frequent coupling stepping carried out at each 100 μs. With Run #2 the intention was to create a statistically unstable, noisy solution that is coupled very frequently to the TH solver inflating the chance that the coupled calculation becomes statistically unstable and observably fails.
For assessing the convergence state a measure has been defined as the L _{2} distance of the distributions in two subsequent time steps:
$$\begin{array}{c}\hfill {L}_{2}\left({X}^{l}{X}^{l1}\right)=\frac{\sqrt{{\sum}_{i}{\left[{X}_{i}^{l}{X}_{i}^{l1}\right]}^{2}}}{\sqrt{N}},\end{array}$$(3)
where l is the time step number and X _{ i } is a certain quantity in the space segment i, with i having the maximum of N. X _{ i } may stand for nodewise power release, temperature, or density. When L _{2} becomes stationary as a function of l, the convergence to a steady state may have been reached. As a result of a stochastic calculation, a fluctuation around the stationary value appears, reaching a steady state is meant within the stochastic error.
The equilibrium state of the neutron population in GD mode can be attributed to an equilibrium state of the neutron flux resulting from an initial source convergence process similar to power iteration for equilibrium fission source distribution. This would reflect also in a steady power evolution as a function of time steps.
Figure 5 shows the initial source convergence. For Run #2, power stabilizes after a few coupling time steps without a visible longterm trend. Fuel temperature and similarly coolant temperature achieve a stable value after 20–30 time steps. Larger sample numbers (and consequently lower relative variances in the power estimates) in Run #1 allow for faster convergence to the initial ground state in terms of the number of time steps, see also Figure 6, showing the data of Figure 5 in an L_{2} sense. Average coolant temperature and power converge to the same value, however, converged fuel temperatures differ. This indicates that stochastic perturbation of the released power has a nonlinear effect on the equilibrium state, triggering an ongoing detailed analysis of THMC coupling, see référence [13]. Lower sample numbers of Run #2 result in stablelooking estimates, devoid of large fluctuations of the averages, with little indication that the convergence may be false.
Fig. 5. Source convergence in dynamic mode, power is normalized to a single starting neutron. 
Fig. 6. L_{2} norm evolution of the GD source convergence for Run #1 (red) and Run #2 (blue). 
A map of the pinwise fuel temperature distribution shows a visible difference in the estimation stability of higher and lower sample numbers, see Figure 7.
Fig. 7. Pinwise fuel temperature distribution at the converged state, for Run #2 (left) and Run #1 (right) accumulated during the corresponding time step. 
The apparently noisy nodewise temperature distribution in Figure 7 hardly suggests that a converged equilibrium state may be reached, or a stable longerterm calculation is feasible. Low sample numbers and short time steps warrant very noisy nodewise power release see Figure 8 that could average in longer time spans, and through that, may provide relaxed, smoother input for the TH calculation with larger coupling time steps. That not being the case, surprisingly, the TH solver tames the noisy input into a globally smooth L_{2} evolution as seen in Figure 9
Fig. 8. Nodewise power estimates for Run #2 (left) and Run #1 (right) accumulated during the corresponding time step. 
Fig. 9. Dynamic simulation average power, fuel, and coolant temperatures, for Run #1 (red) and Run #2 (blue). 
Looking at the time evolution of the integral quantities (see Fig. 9) of the timedependent system after the source convergence and delayed neutron sampling have already been done, power seems to be stable within the statistical fluctuations: smoother for Run #1 (red curve), more undulating for Run #2.
Temperature evolution is very smooth. Coolant temperature seems to be closely matching, but fuel temperatures systematically differ by about 20 degrees. It seems that the statistical poorness of the sampling does not directly translate into the statistical illbeing of the TH integral quantities but may result in an unconverged, false solution.
4.3. Convergence to steady state in kmode
An equilibrium critical state may also be reached by GK calculation. Starting sample number was chosen to be 2^{20}, at each coupling step the neutronics source convergence was run for 100 inactive and 400 active cycles. During the inactive cycles, the estimates of the power release were disregarded but the updated temperature values took effect. Therefore, the new temperature and neutron fields did not go in entwined progress like in the GD mode. Figure 10 shows the alternating convergence to an equilibrium end state. Different random number generator seeds were chosen to show that the large stepwise fluctuations are not caused by randomness and that the converged state is common. The inactive cycles disconnect the previous neutronic state with the next coupling state, hence the alternating deviation from the end state appears. The converged temperature values agree with those found in the GD mode.
Fig. 10. Convergence in static GK mode, evolution of the criticality (k _{eff}), fuel, and coolant temperatures with the coupling steps. 
Nodewise relative error figures per GK step stay below 2% in the middle axial section of the reactor, see Figure 11. Even for the topmost nodes, where a high leakage rate decreases estimation efficiency the nodewise relative standard deviation stays below 10%.
Fig. 11. Relative errors of the nodewise power estimates in percent for the lowest (left), central (middle), and topmost (right) nodes. 
Relevance of these error figures may count if some local anomaly is sought, but stakes are higher if fitness for TH coupling is to be judged. In this respect, it may be more telling to investigate how the error in the power estimation translates into the error in temperatures. A local statistical power anomaly in a millisecondscale time interval is averaged out by the much slower thermalhydraulic processes. Moreover, the directed stream of the coolant flow may transport this statistical error to a different position in the core.
In Figure 12 the fuel temperature standard deviation is plotted against the standard deviation of the power estimates for each pin node. Correlation is not apparent, but fuel temperature errors are generally smaller than the causing statistical error in the power. As power estimation errors increase the amount of statistical inaccuracy transferred seems to decrease.
Fig. 12. Correlation of standard deviation of power and fuel temperature. 
4.4. Comparison of kcode and dynamic convergence
Apart from the globally characteristic number as average fuel and coolant temperatures or average flux, nodewise data were compared, see Figure 13. Dynamic and static neutron fluxes should coincide for a critical reactor, and converged THstates should not differ either.
Fig. 13. Comparison of dynamic and static convergence results: relative difference of power output (left), relative difference of coolant temperatures (middle), and relative difference of fuel temperatures (right). 
Subpercent pinwise differences in all three regarded comparisons indicate that the converged endstates are the same. Apparently random deviations can be observed in the power and the fuel temperature, in the range of a few percent that corresponds to errors in the GK mode. The coolant temperature errors seem to show some weak patterns but the amplitude of the deviations stays way below a percent.
5. Conclusions
Coupling of the direct Monte Carlo code GUARDYAN to the subchannel code SUBCHANFLOW has been achieved and convergence to a steady state was investigated. A literaturesuggested set of TH input settings and high sample numbers resulted in very low statistical errors of the power estimates and stable global measures (L_{2}) of power release, fuel, and coolant temperatures for both static and dynamic convergence. Dynamic mode lowsample simulations provided surprisingly stable global L_{2} measures, correct fuel temperatures, and power release, while coolant temperatures were off, without any indication of the incorrectness of the result. Static convergence showed an alternating, fluctuating L_{2} behavior that did not affect the final stable state.
Further research should target the detailed analysis of the stability, robustness, and unbiasedness of coupled dynamic Monte Carlo and thermal hydraulics.
Conflict of interests
The authors declare that they have no competing interests to report.
Acknowledgments
The authors acknowledge the support and computational resources of the Wigner Scientific Computational Laboratory (WSCLAB) (the former Wigner GPU Laboratory).
Funding
The research reported in this paper is part of project no. BMEEGA02, implemented with the support provided by the Ministry of Innovation and Technology of Hungary from the National Research, Development and Innovation Fund, financed under the TKP2021 funding scheme.
Data availability statement
Data associated with this article cannot be disclosed due to legal reasons.
Author contribution statement
David Legrady: manuscript preparation, research group leader, conceptualization.
Marton Pukler: algorithm testing, modeling.
Balint Panka: reactor geometry modeling.
Gabor Tolnai: simulations, modeling, programming, data analysis, graphs.
References
 B. Molnar, G. Tolnai, D. Legrady, A GPU based direct Monte Carlo simulation of time dependence in nuclear reactors, Ann. Nucl. Energy 132, 46–63 (2019) [CrossRef] [Google Scholar]
 D. Legrady, G. Tolnai, T. Hajas, E. Pazman, T. Parko, I. Pos, Full Core PinLevel VVER440 Simulation of a Rod Drop Experiment with the GPUBased Monte Carlo Code GUARDYAN, Energies 15, 2712 (2022) [CrossRef] [Google Scholar]
 Official webpage of GUARDYAN, http://awing.reak.bme.hu/GUARDYAN [Google Scholar]
 B.L. Sjenitzer, J.E. Hoogenboom, Dynamic Monte Carlo method for nuclear reactor kinetics calculations, Nucl. Sci. Eng. 175, 94–107 (2013) [CrossRef] [Google Scholar]
 M. Faucher, D. Mancusi, A. Zoia, New kinetic simulation capabilities for Tripoli4: Methods and applications, Ann. Nucl. Energy 120, 74–88 (2018) [CrossRef] [Google Scholar]
 J. Leppanen, Development of a dynamic simulation mode in Serpent 2 Monte Carlo code, in Proceedings of the International 426, Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering (M&C 2013), LaGrange Park, Illinois (2013), p. 427 [Google Scholar]
 U. Imke, V.H. Sanchez, Validation of the subchannel code SUBCHANFLOW using the NUPEC PWR Tests (PSBT), Sci. Technol. Nucl. Instal. 12, 465059 (2012) [Google Scholar]
 M.J. Thurgood et al., COBRA/TRAC – A thermalhydraulic code for transient analysis of nuclear reactor vessel and primary coolant systems, in NUREG/CR3046, 1983 [Google Scholar]
 D. Basile, R. Chierici, M. Beghi, E. Salina, E. Brega, COBRAEN, an updated version of the COBRA3C/MIT Code for thermalhydraulic transient analysis of light water reactor fuel assemblies and cores, Report 1010/1 (revised 1.9.99), ENELCRTN Compartimento di Milano (2001) [Google Scholar]
 HighPerformance MonteCarlo Methods for SAFEty Demostration – From proof of concept to realistic Safety Analysis and Industry Applications, Grant Agreement ID 755097. https://cordis.europa.eu/project/rcn/211032/factsheet/en [Google Scholar]
 R.N. Bratton, M. Avramova, K. Ivanov, OECD/NEA benchmark for Uncertainty Analysis in Modeling (UAM) for LWRS – summary and discussion of neutronics cases (Phase I), Nucl. Eng. Technol. 46, 313–342 (2014) [CrossRef] [Google Scholar]
 M. Faucher, Coupling between Monte Carlo neutron transport and thermalhydraulics for the simulation of transients due to reactivity insertions, PhD thesis, University of ParisSaclay, 2019 [Google Scholar]
 T. Hajas, G. Tolnai, D. Légrády, Variance analysis of the coupling of thermalhydraulics and point kinetics with stochastic noise term modeling dynamic Monte Carlo behavior for later use in GUARDYAN, in International Conference on Physics of Reactors 2022, 2022 May 15–20, Pittsburgh [Google Scholar]
Cite this article as: David Legrady, Marton Pukler, Balint Panka, Gabor Tolnai. Coupling the GUARDYAN code to subchanflow, EPJ Nuclear Sci. Technol. 9, 11 (2023)
All Figures
Fig. 1. Power distribution in the core of Paks NPP Unit 4 campaign 32rod drop experiment beginning state (left), halfway (middle), and end state (right). 

In the text 
Fig. 2. Coupling scheme of GUARDYAN Dynamic (GD) mode to thermalhydraulics (TH) solvers. 

In the text 
Fig. 3. Coupling scheme of static calculation (GUARDYAN, GK) and thermalhydraulics (TH). 

In the text 
Fig. 4. TMI Minicore crosssection material models. 

In the text 
Fig. 5. Source convergence in dynamic mode, power is normalized to a single starting neutron. 

In the text 
Fig. 6. L_{2} norm evolution of the GD source convergence for Run #1 (red) and Run #2 (blue). 

In the text 
Fig. 7. Pinwise fuel temperature distribution at the converged state, for Run #2 (left) and Run #1 (right) accumulated during the corresponding time step. 

In the text 
Fig. 8. Nodewise power estimates for Run #2 (left) and Run #1 (right) accumulated during the corresponding time step. 

In the text 
Fig. 9. Dynamic simulation average power, fuel, and coolant temperatures, for Run #1 (red) and Run #2 (blue). 

In the text 
Fig. 10. Convergence in static GK mode, evolution of the criticality (k _{eff}), fuel, and coolant temperatures with the coupling steps. 

In the text 
Fig. 11. Relative errors of the nodewise power estimates in percent for the lowest (left), central (middle), and topmost (right) nodes. 

In the text 
Fig. 12. Correlation of standard deviation of power and fuel temperature. 

In the text 
Fig. 13. Comparison of dynamic and static convergence results: relative difference of power output (left), relative difference of coolant temperatures (middle), and relative difference of fuel temperatures (right). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.