From 21a4773f9276479c6381f5a2a322c68efce7144d Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Thu, 8 Sep 2016 23:09:41 +0100
Subject: [PATCH] Re-factored the whole external potential mechanism to match
what is done for cooling. Renamed 'disk' to 'disc'.
---
doc/Doxyfile.in | 1 +
.../GravityOnly/README | 0
.../GravityOnly/disc-patch.yml} | 8 +-
.../GravityOnly/makeIC.py | 4 +-
examples/DiscPatch/GravityOnly/run.sh | 10 +
.../GravityOnly/test.pro | 2 +-
.../HydroStatic/README | 12 +-
.../HydroStatic/disc-patch-icc.yml} | 8 +-
.../HydroStatic/disc-patch.yml} | 8 +-
.../HydroStatic/dynamic.pro | 0
.../HydroStatic/getGlass.sh | 0
.../HydroStatic/makeIC.py | 8 +-
.../HydroStatic/test.pro | 2 +-
examples/Disk-Patch/GravityOnly/run.sh | 10 -
.../ExternalPointMass/externalPointMass.yml | 2 +-
examples/parameter_example.yml | 12 +-
src/Makefile.am | 4 +-
src/const.h | 5 +-
src/engine.h | 2 +-
src/gravity/Default/gravity.h | 57 +---
src/potential.c | 52 +++
src/potential.h | 54 +++
src/potential/disc_patch/potential.h | 185 ++++++++++
src/potential/isothermal/potential.h | 151 +++++++++
src/potential/none/potential.h | 97 ++++++
src/potential/point_mass/potential.h | 154 +++++++++
src/potentials.c | 130 -------
src/potentials.h | 319 ------------------
src/runner.c | 2 +-
src/swift.h | 2 +-
src/timestep.h | 10 +-
31 files changed, 752 insertions(+), 559 deletions(-)
rename examples/{Disk-Patch => DiscPatch}/GravityOnly/README (100%)
rename examples/{Disk-Patch/GravityOnly/disk-patch.yml => DiscPatch/GravityOnly/disc-patch.yml} (91%)
rename examples/{Disk-Patch => DiscPatch}/GravityOnly/makeIC.py (98%)
create mode 100755 examples/DiscPatch/GravityOnly/run.sh
rename examples/{Disk-Patch => DiscPatch}/GravityOnly/test.pro (99%)
rename examples/{Disk-Patch => DiscPatch}/HydroStatic/README (58%)
rename examples/{Disk-Patch/HydroStatic/disk-patch-icc.yml => DiscPatch/HydroStatic/disc-patch-icc.yml} (91%)
rename examples/{Disk-Patch/HydroStatic/disk-patch.yml => DiscPatch/HydroStatic/disc-patch.yml} (91%)
rename examples/{Disk-Patch => DiscPatch}/HydroStatic/dynamic.pro (100%)
rename examples/{Disk-Patch => DiscPatch}/HydroStatic/getGlass.sh (100%)
rename examples/{Disk-Patch => DiscPatch}/HydroStatic/makeIC.py (98%)
rename examples/{Disk-Patch => DiscPatch}/HydroStatic/test.pro (99%)
delete mode 100755 examples/Disk-Patch/GravityOnly/run.sh
create mode 100644 src/potential.c
create mode 100644 src/potential.h
create mode 100644 src/potential/disc_patch/potential.h
create mode 100644 src/potential/isothermal/potential.h
create mode 100644 src/potential/none/potential.h
create mode 100644 src/potential/point_mass/potential.h
delete mode 100644 src/potentials.c
delete mode 100644 src/potentials.h
diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in
index f444d05000..2a5aeba7d1 100644
--- a/doc/Doxyfile.in
+++ b/doc/Doxyfile.in
@@ -763,6 +763,7 @@ INPUT = @top_srcdir@ @top_srcdir@/src @top_srcdir@/tests @top_
INPUT += @top_srcdir@/src/hydro/Minimal
INPUT += @top_srcdir@/src/gravity/Default
INPUT += @top_srcdir@/src/riemann
+INPUT += @top_srcdir@/src/potential/point_mass
INPUT += @top_srcdir@/src/cooling/const_du
# This tag can be used to specify the character encoding of the source files
diff --git a/examples/Disk-Patch/GravityOnly/README b/examples/DiscPatch/GravityOnly/README
similarity index 100%
rename from examples/Disk-Patch/GravityOnly/README
rename to examples/DiscPatch/GravityOnly/README
diff --git a/examples/Disk-Patch/GravityOnly/disk-patch.yml b/examples/DiscPatch/GravityOnly/disc-patch.yml
similarity index 91%
rename from examples/Disk-Patch/GravityOnly/disk-patch.yml
rename to examples/DiscPatch/GravityOnly/disc-patch.yml
index 78b42e7835..c76e4f6122 100644
--- a/examples/Disk-Patch/GravityOnly/disk-patch.yml
+++ b/examples/DiscPatch/GravityOnly/disc-patch.yml
@@ -19,7 +19,7 @@ Statistics:
# Parameters governing the snapshots
Snapshots:
- basename: Disk-Patch # Common part of the name of output files
+ basename: Disc-Patch # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 8. # Time difference between consecutive outputs (in internal units)
@@ -33,11 +33,11 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
- file_name: Disk-Patch.hdf5 # The file to read
+ file_name: Disc-Patch.hdf5 # The file to read
# External potential parameters
-Disk-PatchPotential:
+DiscPatchPotential:
surface_density: 10.
scale_height: 100.
- z_disk: 300.
+ z_disc: 300.
timestep_mult: 0.03
diff --git a/examples/Disk-Patch/GravityOnly/makeIC.py b/examples/DiscPatch/GravityOnly/makeIC.py
similarity index 98%
rename from examples/Disk-Patch/GravityOnly/makeIC.py
rename to examples/DiscPatch/GravityOnly/makeIC.py
index 702a50ff53..42cd26e235 100644
--- a/examples/Disk-Patch/GravityOnly/makeIC.py
+++ b/examples/DiscPatch/GravityOnly/makeIC.py
@@ -26,7 +26,7 @@ import random
# Generates N particles in a box of [0:BoxSize,0:BoxSize,-2scale_height:2scale_height]
# see Creasey, Theuns & Bower, 2013, for the equations:
-# disk parameters are: surface density sigma
+# disc parameters are: surface density sigma
# scale height b
# density: rho(z) = (sigma/2b) sech^2(z/b)
# isothermal velocity dispersion = <v_z^2? = b pi G sigma
@@ -79,7 +79,7 @@ N = int(sys.argv[1]) # Number of particles
rho = 2. # Density
P = 1. # Pressure
gamma = 5./3. # Gas adiabatic index
-fileName = "Disk-Patch.hdf5"
+fileName = "Disc-Patch.hdf5"
#---------------------------------------------------
diff --git a/examples/DiscPatch/GravityOnly/run.sh b/examples/DiscPatch/GravityOnly/run.sh
new file mode 100755
index 0000000000..9af1011ee6
--- /dev/null
+++ b/examples/DiscPatch/GravityOnly/run.sh
@@ -0,0 +1,10 @@
+#!/bin/bash
+
+# Generate the initial conditions if they are not present.
+if [ ! -e Disc-Patch.hdf5 ]
+then
+ echo "Generating initial conditions for the disc-patch example..."
+ python makeIC.py 1000
+fi
+
+../../swift -g -t 2 disc-patch.yml
diff --git a/examples/Disk-Patch/GravityOnly/test.pro b/examples/DiscPatch/GravityOnly/test.pro
similarity index 99%
rename from examples/Disk-Patch/GravityOnly/test.pro
rename to examples/DiscPatch/GravityOnly/test.pro
index 4bd0d00197..04e0afdf7e 100644
--- a/examples/Disk-Patch/GravityOnly/test.pro
+++ b/examples/DiscPatch/GravityOnly/test.pro
@@ -8,7 +8,7 @@ iplot = 1 ; if iplot = 1, make plot of E/Lz conservation, else, simply compare f
@physunits
indir = './'
-basefile = 'Disk-Patch_'
+basefile = 'Disc-Patch_'
; set properties of potential
uL = phys.pc ; unit of length
diff --git a/examples/Disk-Patch/HydroStatic/README b/examples/DiscPatch/HydroStatic/README
similarity index 58%
rename from examples/Disk-Patch/HydroStatic/README
rename to examples/DiscPatch/HydroStatic/README
index 56578731eb..42853e6b51 100644
--- a/examples/Disk-Patch/HydroStatic/README
+++ b/examples/DiscPatch/HydroStatic/README
@@ -1,4 +1,4 @@
-Generates and evolves a disk-patch, where gas is in hydrostatic
+Generates and evolves a disc-patch, where gas is in hydrostatic
equilibrium with an imposed external gravitational force, using the
equations from Creasey, Theuns & Bower, 2013, MNRAS, Volume 429,
Issue 3, p.1922-1948.
@@ -10,11 +10,11 @@ To generate ICs ready for a scientific run:
2) Generate pre-ICs by running the 'makeIC.py' script.
3) Run SWIFT with an isothermal EoS, no cooling nor feedback, and the
-disk-patch potential switched on and using the parameters from
-'disk-patch-icc.yml'
+disc-patch potential switched on and using the parameters from
+'disc-patch-icc.yml'
4) The ICs are then ready to be run for a science problem. Rename the last
-output to 'Disk-Patch-dynamic.hdf5'. These are now the ICs for the actual test.
+output to 'Disc-Patch-dynamic.hdf5'. These are now the ICs for the actual test.
-When running SWIFT with the parameters from 'disk-patch.yml' and an
-ideal gas EoS on these ICs the disk should stay in equilibrium.
+When running SWIFT with the parameters from 'disc-patch.yml' and an
+ideal gas EoS on these ICs the disc should stay in equilibrium.
diff --git a/examples/Disk-Patch/HydroStatic/disk-patch-icc.yml b/examples/DiscPatch/HydroStatic/disc-patch-icc.yml
similarity index 91%
rename from examples/Disk-Patch/HydroStatic/disk-patch-icc.yml
rename to examples/DiscPatch/HydroStatic/disc-patch-icc.yml
index ebf0467585..6a27016b8a 100644
--- a/examples/Disk-Patch/HydroStatic/disk-patch-icc.yml
+++ b/examples/DiscPatch/HydroStatic/disc-patch-icc.yml
@@ -19,7 +19,7 @@ Statistics:
# Parameters governing the snapshots
Snapshots:
- basename: Disk-Patch # Common part of the name of output files
+ basename: Disc-Patch # Common part of the name of output files
time_first: 0. # Time of the first output (in internal units)
delta_time: 12. # Time difference between consecutive outputs (in internal units)
@@ -33,12 +33,12 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
- file_name: Disk-Patch.hdf5 # The file to read
+ file_name: Disc-Patch.hdf5 # The file to read
# External potential parameters
-Disk-PatchPotential:
+DiscPatchPotential:
surface_density: 10.
scale_height: 100.
- z_disk: 200.
+ z_disc: 200.
timestep_mult: 0.03
growth_time: 5.
diff --git a/examples/Disk-Patch/HydroStatic/disk-patch.yml b/examples/DiscPatch/HydroStatic/disc-patch.yml
similarity index 91%
rename from examples/Disk-Patch/HydroStatic/disk-patch.yml
rename to examples/DiscPatch/HydroStatic/disc-patch.yml
index 55df81a08d..8bd67c5b08 100644
--- a/examples/Disk-Patch/HydroStatic/disk-patch.yml
+++ b/examples/DiscPatch/HydroStatic/disc-patch.yml
@@ -19,7 +19,7 @@ Statistics:
# Parameters governing the snapshots
Snapshots:
- basename: Disk-Patch-dynamic # Common part of the name of output files
+ basename: Disc-Patch-dynamic # Common part of the name of output files
time_first: 968. # Time of the first output (in internal units)
delta_time: 24. # Time difference between consecutive outputs (in internal units)
@@ -33,11 +33,11 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
- file_name: Disk-Patch-dynamic.hdf5 # The file to read
+ file_name: Disc-Patch-dynamic.hdf5 # The file to read
# External potential parameters
-Disk-PatchPotential:
+DiscPatchPotential:
surface_density: 10.
scale_height: 100.
- z_disk: 200.
+ z_disc: 200.
timestep_mult: 0.03
diff --git a/examples/Disk-Patch/HydroStatic/dynamic.pro b/examples/DiscPatch/HydroStatic/dynamic.pro
similarity index 100%
rename from examples/Disk-Patch/HydroStatic/dynamic.pro
rename to examples/DiscPatch/HydroStatic/dynamic.pro
diff --git a/examples/Disk-Patch/HydroStatic/getGlass.sh b/examples/DiscPatch/HydroStatic/getGlass.sh
similarity index 100%
rename from examples/Disk-Patch/HydroStatic/getGlass.sh
rename to examples/DiscPatch/HydroStatic/getGlass.sh
diff --git a/examples/Disk-Patch/HydroStatic/makeIC.py b/examples/DiscPatch/HydroStatic/makeIC.py
similarity index 98%
rename from examples/Disk-Patch/HydroStatic/makeIC.py
rename to examples/DiscPatch/HydroStatic/makeIC.py
index 40b1d1f1ff..48cc578658 100644
--- a/examples/Disk-Patch/HydroStatic/makeIC.py
+++ b/examples/DiscPatch/HydroStatic/makeIC.py
@@ -25,9 +25,9 @@ import math
import random
import matplotlib.pyplot as plt
-# Generates a disk-patch in hydrostatic equilibrium
+# Generates a disc-patch in hydrostatic equilibrium
# see Creasey, Theuns & Bower, 2013, for the equations:
-# disk parameters are: surface density sigma
+# disc parameters are: surface density sigma
# scale height b
# density: rho(z) = (sigma/2b) sech^2(z/b)
# isothermal velocity dispersion = <v_z^2? = b pi G sigma
@@ -79,7 +79,7 @@ Radius = 100. # maximum radius of particles [kpc]
G = const_G
# File
-fileName = "Disk-Patch.hdf5"
+fileName = "Disc-Patch.hdf5"
#---------------------------------------------------
mass = 1
@@ -145,7 +145,7 @@ mass = 0.*h + pmass
entropy_flag = 0
vel = 0 + 0 * pos
-# move centre of disk to middle of box
+# move centre of disc to middle of box
pos[:,:] += boxSize/2
diff --git a/examples/Disk-Patch/HydroStatic/test.pro b/examples/DiscPatch/HydroStatic/test.pro
similarity index 99%
rename from examples/Disk-Patch/HydroStatic/test.pro
rename to examples/DiscPatch/HydroStatic/test.pro
index 31e027e3a3..950aebc65d 100644
--- a/examples/Disk-Patch/HydroStatic/test.pro
+++ b/examples/DiscPatch/HydroStatic/test.pro
@@ -8,7 +8,7 @@ iplot = 1 ; if iplot = 1, make plot of E/Lz conservation, else, simply compare f
@physunits
indir = './'
-basefile = 'Disk-Patch_'
+basefile = 'Disc-Patch_'
; set properties of potential
uL = phys.pc ; unit of length
diff --git a/examples/Disk-Patch/GravityOnly/run.sh b/examples/Disk-Patch/GravityOnly/run.sh
deleted file mode 100755
index a123ad24d7..0000000000
--- a/examples/Disk-Patch/GravityOnly/run.sh
+++ /dev/null
@@ -1,10 +0,0 @@
-#!/bin/bash
-
-# Generate the initial conditions if they are not present.
-if [ ! -e Isothermal.hdf5 ]
-then
- echo "Generating initial conditions for the disk-patch example..."
- python makeIC.py 1000
-fi
-
-../../swift -g -t 2 disk-patch.yml
diff --git a/examples/ExternalPointMass/externalPointMass.yml b/examples/ExternalPointMass/externalPointMass.yml
index d06c165651..ce300b3215 100644
--- a/examples/ExternalPointMass/externalPointMass.yml
+++ b/examples/ExternalPointMass/externalPointMass.yml
@@ -38,7 +38,7 @@ InitialConditions:
shift_z: 50.
# External potential parameters
-PointMass:
+PointMassPotential:
position_x: 50. # location of external point mass in internal units
position_y: 50.
position_z: 50.
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 10bd09a18f..892016d68a 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -66,7 +66,7 @@ DomainDecomposition:
# Parameters related to external potentials --------------------------------------------
# Point mass external potentials
-PointMass:
+PointMassPotential:
position_x: 50. # location of external point mass (internal units)
position_y: 50.
position_z: 50.
@@ -82,12 +82,12 @@ IsothermalPotential:
timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition
# Disk-patch potential parameters
-Disk-PatchPotential:
- surface_density: 10. # Surface density of the disk (internal units)
- scale_height: 100. # Scale height of the disk (internal units)
- z_disk: 200. # Disk height (internal units)
+DiscPatchPotential:
+ surface_density: 10. # Surface density of the disc (internal units)
+ scale_height: 100. # Scale height of the disc (internal units)
+ z_disc: 200. # Disc height (internal units)
timestep_mult: 0.03 # Dimensionless pre-factor for the time-step condition
- growth_time: 5. # (Optional) Time for the disk to grow to its final size (multiple of the dynamical time)
+ growth_time: 5. # (Optional) Time for the disc to grow to its final size (multiple of the dynamical time)
# Parameters related to cooling function ----------------------------------------------
diff --git a/src/Makefile.am b/src/Makefile.am
index 095b11272d..18ea18224a 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -42,7 +42,7 @@ endif
include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
engine.h swift.h serial_io.h timers.h debug.h scheduler.h proxy.h parallel_io.h \
common_io.h single_io.h multipole.h map.h tools.h partition.h clocks.h parser.h \
- physical_constants.h physical_constants_cgs.h potentials.h version.h \
+ physical_constants.h physical_constants_cgs.h potential.h version.h \
hydro_properties.h threadpool.h cooling.h cooling_struct.h
@@ -51,7 +51,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
serial_io.c timers.c debug.c scheduler.c proxy.c parallel_io.c \
units.c common_io.c single_io.c multipole.c version.c map.c \
kernel_hydro.c tools.c part.c partition.c clocks.c parser.c \
- physical_constants.c potentials.c hydro_properties.c \
+ physical_constants.c potential.c hydro_properties.c \
runner_doiact_fft.c threadpool.c cooling.c
# Include files for distribution, not installation.
diff --git a/src/const.h b/src/const.h
index 1082bccc0b..316a13c8eb 100644
--- a/src/const.h
+++ b/src/const.h
@@ -90,9 +90,10 @@
#define const_gravity_eta 0.025f
/* External gravity properties */
-#define EXTERNAL_POTENTIAL_POINTMASS
+#define EXTERNAL_POTENTIAL_NONE
+//#define EXTERNAL_POTENTIAL_POINTMASS
//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
-//#define EXTERNAL_POTENTIAL_DISK_PATCH
+//#define EXTERNAL_POTENTIAL_DISC_PATCH
/* Cooling properties */
#define COOLING_NONE
diff --git a/src/engine.h b/src/engine.h
index e043b7eb76..d36914af61 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -41,7 +41,7 @@
#include "cooling_struct.h"
#include "parser.h"
#include "partition.h"
-#include "potentials.h"
+#include "potential.h"
#include "runner.h"
#include "scheduler.h"
#include "space.h"
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 2415e20ac5..0f13f9557b 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -22,37 +22,7 @@
#include <float.h>
#include "minmax.h"
-#include "potentials.h"
-
-/**
- * @brief Computes the gravity time-step of a given particle due to an external
- *potential.
- *
- * This function only branches towards the potential chosen by the user.
- *
- * @param potential The properties of the external potential.
- * @param phys_const The physical constants in internal units.
- * @param gp Pointer to the g-particle data.
- */
-__attribute__((always_inline)) INLINE static float
-gravity_compute_timestep_external(const struct external_potential* potential,
- const struct phys_const* const phys_const,
- const struct gpart* const gp) {
-
- float dt = FLT_MAX;
-
-#ifdef EXTERNAL_POTENTIAL_POINTMASS
- dt = min(dt, external_gravity_pointmass_timestep(potential, phys_const, gp));
-#endif
-#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
- dt = min(dt, external_gravity_isothermalpotential_timestep(potential,
- phys_const, gp));
-#endif
-#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
- dt = min(dt, external_gravity_disk_patch_timestep(potential, phys_const, gp));
-#endif
- return dt;
-}
+#include "physical_constants.h"
/**
* @brief Computes the gravity time-step of a given particle due to self-gravity
@@ -125,31 +95,6 @@ __attribute__((always_inline)) INLINE static void gravity_end_force(
gp->a_grav[2] *= const_G;
}
-/**
- * @brief Computes the gravitational acceleration induced by external potentials
- *
- * This function only branches towards the potential chosen by the user.
- *
- * @param time The current time in internal units.
- * @param potential The properties of the external potential.
- * @param phys_const The physical constants in internal units.
- * @param gp The particle to act upon.
- */
-__attribute__((always_inline)) INLINE static void external_gravity(
- double time, const struct external_potential* potential,
- const struct phys_const* const phys_const, struct gpart* gp) {
-
-#ifdef EXTERNAL_POTENTIAL_POINTMASS
- external_gravity_pointmass(potential, phys_const, gp);
-#endif
-#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
- external_gravity_isothermalpotential(potential, phys_const, gp);
-#endif
-#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
- external_gravity_disk_patch_potential(time, potential, phys_const, gp);
-#endif
-}
-
/**
* @brief Kick the additional variables
*
diff --git a/src/potential.c b/src/potential.c
new file mode 100644
index 0000000000..5433a05e3e
--- /dev/null
+++ b/src/potential.c
@@ -0,0 +1,52 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * 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/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* This object's header. */
+#include "potential.h"
+
+/**
+ * @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
+ */
+void potential_init(const struct swift_params* parameter_file,
+ const struct phys_const* phys_const,
+ const struct UnitSystem* us,
+ struct external_potential* potential) {
+
+ potential_init_backend(parameter_file, phys_const, us, potential);
+}
+
+/**
+ * @brief Prints the properties of the external potential to stdout.
+ *
+ * @param potential The external potential properties.
+ */
+void potential_print(const struct external_potential* potential) {
+
+ potential_print_backend(potential);
+}
diff --git a/src/potential.h b/src/potential.h
new file mode 100644
index 0000000000..77bd41794a
--- /dev/null
+++ b/src/potential.h
@@ -0,0 +1,54 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * 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_H
+#define SWIFT_POTENTIAL_H
+
+/**
+ * @file src/potential.h
+ * @brief Branches between the different external gravitational potentials.
+ */
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local includes. */
+#include "const.h"
+
+/* Import the right external potential definition */
+#if defined(EXTERNAL_POTENTIAL_NONE)
+#include "./potential/none/potential.h"
+#elif defined(EXTERNAL_POTENTIAL_POINTMASS)
+#include "./potential/point_mass/potential.h"
+#elif defined(EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL)
+#include "./potential/isothermal/potential.h"
+#elif defined(EXTERNAL_POTENTIAL_DISC_PATCH)
+#include "./potential/disc_patch/potential.h"
+#else
+#error "Invalid choice of external potential"
+#endif
+
+/* Now, some generic functions, defined in the source file */
+void potential_init(const struct swift_params* parameter_file,
+ const struct phys_const* phys_const,
+ const struct UnitSystem* us,
+ struct external_potential* potential);
+
+void potential_print(const struct external_potential* potential);
+
+#endif /* SWIFT_POTENTIAL_H */
diff --git a/src/potential/disc_patch/potential.h b/src/potential/disc_patch/potential.h
new file mode 100644
index 0000000000..169e6d1e90
--- /dev/null
+++ b/src/potential/disc_patch/potential.h
@@ -0,0 +1,185 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * 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_DISC_PATCH_H
+#define SWIFT_DISC_PATCH_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <float.h>
+#include <math.h>
+
+/* Local includes. */
+#include "const.h"
+#include "error.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "units.h"
+
+/* External Potential Properties */
+struct external_potential {
+
+ double surface_density;
+ double scale_height;
+ double z_disc;
+ double dynamical_time;
+ double growth_time;
+ double timestep_mult;
+};
+
+/**
+ * @brief Computes the time-step from the acceleration due to a hydrostatic
+ * disc.
+ *
+ * See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
+ *
+ * @param time The current time.
+ * @param potential The properties of the potential.
+ * @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) {
+
+ /* initilize time step to disc dynamical time */
+ const float dt_dyn = potential->dynamical_time;
+ float dt = dt_dyn;
+
+ /* absolute value of height above disc */
+ const float dz = fabs(g->x[2] - potential->z_disc);
+
+ /* vertical cceleration */
+ const float z_accel = 2 * M_PI * phys_const->const_newton_G *
+ potential->surface_density *
+ tanh(dz / potential->scale_height);
+
+ /* demand that dt * velocity < fraction of scale height of disc */
+ float dt1 = FLT_MAX;
+ if (fabs(g->v_full[2]) > 0) {
+ dt1 = potential->scale_height / fabs(g->v_full[2]);
+ if (dt1 < dt) dt = dt1;
+ }
+
+ /* demand that dt^2 * acceleration < fraction of scale height of disc */
+ float dt2 = FLT_MAX;
+ if (fabs(z_accel) > 0) {
+ dt2 = potential->scale_height / fabs(z_accel);
+ if (dt2 < dt * dt) dt = sqrt(dt2);
+ }
+
+ /* demand that dt^3 jerk < fraction of scale height of disc */
+ float dt3 = FLT_MAX;
+ if (abs(g->v_full[2]) > 0) {
+ const float dz_accel_over_dt =
+ 2 * M_PI * phys_const->const_newton_G * potential->surface_density /
+ potential->scale_height / cosh(dz / potential->scale_height) /
+ cosh(dz / potential->scale_height) * fabs(g->v_full[2]);
+
+ dt3 = potential->scale_height / fabs(dz_accel_over_dt);
+ if (dt3 < dt * dt * dt) dt = pow(dt3, 1. / 3.);
+ }
+
+ return potential->timestep_mult * dt;
+}
+
+/**
+ * @brief Computes the gravitational acceleration along z due to a hydrostatic
+ * disc
+ *
+ * See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
+ *
+ * @param time The current time in internal units.
+ * @param potential The properties of the potential.
+ * @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 G_newton = phys_const->const_newton_G;
+ const float dz = g->x[2] - potential->z_disc;
+ const float t_dyn = potential->dynamical_time;
+
+ float reduction_factor = 1.;
+ if (time < potential->growth_time * t_dyn)
+ reduction_factor = time / (potential->growth_time * t_dyn);
+
+ const float z_accel = reduction_factor * 2 * G_newton * M_PI *
+ potential->surface_density *
+ tanh(fabs(dz) / potential->scale_height);
+
+ /* Accelerations. Note that they are multiplied by G later on */
+ if (dz > 0) g->a_grav[2] -= z_accel / G_newton;
+ if (dz < 0) g->a_grav[2] += z_accel / G_newton;
+}
+
+/**
+ * @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(
+ const struct swift_params* parameter_file,
+ const struct phys_const* phys_const, const struct UnitSystem* us,
+ struct external_potential* potential) {
+
+ potential->surface_density =
+ parser_get_param_double(parameter_file, "DiscPatchPotential:surface_density");
+ potential->scale_height =
+ parser_get_param_double(parameter_file, "DiscPatchPotential:scale_height");
+ potential->z_disc =
+ parser_get_param_double(parameter_file, "DiscPatchPotential:z_disc");
+ potential->timestep_mult =
+ parser_get_param_double(parameter_file, "DiscPatchPotential:timestep_mult");
+ potential->growth_time =
+ parser_get_opt_param_double(parameter_file, "DiscPatchPotential:growth_time", 0.);
+ potential->dynamical_time =
+ sqrt(potential->scale_height /
+ (phys_const->const_newton_G * potential->surface_density));
+}
+
+/**
+ * @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 'Disk-patch' with properties surface_density = %e "
+ "disc height= %e scale height = %e timestep multiplier = %e.",
+ potential->surface_density, potential->z_disc, potential->scale_height,
+ potential->timestep_mult);
+
+ if (potential->growth_time > 0.)
+ message("Disc will grow for %f dynamiccal times.", potential->growth_time);
+}
+
+#endif /* SWIFT_DISC_PATCH_H */
diff --git a/src/potential/isothermal/potential.h b/src/potential/isothermal/potential.h
new file mode 100644
index 0000000000..4a43f0d344
--- /dev/null
+++ b/src/potential/isothermal/potential.h
@@ -0,0 +1,151 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * 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_ISOTHERMAL_H
+#define SWIFT_POTENTIAL_ISOTHERMAL_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 "units.h"
+
+/* External Potential Properties */
+struct external_potential {
+
+ double x, y, z;
+ double vrot;
+ float timestep_mult;
+};
+
+/**
+ * @brief Computes the time-step due to the acceleration from an isothermal
+ * potential.
+ *
+ * @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;
+ const float dy = g->x[1] - potential->y;
+ const float dz = g->x[2] - potential->z;
+
+ const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
+ const float drdv =
+ dx * (g->v_full[0]) + dy * (g->v_full[1]) + dz * (g->v_full[2]);
+ const double vrot = potential->vrot;
+
+ const float dota_x =
+ vrot * vrot * rinv2 * (g->v_full[0] - 2 * drdv * dx * rinv2);
+ const float dota_y =
+ vrot * vrot * rinv2 * (g->v_full[1] - 2 * drdv * dy * rinv2);
+ const float dota_z =
+ vrot * vrot * rinv2 * (g->v_full[2] - 2 * drdv * dz * rinv2);
+ const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
+ const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] +
+ g->a_grav[2] * g->a_grav[2];
+
+ return potential->timestep_mult * sqrtf(a_2 / dota_2);
+}
+
+/**
+ * @brief Computes the gravitational acceleration from an isothermal potential.
+ *
+ * Note that the accelerations are multiplied by Newton's G constant
+ * later on.
+ *
+ * @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) {
+
+ const double G_newton = phys_const->const_newton_G;
+ const float dx = g->x[0] - potential->x;
+ const float dy = g->x[1] - potential->y;
+ const float dz = g->x[2] - potential->z;
+ const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
+
+ const double vrot = potential->vrot;
+ const double term = -vrot * vrot * rinv2 / G_newton;
+
+ g->a_grav[0] += term * dx;
+ g->a_grav[1] += term * dy;
+ g->a_grav[2] += term * dz;
+}
+
+/**
+ * @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(
+ const struct swift_params* parameter_file,
+ const struct phys_const* phys_const, const struct UnitSystem* us,
+ struct external_potential* potential) {
+
+ potential->x =
+ parser_get_param_double(parameter_file, "IsothermalPotential:position_x");
+ potential->y =
+ parser_get_param_double(parameter_file, "IsothermalPotential:position_y");
+ potential->z =
+ parser_get_param_double(parameter_file, "IsothermalPotential:position_z");
+ potential->vrot =
+ parser_get_param_double(parameter_file, "IsothermalPotential:vrot");
+ potential->timestep_mult = parser_get_param_float(
+ parameter_file, "IsothermalPotential: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 'Isothermal' with properties are (x,y,z) = (%e, "
+ "%e, %e), vrot = %e "
+ "timestep multiplier = %e.",
+ potential->x, potential->y, potential->z, potential->vrot,
+ potential->timestep_mult);
+}
+
+#endif /* SWIFT_POTENTIAL_ISOTHERMAL_H */
diff --git a/src/potential/none/potential.h b/src/potential/none/potential.h
new file mode 100644
index 0000000000..8b1c3e8415
--- /dev/null
+++ b/src/potential/none/potential.h
@@ -0,0 +1,97 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * 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_NONE_H
+#define SWIFT_POTENTIAL_NONE_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <float.h>
+
+/* Local includes. */
+#include "error.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "units.h"
+
+/**
+ * @brief External Potential Properties
+ */
+struct external_potential {};
+
+/**
+ * @brief Computes the time-step due to the acceleration from nothing
+ *
+ * @param time The current time.
+ * @param potential The properties of the externa potential.
+ * @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) {
+
+ return FLT_MAX;
+}
+
+/**
+ * @brief Computes the gravitational acceleration due to nothing
+ *
+ * We do nothing.
+ *
+ * @param time The current time.
+ * @param potential The proerties of the external potential.
+ * @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) {}
+
+/**
+ * @brief Initialises the external potential properties in the internal system
+ * of units.
+ *
+ * Nothing to do here.
+ *
+ * @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(
+ const struct swift_params* parameter_file,
+ const struct phys_const* phys_const, const struct UnitSystem* us,
+ struct external_potential* potential) {}
+
+/**
+ * @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 'No external potential'.");
+}
+
+#endif /* SWIFT_POTENTIAL_NONE_H */
diff --git a/src/potential/point_mass/potential.h b/src/potential/point_mass/potential.h
new file mode 100644
index 0000000000..7d794347cb
--- /dev/null
+++ b/src/potential/point_mass/potential.h
@@ -0,0 +1,154 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * 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_POINT_MASS_H
+#define SWIFT_POTENTIAL_POINT_MASS_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 "units.h"
+
+/**
+ * @brief External Potential Properties
+ */
+struct external_potential {
+
+ double x, y, z;
+ double mass;
+ float timestep_mult;
+};
+
+/**
+ * @brief Computes the time-step due to the acceleration from a point mass
+ *
+ * We pass in the time for simulations where the potential evolves with time.
+ *
+ * @param time The current time.
+ * @param potential The properties of the externa potential.
+ * @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;
+ const float dx = g->x[0] - potential->x;
+ const float dy = g->x[1] - potential->y;
+ const float dz = g->x[2] - potential->z;
+ const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
+ const float rinv2 = rinv * rinv;
+ const float rinv3 = rinv2 * rinv;
+ const float drdv = (g->x[0] - potential->x) * (g->v_full[0]) +
+ (g->x[1] - potential->y) * (g->v_full[1]) +
+ (g->x[2] - potential->z) * (g->v_full[2]);
+ const float dota_x = G_newton * potential->mass * rinv3 *
+ (-g->v_full[0] + 3.f * rinv2 * drdv * dx);
+ const float dota_y = G_newton * potential->mass * rinv3 *
+ (-g->v_full[1] + 3.f * rinv2 * drdv * dy);
+ const float dota_z = G_newton * potential->mass * rinv3 *
+ (-g->v_full[2] + 3.f * rinv2 * drdv * dz);
+ const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
+ const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] +
+ g->a_grav[2] * g->a_grav[2];
+
+ return potential->timestep_mult * sqrtf(a_2 / dota_2);
+}
+
+/**
+ * @brief Computes the gravitational acceleration of a particle due to a
+ * point mass
+ *
+ * Note that the accelerations are multiplied by Newton's G constant later
+ * on.
+ *
+ * We pass in the time for simulations where the potential evolves with time.
+ *
+ * @param time The current time.
+ * @param potential The proerties of the external potential.
+ * @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;
+ const float dy = g->x[1] - potential->y;
+ const float dz = g->x[2] - potential->z;
+ const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
+ const float rinv3 = rinv * rinv * rinv;
+
+ g->a_grav[0] += -potential->mass * dx * rinv3;
+ g->a_grav[1] += -potential->mass * dy * rinv3;
+ g->a_grav[2] += -potential->mass * dz * rinv3;
+}
+
+/**
+ * @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(
+ const struct swift_params* parameter_file,
+ const struct phys_const* phys_const, const struct UnitSystem* us,
+ struct external_potential* potential) {
+
+ potential->x =
+ parser_get_param_double(parameter_file, "PointMassPotential:position_x");
+ potential->y =
+ parser_get_param_double(parameter_file, "PointMassPotential:position_y");
+ potential->z =
+ parser_get_param_double(parameter_file, "PointMassPotential:position_z");
+ potential->mass =
+ parser_get_param_double(parameter_file, "PointMassPotential:mass");
+ potential->timestep_mult = parser_get_param_float(
+ parameter_file, "PointMassPotential: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 'Point mass' with properties (x,y,z) = (%e, %e, "
+ "%e), M = %e timestep multiplier = %e.",
+ potential->x, potential->y, potential->z, potential->mass,
+ potential->timestep_mult);
+}
+
+#endif /* SWIFT_POTENTIAL_POINT_MASS_H */
diff --git a/src/potentials.c b/src/potentials.c
deleted file mode 100644
index dd7aed8712..0000000000
--- a/src/potentials.c
+++ /dev/null
@@ -1,130 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
- * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- *
- * 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/>.
- *
- ******************************************************************************/
-
-/* Config parameters. */
-#include "../config.h"
-
-/* This object's header. */
-#include "potentials.h"
-
-/**
- * @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
- */
-void potential_init(const struct swift_params* parameter_file,
- const struct phys_const* phys_const,
- const struct UnitSystem* us,
- struct external_potential* potential) {
-
-#ifdef EXTERNAL_POTENTIAL_POINTMASS
-
- potential->point_mass.x =
- parser_get_param_double(parameter_file, "PointMass:position_x");
- potential->point_mass.y =
- parser_get_param_double(parameter_file, "PointMass:position_y");
- potential->point_mass.z =
- parser_get_param_double(parameter_file, "PointMass:position_z");
- potential->point_mass.mass =
- parser_get_param_double(parameter_file, "PointMass:mass");
- potential->point_mass.timestep_mult =
- parser_get_param_float(parameter_file, "PointMass:timestep_mult");
-
-#endif /* EXTERNAL_POTENTIAL_POINTMASS */
-
-#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
-
- potential->isothermal_potential.x =
- parser_get_param_double(parameter_file, "IsothermalPotential:position_x");
- potential->isothermal_potential.y =
- parser_get_param_double(parameter_file, "IsothermalPotential:position_y");
- potential->isothermal_potential.z =
- parser_get_param_double(parameter_file, "IsothermalPotential:position_z");
- potential->isothermal_potential.vrot =
- parser_get_param_double(parameter_file, "IsothermalPotential:vrot");
- potential->isothermal_potential.timestep_mult = parser_get_param_float(
- parameter_file, "IsothermalPotential:timestep_mult");
-
-#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
-#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
- potential->disk_patch_potential.surface_density = parser_get_param_double(
- parameter_file, "Disk-PatchPotential:surface_density");
- potential->disk_patch_potential.scale_height = parser_get_param_double(
- parameter_file, "Disk-PatchPotential:scale_height");
- potential->disk_patch_potential.z_disk =
- parser_get_param_double(parameter_file, "Disk-PatchPotential:z_disk");
- potential->disk_patch_potential.timestep_mult = parser_get_param_double(
- parameter_file, "Disk-PatchPotential:timestep_mult");
- potential->disk_patch_potential.growth_time = parser_get_opt_param_double(
- parameter_file, "Disk-PatchPotential:growth_time", 0.);
- potential->disk_patch_potential.dynamical_time =
- sqrt(potential->disk_patch_potential.scale_height /
- (phys_const->const_newton_G *
- potential->disk_patch_potential.surface_density));
-#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
-}
-
-/**
- * @brief Prints the properties of the external potential to stdout.
- *
- * @param potential The external potential properties.
- */
-void potential_print(const struct external_potential* potential) {
-
-#ifdef EXTERNAL_POTENTIAL_POINTMASS
-
- message(
- "Point mass properties are (x,y,z) = (%e, %e, %e), M = %e timestep "
- "multiplier = %e.",
- potential->point_mass.x, potential->point_mass.y, potential->point_mass.z,
- potential->point_mass.mass, potential->point_mass.timestep_mult);
-
-#endif /* EXTERNAL_POTENTIAL_POINTMASS */
-
-#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
-
- message(
- "Isothermal potential properties are (x,y,z) = (%e, %e, %e), vrot = %e "
- "timestep multiplier = %e.",
- potential->isothermal_potential.x, potential->isothermal_potential.y,
- potential->isothermal_potential.z, potential->isothermal_potential.vrot,
- potential->isothermal_potential.timestep_mult);
-
-#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
-
-#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
-
- message(
- "Disk-patch potential properties are surface_density = %e disk height= "
- "%e scale height = %e timestep multiplier = %e.",
- potential->disk_patch_potential.surface_density,
- potential->disk_patch_potential.z_disk,
- potential->disk_patch_potential.scale_height,
- potential->disk_patch_potential.timestep_mult);
-
- if (potential->disk_patch_potential.growth_time > 0.)
- message("Disk will grow for %f dynamiccal times.",
- potential->disk_patch_potential.growth_time);
-#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
-}
diff --git a/src/potentials.h b/src/potentials.h
deleted file mode 100644
index 74f0fd2856..0000000000
--- a/src/potentials.h
+++ /dev/null
@@ -1,319 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
- * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- *
- * 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_POTENTIALS_H
-#define SWIFT_POTENTIALS_H
-
-/* Config parameters. */
-#include "../config.h"
-
-/* Some standard headers. */
-#include <float.h>
-#include <math.h>
-
-/* Local includes. */
-#include "const.h"
-#include "error.h"
-#include "parser.h"
-#include "part.h"
-#include "physical_constants.h"
-#include "units.h"
-
-/* External Potential Properties */
-struct external_potential {
-
-#ifdef EXTERNAL_POTENTIAL_POINTMASS
- struct {
- double x, y, z;
- double mass;
- float timestep_mult;
- } point_mass;
-#endif
-
-#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
- struct {
- double x, y, z;
- double vrot;
- float timestep_mult;
- } isothermal_potential;
-#endif
-
-#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
- struct {
- double surface_density;
- double scale_height;
- double z_disk;
- double dynamical_time;
- double growth_time;
- double timestep_mult;
- } disk_patch_potential;
-#endif
-};
-
-/* ------------------------------------------------------------------------- */
-
-/* external potential: disk_patch */
-#ifdef EXTERNAL_POTENTIAL_DISK_PATCH
-
-/**
- * @brief Computes the time-step from the acceleration due to a hydrostatic
- * disk.
- *
- * See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
- *
- * @param potential The properties of the potential.
- * @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_disk_patch_timestep(
- const struct external_potential* restrict potential,
- const struct phys_const* restrict phys_const,
- const struct gpart* restrict g) {
-
- /* initilize time step to disk dynamical time */
- const float dt_dyn = potential->disk_patch_potential.dynamical_time;
- float dt = dt_dyn;
-
- /* absolute value of height above disk */
- const float dz = fabs(g->x[2] - potential->disk_patch_potential.z_disk);
-
- /* vertical cceleration */
- const float z_accel = 2 * M_PI * phys_const->const_newton_G *
- potential->disk_patch_potential.surface_density *
- tanh(dz / potential->disk_patch_potential.scale_height);
-
- /* demand that dt * velocity < fraction of scale height of disk */
- float dt1 = FLT_MAX;
- if (fabs(g->v_full[2]) > 0) {
- dt1 = potential->disk_patch_potential.scale_height / fabs(g->v_full[2]);
- if (dt1 < dt) dt = dt1;
- }
-
- /* demand that dt^2 * acceleration < fraction of scale height of disk */
- float dt2 = FLT_MAX;
- if (fabs(z_accel) > 0) {
- dt2 = potential->disk_patch_potential.scale_height / fabs(z_accel);
- if (dt2 < dt * dt) dt = sqrt(dt2);
- }
-
- /* demand that dt^3 jerk < fraction of scale height of disk */
- float dt3 = FLT_MAX;
- if (abs(g->v_full[2]) > 0) {
- const float dz_accel_over_dt =
- 2 * M_PI * phys_const->const_newton_G *
- potential->disk_patch_potential.surface_density /
- potential->disk_patch_potential.scale_height /
- cosh(dz / potential->disk_patch_potential.scale_height) /
- cosh(dz / potential->disk_patch_potential.scale_height) *
- fabs(g->v_full[2]);
-
- dt3 = potential->disk_patch_potential.scale_height / fabs(dz_accel_over_dt);
- if (dt3 < dt * dt * dt) dt = pow(dt3, 1. / 3.);
- }
-
- return potential->disk_patch_potential.timestep_mult * dt;
-}
-
-/**
- * @brief Computes the gravitational acceleration along z due to a hydrostatic
- * disk
- *
- * See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
- *
- * @param time The current time in internal units.
- * @param potential The properties of the potential.
- * @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_disk_patch_potential(
- double time, const struct external_potential* restrict potential,
- const struct phys_const* restrict phys_const, struct gpart* restrict g) {
-
- const float G_newton = phys_const->const_newton_G;
- const float dz = g->x[2] - potential->disk_patch_potential.z_disk;
- const float t_dyn = potential->disk_patch_potential.dynamical_time;
-
- float reduction_factor = 1.;
- if (time < potential->disk_patch_potential.growth_time * t_dyn)
- reduction_factor =
- time / (potential->disk_patch_potential.growth_time * t_dyn);
-
- const float z_accel =
- reduction_factor * 2 * G_newton * M_PI *
- potential->disk_patch_potential.surface_density *
- tanh(fabs(dz) / potential->disk_patch_potential.scale_height);
-
- /* Accelerations. Note that they are multiplied by G later on */
- if (dz > 0) g->a_grav[2] -= z_accel / G_newton;
- if (dz < 0) g->a_grav[2] += z_accel / G_newton;
-}
-#endif /* EXTERNAL_POTENTIAL_DISK_PATCH */
-
-/* ------------------------------------------------------------------------- */
-
-#ifdef EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
-
-/**
- * @brief Computes the time-step due to the acceleration from an isothermal
- * potential.
- *
- * @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_isothermalpotential_timestep(
- 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->isothermal_potential.x;
- const float dy = g->x[1] - potential->isothermal_potential.y;
- const float dz = g->x[2] - potential->isothermal_potential.z;
-
- const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
- const float drdv =
- dx * (g->v_full[0]) + dy * (g->v_full[1]) + dz * (g->v_full[2]);
- const double vrot = potential->isothermal_potential.vrot;
-
- const float dota_x =
- vrot * vrot * rinv2 * (g->v_full[0] - 2 * drdv * dx * rinv2);
- const float dota_y =
- vrot * vrot * rinv2 * (g->v_full[1] - 2 * drdv * dy * rinv2);
- const float dota_z =
- vrot * vrot * rinv2 * (g->v_full[2] - 2 * drdv * dz * rinv2);
- const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
- const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] +
- g->a_grav[2] * g->a_grav[2];
-
- return potential->isothermal_potential.timestep_mult * sqrtf(a_2 / dota_2);
-}
-
-/**
- * @brief Computes the gravitational acceleration from an isothermal potential.
- *
- * Note that the accelerations are multiplied by Newton's G constant
- * later on.
- *
- * @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_isothermalpotential(const struct external_potential* potential,
- const struct phys_const* const phys_const,
- struct gpart* g) {
-
- const float G_newton = phys_const->const_newton_G;
- const float dx = g->x[0] - potential->isothermal_potential.x;
- const float dy = g->x[1] - potential->isothermal_potential.y;
- const float dz = g->x[2] - potential->isothermal_potential.z;
- const float rinv2 = 1.f / (dx * dx + dy * dy + dz * dz);
-
- const double vrot = potential->isothermal_potential.vrot;
- const double term = -vrot * vrot * rinv2 / G_newton;
-
- g->a_grav[0] += term * dx;
- g->a_grav[1] += term * dy;
- g->a_grav[2] += term * dz;
-}
-
-#endif /* EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL */
-
-/* ------------------------------------------------------------------------- */
-
-/* Include exteral pointmass potential */
-#ifdef EXTERNAL_POTENTIAL_POINTMASS
-
-/**
- * @brief Computes the time-step due to the acceleration from a point mass
- *
- * @param potential The properties of the externa potential.
- * @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_pointmass_timestep(
- 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;
- const float dx = g->x[0] - potential->point_mass.x;
- const float dy = g->x[1] - potential->point_mass.y;
- const float dz = g->x[2] - potential->point_mass.z;
- const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
- const float drdv = (g->x[0] - potential->point_mass.x) * (g->v_full[0]) +
- (g->x[1] - potential->point_mass.y) * (g->v_full[1]) +
- (g->x[2] - potential->point_mass.z) * (g->v_full[2]);
- const float dota_x = G_newton * potential->point_mass.mass * rinv * rinv *
- rinv * (-g->v_full[0] + 3.f * rinv * rinv * drdv * dx);
- const float dota_y = G_newton * potential->point_mass.mass * rinv * rinv *
- rinv * (-g->v_full[1] + 3.f * rinv * rinv * drdv * dy);
- const float dota_z = G_newton * potential->point_mass.mass * rinv * rinv *
- rinv * (-g->v_full[2] + 3.f * rinv * rinv * drdv * dz);
- const float dota_2 = dota_x * dota_x + dota_y * dota_y + dota_z * dota_z;
- const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] +
- g->a_grav[2] * g->a_grav[2];
-
- return potential->point_mass.timestep_mult * sqrtf(a_2 / dota_2);
-}
-
-/**
- * @brief Computes the gravitational acceleration of a particle due to a
- * point mass
- *
- * Note that the accelerations are multiplied by Newton's G constant later
- * on.
- *
- * @param potential The proerties of the external potential.
- * @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_pointmass(
- const struct external_potential* restrict potential,
- const struct phys_const* restrict phys_const, struct gpart* restrict g) {
-
- const float dx = g->x[0] - potential->point_mass.x;
- const float dy = g->x[1] - potential->point_mass.y;
- const float dz = g->x[2] - potential->point_mass.z;
- const float rinv = 1.f / sqrtf(dx * dx + dy * dy + dz * dz);
- const float rinv3 = rinv * rinv * rinv;
-
- g->a_grav[0] += -potential->point_mass.mass * dx * rinv3;
- g->a_grav[1] += -potential->point_mass.mass * dy * rinv3;
- g->a_grav[2] += -potential->point_mass.mass * dz * rinv3;
-}
-
-#endif /* EXTERNAL_POTENTIAL_POINTMASS */
-
-/* ------------------------------------------------------------------------- */
-
-/* Now, some generic functions, defined in the source file */
-void potential_init(const struct swift_params* parameter_file,
- const struct phys_const* phys_const,
- const struct UnitSystem* us,
- struct external_potential* potential);
-
-void potential_print(const struct external_potential* potential);
-
-#endif /* SWIFT_POTENTIALS_H */
diff --git a/src/runner.c b/src/runner.c
index 76dc9fd84d..0a08bb5abf 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -152,7 +152,7 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
/* Is this part within the time step? */
if (g->ti_end <= ti_current) {
- external_gravity(time, potential, constants, g);
+ external_gravity_acceleration(time, potential, constants, g);
}
}
diff --git a/src/swift.h b/src/swift.h
index 80a4b1a1f7..3119adfd4b 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -43,7 +43,7 @@
#include "part.h"
#include "partition.h"
#include "physical_constants.h"
-#include "potentials.h"
+#include "potential.h"
#include "queue.h"
#include "runner.h"
#include "scheduler.h"
diff --git a/src/timestep.h b/src/timestep.h
index ca754a439a..599fb4762e 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -68,8 +68,9 @@ __attribute__((always_inline)) INLINE static int get_integer_timestep(
__attribute__((always_inline)) INLINE static int get_gpart_timestep(
const struct gpart *restrict gp, const struct engine *restrict e) {
- const float new_dt_external = gravity_compute_timestep_external(
- e->external_potential, e->physical_constants, gp);
+ const float new_dt_external = external_gravity_timestep(
+ e->time, e->external_potential, e->physical_constants, gp);
+
/* const float new_dt_self = */
/* gravity_compute_timestep_self(e->physical_constants, gp); */
const float new_dt_self = FLT_MAX; // MATTHIEU
@@ -111,8 +112,9 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
float new_dt_grav = FLT_MAX;
if (p->gpart != NULL) {
- const float new_dt_external = gravity_compute_timestep_external(
- e->external_potential, e->physical_constants, p->gpart);
+ const float new_dt_external = external_gravity_timestep(
+ e->time, e->external_potential, e->physical_constants, p->gpart);
+
/* const float new_dt_self = */
/* gravity_compute_timestep_self(e->physical_constants, p->gpart); */
const float new_dt_self = FLT_MAX; // MATTHIEU
--
GitLab