# Multilevel Schwarz preconditioners for singularly perturbed symmetric reaction-diffusion systems

We present robust and highly parallel multilevel non-overlapping Schwarz preconditioners, to solve an interior penalty discontinuous Galerkin finite element discretization of a system of steady state, singularly perturbed reaction-diffusion equations with a singular reaction operator, using a GMRES solver. We provide proofs of convergence for the two-level setting and the multigrid V-cycle as well as numerical results for a wide range of regimes.

## Authors

• 4 publications
• 8 publications
04/29/2020

### Optimization of two-level methods for DG discretizations of reaction-diffusion equations

We analyze and optimize two-level methods applied to a symmetric interio...
04/22/2021

### Efficient discretization and preconditioning of the singularly perturbed Reaction Diffusion problem

We consider the reaction diffusion problem and present efficient ways to...
12/25/2020

09/18/2019

### First-order system least squares finite-elements for singularly perturbed reaction-diffusion equations

We propose a new first-order-system least squares (FOSLS) finite-element...
03/26/2019

### A Multilevel Approach for Trace System in HDG Discretizations

We propose a multilevel approach for trace systems resulting from hybrid...
06/23/2020

### An experimental comparison of a space-time multigrid method with PFASST for a reaction-diffusion problem

We consider two parallel-in-time approaches applied to a (reaction) diff...
03/16/2021

### BDDC Algorithms for Advection-diffusion problems with HDG Discretizations

The balancing domain decomposition methods (BDDC) are originally introdu...
##### 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

In this paper, we present analysis and numerical experiments of two-level Schwarz preconditioners and their multilevel versions for an interior penalty discontinuous Galerkin finite element discretization for a system of reaction-diffusion equations. Our focus is on the singularly perturbed case, where the reaction system has an inertial subspace. Our estimates are robust with respect to the parameters of the system and the experiments confirm efficiency of the method.

Reaction-diffusion systems arise in a variety of physical, chemical and biological contexts. One particular area where these models are widely used is radiation transport, where the reaction-diffusion equation is an approximation of Boltzmann’s linear transport equation that becomes relevant in the so called diffusive regimes, which are characterized by small mean free paths compared to the size of the domain. In these regimes the transport equation is nearly singular and its solution in the interior of the computational domain is close to the solution of a reaction-diffusion equation [16].

We employ the interior penalty discontinuous Galerkin (IP-DG) method to discretize the singularly perturbed reaction-diffusion system in steady state. IP-DG [3, 17, 4, 2, 21] methods are particularly interesting to solve reaction-diffusion equations since they do not produce oscillations near the boundaries in singularly perturbed problems. Using this discretization the reaction operator involves only volume integrals with no coupling between cells. Therefore, all its contributions are included inside the local subspaces of Schwarz methods, which motivates our choice of preconditioner.

We solve the discrete problem with a GMRES solver using multilevel preconditioners with nonoverlapping Schwarz smoothers, effectively solving a full reaction-diffusion problem in each cell (see [15]). Convergence estimates for such methods applied to diffusion problems have been developed in [12]. There, it is assumed that the subdomains defining the decomposition of the fine space are unions of coarse cells. Recently, an extension has been shown covering the case of cell-wise subdomains [10]. However, its application does not extend to quadrilaterals and hexahedra, since the proof uses

nonconforming interpolant and enriching operators for simplices

[7]. We provide an extension for quadrilaterals and hexahedra without such restrictions, which is based on the original proof in [12].

Robust solvers for singularly perturbed reaction-diffusion equations and systems of equations have been studied in [5, 13, 19], some of them use Schwarz methods. But, they all share the assumption of a positive definite reaction system, thus excluding the presence of an inertial subspace. On the other hand, inertial subspaces are important in applications, since they either describe the stationary limit of a reaction system, or in the case of a stiff system with strong separation of time scales, they approximate the slow-changing quantities. Thus, we focus on cases where the coupling matrix is singular. Hence, the reaction-diffusion operator acts as pure diffusion on the inertial subspace, and as reaction dominated on its complement.

Our main results are the proof of the stable decomposition shown in lemma 3.2 to obtain convergence estimates for two-level preconditioners, and the multigrid V-cycle preconditioner estimate in theorem 3.4.

The paper is structured as follows: in section 2 we introduce the continuous problem andthe IP-DG discretization. In section 3 we develop two-level Schwarz and multigrid preconditioners and prove convergence estimates. Finally, in section 4 we demonstrate the efficiency of the proposed methods with experimental results.

## 2. Model problem

We consider the following system of steady state reaction-diffusion equations with a singularly perturbed reaction term

 (1) −∇⋅(ηg∇ug)+1εG∑g′=1(σgg′ug−σg′gug′)=fg in Ω, with g=1…G,

where is the group index identifying each reacting substance, is the diffusion coefficient for each group , is a perturbation parameter defining the relative size of the reaction with respect to the diffusion term, is a convex polyhedral domain in with and is a known source.

The equation is provided with the boundary conditions

 ug=0 on Γ, with g=1…G,

where is the boundary of .

We assume and , for all and there exists such that in . Furthermore, we assume that the reaction matrix is symmetric, in the sense that , .

We introduce the Hilbert spaces

 V=(H10(Ω))G , H=(L2(Ω))G,

where is the standard Sobolev space with zero boundary trace. They are provided with inner products

 (2) (u,v)V=G∑g=1(ηg∇ug,∇vg)L2(Ω) , (u,v)H=G∑g=1(ug,vg)L2(Ω),

and norms

 ∥u∥2V=(u,u)V , ∥u∥2H=(u,v)H.

The weak form of problem (1) is: find such that

 (3) A(u,v)=(f,v)H,

where and the bilinear form is

 (4) A(u,v)=G∑g=1∫Ωηg∇ug⋅∇vgdx+1εG∑g=1G∑g′=1∫Ω(σgg′ugvg−σg′gug′vg)dx=(D∇u,∇v)H+1ε(Σu,v)H=(u,v)V+1ε(Σu,v)H.

The second line uses the vector notation

 u =(u1,…,uG)⊺, v =(v1,…,vG)⊺, D =diag(η1,…,ηG), and Σ= ⎛⎜ ⎜⎝∑g≠1σ1g…−σG1………−σ1G…∑g≠GσGg⎞⎟ ⎟⎠.

According to our assumptions, the reaction matrix is a symmetric, weakly diagonally dominant singular M-matrix111We use the term singular M-matrix, following the terminology in [14, p. 119], to denote a matrix that can be expressed as , where all the elements in are nonnegative, is equal to the maximum of the moduli of the eigenvalues of , and

is an identity matrix.

with zero column and row sum. By the Perron-Frobenius theorem, this implies is singular with rank less than

and by the Geršgorin circle theorem, all eigenvalues are nonnegative.

Physically, the properties of ensure substance conservation and the absence of sinks inside the domain. In a radiation transport context, this implies that the system can have no particle absorption and particles only disappear when they reach the boundary. The presence of absorption would imply all eigenvalues are positive and would be invertible.

Under the assumptions on the parameters of equation (1), the bilinear form is continuous and -coercive relatively to (see [9, §2.6]), i. e. there exist constants such that

 A(u,u)≥γA∥u∥2V , A(u,v)≤CA∥u∥V∥v∥V.

where we remark that even though is independent of , is not.

From Lax-Milgram’s theorem, the variational problem admits a unique solution in .

### 2.1. Discrete problem

We apply a IP-DG discretization to the bilinear form [2]. Let be a subdivision of the domain into quadrilaterals or hexahedra , such that each cell is described by a -linear mapping from the reference cell onto itself. Conformity of the faces of mesh cells is not required, but we assume local quasi-uniformity and shape regularity in the sense that the Jacobians of and their inverses are uniformly bounded.

Let

be the space of tensor product polynomials of degree up to

in each coordinate direction. Then, define the mapped space on the cell as the pull-back of functions under . The vector-valued, discontinuous function space is then defined as

 (5) Vh={v∈H∣∣v|κ∈QGp(κ)}.

Let be the set of all interior faces of the mesh and the set of all boundary faces. Let be two mesh cells with a joint face , and let and be the traces of functions on from and respectively. On face , we define the averaging operator as

 {{u}}=u1+u22.

We introduce the following definition of mesh integrals

 ∫Thudx=∑κ∈Th∫κudx

and integrals over and are defined accordingly. The interior penalty (IP) bilinear form for the scalar Laplacian, as described in [2], reads

 (6) αh(u,v)=∫Th∇u⋅∇vdx−∫FIh∪FBh2({{un}}⋅{{∇v}}+{{∇u}}⋅{{vn}})ds+∫FIh∪FBh4δIPh{{un}}⋅{{vn}}ds

where is the minimum cell diameter adjacent to the face, and . We have replaced the jump operator used in [2] for the equivalent expression: . Coercivity and continuity are proven in [2] under the assumption that is sufficiently large. We will assume in the following that this holds true.

We then define the discrete bilinear form, including the diffusion coefficients follows

 (7) ah(u,v)=∫ThD∇u⋅∇vdx+∫FIh∪FBh4δIPh{{D(un)}}⋅{{vn}}ds−∫FIh∪FBh2({{un}}⋅{{D∇v}}+{{D∇u}}⋅{{vn}})ds.

Under the assumptions made in the previous sections and sufficiently large, is coercive and continuous.

Using (7), our IP-DG formulation for the singularly perturbed reaction diffusion problem reads: Find such that

 (8)

We observe that given the non-negativeness of , the coercivity constant for our problem coincides with the Laplacian case while the continuity constant is now dependent on . In order to obtain a robust solver we precondition the problem to be able to bound the spectral radius of the preconditioned system independently of .

Finally, using a standard basis for the local finite element spaces on each cell and concatenating, we obtain the linear system

 Au=f,

where and are the coefficient vector of the representation of and respectively in terms of the chosen basis.

## 3. Preconditioners

In this section we provide details on our solver and preconditioner choice, as well as the technical tools needed for the numerical analysis of the preconditioned system.

It is known that the convergence of the preconditioned conjugate gradient method for symmetric real operators depends on the condition number of the preconditioned matrix only. Thus, if we find a preconditioner such that the this condition number is independent of and of the parameters of the equation, the number of iterations required for convergence to a certain error is independent of them as well. We will estimate the condition number of the additive Schwarz method by estimating the smallest and largest eigenvalues and as

For the rest of the preconditioners, we will estimate the norm of the error propagation operator of a preconditioned Richardson iteration.

### 3.1. Schwarz preconditioners

We choose Schwarz preconditioners for which there is a well-known framework and theory for symmetric positive definite problems (see [20, 8, 18, 12]). The following sections provide the definitions needed to prove convergence estimates in an abstract formulation.

Let for be Hilbert spaces with norms , where is used to denote the so-called coarse space in a domain decomposition context. For , let

 R⊺j:Vj→Vh

denote prolongation operators for which there holds

 R⊺jVj⊂V , and V=J∑j=0R⊺jVj, for j=0,1,2,…,J.

Here is the range of the linear operator .

Associated with each local space for , we introduce local discrete bilinear forms , defined on , as the restriction of global discrete bilinear form on , with .

For the coarse space we use the rediscretization of the problem on the coarse mesh, namely a bilinear form with a penalty parameter inversely proportional to the diameter of the coarse cells , instead of the inherited coarse space obtained by restriction to .

The abstract theory is based on three main assumptions:

###### Assumption 3.1 (Stable decomposition).

The spaces are said to provide a stable decomposition if there exists a constant such that each admits a decomposition

 v=J∑j=0R⊺jvj,

with such that

 J∑j=0∥vj∥2Aj≤ CV∥v∥2Ah,

where and accordingly.

If , admits a stable decomposition without including the coarse space as follows (see [20, p.49])

 J∑j=1∥vj∥2Aj≤ CV∥v∥2Ah.
###### Assumption 3.2 (Strengthened Cauchy-Schwarz inequality).

There exist constants for such that

 Ah(R⊺ivi,R⊺jvj)≤

We will denote the spectral radius of by .

###### Assumption 3.3 (Local stability).

There exists such that

 Ah(R⊺jvj,R⊺jvj)≤ωAj(vj,vj) ∀vj∈Vj.

We now introduce a set of projection-like operators for . These projection-like operators will serve as the building blocks for the construction of Schwarz methods. For any fixed , define by

 Aj(˜Pjv,wj):= Ah(v,R⊺jwj), ∀wj∈Vj.

We note that the well posedness of the global problem ensures is well defined for . To map the elements of into the global discrete space , we employ the prolongation operator and define the composite operator

 Pj:=R⊺j∘˜Pj, for j=0,1,2,…,J.

Trivially, we have for .

After these preparations, we can write the operator preconditioned with the additive Schwarz method as

To facilitate the comprehension of the method with respect to its implementation, we write the additive operator in a more explicit form. We use the operator notation for the bilinear forms and to obtain the following expression for the local projections

 Aj˜Pjv:= RjAhv, ∀v∈Vh.

Thus,

 ˜Pj= A−1jR⊺jAh, and Pj=R⊺jA−1jRjAh,

While the additive version applies all subspace corrections at once and adds them in the end, the multiplicative version applies them successively. It can be defined easily by the error propagation operator

 Emu=(I−PN)∘(I−PN−1)∘⋯∘(I−P0),

where denotes the identity operator on . Using we define the multiplicative Schwarz preconditioner as

 Pmu=I−Emu,

where denotes the identity operator on .

Finally, we consider the symmetric, hybrid version, which is additive with respect to the subdomain spaces, but applies the coarse grid correction in a multiplicative way:

 (9) Phy=I−(I−N∑i=1Pi)(I−P0)(I−N∑i=1Pi).

Following, we will prove convergence estimates for the operators , and . For we estimate the condition number, for we bound the error operator of a preconditioned Richardson iteration and for we defer the proof to section 3.3, where we study multigrid preconditioners, from which is a special case.

We use the general abstract convergence theory of Schwarz methods given in [20, §2]. We quote the convergence results below.

###### Theorem 3.1.

Let the assumptions 3.1, 3.2 and 3.3 hold, then the following bounds hold for the additive Schwarz preconditioned system

Where and are the smallest and largest eigenvalues of the preconditioned system, respectively.

###### Proof.

See [20, §2.3]. ∎

###### Theorem 3.2.

Let the assumptions 3.1, 3.2 and 3.3 hold, then the following bounds hold for the hybrid Schwarz preconditioned system

 ∣∣Ah([I−P% hy]v,v)∣∣≤cMGAh(v,v), ∀v∈Vh,

where is a constant independent of and .

###### Proof.

We defer this proof to section 3.3, as it is a special case of a multigrid preconditioner and as such its convergence estimate is given in theorem 3.4. ∎

The multiplicative operator is not symmetric and we will consider a simple Richardson iteration applied to the corresponding preconditioned system and provide an upper bound for the norm of the error propagation operator.

###### Theorem 3.3.

Let the assumptions in definitions 3.1, 3.2 and 3.3 hold, then the following bounds hold for the multiplicative Schwarz preconditioned system

 ∥∥Emu∥∥≤1−2−ω(2max{1,ω2}ρ2(Θ)+1)CV≤1
###### Proof.

See [20, §2.3]. ∎

### 3.2. Application to the discrete problem

In this section we define the Schwarz method for the discrete problem in equation (8) and verify that assumptions 3.1, 3.2 and 3.3 apply.

After enumerating the cells for , we choose the local spaces , together with the coarse space , defined on . We remark that we are using a nonoverlapping subdivision order to define the direct decomposition , where is the simple injection. Similarly, for , if and zero otherwise.

Following, we list three standard results from [12] that we need for our proof.

For any , there holds the trace inequality (see [12, Lemma 3.1])

 (10) ∥v∥2H(∂D)≤c[1H∥v∥2H(D)+H∥v∥2V(D)].

Suppose is a convex domain. For any , let be the average value of over . Then we can write a Poincaré inequality as follows (see [12, Lemma 3.2])

 ∥v−¯¯¯v∥H(D)≤c diam(D)∥u∥V(D) on D.

In particular, if

 (11) ∥v−¯¯¯v∥H(D)≤cH∥u∥V(D) on D.

Let , let , , be given (uniquely) by , . Then the following identity holds (see [12, Lemma 3.3])

 (12) ah(v,w)=J∑j=1aj(vj,wj)+I(v,w),

where comprises all terms located outside the block diagonal of the bilinear from , connecting different subdomains.

We then obtain the following interface estimate for cell-wise subdomains

###### Lemma 3.1.

There exists a constant such that

 |I(v,v)|≤c⎡⎣1h2∑K∈Th∥v∥2H(K)+ah(v,v)⎤⎦.
###### Proof.

We extend the result in [12, Lemma 4.3]. The following estimate, from [12, Eq. (4.20)], holds when using cell wise subdomains

where is the -inner product on the faces of cell of the fine mesh.

Using the trace inequality from [12, Eq. (3.9)], we obtain

 |I(v,v)|≤c⎛⎝ah(v,v)+1h∑K∈Th[1h∥v∥2H(K)+h∥∇v∥2H(K)]⎞⎠.

The result follows from observing that . ∎

Finally, we concentrate on a stable decomposition. The convergence theory from [12] requires that the subdomains used for the Schwarz method are at least the same size as the cells in the coarse mesh. Recently, an extension has been published in [10] to include the case of cell-wise subdomains, however, the proof uses nonconforming interpolant and enriching operators for simplices [7].

We achieve a stable decomposition, by a close examination of the proof in [12], which holds for simplices, quadrilaterals, and hexahedra. In particular, it does not require auxiliary spaces with continuity assumptions like Crouzeix-Raviart.

###### Lemma 3.2.

Every admits a decomposition of the form , , which satisfies the bound

 (13) J∑i=0aj(vj,vj)≤CV,Δa(v,v),

with , where and denote the cell diameters used in the fine and coarse meshes respectively.

###### Proof.

Let be the piecewise constant average of on the coarse mesh , let . We decompose in nonoverlapping cell-wise subdomains as follows

 w=J∑j=1vj,

where are uniquely determined. From equation (12) we have

 ah(w,w)= J∑j=1aj(vj,vj)+I(w,w), or equivalently, ah(v−R⊺0v0,v−R⊺0v0)= J∑j=1aj(vj,vj)+I(v−R⊺0v0,v−R⊺0v0),

Reordering and estimating the interface term by its absolute value we obtain

 (14) J∑j=1aj(vj,vj)≤ah(v−R⊺0v0,v−R⊺0v0)+∣∣I(v−R⊺0v0,v−R⊺0v0)∣∣,

using lemma 3.1 we have

 (15) J∑j=1aj(vj,vj)≤c⎛⎝a(v−R⊺0v0,v−R⊺0v0)+1h2∑K∈Th∥v−R⊺0v0∥2H(K)⎞⎠≤c((ah(v,v)1/2+a(R⊺0v0,R⊺0v0)1/2)2+1h2∑D∈TH∥v−R⊺0v0∥2H(D)⎞⎠,

where we used Minkowsky’s inequality and we regrouped the inner products.

We expand the first term and use equation (11) to obtain

 J∑j=1aj(vj,vj)≤ c(ah(v,v)+2ah(v,v)1/2ah(R⊺0v0,R⊺0v0)1/2 +ah(R⊺0v0,R⊺0v0)+H2h2∥v∥2V) ≤ c(2ah(v,v)+2ah(R⊺0v0,R⊺0v0)+H2h2ah(v,v)),

where we used Young’s inequality and coercivity of .

Finally, including the coarse space we achieve

 J∑j=0aj(vj,vj)≤ c(a0(v0,v0)+ah(R⊺0v0,R⊺0v0)+H2h2ah(v,v)).

It remains to bound in a way such that the estimate is independent of the usage of cell-wise or larger subdomains and a constant is achieved as we show below. Since is piecewise constant on , and hence also on ,

 (16) ah(R⊺0v0,R⊺0v0)=δIP∑F∈FIh1h∥R⊺0v+0−R⊺0v−0∥2H(F)+δIP∑F∈FBh1h∥R⊺0v+0∥2H(F),

where we observe

 (17) a0(v0,v0)=hHah(R⊺0v0,R⊺0v0).

Adding and subtracting in equation (16) gives

 ah(R⊺0v0,R⊺0v0)≤cδIP⎛⎜⎝∑F∈FIh1h∥(v−R⊺0v0)+−(v−R⊺0v0)−∥2H(F)+∑F∈FBh1h∥(v−R⊺0v0)+∥2H(F)+∑F∈FIh1h∥v+−v−∥2H(F)+∑F∈FBh1h∥v+∥2H(F)⎞⎟⎠.

The last two terms are obviously bounded by . Also, since is piecewise constant on each , whenever is in the interior of some . Thus,

 ∑F∈FIh1h∥(v−R⊺0v0)+−(v−R⊺0v0)−∥2H(F)+∑F∈FBh1h∥(v−R⊺0v0)+∥2H(F)=∑D∈TH(∑F⊂D∥v+−v−∥H(F)+∑F∈∂D1h∥(v−R⊺0v0)+−(v−R⊺0v0)−∥2H(F)+∑F⊂∂D∈FBh1h∥(v−R⊺0v0)+∥2H(F)⎞⎟⎠≤cah(v,v)+c∑D∈TH1h∥v−R⊺0v0∥2H(∂D).

Now using the trace inequality in equation (10), we obtain

 ∑D∈TH1h∥v−R⊺0v0∥2H(∂D)≤c∑D∈TH1h[1H∥v−R⊺0v0∥2H(D)+H∥v−R⊺0v0∥2V(D)].

Also note that . Hence, applying the approximation result from equation (11) to we obtain

 (18) ah(R⊺0v0,R⊺0v0)≤cHhah(v,v),

therefore, using this result on equation (17) we see that and the result is achieved. ∎

###### Lemma 3.3.

Stable decomposition. The spaces provide a stable decomposition of , with respect to the bilinear form , in the sense of assumption 3.1.

###### Proof.

Let be the stable decomposition constant for the Laplacian, as deduced in lemma 3.2, we then have

 J∑i=0Aj(vj,vj) =J∑i=0{aj(vj,vj)+1ε(Σvj,vi)H} =J∑i=0aj(vj,vj)+1ε(Σv,v)H ≤CV,Δa(v,v)+1ε(Σv,v)H ≤max{CV,Δ,1}Ah(v,v).

It follows that the decomposition for our reaction-diffusion problem is energy stable with where and are the largest and smallest cell diameters respectively. ∎

###### Lemma 3.4.

There exists a strengthened Cauchy-Schwarz inequality in the sense of definition 3.2.

###### Proof.

(See [12, §4.2]). Verifying this inequality consists of obtaining a bound for the spectral radius of the matrix .

That such values exist is a consequence of the Cauchy-Schwarz inequality. The important thing, however, is to obtain a small bound on . To do so, we observe that if the supports of and do not share a face . For the remaining cases, we take . It follows at once from Gershgorin’s circle theorem that

i.e., is bounded by plus the maximum number of adjacent subdomains a given subdomain can have. In practice this number in 2D and in 3D. Even for “unusual” subdomain partitions, this number is not expected to be large. ∎

###### Lemma 3.5 (Local Stability).

There holds

 Ah(R⊺jvj,R⊺jvj)≤ωAj(vj,vj) ∀vj∈Vj,

where for .

###### Proof.

In the case of exact local solvers , in our case the coarse bilinear form uses a penalty parameter depending on the cell diameter of the coarse mesh. Observing the bilinear form (7), we see that for our coarse space bilinear form it holds

 Ah(R⊺0v0,R⊺0v0)≤Hh