PHYSOR 2020: Transition to a Scalable Nuclear Future Cambridge, United Kingdom, March 29th-April 2nd, 2020 ADAPTIVE SOLUTION OF THE NEUTRON DIFFUSION EQUATION WITH HETEROGENEOUS COEFFICIENTS USING THE MIXED FINITE ELEMENT METHOD ON STRUCTURED MESHES

The neutron transport equation can be used to model the physics of the nuclear reactor core. Its solution depends on several variables and requires a lot of high precision computations. One can simplify this model to obtain the SPN equation for a generalized eigenvalue problem. In order to solve this eigenvalue problem, we usually use the inverse power iteration by solving a source problem at each iteration. Classically, this problem can be recast in a mixed variational form, and then discretized by using the RaviartThomas-Nédélec Finite Element. In this article, we focus on the steady-state diffusion equation with heterogeneous coefficients discretized on Cartesian meshes. In this situation, it is expected that the solution has low regularity. Therefore, it is necessary to refine at the singular regions to get better accuracy. The Adaptive Mesh Refinement (AMR) is one of the most effective ways to do that and to improve the computational time. The main ingredient for the refinement techniques is the use of a posteriori error estimates, which gives a rigorous upper bound of the error between the exact and numerical solution. This indicator allows to refine the mesh in the regions where the error is large. In this work, some mesh refinement strategies are proposed on the Cartesian mesh for the source problem. Moreover, we numerically investigate an algorithm which combines the AMR process with the inverse power iteration to handle the generalized eigenvalue problem.


Introduction
Numerical simulations of nuclear reactor core are generally expensive as they require the exact solution of the neutron transport equation. These simulations are computationally expensive since one has to deal with many variables such as the neutron position in space, the motion direction and neutron energy. In this work, we use the multi-group theory for the discretization of the energy variable and we focus on the steady state multi-group neutron diffusion equation. However, we choose to present the results in the form of the multi-group SP N equation since this extension is straightforward. In the time independent case, the multi-group SP N equation corresponds to a generalized eigenvalue problem: where R is a bounded and open subset of R d with d = 1, 2, 3. Let G be the number of energy group and N an odd number which represents for the order of the SP N equation. The number of odd and even moments isN = N +1 2 . The vector φ and p represent the even and odd components of the multi-group flux defined by: ∈ RN , p = (p 1 , . . . , p G ) ∈ R d×N ×G with p g = (p g 1 , p g 3 , . . . , p g N ) ∈ R d×N . The operators T o , T e , M f and H are defined by block matrices: In the above definition, χ g is the fission spectrum and the blocks T g o , T g e , S g →g o , S g →g e , M g f , H g are given by , H g = (δ n,m + δ n,m−1 )N n,m=0 where Σ g t is the total cross section, Σ g →g s is the scattering cross section and the removal macroscopic cross sections defined by Σ g r,n (x) = Σ g t,n (x) − Σ g→g s,n (x). The notation t n stands for the normalization factor defined by t 0 = 1 and t n = 4n 2 −1 nt n−1 for n > 0. The notation ν g is the number of neutron emitted per fission, and Σ g f is the macroscopic fission cross section. The multiplication factor k eff is the inverse of the eigenvalue and characterizes the physical state of the core reactor. In this work, we focus on the development of the project APOLLO3, a shared platform among CEA, FRAMATOME and EDF, which includes different deterministic solvers for the neutron transport equation. Particularly, we are interested in the MINOS solver based on the mixed Raviart-Thomas-Nedelec finite element method and implemented on Cartesian and hexagonal grids for the multi-group SP N equation [1]. Since the removal matrices T o and T e are usually heterogeneous (piecewise constant), the solution of the SP N equation may have some singularities which limit the precision and convergence of the solution [2]. It is well known that the mesh subdivision method is one of the most effective ways to treat this problem [3][4][5]. The purpose of this work is to develop an adaptive solution for the eigenvalue problem. In particular, this article is organized as follows. In Section 2, we present the Adaptive Mesh Refinement (AMR) strategies for the source problem by using a posteriori error estimators. Some marker cell strategies on Cartesian mesh are also detailed in this section. Next, we propose a splitting algorithm for the eigenvalue problem in Section 3. Some numerical test cases are shown in Section 4 to illustrate our purposes. Finally, concluding remarks in Section 5 complete the study.
Proceedings of the PHYSOR 2020, Cambridge, United Kingdom Physics of Reactors Transition to a Scalable Nuclear Future

Adaptive mesh refinement strategies
To be convenient, let us introduce the following function spaces: and the following operators: • : x g n y g n , : x g n · y g n .
Let X h be the discrete space associated to the mixed Raviart-Thomas finite elements. We denote The discrete variational formulation associated to Problem (1) reads: where the bilinear form B(ζ h , ξ h ) and f (ζ h , ξ h ) are respectively defined by The inverse power iteration algorithm is a well known method for the resolution of eigenvalue problem. At iteration m + 1, we deduce (p m+1 , φ m+1 , k m+1 eff ) from (p m , φ m , k m eff ) by solving the following source problem: where the source term is defined by S m = (k m eff ) −1 M f φ m . This iterative process is stopped when the relative error between the current and previous iteration of the flux as well as the multiplication factor k eff are smaller than the given threshold numbers. In this section, we aim to illustrate the adaptive mesh refinement strategies for the source problem (3). In general, the h-refinement is to generate a sequence of mesh T h k from the initial mesh T h 0 by using the AMR strategy which is in general an iterative loop where at each iteration, we consider the four modules: In particular, the module SOLVE is the mixed Raviart-Thomas finite element discretization on the mesh T h k of the source problem (3). In module ESTIMATE, the η K local error indicator on each element is obtained from a posteriori error estimate for the discrete solution. The purpose of the module MARK is to select a set of elements with large error to be refined. In the module REFINE, we use some rule for dividing the marked cells such that the obtained mesh T h k+1 is still conforming. We now detail how we obtain the error indicators for the source problem and some marker cell strategies on Cartesian mesh for the MINOS solver.
The associated discrete variational formulation reads, We now define the norm on X by ζ 2 where the notation dT stands for the diagonal of the matrix T. Finally, let us define the following norm |ζ| + = sup ξ∈X,||ξ|| S ≤1 B(ζ, ξ).
In the mixed dual Raviart-Thomas finite element method, we only have φ h ∈ L h . In particular, the discrete flux is a priori not in V . In order to find a posteriori error estimate for the source problem, it is crucial to introduce a discrete potential reconstruction fluxφ h ∈ V by using the average operator (see [6] and references therein). Next, let ζ and ζ h be respectively the solution of (4) and (5). Letζ h = (p h ,φ h ) be a reconstruction of ζ h in Q × V . For any K ∈ T h , we define respectively the residual estimator, flux estimator, symmetric non-conforming estimator and non-symmetric non-conforming estimator by: Then it stands The proof of this results for the multi-group SP N equation is an extension of the one for the diffusion one group equation proposed in [6,Chapter 8]. For the sake of brevity, details of the numerical analysis for the diffusion equation are given in the companion paper [7].

Refinement strategies on Cartesian mesh
The AMR method is usually applied to the non-structured mesh like triangular grids and the marker strategy is based on the total error indicator on each cell η K = η 2 R,K + η 2 F,K + η 2 N C,K + 4η 2 N C ,K . For example, the classical Dorfler's marking strategy aims to select the minimal cardinality of the subset S h k such that K∈S h k η K ≥ θ K∈T h k η K for a fixed threshold number 0 < θ ≤ 1. In this article, we focus on the solvers which use the structured meshes like the MINOS solver. In this case, to ensure the conformity of the mesh, it is essential to refine the mesh according to the whole line in each direction (x, y, z) which contain the selected cells. As a consequence, it is obvious to Physics of Reactors Transition to a Scalable Nuclear Future see that we use the error indicator of just some selected cells to refine the other cells located in the same line for a given direction. Therefore, it is extremely important to point out some other marker cell strategies based on more information. Instead of using the classical bulk chasing (Dorfler's marking strategy) defined on cell, we modify it to propose the construction of some other error indicators according to the lines of each direction and also on the "cross" value (the total error indicators of all the lines). In order to preserve Cartesian mesh structure, we just simply refine the selected lines by dividing into two lines.

Splitting algorithm for the generalized eigenvalue problem
In order to combine the AMR process with the power inverse iteration, a first approach is the combining algorithm where at each outer iteration, the AMR process is taken into account for the source problem [6,Chapter 8]. In this work, we would like to investigate the splitting algorithm is illustrated in Algorithm 1: The first process is the power inverse iteration for the eigenvalue problem and the second one is the AMR strategy for the source problem. The implementation of the splitting algorithm is easier than the combining one since we can use the solver of the problem (2) as a black box. The stopping criterion of the algorithm is usually defined by the maximum of the error indicators. However, for comparison purposes, if we have the reference value of the multiplication factor k eff , we can define another stopping criterion based on the relative error given by where k h eff is the numerical multiplication factor.

Algorithm 1:
The splitting algorithm Input: The threshold numbers E AMR or E Accuracy and the maximum of the refinement n max Output: The solution (p h k , φ h k , k h k eff ) associated to the mesh T h k .

The one group case
We present here a checkerboard test case for the one group diffusion equation on the domain of the computation R = (0, 100) 2 . The diffusion coefficient D = (3Σ t ) −1 is given as in Figure 1(a) with D = 5 in region and D = 1 in region. The fission cross section Σ f = 1 and the absorption cross section Σ a = Σ t − Σ s = 1. This test aims at investigating the effect of some marker strategies in the discussion above. Therefore, we compare these marker strategies for the splitting algorithm with n max = 1, θ = 0.5 and the lowest order Raviart-Thomas finite element RT 0 . We fix E Accuracy = 10 −5 and the reference eigenvalue is k reff = 0.995194 which is obtained from the computation on the uniform mesh 1000×1000. Figure 1(b) clearly shows that the direction marker is much better than the element and cross marker method since it requires less elements to reach the target accuracy. It was observed that the uniform mesh needs 32400 elements to obtain the same accuracy. Moreover, Figure 1(c)-1(e) illustrate that there are more refinement at singularities with the direction marker strategy which helps to improve the accuracy of the solution.

The multi-group case
We consider here the small FBR core proposed in [8] with 4 energy groups and where the control rod is withdrawn (model 2-case 1). However, we impose zero flux boundary condition on the boundary ∂R. Figure 2 shows the behavior of all error components. As can be seen, the structures of η N C ,K and the total error η K are very similar. Therefore, the non-conforming estimator η N C ,K is the main contributor to the error for this benchmark test case. Figure 3 shows some refined meshes of the refinement process.   Since there is no analytical value for k eff , the stopping criterion is given by E AM R = 0.02. The numerical multiplication factor and the total number of element are given in Table 1.

Conclusion
In this work, we propose an adaptive mesh refinement strategy for the neutron diffusion equation which relies on a posteriori error estimators for the source problem and a splitting approach to handle the generalized eigenvalue problem. Some marker cell strategies on the Cartesian mesh are also investigated and the direction marker method is the best candidate. Future work will be dedicated to the coupling between a posteriori error estimate and non-matching grid domain decomposition method to have locally refinement one each subdomain.