1 / 26

Solving the Poisson Integral for the gravitational potential

Solving the Poisson Integral for the gravitational potential using the convolution theorem. Eduard Vorobyov Institute for Computational Astrophysics. Two avenues for finding the gravitational potential.

upton
Download Presentation

Solving the Poisson Integral for the gravitational potential

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  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 j=3 j=2 j=1 After discretizing, we obtain a set of 3 x 3 linear algebraic equations forunknown potentials (and other boundary values) need to be known from boundary conditions i=0 i=1 i=2 i=3 boundary layer x

  6. Periodic boundary conditions y j=3 j=2 j=1 i=0 i=1 i=2 i=3 boundary layer x

  7. Axis-of-symmetry or equator boundary conditions z Axis of symmetry j=3 j=2 j=1 r equator i=0 i=1 i=2 i=3 boundary layer

  8. Multipole expansion for axisymmetric mass distributions z Axis of symmetry boundary layer j=3 j=2 j=1 r i=1 i=2 i=3 i=4 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, Almust 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 < rB and in the case of Al – over grid zones with r > rB There is no telling how many terms in the above series will be needed!

  11. Axis of symmetry boundary layer z r j=3 j=2 j=1 rB i=1 i=2 i=3 i=4 If we do not take into account the input from grid zones with r > 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 < rB and more formulas for rB < 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 m=N-1 m=2 m=1 m=0 No boundary values involved in the summation! l=0 l=1 l=2 l=N-1 x

  14. Now let’s 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 2 1 3 4 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 simplebecause G ≠ 0

  17. Re-arranging the computational domain to make G periodic N-1 2 1 0 2 1 0 1 2 3 4 5 6 2N-1 -1 -2 -3 -N 3 4

  18. l=-N l=-3 l=-2 l=1 l=0 l=1 l=2 l=N-1 m=N-1 m=2 m=1 m=0 m=-1 m=-2 m=-3 m=-N Let’s assign to some arbitrary values along m=0 2 1 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 > 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, A&A, 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. reduced surface density 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!

More Related