Simulation of a research reactor reactivity transient with deterministic and GPU-assisted Monte Carlo reactor kinetics codes

Reactor kinetic codes are crucial in safety assessment. Validating spatial and high temporal resolution kinetic solvers without thermal feedback is problematic as measurements seldom involve detailed spatial and fine temporal resolution. Benchmarking of deterministic codes thus often resorts to code-to-code comparison against Monte Carlo codes, which can only recently treat direct time dependence. In this paper, we have attempted to compare results from the GUARDYAN directly time-dependent Monte Carlo code and the SEnTRi transient driver developed for the PARTISN deterministic transport code to low power transient measured at the BME Training Reactor. Code-to-measurement comparisons were successful, despite a major uncertainty in the actual timing of the reactivity insertion and withdrawal originating from the instrumentation of the pneumatic rabbit system. Code-to-code comparisons concluded that time dependence was correctly implemented in both GUARDYAN and SEnTRi; furthermore, a hypothetical scenario was set up involving an instantaneous insertion of a negative reactivity into the BME TR core in order to compare spatially and temporally dependent fluxes. The simulations demonstrated the appearance of higher-order modes, and results showed a relatively good match, although fidelity of the comparison could be further improved by reducing the statistical uncertainty of the results provided by GUARDYAN.


Introduction
The computational capacity of the modern high-performance computing systems allows for more and more detailed calculation of reactor kinetics involving spatial effects in high temporal resolution, which are beyond the validity of point-kinetics or quasi-static approximations. Besides the improved deterministic transport approximations (e.g., S N or P L methods), promising developments were also reported recently in the field of Monte Carlo reactor kinetics calculations [1][2][3][4]. Validated spatial-kinetics calculation tools can significantly contribute Focus Point on Advances in the physics and thermohydraulics of nuclear reactors edited by J. Ongena, P. Ravetto, M. Ripani, P Saracco. a e-mail: szieberth@reak.bme.hu (corresponding author) to the improvement in nuclear safety in the field of reactivity initiated transients, severe accidents modeling, etc. Deeper understanding of local power distribution after inadvertent reactivity insertions may help in the determination of appropriate safety limits. Other applications can include advanced simulation of reactor pulses which can be used for several purposes, e.g., testing fuel behavior during reactivity initiated accidents in the TREAT facility [5].
However, notwithstanding the comprehensive number of benchmark scenarios available for criticality calculations, many of them included in [6], it is hard to find an appropriate candidate for benchmarking the transient capabilities of kinetic solvers, especially regarding fast transients and spatial dependence. One reason for this is the extreme difficulties of measurements involving fast reactivity insertion and high-resolution sampling of the neutronics response in space and time. In the lack of measurement data, code-to-code comparison is often used. In such comparisons, Monte Carlo transport codes are often considered as reference due to the continuous energy representation and explicit 3D geometry description. Recent developments in Monte Carlo reactor kinetics may enable to perform transient calculations, but require an extensive validation of the codes.
While there are a number of benchmark problems prepared for deterministic codes with non-continuous energy cross sections, e.g., [4,[7][8][9], existing transient benchmarks on the Monte Carlo field are mostly proposed for coupled neutronics/thermal hydraulics calculations, without explicit time dependence and based on a series of static calculations or the improved quasi-static method as in, e.g., [10][11][12]. The need for defining a benchmark problem for time-dependent Monte Carlo simulations without feedback effects has only been recently realized. Efforts made so far in this respect only involved very simplified models such as [13,14].
The Institute of Nuclear Techniques (NTI) of Budapest University of Technology and Economics (BME) is involved both in the development of Monte Carlo and, in cooperation with Karlsruhe Institute of Technology (KIT), deterministic reactor kinetic methods. In an attempt to provide a validation basis for reactor kinetics codes, reactor physics experiments were performed at the Training Reactor of BME (BME TR). Due to limitations in instrumentation, these experiments cannot depart from the regime where point-kinetics or quasi-static approximations are valid but are well suited to verify that the codes are able to calculate global power in a complicated geometry during reactivity transients. Following a successful validation by measurement, a code-to-code benchmark was defined in the BME TR geometry with the instantaneous insertion of an absorber, where spatial effects can also be observed on a short timescale. This exercise is also considered as a scoping study to define potential transients where spatial effects are significant. Feasibility of the realization of such measurements providing benchmark data for spatial-kinetics codes is under investigation. Section 2 of this paper describes the measurements performed at BME TR for validation purposes. Sections 3 and 4 describe the computational methods and models applied in Monte Carlo and deterministic modeling, respectively. Finally, Sect. 5 presents the results obtained from the experimental benchmark and the code-to-code comparison based on a hypothetical case.

Description of the measurement and the benchmark model
A set of experiments was performed at the Training Reactor of the BME in order to produce reactivity transients for the validation of reactor kinetic codes. The experiments utilized the pneumatic rabbit system of the facility to insert and withdraw a piece of Cadmium absorber to and from the core, while the reactor power was recorded based on ex-core fission chambers. This section describes the BME TR, the performed experiment and its benchmark model. Further details which may be necessary for modeling can be found in technical reports [15,16] available on request from the authors.

The BME Training Reactor
The Training Reactor of the Budapest University of Technology and Economics (BME TR) is a pool-type, water-moderated nuclear reactor. It is located at the campus of the BME and has been operating since 1971, using EK-10-type fuel assemblies with 10% 235 U enrichment. The nominal maximum thermal power of the reactor is 100 kW, and the maximum thermal neutron flux is 2.7 · 10 12 n/cm 2 s. The reactor core is located in a cylinder-shaped tank of 1400 mm diameter and 20 mm thickness, made of Al-alloy, which is filled with demineralized water. The fuel in the reactor is UO 2 in magnesium matrix, and it is placed in cylinder-shaped Al-alloy cladding with a thickness of 1.5 mm. The fuel pins are 10 mm in diameter and 500 mm in active length. At each end of the fuel pins, there are inactive parts with the length of 45 mm made of the Al-alloy of the cladding. The current core of the BME TR (see in Fig. 1) contains a total number of 24 fuel assemblies that contain a total of 369 fuel pins. Besides the 9 regular fuel assemblies with 16 pins in a 17-mm-pitch regular square lattice, 15 fuel assemblies of 5 different designs with missing pins or distorted lattice are used in order to form space for control rods and irradiation positions. The core also contains 41 solid graphite blocks with aluminum cladding and a holder element, which serve as the reflector around the core. Some of the graphite blocks are penetrated by aluminum elements in the cut of the horizontal beam tubes to increase the beam intensity. In the corner assemblies, graphite blocks are cut diagonally, the remaining space is filled by air. The irradiation position of the thermal rabbit system is located in the reflector surrounded by water to maximize the thermal flux. The fast rabbit system is located in the middle of the core replacing a half fuel assembly and is surrounded by air in order to increase the fast flux in the irradiation position.
The numerous irregularities and inhomogeneities in the core make it challenging to model using the homogenization techniques applied for deterministic modeling.

Experimental procedure
During the experiments, the pneumatic rabbit system was used to insert a small absorber into the core, and observe the change in reactor power.
The absorber used in the experiments was a ring made of metallic cadmium. The height of the sample was 40 mm, the outer diameter was 9 mm, and the thickness was 0.5 mm. The ring was placed in a polyethylene casing which is regularly used with the pneumatic rabbit system. The rabbit may be considered as a cylinder of 17 mm outer diameter, 3 mm thickness and 58 mm height-including 4 mm and 5 mm solid caps at its top and bottom, respectively.
At the start of each measurement, the reactor operated at critical state. In order to avoid thermal feedback effects, all experiments were carried out at low power with a maximum of 500 W. The absorber ring was inserted into the core using the 'fast' section (see in Fig. 1) of the pneumatic rabbit system and induced a transient behavior. The bottom of the cadmium ring was 2 cm below the midplane of the core when inserted.
Due to the negative reactivity of the inserted ring, the reactor power started to decrease exponentially in time after a prompt jump. The transient was observed for a few seconds before the absorber was removed, and the reactor returned to a critical state at lower power than at the initial state. The process was repeated for a number of times at different initial power states.  1 Assembly map of the BME TR. The absorbent material was inserted into the 'fast' rabbit system, marked in the figure. Detectors "I1" and "I2" were used during the measurements During the transient, TTL-format signals of CFUL08 fission chambers in ex-core detector tubes "I1" and "I2" were recorded in time-stamp mode and then it was converted into countrate functions with 1 ms time bins. In order to correct for the deadtime effect, the time-interval distribution of the counts was plotted, which showed τ = 150 ns deadtime, and paralyzable deadtime model was used [17].

Benchmark model
The currently available and installed instruments do not make possible to record or to determine the exact position or speed of the rabbits in the pneumatic system. Therefore, some assumptions were made based on the measured power curve in order to make a suitable estimation for the position of the rabbit. Linear functions were fitted over different regions of the normalized count-rate function, and then the intersections were determined. The start and end times of these selected regions are found in Table 1 and the fitting in Fig. 2 for data series III/2 chosen for the present exercise. The intersection of the linear functions was assumed to be the departure or arrival time of the sample from or to the given position which are presented in Table 2. Constant velocity values were assumed and determined based on the traveling time and the traveled distance.   [3] is a GPU-based Monte Carlo code with direct time dependence for applications to nuclear reactor transients. The time evolution of neutron histories is followed directly rather than by a set of quasi-static calculations in consecutive time steps. The GUARDYAN calculation flow has many similarities to the DMC methodology proposed and developed in [1,2] with many of the techniques revised and adjusted to better fit the GPU architecture. These include the time step methodology instead of the classical generation-to-generation tracking, the improved branchless method [18] preventing history branching and thus banishing large particle banks, and the forced decay of precursors ensuring that prompt fission chains are never undersampled. The most important GPU-specific features of the methodology include the complete elimination of neutron banking, the use of Woodcock (delta-) tracking with a supervoxel approach, and the lack of trajectory splitting and Russian Roulette using an importance-based population control between time steps instead. In transient scenarios, perturbations to cross-sectional data, material composition, or geometry must be evidently accounted for in order to model thermal feedback or geometry changes. While the latter could be calculated on-the-fly, considering TH feedback is only possible if neutron histories are synchronized in time. GUARDYAN therefore splits the total time span of the simulation into short intervals. A GPU function termed the kernel function is called periodically at time boundaries to simultaneously process the parts of neutron histories falling in the upcoming time step. Particles are always synchronized at boundaries. Modification of geometry or thermal feedback effects can only be considered at time boundaries. Note that typical time steps for thermal systems in GUARDYAN are in the order of tenth or hundredth of milliseconds, thus feedback/movement mechanisms can be nicely modeled step by step. While this time step methodology is mainly motivated by the change in the physics of the system, it also gives opportunity to observe and interfere with a population that is not dispersed in time or to rearrange computational resources when any changes to the system is made. The initial conditions for a time-dependent simulation are obtained assuming the transient starts from steady state, and thus both neutron and precursor source can be calculated by a static calculation.
While prompt neutrons have an average lifetime in the order of 10 µs (for thermal systems), the expected lifetime of a delayed neutron can range from 0.1 to 100 s. The difference is several orders of magnitudes, and thus analog simulation of delayed neutrons is clearly unfeasible in a time-dependent simulation. In practice, this problem would arise in two ways: First, prompt fission chains would disappear due to the heavy undersampling of delayed neutrons. Second, a delayed neutron produced at a certain time would only participate in the simulation after many generations of prompt neutrons, taking up valuable computational resources in the meantime. Therefore, GUARDYAN treats precursors in a nonanalog way with a MC precursor sample representing a population of precursor atoms by introducing precursor statistical weight. The optimal number of precursors to neutron samples is not trivial to guess, as it changes with the time evolution of the system. In GUARDYAN, a scheme has been developed that allows robust delayed neutron treatment without taking time-dependent importance into account, i.e., the prompt populations and the precursor population adapt automatically even if radical reactivity changes occur.

Data, geometry specification and dependencies
The input logic of GUARDYAN is similar to the logic followed by production-level codes MCNP [19], Serpent [20] or OpenMC [21]. Cells are assumed to be homogeneous and can be defined by second-order bounding surfaces and Boolean operators. Transformations, hexagonal and rectangular lattices, and universes can also be given in an input. Materials are defined by isotopic composition (atomic ratio or mass fraction) and density.
GUARDYAN uses cross-sectional data in the ACE (A Compact ENDF) format, which is generated by NJOY [22] from the ENDF-B-VII.1 cross-sectional library. Other evaluations can also be used to create ACE format for GUARDYAN, e.g., data from JEFF libraries. Isotopic total and majorant cross sections are pre-calculated on a unionized energy grid. Geometry and nuclear data in GUARDYAN are stored using double-precision floating point format. Although single-precision performance of commercial GPU cards highly exceeds that of double precision, in most MC reactor physics applications accuracy of computations is heavily dependent on the precision of nuclear data. Single-precision round-off error could easily lead to the violation of physical laws when sampling tabular distributions in the ACE format, or ignore small details in geometry and yield biased results.
Interactions are modeled similar to MCNP, OpenMC or Serpent; post-collision energy and angle distributions are sampled according to ACE laws. S(α, β) tables are considered for scattering of thermal neutrons on bounded atoms; to deal with unresolved resonances, GUARDYAN uses the probability table method.
Due to GUARDYAN being a novel time-dependent Monte Carlo code, the correct implementation of physics was challenged by comparing GUARDYAN to the general purpose Monte Carlo code MCNP6 [19]. In this code-to-code comparison, differential quantities were analyzed, overall 445 000 data points were compared. In roughly 2000 separate runs in spherical geometry including different isotopic composition and launching neutrons at t = 0 with various starting energies, we tallied the time evolution of energy-dependent flux on the outer surface of a sphere. Only one time step was executed, and tallies were then binned into 9 time bins. Integral comparisons were done on criticality benchmarks. Results compared very well [23] with results of MCNP6.

Modeling the Training Reactor
The GUARDYAN model of the BME Training Reactor followed closely the benchmark model described in [16]; a cross-sectional image of the model geometry is plotted in Fig. 3, consisting of 137 surfaces, 150 cells and 56 isotopes mixed into 13 different materials.
As GUARDYAN does not use correction with the effective multiplication factor to stabilize apparent power level, therefore the reactivity of the system should be set to zero very accurately. For that reason, the automatic control rod of the BME TR had to be adjusted 3 cm higher than specified in the benchmark description. Source convergence was achieved by an initial run starting from a homogeneous 80 cm × 80 cm × 80 cm cube and running the simulation for 0.02s real time but simulating delayed neutrons with instantaneous precursor decay. Given this neutron source, precursors were generated by where C i stands for the concentration of precursors from the family i; β i and λ i are the delayed neutron fraction and decay constant, respectively, Σ f is the fission cross section, v denotes particle velocity, ν stands for the number of fission neutrons and n denotes the density of live neutrons. A further 0.1s real-time calculation was carried out with precursor decay sampled from the realistic exponential distribution to allow for settling of any statistical instabilities. Transient simulations were started from this distribution. Number of particles simulated was about 4 million (2 22 consisting of 2 21 prompt neutrons, 2 20 number of precursors and 2 20 items as a precursor presampling buffer). Running times were about 51 hours per 1s real time on a single Nvidia GeForce GTX 1080 graphics card.

The deterministic reactor kinetics code and model
The PARTISN [24] discrete ordinates neutron transport code developed at Los Alamos National Laboratory (LANL) was used as transport solver for the deterministic transient simulations, which includes a time-dependent option based on the semi-implicit timediscretization scheme without taking into account the delayed neutrons. The code was extended at the Karlsruhe Institute of Technology (KIT) by implementing a theta-scheme, which offers both semi-implicit and fully implicit time-discretization options, and also a treatment of delayed neutrons was included [25]. Furthermore, the use of the quasi-static scheme is now available with PARTISN when it is used with ERANOS [25] and SIMMER codes [26]. However, the simulation of the present transient was performed by a first-order fully implicit time-discretization scheme developed at BME NTI in cooperation with KIT. In the following sections, the computational model, the calculation flow of the SEnTRi transient driver and the PARTISN model of the BME TR are presented.

Numerical method for solving the time-dependent transport equation
The implemented scheme is based on the direct method presented in [27], which is an unconditionally stable first-order fully implicit time-discretization scheme. With the usual notations, the multigroup time-dependent transport equations for angular flux ψ g (r, Ω, t) of energy group g can be written in the following form: 1 v g ∂ ∂t ψ g (r, Ω, t) + Ω∇ψ g (r, Ω, t) + σ t,g (r, t)ψ g (r, Ω, t) = g 4π σ s,g →g (r, Ω → Ω, t)ψ g (r, Ω , t)dΩ where Q g represents the external source, λ i is the decay constant, β i is the delayed neutron fraction and χ i,g is the delayed neutron spectrum in group g for delayed neutron precursor family i. The delayed neutron concentrations C i can be determined as: and the fission source F(r, t) is defined as: where γ is a normalization factor. Since it is assumed that the reactor is in steady-state critical condition at t = 0, the initial conditions can be determined from the eigenvalue equation where γ equals to the eigenvalue: Assuming that ψ g and F change linearly during a time step between t 0 and t 1 = t 0 + t, the time derivative can be expressed as: and C i (r, t 1 ) can also be expressed with C i (r, t 0 ), F(r, t 0 ) and F(r, t 1 ). With these assumptions, Eq.
(2) at t = t 1 can be cast in the following form: where modified total cross section σ t,g (r, t 1 ), fission spectrum χ g (r, t 1 ) and external source Q g (r, Ω, t 1 ) have to be used: Therefore, the time-dependent solution can be obtained from a sequence of time-independent fixed-source problems, where the angular flux, fission source and delayed neutron concentration from the previous time step appear in the external source of the subsequent time step.

The implemented transient driver
The calculational flow of the transient driver based on the direct method implemented in the SEnTRi code of BME NTI is shown in Fig. 4. The PARTISN code is used as transport solver for the simulations. In order to initialize the transient simulation, an eigenvalue calculation is performed for the initial state assuming that the reactor was in steady-state condition until Fig. 4 The workflow of the transient module of the SEnTRi code the beginning of the simulation. After the PARTISN eigenvalue calculation, which solves Eq. (5), the SEnTRi code reads its input and the corresponding binary files: the GEODST geometry description, and the RMFLUX angular flux distribution files. Then it determines the initial delayed neutron precursor concentration based on Eq. (6). According to the previously collected data, the code prepares the following time step by calculating σ t,g (r, t 1 ), χ g (r, t 1 ), and Q g (r, Ω, t 1 ) based on (9). These parameters are passed to the PARTISN via the crosssectional library and the binary external source file FIXSRC, which are updated in each time step by SEnTRi. Afterward, the PARTISN code is called to perform a fixed-source calculation according to Eq. (8), which will create a new RMFLUX file with the current flux distribution. The SEnTRi code reads the RMFLUX file, updates the delayed neutron precursor concentrations and restarts its cycle. The temporal behavior of the angular flux can be obtained from the outputs of the fixed-source calculations. Since PARTISN is capable of reading angular source distribution in spherical harmonics representation, the calculated angular flux of the given time step has to be converted from discrete ordinates to spherical harmonics when used for the external source of the next step. This conversion calls for higher-order angular representation to avoid numerical instabilities.

The deterministic model of the BME Training Reactor
The European Cell COde (ECCO) from the ERANOS [28] program package was used for the preparation of two-group homogenized macroscopic cross section for the BME TR model. The homogenization calculations were performed with a 1986-group ultra-fine JEFF3.1 data library. Separate cell calculations were performed for the regular fuel lattice and for the Vertical and horizontal midplane cross section of the homogenized BME TR model irregular assemblies. Graphite, aluminum and water reflector regions were treated in infinite homogeneous approximation. Special treatment was necessary for the absorbers (control rods and the Cd sample). In the first step, a heterogeneous one-dimensional cylindrical model was created from the actual absorber elements and the neighboring structure, like the guiding tube. This heterogeneous region was surrounded by several layers of fuel and mixture of water, wrapper and cladding material in order to create a representative neutron spectrum. Reactivity worth of the absorber was determined in this model by eigenvalue calculation with the absorber region included and voided. Based on this, an equivalent homogenized composition was created which contained only one homogeneous absorber region and one homogeneous material for the surrounding structures while having the same reactivity worth due to the properly set absorber number density. An example of the heterogeneous and the equivalent models is shown in Fig. 5. The complete core model was built in a finite-differences mesh from the homogenized materials. The final deterministic model is shown in Fig. 6.
Kinetics parameters had to be also determined in order to perform transient calculations. The delayed neutron fraction and the neutron spectra for each family were determined with the help of the ERANOS program package. Average group velocities were calculated in  Fig. 7 Comparison of the different calculation methods for the instantaneously inserted sample 40 groups by the ERANOS and collapsed into two groups using the 40-group spectrum calculated by PARTISN.

Results
Using the above-described models, simulations were performed by both codes for a selected measurement and for hypothetical cases in order to compare the capabilities of the methods.

Estimation of the reactivity worth
Since the given transient is mainly determined by the negative reactivity inserted into the core, calculations were performed for the instantaneous insertion of the absorber in order to determine the reactivity based on the prompt jump approximation and compare the results with point-kinetics. The reactivity worth of the sample was also determined form the difference of static eigenvalue calculations both with PARTISN and MCNP6. A reasonable agreement was found, but PARTISN predicts higher reactivity worth, as shown in Table 3. The comparison of the transient simulations and the point-kinetics estimation based on the reactivity worth estimated by MCNP6 for the instantaneously inserted sample is shown in Fig. 7. In the deterministic transient calculations 10 µs time steps and S 12 approximation were applied. The GUARDYAN follows the point-kinetics within statistics, while the SEnTRi estimates a slightly larger reduction in the power, due to the higher reactivity worth calculated by PARTISN. The simulation of the complete transient was performed using the parameters specified in Sect. 2.3. In the deterministic simulation, 500 µs time step and S 12 approximation were used and the cadmium was inserted and removed from the system in 14 steps of different sizes, as the axial mesh allowed for it. In the GUARDYAN simulations, the Cd absorber was inserted in 10 steps to the system. The global relative power was calculated and compared to the measurements. The results are shown in Figs. 8 and 9. As Fig. 8 shows, both simulations reproduce the global transient behavior with satisfactory agreement. The PARTISN results slightly (in an extent of 1-2%) overestimate the drop in the power due to the overestimated reactivity worth as described in Sect. 5.1. Higher-resolution results are presented for the insertion and removal phases in Fig. 9a, b, where a delay can be observed in the simulations compared to the measurement. This discrepancy can be attributed to the uncertainties concerning the movement of the rabbit and the approximations of the benchmark. For the sake of simplicity in the modeling, constant velocities were assumed for the insertion and withdrawal, while in the reality an acceleration phase must be present. Due to these approximations, the exact reproduction of the measurement results cannot be expected. On the other hand, the two calculations using

Code-to-code comparison
In the above experiments, the point-kinetics could give an almost-perfect estimation for time dependence of the power. In order to test the capabilities of the codes for cases where a considerable change exists in the flux shape during the transient, a hypothetical case was defined in the TR geometry. In this case, a 20 cm long cadmium sample with the polyethylene casing was inserted instantaneously next to the core into the thermal pneumatic rabbit system. The calculated reactivity worth of the sample is shown in Table 4. Fission power was calculated in several fuel assembly quarters (see in Fig. 10) in a 10 cm long section around the midplane to obtain spatial resolution. Since the Monte Carlo code handles detailed geometry and cross sections, such comparison is also useful to test the validity of the homogenized model applied for deterministic calculations. The result of the deterministic calculations is shown in Fig. 11. One can observe that from around 500 µs the curves have the same slopes, which means that the flux shape does not change anymore. On the other hand, in the first part of the curves a significant difference can be observed. Based on these results, the code-to-code comparison was performed for the first 300 µs.
The tallied power changes at the selected locations were in the order of a few percents and had to be locally estimated. Both properties pose a challenge for a Monte Carlo code, especially if direct time dependence is concerned as the statistical uncertainty scales inversely with the number of fissions occurring in a time interval in the specified spatial region. The GPU implementation restricts the number of starting neutrons by the physical memory of the graphics card; therefore, the simulation of the 300 µs long interval was repeatedly simulated starting from the same source distribution but using different random seeds, and the obtained 700 curves were averaged. In the deterministic simulation, 1 µs time steps and S 12 approximation were applied. Due to the small time step applied, the calculation converges only with at least such fine angular resolution.  The results from both simulations are compared in Fig. 12. The figures suggest good agreement as the two estimations agree within the statistics of the Monte Carlo results. One may observe that according to the deterministic calculation, the power in the second fuel region drops only after 40 µs, which is most probably due to the moderator inserted together with the absorber. The Monte Carlo results appear to confirm this finding, although the statistics is not good enough to provide a solid proof.
Despite this, the statistical uncertainty of the GUARDYAN results is considered to be low, since the typical goal of the dynamic Monte Carlo calculations is to produce estimates averaged over 0.1 s within 1% error, in order to make thermal-hydraulic coupling stable. On the other hand, in the simulations presented in Sects. 5.2 and 5.3 the tallies were binned into 1 ms and 1 µs long time bins, respectively. The statistical error of the Monte Carlo results can be manipulated by many parameters, from which the most trivial ones are the number of starting particles and the length of time bins. Using a wider bin would mean the averaging of results over a longer period and statistical uncertainty would drop; however, the required temporal resolution could not be reached. Regarding the number of particles, the uncertainty generally follows a 1 √ N behavior, where N is the number of starting samples. Therefore, an order of magnitude lower statistical uncertainty would require 100 times longer calculation time, which is unfeasible with the computational capacity available.
An important conclusion from this exercise is that the sudden insertion of about 10 cents negative reactivity next to the core causes almost 10% deviation in the flux shape in the nearest fuel pins and less then 1% in the middle of the core. These deviations decay in about 0.5 ms, and after that, the global flux shape in the core does not change and the point-kinetics approximation is fully valid. These observations indicate the required spatial and temporal resolution for an experiment that could provide a validation base for kinetics codes with full spatial dependence.

Conclusion
A benchmarking attempt was made using measurement data from BME TR of two kinetic solvers to compare neutron fluxes with high temporal and spatial resolution. The difficulty in direct measurement of such quantities calls for a method where a possibly approximation-free code may provide flux values with high temporal resolution. The directly time-dependent GUARDYAN Monte Carlo code may provide with such a tool; however, since these code developments only emerged in recent years, GUARDYAN itself needs proper validation. Measurement-to-code comparisons for both GUARDYAN and SEnTRi yielded good agreement, but imprecise timing information of the reactivity insertion left minor discrepancies present. The most important conclusion from the benchmark exercise is that both codes handle properly the global kinetics processes in the very complicated, highly heterogeneous geometry. Since the measured case did not involve the appearance of the higher moments in the flux distribution and high temporal resolution, a code-to-code comparison was performed involving the instantaneous insertion of an absorber into the BME TR core in a hypothetical scenario. The deterministic calculations for the transient showed a strong spatial dependence in the fuel assemblies close to the absorber during the first 0.5 ms of the transient. Results of GUARDYAN and SEnTRi agreed well at different locations regarding the time dependence, but high statistical uncertainties hinder drawing very strong conclusions. The proposed benchmarking methodology could be improved with fast reactivity insertion measurements with precise timing information and also by decreasing Monte Carlo variance.
Funding Open access funding provided by Budapest University of Technology and Economics (BME).
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.