1 Introduction and Motivation
BoSSS^{1}^{1}1 https://github.com/FDYdarmstadt/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 timedependent domains.
One important field of application for the XDG method are moving geometries: Almost all numerical methods for PDEs require the laborintensive and thus expensive (semimanual) creation of numerical grids or meshes. This is particularly cumbersome if timedependent 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 sideproduct, 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 meshmotion techniques like arbitrary LagrangianEulerian (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 crossreference. Since the source code is constantly being developed, the validity of such references would be shortlived.
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 subdomains 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 piecewise Poisson problem
(2) 
with a discontinuous diffusion coefficient
(3) 
Incompressible multiphase flows
In a transient stetting, the phases are usually considered to be timedependent too, i.e. one has the decomposition . Using the same notation as introduced above, we introduce the incompressible twophase NavierStokes equation for material interfaces as
(4) 
with piecewise constant density and viscosity for both phases, i.e.
(5) 
Furthermore, denotes surface tension and denotes the mean curvature of . In a twophase setting, we also assume that the interface does not touches the boundary, i.e. . If this is not the case, a threephase 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 LevelSet:
Note that problem (4) is incomplete without a specification of the interface motion. In order to track the individual domains, a sufficiently smooth LevelSet field is introduced; for sake of simplicity, might also be called ‘LevelSet’. Then, timedependent domains , and the interface can be described as
(6)  
(7)  
(8) 
On , must therefore comply with the LevelSet equation
(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
(10) 
and the mean curvature can be computed by Bonnets formula as
(11) 
If the LevelSet 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
(12) 
Other applications
The XDG framework can be seen as a multipurpose technology. One obvious application is the use as an immersed boundary method (IBM), e.g. in order to circumvent the (labourextensive) 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 NewtonEuler 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 LevelSet.
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 threephase contact line between fluids and and the solid. In order to allow motion of the contact line , the noslip boundary condition has to be relaxed, yielding the socalled NavierSlip boundary condition. It is notable that the XDGimplementation does not need further manipulation of the contact line velocity or contact angle . Furthermore, XDG allows the implementation of the generalized Navierslip boundary condition with
(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 socalled 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 timedependent problems by Chessa and Belytschko Chessa and Belytschko (2004). The first XFEM for twophase flows was presented by Groß and Reusken Groß and Reusken (2007a, b). With the socalled 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 cutcell 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 nodelevel performance. A very important feature is dynamic load balancing, which allows redistribution 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, databasecentric 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:
(14) 
2.1 Extended Discontinuous Galerkin and LevelSet
In order to yield a welldefined interface , the LevelSetfield must be sufficiently smooth. Precisely, we assume to be almosteverywhere in . For certain steadystate 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 (Cutcell mesh).
Timedependent cutcells are given as
(15) 
The set of all cutcells form the timedependent cutcell mesh .
Note that we refer to the ‘original cells’ as background cells, in contrast to the timedependent cutcells . By , resp. a notation for an arbitrary phase is introduced, i.e. can be either or .
XDG is essentially a DG method on cutcells, i.e. one can define the XDG space as follows:
Definition 4 (XDG space).
We define:
(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. outervalue of a field are defined as:
Then, the jump and average value operator are defined as
(17)  
(18) 
For the implementation of an XDGmethod, accurate numerical integration on the cutcells
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 twophase Navier Stokes problem (4) has to be coupled with a second method to compute the evolution of the interface. Since the respective LevelSet 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 NavierStokes equation, with
dependent variables which are not necessarily discretized with the same polynomial degree. We therefore define the degreevector
; and introduce an abbreviation for the function space of test and trial functions,(19) 
Then, the discrete version of some linear problem, for a fixed time formally read as: find so that
(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 lefthandside of Eq. (20). Here, denotes the broken gradient, where differentiation at the jumps on is excluded. The linear form on the righthandside of Eq. (20) is given as
(22) 
The SIP factor is known to scale as , where is a local length scale of the agglomerated cutcell, 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 rulesofthumb 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 NavierStokes Eq. (4)
Due to the nonlinearity in the convectional term, the NavierStokes system is usually split up into the linear Stokes part and the nonlinear convection operator . It is known that, in order to obtain an infsup 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 NavierStokes equation formally reads as: find so that
(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
(24) 
The linearization point may be either set as , which results in a semiimplicit formulation. Alternatively, one can utilize fully implicit approach and iterate over equation (24), so that or employ a Newton method.
On the righthandside, the linear operator performs a lifting from to . This requires that the corresponding part of the righthandside of (24) must be integrated on the old cutcell mesh at time .
The lifting is defined as follows:
Definition 6 (temporal XDG space lifting).
The cutcell mesh and at times and have equal topology, if, and only if for each cutcell one has
On meshes with equal topology, the lifting operator
is uniquely defined requiring polynomial equality on the common domain of old and new cutcell, i.e. by the property
Since is a product of spaces , the lifting naturally extends to . Obviously, it cannot always be ensured that the cutcell 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 threefold: removal of small cutcells, avoiding topology changes in the cutcell mesh for a single timestep 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 nonaggregated 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 subspace of the original space, i.e. .
Temporal discretization and stabilization against small cutcells
As already noted above, in order to discretize temporally evolving systems such as Eq. (24), one has to ensure that the cutcell mesh at time steps and have equal topology in order to obtain a welldefined method. Otherwise, the required lifting operator (see Definition 6) is undefined. For multistep schemes, which involve multiple time steps, the topology has to be equal for all time steps; the same holds for the intermediate steps of RungeKutta schemes, cf. Kummer et al. (2017).
Furthermore, since the interface position is arbitrary, cutcells 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 cutcell mesh, one employs an XDG space on an appropriately agglomerated cutcell mesh. A valid agglomeration map for these purpose must meet the following requirements:

The meshes and have the same topology.

All cutcells 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 cutcells. The cutcell 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
(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
(27) 
Over the course of multiple timesteps, the agglomeration graph typically changes. After each time step is complete, the solution is injected into the nonagglomerated 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 pmultigrid. Aggregation multigrid can be seen as a combination of conventional hmultigrid 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 basisfree 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
(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 cutcells. Since these are defined on the background mesh, in order to be precomputed, one cannot directly apply these aggregation maps onto the cutcell mesh. Therefore, aggregation maps from the background mesh must be mapped onto the cutcell mesh:
Definition 10 (Mapping of an aggregation map onto a cutcell mesh).
Given is an aggregation map on the background mesh . The corresponding aggregation map on the cutcell mesh is formed from all edges by duplicating them for each species, i.e. from edges
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,
where the space on mesh level can be defined as
(29) 
and the aggregation map is the union of the predefined multigrid aggregation sequence (see 28) and the aggregation map to stabilize small cutcells and prevent temporal toplology changes, i.e.
(30) 
3.2 Basis representation and indexing
Up to this point, the XDG method is formulated in a variational, coordinatefree form. In order to notate solver algorithms in a typical matrixnotation, a basis representation of the XDG space is required.
A Basis of
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 NavierStokes, corresponds to the velocity components and corresponds to pressure),

is the species index and

is the DGmode index, where is the dimension of the polynomial space up to degree .
The rowvector of all XDG basis functions, on mesh level is written as
(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.
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
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.
(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
(33) 
Here is the standard basis vector in and
denotes the characteristic function for set
. The matrix provides a reorthonormalization of the basis functions in cutcells and can be obtained e.g. through a Cholesky factorization.Finally, through a combination of multigrid projection matrices and reorthonormalization matrices , one obtains the representation
(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 vectorvalued an advanced indexing notation is required, which is introduced below.
Multiindex mapping
In order to extract submatrices and subvectors 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 multiindex . 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 subsets of , e.g.
Such sets will be used to notate submatrices or subvectors, 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
(35) 
with
(36) 
Through restriction and prolongation matrices, one obtains the relation
(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 righthandside (RHS) vector of the system (35) are denoted as , respectively.
In this section, a series of algorithms is introduced.

A pmultigrid algorithm (Algorithm 11), that operates on a single mesh level, which can be used as a preconditioner for a standard GMRES solver.

The additive Schwarz algorithm (Algorithm 12), which uses pmultigrid 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.
A pmultigrid preconditioner on a single mesh level
For DG or XDG methods, even without aggregation meshes, a pmultigrid seems to be a reasonable idea: At first, a submatrix which corresponds to low order DG modes is extracted. Since the number of degreesoffreedom 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 indexset which restricts the solution to a part of the mesh can be provided.Algorithm 11 (pmultigrid preconditioner).
is computed as follows:
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 righthand sides. Numerical tests hint that usually a singleprecision 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 mediumsized systems.
An additive Schwarz method
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 degreesoffreedom. 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.
Algorithm 12 (Additive Schwarz).
is computed as follows:
Within the scope of this work, the additive Schwarz solver is used as a smoother for the multigrid method. As such, we introduced a minor modification to the original formulation: in the last line of Algorithm 12, the solution is divided by the number of blocks which contain a specific cell. This damping of the approximate solution in the overlapping regions has shown to improve the number of iterations when
Comments
There are no comments yet.