# DISCONTINUOUS GALERKIN DISCRETIZATION AND HP-REFINEMENT FOR THE RESOLUTION OF THE NEUTRON TRANSPORT EQUATION

Raphaele Herbin, Damien Fournier, Romain Le Tellier

► To cite this version:

Raphaele Herbin, Damien Fournier, Romain Le Tellier. DISCONTINUOUS GALERKIN DISCRETIZATION AND HP-REFINEMENT FOR THE RESOLUTION OF THE NEUTRON TRANSPORT EQUATION. SIAM Journal on Scientific Computing, 2013, 35 (2), 10.1137/110844581. hal-00781329

HAL Id: hal-00781329

https://hal.science/hal-00781329

Submitted on 25 Jan 2013

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
DISCONTINUOUS GALERKIN DISCRETIZATION AND
HP-REFINEMENT FOR THE RESOLUTION OF THE NEUTRON
TRANSPORT EQUATION

D. FOURNIER*, R. HERBIN†, AND R. LE TELLIER ‡

**Abstract.** This paper presents a *hp*-refinement method for a first order scalar transport-reaction equation discretized by a discontinuous Galerkin method. First, the theoretical rates of convergence of *h*− and *p*−refinement are recalled and numerically tested. Then, in order to design some meshes, we propose two different estimators of the local error on the spatial domain. These quantities are analysed and compared depending on the regularity of the solution so as to find the best way to lead the refinement process and the best strategy to choose between *h*− and *p*−refinement. Finally, the different possible refinement strategies are compared first on analytical examples and then on realistic applications for neutron transport in a nuclear reactor core.

**Key words.** transport equation, finite volume, discontinuous Galerkin, error estimates, *hp*-refinement, nuclear core

**AMS subject classifications.** 74S10, 65N30, 65L50

**1. Introduction.** Accurate solutions of PDEs require important memory storage and computing time. Mesh adaptation is a way to improve both last points without damaging the convergence. The adaptive strategy is classically driven by an estimation of the local error. This local error estimate must be computable from the approximate solution and some constants related to the problem. It should, of course, approximate the real error as well as possible but be also easy to compute, i.e. have a negligible computation time compared to the CPU time required to solve the problem. For hyperbolic problems discretized by a discontinuous Galerkin (DG) method, error estimators are proposed for example in [2] for a solution belonging to $C^\infty$ and in [10, 9] for less regular functions. The refinement can also be led with respect to the error estimated for a quantity of interest [18], in this case a solution of the dual problem is also required.

With the DG method, once the estimator is chosen, the refinement can be done
according to two different strategies: refining the mesh is called *h*-refinement ([20, 19]), while increasing the degree of the polynomial basis is called *p*-refinement ([7, 3]). The combination of these two approaches (*hp*-refinement) was theoretically introduced in [5] but without a criterion to choose between both strategies. The selection is generally done by estimating the regularity of the solution [23, 4], by using theoretical results of convergence [6, 25], or by using a priori information [24].

In this paper, a *hp*-refinement strategy is proposed in the context of the neutron transport equation. The framework of study is presented in Section 2 and the properties of *h*− and *p*−refinement on this equation are derived in Section 3. Then, error estimators inspired from [2, 12, 10] are derived in Section 4 and analysed numerically in Section 5 on benchmarks for which the regularity of the solution is the one that is expected for realistic applications. A *hp*-algorithm combining both

<footnote>*CEA, DEN, DER/SPRC/LEPh, Cadarache, F-13108 Saint Paul-lez-Durance, France (damien.fournier@cea.fr).</footnote>

<footnote>†Laboratoire d'Analyse et de Topologie de Marseille, Centre de Mathématiques et Informatique (CMI), Université de Provence, Technopôle Château-Gombert, 39, rue F. Joliot Curie, 13453 Marseille Cedex, France.</footnote>

<footnote>‡CEA, DEN, DER/SPRC/LEPh, Cadarache, F-13108 Saint Paul-lez-Durance, France</footnote>
error estimators is deduced, improved by the analysis of Section 3 and compared to other hp-strategies in Section 6.

**2. Framework of study.** The neutron density in a nuclear reactor depends on the energy, direction of motion and spatial position of the particles. The equation is first discretized with respect to the energy and angle [21]; a system of spatial coupled problems has then to be solved on a convex polygonal domain $D \in \mathbb{R}^d$ ($d = 2, 3$) with boundary $\Gamma$. Denoting the inflow boundary $\Gamma_-$ = {$x \in \Gamma \text{ s.t. } \Omega \cdot \nu < 0$} where $\nu$ is the outgoing normal, the spatial neutron transport equation relative to a given energy $E$ and direction $\Omega$ can be written [21] as

$$ \begin{cases} \Omega \cdot \nabla u + cu = g & \text{in } D, \\ u(x) = u_0(x) & \text{on } \Gamma_-, \end{cases} \qquad (2.1) $$

where $u$ is the so-called neutron flux that is the product of the neutron density by the speed, $c$ is related to the probability of interaction between the neutrons and surrounding nuclei, piecewise-constant on $D$, $\Omega$ a given constant vector belonging to the unit sphere and $g$ a source term. In classical applications, $u$, $g$ and $u_0$ belong to $H^{\frac{1}{2}}$ for cores with embedded void regions or $H^{\frac{3}{2}}$ otherwise [27]. This problem is well posed [11, page 231] if $u$ is sought in the space $W = \{w \in L^2(D); \Omega \cdot \nabla u \in L^2(D)\}$ equipped with the norm $\|w\|_W = \|w\|_{L^2(D)} + \|\Omega \cdot \nabla w\|_{L^2(D)}$.

The spatial discretization is ensured by a DG method. Let us denote the mesh discretization by $\mathcal{T}$; in our framework the control volumes are quadrilateral, and $h = \sup \{h_K = \text{diam}(K), K \in \mathcal{T}\}$ denotes the space step. The grid $\mathcal{T}$ is possibly non conforming in order to enable $h$-refinement. The set of edges of each cell $K$ is denoted by $\partial K$ and $\partial K_- = \{x \in \partial K \text{ s.t. } \Omega \cdot \nu < 0\}$ (resp. $\partial K_+ = \{x \in \partial K \text{ s.t. } \Omega \cdot \nu > 0\}$) are the inflow (resp. outflow) edges. The inflow and outflow traces $u_-$ and $u_+$ of a function $u$ on $\partial K$ denoted ) are defined by

$$ u_{\pm}(x \in \partial K) = \lim_{\epsilon \to 0^{\pm}} u(x + \epsilon \Omega). \qquad (2.2) $$

In order to derive the DG scheme, we multiply Eq. (2.1) by a test function $v \in W$ and integrate on a cell $K$; this yields

$$ \int_K (-\Omega \cdot \nabla v + cv) u + \int_{\partial K} \Omega \cdot \nu u|_K v|_K = \int_K gv. \qquad (2.3) $$

Approximate solutions of Eq. (2.1) are computed by using discontinuous piecewise polynomials of degree $p_K$ on each cell $K \in \mathcal{T}$ of the mesh, i.e; functions belonging to the space $V_h^p$, where $p = \max_{K \in \mathcal{T}} p_K$, defined by

$$ V_h^p = \{v_h^p \in L^2(D) \text{ s.t. } v_h^p|_K \in \mathbb{Q}_{p_K}, p_K \le p, \forall K \in \mathcal{T}\}, $$

and

$$ \mathbb{Q}_{p_K} = \left\{ P \text{ s.t. } \exists a_{ij}, P = \prod_{i=1}^{d} \sum_{j=1}^{p_K} a_{ij} x_i^j \right\}. $$

The discrete problem is then obtained by considering Eq. (2.3) with $u_h^p$ and $v_h^p$ belonging to $V_h^p$ instead of $u$ and $v \in W$ (note that $V_h^p \subset W \subset L^2(D)$); we integrate the
resulting equation by parts on each element $K$, and choose, for stability reasons, an upwind conservative numerical flux $F(\Omega, u_h^p)$ for the approximation of $\Omega \cdot \nu u_h^p$, defined by

$$F(\Omega, u_h^p) = \begin{cases} \Omega u_0 & \text{on } \partial K_- \cap \Gamma_-, \\ \Omega u_{h,-}^p & \text{on } \partial K \setminus \Gamma_- \end{cases}. \qquad (2.4)$$

We thus get

$$\int_K u_h^p (cv_h^p - \Omega \cdot \nabla v_h^p) + \int_{\partial K} \nu \cdot F(\Omega, u_h^p) v_h^p|_K = \int_K g v_h^p. \qquad (2.5)$$

Summing Eq. (2.5) over all the control volumes $K$ of the mesh $\mathcal{T}$ leads to the DG scheme:

$$\left. \begin{array}{l} \text{Find } u_h^p \in V_h^p \text{ such that} \\ \qquad B(u_h^p, v_h^p) = L(v_h^p), \quad \forall v_h^p \in V_h^p, \qquad (2.6) \\ \text{where} \\ \qquad B(u_h^p, v_h^p) = \sum_K \left( \int_K (-\Omega \cdot \nabla v_h^p + cv_h^p) u_h^p + \int_{\partial K_- \setminus \Gamma_-} |\Omega \cdot \nu| (v_{h,-}^p - v_{h,+}^p) u_{h,-}^p \right. \\ \qquad \qquad \left. + \int_{\partial K_+ \cap \Gamma_+} |\Omega \cdot \nu| u_{h,-}^p v_{h,-}^p \right) \\ \qquad L(v_h^p) = \sum_K \left( \int_K g v_h^p + \int_{\partial K_- \cap \Gamma_-} |\Omega \cdot \nu| u_0 v_{h,+}^p \right). \end{array} \right\}$$

Integrating by parts, the bilinear form $B$ can also be written as

$$B(u_h^p, v_h^p) = \sum_K \left( \int_K (\Omega \cdot \nabla u_h^p + cu_h^p) v_h^p + \int_{\partial K_- \setminus \Gamma_-} |\Omega \cdot \nu| (u_{h,+}^p - u_{h,-}^p) v_{h,+}^p \right. \\ \left. + \int_{\partial K_- \cap \Gamma_-} |\Omega \cdot \nu| u_0 v_{h,+}^p \right).$$

It can be shown that the (linear) system Eq. (2.6) admits a unique solution $u_h^p \in V_h^p$. In fact, assuming that the solution $u$ of Eq. (2.1) belongs to $H^s$, $s \ge 0$, and denoting $u_h^p$ a solution of Eq. (2.6), Richter [28] and more recently Cockburn [9] proved that there exists $C_p \ge 0$, depending only on $p$ such that

$$\epsilon_{p,h} = \|e_h^p\|_{L^2(\mathcal{D})} = \|u - u_h^p\|_{L^2(\mathcal{D})} \le C_p (h^{\min(p+1,s)}). \qquad (2.7)$$

**3. Comparison of h- and p-refinement strategies.** In this section, an analysis of the behavior of the error depending on $h$ and $p$ is performed to understand the issues of hp-methods. The rate of convergence for uniform h- and p-refinement is derived, then checked numerically. Since the regularity of the solution plays an important role in the rate of convergence, two benchmarks with exact solutions of different regularities are used in the following. The functions $c$ and $g$ defined in Eq. (2.1) are
constant respectively equal to 1 and 0. Two configurations presented in [22] are of particular interest:

<figure><img src="image_2.png" /><figcaption>FIGURE 3.1. Boundary conditions for the MMS0 and MMS1 benchmarks</figcaption></figure>

* The MMS0 benchmark (MMS stands for Method of Manufactured Solutions [22]) generates non continuous solutions (belonging to $H^{1/2}$) which corresponds to the regularity obtained for cores with embedded void regions. It is obtained by imposing a discontinuous flux at incoming boundaries, leading to a discontinuous solution along the characteristic lines.

* The MMS1 benchmark generates continuous solutions but with discontinuous first derivatives (belonging to $H^{3/2}$) by imposing a constant incoming flux. As for the MMS0, the non-regularity is located on the characteristic lines.

<figure><img src="image_3.png" /><figcaption>FIGURE 3.2. Representation of the solution of the MMS0 and MMS1 benchmarks</figcaption></figure>

The boundary conditions and the representation of the solutions are presented in Figures 3.1 and 3.2.

**3.1. Analysis of uniform h-refinement.** The theoretical error estimate given by Eq. (2.7) of [28, 9] was verified numerically in [30] for solutions of escalating degrees of smoothness on triangular meshes. It was shown that the rate of convergence is limited by the regularity except if the mesh is aligned with singularities. In this last case, a $h^{p+1}$ convergence is observed.

In the MMS0 or MMS1 test cases, since $s$ is respectively equal to $1/2$ and $3/2$, Eq. (2.7) implies that there exists $C_p$ independent of the mesh size such that

$$ \epsilon_{p,h} = C_p h^s = \frac{C_p'}{\text{dof}^{s/2}} \qquad (3.1) $$

where 'dof' denotes the number of degrees of freedom:

$$ \text{dof} = \sum_{K \in \mathcal{T}} (p_K + 1)^d, \qquad (3.2) $$

and $p_K$ the degree of the polynomial on $K$.
<figure><img src="image_4.png" /><figcaption>FIGURE 3.3. L²-error for the MMS0 (left) and MMS1 (right) benchmarks with uniform refinement (log-log scale)</figcaption></figure>

For a uniform Cartesian mesh and constant polynomial degree, the number of degrees of freedom is given by

$$ \text{dof} = \text{nbCells} \times (p + 1)^d = \left(\frac{p+1}{h}\right)^d \qquad (3.3) $$

where nbCells is the number of cells $K \in \mathcal{T}$.
Eq. (3.1) implies

$$ \ln(\epsilon_{p,h}) = \ln(C'_{p}) - \frac{s}{2} \ln(\text{dof}) \qquad (3.4) $$

This result can be verified numerically on the MMS0 and MMS1 benchmarks by representing $\epsilon_{p,h}$ as a function of dof as shown in Figure 3.3. The slope is independent of $p$ and is equal to $\frac{1}{4}$ (resp. $\frac{3}{4}$) for the MMS0 (resp. MMS1) case as expected from Eq. (3.4). These numerical results show that the estimate Eq. (2.7) is sharp at least with respect to the power of $h$. Since the power of $h$ is limited by $s$, one may consider that increasing the polynomial basis order is useless to improve the accuracy of the solution. However, it may be useful because of the strong dependency of the parameter $C_p$ on $p$, as we shall see in Section 3.2.

**3.2. Analysis of uniform p-refinement.** Results concerning the rate of convergence of p-refinement methods are less developed particularly in the framework of hyperbolic problems, probably because of the lack of regularity of the exact solutions. The first results were obtained by using stabilization methods [8]. With these methods, the weak form Eq. (2.6) is modified and the discretized scheme becomes:

$$ \tilde{B}(u_h^p, v_h^p) = \sum_K \left( \int_K (\Omega \cdot \nabla u_h^p + c u_h^p) \left( v_h^p + \frac{h_K}{p_K^2} \Omega \cdot \nabla v_h^p \right) \right. \\ \left. + \left( 1 + \frac{h_K}{p_K^2} \right) \int_{\partial K_\perp \setminus \Gamma} |\Omega \cdot n| (u_h^{p,+} - u_h^{p,-}) v_h^{p,+} \right. \\ \left. + \left( 1 + \frac{h_K}{p_K^2} \right) \int_{\partial K_\perp \cap \Gamma} |\Omega \cdot n| u_h^{p,+} v_h^{p,+} \right) \qquad (3.5) $$

$$ \tilde{L}(v_h^p) = \sum_K \left( \int_K g \left( v_h^p + \frac{h_K}{p_K^2} \Omega \cdot \nabla v_h^p \right) + \left( 1 + \frac{h_K}{p_K^2} \right) \int_{\partial K_\perp \cap \Gamma} |\Omega \cdot n| u_0 v_h^p \right) \qquad (3.6) $$
Defining a mesh-dependent norm by

$$ \|v_h^p\|_{hp,\mathcal{D}}^2 = \sum_K \left( \frac{h_K}{p_K^2} \|\Omega \cdot \nabla v_h^p\|_{L^2(K)}^2 + \|v_h^p\|_{L^2(K)}^2 \right. \\ \left. + \int_{\partial K \setminus \Gamma_-} |\Omega \cdot \nu| (v_{h,+}^p - v_{h,-}^p)^2 + \int_{\partial K \cap \Gamma} |\Omega \cdot \nu| (v_h^p)^2 \right) \quad (3.7) $$

Bey et al. [8] proved that, if the solution belongs to $H^s$, $s \ge 1$,

$$ \|e_h^p\|_{hp,\mathcal{D}} = \mathcal{O}\left(\frac{h^{\nu-\frac{1}{2}}}{p^{s-1}}\right) \quad (3.8) $$

where $\nu = \min(p+1, s)$. However, this result is valid only for a quasi-uniform mesh i.e if there exists positive constants $\tau$ and $\sigma$ such that

$$ \frac{h}{h_K} \le \tau \qquad \frac{h_K}{\rho_K} \le \sigma \quad (3.9) $$

where $\rho_K = \sup\{\text{diam}(S_K) \text{ for } S_K \text{ sphere contained in } K\}$.

They noted that this result is not optimal with respect to the $L^2$-norm of the error. The limiting term in this norm is related to the $\Omega \cdot \nabla$ part. Actually, the $\frac{h_K}{p_K^2}$ term was necessary to deal with the slower convergence of the $\Omega \cdot \nabla$ part. This term was bounded using the following lemma [8]:

LEMMA 3.1. Let $w_K$ be a polynomial of degree $p_K$ on a cell $K \in \mathbb{R}^2$ of diameter $h_K$ where $K$ is obtained by an affine transformation of the master element $\hat{K} = [-1, 1] \times [-1, 1]$. Then,

$$ \|\Omega \cdot \nabla w_K\|_{L^2(K)} \le C \frac{p_K^2}{h_K} \|w_K\|_{L^2(K)} \quad (3.10) $$

where $C \in \mathbb{R}$ is a constant depending only on $\Omega$ and on the geometrical parameters $\sigma$ and $\tau$ (see Eq. (3.9)).

Houston et al. [17] found a way to sharply bound the gradient term. They modified the stabilization parameter in Eqs. (3.5, 3.6) and the norm definition in Eq. (3.7) (by removing the terms in bold from these equations). By doing so, they proved that if $s \ge \frac{1}{2}$,

$$ \|e_h^p\|_{hp,\mathcal{D}} = \mathcal{O}\left(\frac{h^{\nu-\frac{1}{2}}}{p^{s-\frac{1}{2}}}\right) \quad (3.11) $$

This bound is optimal with respect to $p$ [17]. They proved [18] that the stabilization is not necessary if one just wants to bound the $L^2$ norm of the solution.

This result is verified on the MMS benchmarks in Figure 3.4. The benchmark problems are discretized with only one control volume and the polynomial order is uniformly increased. The asymptotic result of $O(p^{-s+\frac{1}{2}})$ is checked. For the MMS0 benchmark, $\epsilon_{L^2}$ is constant for large values of $p$; this is expected as the method cannot converge for this test since $\nu = s = \frac{1}{2}$. Similarly, $\epsilon_{L^2}$ converges asymptotically
<figure><img src="image_5.png" /><figcaption>FIGURE 3.4. L²-error versus p for the MMS0 (left) and MMS1 (right) benchmarks</figcaption></figure>

as $O(p^{-1})$ for the MMS1 case. But another phenomenon is observed. There exists a preasymptotic region where the convergence is faster ($O(p^{-s})$). It corresponds to the best interpolation result for a polynomial of degree $p$ proved in [6]:

LEMMA 3.2. Let $K \in \mathcal{T}$ and $u \in H^s$, $s \ge 0$. Then, there exists $C \in \mathbb{R}$ depending only on $\Omega$, $\sigma$ and $\tau$ defined Eq. (3.9), and a function $z \in \mathbb{P}_{p_K}$, $p_K > 0$, such that

$$ \|u - z\|_{L^2(K)} \le C \frac{h_K^{\nu_K}}{p_K^s} \|u\|_{H^s(K)} \quad (3.12) $$

where $\nu_K = \min(p_K + 1, s)$.

This result is a possible explanation for the superconvergence in the preasymptotic region. Thus, increasing $p$ during the first steps can decrease the error even for a discontinuous solution.

3.3. Comparison of the convergence of h- and p-refinement. By using the results of the previous sections, the rate of convergence of the error for uniform refinements can be written as:

$$ \|e_h^p\|_{L^2(\mathcal{D})} = C \frac{h^s}{p^\alpha} \quad (3.13) $$

where $C$ is independent of $h$ and $p$. Since $p+1 > s$ in all the test cases, $\nu$ has been replaced by $s$. The constant $\alpha$ is equal to $s - \frac{1}{2}$ in the asymptotic range but we have seen that it can be equal to $s$ in a pre-asymptotic region. As a measure of the cost of the refinement, a functional $\mathcal{J}$ is defined as:

$$ \left\{ \begin{array}{ll} \mathcal{J}_{p \to p+1,h} = \frac{|\log(\|e_h^{p+1}\|_{L^2(\mathcal{D})}) - \log(\|e_h^p\|_{L^2(\mathcal{D})})|}{\log(\mathrm{dof}_{p+1,h}) - \log(\mathrm{dof}_{p,h})} & \text{for } p\text{-refinement} \\ \\ \mathcal{J}_{p,h \to h/2} = \frac{|\log(\|e_{h/2}^p\|_{L^2(\mathcal{D})}) - \log(\|e_h^p\|_{L^2(\mathcal{D})})|}{\log(\mathrm{dof}_{p,h/2}) - \log(\mathrm{dof}_{p,h})} & \text{for } h\text{-refinement} \end{array} \right. \quad (3.14) $$

The greater this quantity is, the better the strategy is. In particular, $p$-refinement is better than $h$-refinement if

$$ \mathcal{J}_{p \to p+1,h} > \mathcal{J}_{p,h \to h/2} \quad (3.15) $$
Using Eq. (3.13) in Eq. (3.15), $p$-refinement is the best strategy according to this functional if:

$$ \alpha \log \left( \frac{p+1}{p} \right) \geq s \log \left( \frac{p+2}{p+1} \right) \qquad (3.16) $$

If $\alpha = s$, Eq. (3.16) is always satisfied, so $p$-refinement should be always preferred. If $\alpha = s - \frac{1}{2}$, Eq. (3.16) becomes:

$$ s \geq \frac{1}{2} \frac{\log \left( 1 + \frac{1}{p} \right)}{\log \left( 1 + \frac{1}{p(p+2)} \right)} \qquad (3.17) $$

When $s = \frac{1}{2}$ or $\frac{3}{2}$, this relation is never satisfied (once $p$ is at least equal to 2), thus, $h$-refinement is the most interesting. A strategy to choose between both types of refinement with respect to this functional, is to perform $p$-refinement in the pre-asymptotic range then $h$-refinement. As we are in the pre-asymptotic part at least while $p \leq 4$, a basic $hp$-strategy, denoted $hp^B$ in the remainder of the manuscrit, consists in performing $p$-refinement while $p \leq 4$ and $h$-refinement afterwards.

**REMARK 3.3.** *This strategy is based on the global regularity $s$ and could be improved by using an estimation of the local regularity $s_K$. Therefore, performing $p$-refinement at the first steps of the refinement also leads to increase the polynomial order in the regions that are far from the singularities (where $s_K > s$). As the error is mostly located in the regions of singularities (where $s = s_K$), we can reasonably think that the use of $s$ instead of $s_K$ does not modify significantly the rate of convergence.*

**4. Presentation and derivation of the estimators.** In Section 3, we saw that the rate of convergence of the refinement methods depends on the regularity of the solution. We now present two estimators which exhibit different behaviours according to the regularity. The comparison of these estimators is interesting to design a refinement strategy.

**4.1. Radau estimator.** This way of approximating the difference between the exact and computed solution has been proposed in [2] and analysed in [16]. It is obtained by projecting the residual of Eq. (2.1) on a the so called Radau polynomial basis. Note that the analysis is only valid for analytical solutions.

**DEFINITION 4.1.** *The right Radau polynomial of order $k$ is denoted by $R_k$ and defined by*

$$ \left\{ \begin{array}{l} R_0(\xi) = L_0(\xi) \quad -1 \leq \xi \leq 1 \\ R_k(\xi) = L_k(\xi) - L_{k-1}(\xi) \quad k \geq 1, \quad -1 \leq \xi \leq 1, \end{array} \right. \qquad (4.1) $$

where $L_k$ is the Legendre polynomial of order $k$.

Similarly, the left Radau polynomial of order $k$ is

$$ \left\{ \begin{array}{l} R_0(\xi) = L_0(\xi) \quad -1 \leq \xi \leq 1 \\ R_k(\xi) = L_k(\xi) + L_{k-1}(\xi) \quad k \geq 1, \quad -1 \leq \xi \leq 1. \end{array} \right. \qquad (4.2) $$

If the elements $K$ are rectangles (or rectangle parallelepipeds if $d = 3$), there exists a one to one transformation from $K$ to a reference cell $\tilde{K} = [-1, 1]^d$. A variable
expressed in the cell K (resp. $\hat{K}$) is denoted $x = (x_1, ..., x_d)$ (resp. $\hat{x} = (\hat{x}_1, ..., \hat{x}_d)$). The Jacobian matrix associated to this transformation is denoted $J$, $J_{i,j}$ is the coefficient $i, j$ of the matrix and $|J|$ the determinant of this matrix. Using the hat-notation for variables expressed in the reference cell, Theorem 4.2 can be established:

**THEOREM 4.2.** [16] Let $u$ be the solution to Eq. (2.1) and $u_h^p$ the solution to Eq. (2.5). Then, the local error $e_h^p = u - u_h^p$ on each cell $K \in \mathcal{T}$ is given by

$$ \tilde{E}_K^R(x) = \sum_{k=p+1}^{\infty} h^k Q_k(x) \qquad (4.3) $$

with $Q_k \in \mathbb{Q}_{pK}$ and the polynomial of order $p+1$ is given by

$$ Q_{p+1}(x) = \sum_{i=1}^{d} \alpha_{i,K} R_{p+1}(x_i) \quad \text{with } x = (x_1, \dots, x_d), \qquad (4.4) $$

where the coefficients $\alpha_{i,K}$ are uniquely defined for all $i = [1, d]$, and $K \in \mathcal{T}$ by:

$$ \left( 2^d \Omega_i J_{i,i}^{-1} + 2^{d-1} c_{|K} \| R_{p+1} \|_{L^1(\hat{K})}^d |J| \right) \alpha_{i,K} = |J| \int_{\hat{K}} \left( \hat{g} - \Omega J^{-1} \cdot \nabla \hat{u}_h^p - c_{|K} \hat{u}_h^p \right) R_{p+1}(\hat{x}_i) d\hat{V} \quad (4.5) $$

*Proof.* As the exact solution is also solution of the weak form of the equation, the error $e_h^p = u_h^p - u$ satisfies the so-called orthogonality relation $B(e_h^p, v_h^p) = 0 \quad \forall v_h^p \in V_h^p$. A Taylor expansion of $e_h^p$ is then performed to prove the result. The interested reader is refered to [16] for the proof. □

**REMARK 4.3.** Right-Radau polynomials are used if $\Omega_i > 0$ and left ones otherwise.

An error estimator called Radau estimator $E^K$ in the remainder is obtained by truncating the series Eq. (4.3) and keeping only the terms of dominant order. Thus, on each cell $K \in \mathcal{T}$ an approximation of the error $\|e_h^p\|_{L^1(K)}$ is given by

$$ E_K^R = \left\| \sum_{i=1}^{d} \alpha_{i,K} R_{p+1}(x_i) \right\|_{L^1([-1,1]^d)} \qquad (4.6) $$

where the $\alpha_{i,K}$ are defined by Eq. (4.5).

**4.2. Estimator with limited regularity.** As we mentioned above, the Radau estimator is only valid for smooth solutions i.e. where the local regularity $s_K$ is at least bigger than $p_K + 1$. Since the regularity of the solution of problem Eq. (2.1) cannot be higher than $H^{3/2}$ [26], another estimator assuming less regularity has to be found. Error estimates on the space of functions with bounded variations (BV) for hyperbolic problems were first derived within the framework of the finite volume (FV) method [12]. The FV discretized transport equation is written

$$ h_K c u_h |_K + \sum_{\gamma \in \partial K} h_\gamma \nu_\gamma \cdot F_\gamma(\Omega, u_h |_K) = h_K g \quad \forall K \in \mathcal{T} \qquad (4.7) $$
where $F$ is defined by Eq. (2.4).

**THEOREM 4.4.** Let $u$ be the solution to Eq. (2.1) and $u_h$ the solution to Eq. (4.7). Then the error in a ball $B_R(x_0)$ of radius $R$ around the point $x_0$ satisfies:

$$ \|u_h - u\|_{L^1(B_R(x_0))} \le C\sqrt{\eta} $$

where $\tilde{\eta}$ is given by

$$ \tilde{\eta} = \sum_{K \in \mathcal{I}} \tilde{\eta}_K = \sum_{K \in \mathcal{I}} \sum_{\gamma \in \partial K} h_{\gamma} |\Omega \cdot \nu| |u_{h,+} - u_{h,-}| \quad (4.8) $$

with $\mathcal{I} = \{K \in \mathcal{T} | K \cap B_R(x_0) \neq \emptyset\}$.

*Proof.* The error estimate of [12] is obtained for a transient equation discretized by an implicit scheme. In order to adapt the proof to the stationary transport equation with boundary conditions, the main idea is to treat one of the spatial unknowns as the time unknown of the Cauchy problem and the boundary condition as the initial condition. The detailed adaptation of this proof can be found in [14]. □

The generalization to the discontinuous Galerkin discretization was proposed in [10]. In the framework of the transport equation, it may be stated as follows:

**THEOREM 4.5 (DG error estimator).** Let $u$ be the solution of Eq. (2.1) and $u_h^p$ the solution of Eq. (2.5). Then the error in a ball $B_R(x_0)$ of radius $R$ around the point $x_0$ satisfies:

$$ \|u_h^p - u\|_{L^1(B_R(x_0))} \le \eta_h = \sqrt{K_1\eta_1} + \sqrt{K_2\eta_2} $$

where $\eta_i = \sum_{K \in \mathcal{I}} \eta_{i,K}$, for $i = 1, 2$ ($\mathcal{I}$ defined in Theorem 4.4) and the local contributions $\eta_{i,K}$ are given by

$$ \eta_{1,K} = \eta_{1,K}^1 + \eta_{1,K}^2 = \int_K h_K |\Omega \cdot \nabla u_h^p + cu_h^p - g| \\ + \sum_{\gamma \in \partial K} h_{\gamma} \int_{\gamma} |\Omega \cdot \nu| |u_{h,+}^p - u_{h,-}^p| \quad (4.9) $$

$$ \eta_{2,K} = \| \overline{u_h^p} - u_h^p \|_{L^\infty(K)} \int_K |\Omega \cdot \nabla u_h^p + cu_h^p - g| \\ + \sum_{\gamma \in \partial K} \| \overline{u_h^p} - u_h^p \|_{L^\infty(\gamma)} \int_{\gamma} |\Omega \cdot \nu| |u_{h,+}^p - u_{h,-}^p| \quad (4.10) $$

where $\overline{u}_h^p$ denotes the mean value of $u_h^p$ over $K$. We will denote $\eta_i^j = \sum_{K \in \mathcal{I}} \eta_{i,K}^j$.

*Proof.* The interested reader is referred to [10] for the proof in the case of the Cauchy problem and to [14] for its adaptation to the framework of the stationary transport equation. □

**5. Analysis of the estimators.** Pitkäranta [27] proves that the solution of the neutron transport equation only belongs to $H^{\frac{3}{2}}$ and to $H^{\frac{1}{2}}$ if void embedded regions are present in the core configuration. The MMS benchmarks presented in Section 3 reproduce such regularities. The behaviour of the Radau estimator $E^R = (E_K^R)_{K \in \mathcal{T}}$ and the DG $E^{DG} = (E_K^{DG})_{K \in \mathcal{T}}$ on both benchmarks gives some indications about the
refinement strategy. All comparisons are done with respect to the final computational time and the number of degrees of freedom (dof).

Denoting $E_K$ an estimator of the error on the cell $K$, the cell $K$ is refined if

$$E_K \geq a \max_{L \in \mathcal{T}} E_L = \epsilon_r \quad (5.1)$$

where $a$ is a given constant belonging to $[0, 1]$. The choice of $a = 0$ corresponds to the uniform refinement.

**REMARK 5.1.** In the formulation of the neutron transport equation which is angularly discretized by the $S_N$ method [21], Eq. (2.1) is in fact solved for a set of directions $(\Omega_n)_n$. In the following, calculations are done on angular-independent quantities. For example, if $(u_n)_n$ denotes the set of solutions of Eq. (2.1) in the directions $(\Omega_n)_n$, $u$ is defined as $u = \sum_n w_n u_n$ where $w_n$ are the weights associated to the angular quadrature in the direction $\Omega_n$. Similar formulas are used for the different estimators and errors. In the tests presented here, an angular quadrature with 4 directions is used; the discrete angular directions are the diagonals of the domain and all the quadrature weights are equal.

**5.1. Towards a simplified estimator.** In the estimator obtained in Theorem 4.5, four different terms are present. Two of them have to be computed on the cell and the other two on the cell boundaries. To analyse their impact, the residual $\eta_1^1$ and the jump-term $\eta_1^2$ defined by Eqs. (4.9, 4.10) are compared. The differences between $\eta_1$ and $\eta_2$ are also analysed by studying the rate of convergence of $\|\overline{u_h^p} - u_h^p\|_{L^\infty(\mathcal{D})}$ with respect to $h$. The comparison of the convergence rate of each term can be used to determine if some terms can be neglected or to find a strategy to choose between $h$- and $p$-refinement.

<table><caption>TABLE 5.1</caption><thead><tr><th></th><th>MMS0</th><th>MMS1</th></tr></thead><tbody><tr><th>η<sub>1</sub><sup>1</sup></th><td>5/2</td><td>3</td></tr><tr><th>η<sub>1</sub><sup>2</sup></th><td>1/2</td><td>1</td></tr></tbody></table>

Rates of convergence of the error estimator terms $\eta_1^1$ and $\eta_1^2$ with respect to $h$ on the MMS0 and MMS1 benchmarks

Figure 5.1 represents the convergence of $\eta_1^1 = \sum_K \eta_{1,K}^1$ and $\eta_1^2 = \sum_K \eta_{1,K}^2$ as a function of $h$ and $p$ on the MMS0 and MMS1 benchmarks. Independently of the regularity, $\eta_1^1$ is negligible compared to $\eta_1^2$. The result of the polynomial interpolation presented in Table 5.1 shows that $\eta_1^1 \approx h^2\eta_1^2$.

These results seem to indicate that the terms on the edges are dominant compared to the ones on the cells. Therefore, the estimator can be simplified by keeping only $\eta_{1,K}^2$ and $\eta_{2,K}^2$. Let us now compare these last two terms.

Table 5.2 shows that $\|\overline{u_h^p} - u_h^p\|_{L^\infty(\mathcal{D})}$ does not converge when $u$ only belongs to $H^{1/2}(\mathcal{D})$. On the other hand, if $u \in H^{3/2}(\mathcal{D})$, a data fit based on the numerical experience presented in Figure 5.2 leads to

$$\ln(\|\overline{u_h^p} - u_h^p\|_{L^\infty(\mathcal{D})}) \approx -\alpha \ln(\mathrm{dof}) + \beta \quad (5.2)$$
<figure><img src="image_6.png" /><figcaption>FIGURE 5.1. Comparison of $\eta_1^1$ and $\eta_1^2$ as a function of dof and p on the MMS0 (left) and MMS1 (right) benchmarks</figcaption></figure>

with $\alpha \approx 0.5$. It implies that there exists $C_p \ge 0$ only depending on p such that

$$ \|\overline{u_h^p} - u_h^p\|_{L^\infty(D)} \approx C_p h \quad (5.3) $$

REMARK 5.2. Eq. (5.3) comes from a data fit of the MMS1 benchmark; such a relation can be justified by a Taylor expansion.:

$$ u_h^p = \overline{u_h^p} + \mathcal{O}(h) \quad (5.4) $$

Thus, $\eta_2$ is the dominant term compared to $\eta_1$ if $u \in H^{\frac{1}{2}}$ and both terms are equivalent if $u \in H^{\frac{3}{2}}$.

<table><caption>TABLE 5.2
$\|\overline{u_h^p} - u_h^p\|_{L^\infty(D)}$ as a function of h for the MMS0 benchmark</caption><thead><tr><th>h</th><th>2<sup>-1</sup></th><th>2<sup>-2</sup></th><th>2<sup>-3</sup></th><th>2<sup>-4</sup></th><th>2<sup>-5</sup></th><th>2<sup>-6</sup></th><th>2<sup>-7</sup></th></tr></thead><tbody><tr><td>||u<sub>h</sub><sup>p</sup> - u<sub>h</sub><sup>p</sup>||<sub>L<sup>∞</sup>(D)</sub></td><td>0.4316</td><td>0.4546</td><td>0.4661</td><td>0.4718</td><td>0.4746</td><td>0.4760</td><td>0.4767</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>FIGURE 5.2. $\|\overline{u_h^p} - u_h^p\|_{L^\infty(D)}$ as a function of dof for the MMS1 benchmark</figcaption></figure>

<caption>FIGURE 5.3. $L^2$-error on the MMS0 benchmark depending on the kept terms in $EDG$</caption>
<figure><img src="image_8.png" /><figcaption>FIGURE 5.4. $\eta_{1,K}^2$ (left) and $\eta_{1,K}^2 + \eta_{2,K}^2$ (right) on the MMS0 benchmark</figcaption></figure>

Hence, it seems necessary to keep both terms if the solution is discontinuous. Let us then compare the $L^2$-norm of the error with respect to dof by refining using as an error estimator $\eta_{1,K}^2 + \eta_{2,K}^2$ or only $\eta_{2,K}^2$ (Figure 5.3). A cell $K$ is refined in $h$ if it satisfies Eq. (5.1) with $\alpha = 0.5$. Figure 5.3 shows that the rate of convergence is the same with both estimators. Actually, $\|\overline{u_h^p - u_h^p}\|_{L^\infty(K)}$ has the same behaviour as the jumps between the cells. Comparing the values of the estimator by using $\eta_{1,K}^2$ (Figure 5.4 (left)) or $\eta_{1,K}^2 + \eta_{2,K}^2$ (Figure 5.4 (right)) shows that the same cells will be refined if only one term is kept or both. Therefore, since the computation of the infinite norm is costly, it seems reasonable to use only $\eta_{1,K}^2$ for the local error estimator.

The numerical comparison of the different terms contributing to $E^{DG}$ justifies the use of a simplified estimator $E^{FV}$ defined by

$$E^{FV} = \sum_{K \in \mathcal{I}} \eta_{1,K}^2 = \sum_{K \in \mathcal{I}} \sum_{\gamma \in \partial K} h_{\gamma} \int_{\gamma} |\Omega \cdot \nu| |u_{h,+}^p - u_{h,-}^p| \quad (5.5)$$

to lead the refinement process. Note that this quantity is the finite volume estimator defined in Theorem 4.4. This choice is also in agreement with the refinement process proposed in [29] for the neutron transport equation.

**REMARK 5.3.** The quality of an estimator is generally measured by its effectivity index, ratio between the error estimator and the real error in a given norm. Ideally, this ratio should be close to 1 or at least tend to 1 as the mesh size tends to 0 or the polynomial order to infinity. It is the case for the Radau estimator when the solution is smooth [16]. $E^{FV}$ just gives us an upper bound of the real error and thus the effectivity index is not necessarily close to 1 in this case.

**5.2. Towards a hp-refinement strategy.** Without a priori knowledge of the behaviour of the solution in the whole domain, a criterion to choose between h- or p-refinement is generally related to an estimated regularity of the solution [6, 23]. The behaviour of $\frac{1}{h_K} \|\overline{u_h^p - u_h^p}\|_{L^\infty(K)}$ presented in Figure 5.5 for the MMS0 benchmark might give us some information. Indeed, the dominant term is located along the characteristics and does not converge (see Table 5.2). Therefore, a possible criterion for h-refinement versus p-refinement could be to select for h-refinement the cells where this term is greater than a given parameter. Unfortunately, if the solution is
continuous, as in Figure 5.6 for the MMS1 benchmark, the ratio $\frac{\|\overline{u_h^p} - u_h^p\|_{L^\infty(K)}}{h_K}$ does not depend on $h_K$ and does not give additional information to locate irregular regions. Thus, this quantity is not a problem-independent criterion to choose between both types of refinement.

<figure><img src="image_9.png" /><figcaption>FIGURE 5.5. $\frac{\|\overline{u_h^p} - u_h^p\|_{L^\infty(K)}}{h_K}$ with $h = \frac{1}{8}, \frac{1}{16}, \frac{1}{64}$ for the MMS0 benchmark</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 5.6. $\frac{\|\overline{u_h^p} - u_h^p\|_{L^\infty(K)}}{h_K}$ with $h = \frac{1}{8}, \frac{1}{16}, \frac{1}{64}$ for the MMS1 benchmark</figcaption></figure>

Another idea consists in combining both estimators $E^{FV}$ and $E^R$. Actually, if the solution is regular $E^R$ should be a good approximation of the error and its behaviour be similar to that of $E^{FV}$. Both estimators should behave differently in less regular regions. In order to check this conjecture and deduce a hp-criterion, the behaviour of both estimators is studied on the two MMS benchmarks and compared to the local $L^1$-error norm $\epsilon_{L^1(K)}$ between the calculated solution $u_h^p$ and the reference solution $u$ locally defined by:

$$ \epsilon_{L^1(K)} = \int_K |u - u_h^p| dx \qquad (5.6) $$

Starting from a uniform mesh with a constant polynomial basis order equal to 1, the adaptive refinement process is performed with $E^{FV}$. A cell $K$ is refined according to Eq. (5.1) with $a = 0.9$ in order to slowly refine the mesh. At each step, the estimators and the real errors are compared with the values of $E^R$ and $E^{FV}$ normalized by the maximum of each estimator over the whole spatial domain. All estimators and errors are normalized by their maximum, since at each refinement step a cell $K^*$ is refined if $E_{K^*} \ge \alpha \max_{K \in \mathcal{T}} E_K$ so, if

$$ \frac{E_{K^*}}{\max_{K \in \mathcal{T}} E_K} \ge \alpha. $$
Thus, the normalized errors allow to visualize which cells would be refined and to compare the estimators to an optimal strategy (refinement leads by the real error).

**5.2.1. Comparison on the MMS0 Benchmark.** The behaviour of the two estimators and the real error is presented in Figures 5.7 to 5.9 for the initial mesh and after 6 refinement steps in Figures 5.10 to 5.12.

<figure><img src="image_11.png" /><figcaption>FIGURE 5.7. $E^{\text{FV}}$ (64 cells)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 5.8. $E^R$ (64 cells)</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 5.9. $\epsilon_{L^1}$ (64 cells)</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 5.10. $E^{\text{FV}}$ (208 cells)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIGURE 5.11. $E^R$ (208 cells)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 5.12. $\epsilon_{L^1}$ (208 cells)</figcaption></figure>

The error $\epsilon_{L^1}$ is most important in discontinuous regions i.e. along the diagonals (Figures 5.9 and 5.12). The estimators are expected to show the same behaviour and locate the discontinuities in order to refine these parts. The refinement process according to $E^{\text{FV}}$ begins by refining the corners (cells with the larger error in Figure 5.7) and then progresses along the diagonals (Figure 5.10). On the opposite, the behaviour of the refinement process according to $E^R$ (Eq. (4.3)) is not correct. Figures 5.8 and 5.11 show that the most important errors are located just above and below the characteristics but not along these lines. Leading the refinement with such an estimator for a discontinuous solution gives an unadapted final mesh as shown in Figure 5.14 contrary to the one obtained with $E^{\text{FV}}$ (Figure 5.13).

<figure><img src="image_17.png" /><figcaption>FIGURE 5.13. Mesh obtained when refining with $E^{\text{FV}}$ (724 cells)</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>FIGURE 5.14. Mesh obtained when refining with $E^R$ (724 cells)</figcaption></figure>
This first benchmark gives us some indications in order to perform hp-refinement. In the cells where the solution is discontinuous, cells are tagged to be refined by $E^{FV}$ but not by $E^R$. In this case, h-refinement has to be done. If cells are only tagged by $E^R$, e.g. the cells located above and below the diagonals, the relative error (normalized by the maximum) is probably overestimated because of the poor approximation of the maximum. In this case, the cell is not modified at this step.

**5.2.2. Comparison on the MMS1 Benchmark.** The estimators and the real error are compared in Figures 5.15 to 5.20. Let us recall that in this benchmark, the solution is continuous but with a non-continuous first derivative.

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /><figcaption>FIGURE 5.15. $E^{FV}$ (256 cells)</figcaption></figure>

<caption>FIGURE 5.16. $E^R$ (256 cells)</caption>

<figure><img src="image_24.png" /><figcaption>FIGURE 5.17. $\epsilon_{L1}$ (256 cells)</figcaption></figure>

<figure><img src="image_25.png" /></figure>

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /><figcaption>FIGURE 5.18. $E^{FV}$ (460 cells)</figcaption></figure>

<caption>FIGURE 5.19. $E^R$ (460 cells)</caption>

<caption>FIGURE 5.20. $\epsilon_{L1}$ (460 cells)</caption>

Even if it is only the first derivative of the solution that is discontinuous, the maximum error is still located along the characteristic lines (Figures 5.17 and 5.20). The behaviours of the two estimators are closer to one another than in the MMS0 benchmark. Particularly for the conforming mesh, it is difficult to distinguish between the two estimators (Figures 5.15 and 5.16). It means that the solution is sufficiently regular to make $E^R$ and $E^{FV}$ equivalent. In this case, p-refinement could be performed.

Figures 5.18 and 5.19 lead to the same conclusion than MMS0. As $E^R$ underestimates the error in some parts, it could lead to refine some cells where the real error is not important and fail to compute solutions to the desired level of accuracy. Actually, due to non-conformities, the error on the diagonal cells is underestimated, as the error is normalized by the poorly approximated maximum, it leads to an artificial increase of the error in regions where the solution is regular. If the mesh presented in Figure 5.19 is refined according to $E^R$, it leads to select several cells outside the diagonal lines where the real error $\epsilon_{L1}$ is low. On the opposite, a refinement according to $E^{FV}$ chooses the cells on the characteristic lines where $\epsilon_{L1}$ is the most important.

The analysis of these two benchmarks confirms the choice of $E^{FV}$ to lead the refinement process. The estimator $E^R$ can be used to give additional information about the regularity of the solution in a given region when compared to $E^{FV}$. The resulting strategy, denoted by $hp^{2E}$, is summarized in Table 5.3 where $\epsilon_r$ denotes a
given criterion defined by Eq. (5.1): *p*-refinement is performed if $E^R$ and $E^{FV}$ have the same behaviour and *h*-refinement otherwise.

<table><caption>TABLE 5.3
First refinement strategy (hp$^{2E}$-strategy)</caption><thead><tr><th></th><th><i>E</i><sup><i>R</i></sup> &gt; <i>&epsilon;</i><sub><i>r</i></sub></th><th><i>E</i><sup><i>R</i></sup> &lt; <i>&epsilon;</i><sub><i>r</i></sub></th></tr></thead><tbody><tr><th><i>E</i><sup><i>FV</i></sup> &gt; <i>&epsilon;</i><sub><i>r</i></sub></th><td><i>p</i>-refinement</td><td><i>h</i>-refinement</td></tr></tbody></table>

**6. Comparison of hp-refinement methods.** We saw in Section 3 that a *p*-refinement yields a convergence in $O(p^{-s})$ for the first values of *p*, followed by the asymptotic rate in $O(p^{-s+1/2})$. A possible way to take benefit of the preasymptotic region is to modify the strategy proposed in Table 5.3 in the following way: while $p \le 4$, *p*-refinement is still performed on tagged cells; if $p > 4$, the hp$^{2E}$-strategy is used. This new method is denoted hp$^{2E}_{mod}$ and summarized in Table 6.1.

<table><caption>TABLE 6.1
Second refinement strategy (hp$^{2E}_{mod}$-strategy)</caption><thead><tr><th></th><th><i>E</i><sup><i>R</i></sup> &gt; <i>&epsilon;</i><sub><i>r</i></sub></th><th><i>E</i><sup><i>R</i></sup> &lt; <i>&epsilon;</i><sub><i>r</i></sub> and <i>p</i> &le; 4</th><th><i>E</i><sup><i>R</i></sup> &lt; <i>&epsilon;</i><sub><i>r</i></sub> and <i>p</i> &gt; 4</th></tr></thead><tbody><tr><th><i>E</i><sup><i>FV</i></sup> &gt; <i>&epsilon;</i><sub><i>r</i></sub></th><td><i>p</i>-refinement</td><td><i>p</i>-refinement</td><td><i>h</i>-refinement</td></tr></tbody></table>

The strategies presented in Table 5.3 and 6.1 are implemented and tested on the MMS0 and MMS1 benchmarks. The tests on the two benchmarks yield the same conclusions, and therefore only the MMS0 results are presented. The representation of $\epsilon_{L2}$ as a function of dof for the different refinement strategies is given in Figure 6.1.

<figure><img src="image_28.png" /><figcaption>FIGURE 6.1. L$^2$-error versus dof for the MMS0 benchmark with different refinement strategies</figcaption></figure>

Two regions can be distinguished. While dof $\le 2000$, *p*-refinement is the best strategy. As explained in Section 3, it is more efficient to increase *p* and therefore decrease $C_p$ at the beginning than to refine the mesh. Surprisingly, both proposed hp-strategies are less efficient than *p*-refinement. With the hp$^{2E}$-strategy, since the solution is discontinuous, *h*-refinement is often performed. For the MMS0 benchmark, it is the case from the points 0 to 1 (referring to the numbered points in Figure 6.1) and beyond point 2; *p*-refinement is done only between points 1 and 2.
<figure><img src="image_29.png" /><figcaption>FIGURE 6.2. $L^2$-error versus time for the MMS0 benchmark with different refinement strategies</figcaption></figure>

As detailed in Section 3, even if the solution is non-regular, $p$-refinement converges faster for small values of $p$. The improvement obtained with the $hp_{mod}^{2E}$-algorithm is clear. Up to point 1, only $p$-refinement is performed, leading to the best convergence rate. Then the refinement progresses by alternating $h$-refinement (from 1 to 2, 3 to 4, 5 to 6 with the same convergence slope than the exclusive $h$-methods) and $p$-refinement (from 2 to 3 and 4 to 5), leading to the same accuracy as the exclusive $p$-refinement method for dof ≈ 2000. Thus, the order of the polynomial basis and therefore the computational cost has been reduced without damaging the convergence.

The polynomial basis orders and meshes obtained after some refinement steps for both $hp^{2E}$ and $hp_{mod}^{2E}$ strategies are presented in Figure 6.3. The comparison of both estimators allows to find the singular regions and, thus, the diagonal cells are refined with $h$-refinement for the $hp^{2E}$ strategy. Using the $hp_{mod}^{2E}$ strategy allows to perform $p$-refinement before using the $h$-refinement and thus to improve the convergence rate as shown in Figure 6.1.

<figure><img src="image_30.png" /><figcaption>FIGURE 6.3. Polynomial basis orders and meshes for the $hp^{2E}$ and $hp_{mod}^{2E}$ strategies</figcaption></figure>

The analysis becomes different for dof larger than 2000. The $p$-refinement method converges slower, becoming less competitive than all the other strategies. On the opposite, $h$-refinement methods have a constant convergence rate which is
equal to the regularity of the solution. If such strategies are less interesting for a small number of degrees of freedom because of the decay of $C_p$ with $p$, they unconditionally ensure the convergence. A combination of the $h$ and $p$ strategies in the $hp_{mod}^{2E}$-method is really efficient. Actually, we combine the advantages of $p$-refinement at the beginning to decrease the value of $C_p$ and then ensure a constant rate of decay. This analysis remains valid regarding the calculation time (Figure 6.2). Indeed, problems due to large values of $p$ are enhanced and the performance of $hp_{mod}^{2E}$ method is confirmed. The interested reader is referred to [15] for further illustrations of the relationship between time and dof.

A comparison is done on a real core configuration to highlight the interest of $hp$-methods but also to compare the two-estimators based strategy with other $hp$-approaches. The type-parameter method ($hp^{TP}$) proposed in [1] performs $h$-refinement if the ratio between the estimator at the current order $p$ over its value at order ($p-1$) is lower than a parameter $\gamma$ (0.6 in [1]) and $p$-refinement otherwise. This ratio measures the efficiency of $p$-refinement at the previous step. Note that two values are required to begin this strategy, so $p$-refinement is always performed at the first refinement step. The last strategy used for the comparison is a basic strategy presented Section 3.3. A summary of the different methods is proposed in Table 6.2.

<table><caption>TABLE 6.2
*Overview of the different hp-refinement strategies*</caption><thead><tr><th>Strategy name</th><th><i>h</i>-refinement</th><th><i>p</i>-refinement</th></tr></thead><tbody><tr><td>Two-estimators (<i>hp</i><sup>2E</sup>)</td><td><i>E</i><sup><i>R</i></sup>(<i>K</i>) ≤ <i>γ</i></td><td><i>E</i><sup><i>R</i></sup>(<i>K</i>) > <i>γ</i></td></tr><tr><td>Modified two-estimators (<i>hp</i><sub>mod</sub><sup>2E</sup>)</td><td><i>E</i><sup><i>R</i></sup>(<i>K</i>) ≤ <i>γ</i> and <i>p</i><sub>K</sub> > 4</td><td><i>E</i><sup><i>R</i></sup>(<i>K</i>) > <i>γ</i> or <i>p</i><sub>K</sub> ≤ 4</td></tr><tr><td>Type-parameter (<i>hp</i><sup><i>TP</i></sup>)</td><td><sup><i>E</i><sub>p</sub><sup><i>TV</i></sup>(<i>K</i>)</sup>&frasl;<sub><i>E</i><sub>p-1</sub><sup><i>TV</i></sup>(<i>K</i>)</sub> ≤ <i>γ</i></td><td><sup><i>E</i><sub>p</sub><sup><i>TV</i></sup>(<i>K</i>)</sup>&frasl;<sub><i>E</i><sub>p-1</sub><sup><i>TV</i></sup>(<i>K</i>)</sub> > <i>γ</i></td></tr><tr><td>Basic (<i>hp</i><sup><i>B</i></sup>)</td><td><i>p</i><sub>K</sub> > 4</td><td><i>p</i><sub>K</sub> ≤ 4</td></tr></tbody></table>

<figure><img src="image_31.png" /><figcaption>FIGURE 6.4. Presentation of the ZONA2B benchmark and its solution for a given energy</figcaption></figure>

The ZONA2B benchmark is described in [13] and the convergence is studied in [15] with respect to dof and to the computing time for $h$- and $p$-refinement. Figure 6.4 presents the geometry of this benchmark. It contains three different media: FUEL where the neutrons are produced, NASS a reflector part and STELL a shield made of steel. As for the MMS benchmarks, the non-regular parts are located at the interfaces.
<figure><img src="image_32.png" /><figcaption>FIGURE 6.5. $L^2$-error versus dof for the ZONA2B benchmark with different refinement strategies</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>FIGURE 6.6. $L^2$-error versus time for the ZONA2B benchmark with different refinement strategies</figcaption></figure>

Figures 6.5 and 6.6 compare the different hp-strategies presented in Table 6.2. As in the previous analysis for the MMS benchmarks two regions have to be distinguished. In the asymptotic range, the rate of convergence of all the methods is the same. The differences are due to the preasymptotic part. As detailed in the previous sections, performing p-refinement is the best strategy while p is small. It explains the superiority of $hp_{mod}^{2E}$ compared to $hp^{2E}$. The situation is different for $hp^{TP}$. This estimator is based on a direct measure of the efficiency of the p-refinement. As p is efficient at the beginning, $hp^{TP}$ favours p-refinement for small values of p. In realistic benchmarks, this kind of method seems better than an estimation of the regularity. Note however that the rate of convergence is dependent of the value of the parameter $\gamma$. With a good choice of $\gamma$, $hp^{TP}$ and $hp_{mod}^{2E}$ are equivalent. The most surprising aspect is that a similar convergence is even possible with the basic strategy. The importance of the preasymptotic range is at the origin of such a behaviour that leads to the conclusion that on realistic benchmarks, the criterion of selection between both types of refinement does not modify significantly the rate of convergence. Performing p-refinement while $p \le 4$ then h-refinement ($hp^B$) seems a sufficient strategy for the neutron transport equation.

**7. Conclusion.** This article studies some aspects related to the development of hp-refinement methods for the transport equation in neutronics. First, theoretical rates of convergence are recalled showing the importance of the regularity of the solution in the rates of convergence. Two estimators with different behaviour with respect to the regularity are then derived, analysed and compared in order to find a hp-refinement strategy. This method is improved by taking into account the fastest convergence of p-refinement in the preasymptotic range and then compared to other strategies found in the literature. Finally, we show numerically that in the framework of the neutron transport equation, a simple strategy selecting p-refinement up to a given value ($p = 4$ in this study) then h-refinement seems sufficient to obtain an interesting convergence with respect to dof and to the computing time.

**Acknowledgements.** The authors are grateful to the reviewers for their careful reading and for their useful remarks and suggestions.
[1] S. Adjerid, M. Aiffa, and J. E. Flaherty. *Computational Methods for Singularly Perturbed Systems*. AMS, Providence, 1998.

[2] S. Adjerid and T. C. Massey. A posteriori discontinuous finite element error estimation for two-dimensional hyperbolic problems. *Comp. Meth. in Appl. Mech. and Eng.*, 191:5877–5897, 2002.

[3] M. Ainsworth and D. Kay. The approximation theory for the p-version finite element method and application to non-linear elliptic PDEs. *Numer. Math.*, 82:351–388, 1999.

[4] M. Ainsworth and B. Senior. An adaptive refinement strategy for hp-finite element computations. *Applied Numerical Mathematics*, 26:165–178, 1998.

[5] I. Babuška and M. R. Dorr. Error estimates for the combined h and p-version of the finite element method. *Numer. Math.*, 37:257–277, 1981.

[6] I. Babuška and M. Suri. The hp-version of the finite element method with quasiuniform meshes. *M2AN Modél. Math. Anal. Numér.*, 21:199–238, 1987.

[7] I. Babuška, B. A. Szabo, and I. N. Katz. The p version of the finite element method. *SIAM J. Numer. Anal.*, 18:512–545, 1981.

[8] K. S. Bey and J. T. Oden. hp-version discontinuous Galerkin method for hyperbolic conservation laws. *Comput. Methods Appl. Mech. Engrg.*, 133:259–286, 1996.

[9] B. Cockburn, B. Dong, and J. Guzman. Optimal convergence of the original DG method for the transport-reaction equation on special meshes. *SIAM J. Numer. Anal.*, 46(3):1250–1265, 2008.

[10] A. Dedner and M. Ohlberger. *A New hp-Adaptive DG Scheme for Conservation Laws Based on Error Control*. Hyperbolic Problems; Theory, Numerics, Applications, 2008.

[11] A. Ern and J. L. Guermond. *Theory and practice of Finite elements*. Applied Mathematical Sciences 159 Springer, 2004.

[12] R. Eymard, T. Gallouët, M. Ghilani, and R. Herbin. Error estimates for the approximate solutions of a nonlinear hyperbolic equation given by finite volume schemes. *IMA J. Numeric. Anal.*, 18(4):563–594, 1998.

[13] P. J. Finck et al. The CIRANO experimental program in support of advanced fast reactor physics. In *Proc. Int. Conf. PHYSOR'96, Mito, Ibaraki, Japan*, volume 2, pages E66–E75, Sept. 16-20, 1996.

[14] D. Fournier. *Analyse et Développement de Méthodes de Raffinement hp en Espace pour l'Equation de Transport des Neutrons*. PhD thesis, Université de Provence, http://tel.archives-ouvertes.fr/tel-00769546, 2011.

[15] D. Fournier, P. Archier, R. Le Tellier, and C. Suteau. Improvement of neutronic calculations on a MASURCA core using adaptive mesh refinement capabilities. In *International Conference on Mathematics, Computational Methods & Reactor Physics, MEC 2011, Rio de Janeiro, RJ, Brazil*, 2011.

[16] D. Fournier, R. Le Tellier, and C. Suteau. Analysis of an a posteriori error estimator for the transport equation with Sn and discontinuous Galerkin discretizations. *Annals of Nuclear Energy*, 38(2-3):221–231, 2011.

[17] P. Houston, C. Schwab, and E. Süli. Stabilized hp-finite element method for first-order hyperbolic problems. Technical Report NA-98/14, Oxford University Computing Laboratory, 1998.

[18] P. Houston, C. Schwab, and E. Süli. Discontinuous hp-finite element methods for advection-diffusion problems. Technical Report NA-00/15, Oxford University Computing Laboratory, 2000.

[19] C. Johnson and J. Pitkäranta. An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation. *Math. Comp.*, 46:1–26, 1986.

[20] P. Lesaint and P. A. Raviart. *On a Finite Element Method for Solving the Neutron Transport Equation*. Mathematical Aspects of Finite Elements in Partial Differential Equations, C.A. deBoor, ed., Academic Press, New York, 1974.

[21] E. E. Lewis and F. Miller, Jr. *Computational Methods of Neutron Transport*. John Wiley and Sons, Inc., New York, NY, 1984.

[22] C. Lingus. Analytical test cases for neutron and radiation transport codes. In *Proc. Second Conference on Transport Theory*. Los Alamos, New Mexico, US Atomic Energy Commission, CONF-710017, 1971, January 26-29.

[23] C. Mavriplis. Adaptive mesh strategies for the spectral element method. *Comp. Meth. in Appl. Mech. and Eng*, 116:77–86, 1994.

[24] J. M. Melenk. On the robust exponential convergence of hp finite element methods for problems with boundary layers. *IMA J. Numeric. Anal.*, 17:577–601, 1997.

[25] J. M. Melenk and B. I. Wohlmuth. On residual-based a posteriori error estimation in hp-FEM. *Advances in Computational Mathematics*, 15:311–331, 2001.
[26] J. Pitkäranta. On the spatial differencing of the discrete ordinate neutron transport equation. *SIAM J. Numer. Anal.*, 15(5):856-869, 1978.

[27] J. Pitkäranta. Estimates for the derivatives of the solution to the weakly singular fredholm integral equations. *J. Math. Anal.*, 11(6):952-968, 1980.

[28] G. R. Richter. An optimal-order error estimate for the discontinuous Galerkin method. *Math. Comp.*, 50:75-88, 1988.

[29] Y. Wang, W. Bangerth, and J. Ragusa. Three-dimensional h-adaptivity for the multigroup neutron diffusion equations. *Progress in Nuclear Energy*, 51:543-555, 2009.

[30] Y. Wang and J. Ragusa. On the convergence of DGFEM applied to the discrete ordinates transport equation for structured and unstructured triangular meshes. *Nuclear Science and Engineering*, 163:56-72, 2009.