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.