A stabilized cut discontinuous Galerkin framework: II. Hyperbolic problems

07/15/2018 ∙ by Ceren Gürkan, et al. ∙ Umeå universitet 0

We present the second part of a stabilized cut discontinuous Galerkin (cutDG) framework initiated in [1] and develop a novel cutDG method for scalar hyperbolic problems. The domain of interest is embedded into a structured, unfitted background mesh in R^d where the domain boundary can cut through the mesh in an arbitrary fashion. To cope with robustness problems caused by small cut elements, we employ and extend ghost penalty techniques from recently developed, continuous cut finite element methods instead of the cell merging technique commonly used in unfitted discontinuous Galerkin methods, thus allowing for a minimal extension of existing fitted discontinuous Galerkin software to handle unfitted geometries. Identifying a few abstract assumptions on the ghost penalty, we derive geometrically robust a priori error and condition number estimates for the scalar hyperbolic problem which hold irrespective of the particular cut configuration. Possible realizations of suitable ghost penalties are discussed. The theoretical results are corroborated by a number of computational studies for various approximation orders and for two and three-dimensional test problems.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 18

page 19

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

1.1 Background and earlier work

The numerical solution of partial differential equations (PDEs) based on the finite element method requires the generation of high quality meshes to ensure both a proper geometric resolution of the domain features and good approximation properties of the numerical scheme. But even with today’s advanced meshing techniques and continuously growing computer power, the generation of high quality meshes in realistic application problems can be a major bottleneck in the overall simulation pipeline. For instance, most traditional mesh-based numerical solution methods for geological flow and transport models require a series of highly complex preprocessing steps to translate geological image data into conforming domain discretizations which respect intricate geometric structures such as faults and large-scale networks of fractures 

[2]. Similar non-trivial preprocessing steps are necessary when mesh-based domain descriptions are extracted from biomedical image data [3; 4]. The problem of efficient and fast mesh generation becomes even more pronounced when the geometry of the model domain changes drastically in the course of the simulation. For instance, in fluid structure interaction problems [5; 6; 7], multiphase flows [8; 9; 10; 11] or shape optimization problems [12; 13; 14], large or even topological changes in the model geometry may cause even state of the art mesh moving technologies to fail, leaving a costly remeshing as the only resort.

To ease the burden of mesh generation, novel so-called unfitted finite element methods have gained much attention in recent years, see [15] for a recent overview. In unfitted finite element methods, the mesh is exonerated from the task to represent the domain geometry accurately and only used to define proper approximation spaces. The geometry of the model domain is described independently by the means of a separate geometry model, e.g., a CAD or level-set based description, or even another, independently generated tessellation of the geometry boundary. Combined with suitable techniques to properly impose boundary or interface conditions, for instance by means of Lagrange multipliers or Nitsche-type methods [16; 17; 18; 19], complex geometries can be simply embedded into an easy-to-generate background. As the embedded geometry can cut arbitrarily through the background mesh, a main challenge is to devise unfitted finite element methods that are geometrical robustness in the sense that standard a priori error and conditioning number estimates also hold in the unfitted case with constants independent of the particular cut configuration.

One possible solution to achieve geometrical robustness is provided by the approach taken in recent development of the cut finite element method (CutFEM). A distinctive feature of the cut finite element method is that it provides a general, theoretically founded stabilization framework which, roughly speaking, transfers stability and approximation properties from a finite element scheme posed on a standard mesh to its cut finite element counterpart. As a result, a wide range of problem classes has been treated including, e.g., elliptic interface problems [20; 21; 22], Stokes and Navier-Stokes type problems [23; 24; 25; 26; 27; 28; 29; 30; 31], two-phase and fluid-structure interaction problems [32; 9; 33; 34]. As a natural application area, unfitted finite element methods have also been proposed for problems in fractured porous media [35; 36; 37; 38]. Finally, we also mention the finite cell method as another important instance of an unfitted finite element framework with applications to flow and mechanics problems [39; 40; 41], see [42] for a review and references therein.

In addition to the aforementioned unfitted continuous finite element methods, unfitted discontinuous Galerkin (DG) methods have successfully been devised to treat boundary and interface problems on complex and evolving domains [43; 44; 45], including flow problems with moving boundaries and interfaces [46; 47; 48; 49; 50; 51]. In contrast to stabilized continuous cut finite element methods, in unfitted discontinuous Galerkin methods, troublesome small cut elements can be merged with neighbor elements with a large intersection support by simply extending the local shape functions from the large element to the small cut element. As inter-element continuity is enforced only weakly, no additional measures need to be taken to force the modified basis functions to be globally continuous. Consequently, cell merging in unfitted discontinuous Galerkin methods provides an alternative stabilization mechanism to ensure that the discrete systems are well-posed and well-conditioned. For a very recent extension of the cell merging approach to continuous finite elements, we refer to [52; 53]. Thanks to their favorable conservation and stability properties, unfitted discontinuous Galerkin methods remain an attractive alternative to continuous CutFEMs, but some drawbacks are the almost complete absence of numerical analysis except for [54; 55], the implementational labor to track changing matrix sparsity patterns when merging cut elements, and the lack of natural discretization approaches for PDEs defined on surfaces.

So far, most of the theoretically analyzed unfitted finite element schemes have been proposed for elliptic type of problems. Much less attention has been paid to the theoretical development of unfitted finite element methods for advection-dominant or hyberbolic problems, except for [27; 28] considering a continuous interior penalty (CIP) based CutFEM for the Oseen problem, and [56; 57; 58] developing continuous unfitted finite element methods for advection–diffusion equations on surfaces based on the CIP and streamline upwind Petrov Galerkin (SUPG) approach. On the other hand, for the DG community, hyperbolic problems have been of interest from the very beginning, as the first DG method was developed by Reed and Hill [59] to model neutron transport problem. The very first analysis of these methods was presented in [60] and then improved in [61], while Brezzi et al. [62] reformulated and generalized the upwind flux strategy in DG methods for hyperbolic problems by introducing a tunable penalty parameter. The tempting conservation and stability properties, the high locality and particularly the naturally inherited upwind flux term in the bilinear form makes DG methods popular to handle specifically advection dominated problems [63; 64; 65] as well as elliptic ones [66; 67]. As the literature on DG methods for both elliptic and hyberbolic problems is vast, we do not attempt to give a detailed overview but rather refer to the excellent monographs [68; 69], which also includes comprehensive bibliographies.

1.2 Novel contributions and outline of this paper

In this work we continue the development of a novel stabilized cut discontinuous Galerkin (cutDG) framework initiated in [1]. Departing from the elliptic boundary and interface problems considered in [1], we here propose and analyze a cutDG method for scalar hyperbolic equations. In the upcoming work [70], we will combine the presented framework with extension of [71; 72] to introduce cutDGM for mixed-dimensional, coupled problems. The idea of extending stabilization techniques from the continuous CutFEM approach [73] to discontinuous Galerkin based discretizations allows for a minimally invasive extension of existing fitted discontinuous Galerkin software to handle unfitted geometries. Only additional quadrature routines need to be added to handle the numerical integration on cut geometries, and while not being a completely trivial implementation task, we refer to the numerous quadrature algorithms capable of higher order geometry approximation  [74; 45; 75; 76; 77] which have been proposed in recent years. With a suitable choice of the ghost penalty, the sparsity pattern associated with the final system matrix does not require any manipulation and is identical to its fitted DG counterpart.

We start by briefly recalling the advection-reaction model problem and its corresponding weak formulation in Section 2, followed by the presentation of the stabilized cut discontinuous Galerkin method in Section 3 which includes an additional abstract ghost penalty. A main feature of the presented numerical analysis is that we identify a number of abstract assumptions on the ghost penalty to prove geometrically robust optimal a priori error and condition number estimates which hold independent of the particular cut configuration. To prepare the a priori error analysis, Section 4 collects a number of useful inequalities and explains the construction of an unfitted but stable projection operator. In Section 5, we derive an inf-sup condition in a ghost penalty enhanced streamline-diffusion type norm. A key observation is that the classical argument in [62] based on boundedness on orthogonal subscales does not work as the local orthogonality properties of the unfitted projection is perturbed due to the mismatch of the physical domain and the active mesh where the discontinuous ansatz functions are defined. In Section 6, the results from the previous sections are combined to establish an optimal a priori error estimate for the proposed stabilized cutDG method with constants independent of the particular cut configuration. We continue our abstract analysis in Section 7 and prove that condition number scales like , again with geometrically robust constants. Our theoretical investigation concludes with discussing a number of ghost penalty realizations in Section 8. Finally, we corroborate our theoretical analysis with numerical examples in Section 9 where we study both the convergence properties and geometrical robustness of the proposed cutDG method.

1.3 Basic notation

Throughout this work, , denotes an open and bounded domain with piecewise smooth boundary , while denotes its topological boundary. For and , , let be the standard Sobolev spaces consisting of those -valued functions defined on which possess -integrable weak derivatives up to order . Their associated norms are denoted by . As usual, we write and and for the associated inner product and norm. If unmistakable, we occasionally write and for the inner products and norms associated with , with being a measurable subset of . Any norm used in this work which involves a collection of geometric entities should be understood as broken norm defined by whenever is well-defined, with a similar convention for scalar products . Any set operations involving are also understood as element-wise operations, e.g., and allowing for compact short-hand notation such as and . Finally, throughout this work, we use the notation for for some generic constant (even for ) which varies with the context but is always independent of the mesh size and the position of relative to the background , but may be dependent on the dimensions , the polynomial degree of the finite element functions, the shape regularity of the mesh and the curvature of .

2 Model problem

For a given vector field

and a scalar function we consider the advection-reaction problem of the form

(2.1a)
(2.1b)

where describes the boundary value on the inflow boundary of defined by

(2.2)
Here, denotes the outer normal associated with . Correspondingly, the outflow and characteristic boundary are defined by, respectively,
(2.3)
(2.4)

As usual in this setting, we assume that

(2.5)

The corresponding weak formulation of scalar hyperbolic problem (2.1) is to seek such that

(2.6)

with the bilinear form and linear form given by

(2.7)

3 A cut discontinuous Galerkin method for the advection-reaction problem

Let be a background mesh consisting of -dimensional, shape-regular (closed) simplices covering . As usual, we introduce the local mesh size and the global mesh size . For we define the so-called active mesh

(3.1)

and its submesh consisting of all elements cut by the boundary,

(3.2)

Note that since the elements are closed by definition, the active mesh still covers . The set of interior faces in the active background mesh is given by

(3.3)
Figure 3.1: Computational domains for the boundary value problem (2.1). (Left) Physical domain with boundary and outer normal . (Right) Background mesh and active mesh used to define the approximation space. Dashed faces are used to define face based ghost penalties discussed in Section 8.

To keep the technical details the forthcoming numerical analysis at a moderate level, we make two reasonable geometric assumptions on and .

Assumption G1

The mesh is quasi-uniform.

Assumption G2

For there is an element in with a “fat” intersection111The constant is typically user-defined. such that

(3.4)

for some mesh independent . Here, denotes the set of elements sharing at least one node with .

Quasi-uniformity of is assumed mostly for notational convenience. Except for the condition number estimates, all derived estimates can be easily localized to element or patch-wise estimates. The “fat” intersection property guarantees that is reasonably resolved by the active mesh . It is automatically satisfied if, e.g., and is small enough, that is; if where are the principal curvatures of , see [22]. On the active mesh , we define the discrete function space as the broken polynomial space of order ,

(3.5)

To formulate our cut discontinuous Galerkin discretization of the weak problem (2.6), we recall the usual definitions of the averages

(3.6)
(3.7)

and the jump across an interior facet ,

(3.8)

Here, for some chosen unit facet normal on . To keep the notation at a moderate level, we will from hereon drop any subscripts indicating whether a normal belongs to a face or to the boundary as it will be clear from context. With these definitions in place, our cut discontinuous Galerkin method for advection-reaction problem (2.1) can be formulated as follows: find such that

(3.9)

The discrete bilinear forms and represent the classical discrete, discontinuous Galerkin counterparts of (2.6) defined on mesh elements and faces restricted to the physical domain . More precisely, for , we let

(3.10)
(3.11)

The main difference to the classical discontinuous Galerkin method for Problem (2.1), cf. [62; 68], is the appearance of an additional stabilization form , also known as ghost penalty. The role of the ghost penalty is to ensure that the favorable stability and approximation properties of the classical upwind stabilized discontinuous Galerkin method carry over to the unfitted mesh scenario. The precise requirements on will be formulated as a result of the theoretical analysis performed in the remaining sections, but we already point out that since Problem (2.1) involves an advection and a reaction term, it will be natural to assume that can be decomposed in a reaction related ghost penalty and an advection related ghost penalty , both of which we assume to be symmetric,

(3.12)

4 Norms, useful inequalities and approximation properties

Following the presentation in [68], we introduce a reference velocity and a reference time ,

(4.1)

and define the central flux, upwind and the scaled streamline diffusion norm by setting

(4.2)
(4.3)
(4.4)

respectively. Thanks to the assumed quasi-uniformity, the global scaling factor is given by

(4.5)

For each norm with , its ghost penalty enhanced version is given by

(4.6)

To simplify a number of estimates, we will also use a slightly stronger norm than , namely

(4.7)

Finally, we assume as in [68] that the mesh  is sufficiently resolved such that the inequality

(4.8)

is satisfied, which means that Problem (2.1) can be considered as advection dominant and that the advective velocity field is sufficiently resolved in the sense that the following estimates hold,

(4.9)

Before we turn to the stability and a prior error analysis in Section 5, we briefly review some useful inequalities needed later and explain how a suitable approximation operator can be constructed on the active mesh . Recall that for , the local trace inequalities of the form

(4.10)
(4.11)

hold, see [78] for a proof of the second one. For , let be the -th total derivative, then the following inverse inequalities hold,

(4.12)
(4.13)
(4.14)

see [78] for a proof of the last one. Next, to define a suitable approximation operator, we depart from the -orthogonal projection which for and satisfies the error estimates

(4.15)
(4.16)

whenever . Now to lift a function to

, where we for the moment use the notation

, we recall that for Sobolev spaces , , there exists a bounded extension operator satisfying

(4.17)

for , see [79] for a proof. We can now define an unfitted projection variant by setting

(4.18)

Note that this -projection is slightly perturbed in the sense that it is not orthogonal on but rather on . Combining the local approximation properties of with the stability of the extension operator , we see immediately that satisfies the global error estimates

(4.19)
(4.20)
(4.21)

To assure that the consistency error caused by does not affect the convergence order, we require that the ghost penalties and are weakly consistent in the following sense.

Assumption HP1 (Weak consistency estimate)

For and , the semi-norms and satisfy the estimates

(4.22)
(4.23)

With these assumptions, the approximation error of the unfitted projection can be quantified with respect to the norm.

Proposition 4.1

Let and assume that . Then for , the approximation error of satisfies

(4.24)

Set , then by definition

(4.25)
(4.26)

Thanks to Assumption HP1, we only need to estimate . Using (4.19), is bounded by

(4.27)

while the last three terms can be bounded by combining (4.19)–(4.21) leading to

(4.28)

5 Stability analysis

The goal of this section is show that the stabilized cutDG method (3.9) satisfies a discrete inf-sup condition with respect to the -norm, similar to the results presented in [80; 68] for the corresponding classical fitted discontinuous Galerkin method. In the unfitted case, one main challenge is that the proof of the inf-sup condition involves various inverse estimates which in turn forces us to gain control over the considered norm on the entire active mesh.

We start our theoretical analysis by recalling that a partial integration of the element contributions in bilinear form (3.10) leads to a mathematically equivalent expression for , namely,

(5.1)

By taking the average of (3.10) and (5.1), we obtain the well-known decomposition of

into a symmetric and skew-symmetric part,

(5.2)
where
(5.3)
(5.4)
(5.5)

Consequently, the total bilinear form is discretely coercive with respect to :

Lemma 5.1

For it holds

(5.6)

As in the classical fitted mesh case, the proof follows simply from observing that

(5.7)
(5.8)

noting that thanks to assumption (2.5) and the definition of , cf. (4.1). ∎ We point out that for the advection-reaction problem, the ghost penalty  is not needed to establish discrete coercivity in the  norm, in contrast to the cutDGM for the Poisson problem presented in [1]. Nevertheless, we will see that is required to prove an inf-sup condition for the discrete problem (3.9) using the stronger  norm.

Following the derivation of the inf-sup condition in the fitted mesh case, see for instance [81], a proper test function needs to be constructed which establishes control over the streamline derivative. To construct such a suitable test function, we assume the existence of discrete vector field satisfying

(5.9)

Note that for such a discrete vector field can always be constructed. To establish the inf-sup result, the ghost penalties and for the discretized advection-reaction problem are supposed to satisfy the following assumptions.

Assumption HP2

The ghost penalty extends the norm from the physical domain to the active mesh in the sense that for , it holds that

(5.10)
Assumption HP3

The ghost penalty extends the streamline diffusion semi-norm from the physical domain to the active mesh in the sense that for , it holds that

(5.11)

At first it might be more natural to assume that

(5.12)

is valid for instead of the more convoluted estimate (5.11), but the forthcoming numerical analysis as well as the actual design of will require to frequently pass between certain locally constructed discrete vector fields and the original vector field . This is also the content of the next lemma.

Lemma 5.2

Let and let be an element patch with . Let be a patch-wide defined velocity field satisfying

(5.13)

Then

(5.14)

Before we turn to proof we note an immediate corollary of Lemma 5.2 and Assumption HP3.

Corollary 5.3

Let be a collection of patches and assume that the number of patch overlaps is uniformly bounded. Let and satisfy the assumptions in Lemma 5.2. Then

(5.15)

[Lemma 5.2] Thanks to the assumptions (5.13) and the fact that , we have the following chain of estimates,

(5.16)

Thus an application of the Cauchy-Schwarz inequality and the inequality (4.12) to yields

(5.17)

∎ Note that thanks to assumption and the existence of an extension operator , cf. Section 4, a patch-wise defined velocity field satisfying Lemma 5.2, Eq. (5.13) can be always constructed, e.g., by simply taking the value of at some point in the patch.

Next, we state and prove the main result of this section, showing that the cut discontinuous Galerkin method (3.9) for the advection-reaction problem (2.1) is inf-sup stable with respect to the norm. The main challenge here is to show that this result holds with a stability constant which is independent of the particular cut configuration.

Theorem 5.4

Let be the bilinear form defined by (3.9). Then for it holds that

(5.18)

with the hidden constant independent of the particular cut configuration.

We follow the proof for the fitted DG method, see for instance [68], to construct a suitable test function for given such that (5.18) holds.
Step 1. Setting , we gain control over the , thanks to the coercivity result (5.6),

(5.19)

Step 2. Next, set and note that since . Then

(5.20)
(5.21)

Term can be bounded by successively applying a Cauchy-Schwarz inequality and (5.14),

(5.22)

To estimate the remaining terms , let us for the moment assume that the stability estimate

(5.23)

holds. Then one can easily see that

(5.24)

Combing (5.22) and (5.24), we can estimate the right-hand side of (5.21) further by employing a Young inequality of the form yielding

(5.25)
(5.26)
(5.27)
(5.28)

for some constant and . This gives us the desired control over the streamline derivative.
Step 3. To construct a suitable test function for given , we set now . Thanks to stability estimate (5.23) we have and thus combining (5.19) and (5.28) leads us to

(5.29)

for some constant and small enough. Dividing by and taking the supremum over proves (5.18). To complete the proof, it remains to establish the stability bound (5.23).
Estimate (5.23). We start by unwinding the definition of ,

(5.30)
(5.31)

The main tool to estimate is inequality (5.15) with , and . Inserting the definition of and the fact that thanks to assumption (4.8) yields

(5.32)

The second term can be dealt with by recalling the definition of , applying the inverse estimate (4.12) and subsequently moving to the streamline diffusion norm via (5.15),

(5.33)
(5.34)

Next, invoking the inverse trace estimates (4.14) followed by an application of (5.15), we see that

(5.35)
and similarly for IV,
(5.36)

After using (5.11), the remaining last term can be estimated by proceeding as for and :

(5.37)

∎ We conclude this section with a couple of remarks, elucidating the role of the ghost penalties and