Predictor-Corrector Domain Decomposition Algorithm for Parabolic Problems on Graphs

In this paper, we present a predictor-corrector type algorithm for solution of linear parabolic problems on graph structure. The graph decomposition is done by dividing some edges and therefore we get a set of problems on sub-graphs, which can be solved efficiently in parallel. The convergence analysis is done by using the energy estimates. It is proved that the proposed finite-difference scheme is unconditionally stable but the predictor step error gives only conditional approximation. In the second part of the paper it is shown that the presented algorithm can be written as Douglas type scheme, based on the domain decomposition method. For a simple case of one dimensional parabolic problem, the convergence analysis is done by using results from [P. Vabishchevich. A substracturing domain decomposition scheme for unsteady problems. Comp. Meth. Appl. Math. 11(2):241-268, 2011]. The optimality of asymptotical error estimates is investigated. Results of computational experiments are presented.


Introduction
The predictor-corrector splitting and domain decomposition methods are widely used to solve elliptic and parabolic PDEs in multidimensional domains [3,7,8,15]. We also note that new developments of these ideas are presented in [13,14], where non-standard domain decomposition methods with overlapping and non-overlapping subdomains are investigated.
Domain Decomposition and predictor-corrector techniques can be applied for problems on graph structures (see [1,2,9,11] for related discussions). In [1], it is proposed to divide the graph by decomposing the set of graph vertexes. The stability analysis is done in the maximum norm, but only conditional stability can be proved by this method and no influence of the corrector step on the stability is obtained.
In this paper we consider predictor-corrector type finite difference schemes for solution of linear parabolic problems on graph domains. The graph decomposition is done by dividing some edges of the graph, therefore a set of problems on sub-graphs is obtained and these subproblems can be solved efficiently in parallel. The main goal is to investigate the stability of the obtained discrete algorithms with respect to the approximation errors introduced at the predictor step of the algorithm.
The rest of the paper is organized as follows. In Section 2, the mathematical model is formulated, and the predictor-corrector finite difference scheme is constructed. The stability and convergence analysis is presented in Section 3. By using the energy method it is proved that the finite difference scheme is unconditionally stable but the predictor step error gives only conditional approximation. In Section 4, it is shown that the presented algorithm can be written as Douglas type scheme. The decomposition operators are based on a partition of unity over the computational domain. For a simple case of one dimensional parabolic problems, the convergence analysis is done by using results from [13]. The optimality of asymptotical error estimates is investigated. Results of computational experiments are presented. Finally, some conclusions are given in Section 5.

Formulation of the Problem
In this section we present a formulation of the mathematical model describing linear parabolic problems on graph structures.

Mathematical model
Let P = {p j , j = 1, . . . , J} be a set of branching points or vertexes. Some of these points are joined by individual directed edges making a set of edges E (or a set of ordered pairs of vertexes). Each edge has a direction and these directions can be selected by using any convenient scheme: where p ks is a starting point of edge e k , and p k f is an end point of the same edge. The graph is connected, i.e. there is a path from any vertex to any other vertex in the graph, including termination points (in definition of the connectivity, the direction of edges are not taken into account).
Let l k be a length of e k and let x be a distance along the interval (0, l k ), which defines the given edge in a geometrical domain. For each point p j ∈ P we denote by N ± (p j ) the sets of edges, having x = 0 or x = l k as an end point of edge at p j : N + (p j ) = {e k : e k = (p js , p j ) ∈ E, s = 1, . . . , S j }, On the graph (E, P ) we consider a system of parabolic linear problems for functions {u k (x, t)}: Let us divide all branch points into two sets P = T ∪ P 1 . Point p j ∈ P is assigned to a set of termination points T if there exists only one edge e kj , which terminates at p j . Let x = b kj be a coordinate of the end of interval (0, l kj ), which is joined to p j , i.e. b kj = 0 or b kj = l kj . The following boundary conditions are formulated: where µ j are given functions. At the branch points p j ∈ P 1 the fluxes of the solution are conserved [6]: At all vertexes of the graph the continuity constraints are satisfied u m (p j , t) = u k (p j , t), ∀p j ∈ P 1 , e m , e k ∈ N ± (p j ).

The fully implicit approximation
We start from the fully implicit approximation. The governing equations (2.1)-(2.6) are discretized by the Finite Volume Method. On each edge e k , k = 1, . . . , K we define a discrete uniform grid here U k,n j approximates the solution u k (x j , t n ) on edge e k . The backward time finite difference, the forward and backward space finite differences are defined as: We approximate differential equations (2.1) by the following finite difference equations (see, also [1,2,9]): Boundary conditions are defined at the termination points The flux balance equations (2.4) at the branch points are approximated by discrete conservation equations (see [1,2]): At all branch points the continuity constraints are satisfied Discrete problem (2.7)-(2.10) has a unique solution, since for each t n we solve a system of linear equations with a M -matrix. The solution of finite volume scheme (2.7)-(2.10) can be computed efficiently by using a modified factorization algorithm (see [1]).

Predictor-corrector type parallel algorithm
On the basis of the sequential discrete algorithm (2.7)-(2.10), the parallel algorithm is developed using domain decomposition method. It is based on the predictor-corrector splitting technique and it is used in [2] for numerical approximation of 1D linear parabolic problem. Similar techniques also can be applied for problems on graph structures (see [9,11] for related discussions). We also mention a recent paper [12], where the domain decomposition method is used to parallelize a solver for simulation of flows in oil filters.
Graph decomposition step. The load balancing problem should be solved during the implementation of this step. First, it is aimed to guarantee that each processor has about the same number of mesh grid points, since this number defines the computational complexity for all parts of the discrete algorithm.
Due to the stencil of discretization, the computational domains of different processors are overlapping. The information belonging to the overlapped regions should be exchanged among processors. The time of data exchanges is contributing to the additional costs of the parallel algorithm. Thus, a second goal of defining the optimal data mapping is to minimize the overlapping regions. Metis tool [5] is applied to distribute the graph with weighted edges, where the weights are taken equal to the number of mesh points on the given edge.
Let us denote a set of edges, which connect vertexes belonging to different subdomains by D = e s ∈ E, s = 1, . . . , S . Then each grid ω h (s) on these edges is divided into two subgrids by midpoints Predictor step. First, we compute in parallel new values of the solution at the splitting points x js ∈ B. The explicit Euler approximation is used to discretize the differential equation (2.1):

Convergence Analysis
Our goal is to investigate the stability and the accuracy of the discrete predictor-corrector algorithm. We define subgrids Let denote the error functions of the discrete solution as Z s,n js = U s,n js − u s (x js , t n ), s = 1, . . . , S.
By putting them into the finite-difference scheme we get a discrete problem for the error functions: where functions ψ k,n , ψ s,n and ϕ n define the truncation errors. At all branch points the error functions satisfy the continuity conditions Let us assume that the initial data, coefficients and the exact solution of problem (2.1)-(2.6) are smooth enough. By using Taylor's formula the truncation errors of the finite difference scheme (2.11)-(2.12) can be estimated as [2]: |ψ s,n js − ψ s,n js | Cτ, ∀x js ∈ B. Our stability analysis essentially uses results of [2,11]. Let For any V ∈ V h we define the L 2 , L ∞ -norms and energy seminorm: A simple generalization of the well known embedding inequality is valid: for any V ∈ V h we have Lemma 1. Let Z k,n satisfy problem (3.1) -(3.6). Then it holds that where J 1 = |P 1 | and the functional E n is defined as Proof. Multiplying equations (3.1) by 2h k τ ∂tZ k,n j , x k j ∈ w h (k), k = 1, . . . , K and equations (3.4) by 2h s τ ∂tZ s,n js±1 , x s ∈ B and adding the resulting equalities, we get The well-known Green summation formula states: for any mesh functions V, W ∈ V h the following equality is valid Taking into account boundary conditions at the termination vertexes, the fluxes at the remaining p j ∈ P 1 vertexes can be grouped in the following way From (3.10), applying equality 2Z n = Z n + Z n−1 + τ ∂tZ n and equation ( ϕ n (p j ) ∂tZ n (p j ).
Using the Schwarz inequality and the ε-inequality, we have Similarly we get the estimate It remains to estimate the error introduced at the predictor step. The analysis is similar to one presented in [2]. From Applying the Schwarz inequality and the ε-inequality, we get the following estimates: By repeated application, this yields the desired result (3.9). Now, using the stability estimate given in Lemma 1 and the estimates of truncation errors (3.7) we get the following convergence result Theorem 1. Let U be the solution of the finite difference scheme (2.11)-(2.12) and u(x, t) be the solution of the differential problem (2.1)-(2.6). Then It follows from Lemma 1 and (3.15) that finite difference scheme (2.11)-(2.12) is unconditionally A-stable, but the approximation of the scheme in the related norm is only conditional.
It is shown in [2] for a linear 1D parabolic problem that the global error of the discrete predictor-corrector scheme solution depends on the prediction error as O(τ 2 /h). The L ∞ norm is used in the analysis.

Domain Decomposition Method
In this section we compare the obtained above results with error estimates presented in [13]. Here the stability and accuracy of the predictor-corrector type schemes is investigated by using the Domain Decomposition (DD) method, which is based on additive splitting schemes.
For simplicity of notations, here we restrict to 1D parabolic problem (see, also [2]): We note that equation (4.1) includes a linear sink term q(x)u(x). The domainΩ = [0, 1] is divided into two subdomains (a generalization to the case of K subdomains is quite straightforward) In accordance with this domain decomposition, the space grid is decomposed Note, that here we use a global numbering of grid points.
• Predictor step. The explicit Euler approximation is used to compute the solution at the boundary point of two subdomains • Domain decomposition step. Solutions on each subdomain are computed in parallel using the implicit finite difference scheme (4.3), the predicted valueŨ n j1 is used as the interface boundary condition: • Corrector step. Third, using the implicit finite difference scheme (4.3) and taking the solution U n , computed at the second step, we update the value of the solution at the boundary point of subdomains: Following [13], we can write this predictor-corrector scheme as an additive splitting scheme. With the decomposition of the discrete domain ω h we associate a corresponding additive representation of the identity operator I: 2 j=1 χ j = I, χ j 0, j = 1, 2.
Then the following discrete operators are defined: It is easy to check, that the operator A h is self-adjoint (with respect to a standard inner product) and positive definite The operators A h,α do not inherit these basic properties, therefore a construction and analysis of splitting schemes is quite non-trivial task. The predictor-corrector scheme (4.2)-(4.4) can be written as the Douglas scheme (see, [13]) Excluding the intermediate solution U n−1/2 , we get the equivalent factorized scheme The corresponding problem for the error Z n j = U n j − u(x j , t n ) is defined by where ψ n is the approximation error The following convergence estimate is proved in [13] (for more results on stability of the Douglas scheme see [4,10]): where the energy norm is defined as Z A h = (A h Z, Z) 1/2 . We split the approximation error into two parts ψ n j = ψ n 1,j + ψ n 2,j . The first term defines the approximation error of the backward Euler scheme ψ n 1,j = f n j − ∂τ u n j − A h u n j , x j ∈ ω h , the second term defines the splitting error (i.e., the predictor-corrector error): We consider both terms in detail. Note, that in [13] the analysis of the first term is skipped on the basis that it defines a standard error of the classical Euler scheme with weights. Let us consider the operator I + τ A h,1 , it is easy to see that for τ = Ch 2 the estimate I + τ A h,1 ≈ cI is valid. Thus the influence of the diffusion is very weak in this case. As an example we investigate a discrete parabolic problem with a constant diffusion coefficient d(x) = 1 and q(x) = 0 when the exact solution of the differential problem is u = exp(x − t). Then the approximation error is equal to ψ n 1,j = (t + h 2 )u(x j , t n ) + O(τ 2 + h 4 ), x j ∈ ω h . For τ = Ch 2 the energy norm of ψ n 1 can be estimated as Thus due to the need to estimate the approximation error in the strong energy norm the convergence rate is reduced and only the conditional approximation is valid.
In Table 1 values of (I +τ A h,1 ) −1 ψ n 1 A h are presented for different discrete time and space steps. Numerical experiments are done for two relations τ = h and τ = 50h 2 of discrete steps. The experimental convergence rates p are defined as (I + τ A h,1 ) −1 ψ n 1 A h = Ch p . It is proved in [13], that the splitting error can be estimated as (4.10) Table 2. Theoretical analysis of the domain decomposition error ψ n 2 .
Now, in Table 2 we collect all three estimates obtained for the splitting error ψ n 2 , these estimates are applied for different relations on τ and h.
We have solved numerically 1D parabolic problem with a constant diffusion coefficient d(x) = 1 and q(x) = 0. The exact solution of the differential problem is u = tx(1 − x). In this case ψ n 1 = 0, and only domain decomposition error ψ n 2 is included in the analysis. Since, the behavior of the error Z n is the same in all tested norms, we present in Table 3 results only for the energy norm Z n E at t = 1. The results of numerical experiments agree well with the theoretical results presented above.

Conclusions
In this work, we have presented a parallel predictor-corrector type algorithm for solution of linear one-dimensional parabolic problems on graphs. By using the energy estimates it is proved that the predictor-corrector algorithm is unconditionally stable. The relation between space and time steps is still required in order to get the convergence of the discrete solution, since only a conditional approximation is obtained due to the truncation error introduced at the prediction step.
Applying results from [13] and using the equivalence of the predictorcorrector scheme and the Douglas type scheme, it is possible to get new convergence estimates. The asymptotic optimality of different theoretical accuracy estimates is compared with results of computational experiments.
An interesting goal of future work is to compare the accuracy of fully implicit and predictor-corrector splitting schemes (including different strategies for graph decomposition).