A phase-space discontinuous Galerkin scheme for the radiative transfer equation in slab geometry

: We derive and analyze a symmetric interior penalty discontinuous Galerkin scheme for the approximation of the second-order form of the radiative transfer equation in slab geometry. Using appropriate trace lemmas, the analysis can be carried out as for more standard elliptic problems. Supporting examples show the accuracy and stability of the method also numerically, for different polynomial degrees. For discretization, we employ quad-tree grids, which allow for local refinement in phase-space, and we show exemplary that adaptive methods can efficiently approximate discontinuous solutions. We investigate the behavior of hierarchical error estimators and error estimators based on local averaging.


Introduction
We consider the numerical solution of the radiative transfer equation in slab geometry, which has several applications such as atmospheric science [27], oceanography [5], pharmaceutical powders [9] or solid state lighting [38]; see also [10] for a recent introduction.
Due to the product structure of Ω, it seems natural to use separate discretization techniques for the spatial variable  and the angular variable .This is for instance done in the spherical harmonics method, in which a truncated Legendre polynomial expansion is employed to discretize  [18].The resulting coupled system of Legendre moments, which still depend on , is then discretized for instance by finite differences or finite elements [18].Another class of approximations consists of discrete ordinates methods which perform a collocation in  and the integral in eq. ( 7) is approximated by a quadrature rule [18].The resulting system of transport equations is then discretized for instance by finite differences [18] or discontinuous Galerkin methods [26,24], and also spatially adaptive schemes have been used [41].
A major drawback of the independent discretization of the two variables  and  is that a local refinement in phase-space is not possible.Such local refinement is generally necessary to achieve optimal schemes.For instance, the solution can be non-smooth in the two points (, ) = (0, 0) and (, ) = (, 0), which are exactly the two points separating the inflow from the outflow boundary.Although certain tensor-product grids can resolve this geometric singularity for the slab geometry, such as double Legendre expansions [18], they fail to do so for generic multi-dimensional situations.Moreover, local singularities of the solution due to singularities of the optical parameters or the source terms can in general not be resolved with optimal complexity.
Phase-space discretizations have been used successfully for radiative transfer in several applications, see, e.g., [15,35,36,37] for slab geometry, [32] for geometries with spherical symmetries, or [21,33] for more general geometries.Let us also refer to [31] for a phase-space discontinuous Galerkin method for the nonlinear Boltzmann equation.A non-tensor product discretization that combines ideas of discrete ordinates to discretize the angular variable with a discontinuous Petrov-Galerkin method to discretize the spatial variable has been developed in [13].
In this work, we aim to develop a numerical method for ( 7)-( 8) that allows for local mesh refinement in phase-space and that allows for a relatively simple analysis and implementation.To accomplish this, we base our discretization on a partition of Ω such that each element in that partition is the Cartesian product of two intervals.Local approximations are then constructed from products of polynomials defined on the respective intervals.In order to easily handle hanging nodes, which such partitions generally contain, we use globally discontinuous approximations.In case the resulting linear systems are very large, iterative solution techniques with small additional memory requirements may be employed for their numerical solution, such as the conjugate gradient method, which, however, requires the linear system to be symmetric positive definite.Therefore, we employ a symmetric interior penalty discontinuous Galerkin formulation.Besides the proper treatment of traces, which requires the inclusion of a weight function in our case, the analysis of the overall scheme is along the standard steps for the analysis of discontinuous Galerkin methods [16].As a result, we obtain a scheme that enjoys an abstract quasi-best approximation property in a mesh-dependent energy norm.Our choice of meshes also allows to explicitly estimate the constants in auxiliary tools, such as inverse estimates and discrete trace inequalities.As a result, we can give an explicit lower bound on the penalty parameter required for discrete stability.This lower bound for the penalty parameter depends only on the polynomial degree for the approximation in the -variable and is relatively simple to compute; see [20] for the estimation of the penalty parameter in the context of standard elliptic problems.Our theoretical results about accuracy and stability of the method are confirmed by numerical examples, which show optimal convergence rates for different polynomial degrees assuming sufficient regularity of the solution.Moreover, we show that adaptively refined grids are able to efficiently construct approximations to non-smooth solutions.
For the local adaptation of the grid we investigate several error estimators.First, we consider two hierarchical error estimators, which either use polynomials of higher degree or the discrete solution on a uniformly refined mesh, respectively.Such estimators have been investigated in the elliptic context, e.g., in [7,30].Our numerical results show that these error indicators can be used to refine the mesh towards the singularity of the solution.A drawback of these estimators is that an additional global problem has to be solved in every step.Since the solutions to ( 7)-( 8) can be discontinuous in , the proofs developed for elliptic equations to show that the global estimator is equivalent to a locally computable quantity, see, e.g., [30], do not apply.To overcome the computational complexity of building estimators that require to solve a global problem, we propose an a posteriori estimator based on a local averaging procedure.This cheap estimator shows a similar performance compared to the more expensive hierarchical ones mentioned before.
The outline of the rest of the manuscript is as follows.In Section 2 we introduce notation and collect technical tools, such as trace theorems.In Section 3 we derive and analyze the discontinuous Galerkin scheme.Section 4 presents numerical examples confirming the theoretical results of Section 3. Section 5 shows that our scheme works well with adaptively refined grids.We introduce here two hierarchical error estimators and one based on local post-processing.The paper closes with some conclusions in Section 6.

Preliminaries
We denote by  2 (Ω) the usual Hilbert space of square integrable functions and denote the corresponding inner product by Furthermore, we introduce the Hilbert space which consists of square integrable functions for which the weighted derivative is also square integrable; see [2, Section 2.2].We endow the space  with the graph norm To treat the boundary condition eq. ( 8), let us introduce the following inner product
If  ∈  2 (Ω), testing eq. ( 10) with functions in  ∞ 0 (Ω) shows that      has a weak   -derivative in  2 (Ω) and eq.( 7) holds a.e. in Ω.In particular,      ∈  and      has a trace.For  ∈  , an integration by parts in eq. ( 10) then shows that Since the trace operator is surjective from  to  2 (Γ; ) [2, Theorem 2.9], it follows that eq. ( 8) holds in  2 (Γ; ).We denote the space of solutions with data  ∈  2 (Ω) and  ∈  2 (Γ; ) by 3 Discontinuous Galerkin scheme In the following we will derive the numerical scheme to approximate solutions to eq. ( 10).After introducing a suitable partition of Ω using quad-tree grids and corresponding broken polynomial spaces, we can essentially follow the standard procedure for elliptic problems, cf.[16].One notable difference is that we need to incorporate the weight function  on the faces.

Mesh and broken polynomial spaces
In order to simplify the presentation, and subsequently the implementation, we consider quad-tree meshes [23] as follows.Let  be a partition of Ω such that   is constant on each element  ∈  , and that Next, let us introduce some standard notation.Denote P  the space of polynomials of one real variable of degree  ≥ 0, and let the broken polynomial space  ℎ be denoted by with   ,   ≥ 0. Here, P +1 ⊗P  denotes the tensor product of P +1 and P  .Moreover, let  (ℎ) :=  + ℎ .By ℱ  ℎ  we denote the set of interior vertical faces, that is for any  ∈ ℱ  ℎ  there exist two disjoint elements we define the jump and the average of  ∈  ℎ by In order to take into account local variations in the mesh size and diffusion coefficient 1/  , we furthermore define the dimensionless quantity where ℎ  ,  ∈ {1, 2}, denotes the local mesh size of the element   in -direction.For an interior face , which is shared by two elements    ∈  ,  = 1, 2, as above, let us introduce the sub-elements We note that the inclusion in eq. ( 16) can be strict in the case of hanging nodes, see for instance Figure 1.
Proof.Using Lemma 1 we have that Using eq. ( 17), we estimate the weighted derivative term as follows Using that ℎ  =    −    and  ≤ 1, we thus obtain that which concludes the proof.
Remark 1.The value of   () of the inverse inequality in eq. ( 17) can be computed by solving a small eigenvalue problem of dimension  + 1, which is obtained by transforming eq. ( 17) to the unit interval.In fact,   () is the maximal eigenvalue of where of the space of polynomials of degree at most  on the unit interval.Explicit bounds for   (), which are optimal for  = 1, 2, are given in [12].

Derivation of the DG scheme
In order to extend the bilinear form defined in eq. ( 11) to the broken space  ℎ , we denote with  ℎ  the broken derivative operator such that In view of eq. ( 10), let us then introduce the bilinear form which is defined on  (ℎ).Note that   and   ℎ coincide on  .In order to obtain a consistent bilinear form,   ℎ needs to be modified.We follow [16,Chapter 4] to determine the required modification.Choosing  ∈  *ℎ :=  * +  ℎ and  ℎ ∈  ℎ , integration by parts in  shows that where we used the identity with  , defined in eq. ( 15) and with positive penalty parameter   > 0, which will be specified below.
Since  = 0 on any  ∈ ℱ  ℎ  and  ∈  , it follows that  ℎ is consistent, i.e., for  ∈  * it holds that The discrete variational problem is formulated as follows: Find  ℎ ∈  ℎ such that

Analysis
For the analysis of ( 20), let us introduce mesh-dependent norms In order to show discrete stability and boundedness of  ℎ , we will use the following auxiliary lemma.
Proof.By definition of the average, we have that where  1 ,  2 and  1  ,  2  denote the restrictions of  and   to  1  and  2  , respectively.To estimate the first integral on the right-hand side, we employ the Cauchy-Schwarz inequality to obtain where we used Lemma 2 applied to      1 , which is a piecewise polynomial of degree   in .A similar estimate holds for the second integral.Using the Cauchy-Schwarz inequality, we then obtain that which, in view of eq. ( 15), concludes the proof.
The auxiliary lemma allows to bound the consistency terms in  ℎ , which gives discrete stability of  ℎ .
Lemma 4 (Discrete stability).For any  ∈  ℎ it holds that Proof.Let  ℎ ∈  ℎ , and consider Using Lemma 3, and the fact that each sub-element    touches at most two interior vertical faces, an application of the Cauchy-Schwarz yields for any  > 0, Hence, by choosing  = 1/2, from which we obtain the assertion.
Proof.The space  ℎ is finite-dimensional.Hence, Lemma 4 implies the assertion.
To proceed with an abstract error estimate, we need the following boundedness result.
Proof.We have that The first two terms can be estimated using the Cauchy-Schwarz inequality as follows For the third term, we use Lemma 3 to obtain ‖  ‖  2 ( ;) .
To separate the terms that include  and , respectively, we apply the Cauchy-Schwarz inequality once more and use again that each sub-element    touches at most two interior faces, to arrive at Remark 4. In view of Remark 1, the value of   (  ) can be computed explicitly once   is known.Hence, we can give an explicit value for the penalization parameter   such that the discontinuous Galerkin scheme eq. ( 20) is well-posed and the error enjoys the bound given in Theorem 2. We note that we choose here   to be the same for all interior faces.Moreover, the choice of   is independent of the partition  and the mean-free path 1/  ; while the dependence on the mean-free path is explicit through  , , which might be exploited if the behavior of the scheme is investigated in the diffusion limit where the mean-free path tends to zero.Let us refer to [24] for a detailed discussion about issues of DG schemes for radiative transfer in the diffusion limit.
Remark 5. Instead of using the symmetric bilinear form   ℎ to define  ℎ in eq. ( 18), we may use the more general bilinear form with parameter  ∈ [−1, 1], cf.[14].The choice  = 1 leads to   ℎ , while the choices  = 0 or  = −1 yield incomplete interior penalty and the non-symmetric interior penalty discontinuous Galerkin schemes, respectively, see also [28,42] for the case  = −1.We note that for  = −1, the terms involving the face integral vanish for  = , and hence coercivity can be proven straight-forward.However, since symmetry is lost, improved  2 -convergence rates for sufficiently smooth solutions do not hold in general, see [4] and Table 3 below.Moreover, the numerical solution of large non-symmetric linear systems can be more difficult than in the symmetric case.We mention that the results in this section can be extended with minor modifications to the general case −1 ≤  ≤ 1.

Numerical examples
In the following we confirm the theoretical statements about stability and convergence of Section 3 numerically.Let   = 1/2 and   = 1 and the width of the slab be  = 1.We then define the source terms  and  in ( 7)-( 8) such that the exact solution is given by the following function Here,  {>1/2} () denotes the indicator function of the interval (1/2, 1), i.e.,  is discontinuous in  = 1/2, but note that  ∈  * .We compute the DG solution  ℎ of ( 20) on a sequence of uniformly refined meshes, initially consisting of 16 elements, see Figure 1.Hence, the discontinuity in  is resolved by the mesh.
For our computations we use the spaces  ℎ with   and   in eq. ( 14), that is piecewise polynomials of degree   in  and piecewise polynomials of degree   + 1 in .The value of   (  ) of the inverse inequality in eq. ( 17) is computed numerically by solving a small eigenvalue problem of dimension   + 1, see Remark 1.
For the numerical solution of the resulting linear systems, we a usual fixed-point iteration [1]: Introducing the auxiliary bilinear form  ℎ (, ) =  ℎ (, ) − (   , ), the fixed-point iteration maps   ℎ to  +1 ℎ by solving The fixed-point iteration converges linearly with a rate   /  [1], which is bounded by 1/2 in this example.The iteration is stopped as soon as ‖ +1 ℎ −   ℎ ‖  2 (Ω) < 10 −10 .For acceleration of the source iteration by preconditioning see also [1,39].The matrix representation of  ℎ has a block structure for the uniformly refined meshes considered in this section, and its inverse can be applied efficiently via LU factorization.
Table 1 shows the  ℎ norm of the error  −  ℎ between the exact and the numerical solution for  =   =   .For fixed , we observe a convergence rate of  + 1 under mesh refinement, which is expected from the smoothness of  per element and Remark 3. In particular, inspecting Table 1 by rows, we notice linear convergence for  = 0, quadratic convergence for  = 1, and so on.
Tab. 1: Error ‖ −  ℎ ‖  ℎ for different local polynomial degrees with  =  = , see eq. ( 14), and uniformly refined meshes with  elements and solution  defined in eq.(22).Since the coefficients are smooth, we may expect higher order convergence in the  2 -norm for the symmetric formulation if   =   + 1, see Remark 5.In Table 2 and Table 3 we compare the symmetric interior penalty method ( = 1) with its non-symmetric counterpart ( = −1), with  introduced in Remark 5. Table 2 shows that, for fixed   , the  2 -error decays upon mesh refinement at an improved rate of (ℎ +1 ) for the symmetric interior penalty method.This improved convergence rate can also be observed for the non-symmetric interior penalty method if the employed polynomial degrees are odd, while the suboptimal rate (ℎ  ) can be observed if the used polynomial degrees are even, cf.[29] for a similar observation on the convergence rates for the unsymmetric interior penalty method in the context of non-stationary convection diffusion problems.

Adaptivity
In this section we show, by examples, that hierarchical error estimators, see, e.g., [30] for the elliptic case, as well as estimators based on averaging the approximate solutions are a possible choice to adaptively construct optimal partitions  of Ω to approximate non-smooth solutions to eq. ( 7).In particular, we show how adaptive mesh refinement is beneficial if the discontinuity of the solution is not resolved by the mesh.
To highlight the dependency on the partition  of Ω and on the polynomial degree, we will write   , instead of  ℎ for the solution of the discontinuous Galerkin scheme (20).Similarly, assuming  :=   =   , we write   , for the corresponding approximation space instead of  ℎ , see eq. ( 14).Let  ′ be another partition of Ω such that  ⊂  ′ , and let  ′ ≥ .Denoting ‖ • ‖ some norm defined on  +   , +   ′ , ′ , and supposing the saturation assumption, which has been used, e.g., also in [7], for some universal constant  < 1, we obtain the equivalence between the approximation error and the estimator  :=   ′ , ′ −   , , i.e., ( For a justification of the saturation assumption in the context of elliptic problems we refer to [11].In the following numerical experiments, we use the norm ‖•‖  , defined as to investigate the behaviour of two hierarchical error indicators for different test cases.The local error contributions are then given by where  ∈  .The mesh is then refined by a Dörfler marking strategy [17,43], that is all elements in the set  ⊂  are refined, where  ⊂  is the set of smallest cardinality such that where 0 <  ≤ 1 is the bulk-chasing parameter.Differently from the previous section, we assume   = 1 and   = 0.Moreover, we consider two different manufactured solution  1 and  2 given by The choice of 1/ √ 2 in the indicator function ensures that the corresponding line discontinuity of  2 is never resolved by our mesh.Furthermore,    1 is bounded and vanishes in (0, 0).In particular, we note that  1 ,  2 ∈  * .In the following we report on numerical examples using the Dörfler parameter  := 0.75.We note that we obtained similar results for the choice  = 0.3.

Hierarchical p-error estimator
Setting  ′ :=  and  ′ :=  + 1, the hierarchical -error estimator is defined as We note that   := 1/2 +   ( + 1) is used for the stabilization parameter to obtain both numerical solutions   , and   ,+1 .Figure 2 and Figure 3 show the convergence rates for adaptively refined meshes using the   indicator, for different values of the polynomial degree .We observe that for the manufactured solution  1 , which has a point singularity in the origin, the indicator follows tightly the curve of the actual error.Moreover, the error decays at the optimal rate 1/ √  +1 , with  denoting the number of degrees of freedom in   , , also shown for comparison.
For the manufactured solution  2 with line discontinuity defined in eq. ( 27), the convergence behavior is different.For  = 0, the error and the error estimator stay rather close, and follow the curve for the optimal rate.For  ≥ 1 the rate is sub-optimal, which is expected from a counting argument.Moreover, also the error estimator is not as close to the true error anymore, compared to the test case with  1 .

Hierarchical h-error estimator
Using once again the test cases eq. ( 27) and eq.( 26), we now keep  ′ :=  fixed and we construct  ′ by uniform refinement of  , i.e., every element in  is subdivided in 4 new elements with halved edge length.The error estimator is now Some comments are due for the computation of   ′ , , for which we use, as in definition (18) but with  replaced by  ′ , the bilinear form  ℎ ′ :   ′ , ×   ′ , → R. Since   , is a subspace of   ′ , , we require, similar to [30], that  ℎ is the restriction of  ℎ ′ to   , , in the sense that Fig. 2: Test case eq. ( 26) with singularity in (0, 0).Broken  1 norm of the approximation error and of the -estimator plotted against the theoretical optimal rate, for different values of the starting polynomial degree  = 0, 1, 2, 3, in a double logarithmic scale.26) with singularity in (0, 0).Broken  1 norm of the approximation error and of the ℎ-estimator plotted against the theoretical optimal rate, for different values of the polynomial degree  = 0, 1, 2, 3, in a double logarithmic scale.
Comparing the penalty terms in  ℎ and  ℎ ′ shows that the above restriction is fulfilled if   = 2 ′  .Since we need to ensure discrete stability of both bilinear forms, we choose  ′  := 1/2 +   ().Figure 4 and Figure 5 show the optimality of the estimator for the manufactured solution  1 with point singularity defined in eq. ( 26), and sub-optimality in the case of discontinuous exact solutions eq. ( 27), except for  = 0, where a similar comment as for the -hierarchical estimator applies.We note that in all cases, the estimator is close to the actual error.Moreover, as shown in Figure 6 we observe that the estimator is able to detect the line discontinuity present in  2 .We note that the condition in eq. ( 30) is not essential for the results shown in this section.In fact, similar results can be obtained by using   =   ′ = 1/2 +   ().The condition eq. ( 30) is required in the next section.

Error estimator based on the solution of local problems
Since the computation of the global error estimators  presented in Section 5.1 and Section 5.2 is in general expensive, we investigate also an error estimator based on the solution of local problems, as presented in [22,30] for corresponding elliptic problems.In this approach, the computed solution   , is understood as the coarse-mesh approximation to some function, here   ′ , .Instead of using  =   , −   ′ , as before, local approximations   are computed element-wise.27) with line discontinuity.Broken  1 norm of the approximation error and of the ℎ-estimator plotted against the theoretical optimal rate, for different values of the polynomial degree  = 0, 1, 2, 3, in a double logarithmic scale.
At this point we observe that ( 20), ( 30) and ( 34) imply, for all  ∈  ℎ , Eventually we introduce the functions {  ∈   ′ , ()| ∈  ℎ } as solutions to the local problems Each   can be computed independently of each other.The function  = ∑︀ ∈   ∈   ′ , then may serve as an approximation to the estimator  =   ′ , −   , on   ′ , , as in [30] for elliptic problems.Following [30,Theorem 4.1], we can prove a lower bound for  in terms of the local error estimator , i.e.Lemma 6.We have that Proof.We first rewrite (36) in terms of the estimator: Plugging  =   ∈   ′ , () into the previous equation and recalling that  = ∑︀    , we have where we used Corollary 1 in the last step.Coercivity of   ℎ ′ , expressed in (32), and eq.( 31) entail Combining ( 39) and ( 40) we have eventually which concludes the proof.
Due to the lack of appropriate interpolation operators for functions in the space  , one cannot adapt the proofs given in [30] in a straight-forward fashion to show a bound of  in terms of the local contributions .
In fact, some preliminary numerical tests, based on the broken  1 -norm eq. ( 24), suggest that the desired equivalence of  and  might not be true.In a similar spirit, our preliminary numerical tests indicate that standard residual error estimators are either not reliable or efficient, which, again, may be explained by the lack of suitable interpolation error estimates required to obtain the correct scaling in terms of the local mesh size of the different local contributions, cf., [16,Section 5.6] or [3,43].Therefore, we investigate in the next section another error estimator based on local averages.

Error estimator based on averaging the approximate solution
In the context of a posteriori error estimation and adaptive mesh refinement, ZZ-error estimators named after Zienkiewicz and Zhu [44] are widely used in practice.Compared to the previously mentioned hierarchical error estimators, their major advantage is the fact that no further mesh nor a further solution is required.We consider the case   =   = 0.In order to obtain a reliable error estimator, one simply takes a discontinuous  ℎ ∈  ℎ and approximates it by some continuous piecewise linear polynomial ũℎ by a post-processing step.
In the presence of a geometrically conforming triangulation, such a continuous piecewise polynomials ũℎ can be described as a linear combination of the well-known Lagrange nodal basis functions.However, our approximation involves hanging nodes and we therefore restrict the construction to the set of regular nodes  ℎ , i.e.
If a regular node  ∈  ℎ is shared by four quadrilaterals  1 , ...,  4 of the same area, the idea is to set the value of a continuous polynomial to where the local contributions are used to refine the mesh using Dörfler marking as described above.Figure 7 shows the convergence rates for adaptively refined meshes using the averaging indicator for the test eq.( 27).The indicator behaves correctly and replicates the curve of the actual error.These curves have the same slope as the optimal rate 1/ √  curve, with  number of elements in the quad-tree mesh, also shown for comparison.
In comparison to the hierarchical estimators, cf. Figure 3 and Figure 5, the averaging error estimator follows the actual error curve more closely.

Conclusions and discussion
We developed and analyzed a discontinuous Galerkin approximation for the radiative transfer equation in slab geometry.The use of quad-tree grids allowed for a relatively simple analysis with similar arguments as for more standard elliptic problems.While such grids allow for local mesh refinement in phase-space, the implementation of the numerical scheme is straightforward.For sufficiently regular solutions, we showed optimal rates of convergence.We showed by example that non-smooth solutions can be approximated well by adaptively refined grids.The ability to easily adapt the computational mesh can also be useful when complicated geometries must be resolved, which may occur in higher-dimensional situations.Also more general elements could be employed at the expense of a more complicated notation and analysis; we leave this to future research.In order to automate the mesh adaptation procedure, an error estimator is required.We investigated numerically hierarchical error estimators and estimators based on local averaging in a post-processing step.All three estimators closely follow the actual error, and, in the case of point singularities, they can be used to obtain optimal convergence rates.We note that the hierarchical error estimators require to solve global problems, and it is left for future research to investigate whether a localization is possible.Upper bounds for the error can be derived for consistent approximations using duality theory [25].Rigorous a posteriori error estimation has also been done using discontinuous Petrov-Galerkin discretizations [13].We leave it to future research to analyze the error estimators for the discontinuous Galerkin scheme considered here and to generalize the present method to a corresponding ℎ −  version, where the polynomial degrees can be varied independently over the elements.
While the solution of the linear systems for uniformly refined grids can be implemented using the established preconditioned iterative solvers [1,39], the structure of the linear systems for adaptively refined grids is more complex because the equations do not fully decouple in ; compare the situations in Figure 1.One possible direction is to develop nested solvers, or to adapt the methodology of [40].We leave this for future research.
Another direction of future research entails the regularity of the right-hand side  in eq. ( 7) and  in eq. ( 8).If  and  define only an element in the dual space of  , see eq. ( 10), then the flux  −1     / ∈  in general, and thus the flux may not have a trace.In this low regularity regime, the analysis of Section 4 cannot be carried out.A possible remedy might be to use a lifting operator to replace the face integral by integrals over Ω, see, e.g., [16, p. 138] or [34].

Fig. 3 : 2 .
Fig. 3: Test case eq.(27) with discontinuity along the line  = 1/ √ 2. Broken  1 norm of the approximation error and of the -estimator plotted against the theoretical optimal rate, for different values of the starting polynomial degree  = 0, 1, 2, 3, in a double logarithmic scale.

Fig. 4 :
Fig. 4: Test case eq.(26) with singularity in (0, 0).Broken  1 norm of the approximation error and of the ℎ-estimator plotted against the theoretical optimal rate, for different values of the polynomial degree  = 0, 1, 2, 3, in a double logarithmic scale.

Fig. 5 :
Fig. 5: Test case eq.(27) with line discontinuity.Broken  1 norm of the approximation error and of the ℎ-estimator plotted against the theoretical optimal rate, for different values of the polynomial degree  = 0, 1, 2, 3, in a double logarithmic scale.

Fig. 7 :
Fig.7: Non-smooth test case eq.(27).Top Locally refined mesh with the average error estimator after 6 (left) and 9 (right) refinements.Bottom: Convergence of the DG solution and the averaging error estimator on adaptively refined grids as well as the optimal rate 1/ √  (light dotted line) in a double logarithmic scale.The dashed line with o shows the behavior of the  2 -error using the averaging error estimator for grid adaptation.The solid line with x shows the corresponding values of the averaging error estimator.For comparison, we show convergence of the  2 -error (dash-dotted with x), where the grid adaptation is based on the  2 -error itself.The dotted line with o shows the values of the corresponding averaging error estimator on that grid.

2 ,
which concludes the proof as   +   ≥ 3/2.Before continuing, an inspection of the previous proof shows that we have the following corollary stating boundedness of  ℎ on  ℎ .
Taking into account the possibility of quadrilaterals of different area, for a node  ∈  ℎ , we define by   the union of all elements of  ∈  sharing the vertex .