# A subcell-enriched Galerkin method for advection problems

In this work, we introduce a generalization of the enriched Galerkin (EG) method. The key feature of our scheme is an adaptive two-mesh approach that, in addition to the standard enrichment of a conforming finite element discretization via discontinuous degrees of freedom, allows to subdivide selected (e.g. troubled) mesh cells in a non-conforming fashion and to use further discontinuous enrichment on this finer submesh. We prove stability and sharp a priori error estimates for a linear advection equation by using a specially tailored projection and conducting some parts of a standard convergence analysis for both meshes. By allowing an arbitrary degree of enrichment on both, the coarse and the fine mesh (also including the case of no enrichment), our analysis technique is very general in the sense that our results cover the range from the standard continuous finite element method to the standard discontinuous Galerkin (DG) method with (or without) local subcell enrichment. Numerical experiments confirm our analytical results and indicate good robustness of the proposed method.

## Authors

• 9 publications
• 3 publications
• 3 publications
• ### Staggered DG method with small edges for Darcy flows in fractured porous media

In this paper, we present and analyze a staggered discontinuous Galerkin...
05/22/2020 ∙ by Lina Zhao, et al. ∙ 0

• ### A posteriori goal-oriented bounds for the Poisson problem using potential and equilibrated flux reconstructions: application to the hybridizable discontinuous Galerkin method

We present a general framework to compute upper and lower bounds for lin...
06/21/2021 ∙ by Nuria Pares, et al. ∙ 0

• ### Fast formation and assembly of isogeometric Galerkin matrices for trimmed patches

This work explores the application of the fast assembly and formation st...
01/20/2021 ∙ by Benjamin Marussig, et al. ∙ 0

• ### On Approximating Discontinuous Solutions of PDEs by Adaptive Finite Elements

For singularly perturbed problems with a small diffusion, when the trans...
07/08/2019 ∙ by Shun Zhang, et al. ∙ 0

• ### Adaptive non-hierarchical Galerkin methods for parabolic problems with application to moving mesh and virtual element methods

We present a posteriori error estimates for inconsistent and non-hierarc...
05/12/2020 ∙ by Andrea Cangiani, et al. ∙ 0

• ### Low-memory, discrete ordinates, discontinuous Galerkin methods for radiative transport

The discrete ordinates discontinuous Galerkin (S_N-DG) method is a well-...
07/01/2019 ∙ by Zheng Sun, et al. ∙ 0

• ### Mesh sampling and weighting for the hyperreduction of nonlinear Petrov-Galerkin reduced-order models with local reduced-order bases

The energy-conserving sampling and weighting (ECSW) method is a hyperred...
08/06/2020 ∙ by Sebastian Grimberg, et al. ∙ 0

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1. Introduction

The main idea of the enriched Galerkin (EG) method is to extend the approximation space of the continuous finite elements by including some element-local discontinuous functions and to utilize a solution procedure similar to that of the discontinuous Galerkin (DG) method (Riemann solvers, edge fluxes, …). The latter feature makes the EG schemes fundamentally different from the XFEM methods that frequently also rely on local approximation space enrichments. The resulting discretization is locally conservative and robust but, in multidimensions, has substantially fewer degrees of freedom than a DG method of the same order.

In [RuppL2020], the EG methods were re-cast as a generalization of the classical finite elements, i.e. continuous Galerkin (CG) methods by considering the EG space as a combination of arbitrary continuous and discontinuous Galerkin (DG) test and trial spaces. However, the original EG scheme proposed in [Becker2003] for the advection equation was a combination of lowest order finite elements and finite volumes discretized using the DG framework. This methodology was further developed and investigated by Wheeler, Lee, and coworkers, who also considered higher order enriched CG methods and a wider range of applications [lee2016locally, sun2009locally, LEE2019, lee2018phase, lee2017adaptive, lee2018enriched, choo2018enriched, Kadeethum2019, kadeethum2019comparison]. The analysis of EG method in [RuppL2020] used a special EG-type projection and was limited to elliptic and parabolic problems. Nonetheless, it paved the way to the analysis for hyperbolic equations in this work.

Similarly to CG approximations, EG methods for hyperbolic equations may develop spurious oscillations. Kuzmin et al. [KuzminHR2020] proposed several algebraic flux correction schemes to ensure the validity of local maximum principles. Limiting techniques of this kind have also been successfully applied to CG [Kuzmin2012, Lohmann2019] and DG [FrankRK2020] discretizations. The use of localized subcell limiters was found to be essential in extensions to high-order Bernstein finite elements [HajdukEtAl2020a, HajdukEtAl2020b, KuzminMQL2020, Lohmann2019]. An -adaptive approach to subcell limiting was introduced in [KuzminEtAl2019]. Using continuous blending functions, a high-order finite element approximation on a large macrocell was combined with a bound-preserving piecewise (multi-)linear subcell approximation.

Another well-known class of methods relying on subcell limiting to suppress spurious oscillations has been introduced in [Dumbser2014] and generalized to unstructured meshes in [Dumbser2016]. These techniques are based on the ADER-DG schemes proposed in [Dumbser2008] and possess a very attractive capability to detect high- and low-regularity solution behavior. The underlying a posteriori limiting strategy was inspired by the Multi-dimensional Optimal Order Detection (MOOD) approach originally developed for finite volumes. In the context of the ADER-DG methods, physical and numerical admissibility conditions are enforced by, first, advancing the solution in time using a high-order DG method on the coarse mesh, and, for troubled cells, repeating the last time step locally via a low-order DG (i.e. finite volume) method on the submesh.

Our subcell EG method has the potential to further customize the local approximation space by supporting the whole range of local polynomial orders on both, the coarse and the fine (subcell) mesh. This feature of our approach makes it possible to combine popular - and -adaptivity techniques with the two-mesh approach , while exploiting its intrinsic ability to assess the local solution regularity.

The main purpose of this work is to present a stability and a priori error analysis for the subcell-enriched EG method for the linear advection equation and to demonstrate the performance of the new scheme using some test problems. As in [RuppL2020] this analysis is conducted in a unified framework that covers the CG, DG, and EG (with and without subcell enrichment) discretizations. The implementation of the new numerical scheme was carried out in our FESTUNG framework [FrankRAK2015, ReuterAWFK2016, JaustRASK2018, ReuterHRFAK2019, ReuterRAFK2020] based on our EG scheme for the shallow-water equations [HauckAFHR2020].

### 1.1. Model problem

We consider a non-stationary advection equation on a bounded Lipschitz domain (with ). The precise formulation of the linear hyperbolic problem to be solved is as follows:

 ∂tu+∇⋅(a(t,x)u)=f(t,x) in (0,T)×Ω, (1.1)

for a given velocity field and a right-hand side function . Additionally, initial data is prescribed, and we denote by the outward unit normal to . Furthermore, we assume that the inflow boundary

 Γ− \coloneqq {x∈∂Ω:a(t,x)⋅νΩ<0}

is independent of time and disjointly subdivided into Dirichlet and flux boundaries (this subdivision is also assumed to be independent of time), i.e.

 u=uD on ΓD and |a⋅νΩ|u=gF on ΓF,uD∈L2(0,T;H1/2(ΓD)),gF∈L2(0,T;H1/2(ΓF)).

For the sake of simplicity, we assume that there exists such that .

### 1.2. Structure of the manuscript

The remainder of this manuscript is structured as follows: In Section 2, we introduce the enriched Galerkin method with local subcell enrichment for advection equations. Section 3 investigates the energy stability of the new scheme, while its a priori convergence is proved in Section 4 and verified numerically in Section 5. A short conclusions section wraps up the article.

## 2. The enriched Galerkin finite element method

### 2.1. Basic definitions and notations

In the following, denotes a successively refined family of ( is the number of elements) of -dimensional non-overlapping partitions of (see [PietroErn, Def. 1.12]) that is assumed to be regular (in the sense of [PietroErn, Def. 1.38]) and geometrically conformal (in the sense of [ErnGuermond, Def. 1.55]). For the sake of simplicity, we assume that consists of simplices and/or quadrilaterals/hexahedrons.

Furthermore, denotes a mesh of which some elements have been refined (Fig. 2.1 (middle)). The mesh can be geometrically non-conformal. By construction, can be embedded into a regular and conformal mesh (Fig. 2.1 (right)), which contains the elements added during the refinement process. Hence, we can write as disjoint union

denoting the subsets of unrefined and refined elements, respectively. Writing for the set of faces we define the skeleton of as

 Σ\coloneqq⋃F∈F(TH|h)F.

We write for the diameter of ; furthermore, parameter refers to the maximum diameter of an element of a mesh, i.e., . If without an index is evaluated on a face, a unit normal with respect to the face is arbitrarily chosen.

The double mesh sequence is called weakly quasi-uniform if there exists a constant such that for all , all , , and all , we have

 hH|ˆK := max{hK:K∈Sh% with K⊂ˆK} ≤ ρmin{hK:K∈Sh with K⊂ˆK}.

To simplify notation we set for .

The test and trial spaces for our EG method utilize the broken polynomial spaces of order on some mesh . They are denoted by and consist of element-wise polynomials of degree at most

(simplices) or tensor-product polynomials of degree at most

in each spatial coordinate (quadrilaterals/hexahedrons) without any continuity constraints. Thus,

 Vkℓ,m\coloneqq(Pk(TH)∩C(Ω))+Pℓ(TH)+Pm(TH|h)

for , . Here, , and one can observe that is the standard continuous finite element space. Obviously .

In this work, we utilize several types of projection/interpolation operators denoted as follows:

• and are the -projections into the spaces and , respectively.

• is the standard interpolation operator for finite element space .

• is the mapping used to project the initial data into proposed in [RuppL2020] and given by

 π:L2(Ω)∩C(Ω)→Vkℓ,m,πu\coloneqqIkHu+Πℓ,mH|h(u−IkHu). (2.1)

### 2.2. Semi-discrete formulation

The semi-discrete EG formulation of the problem can be constructed by using the standard DG bilinear and linear forms for the advection equation on . The bilinear form uses the notion of averages and jumps , which for with are defined as

 {[g]}\coloneqq12(g|K−+g|K+),[[g]]\coloneqqg|K−νK−+g|K+νK+,

for a scalar that is element-wise smooth enough to have traces. On , this definition is modified as follows:

 {[g]}\coloneqqg,[[g]]\coloneqqgνΩ,

Hence, the jump turns a scalar into a vector. Also note the following property of jumps used in our analysis

 [[g2]]=2{[g]}[[g]].

Given a velocity field , we define the upwind value of as

 (g)↑a\coloneqq{[g]}a+sign(a⋅ν)2[[g]]⋅νa,

where is the standard signum function.

Using this notation, we can formulate our semi-discrete problem

 a(U,φ)=b(φ)

with trial function and test function from for almost every and , where

 a(U,φ)\coloneqq ∫Ω∂tUφdx−∑K∈TH|h∫KUa⋅∇φdx+∫Σ∖Γ−(U)↑a⋅[[φ]]dσ, b(φ)\coloneqq ∫Ωfφdx+∫ΓDuD|a⋅νΩ|φdσ+∫ΓFgFφdσ.

Note that is a standard DG bilinear form; its consistency implies that the EG bilinear form is also consistent since .

## 3. Stability analysis

The stability of the method can be obtained exactly as the stability of the DG methods. Thus,

###### Theorem 3.1.

The EG solution is stable.

###### Proof.

We test with , use the identity and integrate by parts to obtain

 12∂t∥U ∥2L2(Ω)+12∫Ω(∇⋅a)U2dx− 12∫Σ[[U2]]⋅adσ+∫Σ∖Γ−(U)↑a⋅[[U]]dσ=12∫Σ|a⋅νΩ|[[U]]2dσ = ∫ΩfUdx+∫ΓDuD% |a⋅νΩ|Udσ+∫ΓFgFUdσ ≤ ∫ΩfUdx+12∫ΓD[|a⋅νΩ|u2D+|a⋅νΩ|[[U]]2]dσ+12∫ΓF[g2F|a⋅νΩ|+|a⋅νΩ|[[U]]2]dσ,

where the last inequality follows from the Young’s and Cauchy–Schwarz inequalities and uses the assumption on . This directly implies the -stability without exponential growth of constants if and after integrating with respect to time. Otherwise Grönwall’s, Young’s, and Cauchy–Schwarz inequalities give the result (after moving to the right-hand side). ∎

## 4. Error analysis

For the error analysis, we need some auxiliary results:

###### Lemma 4.1.

The operator of (2.1) is an orthogonal projection into with respect to the -inner product, i.e.,

 ∫Ω(πu−u)φdx=0∀φ∈Pℓ(TH)+Pm(TH|h). (4.1)
###### Proof.

Result follows directly from the fact that and the -orthogonality of . ∎

###### Lemma 4.2 (Best approximation property of Πℓ,mH|h).

For all , , and all ,

 ∥Πℓ,mH|hg−g∥L2(K)≤∥φ−g∥L2(K) (4.2)
###### Proof.

Follows directly from the -orthogonality of and the possibility to localize the projection to all . ∎

###### Lemma 4.3 (Inverse inequality).

Let be a regular mesh sequence. There exists a constant such that for all , all , and all

 |φ|H1(K)≤C h−1K∥φ∥L2(K). (4.3)
###### Proof.

This is [PietroErn, Lem. 1.44]. ∎

###### Lemma 4.4 (Discrete trace inequality).

Let be a regular mesh sequence. There exists a constant such that for all , all , all , and all with

 ∥φ∥L2(F)≤C h−1/2K∥φ∥L2(K). (4.4)
###### Proof.

This is [PietroErn, Lem. 1.46]. ∎

###### Lemma 4.5 (Continuous trace inequality).

Let be a regular mesh sequence. There exists a constant such that for all , all , all , and all with

 ∥v∥2L2(F)≤C (|v|H1(K)+h−1K∥v∥L2(K))∥v∥L2(K). (4.5)
###### Proof.

This is [PietroErn, Lem. 1.49]. ∎

###### Lemma 4.6 (Approximation property).

Let be a regular mesh sequence. There exists a constant such that for all , all , and all

 |ΠkHv−v|Hm(K) ≤C hk+1−mK|v|Hk+1(K)for k∈N∪{0}, (4.6) |IkHv−v|Hm(K) ≤C hk+1−mK|v|Hk+1(K)for k∈N∖{0}. (4.7)
###### Proof.

The first inequality is [KnabnerAngermann, Theo. 3.29]. The second inequality is [PietroErn, Lem. 1.58]. ∎

###### Lemma 4.7.

Let be a regular mesh sequence. There exists a constant such that for all , all , and all

 |Π0Hv∥L∞(K) ≤∥v∥L∞(K), (4.8) ∥Π0Hv−v∥L∞(K) ≤hK|v|W1,∞(K). (4.9)
###### Proof.

The first inequality is the observation that is the element-wise mean of which needs to be smaller than or equal to its essential maximum. The second inequality is a simple combination of [KnabnerAngermann, Theo. 3.24 & 3.26]. ∎

Next, we formulate our main result.

###### Theorem 4.8.

Let be a weakly quasi-uniform mesh (double) sequence, and let , . Then, the EG approximation converges in to the analytical solution , i.e., there exists independent of and such that

 ∥u−U∥2L∞(0,T;L2(Ω))≤C∑K∈THh2˜mH|KH2(k−m)K|u|2H1(0,T;Hk+1(K))

with for simplicial meshes and or general meshes and (i.e. in the case of DG). Otherwise, .

###### Proof.

Defining

 eu\coloneqqU−πu and θu\coloneqqπu−u,

we have due to the consistency and since that

 a(eu+θu,eu)=0 for almost every t∈(0,T).

This can be rewritten as

 12∂t∥eu∥2L2(Ω)+ 12∫Σ|a⋅ν|[[eu]]2dσ=−12∫Ω(∇⋅a)e2udx\eqqcolonΞ1−∫Σ∖Γ−(θu)↑a⋅[[eu]]dσ\eqqcolonΞ2−∫Ω∂tθueudx\eqqcolonΞ3 +∑K∈TH|h∫Kθu(a−Π0H|ha)⋅∇eudx\eqqcolonΞ4+∑K∈TH|h∫KθuΠ0H|ha⋅∇eudx\eqqcolonΞ5.

We can immediately deduce that:

• If then holds, and this term can be moved to the left hand side and integrated into the energy norm. This is consistent with the continuous case, when the mass sinks lead to an increase in the stability.

• If is element-wise constant, then .

• If the mesh is simplicial, and the globally continuous polynomials are from the space , then . Using (4.1) yields , provided that .

• If the mesh is quadrilateral, and the globally continuous polynomials are from the space , then . Using (4.1) yields provided that , i.e., in the case of DG.

Next, we estimate terms :

 |Ξ1| ≤ 12 ∥∇⋅a∥L∞(Ω)∥eu∥2L2(Ω), |Ξ2| ≤ ∫Σ∖Γ−(θu)↑2dσ+14∫Σ|a⋅ν|[[eu]]2dσ ≤ C ∑K∈TH|h∫∂Kθ2udσ+14∫Σ|a⋅ν|[[eu]]2dσ C∑K∈TH|h(|θu|H1(K)+h−1K∥θu∥L2(K))∥θu∥L2(K)+14∫Σ|a⋅ν|[[eu]]2dσ ≤ C ∑K∈TH|hhK|θu|2H1(K)+ C∑K∈TH|hh−1K∥θu∥2L2(K)+14∫Σ|a⋅ν|[[eu]]2dσ, |Ξ3| ≤ 12 ∥∂tθu∥L2(Ω)+12 ∥eu∥L2(Ω), |Ξ4| ≤ ∑K∈TH|h∥θu∥L2(K)∥a−Π0H|ha∥L∞(K)∥∇eu∥L2(K) C ∑K∈TH|h∥θu∥L2(K)hK|a|W1,∞(K)h−1K∥eu∥L2(K) ≤ |a|2W1,∞(Ω)∑K∈TH|h∥θu∥2L2(K)+C∑K∈TH|h∥eu∥2L2(K), |Ξ5| ≤ ∑K∈TH|h∥θu∥L2(K)∥Π0H|ha∥L∞(K)∥∇eu∥L2(K) C ∑K∈TH|h∥θu∥L2(K)∥a∥L∞(K)h−1K∥eu∥L2(K) ≤ ∥a∥2L∞(Ω)∑K∈TH|hh−2K∥θu∥2L2(K)+C∑K∈TH|h∥eu∥2L2(K).

This would give the desired result (after applying Grönwall’s inequality – if needed) provided that we could find good estimates for the terms involving , and . Note that only the norm enters the exponential term in the Grönwall estimate.

We consider the cases and separately. In the first case, we can estimate

 =∑ˆK∈TH∖SHhrH|ˆK∥Πℓ,mH|h(u−IkHu)−(u−IkHu)∥2L2(ˆK) ∑ˆK∈TH∖SH hrH|ˆK∥ΠmH|h(u−IkHu)−(u−IkHu)∥2L2(ˆK) =ˆK∈TH∖SH∑Sh∋K⊂ˆKhrH|ˆK∥ΠmH|h(u−IkHu)−(u−IkHu)∥2L2(K) C ∑ˆK∈TH∖SHh2m+2+rH|ˆK |u−IkHu|2Hm+1(ˆK) C ∑ˆK∈TH∖SHh2m+2+rH|ˆKH2k−2mˆK |u|2Hk+1(ˆK).

In the second case, we obtain for using the same arguments

 ∑K∈SHHrK∥πu−u∥2L2(K)≤ C ∑K∈SHH2k+2+rK|u|2Hk+1(K).

The estimate for is conducted analogously. Here, the projection is used to obtain

 ˆK∈TH∖SH∑Sh∋K⊂ˆKhrK|πu−u|2H1(K) C ˆK∈TH∖SH∑Sh∋K⊂ˆK(hrK|π⋆u−u|2H1(K)+hr−2K∥πu−π⋆u∥2L2(K)),

which gives the needed estimate after inserting into the second summand and redoing the aforementioned arguments. Collecting all terms gives the result. ∎

###### Remark 4.9.

This result is not optimal, since it uses high regularity of the temporal derivative. However, in the case of DG, i.e. and , the proof can be streamlined by replacing (and ) by —this also implies that the initial data is constructed using an orthogonal projection with respect to the -norm. Here, also the distinction between simplices and quadrilaterals/hexahedrons becomes unnecessary, and the polynomial approximation spaces may all be of type. This results in and yields the optimal estimate

 ∥u−U∥2L∞(0,T;L2(Ω)) ≤ C∑K∈THh2m+1H|KH2(k−m)K|u|2L2(0,T;Hk+1(K)) ≤ CH2k+1|u|2L2(0,T;Hk+1(Ω)),

where is only assumed to be an element of .

## 5. Numerical results

### 5.1. Analytical convergence test

In order to verify the convergence of the numerical schemes, we use the method of manufactured solution. On the domain and the time interval , we define the analytical solution and velocity filed by

The right-hand side of the problem is chosen so that and satisfy (1.1). We prescribe Dirichlet boundary conditions on the inflow boundary, i.e., , and use and .

Let and denote the refinement levels for the meshes with element sizes and , respectively. The initial mesh (==1) consisting of four triangles is obtained by diagonally subdividing ; finer meshes are produced by connecting the edge midpoints of every triangle. As temporal discretization, we use an explicit SSP Runge–Kutta method with stages.

We utilize the EG method with polynomial orders and on the coarse grid (of refinement level ) enriched by the DG method of order at most on the fine grid (of refinement level ).

Our implementation currently supports the approximation orders up to two. This yields four possible combinations of , and . In Table 5.1, the -th entry corresponds to the -error at time using the refinement levels and .

We observe that our analytical convergence rates are confirmed by the numerical tests; however, somewhat better convergence (by an order of ca. ) is apparent. This is a well-known phenomenon also experienced in numerical experiments for the DG method on regular meshes. Moreover, the first subcell refinement step has the tendency to show a deteriorated rate of convergence – presumably due to an increasing constant when switching between the two branches in the proof of Theorem 4.8 (discriminating between locally refined and not locally refined elements).

In Fig. 5.1, one can see the respective convergence plots for different local refinement strategies. Note that the solutions for and are very similar, their error plots in Fig. 5.1 lie on top of each other. The error plots for the local refinements with (dashed lines) and (solid lines) are shown in Fig. 5.1 (left). We observe that the slopes of the error plots match (or exceed by ca. ) the convergence rates in Theorem 4.8. In Fig. 5.1 (right), the convergence for fixed and successively refined is shown. In line with Theorem 4.8, we observe order of convergence one for numerical methods with and order of convergence two for .

### 5.2. Solid body rotation

As the next benchmark problem, we use solid body rotation test proposed by LeVeque [LeVeque1996]. It consists of a slotted cylinder, a sharp cone, and a smooth hump (see Fig. 5.2 (right)) that are placed in a square domain and transported by a time-independent velocity field

 a(t,x1,x2)=(0.5−x2x1−0.5)

in a counterclockwise rotation about . Using and , we choose the following initial data

 u0(x)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩1if(x1−0.5)2+(x2−0.75)2≤r∧(x1≤0.475∨x1≥0.525∨x2≥0.85)1−G(x,(0.50.25))if(x1−0.5)2+(x2−0.25)2≤r14(1+cos(πG(x,(0.250.5))))if(x1−0.25)2+(x2−0.5)2≤r0otherwise

At the inlet , we prescribe the Dirichlet boundary condition . The right-hand side of the advection equation is given by . In order to obtain a discrete initial condition preserving the bounds of the analytical solution (), we define using the -projection into the space of piecewise constant functions instead of our special EG projection operator .