close

Вход

Забыли?

вход по аккаунту

?

321

код для вставкиСкачать
INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
ITERATIVE SOLUTION OF SYSTEMS OF EQUATIONS IN
THE DUAL RECIPROCITY BOUNDARY ELEMENT
METHOD FOR THE DIFFUSION EQUATION
V. BULGAKOV1,s, B. S[ ARLER2,* AND G. KUHN1
1¸ehrstuhl f ür ¹echnische Mechanik, ºniversität Erlangen-Nürnberg, Egerlandstr. 5, D-91058 Erlangen, Germany
2¸aboratory for Fluid Dynamics and ¹hermodynamics, Faculty of Mechanical Engineering, ºniversity of ¸jubljana,
As\ kerc\ eva 6, SI-1000 ¸jubljana, Slovenia
ABSTRACT
In this paper the diffusion equation is solved in two-dimensional geometry by the dual reciprocity boundary
element method (DRBEM). It is structured by fully implicit discretization over time and by weighting with
the fundamental solution of the Laplace equation. The resulting domain integral of the diffusive term is
transformed into two boundary integrals by using Green’s second identity, and the domain integral of the
transience term is converted into a finite series of boundary integrals by using dual reciprocity interpolation
based on scaled augmented thin plate spline global approximation functions. Straight line geometry and
constant field shape functions for boundary discretization are employed. The described procedure results in
systems of equations with fully populated unsymmetric matrices. In the case of solving large problems, the
solution of these systems by direct methods may be very time consuming. The present study investigates the
possibility of using iterative methods for solving these systems of equations. It was demonstrated that
Krylov-type methods like CGS and GMRES with simple Jacobi preconditioning appeared to be efficient
and robust with respect to the problem size and time step magnitude. This paper can be considered as
a logical starting point for research of iterative solutions to DRBEM systems of equations.( 1998
John Wiley & Sons, Ltd.
KEY WORDS: diffusion equation; boundary elements; dual reciprocity; iterative methods
INTRODUCTION
The boundary element method is a weighted residual method for solving partial differential
equations, characterized by choosing an appropriate fundamental solution as a weighting
function and by using the generalized Green’s formula for complete transfer of one or more
partial differential operators on the weighting function. The main comparative advantage of the
*Correspondence to: B. S[ arler, Faculty of Mechanical Engineering, University of Ljubljana, As\ kerc\ eva 6, SI-1000
Ljubljana, Slovenia. E-mail: bozidar.sarler@fs.uni-lj.si
sOn leave from: Department of Applied Mathematics, Moscow State University of Civil Engineering, 129337 Moscow,
Yaroslavskoe Shosse 26, Russia
Contract/grant sponsor: Alexander von Humboldt Foundation
Contract/grant sponsor: Slovenian Ministry of Science and Technology
Contract/grant sponsor: International Bureau of BMBF
CCC 0029—5981/98/040713—20$17.50
( 1998 John Wiley & Sons, Ltd.
Received 11 March 1997
Revised 19 December 1997
714
V. BULGAKOV, B. S[ ARLER AND G. KUHN
boundary element method over other discrete approximative methods is demonstrated in cases
where all of the resultant domain integrals could be represented by the boundary integrals. This
turns out to be possible only for some partial differential equations. When dealing with the
boundary element method for the transport equation structured by weighting with the fundamental solution of the Laplace equation, domain integrals appear at least from the transience,
convective, and source terms.
The dual reciprocity boundary element method (DRBEM)1 represents one of the possibilities
for transforming the resultant domain integrals into a finite series of boundary integrals. The key
point of the DRBEM is approximation of the field in the domain by a set of global approximation
functions and subsequent representation of the domain integrals of these global approximation
functions by the boundary integrals.
The method was first applied to simple diffusion-governed linear2 and non-linear problems,3
then to steady convective—diffusive problems,4 and finally successfully used in a completely
non-linear transient convective—diffusive context including phase-change effects and non-linearities arising simultaneously from material properties and boundary conditions.5
Despite the proved practical applicability of the method, important theoretical questions
remain. For example, which global interpolation functions are suitable for representation of the
fields in the domain, and where should the collocation points of these global interpolation
functions be put. Up to present, both questions have been addressed mostly from the heuristic
point of view. The most popular global interpolation function used in the majority of the
DRBEM calculations so far in the field of transport phenomena is 1#r , with r representing the
n
n
Euclidean distance between field point p and collocation point p . The convergence of these
n
functions was studied both numerically and in a more formal way in Reference 6. Up to now, only
ad hoc distributed collocation points have been used, since an appropriate error estimator does
not exist. A typical example of a completely non-uniform DRBEM collocation point mesh in two
dimensions as well as in axisymmetry can be found in Reference 7.
Some mathematically substantiated answers to the first question have been recently rediscovered in works concerning the theory of radial basis functions.8 It has been demonstrated9
that the use of thin plate spline radial basis functions gives an approximation which minimizes its
curvature. The first theoretical investigations of the convergence of these functions were carried
out in Reference 10.
The global interpolation of the fields over a domain introduces unknowns in the domain in
addition to the unknowns at the boundary. The number of these additional unknowns usually
exceeds the number of boundary unknowns. The systems of algebraic equations resulting from
DRBEM are thus large, fully populated and unsymmetric. A third important, but not at all
sufficiently investigated element of the DRBEM method, is the possible iterative solution of
associated systems of algebraic equations. This issue is of utmost importance when solving
large-scale problems.
In the last decade a considerable amount of effort has been put into the application of iterative
methods to solve ‘standard’ (without interior points) BEM systems of equations. We quote here
some papers devoted to this question. In References 11 and 12 iterative methods of Krylov type
were applied to BEM systems. In Reference 13 an attempt was made to use the hierarchical-basisbased approximation and preconditioner in 2D BEM similar to that of finite elements. In
Reference 14 a new aggregation-based transformation method was proposed to construct efficient
sparse preconditioners for an iterative solution of BEM. In Reference 15 iterative solution
procedures with ILU preconditioning derived from substructuring were applied to Navier—Stokes
( 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
DUAL RECIPROCITY BOUNDARY ELEMENT METHOD
715
equations. In Reference 16 a transformation technique was considered which leads to diagonally
dominant BEM matrices. Reference 17 is devoted to extensive investigation of Krylov iterative
techniques applied to elasticity and heat conduction problems. In Reference 18 the question of
iterative treatment of BEM equations for practical applications and some views on the theoretical
aspects of using Krylov subspace methods in practical situations are considered. Reference 19 is
devoted to the application of rapidly developing Wavelet transform methods to 2-D BEM
problems.
On the other hand, publications which consider the iterative solution of DRBEM systems of
equations are not known to the authors.
Regarding the potential use of the dual reciprocity method in the context of large-scale
industrial transport phenomena calculations, for example solidification simulations of semicontinuously casted aluminium ingots,7 the lack of research on the iterative solution of respective
systems of algebraic equations together with the lack of error estimators present two principal
points which definitely require a substantial amount of additional research. These two points are
recognized as major limiting factors in the state of the art of the DRBEM.
In this study the diffusion equation in two dimensions is considered as a logical starting point
of the research of iterative solutions of DRBEM systems of equations. Its final goal would be an
iterative solution of the systems of algebraic equations arising in the DRBEM solution of the
complete non-linear transport equation.
GOVERNING EQUATIONS
Consider a connected fixed domain ) with boundary !. We seek the solution of the dimensionless
diffusion equation
L
¹"+2¹
Lt
(1)
subject to the initial condition
¹(p, t )"¹ ,
0
0
p3)
(2)
p3!D, t )t)t #*t
0
0
(3)
the Dirichlet boundary condition
¹(p, t)"¹ ,
!
and the Neuman boundary condition
L¹
(p, t)"F ; p3!N, t )t)t #*t
(4)
!
0
0
Ln
!
where p is a position vector in a two-dimensional Cartesian frame, t time, t initial time, *t
0
positive time increment, and n normal vector on the boundary !. !D and !N stand for not
!
necessarily connected parts of the boundary ! where the Dirichlet and Neuman boundary
conditions are prescribed, respectively
!"!DX!N
(5)
¹ , ¹ , and F are known functions. The solution of equation (1) according to the initial
0 !
!
conditions (2) and boundary conditions (3), (4) is the field ¹(p, t #*t); p3)X!N .
0
( 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
716
V. BULGAKOV, B. S[ ARLER AND G. KUHN
DISCRETIZATION
Discretization over time
The governing equation is discretized with respect to time by using the fully implicit scheme
¹j`1!¹j
"+2¹j`1
*t
(6)
where index j#1 denotes the time level t #*t and index j stands for the time level t .
0
0
Discretization over space
By weighting the time-discretized governing equation by the fundamental solution ¹* of the
Laplace equation over the domain ) one arrives at the following domain integral equation
P)
P)+2¹j`1¹* d)
¹j`1!¹j
¹* d)"
*t
(7)
Equation (7) has to be properly transformed and discretized in order to obtain the DRBEM
solution procedure. The following two integral types have to be considered accordingly;
P)F¹* d)
(8)
P)+2F¹* d)
(9)
with F standing for the scalar-valued function of position. The integral type (8) is representative
for the left-hand side of equation (7), and the integral type (9) is representative for the integral on
the right-hand side of equation (7). We focus on two-dimensional situations, e.g.
1
r
¹*" log 0 ,
2n
r
r2"r ) r,
r"p!s
(10)
where s stands for the source point position and r for the reference radius. The integral (9) could
0
be straightforwardly transformed from the domain type into the boundary type by using Green’s
second identity
LF
L¹*
¹* d!! F
d!!c*F(s)
(11)
Ln
Ln
)
! !
!
!
The simplest discretization is used to discretize the two boundary integrals on the right-hand side
of equation (11). Boundary geometry is approximated by N straight line segments, and spatial
!
variation of the fields on each of the boundary segments is represented by constant interpolation
functions with gridpoints coinciding with the geometrical centres of the straight line segments.
The Einstein summation convention is used in this text, i.e. any index which is repeated twice in
a product is summed up. Underlined indices are not summed up. This gives
P
P
+2F¹* d)"
LF
L¹*
P
LF
P! Ln! ¹*l d!!P!F Ln!l d!!c*l6 Fl+Glkdki Ln!i !HlkdkiFi!c*l6 dli Fi
( 1998 John Wiley & Sons, Ltd.
(12)
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
DUAL RECIPROCITY BOUNDARY ELEMENT METHOD
717
where k"1, 2, . . . , N and i, l"1, 2, . . . , N. N"N #N is the total number of points in
!
!
)
which the solution is calculated. The first N points p coincide with the boundary gridpoints and
!
n
the last N points are arbitrarily distributed over the domain ). Index l stands for s "p .
)
l
l
d represents the Kronecker symbol. Matrix elements G and H are defined as follows:
lk
lk
L¹*
l d!
G " ¹* d!, H "
(13)
l
lk
lk
!
! Ln
k
k
!
where ! represents the kth boundary segment, and c* is equal to
l
k
c*"0, s N)X!; c*"1 , s 3!; c*"1, s 3)
(14)
2
l
l
l
l
l
l
The integral (8) is transformed by considering the approximation of the spatial variation of the
fields in ) by the global interpolation functions of the form
P
P
F(p)+t (p) 1 , u"1, 2, . . . , N#3
(15)
u
u
Two-dimensional scaled augmented thin plate splines are used in this work, due to the theoretical
considerations described in Reference 9:
t (p)"r2 log r , n"1, 2, . . . , N
n
n
n
t
(p)"p !p0
x
N`1
x
t
(p)"p !p0
N`2
y
y
t
(p)"1
N`3
(16)
r2"(p!p ) ) (p!p ), p "p ) i , p "p ) i
x
x
y
y
n
n6
n6
where i and i represent the Cartesian base vectors. Scaling constants p0 and p0 stand for the
y
x
x
y
mean co-ordinates of the domain ). Coefficients 1 are calculated by constructing a system of
u
N#3 algebraic equations
W1"F
(17)
The vectors are 1"(1 , 1 , . . . , 1
)T and F"(F , F , . . . , F , 0, 0, 0)T. The first
1 2
N`3
1 2
N
v"1, 2, . . . , N rows of matrix W are of the form (t , t , . . . , t
), and the last three rows
v1 v2
vN`3
v"N#1, N#2, N#3 are of the form (t , t , 2 , t , 0, 0, 0), where the notation has been
1v 2v
Nv
shortened to F ,F(p ), t ,t (p ). Coefficients 1 follows by inverting the system (17)
n
n nu
u n
u
1"W~1F
(18)
Consequently, integral (8) can be written in a compact dual reciprocity form (k"1, 2, . . . , N ,
!
i, l"1, 2, . . . , N, u"1, 2, . . . , N#3)
P)F¹*l d)+P)tu(~1ui Fi¹*l d)"(lu(~1ui Fi
(19)
with
P)tu¹*l d)
( ,
lu
( 1998 John Wiley & Sons, Ltd.
(20)
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
718
V. BULGAKOV, B. S[ ARLER AND G. KUHN
The integral ( is calculated by defining the harmonic functions tK
lu
u
+2tK (p)"t (p)
u
u
which allows us to write (see equation (11))
(21)
LtK
L¹*
u ¹* d!! tK
d!!c*tK (s)
(22)
u
u
Ln
Ln
! !
!
!
Upper integrals are numerically evaluated by using the same discretization strategy that leads to
expression (12)
P
( "
lu
P
LtK
iu !Ht d tK !c*d tK
(23)
+2tK ¹* d) + Gt d
lk ki Ln
l6 li iu
u l
lk ki iu
)
!
However, matrices Gt and Ht could differ from matrices G and H because the boundary
subdivision that corresponds to calculation of the fields on the boundary and integrals ( could
lu
differ in general. Let us denote the number of boundary segments leading to calculation of
matrices Gt and Ht with N t . Therefore, index k in equation (23) runs as k"1, 2, . . . , Nt .
!
!
The adjacent harmonic functions tK for the thin plate splines (16) are
u
P
tK " 1 r4 log r ! 1 r4
n 16 n
n 32 n
tK
"1 (p !p0)3
x
N`1 6 x
tK
"1 (p !p0)3
N`2 6 y
y
tK
"1 (p !p0)2#1 (p !p0)2
(24)
4 y
N`3 4 x
x
y
After transforming and discretizing the domain integrals in equation (7) as described, the
following compact form is obtained:
with the coefficients
L¹j`1
C¹j`1¹j`1#CFj`1
"C¹j¹j
i
li
li
li i
Ln
!i
(25)
1
C¹j`1" ( (~1#H d #c*d
lk ki
l6 li
li
*t lu ui
1
C¹j" ( (~1
li
*t lu ui
CFj`1"!G d
(26)
li
lk ki
The terms in the upper coefficients that involve ( came from the dual reciprocity transformation
lu
of the transience term. Others correspond to the diffusive term.
Application of the boundary conditions (3), (4) yields the following N]N system of algebraic
equations:
Ax"b
( 1998 John Wiley & Sons, Ltd.
(27)
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
DUAL RECIPROCITY BOUNDARY ELEMENT METHOD
719
with the system matrix A in the form
The vector of unknowns x is
A "f C¹j`1#(1!f ) CFj`1
i6 li
li
i6 li
L¹j`1
, i"1, 2, . . . , N
x "f ¹j`1#(1!f )
!
i
i6 i
i6 Ln
!i
x "¹j`1, i"N
,. . .,N
i
!`1
i
and the adjacent right-hand side vector b is
(28)
(29)
!f CFj`1 Fj`1
(30)
b "C¹j¹j!(1!f ) C¹j`1¹j`1
!
!
i
li i
i6
li
i6 li
l
i
The boundary conditions indicator f is defined 0 for p 3!DX) and equal to 1 for p 3!N.
i
i
i
ITERATIVE METHODS
Now we concentrate on solving equation systems (27) arising in each time step of the described
DRBEM procedure. They are non-symmetric and fully populated and therefore special care in
applying iterative methods is necessary. Preconditioned CGS (Conjugate Gradient Squared)20,21
and GMRES (General Minimal Residual)22 algorithms have already proved to be efficient when
applied to classical direct BEM problems.11,17,18 A short description of these methods is given in
the following subsection, followed by an evaluation of the theoretical and numerical observations
of their behaviour in the DRBEM context.
Preconditioned CGS and GMRES algorithms
A preconditioned CGS algorithm was implemented in the following schematic form with
B standing for the left preconditioner:
SE¹ r0"b, u0"g~1"h0"0, b "0
0
¸OOP k"0, 1, 2, . . .
x"B~1rk
o "(r0)Tx
k
if k'0 b "o /o
k
k k~1
x"x#b hk
k
gk"x#b (b gk~1#hk)
k k
z"B~1Agk
p "(r0)Tz
k
a"o /p
k k
hk`1"x!a z
k
uk`1"uk#a (x#hk`1)
k
rk`1"rk!a A(x#hk`1)
k
¹est for convergence:
[(rk`1)T rk`1]1@2/[(r0)T r0]1@2)e
END ¸OOP
( 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
720
V. BULGAKOV, B. S[ ARLER AND G. KUHN
We use the left-preconditioned GMRES algorithm. In spite of the fact that our numerical
examples exhibited a very fast convergence and therefore no restart in the GMRES method was
required, we give a description of a restarted version of GMRES denoted by GMRES(m) which
may be important for other dual reciprocity applications. Compared to the algorithm given in
Reference 22 our algorithm was modified in order to take into account preconditioning and
a suitable process termination:
1. S¹AR¹: Choose x , Compute r "B~1(b!Ax ), b"Er E, v "r /b
0
0
0
0
1
0
2. Iterations:
¸OOP j"1, 2, . . . , m Do:
Compute w"B~1Av
j
¸OOP i"1, . . . , j Do:
h "(w, v )
i,j
i
w"w!h v
i,j i
END ¸OOP
Compute h
"EwE and v "w/h
j`1,j
j`1
j`1,j
Add to QR decomposition of H
j
Check for ¹ermination/Restart: (see below)
If satisfied consider j as m, go to 3
END ¸OOP
3. Define V "[v , . . . , v ]
m
1
m
4. Compute y "argmin Ebe !H yE
m
y 1
m
5. Compute x "x #V y
m
0
m m
6. If EAx !bE/EbE(e then S¹OP
m
else x "x go to 1
0
m
The QR factorization of H needed for defining y from the least squares problem
j
min Ebe !H yE is updated progressively. It is shown in Reference 22 that this function is equal
y 1
k
to Eg !R yE, where g "Q be , and that the residual Euclidean norm is equal to the (k#1)th
k
k
k
k 1
component of g . This is very convenient since x need not be calculated for obtaining the residual
k
k
norm necessary for the convergence check. However, in our case this is the norm of the residual
vector of the preconditioned system. Therefore, if this norm satisfies a given tolerance we go to
step 3 of the algorithm in order to calculate the true residual norm. If the prescribed accuracy is
not achieved we divide the tolerance for the preconditioned residual vector (which originally
equals e) by 10 and go to the first step. Usually after the next termination upon achieving accuracy
in step 2, global convergence is achieved.
The restart can also occur after performing m iterations of step 2 without the above mentioned
termination.
Theoretical observations
It is known that both the CGS and GMRES methods terminate at most when the number of
iterations becomes equal to the problem size in case of an absence of roundoff errors. This is of
course not acceptable in real computations. In order to accelerate convergence we have to require
that the matrix A has some special properties.
( 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
DUAL RECIPROCITY BOUNDARY ELEMENT METHOD
721
As was mentioned before, the unknown values of the discrete problem are split into two parts,
numbered in such a way that the first N values correspond to the boundary grid points and the
!
last N values correspond to the domain grid points. Therefore, the whole matrix A can be split
)
into four blocks
C
D
A
A
!!
!)
(31)
A
A
)!
))
According to our DRBEM discretization scheme, the block A is related to the standard BE
!!
unknowns and the remaining blocks result from the dual reciprocity method. A lack of theory for
predicting the convergence of iterative methods for standard BE systems, namely with matrices
like A , is doubled for the whole matrix A. We therefore investigate the properties of our
!!
matrices on the basis of numerical experiments.
For symmetric positive-definite problems, when
A"
AT"A, (Au, u)'0,
∀u3Rn, uO0
(32)
a spectral condition number j (A)/j (A) is often a useful characteristic for the performance of
.!9
.*/
iterative methods. Being equal to the ratio of maximal and minimal matrix eigenvalues this
characteristic is not available for nonsymmetric matrices in which eigenvalues can be negative
and complex. However, the distribution of matrix eigenvalues is still very informative and useful.
It is very important for the convergence of CGS and GMRES methods that the matrix is
positive definite
(Au, u)'0,
∀u3Rn, uO0,
(33)
or ‘nearly’ positive definite. From the decomposition of the matrix into Hermitian and SkewHermitian parts (or symmetric and skew-symmetric for real matrices)
A"H#iS H"1/2 (A#AH) and S"1/2i(A!AH)
(34)
and the fact that our matrices are real it follows that
(Au, u)"(Hu, u)
(35)
for real u. It is known that in this case (Au, u)*j (H) (u, u) and therefore our matrix is positive
.*/
definite when its symmetric (Hermitian) part is positive definite. For the GMRES (m) method, e.g.
with positive-definite matrices the following estimate22 exists for the residual vectors
ErmE)(1!a/p)m@2Er0E
(36)
where a"(j (H))2 and p"j (ATA) which proves the convergence of GMRES(m). Taking
.*/
.!9
into account that
j (H))Re(j(A)))j (H)
(37)
.*/
.!9
we can state that if H is positive definite then the real parts of the eigenvalues of matrix A are
positive. The reverse statement is unfortunately not true. However, it is still very important that
all or almost all Re (j(A)) are positive. In this case the CGS method gives the contraction for the
residual vector rk"b!Axk by means of its matrix polynomial / (A)20
k
(38)
rk"/2(A) r0
k
( 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
722
V. BULGAKOV, B. S[ ARLER AND G. KUHN
which is calculated recursively. This contraction usually has asymptotic character, i.e. after some
m iterations one can observe that
ErmE(Er0E
(39)
Another known aspect important for iterative methods convergence is eigenvalue clustering.
‘Good clustering’ means that all eigenvalues fall into a small number of size classes. In this case it
can be expected that the CGS method needs a small number of search directions in order to
suppress the residual vector and the GMRES method needs a small number of vectors of the
Krylov subspace to construct a representative basis. For a matrix scaled row-wise by its diagonal,
which corresponds to the left diagonal preconditioning D~1A applied in the subsequent study it
follows trivially from the Gershgorin theorem23 that
&i D j(D~1A)!1D)o
i
n
where o " + Da /a D
i j"1; jOi ij ii
(40)
This means that for any eigenvalue of D~1A there is a matrix row for which o is the radius of the
i
circle on the complex plane centred at 1 where this eigenvalue is located. It is clear that when the
matrix has a strong diagonal predominance (as an extreme case we can take a unity matrix) then
it is well clustered around 1. When only some single off-diagonal terms are large and the other
two are small compared to 1 then in agreement with (40) we can expect that the majority of
eigenvalues will be still well clustered around 1. In order to get some confirmation of this fact
consider three 4]4 matrices
1
0
0
0
0)1
1
0)1
0
0
0)2
1
0)2
0
0
0)1
1
,
1
100
0
0
0)1
1
0)1
0
0
0)2
1
0)2
0
0
0)1
1
,
1
!100
0
0
0)1
1
0)1
0
0
0)2
1
0)2
0
0
0)1
1
(41)
The first one has a good diagonal predominance. Its eigenvalues are
0)80,
0)99,
1)0,
1)20
They are well clustered around 1. We change this matrix by setting a large off-diagonal entry
equal to 100. The eigenvalues of this matrix are
!2)16,
0)85,
1)14,
4)16
We can see that only two eigenvalues are strongly influenced by this matrix change while the
other two are still well clustered near 1. The third matrix is changed in such a way that we insert
!100 instead of 100. It has the following spectrum:
1#i]3)16,
1!i ] 3)16,
0)85,
1)14
which now has two complex eigenvalues and two other equal to those of the previous matrix. One
can observe that the real part of the complex eigenvalues as well as the real eigenvalues are again
well clustered around 1.
( 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
DUAL RECIPROCITY BOUNDARY ELEMENT METHOD
723
NUMERICAL EXPERIMENTS
Numerical implementation
The DRBEM as well as the equation solvers have been implemented as FORTRAN 77 codes
using double precision accuracy. The regular integrals in the DRBEM method have been
calculated using the eight-point Gaussian quadrature, the weakly singular integrals analytically,
and the strongly singular integrals by the thermal equilibrium technique.
The problems reported herein were solved on a SUN Sparc SUNW Ultra-1, SPARC-SUN4U
workstation. In preconditioned CGS and GMRES methods the diagonal (Jacobi) preconditioners have been used. Calculation of matrix eigenvalues as well as the LU matrix decomposition
were performed by means of programs from the FORTRAN-77 standard subroutine library.24
Calculation of 4]4 matrix eigenvalues, matrix surface plots and eigenvalue distribution graphs
have been performed using the MAPLE package.25 Convergence tolerances for both CGS and
GMRES methods used in the numerical study corresponded to the 106 reduction in the residual
vector norm, namely
ErkE/EbE"EAxk!bE/EbE)e"10~6
(42)
¹est cases. To demonstrate the equation solver performance two geometrically different
arrangements were considered. The domain ) in the first arrangement represents a square region
x~(p (x`, y~(p (y`, y~"!0)5, y`"0)5, x~"!0)5, x`"0)5. The initial condix
y
tion at t"t is ¹ "1. The boundaries x"x~ and y"y~ are subject to the Dirichlet boundary
0
0
condition ¹ "0; t*t and the boundaries x"x` and y"y` are subject to the Neumann
!
0
boundary condition F "0; t*t . The described Dirichlet jump problem is solved with five
!
0
different regular meshes differing in space discretization and timestep size. The boundary element
mesh for calculating the integrals of the global interpolation functions (equation (23)) is made by
subdivision of the boundary elements for calculation of the fields into 11 equidistant subelements, i.e. Nt"11 N . The mesh characteristics of the first mesh arrangement are summarized
!
!
in Table I (1)1—1)5) and the schematic of mesh 1)1 is depicted in Figure 1(a).
The domain ) is in the second arrangement represented by the four-pointed star-shaped region
depicted in Figure 1(b). The arrangement of the mesh used in test case 2)1 is similar to the mesh
from Figure 1(b) (with N "40, N "21) except that it is 20 times more dense. The dense mesh
!
)
characteristics are in Table I (2)1) (In the star example, N( has been set equal to N ).
!
!
Table I. Discretization in test cases
N
!
Nt
!
N
)
N
*t
1)1
1)2
1)3
1)4
1)5
40
80
80
80
120
440
880
880
880
1320
81
361
361
361
841
121
441
441
441
961
0)010
0)100
0)010
0)001
0)010
2)1
800
800
401
1201
0)010
Case
( 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
724
V. BULGAKOV, B. S[ ARLER AND G. KUHN
The Dirichlet boundary conditions are applied at the star boundaries in the first quadrant, and
the Neumann boundary conditions apply on the remaining boundary. The same magnitude of
the initial (¹ ) and boundary conditions forcing functions (¹ ) and (F ) have been used for the
0
!
!
star case as for the square case.
Discussion. Figures 2—7 correspond to the rough mesh 1)1 from Table I. The small example
with mesh 1)1 allows an appropriate graphical representation of the matrix surface plots. The
solution of the algebraic equation system is considered at the first timestep and the results are
summarized in Table II. Figures 2 and 3 demonstrate the surface plot of matrices A and D~1A,
respectively. From Figure 2 it can be observed that the block A , corresponding to the internal
))
points, has a strong diagonal dominance. However, block A has small entries on the main
!!
diagonal which after scaling by the left diagonal preconditioning produce some large off-diagonal
entries. Nevertheless, the number of these large entries is sufficiently small and the majority of
Figure 1(a). Schematic of the dual reciprocity mesh with respect to the first geometrical arrangement. Boundary elements
are denoted with f, domain points with C and the borders between boundary elements with ]
( 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
DUAL RECIPROCITY BOUNDARY ELEMENT METHOD
725
Figure 1(b). Schematic of the dual reciprocity mesh with respect to the second geometrical arrangement. Same symbols as
in Figure 1(a)
matrix coefficients are small compared to 1. According to our theoretical considerations we can
expect that the scaled matrix maintains good clustering around 1. Figures 4 and 5 demonstrate
the distribution of the real part of the eigenvalues of matrices A and D~1A respectively. One can
observe that D~1 scaling improves the matrix in such a way that its eigenvalues are well clustered
and positive. Both these properties provide a fast convergence for iterative methods which is
clearly seen from Table II. Test provide a fast convergence for iterative methods which is clearly
seen from Table II. Test cases 1)2—1)5 demonstrate the same problem with finer space discretization and finer meshes.
One can observe from Table II that by reducing the timestep by a factor of 10 and then again
by a further factor of 10 causes only a small deterioration of convergence (cases 1)2—1)4). The
gradient of the solution in the vicinity of the boundary approaches infinity as the timestep tends
to zero. This explains why the number of iterations increases with the reduction of the timestep.
( 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
726
V. BULGAKOV, B. S[ ARLER AND G. KUHN
Figure 2. Original matrix surface plot for the problem from test case 1)1
Figure 3. Scaled (preconditioned) matrix surface plot for the problem from test case 1)1
( 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
DUAL RECIPROCITY BOUNDARY ELEMENT METHOD
727
Figure 4. Distribution of real parts of eigenvalues of the original matrix from test case 1)1
Figure 5. Distribution of real parts of eigenvalues of the scaled matrix from the test case 1)1
( 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
728
V. BULGAKOV, B. S[ ARLER AND G. KUHN
Figure 6. Distribution of real parts of eigenvalues of the original matrix from the test case 1)5
Figure 7. Distribution of real parts of eigenvalues of the scaled matrix from the test case 1)5
( 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
729
DUAL RECIPROCITY BOUNDARY ELEMENT METHOD
One can observe from Table II that the CGS method appears to be more sensitive to the mesh
size than the GMRES method (cases 1)1, 1)3, 1)5). The role of the preconditioner is more
important for the CGS method than for the GMRES method.
Table II shows that in spite of the fact that the problem size has increased considerably in cases
1)2—1)5 compared to case 1)1, the convergence rate is on the same order of magnitude. The
explanation of this fact is straightforward if one resorts to the eigenvalue distribution depicted in
Figures 6 and 7 for the original and preconditioned matrices for case 1)5. After preconditioning all
the real parts of the eigenvalues are positive and well clustered. In addition, most of them are at
the same clustering level associated with 1. The distribution of eigenvalues in case 1)5 is similar to
that of case 1)1.
For the problem with 961 unknowns (case 1)5) we compare the time required for the direct and
iterative solves. For this comparison the programs LUDCMP for the LU decomposition and
LUBKSB for the back substitution which are available in Reference 24 are used for direct
method. The preconditioned CGS and GMRES methods in our equation solver are used for
iterative method. A speed-up by a factor of 19 relative to the direct method has been achieved
with the CGS method and by a factor of 37 with the GMRES method as seen in Table III.
Further numerical investigation was carried out to study the performance of iterative methods
in the context of the time evolution process for both the square domain problem with 441
unknowns and the star domain problem with 1201 unknowns. Tables IV and V demonstrate the
convergence of preconditioned GMRES and CGS methods at successive time steps 1—3 and
10—12. As an initial guess vector we take both a conventional zero vector and one which
corresponds to the solution obtained from the previous time step. The later is taken from the
initial conditions in the first time step.
As expected we obtained some benefit from using solutions from the previous time step as
initial guess vectors and this can be recommended in practical applications. In this case the
starting and subsequent Krylov subspace vectors calculated progressively in both the CGS and
Table II. Number of iterations of CGS and GMRES methods with and without preconditioning for test
cases 1)1—1)5. The number of boundary and internal grid points as well as achieved accuracy are shown in
parentheses
Iterative method
CGS
Case
GMRES
Pr.Size
Timestep
Precond.
No precond.
Precond.
No precond.
1)1
121
(40#81)
0)01
10
(2)7E!7)
75
(6)6E!7)
19
(1)0E!7)
33
(7)1E!7)
1)2
441
(80#361)
0)1
21
(3)1E!7)
27
(5)0E!7)
20
(1)9E!7)
30
(7)5E!7)
1)3
441
0)01
24
(8)7E!7)
110
(3)3E!8)
23
(2)5E!7)
43
(3)8E!7)
1)4
441
0)001
44
(2)9E!8)
No convergence
—
28
(5)5E!7)
60
(9)6E!7)
1)5
961
(120#841)
0)01
39
(3)0E!7)
116
(7)1E!7)
24
(1)2E!7)
44
(7)6E!7)
( 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
730
V. BULGAKOV, B. S[ ARLER AND G. KUHN
Table III. Comparison of time requirement for solving the system of
equations in test case 1)5 with 961 grid points
Direct LU
Precond. CGS
Precond. GMRES
112
6
3
Time (s)
Table IV. The number of iterations of preconditioned CGS and GMRES
methods required for solving equation systems for time steps 1, 2, 3 and 10, 11,
12 with zero initial guess (ZG) and that provided from the previous time step
(PG). Case 1)3
Time
CGS (ZG)
CGS (PG)
GMRES (ZG)
GMRES (PG)
0)01
0)02
0)03
33
12
9
30
10
10
23
22
21
23
22
21
0)10
0)11
0)12
9
9
9
8
7
6
21
21
21
9
9
14
Table V. The number of iterations of preconditioned CGS and GMRES
methods required for solving equation systems for time steps 1, 2, 3 and 10, 11,
12 with zero initial guess (ZG) and that provided from the previous time step
(PG). Case 2)1
Time
CGS (ZG)
CGS (PG)
GMRES (ZG)
GMRES (PG)
0)01
0)02
0)03
83
67
68
64
65
41
48
46
46
48
37
43
0)10
0)11
0)12
47
68
57
47
37
55
48
47
47
41
41
41
GMRES methods inherit more information about the problem, which is expected to give faster
convergence. It is also important to note that in spite of the complex geometry of the star-shaped
example, the closeness of grid points at sharp corners, same discretization for calculating matrices
G, H and Gt, Ht, and the larger size of the problem, the rate of convergence remains fast and
stable with the increasing time and does not increase substantially compared to square domain
examples.
CONCLUDING REMARKS
The systems of equations produced in each time step of the fully implicit scheme of the boundary
element dual reciprocity method for the diffusion equation with scaled augmented thin plate
( 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
DUAL RECIPROCITY BOUNDARY ELEMENT METHOD
731
spline global interpolation functions were investigated. It was demonstrated that the matrices of
these systems have properties which allow the efficient use of Krylov iterative solvers like the
CGS and GMRES methods with Jacobi preconditioning. The GMRES method appeared to be
more efficient than the CGS method. The same situation is usually observed when solving BEM
problems with only boundary discretization (cf. Reference 18). A simple Jacobi preconditioner
transforms the matrices in such a way that all the real parts of the eigenvalues become positive
and better clustered. Comparison of convergence for the original and preconditioned matrices
clearly demonstrates the positive influence of these two properties on the rate of convergence. The
efficiency of the iterative solution does not drop when the size of the problem is increased. The
iterative solves appeared to be robust with respect to a wide range of time-stepping variations and
clearly far superior to direct methods.
It is also important to speed up the procedure of solving system (17). This will be the subject of
further investigation.
The present study should be considered as a first step in the iterative solving of DRBEM
algebraic equation systems, logically focused on the diffusion equation sub-problem only. The
upgrades of the present work will first consider the steady-state convective—diffusive problems,
and subsequently the full non-linear scalar transport equation.
ACKNOWLEDGEMENTS
The work of the first author has been sponsored by the Alexander von Humboldt Foundation.
This support is gratefully acknowledged. The second and third authors would like to acknowledge the Slovenian Ministry of Science and Technology and the International Bureau of BMBF,
Germany for support in the framework of the bilateral Slovene-German project ‘Solidification
Modelling by the Boundary Element Method’.
REFERENCES
1. P. W. Partridge, C. A. Brebbia and L. C. Wrobel, ¹he Dual Reciprocity Boundary Element Method, Elsevier, London,
1992.
2. L. C. Wrobel, C. A. Brebbia and D. Nardini, ‘The dual reciprocity boundary element formulation for transient heat
conduction’, in A. Sa da Costa et al. (eds.), Advances in ¼ater Resources »I, Springer, Berlin, 1986.
3. L. C. Wrobel and C. A. Brebbia, ‘The dual reciprocity boundary element formulation for nonlinear diffusion
problems’, Comput. Meth. Appl. Mech. Engng., 65, 147—164 (1987).
4. L. C. Wrobel and D. B. DeFigueredo, ‘A dual reciprocity boundary element formulation for convection-diffusion
problems with variable velocity fields’, Engng. Anal., 8, 312—319 (1991).
5. B. S[ arler, ‘Boundary integral formulation of general source based method for convective-diffusive solid liquid phase
change problems’, in C. A. Brebbia et al. (eds.), Boundary Elements X»III, CMP, Southampton, 1996, pp. 551—560.
6. T. Yamada, L. C. Wrobel and H. Power, ‘On the convergence of the dual reciprocity boundary element method’,
Engng. Anal., 13, 291—298 (1994).
7. B. S[ arler, ‘The dual reciprocity boundary element method solution of the continuous casting problem’, in J. I. Frankel,
C. A. Brebbia and M. A. H. Aliabadi (eds.), BE¹ECH-XII, CMP, Boston, 1997, pp. 309—320.
8. M. A. Golberg and C. S. Chen, ‘The theory of radial basis functions applied to the BEM for inhomogenous partial
differential equations’, Boundary Elements Comm., 5, 57—61 (1994).
9. J. Duchon, ‘Splines minimizing rotation invariant seminorms in Sobolev spaces’, in W. Schemp and K. Zeller (eds.),
Constructive ¹heory of Functions of Several »ariables, Lecture Notes in Math., vol. 571, Springer, Berlin, 1977, pp.
85—110.
10. M. D. Buhmann, ‘Multivariate cardinal interpolation with radial basis functions’, Constr. Approx., 6, 225—255 (1990).
11. W. J. Mansur, F. C. Araujo and J. E. B. Malaghini, ‘Solution of BEM systems of equations via iterative techniques’,
Int. J. Numer. Meth. Engng., 33, 1823—1844 (1992).
12. L. P. S. Barra, A. L. G. A. Coutinho, W. J. Mansur and J. C. F. Telles, ‘Iterative solution of BEM equations by
GMRES algorithm’, Comput. Struct., 44, 1249—1253 (1992).
( 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
732
V. BULGAKOV, B. S[ ARLER AND G. KUHN
13. L. P. S. Barra, A. L. G. A. Coutinho, J. C. F. Telles and W. J. Mansur, ‘Multi-level hierarhical preconditioners for
boundary element systems’, Engng. Anal. Boundary Elements, 12, 103—109 (1993).
14. V. E. Bulgakov, R. A. Bialecki and G. Kuhn, ‘Coarse division transform based preconditioner for boundary element
problems’, Int. J. Numer. Meth. Engng., 38, 2115—2129 (1995).
15. M. Hribers\ ek and L. S[ kerget, ‘Iterative methods in solving Navier—Stokes equations by the boundary element
method’, Int. J. Numer. Meth. Engng., 39, 115—139 (1996).
16. T. J. Urekew and J. J. Rencis, ‘The importance of diagonal dominance in the iterative solution of equations generated
from the boundary element method’, Int. J. Numer. Meth. Engng., 36, 3509—3527 (1993).
17. K. Guru Prasad, J. H. Kane, D. E. Keyes and C. Balakrishna, ‘Preconditioned Krylow solvers for BEA’, Int. J. Numer.
Meth. Engng., 37, 1651—1672 (1994).
18. V. E. Bulgakov, M. Merkel, R. A. Bialecki and G. Kuhn, ‘On the iterative solution of large-scale BEM industrial
problems’, submitted.
19. K. Guru Prasad, D. E. Keyes and J. H. Kane, ‘Generalized wavelet bases for boundary element matrices’, submitted.
20. P. Sonneveld, ‘CGS, A fast Lanczos-type solver for nonsymmetric linear systems’, SIAM J. Sci. Statist. Comput., 10,
36—52 (1989).
21. P. Sonneveld, P. Wesseling and P. M. de Zeeuw, ‘Multigrid and conjugate gradient methods as convergence
acceleration techniques’, in D. J. Paddon and H. Holstein (eds.), Multigrid Methods for Integral and Differential
Equations, Clarendon Press, Oxford, 1985, pp. 117—167.
22. Y. Saad and M. H. Schultz, ‘GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear
systems’, SIAM J. Scient. Statist. Comput., 7, 856—869 (1986).
23. P. S. Shields, ¸inear Algebra, Addison-Wesley, Reading, MA, 1964.
24. W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, Numerical Recipes in Fortran, Cambridge
University Press, Cambridge, 1992.
25. B. W. Char et al., MAP¸E 5, ¸ibrary Reference Manual, Springer, Berlin, 1991.
( 1998 John Wiley & Sons, Ltd.
Int. J. Numer. Meth. Engng. 43, 713—732 (1998)
Документ
Категория
Без категории
Просмотров
3
Размер файла
233 Кб
Теги
321
1/--страниц
Пожаловаться на содержимое документа