ADAPTIVE TIME-STEP CONTROL FOR THE MODAL METHOD TO INTEGRATE THE MULTIGROUP NEUTRON DIFFUSION EQUATION

The distribution of the power inside a reactor core can be described by the time dependent multigroup neutron diffusion equation. One of the approaches to integrate this timedependent equation is the modal method, that assumes that the solution can be described by the sum of amplitude function multiplied by shape functions of modes. These shape functions can be computed by solving a λ-modes problems. The modal method has a great interest when the distribution of the power cannot be well approximated by only one shape function, mainly, when local perturbations are applied during the transient. Usually, the shape functions of the modal methods are updated for the time-dependent equations with a constant time-step size to obtain accurate results. In this work, we propose a modal methodology with an adaptive control time-step to update the eigenfunctions associated with the modes. This algorithm improves efficiency because of time is not spent solving the systems to a level of accuracy beyond relevance and reduces the step size if they detect a numerical instability. Step size controllers require an error estimation. Different error estimations are considered and analyzed in a benchmark problem with a out of phase local perturbation.


INTRODUCTION
The distribution of the neutrons inside a reactor core along time can be described by the time dependent multigroup neutron diffusion equation. This equation depends on the position and the time. A high order finite element method (FEM) is considered to make the spatial discretization of the neutron diffusion equation to get a system of semi-ordinary differential equations [1] V −1 dΦ dt where Φ is the algebraic neutron flux, C is the algebraic concentration of precursors, V is the matrix with the velocities in their diagonal, L is the discretized leakage operator, M the discretized production operator. The matrix M 1 corresponds to the discretization of the first row of the differential production operator. Finally, if I is the identity matrix, the matrix X is equal to The FEM has been implemented by using the open source library Deal.ii [2]. This system of ODE is, in general, stiff. Several approaches have been studied to integrate this time-dependent equation. For instance, one can use the backward differential method that is based on the implicit Euler method. Other one is the quasi-static method, that decomposes the solution into an amplitude function (which have a fast change over time) and a shape function (which change slowly over time). Other strategy is the modal method, that assumes that the solution can be described by the sum of several amplitude functions multiplied by shape functions computed by solving the λ-modes problem Associated to this spatial problem, we can define the adjoint problem associated where the eigenvalues of both problems are equal. The adjoint functions obtained, ψ † l , l = 1, . . . , q satisfy the biorthogonality condition ψ † l , M ψ m = ψ † m , M ψ m δ l,m , l, m = 1, . . . , q.
The modal approach has a great interest when the distribution of the power cannot be well approximated by only one shape function, mainly, when local perturbations are applied during the transient. Usually, the shape functions of the modal methods are updated for the time-dependent equations with a constant time-step size to obtain accurate results [3]. Recent time schemes have been developed to select automatically an appropriate step size from a given tolerance. These algorithms have the advantage that they do no wast time for solving the systems to a level of accuracy beyond relevance and reduce the step size if they detect a numerical instability. In addition, they remove the necessity of selecting an adequate time step before the simulation has been run.
In this work, we propose an adaptive control time-step for the modal methodology to integrate the time-dependent neutron diffusion equation.

THE UPDATED MODAL METHOD
The modal methodology supposes that Φ( r, t) admits the following expansion where n m (t) are the amplitude coefficients and ψ m ( r) the shape functions obtained of solving the problem (3). It also expresses the matrices L and M of the problem (3) as where L 0 and M 0 are the matrices at t = 0 and in critical state (M 0 is divided by k ef f = λ 1 ). First, the expansion (4) and the assumption (5) is substituted into Equation (1) to get Now we use the definition of the λ-modes problem (3) and then, the resulting system is multiplied (by the left) by the adjoint modes ψ † l , l = 1, . . . , q, to obtain the system of ODE's d dt where The initial conditions values are that are obtained from the equations in the critical state. The ODE is solved with the open source library SUNDIALS, in particular with the CVODE module [4]. The time-step to solve this system is selected by an internal control step implemented in the library.

Updated modal methodology
In realistic transient computations, the flux can suffer extremely spatial variations. Sometimes, obtaining good approximations implies a high number of modes that means high computational cost [5]. As a solution it was proposed a modal methodology where the modes are updated in a certain time interval as we described in the following. Therefore, the time domain is divided into , the neutron diffusion equation can be integrated by using the λ-modes associated with the problem where L i and M i are the matrices associated with the nuclear reactor at time t i .
The differential equations that are needed to integrate at each time-step have the same structure than the problem without updating (Equation (7)). However, the initial conditions for n i,λ m at time t i must be defined to guarantee the continuity of the solution in the limits of the intervals. In the interval and the value of Φ(t i ) from the solution computed in the previous modal step. Likewise, to compute . The concentration of precursors at time t i can be computed as where The modes are computed with the Block inverse-free preconditioned Arnoldi method (BIFPAM). The BIFPAM is a block method based on the development of the Krylov subspaces associated with the residual error (see more details in [6]). Moreover, the computation of these modes and the modal methodology have been implemented using a matrix-free strategy, that avoids the assembly of the matrices and saves a lot of computer memory for the storage.

ADAPTIVE TIME-STEP CONTROL
The modes can be updated with fix time-step. This implies the necessity to select a time-step previously (that leads to obtain results with unpredictable errors). On the other hand, if we set small time-steps to obtain accurate approximations, the computational cost also increases. For this reason, we propose an adaptive control of the time-step. We need to define an error estimation and a suitable constraint to select the time-step based on the error estimation.
Estimation error. The first estimation is based on the difference between eigenfunctions. One can compute the modes in the next time-step to predict how the total flux will change. The (modal difference error) is defined as The second approach is based on the residual error that appear when the actual modes are substituted on the problem in other time step. The modal residual error is defined as Finally, we assume that the flux along the time will change depending on the variation in the cross-sections. The cross-section perturbation error is defined as where c is a cell of the spatial discretization.
Control Algorithm. We have studied two strategies based on the errors in the previous step. The first one, the Banded Error Control, changes the time-step in a fixed way The values min le and max le are dependent on the transient studied. The second option, the Adaptive Error Control, is based on the control algorithms defined for differential methods [7]. In particular, we obtain the step ∆t i as

NUMERICAL RESULTS
The performance of the adaptive updated modal method is tested in a transient defined from the Langenbuch reactor [8]. A top view of this reactor is represented in Figure 1. The transient is defined from an out of phase perturbation in the fission cross-sections of material 1 with striped pattern. The perturbation applied to the initial cross-section is, in P 1 and P 2 cells, respectively δΣ P 1 f,g (t) = 10 −3 sin(2πt), δΣ P 2 f,g (t) = 10 −3 sin(2πt + π) g = 1, 2.
The diffusion equation for all cases has been discretized by using polynomials of degree 3 in the FEM. The methodology proposed is compared with the backward differential method (BKM) [1]. The CPU time to obtain the solution with the BKM has been 232 min. This transient cannot be solved by using only one mode in the modal expansion and high number of modes are needed to used to obtain accurate approximations without updating [5]. The evolution of the global power computed with the BKM and with the updated modal method with some fix time-step is represented in Figure  2. It is observed that large errors are produced when the perturbations reach their maximums.
We analyze the local error (difference between the power obtained with the BKM and the modal method) for some settings of the updated modal method computed with the matrix-free method. Without updating the number of modes required to obtain accurate results and the CPU Time are very large [5]. Table 1 shows that small time-steps in the updating gives more accurate approximations. In contrast, the CPU times are also higher. On the other hand, if the number of modes is compared, one can observe that better approximations are obtained when the number of modes is higher. In this case, a selection of 4 modes and time-step equal to 0.1 s is the most effective option. However, these parameters can change for other reactor computations. Figure 3 represents the local error along the time. It shows large errors near to the extremes of the power (maximums and minimums) and before to the modal updating. More distributed errors are obtained when the time-step is smaller.
To test the adaptive modal method, Table 2 shows the local errors and CPU times obtained by setting the different error estimations and control errors with 4 eigenvalues. One can deduced that the modal difference error (md) is not very efficient because it needs to compute the modes to estimate the error that is very expensive. Regarding the other error estimations, there are not big differences between them. If the type of control error is compared, the adaptive control gives similar approximations than the fixed control. Figure 4 shows the evolution of the ∆t computed with the adaptive error control.

CONCLUSIONS
In this work, we have studied an adaptive modal method to integrate the time-dependent neutron diffusion equation obtained from applying an out-of-phase perturbation. The numerical results have shown that the modal expansion has to be updated with a small time-step and a considerable number of modes to obtain accurate results. The implementation of the adaptive control time-step allows not to set the time-step previously and it is adapted depending on an error estimation in the modal expansion. The results prove that the adaptive control decreases the errors in similar CPU times. Moreover, it has been observed that the modal method is an efficient way to integrate the time dependent equation because can be implemented with a matrix-free strategy saving up memory resources and computational times. In future works, this methodology will studied for other approximations of the neutron transport theory where the computational cost is very high.