Gravity Unit
Introduction
The Gravity
unit supplied with FlashX computes gravitational source
terms for the code. These source terms can take the form of the
gravitational potential \(\phi({\bf x})\) or the gravitational
acceleration \({\bf g}({\bf x})\),
The gravitational field can be externally imposed or selfconsistently computed from the gas density via the Poisson equation,
where \(G\) is Newton’s gravitational constant. In the latter case, either periodic or isolated boundary conditions can be applied.
Externally Applied Fields
The FlashX distribution includes three externally applied gravitational
fields, along with a placeholder module for you to create your own. Each
provides the acceleration vector \({\bf g}({\bf x})\) directly,
without using the gravitational potential \(\phi({\bf x})\) (with
the exception of UserDefined
, see below).
When building an application that uses an external, timeindependent
Gravity
implementation, no additional storage in unk
for holding
gravitational potential or accelerations is needed or defined.
Constant Gravitational Field
This implementation creates a spatially and temporally constant field
parallel to one of the coordinate axes. The magnitude and direction of
the field can be set at runtime. This unit is called
Gravity/GravityMain/Constant
.
Planeparallel Gravitational field
This PlanePar
version implements a timeconstant gravitational field
that is parallel to one of the coordinate axes and falls off with the
square of the distance from a fixed location. The field is assumed to be
generated by a point mass or by a spherically symmetric mass
distribution. A finite softening length may optionally be applied.
This type of gravitational field is useful when the computational domain
is large enough in the direction radial to the field source that the
field is not approximately constant, but the domain’s dimension
perpendicular to the radial direction is small compared to the distance
to the source. In this case the angular variation of the field direction
may be ignored. The PlanePar
field is cheaper to compute than the
PointMass
field described below, since no fractional powers of the
distance are required. The acceleration vector is parallel to one of the
coordinate axes, and its magnitude drops off with distance along that
axis as the inverse distance squared. Its magnitude and direction are
independent of the other two coordinates.
Gravitational Field of a Point Mass
This PointMass
implementation describes the gravitational field due
to a point mass at a fixed location. A finite softening length may
optionally be applied. The acceleration falls off with the square of the
distance from a given point. The acceleration vector is everywhere
directed toward this point.
UserDefined Gravitational Field
The UserDefined
implementation is a placeholder module for the user
to create their own external gravitational field. All of the subroutines
in this module are stubs, and the user may copy these stubs to their
setup directory to write their own implementation, either by specifying
the gravitational acceleration directly or by specifying the
gravitational potential and taking its gradient. If your userdefined
gravitational field is timevarying, you may also want to set
PPDEFINE flashx_GRAVITY_TIMEDEP
in your setup’s Config file.
Selfgravity
The selfconsistent gravity algorithm supplied with FlashX computes the Newtonian gravitational field produced by the matter. The produced potential function satisfies Poisson’s equation [Eqn:Poisson]. This unit’s implementation can also return the acceleration field \({\bf g}({\bf x})\) computed by finitedifferencing the potential using the expressions
In order to preserve the secondorder accuracy of these expressions at jumps in grid refinement, it is important to use quadratic interpolants when filling guard cells at such locations. Otherwise, the truncation error of the interpolants will produce unphysical forces at these block boundaries.
Two algorithms are available for solving the Poisson equations:
Gravity/GravityMain/Multipole
and Gravity/GravityMain/Multigrid
.
The initialization routines for these algorithms are contained in the
Gravity
unit, but the actual implementations are contained below the
Grid
unit due to code architecture constraints.
The multipolebased solver described in for self gravity is appropriate for spherical or nearlyspherical mass distributions with isolated boundary conditions. For nonspherical mass distributions higher order moments of the solver must be used. Note that storage and CPU costs scale roughly as the square of number of moments used, so it is best to use this solver only for nearly spherical matter distributions.
The multigrid solver described in is appropriate for general mass distributions and can solve problems with more general boundary conditions.
The tree solver described in [Sec:GridSolversBHTree] is appropriate for general mass distributions and can solve problems with both isolated and periodic boundary conditions set independently in individual directions.
Coupling Gravity with Hydrodynamics
The gravitational field couples to the Euler equations only through the momentum and energy equations. If we define the total energy density as
where \(\epsilon\) is the specific internal energy, then the gravitational source terms for the momentum and energy equations are \(\rho{\bf g}\) and \(\rho{\bf v}\cdot{\bf g}\), respectively. Because of the variety of ways in which different hydrodynamics schemes treat these source terms, the gravity module only supplies the potential \(\phi\) and acceleration \({\bf g}\), leaving the implementation of the fluid coupling to the hydrodynamics module. Finitedifference and finitevolume hydrodynamic schemes apply the source terms in their advection steps, sometimes at multiple intermediate timesteps and sometimes using staggered meshes for vector quantities like \({\bf v}\) and \({\bf g}\).
For example, the PPM algorithm supplied with FlashX uses the following update steps to obtain the momentum and energy in cell \(i\) at timestep \(n+1\)
Here \(g_i^{n+1}\) is obtained by extrapolation from \(\phi_i^{n1}\) and \(\phi_i^n\). The \({\tt Poisson}\) gravity implementation supplies a mesh variable to contain the potential from the previous timestep; future releases of FlashX may permit the storage of several time levels of this quantity for hydrodynamics algorithms that require more steps. Currently, \({\bf g}\) is computed at cell centers.
Note that finitevolume schemes do not retain explicit conservation of momentum and energy when gravity source terms are added. Godunov schemes such as PPM, require an additional step in order to preserve secondorder time accuracy. The gravitational acceleration component \(g_i\) is fitted by interpolants along with the other state variables, and these interpolants are used to construct characteristicaveraged values of \({\bf g}\) in each cell. The velocity states \(v_{L,i+1/2}\) and \(v_{R,i+1/2}\), which are used as inputs to the Riemann problem solver, are then corrected to account for the acceleration using the following expressions
Here \(g^\pm_{X,i+1/2}\) is the acceleration averaged using the interpolant on the \(X\) side of the interface (\(X=L,R\)) for \(v\pm c\) characteristics, which bring material to the interface between cells \(i\) and \(i+1\) during the timestep.
Tree Gravity
The Tree
implementation of the gravity unit in
physics/Gravity/GravityMain/Poisson/BHTree
is meant to be used
together with the tree solver implementation
Grid/GridSolvers/BHTree/Wunsch
. It either calculates the
gravitational potential field which is subsequently differentiated in
subroutine physics/Gravity/Gravity_accelOneRow
to obtain the
gravitational acceleration, or it calculates the gravitational
acceleration directly. The latter approach is more accurate, because the
error due to numerical differentiation is avoided, however, it consumes
more memory for storing three components of the gravitational
acceleration. The direct acceleration calculation can be switched on by
specifying bhtreeAcc=1
as a command line argument of the setup
script.
The gravity unit provides subroutines for building and walking the tree called by the tree solver. In this version, only monopole moments (node masses) are used for the potential/acceleration calculation. It also defines new multipole acceptance criteria (MACs) that estimate the error in gravitational acceleration of a contribution of a single node to the potential (hereafter partial error) much better than purely geometrical MAC defined in the tree solver. They are: (1) the approximate partial error (APE), and (2) the maximum partial error (MPE). The first one is based on an assumption that the partial error is proportional to the multipole moment of the node. The node is accepted for calculation of the potential if
where \(D\) is distance between the pointofcalculation and the
node, \(M\) is the node mass, \(S_\mathrm{node}\) is the node
size, \(m\) is a degree of the multipole approximation and
\(\Delta a_\mathrm{p,APE}\) is the requested maximum error in
acceleration (controlled by runtime parameter Gravity/grv_bhAccErr
).
Since only monopole moments are used for the potential calculation, the
most reasonable choice of \(m\) seems to be \(m=2\). This MAC is
similar to the one used in Gadget2 (see Springel, 2005, MNRAS, 364,
1105).
The second MAC (maximum partial error, MPE) calculates the error in acceleration of a single node contribution \(\Delta a_\mathrm{p,MPE}\) according to formula 9 from Salmon&Warren (1994; see this paper for details):
where \(B_n = \Sigma_i m_i \mathbf{r}_i\mathbf{r}_0^n\) where
\(m_i\) and \(\mathbf{r}_i\) are masses and positions of
individual grid cells within the node and \(\mathbf{r}_0\) is the
node mass center position. Moment \(B_2\) can be easily determined
during the tree build, moment \(B_3\) can be estimated as
\(B_3^2 \ge
B_2^3/M\). The maximum allowed partial error in gravitational
acceleration is controlled by runtime parameters
Gravity/grv_bhAccErr
and Gravity/grv_bhUseRelAccErr
(see
1.4.1).
During the tree walk, subroutine
physics/Gravity/Gravity_bhNodeContrib
adds contributions of tree
nodes to the gravitational potential or acceleration fields. In case of
the potential, the contribution is
if isolated
boundary conditions are used, or
if periodic periodic
or mixed
boundary conditions are used. In
case of the acceleration, the contributions are
for isolated
boundary conditions, or
for periodic
or mixed
boundary conditions. In In the above
formulae, \(G\) is the constant of gravity, \(M\) is the node
mass, \(\vec{r}\) is the position vector between
pointofcalculation and the node mass center and
\(f_\mathrm{EF,\Phi}\) and \(f_\mathrm{EF,a}\) are the Ewald
fields for the potential and the acceleration (see below).
Boundary conditions can be isolated or periodic, set independently for each direction. If they are periodic at least in one direction, the Ewald method is used for the potential calculation (Ewald, P. P., 1921, Ann. Phys. 64, 253). The original Ewald method is an efficient method for computing gravitational field for problems with periodic boundary conditions in three directions. Ewald speeded up evaluation of the gravitational potential by splitting it into two parts, \(Gm/r=Gm \, \mathrm{erf}(\alpha r)/r + Gm \, \mathrm{erfc}(\alpha r)/r\) (\(\alpha\) is an arbitrary constant) and then by applying Poisson summation formula on erfc terms, gravitational field at position \(\vec{r}\) can be written in the form
the first sum runs over whole computational domain, where at position \(\vec{r}_a\) is mass \(m_a\). Second sum runs over all neigbouring computational domains, which are at positions \(\vec{l}_{i_1,i_2,i_3}\) and \(A_S(\vec{r},\vec{r}_a,\vec{l}_{i_1,i_2,i_3})\) and \(A_L(\vec{r},\vec{r}_a,\vec{l}_{i_1,i_2,i_3})\) are short and long–range contributions, respectively. It is sufficient to take into account only few terms in eq. [ewald_sum]. The Ewald field for the acceleration, \(f_\mathrm{EF,a}\), is obtained using a similar decomposition. We modified Ewald method for problems with periodic boundary conditions in two directions and isolated boundary conditions in the third direction and for problems with periodic boundary conditions in one direction and isolated in two directions.
The gravity unit allows also to use a static external gravitational
field read from file. In this version, the field can be either
spherically symmetric or planar being only a function of the
zcoordinate. The external field file is a text file containing three
columns of numbers representing the coordinate, the potential and the
acceleration. The coordinate is the radial distance or zdistance from
the center of the external field given by runtime parameters. The
external field if mapped to a grid using a linear interpolation each
time the gravitational acceleration is calculated (in subroutine
physics/Gravity/Gravity_accelOneRow
).
Usage
To include the effects of gravity in your FlashX executable, include the option
withunit=physics/Gravity
on your command line when you configure the code with setup
. The
default implementation is Constant
, which can be overridden by
including the entire path to the specific implementation in the command
line or Config
file. The other available implementations are
Gravity/GravityMain/Planepar
, Gravity/GravityMain/Pointmass
and
Gravity/GravityMain/Poisson
. The Gravity
unit provides accessor
functions to get gravitational acceleration and potential. However, none
of the external field implementations of Section explicitly compute the
potential, hence they inherit the null implementation from the API for
accessing potential. The gravitation acceleration can be obtained either
on the whole domain, a single block or a single row at a time.
When building an application that solves the Possion equation for the
gravitational potential, additional storage is needed in unk
for
holding the last, as well as (usually) the previous, gravitational
potential field; and, depending on the Poisson solver used, additional
variables may be needed. The variables GPOT_VAR
and GPOT_VAR
,
and others as needed, will be automatically defined in Simulation.h
in those cases. See physics/Gravity/Gravity_potentialListOfBlocks
for more information.
Tree Gravity Unit Usage
Calculation of gravitational potential can be enabled by compiling in
this unit and setting the runtime parameter Gravity/useGravity
true.
The constant of gravity can be set independently by runtime parameter
Gravity/grv_bhNewton
; if it is not positive, the constant Newton
from the FlashX PhysicalConstants
database is used. If parameters
Grid/gr_bhPhysMACTW
or Grid/gr_bhPhysMACComm
are set, the
gravity unit MAC is used and it can be chosen by setting
Gravity/grv_bhMAC
to either ApproxPartialErr
or
MaxPartialErr
. If the first one is used, the order of the multipole
approximation is given by Gravity/grv_bhMPDegree
.
The maximum allowed partial error in gravitational acceleration is set
with the runtime parameter Gravity/grv_bhAccErr
. It has either the
meaning of an error in absolute acceleration or in relative acceleration
normalized by the acceleration from the previous timestep. The latter
is used if Gravity/grv_bhUseRelAccErr
is set to True, and in this
case the first call of the tree solver calculates the potential using
purely geometrical MAC (because the acceleration from the previous
timestep does not exist).
Boundary conditions are set by the runtime parameter
Gravity/grav_boundary_type
and they can be isolated
,
periodic
or mixed
. In the case of mixed boundary conditions,
runtime parameters Gravity/grav_boundary_type_x
,
Gravity/grav_boundary_type_y
and Gravity/grav_boundary_type_z
specify along which coordinate boundary conditions are periodic and
isolated (possible values are periodic
or isolated
). Arbitrary
combination of these values is permitted, thus suitable for problems
with planar resp. linear symmetry. It should work for computational
domain with arbitrary dimensions. The highest accuracy is reached with
blocks of cubic physical dimensions.
If runtime parameter Gravity/grav_boundary_type
is periodic
or
mixed
, then the Ewald field for appropriate symmetry is calculated
at the beginning of the simulation. Parameter
Gravity/grv_bhEwaldSeriesN
controls the range of indices
\(i_1,i_2,i_3\) in (eq. [ewald_sum]). There are two
implementations of the Ewald method: the new one (default) requires less
memory and it should be faster and of comparable accuracy as the old
one. The default implementation computes Ewald field minus the singular
\(1/r\) term and its partial derivatives on a single cubic grid, and
the Ewald field is then approximated by the first order Taylor formula.
Parameter Gravity/grv_bhEwaldNPer
controls number of grid points in
the \(x\) direction in the case of periodic
or in periodic
direction(s) in the case of mixed
boundary conditions. Since an
elongated computational domain is often desired when
Gravity/grav_boundary_type
is mixed
, the cubic grid would lead
to a huge field of data. In this case, the amount of necessary grid
points is reduced by using an analytical estimate to the Ewald field
sufficiently far away of the symmetry plane or axis.
The old implementation (from Flash4.2) is still present and is enabled
by adding bhtreeEwaldV42=1
on the setup command line. The Ewald
field is then stored in a nested set of grids, the first of them
corresponds in size to full computational domain, and each following
grid is half the size (in each direction) of the previous grid. Number
of nested grids is controlled by runtime parameter
Gravity/grv_bhEwaldNRefV42
. If Gravity/grv_bhEwaldNRefV42
is too
low to cover origin (where is the Ewald field discontinuous), then the
run is terminated. Each grid is composed of
Gravity/grv_bhEwaldFieldNxV42
\(\times\)
Gravity/grv_bhEwaldFieldNyV42
\(\times\)
Gravity/grv_bhEwaldFieldNzV42
points. When evaluation of the Ewald
Field at particular point is needed at any time during a run, the field
value is found by interpolation in a suitable level of the grid. Linear
or semiquadratic interpolation can be chosen by runtime parameter
Gravity/grv_bhLinearInterpolOnlyV42
(option true
corresponds to
linear interpolation). Semiquadratic interpolation is recommended only
in the case when there are periodic boundary conditions in two
directions.
The external gravitational field can be switched on by setting
Gravity/grv_useExternalPotential
true. The parameter
Gravity/grv_bhExtrnPotFile
gives the name of the file with the
external potential and Gravity/grv_bhExtrnPotType
specifies the
field symmetry: spherical
for the spherical symmetry and planez
for the planar symmetry with field being a function of the zcoordinate.
Parameters Gravity/grv_bhExtrnPotCenterY
,
Gravity/grv_bhExtrnPotCenterX
and Gravity/grv_bhExtrnPotCenterZ
specify the position (in the simulation coordinate system) of the
external field origin (the point where the radial or zcoordinate is
zero).
calculation.
Variable
Type
Default
Description
Gravit y/grv_bhNewton
real
1.0
constant of gravity; if \(<\) 0, it is obtained
from internal physical constants database
Gra vity/grv_bhMAC
string
“A pproxPartialErr”
MAC, other option: “MaxPartialErr”
Gravity/ grv_bhMPDegree
integer
2
degree of multipole in error estimate in APE MAC
Gravity/grv_ bhUseRelAccErr
logical
.false.
if .true., grv_bhAccErr has meaning of
relative error, otherwise absolute
Gravit y/grv_bhAccErr
real
0.1
maximum allowed error in gravitational
acceleration
conditions.
Variable
Type
Default
Description
Gravity/grav _boundary_type
string
“isolated”
or “periodic” or “mixed”
Gravity/grav_b oundary_type_x
string
“isolated”
or “periodic”
Gravity/grav_b oundary_type_y
string
“isolated”
or “periodic”
Gravity/grav_b oundary_type_z
string
“isolated”
or “periodic”
Gra vity/grv_bhEwald AlwaysGenerate
boolean
true
whether Ewald field should be regenerated
Gravity/grv_ bhEwaldSeriesN
integer
\(10\)
number of terms in the Ewald series
Gravity/g rv_bhEwaldNPer
integer
32
number of points+1 of the Taylor expansion
Gravity/gr v_bhEwaldFName
string
“ewald_coeffs”
file with coefficients of the Ewald field Taylor expansion
` Gravity/grv_bhE waldFieldNxV42`
integer
\(32\)
size of the Ewald field grid in xdirection
` Gravity/grv_bhE waldFieldNyV42`
integer
\(32\)
size of the Ewald field grid in ydirection
` Gravity/grv_bhE waldFieldNzV42`
integer
\(32\)
size of the Ewald field grid in zdirection
Gravity/grv_ bhEwaldNRefV42
integer
1
number of refinement levels (nested grids) for the Ewald
field; if \(<\) 0, determined automatically
Gravi ty/grv_bhLinearI nterpolOnlyV42
logical
.true.
if .false., semiquadratic interpolation is used for
interpolation in the Ewald field
`` Gravity/grv_bhEw aldFNameAccV42``
string
” ewald_field_acc”
file with the Ewald field for acceleration
`` Gravity/grv_bhEw aldFNamePotV42``
string
” ewald_field_pot”
file with coefficients of the Ewald field for potential
Tree gravity unit parameters controlling the external gravitational field.
Variable 
Type 
Default 
Description 


logical 
.false. 
whether to use external field 

logical 
.true. 
whether to use gravitational field calculated by the 
tree solver 


string 
“externa l_potential.dat” 
file containing the external gravitational field 

string 
“planez” 
type of the external field: planar or spherical symmetry 
` Gravity/grv_bhE xtrnPotCenterX` 
real 
0.0 
xcoordinate of the center of the external field 
` Gravity/grv_bhE xtrnPotCenterY` 
real 
0.0 
ycoordinate of the center of the external field 
` Gravity/grv_bhE xtrnPotCenterZ` 
real 
0.0 
zcoordinate of the center of the external field 
Unit Tests
There are two unit tests for the gravity unit. Poisson3
is
essentially the Maclaurin spheroid problem described in . Because an
analytical solution exists, the accuracy of the gravitational solver can
be quantified. The second test, Poisson3_active
is a modification of
Poisson3
to test the mapping of particles in
Grid/Grid_mapParticlesToMesh
. Some of the mesh density is
redistributed onto particles, and the particles are then mapped back to
the mesh, using the analytical solution to verify completeness. This
test is similar to the simulation PoisParticles
discussed in .
PoisParticles
is based on the HuangGreengard Poisson gravity test
described in .