Commit c808c1bf authored by Folkert Nobels's avatar Folkert Nobels Committed by Matthieu Schaller

Add Hernquist and NFW potential to SWIFT

parent 56a38861
......@@ -1518,7 +1518,7 @@ esac
# External potential
AC_ARG_WITH([ext-potential],
[AS_HELP_STRING([--with-ext-potential=<pot>],
[external potential @<:@none, point-mass, point-mass-ring, point-mass-softened, isothermal, softened-isothermal, disc-patch, sine-wave, default: none@:>@]
[external potential @<:@none, point-mass, point-mass-ring, point-mass-softened, isothermal, softened-isothermal, nfw, hernquist, disc-patch, sine-wave, default: none@:>@]
)],
[with_potential="$withval"],
[with_potential="none"]
......@@ -1533,6 +1533,12 @@ case "$with_potential" in
isothermal)
AC_DEFINE([EXTERNAL_POTENTIAL_ISOTHERMAL], [1], [Isothermal external potential])
;;
hernquist)
AC_DEFINE([EXTERNAL_POTENTIAL_HERNQUIST], [1], [Hernquist external potential])
;;
nfw)
AC_DEFINE([EXTERNAL_POTENTIAL_NFW], [1], [Navarro-Frenk-White external potential])
;;
disc-patch)
AC_DEFINE([EXTERNAL_POTENTIAL_DISC_PATCH], [1], [Disc-patch external potential])
;;
......
.. External potentials in SWIFT
Folkert Nobels, 25th October 2018
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:
1. No potential (none)
2. Point mass potential (point-mass): classical point mass, can be placed at
a position with a mass.
3. Plummer potential (point-mass-softened): in the code a softended point mass
corresponds to a Plummer potential, can be placed at a position with a mass.
4. Isothermal potential (isothermal): An isothermal potential which corresponds
to a density profile which is :math:`\propto r^{-2}` and a potential which is
logarithmic. This potential has as free parameters the rotation velocity
and the position.
5. Hernquist potential (hernquist): A potential that is given by the Hernquist
potential:
:math:`\Phi(r) = - \frac{GM}{r+a}.`
The free paramters of Hernquist potential are mass, scale length,
and softening. The potential can be set at any position in the box.
6. NFW potential (nfw): The most used potential to describe dark matter halos, the
potential is given by:
:math:`\Phi(r) = - \frac{4\pi G \rho_0 R_s^3}{r} \ln \left( 1+
\frac{r}{R_s} \right).`
This potential has as free paramters the concentration of the DM halo, the
virial mass (:math:`M_{200}`) and the critical density.
7. Sine wave (sine-wave)
8. Point mass ring (point-mass-ring)
9. Disc Patch (disc-patch)
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.
......@@ -21,6 +21,7 @@ difference is the parameter file that will need to be adapted for SWIFT.
HydroSchemes/index
Cooling/index
EquationOfState/index
ExternalPotentials/index
NewOption/index
Task/index
VELOCIraptorInterface/index
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.988e+33 # Grams
UnitLength_in_cgs: 3.086e+21 # Centimeters
UnitVelocity_in_cgs: 1e5 # Centimeters per second
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters governing the time integration (Set dt_min and dt_max to the same value for a fixed time-step run.)
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 2.0 # The end time of the simulation (in internal units).
dt_min: 1e-10 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e0 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the snapshots
Snapshots:
basename: output # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 1e-3 # Time difference between consecutive outputs (in internal units)
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e0 # Time between statistics output
# Parameters related to the initial conditions
InitialConditions:
file_name: circularorbitshernquist.hdf5 # The file to read
periodic: 0
# Hernquist potential parameters
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: 2e12 # Mass of the Hernquist potential
scalelength: 10.0 # Scale length of the potential
timestep_mult: 0.005 # Dimensionless pre-factor for the time-step condition
epsilon: 0.1 # Softening size (internal units)
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2018 Folkert Nobels (nobels@strw.leidenuniv.nl)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
from galpy.potential import NFWPotential
from galpy.orbit import Orbit
from galpy.util import bovy_conversion
import numpy as np
import matplotlib.pyplot as plt
from astropy import units
import h5py as h5
C = 8.0
M_200 = 2.0
N_PARTICLES = 3
print("Initial conditions written to 'test_nfw.hdf5'")
pos = np.zeros((3, 3))
pos[0, 2] = 50.0
pos[1, 2] = 10.0
pos[2, 2] = 2.0
pos = pos + 500.0
vel = np.zeros((3, 3))
vel[0, 1] = 348.0
vel[1, 1] = 466.9
vel[2, 1] = 348.0
ids = np.array([1.0, 2.0, 3.0])
mass = np.array([1.0, 1.0, 1.0])
# File
file = h5.File("circularorbitshernquist.hdf5", "w")
# Units
grp = file.create_group("/Units")
grp.attrs["Unit length in cgs (U_L)"] = 3.086e21
grp.attrs["Unit mass in cgs (U_M)"] = 1.988e33
grp.attrs["Unit time in cgs (U_t)"] = 3.086e16
grp.attrs["Unit current in cgs (U_I)"] = 1.0
grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
# Header
grp = file.create_group("/Header")
grp.attrs["BoxSize"] = 1000.0
grp.attrs["NumPart_Total"] = [0, N_PARTICLES, 0, 0, 0, 0]
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
grp.attrs["NumPart_ThisFile"] = [0, N_PARTICLES, 0, 0, 0, 0]
grp.attrs["Time"] = 0.0
grp.attrs["NumFilesPerSnapshot"] = 1
grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
grp.attrs["Dimension"] = 3
# Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = 1
# Particle group
grp1 = file.create_group("/PartType1")
ds = grp1.create_dataset("Velocities", (N_PARTICLES, 3), "f", data=vel)
ds = grp1.create_dataset("Masses", (N_PARTICLES,), "f", data=mass)
ds = grp1.create_dataset("ParticleIDs", (N_PARTICLES,), "L", data=ids)
ds = grp1.create_dataset("Coordinates", (N_PARTICLES, 3), "d", data=pos)
file.close()
#!/usr/bin/env python
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2018 Folkert Nobels (nobels@strw.leidenuniv.nl)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
import numpy as np
import h5py
import matplotlib.pyplot as plt
from scipy.integrate import odeint
t = np.linspace(0, 40, 1e5)
y0 = [0, 10]
a = 30.0
G = 4.300927e-06
M = 1e15
GM = G * M
lengthrun = 2001
numbpar = 3
radius = np.zeros((numbpar, lengthrun))
xx = np.zeros((numbpar, lengthrun))
yy = np.zeros((numbpar, lengthrun))
zz = np.zeros((numbpar, lengthrun))
time = np.zeros(lengthrun)
for i in range(0, lengthrun):
Data = h5py.File("output_%04d.hdf5" % i, "r")
header = Data["Header"]
time[i] = header.attrs["Time"]
particles = Data["PartType1"]
positions = particles["Coordinates"]
xx[:, i] = positions[:, 0] - 500.0
yy[:, i] = positions[:, 1] - 500.0
zz[:, i] = positions[:, 2] - 500.0
col = ["b", "r", "c", "y", "k"]
print(np.shape(xx), np.shape(yy), np.shape(zz))
for i in range(0, numbpar):
plt.plot(xx[i, :], yy[i, :], col[i])
plt.ylabel("y (kpc)")
plt.xlabel("x (kpc)")
plt.savefig("xyplot.png")
plt.close()
for i in range(0, numbpar):
plt.plot(xx[i, :], zz[i, :], col[i])
plt.ylabel("z (kpc)")
plt.xlabel("x (kpc)")
plt.savefig("xzplot.png")
plt.close()
for i in range(0, numbpar):
plt.plot(yy[i, :], zz[i, :], col[i])
plt.ylabel("z (kpc)")
plt.xlabel("y (kpc)")
plt.savefig("yzplot.png")
plt.close()
#!/bin/bash
if [ ! -e circularorbitshernquist.hdf5 ]
then
echo "Generate initial conditions for circular orbits"
if command -v python3 &>/dev/null; then
python3 makeIC.py
else
python makeIC.py
fi
fi
# self gravity G, external potential g, hydro s, threads t and high verbosity v
../swift -g -t 6 hernquistcirc.yml 2>&1 | tee output.log
echo "Save plots of the circular orbits"
if command -v python3 &>/dev/null; then
python3 plotprog.py
else
python plotprog.py
fi
This example generates 5 particles at radii of 10, 20, 30, 40 and 50 kpc
without velocitiy and follows the evolution of these particles in an Hernquist
potential as they are free falling.
# Define the system of units to use internally.
InternalUnitSystem:
UnitMass_in_cgs: 1.98848e33 # M_sun
UnitLength_in_cgs: 3.08567758e21 # kpc
UnitVelocity_in_cgs: 1e5 # km/s
UnitCurrent_in_cgs: 1 # Amperes
UnitTemp_in_cgs: 1 # Kelvin
# Parameters governing the time integration
TimeIntegration:
time_begin: 0. # The starting time of the simulation (in internal units).
time_end: 40. # The end time of the simulation (in internal units).
dt_min: 9e-10 # The minimal time-step size of the simulation (in internal units).
dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-3 # Time between statistics output
# Parameters governing the snapshots
Snapshots:
basename: hernquist # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 0.02 # Time difference between consecutive outputs (in internal units)
# Parameters related to the initial conditions
InitialConditions:
file_name: Hernquist.hdf5 # The file to read
periodic: 1
shift: [200.,200.,200.] # Shift all particles to be in the potential
# External potential parameters
HernquistPotential:
useabspos: 0 # Whether to use absolute position (1) or relative potential to centre of box (0)
position: [0.,0.,0.]
mass: 1e9
scalelength: 1.0
timestep_mult: 0.01 # controls time step
epsilon: 2.0 # No softening at the centre of the halo
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2018 Folkert Nobels (nobels@strw.leidenuniv.nl)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
##############################################################################
import h5py
import sys
import numpy
import math
import random
import numpy as np
# Generates N particles in a spherical distribution centred on [0,0,0], to be moved in an isothermal potential
# usage: python makeIC.py 1000 0 : generate 1000 particles on circular orbits
# python makeIC.py 1000 1 : generate 1000 particles with Lz/L uniform in [0,1]
# all particles move in the xy plane, and start at y=0
# physical constants in cgs
NEWTON_GRAVITY_CGS = 6.67408e-8
SOLAR_MASS_IN_CGS = 1.98848e33
PARSEC_IN_CGS = 3.08567758e18
YEAR_IN_CGS = 3.15569252e7
# choice of units
const_unit_length_in_cgs = 1000 * PARSEC_IN_CGS
const_unit_mass_in_cgs = SOLAR_MASS_IN_CGS
const_unit_velocity_in_cgs = 1e5
# Properties of the Hernquist potential
Mass = 1e15
scaleLength = 30.0 # kpc
# derived units
const_unit_time_in_cgs = const_unit_length_in_cgs / const_unit_velocity_in_cgs
const_G = (
NEWTON_GRAVITY_CGS
* const_unit_mass_in_cgs
* const_unit_time_in_cgs
* const_unit_time_in_cgs
/ (const_unit_length_in_cgs * const_unit_length_in_cgs * const_unit_length_in_cgs)
)
print("G=", const_G)
def hernquistcircvel(r, M, a):
""" Function that calculates the circular velocity in a
Hernquist potential.
@param r: radius from centre of potential
@param M: mass of the Hernquist potential
@param a: Scale length of the potential
@return: circular velocity
"""
return (const_G * M * r) ** 0.5 / (r + a)
# Parameters
periodic = 1 # 1 For periodic box
boxSize = 400.0 # [kpc]
Radius = 100.0 # maximum radius of particles [kpc]
G = const_G
N = 5
L = N ** (1.0 / 3.0)
fileName = "Hernquist.hdf5"
# ---------------------------------------------------
numPart = N
mass = 1
# --------------------------------------------------
# File
file = h5py.File(fileName, "w")
# Units
grp = file.create_group("/Units")
grp.attrs["Unit length in cgs (U_L)"] = const_unit_length_in_cgs
grp.attrs["Unit mass in cgs (U_M)"] = const_unit_mass_in_cgs
grp.attrs["Unit time in cgs (U_t)"] = (
const_unit_length_in_cgs / const_unit_velocity_in_cgs
)
grp.attrs["Unit current in cgs (U_I)"] = 1.0
grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
# Header
grp = file.create_group("/Header")
grp.attrs["BoxSize"] = boxSize
grp.attrs["NumPart_Total"] = [0, numPart, 0, 0, 0, 0]
grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
grp.attrs["NumPart_ThisFile"] = [0, numPart, 0, 0, 0, 0]
grp.attrs["Time"] = 0.0
grp.attrs["NumFilesPerSnapshot"] = 1
grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
grp.attrs["Dimension"] = 3
# set seed for random number
numpy.random.seed(1234)
# Particle group
grp1 = file.create_group("/PartType1")
# generate particle positions
# radius = Radius * (numpy.random.rand(N))**(1./3.) + 10.
radius = np.zeros(N)
radius[0] = 10
radius[1] = 20
radius[2] = 30
radius[3] = 40
radius[4] = 50
# this part is not even used:
# ctheta = -1. + 2 * numpy.random.rand(N)
# stheta = numpy.sqrt(1.-ctheta**2)
# phi = 2 * math.pi * numpy.random.rand(N)
# end
r = numpy.zeros((numPart, 3))
r[:, 0] = radius
# import matplotlib.pyplot as plt
# plt.plot(r[:,0],'.')
# plt.show()
# print('Mass = ', Mass)
# print('radius = ', radius)
# print('scaleLength = ',scaleLength)
#
v = numpy.zeros((numPart, 3))
# v[:,0] = hernquistcircvel(radius,Mass,scaleLength)
omega = v[:, 0] / radius
period = 2.0 * math.pi / omega
print("period = minimum = ", min(period), " maximum = ", max(period))
print("Circular velocity = minimum =", min(v[:, 0]), " maximum = ", max(v[:, 0]))
omegav = omega
v[:, 0] = -omegav * r[:, 1]
v[:, 1] = omegav * r[:, 0]
ds = grp1.create_dataset("Velocities", (numPart, 3), "f", data=v)
m = numpy.full((numPart,), mass, dtype="f")
ds = grp1.create_dataset("Masses", (numPart,), "f", data=m)
ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False)
ds = grp1.create_dataset("ParticleIDs", (numPart,), "L", data=ids)
ds = grp1.create_dataset("Coordinates", (numPart, 3), "d", data=r)
file.close()
#!/usr/bin/env python
###############################################################################
# This file is part of SWIFT.
# Copyright (c) 2018 Folkert Nobels (nobels@strw.leidenuniv.nl)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
import numpy as np
import h5py
import matplotlib.pyplot as plt
from scipy.integrate import odeint
lengthrun = 2001
numbpar = 5
radius = np.zeros((numbpar, lengthrun))
time = np.zeros(lengthrun)
for i in range(0, lengthrun):
Data = h5py.File("hernquist_%04d.hdf5" % i, "r")
header = Data["Header"]
time[i] = header.attrs["Time"]
particles = Data["PartType1"]
positions = particles["Coordinates"]
radius[:, i] = positions[:, 0] - 200.0
col = ["b", "r", "c", "y", "k"]
for i in range(0, numbpar):
plt.plot(time, radius[i, :], col[i])
plt.axhline(np.max(radius[i, :]), color=col[i], linestyle="--")
plt.axhline(-np.max(radius[i, :]), color=col[i], linestyle="--")
plt.ylabel("Radial distance (kpc)")
plt.xlabel("Simulation time (internal units)")
plt.savefig("radial_infall.png")
plt.close()
#!/bin/bash
# Generate the initial conditions if they are not present.
if [ ! -e Hernquist.hdf5 ]
then
echo "Generate initial conditions for radial orbits"
if command -v python3 &>/dev/null; then
python3 makeIC.py
else
python makeIC.py
fi
fi
rm -rf hernquist_*.hdf5
../swift -g -t 1 hernquist.yml 2>&1 | tee output.log
echo "Make plots of the radially free falling particles"
if command -v python3 &>/dev/null; then
python3 plotprog.py
else
python plotprog.py
fi
......@@ -86,7 +86,7 @@ for i in range(402):
time_snap[i] = f["Header"].attrs["Time"]
E_kin_snap[i] = np.sum(0.5 * mass * (vel_x[:]**2 + vel_y[:]**2 + vel_z[:]**2))
E_pot_snap[i] = np.sum(-mass * Vrot**2 * log(r))
E_pot_snap[i] = np.sum(mass * Vrot**2 * log(r))
E_tot_snap[i] = E_kin_snap[i] + E_pot_snap[i]
Lz_snap[i] = np.sum(Lz)
......
......@@ -31,6 +31,8 @@ InitialConditions:
# External potential parameters
IsothermalPotential:
useabspos: 0 # Whether to use absolute position (1) or relative potential to centre of box (0)
position: [0.,0.,0.]
vrot: 200. # rotation speed of isothermal potential in internal units
timestep_mult: 0.01 # controls time step
epsilon: 0. # No softening at the centre of the halo
This just provides a test that the NFW potential is giving the correct orbit
for an elliptical orbit as calculated by Jo Bovy's galpy package. If
galpy is not installed on your system you can install it by using:
pp install galpy --user
################################################################################
# This file is part of SWIFT.
# Copyright (c) 2018 Ashley Kelly ()
# Folkert Nobels (nobels@strw.leidenuniv.nl)
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as published
# by the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
################################################################################
import numpy as np
import matplotlib.pyplot as plt
from astropy import units
import h5py as h5
C = 8.0
M_200 = 2.0
N_PARTICLES = 1
print("\nInitial conditions written to 'test_nfw.hdf5'")
pos = np.array([8.0, 0.0, 0.0]) + 500.0
vel = np.array([0.0, 240.0, 5.0])
ids = np.array([1.0])
mass = np.array([1.0])
# File
file = h5.File("test_nfw.hdf5", "w")
# Units
grp = file.create_group("/Units")
grp.attrs["Unit length in cgs (U_L)"] = 3.086e21
grp.attrs["Unit mass in cgs (U_M)"] = 1.988e33
grp.attrs["Unit time in cgs (U_t)"] = 3.086e16
grp.attrs["Unit current in cgs (U_I)"] = 1.0
grp.attrs["Unit temperature in cgs (U_T)"] = 1.0
# Header
grp = file.create_group("/Header")