Novel Finite Difference Discretization of Interface Boundary Conditions for Stablized Explicit-Implicit Domain Decomposition Methods

Stabilized explicit-implicit domain decomposition is a group of methods for solving time-dependent partial difference equations of the parabolic type on parallel computers. They are efficient, stable, and highly parallel, but suffer from a restriction that the interface boundaries must not intersect inside the domain. Various techniques have been proposed to handle this restriction. In this paper, we present finite difference schemes for discretizing the equation spatially, which is of high simplicity, easy to implement, attains second-order spatial accuracy, and allows interface boundaries to intersect inside the domain.

In parallel implementation of corrected EIDD methods, the correction step is communicationally expensive to be parallelized when the interior boundaries cross into each other inside the domain, e.g. as in Figure 1(a).While for some problems [10], it causes no trouble to partition a domain with no intersecting interior boundaries, in many cases corrected EIDD methods suffer from low accuracy when partitioned into a large number of narrow strip subdomains when a large number of processors is used [34].To address the problem of no crossover for interior boundaries, Shi and Liao [29] introduced in 2006 zigzag interior (ZI) boundaries so that in the implicit correction, spatial discretization does not result in coupling of all grip points on the interior boundaries into one single equation.In 2009, Liao, Shi, and Sun [27] developed composite interior boundaries by replacing the ZI boundaries in [29] with straight-line interior boundaries at locations neighboring intersection points of interior boundaries, leading to improved programming simplicity for the treatment of interior boundaries than the ZI boundaries.Zhu, Yuan, and Du [30] [31] used a different technique to handle the crossover of interior boundaries.Jun and Mai [20] [39] used special treatment for the implicit discretization at points neighboring intersection points while maintaining unconditional stability.The interface boundary treatment introduced by Jun and Mai for their modified implicit prediction method [20] [39] can also be used to solve the intersecting interior boundary problem for corrected EIDD methods.Zhuang and Sun [35] and Wang, Wu, and Zhuang [37] tackled disadvantages of no-crossover interface boundaries by using a data partition different from the domain partition, where the domain is partitioned with no crossover interface boundaries as in Figure 1(b) but the data of each subdoman is further partitioned into multiple data subsets like in Figure 1(a) for distribution to different processors.
To allow crossover of interface boundaries, in this paper we propose new finite difference schemes for interface boundary conditions at intersecting points of crossover interface boundaries.The technical motivation of this finite difference approximation is given in Section 2, which describes a Stabilized Explicit-Implicit Domain Decomposition Methods and the problems of the interface boundary condition treatment when standard finite different approximation is used.Section 3 describes the new finite different schemes for interface boundary conditions.In Section 4 we present numerical tests, and Section 5 gives the concluding remarks.

The Stabilized EIDD Method and Need for Interface Boundary Condition Treatment
Stabilized EIDD (SEIDD) methods [34] [35] and the more general corrected EIDD (CEIDD) methods [15] [27] [28] [29] [30]- [35] are operator splitting time discretization methods for time dependent partial differential equations, where operator splitting is domain decomposition-based.A common feature of the SEIDD and CEIDD methods that call for a communication-efficiency-targeting treatment of the interface boundary conditions is the stabilization or correction of the interface boundary conditions by an implicit time discretization scheme.To see how this technical issue arises, for reading convenience, we give the description of a SEIDD method below.To that end, we first list some notations.To numerically solve problem (1 -2), we choose a discrete spatial grid Ω ℎ with mesh size h, and discretize Equations ( 1) and ( 2) spatially into ( ) ( ) where A is the discrete approximation of the spatial operator on the right hand side of Equation (1).Our description will be based on this spatially discrete form of the equation.For a domain partitioned as in Figure 1(b) or Figure 1(c), let B be the set of grid points on interface boundaries.With k u denoting the numerical solution of the k -th time step, the SEIDD method for computing the solution 1 k u + at the (k + 1)-th time step from the current k-th time step is given below.

A SEIDD Method
1) Compute the interface boundary condition using the explicit forward Euler scheme 1 ( ) + ∆ on inter face boundaries points B , where I is the identity matrix.
2) Using the interface boundary conditions computed at step 1 together with exterior boundary conditions, compute the solution on the subdomains using the implicit backward Euler scheme ( ) 3) Throw away the interface boundary condition computed at step 1.Using solution data 1 k u + on nearby sub- domain as boundary conditions, re-compute interface boundary condition on interface boundary with the backward Euler ( ) A SEIDD method or a CEIDD method uses an implicit scheme, e.g. the backward Euler, to implicitly recompute solution 1 k u + on the interface boundary .When the domain is partitioned with no intersection of interface boundaries, like in

The New Finite Difference Approximation of the Interface Boundary Conditions
With the discussion in Section 2, it is desirable that the domain be partitioned as in Figure 1(a).To handle such partitioned domains, we propose finite difference schemes for approximating interface boundary conditions at intersecting points.Our presentation of the finite difference schemes will be based on two-dimensional problems, i.e.Equation (1) has two independent variables, and Equation (1) has the form x y x x y u t x y a x y u a x y u c x y u d x y u s x y u t For the discrete domain, we assume that uniform mesh size is used with mesh size ℎ.Let ( , ) i j x y be an intersecting point of two interface boundaries as in Figure 2, and let ) i j x y ( , ) ( , ) i j x y be four neighboring grid points on the subdomains.For any grid point ( , ) m n x y , the notation , ( ) when no confusion arises, is used to denote the solution at grid point ( , ) m n x y With these notations, the finite difference schemes for approximating the differential operators on the right hand side of Equation (4) at ( , ) i j x y are given below.
( ) ( ) ( ) If ( , ) i j x y is not an intersecting point of interface boundaries, the conventional central finite difference is use to approximate the right hand of Equation ( 4), i.e., ( ) ( ) , , It is known that central finite difference schemes (10) and (11) have second-order accuracy with ( )

3
O h er- rors.It is also easy to verify using Taylor expansion that schemes (6) and ( 7) have second-order accuracy.Also, comparing schemes ( 6) and (10), one can see that scheme (6) is the average of scheme (10) used at the two grid points 1 ( , ) i j x y + and 1 ( , ) i j x y − , and then the second order accuracy of scheme ( 6) follows from the fact that ( ) ( ) for any smooth function f .Scheme (5) also has second-order accuracy.To see that, one can verify that under the change of variables With Equality (12) under the change of variables (13), one can easily prove that scheme (5) has second-order accuracy.Actually, Scheme (5) is the application of schemes ( 8) and ( 9) along the two new variables w and v , which are actually representing the two diagonal lines y x = and y x = − in the original coordinate space of variables x and y .These discussions show that the new finite difference schemes (5-7) have second-order accuracy, which is stated in the following Theorem.
Theorem.Finite difference schemes (5-7) for approximating the differential operators on the right hand side of Equation ( 4) at intersecting point ( , ) i j x y of interface boundaries have second-order accuracy with errors of order ( ) The theorem above means that finite different Schemes (5 -7) have good accuracy, as good as the conventional schemes (8 -11).Now, let us look at another feature of Schemes (5 -7) in term of its impact on the SEIDD methods.From the formulas of schemes (5-7), one can see that the solution  is used at points ( , ) x y is to approximate the right hand side of Equation (4).Since ( , ) i j x y is an intersecting point of interface boundaries, points x y x y + + are in subdomains.So for the stabilization step of the SEIDD method, i.e.
Step 3 of the SEIDD method, Schemes (5-7) use the solution at points x y x y and 1, 1 i j u + + , components of the solution on the subdomains which have been computed in Step 2. But Schemes (5-7) do not use the solution on any point on the interface boundaries except point ( , ) i j x y .Thus, the discrete equation of the stabilization step of a SEIDD method that involves solution at point ( , ) i j x y does not involve any other grid points on the interface boundaries, leading to the decoupling of discrete equation on an intersecting point of interface boundaries from all other discrete equations on other interface boundaries.Such decoupling enables the efficient parallel processing of the solution process without all-to-all communications for the stabilization step.
In addition to decoupling the discrete equations on interface boundaries, another feature of Schemes (5-7) is that the scheme is very simple to implement into simulation code and does not require sophisticated programming techniques, which is very helpful for simulation code development and hence lower the chance of code bugs due to its code implementation simplicity.

Numerical Experiments
To experimentally examine the performance of finite difference schemes (5-7), we choose two testing problems on the same spatial domain square domain . The two problems are 1.
( ) Uniform spatial grid of 2048 2048 × was chosen so that the x-direction and y-direction mesh size is 2 / 2048 h π = .The simulation time interval is [ ] 0,1 and the total simulation time steps is 4000 with a time step size Δ 1 / 4000 t = .To test the proposed finite difference schemes (5-7), the spatial domain is divided into p p × equal-size square sub-subdomains as in Figure 3, with 2 p ranging from 1 to 256, and each square sub-subdomain is assigned to a different processor.
On intersecting points of interface boundaries, Schemes (5-7) was used to discretize the equation.On other points, including points in subdomains and on interface boundaries, standard finite difference schemes (8)(9)(10)(11) were used.With the spatial discretization, the equation on each subdomain is solved by a modified symmetric successive over-relaxation (QSSOR) tailored to non-symmetric matrices [40] with 30 QSSOR iterations.The equations on the interface boundaries in Step 3 (except equations at the intersecting points) are tridiagonal systems and are solved by a tridiagonal solver.The computation on Step 1 is a matrix-vector multiplication since it uses the forward Euler method.We measured the errors of the computed solution in  ∞ -norm, i.e. the maximum errors at time t = 1, and listed the errors in Table 1 for the two problems for the indicated number of subdomains.
To compare the result of using Schemes (5-7) with the result of using the conventional central finite difference, we also used schemes (8)(9)(10)(11) for discretizing Equation (4).Since using Schemes (8)(9)(10)(11) would couple all grid points on the interface boundaries together into a discrete equation on all interface boundary points, the domain is partitioned as in Figure 1(b) or Figure 1(c) with no intersecting points of interface boundaries.We measured the errors of the computed solution in  ∞ -norm at time t = 1, and listed the errors in Table 2 for the two problems for the indicated number of subdomains.The test data show that when the number of subdomains reaches 64 or larger, the new finite difference scheme produces higher accuracy.

Concluding Remarks
Stabilized explicit-implicit domain decomposition provides efficient, stable, and highly parallel methods for solving time-dependent partial difference equations of the parabolic type on parallel computers, but suffer from a restriction that the interface boundaries must not intersect inside the domain.In this paper, we present a finite difference scheme for discretizing the equation spatially, which is of high simplicity, easy to implement, attains second-order spatial discretization accuracy as does the conventional central finite difference, and allows interface boundaries to intersect inside the domain.equal-size subdomains, where P is the number of processors.The second column under m × n indicates the discrete grid size of each subdomain.When P = 1, there is Scheme 5 is not used since there is no interface boundary.

Figure 1 .
Figure 1.Different ways of domain partitions.

Figure 1 (
b) or Figure 1(c), the implicit re-computation of interface boundary condition on different interface boundaries can be executed by processors independently and hence in parallel when different interface boundaries are assigned to different processors.But when a domain is partitioned as in Figure 1(a), conventional finite difference approximation of Equation (1) on the interface boundaries would generate a discrete equation coupling all grid points on the interface boundaries, and since these interface boundary grid points are distributed on different processors, to solve a discrete equation involving all grid points of the interface boundaries would require expensive all-to-all communication.Domain partition with no crossover interface boundaries like Figure 1(b) or Figure 1(c) will not have this problem but would require the domain be decomposed into many long and narrow subdomains as in Figure 1(c) which has nine subdomains, the same number of subdomains as in Figurte 1(a).

Figure 2 .
Figure 2.An intersection point of interface boundaries.

Figure 3 .
Figure 3. Illustration of the domain partition.

Table 1 .
Maximum errors with the new finite difference scheme.

Table 2 .
Maximum errors with the standard finite difference Note that when P = 1, the method and the data are exactly same as those in Table1since no interface boundaries.