External Potentials
SWIFT can be run with an external potential on this page we will summarize the current potentials which can be run with SWIFT and how to implement your own potential in SWIFT.
Implemented External Potentials
Currently there are several potentials implemented in SWIFT. On this page we
give a short overview of the potentials that are implemented in the code. They
are all switched on at configuration time via the argument
--with-ext-potential=XYZ and get their own independant section in the
parameter file. The name of the potential to pass to the configuration script is
given in the parenthesis of the titles below; for instance
--with-ext-potential=herquist.
1. No potential (none)
This is the default setup. No external potential is used. There are no parameters associated with this model.
2. Constant acceleration (constant)
This “potential” just applies a constant acceleration vector at every position \(\vec{x}\) in the simulation volume. This is very common in idealised test cases like the Rayleigh-Taylor instability or engineering applications relying on Earth’s constant gravity field.
\(\phi(\vec{x}) = -\vec{g} \cdot \vec{x}\)
\(\vec{a}(\vec{x}) = \vec{g}\)
The only parameter of this model is the vector \(\vec{g}\) given in cgs units:
ConstantPotential:
g_cgs: [0., 0., -9,81] # Earth acceleration along z-axis (cgs units)
3. Point mass potential (point-mass)
A simple single point mass that can be placed at any position \(\vec{p}\). The position can be specified either absolutely or with respect to the centre of the simulation volume.
\(\phi(\vec{x}) = -\displaystyle\frac{G_{\rm N}m}{|r|}\)
\(\vec{a}(\vec{x}) = -\displaystyle\frac{G_{\rm N}m}{|r|^3}\vec{r}\),
where \(\vec{r} = \vec{x} - \vec{p}\).
The code also imposes an extra limit of the size of the particles’ time-step. The time-step has to be shorter than \(\Delta t_{pot} = C |\vec{a}(\vec{x})| / |\dot{\vec{a}}(\vec{x})|\). This criterion is designed to make sure the changes in accelerations remain small. The jerk (\(\dot{\vec{a}}\equiv\frac{d}{dt}\vec{a}\)) is calculated from the positions of the particles and only includes the contribution of the external potential itself. The other criteria (CFL, self-gravity, …) are applied on top of this criterion. The dimensionless constant \(C\) defaults to FLT_MAX if left unspecified, effectively making SWIFT run without this time-step criterion.
The parameters of the model are:
PointMassPotential:
position: [3., 4., 5.] # Location of the point mass (internal units)
useabspos: 1 # Use absolute positions (0 for relative to centre)
mass: 5.972e24 # Mass of the point (internal units)
timestep_mult: 0.1 # (Optional) The dimensionless constant C in the time-step condition
4. Plummer potential (point-mass-softened)
A single point mass with a fixed softening length using a Plummer shape that can be placed at any position \(\vec{p}\). The position can be specified either absolutely or with respect to the centre of the simulation volume.
\(\phi(\vec{x}) = -\displaystyle\frac{G_{\rm N}m}{\sqrt{|r|^2 + \epsilon^2}}\)
\(\vec{a}(\vec{x}) = -\displaystyle\frac{G_{\rm N}m}{(|r|^2 + \epsilon^2)^{3/2}}\vec{r}\),
where \(\vec{r} = \vec{x} - \vec{p}\).
The code also imposes an extra limit of the size of the particles’ time-step. The time-step has to be shorter than \(\Delta t_{pot} = C |\vec{a}(\vec{x})| / |\dot{\vec{a}}(\vec{x})|\). This criterion is designed to make sure the changes in accelerations remain small. The jerk (\(\dot{\vec{a}}\equiv\frac{d}{dt}\vec{a}\)) is calculated from the positions of the particles and only includes the contribution of the external potential itself. The other criteria (CFL, self-gravity, …) are applied on top of this criterion. The dimensionless constant \(C\) defaults to FLT_MAX if left unspecified, effectively making SWIFT run without this time-step criterion.
The parameters of the model are:
PointMassPotential:
position: [3., 4., 5.] # Location of the point mass (internal units)
useabspos: 1 # Use absolute positions (0 for relative to centre)
mass: 5.972e24 # Mass of the point (internal units)
softening: 0.01 # Softening length (internal units)
timestep_mult: 0.1 # (Optional) The dimensionless constant C in the time-step condition
5. Isothermal potential (isothermal)
An isothermal potential which corresponds to a density profile which is \(\propto r^{-2}\) and a potential which is logarithmic. This potential is entirely set by its centre \(\vec{p}\) and the (radially-constant) rotation velocity of the potential \(v_{\rm rot}\). The centre of the potential is softened.
\(\phi(\vec{x}) = -\displaystyle\frac{1}{4\pi G_{\rm N}}\log(\sqrt{|r|^2 + \epsilon^2})\)
\(\vec{a}(\vec{x}) = -\displaystyle\frac{v_{\rm rot}^2} {|r|^2 + \epsilon^2}\),
where \(\vec{r} = \vec{x} - \vec{p}\).
The code also imposes an extra limit of the size of the particles’ time-step. The time-step has to be shorter than \(\Delta t_{pot} = C |\vec{a}(\vec{x})| / |\dot{\vec{a}}(\vec{x})|\). This criterion is designed to make sure the changes in accelerations remain small. The jerk (\(\dot{\vec{a}}\equiv\frac{d}{dt}\vec{a}\)) is calculated from the positions of the particles and only includes the contribution of the external potential itself. The other criteria (CFL, self-gravity, …) are applied on top of this criterion. The dimensionless constant \(C\) defaults to FLT_MAX if left unspecified, effectively making SWIFT run without this time-step criterion.
The parameters of the model are:
IsothermalPotential:
position: [3., 4., 5.] # Location of the centre of the profile (internal units)
useabspos: 1 # Use absolute positions (0 for relative to centre)
vrot: 200 # Rotation velocity of the profile (internal units)
softening: 0.01 # Softening length (internal units)
timestep_mult: 0.1 # (Optional) The dimensionless constant C in the time-step condition
6. Hernquist potential (hernquist)
We can set up a potential as given by Hernquist (1990):
\(\Phi(r) = - \frac{G_{\rm N}M}{r+a},\)
where \(M\) is the total Hernquist mass and :math: a is the Hernquist- equivalent scale radius of the potential. The potential can be set at any position in the box. It adds an additional time step constraint that limits the time step to a maximum of a specified fraction of the circular orbital time at the current radius of the particle. The other criteria (CFL, self-gravity, …) are applied on top of this criterion. For example, a fraction of 0.01 means that an orbital period would be resolved by 100 time steps.
In the most basic version, the Hernquist potential can be run by providing only the Hernquist mass, scale length, softening length and fraction of the orbital time for the time stepping. In this case the model parameters may look something like:
HernquistPotential:
useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
position: [0.,0.,0.] # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
mass: 1e12 # Mass of the Hernquist potential
scalelength: 2.0 # scale length a (internal units)
timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, determines the fraction of the orbital time we use to do the time integration; fraction of 0.01 means we resolve an orbit with 100 timesteps
epsilon: 0.2 # Softening size (internal units)
Besides the basic version, it is also possible to run the Hernquist potential for idealised disk galaxies that follow the approach of Springel, Di Matteo & Hernquist (2005). This potential, however, uses a corrected value of the formulation that improves the match with the NFW profile (below) with the same M200 (Nobels+ in prep). Contrary to the above (idealizeddisk: 0 setup), the idealised disk setup runs by specifying one out of \(M_{200}\), \(V_{200}\), or \(R_{200}\), plus concentration and reduced Hubble constant.
In this case, we don’t provide the ‘mass’ and ‘scalelength’ parameters, but ‘M200’ (or ‘V200’, or ‘R200’) and ‘concentration’ \(c\), as well as reduced Hubble constant \(h\) to define the potential. The parameters of the model may look something like:
HernquistPotential:
useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
position: [0.,0.,0.] # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
idealizeddisk: 1 # Run with an idealized galaxy disk
M200: 137.0 # M200 of the galaxy disk
h: 0.704 # reduced Hubble constant (value does not specify the used units!)
concentration: 9.0 # concentration of the Halo
diskfraction: 0.040 # Disk mass fraction
bulgefraction: 0.0 # Bulge mass fraction
timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, determines the fraction of the orbital time we use to do the time integration; fraction of 0.01 means we resolve an orbit with 100 timesteps
epsilon: 0.2 # Softening size (internal units)
The user should specify one out of ‘M200’, ‘V200’, or ‘R200’ to define the potential. The reduced Hubble constant is then used to determine the other two. Then, the scale radius is calculated as \(R_s = R_{200}/c\), where \(c\) is the concentration, and the Hernquist-equivalent scale-length is calculated as:
\(a = \frac{b+\sqrt{b}}{1-b} R_{200}\),
where \(b = \frac{2}{c^2}(\ln(1+c) - \frac{c}{1+c})\).
Two examples using the Hernquist potential can be found in swiftsim/examples/GravityTests/.
The Hernquist_radialinfall example puts 5 particles with zero velocity in a Hernquist
potential, resulting in radial orbits. The Hernquist_circularorbit``example puts three
particles on a circular orbit in a Hernquist potential, one at the inner region, one at the
maximal circular velocity, and one in the outer region. To run these examples, SWIFT must
be configured with the flag ``--with-ext-potential=hernquist, or hernquist-sdmh05
(see below).
7. Hernquist SDMH05 potential (hernquist-sdmh05)
This is the same potential as Hernquist with the difference being the way that the idealised disk part is calculated. This potential uses exactly the same approach as Springel, Di Matteo & Hernquist (2005), which means that ICs generated with the original MakeNewDisk code can be used with this potential. Contrary to the updated Hernquist potential (above), it is not possible to have an identically matched NFW potential in this case.
The parameters needed for this potential are the same set of variables as above, i.e. ‘mass’ and ‘scalelength’ when we don’t use the idealised disk, and ‘concentration’ plus one out of ‘M200’, ‘V200’, or ‘R200’ if we do. As one of the three is provided, the reduced Hubble constant is used to calculate the other two. Then, the scale radius is calculated using the NFW definition, \(R_s = R_{200}/c\), and the Hernquist- equivalent scale length is given by
\(a = R_s \sqrt{2(\ln(1+c) - \frac{c}{1+c})}.\)
Runs that provide e.g. M200 and c (using idealised disk) are thus equivalent to providing mass and scale length if calculated by the above equation (without idealized disk).
9. NFW poential + Miyamoto-Nagai potential (nfw-mn)
This includes and NFW potential (identical to nfw) plus an axisymmetric Miyamoto-Nagai potential. The Miyamoto-Nagai potential is given by:
\(\Phi(R,z) = - \frac{G_{\rm N} M_{d}}{\sqrt{R^2 + \left ( R_d + \sqrt{z^2 + Z_d^2} \right )^2}}\),
where \(R^2 = x^2 + y^2\) is the projected radius and \(M_d\), \(R_d\), \(Z_d\) are the mass, scalelength and scaleheight of the disk (in internal units), respectively.
The parameters of the model are:
NFW_MNPotential:
useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
position: [0.,0.,0.] # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
M_200: 137.0 # M200 of the halo in internal units
critical_density: 123.4 # Critical density of the universe in internal units
concentration: 9.0 # concentration of the NFW halo
Mdisk: 3.0 # Mass of the disk in internal units
Rdisk: 3.0 # Disk scale-length in internal units
Zdisk: 3.0 # Disk scale-height in internal units
timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
10. Sine wave (sine-wave)
This “potential” is designed for specific idealised tests. It is a basic (not really physical!) sinusoid wave along the x-axis with a unit wavelength and a tunable amplitude and growth time.
\(\phi(\vec{x}) = A \cos\left(2 \pi x_{0}\right) / 2\pi\)
\(a_0(\vec{x}) = A \sin\left(2 \pi x_{0}\right)\),
where the 0 subscript indicates the x-component. The y- and z-components are zero.
The amplitude \(A\) can be growing linearly to its maximum over a fixed time-scale.
Optionally, a constant maximail time-step size can be used with this potential.
The parameters of the model are:
SineWavePotential:
amplitude: 2.5 # The amplitude of the wave (internal units)
growth_time: 1.2 # Time for the amplitude to grow to its final value (internal units)
timestep_limit: 5.0 # (Optional) The max time-step of the particles (internal units)
11. Disc Patch (disc-patch)
A potential corresponding to a vertical slice through a galactic disk. This follows the definitions of Creasey, Theuns & Bower (2013) equations (16) and (17). The potential is implemented along the x-axis.
12. MWPotential2014 (MWPotential2014)
This potential is based on galpy’s MWPotential2014 from Jo Bovy (2015) and consists in a NFW potential for the halo, an axisymmetric Miyamoto-Nagai potential for the disk and a bulge modeled by a power spherical law with exponential cut-off. The bulge is given by the density:
\(\rho(r) = A \left( \frac{r_1}{r} \right)^\alpha \exp \left( - \frac{r^2}{r_c^2} \right)\),
where \(A\) is an amplitude, \(r_1\) is a reference radius for amplitude, \(\alpha\) is the inner power and \(r_c\) is the cut-off radius.
The resulting potential is:
\(\Phi_{\mathrm{MW}}(R, z) = f_1 \Phi_{\mathrm{NFW}} + f_2 \Phi_{\mathrm{MN}} + f_3 \Phi_{\text{bulge}}\),
where \(R^2 = x^2 + y^2\) is the projected radius and \(f_1\), \(f_2\) and \(f_3\) are three coefficients that adjust the strength of each individual component.
The parameters of the model are:
MWPotential2014Potential:
useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
position: [0.,0.,0.] # Location of centre of potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
timestep_mult: 0.005 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
epsilon: 0.001 # Softening size (internal units)
concentration: 9.823403437774843 # concentration of the Halo
M_200_Msun: 147.41031542774076e10 # M200 of the galaxy disk (in M_sun)
H: 1.2778254614201471 # Hubble constant in units of km/s/Mpc
Mdisk_Msun: 6.8e10 # Mass of the disk (in M_sun)
Rdisk_kpc: 3.0 # Effective radius of the disk (in kpc)
Zdisk_kpc: 0.280 # Scale-height of the disk (in kpc)
amplitude: 1.0 # Amplitude of the bulge
r_1_kpc: 1.0 # Reference radius for amplitude of the bulge (in kpc)
alpha: 1.8 # Exponent of the power law of the bulge
r_c_kpc: 1.9 # Cut-off radius of the bulge (in kpc)
potential_factors: [0.4367419745056084, 1.002641971008805, 0.022264787598364262] #Coefficients that adjust the strength of the halo (1st component), the disk (2nd component) and the bulge (3rd component)
Note that the default value of the “Hubble constant” here seems odd. As it enters multiplicatively with the \(f_1\) term, the absolute normalisation is actually not important.
How to implement your own potential
The first step in implementing your own potential is making a directory of your
potential in the src/potential folder and creating a file in the folder
called potential.h.
Configuring the potential
To get started you can copy a potential.h file from an already implemented
potential. In this potential the header guards (e.g. #IFDEF <>) need to be
changed to the specific potential and the struct and
potential_init_backend need to be changed such that it uses your potential
and reads the correct potential from the parameter file during running the
program.
Add the potential to the potential.h file in the src directory such that
the program knows that it is possible to run with this potential.
Furthermore during the configuration of the code it also needs to be clear for
the program that the code can be configured to run with the different
potentials. This means that the configure.ac file needs to be changed.
This can be done to add an other case in the potential:
case "$with_potential" in
none)
AC_DEFINE([EXTERNAL_POTENTIAL_NONE], [1], [No external potential])
;;
newpotential)
AC_DEFINE([EXTERNAL_POTENTIAL_NEWPOTENTIAL], [1], [New external potential])
;;
After this change it is possible to configure the code to use your new potential.