1. Introduction
This study is devoted to the development of computationally efficient methods for numerically solving multiphysics and multiscale problems. A representative set of such examples is discussed in the review article [
1]. Such applications are also considered in [
2,
3,
4]. More specifically, our attention is focused on a class of problems that involve coupling at interfaces, which are manifolds of lower dimensions. Such mathematical models are often described by saddle-point systems of differential equations. For example, this is the case when either a mixed formulation is used to ensure better accuracy and robustness or when Lagrangian multipliers are introduced to interpret the coupling at the interface. Thus, interface conditions in fractional Sobolev spaces naturally appear. For more related details, we can refer to [
5,
6,
7].
Aimed at presumably more realistic applications, we are interested in the numerical solutions of multidimensional boundary value problems in bounded domain with general geometry . We assume that the geometry of the interface is general, as well. After discretization by a certain mesh method (such as finite differences, finite elements, finite volumes, etc.) on a properly generated unstructured mesh (with possible mesh refinements), we obtain a system of linear algebraic equations with a sparse unstructured matrix . Under the above assumptions, the discrete problems are apparently very large scale. For such problems, the advantages of iterative solution methods (solvers) are well expressed. Note that the condition number of the corresponding matrix can be very large. It could be due not only to the number of degrees of freedom but also to the structure of the mesh. For the multiscale and multiphysical problems under consideration, the media can be very highly heterogeneous and the anisotropy ratio can be very large. Thus, looking for solvers of optimal computational complexity gives rise to the requirement for efficient and robust preconditioning.
The block-diagonal preconditioning is among the most efficient approaches for numerically solving the studied saddle-point problems. Let us consider the system of differential equations
written in the block-operator form
where the first
rows correspond to the volume equations while the last one represents the interface condition.
Definition 1. The block-diagonal preconditioners have zero off-diagonal blocks and nonsingular diagonal blocks. They are also known as additive preconditioners.
In this study, we analyze two types of block preconditioners that have optimal convergence rate, that is, the number of iterations is
. Following the notation from Equation (
1), they can be written as follows:
and
where
, respectively
. The first
blocks of
and
are selfadjoint positive definite differential operators,
is the identity operator at the interface
, and
stands for the Laplacian at
.
At discrete level, the implementation of the first blocks of requires solving systems with sparse symmetric and positive definite (SPD) matrices. The same holds true for the first blocks of . For all of them, a standard PCG (e.g., AMG) solver of optimal computational complexity can be applied. The matrices, corresponding to the last blocks, are dense.
Thus, for the application of the preconditioner , the solution of a system with a matrix in the form is required, while the implementation of leads to matrix vector multiplication with . Here, stands for the matrix corresponding to the discretization of . In both cases, numerical computations with a fractional power in the unit interval of an SPD matrix are involved—either via solving a linear system of equations or via a matrix-vector multiplication. As we will see later on, although the implementation of the last block of is inverse free, it is not less challenging (but even more complicated) than the case of .
In this context, multigrid methods for discrete fractional Sobolev spaces are developed in [
5]. The abstract additive multilevel framework is adapted there. The constructed smoother involves small
blocks in the form
, where
is the graph-degree of the mesh node
k.
As an alternative, we propose to use the best uniform rational approximation (BURA) method for
[
8,
9,
10,
11]. The aim is to develop a robust preconditioning algorithm of optimal computational complexity. As correctly noted in [
5], a disadvantage of the method introduced in 2018 in [
11] was that the accuracy depended on the condition number of the matrix
. This drawback has been overcome in the improved BURA approximation developed in [
10]; see also the survey paper [
8]. It reads as
, where
is the best uniform rational approximation of degree
k of
on
and
is the smallest eigenvalue of
. The
error of the BURA approximation decreases with
k, satisfying the asymptotic estimate
. Further improvement of the computational efficiency of the method is presented in [
9], where problem-specific model reduction techniques are utilized.
The structure of the paper is as follows. Two examples of coupled multiphysics problems, leading to the class of discussed block-diagonal preconditioners, are presented in
Section 2.
Section 3 is devoted to the construction and the basic properties of the BURA method. Then, the accuracy of the BURA-based preconditioners is analyzed in
Section 4, including the derived estimates for the corresponding condition number. The theoretical results are supported by representative numerical experiments. The implementation issues are discussed in
Section 5, ending with the proof of optimal computational complexity of the studied methods and algorithms. The motivation for the choice of numerical tests presented in
Section 6 is twofold. On one hand, they illustrate the accuracy of the derived theoretical estimates, while on the other, they provide an additional practical insight for the efficient usage of BURA preconditioners with a small degree
k. Finally, we give some conclusions in
Section 7.
2. Monolithically Coupled Multiphysics Problems: Two Examples
In this section, we present two examples of related multiphysics problems that lead to block-diagonal preconditioners in the form (
2) and (
3). These preconditioners are already known and they are important for the analysis, presented in this paper as real world examples, where fractional elliptic problems are formulated. In this paper, we do not analyze the specific properties of the approximations of multiphysics equations or the influence of the different blocks of the preconditioner on the total condition number of the preconditioned system. The main goal of our study is to obtain detailed and asymptotically exact estimates of the accuracy of the BURA method. To simplify notation, throughout the paper we will use
N for both the size of an SPD matrix and the number of unknowns of a coupled discrete problem. From the context, it will be clear to which of the two possibilities
N refers to. Similarly, we denote by
various SPD matrices and again the context will clarify the particular object of interest.
2.1. Example 1: Coupling of Darcy and Stokes Flows
Let
be a closed domain in
, which is a disjoint union of two subdomains
and
, that is,
.
Figure 1 shows two examples of model problems with simple geometry and more realistic problem with a more general geometry of the interface.
The exchange through the interface
between a Stokes flow in
and a Darcy flow in
is modeled by the system of equations in the form [
7,
12,
13]:
Here, , are the velocities and pressure for the Stokes problem in , , are the velocities and pressure for the Darcy problem in , K and are the hydraulic conductivity, respectively, the fluid viscosity, while the coefficient D corresponds to the Beavers–Joseph–Saffman condition. Finally, and are unit vectors that are normal and tangential to the interface . We further assume that the problem is equipped with correct boundary conditions.
Similarly to (
1), the system of differential Equation (
4) can be written in the block-operator form
, where
The block sparsity structure is shown here. Then, after discretization by some mesh method (such as finite elements, finite differences, or finite volumes) we will obtain a linear system, where all blocks that correspond to nonzero operators in (
5) will be also sparse.
There are two groups of methods for preconditioned iterative solution of coupled multiphysics problems: monolithic and based on some splitting [
5,
14,
15]. For example, splitting techniques include operational splitting (by physical processes), domain decomposition (geometrical splitting), and block partitioning (algebraic splitting). The methods considered in this paper belong to the second group, following the motivation from [
13]. Additionally, our study is supported by some recent achievements in the solution methods for fractional diffusion problems, or more generally for linear algebraic systems with fractional power of sparse SPD matrices.
The case
is considered in [
13]. It should be noted that this simplified problem is representative enough to illustrate why the norms studied in [
12] are not sufficient for the proposed preconditioning approach. In [
16], it is shown that the block-diagonal preconditioner
provides an optimal condition number estimate. Here,
is a tangential trace operator, while
stands for the identity operator. The first four blocks of (
6) are standard components for preconditioning Darcy and Stokes problems. Thus, after discretization, the matrices corresponding to the first two blocks are sparse SPD, to which some Preconditioned Conjugate Gradients (PCG) solver (e.g., AMG) of optimal computational complexity can be applied. Note that the local operator
only increases the diagonal dominance in the first block. The last block leads to a system with the dense matrix
.
This example raises the question of the efficient solution/preconditioning of systems with matrices in the form , is sparse SPD, . Until recently, this was a challenge.
2.2. Example 2: Coupling a 1D Structure Embedded in
Let
be a bounded domain, while
represents a 1D manifold (structure) inside
, see
Figure 2. The following trace coupled problem is considered [
17]:
Here u, v, p are unknowns, the 1D manifold is parameterized in terms of s, the Laplacian on is understood in terms of the Laplace–Beltrami operator, is a suitable trace operator, and is a Dirac measure such that for a continuous function w.
At the same time, it is worth noting that the equations used in the real-life simulations are more complicated than (
7); thus, the above example can be considered as a relevant model problem.
As it can be seen in [
17], the standard finite element methods provide accurate error estimates in the norms of the introduced fractional order Sobolev spaces. The saddle-point matrix corresponding to the obtained system of linear algebraic equations is sparse. Arguments, similar to those discussed after (
4) in the case of the Darcy–Stokes problem, are applicable here. However, standard techniques are not directly applicable to effective preconditioning. In what follows, we will consider the block-diagonal preconditioner
introduced in [
6]. The choice of the fractional power in the third block of (
8) is analyzed there, observing that values
yield bounded condition numbers. The presented additional stability analysis has led to the preferred value
. Similar arguments to those concerning the implementation of the preconditioner
are applicable to
.
The efficient approximate matrix vector multiplication with the matrix , where is a sparse SPD matrix and , is the challenging question raised by this example.
2.3. Implementation of and : The Key Question
The key questions in implementing the block-diagonal preconditioners
and
are raised by the last blocks. They concern the fractional powers of SPD matrices defined through the spectral decomposition
where the columns of
are the orthonormal eigenvectors of
and
stands for the diagonal matrix of the related eigenvalues. Under these assumptions,
, the matrix
is unitary,
, which leads directly to (
9). Therefore
The Formula (
10) could be used in practical computations if the
N eigenvectors and eigenvalues are explicitly known and if the matrix vector multiplication with
is equivalent to a fast Fourier transform. In such cases, the computational complexity is almost linear,
. However, when talking about fractional Laplacian on
, this limits the applications to interfaces with simple geometry (like the left sides of
Figure 1 and
Figure 2), uniform meshes, and to the lowest order finite element approximations.
The multigrid methods proposed and analyzed in [
5] are aimed at handling general geometries of the interface
. As noted in the Introduction, the smoother includes non-local (dense) but small
blocks
, where
is the graph-degree of the mesh node
k. Under certain regularity assumptions and fixed number of the multigrid levels
J, the condition number
is uniformly bounded for each given
. More precisely,
or
, depending on the regularity assumptions.
The case of negative power
corresponds to the preconditioner (
3), and more specifically to (
8). The implementation of this kind of preconditioner requires a matrix vector multiplication with
. For this multiplication, one can use the representation
thus reducing the problem to numerically solving linear system with
,
. Following (
11), the second variant of the multigrid method in [
5] is adapted to handle negative powers. The condition number of the proposed preconditioner depends on
N, including the case of a fixed number of multigrid levels
J.
In the following sections, we consider the application of the BURA method in the general framework of the preconditioners
and
(Equations (
2) and (
3)) and, respectively,
and
(Equations (
6) and (
8)). Our study is aimed at developing uniformly bounded condition number estimates, which are robust in terms of the geometry of the
interface, as well as in terms of regularity of the solution. The ultimate goal is to obtain preconditioned solution methods of optimal computational complexity.
3. The BURA Method
The fractional diffusion problems are non-local and their numerical solution is costly. The survey paper [
8] tracks some of the last years’ achievements on this topic. They include a variety of methods based on: (i) extension of the problem from
to a local elliptic problem in
, (ii) reformulation as a pseudo-parabolic problem in the cylinder
, (iii) methods based on approximation of the Dunford–Taylor integral representation of the solution, and (iv) methods based on the best uniform rational approximation (BURA) of
on
. Applying rational approximations is one of the most successful approaches to efficiently solve non-local problems. For example, following the reformulation to a pseudo-parabolic problem, an efficient method based on diagonal Padé approximation to
is developed in [
18]. Recently, a unified view to the methods based on the approaches (i–iv) was presented in [
19], showing that all of them can be interpreted as generated by some rational approximation of
on
. In this sense, BURA methods are asymptotically the best.
The BURA method is originally introduced in [
11] where the solution of the linear algebraic system
is written in the form
with
, with
s a (presumably small) positive integer. The proposed approximate solution method is based on the best uniform rational approximation (BURA) of the function
for
. Quite promising numerical results for
, accompanied by representative comparison analysis with other known methods, were presented in [
11]. However, the convergence rate of the first BURA method was not robust with respect to the condition number of
, which is correctly noted in [
5]. This shortcoming is completely overcome in [
10]. Several aspects of the step-by-step development of the BURA method can be traced in the survey paper [
8]. Note that in [
19], it was shown that the available numerical solvers of (
12) are based on uniform rational approximation of
,
. Therefore, methods based on best uniform rational approximation are expected to be the most efficient. From now on, when we talk about BURA, we will have in mind the following variant of the method proposed and analyzed in [
10].
Let
be the best uniform rational approximation of
in
in the class of rational functions
,
and
are polynomials of degree
k. The case of rational approximations with equal degrees of the polynomials
P and
Q is known as the diagonal sequence of the Walsh table. The minimizer
satisfies the following sharp exponential error estimate with respect to the degree
k [
20]:
Then, the BURA approximation of
is defined as
where
is the smallest eigenvalue of the SPD matrix
.
The Padé approximation is also defined in the class of rational functions. However, while the minimizer is implicitly defined, the Padé approximation is a Hermite interpolation. The advantage of the Padé approximation is that it is easier to compute.
Let be obtained after finite difference discretization of a given elliptic operator. In the context of this paper, we can assume that corresponds to the Laplacian with homogeneous Dirichlet boundary conditions. Then, , which follows, e.g., from the Poincaré–Steklov inequality. The same is true if the finite element method (FEM) is applied to discretize the Laplacian. In this case, where and are the stiffens and mass matrices, respectively. Then, the SPD property of the FEM matrix is understood in terms of the inner product induced by the mass matrix .
Let us denote by
the BURA approximation of the solution
of the fractional diffusion linear system (
12). Then, applying spectral basis representation of
,
and of
, and the Parseval’s formula, we obtain the error estimate (see, for more details [
10])
The estimate (
14) clearly shows the exponential rate with respect to the degree
k of the rational approximation. Moreover, its dependence on the fractional power is explicit. Smaller
means lower accuracy. Since
is uniformly bounded from below, the estimate is robust with respect to the condition number
. The estimate (
14) is very sharp. The numerical tests conducted in [
9] illustrate how close the computed errors are to the theoretical estimate, when the latter is greater than the double-precision accuracy.
The properties of the single roots and the related interlacing poles of
are well studied [
21]. Both the additive and the multiplicative representation of the rational function
can be used in the implementation of the BURA method. The corresponding matrix valued expressions read as
and
where
and
. A numerically stable approach for the computation of these coefficients is developed in [
22]. Both the additive (
15) and the multiplicative (
16) algorithms reduce the problem to solving
k systems with sparse SPD matrices
, which are positive diagonal perturbations of
. This leads to the understanding of the degree
k as a universal measure for algorithmic efficiency, thus giving rise to the following asymptotic estimate of computational efficiency of the BURA methods
The estimate (
17) holds true if a method with optimal computational complexity is used to solve the auxiliary linear systems that appear in (
15) or (
16). In the multidimensional case, this means that such a preconditioned iterative solver is applied. In the numerical tests presented in our previous papers [
9,
11,
23], the Boomer-AMG implementation from HYPRE [
14] of the algebraic multigrid (AMG) preconditioner in the PCG framework is used for this purpose.
The computational complexity depends on the target accuracy of the numerical solution. Assume that the matrix
corresponds to FEM discretization of the fractional diffusion problem. The rigorous error analysis is a task beyond the standard theory of FEM. Such error bounds are derived in [
24] under certain assumptions that are an interplay among data regularity, regularity pickup of the differential operator (e.g., the Laplacian), and the fractional power
. If linear finite elements are used and if
, the
error is essentially
The case where the finite difference method (FDM) is applied for discretization of the fractional diffusion problem is analyzed in [
10]. The presented error estimates are based on certain equivalence between the lumped mass FEM approximation and the
-stencil FDM approximation.
The further analysis integrates the FEM/FDM error estimate (
18) in the complexity estimate (
17). Assuming the mesh is quasi uniform,
, with
, and therefore
Then, for a given mesh parameter
h, the degree of BURA approximation
k is chosen such that the contributions of the discretization error (
18) and the BURA error (
14) to the total error are balanced. In this way, and taking into account (
17), the following almost optimal computational complexity
is derived in [
10].
More recently, the so-called reduced additive BURA-AR and multiplicative BURA-MR methods were introduced in [
9]. This study follows the observation that the matrices of part of the auxiliary systems possess very different properties. An important observation is that under certain assumptions, the computational complexity of the reduced BURA can be improved to
.
7. Concluding Remarks
The fractional diffusion operators are non-local, and the related matrices are dense. Furthermore, when using the spectral definition, the matrix is generally non-computable. In such a situation, we are not even able to compute the residual, which should directly rule out the idea of applying any iterative solution methods to systems in the form .
Why, then, do we discuss preconditioning methods for
? For the considered coupled problems such preconditioning makes sense, as the corresponding systems have sparse matrices (so computing of the residue is not a problem), while the fractional power matrices appear only in the block-diagonal preconditioners (
6) and (
8).
Our study is motivated by the results of the recently published papers [
5,
13,
17], where the developed multigrid methods for discrete fractional Sobolev spaces are applied to the Darcy–Stokes and 3D–1D coupled problems. Alternatively, we have proposed a fully algebraic preconditioning approach.
For
, the new preconditioner
relies directly on the best uniform rational approximation (BURA) of degree
k of
. An important further contribution is the construction of
for
leading to the equality
This allowed us to develop a unified approach to the relative condition number estimates (
21) and (
24), applied later in Theorem 1 and Theorem 2, where the optimal computational complexity for both considered coupled problems is proved. The numerical study of
and
is also important for the value of the paper. As shown, a very small degree of
k is practically sufficient to ensure high quality BURA preconditioning.
In general, it is clear that a required accuracy can be obtained by applying different known methods for spectral fractional diffusion problems. Here, we apply the BURA method, motivated by its advantages, as proven in [
19] (see the discussion at the beginning of
Section 3).
For the simplest case of
, the comparison with the multigrid preconditioners from [
5] illustrates the competitiveness of the newly proposed
and
. For most of the BURA preconditioners, even stronger numerical advantages can be expected when more complex coupled problems are considered, due to the more general scope of applicability of the presented results. Among the other advantages, we note the following: BURA preconditioners are purely algebraic, so regularity assumptions are not made in the theoretical analysis. The theoretical estimates (
21) and (
24), as well as those shown in
Table 1, are exactly the same for every type of geometry of
. The BURA approach and the obtained condition number estimates are applicable to a wider class of problems, including, for example, coupling of a 2D manifold embedded in
.