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).

Description of why we want 2.5D self gravity

Big picture of Hypre /Astrobear connection

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)

This means that we need to modify the elements StencilValues(1) and StencilValues(2) for Cylindrical self-gravity, as ticket 150 shows. This means setting,

StencilValues(1)=1- \frac{\Delta r}{2r}

StencilValues(2)=1+ \frac{\Delta r}{2r}

Poisson Matrix Set Box Values routine

This routine is responsible for making changes to the stencil coefficients, for say, ghost zones, or differing geometry. This is where the mods to the stencil for 2.5D geometry will take place.

Since stencil values is a module wide variable, we should not modify it itself in any routine. Instead, we define a new array, tempStencilValues, set it = StencilValues, and then modify elements of this new array when necessary.

The routine begins by allocating an array to hold all of the coefficients for all of the stencil's cells for each of the cells in the patch. This array is aptly called 'matrixvalues' and is mx*my*mz*num_stencil_values long (where mx,y,z are the number of cells in each direction of the patch, and num_stencil_values is self explanatory).

The main body of the routine loops through all of the cells in the mesh, starting with each z, y, and then x. Since we are only concerned with modifying the 2D geometry for cylindrical, we really are looping through y first and then x. This means that for each y value (can think of it as a 'row' in the 2D mesh), we loop through each cell in x on that row, and calculate any modifications to the stencil that need to happen. After these mods are made, the 5 (in this special case) stencil values are fed into the 5 matrix values slots for that particular cell of the grid. And the loops continues. These values are then sent off to Hypre for solution.

2.5D Mod

What we then need to add is a modification in these loops for cylindrical. As stated above, this modification is going to need both the dx of the patch, and the x position of the cell (also patch dependent - because x(dx)). Instead of calculating for each iteration of the loop, we can try and save some computational power and put the calculation of x and dx before the loop.

Dx for the patch is easy, that is just,

dx = levels(Info\%level)\% dx

To get all of the x_i's of the cells in the patch we can allocate an array 'x_i' that has length = number of cells along x axis of patch,

Allocate ~x(Info\%mGlobal(1,1): Info\%mGlobal(1,2))

(For a refresher page describing patch specific quantities, refer here)

Now to populate this array, we loop through:

Do ~i = Info\%mGlobal(1,1), Info\%mGlobal(1,2)

~x(i) = GxBounds(1,1) + (i - half)*dx

End ~Do

This starts at the global lower left hand corner of the entire box (gXbounds(1,1)) and marches over i*dx_patch spaces (subtracting off the .5dx to get to the cell center).

Another way to do this would have been to use the lower x position of the patch itself, say:

Once we have these 2 things calculated, we can then add in the meat and potatoes:

If ~(iCylindrical > 0)~Then

~tempStencilValues(1) = 1-dx/2x(i)
tempStencilValues(2) = 1+dx/2x(i)

End~ If

Then set the matrixvalues for the stencil of the ith cell, by doing

matrixvalues(m:m+4) = tempStencilValues(1:5)

where m increases by 5 for each cell in the mesh. For example,

matrixvalues(:)=(-4,1,1,1,1, -4, 1, 1, 1, 1, -4, 1, 1, 1, 1, ...., -4, 1, 1, 1, 1)

for regular 2D Cartesian geometry.

This means that matrixvalues(1) is the center cell stencil coefficient of the 1st cell in the mesh, matrix(6) is the center cell stencil coefficient of the 2nd cell in the mesh, and so on.

Boundary Conditions

Vector Fixes

Force Fixes