4/19/2004 CS267, Yelick 1CS 267: Applications of Parallel ComputersTree-Structured Codes for N-Body SimulationsKathy Yelickhttp://www.cs.berkeley.edu/~yelick/cs2674/19/2004 CS267, Yelick 2Motivation• Particle methods are used for a variety of applications• Astrophysics• The particles are stars or galaxies• The force is gravity• Particle physics• The particles are ions, electrons, etc.• The force is due to Coulomb’s Law• Molecular dynamics• The particles are atoms or • The forces is electrostatic• Vortex methods in fluid dynamics• Particles are blobs of fluid4/19/2004 CS267, Yelick 3Outline• Motivation• Obvious algorithm on N bodies takes O(N2) work• How to reduce the number of particles in the force sum• We must settle for an approximate answer (say 2 decimal digits, or perhaps 16 …)• Basic Data Structures: Quad Trees and Oct Trees• The Barnes-Hut Algorithm (BH)• An O(N log N) approximate algorithm for the N-Body problem• The Fast Multipole Method (FMM)• An O(N) approximate algorithm for the N-Body problem• Parallelizing BH, FMM and other algorithms• Example applications• Alternative approach: lots of hardware4/19/2004 CS267, Yelick 4Particle Simulation• Let f(i) be the force on particle i• f(i) = external_force + nearest_neighbor_force + NBody_forcet = 0while t < t_finalfor i = 1 to n … n = number of particlescompute f(i) = force on particle ifor i = 1 to nmove particle i under force f(i) for time dt … using F=macompute interesting properties of particles (energy, etc.)t = t + dtend while4/19/2004 CS267, Yelick 5Particle Simulationf(i) = external_force + nearest_neighbor_force + NBody_force• External_force (e.g., current) is usually embarrassingly parallel, O(N) • Nearest_neighbor_force is with a few neighbors, so still O(N)• N-Body_force (gravity or electrostatics) requires all-to-all interactions•f(i) = Σk!=if(i,k) … f(i,k) = force on i from k•f(i,k) = c*v/||v||3 in 3 dimensions or•f(i,k) = c*v/||v||2 in 2 dimensions• c = product of masses or charges and appropriate constants• v = vector from particle i to particle k, ||v|| = length of v •Obvious algorithm costs O(N2), but we can do better...4/19/2004 CS267, Yelick 6Reducing the Number of Particles in the Sum• Previous divide and conquer algorithms use same intuition• Consider computing force on earth due to all celestial bodies• Look at night sky, # terms in force sum >= number of visible stars• One “star” is really the Andromeda galaxy, which is billions of stars• A lot of work if we compute this per star … • OK to approximate all stars in Andromeda by a single point at its center of mass (CM) with same total mass• D = size of box containing Andromeda , r = distance of CM to Earth• Require that D/r be “small enough”• Idea not new: Newton approximated earth and falling apple by CMsEarthr = distance to center of massx = location of center of massAndromedaDDrx4/19/2004 CS267, Yelick 7What is new: Using points at CM Recursively• From Andromeda’s point of view, Milky Way is also a point mass• Within Andromeda, picture repeats itself• As long as D1/r1 is small enough, stars inside smaller box can be replaced by their CM to compute the force on Vulcan• Boxes nest in boxes recursivelyEarthAndromedaDDrxVulcanD1D1Replacing clusters by their Centers of Mass Recursivelyr14/19/2004 CS267, Yelick 8Quad Trees• Data structure to subdivide the plane• Nodes can contain coordinates of center of box, side length• Eventually also coordinates of CM, total mass, etc.• In a complete quad tree, each nonleaf node has 4 children4/19/2004 CS267, Yelick 9Oct Trees• Similar Data Structure to subdivide 3D space• Analogous to 2D Quad tree--each cube is divided into 8 sub-cubesTwo Levels of an OctTree4/19/2004 CS267, Yelick 10Using Quad Trees and Oct Trees• All our algorithms begin by constructing a tree to hold all the particles• Interesting cases have non-uniform particle distribution• In a complete tree (full at lowest level), most nodes would be empty, a waste of space and time• Adaptive Quad (Oct) Tree only subdivides space where particles are located • More compact and efficient computationally, but harder to program4/19/2004 CS267, Yelick 11Example of an Adaptive Quad TreeChild nodes enumerated counterclockwise from SW cornerEmpty ones excludedAdaptive quad tree where no space contains more than 1 particle4/19/2004 CS267, Yelick 12Adaptive Quad Tree Algorithmclass QuadTreebuild (particles)QuadTree t = new QuadTree();for each j in particles t.insert(j) … loop over all N particlest.insert(j) … insert particle j in QuadTreeend buildinsert(j) … Try to insert particle j at node n in QuadTreeif this node is empty … emptyadd j as the (only) particle in this nodeif this is an internal node (has 4 children) … internal determine which child c contains particle jc.insert(j)else (this quadtree contains 1 particle) … leafadd n’s 4 children to the QuadTreelet c be the child of the particle k already herec.insert(k) let c be the child of n containing jc.insert(j)end insert4/19/2004 CS267, Yelick 13Cost of Adaptive QuadTree ConstructionCost <= N * maximum cost of QuadTreeInsert = O(N * maximum depth of QuadTree)1. Uniform distribution =>depth of QuadTree = O(log N), so Cost = O( N log N)2. Arbitrary distribution => depth of Quad Tree = O(b), where b= #bits in particle coordinates, so Cost = O(bN)4/19/2004 CS267, Yelick 14Barnes-Hut Algorithm• “A Hierarchical O(n log n) force calculation algorithm”, J. Barnes and P. Hut, Nature, v. 324 (1986), plus many later papers.• Good for low accuracy calculations:RMS error = (Σk|| approx f(k) - true f(k) ||2/ || true f(k) ||2/N)1/2~ 1%(other measures better if some true f(k) ~ 0)• High Level Algorithm (in 2D, for simplicity)1) Build the QuadTree using QuadTree.build… already described, cost = O( N log N) or O(b N)2) For each node = subsquare in the QuadTree, compute the CM and total mass (TM) of all the particles it contains… “post order traversal” of QuadTree, cost = O(N log N) or O(b N)3) For each particle, traverse the QuadTree to compute the force on it, using the CM and TM of “distant” subsquares… core of algorithm … cost depends on accuracy desired but still O(N log N) or O(bN)4/19/2004 CS267, Yelick 15Step 2: Compute CM and TM of Each NodeCost = O(# nodes in QuadTree) = O( N )… Compute the CM =
View Full Document
Unlocking...