This page mirrors docs/wiki/Multigrid-Self-Gravity.md in the repository, which is the authoritative version.

Adding a <gravity> block turns on Newtonian self-gravity. A geometric multigrid solver computes the gravitational potential from the gas density on the whole MeshBlock hierarchy. The force and its work are then added to the hydro or MHD update as a source term. The module works on uniform meshes, SMR and AMR, with periodic, isolated (multipole) or fixed-potential boundaries. For hydro runs it also works with sink particles and with local adaptive time stepping (LAT).

Use it for self-gravitating gas: a star that holds itself together, a collapsing cloud, or a stream that fragments. Do not use it for relativistic runs (it is Newtonian and has no GR version), for curvilinear coordinates (anisotropic Cartesian cells are supported), or on 1D and 2D meshes (the solver is 3D only). The cost is dominated by the Poisson solve, not by the coupling. Under LAT the potential is frozen over a window of fine steps and is solved only between windows.

An isolated problem (a star in vacuum) with multipole boundaries:

<gravity>
four_pi_G  = 1.0        # 4 pi G in code units, so G = 1/(4 pi)
mg_bc      = multipole  # isolated system
threshold  = 1.0e-6     # stop when the L2 defect is below this
niteration = 20         # but do at most 20 defect checks

Everything is set in <gravity>. There is no <multigrid> block, and a deck that has one aborts. Under LAT also set solve_dt (see LAT). The keys you will use most:

The solver finds the potential $\phi$ from the Poisson equation $\nabla^2\phi = 4\pi G\,\rho$, where $\rho$ is the gas density, and the gas feels the acceleration $-\nabla\phi$. The only number the solver reads is the product $4\pi G$ in your code units, which is four_pi_G. There is no separate input for $G$. For $G = 1$ set four_pi_G = 12.566370614359172. The value four_pi_G = 1.0 used in the TDE decks means $G = 1/(4\pi)$.

Two other places need $G$ itself, and both compute it as $G$ = four_pi_G/(4 pi). One is the analytic black-hole potential of the tde_external problem generator. The other is the sink particles, but only when the deck has no <units> block. With a <units> block the sinks take $G$ from the unit system instead, and nothing checks that this agrees with four_pi_G. If you use <units>, set four_pi_G to $4\pi$ times that $G$ yourself.

All keys are in <gravity> unless the table says otherwise. Defaults are the values in the source.

The defaults are right for most runs. These keys trade extra local work for fewer halo exchanges on several ranks.

The solver finds the potential $\phi$ that satisfies $\nabla^2\phi = 4\pi G\,\rho$. The density is the conserved density of whichever fluid is active, multiplied by the coupling weight of rho_grav_min.

fine grid   smooth --restrict--> smooth --restrict--> ... coarsest: solve
fine grid   smooth <-prolong---- smooth <-prolong---- ... (correct upwards)

• Multigrid. A smoothing sweep (red-black Gauss-Seidel with over-relaxation) removes short-wavelength error quickly but long-wavelength error slowly. So the solver restricts the residual to a coarser grid, where those long wavelengths are short, smooths there, and interpolates the correction back up. One trip down and up is a V-cycle. The coarsest grid is relaxed many times. A full-multigrid (FMG) start builds a good first guess by solving on the coarsest level and working up. Later solves start from the previous potential.
• Refinement. The levels between the root grid and the MeshBlocks that AMR or SMR creates are handled on the host as small 2x2x2 groups of cells (octets).
• Stopping. With a negative threshold the solver does exactly niteration V-cycles. With a positive one it checks the L2 defect after every defect_check_interval V-cycles and stops when it is below threshold, or after niteration checks (160 if niteration is not positive). With threshold = 0 an automatic controller runs a first burst of V-cycles (niteration of them if set, else fmg_ncycle + 1) and then extra cycles, until the defect is below $10^{-6}$ times the L2 norm of the source (never less than $10^{-12}$), until auto_max_extra_cycles extra cycles have run, or until the defect stops improving. The $10^{-6}$ is fixed in the code and is not an input key. A threshold relative to the source carries over between problems, while an absolute threshold depends on your units.
• Boundaries. zerograd is a zero-gradient condition and zerofixed sets the potential to zero on the face. multipole sets the potential on the face from a multipole expansion of the source (up to the quadrupole for mporder = 2, up to the hexadecapole for 4) about its centre of mass.
• Coupling. The force enters as a source term in the form of Mullen, Hanawa and Gammie (2020). The momentum changes by the potential gradient. The energy changes by the same potential differences multiplied by the Godunov mass flux through each cell face. Energy conservation is therefore limited by the multigrid residual and by the age of the potential. The dual-energy auxiliary variable is not changed by this term.
• When it solves. By default the potential is solved at every Runge-Kutta stage, because the stage density changes. With solve_every above 1 or with solve_dt, it is solved only on the last stage of a due cycle and reused in between. The first solve of a run, and the solves after a regrid, always run. With solve_dt the due times sit on the grid 0, solve_dt, 2 solve_dt, ... of code time, so a restarted run solves on the same cycles as an uninterrupted one.
• Under LAT. The potential is frozen inside a window, and no Poisson solve runs inside it, so the source density is never loaded from blocks at different local times. The driver solves before a window, and again at a window end when the run stops there (with lat_time_centered_work it solves at every window end). Before a window it first shortens the window to at most solve_dt. If the potential is invalid or due, or the window would run past the next due time, it solves first instead of cutting the window.

• Start with threshold = 0 (all the shipped TDE decks use it) for an automatic stopping rule. Use a fixed niteration (with threshold negative) for a reproducible cost, or a positive threshold for a controlled accuracy.
• Leave full_multigrid = true. It matters only for the first solve and for restarts.
• With SMR or AMR the defaults already switch to npostsmooth = 2 and tricubic prolongation.
• Leave the parallel tuning keys at their defaults unless you are tuning performance.
• To keep a vacuum background out of the gravity, set rho_grav_min above the background density. The TDE decks use hydro/dfloor = 1e-10 and rho_grav_min = 1e-9, so the coupling rises from 0 to 1 over that decade.

• Use mg_bc = multipole for an isolated system and a periodic mesh for a box.
• With all-Neumann boundaries (the default zerograd) the mean of the source is subtracted, so the equation solved is the Jeans swindle, not the isolated problem. The driver only warns.
• mporder must be 2 or 4. With auto_mporigin = false give the three mporigin_x* keys.

• Set gravity/solve_dt to the interval at which the potential should be refreshed. The potential is frozen inside a LAT window, so this key also caps the window length. A smaller value gives a fresher potential and shorter windows.
• LAT is available for hydro only in this release, so LAT with self-gravity is for hydro runs.
• Local Adaptive Time Stepping explains how a window works. Its list of supported and refused modules applies here too. The conditions that self-gravity adds are in Requirements and refusals.

Sink particles are deliberately absent from the Poisson source. Their potential enters through a separate source term that uses the same softened potential as the sink-to-gas reaction, so the force pair is equal and opposite cell by cell. <sink_particles>/create = true requires a <gravity> block, and sinks require <hydro>. mask_radius is the deck-level way to drop the Poisson source outside a radius.

The potential is not written to restart files, so the first solve after a restart is a cold FMG solve. A restarted run therefore agrees with the uninterrupted run only to the solver tolerance, not bit for bit. With a solve_every or solve_dt cadence and no LAT, the first stage after a restart also uses a freshly solved potential, where the uninterrupted run would still use an older one. Under LAT the potential is solved at the first window start. LAT for MHD is not in this release, and the LAT gravity logic reads the hydro source terms only.

Each of these makes the run exit with an error.

Mesh and setup:

• The mesh is not 3D, or the MeshBlocks are not logically cubic or not a power of two in size.
• The deck has a <multigrid> block, which is not a valid block name.
• There is neither a hydro nor an MHD module.
• <sink_particles>/create = true without a <gravity> block.
• hydro_srcterms/self_gravity = true without a <gravity> block.

Gravity keys:

• four_pi_G is not set, is 0, or is negative.
• Neither threshold nor niteration is set, or niteration = 0 with a negative threshold.
• auto_target_defect is set. It is no longer a key.
• mg_bc is periodic or an unknown word, or a face is periodic for gravity but not for the mesh.
• mporder is not 2 or 4, or both auto_mporigin and nodipole are true.
• mg_nghost is above 1 on a refined mesh.
• show_defect is not true, false or an integer, or solve_dt is negative.

Time stepping and LAT:

• time/lat with self-gravity and no positive gravity/solve_dt. The potential is frozen inside a window, and this key bounds the window length.
• time/lat_same_level = true with lat_neighbor_limiter = face, with self-gravity or with the analytic black hole. The limiter must be all or hybrid.
• time/lat with an <mhd> block, because LAT for MHD is not in this release. LAT also refuses other modules, such as an ion-neutral pack (see the LAT page linked above).
• lat_boundary_flux_diagnostics = true without lat_time_centered_work = true.
• lat_time_centered_work = true unless all of these hold: time/lat = true, rho_grav_min is at most hydro/dfloor, there is no positive mask_radius, there is no <sink_particles> block, and either solve_dt > 0 or solve_every = 1. At run time it also needs hydro with an energy equation (not isothermal), because the correction is added to the gas energy.

Two root-grid conditions only warn. A root grid that cannot be coarsened to one cell prints a warning, and its coarsest level is then solved iteratively, which is less efficient. A root grid with more than 100 coarsest cells prints a warning about the cost of the coarsest solve. Separately, auto_max_extra_cycles and warm_final_niter are silently ignored unless threshold = 0.

• Set show_defect = 1 for the final defect of each solve, or 2 for every iteration.
• show_timing = true prints the time per phase.
• The grav_phi output variable writes the gravitational potential (the solved potential plus the analytic black-hole potential if enabled). Before the first solve the solved part is zero. It is also added automatically to the hydro_u, hydro_w, mhd_u and mhd_w output bundles.
• The regression tests are in tst/test_suite/multigrid/. From tst/, run python run_test_suite.py --cpu --test test_suite/multigrid, or use --mpicpu or --gpu. They cover a two-sphere Poisson problem, Jeans waves and 64^3 hydro and MHD solves, with checks on convergence rate and on independence from MeshBlock size and rank count. Further gravity tests are in tst/test_suite/nr/.

• tde_external uses the solved potential to hold the star together, and applies the black hole as an analytic potential in the same source-term kernel (problem/external_bh_gravity_source).
• The decks that switch self-gravity on are the five inputs/TDE_examples/* decks and inputs/sink_particles/creation.athinput.

• Mullen, Hanawa and Gammie (2020): the flux-consistent gravitational-work source form.
• Local Adaptive Time Stepping: how a window works and how LAT and gravity stay synchronized.