Self Gravity with FFT - PrincetonUniversity/athena GitHub Wiki

A new Gravity class is constructed as a member of MeshBlock. The Gravity class contains the gravitational potential array phi. For the FFT gravity solver, the boundary condition is also separately set for phi.

When self-gravity is enabled, phi is automatically included in any Outputs that include either the prim or cons variables. Alternatively, it can be explicitly included in the variable field as phi.

A new FFTGravityDriver class is constructed and solves Poisson's equation using FFT assuming periodic BCs in all directions. The FFTGravityDriver::Solve() function is called in the main loop at every time-integrator stage, and it performs the following steps:

1. load density from MeshBlocks to the FFT input array
2. execute forward FFT
3. multiply kernel (-k^-2; SetNormFactor() takes the four_pi_G factor into account)
4. execute backward FFT
5. retrieve the real part of the FFT output, and store in the gravitational potential array phi
6. call gravity tasklist to send/receive boundary values

The gravitational constant (four_pi_G) should be set in Mesh::InitMeshData using SetGravityConstant or SetFourPiG function.

The background mean density (grav_mean_rho) should be set in Mesh::InitMeshData using SetMeanDensity function. For problems with the Jeans' swindle (= fully periodic BCs), this should be the initial mean density (see athena/src/pgen/jeans.cpp; see also #155). For problems with other BCs (not implemented yet), this should be zero.

Since we have switched back to the source term implementation, we do not need to set the mean density any more. A pgen file using SetMeanDensity has to be updated as the function is deleted.

> cd tst/regression/
> ./run_tests.py grav

> ./configure.py --prob=poisson -fft --grav=fft
> ./athena -i ../inputs/hydro/athinput.poisson 

Updates to momenta and total energy due to self-gravity are added via source terms, not via the divergence of a gravitational stress tensor or gravitational energy "flux", respectively. Source terms are added following the prescription detailed in MHG2020, with the exception that the energy source term does not use the time-averaged gravitational potential. With this simplification, total energy is not conserved to round-off error. Momentum conservation is determined by the accuracy of the gravitational potential.