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 .
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 ). Convergence estimates for such methods applied to diffusion problems have been developed in . 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 . However, its application does not extend to quadrilaterals and hexahedra, since the proof uses
nonconforming interpolant and enriching operators for simplices. We provide an extension for quadrilaterals and hexahedra without such restrictions, which is based on the original proof in .
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
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
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
where is the standard Sobolev space with zero boundary trace. They are provided with inner products
The weak form of problem (1) is: find such that
where and the bilinear form is
The second line uses the vector notation
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.
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.
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 . 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.
be the space of tensor product polynomials of degree up toin 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
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
We introduce the following definition of mesh integrals
and integrals over and are defined accordingly. The interior penalty (IP) bilinear form for the scalar Laplacian, as described in , reads
where is the minimum cell diameter adjacent to the face, and . We have replaced the jump operator used in  for the equivalent expression: . Coercivity and continuity are proven in  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
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
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
where and are the coefficient vector of the representation of and respectively in terms of the chosen basis.
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
denote prolongation operators for which there holds
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
with such that
where and accordingly.
If , admits a stable decomposition without including the coarse space as follows (see [20, p.49])
Assumption 3.2 (Strengthened Cauchy-Schwarz inequality).
There exist constants for such that
We will denote the spectral radius of by .
Assumption 3.3 (Local stability).
There exists such that
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
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
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
Finally, our additive Schwarz preconditioned system reads
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
where denotes the identity operator on . Using we define the multiplicative Schwarz preconditioner as
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:
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.
Where and are the smallest and largest eigenvalues of the preconditioned system, respectively.
See [20, §2.3]. ∎
where is a constant independent of and .
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.
3.2. Application to the discrete problem
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  that we need for our proof.
For any , there holds the trace inequality (see [12, Lemma 3.1])
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])
In particular, if
Let , let , , be given (uniquely) by , . Then the following identity holds (see [12, Lemma 3.3])
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
There exists a constant such that
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
The result follows from observing that . ∎
Finally, we concentrate on a stable decomposition. The convergence theory from  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  to include the case of cell-wise subdomains, however, the proof uses nonconforming interpolant and enriching operators for simplices .
We achieve a stable decomposition, by a close examination of the proof in , which holds for simplices, quadrilaterals, and hexahedra. In particular, it does not require auxiliary spaces with continuity assumptions like Crouzeix-Raviart.
Every admits a decomposition of the form , , which satisfies the bound
with , where and denote the cell diameters used in the fine and coarse meshes respectively.
Let be the piecewise constant average of on the coarse mesh , let . We decompose in nonoverlapping cell-wise subdomains as follows
where are uniquely determined. From equation (12) we have
Reordering and estimating the interface term by its absolute value we obtain
using lemma 3.1 we have
where we used Minkowsky’s inequality and we regrouped the inner products.
We expand the first term and use equation (11) to obtain
where we used Young’s inequality and coercivity of .
Finally, including the coarse space we achieve
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 ,
where we observe
Adding and subtracting in equation (16) gives
The last two terms are obviously bounded by . Also, since is piecewise constant on each , whenever is in the interior of some . Thus,
Now using the trace inequality in equation (10), we obtain
Also note that . Hence, applying the approximation result from equation (11) to we obtain
therefore, using this result on equation (17) we see that and the result is achieved. ∎
Stable decomposition. The spaces provide a stable decomposition of , with respect to the bilinear form , in the sense of assumption 3.1.
Let be the stable decomposition constant for the Laplacian, as deduced in lemma 3.2, we then have
It follows that the decomposition for our reaction-diffusion problem is energy stable with where and are the largest and smallest cell diameters respectively. ∎
There exists a strengthened Cauchy-Schwarz inequality in the sense of definition 3.2.
(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).
where for .
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