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, with appropriate boundary conditions that closes the system of discretized equations on a mesh. These boundary conditions are either periodic, or specifying the potential or its derivative on the boundary. In astrobear, we use the potential to solve for the gravitational forces in the fluid.

Multidimensions

Poisson's equation can be simplified (i.e. written with last terms) under various symmetries in the problem setup. Some examples are:

2D-symmetry— Can drop the z-derivative term in the PDE. This modified equation is appropriate for an infinite slab. The variables in the PDE are x & y.

2.5D-symmetry— Use cylindrical geometry to rewrite the PDE, and drop the theta derivative terms. This is appropriate for cylindrical symmetry. The variables are r & z.

If the simulations have these symmetries, then instead of solving Poisson's full 3D equation, we can solve a modified, simpler equation.

2.5D Geometry

2.5D means cylindrical symmetry, i.e. the xy plane now represents the rz plane — there is azimuthal symmetry in phi.

The hyperbolic solvers for the Euler equations have been modified (i.e. include 1/r terms now), however, up until this fix in the code, any calculations with gravity used regular 2D gravity, which would introduce errors. Here is an example:

Results

Check out my astrobear blog to see results.

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, since cylindrical fixes amount to modifying the 2D stencil coefficients.

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 (and 2.5D) 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
stencil values(m)	-4	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. Originally it was written to fix coefficients for boundary/ghost zones. This makes it a natural place to modify coefficients for differing geometry as well; this is where the fixes for 2.5D gravity will go.

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 patch, 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 a 2D matrix), 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 the case of 2D) stencil values are fed into the 5 matrix value slots for that particular cell of the grid. And the loops continues. These values are 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 positions 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 parameters calculated, we can then add in the fix for 2.5D:

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 mth cell, by doing

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

Note, for 2/2.5D, 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

Testing

Added a print statement for tempStencilValues