Class 20. N -body Techniques, Part 3The PM Method, ContinuedThere are several distinct steps in PM process:1. Assign particles to mesh to compute ρi.2. Get boundary conditions for Φ (Φ0and ΦN+1).3. Solve discretized version of Poisson’s equation.4. Compute F from discretized version of force equation.Step 1: Assigning particles to meshDiscuss two schemes here:1. Nearest Grid Point (NGP) scheme:• Assign entire mass of particle to grid zone that contains it.• E.g., discretize space into N zones in x-dimension:Set ρi= nim/∆x, where ni= number of particles in cell i (equal mass).• Leads to a very coarse distribution of ρi:2. Charge-In-Cell (CIC) or Particle-In-Cell (PIC):• Assign a “shap e” or “cloud” to particle.• Assume a distribution of ρ inside this shape.• Then distribute mass to zones according to overlap.• E.g., assume particle has top-hat ρ distribution, width w, height ρ0= m/w:1Then (in 1-D),R∞−∞ρ(x) dx = m. Distribute mass of particle according to overlap:Leads to smoother ρi.• Can adopt more complex shapes for density. E.g.,TriangleGaussianetc.Higher-order “shapes” introduce ringing into system.Step 2: Boundary conditions• Given ρi, i = 1, ..., N, need a boundary value for Φ, i.e., need Φ0and ΦN+1.• Often can use periodicBC, i.e., Φ0= ΦN, ΦN+1= Φ1. Appropriate for, e.g., cosmologysimulations.• Otherwise, standard to use multipole expansion (e.g., Jackson 1975) to compute po-tential on boundary due to mass in each cell.– Often, first (monopole) term is good enough:ΦB(r) = −NXi=1Gmi|r − ri|.– See Binney & Tremaine Eq. 2-122 for full series (involves spherical harmonics).Step 3: Solve Poisson’s equation• Can see that discretized equationΦi−1− 2Φi+ Φi+1(∆x)2= 4πGρileads to tri-diagonal (tri-di) matrix:−2 11 −2 11 −2............1 −2 11 −2Φ1Φ2Φ3......ΦN=4πGρ1(∆x)2− Φ04πGρ2(∆x)24πGρ3(∆x)2......4πGρN(∆x)2− ΦN+1.2• There is an extremely efficient algorithm for solving tri-di systems.– Write discretized system as:aiΦi−1+ biΦi+ ciΦi+1= di.– Then forward elimination gives (Hockney & Eastwood, p. 185):1w1=c1b1wi=cibi− aiwi−1,(i = 2, 3, ..., N − 1), and,g1=d1b1gi=di− aigi−1bi− aiwi−1.– Backsubstitution:ΦN= gNΦi= gi− wiΦi+1,with i = N − 1, N − 2, ..., 1.– If a, b, c constant, can precompute wiand 1/(bi− aiwi−1).– If a = 1, b = −2, c = 1, only need 4N operations.– For periodic BC, even simpler method possible (Hockney & Eastmood, p. 35).Step 4: Force interpolation• Once potential is given, must compute force (per unit mass) from F = −∇Φ.• In 1-D, F = −∂Φ/∂x ⇒ FDE Fi+1/2= −(Φi+1− Φi)/∆x.– Forces centered at cell boundaries:• Must interpolate forces to particle positions.• Linear interpolat io n simplest. For each particle, position xi−1/2< x < xi+1/2, compute:F(x) = Fi−1/2+x − xi−1/2∆xFi+1/2− Fi−1/2.• Higher-order interpolation used in conjunction with higher-order charge-assignmentschemes.1Also see tridag() (NRiC §2.4).3We now have ingredients necessary for a 1-D PM code:1. Particle assignment;2. Boundary conditions;3. Solve Poisson’s equation;4. Force interpolation.Result is F for every particle.Generalizing to 3-D• Generalizing to 3-D is straightforward:1. Particle a ssignment: use NGP; or for PIC, particle is sphere.2. BCs: periodic, or use 3-D multipole expansion.3. Solve Poisson’s equation in 3-D (see below).4. Interpolate F in 3-D (easy).• Poisson’s equation in 3-D:∂2Φ∂x2+∂2Φ∂y2+∂2Φ∂z2= 4πGρ.• Discretize Φ in 3-D:Φ(x, y, z) → Φi,j,k,ρ(x, y, z) → ρi,j,k.• FDE becomes:Φi−1,j,k− 2Φi,j,k+ Φi+1,j,k(∆x)2+Φi,j−1,k− 2Φi,j,k+ Φi,j+1,k(∆y)2+Φi,j,k−1− 2Φi,j,k+ Φi,j,k+1(∆z)2= 4πGρi,j,k.• Can be written in matrix form:aiΦi,j,k−1+ biΦi,j−1,k+ ciΦi−1,j,k+ diΦi,j,k+ eiΦi+1,j,k+ fiΦi,j+1,k+ giΦi,j,k+1= hi,where i = 1, ..., Nx, j = 1, ..., Ny, k = 1, ..., Nzandci= ei= 1/(∆x)2di= −2(1/∆x)2+ (1/∆y)2+ (1/∆z)2 bi= fi= 1/(∆y)2hi= 4πGρi,j,k(modulo BCs)ai= gi= 1/(∆z)24• Leads to very large sparse banded matrix:d e f gc · · · ·· · · · g· · · f· · · ·b · d e ·· c d · f· · · ·b · · ·a · · · ·· · · · ea b c d– Dimension is (NxNyNz) × (NxNyNz)!– =⇒ even very small problem (203) → large matrix 8000 × 8000.– “Reasonable” sized problem (1003) → 106× 106matrix!– Clearly need efficient ways to solve matrix:1. Relaxation schemes — guess solution, then relax (Cf. NRiC §19.5–19.6).E.g., “Successive Over-Relaxation” (SOR), “Alternating-Direction Implicit”(ADI), multi-grid (use exact solution on coarse grid as initial guess foriterative solution on fine grid), etc.2. Sparse banded solvers, e.g., conjugate gradient method (NRiC, §2.7).3. Fourier methods — solution of FDE in Fourier space is very simple, then caninverse Fourier t r ansform solution back to real space (NRiC §19.4).∗ Very powerful, but requires periodic BCs.Summary: PM Method• What is advantage of PM code?– Force solving scales a s O(Ng), where Ng= number of mesh grid points.– Leapfrog scales as O(Np), where Np= number of particles.– Work associated with leapfrog  solving Poisson’s equation.∴ can afford very large Np, e.g., Np106–8with Ng∼ 104–6.– Not good fo r correlated systems (in which 2-body encounters important) but greatfor uncorrelated systems (where it takes the place of