Self Gravity in astrobear

The governing equation is Poisson's equation for gravity, which for a given mass distribution can be solved for the gravitational potential. In astrobear, we use the potential to solve for the gravitational forces in the fluid. The equation we use for the resultant force is either in conservative or non-conservative form (want to read more into the numerical methods here and link to these pages).

Laplace Operator

In cartesian coordinates, the Laplacian contains only simple derivatives, i.e. does not contain any functions of position as it does in cylindrical or spherical coordinates:

in Cartesian coordinates (x,y,z):

\nabla^2 = \frac{\partial^2f}{\partial x^2} + \frac{\partial^2f}{\partial y^2} + \frac{\partial^2f}{\partial z^2}

in Cylindrical coordinates (\rho, \phi, z):

\nabla^2 = \frac{1}{\rho} \frac{\partial}{\partial \rho}(\rho \frac{\partial f}{\partial \rho}) + \frac{1}{\rho ^2}\frac{\partial^2 f }{\partial \phi^2} + \frac{\partial^2f}{\partial z^2}

in Spherical (r, \theta, \phi) coordinates:

\nabla^2 = \frac{1}{r ^2} \frac{\partial}{\partial r}(r^2 \frac{\partial f}{\partial r}) + \frac{1}{r ^2 \sin \theta}\frac{\partial}{\partial \theta}(\sin \theta \frac{\partial f}{\partial \theta}) + \frac{1}{r^2 \sin^2 \theta} \frac{\partial^2 f}{\partial \phi ^2}

Currently astrobear is configured for gravity in Cartesian coordinates. I will be modifying the code so that it can solve for gravity in 2.5D (aka cylindrical, axisymmetric symmetry) and 1D spherical geometry (aka spherical coordinates with polar and azimuthal symmetry). With these symmetries, the Laplacian becomes:

in 2.5 d:

\nabla^2 = \frac{1}{\rho} \frac{\partial}{\partial \rho}(\rho \frac{\partial f}{\partial \rho}) + \frac{\partial^2f}{\partial z^2}

in 1d, spherical:

\nabla^2 = \frac{1}{r ^2} \frac{\partial}{\partial r}(r^2 \frac{\partial f}{\partial r})

The game then becomes expanding out these derivatives, discretizing the Laplacian, and putting the resulting system of poisson equations for the mesh in matrix form. This matrix will differ from the cartesian version in its coefficients. So, to modify the existing code, I will add different coefficients to the matrix for different geometries. These matrices are then sent off to hypre to be solved for phi.

The method in Cartesian coordinates

• What is our stencil? Eg. for 2D are a 5 point stencil, or a 9 point stencil.
• How do we order the matrix?
• What the matrix looks like

Expanding out derivatives for 2.5d and 1d spherical

Modified matrices

Relevent regions of the code

Tests for 2.5d and 1d gravity

Downstream algoritms, gravitational force solvers

Erica's latest updates:

Poisson Init routine

The first chunk of code in poisson.f90, 'PoissonInit', goes through and checks that the boundary conditions are set up correctly for the elliptic solver, and then sets up the 'Hypre' arrays (those arrays which are sent off to Hypre for solving Poisson's equation). In what follows, we'll focus on the ndim=2 case for illustration of the code.

There are 2 main arrays in this section of the code: 1. 'offsets', and 2. 'stencil values'. (ientries seems less important).

'Offsets' is an (ndim x 0:2*ndim) matrix; for 2D, offsets is (2x5). Let's call these indices (i,m). The purpose of offsets is to set the stencil geometry for Hypre. Recall, each stencil is comprised of cells, and each cell has a stencil. The first index (i) of the matrix denotes the dimension of the cell in the stencil. The second (m) indexes the cell in the stencil. m=0 denotes the center cell of the stencil, m=1 the 1st cell in the stencil, m=2 the second, and so on until you reach the final number of cells in the stencil (for the 2D case this is m=4 as it is a 5-point stencil and m starts at 0).

Offsets (:,0) = 0 is the starting line in the routine. It states that in any dimension, the center cell itself has a '0 offset'.

When all is said and done, the array's values are given as:

Offsets (i,m) =

i\m	0	1	2	3	4
1	0	-1	1	0	0
2	0	0	0	-1	1

-1 corresponds to left/down, and +1 to right/up. So this states the 1st cell in the stencil (m=1) has a position to the left of the center cell in the i=1 (x) direction, and the 2nd cell is to the right. The 3rd and 4th cell has a 0 offset in x (because they are above and below the center cell). This is why their offsets are listed in the 2nd row of the table. Again, the center cell itself (m=0) has 0 offsets from itself in both directions.

The stencil values are also set here for regular, Cartesian geometry. They need to be modified later in the code for use with 2.5 geometry.

Stencil values is an (0:2*ndim) array. For 2D then, stencil values is 5 elements long. It is initialized with the following values which are the coefficients for the Poisson matrix equation described above for Cartesian geometry:

Stencil values (m) =

m	0	1	2	3	4	5
stencil values(m)	-4	1	1	1	1	1

Together with the offsets array, this means that the values of stencil values corresponds to the cells:

*StencilValues = (center, left, right, down, up)*