Download Particle Mesh Method: Solving Poisson's Equation in Gravitational Systems - Prof. John Wal and more Study notes Computer Science in PDF only on Docsity!
CSI 761- Fall 2008 - Lecture 4 Notes
September 22, 2008
1 Summary
Syllabus topics:
- Finite Difference Methods
- Boundary Conditions
- Linear Solvers
Educational Objectives: At the end of the lecture and homework, students will be able to
- Understand the connections between the integral and differential forms used to solve gravita- tional problems
- Be able to put Poisson’s equation into a finite difference form
Lesson Plan:
- review homework assignments
- talk about octave implementation
2 Gravity
Some useful constants
- G = 6. 67 × 10 −^8 cm^3 gm−^1 s−^2
- M⊙ = 2 × 1033 gm
- 1 AU = 1. 5 × 1013 cm
- 1 year = 3. 1 × 108 seconds
- 1 pc = 3 × 1018 cm
F^ ~ = −GM m r^2
rˆ (1)
U = −
GM m r
F~ = m~a (3)
~a =
d~v dt
~v =
d~r dt
~r = (x, y, z) (6)
3 An Overview of the Particle Mesh Method
3.1 Introduction
The particle mesh method works allowing particles to freely move, but calculating the potentials and the forces on them using a grid. In a direct particle-particle code, we need to calculate the forces on each particle using
F^ ~i =
∑^ N
j=1,j 6 =i
GMiMj r ij^2
r ˆij
where N is the total number of particles in the system. For a particle mesh code, in principle we have to do a calculation for each grid cell k
F~k =
∑^ M
l=1l 6 =k
GMkMl r kl^2
ˆrkl
Where M is the total number of cells in the grid. The real advantage in the particle mesh method comes in for several reasons-
- In general, there are many particles in a single cell. Thus, there are fewer calculations to complete compared to a particle-particle code.
- Because we are using a grid, there are rapid potential solvers that can be used to decrease the computational time to do the potential and force calculations.
The second reason, namely the efficiencies of using a regular grid for the calculation can be substantial.
where δ is the Dirac delta function such that
δ(x) =
{ 0 , if x 6 = 0 ∞, if x = 0
and
∫ (^) ∞
−∞
δ(x)dx = 1
5 Gravity
Given a vector field of the form
F^ ~ = −GM m r^2
rˆ
We can use the divergence theorem ∫
V
( ∇ · F~
) dτ =
∮
S
F^ ~ · d~a ∫
V
( ∇ · F~
) dτ = −
∮
S
GM m r^2
rˆ · d~a
Since
d~a = r^2 sin θdθ dφ rˆ
On the right side, we have
∮
S
GM m r^2
ˆr · d~a = −GM m
∮
S
ˆr r^2
· d~a
= −GM m
∮
S
ˆr r^2
r^2 · ˆr sin θdθ dφ
= −GM m
∮
S
sin θdθ dφ = − 4 πGM m
The mass of our test particle m is fixed, as is G. However, the total mass M is a function of the density ρ(~r) and the volume
M =
∫
V
ρdτ
Thus, we can write ∫
V
( ∇ · F~
) dτ = − 4 πGm
∫
V
ρdτ
Remembering our definitions for gravitational potential
F^ ~ = −∇U
∫
V
(∇ · ∇U ) dτ = − 4 πGm
∫
V
ρdτ ∫
V
( ∇^2 U
) dτ = 4 πGm
∫
V
ρdτ
Simplifying this, we have Poisson’s equation for gravitational systems
∇^2 U = 4 πGm ρ(~r)
Note that this is identical in form to the electromagnetic Poisson equation of
∇^2 U = −
ǫ 0
ρ(~r)
Instead of calculating forces or potentials using integration, we can instead solve this as a system of partial differential equations.
5.1 Poission equation check
Let
ρ(~r) = M δ(~r)
∇^2 U = 4 πGm M δ(r)
Integrating this over a volume V ∫
V
∇^2 U dτ = 4 πGm M
∫
V
δ(r)dτ ∫
V
∇^2 U dτ = 4 πGm M
Using the spherical form of the Laplacian
∇^2 U =
r^2
∂r
( r^2
∂U
∂r
)
∫
V
r^2
∂r
( r^2
∂V
∂r
) dτ = 4 πGm M
The spherical volume integration element is
dτ = r^2 dr sinθdθ dφ
or for spherically symmetric systems
dτ = 4πr^2 dr
Let’s examine the one dimensional case with cell spacing h. The cell centers are designated as
{· · · , xi− 2 , xi− 1 , xi, xi+1, xi+2, · · ·}
Let’s assume the particle is between
xi ≥ x ≥ xi+ Then the mass will be divided between cell i and i + 1 αi = (x − xi)
ρi =
M
h
αi
ρi+1 =
M
h
(1 − αi)
Let’s assume we have a two dimensional problem with cell centers set as {· · · , xi− 2 , xi− 1 , xi, xi+1, xi+2, · · ·} {· · · , yj− 2 , yj− 1 , yj , yj+1, yj+2, · · ·}
where the spacing is h and k in the x and y dimensions respectively. In two dimensions, we can generalize this to αi = (x − xi) αj = (y − yj ) Then we can write
ρi,j =
M
hk
αiαj
ρi+1,j =
M
hk
(1 − αi) αj
ρi,j+1 =
M
hk
αi (1 − αj )
ρi+1,j+1 =
M
hk
(1 − αi) (1 − αj )
6.2 Interpolating Forces from the Grid
Once we solve for the potential on the grid, we can then solve for the forces using the finite difference approximation for the derivative
F^ ~ = −∇U
F~ = −
[ ∂U ∂x
xˆ −
∂U
∂y
yˆ
]
Using the center difference approximation
F^ ~ij = −
[( Ui+1,j − Ui− 1 ,j h
) x ˆ +
( Ui,j+1 − Ui,j− 1 k
) y ˆ
]
This gives the force at the center of the grid cell at i, j. We still need to interpolate the force back to the particles.
6.2.1 Nearest Grid Point Interpolation
Just as in the particle mapping, the nearest grid point interpolation uses the center of the cell as the reference for the force. In short, for particle k located in the grid cell centered on i, j, we have
F^ ~k = F~i,j
6.3 Cloud in Cell Interpolation
As with the cloud in the cell particle mapping, we can use linear interpolation to calculate the forces on the particles. Assuming the particle k is cell i, j and x ≥ xi and y ≥ yj We use the definitions of
αi = (x − xi) αj = (y − yj )
We use the same geometric weighting that was used to distribute the particle’s mass to interpolate the forces back to the particle. We can define some constants-
k 1 = (1 − αi)(1 − αj ) k 2 = αi(1 − αj ) k 3 = (1 − αi)αj k 4 = αiαj
Then we can write
F^ ~ = k 1 F~i,j + k 2 F~i+1,j + k 3 F~i,j+1 + k 4 F~i+1,j+
7.1.1 Characteristic Equations
We classify PDE’s on the basis of their characteristic equations. Consider
auxx + buxy + cuyy = f
Under what conditions does u, ux, and uy uniquely determine uxx, uxy, uyy such that our equation is satisfied? To understand this, we need to add two identities using the definition of derivatives of u with respect to x and y.
d(ux) = uxxdx + uxydy d(uy) = uxydx + uyydy
We examine the case and assume that multiplicity of solutions are possible. Using the previous three equations in matrix form
a b c dx dy 0 0 dx dy
uxx uxy uyy
=
f d(ux) d(uy)
The solution for the PDE exists and is unique unless the determinant of the coefficent matrix vanishes. If multiple solutions do exist, the determinant is zero. The determinant of
a b c dx dy 0 0 dx dy
is
a(dy)^2 − bdydx + c(dx)^2 = 0
This is clearly a quadratic equation. The behavior of this quadratic is closely related to the character of the the PDE.
- b^2 − 4 ac > 0 - equation is hyperbolic
- b^2 − 4 ac = 0 - equation is parabolic
- b^2 − 4 ac < 0 - equation is elliptic
This type of analysis is a bit tricky, but can be done for any PDE.
7.1.2 Classification of PDEs
Normally, we consider differential equations to have three distinct classes based on their mathematical behavior.
∂u ∂t
= α
∂^2 u ∂x^2
∂^2 u ∂y^2
∂^2 u ∂z^2
α
∂^2 u ∂x^2
∂^2 u ∂y^2
∂^2 u ∂z^2
= g
∂^2 u ∂t^2
= α
∂^2 u ∂x^2
∂^2 u ∂y^2
∂^2 u ∂z^2
These classifications come from the analysis of the characteristic equations associated with each type of PDE.
7.1.3 Well Posed Problems
A problem is “well posed” if it
- a solution exists
- the solution is unique
- the unique solution depends on auxiliary data
Problems which are not well posed are difficult to solve, but may be of interest in some cases.
7.1.4 Separation of Variables
The basic way we approach analytic PDE solutions is separation of variables. In short, we assume the partial derivatives are independent of each other. The heat equation is:
∂u ∂t
= β^2
∂^2 u ∂x^2 If we assume a solution of the form
u(x, t) = X(x)T (t)
then
β^2 X′′T = XT ′
which gives us
X′′ X
β^2
T ′
T
= σ
where σ is the separation constant From
X′′ X
β^2
T ′
T
= σ
8 Finding the Forces on the Grid
In the Particle-Mesh method, you calculate the potentials on the grid. However, we need to find the forces so we can move the particles. To find the forces, we use a centered finite difference scheme.
F^ ~ = −∇U
F~ = −
( ∂U ∂x
ˆx +
∂U
∂y
yˆ
)
Using the centered difference scheme described earlier, we can write
F^ ~ = −∇U
F~ = −
([ U (xi+1,j ) − U (xi− 1 ,j ) 2 h
] ˆx +
[ U (xi,j+1) − U (xi,j− 1 ) 2 k
] y ˆ
)
where h is the spacing the x direction and k is the cell spacing in the y direction. It is often simplier to write this as
F^ ~ = −
([ U
i+1,j −^ Ui− 1 ,j 2 h
] ˆx +
[ U
i,j+1 −^ Ui,j− 1 2 k
] y ˆ
)
Where we have used
Ui,j = U (xi,j )
