A Fast Multipole Accelerated Singular Boundary Method for Potential Problems

W.ChenC.J.Liuand Y.Gu

A Fast Multipole Accelerated Singular Boundary Method for Potential Problems

W.Chen1,2,C.J.Liu1and Y.Gu2,3

The singular boundary method(SBM)is a recently-developed meshless boundary collocation method.This method overcomes the well-known fictitious boundary issue associated with the method of fundamental solutions(MFS)while remaining the merits of the later of being truly meshless,integral-free,and easy-to-program.Similar to the MFS,this method,however,produces dense and unsymmetrical coef ficient matrix,which although much smaller in size compared withdomaindiscretizationmethods,requiresO(N2)operationsintheiterativesolution of the resulting algebraic system of equations.To remedy this bottleneck problem for its application to large-scale problems,this paper makes the first attempt to develop a fast multipole SBM(FM-SBM)formulation for two-dimensional(2D)potential problems.The proposed strategy can solve large-scale problems with several millions boundary discretization nodes on a desktop computer.

Singular boundary method,method of fundamental solutions,fast multipole method,meshless boundary collocation method,potential problem.

1 Introduction

The singular boundary method(SBM)[Chen and Wang(2010)]is a recentlydeveloped meshless boundary collocation method[Atluri et al.(2006);Dong and Atluri(2011),(2012);Zhang et al.(2002);Zhang et al.(2003)]and has since been successfully applied to a variety of physical problems,such as potential[Chen and Gu(2012);Gu and Chen(2013);Gu et al.(2012a)]and elastostatic problems[Gu et al.(2014);Gu et al.(2012b);Gu et al.(2011)].This method inherits the merits of the method of fundamental solutions(MFS)[Chen(1995);Fairweather andKarageorghis(1998);Sarler(2009);Sarler and Liu(2013);Young et al.(2007)]and is truly free of mesh and integration.The SBM avoids the perplexing fictitious boundary issue associated with the traditional MFS.However,the application of this method has so far been limited to model relatively small scale problems.This is mainly because the SBM produces dense interpolation matrix that,although smaller in size compared with the domain discretization methods such as the FEM,requires O(N2)operations in the iterative solution of its resulting algebraic equations.The SBM solution of large-scale problems is computationally expensive.

Inspired by the successful application of the fast multipole method(FMM)[Greengard and Rokhlin(1987);Rokhlin(1985)]in signi ficantly reducing computing costs of the desnse matrix equations from boundary integraltion methods(BIEs)[Liu and Nishimura(2006);Liu et al.(2005);Nishimura(2002);Wang and Yao(2004)],this study makes the first attempt to apply the FMM to accelerating the solutions of the SBM discretization equations.Applying the FMM to the SBM can potentially reduce both the operation count and memory requirement to O(N).This will be a signi ficant improvement to the SBM in extending its applications.

This paper demonstrates that the proposed FM-SBM strategy is capable of solving large-scale problems.The FM-SBM formulations and its numerical implementation will be discussed in details.Numerical examples with up to 4,000,000 collocation points are successfully solved on a Core II PC using the developed FM-SBM code.

A brief outline of the rest of this paper is as follows:Section 2 reviews the conventional SBM formulations for potential problems.Section 3 provides detailed formulations and algorithms used in the FM-SBM,followed in Section 4 two benchmark numerical examples are well studied to illustrate the accuracy and ef ficiency of the proposed strategy.Finally,the conclusions and remarks are provided in Section 5.

2 Formulations of the conventional SBM

Without loss of generality,we consider the following Laplace equation governing potential problems in a 2D domain Ω(see Fig.1)

Figure 1:Domain Ω and its boundary S.

whereudenotes the potential field,Ω is a bounded domain in R2with boundaryS=S1+S2which is assumed to be piecewise smooth,n presents the outward normal,and the barred quantities indicate given values along the boundary.

In the SBM,the collocation and source points are coincident and are both placed on the physical boundary,without requiring a fictitious boundary as compared with the traditional MFS.The SBM interpolation formula for 2D potential problems is given by[Chen and Gu(2012);Gu et al.(2012a)]

where xiis theithcollocation point,sjthejthsource point,αjthejthunknown intensity of the distributed source at sj,Nthe number of the collocation points,and

the fundamental solutions for 2D potential problems.

In the above Eqs.(4)and(5),uiiandqiiare de fined as origin intensity factors,i.e.,the diagonal elements of the SBM interpolation matrix.The fundamental assumption of the SBM is the existence of the origin intensity factor upon the singularity of the coincident source-collocation point for mathematically well-posed problems.Our numerical experiments illustrate that the origin intensity factors do exist and have finite values,only depending on the distribution of boundary nodes and their respective boundary conditions.We refer readers to Refs.[Gu et al.(2014);Gu et al.(2011)]for further detailed derivation of the aforementioned origin intensity factors.

3 The fast multipole method for SBM

The key idea of the FMM is to convert the point-to-point interactions to cell-to-cell interactions,in which the procedure is ef ficiently carried out with tree structures.An iterative solver,such as GMRES,is employed to reduce the computational complexity of the matrix-vector multiplication from O(N2)to O(N).Using the FMM,both the solution time and the memory requirements can be reduced to order of O(N)[Rokhlin(1985)].

This section will present a FMM formulation for the SBM,which is similar to that for the BEM,except for the calculation of the kernel expansions where integrals are replaced by summations.For completeness,all the required formulas are provided.The detailed derivation of some of these formulas can be found in Refs.[Liu and Nishimura(2006);Liu et al.(2005)].

For convenience,complex notations are used to represent the collocation point x and the source point s by z0=x0+iy0and z=x1+iy1,respectively(see Fig.2).The fundamental solution(6)can then be rewritten as

in which the real part gives the values of fundamental solution.And the kernels A and B shown blow related to the G kernel,can be rewritten as

with n=n1+in2,in complex notation.

3.1 Multipole moments

Suppose that zcis a point close to source point z,that is,inequalityis valid(see Fig.2).Using the subtracting and adding-back technique,

Figure 2:Main steps of FMM.

Eq.(8)can be expressed as

Applying the Taylor series expansion to the second log term of Eq.(11),we get the following equation

Considering a group of source points,when the collocation point z0is suf ficiently far away from this group,the sum∑jG(z0,zj)αjon the right-hand side of Eq.(4)can be obtained by the following multipole expansion and only need to be computed once.

In addition,for the sum on the right-hand side of Eq.(5),we can also derive the following multipole expansion

by using the same moments given in Eq.(14).

In order to calculate the diagonal coef ficients of Neumann boundary Eq.(5),the sumshould be calculated.By using Eqs.(10)and(12),we can obtain the following multipole equation

are another group of moments about zc.

3.2 M2M translation

Suppose that the point zcis moved to a new location zc′(see Fig.2),we obtain the following equations

When applying the binomial formula and rearranging the terms to Eqs.(18)and(19),respectively,we can further get

Note that Eqs.(20)and(21)translate the moments centered at zcto the moments centered at zc′,which are called M2M translation.

3.3 Local expansions and M2L translation

SupposezLisapointclosetothecollocationpointz0,thatis,theinequalityis valid(see Fig.2).Using the multipole expansions Eqs.(13)and(16),we can get

respectively.

Expanding f(z0)and g(z0)about zLusing the Taylor series expansion,we can derive the following local expansions

which are called M2L translation.

Using Eq.(26),the local expansion for the A kernel sum can be naturally obtained by

3.4 L2L translation

Suppose the point for local expansion is moved from zLto zL′(see Fig.2),using Eq.(24),we obtain the following local equation

Applying the binomial formula and the relationto Eq.(29),we can derive the following equation

where the new coef ficients are given by the following L2L translation

Note that Eqs.(24)and(25)have the same form,and the L2L translations for B kernel sum can be easily obtained by replacingrespectively.

It is noted that the M2M and L2L translations for A kernel sum and G kernel sum remain the same.

3.5 FMM formulation for the SBM

The algorithms or procedures for the FM-SBM are similar to those for the FMBEM[Liu et al.(2005);Nishimura(2002)].The main steps in the FM-SBM can be described as follows:

Step 1.Discretization.Discretize the boundary S in the same manner as in the conventional SBM,that is,place a certain number of collocation points,which are also the source points at the same time,on S as shown in Fig.1.

Step 2.Construction of quad-tree structure.Consider a square that covers all the collocation points.Call this square the cell of level 0(see Fig.3).Now take a cell(a parent cell)of levell(l≥0)and divide it into four equal sub squares whose edge length is half of that of the parent cell and call any of them a cell(a child cell)of levell+1 if some collocation points belong to this square(see Fig.3).Stop dividing a cell if the number of points in that cell is less than a given number as one in the example is shown in Fig.3.A childless cell is called a leaf.

Figure 3:Construction of a tree structure.

Step 3.Upward pass.First we compute the multipole moments associated with leaves using Eqs.(14)and(17)for up topterms.For a parent cell,we compute the multipole moments by adding all the moments of its children using the M2M translations Eqs.(20)and(21),in whichzc′is the centroid of the parent cell andzcthe centroid of a child cell(see Fig.2).This procedure is repeated forl≥2 tracing the tree structure upward.

Step 4.Downward pass.Two cells are called"adjacent cells at levell"if these cells are both of the levelland share at least one common vertex.Two cells are"well-separated at levell"if they are not adjacent at levellbut their parent cells are adjacent at levell−1.The list of all the well-separated cells from a levellcellCis called the interaction list ofC.Now we compute the local expansion associated with a cellC.The local expansion associated withCis the sum of the contributions due to the collocation points in cells of the interaction list ofCand the contributions due to collocations points in cells which are not adjacent to the parent of cellC.The former is calculated by using M2L translations Eqs.(26)and(27),with moments associated with cells in the interaction list,and the latter by using L2L translations which shifts the local expansions from the centroid of C,s parent to the centroid of cellC(see Fig.2).We repeat these procedures starting from level 2 and tracing the tree structure downward to all leaves.For a cell at level 2,the coef ficients of the local expansions are only involved M2L translation.

Step 5.Evaluation of the sum in Eq.(4)or(5).LetC be a leaf to which the collocation point z0belongs,and the values of αjbe given.We calculate the contributions from source points in leafC and its adjacent cells directly as in conventional SBM.We also compute contributions from all other cells by using the coef ficients of the local expansions associated withC which have been computed in Step 4.

Step 6.First,we let αjequals 1 to compute the diagonal coef ficients of the coefficient matrix in Neumann boundary equation(5)using Step 5.Then,we update the unknown vector in the system produced by using Eqs.(4)or(5),and go back to Step 3 for the matrix and unknown vector multiplication until the solution converges within a given residue using generalized minimal residual(GMRES)solver.

3.6 Preconditioners for the FM-SBM

Similar to FM-BEM,a good pre-conditioner for FM-SBM is crucial for its convergence and computing ef ficiency.Since the conditioning number of the SBM is almost the same as that of the BEM,the present FM-SBM directly utilizes the block diagonal pre-conditioner proposed for FM-BEM[Wang and Yao(2004)].The preconditioner is based on leaves of the hierarchical tree,which is formed on each leaf using direct evaluation of the fundamental solution with the collocation points within that leaf and their corresponding source points.This pre-conditioner will be inef ficient if the number of collocation points in a leaf is large.

4 Numerical examples

To verify the feasibility and ef ficiency of the proposed FM-SBM,two numerical examples are examined in which the results are compared with the corresponding exact solutions and those obtained using conventional SBM.To assess the accuracy of the FM-SBM,the numerical results will be compared with exact solutions in terms of the relative error de fined by

Figure 4:An amoeba-like domain with ten holes.

4.1 An amoeba-like domain with ten holes

First,we consider a simple heat condition problem in an amoeba-like domain with ten inner holes(see Fig.4),whose outer boundary is de fined as

and the ten inner holes are randomly distributed,whose radius are all equal to 0.2.In this case,the analytical solution is given by

The Dirichlet boundary conditions are speci fied on the outer boundary while the Numann boundary conditions are speci fied on the rest of the boundaries.

We discretize the inner and outer boundaries with the same number of boundary points,and every inner holes boundary has the same points.The conventional SBM code uses direct solver Gauss elimination for solving the linear system equations.For FM-SBM,the numbers of the terms for both moments and local expansions were set to 20,and the maximum number of collocation points allowed in a leaf of the tree structure was prescribed to 400.The residue for iteration is taken as 1×10−6,and the initial guess is zero.

With the number of collocation points varying from 2,000 to 360,000,the relative error of u and q by either the FM-SBM or the conventional SBM are listed in Table 1.It can be observed from this table that the FM-SBM performs as well asthe conventional SBM when the number of collocation points less than 8,000.As the number of boundary collocation point further increases,the conventional SBM becomes infeasible due to explosive increase of memory demands.On the other hand,the proposed FM-SBM can continue to provide accurate results even asNreaches 360,000 without any dif ficulty.

Table 1:Relative errors of computed potential and normal derivative.

The comparisons of CPU time and memory requirement between FM-SBM and conventional SBM are illustrated in Figs.5 and 6,respectively.To show the order of complexity of FM-SBM more clearly,logarithmic scale was used in Figs.5 and 6.The vertical axis of Fig.5 denotes the 10-base logarithmic of CPU time,simplified aslog(t),and the horizontal axis represents the 10-base logarithmic of number of collocation points,simpli fied aslog(N).Similarly,the vertical axis of Fig.6 denotes the 10-base logarithmic of memory requirement,simpli fied aslog(m),and the horizontal axis represents the 10-base logarithmic of number of collocation points,simpli fied aslog(N).These figures show a signi ficant edge of the present FM-SBM in the savings of CPU time and memory,which are observed both approximately proportional to the number of collocation points.The results verify that both multiplication operations and storage costs of FM-SBM areO(N).This shows that the FM-SBM has great promise to handle large-scale problems.

Fig.7 show the iteration number and maximum number of boundary points in a leaf node vary with the number of collocation points,respectively,as the residue for iteration and the maximum of points allowed in a leaf of the tree structure were fixed.It can be concluded that,in general,the iteration number and the maximum number of boundary increase with the number of collocation points.

The relationship between the iteration number and residue of iteration is given in

Figure 5:Comparison of the CPU time for the problem of the amoeba-like domain with ten holes.

Figure 6:Comparison of the memory requirement for the problem of the amoebalike domain with ten holes.

Figure 7:Relationships of iteration number and maximum number of points in leaves with the number of collocation points.

Figure 8:Iteration number versus residue of iteration under different scale of problems.

Fig.8.From this figure,we can obtain that,when the residue of iteration gets less,FM-SBM needs more iteration steps to get the resolution.And,this tendency doesn’t matter with the scale of the problem.We can also get that,for the same residue of iteration,when the scale becomes larger,FM-SBM will need more iteration steps.

Table 2:Relative errors of computed potential and normal derivative.

In Table 2,Nleafand Nleveldenote the number of collocation points in a leaf and the number of levels of the tree structure,respectively.Table 2 shows that,when the number of collocation points increases,the number of levels decreases.As we also know,when the number of collocation points increases,it will need more CPU time to calculate the interaction of nodes in a leaf.So,as illustrated in Table 2,there will be a better value of the number of collocation points which is 400 in this example to get less CPU time.

4.2 A square plate with numerous circular holes

Here,the thermal conductivity of a square plate with uniformly distributed circular holes,or a perforated plate,is studied using the FM-SBM(see Fig.9).Several models of the plate,with increasing number of holes,are considered.Each model of the plate contains a total of 2m×2m holes,with m=2,3,...,28,30.The radii of the holes are the same and equal to 4.0/(2m+2).Each hole is discretized with 100 collocation points and the outer boundary of the plate is discretized with 400mcollocation points.

In this example,the analytical solution is given by

The Dirichlet boundary conditions are speci fied on the outer boundary,and the Numann boundary condition are speci fied on the remaining boundaries.For FMSBM,the number of terms equals 20 for multipole and local expansions,and the maximum of collocation points allowed in a leaf is set to 400.The residue for iteration is prescribed to 1×10−6,and the initial guess is zero.

Figure 9:A model of the square plate with 400 uniformly distributed circular holes

Table 3:Relative error of computed potential and normal derivative for the square plate with circular holes.

Table 3 illustrates the relative errors ofuandqby both the FM-SBM and conventional SBM.As shown in this table,the conventional SBM can only provide accurate numerical solution when the number of circular holes less than 100,after which the method becomes infeasible due to memory limitations.In sharp contrast,the proposed FM-SBM is still workable when the number of circular holes is as large as 3,600,with the largest relative error less than 3E-5.

Figure 10:Comparison of the CPU time for the square plate with circular holes.

Figure 11:Comparison of the memory for the square plate with circular holes.

The relationships of iteration number and maximum number of points in leaves with the number of collocation points,is shown in Fig.12.From this figure,we can conclude that the iteration number and the maximum number of boundary increase with the number of collocation points,generally.

Figure 12:Relationships of iteration number and maximum number of points in leaves with the number of collocation points.

5 Conclusions

In this paper,a fast multipole accelerated singular boundary method is presented for 2D potential problems.Complete formulations and implementation details for FM-SBM are provided.Two numerical examples demonstrate the ef ficiencies of the present strategy for solving large-scale problems.

The convergence issue in the FM-SBM requires further investigations,and a better pre-conditioner is required to devise.Furthermore,the algorithms of the FM-SBM can also be improved.In this work,the calculation of the potential and its derivatives were done by direct calculation using SBM interpolation formula,which may further be speeded up by applying FMM,when the number of data points is comparable to that of collocation points.All these issues merit further investigations based on this study.

Acknowledgement:The work described in this paper was supported by the National Basic Research Program of China(973 Project No.2010CB832702),the National Science Funds for Distinguished Young Scholars of China(No.11125208),andtheNationalScienceFundsofChina(No.11402075,11302069),andtheChina Postdoctoral Science Foundation(No.2015M570572).

Atluri,S.N.;Liu,H.T.;Han,Z.D.(2006):Meshless local Petrov-Galerkin(MLPG)mixed collocation method for elasticity problems.CMES:Computer Modeling in Engineering&Sciences,vol.4,no.3,pp.141-152.

Chen,C.S.(1995):The method of fundamental solutions for non-linear thermal explosions.Commun.Numer.Methods.Eng.,vol.11,no.8,pp.675-681.

Chen,W.;Gu,Y.(2012):An improved formulation of singular boundary method.Adv.Appl.Math.Mech.,vol.4,no.5,pp.543-558.

Chen,W.;Wang,F.Z.(2010):A method of fundamental solutions without fictitious boundary.Eng.Anal.Bound.Elem.,vol.34,no.5,pp.530-532.

Dong,L.;Atluri,S.N.(2011):Development of T-Trefftz four-node quadrilateral and Voronoi Cell Finite Elements for macro-µmechanical modeling of solids.CMES:Computer Modeling in Engineering&Sciences,vol.81,no.1,pp.69-118.

Dong,L.;Atluri,S.N.(2012):A simple multi-source-point Trefftz method for solvingdirect/InverseSHMproblemsofplaneelasticityinarbitrarymultiply-connected domains.CMES:Computer Modeling in Engineering&Sciences,vol.85,no.1,pp.1-44.

Fairweather,G.;Karageorghis,A.(1998):The method of fundamental solutions for elliptic boundary value problems.Adv.Comput.Math.,vol.9,no.1,pp.69-95.

Greengard,L.;Rokhlin,V.(1987):A fast algorithm for particle simulations.J.Comput.Phys.,vol.73,no.2,pp.325-348.

Gu,Y.;Chen,W.(2013):In finite domain potential problems by a new formulation of singular boundary method.Appl.Math.Model.,vol.37,no.4,pp.1638-1651.

Gu,Y.;Chen,W.;He,X.-Q.(2012a):Singular boundary method for steady-state heat conduction in three dimensional general anisotropic media.Int.J.Heat.Mass.Transf.,vol.55,no.17–18,pp.4837-4848.

Gu,Y.;Chen,W.;He,X.(2014):Improved singular boundary method for elasticity problems.Comput.Struct.,vol.135,no.0,pp.73-82.

Gu,Y.;Chen,W.;He,X.Q.(2012b):Domain-Decomposition Singular Boundary Method for Stress Analysis in Multi-Layered Elastic Materials.CMC:Computers Materials&Continua,vol.29,no.2,pp.129-154.

Gu,Y.;Chen,W.;Zhang,C.-Z.(2011):Singular boundary method for solving plane strain elastostatic problems.Int.J.Solids.Struct.,vol.48,no.18,pp.2549-2556.

Liu,Y.J.;Nishimura,N.(2006):The fast multipole boundary element method for potential problems:A tutorial.Engineering Analysis with Boundary Elements,vol.30,no.5,pp.371-381.

Liu,Y.J.;Nishimura,N.;Yao,Z.H.(2005):A fast multipole accelerated method of fundamental solutions for potential problems.Engineering Analysis with Boundary Elements,vol.29,no.11,pp.1016-1024.

Nishimura,N.(2002):Fastmultipoleacceleratedboundaryintegralequationmethods.Applied Mechanics Reviews,vol.55,no.4,pp.299-324.

Rokhlin,V.(1985):Rapid solution of integral equations of classical potential theory.Journal of Computational Physics,vol.60,no.2,pp.187-207.

Sarler,B.(2009):Solution of potential flow problems by the modi fied method of fundamental solutions:Formulations with the single layer and the double layer fundamental solutions.Eng.Anal.Bound.Elem.,vol.33,no.12,pp.1374-1382.

Sarler,B.;Liu,Q.G.(2013):Non-singular method of fundamental solutions for two-dimensional isotropic elasticity problems.CMES:Computer Modeling in Engineering&Sciences,vol.91,no.4,pp.235-266.

Wang,H.T.;Yao,Z.H.(2004):Application of a new fast multipole BEM for simulation of 2D elastic solid with large number of inclusions.Acta Mechanica Sinica,vol.20,no.6,pp.613-622.

Young,D.L.;Chen,K.H.;Chen,J.T.;Kao,J.H.(2007):A modi fied method of fundamental solutions with source on the boundary for solving Laplace equations with circular and arbitrary domains.CMES:Computer Modeling in Engineering&Sciences,vol.19,no.3,pp.197-221.

Zhang,J.M.;Yao,Z.H.;Li,H.(2002):A hybrid boundary node method.Int.J.Numer.Methods.Eng.,vol.53,no.4,pp.751-763.

Zhang,J.M.;Yao,Z.H.;Tanaka,M.(2003):The meshless regular hybrid boundary node method for 2D linear elasticity.Eng.Anal.Bound.Elem.,vol.27,no.3,pp.259-268.

1State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering,International Center for Simulation Software in Engineering and Sciences,College of Mechanics and Materials,Hohai University,Nanjing 210098,China.

2Corresponding author.E-mail:chenwen@hhu.edu.cn;guyan1913@163.com

3College of Mathematics,Qingdao University,Qingdao 266071,PR China.