A POSTERIORI ERROR ESTIMATES FOR MIXED FINITE ELEMENT DISCRETIZATIONS OF THE NEUTRON DIFFUSION EQUATIONS

. We analyse a posteriori error estimates for the discretization of the neutron diffusion equations with mixed finite elements. We provide guaranteed and locally efficient estimators on a base block equation, the one-group neutron diffusion equation. We pay particular attention to AMR strategies on Cartesian meshes, since such structures are common for nuclear reactor core applications. We exhibit a robust marker strategy for this specific constraint, the direction marker strategy.


Introduction
The diffusion equation can model different physical phenomena, for instance Darcy's law, Fick's law or the neutron diffusion.Among models that are used in the nuclear industry, the multigroup neutron diffusion equation plays a central role [1].The base block is the one-group neutron diffusion equation.In [2,3], the first author and co-authors carried out the numerical analysis of this one-group neutron diffusion equation with a source term, discretized with mixed finite elements.The analysis included in particular the case of low-regularity solutions.A priori estimates were derived in the process.A natural question is then the a posteriori analysis of the method, to further optimize the cost of the numerical method.This is the main topic we address in this paper.
Nuclear reactor cores often have a Cartesian geometry.Indeed, in the models, the base brick, which is called a cell, is a rectangular cuboid of R 3 .The global layout is a set of cells that are distributed on a 3D grid, so that the global domain of the reactor core is represented by a rectangular cuboid of R 3 .Each cell can be made of fuel, absorbing or reflector material.To account for the different materials, the coefficients in the models are piecewise polynomials (possibly piecewise constant) with respect to the position, i.e. their restriction to each cell is a polynomial [1,11,12].In practice the coefficients characterizing the materials may differ from one cell to another by a factor of order 10 or more.
The outline is as follows.In Sections 2 and 3, we introduce some notations and our model problem.Then in Section 4, we recall how it can be solved in a mixed setting.To that aim we build the standard equivalent variational formulation, and provide the existing a priori numerical analysis results that allow one to compare the discrete solution to the exact one.For the discretization, we choose the well-known Raviart-Thomas-Nédélec finite element RTN  , where  ≥ 0 denotes the order.
In Section 5, we propose the a posteriori analysis of the model.We begin by the reconstruction of the solution (via post-processing), which can be devised in at least two ways: one is specific to the lowest-order, and the second one can be applied to any order.We also mention an averaging approach for the reconstruction.In Section 6, we propose some numerical experiments to compare the resulting strategies.For that, we focus on a specific discretization, based on Cartesian meshes.This kind of discretization is of particular importance for nuclear core simulations.

Notations
We choose the same notations as in [2,3].Throughout the paper,  is used to denote a generic positive constant which is independent of the mesh size, the mesh and the quantities/fields of interest.We also use the shorthand notation   for the inequality  ≤ , where  and  are two scalar quantities, and  is a generic constant.
If moreover the boundary  is Lipschitz, n denotes the unit outward normal vector field to .Finally, it is assumed that the reader is familiar with vector-valued function spaces related to the diffusion equation, such as H(div ; ), H 0 (div ; ) etc.
Specifically, we let Ω be a bounded, connected and open subset of R  for  = 2, 3, having a Lipschitz boundary which is piecewise smooth.We split Ω into  connected open disjoints parts {Ω  } 1≤≤ with Lipschitz, piecewise smooth boundaries: Ω = ∪ 1≤≤ Ω  and the set {Ω  } 1≤≤ is called a partition of Ω.For a field  defined over Ω, we shall use the notations   =  |Ω , for 1 ≤  ≤  .

The model
Given a source term   ∈  2 (Ω), we consider the following neutron diffusion equation, with vanishing Dirichlet boundary condition.In its primal form, it is written: {︂ Find  ∈  1 0 (Ω) such that −div D grad  + Σ   =   in Ω, (3.1) where , D, and Σ  denote respectively the neutron flux, the diffusion coefficient and the macroscopic absorption cross section.Finally,   denotes the fission source.When solving the neutron diffusion equation, D is scalarvalued.We choose to consider more generally that D is a (symmetric) tensor-valued coefficient.The coefficients defining Problem (3.1) satisfy the assumptions: a.e. in Ω. (3.2) Classically, Problem (3.1) is equivalent to the following variational formulation: Under the assumptions (3.2) on the coefficients, the primal problem (3.1) is well-posed, in the sense that for all   ∈  2 (Ω), there exists one and only one solution  ∈  1 0 (Ω) that solves (3.1), with the bound ‖‖ 1,Ω ‖  ‖ 0,Ω .Provided that the coefficient D is piecewise smooth, the solution has extra smoothness (see e.g.Prop. 1 in [2]).Instead of imposing a Dirichlet boundary condition on Ω, one can consider a Neumann or Fourier boundary condition    + (D grad ) • n = 0, with   ≥ 0. Results are similar.Throughout the paper, we add remarks on the extension to the situation where Σ  ≥ 0 may vanish.In particular, the analysis we propose covers both the pure diffusion case, and the diffusion-reaction case.

Variational formulation and discretization
Let us introduce the function spaces: From now on, we use the notations:  = (p, ) and  = (q, ).
One may prove that the mixed formulation (4.5) is well-posed, see Theorem 4.4 in [3].As a matter of fact, the result is obtained by proving an inf-sup condition in  , which we recall here.

Discretization and a priori error analysis
We study conforming discretizations of (4.5).Let ( ℎ ) ℎ be a family of meshes, made for instance of simplices, or of rectangles ( = 2), resp.cuboids ( = 3), indexed by a parameter ℎ equal to the largest diameter of elements of a given mesh.We introduce discrete, finite-dimensional, spaces indexed by ℎ as follows: The conforming discretization of the variational formulation (4.5) is then: Following Definition 2.14 in [13], we assume that (Q ℎ ) ℎ , resp.( ℎ ) ℎ have the approximability property in the sense that ∀q ∈ H(div , Ω), lim ℎ→0 We also impose that the space  0 ℎ of piecewise constant fields on the mesh is included in  ℎ , and that div Q ℎ ⊂  ℎ .We finally define: Remark 4.1.At some point, the discrete spaces are considered locally, i.e. restricted to one element of the mesh.So, one introduces the local spaces Q ℎ (),  ℎ (),  ℎ () for every  ∈  ℎ .
Provided the above conditions are fulfilled, one may derive a uniform discrete inf-sup condition under the same assumptions as in Theorem 4.1 (cf.[3], Thm.4.5).
Then the bilinear form  fulfills a uniform discrete inf-sup condition in  ℎ .
The classical a priori error analysis follows.Let  ℎ = (p ℎ ,  ℎ ) be the solution to (4.7).
Corollary 4.1.Under the assumptions of Theorem 4.2, there holds: Explicit a priori error estimates may be derived, see e.g.[3].
In this paper, we focus on the Raviart-Thomas-Nédélec (RTN) Finite Element [14,15].For simplicial meshes, that is meshes made of simplices, the finite element spaces RTN  can be described as follows, where  ≥ 0 is the order of the discretization for the scalar fields of  ℎ , see e.g.[16].
The boundary of a simplex  ∈  ℎ is made of the union of ( − 1)-simplices, called facets from now on, and denoted by (   ) 1≤≤+1 .We let P  () be the space of polynomials of maximal degree  on , resp.P  (   ) the space of polynomials of maximal degree  on    .The definition is Observe that for all q ∈ RTN  (), for all . The definitions of the finite element spaces RTN  are then For rectangular or Cartesian meshes, a description of the Raviart-Thomas-Nédélec (RTN) finite element spaces can be found for instance in Section 4.2 of [12].We consider those meshes explicitly for the numerical examples, see Section 6.

A POSTERIORI studies for a mixed finite element discretization
To develop the study of a posteriori estimates, we use the so-called reconstruction of the discrete solution  ℎ .In what follows, we denote by ζℎ := ζℎ ( ℎ ) a reconstruction, and by  := ( ζℎ ) an estimator.Classically, our aim is to obtain reliable and efficient estimators for the reconstructed error  − ζℎ , meaning that: where C and c are generic constants, and ‖ • ‖ is some norm to measure the error.To that aim, we consider that  =  1 0 (Ω), the original space of solutions, see (3.1), is the default space of scalar reconstructed fields.We also introduce the broken spaces A first approach has been suggested in [17], Chapter 8.The reconstruction ζℎ = (p ℎ , φℎ ) is defined as Remark 5.1.For other boundary conditions, i.e. for a Neumann or Fourier boundary condition, the default space  of scalar reconstructed fields would be equal to  1 (Ω).
In Section 5.1, we recall some reconstruction approaches for RTN finite element spaces.Section 5.2 is devoted to the derivation of a posteriori estimates.

Reconstruction of the discrete solution
In this section, we investigate some approaches to devise a reconstruction of the discrete solution (p ℎ ,  ℎ ), here obtained with the RTN  finite element discretization, for  ≥ 0.
For illustrative purposes, we consider simplicial meshes (see Rem. 5.2).Let us introduce some further notations, given such a mesh  ℎ .The set of facets of  ℎ is denoted ℱ ℎ , and it is split as ℱ ℎ = ℱ  ℎ ∪ ℱ  ℎ , with ℱ  ℎ (resp.ℱ  ℎ ) being the set of boundary facets (resp.interior facets).We denote by P  ( ℎ ) the space of piecewise polynomials of maximal degree  on each simplex  ∈  ℎ .We let   ℎ be the set of interpolation points (or nodes) where the degrees of freedom of the  -conforming Lagrange Finite Element space of order  are defined.And, for a node  ∈   ℎ , we denote by   the set of simplices  such that  ∈ .We recall the definition of the (original) Oswald interpolation operator [18] ℐ Os : P  ( ℎ ) → P  ( ℎ ) ∩  such that A second, modified Oswald operator is defined in [10] as follows.Let denotes the jump of  ℎ on the facet  ∈ ℱ  ℎ shared by elements  1 and  2 and n 1,2 is the unit outer normal of the mesh element  1,2 ∈  ℎ .Then, the modified Oswald operator 1 Above, the values (( ℎ ,   ))  ∈ℱ ℎ at the barycenters of the facets are then defined so that the mean value of ℐ MO ( ℎ ) on every facet is equal to the mean value of  ℎ on the same facet.
Remark 5.2.Observe that the results presented in this section can be extended to the case of rectangular or cuboid meshes [7].

Post-processing approaches
In order to recover the relation p = −D grad  at the discrete level, some post-processing techniques have been introduced for mixed finite element method [7,10].The first one is specific to a discretization with the RTN 0 finite element, whereas the second one can be applied to any discretization with a RTN  finite element, i.e.  can be any integer, possibly equal to 0.
For  = 0, the author in [10] chooses one post-processed scalar variable Problems (5.2) are local and independent on each element  ∈  ℎ .We define the RTN 0 post-processing by ℐ Os ∘ ℐ  .The reconstruction associated to the RTN 0 post-processing is On the other hand, for  ≥ 1, there exists no solution to Problem (5.2).We present here the approach proposed in [19], valid for  ≥ 0. It is shown there that the solution to (4.7),  ℎ = (p ℎ ,  ℎ ) ∈  ℎ , is also equal to the first argument of the solution of a hybrid formulation (see (5.4) below), where the constraint on the continuity of the normal trace of p ℎ is relaxed.Let be the space of the Lagrange multipliers and let Xℎ = Π ∈ ℎ  ℎ () be the unconstrained approximation space with the RTN  local finite element spaces.By definition,  ℎ is a strict subset of Xℎ .The hybrid formulation is: Let Π  ℎ : Xℎ × Λ ℎ →  ℎ be the projection onto an appropriate space 2 such that, given Remark 5.3.In the situation where Σ  ≥ 0 may vanish, the projection ̂︀ where for all  ∈  ℎ , (5.5) The space  ℎ is defined here as where P+2 ( ) denotes the  2 ( )-orthogonal complement of P +1 ( ) in P +2 ( ) for any facet  ∈ ℱ ℎ .We refer to [19] for the definition of ad hoc finite-dimensional spaces  ℎ for various families and types of elements.
Finally, we refer to [20] for an application of this technique in the field of neutronics.The RTN post-processing is defined here by The reconstruction associated to the RTN post-processing is (5.6)

Adding bubbles functions
This section details a possible correction of a reconstruction ζℎ = (p ℎ , φℎ ) to enforce the conservation of local averages, such as (5.14) below.It consists in adding bubble functions ( [21], Sect.3.2.2).The resulting reconstruction with bubble correction is defined as where   is the bubble function on  ∈  ℎ defined as the product of the barycentric coordinates of .

A posteriori error estimates
We now detail the derivation of a posteriori estimates.We define The definition is extended to piecewise smooth fields on  ℎ by replacing ∫︀ Ω by ∑︀ ∈ ℎ ∫︀  .Given  ∈  ℎ , we also define   0 the  2 ()-orthogonal projection on the space  0 ℎ (), and

•
In order to state the estimates, at some point we will use the following assumptions.Finally, we recall that there exists  , > 0, the so-called Poincaré constant (see e.g.Eq. (2.1) in [10]), such that, for all ℎ, for all  ∈  ℎ and for all  ∈  1 (), it holds that Note that  , = 1  in the case where the considered mesh elements are convex; cf.[22,23].In Section 5.2.1, we recall a classical a posteriori error framework in the primal setting (unknown ), where the error is measured in  1 ( ℎ ) norm.We propose two alternatives in the mixed setting (unknown (p, )): for the first one we measure the error with respect to the ℋ norm, while for the second one we use a weighted H(div ;  ℎ ) ×  2 (Ω) norm.Both approaches are respectively developed in Sections 5.2.2 and 5.2.3.

Estimates in 𝐻 1 norm
In this section, we aim to briefly recall the a posteriori error framework introduced in [10].The energy norm associated to the primal form is .
Theorem 5.1.Let  and  ℎ be respectively the solution to (4.5) and (4.7) with RTN 0 finite elements, and let For all  ∈  ℎ , we define the residual estimator and the nonconformity estimator .
Theorem 5.2 (Local efficiency of the a posteriori error estimators).Let  and  ℎ be respectively the solution to (4.5) and (4.7) with RTN 0 finite elements, and let ̂︀  ℎ = ℐ  ( ℎ ).Under Assumption 5.1, there holds on every with the constant C depending only on the polynomial degree  of   , the space dimension , and the shaperegularity parameter   = ||/ℎ   .

Estimates in ℋ norm
In this section, we use the broken norm associated to the bilinear form   , i.e.
For any  ∈  ℎ , we define the residual estimator where and   = min and the flux estimator (5.17) One has the reliability estimate (5.18) Proof.Using Lemma 5.1, we have, for all  ∈ ℋ, For any  ∈  , let  = (−D grad , ) ∈ ℋ.One has, by symmetry of D and according to (3.3), that Owing to the fact that  −  ∈  and the symmetry of D, we can integrate by parts the second term We also obtain where we used (5.14) and (4.7) in the first line, and we applied the Poincaré inequality (5.8) in the second line.Collecting (5.20) and (5.21) and using the definition of  , , we find Next, using (5.19) and the Cauchy-Schwarz inequality, we get We conclude the proof by choosing  = (−D grad φℎ , φℎ )∈ ℋ and using the definition of  , .
Remark 5.4.We recall the illuminating Prager-Synge theorem [24] which states that, given  ∈  the weak solution of (3.3), φℎ ∈  and p ℎ ∈ H(div , Ω) with div p ℎ + Σ  φℎ =   arbitrary, then one has the equality This result yields that, if one chooses a reconstruction that is H(div , Ω) ×  -conforming, one can derive simultaneous reliability and efficiency in a straightforward manner.
Theorem 5.4.Under the assumptions of Theorem 5.3, one has the reliability estimate Proof.Recall that φℎ ∈  .Since  solves (3.3), one has that .
Then, we have To reach the last line, we integrated by parts the second term.By the symmetry of D and the Cauchy-Schwarz inequality, we find for the second term (5.23) On the other hand, for the first term, one has by the Cauchy-Schwarz inequality We may also use (5.14) to write By definition (see (5.9)), we know that ‖Σ (5.26) And, with the help of (5.23), we get , which leads to the conclusion (5.22) by applying the Cauchy-Schwarz inequality one last time.
Theorem 5.5 (Local efficiency of the a posteriori error estimators).Let Assumption 5.2 be fulfilled.For  ∈  ℎ , let  , and  , be the residual and flux estimators respectively given by (5.15), and (5.17).In addition, we suppose that φℎ is piecewise polynomial on  ℎ .The following estimates hold true where c and C are constants which depend only on the polynomial degree of   , Σ  and φℎ , , and on the shape-regularity parameter   .
Proof.The proof follows that given in [10], Lemma 7.6.Let   be the bubble function on , given as the product of the  + 1 linear functions that take the value 1 at one vertex of  and vanish at the other vertices.
Let   = (  − div p ℎ − Σ  ̃︀  ℎ ).Note that   is a polynomial in , because each term appearing in its definition is a polynomial (thanks to Assumption 5.2 for   and Σ  ).Then the equivalence of norms on finite-dimensional spaces, the definition of   and the inverse inequality (cf., e.g.[25], Thm.3.2.6)respectively give ) ) with the constants  and  depending only on the polynomial degree of   , Σ  and φℎ , , and   .Let  , = (0,     ) in , and 0 elsewhere: by construction,  , ∈  .Then we have, by the definition of the bilinear form  and of , On the other hand, where we integrated by parts the second term in the second line.Combining (5.29)-(5.32),one comes to Considering the definition (5.15) of  , , we infer that , we have   = 1 and Otherwise, we have This concludes the proof of (5.27).We now proceed with the triangle inequality for the second estimate, Considering the definition of  , by (5.17) concludes the proof.
Remark 5.5.Assume in addition in Theorem 5.5 that there exists a constant  > 0, such that min ∈ ℎ   ≥ , for all ℎ > 0.Then, the constants c and C do not depend on   (but on ).
Remark 5.6.The results of this section extend with the same arguments to the situation where Σ  ≥ 0 may vanish.Under the assumptions of Theorem 5.3, one has the reliability estimates where the residual estimator  , =  ,   with and Under the assumptions of Theorem 5.5, one has the efficiency3 estimate Remark 5.7.Note that the estimators are robust with respect to the interplay of the sizes of D and Σ  .We refer to [26,27] for similar works on this issue.
For  ∈  ℎ , we introduce (, ). (5.34) Lemma 5.2.Let  and  ℎ be respectively the solution to (4.5) and (4.7).Let ζℎ = (p ℎ , φℎ )∈ Q ℎ ×  be a reconstruction of  ℎ .We have for all = (q, ) ∈  , Proof.Let  be in  .According to (4.5), we have Owing to the fact that φℎ is in  , we can integrate by part the last integral: This concludes the proof.
Theorem 5.7 (Local efficiency of the a posteriori error estimators).Let Assumption 5.2 be fulfilled.For  ∈  ℎ , let  , and  , be the residual and flux estimators respectively given by (5.16), and (5.17).In addition, we suppose that φℎ is piecewise polynomial on  ℎ .The following estimates hold true where c and C are constants which depend only on the polynomial degree of   , D, Σ  and φℎ , , and the shape-regularity parameter   .
Since the support of  , is equal to , one has actually  , ∈   .So, by definition (5.34) of the strengthened Combining (5.29), (5.30) and ( 5.40), one comes to Using the definition of  , by (5.16) concludes the proof of (5.38): .
We now proceed similarly for the second estimate.Let us denote q  = (D −1 p ℎ + grad φℎ ) on a given  ∈  ℎ .Note that q  is a polynomial in  (thanks to Assumption 5.2 for D −1 ).Then the equivalence of norms on finite-dimensional spaces, the definition of   and the inverse inequality (cf., e.g.[25], Thm.3.2.6)give ‖q  ‖ 2 0, ≤ (q  ,   q  ) 0, , (5.41) with the constants  and  depending only on the polynomial degree of D −1 and φℎ , , and   .Let  , = (  q  , 0) in , and 0 elsewhere.We observe that   q  is smooth in  (a closed subset of R  ), and moreover that (  q  ) | = 0 thanks to the definition of   .Hence,  , ∈   .According to Lemma 5.2 where we used the inverse inequality (5.43) to reach the last line.Combining (5.41), (5.42) and (5.44), one comes to Considering the definition of  , by (5.17) concludes the proof.
Remark 5.8.Assume in addition in Theorem 5.7 that there exists a constant  > 0, such that min ∈ ℎ   ≥ , for all ℎ > 0.Then, the constants c and C do not depend on   (but on ).
Remark 5.9.Similarly to Remark 5.6, the results of this section extend with the same arguments to the situation where Σ  ≥ 0 may vanish if one slightly modifies the definition of the norms by where Σ ⋆ is defined in (5.5).
Let us define for all  ∈  ℎ ,

•
Under the assumptions of Theorem 5.6, one has the reliability estimate , where the residual estimator becomes  ,⋆, =  ,⋆,  ⋆, with Under the assumptions of Theorem 5.7, one has the efficiency estimates

Numerical results
This section is devoted to the numerical experiments we performed on adaptive mesh refinement strategies (AMR).In fact, the AMR strategy can be classified into several categories: the ℎ-refinement (mesh subdivision), which amounts to refining the mesh where large errors occur [28]; the -refinement (local high order approximation), which increases the order of the polynomial functions [29], or the -refinement (moving mesh) that moves the nodes of the mesh to increase the mesh density [30], in the regions of interest where large variations of the solution occur.The above strategies can be mixed, such as ℎ-refinement [31,32] and ℎ-refinement [33].
We are interested in the case of heterogeneous coefficients which may induce some singularities in the solution of Problem (3.1), that is a loss of regularity of the solution due to the discontinuities in the data.Therefore, we focus on mesh subdivision strategy in this section.
The performance of adaptive mesh refinement is assessed with respect to various criteria such as the error estimator, the marker strategy and the threshold parameter.We recall that the context of our applications is modelling nuclear reactor cores, in particular geometries composed of rectangular cuboids of R 3 .This is the reason why the discretization in this section is applied on Cartesian meshes.We recall that we refer to [19] for the definition of the corresponding finite-dimensional spaces  ℎ , see footnote 2 page 7.
Mesh subdivision strategies are introduced in Section 6.1.Section 6.2 presents the set of test cases considered throughout the whole Section 6. Section 6.3 focuses on the marker strategies.The sensitivity with respect to the threshold parameter is investigated in Section 6.4.Section 6.5 compares various reconstruction approaches.Section 6.6 examines different error estimators.

An adaptive mesh refinement strategy
We recall in this section the ℎ-refinement approach.From the initial mesh  ℎ0 , the AMR strategy generates a sequence of meshes  ℎ  .This strategy corresponds in general to an iterative loop where at each iteration, we consider the following four modules: Assuming that the mesh  ℎ  is computed, module Solve indeed corresponds to solving Problem (4.7).The output of the module Estimate is (  ) ∈ ℎ  where   is an a posteriori error estimator.The stopping criterion of the algorithm is given by max where  AMR > 0. Module Mark returns the set of the marked cells based on the error estimators (  ) ∈ ℎ  .
In other words, this module selects a set of elements to be refined.To be convenient, for  ⊂  ℎ  , let us denote () = ( ∑︀ ∈  2  ) 1/2 .For a fixed threshold parameter 0 <  ≤ 1, the classical bulk-chasing criterion (Dörfler's marking strategy [34]) is to select the (smallest) set of elements such that () ≥ ( ℎ  ).( Lastly, module Refine performs the mesh refinement according to the selected mesh elements.This strategy is generic and can be applied to any kind of mesh. In order to preserve the Cartesian structure of the mesh at each iteration, it is essential to refine the mesh according to lines, along the directions (e  ) =1, , which contain at least one of the selected cells.As a consequence, it is obvious to see that the above cell marker strategy is extremely costly since we use the error indicator of some selected cells to refine the other cells located in the same line, for a given direction (e  ) =1, .Due to this drawback, it is relevant to consider some other marker cell strategies.Therefore, instead of using the classical bulk-chasing criterion defined on a single cell, we introduce: -aggregate error indicators according to a line containing the cell; -the sum of aggregate error indicators, taken on all lines containing the cell.
We call them respectively the direction marker and cross marker method.
The direction marker method consists in selecting for each direction e  ,  = 1, . . ., , the smallest set of lines where 0 <   ≤ 1 is a fixed threshold parameter.The resulting selected set is ∪ =1,   .
The cross marker method corresponds to selecting for each  ∈  ℎ , the smallest set of elements  ⊂  ℎ  such that ∑︁ ∈ (  ) ≥   ( ℎ  ), ( where 0 <   ≤ 1 is a fixed threshold parameter and   is the union over all directions (e  ) =1, of the lines containing .The resulting selected set is ∪ ∈   .Interestingly, performing the mesh refinement is straightforward with both the direction marker and cross marker methods.In addition, they both preserve the Cartesian structure of the mesh.

Setting of the test cases
This section is devoted to the definition of the test cases considered, namely the Dauge test case, the Checkerboard test case, the Center test case and the Rotation Center test cases.In the following test cases, we perform the numerical simulations on the domain Ω = (0, 1) 2 .We consider here a simple source term given by   = 1.Moreover, we assume that the diffusion coefficient D is a scalar, piecewise constant given by Figure 1.We also set Σ  = 1.
In the following, the initial mesh of the AMR strategy is chosen to be uniform in all directions.The mesh size of the initial mesh of the Dauge test case, the Checkerboard test case, the Center test case and the Rotation Center test cases are respectively equal to 1/4, 1/8, 1/6 and 1/8.

Influence of the marker cell strategy
We now study the influence of the marker cell strategy on the AMR approach for our set of test cases.In this section, the Dauge test case, the Center test case and the Checkerboard test case are performed with RTN 0 finite elements, while the Rotation Center test cases are performed with RTN 1 finite elements.The AMR strategies are based on the error estimator introduced in (5.36) and the average reconstruction (5.1).
The Dauge test case is a singular toy problem (see also in [2,17,35] and references therein for more details).In this test case, the singularity is located at (0.5, 0.5) and we expect refinement in this region.Adaptive mesh refinement is performed with a stopping criterion equal to  AMR = 2 × 10 −3 .Figure 2 shows that mesh refinement is more located near the singularity for the direction marker strategy than the other strategies.Moreover, Figure 3 shows that the direction marker needs three to four times less mesh elements than the other strategies.All the other test cases yield the same conclusions.
So, from now on, the adaptive mesh refinement is always performed with the direction marker method.

Sensitivity with respect to the threshold parameter
In this section, we evaluate the sensitivity with respect to the threshold parameter   defined in (6.2) on our set of test cases.For unstructured meshes like triangular meshes, the typical value for the threshold parameter   is 0.5.However, the choice of an optimal value for the threshold parameter   remains a difficult question.Therefore, we numerically investigate the optimal value of the threshold parameter.
The stopping criterion of the Checkerboard test case is  AMR = 5×10 −3 .For the other test cases, the stopping criterion is set to  AMR = 2 × 10 −3 .
For the sake of brevity, we only show Figure 4 that indicates that the optimal value of   for the Dauge test case is around 0.35. Figure 5 shows the numerical flux on the refined mesh with an optimal value of   for the different test cases.

Influence of the reconstruction
In this section, we investigate the influence of the reconstruction on the error estimator defined in (5.36).To this aim, we compare the reconstruction approaches defined in Section 5.1 on the Dauge test case.
First, the stopping criterion is fixed at  AMR = 1.5 × 10 −3 .As can be seen in Figure 6, the average reconstruction and RTN post-processing need more elements to reach the stopping criterion than the RTN 0 post-processing.It is related to the fact that the flux estimator is the dominant contribution of the total estimator.Figure 7 shows the numerical flux on the refined mesh for the different reconstructions.Second, we modify the stopping criterion.Now, the stopping criterion is based on the  2 error with respect to a reference solution.That is to say, the algorithm is stopped when    where  ref is a reference solution computed on a fine mesh.We fix  rel = 2 × 10 −2 .Figure 8 shows that the RTN 0 post-processing and RTN post-processing give similar AMR strategies and that the resulting mesh have fewer elements than with the average reconstruction.

Comparison of the error estimators
In this section, we perform the Dauge test case with a stopping criterion on the relative  2 error (6.4) with respect to a reference solution at  rel = 2×10 −2 .To be convenient, let Estimator 1, Estimator 2, Estimator 3 and Estimator 4 respectively stand for the error estimator defined in ( [17], Theorem 8.4), (5.36), (5.18) and (5.11).For the sake of completeness, we recall here the different estimators for all  ∈  ℎ , We apply the reconstruction associated to the RTN 0 post-processing with bubble correction (5.3)-(5.7) to Error estimators 1, 2 and 3 (we refer to Thm. 5.1 for Estimator 4).As can be seen in Figure 9, we obtain similar meshes for the various AMR strategies using Estimators 1, 2 and 4. On the other hand, there is more refinement near the discontinuities for Estimator 3. Also, the relative  2 error on the neutron flux are similar for the different AMR strategies according to Figure 10.

Conclusion
In this manuscript, we derive a posteriori estimates associated to different norms for the numerical solution of the neutron diffusion equation in mixed form.
We discuss the approach presented in [17], Chapter 8.Although reliability can be proven, it remains difficult to achieve local efficiency of the estimators.We address this issue by proposing a posteriori estimators that are both reliable and locally efficient.We also propose two norms to measure the errors.
Regarding the numerical aspects, we focus on Cartesian meshes, since such structures are relevant in nuclear core applications, and outline a robust marker strategy for this specific constraint, the direction marker strategy.We observe numerically that the AMR strategy is sensitive to the choice of the threshold parameter.We compare various a posteriori estimators under different criteria.We show that the choice of the reconstruction has a strong influence on the AMR strategy.The post-processing approaches are shown to be more efficient than the average reconstruction.In the case of the lowest-order Raviart-Thomas-Nédélec finite element, we observe that the RTN 0 post-processing gives a more accurate reconstruction compared to the RTN post-processing.Also, we compare the different estimators with the same choice of reconstruction.And we note that, if the stopping criterion is based on the  2 error with respect to a reference solution, the various refinement strategies yield similar results.
In a companion paper, we will consider more general models or settings.First, the extension of our a posteriori estimators to a Domain Decomposition Method, the so-called DD+ 2 jumps method.Then, we will study a more general model, widely used for nuclear core simulations, the multigroup diffusion problem, for which we will also provide a posteriori estimators.In both cases, these will be proven to be reliable and locally efficient.

Figure 2 .
Figure 2. Dauge test case: the numerical flux on refined meshes for different marker strategies with RTN 0 .(A) Cell marker.(B) Direction marker.(C) Cross marker.

Figure 3 .
Figure 3. Dauge test case: evolution of the number of elements and the maximum of the total error estimator for different marker cell strategies.(A) Evolution of the number of elements.(B) Evolution of the maximum of total error estimator.

Figure 4 .
Figure 4. Dauge test case with varying thresholds.(A) Evolution of the number of elements.(B) Evolution of the maximum of total error estimator.

Figure 6 .
Figure 6.Evolution of the number of elements and the maximum of the total estimator by using different reconstruction methods: namely the average reconstruction (5.1), the reconstruction associated to the RTN post-processing (5.6) and the reconstruction associated to the RTN 0 post-processing (5.3).(A) Number of elements.(B) Maximum of the total estimator.(C) Maximum of the residual estimator.

Figure 8 .
Figure 8. Evolution of the number of elements and the relative  2 error (6.4) for different reconstruction methods : namely the average reconstruction (5.1), the reconstruction associated to the RTN post-processing (5.6) and the reconstruction associated to the RTN 0 postprocessing (5.3).(A) Number of elements.(B) Relative  2 error.

Figure 10 .
Figure 10.Evolution of the number of elements and the relative  2 error (6.4) for different error estimators.(A) Number of elements.(B) Relative  2 error.