# Symmetric Interior Penalty Discontinuous Galerkin Discretisations
and Block Preconditioning for Heterogeneous Stokes Flow^{†}^{†}thanks: Author
S. M. Schnepp acknowledges financial support from the Swiss University Conference and the Swiss Council of Federal Institutes of Technology through the
Platform for Advanced Scientific Computing (PASC) program.
Author D. E. Charrier acknowledges financial support from the
European Unionâs Horizon 2020 research and innovation programme under grant
agreement No 671698.

###### Abstract

Provable stable arbitrary order symmetric interior penalty discontinuous Galerkin (SIP) discretisations of heterogeneous, incompressible Stokes flow utilising – elements and hierarchical Legendre basis polynomials are developed and investigated. For solving the resulting linear system, a block preconditioned iterative method is proposed. The nested viscous problem is solved by a -multilevel preconditioned Krylov subspace method. For the -coarsening, a twolevel method utilising element-block Jacobi preconditioned iterations as a smoother is employed. Piecewise bilinear () and piecewise constant () -coarse spaces are considered. Finally, Galerkin -coarsening is proposed and investigated for the two -coarse spaces considered. Through a number of numerical experiments, we demonstrate that utilising the coarse space results in the most robust -multigrid method for heterogeneous Stokes flow. Using this coarse space we observe that the convergence of the overall Stokes solver appears to be robust with respect to the jump in the viscosity and only mildly depending on the polynomial order . It is demonstrated and supported by theoretical results that the convergence of the SIP discretisations and the iterative methods rely on a sharp choice of the penalty parameter based on local values of the viscosity.

Key words. heterogeneous Stokes flow, variable viscosity, incompressible flow, block preconditioners, DG, SIP, Galerkin multigrid, geodynamics,

AMS subject classifications. 76D07,65M55,65N30,65N12

## 1 Introduction

### 1.1 Background and Motivations

Earth exhibits a diverse range of unique geological processes: mountain building, subduction and continental rifting, earthquakes and volcanism. These phenomena are the result of multi-phase, history-dependent, large-deformation processes spanning million year time scales.

Computational models provide a viable technique to study the evolution in both space and time of geological processes. A prototypical continuum description of the behaviour of rocks is stationary, incompressible Stokes flow with Boussinesq approximation [44, 40]:

\hb@xt@.01(1.1) |

where , , is the velocity, pressure and strain rate, respectively, is the temperature, is the material composition, the effective viscosity, is the reference density at reference temperature , and is the gravity vector. The conservation of momentum and mass for the creeping fluid is coupled with the conservation of energy equation:

\hb@xt@.01(1.2) |

where is the specific heat at constant pressure, the conductivity, and the external heat source; and the evolution of the composition:

\hb@xt@.01(1.3) |

From high-pressure and temperature laboratory experiments of minerals, it is known that rocks exhibit thermally activated creep and follow an Arrhenius type law [37, 26]:

\hb@xt@.01(1.4) |

where is a compositional dependent experimentally determined constant, is the second invariant of the strain rate tensor, are the activation energy and activation volume, is the universal gas constant and is the power-law exponent. To facilitate brittle behaviour at low temperature, the ductile creep laws are augmented with a plasticity model (e.g. Drucker-Prager [37]).

When such a composite flow law is applied to geodynamics scenarios, the effective viscosity () is highly heterogeneous. At depths km, ductile behaviour dominates and the viscosity profile can to first order be characterised by a smooth, exponential function. Above, due to material failure or compositional variations associated with crustal layers, the viscosity profile will be discontinuous and can possess jumps on the order of – Pa s. Realistic forward models of both mantle- and crust-scale simulations are adversely affected by the degree of heterogeneity within the viscosity due to both accuracy concerns associated with the particular spatial discretisation used, and a lack of solver (linear and nonlinear) robustness.

### 1.2 Related work

Geodynamics forward models of incompressible Stokes which permit highly heterogeneous viscosity structures have traditionally utilised finite difference (FD) (e.g. [48, 49]), finite volume (e.g. [23, 41, 42]), or finite element (FE) (e.g. [35, 21, 33, 36, 9, 28, 31]) spatial discretisations. The relative merits of FD and FE methods for geodynamics applications can be broadly summarised as follows:

Staggered grid FD methods are “cheap” (few non-zeros in the stencil), the general implementation is rather straight-forward. However, geometric flexibility is limited, and boundary condition imposition is non-trivial. Introducing new physics may further require the development of a modified stencil; see e.g. [22, 18]. In the context of nonlinear problems, Newton linearisation causes stencil growth and thus increases the overall cost of the discretisation.

Whilst being more expensive than FD methods (on the same grid), inf-sup stable FE methods permit geometric versatility and natural boundary conditions are trivial to impose. Spatial adaptivity () can be readily introduced without requiring redeveloping the underlying numerical method (e.g. [30, 13, 27]). Newton linearisation does not cause the equivalent of stencil growth. For both FD and FE discretisations, robust multi-level preconditioners suitable for highly heterogeneous viscosity structures exist [42, 13, 27, 31].

Inf-sup stable discontinuous Galerkin (DG) methods of the interior penalty type for the Stokes equations can be constructed by using the tensor product element pairs –, and – [43, 38], as well as the -conforming Raviart-Thomas, Brezzi-Douglas-Marini, and Brezzi-Douglas-Fortin-Marini kind element pairs; see [47, 14, 3, 25].

The fact that adaptivity in space and approximation order ( or “” in the preconditioner context) is realised with comparably less effort than for other discretisations makes DG very appealing for geodynamics applications given the very nature of mantle-lithosphere-crust systems. The inf-sup constants of stable DG discretisations are further not sensitive to the element aspect ratio which is highly desirable, for example, within a crustal scale model with a domain spanning km where the domain itself possesses a high aspect ratio, or if anisotropic refinement was employed [39, Theorem 9].

Regarding the solution of the equation system arising from symmetric interior penalty DG (SIP) [1, 2] based discretisations of incompressible Stokes flow, we note that for -conforming discretisations, efficient preconditioners have been introduced very recently [3] [25]. Recent advances in developing efficient and robust solvers for interior penalty DG discretisations of second order elliptic problems with heterogeneous coefficients involve the algebraic multigrid preconditioner proposed in [8, 7], as well as the twolevel methods proposed in [15] and [45, 46].

To the best of our knowledge, employing DG methods for heterogeneous Stokes flow problems in geodynamics was so far only considered in [29]. Preconditioning was not discussed there.

### 1.3 Contributions

We examine the applicability of using mixed SIP based Stokes discretisations for studying heterogeneous, incompressible Stokes flow problems associated with prototype problems arising in geodynamics. Through comparison with an analytic solution employing a discontinuous viscosity structure, we numerically demonstrate that the discretisation yields optimal order of accuracy for – elements for .

Our main contribution is the development of a preconditioned iterative method for the discrete saddle point system resulting from the Stokes discretisation. To this end, we follow a block preconditioning approach; e.g. [20]. The nested viscous problem is solved by a -multilevel preconditioned Krylov subspace method. For the -coarsening, a twolevel method utilising element-block Jacobi preconditioned iterations as a smoother is employed. Two -coarse spaces are considered: the space of element-wise constants and the space of continuous, element-wise bilinear functions Through numerical experiments with heterogeneous viscosity, we demonstrate that the variant utilising the element-wise bilinear coarse space has a convergence rate which is independent of the number of elements, largely insensitive to the jump in viscosity, and only weakly dependent on the polynomial order.

The heterogeneous nature of the viscosity in geodynamics applications requires a careful choice of the SIP penalty parameters. We provide a brief analysis of the influence of the penalty parameters on the discretisation error as well as on the quality of the element-block Jacobi smoothers.

### 1.4 Limitations

We restrict ourselves to linear problems with element-wise constant viscosity distributions in this study.

## 2 Governing equations

Neglecting nonlinearities and the effect of temperature on the viscosity, we restrict ourselves to a model that is solely depending on the material composition of the rocks in the Earth’s mantle:

\hb@xt@.01(2.1a) | |||||

\hb@xt@.01(2.1b) | |||||

where is a velocity field, is a pressure, denotes the viscosity, denotes a volumetric force, denotes the material composition with , and denotes the (linearised) strain rate tensor with , . Further, , denotes a rectangular domain with boundary consisting of a Neumann part , and a Navier part . We require the velocity and pressure to satisfy homogeneous Neumann boundary conditions, | |||||

\hb@xt@.01(2.1c) | |||||

as well as homogeneous Navier boundary conditions, | |||||

\hb@xt@.01(2.1d) | |||||

where denotes the outward normal to the boundary, and denotes a vector belonging to the tangential space of . In case the Neumann boundary is empty, we ensure uniqueness of the pressure solution by enforcing the following constraint: | |||||

\hb@xt@.01(2.1e) |

We further require that the domain is restrained against rigid motions (A1).

Let us denote by and the usual Sobolev spaces and by and their norms.

Let us assume that the viscosity is bounded according to

\hb@xt@.01(2.2) |

(AS1), and that (AS2), then it is well-known that problem (LABEL:eq:governing_equations:momentum_balance) – (LABEL:eq:governing_equations:navier_dirichlet_bc) s.t. to the other named constraints admits a unique weak solution; see e.g. [11].

## 3 Computational grid and trace operators

Let be a regular Cartesian grid on (Ah1). We refer to the disjoint open sets as elements and denote their diameter by . The number of elements is denoted by . Finally, denotes the outward normal unit vector to the element boundary .

An interior face of is the dimensional intersection , where and are two adjacent elements of . Similarly, a boundary face of is the dimensional intersection which consists of entire faces of . We denote by the union of all interior faces of , by , and the union of all boundary faces belonging to the Neumann part, and the Navier part of the boundary, respectively, and set . Here and in the following, we refer generically to a “face” although we consider only two-dimensional problems in this paper.

Let ; we denote by the space of real-valued functions such that the function and its weak derivatives up to order are measurable and square integrable in . We will denote the norms on all three spaces , , and by the symbol . We further introduce the broken Sobolev space

\hb@xt@.01(3.1) |

Let us denote the norms of , , and by the symbol .

Let , and either belong to , , or . Let be an interior face shared by the elements and . Let and denote the traces of and on from the interior of , respectively. Further, let denote the outward normal unit vector to the boundary . We define the mean value and the jump at by

\hb@xt@.01(3.2) |

Let denote a vector-valued function in , and denote its traces on from the interior of . We define the jumps and at by

\hb@xt@.01(3.3) |

where “” denotes the dyadic product.

## 4 Discretisation of the Stokes problem

Let us introduce , , and define:

\hb@xt@.01(4.1) |

For simplicity, we assume here that the viscosity is element-wise constant (ASh1). Let us additionally define the constants

\hb@xt@.01(4.2) |

where the constants , , stem from the following discrete trace inequality [24]:

###### Lemma 4.1 (Discrete trace inequality)

Let be an affine quadrilateral, and let be an edge belonging to the boundary of . Then, it holds that

\hb@xt@.01(4.3) | ||||

with the trace inequality constant | ||||

\hb@xt@.01(4.4) |

We approximate velocity and pressure in the discontinuous finite element spaces

\hb@xt@.01(4.5) | ||||

\hb@xt@.01(4.6) |

where is the space of polynomials of maximum degree in each variable on the mesh cell .

As approximation to (LABEL:eq:governing_equations:momentum_balance) – (LABEL:eq:governing_equations:pressure), we then consider the problem of finding and such that:

\hb@xt@.01(4.7) | |||||

\hb@xt@.01(4.8) |

where we use a SIP form , and a form similar to the one used in [43]:

\hb@xt@.01(4.9) | ||||

and

\hb@xt@.01(4.10) | ||||

\hb@xt@.01(4.11) |

with , , and the face-wise penalties , . The parameters are the so-called penalty or stability parameters that must be chosen sufficiently large (to be specified below) to guarantee that the bilinear form is coercive on the discrete space .

Consistency of the discrete variational problem can be shown by following the proof of [43, Lemma 7.5.].

### 4.1 Stability of the Stokes discretisation

For the analysis of the Stokes discretisation, it is necessary to introduce the functionals

\hb@xt@.01(4.12) | ||||

\hb@xt@.01(4.13) | ||||

\hb@xt@.01(4.14) |

with and . By definition, is a norm on . Under assumption LABEL:itm:governing_equations:assumption_rigid_motions, we assume that is a norm on for the considered boundary conditions. (ASh2).

Stability of the Stokes discretisation follows from the discrete inf-sup stability of , the discrete coercivity of , and Brezzi’s lemma:

###### Lemma 4.2 (Discrete inf-sup stability)

Let the assumptions LABEL:itm:grid:assumption_grid be satisfied and let , . Then, it holds that

\hb@xt@.01(4.15) |

with independent of and .

We rely on a discrete Korn inequality to show coercivity of the bilinear form .

###### Lemma 4.3 (Discrete Korn inequality,[10])

Let the assumptions LABEL:itm:governing_equations:assumption_rigid_motions and LABEL:itm:grid:assumption_grid be satisfied. Then, it holds for all that

\hb@xt@.01(4.16) |

with a constant independent of .

Proof. This is a slight modification to the result in [10] relying on assumption LABEL:itm:governing_equations:assumption_rigid_motions as well as on the affinity and shape regularity of the grid elements LABEL:itm:grid:assumption_grid.

###### Lemma 4.4 (Discrete coercivity)

Let the assumptions LABEL:itm:governing_equations:assumption_rigid_motions, LABEL:itm:grid:assumption_grid, and LABEL:itm:stokes_discretisation:assumption_eta – LABEL:itm:stokes_discretisation:assumption_norm_h hold. Assume that the penalty parameters are chosen according to

\hb@xt@.01(4.17) |

with from (LABEL:eq:stokes_discretisation:eta_face), and with denoting the number of faces of an element. Then, it holds that

\hb@xt@.01(4.18) |

with a constant independent of and . Here, is from assumption LABEL:itm:stokes:assumption_eta_bounds. The parameter is a small value that is necessary to use the discrete Korn inequality.

Proof. The proof is based on the most part on the stability analysis for the element-wise penalty approach presented in [16, Section 3.3.3.]. Let , ; we apply a Young’s inequality to the second and third term of (LABEL:eq:stokes_discretisation:stiffness_form). Thus, we obtain

\hb@xt@.01(4.19) |

where , . As a next step, we will bound the second and third term of (LABEL:eq:stokes_discretisation:coercivity_proof_1) below. We have for interior faces ,

\hb@xt@.01(4.20) | ||||

In step (I), we have applied the trace inequality (LABEL:eq:stokes_discretisation:discrete_trace) from lemma (LABEL:thm:stokes_discretisation:discrete_trace), and further have introduced constant defined as in (LABEL:eq:stokes_discretisation:const_trace_face), as well as constant defined as in (LABEL:eq:stokes_discretisation:eta_face). Analogously, we obtain for boundary faces ,

\hb@xt@.01(4.21) |

Inserting (LABEL:eq:stokes_discretisation:coercivity_proof_2) – (LABEL:eq:stokes_discretisation:coercivity_proof_3b) in (LABEL:eq:stokes_discretisation:coercivity_proof_1), leads to

\hb@xt@.01(4.22) | ||||

We see that this expression is positive for any if

\hb@xt@.01(4.23) | ||||

\hb@xt@.01(4.24) |

where denotes the number of faces of a quadrilateral element. In the following, we choose the penalty parameters according to (LABEL:eq:stokes_discretisation:sipg_penalty). We add a small value once/twice to each penalty parameter in order to use the discrete Korn inequality (LABEL:eq:stokes_discretisation:discrete_korn) from lemma LABEL:thm:stokes_discretisation:discrete_korn in step (I) of the next derivations. Thus, if we choose the penalty value according to (LABEL:eq:stokes_discretisation:sipg_penalty), we obtain from (LABEL:eq:stokes_discretisation:coercivity_proof_4) that

\hb@xt@.01(4.25) | ||||

where is from assumption LABEL:itm:stokes:assumption_eta_bounds, and is from the discrete Korn inequality (LABEL:eq:stokes_discretisation:discrete_korn) from lemma LABEL:thm:stokes_discretisation:discrete_korn. In step (II), we have used that due to the affinity and shape regularity of the grid elements LABEL:itm:grid:assumption_grid, it holds that for and , with denoting the shape regularity constant. We have further bounded below by . The constant is assigned a different value in every step and independent of and .

We typically choose the penalty parameters as the lower bound since estimate (LABEL:eq:stokes_discretisation:sipg_penalty) is not totally sharp due to the utilised inequalities. Even smaller values can be chosen in practice [24]. We note that the parameter is set to zero in our computations. We further remark that estimate (LABEL:eq:stokes_discretisation:sipg_penalty) assumes piecewise constant viscosity distributions. For piecewise polynomial viscosity distributions, it is necessary to replace by in the discrete trace inequality constant.

### 4.2 A-priori error estimates for the Stokes discretisation

Let us state the continuity properties of the forms and :

###### Lemma 4.5 (Continuity of and )

Under assumptions LABEL:itm:governing_equations:assumption_rigid_motions – LABEL:itm:stokes:assumption_force, LABEL:itm:grid:assumption_grid, and LABEL:itm:stokes_discretisation:assumption_eta, bilinear forms and are continuous in the sense that

\hb@xt@.01(4.26) | |||||

\hb@xt@.01(4.27) |

with and . Here, is from assumption LABEL:itm:stokes:assumption_eta_bounds, and denotes the penalty value on face . Both constants and are independent of and .

Proof. The proof of both inequalities follows from standard inequalities. See e.g. the proofs of [43, Lemma 7.1] and [43, Lemma 7.2] for more details.

Further note that the norms and are equivalent on the discrete spaces and . This is also the case for the norms and .

Using the consistency of the discretisation, the discrete inf-sup stability of , the discrete coercivity of , the continuity of both bilinear forms together with the discrete equivalency of norms, as well as suitable element-wise -interpolants, we can derive the following a-priori error estimate:

###### Lemma 4.6

Let the assumptions. LABEL:itm:governing_equations:assumption_rigid_motions, LABEL:itm:stokes:assumption_eta_bounds – LABEL:itm:stokes:assumption_force, LABEL:itm:grid:assumption_grid, and LABEL:itm:stokes_discretisation:assumption_eta hold. Assume that the weak solution to (LABEL:eq:governing_equations:momentum_balance) – (LABEL:eq:governing_equations:pressure) belongs to , , with and . Further, let and denote the discrete solution to problem (LABEL:eq:stokes_discretisation:discrete_momentum_balance) – (LABEL:eq:stokes_discretisation:discrete_incompressible_constraint). Then, it holds hat

\hb@xt@.01(4.28) |

with and . The constant is independent of and but depends on the shape regularity of the grid elements. The constant is from lemma LABEL:thm:stokes_discretisation_error:continuity, and is from lemma LABEL:thm:stokes_discretisation:discrete_coercivity. Note that this estimate holds up to the discrete inf-sup constant from lemma LABEL:thm:stokes_discretisation:discrete_inf_sup that depends on .

Proof. The proof requires arguments analogous to the ones used for the proofs of [43, Lemma 8.1] and [43, Lemma 8.2].

Utilising a duality argument, additional requirements on the regularity of the solution to the continuous variational problem (and on the solution to the adjoint problem), as well as -interpolants for the whole domain , we can derive the following error estimate for the discrete velocity solution:

###### Lemma 4.7 ( error)

Let the assumptions LABEL:itm:governing_equations:assumption_rigid_motions – LABEL:itm:stokes:assumption_force, LABEL:itm:grid:assumption_grid, and LABEL:itm:stokes_discretisation:assumption_eta hold. Assume further that the weak solution to (LABEL:eq:governing_equations:momentum_balance) – (LABEL:eq:governing_equations:pressure) belongs to , with . Choose the penalty values according to (LABEL:eq:stokes_discretisation:sipg_penalty) in lemma LABEL:thm:stokes_discretisation:discrete_coercivity. Let and denote the discrete solution to problem (LABEL:eq:stokes_discretisation:discrete_momentum_balance) – (LABEL:eq:stokes_discretisation:discrete_incompressible_constraint). Then, it holds that

\hb@xt@.01(4.29) |

where the constant is independent of and but depends on the shape regularity of the grid. See Lemma LABEL:thm:stokes_discretisation_error:error_norm_h for a definition of the remaining constants.

###### Remark 4.8

We emphasise that the discretisation errors depend on the size of the largest penalty parameter as well as on the viscosity contrast by means of the constants and .

## 5 Preconditioning the Stokes system

The discrete variational problem (LABEL:eq:stokes_discretisation:discrete_momentum_balance) – (LABEL:eq:stokes_discretisation:discrete_incompressible_constraint) is equivalent to the saddle point problem

\hb@xt@.01(5.1) |

We solve (LABEL:eq:discretestokes) using a right preconditioned Krylov method, with an upper block triangular preconditioner of the form:

\hb@xt@.01(5.2) |

where is the pressure Schur complement.

Our implementation employs tensor products of pairwise orthogonal Legendre polynomials as basis functions for the pressure and the velocity components. The zero pressure average is not build into the pressure basis functions. In case the Neumann boundary is empty, we thus solve a singular system. As we typically use Krylov methods which are mathematically equivalent to GMRES, we then require that the right-hand side is consistent (e.g., we remove the constant pressure null space) [12, Theorem 2.4.].

Noting that both is symmetric positive definite (stemming from the SIP formulation) and is symmetric positive definite (stemming from the inf-sup stability), this choice for will result in convergence in at most two iterations in exact arithmetic [19]. Whilst optimal (in the sense of iteration counts), the definition of