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)

1/--страниц