Solving the Poisson Integral

1
Solving the Poisson Integral
for the gravitational potential
using the convolution theorem
Eduard Vorobyov Institute for
Computational Astrophysics
2

• Two avenues for finding the gravitational
  potential.
• Solution techniques for the Poisson equation.
• The problem of boundary conditions.
• The Poisson integral for the gravitational
  potential.
• The convolution theorem.
• Applications of the convolution theorem for
  solving
• the Poisson integral.
• FFT libraries.

3
Although gravity is omnipresent in the Universe,
its effect is often simplified or even neglected.
However, there are situations when an accurate
calculation of gravity is a necessity (galaxy
formation and evolution, star and planet
formation, circumstellar disk dynamics, etc.)
Gravitational potential
Gravity force per unit mass
How to calculate the gravitational potential?
Solving for the Poisson equation (grid-based
hydrocodes)
Solving for the Poisson integral (SPH and
N-body codes)
4
Solution procedure for the Poisson equation
Discretization (equidistant or non-equidistant
grid)

Boundary conditions (periodic or non-periodic)

• Equidistant grid and periodic BC
• FFT

• Non-equidistant grid or non-periodic BC
• ADI method (3D with axial symmetry)
• SOR (slow in full 3D)
• Multigrid methods (fast only on Cartesian
• geometry?)

5
The Poisson equation, when discretized, needs a
knowledge of the boundary values!
For a simple case of 2D Cartesian equidistant
mesh we obtain a five-zone molecule
y






j3 j2 j1
After discretizing, we obtain a set of 3 x 3
linear algebraic equations for unknown
potentials
(and other boundary values)
need to be known from boundary conditions
i0 i1 i2 i3
boundary layer
x
6
Periodic boundary conditions
y






j3 j2 j1
i0 i1 i2 i3
boundary layer
x
7
Axis-of-symmetry or equator boundary conditions
z
Axis of symmetry






j3 j2 j1
r
equator
i0 i1 i2 i3
boundary layer
8
Multipole expansion for axisymmetric mass
distributions
z
Axis of symmetry
boundary layer






j3 j2 j1
r
i1 i2 i3 i4
Laplace equation in spherical coordinates (r, q)
9
The method of separation of variables (Jackson
1975)
-- Legendre polynomials
(can be pre-computed and stored)
So far, we have not specified the location of our
boundary with respect
to the computational domain
In the case of the outer boundary, when ALL mass
is confined within radius rB, Al must go to zero
for the potential to have a finite value at rB?
inf In the opposite case of the inner boundary,
Bl 0.
10
Bl and Al are the so-called interior and exterior
multipole moments
In the case of Bl, the integration (summation) is
performed over ALL grid zones with r lt rB and in
the case of Al over grid zones with r gt rB
There is no telling how many terms in the above
series will be needed!
11
Axis of symmetry
boundary layer
z








r
j3 j2 j1
rB
i1 i2 i3 i4
If we do not take into account the input from
grid zones with r gt rB, the series may diverge!
12
Multipole expansion for non-axisymmetric mass
distributions
Ylm are the spherical harmonic functions
(array of 4 variables!) and Blm are multipole
moments
The integration (summation) is performed over
grid zones with r lt rB
and more formulas for rB lt r ..
The fully 3D case is a lot more complicated than
2.5D case and it takes substantial
computational resources .
See Cohl Tohline (ApJ 1999) Binney Tremaine,
Galactic Dynamics
13
Finding the gravitational potential using the
Poisson integral
For a simple 2D Cartesian grid
M(xl,ym) is the mass contained in
grid zone (l,m)
y




mN-1 m2 m1 m0
No boundary values involved in the summation!
l0 l1 l2 lN-1
x
14
Now lets assume that our computational grid is
equidistant.
Gravitational potential in zone (l,m) created by
unit mass located in zone
Direct summation takes N2 operations, where N is
the total number of grid zones
A much faster way for evaluating the double sum
is to use the convolution theorem
15
The convolution theorem
B and C are periodic with a period of 2N
This sum can be calculated using the following
three steps
direct Fourier transform
product of Fourier transforms
inverse Fourier Transform
16
Convolution sum
Our gravitational potential
Doubling the computational domain
M and G are in general non-periodic and we have
to make them periodic








N-1 2 1 0








-1 -2 -3 -N
We may require M be periodic with a period of 2N
because M 0 in zones 2,3,4. With G it is not
that simple because G ? 0
17
Re-arranging the computational domain to make G
periodic








N-1 2 1 0
1
2








-1 -2 -3 -N
0 1 2 3 4
5 6 2N-1
3
4
18
Lets assign to some arbitrary
values along m0








1
2
40 50 60 70
80 90 100 90








3
4
19
Singularity of the Green function
However, it is possible to calculate the
contribution of the material in the (l,m)th cell
to the potential in the same cell by assuming
constant surface density within the cell and
integrating over the cell area
within an individual cell (l,m)
20
defining
and noticing that
after a few pages of algebra
21
3D Cartesian coordinates
The extension to 3D Cartesian coordinates is
straightforward
has to be taken numerically .
Problem the convolution method takes a lot of
memory in 3D due to doubling of the
computational grid. Some remedy see Hockney and
Eastwood, Computer simulations using particles.
The Fourier transform of the Green function has
to be taken only once if the grid is not
arbitrarily varying during simulations. This
leaves us with 2 FFTs each taking 2 N log2N
operations where N is the total number of grid
zones. The direct summation takes N2 operations
and we have a speedup for N gt 16.
22
3D cylindrical coordinates
If j and z coordinates are discretized evenly,
the sums over and are a convolution,
but the sum over is not, irrespective of
the discretization!
The procedure is to rearrange the triple sum and
take the inner two sums for each and every
cylindrical layer using the convolution
theorem (thus finding the gravitational potential
of the layer) and then perform a direct summation
over all cylindrical layers.
constants
Slower, but takes less memory since doubling is
needed in z-direction only. See more details in

Pfenniger Friedly, AA, 1993
23
2D polar coordinates
Logarithmically spaced grid in r-direction
Simulations of galactic and stellar disk dynamics
require high resolution in the inner regions,
while a lower resolution may be
sufficient in the outer regions
24
We introduce a new radial coordinate
reduced potential
See more in Binney Tremaine, Galactic Dynamics,
pp 96-97.
25
FFT libraries

• ACML (AMD architecture, OpenMP parallelized,
  free)
• ICML (Intel architecture, commercial)
• MKL (Intel architecture, commercial)
• FFTW (MPI parallelized, OpenMP parallelized?,
  free)

26
Thank you!