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
where is the distributional gradient of the measurable function .
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 ofand , 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
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,
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
is the identity matrix, vectors, and time-grid propagators are assumed not to be connected with any specific forms.
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.
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 .
(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.
By exploiting the closed form
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 ofthrough the minimum nonzero singular value of with the advantage of its simplicity.
Based on the simultaneous diagonalizability from Stage S1:
consisting of diagonal blocks,
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