May 23, 2005Lecture 24CME342/AA220 /CS238 - Parallel Methods in Numerical AnalysisFast methods for N-body problem IIParticle simulationst = 0while t < t_final for i = 1 to n … n = number of particles compute f(i) = force on particle i for i = 1 to n move particle i under force f(i) for time dt … using F=ma compute interesting properties of particles (energy, etc.) t = t + dtend whileforce =externalF+neighborF+farfieldFO( N) O( N)O( N2)Using a classical approach, this will require O(N^2) operations,but with a clever divide and conquer approach we can solve theproblem in O(N logN) or O(N) operationsFast hierarchical methods•Use of adaptive quadtree or octree to subdivide the space•Barnes-Hut algorithm:•can solve the problem in O(N logN) •not very accurate •use the center of mass and total mass•Fast Multipole method :•can solve the problem in O(N) •accurate•use more information than BH•more difficult to implementAdaptive quad tree•All our algorithms begin by constructing a tree to hold all the particles•Adaptive Quad (Oct) Tree only subdivides space where particles are locatedChild nodes enumerated counterclockwisefrom SW corner, empty ones excludedFast Multipole method• “A fast algorithm for particle simulation”, Greengardand Rokhlin, J. Comp. Phys. V. 73, 1987, manylater papers• Differences from Barnes-Hut– FMM computes the potential at every point, not just the force– FMM uses more information in each box than the CM and TM, so it isboth more accurate and more expensive– To compensate, FMM accesses a fixed set of boxes at every level,independent of D/r– BH uses fixed information (CM and TM) in every box, but # boxesincreases with accuracy. FMM uses a fixed # boxes, but the amount ofinformation per box increases with accuracy.Fast Multipole method• Use the potential instead of the force• FMM uses two kinds of expansions– Outer expansions represent potential outside node due toparticles inside, analogous to (CM,TM)– Inner expansions represent potential inside a node due toparticles outside;• Computing these for every leaf node is thecomputational goal of FMMOuter expansion• Outer(n) = (M, 1 , 2 , … , r, zn )– Stores data for evaluating potential (z) outside node n dueto particles inside n• Will be computed for each node in QuadTree• Cost of evaluating (z) is O( p ), independent of thenumber of particles inside n• Outer expansion needed to compute the potential away from ndue to particles inside n. We also want to compute potentialinside n due to particles away from n1(z) =1M log(z  z1) + 1( j)(z  z1)jj =1pInner expansion• Inner(n) = ( 1 , 2 , … , r , zn )– Stores data for evaluating potential (z) inside node n dueto particles outside n• Computing this for every leaf node is thecomputational goal of FMMc(z) = (i)(z  zc)ii=0pConverting Outer ExpansionsWe want to compute Outer(n) cheaply from Outer( c ) for all children of nHow to combine outer expansions around different points?– First step: make them expansions around same pointn1 is a subsquare of n2,zk = center(nk) for k=1,2Outer(n1) expansion accurate outsideblue dashed square, so also accurateoutside black dashed squareSo there is an Outer(n2) expansionwith different ak and center z2 whichrepresents the same potential asOuter(n1) outside the black dashed boxConverting outer expansionsGivenSolve for M2 and 2 inGet M2 = M1 and each 2 is a linear combination of the 1multiply r-vector of 1 values by a fixed r-by-r matrix to get 2( M2 , 12 , … , r2 ) = Outer_shift( Outer(n1) , z2 ) = desired Outer( n2 )1(z) =1M log(z  z1) + 1( j)(z  z1)jj =1p1(z) 2(z) =2M log(z z2) +2( j)(z z2)jj=1pConverting inner expansionsHow to combine inner expansions around different points?– First step: make them expansions around same pointn1 is a subsquare of n2,zk = center(nk) for k=1,2Inner(n2) expansion accurate insideblack dashed square, so also accurateinside blue dashed squareSo there is an Inner(n1) expansionwith different k and center z1 whichrepresents the same potential asInside(n2) inside the blue dashed boxConverting Inner ExpansionsGivenSolve for 2 such thatEach 2 is a linear combination of the 1multiply p-vector of 1 values by a fixed p-by-p matrix to get 2(12 , … , r2 ) = Inner_shift( Inner(n1) , z2 ) = desired Inner( n2 )1(z) = (i)(z  z1)ii=0p 2(i)(z  z2)ii=0p= 1(i)(z  z1)ii= 0pConverting Outer to InnerWe need to show how to convert Outer(n3) to Inner(n4)when n4 is far enough from n3Inner(n4) = ( 0 , 1 , … , r , z4 )Outer(n3)=(M, 1 , 2 , … , r , z3 )Once again is a weighted sumInner(n4)=Convert(Outer(n3),center(n4))3Mlog( z  z3) + 3( j)(z  z3)j=j=1p4(i)(z  z4)ii=1pTop Level Description of FMM1) Build the QuadTree2) Call Build_Outer(root), to compute outer expansions of each node n in the QuadTree … Traverse QuadTree from bottom to top, … combining outer expansions of children … to get out outer expansion of parent3) Call Build_ Inner(root), to compute inner expansions of each node n in the QuadTree… Traverse QuadTree from top to bottom, … converting outer to inner expansions … and combining them4) For each leaf node n, add contributions of nearest particles directly into Inner(n) … final Inner(n) is desired output: expansion for potential at each point due to all particlesCompute Outer(n) for each node… Compute Outer(n) for each node of the QuadTreeouter = Build_Outer( root )function ( M, a1,…,ar , zn) = Build_Outer( n ) … compute outer expansion of node n if n if a leaf … it contains 1 (or a few) particles compute and return Outer(n) = ( M, a1,…,ar , zn) directly from its definition as a sum else … “post order traversal”: process parent after all children Outer(n) = 0 for all children c(k) of n … k = 1,2,3,4 Outer( c(k) ) = Build_Outer( c(k) ) Outer(n) = Outer(n) + Outer_shift( Outer(c(k)) , center(n)) … just add component by component endfor return Outer(n)end ifCost = O(# nodes in QuadTree) = O( N ) same as for