Next Article in Journal
The Smallest Singular Value Anomaly and the Condition Number Anomaly
Previous Article in Journal
The Local Antimagic Total Chromatic Number of Some Wheel-Related Graphs
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

A Comparison of Parallel Algorithms for Numerical Solution of Parabolic Problems with Fractional Power Elliptic Operators

1
Department of Mathematical Modelling, Faculty of Fundamental Sciences, Vilnius Gediminas Technical University, Saulėtekio al. 11, LT-10223 Vilnius, Lithuania
2
Kaunas Faculty, Institute of Social Sciences and Applied Informatics, Vilnius University, Muitinės St 8, LT-44280 Kaunas, Lithuania
*
Author to whom correspondence should be addressed.
Axioms 2022, 11(3), 98; https://doi.org/10.3390/axioms11030098
Submission received: 28 January 2022 / Revised: 20 February 2022 / Accepted: 22 February 2022 / Published: 25 February 2022

Abstract

:
In this article we construct parallel solvers analyze the efficiency and accuracy of general parallel solvers for three dimensional parabolic problems with the fractional power of elliptic operators. The proposed discrete method are targeted for general non-constant elliptic operators, the second motivation for the usage of such schemes arises when non-uniform space meshes are essential. Parallel solvers are required to solve the obtained large size systems of linear equations. The detailed scalability analysis is done in order to compare the efficiency of prposed parallel algorithms. Results of computational experiments are presented and analyzed.

1. Introduction

Currently, fractional differential equations are used to simulate non-standard diffusive processes when the diffusion processes can’t be desribed by classical mathematical models. We mention applications of chemical and contaminant transport in heterogeneous aquifer [1], recent models in physics [2], biomedicine [3], and optics and image applications [4]. More examples are given in [5,6].
It is well-known that the fractional power of an elliptic operator A h α can be defined in a non-unique way. We restrict to the spectral definition, and detail definitions are presented in Section 2. It is interesting to stress, that the classical spectral algorithm and some modifications can be used as computational tools to solve parabolic type problems with nonlocal operator A h α , making the complexity of the nonlocal algorithms the same as for spectral methods targeted to solve parabolic problems with standard elliptic operators. However, clearly, this strategy is computationally efficient only if the differential problem is solved in a rectangular domain, if the complete set of eigenfunctions of operator A h are known in advance, and if the FFT technique can be applied.
Another general approach is based on the Lanczos method. The main benefit of this algorithm is that the required matrix-function-vector products are computed without ever forming the dense matrix A h α . A preconditioned Lanczos method can be used to accelerate the convergence of this method (see [7] and references given therein).
Thus, alternative approaches should be used in the case of non-uniform space meshes and/or elliptic operators with non-constant coefficients. The main idea of such techniques is to transform non-local problems to some local classical differential problems. A very good review on these methods is given in [8].
At the end of this short review we note one more important approach. There exist very interesting analytical methods that can solve the same equations in the form of a series on the root vectors of the right-hand [9,10]. They look very promising for the development of not only theoretical results on the existence and uniqueness of solutions, but also for the numerical solution of real world applied problems. A practical comparison of these approaches for well-selected benchmark problems is an important task.
The rest of the paper is organized in the following way. In Section 2, the problem is formulated. The non-stationary parabolic equation with a fractional power elliptic operator is formulated in the 3D parallelepiped, and the initial and boundary conditions are specified. As mentioned above, we use the spectral definition of the fractional power of elliptic operators.
Implicit approximations of parabolic problems with fractional powers of discrete elliptic operators are defined in Section 3. General stability results are presented in the case of a linear source term. Two approaches are proposed on how to construct efficient discrete algorithms. The AAA algorithm and the splitting type scheme are used to define discrete solvers for the non-local time dependent problem. At the end of this section, parallel versions of all algorithms are described and the motivation of the selected approach is discussed.
In Section 4, results of numerical experiments are presented. Some final conclusions are done in Section 5.

2. Problem Formulation

Let Ω be some open and bounded domain Ω R d , α ( 0 , 1 ) . Define a self-adjoint elliptic diffusion operator:
A u = div ( K u ) in Ω
with K ( x ) R d × d symmetric and uniformly positive definite. Operator A is supplemented with either a homogeneous Dirichlet boundary condition on Ω or a homogeneous Neumann one.
For (1), we define the bilinear form
a ( u , v ) = ( K u , v ) ,
where ( · , · ) denotes the L 2 -inner product in Ω and the norm is defined as u = ( u , u ) 1 / 2 .
Let V h H 0 1 ( Ω ) denote a finite-dimensional space of functions that satisfy the homogeneous Dirichlet boundary conditions and we assume that this space is spanned by piecewise linear basis { φ j , j = 1 , , J } . Then, the discrete elliptic operator A h is defined as
( A h U , V ) = a ( U , V ) , U , V V h .
As was mentioned above the spectral definition of fractional power of the operator A h is used in our paper [8,11]. Under standard assumptions on the diffusion coefficients K and the boundary Ω , discrete operator A h admits a system of eigenfunctions Φ j with corresponding eigenvalues λ j > 0 such that
A h Φ j = λ j Φ j , j = 1 , , J .
Then, a fractional power of A h is defined as
A h α U = j = 1 J λ j α ( U , Φ j ) Φ j .
It is easy to show that the nonlocal operator A h α is also self-adjoint and positive definite:
A h α = ( A h α ) * , 0 < λ 1 α I A h α λ J α I .
We solve a Cauchy problem for the evolutionary first-order equation:
d U d t + A h α U = F ( X , t , U ) , 0 < t T , X Ω U ( 0 ) = U 0 , U 0 V h ,
where F represents an associated reaction term.

3. Implicit Approximations of Parabolic Problems with Fractional Powers of Discrete Elliptic Operators

In order to construct a discrete approximation of the differential problem, we define a non-uniform time grid
ω t = { t n : t n = t n 1 + τ n 1 , n = 0 , , N , t 0 = 0 , t N = T } .
Let U n be a numerical approximation to the exact solution U ( t n ) of problem (5). We define the averaging operator
U n + σ = σ U n + 1 + ( 1 σ ) U n , 0 σ 1 .
All parallel algorithms are defined for a general case of σ .
Let us consider the implicit scheme
U n + 1 U n τ n + A h α U n + σ = F ( X , t n + σ , U n + σ ) , n = 0 , , N 1 , U 0 = U 0 , U 0 V h .
It is clear that for σ = 0.5 the symmetrical scheme has the second order of approximation. First, in this paper we always assume that the given problem has a sufficiently smooth solution U ( t ) . Then, applying Taylor series we can compute the approximation (truncation) error
T ( τ ) = τ 2 3 | 3 t 3 U ( t ¯ ) | + τ 2 2 | A h α 2 t 2 U ( t ¯ ) | , t n < t ¯ < t n + 1 .
Second, the discrete fractional power operator A h α is bounded for sufficiently smooth functions and 0 < α < 1 .
We note that the investigation of the approximation accuracy in space is more complicated and results can depend essentially on the fractional order parameter α . These questions are investigated in detail in many papers dealing with discrete approximations of the stationary problems with a fractional power of elliptic operators. We can recommend [12], where h p -adaptive methods are applied (see also references given therein).
For the analysis of time dependent errors we demand a sufficient regularity (in time) of the right-hand side along with compatibility conditions for the initial condition datum. The source functions should satisfy the standard regularity (in time) conditions similar to classical parabolic problem cases. The regularity of F in space is a more complicated question and the appropriate fractional Sobolev spaces are specified in papers dealing with non-stationary problems (see [13]).
We also present results on the stability of the difference scheme (6) in the case of a linear source term F ( x , t ) (the proof is similar to one presented in [6]).
Lemma 1. 
If the weight parameter σ 0.5 , then the difference scheme (6) is unconditionally stable
U n + 1 U n + τ n F ( t n + σ ) , n = 0 , , N 1 .
We repeat our note that, if the direct application of the spectral formula (3) is possible and the FFT algorithm is applicable, then the nonlocal parabolic problem (6) can be solved efficiently at each time layer (see [5,14]). Still, this approach is valid only in very special cases.
In general case, the implementation of the discrete scheme (6) requires one to use some approximations of the operators with fractional powers of discrete elliptic operators A h α . Let us write the discrete scheme in the following form
U n + σ U n σ τ n + A h α U n + σ = F ( X , t n + σ , U n + σ ) .
One possibility to linearize the nonlinear problem (8) is to use the following fixed point iteration method
U n + σ , k U n σ τ n + A h α U n + σ , k = F ( X , t n + σ , U n + σ , k 1 ) , k = 1 , , K ,
where U n + σ , k = σ U n + 1 , k + ( 1 σ ) U n . If the iteration number is restricted to K = 2 , then we get the predictor-corrector method, which is sufficient to guarantee the second order of convergence for the symmetrical scheme with σ = 0.5 .
In this paper we consider 3D elliptic operators. It is well known that for the solution of parabolic problems with standard multidimensional elliptic operators, the most efficient approaches are based on splitting algorithms. In the case of fractional powers of elliptic operators such algorithms are still not constructed and therefore this topic can be formulated as a new important open problem in numerical mathematics.
Here we restrict to general but less efficient algorithms. At each iteration, the solution U n + σ , k is obtained by solving a non-local problem
I + σ τ n A h α U n + σ , k = U n + σ τ n F ( X , t n + σ , U n + σ , k 1 ) .
The first approach to construct an efficient discrete algorithm is based on the idea to approximate the non-local operator ( I + σ τ n A h α ) 1 by some linear local operator B h
I + σ τ n A h α 1 B h .
Then, we compute an approximate solution as
U ˜ n + σ , k = B h U ˜ n + σ τ n F X , t n + σ , U ˜ n + σ , k 1 .
In the next lemma, simple sufficient stability conditions of scheme (11) are formulated in the case of the linear source function F ( X , t n + σ ) (the proof is obtained by modifying the analysis given in [6] for one dimensional operators).
Lemma 2. 
Let us assume that
B h 1 = I + σ τ n C h , C h = C h * > 0 .
Then, for σ 0.5 the difference scheme (11) is unconditionally stable
U ˜ n + 1 U ˜ n + τ n F ( t n + σ ) , n = 0 , , N 1 .
In the second approach, we separately resolve the non-local diffusion and non-linear source terms:
U n + σ , 0 = U n , A h α W n + σ , k = F ( X , t n + σ , U n + σ , k 1 ) , k = 1 , , K , V n , k = U n W n + σ , k ,
V n + σ , k V n , k σ τ n + A h α V n + σ , k = 0 , U n + σ , k = V n + σ , k + W n + σ , k , U n + σ = U n + σ , K .
The main goal of this modification is to approximate more accurately the impact of the non-linear source function, since efficient and accurate solvers are developed for approximation of the fractional powers of stationary diffusion operators. One drawback of this approach is that two subproblems (13) and (14) for different non-local operators should be solved at each iteration. Therefore, it is interesting to investigate if the improvement in the accuracy of approximation compensates for the increase of computational costs of this modified algorithm. It is clear that the standard algorithm of symmetrical splitting with respect to different physical processes can also be used
d U ˜ d t = F ( X , t n + 1 2 , U ˜ ) , U ˜ ( t n ) = U n , t n < t t n + 1 2 , U n + 1 3 = U ˜ ( t n + 1 2 ) ,
U n + 2 3 U n + 1 3 σ τ n + A h α U n + 2 3 + U n + 1 3 2 = 0 ,
d U ˜ d t = F ( X , t n + 1 2 , U ˜ ) , U ˜ ( t n + 1 2 ) = U n + 2 3 , t n + 1 2 < t t n + 1 , U n + 1 = U ˜ ( t n + 1 ) .
This algorithm has two important properties. First, only one subproblem with non-local operators is solved for each time step. Second, the dynamics of the non-linear interaction can be resolved by using specialized solvers targeted for specific non-linear functions, in many cases this part of the problem can even be solved exactly. Therefore this splitting algorithm is a competitive alternative to the other proposed algorithms.
Next we consider an efficient approach for how to construct operators B h to solve the implicit non-local systems (10) and (13) and (14).

3.1. The AAA Algorithm

The construction of B h is based on the so-called AAA algorithm proposed in [15]. We will apply this general algorithm for implementation of the discrete scheme (10). Let us consider a function
f ( z ) = 1 1 + σ τ n z α , z > 0 .
Then, a rational approximation of f ( z ) is given in partial fraction decomposition form [15]
r m ( z ) = N m ( z ) D m ( z ) = c 0 + j = 1 m c j z d j .
The approximate solution of the scheme (10) is computed as
U ˜ n + σ , k = c 0 I + j = 1 m c j A h d j I 1 U ˜ n + σ τ n F ( X , t n + σ , U ˜ n + σ , k 1 ) .
In order to apply the AAA algorithm for the implementation of the discrete scheme (13) and (14), let us consider a function
f 2 ( z ) = z α , z > 0
and compute a rational approximation of f 2 ( z ) given in partial fraction decomposition form
r ˜ m 2 ( z ) = N m 2 ( z ) D m 2 ( z ) = c ˜ 0 + j = 1 m 2 c ˜ j z d ˜ j .
The approximate solution of the scheme (13) and (14) is computed as
U ˜ n + σ , 0 = U ˜ n , W ˜ n + σ , k = c ˜ 0 I + j = 1 m 2 c ˜ j A h d ˜ j I 1 F ( X , t n + σ , U ˜ n + σ , k 1 ) , V ˜ n = U ˜ n W ˜ n + σ , k ,
V ˜ n + σ , k = c 0 I + j = 1 m c j A h d j I 1 V ˜ n , U ˜ n + σ , k = V ˜ n + σ , k + W ˜ n + σ , k , k = 1 , , K , U ˜ n + 1 = 2 U ˜ n + σ , K U ˜ n .

Stability Analysis

The approximation of the AAA algorithm is based on a rational approximation of function f ( z ) . Let us consider the solution of Equation (21) for F = 0 , and represent it as a linear combination of the eigenvectors of operator A h
U ˜ n = j = 1 J c j n Φ j .
Substituting this expression into (21) we get the equations for coefficients c j :
σ c j n + 1 + ( 1 σ ) c n = r m ( λ j ) c j , j = 1 , , J .
We assume that z L λ j z R . Then the stability of scheme (19) is guaranteed if the following condition is satisfied
max z L z z R 1 σ | r m ( z ) ( 1 σ ) | 1 .
As in [6], we investigate this condition experimentally. The results of numerical experiments confirmed that for fractional powers α = 0.25 , 0.5 , 0.75 , m = 4 , , 15 , time steps τ = 10 l , l = 1 , , 5 , σ = 0.5 and z L , z R defined by the minimal and maximal eigenvalues of the operator A h , for N = 128 , 256 , 512 , 1024 , the stability condition (23) was satisfied unconditionally. Such stability tests are recommended before starting large scale computations for simulations of real world applications.

3.2. Symmetrical Splitting Scheme

By using the same AAA algorithm, we get an implementation of the symmetrical splitting scheme
d U ˜ d t = F ( X , t n + 1 2 , U ˜ ) , U ˜ ( t n ) = U ˜ n , t n < t t n + 1 2 , U ˜ n + 1 3 = U ˜ ( t n + 1 2 ) ,
U ˜ n + 2 3 = 2 c 0 I + j = 1 m c j A h d j I 1 U ˜ n + 1 3 U ˜ n + 1 3 ,
d U ˜ d t = F ( X , t n + 1 2 , U ˜ ) , U ˜ ( t n + 1 2 ) = U ˜ n + 2 3 , t n + 1 2 < t t n + 1 , U ˜ n + 1 = U ˜ ( t n + 1 ) .
The non-linear problems (24) and (26) are solved by applying an ODE solver optimal for a given type of non-linearities.

3.3. Additive Splitting Scheme

This scheme is constructed using techniques from [16] where BURA type discrete approximations are constructed for stationary problems with fractional powers of elliptic operators. We note that similar splitting schemes were used also in [17].
First we rewrite the linearized Equation (9) in the following form
U n + σ , k U n σ τ n + A h β A h U n + σ , k = F ( X , t n + σ , U n + σ , k 1 ) , k = 1 , , K ,
where β = 1 α . Next, the AAA method is used to construct a rational approximation of fractional power operator A h β
A h β R m ( β ) = c 0 I + j = 1 m c j A h d j I 1 .
Let us write the obtained discrete scheme as
U n + σ , k U n σ τ n + j = 0 m c j D h , j 1 A h U n + σ , k = F ( X , t n + σ , U n + σ , k 1 ) , k = 1 , , K ,
where
D h , 0 = I , D h , j = A h d j I , j = 1 , , m .
Here we use the same notation U n for a new discrete solution.
Since all discrete operators commute, the solution U n + 1 , k is computed by using the classical splitting technique in m + 1 steps:
U ˜ n 1 / m , k = U ˜ n , D h , j U ˜ n + j / m , k U ˜ n + ( j 1 ) / m , k τ n + c j A h σ U ˜ n + j / m , k + ( 1 σ ) U ˜ n + ( j 1 ) / m , k = δ 0 j D h , j F ( X , t n + σ , U ˜ n + σ , k 1 ) , j = 0 , , m .
The splitting technique introduces additional approximation errors and therefore this algorithm can require to use smaller time steps. This scheme is unconditionally stable but again the transfer operator is different than in the Crank–Nicolson method.
As was done for the AAA algorithm, we consider a modification of this scheme when the influence of the non-homogeneous source term is computed separately by solving Equation (21):
U ˜ n + σ , 0 = U ˜ n , W ˜ n + σ , k = c ˜ 0 I + j = 1 m 2 c ˜ j A h d ˜ j I 1 F ( X , t n + σ , U ˜ n + σ , k 1 ) , V n 1 / m , k = U ˜ n W ˜ n + σ , k ,
D h , j V ˜ n + j / m , k V n + ( j 1 ) / m , k σ τ n + c j A h V ˜ n + j / m , k = 0 , j = 0 , , m , V n + j / m , k = 2 V ˜ n + j / m , k V n + ( j 1 ) / m , k , j = 0 , , m , U ˜ n + σ , k = V ˜ n + 1 , k + W ˜ n + σ , k , k = 1 , , K , U ˜ n + 1 = V n + 1 , K + W ˜ n + σ , K .

3.4. Parallel Algorithms

All parallel solvers are constructed by using the same approach, when parallelization is done at the level of solvers for solution of the obtained systems of linear equations. Thus we can use efficient parallel solvers developed for such problems.
Our aim is to investigate how the efficiency of these solvers depend on specific properties of discrete elliptic operators ( A d j I ) and parameters d j . The values of these parameters depend on the selected non-local operators and on the number m in approximations (19) and (20).
We note one important difference among the developed parallel algorithms. At each time step, the AAA algorithms (19), (21), and (22), and symmetrical splitting schemes (24)–(26) define m independent systems of linear equation, and these problems can be solved in parallel. Thus, an additional parallelization level is obtained and the distribution of tasks can be optimized by using general parallelization templates (see [18]).
While for the additive splitting scheme (27) all auxiliary problems should be solved sequentially, one after another. For this method, parallelization can be done only on the level of a solution of linear systems.
It is important to note that the total number of non-local problems solved at one time step is also different for specific schemes. Two problems are solved for modified AAA and additive splitting schemes, while one problem is sufficient for the remaining schemes. In the case of non-linear source functions, the predictor-corrector iterations are done and the number of non-local problems solved at one time step is increased twice. The only exception is the symmetrical splitting scheme, where again it is sufficient to solve only one non-local problem per time step.
As our tool for development parallel versions of the constructed finite difference schemes, we use AGMG multigrid solver [19]. The main reason to base the construction of parallel solvers on some state of the art tool is that general mathematical libraries are adjusted and optimized to new computer architectures by library developers. At the same time most tools have a sufficiently large set of parameters to fit these libraries to specific needs arising when solving new classes of problems. We also note that in [20], an extended analysis of parallel solvers is done for the solution of stationary fractional power elliptic problems. The hyper AMG solver BoomerAMG was used in [20] to construct parallel versions of the finite difference schemes for a 7-point 3D Laplace problem.
For comparison of different approaches we have implemented and investigated the scalability of symmetrical splitting scheme when the spectral solver based on FFT is used. The parallelization is done by using the well-known library FFTW [21].

4. Numerical Experiments

In this section we solve three test problems. The first one defines a nonlinear problem and is used to analyze the accuracy of different schemes. The second problem has a non-smooth solution and the scalability analysis of the parallel algorithms is done on its basis. The last system of coupled 3D nonlinear parabolic equations enables to analyze the influence of the fractional order parameter on complicated pattern formation. In addition, this problem is used to investigate the parallel efficiency of the symmetric splitting scheme when the non-local operators are resolved by using an FFT algorithm.

4.1. Convergence Tests

In order to investigate the convergence accuracy of the proposed discrete schemes and integration methods, we consider a linear fractional heat equation. Its implementation still requires to apply the predictor-corrector iterative algorithm.
Example 1. 
In this example we solve a linear fractional heat equation with a linear source term in three-dimensional space Ω = [ 0 , 1 ] × [ 0 , 1 ] × [ 0 , 1 ] (for a similar two-dimensional example, see [5]):
u t + ( Δ ) α u = F ( x , y , z , t , u ) , ( x , y , z ) Ω , F ( x , y , z , t , u ) = t 2 α i = 1 8 β i λ i α v i + 2 α t 2 α 1 + t 2 α sin 3 ( π x ) sin 3 ( π y ) sin 3 ( π z ) u , u ( x , y , z , t ) = 0 , ( x , y , z ) Ω ,
where v 1 = sin ( π x ) sin ( π y ) sin ( π z ) , λ 1 = 3 π 2 , β 1 = 27 / 64 , v 2 = sin ( π x ) sin ( π y ) sin ( 3 π z ) , λ 2 = 11 π 2 , β 2 = 9 / 64 , v 3 = sin ( π x ) sin ( 3 π y ) sin ( π z ) , λ 3 = 11 π 2 , β 3 = 9 / 64 , v 4 = sin ( π x ) sin ( 3 π y ) sin ( 3 π z ) , λ 4 = 19 π 2 , β 4 = 3 / 64 , v 5 = sin ( 3 π x ) sin ( π y ) sin ( π z ) , λ 5 = 11 π 2 , β 5 = 9 / 64 ,, v 6 = sin ( 3 π x ) sin ( π y ) sin ( 3 π z ) , λ 6 = 19 π 2 , β 6 = 3 / 64 , v 7 = sin ( 3 π x ) sin ( 3 π y ) sin ( π z ) , λ 7 = 19 π 2 , β 7 = 3 / 64 , v 8 = sin ( 3 π x ) sin ( 3 π y ) sin ( 3 π z ) , λ 8 = 27 π 2 , β 8 = 1 / 64 .
Then this problem has the exact solution
u ( x , y , t ) = t 2 α sin 3 ( π x ) sin 3 ( π y ) sin 3 ( π z ) .
We define uniform space meshes in each direction Ω ¯ h = ω ¯ x × ω ¯ y × ω ¯ z :
ω ¯ x = { x i : x i = i h x , i = 0 , , J x , h x = 1 / J x } , ω ¯ y = { y j : y j = j h y , j = 0 , , J y , h y = 1 / J y } , ω ¯ z = { z j : z k = k h z , k = 0 , , J z , h z = 1 / J z }
and approximate the differential diffusion operator by using the standard central differences
A h U = ( U i + 1 , j , k 2 U i j k + U i 1 , j , k h x 2 + U i , j + 1 , k 2 U i j k + U i , j 1 , k h y 2 + U i , j , k + 1 2 U i j k + U i , j , k 1 h z 2 ) .
The predictor-corrector method takes into account the dependence of the source term on the solution. The discrete solution is computed using the spectral method, since only eight spectral modes should be resolved. The obtained results give reference solutions for the estimation of the accuracy of more general integration methods.
We have investigated the space and time discretization errors separately. First, we have fixed α = 0.5 , J x = J y = J z = 1024 and solved the discrete problem till T = 1 with different sizes of time step τ = 0.1 , 0.05 , 0.025 . The computed maximum errors
e = max ( x i , y j , z k ) Ω h | U i j k U ˜ i j k |
confirm the second order accuracy of the approximation:
e ( 0.1 ) = 2.1089 × 10 4 , e ( 0.05 ) = 5.8807 × 10 5 , e ( 0.025 ) = 1.4636 × 10 5 .
Next, we have fixed τ = 0.002 and solved the same problem for different space meshes J x = J y = J z . The results of computational experiments again confirmed the second order accuracy of the approximation:
e ( 16 ) = 6.659 × 10 3 , e ( 32 ) = 1.658 × 10 3 , e ( 64 ) = 4.141 × 10 4 , e ( 128 ) = 1.0341 × 10 4 .
Note, that in [5] the accuracy of the 2D solver is tested only for fractional order parameters α = 1 , 0.85 , 0.75 . It is obvious that the smoothness of the solution changes for smaller values of the fractional order and for α = 0.25 the time derivative of the solution is singular at the initial time moment t = 0 .
In order to test the accuracy of the finite difference scheme we have fixed α = 0.25 , J x = J y = J z = 1024 and solved the discrete problem till T = 1 with different sizes of time step τ :
e 1 20 = 1.4157 × 10 3 , e 1 40 = 1.0443 × 10 3 , e 1 160 = 5.3554 × 10 4 , e 1 640 = 2.6897 × 10 4 .
The computed maximum errors show a reduced convergence rate O ( τ 0.5 ) of the discrete solution. Therefore, some adaptive (non-uniform) time step should be used to resolve this singularity.
It is clear that such asymptotical dynamics of the error depends mainly on the singularity of the exact solution, but not on the fractional power parameter α . In order to illustrate this statement we solved the parabolic problem (30) when α = 0.25 , but the source function F is selected such that the exact solution is defined as
u ( x , y , t ) = t sin 3 ( π x ) sin 3 ( π y ) sin 3 ( π z ) .
The errors of the discrete solution in the maximum norm are the following
e 1 10 = 5.7916 × 10 4 , e 1 20 = 1.5065 × 10 4 , e 1 40 = 3.8083 × 10 5 , e 1 80 = 9.2027 × 10 6 .
It follows from the computational experiments that the second order convergence rate is restored for the discrete scheme also for α = 0.25 .
Next we present results of computational experiments when different algorithms are used to approximate non-local fractional power discrete elliptic operators. In all experiments the space mesh size is fixed to J x = J y = J z = 256 and the test problem is solved till T = 1 with the time step τ = 1 256 .

The AAA Algorithm

First we consider the AAA algorithm (19). In order to apply the AAA approximation method. We sample function ( 1 + 0.5 τ z α ) 1 with M = 50,000 points over the exact interval of eigenvalues of A h operator [ λ h m i n , λ h m a x ] . The resulting errors for three values of fractional order parameters α = 0.75 , 0.5 , 0.25 and various orders m of the rational functions r m ( z ) are presented in Table 1. For comparison we also present the error of the discrete solution of the symmetric Crank–Nicolson finite difference scheme, when at each time step the non-local problem is solved by using FFT method.
It follows from the presented results, that a small number m 10 is sufficient to reach the optimal error levels. These values of m do not depend on the fractional parameter α .
Next we consider the modified AAA algorithms (21) and (22). The resulting errors for α = 0.75 , 0.5 , 0.25 and various m of the rational functions r m ( z ) are presented in Table 2.
It follows from the presented results, that a small number m 10 is again sufficient to reach the optimal error levels. The accuracy of the modified AAA algorithm is a little better for most cases of parameters.
Next we consider the modified additive splitting algorithm (28) and (29). The resulting errors for α = 0.75 , 0.5 , 0.25 and various orders m of the rational functions r m ( z ) are presented in Table 3.
The accuracy of the modified additive splitting scheme (28) and (29) is better than for the modified AAA algorithm for fractional order parameters α = 0.5 and α = 0.75 .

4.2. Scalability of the Parallel Algorithms

Example 2. 
In this example we solve the following 3D linear problem
u t + ( Δ ) α u = F ( x , y , z ) , ( x , y , z ) Ω = ( 0 , 1 ) × ( 0 , 1 ) × ( 0 , 1 ) , u ( x , y , z , 0 ) = 0 , ( x , y , z ) Ω , u ( x , y , z , t ) = 0 , ( x , y , z ) Ω , t > 0 ,
where the right-hand-side F is the well-known checkerboard function:
F ( x , y , z ) = 1 , i f ( x 0.5 ) ( y 0.5 ) ( z 0.5 ) > 0 ; 0 , o t h e r w i s e .
We note that for stationary problems, similar test problems (mostly in 2D domains) with checkerboard right-hand-side functions are used to demonstrate the convergence of the selected numerical methods in many papers (see [11,16,20] and references given therein). This type of source functions, but for a 1D case, was also used in [6] for testing the accuracy of different integration methods for parabolic problems with fractional power elliptic operators.
We note that all constructed finite difference algorithms have the same structure, when at each time step a set of systems of linear equations is solved, where each system approximates classical diffusion operators perturbed by adding scaled identity operators.
As explained above, the AAA algorithms (19), (21), and (22), and the symmetrical splitting scheme (24)–(26) define sets of independent systems of linear equations, and these problems can be solved in parallel. While for the additive splitting scheme (27) all auxiliary problems should be solved sequentially, one after another.
In this paper we restrict ourselves to the scalability analysis of parallel algorithms when all processes are used to solve one system of linear equations. This strategy has two advantages: first, it is possible to solve larger systems of equations (i.e., to get more accurate approximations); second, data communication costs are reduced essentially.
In all computational experiments we restrict to the strong scalability analysis when the size of the discrete problem is fixed to the 201 × 201 × 201 grid.
Parallel numerical tests in this work were performed on the computer cluster “HPC Vanagas” at the High Performance Computing Laboratory of Vilnius Gediminas technical university. We have used up to 16 cores of Intel® Xeon® processor E5-2630 with 20 cores ( 2.20 GHz) and 32 GB DDR4 of RAM. The parallel algorithms are implemented using the M P I library.
In order to estimate the scalability of parallel algorithms for the solution of parabolic problems with fractional order elliptic operators as a benchmark, we use results obtained for a classical Crank–Nicolson finite difference scheme. The parabolic Poisson problem with a 7-point 3D discrete Laplace operator on a J × J × J grid is solved and the AGMG multigrid solver is used to solve systems of linear equations with a banded matrix.
The discrete problem is solved on the J x = J y = J z = 201 grid and the time step is fixed to τ = 0.5 . Five time steps are computed in all the experiments of this example.
The results of the computational experiments are presented in Table 4. Here S p = T 1 / T p is the speedup of the parallel algorithm, E p = S p / p is the efficiency number, and T p is the parallel run-time on p cores.
The presented results define the reference estimates of the AGMG solver on this parallel computer. They will be used to analyze the strong scalability of parallel solvers constructed to solve problems with fractional power elliptic operators. A small reduction in the efficiency of the parallel solver for p = 2 processes can be explained as a possible switching of the system from two processes running on two physical cores to two virtual processes running on one physical core (so-called multi-threading).

Strong Scalability Analysis of the Parallel AAA Algorithm

Computational experiments are done for two fractional order parameters α = 0.5 , 0.25 and different numbers of cores p = 1 , 2 , 4 , 8 , 16 . The order of AAA algorithm is fixed to m = 10 . Results of the experiments are presented in Table 5.
The obtained results show that the run-time T p is decreased in comparison with results in Table 4. The improvement is due to the faster convergence of the AGMG solver when a positive operator ( d j ) I is added to the elliptic operator A h (all coefficients d j are negative for the constructed AAA approximation of the function f ( z ) = ( 1 + 0.5 τ z α ) 1 ). The efficiency numbers E p of the parallel AAA discrete scheme are slightly improved, this fact confirms the robustness of the AGMG solver.
Next we estimated the influence of the order of AAA algorithm m, since for different values of m we get different sets of diagonalization coefficients d j . In the following experiments m is fixed to m = 5 and fractional order parameter α = 0.5 . The results of the computational experiments are presented in Table 6.
As expected, the run-time T p is reduced and even more than twice, but the efficiency of the parallel algorithm is changed only slightly.
Next we investigated the influence of the time step τ size. The order of AAA algorithm is taken m = 10 , the fractional parameter α = 0.5 , and the time step is increased to τ = 0.025 . Five time steps are computed, as in previous experiments. The results of the computational experiments are presented in Table 7.
We see that the run-time is increased when the time step size is increased, since the AGMG solver converges slower. However, the efficiency of the parallel algorithm changes only very slightly.
It follows from the presented results of computational experiments that the efficiency of the parallel AAA algorithm is robust with respect to fractional power parameter, the order of AAA algorithm parameter m and discrete time step τ variations.
Example 3. 
In this example, we consider the 3D space-fractional Gray–Scott model (see [14] for 2D model). It is described by a system of coupled 3D diffusion-reaction equations
u t + K u ( Δ ) α u = ( 1 u ) v 2 F u , ( x , y , z ) Ω = ( 0 , 1 ) × ( 0 , 1 ) × ( 0 , 1 ) , v t + K v ( Δ ) α v = ( 1 u ) v 2 ( F + λ ) v ,
where the diffusion rates K u = 2 × 10 5 , K v = 1 × 10 5 , the dimensionless feed rate F = 0.03 , and the dimensionless decay rate constant of the second reaction λ = 0.061 .
The homogeneous Direchlet boundary conditions are defined on the boundary Ω , therefore a new function u = 1 u ˜ is used in the mathematical model (34) (compare it with [14]). The initial conditions are given by
u ( x , y , z , 0 ) = 0.5 , v ( x , y , z , 0 ) = 0.25 , ( x , y , z ) Ω 1 , u ( x , y , z , 0 ) = 0 , v ( x , y , z , 0 ) = 0 , ( x , y , z ) Ω \ Ω 1 ,
where Ω 1 = { ( x , y , z ) : ( x 0.5 ) 2 + ( y 0.5 ) 2 + ( z 0.5 ) 2 0 . 04 2 } .
It is well-known that the solution of the Gray-Scott model exhibits a complicated set of patterns and the process of pattern formation is far from trivial. First, we would like to note that the specific trend of patterns depends essentially on the discrete size of the model (i.e., the space mesh size), therefore, we here fix the mesh to J x = J y = J z = 256 . Our main goal is to show the dependents of the solution on the fractional power parameter α (i.e., to evaluate the effects of the super-diffusion). A second aim is to investigate the efficiency of the parallel symmetrical splitting algorithm for solving this complicated system of three dimensional non-linear and non-local parabolic equations. Figure 1 presents the evolution of the numerical solution v for three fractional orders α = 1 , 0.85 , 0.75 . The contours of v are shown at times t = 2000 , 3000 from top to bottom.
We see from the presented results that the way of the replication in pattern formation depends essentially on parameter α .
We also investigated the scalability of the parallel symmetrical splitting algorithm (15)–(15). The non-local fractional power operator is approximated by using the spectral method. The FFTW library [21] is used to solve the obtained discrete problems. The results of the computational experiments are presented in Table 8. The fractional order parameter α = 0.85 and the four time steps are computed.
It follows from the presented results that for the given parallel computer, the efficiency of the spectral FFTW solver is better than the efficiency of the GAMG solver. Still, the latter solver can be used to implement the constructed general finite difference schemes, when the spectral method is not applicable. We note that computations of images presented in Figure 1 required up to 3 h of run-time on 20 cores system in order to solve one case of the problem.

5. Conclusions

In this paper, we propose efficient parallel finite difference schemes for computing the solutions of 3D parabolic problems with a fractional power of elliptic operators and non-linear source terms. These algorithms are targeted to solve problems with general elliptic operators in bounded domains. Discrete schemes include approximations of non-local operators by using the AAA algorithm. Three different approaches are used to construct efficient implicit in time schemes. The stability of schemes and convergence rates of the discrete solution are proved.
The main aim is to investigate the scalability of parallel solvers constructed for each discrete scheme. The general AGMG multigrid solver is selected as a basic solver and it is shown that the obtained parallel algorithms are robust to changes of parameters of the differential problem and discrete schemes.
The results of the extended experiments are reported to illustrate the accuracy and efficiency of the constructed algorithms. In particular, applications to non-stationary fractional Poisson-like equations and systems of such equations are given. In a forthcoming paper, we will investigate some realistic applications for the simulation of smart biosensors, when fractional power diffusion operators describe more accurately the diffusion processes.

Author Contributions

R.Č. (Raimondas Čiegis): conceptualization; methodology; writing—original draft preparation; numerical algorithms; R.Č. (Remigijus Čiegis): analysis; validation; writing—editing; I.D.: software; computational experiments, analysis of results. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Not applicable.

Acknowledgments

The authors are thankful to the reviewers for their remarks, which helped to improve the presentation of results.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Baeumer, B.; Benson, D.A.; Meershaert, M.M.; Wheatcraft, S.W. Subordinated advection-dispersion equation for contaminant transport. Water Resour. Res. 2001, 37, 1543–1550. [Google Scholar] [CrossRef]
  2. Metzler, R.; Klafter, J. The restaurant at the end of the random walk: Recent developments in the description of anomalous transport by fractional dynamics. J. Phys. A Math. Gen. 2004, 37, R161–R208. [Google Scholar] [CrossRef]
  3. Bueno-Orovio, A.; Kay, D.; Grau, V.; Rodriguez, B.; Burrage, K. Fractional diffusion models of cardiac electrical propagation: Role of structural heterogeneity in dispersion of repolarization. J. R. Soc. Interface 2014, 11, 20140352. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  4. Harizanov, S.; Margenov, S.; Marinov, P.; Vutov, Y. Volume constrained 2-phase segmentation method utilizing a linear system solver based on the best uniform polynomial approximation of x−1/2. J. Comput. Appl. Math. 2017, 310, 115–128. [Google Scholar] [CrossRef]
  5. Lee, H.G. A second-order operator splitting Fourier spectral method for fractional-in-space reaction–diffusion equations. J. Comput. Appl. Math. 2018, 333, 395–403. [Google Scholar] [CrossRef]
  6. Čiegis, R.; Čiegis, R.; Dapšys, I. A Comparison of discrete schemes for numerical solution of parabolic problems with fractional power elliptic operators. Mathematics 2021, 9, 1344. [Google Scholar] [CrossRef]
  7. Yang, Q.; Turner, I.; Moroney, T.; Liu, F. A finite volume scheme with preconditioned Lanczos method for two-dimensional space-fractional reaction-diffusion equations. Appl. Math. Model. 2014, 38, 3755–3762. [Google Scholar] [CrossRef] [Green Version]
  8. Hofreither, C. A unified view of some numerical methods for fractional diffusion. Comput. Math. Appl. 2020, 80, 332–350. [Google Scholar] [CrossRef]
  9. Kukushkin, M.V. Natural lacunae method and Schatten-von Neumann classes of the convergence exponent. arXiv 2021, arXiv:2112.10396v2. [Google Scholar]
  10. Kukushkin, M.V. Evolution equations in Hilbert spaces via the lacunae method. arXiv 2022, arXiv:2202.07338. [Google Scholar]
  11. Bonito, A.; Pasciak, J.E. Numerical approximation of fractional powers of elliptic operators. Math. Comput. 2015, 84, 2083–2110. [Google Scholar] [CrossRef] [Green Version]
  12. Melenk, J.M.; Rieder, A. hp-FEM for the fractional heat equation. IMA J. Numer. Anal. 2021, 41, 412–454. [Google Scholar] [CrossRef]
  13. Nochetto, R.H.; Otarola, E.; Salgado, A.J. A PDE approach to space-time fractiobal parabolic problems. SIAM J. Numer. Anal. 2016, 54, 848–873. [Google Scholar] [CrossRef] [Green Version]
  14. Zhang, H.; Jianga, X.; Zeng, F.; Karniadakis, G.E. A stabilized semi-implicit Fourier spectral method for nonlinear space-fractional reaction-diffusion equations. J. Comput. Phys. 2020, 405, 109141. [Google Scholar] [CrossRef] [Green Version]
  15. Nakatsukasa, Y.; Sete, O.; Trefethen, L.N. The AAA algorithm for rational approximation. SIAM J. Sci. Comput. 2018, 40, A1494–A1522. [Google Scholar] [CrossRef]
  16. Harizanov, S.; Lazarov, R.; Margenov, S.; Marinov, P.; Vutov, Y. Optimal solvers for linear systems with fractional powers of sparse SPD matrices. Numer. Linear Algebra Appl. 2018, 25, e2167. [Google Scholar] [CrossRef]
  17. Vabishchevich, P.N. Splitting schemes for non-stationary problems with a rational approximation for fractional powers of the operator. Appl. Numer. Math. 2021, 165, 414–430. [Google Scholar] [CrossRef]
  18. Kriauzienė, R.; Bugajev, A.; Čiegis, R. A three-level parallelisation scheme and application to the Nelder-Mead algorithm. Math. Model. Anal. 2020, 25, 584–607. [Google Scholar] [CrossRef]
  19. Notay, Y. Aggregation-based algebraic multigrid for convection-diffusion equations. SIAM J. Sci. Comput. 2012, 34, A2288–A2316. [Google Scholar] [CrossRef] [Green Version]
  20. Čiegis, R.; Starikovičius, V.; Margenov, S.; Kriauzienė, R. Parallel solvers for fractional power diffusion problems. Concurr. Comput. Pract. Exp. 2017, 29, e4216. [Google Scholar] [CrossRef]
  21. Johnson, S.G.; Frigo, M. A modified split-radix FFT with fewer arithmetic operations. IEEE Trans. Signal Process. 2007, 55, 111–119. [Google Scholar] [CrossRef]
Figure 1. Dynamics of the v wave in the Gray-Scott model at different times t = 2000 , 3000 . The v contours are shown from top to bottom for three fractional orders: (a) α = 1 , (b) α = 0.85 , (c) α = 0.75 .
Figure 1. Dynamics of the v wave in the Gray-Scott model at different times t = 2000 , 3000 . The v contours are shown from top to bottom for three fractional orders: (a) α = 1 , (b) α = 0.85 , (c) α = 0.75 .
Axioms 11 00098 g001
Table 1. The AAA algorithm: errors of the solution of the discrete scheme (19) for α = 0.75 , 0.5 , 0.25 and a sequence of parameters m. The column FFT defines errors of the spectral solver.
Table 1. The AAA algorithm: errors of the solution of the discrete scheme (19) for α = 0.75 , 0.5 , 0.25 and a sequence of parameters m. The column FFT defines errors of the spectral solver.
α m = 5 m = 7 m = 10 m = 12 m = 15 FFT
0.75 6.62 × 10 4 1.60 × 10 5 4.20 × 10 5 4.13 × 10 5 4.57 × 10 5 4.15 × 10 5
0.5 6.92 × 10 4 8.85 × 10 5 2.58 × 10 5 2.39 × 10 5 3.10 × 10 5 2.55 × 10 5
0.25 1.27 × 10 3 5.09 × 10 4 4.14 × 10 4 4.17 × 10 4 4.04 × 10 4 4.15 × 10 4
Table 2. The modified AAA algorithm: errors of the solution of the discrete scheme (21) and (22) for α = 0.75 , 0.5 , 0.25 and a sequence of parameters m.
Table 2. The modified AAA algorithm: errors of the solution of the discrete scheme (21) and (22) for α = 0.75 , 0.5 , 0.25 and a sequence of parameters m.
α m = 5 m = 7 m = 10 m = 12 m = 15
0.75 9.960 × 10 5 4.430 × 10 5 4.119 × 10 5 4.175 × 10 5 4.102 × 10 5
0.5 1.982 × 10 4 2.536 × 10 5 2.543 × 10 5 2.582 × 10 5 2.421 × 10 5
0.25 2.351 × 10 4 3.970 × 10 4 4.150 × 10 4 4.140 × 10 4 4.189 × 10 4
Table 3. The modified additive splitting algorithm: errors of the solution of the discrete scheme (28) and (29) for α = 0.75 , 0.5 , 0.25 and a sequence of parameters m.
Table 3. The modified additive splitting algorithm: errors of the solution of the discrete scheme (28) and (29) for α = 0.75 , 0.5 , 0.25 and a sequence of parameters m.
α m = 5 m = 7 m = 10 m = 12 m = 15
0.75 3.144 × 10 5 7.124 × 10 6 8.155 × 10 6 7.866 × 10 6 8.189 × 10 6
0.5 1.338 × 10 4 2.299 × 10 5 1.569 × 10 5 1.561 × 10 5 1.553 × 10 5
0.25 2.930 × 10 4 4.153 × 10 4 4.163 × 10 4 4.163 × 10 4 4.164 × 10 4
Table 4. The Crank–Nicolson discrete scheme: the values of the speedup S p and efficiency E p numbers for a sequence of cores p.
Table 4. The Crank–Nicolson discrete scheme: the values of the speedup S p and efficiency E p numbers for a sequence of cores p.
p = 1 p = 2 p = 4 p = 8 p = 16
T p 344.9193.591.7163.6856.65
S p 1.01.7823.7615.4166.088
E p 0.8910.9400.6780.381
Table 5. The AAA algorithm (19): the values of the speedup S p and efficiency E p numbers for α = 0.5 , 0.25 and different numbers of cores. The order of AAA algorithm is fixed to m = 10 .
Table 5. The AAA algorithm (19): the values of the speedup S p and efficiency E p numbers for α = 0.5 , 0.25 and different numbers of cores. The order of AAA algorithm is fixed to m = 10 .
p = 1 p = 2 p = 4 p = 8 p = 16
α = 0.5
T p 267.9149.172.0448.4944.24
S p 1.01.7973.7195.5256.056
E p 0.8990.9300.6900.378
α = 0.25
T p 267.6149.271.6348.6544.05
S p 1.01.7973.7365.5016.075
E p 0.8990.9340.6880.380
Table 6. The AAA algorithm (19): the values of the speedup S p and efficiency E p numbers for α = 0.5 and different numbers of cores. The order of the AAA algorithm is fixed to m = 5 .
Table 6. The AAA algorithm (19): the values of the speedup S p and efficiency E p numbers for α = 0.5 and different numbers of cores. The order of the AAA algorithm is fixed to m = 5 .
p = 1 p = 2 p = 4 p = 8 p = 16
T p 107.161.2729.7620.1518.14
S p 1.01.7483.5995.3155.904
E p 0.8740.9000.6520.369
Table 7. The AAA algorithm (19): the values of the speedup S p and efficiency E p numbers for τ = 0.025 , α = 0.5 and different numbers of cores. The order of the AAA algorithm is fixed to m = 10 .
Table 7. The AAA algorithm (19): the values of the speedup S p and efficiency E p numbers for τ = 0.025 , α = 0.5 and different numbers of cores. The order of the AAA algorithm is fixed to m = 10 .
p = 1 p = 2 p = 4 p = 8 p = 16
T p 293.9163.9679.3454.249.61
S p 1.01.7933.7045.4225.924
E p 0.8960.9260.6780.370
Table 8. The spectral symmetrical splitting algorithm: the values of the speedup S p and efficiency E p numbers for α = 0.85 and different numbers of cores.
Table 8. The spectral symmetrical splitting algorithm: the values of the speedup S p and efficiency E p numbers for α = 0.85 and different numbers of cores.
p = 1 p = 2 p = 4 p = 8 p = 16
T p 173.190.1246.3222.1512.84
S p 1.01.9213.7377.81513.481
E p 0.9600.9340.9770.843
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Čiegis, R.; Dapšys, I.; Čiegis, R. A Comparison of Parallel Algorithms for Numerical Solution of Parabolic Problems with Fractional Power Elliptic Operators. Axioms 2022, 11, 98. https://doi.org/10.3390/axioms11030098

AMA Style

Čiegis R, Dapšys I, Čiegis R. A Comparison of Parallel Algorithms for Numerical Solution of Parabolic Problems with Fractional Power Elliptic Operators. Axioms. 2022; 11(3):98. https://doi.org/10.3390/axioms11030098

Chicago/Turabian Style

Čiegis, Raimondas, Ignas Dapšys, and Remigijus Čiegis. 2022. "A Comparison of Parallel Algorithms for Numerical Solution of Parabolic Problems with Fractional Power Elliptic Operators" Axioms 11, no. 3: 98. https://doi.org/10.3390/axioms11030098

APA Style

Čiegis, R., Dapšys, I., & Čiegis, R. (2022). A Comparison of Parallel Algorithms for Numerical Solution of Parabolic Problems with Fractional Power Elliptic Operators. Axioms, 11(3), 98. https://doi.org/10.3390/axioms11030098

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop