From c808c1bf13b986be2d5a11a866dafe580b84d98a Mon Sep 17 00:00:00 2001
From: Folkert Nobels <nobels@strw.leidenuniv.nl>
Date: Tue, 6 Nov 2018 09:08:31 +0000
Subject: [PATCH] Add Hernquist and NFW potential to SWIFT
---
configure.ac | 8 +-
doc/RTD/source/ExternalPotentials/index.rst | 80 ++++++
doc/RTD/source/index.rst | 1 +
.../Hernquist_circularorbit/hernquistcirc.yml | 38 +++
examples/Hernquist_circularorbit/makeIC.py | 81 ++++++
examples/Hernquist_circularorbit/plotprog.py | 77 ++++++
examples/Hernquist_circularorbit/run.sh | 23 ++
examples/Hernquist_radialinfall/README | 3 +
examples/Hernquist_radialinfall/hernquist.yml | 39 +++
examples/Hernquist_radialinfall/makeIC.py | 167 +++++++++++
examples/Hernquist_radialinfall/plotprog.py | 50 ++++
examples/Hernquist_radialinfall/run.sh | 24 ++
examples/IsothermalPotential/energy_plot.py | 2 +-
examples/IsothermalPotential/isothermal.yml | 2 +
examples/NFW_Halo/README | 5 +
examples/NFW_Halo/makeIC.py | 75 +++++
examples/NFW_Halo/makePlots.py | 73 +++++
examples/NFW_Halo/run.sh | 20 ++
examples/NFW_Halo/test.yml | 41 +++
examples/parameter_example.yml | 20 ++
src/potential.h | 4 +
src/potential/hernquist/potential.h | 223 +++++++++++++++
src/potential/isothermal/potential.h | 16 +-
src/potential/nfw/potential.h | 260 ++++++++++++++++++
src/potential/point_mass/potential.h | 15 +-
src/potential/point_mass_softened/potential.h | 13 +
26 files changed, 1353 insertions(+), 7 deletions(-)
create mode 100644 doc/RTD/source/ExternalPotentials/index.rst
create mode 100755 examples/Hernquist_circularorbit/hernquistcirc.yml
create mode 100755 examples/Hernquist_circularorbit/makeIC.py
create mode 100755 examples/Hernquist_circularorbit/plotprog.py
create mode 100755 examples/Hernquist_circularorbit/run.sh
create mode 100644 examples/Hernquist_radialinfall/README
create mode 100644 examples/Hernquist_radialinfall/hernquist.yml
create mode 100644 examples/Hernquist_radialinfall/makeIC.py
create mode 100755 examples/Hernquist_radialinfall/plotprog.py
create mode 100755 examples/Hernquist_radialinfall/run.sh
create mode 100755 examples/NFW_Halo/README
create mode 100755 examples/NFW_Halo/makeIC.py
create mode 100755 examples/NFW_Halo/makePlots.py
create mode 100755 examples/NFW_Halo/run.sh
create mode 100755 examples/NFW_Halo/test.yml
create mode 100644 src/potential/hernquist/potential.h
create mode 100644 src/potential/nfw/potential.h
diff --git a/configure.ac b/configure.ac
index 345f732742..ab2ed97305 100644
--- a/configure.ac
+++ b/configure.ac
@@ -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])
;;
diff --git a/doc/RTD/source/ExternalPotentials/index.rst b/doc/RTD/source/ExternalPotentials/index.rst
new file mode 100644
index 0000000000..e313816246
--- /dev/null
+++ b/doc/RTD/source/ExternalPotentials/index.rst
@@ -0,0 +1,80 @@
+.. 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.
+
diff --git a/doc/RTD/source/index.rst b/doc/RTD/source/index.rst
index 343713f5f8..d148398c1b 100644
--- a/doc/RTD/source/index.rst
+++ b/doc/RTD/source/index.rst
@@ -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
diff --git a/examples/Hernquist_circularorbit/hernquistcirc.yml b/examples/Hernquist_circularorbit/hernquistcirc.yml
new file mode 100755
index 0000000000..5e81d18000
--- /dev/null
+++ b/examples/Hernquist_circularorbit/hernquistcirc.yml
@@ -0,0 +1,38 @@
+# 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)
diff --git a/examples/Hernquist_circularorbit/makeIC.py b/examples/Hernquist_circularorbit/makeIC.py
new file mode 100755
index 0000000000..474450f0e2
--- /dev/null
+++ b/examples/Hernquist_circularorbit/makeIC.py
@@ -0,0 +1,81 @@
+###############################################################################
+# 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()
diff --git a/examples/Hernquist_circularorbit/plotprog.py b/examples/Hernquist_circularorbit/plotprog.py
new file mode 100755
index 0000000000..849ae4e365
--- /dev/null
+++ b/examples/Hernquist_circularorbit/plotprog.py
@@ -0,0 +1,77 @@
+#!/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()
diff --git a/examples/Hernquist_circularorbit/run.sh b/examples/Hernquist_circularorbit/run.sh
new file mode 100755
index 0000000000..a2fedd0914
--- /dev/null
+++ b/examples/Hernquist_circularorbit/run.sh
@@ -0,0 +1,23 @@
+#!/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
diff --git a/examples/Hernquist_radialinfall/README b/examples/Hernquist_radialinfall/README
new file mode 100644
index 0000000000..be22a1a11a
--- /dev/null
+++ b/examples/Hernquist_radialinfall/README
@@ -0,0 +1,3 @@
+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.
diff --git a/examples/Hernquist_radialinfall/hernquist.yml b/examples/Hernquist_radialinfall/hernquist.yml
new file mode 100644
index 0000000000..adea54ed9a
--- /dev/null
+++ b/examples/Hernquist_radialinfall/hernquist.yml
@@ -0,0 +1,39 @@
+# 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
diff --git a/examples/Hernquist_radialinfall/makeIC.py b/examples/Hernquist_radialinfall/makeIC.py
new file mode 100644
index 0000000000..567e15a953
--- /dev/null
+++ b/examples/Hernquist_radialinfall/makeIC.py
@@ -0,0 +1,167 @@
+###############################################################################
+# 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()
diff --git a/examples/Hernquist_radialinfall/plotprog.py b/examples/Hernquist_radialinfall/plotprog.py
new file mode 100755
index 0000000000..d8de00a6b6
--- /dev/null
+++ b/examples/Hernquist_radialinfall/plotprog.py
@@ -0,0 +1,50 @@
+#!/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()
diff --git a/examples/Hernquist_radialinfall/run.sh b/examples/Hernquist_radialinfall/run.sh
new file mode 100755
index 0000000000..b2bcd99f9e
--- /dev/null
+++ b/examples/Hernquist_radialinfall/run.sh
@@ -0,0 +1,24 @@
+#!/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
diff --git a/examples/IsothermalPotential/energy_plot.py b/examples/IsothermalPotential/energy_plot.py
index dab30715fb..d157e4233c 100644
--- a/examples/IsothermalPotential/energy_plot.py
+++ b/examples/IsothermalPotential/energy_plot.py
@@ -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)
diff --git a/examples/IsothermalPotential/isothermal.yml b/examples/IsothermalPotential/isothermal.yml
index f9cae4d3b4..4f8d98a1f7 100644
--- a/examples/IsothermalPotential/isothermal.yml
+++ b/examples/IsothermalPotential/isothermal.yml
@@ -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
diff --git a/examples/NFW_Halo/README b/examples/NFW_Halo/README
new file mode 100755
index 0000000000..059d35c9a9
--- /dev/null
+++ b/examples/NFW_Halo/README
@@ -0,0 +1,5 @@
+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
+
diff --git a/examples/NFW_Halo/makeIC.py b/examples/NFW_Halo/makeIC.py
new file mode 100755
index 0000000000..68d8108f84
--- /dev/null
+++ b/examples/NFW_Halo/makeIC.py
@@ -0,0 +1,75 @@
+################################################################################
+# 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")
+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()
diff --git a/examples/NFW_Halo/makePlots.py b/examples/NFW_Halo/makePlots.py
new file mode 100755
index 0000000000..5e6f24d7a7
--- /dev/null
+++ b/examples/NFW_Halo/makePlots.py
@@ -0,0 +1,73 @@
+################################################################################
+# 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/>.
+#
+################################################################################
+from galpy.potential import NFWPotential
+from galpy.orbit import Orbit
+import numpy as np
+import matplotlib.pyplot as plt
+from astropy import units
+import h5py as h5
+
+C = 8.0
+M_200 = 2.0
+
+
+def read_data():
+ R = np.array([])
+ z = np.array([])
+ for frame in range(0, 599, 1):
+ try:
+ sim = h5.File("output_%04d.hdf5" % frame, "r")
+ except IOError:
+ break
+
+ boxSize = sim["/Header"].attrs["BoxSize"][0]
+ pos = sim["/PartType1/Coordinates"][:, :] - boxSize / 2.0
+ R = np.append(R, np.sqrt(pos[0, 0] ** 2 + pos[0, 1] ** 2))
+ z = np.append(z, pos[0, 2])
+ return (R, z)
+
+
+def galpy_nfw_orbit():
+ # Setting up the potential
+ nfw = NFWPotential(conc=C, mvir=M_200, H=70.0, wrtcrit=True, overdens=200)
+ nfw.turn_physical_on()
+ vxvv = [
+ 8.0 * units.kpc,
+ 0.0 * units.km / units.s,
+ 240.0 * units.km / units.s,
+ 0.0 * units.pc,
+ 5.0 * units.km / units.s,
+ ]
+
+ # Calculating the orbit
+ ts = np.linspace(0.0, 0.58, 1000) * units.Gyr
+ o = Orbit(vxvv=vxvv)
+ o.integrate(ts, nfw, method="odeint")
+
+ return o
+
+
+o = galpy_nfw_orbit()
+(R, z) = read_data()
+
+o.plot()
+plt.scatter(R, z, s=1, color="black", marker="x")
+plt.savefig("comparison.png")
+plt.close()
diff --git a/examples/NFW_Halo/run.sh b/examples/NFW_Halo/run.sh
new file mode 100755
index 0000000000..5501f09ef4
--- /dev/null
+++ b/examples/NFW_Halo/run.sh
@@ -0,0 +1,20 @@
+#!/bin/bash
+
+if [ ! -e test_nfw.hdf5 ]
+then
+ echo "Generate initial conditions for NFW example"
+ 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 test.yml 2>&1 | tee output.log
+
+if command -v python3 &>/dev/null; then
+ python3 makePlots.py
+else
+ python makePlots.py
+fi
diff --git a/examples/NFW_Halo/test.yml b/examples/NFW_Halo/test.yml
new file mode 100755
index 0000000000..73831af307
--- /dev/null
+++ b/examples/NFW_Halo/test.yml
@@ -0,0 +1,41 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.988e+33 # Solar mass
+ UnitLength_in_cgs: 3.086e+21 # kpc
+ UnitVelocity_in_cgs: 1e5 # km / s
+ 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: 0.6 # The end time of the simulation (in internal units).
+ dt_min: 1e-8 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-1 # 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: 1e-3 # Time between statistics output
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: test_nfw.hdf5 # The file to read
+ shift_x: 0. # (Optional) A shift to apply to all particles read from the ICs (in internal units).
+ shift_y: 0.
+ shift_z: 0.
+ periodic: 0
+
+# Isothermal potential parameters
+NFWPotential:
+ useabspos: 0
+ position: [0.0,0.0,0.0] # Location of centre of potential with respect to centre of the box (internal units)
+ concentration: 8.
+ M_200: 2.0e+12 # Virial mass (internal units)
+ critical_density: 140 # Critical density (internal units)
+ timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 64e31b97d0..63a61af3e7 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -175,6 +175,7 @@ EoS:
# Point mass external potentials
PointMassPotential:
+ useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
position: [50.,50.0,50.] # location of external point mass (internal units)
mass: 1e10 # mass of external point mass (internal units)
timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition
@@ -182,10 +183,29 @@ PointMassPotential:
# Isothermal potential parameters
IsothermalPotential:
+ useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
position: [100.,100.,100.] # Location of centre of isothermal potential with respect to centre of the box (internal units)
vrot: 200. # Rotation speed of isothermal potential (internal units)
timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition
epsilon: 0.1 # Softening size (internal units)
+
+# Hernquist potential parameters
+HernquistPotential:
+ useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
+ position: [100.,100.,100.] # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
+ mass: 1e10 # Mass of the Hernquist potential
+ scalelength: 10.0 # Scale length of the potential
+ 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
+ epsilon: 0.1 # Softening size (internal units)
+
+# Isothermal potential parameters
+NFWPotential:
+ useabspos: 0
+ position: [0.0,0.0,0.0] # Location of centre of isothermal potential with respect to centre of the box (internal units) if useabspos=0 otherwise with respect to the 0,0,0, coordinates.
+ concentration: 8. # Concentration of the halo
+ M_200: 2.0e+12 # Mass of the halo (M_200 in internal units)
+ critical_density: 127.4 # Critical density (internal units).
+ timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, basically determines fraction of orbital time we need to do an integration step
# Disk-patch potential parameters
DiscPatchPotential:
diff --git a/src/potential.h b/src/potential.h
index 814b83c691..59567fe922 100644
--- a/src/potential.h
+++ b/src/potential.h
@@ -34,6 +34,10 @@
#include "./potential/point_mass/potential.h"
#elif defined(EXTERNAL_POTENTIAL_ISOTHERMAL)
#include "./potential/isothermal/potential.h"
+#elif defined(EXTERNAL_POTENTIAL_HERNQUIST)
+#include "./potential/hernquist/potential.h"
+#elif defined(EXTERNAL_POTENTIAL_NFW)
+#include "./potential/nfw/potential.h"
#elif defined(EXTERNAL_POTENTIAL_DISC_PATCH)
#include "./potential/disc_patch/potential.h"
#elif defined(EXTERNAL_POTENTIAL_SINE_WAVE)
diff --git a/src/potential/hernquist/potential.h b/src/potential/hernquist/potential.h
new file mode 100644
index 0000000000..d0ec339d91
--- /dev/null
+++ b/src/potential/hernquist/potential.h
@@ -0,0 +1,223 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_POTENTIAL_HERNQUIST_H
+#define SWIFT_POTENTIAL_HERNQUIST_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local includes. */
+#include "error.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "space.h"
+#include "units.h"
+
+/**
+ * @brief External Potential Properties - Hernquist potential
+ */
+struct external_potential {
+
+ /*! Position of the centre of potential */
+ double x[3];
+
+ /*! Mass of the halo */
+ double mass;
+
+ /*! Scale length (often as a, to prevent confusion with the cosmological
+ * scale-factor we use al) */
+ double al;
+
+ /*! Square of the softening length. Acceleration tends to zero within this
+ * distance from the origin */
+ double epsilon2;
+
+ /* Minimum timestep of the potential given by the timestep multiple
+ * times the orbital time at the softening length */
+ double mintime;
+
+ /*! Time-step condition pre-factor, is multiplied times the circular orbital
+ * time to get the time steps */
+ double timestep_mult;
+};
+
+/**
+ * @brief Computes the time-step in a Hernquist potential based on a
+ * fraction of the circular orbital time
+ *
+ * @param time The current time.
+ * @param potential The #external_potential used in the run.
+ * @param phys_const The physical constants in internal units.
+ * @param g Pointer to the g-particle data.
+ */
+__attribute__((always_inline)) INLINE static float external_gravity_timestep(
+ double time, const struct external_potential* restrict potential,
+ const struct phys_const* restrict phys_const,
+ const struct gpart* restrict g) {
+
+ const float G_newton = phys_const->const_newton_G;
+
+ /* Calculate the relative potential with respect to the centre of the
+ * potential */
+ const float dx = g->x[0] - potential->x[0];
+ const float dy = g->x[1] - potential->x[1];
+ const float dz = g->x[2] - potential->x[2];
+
+ /* calculate the radius */
+ const float r = sqrtf(dx * dx + dy * dy + dz * dz + potential->epsilon2);
+ const float sqrtgm_inv = 1.f / sqrtf(G_newton * potential->mass);
+
+ /* Calculate the circular orbital period */
+ const float period = 2.f * M_PI * sqrtf(r) * potential->al *
+ (1 + r / potential->al) * sqrtgm_inv;
+
+ /* Time-step as a fraction of the cirecular orbital time */
+ const float time_step = potential->timestep_mult * period;
+
+ return max(time_step, potential->mintime);
+}
+
+/**
+ * @brief Computes the gravitational acceleration from an Hernquist potential.
+ *
+ * Note that the accelerations are multiplied by Newton's G constant
+ * later on.
+ *
+ * a_x = - GM / (a+r)^2 * x/r
+ * a_y = - GM / (a+r)^2 * y/r
+ * a_z = - GM / (a+r)^2 * z/r
+ *
+ * @param time The current time.
+ * @param potential The #external_potential used in the run.
+ * @param phys_const The physical constants in internal units.
+ * @param g Pointer to the g-particle data.
+ */
+__attribute__((always_inline)) INLINE static void external_gravity_acceleration(
+ double time, const struct external_potential* potential,
+ const struct phys_const* const phys_const, struct gpart* g) {
+
+ /* Determine the position relative to the centre of the potential */
+ const float dx = g->x[0] - potential->x[0];
+ const float dy = g->x[1] - potential->x[1];
+ const float dz = g->x[2] - potential->x[2];
+
+ /* Calculate the acceleration */
+ const float r = sqrtf(dx * dx + dy * dy + dz * dz + potential->epsilon2);
+ const float r_plus_a_inv = 1.f / (r + potential->al);
+ const float r_plus_a_inv2 = r_plus_a_inv * r_plus_a_inv;
+ const float term = -potential->mass * r_plus_a_inv2 / r;
+
+ g->a_grav[0] += term * dx;
+ g->a_grav[1] += term * dy;
+ g->a_grav[2] += term * dz;
+}
+
+/**
+ * @brief Computes the gravitational potential energy of a particle in an
+ * Hernquist potential.
+ *
+ * phi = - GM/(r+a)
+ *
+ * @param time The current time (unused here).
+ * @param potential The #external_potential used in the run.
+ * @param phys_const Physical constants in internal units.
+ * @param g Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+external_gravity_get_potential_energy(
+ double time, const struct external_potential* potential,
+ const struct phys_const* const phys_const, const struct gpart* g) {
+
+ const float dx = g->x[0] - potential->x[0];
+ const float dy = g->x[1] - potential->x[1];
+ const float dz = g->x[2] - potential->x[2];
+ const float r = sqrtf(dx * dx + dy * dy + dz * dz);
+ const float r_plus_alinv = 1.f / (r + potential->al);
+ return -phys_const->const_newton_G * potential->mass * r_plus_alinv;
+}
+
+/**
+ * @brief Initialises the external potential properties in the internal system
+ * of units.
+ *
+ * @param parameter_file The parsed parameter file
+ * @param phys_const Physical constants in internal units
+ * @param us The current internal system of units
+ * @param potential The external potential properties to initialize
+ */
+static INLINE void potential_init_backend(
+ struct swift_params* parameter_file, const struct phys_const* phys_const,
+ const struct unit_system* us, const struct space* s,
+ struct external_potential* potential) {
+
+ /* Read in the position of the centre of potential */
+ parser_get_param_double_array(parameter_file, "HernquistPotential:position",
+ 3, potential->x);
+
+ /* Is the position absolute or relative to the centre of the box? */
+ const int useabspos =
+ parser_get_param_int(parameter_file, "HernquistPotential:useabspos");
+
+ if (!useabspos) {
+ potential->x[0] += s->dim[0] / 2.;
+ potential->x[1] += s->dim[1] / 2.;
+ potential->x[2] += s->dim[2] / 2.;
+ }
+
+ /* Read the other parameters of the model */
+ potential->mass =
+ parser_get_param_double(parameter_file, "HernquistPotential:mass");
+ potential->al =
+ parser_get_param_double(parameter_file, "HernquistPotential:scalelength");
+ potential->timestep_mult = parser_get_param_float(
+ parameter_file, "HernquistPotential:timestep_mult");
+ const float epsilon =
+ parser_get_param_double(parameter_file, "HernquistPotential:epsilon");
+ potential->epsilon2 = epsilon * epsilon;
+
+ /* Compute the minimal time-step. */
+ /* This is the circular orbital time at the softened radius */
+ const float sqrtgm = sqrtf(phys_const->const_newton_G * potential->mass);
+ potential->mintime = 2.f * sqrtf(epsilon) * potential->al * M_PI *
+ (1 + epsilon / potential->al) / sqrtgm *
+ potential->timestep_mult;
+}
+
+/**
+ * @brief prints the properties of the external potential to stdout.
+ *
+ * @param potential the external potential properties.
+ */
+static inline void potential_print_backend(
+ const struct external_potential* potential) {
+
+ message(
+ "external potential is 'hernquist' with properties are (x,y,z) = (%e, "
+ "%e, %e), mass = %e "
+ "scale length = %e , minimum time = %e "
+ "timestep multiplier = %e",
+ potential->x[0], potential->x[1], potential->x[2], potential->mass,
+ potential->al, potential->mintime, potential->timestep_mult);
+}
+
+#endif /* SWIFT_POTENTIAL_HERNQUIST_H */
diff --git a/src/potential/isothermal/potential.h b/src/potential/isothermal/potential.h
index b5f8d7c397..160372210e 100644
--- a/src/potential/isothermal/potential.h
+++ b/src/potential/isothermal/potential.h
@@ -148,7 +148,7 @@ external_gravity_get_potential_energy(
const float dy = g->x[1] - potential->x[1];
const float dz = g->x[2] - potential->x[2];
- return -0.5f * potential->vrot * potential->vrot *
+ return 0.5f * potential->vrot * potential->vrot *
logf(dx * dx + dy * dy + dz * dz + potential->epsilon2);
}
@@ -166,11 +166,19 @@ static INLINE void potential_init_backend(
const struct unit_system* us, const struct space* s,
struct external_potential* potential) {
+ /* Read in the position of the centre of potential */
parser_get_param_double_array(parameter_file, "IsothermalPotential:position",
3, potential->x);
- potential->x[0] += s->dim[0] / 2.;
- potential->x[1] += s->dim[1] / 2.;
- potential->x[2] += s->dim[2] / 2.;
+
+ /* Is the position absolute or relative to the centre of the box? */
+ const int useabspos =
+ parser_get_param_int(parameter_file, "IsothermalPotential:useabspos");
+
+ if (!useabspos) {
+ potential->x[0] += s->dim[0] / 2.;
+ potential->x[1] += s->dim[1] / 2.;
+ potential->x[2] += s->dim[2] / 2.;
+ }
potential->vrot =
parser_get_param_double(parameter_file, "IsothermalPotential:vrot");
diff --git a/src/potential/nfw/potential.h b/src/potential/nfw/potential.h
new file mode 100644
index 0000000000..28bafd439a
--- /dev/null
+++ b/src/potential/nfw/potential.h
@@ -0,0 +1,260 @@
+/*******************************************************************************
+ * 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_POTENTIAL_NFW_H
+#define SWIFT_POTENTIAL_NFW_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <float.h>
+#include <math.h>
+
+/* Local includes. */
+#include "error.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "space.h"
+#include "units.h"
+
+/**
+ * @brief External Potential Properties - NFW Potential
+ rho(r) = rho_0 / ( (r/R_s)*(1+r/R_s)^2 )
+
+ We however parameterise this in terms of c and virial_mass
+ */
+struct external_potential {
+
+ /*! Position of the centre of potential */
+ double x[3];
+
+ /*! The scale radius of the NFW potential */
+ double r_s;
+
+ /*! The pre-factor \f$ 4 \pi G \rho_0 \r_s^3 \f$ */
+ double pre_factor;
+
+ /*! The critical density of the universe */
+ double rho_c;
+
+ /*! The concentration parameter */
+ double c_200;
+
+ /*! The virial mass */
+ double M_200;
+
+ /*! Time-step condition pre_factor, this factor is used to multiply times the
+ * orbital time, so in the case of 0.01 we take 1% of the orbital time as
+ * the time integration steps */
+ double timestep_mult;
+
+ /*! Minimum time step based on the orbital time at the softening times
+ * the timestep_mult */
+ double mintime;
+
+ /*! Common log term \f$ \ln(1+c_{200}) - \frac{c_{200}}{1 + c_{200}} \f$ */
+ double log_c200_term;
+
+ /*! Softening length */
+ double eps;
+};
+
+/**
+ * @brief Computes the time-step due to the acceleration from the NFW potential
+ * as a fraction (timestep_mult) of the circular orbital time of that
+ * particle.
+ *
+ * @param time The current time.
+ * @param potential The #external_potential used in the run.
+ * @param phys_const The physical constants in internal units.
+ * @param g Pointer to the g-particle data.
+ */
+__attribute__((always_inline)) INLINE static float external_gravity_timestep(
+ double time, const struct external_potential* restrict potential,
+ const struct phys_const* restrict phys_const,
+ const struct gpart* restrict g) {
+
+ const float dx = g->x[0] - potential->x[0];
+ const float dy = g->x[1] - potential->x[1];
+ const float dz = g->x[2] - potential->x[2];
+
+ const float r =
+ sqrtf(dx * dx + dy * dy + dz * dz + potential->eps * potential->eps);
+
+ const float mr = potential->M_200 *
+ (logf(1.f + r / potential->r_s) - r / (r + potential->r_s)) /
+ potential->log_c200_term;
+
+ const float period =
+ 2 * M_PI * r * sqrtf(r / (phys_const->const_newton_G * mr));
+
+ /* Time-step as a fraction of the circular period */
+ const float time_step = potential->timestep_mult * period;
+
+ return max(time_step, potential->mintime);
+}
+
+/**
+ * @brief Computes the gravitational acceleration from an NFW Halo potential.
+ *
+ * Note that the accelerations are multiplied by Newton's G constant
+ * later on.
+ *
+ * a_x = 4 pi \rho_0 r_s^3 ( 1/((r+rs)*r^2) - log(1+r/rs)/r^3) * x
+ * a_y = 4 pi \rho_0 r_s^3 ( 1/((r+rs)*r^2) - log(1+r/rs)/r^3) * y
+ * a_z = 4 pi \rho_0 r_s^3 ( 1/((r+rs)*r^2) - log(1+r/rs)/r^3) * z
+ *
+ * @param time The current time.
+ * @param potential The #external_potential used in the run.
+ * @param phys_const The physical constants in internal units.
+ * @param g Pointer to the g-particle data.
+ */
+__attribute__((always_inline)) INLINE static void external_gravity_acceleration(
+ double time, const struct external_potential* restrict potential,
+ const struct phys_const* restrict phys_const, struct gpart* restrict g) {
+
+ const float dx = g->x[0] - potential->x[0];
+ const float dy = g->x[1] - potential->x[1];
+ const float dz = g->x[2] - potential->x[2];
+
+ const float r =
+ sqrtf(dx * dx + dy * dy + dz * dz + potential->eps * potential->eps);
+ const float term1 = potential->pre_factor;
+ const float term2 = (1.0f / ((r + potential->r_s) * r * r) -
+ logf(1.0f + r / potential->r_s) / (r * r * r));
+
+ g->a_grav[0] += term1 * term2 * dx;
+ g->a_grav[1] += term1 * term2 * dy;
+ g->a_grav[2] += term1 * term2 * dz;
+}
+
+/**
+ * @brief Computes the gravitational potential energy of a particle in an
+ * NFW potential.
+ *
+ * phi = -4 * pi * G * rho_0 * r_s^3 * ln(1+r/r_s)
+ *
+ * @param time The current time (unused here).
+ * @param potential The #external_potential used in the run.
+ * @param phys_const Physical constants in internal units.
+ * @param g Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+external_gravity_get_potential_energy(
+ double time, const struct external_potential* potential,
+ const struct phys_const* const phys_const, const struct gpart* g) {
+
+ const float dx = g->x[0] - potential->x[0];
+ const float dy = g->x[1] - potential->x[1];
+ const float dz = g->x[2] - potential->x[2];
+
+ const float r =
+ sqrtf(dx * dx + dy * dy + dz * dz + potential->eps * potential->eps);
+ const float term1 = -potential->pre_factor / r;
+ const float term2 = logf(1.0f + r / potential->r_s);
+
+ return term1 * term2;
+}
+
+/**
+ * @brief Initialises the external potential properties in the internal system
+ * of units.
+ *
+ * @param parameter_file The parsed parameter file
+ * @param phys_const Physical constants in internal units
+ * @param us The current internal system of units
+ * @param potential The external potential properties to initialize
+ */
+static INLINE void potential_init_backend(
+ struct swift_params* parameter_file, const struct phys_const* phys_const,
+ const struct unit_system* us, const struct space* s,
+ struct external_potential* potential) {
+
+ /* Read in the position of the centre of potential */
+ parser_get_param_double_array(parameter_file, "NFWPotential:position", 3,
+ potential->x);
+
+ /* Is the position absolute or relative to the centre of the box? */
+ const int useabspos =
+ parser_get_param_int(parameter_file, "NFWPotential:useabspos");
+
+ if (!useabspos) {
+ potential->x[0] += s->dim[0] / 2.;
+ potential->x[1] += s->dim[1] / 2.;
+ potential->x[2] += s->dim[2] / 2.;
+ }
+
+ /* Read the other parameters of the model */
+ potential->timestep_mult =
+ parser_get_param_double(parameter_file, "NFWPotential:timestep_mult");
+ potential->c_200 =
+ parser_get_param_double(parameter_file, "NFWPotential:concentration");
+ potential->M_200 =
+ parser_get_param_double(parameter_file, "NFWPotential:M_200");
+ potential->rho_c =
+ parser_get_param_double(parameter_file, "NFWPotential:critical_density");
+ potential->eps = 0.05;
+
+ /* Compute R_200 */
+ const double R_200 =
+ cbrtf(3.0 * potential->M_200 / (4. * M_PI * 200.0 * potential->rho_c));
+
+ /* NFW scale-radius */
+ potential->r_s = R_200 / potential->c_200;
+ const double r_s3 = potential->r_s * potential->r_s * potential->r_s;
+
+ /* Log(c_200) term appearing in many expressions */
+ potential->log_c200_term =
+ log(1. + potential->c_200) - potential->c_200 / (1. + potential->c_200);
+
+ const double rho_0 =
+ potential->M_200 / (4.f * M_PI * r_s3 * potential->log_c200_term);
+
+ /* Pre-factor for the accelerations (note G is multiplied in later on) */
+ potential->pre_factor = 4.0f * M_PI * rho_0 * r_s3;
+
+ /* Compute the orbital time at the softening radius */
+ const double sqrtgm = sqrt(phys_const->const_newton_G * potential->M_200);
+ const double epslnthing = log(1.f + potential->eps / potential->r_s) -
+ potential->eps / (potential->eps + potential->r_s);
+
+ potential->mintime = 2. * M_PI * potential->eps * sqrtf(potential->eps) *
+ sqrtf(potential->log_c200_term / epslnthing) / sqrtgm *
+ potential->timestep_mult;
+}
+
+/**
+ * @brief Prints the properties of the external potential to stdout.
+ *
+ * @param potential The external potential properties.
+ */
+static INLINE void potential_print_backend(
+ const struct external_potential* potential) {
+
+ message(
+ "External potential is 'NFW' with properties are (x,y,z) = (%e, "
+ "%e, %e), scale radius = %e "
+ "timestep multiplier = %e, mintime = %e",
+ potential->x[0], potential->x[1], potential->x[2], potential->r_s,
+ potential->timestep_mult, potential->mintime);
+}
+
+#endif /* SWIFT_POTENTIAL_NFW_H */
diff --git a/src/potential/point_mass/potential.h b/src/potential/point_mass/potential.h
index f9d56a1ff1..5ae03f8637 100644
--- a/src/potential/point_mass/potential.h
+++ b/src/potential/point_mass/potential.h
@@ -137,7 +137,7 @@ external_gravity_get_potential_energy(
const float dx = g->x[0] - potential->x[0];
const float dy = g->x[1] - potential->x[1];
const float dz = g->x[2] - potential->x[2];
- const float rinv = 1. / sqrtf(dx * dx + dy * dy + dz * dz);
+ const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
return -phys_const->const_newton_G * potential->mass * rinv;
}
@@ -156,8 +156,21 @@ static INLINE void potential_init_backend(
const struct unit_system* us, const struct space* s,
struct external_potential* potential) {
+ /* Read in the position of the centre of potential */
parser_get_param_double_array(parameter_file, "PointMassPotential:position",
3, potential->x);
+
+ /* Is the position absolute or relative to the centre of the box? */
+ const int useabspos =
+ parser_get_param_int(parameter_file, "PointMassPotential:useabspos");
+
+ if (!useabspos) {
+ potential->x[0] += s->dim[0] / 2.;
+ potential->x[1] += s->dim[1] / 2.;
+ potential->x[2] += s->dim[2] / 2.;
+ }
+
+ /* Read the other parameters of the model */
potential->mass =
parser_get_param_double(parameter_file, "PointMassPotential:mass");
potential->timestep_mult = parser_get_param_float(
diff --git a/src/potential/point_mass_softened/potential.h b/src/potential/point_mass_softened/potential.h
index 0e35e7bb98..050bc1a00c 100644
--- a/src/potential/point_mass_softened/potential.h
+++ b/src/potential/point_mass_softened/potential.h
@@ -183,8 +183,21 @@ static INLINE void potential_init_backend(
const struct unit_system* us, const struct space* s,
struct external_potential* potential) {
+ /* Read in the position of the centre of potential */
parser_get_param_double_array(parameter_file, "PointMassPotential:position",
3, potential->x);
+
+ /* Is the position absolute or relative to the centre of the box? */
+ const int useabspos =
+ parser_get_param_int(parameter_file, "PointMassPotential:useabspos");
+
+ if (!useabspos) {
+ potential->x[0] += s->dim[0] / 2.;
+ potential->x[1] += s->dim[1] / 2.;
+ potential->x[2] += s->dim[2] / 2.;
+ }
+
+ /* Read the other parameters of the model */
potential->mass =
parser_get_param_double(parameter_file, "PointMassPotential:mass");
potential->timestep_mult = parser_get_param_float(
--
GitLab