# BoSSS: a package for multigrid extended discontinuous Galerkin methods

The software package BoSSS serves the discretization of (steady-state or time-dependent) partial differential equations with discontinuous coefficients and/or time-dependent domains by means of an eXtended Discontinuous Galerkin (XDG, resp. DG) method, aka. cut-cell DG, aka. unfitted DG. This work consists of two major parts: First, the XDG method is introduced and a formal notation is developed, which captures important numerical details such as cell-agglomeration and a multigrid framework. In the second part, iterative solvers for extended DG systems are presented and their performance is evaluated.

## Authors

• 3 publications
• 2 publications
• 1 publication
• ### A postprocessing technique for a discontinuous Galerkin discretization of time-dependent Maxwell's equations

We present a novel postprocessing technique for a discontinuous Galerkin...
10/03/2020 ∙ by G. Nehmetallah, et al. ∙ 0

• ### A mixed discontinuous Galerkin method without interior penalty for time-dependent fourth order problems

A novel discontinuous Galerkin (DG) method is developed to solve time-de...
09/30/2019 ∙ by Hailiang Liu, et al. ∙ 0

• ### Error analysis of Runge–Kutta discontinuous Galerkin methods for linear time-dependent partial differential equations

In this paper, we present error estimates of fully discrete Runge–Kutta ...
01/03/2020 ∙ by Zheng Sun, et al. ∙ 0

• ### Multiphysics Simulation of Plasmonic Photoconductive Antennas using Discontinuous Galerkin Methods

Plasmonic nanostructures significantly improve the performance of photoc...
12/08/2019 ∙ by Liang Chen, et al. ∙ 0

• ### Positivity-preserving third order DG schemes for Poisson–Nernst–Planck equations

In this paper, we design and analyze third order positivity-preserving d...
01/29/2021 ∙ by Hailiang Liu, et al. ∙ 0

• ### A matrix-free Discontinuous Galerkin method for the time dependent Maxwell equations in unbounded domains

A Discontinuous Galerkin (DG) Finite Element Method (FEM) approach for t...
02/20/2020 ∙ by Bernard Kapidani, et al. ∙ 0

• ### Multiphysics Modeling of Plasmonic Photoconductive Devices using Discontinuous Galerkin Methods

Plasmonic nanostructures can significantly improve the performance of ph...
12/08/2019 ∙ by Liang Chen, 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 and Motivation

BoSSS (Bounded Support Spectral Solver) is a flexible framework for the development, evaluation and application of numerical discretization schemes for partial differential equations (PDEs) based on the eXtended Discontinuous Galerkin (XDG) method. The main focus of the XDG method are problems with dynamic interfaces and/or time-dependent domains.

One important field of application for the XDG method are moving geometries: Almost all numerical methods for PDEs require the labor-intensive and thus expensive (semi-manual) creation of numerical grids or meshes. This is particularly cumbersome if time-dependent domains are of interest, e.g. flows with moving boundaries, as can be found when dealing with a moving geometry (e.g. rotors, pistons but also flexible structures such as heart valves). As a side-product, one could use XDG in order to embed a static geometry, which would be cumbersome to mesh with higher order elements, on a Cartesian background mesh.

A second field of application are multiphase flows, like mixtures of oil and water: due to topology changes, e.g. merging or breakup of droplets, such scenarios would be very difficult to handle with some conformal meshing or mesh-motion techniques like arbitrary Lagrangian-Eulerian (ALE) meshes.

#### Purpose and scope of this work

All definitions and algorithms found in this paper have their direct counterpart in the source code. However, this work is not intended to be a software manual – in such, the link to the corresponding place in the source code would have to be established with a suitable kind of cross-reference. Since the source code is constantly being developed, the validity of such references would be short-lived.

Instead, this work is intended to serve interested readers as an axiomatic description of the extended discontinuous Galerkin framework in BoSSS, in order to assess its capabilities and limitations. We therefore emphasize on issues which are not discussed in publications yet, in specific, details on multigrid solvers and their performance will be presented.

### 1.1 Prototype Problems

#### Poisson with a jump in coefficients

Let () be some polygonal domain, which is partitioned into two disjoint but adjacent sub-domains and . These might be referred to as phases, since they represent e.g. the liquid and gaseous phase in multiphase flow simulations. We assume that the interface is a - dimensional manifold which is at least almost everywhere. Therefore, one can define an oriented normal field on ; w.l.o.g., we assume that “points from to , i.e. is, on equal to the outer of and an inner normal of . Then, for any property which is continuous in , one can define the jump operator

 (1)

for . Then, we can define the piece-wise Poisson problem

 (2)

with a discontinuous diffusion coefficient

 μ(→x)={μAfor →x∈A,μBfor →x∈B. (3)

#### Incompressible multiphase flows

In a transient stetting, the phases are usually considered to be time-dependent too, i.e. one has the decomposition . Using the same notation as introduced above, we introduce the incompressible two-phase Navier-Stokes equation for material interfaces as

 (4)

with piece-wise constant density and viscosity for both phases, i.e.

 ρ(→x)={ρAfor →x∈AρBfor →x∈Bandμ(→x)={μAfor →x∈AμBfor →x∈B. (5)

Furthermore, denotes surface tension and denotes the mean curvature of . In a two-phase setting, we also assume that the interface does not touches the boundary, i.e. . If this is not the case, a three-phase contact line occurs, where two fluids and the solid boundary met at a – dimensional manifold. This usually requires more sophisticated boundary conditions (see below).

##### Interface tracking and Level-Set:

Note that problem (4) is incomplete without a specification of the interface motion. In order to track the individual domains, a sufficiently smooth Level-Set field is introduced; for sake of simplicity, might also be called ‘Level-Set’. Then, time-dependent domains , and the interface can be described as

 A(t):= {→x∈Ω; φ(t,→x)<0}, (6) I(t):= {→x∈Ω; φ(t,→x)=0}, (7) B(t):= {→x∈Ω; φ(t,→x)>0}. (8)

On , must therefore comply with the Level-Set equation

 ∂φt+∇φ⋅→u=0, (9)

which states that the interface speed in normal direction is equal to the flow velocity in normal direction. From the , the normal can be computed as

 →nI=∇φ/|∇φ|2 (10)

and the mean curvature can be computed by Bonnets formula as

 κ=div(∇φ|∇φ|2). (11)

If the Level-Set is prescribed, e.g. in the case of some externally forced motion, one can infer the interface speed in normal direction from Eq. (9) as

 s=−∂tφ|∇φ|. (12)

#### Other applications

The XDG framework can be seen as a multi-purpose technology. One obvious application is the use as an immersed boundary method (IBM), e.g. in order to circumvent the (labour-extensive) meshing of domains with complex geometrical details or to represent domains with moving parts. Boundary conditions are weakly enforced at the interface in the same formulation as it would be for boundary fitted methods. The fluid domain is represented by , whereas just describes a void region. This was demonstrated for e.g. for moving body flows where the motion of solid particles is characterized by the Newton-Euler equations in the work of Krause and Kummer Krause and Kummer (2017).

Immersed boundaries can also be used in the context of compressible flows. For inviscid flows around immersed boundaries, such as cylinders and airfoils, Müller et. al. Müller et al. (2017) demonstrated the expected high order of convergence for such geometries using an embedded boundary described by a static Level-Set.

A further issue in the scope of multiphase flows (4) is the interaction of the interface with the domain boundary, i.e. the case : this is referred to as the three-phase contact line between fluids and and the solid. In order to allow motion of the contact line , the no-slip boundary condition has to be relaxed, yielding the so-called Navier-Slip boundary condition. It is notable that the XDG-implementation does not need further manipulation of the contact line velocity or contact angle . Furthermore, XDG allows the implementation of the generalized Navier-slip boundary condition with

 σ(cosΘstat−cosΘ)→nL=βL(→u⋅→nL)→nL, (13)

localized at the contact line , which is a – dimensional manifold.

However, for sake of simplicity and moreover, compactness of the presentation, in this work most issues of the XDG method will be discussed with respect to the Poisson problem (2).

### 1.2 Historical Development and State of the Art

The origins of discontinuous Galerkin (DG) methods can be tracked back to the work of Reed and Hill Reed and Hill (1973). The name ‘Discontinuous Galerkin’ was mainly established through the works of Cockburn and Shu Cockburn and Shu (1991)

, although similar methods, based on broken polynomial approximation, ideas were already established earlier: The probably most popular discretization for the Laplace operator, the so-called symmetric interior penalty (SIP) method was proposed by Arnold

Arnold (1982) about nine years earlier.

The idea of adapting a finite element method (FEM) to allow jumps in parameters can be traced back to the 1970s, where Poisson problems with a jump in the diffusion parameter along a smooth interface were investigated by Babus̆ka Babus̆ka (1970) as well as by Barrett and Elliott Barrett and Elliott (1987). These works mainly relied on isoparametric elements, fitted to the location of the interface at which the discontinuity occurs. Thus, in a transient setting complex motion of the interface is quite difficult to address.

The continuous extended finite element method (XFEM), was presented by Moës et al. Moës et al. (1999) to simulate cracking phenomena in solid mechanics. Those ideas were extended to time-dependent problems by Chessa and Belytschko Chessa and Belytschko (2004). The first XFEM for two-phase flows was presented by Groß and Reusken Groß and Reusken (2007a, b). With the so-called intrinsic XFEM presented by Fries Fries (2009), it became possible to represent also kinks in the velocity field. This work was later extended by Cheng and Fries Cheng and Fries (2012) and further by Sauerland and Fries Sauerland and Fries (2013). The first extended DG (XDG) method (also called unfitted DG or cut-cell DG) was presented by Bastian and Engwer Bastian and Engwer (2009) in order to model flows in porous media. Later, those approaches were extended to multiphase flows by Heimann et. al. Heimann et al. (2013).

### 1.3 The BoSSS code

The development of BoSSS has been initiated at the Chair of Fluid Dynamics in 2008. Since 2017, it is publicly available under the Apache License. The very first motivation for BoSSS is to have a suitable code base for the development of XDG methods. Very early, it was decided to investigate both, compressible as well as incompressible flows. The highlights of the code package are:

• XDG and support for flows with dynamic interfaces, which is the main topic of this paper.

• Suitability for High Performance Computing (HPC): All production algorithms in BoSSS are implemented MPI parallel. Furthermore, BoSSS provides instrumentation output in order to analyze and optimize parallel scaling as well as node-level performance. A very important feature is dynamic load balancing, which allows re-distribution of the computational mesh when local processor load changes, e.g. due to motion of the fluid interface.

• Rapid prototyping capability: New models resp. equations can be added in a notation that is close to the usual presentation of DG methods in textbooks, with low development effort. Technical issues, e.g. the handling of the numerical mesh, are provided by the software library.

• Sophisticated workflow and data management facilities: In order to organize and analyze e.g. large parameter studies, database-centric workflow tools were developed.

## 2 Discontinuous and Extended Discontinuous Galerkin Methods

For sake of completeness, and introducing the notation, we briefly summarize the DG method; The following definitions are fairly standard and can be found in similar form in many textbooks Di Pietro and Ern (2011); Hesthaven and Warburton (2008).

###### Definition 1 (basic notations).

We define:

• a numerical mesh/grid for is a set . The cells are simply connected, cover the whole domain (), but do not overlap ( for ).

• the skeleton of the background mesh: . Furthermore, internal edges: ;

• a normal field on . On , it represents an outer normal, i.e., on , . By , we denote a normal field that is equal to on and equal to on ;

For sake of completeness, we notate the approximation space of the ‘standard’ DG method:

###### Definition 2 (DG space).

The broken polynomial space of total degree is defined as:

 Pk(Kh):={f∈L2(Ω); ∀ K∈Kh:f|K is polynomialand deg(f|K)≤k}; (14)

### 2.1 Extended Discontinuous Galerkin and Level-Set

In order to yield a well-defined interface , the Level-Set-field must be sufficiently smooth. Precisely, we assume to be almost-everywhere in . For certain steady-state problems with simple interface shapes, one can represent by an explicit formula that can e.g. be inserted into to code directly. In more complicated cases, or for temporally evolving problems, may be itself a continuous broken polynomial on the background mesh, i.e. . In both cases, one can infer an interface normal (Eq. 10) as well as the interface speed (Eq 12) from .

###### Definition 3 (Cut-cell mesh).

Time-dependent cut-cells are given as

 Kj,A(t):=Kj∩A(t),Kj,B(t):=Kj∩B(t). (15)

The set of all cut-cells form the time-dependent cut-cell mesh .

Note that we refer to the ‘original cells’ as background cells, in contrast to the time-dependent cut-cells . By , resp. a notation for an arbitrary phase is introduced, i.e. can be either or .

XDG is essentially a DG method on cut-cells, i.e. one can define the XDG space as follows:

###### Definition 4 (XDG space).

We define:

 PXk(Kh,t):=Pk(KXh(t)) (16)

Most DG methods are written in terms of jump- and average operators, as already defined on , see Eq. (1). This notation is extended onto the skeleton of the mesh :

###### Definition 5 (inner and outer value, jump and average operator).

At the mesh skeleton, the inner- resp. outer-value of a field are defined as:

 uin(→x):=limξ↘0u(→x−ξ→nI,Γ)for →x∈Γ∪Iuout(→x):=limξ↘0u(→x+ξ→nI,Γ)for →x∈Γint∪I

Then, the jump and average value operator are defined as

 := {uin−uout% on Γintuinon ∂Ω, (17) {{u}} := {(uin+uout)/2on Γintuinon ∂Ω. (18)

For the implementation of an XDG-method, accurate numerical integration on the cut-cells

is required. In BoSSS the user can either select the Hierarchical Moment Fitting (HMF) procedure, developed by Müller et. al.

Müller et al. (2013), or, alternatively, a method proposed by Saye Saye (2015). While the HMF supports all types of cells (triangles, quadrilaterals, tetra- and hexahedrons) the method of Saye is generally faster, but restricted to quadrilaterals and hexahedrons.

#### Ensuring continuity of φ, computation of normals and curvature

Some XDG method for the two-phase Navier Stokes problem (4) has to be coupled with a second method to compute the evolution of the interface. Since the respective Level-Set equation (9) is of hyperbolic type, a conventional DG approximation seems a reasonable choice. However, such a method will typically yield a discontinuous field . In such cases, needs to be projected to a continuous space before it can be used to setup the XDG method. Furthermore, the curvature required in problem (4), see also Eq. (11), needs to be handled with care. For a detailed presentation of the filtering procedures used in BoSSS, we refer to the work of Kummer and Warburton Kummer and Warburton (2016). For sake of simplicity, in this work we assume to be sufficiently smooth with respect to space and time.

### 2.2 Variational formulations

Now, problems (2) and (4) can be discretized in the XDG space. In general, we are interested in systems like the Navier-Stokes equation, with

dependent variables which are not necessarily discretized with the same polynomial degree. We therefore define the degree-vector

; and introduce an abbreviation for the function space of test and trial functions,

 VXk––(t):=Dv∏γ=1P% Xkγ(Kh,t). (19)

Then, the discrete version of some linear problem, for a fixed time formally read as: find so that

 a(U,V)=b(V)∀V∈VXk––(t). (20)

#### Discrete variational formulation for Poisson Eq. (2)

The variational formulation of the symmetric interior penalty method, originally proposed by Arnold Arnold (1982), reads as

 (21)

for the left-hand-side of Eq. (20). Here, denotes the broken gradient, where differentiation at the jumps on is excluded. The linear form on the right-hand-side of Eq. (20) is given as

 b(v):=∫Ωfv dV−∮ΓDiriμg% Diri(∇hv⋅→n∂Ω−ηv) dS+∮ΓNeuμgNeuv dS (22)

The SIP factor is known to scale as , where is a local length scale of the agglomerated cut-cell, see section 2.3. For certain specific cell shapes, explicit formulas for can be given, a comprehensive overview is given in the thesis of Hillewaert Hillewaert (2013). In the case of XDG methods whith arbitrary cell shapes, some rules-of-thumb can be used, given that a sufficiently large multiplicative constant is used. Alternatively, a symmetric weighted interior penalty (SWIP) form, as presented in the textbook of Di Pietro & Ern Di Pietro and Ern (2011) might be used.

#### Discrete variational formulation for Navier-Stokes Eq. (4)

Due to the non-linearity in the convectional term, the Navier-Stokes system is usually split up into the linear Stokes part and the nonlinear convection operator . It is known that, in order to obtain an inf-sup stable discretization, the DG polynomial degree of the pressure has to be one lower than for velocity, i.e. velocity is discretized with degree and pressure with degree . ()Note this is proven only for special cases, see the work of Girault et. al. Girault et al. (2005).) In this notational framework, for spatial dimension , we write and . The variational formulation of the Navier-Stokes equation formally reads as: find so that

 ∫Ω∂tU⋅V dV+c(U,U,V)+a(U,V)=b(V)∀V∈VXk––(t). (23)

A complete specification of the involved forms would be too lengthy here; hence, we refer to the works of Heimann et. al. Heimann et al. (2013) and Kummer Kummer (2017). Obviously, (i) this equation still requires a temporal discretization and (ii) a nonlinear solver.

For sake of simplicity, we assume an implicit Euler discretization in time, i.e. , where denotes the known value from previous timestep , and denotes the unknown value at the new timestep , i.e. . We also fix the first argument of . Then, scheme (23) reduces to a form which is formally equivalent to scheme (20), namely: find so that

 ∫ΩU1⋅V dV+c(U0,U1,V)+a(U1,V)=b(V)+∫ΩU0⋅Lt0t1V dV∀V∈VXk––(t1). (24)

The linearization point may be either set as , which results in a semi-implicit formulation. Alternatively, one can utilize fully implicit approach and iterate over equation (24), so that or employ a Newton method.

On the right-hand-side, the linear operator performs a lifting from to . This requires that the corresponding part of the right-hand-side of (24) must be integrated on the old cut-cell mesh at time .

The lifting is defined as follows:

###### Definition 6 (temporal XDG space lifting).

The cut-cell mesh and at times and have equal topology, if, and only if for each cut-cell one has

 |Kj,s(t0)|>0 ⇒ |Kj,s(t0)∩Kj,s(t1)|>0% and |Kj,s(t1)|>0 ⇒ |Kj,s(t1)∩Kj,s(t0)|>0.

On meshes with equal topology, the lifting operator

 PXk(Kh,t1)∋u1↦Lt0t1u1=:u0∈PXk(Kh,t0)

is uniquely defined requiring polynomial equality on the common domain of old and new cut-cell, i.e. by the property

 u0∣∣Kj,s(t0)∩Kj,s(t1)=u1∣∣Kj,s(t0)∩Kj,s(t1)∀j,s.

Since is a product of spaces , the lifting naturally extends to . Obviously, it cannot always be ensured that the cut-cell meshes for and have equal topology. However, this important property can be achieved through cell agglomeration, which is addressed in the next section.

### 2.3 Cell agglomeration

The motivation for aggregation/agglomeration meshes is three-fold: removal of small cut-cells, avoiding topology changes in the cut-cell mesh for a single time-step and formulation of multigrid methods without the usual hierarchy of meshes, cf. Section 3 and Figure 1).

Formally, the aggregation mesh is introduced by means of graph theory:

###### Definition 7 (Graph of a numerical mesh).

Let be a numerical mesh; For , the set is called a logical edge if, and only if . Furthermore, let be the set of all logical edges. Then, the pair forms an undirected graph in the usual sense.

###### Definition 8 (Aggregation maps and meshes).

Let be an aggregation map and be the nodes of a connected component of . Note that might consist of only a single element, i.e. an isolated node, which is called an non-aggregated cell with respect to . The aggregate cell is defined as the union of all cells, i.e. (Rem.: in order to ensure that the aggregate cell is again a simply connected, open set, one has to take the closure of each cell first and then subtract the boundary, therefore we define a modified union as .)

For the aggregation mesh is the set of all aggregate cells which can be formed w.r.t. .

Based upon the aggregation mesh, an aggregated XDG space can be defined:

###### Definition 9 (Aggregated XDG space).

For some agglomeration map we define the agglomerated XDG space as:

 (25)

Obviously, the agglomerated XDG space is a sub-space of the original space, i.e. .

#### Temporal discretization and stabilization against small cut-cells

As already noted above, in order to discretize temporally evolving systems such as Eq. (24), one has to ensure that the cut-cell mesh at time steps and have equal topology in order to obtain a well-defined method. Otherwise, the required lifting operator (see Definition 6) is undefined. For multi-step schemes, which involve multiple time steps, the topology has to be equal for all time steps; the same holds for the intermediate steps of Runge-Kutta schemes, cf. Kummer et al. (2017).

Furthermore, since the interface position is arbitrary, cut-cells can be arbitrarily small, i.e. its volume fraction w.r.t. the background cell can be small. This leads e.g. to large penalty parameters in the SIP form (21) which is known to cause undesirably high condition numbers of the discretized system.

Therefore, instead of solving the variational system on the space , which is induced by the cut-cell mesh, one employs an XDG space on an appropriately agglomerated cut-cell mesh. A valid agglomeration map for these purpose must meet the following requirements:

• The meshes and have the same topology.

• All cut-cells with a volume fraction are agglomerated

• There is no agglomeration across species, i.e. there exists no edge in .

The formulation of an algorithm that constructs an agglomeration map which fulfills the properties noted above is left to the reader. It may consist of a loop over all cut-cells. The cut-cell must be agglomerated if it is a new cell (i.e. and ) or a vanished cell (i.e. and ) or if its volume fraction is below the threshold . A decent value for lies in the range of 0.1 to 0.3, cf. Kummer (2017); Müller et al. (2017). In our implementation, such a cell is agglomerated to its largest neighbor cell in the same species.

#### The final system

Instead of solving the generic variational system (20) on the space , the aggregated space

 VX,α,Δtk––:=Dv∏γ=1PXkAα,t0,t1(Kh,t1) (26)

is used for discretization. In comparison to (19), the dependence on time is dropped since w.l.o.g. all temporally evolving systems are solved for the ‘new’ time step .

Hence, the final discretization reads as: find so that

 a(U,V)=b(V)∀V∈VX,α,Δtk––. (27)

Over the course of multiple time-steps, the agglomeration graph typically changes. After each time step is complete, the solution is injected into the non-agglomerated space. For the next time step it is projected back (in an -sense) onto the (potentially different) agglomerated space to serve as an initial value.

## 3 Multigrid solvers

The remaining of this paper is dedicated to the solution of the linear system (27). The solvers presented are use a combination of aggregation- and p-multigrid. Aggregation multigrid can be seen as a combination of conventional h-multigrid and algebraic multigrid methods: there is still an underlying mesh of polyhedral cells. Due to these cells, the flexibility is comparable to algebraic multigrid, see Figure 1.

This section is organized as follows: first, the aggregation multigrid framework will be unified with the XDG method established so far (section 3.1). Next, the construction of a basis for the nested sequence of approximation spaces will be laid out, in order tor transfer the basis-free formulation (27) into a matrix form. Finally, the specific combination of multigrid algorithms are presented (3.3).

### 3.1 Aggregation multigrid for XDG

The starting point of the aggregation multigrid is a sequence of aggregation maps

 ∅=A1⊂A2⊂…⊂AΛ⊂Edg(Kh) (28)

on the background mesh. Note that the injection/projection operator between aggregation grid are quite expensive to compute. Therefore, it is beneficial to compute them initially and only update when necessary in cut-cells. Since these are defined on the background mesh, in order to be precomputed, one cannot directly apply these aggregation maps onto the cut-cell mesh. Therefore, aggregation maps from the background mesh must be mapped onto the cut-cell mesh:

###### Definition 10 (Mapping of an aggregation map onto a cut-cell mesh).

Given is an aggregation map on the background mesh . The corresponding aggregation map on the cut-cell mesh is formed from all edges by duplicating them for each species, i.e. from edges

 {Kj,s(t),Kl,s(t)}∀s∈{A,B}, if |Kj,s(t)|>0 and |Kl,s(t)|>0.

Note that this construction avoids aggregation across the interface , as illustrated in Figure 2.

Through the mapping of an aggregation map , a sequence of aggregated XDG spaces is induced,

 VX,α,Δtk––=:VX,1k–– R≥ VX,2k–– R≥ … R≥ VX,Λk––

where the space on mesh level can be defined as

 VX,λk––:=Dv∏γ=1PXkA′λ(Kh,t1) (29)

and the aggregation map is the union of the predefined multigrid aggregation sequence (see 28) and the aggregation map to stabilize small cut-cells and prevent temporal toplology changes, i.e.

 A′λ:=(Aλ)X(t1)∪Aα,t1,t0. (30)

### 3.2 Basis representation and indexing

Up to this point, the XDG method is formulated in a variational, coordinate-free form. In order to notate solver algorithms in a typical matrix-notation, a basis representation of the XDG space is required.

#### A Basis of VX,λk––

The elements of the basis are written as , with the following index conventions:

• is the index related to the background cell, where for aggregation cells one picks the minimum cell index of all aggregated background cells on mesh level as a representative,

• is the variable index (e.g. for Navier-Stokes, corresponds to the velocity components and corresponds to pressure),

• is the species index and

• is the DG-mode index, where is the dimension of the polynomial space up to degree .

The row-vector of all XDG basis functions, on mesh level is written as

 (→ΦX,λj,γ,s,n)j,γ,s,n=:→Φ––X,λ. (31)

(Here, we skik the specification of all valid combinations of for sake of compactness.) In the implementation, this basis is constructed from a basis on the space , with in the following way: first, a basis of the aggregated space is created. Since , one can express this basis in terms of the original basis, i.e.

 Φ––λ=Φ––⋅Q––––λ,

with a suitable matrix . It can be derived from an Ansatz which projects a polynomial basis on the bounding box of an aggregation cell onto the background cells, for details see ref. Kummer (2017).

The matrix obviously is a prolongation operator. If both, and are orthonormal, the mapping

 Pk(Kh)∋Φ––⋅u––↦Φ––λ⋅((Q––––λ)Tu––)∈Pk(Agg(Kh,Aλ))

is a projector in the -sense.

For computational efficiency and in order to save memory, the matrix should not be stored. Instead, is expressed in terms of , i.e.

 Φ––λ+1=Φ––λ⋅R––––λ. (32)

For orthonormal bases, is the prolongation matrix from the coarse to fine mesh, while is a restriction, resp. projection matrix from fine to coarse mesh. For a specific aggregation cell, it can be computed by projection of polynomials onto one representative background cell for each part of the aggregation cell.

Second, the basis of the XDG space is constructed: the basis elements of are expressed as

 →ΦX,λj,γ,s,n(→x):=→eγNkγ∑m=1Φλj,m(→x)χs(t)(→x) Sλmn. (33)

Here is the standard basis vector in and

denotes the characteristic function for set

. The matrix provides a re-orthonormalization of the basis functions in cut-cells and can be obtained e.g. through a Cholesky factorization.

Finally, through a combination of multigrid projection matrices and re-orthonormalization matrices , one obtains the representation

 →Φ––X,λ+1=→Φ––X,λ⋅R––––X,λ. (34)

Formally, is a matrix product of and . The notation of its exact shape is rather technical and therefore skipped. Mainly, since e.g. are vector-valued an advanced indexing notation is required, which is introduced below.

#### Multi-index mapping

In order to extract sub-matrices and sub-vectors which correspond to certain cells, variables and DG modes, a sophisticated index notation is required. A single basis element of can be associated with a multi-index . One may think of as a bijection between all valid combinations of , , and and the set . We use a notation where the mapping is employed to select sub-sets of , e.g.

 m(j,−,−,≤Nk––):={m(j,dv,s,n); ∀dv, ∀s present in cell j, ∀n≤Nkdv}.

Such sets will be used to notate sub-matrices or sub-vectors, similar to the typical MATLAB notation.

#### Agglomeration Algebra

Given that a basis is established, the generic system (27) which searches for a solution can be transfered to an equivalent matrix formulation

 M––––––λ U––=b–λ (35)

with

 Mλm(j,dv,s,n) m(l,ev,r,m)=a(→ΦX,λl,ev,r,m,→ΦX,λj,dv,s,n) and bλm(j,dv,s,n)=b(→ΦX,λl,ev,r,m). (36)

Through restriction and prolongation matrices, one obtains the relation

 M––––––λ+1=(R––––X,λ)T M––––––λ R––––X,λ and b–λ+1=(R––––X,λ)T R––λ, (37)

as usual in multigrid methods.

### 3.3 Preconditioners and Solvers

Within this section the discussion is focused on a single mesh level , hence the level index is dropped. On this grid level, let be the vector space dimension, i.e. . The (approximate) solution and right-hand-side (RHS) vector of the system (35) are denoted as , respectively.

In this section, a series of algorithms is introduced.

• A p-multigrid algorithm (Algorithm 11), that operates on a single mesh level, which can be used as a pre-conditioner for a standard GMRES solver.

• The additive Schwarz algorithm (Algorithm 12), which uses p-multigrid as a block solver. It contains a minor modification to its original form, which we found helpful in decreasing the number of iterations and is therefor also presented here.

• Finally, the orthonormalization multigrid algorithm (Algorithm 14) which employs the additive Schwarz as a smoother and a residual minimization (Algorithm 13) to ensure non-increasing residual norms.

#### A p-multigrid pre-conditioner on a single mesh level

For DG or XDG methods, even without aggregation meshes, a p-multigrid seems to be a reasonable idea: At first, a sub-matrix which corresponds to low order DG modes is extracted. Since the number of degrees-of-freedom for this is typically low, a sparse direct solver can be used. The degrees of freedom which correspond to higher order DG modes are then solved locally in each cell. Since Algorithm

11 will also be used as a block solver for the additive Schwarz method (Algorithm 12), an optional index-set which restricts the solution to a part of the mesh can be provided.

###### Algorithm 11 (p-multigrid pre-conditioner).

is computed as follows:

A right-hand-side , a vector of DG polynomial degrees which separates low from high order modes, for each variable. Optionally, an index-set of sub-cells, denoting the block to solve.
An approximate solution
initialize approximate solution
Indices for low-order modes in cells
extract low-order system
extract low-order RHS
Solve: usually using a sparse direct solver
accumulate low-order solution
for all  do loop over cells…
Indices for high-order modes in cell
extract high-order system in cell
extract high-order RHS in cell
Solve: using a dense direct solver
accumulate local high-order solution
end for

Since a preconditioner typically has to be applied multiple time, it is essential for performance to use a sparse direct solver for the low order system which is able to store the factorization and apply multiple right-hand sides. Numerical tests hint that usually a single-precision solvers is sufficient as long as the residuals are computed in double precision. The PMG aglorithm presented above can directly be used as a preconditioner, e.g. for GMRES. However, since only two multigrid levels are used, its application typically is only reasonable up to medium-sized systems.

In the scope of this work, is also used as a block solver for an additive Schwarz method. For a general discussion on Schwarz methods, we refer to the review article of Gander Gander (2008). In our implementation the blocking is determined on the level of cells, i.e. on the basis of the mesh . The METIS Karypis and Kumar (1998) software library is used to partition ; the number of partitions is determined so that a single partition contains roughly 10,000 degrees-of-freedom. After the partitioning by METIS, each partition is enlarged by its neighbor cells to generate the overlap layer that is typically used with Schwarz methods.