© C. Montecchio et al., Published by EDP Sciences, 2025

1. Introduction

Monte Carlo methods are traditionally considered the gold standard for stationary neutron-transport problems in reactor physics applications. In recent years, they have been also successfully extended to reactor kinetics (including delayed neutron precursors contributions), a domain previously handled exclusively by deterministic solvers [1–9]. This improvement has been made possible thanks to the increasing availability of computational resources, combined with advancements in variance-reduction techniques [10, 11], prompted by the pioneering work by Sjenitzer and Hoogenboom [12–17]. The development of efficient Monte Carlo methods for time-dependent problems is essential for the numerical validation of deterministic kinetic solvers, which rely on the discretization of the phase space and thus involve model-dependent approximations; in this respect, kinetic simulation capabilities have been implemented in the TRIPOLI-4^® Monte Carlo code [18, 19], developed at CEA, including several variance-reduction and population-control methods [3, 10], and extensive benchmark comparisons have been carried out [2, 4, 20].

Despite these advances, kinetic Monte Carlo methods are still in their infancy and face some serious limitations in terms of Figure of Merit (FoM), since they typically require extremely long computer times in order to achieve acceptable statistical accuracy on the sought response functions. These considerations severely hinder the simulation of transients involving large reactivity insertions or withdrawals. One of the toughest challenges in time-dependent Monte Carlo simulations arises from the imbalance between neutron and precursor population sizes, which stems from their largely different respective time-scales [13]. Indeed, the characteristic time-scale of the precursor population (i.e. the average decay time $1/\overline{\lambda }$) is much longer than the typical neutron time-scale (i.e. the mean generation time Λ). This is mirrored in the ratio $\beta /(\overline{\lambda }\mathrm{\Lambda })$ between the two population sizes at critical conditions: for a typical PWR, a snapshot of the core at equilibrium would show 10⁴ more precursors than neutrons; the ratio would exceed 10⁵ for a typical fast reactor [13]. This substantial difference, possibly coupled to the unbounded growth or disappearance of the particles during a transient, imposes significant constraints on kinetic Monte Carlo simulations in terms of memory burden and computer time.

Extreme care is therefore required to manage both the population size and its evolution while addressing the neutron-precursor disparity. This is typically accomplished using population-control methods, which keep the overall particle population within reasonable bounds while optimizing the allocation of computational resources [10]. Building on previous investigations concerning the use of the combing strategy [21], we have recently implemented a preliminary version of an adaptive population control algorithm in TRIPOLI-4 [22]. In this work, we considerably extend our findings and show how combing can be made highly efficient by using an adjoint-based importance-sampling scheme: a thorough numerical assessment of the importance-dependent strategy on a collection of benchmark configurations will substantiate our claims.

This manuscript is organized as follows: in Section 2 we briefly recall the key features of the legacy implementation of the kinetic algorithm in TRIPOLI-4, focusing primarily on the population control methods. In Section 3 we introduce an importance-dependent adaptive population control strategy, which considerably enhances the effectiveness of the method proposed in Ref. [22]. Then, in Section 4 we examine the FoM of the legacy and importance-dependent algorithms relying on the Flattop-Pu [23] and on the SPERT III E-core benchmark [24]. Further analyses are carried out in order to probe the performance of the population-control methods for scenarios corresponding to large reactivity and long transients, relying on the Flattop-Pu and on the CROCUS benchmark [25]. Conclusions are finally drawn in Section 5.

2. Overview of the kinetic capabilities of TRIPOLI-4

We begin by briefly recalling the key features of the legacy algorithm currently implemented in TRIPOLI-4. For an exhaustive description of the kinetic capabilities of the code, the reader is referred to the original work by Faucher et al. [3].

In TRIPOLI-4, the kinetic simulation scheme is divided into two main stages, as illustrated in Figure 1. The first stage is used to sample the critical distribution of neutrons and precursors, if needed. It consists of a power iteration over a number of cycles, called inactive cycles, because they do not contribute to the final tallies. Their number must be large enough to ensure convergence of the fission source. The power iteration is performed on the equilibrium configuration of the simulated system; once convergence is attained, the neutron and precursor populations are sampled within a further power iteration cycle. These populations then become the initial condition for the second stage, where the particles are simulated in time-dependent conditions, over a user-defined time mesh used for tallying, and the geometry and the material compositions of the system are updated according to the chosen kinetic scenario. This two-stage scheme is repeated for a sufficiently large number of batches, and the final result reported by TRIPOLI-4 is obtained by averaging over them. In principle, each batch should be an independent calculation including both a source-convergence and a kinetic phase, but this is computationally demanding. To reduce this cost, the fission neutrons produced in the last inactive cycle of one batch are reused to initialize the next. This approach reduces overhead but introduces correlations between batches, which are mitigated by adding a few additional de-correlation inactive cycles. In parallel runs, each task executes the scheme shown in Figure 1.

Fig. 1. Sketch of the current implementation of the kinetic algorithm in TRIPOLI-4, taken from Ref. [3]. In the kinetic simulation (green) neutrons and precursors are sampled during the last iteration of the power iteration (white) and their subsequent time evolution is followed on a time grid, during which forced precursor decay and population control algorithms are applied.

Due to the combined effect of nuclear data libraries and technological data, the effective multiplication factor k_eff of a simulated system may deviate from unity for a configuration that is experimentally observed to be critical. Correspondingly, the associated fundamental eigenmode φ_k will also deviate with respect to the actual stationary flux. When performing a kinetic simulation starting from equilibrium conditions, it is therefore very important to correct this bias, otherwise we would observe a departure from the stationary condition for the neutron and precursor population of the supposedly critical system. For this purpose, in TRIPOLI-4 we normalize the average number of neutrons generated by fission by k_eff, forcing the system to be critical.

2.1. Non-analog precursor sampling

Precursors are treated as a distinct particle. Currently, precursors retain the spatial coordinates of the neutrons from which they are created, i.e. their migration is neglected. The large difference in the time-scales of neutrons and precursors, which would lead to intractable issues in a naive time-dependent simulation, is taken into account by suitably altering the sampling laws of precursors, and correspondingly adjusting their statistical weights in order to preserve the unbiasedness of the tallies, based on a strategy proposed by Legrady and Hoogenboom [12]. Precursors belonging to different families are grouped together in a single representative precursor particle, which reduces the variance due to the dispersion of the decay constants λ_j. Furthermore, the decay time distribution of such representative precursor is forced to be uniform within every simulated time step [t_i, t_i + 1], namely,

$$\begin{array}{c}\hfill p_{\text{decay}}(t)=\frac{1}{t_{i+1}-t_i}\chi (t_i,t_{i+1}),\end{array}$$(1)

where χ(a, b) is the marker function of the interval [a, b], and i is the index of the time step. This enhances the emission of a delayed neutron in each time interval, compared to the analog exponential distribution of the decay times. To preserve an unbiased Monte Carlo game, the statistical weight of the delayed neutrons emitted upon precursor decay must be modified accordingly, i.e.

$$\begin{array}{c}\hfill W_{\text{delayed}}(t|t_0)=(t_{i+1}-t_i)\sum _j\frac{{\lambda }_j{\beta }_j}{\beta }{e}^{-{\lambda }_j(t-t_0)}{w}_C,\end{array}$$(2)

where w_C is the initial weight of the emitting precursor, i.e. the weight of the neutron from which the precursor was generated, β_j is the fraction associated with the delayed neutron precursor family j (with β = ∑_jβ_j), and t₀ is the precursor creation time. Upon decay, the precursor is not killed but instead added to a buffer of particles to be followed in the next time step i + 1, with a weight equal to:

$$\begin{array}{c}\hfill W(t_{i+1}|t_0)=w_C\sum _j\frac{{\beta }_j}{\beta }{e}^{-{\lambda }_j(t_{i+1}-t_0)}.\end{array}$$(3)

2.2. Population control algorithms

If the particles are left free to evolve during a kinetic transient, the population size over time would either grow unbounded if the system is supercritical, or shrink to zero if the system is subcritical, both conditions being unacceptable for the Monte Carlo simulation. Even a critical system, where the population stays constant on average, would lead in practice to intractable conditions for the simulation, since the variance would increase linearly in time because of the fluctuations induced by the fission chains [26]. Furthermore, the significant disparity between neutron and precursor time-scales creates an imbalance among the two populations, favoring precursors over neutrons; since tallies in the simulation are carried by neutrons, this situation is again detrimental to the overall variance of the sought responses.

The set of approaches designed to tackle these issues is collectively referred to as population control strategies. The goal of such techniques is to maintain a suitable number of particles to be simulated while also addressing the pronounced disparity between neutrons and precursors, thus ensuring acceptable statistics. In TRIPOLI-4, population control strategies are applied on the same time grid as forced decay, minimizing the need to define new grids. For convenience, we can distinguish two main procedures in a population control algorithm:

1. The scaling of the neutron and precursor populations according to an importance-sampling scheme;

2. The application of either combing to both neutrons and precursors, or Russian roulette specifically to precursors to control the actual population size at the end of each time step.

The role of the first procedure is to promote the allocation of computer resources for the neutron population in the transport simulation, in view of the previous considerations about tallies being sensitive to neutrons only. In this work we have extensively modified this procedure with respect to the legacy implementation in TRIPOLI-4 described in Ref. [3]: a thorough description of the novel implementation is provided in Section 3.

The role of the second procedure is to enforce the control of the population sizes, in order to restore a convenient amount of particles for the subsequent time step, based on the particles that survived until the end of the current time step.

If combing is chosen, the number of particles is set to a target value M, while their total weight and the proportion of respective weights (in the case of multiple particle species) are preserved. The new M particles are sampled among the K current particles by using two equally spaced weight intervals, one for neutrons and one for precursors, over the total available weight. It has been shown that this algorithm is unbiased, provided that a random offset is imposed when selecting the first particle [21].

If Russian roulette for precursors is chosen, then the population size of precursors is controlled through the comparison with a reference weight: if the weight of the particle is lower than a predefined threshold, it is subjected to a random survival test. The roulette is applied on the basis of the expected delayed neutron weight w_delayed that can be emitted during the following time step [t_i, t_i + 1], which follows from [27]

$$\begin{array}{cc}& w_{\text{delayed}}(t_i,t_{i+1}\mid t_0)\hfill \\ \hfill & =w_C\sum _j\frac{{\beta }_j}{\beta }(e^{-{\lambda }_j(t_i-t_0)}-e^{-{\lambda }_j(t_{i+1}-t_0)}).\hfill \end{array}$$(4)

This quantity is compared to a threshold value, which in TRIPOLI-4 is set by default to 0.8. A detailed comparison of the performance of combing versus Russian roulette for population control will be provided in Section 4.

In addition to either combing or Russian roulette at the end of each time step, in TRIPOLI-4 Russian roulette and splitting are also applied to neutrons within each time step after every collision, as a standard variance-reduction technique.

2.3. Legacy importance-sampling scheme for neutrons and precursors

The kinetic algorithm of TRIPOLI-4 features an importance-sampling scheme, whose legacy implementation was described in Refs. [3] and [10]. In a recent work, we have introduced a preliminary version of an enhanced importance-sampling scheme for kinetics calculations [22], whose comprehensive description is provided in Section 3. The key motivation behind such importance-sampling scheme is to adjust the weights of particles and tame the imbalance between neutron and precursors, prioritizing neutrons in memory allocation. For the sake of clarity, we stress that the notation used in this paper differs from the one adopted in Faucher et al. [3, 10]: we use 1/I_n instead of I_n, and 1/I_c instead of I_c, respectively. The importance ratio R = I_n/I_c is still defined as the ratio between the neutron importance and the precursor importance, but the values of the absolute importance are altered. The reason for this change, which was suggested in Ref. [28], is that the new definition yields importance ratios significantly larger than one for short time steps, which aligns better with the intuitive understanding that neutrons play a more significant role than precursors on short time-scales. Moreover, this definition is more consistent with the interpretation of importance as the expected contribution to the detector response and as the solution to the adjoint Boltzmann equation.

We have then to distinguish physical and simulation weights. The statistical weights assigned to neutrons and precursors in the simulation, denoted respectively as w_n and w_c, are referred to as physical weights. However, to optimize the computational handling of neutron and precursor populations, we replace them with the simulation weights, ${\tilde{w}}_n$ and ${\tilde{w}}_c$. In the legacy algorithm, the simulation weights ${\tilde{w}}_n$ and ${\tilde{w}}_c$ are derived from the physical weights w_n and w_c through the following relationships:

$$\begin{array}{cc}\hfill {\tilde{w}}_n& =w_n·I_n\hfill \end{array}$$(5a)

$$\begin{array}{cc}\hfill {\tilde{w}}_c& =w_c·I_c.\hfill \end{array}$$(5b)

More specifically, the legacy algorithm is structured so that I_c = 1 and I_n = R = I_n/I_c, meaning that only the neutron weights are effectively corrected, while the precursor weights remain unmodified. The sum of the physical weights is required to be equal to the sum of the simulation weights, ensuring the conservation law:

$$\begin{array}{c}\hfill w_n+w_c={\tilde{w}}_n+{\tilde{w}}_c.\end{array}$$(6)

Using the constraint in equation (6) and referring to the definitions of ${\tilde{w}}_n$ and ${\tilde{w}}_c$ in equation (6), we can derive the explicit expressions of the neutron importance I_n and of the precursor importance I_c as a function of the importance ratio R, namely:

$$\begin{array}{cc}\hfill I_n& =R·\frac{w_{n,\mathrm{tot}}+w_{c,\mathrm{tot}}}{R·w_{n,\mathrm{tot}}+w_{c,\mathrm{tot}}}\hfill \end{array}$$(7a)

$$\begin{array}{cc}\hfill I_c& =\frac{w_{n,\mathrm{tot}}+w_{c,\mathrm{tot}}}{R·w_{n,\mathrm{tot}}+w_{c,\mathrm{tot}}}\hfill \end{array}$$(7b)

$$\begin{array}{cc}\hfill R& =\frac{I_n}{I_c}·\hfill \end{array}$$(7c)

The subscript tot denotes the sum of all particle weights. The importance ratio R represents a degree of freedom in this importance-sampling scheme. In the legacy version of TRIPOLI-4, R is a user-defined parameter. When a list of importance ratios is specified, along with the corresponding time values at which changes must be applied, the algorithm adjusts the weights at the beginning of each time step by satisfying two conditions: (1) the equivalence of the sum of the physical weights from time step i to the subsequent time step i + 1 and (2) the preservation of the physical weights for both neutrons and precursors, namely

$$\begin{array}{cc}\hfill w_{n,i}+w_{c,i}& =w_{n,i+1}+w_{c,i+1}\hfill \end{array}$$(8a)

$$\begin{array}{cc}\hfill \frac{{\tilde{w}}_{n,i+1}}{I_{n,i+1}}& =\frac{{\tilde{w}}_{n,i}}{I_{n,i}}\hfill \end{array}$$(8b)

$$\begin{array}{cc}\hfill \frac{{\tilde{w}}_{c,i+1}}{I_{c,i+1}}& =\frac{{\tilde{w}}_{c,i}}{I_{c,i}}·\hfill \end{array}$$(8c)

In the legacy algorithm, the transformation of physical weights into simulation weights are carried out at the end of the inactive cycles of the power iteration, for the particles that have been sampled during the last iteration of the critical phase and represent the initial condition for the kinetic simulation [3].

During the simulation of the time steps, each precursor created by a neutron is assigned a weight equal to the neutron simulation weight divided by the importance ratio, namely ${\tilde{w}}_{c,i}={\tilde{w}}_{n,i}·I_c/I_n$. Conversely, the weight of a precursor decaying into a neutron is multiplied by the importance ratio, namely ${\tilde{w}}_{n,i}={\tilde{w}}_{c,i}·I_n/I_c$. If the importance ratio changes from one step to the next, the weights of neutrons and precursors are adjusted according to equations (8b) and (8c), respectively. Finally, during the scoring phase, any tally involving neutron weights is divided by the neutron importance to ensure unbiased results. A graphical representation of this scheme, considering two time steps for which two different importance ratios are defined, is shown in Figure 2. If only one scalar value is provided for R across the whole time grid, the conversion from the physical weights to the simulation weights (as per Eq. (5)) is performed only once, at the beginning of the kinetic simulation, based on the particle weights inherited from the critical source (Eqs. (7) and (7b)). Selecting a global value for R for the entire transient might be a poor choice, especially for cases evolving over long time-scales and/or with large reactivity excursions, where the weights of the particles can change significantly. As a result, this approach is not particularly effective to deal with such scenarios.

Fig. 2. Schematic representation of the importance sampling strategy implemented as a variance reduction technique in the kinetic algorithm in TRIPOLI-4. The scheme is applied to a kinetic simulation with two time steps, considering two different values Ri and Ri + 1 for the importance ratio, one for each time step.

Despite its apparent simplicity and its easy application, the legacy implementation of the importance-sampling scheme suffers from two major drawbacks.

First, a comprehensive investigation has shown that there exists an optimal value of the importance ratio R, which maximises the simulation performance. However, this value is system-dependent and varies with the reactivity introduced into the system and with the simulation time steps used for applying the population control algorithm [10, 20]. A non-homogeneous time mesh with varying time step sizes might be extremely convenient in many practical problems. Indeed, in the simulation of a critical system a very fine grid for the population control is not needed, since the population is expected to remain stable. On the contrary, when the reactivity excursion is very fast (e.g. during a prompt jump or drop) one might be interested in controlling the population more often; in such case, finer grids might be necessary. Thus, determining the optimal value of R for each time step, while also accounting for reactivity changes, is a challenging task that often relies on a trial-and-error approach [3]. Optimizing algorithm performance would require extensive verification to establish the most suitable set of importance ratios across the entire simulation domain.

Second, more subtly, the population importance values I_n and I_c used in the algorithm are stochastic quantities. In equation (5), the weights are multiplied by a stochastic value, itself depending on the particle weights; thus, equation (5) is a non-linear operation, and has been pointed out to possibly lead to a bias [22], similarly to what happens during the normalization procedure used in the regular power iteration [29, 30]. In particular, the expected value $\mathbb{E}[f({\tilde{w}}_n)]$ of any estimator involving a nonlinear function f will be different from the true value $f(\mathbb{E}[{\tilde{w}}_n])$.

3. A refined adjoint-based importance-sampling scheme

In order to address the aforementioned issues, we have implemented an improved version of the importance-sampling scheme. The conversion from physical weights to simulation weights is performed according to equation (5), as in the legacy version. The major change in the algorithm concerns the elements appearing in the definition of the importance functions, shown in equation (7), which are required for the weight conversion. Instead of estimating these quantities based on the neutron and precursor weights obtained along the Monte Carlo simulation, we rely on the system of deterministic one-family point-kinetics equations

$$\begin{array}{cc}& \frac{dn(t)}{\mathit{dt}}=\frac{{\rho }_k(t)-{\beta }_{\text{eff}}}{{\mathrm{\Lambda }}_{\text{eff}}}·n(t)+\overline{\lambda }·c(t)\hfill \end{array}$$(9a)

$$\begin{array}{cc}& \frac{dc(t)}{\mathit{dt}}=\frac{{\beta }_{\text{eff}}}{{\mathrm{\Lambda }}_{\text{eff}}}·n(t)-\overline{\lambda }·c(t),\hfill \end{array}$$(9b)

where Λ_eff is the effective mean generation time, β_eff the effective delayed neutron fraction, $\overline{\lambda }$ the average precursor decay time and ρ_k the static reactivity. The parameters Λ_eff, β_eff, and $\overline{\lambda }$, which are required in equations (9), are estimated by TRIPOLI-4 through an eigenvalue calculation using the IFP method [31]. This calculation must be performed prior to the kinetic simulation. It is implicitly assumed that these values are only weakly affected by the system evolution during the transient, so that it is safe to consider them to be time-independent. On the contrary, the system reactivity is estimated using the definition of static reactivity ρ_k = 1/k_initial – 1/k_perturbed, running a series of eigenvalue calculations for each time-step configuration where the reactivity changes. It is worth noting that reactivity values can be obtained from power iteration calculations, which are much less computationally expensive than the target kinetic simulation¹. The analytical solution of the point-kinetics equations (9) is provided in the Appendix. The decision to collapse all precursor families into a single effective family is driven by the need of a tractable analytical solution. The investigation of the impact of considering the full system of equations with multiple precursor families is left for future work.

Based on the solution of equations (9), we can determine the importance values for each time step, namely

$$\begin{array}{cc}& I_n(t_i)=R(t_i)·\frac{n(t_i)+c(t_i)}{R(t_i)·n(t_i)+c(t_i)}\hfill \end{array}$$(10a)

$$\begin{array}{cc}& I_c(t_i)=\frac{n(t_i)+c(t_i)}{R(t_i)·n(t_i)+c(t_i)}·\hfill \end{array}$$(10b)

These expressions are equivalent to equations (7a) and (7b), the only difference being that the physical weights of the particles are replaced by the deterministic estimates n(t_i) and c(t_i), which are evaluated at the end of each time step. While equations (9) provide a highly simplified and approximate representation of the kinetic response of the system, we will later show that they are sufficient in this context to capture the overall time evolution and contribute to an effective correction of particles weights.

In the improved version of the importance-sampling scheme, the importance values obtained from equations (10) are non-stochastic, which suppresses the emergence of the aforementioned bias due to the non-linear operation involving statistical weight multiplications. This is a remarkable advantage over the legacy implementation. Furthermore, the proposed approach naturally takes into account the evolution of the importance values over time, which allows for a built-in dynamic adjustment of physical weights to simulation weights, a feature which is not achievable with the legacy algorithm, since the definition of the importance values given in equations (7a) and (7b) remains constant over time and changes only if the user provides a different value of R.

Equations (10) depend on the importance ratio R(t_i), which should be assigned for each time step. In the improved version of the importance-sampling scheme, we have implemented a feature that enables to estimate optimal values of R(t_i) automatically, in order to alleviate the difficulties of the legacy algorithm mentioned in Section 2.3. The key idea is to leverage the equivalence between the importance function and the solution of the adjoint time-dependent transport equations [32]: the adjoint neutron density n^† and the adjoint precursor density c^† physically represent the importance of a neutron (respectively, precursor) with respect to a detector. We define the detector response as the integral over the time interval [t_i, t_i + 1] of the product between the adjoint neutron source Q^†, interpreted as the detector response function, and the neutron density n. Here the adjoint source Q^† is equal to 1 between t_i and t_i + 1, and zero elsewhere. It is assumed that precursors do not directly contribute to the detector response. The importance ratio can be thus derived by taking the ratio between the solutions of the adjoint transport equations, i.e., R = n^†/c^†. By analogy with the strategy used to infer the forward neutron and precursor densities n and c, we assume that the quantities n^† and c^† are determined by solving the set of adjoint point-kinetics equations, namely

$$\begin{array}{cc}& -\frac{d{n}^{†}(t)}{\mathit{dt}}=\frac{{\rho }_k(t)-{\beta }_{\text{eff}}}{{\mathrm{\Lambda }}_{\text{eff}}}{n}^{†}(t)+\frac{{\beta }_{\text{eff}}}{{\mathrm{\Lambda }}_{\text{eff}}}{c}^{†}(t)+Q^{†}(t)\hfill \end{array}$$(11a)

$$\begin{array}{cc}& -\frac{d{c}^{†}(t)}{\mathit{dt}}=\overline{\lambda }{n}^{†}(t)-\overline{\lambda }{c}^{†}(t).\hfill \end{array}$$(11b)

The exact solutions of equations (11) are recalled in the Appendix; for a thorough discussion on the role of the adjoint flux in zero-variance schemes, see e.g. Ref. [11]. The evolution of n^† and c^† with respect to time depends on the physical features of the simulated system (i.e., Λ_eff, β_eff and $\overline{\lambda }$), on the reactivity ρ_k, and on the simulation parameters (i.e., the time step Δt for applying the population control strategy). The solution of the system of equations (11) provides the values n^†(t_i) and c^†(t_i), whose ratio at each time step i serves as an estimate of the importance ratio R(t_i)=n^†(t_i)/c^†(t_i).

At the initial time, assuming a critical reactor condition, the precursor-to-neutron population ratio is given by ${\beta }_{\text{eff}}/({\mathrm{\Lambda }}_{\text{eff}}\overline{\lambda })$, which is a very large number for typical thermal and fast systems. Correspondingly, the ratio between the neutron and precursor importance is expected to be very large, reflecting the imbalance in the respective particle contributions to the response, i.e., the integral of the neutron density n over [t_i, t_i + 1]. When a reactivity change is induced in the system, the importance of neutrons and precursors adjusts accordingly, depending on the sign of the reactivity change. For a positive insertion, the net population of both neutrons and precursors increases, reducing the relative importance of a single particle in contributing to the response, since the overall larger population makes the individual contribution less meaningful. Conversely, for a negative reactivity insertion, the population decreases, and the importance of individual particles increases as their contributions will have a larger impact. Furthermore, if the response involves the simulation of very small time steps, the importance of a neutron will be significantly higher than in simulations with larger time steps.

The time dependence of the importance ratio R(t_i) computed according to this scheme ultimately reduces to a dependence on the size of the time step [t_i, t_i + 1]. Indeed, as easily inferred from equation (6) in the Appendix, the solutions of equations (11) are evaluated at t = t_i, meaning that they depend only on the time-step size t_i + 1 − t_i. In the improved version of the importance-sampling scheme, it is still possible to override the importance ratio determined from solving the adjoint point-kinetics equations (11) and impose instead a constant user-provided value. However, even when a constant value of R is used, the new method differs from the legacy approach, as the importance values in the new algorithm depend on the neutron and precursor concentrations computed using the forward point-kinetics equations (10), whereas in the legacy algorithm they are based on the weights given in equation (7). Therefore, the improved importance-sampling is better tailored to follow the kinetic evolution of the simulated scenario and is expected to perform better than the legacy algorithm.

In conclusion, the main features of the new importance-sampling scheme calculation are:

• (i)

  The definition of the importance values depends on the importance ratio and on estimates coming from the solution of the point-kinetics equations. This dependence is particularly advantageous as it enables a step-by-step tracking of the system evolution. Consequently, it provides a finely tailored update of the importance values, taking into account the changing sizes of the neutron and precursor populations;

• (ii)

  The user is not required to carry out a trial-and-error procedure to determine the best importance ratio for the system: a set of values is automatically determined based on the solution of the adjoint point-kinetics equations;

• (iii)

  The importance ratio depends on specific parameters of the system and varies according to changes in the time-step size or in the reactivity value.

4. Testing the new kinetic capabilities of TRIPOLI-4

When comparing the effectiveness of various Monte Carlo simulations, a commonly used metric is the figure of merit (FoM), namely

$$\begin{array}{c}\hfill \text{FoM}=\frac{1}{{\sigma }^{2}T},\end{array}$$(12)

where σ² is the estimated variance for a given observable and T is the duration of the simulation. In this section we investigate the population control strategy proposed in the previous section using the FoM as a comparative metric. Specifically, we evaluate the performance of the two strategies described in Section 2.2, which are combing and Russian roulette on precursors, when combined with the importance-sampling scheme, exploiting different values of the importance ratio. All FoM calculations presented hereafter are conducted using alternatively the legacy and the refined algorithm and they are executed as single-processor simulations, in order to avoid potential distortions associated with parallelism. We emphasize that our analysis is strictly limited to configurations that are within reach of the legacy algorithm; otherwise, a meaningful comparison would not be possible. All the simulation results presented in this section are based on the ENDF/B-VIII.0 nuclear data library [33].

4.1. A critical system: Flattop-Pu

We start by considering the PU-MET-FAST-006 (Flattop-Pu) benchmark [23], which has been previously used for the assessment of the preliminary implementation of the revised importance sampling algorithm in TRIPOLI-4 [22]. This benchmark consists of a plutonium sphere surrounded by a natural uranium reflector, and is characterized by a fast spectrum. For illustration, see Figure 3; the kinetics parameters of Flattop-Pu are provided in Table 1. Despite its apparent simplicity, the Flattop-Pu benchmark provides an extremely challenging scenario for kinetics [3], because of the very low value of the mean generation time, Λ_eff ∼ 10 ns, which induces a huge imbalance between the neutron and precursor population sizes at equilibrium.

Fig. 3. Graphical representation of the Flattop-Pu benchmark in its critical configuration.

Our first FoM comparison is based on the following reference simulation: the kinetic evolution is performed over a single time step, during which the system is maintained in its critical state². The observable chosen for the FoM analysis is the ensemble-averaged time integral of the energy and volume-integrated neutron flux. We run this simulation using the legacy version of TRIPOLI-4 and the refined algorithm. To evaluate the optimal population control strategy, we investigate the following schemes:

• (i)

  Russian roulette combined with the importance-sampling scheme, exploiting one by one six different values of importance ratio, namely R = {1, 10³, 10⁶, 10⁹, 10¹², 10¹⁵}. For the new algorithm, we test also the value of R computed solving the adjoint point-kinetics equations, as described in Section 3.

• (ii)

  Combing combined with the importance-sampling scheme, exploiting one by one the same six values of importance ratio as listed in (i) and, for the new algorithm, the values computed solving the adjoint point-kinetics equations.

All these simulations are repeated for three different sizes of the simulated time step, namely Δt = {10⁻⁸, 10⁻⁹, 10⁻¹⁰} s. The values of Δt have been selected based on the value of Λ_eff. Additionally, since this analysis encompasses simulations over a wide range of importance ratio values, using larger time steps would be impractical due to the increased computation time.

For each kinetic simulation, the number of particles and cycles is adjusted to ensure that the estimated standard error on the total neutron flux remains below 20%. When using lower importance ratio values, simulations require a larger number of histories to achieve the target accuracy, due to the very large imbalance between neutrons and precursors. In this kinetic scenario, involving only a single time step over which the system is critical, the differences between the legacy and the new algorithm are not expected to be significant for the same value of importance ratio. This occurs because here the only discrepancy between the importance-sampling schemes lies in how the importance values are computed at the initial time to obtain the simulation weights for the single time step: the legacy algorithm will use weight estimates from the power iteration, whereas the improved algorithm will use the point-kinetics equations. Consequently, this investigation of a critical system cannot fully showcase the improvements of the revised importance-sampling scheme in dynamically adjusting the importance values according to the system evolution; nonetheless, the results illustrated in the following clearly demonstrate the impact of the chosen population control strategy on the FoM and confirm the existence of an optimal range of values for the importance ratio.

In Figure 4 we show the FoM gain associated with different population control strategies and importance ratios used in the importance-sampling scheme, for three different sizes of the time step. For each time step, the FoM gain is always measured relative to a reference simulation in which Russian roulette is employed, and the importance-sampling scheme is not applied (R = 1). Notably, the FoM gain shows a strong dependence on the selected value of importance ratio R, whereas no significant differences are observed between using combing over Russian roulette. The observed behavior of the FoM gain as a function of the importance ratio, characterized by a steep initial increase followed by a plateau, is consistent across all tested time step sizes and reflects the fact that the simulation is sub-optimal for low values of R. This is due to the disparity in the weight values between precursors and neutrons, leading to an under-sampling of neutrons. Precursors decay within the time step based on the forced decay scheme described in Section 2.1: without an appropriate importance adjustment, they will produce delayed neutrons with extremely low weights. Most of these neutrons are subsequently eliminated through the Russian roulette that is enforced after each collision in TRIPOLI-4. Therefore, choosing an importance ratio that increases the weight of neutrons are crucial and clearly improves the performance of the simulation by various orders of magnitude. The FoM gain associated with the importance ratio relying on the adjoint point-kinetics equations are located within the plateau region, indicating that the computed importance ratio provides a optimal FoM gain.

Fig. 4. Analysis of the FoM gain associated to different population control strategies and choices of the population importance ratio for the kinetic simulation of a single time step of the Flattop-Pu benchmark in its critical configuration. The comparison is repeated for three different choices of the time step size: Δt = 10−8 s, Δt = 10−9 s and Δt = 10−10 s. For each time step size all the computed FoM values are normalized by the FoM value of a reference simulation where Russian roulette is employed and no importance-sampling scheme is applied (R = 1). The results of the legacy algorithm are displayed in red, while those of the new algorithm are shown in black. The circle in the curve of the new algorithm represents the value of FoM gain associated with the importance ratio based on the solution of the adjoint point-kinetics equations.

Previous studies have shown that the FoM gain also depends on the time step size [10, 20]. Specifically, it has been observed that shorter time steps allow for a higher maximum gain achievable through the importance-sampling scheme, for both combing and Russian roulette. However, in the simulations conducted in this work, the observed differences in FoM gain are relatively small, suggesting that the previously identified trend is not strongly pronounced for our choices of time step sizes. Overall, as briefly outlined above, the computational efficiency trends of the legacy and improved implementations are rather consistent; the importance-sampling scheme has a significant impact on improving the FoM, and the value of the importance ratio derived from the solution of the adjoint point-kinetics equations is highly effective.

4.2. Analysis of a reactivity insertion in the Flattop-Pu benchmark

The next step of our analysis addresses more challenging scenarios for kinetic simulations, where the system transitions from a critical state to a supercritical or sub-critical state. For this purpose, we consider a positive reactivity insertion into the initially critical Flattop-Pu benchmark, inspired by previous investigations in Ref. [22]. After 0.5 μs we apply a step-wise modification of the volume of the inner plutonium sphere, from an unperturbed radius R₁ = 4.5332 cm to a perturbed radius R₁ = 4.5650 cm, by keeping its density fixed, which results in a reactivity insertion of approximately $2.2. The radius of the external sphere remains constant. The inner sphere is then brought back to its original size after 2.5 μs, which restores criticality. Based on estimates derived from point kinetics, using the parameters shown in Table 1, we find that this transient causes the flux to increase to approximately 2.7 times its initial value.

We assess the FoM across six distinct time intervals of different length. In the first step ([0 − 0.5] μs), the system remains critical. During the second ([0.5 − 2] μs), third ([2 − 2.75] μs), and fourth ([2.75 − 3] μs) steps, the system is kept in the supercritical state described above. Finally, in the fifth ([3 − 3.25] μs) and last step ([3.25 − 4] μs), the critical configuration is restored. Time steps of different size have been chosen in order to increase the frequency of population control during the supercritical excursion. For the performance assessment, we estimate the FoM(t_i) value for each time step t_i, according to equation (12). We rely on the same observable as in Section 4.1 (i.e. the time, energy and volume-integrated neutron flux) tallied for the specific time step. Regarding the simulation time, we consider the cumulative time required to perform the dynamic simulation up to and including the time step for which the FoM is assessed. The simulations are performed alternatively with both importance-sampling schemes, combined with combing and Russian roulette. The importance ratio is either computed by solving the adjoint point-kinetics equations, or imposed by choosing values in R = {1, 10³, 10⁶, 10⁹, 10¹²}.

The point-kinetics evolution of the kinetic scenario is illustrated in Figure 5, using the parameters provided in Table 1. We also display the evolution of the importance ratio computed using the improved importance-sampling scheme described in Section 3, along with the evolution of the neutron and precursor importance. The importance ratio depending on the adjoint point-kinetics equations changes whenever there is a variation either in the time step size or the static reactivity value, which leads to an adaptive time-dependent control of the neutron population through the factor I_n. In Figure 6 we present the evolution of the FoM gain as a function of the importance ratio across the six time steps analyzed, combined with either combing (Fig. 6i) or Russian roulette (Fig. 6ii). The improved algorithm is compared to the legacy implementation. For each of the six time steps, the FoM gain is assessed with respect to a reference simulation using the importance ratio R = 10³, combined with either combing (Fig. 6i) or Russian roulette (Fig. 6ii). The case of R = 1 (i.e. no importance-sampling) is not used for our comparisons, since numerical tests have shown that for this configuration it is extremely difficult to achieve statistical uncertainties below 20% for all time steps, due to the strong under-sampling of the neutron population.

This analysis suggests that variance-reduction techniques are not only beneficial for improving statistics but are in practice mandatory for obtaining meaningful average results in the first place. When combing is employed, the trend across the time steps shows a decrease in the maximum achievable FoM gain. In contrast, with Russian roulette, the FoM gain remains relatively constant across time steps, except for the very first one, where it reaches a slightly higher value. In both cases, the use of the importance-sampling scheme is crucial for optimizing the simulation performance. In general, after a sharp increase, a plateau of the FoM is observed, extending over a wide range of values of importance ratio; beyond the plateau, for high values of the importance ratio, a steep decrease occurs. This behaviour is likely due to an extreme increase of the neutron weight, causing the sampled neutrons to undergo excessive splitting, which in turn leads to strong correlations that reduce the FoM [20, 28]. Overall, for this benchmark configuration the performance of the improved algorithm is comparable to that of the legacy implementation, except for the first time step, where the improved algorithm performs better in the plateau region. The FoM gain associated with the importance ratio value obtained from the adjoint point-kinetics equations falls within the plateau region, confirming also in this case its optimality.

Fig. 5. Left: evolution of the neutron flux in the Flattop-Pu benchmark after a positive reactivity insertion of $2.2, according to the solution of the point-kinetics equations with one family of precursors. Middle: Evolution of the importance ratio for every time-step used for the kinetic simulation. The population control strategy is performed on a grid that becomes finer as the neutron flux increases. The values of In and Ic are calculated according to equations (6); consequently, they vary with the simulated time-step length and with the considered reactivity. Right: Evolution of the importance values (In and Ic) for the same simulation.

Fig. 6. Evolution of the FoM gain for every time step simulated in the kinetic scenario presented in Figure 5. The displayed FoM gain values for every time step are computed choosing combing (top) or Russian roulette on precursors (bottom) as the population control strategy. In both cases we investigate a broad range of importance ratio values, spanning from 103 to 1012. The circle represents the FoM gain obtained using the value of R determined from the solution of the adjoint point-kinetics equations.

4.3. Analysis of the SPERT III E-core benchmark

In order to further broaden our analysis, we consider a scenario corresponding to a negative reactivity insertion, using the Special Power Excursion Reactor Test III (SPERT III) E-core benchmark [34–36], based on a small pressurized-water research reactor built and operated in the United States in the 1960s to probe the effects of power excursions on fuel and structural materials. The core contains 60 assemblies, including 48 fuel assemblies with 25 (5 × 5) fuel rods, 4 assemblies with 16 (4 × 4) fuel rods, and 4 pairs of borated steel control rods. The SPERT III E-core benchmark was previously used in a series of investigations carried out with TRIPOLI-4, to which we refer for a thorough description of the reactor components [4, 24, 37–39]. The corresponding model is illustrated in Figure 7.

Fig. 7. Radial and axial cut of the SPERT III E-core benchmark, as modelled with TRIPOLI-4.

Here we will consider a shut-down scenario: starting from a critical (cold zero power) configuration, all the control rods are inserted into the core, leading to an instantaneous reactivity withdrawal of ∼$19. Following the control rod insertion, the neutron population first undergoes a prompt drop and then a slower decrease driven by the precursor decay. In view of the position of the control rods in the core, the reactivity insertion is radially symmetric. This kind of kinetic simulation has been previously investigated using the legacy kinetic algorithm in TRIPOLI-4 and a constant importance ratio [3], and it is thus interesting to revisit it using the improved algorithm, with a customized time grid. Specifically, the time grid is divided as follows:

$$\begin{array}{c}\hfill \mathrm{\Delta }t=\{\begin{array}{cc}100\mathrm{\mu }\mathrm{s},\hfill & \text{for}0\mathrm{s}\le t<1\mathrm{ms},\hfill \\ 1\mathrm{ms},\hfill & \text{for}1\mathrm{ms}\le t<0.011\mathrm{s},\hfill \\ 10\mathrm{ms},\hfill & \text{for}0.011\mathrm{s}\le t<0.111\mathrm{s},\hfill \\ 0.1\mathrm{s},\hfill & \text{for}0.111\mathrm{s}\le t<1.111\mathrm{s},\hfill \\ 1\mathrm{s},\hfill & \text{for}1.111\mathrm{s}\le t<10.111\mathrm{s}.\hfill \end{array}\end{array}$$(13)

This choice is motivated by the fact we want to ensure a finer time grid during the early phases of the simulation, allowing frequent application of population control to maintain a roughly constant and convenient number of particles during the rapid transient. At later times, the grid becomes coarser as the decrease of the flux is regulated by slower time-scale of the delayed neutrons precursors.

To test the simulation capabilities of TRIPOLI-4 on the transient described above, we rely on the importance-sampling scheme with the adaptive calculation of R implemented in the improved algorithm. In this scenario, the improved algorithm automatically adapts the importance ratio according to the solution of the adjoint point-kinetics equations, reflecting the fact that in the first part of the grid (with shorter time steps) the neutrons are more important than the precursors, while in the second part (with longer time steps) the contribution from the precursors is the most relevant.

For comparison, the legacy algorithm offers two approaches for handling the two time-scales of this scenario. We can either use a single constant value for the importance ratio R, or assign a list of importance values for each time step. The use of a constant value is immediately faced with a serious issue: choosing a large value of R to boost the number of neutrons in the first phase of the transient (where time steps are very short and the contribution from precursor decay is not significant) will allocate excessive computational resources to neutrons in the second phase, potentially leading to a catastrophic memory burden. Conversely, if a low value of R is chosen in view of optimizing the second phase of the transient, then in the first phase neutrons will be significantly under-sampled, leading to unacceptable statistical uncertainty. In order to overcome the problems posed by using a constant R for the simulation of a long transient in SPERT III, in Ref. [3] two different values of importance ratio were determined through trial and error, one to deal with the prompt drop and one for the rest of the transient. For reference, we will reproduce this configuration in our investigation, setting R = 10⁶ from t = 0 s to t = 0.001 s, and R = 10² thereafter.

The values of the importance ratio computed with the improved algorithm for this scenario are shown in Figure 8, as a function of all the time steps sizes listed in equation (13). The guess values used in the legacy algorithm are also displayed. It is apparent that the value of R changes by several orders of magnitude over the transient, significantly more than what was observed for the positive transient described in Section 4.2, where the ratio varies only by a factor of two. This suggests that the adaptive calculation of the importance ratio performed by the improved algorithm might be particularly efficient for this scenario.

Fig. 8. Importance ratio as a function of the time-step size used in the dynamic simulation. In the new algorithm these values are computed by the adjoint point-kinetics scheme, whereas in the legacy algorithm they are provided by the user.

The comparison between the improved and the legacy algorithm is shown in Figure 9, using the neutron flux as the observable and the FoM(t_i) values, computed at each time step t_i, as metrics. We consider two cases: combing (Fig. 9i) and Russian roulette (Fig. 9ii). In both cases simulations are performed using 100 cycles and 500 neutrons per cycle.

Fig. 9. A comparison of the neutron flux (left) and the standard error of the neutron flux (middle) is presented between the legacy algorithm, using two different values of R, and the new algorithm, where R is adapted at each change in time step size according to the adjoint-based scheme. The figure on the right displays the FoM ratio between the new and legacy algorithms. This analysis is provided for both Russian roulette (top) and combing (bottom).

Our results confirm that the adaptive procedure for the importance ratio calculation provides the most efficient strategy for the FoM throughout the whole simulation. This trend is observed regardless of the choice of the population control algorithm. However, while the improved algorithm demonstrates comparable results and uncertainties below 5% in the flux score when using either combing or Russian roulette, the legacy algorithm appears to yield more accurate results when combing is applied, with a maximum uncertainty of approximately 7% compared to 22%. In conclusion, the improved algorithm proves to be a highly effective variance-reduction technique, particularly when combined with the adjoint-based calculation of the importance ratio.

4.4. Handling strong reactivity excursions

After evaluating the FoM achieved by the improved algorithm in TRIPOLI-4, our goal is to further demonstrate its practical capabilities. In particular, we would like to assess its impact on the simulation of strong flux excursions and long transients, which are two of the major challenges for kinetic Monte Carlo simulations. For this purpose, we reconsider the Flattop-Pu benchmark addressed in Section 4.1. Reactivity is again inserted by varying the inner Pu sphere radius. After maintaining a critical state for 1 μs, we inject $2.2: in four distinct scenarios, the system is kept in supercritical conditions for (i) 24 μs, (ii) 34 μs, (iii) 44 μs, and (iv) 51 μs, respectively. Adjusting the time duration over which the system remains supercritical, we can generate neutron flux excursions of different amplitudes, spanning several orders of magnitude. According to the point-kinetics equations, the neutron population is expected to increase by a factor of the order of 10³ in the first case, 10⁴ in the second case, 10⁵ in the third case, and 10⁶ in the fourth case. At the end of each reactivity insertion, the system is then brought back to criticality, and the neutron population is tracked over a total simulation time of 100 μs.

In the first, second and third cases, the time step size used for the application of the population control algorithm is kept equal to Δt = 1/3 μs, while in the last case a finer grid with Δt = 0.25 μs is employed, given the extremely steep increase of the flux during the supercritical regime. All cases were simulated with 1000 cycles and 10000 neutrons per cycle. Our simulation results are shown in Figure 10, where the flux obtained with the TRIPOLI-4 kinetic simulation is compared to the solution of the point-kinetics equations, for reference. The upper and lower bounds of the error band on the point-kinetics curve correspond to solutions obtained with ρ_k = 1/k_critical − 1/k_perturbed varied within its uncertainty, which in turn derives from the statistical uncertainty of the effective multiplication factor estimated through the power iteration algorithm in TRIPOLI-4.

Fig. 10. Neutron flux obtained on the Flattop-Pu benchmark with the new kinetic algorithm, juxtaposed to the evolution predicted by the point-kinetics equations. The transient is initiated by enlarging the internal Pu sphere, and it consists of a positive reactivity insertion of $2.2 on a system initially kept in its critical state for 1 μs. Different durations of the super-critical state are investigated, namely (i) 24 μs, (ii) 34 μs, (iii) 44 μs, (iv) 51 μs. The evolution of the importance values and the importance ratio used in the new importance-sampling scheme is also presented, demonstrating the highly adaptive and dynamic nature of this new variance reduction strategy.

In the same figure, we also present the evolution of the importance ratio and the neutron and precursor importance values, normalized to their initial values. In all the tested cases, the importance ratio remains relatively stable, as there are no significant changes in the time step size or reactivity. In contrast, the evolution of the neutron and precursor importance values spans several orders of magnitude, following the duration of the supercritical state. The significant variation in the importance values, despite the importance ratio remaining nearly constant, arises from the fact that these values depend on the neutron and weight proportions, which in the improved algorithm are continuously updated based on the direct point-kinetics equations (10a) and (10b). This is in contrast to the legacy algorithm, where these proportions are exclusively based on the initial weights (Eqs. 7a and 7b) and are not updated during the simulation. The computational cost of these simulations, expressed in kh_CPU³, is presented in Table 2.

The newly developed variance-reduction technique is ideally suited for this challenging scenario, as witnessed by the evolution of the importance values to closely follow the transient. To illustrate the impact of the features added in the new algorithm, Figure 11 presents the results obtained using the legacy kinetic algorithm for one of the cases just investigated, specifically the scenario where the system remains in a supercritical state for 34 μs. In this configuration, the importance ratio and importance values are not dynamically adjusted during the simulation but remain constant, resulting in a less effective variance reduction and a significant underestimation of the neutron flux.

Fig. 11. Neutron flux obtained on the Flattop-Pu benchmark with the legacy kinetic algorithm, juxtaposed to the evolution predicted by the point-kinetics equations. The transient is initiated by enlarging the internal Pu sphere, and it consists of the maintenance of a positive reactivity insertion of $2.2 for 34 μs on a system initially kept in its critical state for 1 μs. The importance values and importance ratio used in the legacy importance-sampling scheme in this case remain constant, as the legacy algorithm does not support their dynamic adaptation. As a consequence of this less tailored variance reduction technique, the kinetic results from TRIPOLI-4 significantly underestimate the system’s evolution.

4.5. Handling long transients

To further substantiate our analysis and probe the behaviour of the improved algorithm with respect to scenarios involving long transients, we consider finally the CROCUS benchmark [25], based on a thermal two-region uranium-fueled and light-water-moderated zero-power reactor located at the École Polytechnique Fédérale de Lausanne (EPFL) in Switzerland. The CROCUS core is loaded with two types of fuel rods, using respectively uranium oxide (UO₂) at 1.806%_wt²³⁵U enrichment, and uranium metal at 0.947%_wt²³⁵U enrichment. The experimental configurations that we have used for our Monte Carlo simulations are documented in the international OECD/NEA benchmark collection [25], and have been previously modelled with the TRIPOLI-4 code [40, 41] (see Fig. 12).

Fig. 12. Radial cut of the TRIPOLI-4 model of the CROCUS core.

In this work we consider three transient scenarios, where the system departs from a critical configuration (H₁) and attains three different supercritical configurations, called H₂, H₃, and H₄, respectively, by varying the water level inside the core tank. Therefore, the reactivity insertion for these scenarios are again radially symmetric. The water levels, along with the corresponding experimentally measured inverse reactor period α, are shown in Table 3. In our simulations, the system is initially critical for 10 s, then suddenly becomes supercritical by changing the water level and stays supercritical for 1 minute; finally, criticality is restored and kept for extra 10 s. We used 3000 neutrons per cycle for a total of 150 cycles. Population control is performed on a time grid with Δt = 0.5 s.

4.6. Comparison between the results obtained with the kinetic algorithm and point kinetics

The neutron fluxes estimated by TRIPOLI-4 are shown in Figure 14 and compared to the solution of the point-kinetics equations, for reference. As in Section 4.4, the uncertainty on the point kinetics curve is derived from the uncertainty of the effective multiplication factors used to compute the inserted static reactivity. Exclusively for this benchmark study, the JEFF-3.3 nuclear data library [42] has been tested alongside ENDF/B-VIII.0 to evaluate the impact of nuclear data on the results. Overall, a good agreement is observed between the results of the kinetics simulations in TRIPOLI-4 and point-kinetics for the configurations H₃ and H₄, for both libraries. This is not surprising, given that the CROCUS reactor is a highly compact system where spatial effects are minimal, particularly when reactivity is introduced homogeneously via water level control in the tank. As a general trend, the flux estimated using the ENDF/B-VIII.0 nuclear data library is higher than that estimated using JEFF-3.3, consistent with the fact that the injected reactivity is estimated to be larger when using the former library. This trend is also reflected in the time required to complete the simulation, measured in kh_CPU, for the three H₂, H₃, H₄ supercritical states, which are consistently higher for simulations using ENDF/B-VIII.0, as shown in Table 4.

Fig. 13. Evolution of the CROCUS system from a critical state (H1) to a super critical configuration (H2), which is kept for one minute before going back to criticality. The results of the kinetic simulation with TRIPOLI-4 are compared with the behavior predicted by the point-kinetics equation, using the static reactivity ρk = 1/kcritical − 1/kperturbed and the k-adjoint weighted kinetics parameters (green) but also the dynamic reactivity ρα and the α-adjoint weighted kinetics parameters (red).

Fig. 14. Transient simulation, going from the critical configuration called H1 to different supercritical states (H2, H3, H4) by adjusting the water level in the tank. The results of the kinetic simulation in TRIPOLI-4 are compared with the solution of the point-kinetics equation, using the adjoint-weighted kinetics parameters. Results have been obtained both with the ENDF/B-VIII.0 (left) and the JEFF-3.3 (right) nuclear data library.

The H₂ configuration shows the largest discrepancy between the point-kinetics and TRIPOLI-4 kinetics results, as illustrated in Figure 14i. Indeed, in this case the point-kinetics equations seem to underestimate the flux in the system compared to kinetic simulations. Stimulated by this puzzling behaviour, we resorted to a complementary analysis: instead of relying on the static reactivity ρ_k and adjoint k-weighted kinetics parameters, we solved the point-kinetics equations using the dynamic reactivity ρ_α, defined by the equation

$$\begin{array}{c}\hfill {\rho }_{\alpha }=\alpha ({\mathrm{\Lambda }}_{\alpha }+\sum _{j=1}^6\frac{{\beta }_{\alpha ,j}}{\alpha +{\lambda }_j}),\end{array}$$(14)

where α is the dominant time eigenvalue, and the adjoint α-weighted kinetics parameters Λ_α and β_α, j, which can be computed using TRIPOLI-4 [43–46]. The corresponding point-kinetics results are illustrated in Figure 13 and align more closely with the behavior predicted by the kinetic simulation for both nuclear data libraries. For the H₃ and H₄ configurations, the good agreement already observed with the standard point-kinetics approach is also maintained with the dynamic reactivity ρ_α and the α-weighted kinetics parameters (not shown here). These considerations clearly demand further investigations, since a systematic comparison between k-point kinetics and α-point kinetics could provide valuable insights, especially in view of the possibility of using kinetic simulations as a reference, together with experimental results (when available).

5. Conclusions

In this manuscript we have presented the development of an improved importance-sampling scheme for time-dependent simulations in TRIPOLI-4, in order to effectively tame the imbalance between the neutron and precursor populations. In the first part of this work we have provided a general overview of this scheme, highlighting the differences with respect to the legacy implementation. We have then discussed a thorough analysis of the algorithm over a set of benchmark configurations and various reactivity insertion and withdrawal scenarios, using the FoM as a metric for the performance evaluation. In the last portion of the paper we have showcased the enhanced kinetic capabilities enabled by the new importance-sampling scheme, demonstrating the simulation of longer and more steep transients previously unattainable with the legacy version of the algorithm. Our findings contribute thus to establishing Monte Carlo methods as a reference simulation tool for the analysis of non-stationary configurations in reactor physics applications. A possible shortcoming of our findings is that we only considered scenarios where the cores are tightly coupled and the reactivity perturbation is radially symmetric, which improves the agreement between point kinetics and full three-dimensional kinetics; this might not be the case for more general scenarios.

Future work will concern the validation of the improved time-dependent capabilities of TRIPOLI-4 against experimental data stemming from experimental campaigns that are currently under way in the CROCUS reactor, encompassing spatially symmetric and asymmetric reactivity insertions and withdrawals.

Acknowledgments

TRIPOLI-4^® is a registered trademark of CEA.

Funding

This work was partially funded by EDF.

Conflicts of interest

The authors declare no conflict of interest.

Data availability statement

The datasets generated and/or analysed during the current study are not publicly available due to them being generated using export-controlled codes, but are available from the corresponding author on reasonable request.

Author contribution statement

Cecilia Montecchio implemented the new version of the importance-sampling scheme, performed the simulations and wrote the manuscript. Davide Mancusi and Andrea Zoia supervised the methodologies, verified the algorithm implementation, and wrote the manuscript. Vincent Lamirand and Wilfried Monange reviewed the article and co-supervised the work.

References

Appendix

Exact solutions of the forward and adjoint point-kinetics equations

We report the exact solution of the forward point kinetics equations, considering only one family of precursors:

$$\begin{array}{cc}\hfill n(t)& =n_0(\frac{exp({\alpha }_{1}t)-exp({\alpha }_{2}t)}{{\mathrm{\Lambda }}_{\text{eff}}({\alpha }_1-{\alpha }_2)}\hfill \\ \hfill & +\frac{{\alpha }_{2}exp({\alpha }_{1}t)-{\alpha }_{1}exp({\alpha }_{2}t)}{{\alpha }_2-{\alpha }_1})\hfill \end{array}$$(A.1a)

$$\begin{array}{c}\hfill c(t)=\frac{{\beta }_{\text{eff}}}{\lambda {\mathrm{\Lambda }}_{\text{eff}}}\left(\frac{{\alpha }_{2}exp({\alpha }_{1}t)-{\alpha }_{1}exp({\alpha }_{2}t)}{{\mathrm{\Lambda }}_{\text{eff}}({\alpha }_2-{\alpha }_1)}\right).\end{array}$$(A.1b)

We also provide the exact solution of the adjoint point kinetics equations evaluated at t_i:

$$\begin{array}{cc}\hfill n^{†}(t_i)=v& (\frac{\overline{\lambda }{\alpha }_2(e^{{\alpha }_{1(t_{i+1}-t_i)}}-1)}{{\alpha }_1-{\alpha }_2}+\frac{\overline{\lambda }{\alpha }_1(1-e^{{\alpha }_{2(t_{i+1}-t_i)}})}{{\alpha }_1-{\alpha }_2}\hfill \\ \hfill & +\frac{{\alpha }_1{\alpha }_2(e^{{\alpha }_{1(t_{i+1}-t_i)}}-e^{{\alpha }_{2(t_{i+1}-t_i)}})}{{\alpha }_1{\alpha }_2({\alpha }_1-{\alpha }_2)})H(t_{i+1}-t_i)\hfill \end{array}$$(A.2a)

$$\begin{array}{cc}\hfill c^{†(t_i)}=v\overline{\lambda }& (\frac{{\alpha }_1(1-e^{{\alpha }_{2(t_{i+1}-t_i)}})}{{\alpha }_1{\alpha }_2({\alpha }_1-{\alpha }_2)}+\frac{{\alpha }_2(e^{{\alpha }_{1(t_{i+1}-t_i)}}-1)}{{\alpha }_1{\alpha }_2({\alpha }_1-{\alpha }_2)}),\hfill \end{array}$$(A.2b)

where v is the neutron velocity, which cancels out when taking the ratio n^†(t_i)/c^†(t_i), and H is the Heaviside step function. Here α₁ and α₂ are the roots of the characteristic polynomial associated with the Nordheim equation, namely:

$$\begin{array}{c}\hfill {\alpha }_{1,2}=\frac{({\rho }_k-{\beta }_{\text{eff}}-{\mathrm{\Lambda }}_{\text{eff}}\overline{\lambda })+\sqrt{{({\mathrm{\Lambda }}_{\text{eff}}\overline{\lambda }-{\rho }_k\pm {\beta }_{\text{eff}})}^2+4{\rho }_k\overline{\lambda }{\mathrm{\Lambda }}_{\text{eff}}}}{2{\mathrm{\Lambda }}_{\text{eff}}}\end{array}$$(A.3)

Among all the parameters appearing in equation (A.3), only the reactivity ρ_k changes over the time grid, mirroring the changes in the simulated configuration. All the other parameters are assumed to stay constant during the transient.

All Tables

All Figures

	Fig. 1. Sketch of the current implementation of the kinetic algorithm in TRIPOLI-4, taken from Ref. [3]. In the kinetic simulation (green) neutrons and precursors are sampled during the last iteration of the power iteration (white) and their subsequent time evolution is followed on a time grid, during which forced precursor decay and population control algorithms are applied.
In the text

	Fig. 2. Schematic representation of the importance sampling strategy implemented as a variance reduction technique in the kinetic algorithm in TRIPOLI-4. The scheme is applied to a kinetic simulation with two time steps, considering two different values Ri and Ri + 1 for the importance ratio, one for each time step.
In the text

	Fig. 3. Graphical representation of the Flattop-Pu benchmark in its critical configuration.
In the text

Fig. 4. Analysis of the FoM gain associated to different population control strategies and choices of the population importance ratio for the kinetic simulation of a single time step of the Flattop-Pu benchmark in its critical configuration. The comparison is repeated for three different choices of the time step size: Δt = 10−8 s, Δt = 10−9 s and Δt = 10−10 s. For each time step size all the computed FoM values are normalized by the FoM value of a reference simulation where Russian roulette is employed and no importance-sampling scheme is applied (R = 1). The results of the legacy algorithm are displayed in red, while those of the new algorithm are shown in black. The circle in the curve of the new algorithm represents the value of FoM gain associated with the importance ratio based on the solution of the adjoint point-kinetics equations.

In the text

Fig. 5. Left: evolution of the neutron flux in the Flattop-Pu benchmark after a positive reactivity insertion of $2.2, according to the solution of the point-kinetics equations with one family of precursors. Middle: Evolution of the importance ratio for every time-step used for the kinetic simulation. The population control strategy is performed on a grid that becomes finer as the neutron flux increases. The values of In and Ic are calculated according to equations (6); consequently, they vary with the simulated time-step length and with the considered reactivity. Right: Evolution of the importance values (In and Ic) for the same simulation.

In the text

Fig. 6. Evolution of the FoM gain for every time step simulated in the kinetic scenario presented in Figure 5. The displayed FoM gain values for every time step are computed choosing combing (top) or Russian roulette on precursors (bottom) as the population control strategy. In both cases we investigate a broad range of importance ratio values, spanning from 103 to 1012. The circle represents the FoM gain obtained using the value of R determined from the solution of the adjoint point-kinetics equations.

In the text

	Fig. 7. Radial and axial cut of the SPERT III E-core benchmark, as modelled with TRIPOLI-4.
In the text

	Fig. 8. Importance ratio as a function of the time-step size used in the dynamic simulation. In the new algorithm these values are computed by the adjoint point-kinetics scheme, whereas in the legacy algorithm they are provided by the user.
In the text

Fig. 9. A comparison of the neutron flux (left) and the standard error of the neutron flux (middle) is presented between the legacy algorithm, using two different values of R, and the new algorithm, where R is adapted at each change in time step size according to the adjoint-based scheme. The figure on the right displays the FoM ratio between the new and legacy algorithms. This analysis is provided for both Russian roulette (top) and combing (bottom).

In the text

Fig. 10. Neutron flux obtained on the Flattop-Pu benchmark with the new kinetic algorithm, juxtaposed to the evolution predicted by the point-kinetics equations. The transient is initiated by enlarging the internal Pu sphere, and it consists of a positive reactivity insertion of $2.2 on a system initially kept in its critical state for 1 μs. Different durations of the super-critical state are investigated, namely (i) 24 μs, (ii) 34 μs, (iii) 44 μs, (iv) 51 μs. The evolution of the importance values and the importance ratio used in the new importance-sampling scheme is also presented, demonstrating the highly adaptive and dynamic nature of this new variance reduction strategy.

In the text

Fig. 11. Neutron flux obtained on the Flattop-Pu benchmark with the legacy kinetic algorithm, juxtaposed to the evolution predicted by the point-kinetics equations. The transient is initiated by enlarging the internal Pu sphere, and it consists of the maintenance of a positive reactivity insertion of $2.2 for 34 μs on a system initially kept in its critical state for 1 μs. The importance values and importance ratio used in the legacy importance-sampling scheme in this case remain constant, as the legacy algorithm does not support their dynamic adaptation. As a consequence of this less tailored variance reduction technique, the kinetic results from TRIPOLI-4 significantly underestimate the system’s evolution.

In the text

	Fig. 12. Radial cut of the TRIPOLI-4 model of the CROCUS core.
In the text

Fig. 13. Evolution of the CROCUS system from a critical state (H1) to a super critical configuration (H2), which is kept for one minute before going back to criticality. The results of the kinetic simulation with TRIPOLI-4 are compared with the behavior predicted by the point-kinetics equation, using the static reactivity ρk = 1/kcritical − 1/kperturbed and the k-adjoint weighted kinetics parameters (green) but also the dynamic reactivity ρα and the α-adjoint weighted kinetics parameters (red).

In the text

Fig. 14. Transient simulation, going from the critical configuration called H1 to different supercritical states (H2, H3, H4) by adjusting the water level in the tank. The results of the kinetic simulation in TRIPOLI-4 are compared with the solution of the point-kinetics equation, using the adjoint-weighted kinetics parameters. Results have been obtained both with the ENDF/B-VIII.0 (left) and the JEFF-3.3 (right) nuclear data library.

In the text