Log In Sign Up

A parallel-in-time multigrid solver with a new two-level convergence for two-dimensional unsteady fractional Laplacian problems

The multigrid-reduction-in-time (MGRIT) technique has proven to be successful in achieving higher run-time speedup by exploiting parallelism in time. The goal of this article is to develop and analyze an MGRIT algorithm, using FCF-relaxation with time-dependent time-grid propagators, to seek the finite element approximations of two-dimensional unsteady fractional Laplacian problems. Motivated by [B. S. Southworth, SIAM J. Matrix Anal. Appl. 40 (2019), pp. 564-608], we provide a new temporal eigenvalue approximation property and then deduce a generalized two-level convergence theory which removes the previous unitary diagonalization assumption on the fine and coarse time-grid propagators. Numerical computations are included to confirm theoretical predictions and demonstrate the sharpness of the derived convergence upper bound.


page 1

page 2

page 3

page 4


Finite element approximation of fractional Neumann problems

In this paper we consider approximations of Neumann problems for the int...

Two-level error estimation for the integral fractional Laplacian

For the singular integral definition of the fractional Laplacian, we con...

On the approximation of dispersive electromagnetic eigenvalue problems in 2D

We consider time-harmonic electromagnetic wave equations in composites o...

Asynchronous Truncated Multigrid-reduction-in-time (AT-MGRIT)

In this paper, we present the new "asynchronous truncated multigrid-redu...

Fractional Laplacian - Quadrature rules for singular double integrals in 3D

In this article, quadrature rules for the efficient computation of the s...

On the usefulness of lattice approximations for fractional Gaussian fields

Fractional Gaussian fields provide a rich class of spatial models and ha...

1 Introduction

Time-dependent problems involving fractional Laplacian operator are of particular interest in fractional calculus p-01 ; w-02 ; m-02 ; c-08 . The fractional Laplacian, which is the infinitesimal generator of a standard isotropic -stable Lévy flight s-01 , has been used in a broad range of applications including quantum mechanics l-03 , finance c-06 , elasticity d-01 , porous medium p-02 and stochastic dynamics g-04 . Further details on this topic and equivalent characterizations can be found in the review article l-01 . We are concerned in this paper with a non-intrusive parallel time integration with multigrid to a class of two-dimensional unsteady fractional Laplacian problems of the form


where , the fractional order , solution domain is an open and bounded subset with Lipschitz continuous boundary .

One of the main challenges of problem (1.1) lies in the nonlocality of the fractional Laplacian. Numerous numerical approximations have been proposed to treat the fractional Laplacian on bounded domains in recent years, such as finite difference h-01 ; s-02 ; w-08 ; d-03 ; h-04 ; f-03 ; w-05 ; w-07 ; l-06 ; w-09 , finite element (FE) k-01 ; w-10 ; n-01 ; g-05 ; c-01 ; a-02 ; a-05 ; j-01 , finite volume y-02 ; c-07 ; g-08 , discontinuous Galerkin x-02 ; l-05 ; z-05 ; a-04 and spectral (element) b-01 ; m-04 ; d-04 ; z-01 ; h-02 ; w-06 ; z-02 ; w-11 ; t-02 ; x-03 methods. Utilizing a similar technique to the fractional Laplacian in c-05 , we can obtain the Caffarelli-Silvestre extension, which is reformulated the original problem posed over on the semi-infinite cylinder with dynamic boundary condition:


where the parameter , the lateral boundary , the positive normalization constant and the conormal exterior are respectively defined by

A widely used approach to solve nonlocal problem (1.1) is based on the solution of local problem (1.2), however, incorporating one more dimension in space n-01 ; c-01 . Specifically, exploiting the intermediate function , we immediately get the desired solution subject to the expression for all and . It is worthwhile to point out that the anisotropic spatial mesh must be used to capture the singular behavior of as d-05 , and Nochetto et al. have proved the existence and uniqueness for solutions of problems (1.1) and (1.2) (see (n-02, , Theorem 6)) in a purely theoretical manner.

Moving from two-dimensional problems to three-dimensional problems can be difficult for solvers. Chen et al. developed uniformly convergent multilevel methods to tackle steady fractional Laplacian problems c-01 , which is an important constructive study in the literature. However, a more severe challenge to be faced with arises from the facts that (i) the higher the spatial resolution signifies a comparative or greater scaling in the temporal direction, (ii) future computing speed must rely on the increased concurrency provided by not faster but more processors. The practical consequence is that solution algorithms limited to spatial parallelism for problems with evolutionary behavior shall entail long overall computation time. Therefore, solution algorithms achieving parallelism in time have been of especially high demand. Although time is sequential in essence, parallel-in-time methods are not new and may be dated back to n-03 ; m-05 . A recent review of the literature can be found in g-07 . At present, parareal in time l-02 and multigrid-reduction-in-time (MGRIT) f-02 are two of the most popular choices. Wu and Zhou analyzed convergence properties of the parareal algorithm for space fractional diffusion problems w-01 ; w-04 and time fractional equations with local time-integrators w-03 in constant time-steps. Note that Gander and Vandewalle interpreted parareal as a two-level multigrid full approximation scheme g-02 , whereas Falgout et al. identified parareal with an optimal two-level MGRIT algorithm and further pointed out that the obvious generalization of parareal to multiple levels couldn’t produce an optimal algorithm f-02 . In addition, the concurrency of parareal is still limited as its relatively large sequential coarse-grid solve. Conversely, MGRIT, a truly multilevel algorithm implemented in the open-source package XBraid x-01 , has been proven to be rather effective and analyzed in two-level d-02 and multilevel h-03 settings, for integer order basic parabolic and hyperbolic problems using L-stable Runge-Kutta schemes, with the limitation on analysis that they only consider the same fine time-grid propagators. Extensions on the two-level convergence analysis to the case of nonuniform temporal meshes are provided in y-01 ; s-03 . However, we should point out that only MGRIT using F-relaxation was abstractly discussed in s-03 and the unitary diagonalization assumption must be enforced in y-01 . Moreover, from the survey of references, there are no calculations taking into account of MGRIT to the fully implicit FE discretization of problem (1.1).

The arrangement of this work proceeds as follows. Section 2 introduces the functional framework and outlines the fully implicit discrete FE scheme. In Section 3, we present the MGRIT V(1,0)-cycle algorithm using FCF-relaxation followed by a generalized two-level convergence analysis without the assumption of unitary diagonalization. Section 4 conducts numerical experiments to illustrate the optimal order of convergence and the computational error-norms decay rate of the fully discrete FE scheme, the robustness of Chen’s multigrid with line smoother with respect to the fractional order and the spatial mesh size, as well as the sharpness of the derived convergence upper bound on the two-level MGRIT algorithm. Finally, in Section 5, we draw some conclusions and discuss future work.

2 Mathematical preliminaries

This section is devoted to introducing some mathematical preliminaries, including the variational formulation of the local problem (1.2) and its fully implicit discrete FE scheme.

2.1 The variational formulation

Let be open, and denote the weighted Lebesgue and Sobolev spaces g-06

with norms

where is the distributional gradient of the measurable function .

Since , weight belongs to the Muckenhoupt class m-03 , which implies that is Hilbert. In order to settle problem (1.2), we define

Thus, given a forcing term and an initial datum , the function solves problem (1.1) if and only if its harmonic extension solves problem (1.2), where


Recall that the trace operator onto satisfies (see (n-01, , Proposition 2.5))

we get the variational formulation of problem (1.2): to find subject to and


where is the duality pairing between and .

2.2 The fully implicit discrete FE scheme

Note that problem (1.2) can’t be approximated except that is truncated to for a suitable with and an additional homogeneous Dirichlet on , whose theoretical foundation is the exponential decay of in the extended direction on the basis of the length , see (n-02, , Proposition 22).

The first step toward the discretization is to obtain a semi-discrete approximation of the truncated problem of the original problem (2.1). To do so, we associate a (possibly nonuniform) mesh by points and time step sizes for . The backward Euler discretization is applied: to seek subject to and

for all , where . Obviously, the sequence

is a series of piecewise constant approximations to the solution of problem (1.1).

Next, we turn to discuss the spatial discretization. Given a quasi-uniform triangulation over and a graded partition of by the transition points


as well as segments for , we construct the mesh over

as the tensor product of

and , and define

and the weighted elliptic projector such that

is valid for any function , where denotes the space of all polynomials of at most one degree.

Finally, we arrive at the fully discrete FE scheme: with the initial guess , to calculate satisfying the constraint


As before, approximate solutions of problem (1.1) are given by

In order to get the matrix formulation of (2.3), it is necessary to introduce


where is the number of interior vertices in . It is straightforward to verify that and can be respectively written as

Choosing the test function in (2.3), we deduce

by Green’s first identity and the fact that , and one further expansion via separation of variables

where entries

Thus, the system of linear equations obtained from (2.3) is of the form



  • matrixes , , and ,

  • represents a zero matrix of suitable size,

    denotes the Kronecker product,

  • with for to approximate the initial ,

  • is the auxiliary vector produced by the other points in the extended direction,

  • with .

It is trivial to see that the coefficient matrix of (2.4) is symmetric positive definite (SPD). Therefore, the multigrid with line smoother proposed in c-01 is a very suitable method and is employed in our implementation.

3 MGRIT formulation and a generalized two-level convergence upper bound

The main focus of this section is the general MGRIT algorithm using FCF-relaxation, followed by its two-level convergence analysis for the general case where the time-grid propagators don’t have to be the same and need not be simultaneously diagonalizable with a unitary matrix.

3.1 The general MGRIT algorithm using FCF-relaxation

Consider the following block unit lower bidiagonal system


where ,

is the identity matrix, vectors

, and time-grid propagators are assumed not to be connected with any specific forms.

Figure 1: Schematic view of the CF-splitting by a factor of four.

Various components are required in order to solve the global space-time problem (3.1) through the MGRIT algorithm based on the multigrid reduction technique r-01 . Define the coarse temporal mesh with nodes and time step sizes for some coarsening factor , where . Then all are disconnected C-points representing those variables in the coarse level and the others are F-points grouped adjacent points into a bunch. The particular case when is illustrated in Fig. 1. It should be mentioned that is still viewed to be nonuniform. Assume for simplicity of presentation here that F-points are ordered first. MGRIT utilizes the coarse injection as the restriction operator, i.e., only the residual at C-points is computed, where is the identity matrix. The FCF-relaxation (i.e., an initial F-relaxation followed by a C-relaxation and then a second F-relaxation) is often the most preferable choice to produce optimal and scalable multilevel iterations f-02 , where F-relaxation refers to updating variables at F-points in turn yet independent of other F-intervals, and C-relaxation stands for performing (2.4

) to propagate the solution at a C-point. It is apparent that these two processes are both highly parallel. The ideal interpolation is used to provide the coarse-grid correction, which is implemented at the cost of updating C-variables by the coarse approximations immediately followed by an F-relaxation over the time-grid propagators

. It is worth mentioning that the resulting Petrov-Galerkin triple matrix product is exactly the Schur complement


where with , positive integer and nonnegative integer . However, this coarse-grid problem would be as expensive as the original fine-grid problem. MGRIT approximates the coarse-grid operator with


where is an approximate coarse-grid time-stepping operator to . The most natural choice of is to re-discretize the time-dependent linear continuous problem on the coarse temporal mesh. The final stage is the application of the above processes recursively, with possibly different coarsening factors at different time levels, to construct a multilevel hierarchy under certain stopping criteria, e.g., by setting the maximum number of time levels and the minimum coarsest temporal grid size. Below is a sketch of the resulting MGRIT V(1,0)-cycle algorithm using FCF-relaxation in the recursive fashion, where is the number of time levels determined by the user supplied stopping procedure, , , matrices , and () respectively correspond to the -th coarse-scale time re-discretization, the mentioned above restriction and interpolation.

Algorithm 1

MGRIT V(1,0)-cycle algorithm using FCF-relaxation:

Step 1

Apply FCF-relaxation to with as the initial guess.

Step 2

Inject the residual to the coarse time level .

Step 3

If , then solve ; Else compute .

Step 4

Perform the coarse-grid correction using the ideal interpolation .

Remark 1

(Reduction of the computational overhead) Here, we point out that Step 4 of Algorithm 1 is executed by just correcting C-variables, and updating F-variables only when the Euclidean norm of the residual is small enough. The underlying reason is that corrections at F-points are equivalent to an F-relaxation, which will be performed in the subsequent MGRIT iteration.

3.2 A generalized two-level convergence analysis

Setting , Algorithm 1 in this case reduces to the well-known two-level scheme, whose residual propagation operator for a single iteration, as in the derivations of y-01 ; s-03 , can be expressed as in an FC-partitioning of :


where the matrix used in forming the “ideal restriction” operator

Below, we wish to find an upper bound for a certain norm of the residual propagation .


Assume that time-stepping operators and are both strongly stable l-04 (i.e., the discrete -norms and ) and simultaneously diagonalizable by the matrix . Note that the assumption on simultaneous diagonalization will be tenable when the same spatial discretization is exploited on fine and coarse temporal meshes. Denote by and the sets of all complex eigenvalues of and , respectively.


Through simple algebraic calculations, we arrive at an estimate of the second term in (


From the relation for in Stage S1, we have

and hence .


By exploiting the closed form

the other two matrix forms appeared in are obtained from block lower bidiagonal structures (3.2)-(3.3):


As a result, the expression of can be rewritten as


which can be reformulated as the under-mentioned triple matrix product

Then, under the assumption that is invertible for , using the relation , and recalling that the Moore-Penrose conditions are used to define the pseudoinverse of any matrix :

the unique pseudoinverse of can be established in a straightforward way


by noticing the following two important results

It is easy to see that the use of

allows for a simpler technical path to analyze the maximum singular value of

through the minimum nonzero singular value of with the advantage of its simplicity.


Based on the simultaneous diagonalizability from Stage S1:

we have


  • consisting of diagonal blocks,

  • and are obtained from replacing , by , in (3.5) and (3.6), respectively,

  • denotes the orthogonal permutation matrix which makes block diagonal,

  • and represent and , respectively, evaluated at () and ().

We need to determine the nonzero singular value . For this purpose, deleting the last two zero-rows and zero-columns of which correspond to two zero eigenvalues, the remaining algebraic manipulation is to pursue after the square root of the smallest nonzero eigenvalue of the complex Hermitian tridiagonal matrix