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 elementlocal 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 recast 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 EGtype 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 highorder Bernstein finite elements [HajdukEtAl2020a, HajdukEtAl2020b, KuzminMQL2020, Lohmann2019]. An adaptive approach to subcell limiting was introduced in [KuzminEtAl2019]. Using continuous blending functions, a highorder finite element approximation on a large macrocell was combined with a boundpreserving piecewise (multi)linear subcell approximation.
Another wellknown 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 ADERDG schemes proposed in [Dumbser2008] and possess a very attractive capability to detect high and lowregularity solution behavior. The underlying a posteriori limiting strategy was inspired by the Multidimensional Optimal Order Detection (MOOD) approach originally developed for finite volumes. In the context of the ADERDG methods, physical and numerical admissibility conditions are enforced by, first, advancing the solution in time using a highorder DG method on the coarse mesh, and, for troubled cells, repeating the last time step locally via a loworder 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 twomesh 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 subcellenriched 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^{1}^{1}1https://github.com/FESTUNG framework [FrankRAK2015, ReuterAWFK2016, JaustRASK2018, ReuterHRFAK2019, ReuterRAFK2020] based on our EG scheme for the shallowwater equations [HauckAFHR2020].
1.1. Model problem
We consider a nonstationary advection equation on a bounded Lipschitz domain (with ). The precise formulation of the linear hyperbolic problem to be solved is as follows:
(1.1) 
for a given velocity field and a righthand side function . Additionally, initial data is prescribed, and we denote by the outward unit normal to . Furthermore, we assume that the inflow boundary
is independent of time and disjointly subdivided into Dirichlet and flux boundaries (this subdivision is also assumed to be independent of time), i.e.
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 nonoverlapping 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 nonconformal. 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
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 quasiuniform if there exists a constant such that for all , all , , and all , we have
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 elementwise polynomials of degree at most
(simplices) or tensorproduct polynomials of degree at most
in each spatial coordinate (quadrilaterals/hexahedrons) without any continuity constraints. Thus,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
(2.1)
2.2. Semidiscrete formulation
The semidiscrete 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
for a scalar that is elementwise smooth enough to have traces. On , this definition is modified as follows:
Hence, the jump turns a scalar into a vector. Also note the following property of jumps used in our analysis
Given a velocity field , we define the upwind value of as
where is the standard signum function.
Using this notation, we can formulate our semidiscrete problem
with trial function and test function from for almost every and , where
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
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 righthand 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.,
(4.1) 
Proof.
Result follows directly from the fact that and the orthogonality of . ∎
Lemma 4.2 (Best approximation property of ).
For all , , and all ,
(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
(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
(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
(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
(4.6)  
(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
(4.8)  
(4.9) 
Proof.
The first inequality is the observation that is the elementwise 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 quasiuniform mesh (double) sequence, and let , . Then, the EG approximation converges in to the analytical solution , i.e., there exists independent of and such that
with for simplicial meshes and or general meshes and (i.e. in the case of DG). Otherwise, .
Proof.
Defining
we have due to the consistency and since that
This can be rewritten as
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 elementwise 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 :
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
In the second case, we obtain for using the same arguments
The estimate for is conducted analogously. Here, the projection is used to obtain
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
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 righthand 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 .
space  , ,  , ,  

1  2  3  4  5  1  2  3  4  5  
1  5.85E01  —  —  —  —  4.48E01  —  —  —  — 
2  5.42E01  2.85E01  —  —  —  2.95E01  7.50E02  —  —  — 
3  2.95E01  1.95E01  7.24E02  —  —  2.21E01  7.39E02  1.06E02  —  — 
4  1.60E01  1.10E01  7.37E02  1.88E02  —  1.23E01  6.79E02  9.67E03  1.46E03  — 
5  8.40E02  5.78E02  4.18E02  1.99E02  4.80E03  6.46E02  3.93E02  8.73E03  1.27E03  2.07E04 
6  4.33E02  2.98E02  2.21E02  1.12E02  5.16E03  3.33E02  2.10E02  5.03E03  1.17E03  1.59E04 
7  2.21E02  1.53E02  1.15E02  5.90E03  2.89E03  1.70E02  1.10E02  2.68E03  6.72E04  1.49E04 
space  , ,  , ,  
1  2  3  4  5  1  2  3  4  5  
1  4.47E01  —  —  —  —  4.47E01  —  —  —  — 
2  2.95E01  7.11E02  —  —  —  1.72E01  7.11E02  —  —  — 
3  2.21E01  7.34E02  9.80E03  —  —  6.62E02  5.96E02  9.80E03  —  — 
4  1.23E01  6.76E02  9.59E03  1.29E03  —  1.76E02  1.63E02  7.57E03  1.29E03  — 
5  6.46E02  3.91E02  8.69E03  1.26E03  1.63E04  4.51E03  4.24E03  2.10E03  1.00E03  1.63E04 
6  3.33E02  2.09E02  5.00E03  1.17E03  1.59E04  1.15E03  1.09E03  5.41E04  2.79E04  1.27E04 
7  1.70E02  1.09E02  2.66E03  6.70E04  1.49E04  3.29E04  2.77E04  1.37E04  7.18E05  3.56E05 
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 wellknown 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 timeindependent velocity field
in a counterclockwise rotation about . Using and , we choose the following initial data
At the inlet , we prescribe the Dirichlet boundary condition . The righthand 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 .
a) , err: 1.28E01  b) , err: 1.15E01 
[draft = false, width=0.33trim=100 20 100 50,clip]pictures/solid_body_4_4_tecplot  [draft = false, width=0.33trim=100 20 100 50,clip]pictures/solid_body_4_5_tecplot 
c) , 
Comments
There are no comments yet.