A Parallel Grid Modification and Domain Decomposition Algorithm for Local Phenomena Capturing and Load Balancing.
Lion's nonoverlapping Schwarz domain decomposition method based on a finite difference discretization is applied to problems with fronts or layers. For the purpose of getting accurate approximation of the solution by solving small linear systems, grid refinement is made on subdomains that contain fronts and layers and uniform coarse grids are applied on subdomains in which the solution changes slowly and smoothly. In order to balance loads among different processors, we employ small subdomains with fine grids for rapidly-changing-solution areas, and big subdomains with coarse grids for slowly-changing-solution areas. Numerical implementations in the SPMD mode on an nCUBE2 machine are conducted to show the efficiency and accuracy of the method.
Introduction
Grid refinement methods have been proved to be essential and efficient in solving large-scale problems
with localized phenomena, such as boundary layers or wave fronts. For many engineering
problems, however, this still leads to large-size linear systems of algebraic equations, which can not
be solved easily even on today's largest computing systems.
Domain decomposition methods have been extensively investigated recently since they provide
a mechanism for dividing a large problem into a collection of smaller problems; each of them can
be solved sequentially or in parallel, on a workstation or a processor with moderate computing
capability.
In this work, we investigate the possibility of combining domain decomposition methods and
grid refinement techniques. We first divide the original physical domain into a collection of subdo-
mains, and then apply fine grids in subdomains that contain local phenomena and coarse grids in
This paper appeared in: Journal of Scientific Computing, 12(1997), 99-117.
y Department of Mathematics, Wayne State University, Detroit, Michigan 48202. Email: yang@math.wayne.edu,
dyang@na-net.ornl.gov
subdomains in which the solution changes slowly. For the purpose of load balancing among different
processors, when solved on a parallel computer or a distributed computing system, the size of the
subdomains that contain local phenomena should be kept small, while the subdomains in which the
solution changes smoothly should be large. In our implementation, we try to keep approximately
the same number of degrees of freedom for each subdomain.
Our results are experimental and for one-dimensional problems. Multi-dimensional problems
can be considered similarly but with more complexity. In section 2, we will describe the domain
decomposition method with grid modification, and in section 3, we will conduct some numerical
experiments.
domain decomposition method with local grid refine-
ment
In this section, we describe a domain decomposition method which will be combined with grid
refinement techniques.
2.1 The differential problem
Consider the model problem in 1-D:
@x
@x
@x
(1)
@-
(2)
@-
where - is the outward normal, and - are given constants.
Let the
decomposed into K nonoverlapping subdomains satisfying
Then the system (1)-(3) is equivalent to the following split problem:
@x
@x
@x
@-m
@-
@-
is the outward normal to the boundary of
subdomain\Omega k . In [5], Lions proposed a non-overlapping
Schwarz algorithm which can be heuristically stated as follows. Choose the initial
guess
k satisfying the given boundary condition (2)-(3). For the sequence u n
k is
constructed such that
\Gammaq @u n+1
@-m
@-
where - km is a parameter which can be chosen to speed up the convergence, and L is the given
elliptic operator.
Despr'es [1], Douglas et al [2], and Kim [4] have considered discretizations of this procedure
by the hybridized mixed finite element method and the finite difference method and applied the
procedure to wave equations. Mu and Rice [7] considered collaborating PDE solvers in the context
of domain decompositions.
2.2 The finite difference discretization
Assume that uniform grid with size h k is used in
subdomain\Omega k and that there are d k grid points
in
\Omega k with the i-th grid point being denoted by x k;i , where
the interface point between
and\Omega k+1 is x k;d k (= x k+1;1 1: The
endpoints
of\Omega are x We denote, for example,
In the following, we describe two schemes based on finite difference discretization, which differ
from each other in the way of discretizing the convection term.
Scheme 1:
For
using the exterior point x we have the second
order finite difference scheme and the first order interface condition at x k;1 (see figure 1)
\Gammaq k;1
k\Gamma1;d
k\Gamma1;d
Note that q
only if q is continuous at x k;1 (= x k\Gamma1;d
x
x
Thus (16) represents a general form of continuity of function value and flux.
x k\Gamma1;1
a
x k\Gamma1;2
a a
x k\Gamma1;d
\Omega
x k;1

Figure

1: Grid points in
subdomain\Omega k , and interface between
and\Omega
Eliminating
k;0 from (15) and (16) we have the equation for the first point in
subdomain\Omega k .
For interior points of
subdomain\Omega k , we apply the second order finite difference scheme. For the
last point of
subdomain\Omega k , we have
Similarly, in (17)-(18) the point x k;d k is an imaginary one. The value u n
k;d k +1 in
(17)-(18) will be eliminated, under the assumption that q k;d k
Thus we have the following linear system for
k\Gamma1;d
where for the first
the equation (19) should be changed to
and for the last
subdomain\Omega K , the equation (21) should be changed to
K;d K
Here we assumed that h k has been chosen such that
It is easy to see that when the sequence u n
k;j converges as n !1, the limit, denoted by u k;j , will
satisfy the standard global finite difference scheme for problem (1)-(3), and thus u
O(h) in the L 1 norm, where g.
Scheme 2:
In this scheme we employ the first order finite difference approximation to the convection term
on each interface. To be more specific, applying the first order finite difference to the convection
term in (15), we have
Similarly we replace (17) by
Then the equation (19) should be changed to
k\Gamma1;d
k\Gamma1;d
And the equation (21) should be changed to
All other equations remain unchanged. Thus Scheme 2 consists of equations (27), (20), (28), (22),
and (23). For this scheme, the assumption (24) is not required.
2.3 The matrix form
We now rewrite Scheme 1 into a matrix form which may be helpful when implementing the scheme.
Define the tridiagonal matrices,
. ff k;d k
ff K;d
and define the block tridiagonal matrix
where
Then the system (19)-(21) can be rewritten into
K;d K
Scheme 2 can also be put in the matrix form (39), but with different definitions of G and H.
2.4 Convergence analysis and the relaxation parameter
Let ae(A) denote the spectral radius of a matrix A. The relaxation parameters - k;m will be chosen
such that ae(G Let the elements of D \Gamma1
k be denoted by y k
i;j , that is, denote D \Gamma1
(y
We will show that such - k;m can be obtained by the following relation:
y
if we require that -
Note that the formula (41) is valid for both Scheme 1 and Scheme 2. Following Tang [8], we
give a proof of (41) in case of three subdomains, that is, deleting all the columns except
the eight nonzero ones and corresponding rows of G \Gamma1 H, we see that ae(G
y (2)
y (2)
y (2)
y (2)
Let
Then
y (2)
z
z
y (2)
where
z
z
z
By virtue of (43) we have that ae(G
z
z
In order to let ae(G choose z 0: In view of (35)-(38) for Scheme 1
and similar formulas for Scheme 2, we have
y (1)
y (1)
y (2)
y (2)
This finishes the proof of (41) in case
In order to find the optimal parameters by (41), we need to know the ratio of the two elements
of the inverse of the matrix D k . Assume that - k\Gamma1;k has been found, the ratio can be determined
from the relation
y
y
Indeed, an incomplete (omitting the last step) LU-factorization to D k leads to : D
L has diagonal elements 1, and U has elements U (k)
Now the (d k \Gamma 1)\Gammath row
of (46) gives us that
U
y
In view of (41) we have
U
Thus optimal parameters - k;k+1 can be found recursively from (47) by incomplete LU-decompositions
so that the iterative matrix of the system (39) has spectral radius 0. This also gives the convergence
for both Scheme 1 and Scheme 2. Note that the matrix D k is tridiagonal. Thus the ratio
U
can be obtained by a single (not nested) loop (with less than 3d k floating point
operations), and an incomplete LU decomposition is not needed. In our implementation, each processor
(except the last one) computes the parameter at the right-hand boundary of the subdomain.
For constant coefficients, Tang [8] used (41) to find optimal parameters by directly computing the
inverse of the matrices D k . Note that Tang's analysis in case of minimum subdomain overlapping
is equivalent to the nonoverlapping case presented here. Kim [4] applied the procedure (47) to
wave problems. The analysis here extends Tang and Kim's to nonuniform grids with an explicit
treatment of discontinuous diffusion coefficients.
3 Numerical experiments
In this section, we implement our grid refinement and domain decomposition strategy on the
nCUBE2, a MIMD parallel computer with distributed memory, located in Department of Computer
Science at Purdue. It is observed that Scheme 2 gives a little more accuracy than Scheme 1
and does not have the restriction as described in (24). Thus we just report our results for Scheme
2. It should be noted that Scheme 1 gives second order approximation to the convection term as
well the diffusion term and provides a natural approach to discretizing the differential equation.
3.1 Hypercube machines
We now describe briefly hypercube machines and their properties that have influence on our imple-
mentation. By definition, an r-dimensional hypercube machine has nodes and r2 r\Gamma1 edges.
Each node corresponds to an r-bit binary string, and two nodes are linked with an edge if and
only if their binary strings differ in precisely one bit. See Figure 2. As a consequence, each node
is incident to r = log N other nodes, one for each bit position. In addition to a simple recursive
structure, the hypercube also has low diameter (log N) and high bisection width (N=2).
At each iteration of our algorithm, each node computes an error between the current solution
and the solution at the previous iteration on its subdomain. The maximum among the errors on all
the nodes need be found to decide if the iterative process should be stopped. In our implementation,
we first let each of the nodes whose binary string have the form i its value
to the node whose binary string is then let the node i 1
compute and keep the maximum of the two values. After this operation, we have only half of the
values to be dealt with, stored in nodes By the recursive structure of the
hypercube machine, we need to continue only log N steps to get the maximum value among all the
subdomain errors. The maximum is thus finally found on node 0. If this maximum is less than a
prescribed small number, node 0 then sends a message to every other node to stop. Otherwise the
iterative process continues until, at some later iteration, the maximum error is small.

Figure

2: Examples of interprocessor connection and data network for hypercube machines.
For example, consider a 2-dimensional hypercube, corresponding to Figure 2. At
the first step, nodes 10 and 11 send their values to nodes 00 and 01, respectively. Then node 01
computes and stores the maximum value among nodes 01 and 11, and node 00 computes and stores
the maximum value among nodes 00 and 10. At the second step, node 01 sends its updated value
to node 00, and node 00 computes and stores the maximum. This procedure can be efficiently
implemented using bit masks provided by the C programming language.
3.2 Parallel implementation
We will adopt the following stopping criterion for the iterative procedure:
k1 denotes the discrete L 1 norm. However, subdomain problems are solved by Gauss
elimination with banded storage. In all the computations below, we choose the initial guess to be
zero. When the stopping criterion (48) is satisfied, the procedure is stopped and the relative errors
in the L 1 norm between the iterative solution and the true solution are computed.
To test the accuracy of the method, consider
Example 1:
du
dx
The exact solution to the problem is
where the dimensionless quantity is the Peclet number. When B ?? 1, u(x) exhibits
a boundary layer of thickness O(B 1. In this example, we keep
different values for q.

Table

1: Numbers of iterations and L 1 norms of the errors between the true solution and the domain
decomposition solution for Example 1 with different Peclet numbers and grid sizes, implemented
on processors with subdomains.
Number of iterations 21 20 20 20
of GMDD 3.9E-2 4.7E-2 3.2E-2 3.4E-2
of FDM1 8.7E-2 2.5E-1 4.5E-1 5.8E-1
of FDM2 3.4E-2 3.5E-2 3.5E-2 3.5E-2
The domain is decomposed into 16 subdomains. Since the boundary layer is near the point
apply coarse grid h 1for the left 10 subdomains, fine grid h
is the
Peclet number, for the rightmost 4 subdomains, and ffor the middle two (the 11-th and
12-th, if counting from the left) subdomains. The sizes of the subdomains are arranged such that
each subdomain contains approximately the same number of degrees of freedom. Let H c ; H f ; and
Hm be the size of coarse, fine, and medium grid subdomains, respectively. Then
Solving the system we obtain the size of each subdomain. We denote this Grid Modification and
Domain Decomposition method by GMDD. To compare its accuracy, we also employ the second
order finite difference method without domain decomposition with
1. (denoted by FDM1) uniform grid such that the total number of unknowns is equal to that of
GMDD.
2. (denoted by FDM2) uniform fine grid with . The total number of unknowns for
FDM2 is much bigger than that for GMDD.
The results in Table 1 show that the domain decomposition solution is more accurate than the
global method with the same number of degrees of freedom, and as accurate as the global method
with a much larger number of degrees of freedom. For example, for the our grid
modification and domain decomposition method has only 80 degrees of freedom while the global
finite difference method with fine grid (FDM2 in Table 1) has 600 degrees of freedom. But the
accuracies of the two methods are about the same.
We now consider the following non-selfadjoint problem with variable coefficients.
Example 2:
@x
@x
@x

Table

2: Numerical results in the L 1 norm for Example 2 with different numbers of processors;
each processor takes care of one subdomain.
Number of processors 8
Number of iterations
Errors in L 1 norm 7.21E-2 6.26E-2 6.27E-2 5.19E-2

Table

3: Numerical results for Example 3 with discontinuous diffusion coefficient and different grid
sizes, in the case of 16 subdomains.
Number of iterations 78 79 79
The function f on the right hand side is chosen such that the exact solution to the problem is
The true solution has a boundary layer near the right boundary. The domain is decomposed into
different numbers of subdomains, with the right most three subdomains being discretized with grid
size 1, and the rest of the subdomains having grid size 1. The number of subdomains is equal
to the number of the processors in the implementation. The results for this example with different
numbers of processors are shown in Table 2. From Table 2 we know that the number of iterations
required for the procedure to stop is bigger than for problems with constant coefficients, see Table
1. It is interesting to observe that the accuracy improves as the number of subdomains increases.
This agrees with our computing experience in [9].
Then, we see how well this algorithm performs with discontinuous diffusion coefficients and
consider
Example 3:
@x
@x
@x
where
and the exact solution is chosen as:

Table

4: Numerical results for Example 4 with different numbers of processors; each processor takes
care of one subdomain.
Number of processors 8
Number of iterations 29 55 104 197
Errors in L 1 norm 2.99E-2 2.39E-2 2.37E-2 2.32E-2
We decompose the domain into 16 subdomains. The accuracy and number of iterations with
different grid sizes for Example 3 are shown in Table 3. From Table 3 we see that the number
of iterations and the errors are larger for discontinuous diffusion coefficients than for continuous
coefficients.
Finally, we consider another convection-diffusion problem
Example 4:
@x
@x
@x
The function f on the right hand side is chosen such that the exact solution to the problem is
The true solution has a boundary layer near the right boundary. The domain is decomposed into
different numbers of subdomains, with one third of the subdomains on the right hand side being
discretized with grid size h f = 1, two thirds of the subdomains on the left hand side having grid
size h c = 1, and one middle subdomain between them having grid size (h c )=2. The length
of each subdomain is chosen as in Example 1 such that each subdomain has approximately the
same number of degrees of freedom. The number of subdomains is still equal to the number of the
processors in the implementation. Numerical results for this example are shown in Table 4.
Concluding Remarks
We have implemented a non-overlapping Schwarz domain decomposition method with grid modification
for elliptic problems involving localized phenomena, such as fronts or layers. In order to
capture local phenomena and save computational work, we apply fine grids in subdomains that
contain fronts or layers, and coarse grids in other subdomains. However, when implemented on
parallel computing systems, the subdomain problems must have approximately the same computational
complexity so that loads on different processors or workstations can be balanced. This is
important for synchronization and achievement of a good speedup.
Though our implementation is for one-dimensional and steady problems, the methodology applies
to time-dependent and multi-dimensional problems. Time-dependent problems are usually
first discretized in time to elliptic problems at each time step, and then domain decomposition
algorithms can be applied. However, for advection-dominated transport problems, discretization in
space and time should be coordinated. For example, unwinding techniques can be incorporated in
the discretization process.
Nonoverlapping Schwarz domain decomposition methods have recently received a lot of atten-
tion, since its efficiency and elegant simplicity in implementation, great savings in computer storage
(in 3D, even small overlapping of subdomains could cause a lot more storage), and direct applicability
to transmission problems. In this direction, several other types of domain decomposition have
been considered [6, 11, 10, 12]. Their modifications can apply directly to domain decompositions
with cross points and with long and narrow subdomains. In [3], we successfully implemented a variant
of the methods [12] for selfadjoint and non-selfadjoint elliptic partial differential equations with
variable coefficients and full diffusion tensor, in the case of 100 subdomains with 81 cross points in
2-D. The variant also works well for long and narrow subdomains with length of a subdomain being
times as large as the width of the subdomain. In [9], the author considered a dynamic domain
decomposition method based on finite element discretization for two-dimensional problems. The
dynamic change of domain decompositions can provide a mechanism for capturing moving layers
and fronts, and for load balancing on different processors.

Acknowledgement

: The author would like to thank Professor John R. Rice for his valuable advice
and for providing the computing facility in Department of Computer Science at Purdue. He also
feels grateful to Professor Jim Douglas, Jr. and Dr. S. Kim for helpful discussions. The referee's
suggestions also improved the quality of the paper.



--R

Domain decomposition method and the Helmholz problem


Finite difference domain decomposition procedures for solving scalar waves
On the Schwarz alternating method III: a variant for nonoverlapping subdomains
A relaxation procedure for domain decomposition methods using finite elements
Modeling with collaborating PDE solvers: theory and practice

Different domain decompositions at different times for capturing moving local phenomena
A parallel iterative nonoverlapping domain decomposition procedure for elliptic problems
A parallel iterative domain decomposition algorithm for elliptic problems
A parallel iterative nonoverlapping domain decomposition method for elliptic interface problems
--TR
A relaxation procedure for domain decomposition methods using finite elements
Generalized Schwarz splittings
Different domain decompositions at different times for capturing moving local phenomena

--CTR
Daoqi Yang, Finite elements for elliptic problems with wild coefficients, Computational science, mathematics and software, Purdue University Press, West Lafayette, IN, 2002
J. R. Rice , P. Tsompanopoulou , E. Vavalis, Fine tuning interface relaxation methods for elliptic differential equations, Applied Numerical Mathematics, v.43 n.4, p.459-481, December 2002
finite difference method;grid modification;parallel computing;domain decomposition method