diff --git a/AUTHORS b/AUTHORS
index c822300c22885a05b42d58a51cc86af9da410429..6f283405b69a7d3a5397916f0a3afa7f4fb54a4a 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -10,3 +10,4 @@ Tom Theuns tom.theuns@durham.ac.uk
Richard G. Bower r.g.bower@durham.ac.uk
Stefan Arridge stefan.arridge@durham.ac.uk
Massimiliano Culpo massimiliano.culpo@googlemail.com
+Yves Revaz yves.revaz@epfl.ch
diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in
index 0658dd27e043b11c7480eaf8b01c0c40809e215c..2a5aeba7d1db0b1e1e56a9a6eed3059aba6a09ff 100644
--- a/doc/Doxyfile.in
+++ b/doc/Doxyfile.in
@@ -763,7 +763,8 @@ 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/cooling/const_lambda
+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
# that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
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 78b42e78356f83e80eee8e7f5f91ad7dcf90c37f..c76e4f612250d180f2ba2fccd0c6209878173433 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 702a50ff53b73d004ff36be8049823515675cccf..42cd26e235deb17a899a65050ef5caa9c832c59c 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 0000000000000000000000000000000000000000..9af1011ee653253f0d1b2cd26db0ac13cf11adc0
--- /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 4bd0d001975d80b6729cf2ef7b95f81da5bc4fe8..04e0afdf7e6d2b4f0122a3d7d1bd1084539c405e 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 56578731eb16a27febb3627524956b4e38b6428f..42853e6b51983f2868528202adec3fc829c2ddbc 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 ebf04675852a7663119ed1ecfd33a05da6b7bb15..6a27016b8a3f484b7c1c9b74594073d5f28efe90 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 55df81a08d16c6a4f39aa5e9e9205101dedaa3a9..8bd67c5b08de82bb6a3d47ccf3419f85e3e5c6b1 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 40b1d1f1ff9e08dae0c4b0b1539ca773c93009b4..48cc578658a9520477d40bf504e3eb7c3c8e5164 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 31e027e3a308b04f1c59222e1a339786857061ac..950aebc65d7d34cd7aaeb2368734e5492902a912 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 a123ad24d7ca34105c22f5f31e75c688c681288f..0000000000000000000000000000000000000000
--- 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 d06c165651ce8f33692d1512ddd8fdae80ffb556..ce300b32157361e7860d201c186823471a179c0a 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/MultiTypes/multiTypes.yml b/examples/MultiTypes/multiTypes.yml
index a9e4e46744a703d45ce6344cdadc116a1d0f306d..51a6d2b478681e2e1c61e199f758e35c507ec195 100644
--- a/examples/MultiTypes/multiTypes.yml
+++ b/examples/MultiTypes/multiTypes.yml
@@ -35,7 +35,7 @@ InitialConditions:
file_name: ./multiTypes.hdf5 # The file to read
# 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/UniformBox_3D/uniformBox.yml b/examples/UniformBox_3D/uniformBox.yml
index 7c9c74e1342bffb939131a265188cae269cc773f..8aaa802b64de46244f7066bce00f342cad8c5ef0 100644
--- a/examples/UniformBox_3D/uniformBox.yml
+++ b/examples/UniformBox_3D/uniformBox.yml
@@ -33,11 +33,3 @@ SPH:
# Parameters related to the initial conditions
InitialConditions:
file_name: ./uniformBox.hdf5 # The file to read
-
-
- # External potential parameters
-PointMass:
- position_x: 50. # location of external point mass in internal units
- position_y: 50.
- position_z: 50.
- mass: 1e10 # mass of external point mass in internal units
diff --git a/examples/main.c b/examples/main.c
index 084a34723a7b1c3d6e075a67936c14b2380e69ad..895bb9003e04a4723b13ea1eaa360c5aa44e7228 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -439,9 +439,9 @@ int main(int argc, char *argv[]) {
if (with_external_gravity && myrank == 0) potential_print(&potential);
/* Initialise the cooling function properties */
- struct cooling_data cooling;
- if (with_cooling) cooling_init(params, &us, &prog_const, &cooling);
- if (with_cooling && myrank == 0) cooling_print(&cooling);
+ struct cooling_function_data cooling_func;
+ if (with_cooling) cooling_init(params, &us, &prog_const, &cooling_func);
+ if (with_cooling && myrank == 0) cooling_print(&cooling_func);
/* Construct the engine policy */
int engine_policies = ENGINE_POLICY | engine_policy_steal;
@@ -457,7 +457,7 @@ int main(int argc, char *argv[]) {
struct engine e;
engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, with_aff,
engine_policies, talking, &us, &prog_const, &hydro_properties,
- &potential, &cooling);
+ &potential, &cooling_func);
if (myrank == 0) {
clocks_gettime(&toc);
message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 10bd09a18f7e8099e4db559e76957d19d8f90164..be2f512613edffd72eaf03689bccee0dd755726b 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. # Position of the disc along the z-axis (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 2343ab99ffd90a27e588344c1fae4f1491b4625e..f7cb52ba40a34269173ab8c5019d3011b0c34b61 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -42,8 +42,8 @@ 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 \
- hydro_properties.h threadpool.h cooling.h
+ physical_constants.h physical_constants_cgs.h potential.h version.h \
+ hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h
# Common source files
@@ -51,11 +51,11 @@ 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.
-nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
+nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
kernel_long_gravity.h vector.h runner_doiact.h runner_doiact_grav.h runner_doiact_fft.h \
units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
dimension.h equation_of_state.h \
@@ -71,9 +71,13 @@ nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_h
hydro/Gadget2/hydro_debug.h hydro/Gadget2/hydro_part.h \
hydro/Gizmo/hydro.h hydro/Gizmo/hydro_iact.h hydro/Gizmo/hydro_io.h \
hydro/Gizmo/hydro_debug.h hydro/Gizmo/hydro_part.h \
- riemann.h riemann/riemann_hllc.h riemann/riemann_trrs.h \
+ riemann/riemann_hllc.h riemann/riemann_trrs.h \
riemann/riemann_exact.h riemann/riemann_vacuum.h \
- cooling/const_du/cooling.h cooling/const_lambda/cooling.h
+ potential/none/potential.h potential/point_mass/potential.h \
+ potential/isothermal/potential.h potential/disc_patch/potential.h \
+ cooling/none/cooling.h cooling/none/cooling_struct.h \
+ cooling/const_du/cooling.h cooling/const_du/cooling_struct.h \
+ cooling/const_lambda/cooling.h cooling/const_lambda/cooling_struct.h
# Sources and flags for regular library
diff --git a/src/align.h b/src/align.h
new file mode 100644
index 0000000000000000000000000000000000000000..84e2909c0866c18f0f8378df9d0efc8d0f6545b5
--- /dev/null
+++ b/src/align.h
@@ -0,0 +1,27 @@
+/*******************************************************************************
+ * 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_ALIGN_H
+#define SWIFT_ALIGN_H
+
+/**
+ * @brief Defines alignment of structures
+ */
+#define SWIFT_STRUCT_ALIGN __attribute__((aligned(32)))
+
+#endif /* SWIFT_ALIGN_H */
diff --git a/src/cell.h b/src/cell.h
index fd836206241c55cd6b9da2e29157c11a14c5a892..b78cc0a8f842770f60777e3986616a175d2f33ca 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -30,6 +30,7 @@
#include <stddef.h>
/* Local includes. */
+#include "align.h"
#include "lock.h"
#include "multipole.h"
#include "part.h"
@@ -44,7 +45,7 @@ struct space;
* The maximum was lowered by a further factor of 2 to be on the safe side.*/
#define cell_max_tag (1 << 29)
-#define cell_align 32
+#define cell_align 128
/* Global variables. */
extern int cell_next_tag;
@@ -74,7 +75,8 @@ struct pcell {
/* Relative indices of the cell's progeny. */
int progeny[8];
-};
+
+} SWIFT_STRUCT_ALIGN;
/* Structure to store the data of a single cell. */
struct cell {
@@ -182,7 +184,7 @@ struct cell {
double mom[3], ang_mom[3];
/* Mass, potential, internal and kinetic energy of particles in this cell. */
- double mass, e_pot, e_int, e_kin, entropy;
+ double mass, e_pot, e_int, e_kin, e_rad, entropy;
/* Number of particles updated in this cell. */
int updated, g_updated;
@@ -208,7 +210,7 @@ struct cell {
#endif
-} __attribute__((aligned(cell_align)));
+} SWIFT_STRUCT_ALIGN;
/* Convert cell location to ID. */
#define cell_getid(cdim, i, j, k) \
diff --git a/src/const.h b/src/const.h
index 552d49c8f1e4dd9f4fa2855e5fed548acd5c0b3f..316a13c8ebba4785e8e1a66b8307f8681697c938 100644
--- a/src/const.h
+++ b/src/const.h
@@ -90,13 +90,15 @@
#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
//#define COOLING_CONST_DU
-#define COOLING_CONST_LAMBDA
+//#define COOLING_CONST_LAMBDA
//#define COOLING_GRACKLE
/* Are we debugging ? */
diff --git a/src/cooling.c b/src/cooling.c
index 102adba9d521dbe53caeff3c4409292bb64a29b2..e0208dbb591445d0877ef1e703d6e8cf349ddfd6 100644
--- a/src/cooling.c
+++ b/src/cooling.c
@@ -36,7 +36,7 @@
void cooling_init(const struct swift_params* parameter_file,
const struct UnitSystem* us,
const struct phys_const* phys_const,
- struct cooling_data* cooling) {
+ struct cooling_function_data* cooling) {
cooling_init_backend(parameter_file, us, phys_const, cooling);
}
@@ -48,7 +48,7 @@ void cooling_init(const struct swift_params* parameter_file,
*
* @param cooling The properties of the cooling function.
*/
-void cooling_print(const struct cooling_data* cooling) {
+void cooling_print(const struct cooling_function_data* cooling) {
cooling_print_backend(cooling);
}
diff --git a/src/cooling.h b/src/cooling.h
index 034ee2329d91b17b875932bc4eb40d0d80a05111..1b326f6dc4fdf796dd1587e73e9b591f0f500ccb 100644
--- a/src/cooling.h
+++ b/src/cooling.h
@@ -31,7 +31,9 @@
#include "const.h"
/* Import the right cooling definition */
-#if defined(COOLING_CONST_DU)
+#if defined(COOLING_NONE)
+#include "./cooling/none/cooling.h"
+#elif defined(COOLING_CONST_DU)
#include "./cooling/const_du/cooling.h"
#elif defined(COOLING_CONST_LAMBDA)
#include "./cooling/const_lambda/cooling.h"
@@ -45,8 +47,8 @@
void cooling_init(const struct swift_params* parameter_file,
const struct UnitSystem* us,
const struct phys_const* phys_const,
- struct cooling_data* cooling);
+ struct cooling_function_data* cooling);
-void cooling_print(const struct cooling_data* cooling);
+void cooling_print(const struct cooling_function_data* cooling);
#endif /* SWIFT_COOLING_H */
diff --git a/src/cooling/const_du/cooling.h b/src/cooling/const_du/cooling.h
index 634723f5a5a9af56aefc5fee1d409f583b243eba..b25980ff2269ca9ea176edcc2a3c771647819133 100644
--- a/src/cooling/const_du/cooling.h
+++ b/src/cooling/const_du/cooling.h
@@ -22,11 +22,12 @@
#define SWIFT_COOLING_CONST_DU_H
/**
- * @file src/cooling/const/cooling.h
+ * @file src/cooling/const_du/cooling.h
* @brief Routines related to the "constant cooling" cooling function.
*
* This is the simplest possible cooling function. A constant cooling rate with
- * a minimal energy floor is applied.
+ * a minimal energy floor is applied. Should be used as a template for more
+ * realistic functions.
*/
/* Some standard headers. */
@@ -34,6 +35,7 @@
/* Local includes. */
#include "const.h"
+#include "cooling_struct.h"
#include "error.h"
#include "hydro.h"
#include "parser.h"
@@ -41,33 +43,24 @@
#include "physical_constants.h"
#include "units.h"
-/**
- * @brief Properties of the cooling function.
- */
-struct cooling_data {
-
- /*! Cooling rate in internal units. du_dt = -cooling_rate */
- float cooling_rate;
-
- /*! Minimally allowed internal energy of the particles */
- float min_energy;
-
- /*! Constant multiplication factor for time-step criterion */
- float cooling_tstep_mult;
-};
-
/**
* @brief Apply the cooling function to a particle.
*
+ * In this simple example we just apply the const cooling rate
+ * and check that we don't go below the given floor.
+ *
* @param phys_const The physical constants in internal units.
* @param us The internal system of units.
- * @param cooling The #cooling_data used in the run.
+ * @param cooling The #cooling_function_data used in the run.
* @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
* @param dt The time-step of this particle.
*/
__attribute__((always_inline)) INLINE static void cooling_cool_part(
- const struct phys_const* const phys_const, const struct UnitSystem* us,
- const struct cooling_data* cooling, struct part* p, float dt) {
+ const struct phys_const* restrict phys_const,
+ const struct UnitSystem* restrict us,
+ const struct cooling_function_data* restrict cooling,
+ struct part* restrict p, struct xpart* restrict xp, float dt) {
/* Get current internal energy (dt=0) */
const float u_old = hydro_get_internal_energy(p, 0.f);
@@ -86,27 +79,69 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
/* Update the internal energy */
hydro_set_internal_energy(p, u_new);
+
+ /* Store the radiated energy */
+ xp->cooling_data.radiated_energy += hydro_get_mass(p) * (u_old - u_new);
}
/**
* @brief Computes the cooling time-step.
*
- * @param cooling The #cooling_data used in the run.
+ * In this simple example, we return \f$ \alpha \frac{u}{du/dt} \f$.
+ * This is used to compute the time-step of the particle. Cooling functions
+ * that are solved implicitly can simply return FLT_MAX here.
+ *
+ * @param cooling The #cooling_function_data used in the run.
* @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
* @param p Pointer to the particle data.
*/
-__attribute__((always_inline)) INLINE static double cooling_timestep(
- const struct cooling_data* cooling,
- const struct phys_const* const phys_const, const struct part* const p) {
+__attribute__((always_inline)) INLINE static float cooling_timestep(
+ const struct cooling_function_data* restrict cooling,
+ const struct phys_const* restrict phys_const,
+ const struct UnitSystem* restrict us, const struct part* restrict p) {
const float cooling_rate = cooling->cooling_rate;
- const float internal_energy =
- hydro_get_internal_energy(p, 0); // dt = 0 because using current entropy
+ const float internal_energy = hydro_get_internal_energy(p, 0);
return cooling->cooling_tstep_mult * internal_energy / cooling_rate;
}
/**
- * @brief Initialises the cooling properties.
+ * @brief Sets the cooling properties of the (x-)particles to a valid start
+ * state.
+ *
+ * In this case, we set the total radiated energy to 0. Note that the particle
+ * structure is just passed in for cases where information needs to be read
+ * from there.
+ *
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static void cooling_init_part(
+ const struct part* restrict p, struct xpart* restrict xp) {
+
+ xp->cooling_data.radiated_energy = 0.f;
+}
+
+/**
+ * @brief Returns the total radiated energy by this particle.
+ *
+ * In this simple example we jsut return the quantity accumulated in the
+ * #cooling_xpart_data of this particle.
+ *
+ * @param xp The extended particle data
+ */
+__attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
+ const struct xpart* restrict xp) {
+
+ return xp->cooling_data.radiated_energy;
+}
+
+/**
+ * @brief Initialises the cooling function properties from the parameter file
+ *
+ * In this example, we just read in the values from the YAML file without
+ * doing any conversions or multiplying any constants in.
*
* @param parameter_file The parsed parameter file.
* @param us The current internal system of units.
@@ -115,7 +150,8 @@ __attribute__((always_inline)) INLINE static double cooling_timestep(
*/
static INLINE void cooling_init_backend(
const struct swift_params* parameter_file, const struct UnitSystem* us,
- const struct phys_const* phys_const, struct cooling_data* cooling) {
+ const struct phys_const* phys_const,
+ struct cooling_function_data* cooling) {
cooling->cooling_rate =
parser_get_param_double(parameter_file, "ConstCooling:cooling_rate");
@@ -130,9 +166,10 @@ static INLINE void cooling_init_backend(
*
* @param cooling The properties of the cooling function.
*/
-static INLINE void cooling_print_backend(const struct cooling_data* cooling) {
+static INLINE void cooling_print_backend(
+ const struct cooling_function_data* cooling) {
- message("Cooling function is 'Constant cooling' with rate %f and floor %f",
+ message("Cooling function is 'Constant cooling' with rate %f and floor %f.",
cooling->cooling_rate, cooling->min_energy);
}
diff --git a/src/cooling/const_du/cooling_struct.h b/src/cooling/const_du/cooling_struct.h
new file mode 100644
index 0000000000000000000000000000000000000000..cc00b001cf6b576266de02dac885f87d089bd8e4
--- /dev/null
+++ b/src/cooling/const_du/cooling_struct.h
@@ -0,0 +1,60 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Richard Bower (r.g.bower@durham.ac.uk)
+ * Stefan Arridge (stefan.arridge@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_COOLING_STRUCT_CONST_DU_H
+#define SWIFT_COOLING_STRUCT_CONST_DU_H
+
+/**
+ * @file src/cooling/const_du/cooling_struct.h
+ * @brief Structure related to the "constant cooling" cooling function.
+ *
+ * This is the simplest possible cooling function. A constant cooling rate with
+ * a minimal energy floor is applied. Should be used as a template for more
+ * realistic functions.
+ */
+
+/**
+ * @brief Properties of the cooling function.
+ */
+struct cooling_function_data {
+
+ /*! Cooling rate in internal units. du_dt = -cooling_rate */
+ float cooling_rate;
+
+ /*! Minimally allowed internal energy of the particles */
+ float min_energy;
+
+ /*! Constant multiplication factor for time-step criterion */
+ float cooling_tstep_mult;
+};
+
+/**
+ * @brief Properties of the cooling stored in the particle data.
+ *
+ * This is used to carry properties such as the total amount of
+ * energy radiated away.
+ */
+struct cooling_xpart_data {
+
+ /*! Energy radiated away by this particle since the start of the run */
+ float radiated_energy;
+};
+
+#endif /* SWIFT_COOLING_STRUCT_CONST_DU_H */
diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h
index 3fd583b55a31534996e23e4dce71bcae9394d7de..11cf2cae51f1ab2646d1391d2164c399c77a7bba 100644
--- a/src/cooling/const_lambda/cooling.h
+++ b/src/cooling/const_lambda/cooling.h
@@ -35,40 +35,17 @@
#include "physical_constants.h"
#include "units.h"
-/* Cooling Properties */
-struct cooling_data {
-
- /*! Cooling rate in cgs units. Defined by 'rho * du/dt = -lambda * n_H^2'*/
- float lambda;
-
- /*! Minimum temperature (in Kelvin) for all gas particles*/
- float min_temperature;
-
- /*! Fraction of gas mass that is Hydrogen. Used to calculate n_H*/
- float hydrogen_mass_abundance;
-
- /* 'mu', used to convert min_temperature to min_internal energy*/
- float mean_molecular_weight;
-
- /*! Minimally allowed internal energy of the particles */
- float min_energy;
- float min_energy_cgs;
-
- /*! Constant multiplication factor for time-step criterion */
- float cooling_tstep_mult;
-};
-
/**
* @brief Calculates du/dt in code units for a particle.
*
* @param phys_const The physical constants in internal units.
* @param us The internal system of units.
- * @param cooling The #cooling_data used in the run.
+ * @param cooling The #cooling_function_data used in the run.
* @param p Pointer to the particle data..
*/
__attribute__((always_inline)) INLINE static float cooling_rate(
const struct phys_const* const phys_const, const struct UnitSystem* us,
- const struct cooling_data* cooling, const struct part* p) {
+ const struct cooling_function_data* cooling, const struct part* p) {
/* Get particle properties */
/* Density */
@@ -101,13 +78,15 @@ __attribute__((always_inline)) INLINE static float cooling_rate(
*
* @param phys_const The physical constants in internal units.
* @param us The internal system of units.
- * @param cooling The #cooling_data used in the run.
+ * @param cooling The #cooling_function_data used in the run.
* @param p Pointer to the particle data.
* @param dt The time-step of this particle.
*/
__attribute__((always_inline)) INLINE static void cooling_cool_part(
- const struct phys_const* const phys_const, const struct UnitSystem* us,
- const struct cooling_data* cooling, struct part* p, float dt) {
+ const struct phys_const* restrict phys_const,
+ const struct UnitSystem* restrict us,
+ const struct cooling_function_data* restrict cooling,
+ struct part* restrict p, struct xpart* restrict xp, float dt) {
/* Get current internal energy (dt=0) */
const float u_old = hydro_get_internal_energy(p, 0.f);
@@ -135,20 +114,23 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
* " */
/* "%g, du_dt*dt = %g, u_old + du_dt*dt = %g, u_new = %g\n", */
/* u_old, du_dt, dt, du_dt * dt, u_new, u_new_test); */
+
+ /* Store the radiated energy */
+ xp->cooling_data.radiated_energy += hydro_get_mass(p) * (u_old - u_new);
}
/**
* @brief Computes the time-step due to cooling
*
- * @param cooling The #cooling_data used in the run.
+ * @param cooling The #cooling_function_data used in the run.
* @param phys_const The physical constants in internal units.
* @param us The internal system of units.
* @param p Pointer to the particle data.
*/
__attribute__((always_inline)) INLINE static float cooling_timestep(
- const struct cooling_data* cooling,
- const struct phys_const* const phys_const, const struct UnitSystem* us,
- const struct part* const p) {
+ const struct cooling_function_data* restrict cooling,
+ const struct phys_const* restrict phys_const,
+ const struct UnitSystem* restrict us, const struct part* restrict p) {
/* Get du_dt */
const float du_dt = cooling_rate(phys_const, us, cooling, p);
@@ -159,6 +141,30 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
return u / du_dt;
}
+/**
+ * @brief Sets the cooling properties of the (x-)particles to a valid start
+ * state.
+ *
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static void cooling_init_part(
+ const struct part* restrict p, struct xpart* restrict xp) {
+
+ xp->cooling_data.radiated_energy = 0.f;
+}
+
+/**
+ * @brief Returns the total radiated energy by this particle.
+ *
+ * @param xp The extended particle data
+ */
+__attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
+ const struct xpart* restrict xp) {
+
+ return xp->cooling_data.radiated_energy;
+}
+
/**
* @brief Initialises the cooling properties.
*
@@ -169,7 +175,8 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
*/
static INLINE void cooling_init_backend(
const struct swift_params* parameter_file, const struct UnitSystem* us,
- const struct phys_const* phys_const, struct cooling_data* cooling) {
+ const struct phys_const* phys_const,
+ struct cooling_function_data* cooling) {
cooling->lambda =
parser_get_param_double(parameter_file, "LambdaCooling:lambda");
@@ -199,7 +206,8 @@ static INLINE void cooling_init_backend(
*
* @param cooling The properties of the cooling function.
*/
-static INLINE void cooling_print_backend(const struct cooling_data* cooling) {
+static INLINE void cooling_print_backend(
+ const struct cooling_function_data* cooling) {
message(
"Cooling function is 'Constant lambda' with "
diff --git a/src/cooling/const_lambda/cooling_struct.h b/src/cooling/const_lambda/cooling_struct.h
new file mode 100644
index 0000000000000000000000000000000000000000..27c5df16bffbe7d165237d201ca68ea4ba89dd73
--- /dev/null
+++ b/src/cooling/const_lambda/cooling_struct.h
@@ -0,0 +1,60 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Richard Bower (r.g.bower@durham.ac.uk)
+ * Stefan Arridge (stefan.arridge@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_COOLING_STRUCT_CONST_LAMBDA_H
+#define SWIFT_COOLING_STRUCT_CONST_LAMBDA_H
+
+/**
+ * @brief Properties of the cooling function.
+ */
+struct cooling_function_data {
+
+ /*! Cooling rate in cgs units. Defined by 'rho * du/dt = -lambda * n_H^2'*/
+ float lambda;
+
+ /*! Minimum temperature (in Kelvin) for all gas particles*/
+ float min_temperature;
+
+ /*! Fraction of gas mass that is Hydrogen. Used to calculate n_H*/
+ float hydrogen_mass_abundance;
+
+ /* 'mu', used to convert min_temperature to min_internal energy*/
+ float mean_molecular_weight;
+
+ /*! Minimally allowed internal energy of the particles */
+ float min_energy;
+ float min_energy_cgs;
+
+ /*! Constant multiplication factor for time-step criterion */
+ float cooling_tstep_mult;
+};
+
+/**
+ * @brief Properties of the cooling stored in the particle data.
+ */
+struct cooling_xpart_data {
+
+ /*! Energy radiated away by this particle since the start of the run */
+ float radiated_energy;
+};
+
+#endif /* SWIFT_COOLING_STRUCT_CONST_LAMBDA_H */
diff --git a/src/cooling/none/cooling.h b/src/cooling/none/cooling.h
new file mode 100644
index 0000000000000000000000000000000000000000..0461100dc11e7ffbb4616766923142442b4ac943
--- /dev/null
+++ b/src/cooling/none/cooling.h
@@ -0,0 +1,125 @@
+/*******************************************************************************
+ * 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_COOLING_NONE_H
+#define SWIFT_COOLING_NONE_H
+
+/**
+ * @file src/cooling/none/cooling.h
+ * @brief Empty infrastructure for the cases without cooling function
+ */
+
+/* Some standard headers. */
+#include <float.h>
+#include <math.h>
+
+/* Local includes. */
+#include "error.h"
+#include "hydro.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "units.h"
+
+/**
+ * @brief Apply the cooling function to a particle.
+ *
+ * We do nothing.
+ *
+ * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
+ * @param cooling The #cooling_function_data used in the run.
+ * @param p Pointer to the particle data.
+ * @param dt The time-step of this particle.
+ */
+__attribute__((always_inline)) INLINE static void cooling_cool_part(
+ const struct phys_const* restrict phys_const,
+ const struct UnitSystem* restrict us,
+ const struct cooling_function_data* restrict cooling,
+ struct part* restrict p, struct xpart* restrict xp, float dt) {}
+
+/**
+ * @brief Computes the cooling time-step.
+ *
+ * We return FLT_MAX so as to impose no limit on the time-step.
+ *
+ * @param cooling The #cooling_function_data used in the run.
+ * @param phys_const The physical constants in internal units.
+ * @param us The internal system of units.
+ * @param p Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float cooling_timestep(
+ const struct cooling_function_data* restrict cooling,
+ const struct phys_const* restrict phys_const,
+ const struct UnitSystem* restrict us, const struct part* restrict p) {
+
+ return FLT_MAX;
+}
+
+/**
+ * @brief Sets the cooling properties of the (x-)particles to a valid start
+ * state.
+ *
+ * Nothing to do here.
+ *
+ * @param p Pointer to the particle data.
+ * @param xp Pointer to the extended particle data.
+ */
+__attribute__((always_inline)) INLINE static void cooling_init_part(
+ const struct part* restrict p, struct xpart* restrict xp) {}
+
+/**
+ * @brief Returns the total radiated energy by this particle.
+ *
+ * No cooling, so return 0.
+ *
+ * @param xp The extended particle data
+ */
+__attribute__((always_inline)) INLINE static float cooling_get_radiated_energy(
+ const struct xpart* restrict xp) {
+
+ return 0.f;
+}
+
+/**
+ * @brief Initialises the cooling properties.
+ *
+ * Nothing to do here.
+ *
+ * @param parameter_file The parsed parameter file.
+ * @param us The current internal system of units.
+ * @param phys_const The physical constants in internal units.
+ * @param cooling The cooling properties to initialize
+ */
+static INLINE void cooling_init_backend(
+ const struct swift_params* parameter_file, const struct UnitSystem* us,
+ const struct phys_const* phys_const,
+ struct cooling_function_data* cooling) {}
+
+/**
+ * @brief Prints the properties of the cooling model to stdout.
+ *
+ * @param cooling The properties of the cooling function.
+ */
+static INLINE void cooling_print_backend(
+ const struct cooling_function_data* cooling) {
+
+ message("Cooling function is 'No cooling'.");
+}
+
+#endif /* SWIFT_COOLING_NONE_H */
diff --git a/src/cooling/none/cooling_struct.h b/src/cooling/none/cooling_struct.h
new file mode 100644
index 0000000000000000000000000000000000000000..a08530c44d7405df934136f2861f84ba619d2595
--- /dev/null
+++ b/src/cooling/none/cooling_struct.h
@@ -0,0 +1,37 @@
+/*******************************************************************************
+ * 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_COOLING_STRUCT_NONE_H
+#define SWIFT_COOLING_STRUCT_NONE_H
+
+/**
+ * @file src/cooling/none/cooling_struct.h
+ * @brief Empty infrastructure for the cases without cooling function
+ */
+
+/**
+ * @brief Properties of the cooling function.
+ */
+struct cooling_function_data {};
+
+/**
+ * @brief Properties of the cooling stored in the particle data
+ */
+struct cooling_xpart_data {};
+
+#endif /* SWIFT_COOLING_STRUCT_NONE_H */
diff --git a/src/cooling_struct.h b/src/cooling_struct.h
new file mode 100644
index 0000000000000000000000000000000000000000..0c567788423ae39507864de8b4a687eeed358cb6
--- /dev/null
+++ b/src/cooling_struct.h
@@ -0,0 +1,46 @@
+/*******************************************************************************
+ * 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_COOLING_STRUCT_H
+#define SWIFT_COOLING_STRUCT_H
+
+/**
+ * @file src/cooling_struct.h
+ * @brief Branches between the different cooling functions.
+ */
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local headers. */
+#include "const.h"
+
+/* Import the right cooling definition */
+#if defined(COOLING_NONE)
+#include "./cooling/none/cooling_struct.h"
+#elif defined(COOLING_CONST_DU)
+#include "./cooling/const_du/cooling_struct.h"
+#elif defined(COOLING_CONST_LAMBDA)
+#include "./cooling/const_lambda/cooling_struct.h"
+#elif defined(COOLING_GRACKLE)
+#include "./cooling/grackle/cooling_struct.h"
+#else
+#error "Invalid choice of cooling function."
+#endif
+
+#endif /* SWIFT_COOLING_STRUCT_H */
diff --git a/src/engine.c b/src/engine.c
index 694528ed57d5692cbb00e7d30ba1364cfdc6dedd..6fed485fa52d6ac7d5346a17122dae9254460dfb 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -2418,7 +2418,8 @@ void engine_print_stats(struct engine *e) {
const ticks tic = getticks();
const struct space *s = e->s;
- double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, entropy = 0.0, mass = 0.0;
+ double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, e_rad = 0.0;
+ double entropy = 0.0, mass = 0.0;
double mom[3] = {0.0, 0.0, 0.0}, ang_mom[3] = {0.0, 0.0, 0.0};
/* Collect the cell data. */
@@ -2429,6 +2430,7 @@ void engine_print_stats(struct engine *e) {
e_kin += c->e_kin;
e_int += c->e_int;
e_pot += c->e_pot;
+ e_rad += c->e_rad;
entropy += c->entropy;
mom[0] += c->mom[0];
mom[1] += c->mom[1];
@@ -2441,33 +2443,35 @@ void engine_print_stats(struct engine *e) {
/* Aggregate the data from the different nodes. */
#ifdef WITH_MPI
{
- double in[11] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
- double out[11];
+ double in[12] = {0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.};
+ double out[12];
out[0] = e_kin;
out[1] = e_int;
out[2] = e_pot;
- out[3] = mom[0];
- out[4] = mom[1];
- out[5] = mom[2];
- out[6] = ang_mom[0];
- out[7] = ang_mom[1];
- out[8] = ang_mom[2];
- out[9] = mass;
- out[10] = entropy;
+ out[3] = e_rad;
+ out[4] = mom[0];
+ out[5] = mom[1];
+ out[6] = mom[2];
+ out[7] = ang_mom[0];
+ out[8] = ang_mom[1];
+ out[9] = ang_mom[2];
+ out[10] = mass;
+ out[11] = entropy;
if (MPI_Reduce(out, in, 11, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD) !=
MPI_SUCCESS)
error("Failed to aggregate stats.");
e_kin = out[0];
e_int = out[1];
e_pot = out[2];
- mom[0] = out[3];
- mom[1] = out[4];
- mom[2] = out[5];
- ang_mom[0] = out[6];
- ang_mom[1] = out[7];
- ang_mom[2] = out[8];
- mass = out[9];
- entropy = out[10];
+ e_rad = out[3];
+ mom[0] = out[4];
+ mom[1] = out[5];
+ mom[2] = out[6];
+ ang_mom[0] = out[7];
+ ang_mom[1] = out[8];
+ ang_mom[2] = out[9];
+ mass = out[10];
+ entropy = out[11];
}
#endif
@@ -2475,11 +2479,11 @@ void engine_print_stats(struct engine *e) {
/* Print info */
if (e->nodeID == 0) {
- fprintf(
- e->file_stats,
- " %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e\n",
- e->time, mass, e_tot, e_kin, e_int, e_pot, entropy, mom[0], mom[1],
- mom[2], ang_mom[0], ang_mom[1], ang_mom[2]);
+ fprintf(e->file_stats,
+ " %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e %14e "
+ "%14e\n",
+ e->time, mass, e_tot, e_kin, e_int, e_pot, e_rad, entropy, mom[0],
+ mom[1], mom[2], ang_mom[0], ang_mom[1], ang_mom[2]);
fflush(e->file_stats);
}
@@ -3101,7 +3105,7 @@ void engine_unpin() {
* @param physical_constants The #phys_const used for this run.
* @param hydro The #hydro_props used for this run.
* @param potential The properties of the external potential.
- * @param cooling The properties of the cooling function.
+ * @param cooling_func The properties of the cooling function.
*/
void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, int nr_nodes, int nodeID,
@@ -3110,7 +3114,7 @@ void engine_init(struct engine *e, struct space *s,
const struct phys_const *physical_constants,
const struct hydro_props *hydro,
const struct external_potential *potential,
- const struct cooling_data *cooling) {
+ const struct cooling_function_data *cooling_func) {
/* Clean-up everything */
bzero(e, sizeof(struct engine));
@@ -3162,7 +3166,7 @@ void engine_init(struct engine *e, struct space *s,
e->physical_constants = physical_constants;
e->hydro_properties = hydro;
e->external_potential = potential;
- e->cooling_data = cooling;
+ e->cooling_func = cooling_func;
e->parameter_file = params;
engine_rank = nodeID;
@@ -3313,11 +3317,11 @@ void engine_init(struct engine *e, struct space *s,
engine_default_energy_file_name);
sprintf(energyfileName + strlen(energyfileName), ".txt");
e->file_stats = fopen(energyfileName, "w");
- fprintf(
- e->file_stats,
- "#%14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s\n",
- "Time", "Mass", "E_tot", "E_kin", "E_int", "E_pot", "Entropy", "p_x",
- "p_y", "p_z", "ang_x", "ang_y", "ang_z");
+ fprintf(e->file_stats,
+ "#%14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s %14s "
+ "%14s\n",
+ "Time", "Mass", "E_tot", "E_kin", "E_int", "E_pot", "E_radcool",
+ "Entropy", "p_x", "p_y", "p_z", "ang_x", "ang_y", "ang_z");
fflush(e->file_stats);
char timestepsfileName[200] = "";
diff --git a/src/engine.h b/src/engine.h
index c35163a44e1c6fa4e399fb421bfbfd0224b9486c..d36914af611723bc4d496d5c2dc68050eea6ffe6 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -38,10 +38,10 @@
/* Includes. */
#include "clocks.h"
-#include "cooling.h"
+#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"
@@ -205,7 +205,7 @@ struct engine {
const struct external_potential *external_potential;
/* Properties of the cooling scheme */
- const struct cooling_data *cooling_data;
+ const struct cooling_function_data *cooling_func;
/* The (parsed) parameter file */
const struct swift_params *parameter_file;
@@ -223,7 +223,7 @@ void engine_init(struct engine *e, struct space *s,
const struct phys_const *physical_constants,
const struct hydro_props *hydro,
const struct external_potential *potential,
- const struct cooling_data *cooling);
+ const struct cooling_function_data *cooling);
void engine_launch(struct engine *e, int nr_runners, unsigned int mask,
unsigned int submask);
void engine_prepare(struct engine *e, int nodrift);
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 2415e20ac5eb68f1b773b990bab232707166903c..9e0ca81edff06b8a32afb185f24a88b41dc87da7 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -22,47 +22,14 @@
#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;
-}
/**
* @brief Computes the gravity time-step of a given particle due to self-gravity
*
- * @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_self(const struct phys_const* const phys_const,
- const struct gpart* const gp) {
+gravity_compute_timestep_self(const struct gpart* const gp) {
const float ac2 = gp->a_grav[0] * gp->a_grav[0] +
gp->a_grav[1] * gp->a_grav[1] +
@@ -125,31 +92,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/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index 1850ff0a1644d3593f78f150646eae8b2f074e1e..f06e65e5b30ebcd609c0c6204de33da17b770add 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -53,6 +53,6 @@ struct gpart {
which this gpart is linked. */
long long id_or_neg_offset;
-} __attribute__((aligned(gpart_align)));
+} SWIFT_STRUCT_ALIGN;
#endif /* SWIFT_DEFAULT_GRAVITY_PART_H */
diff --git a/src/hydro/Default/hydro_part.h b/src/hydro/Default/hydro_part.h
index f42c3dc886ae1ab8f472ffdf5ff508f6735d1bb1..c7464bcf338b1c5b81ffa91d92264c2bd35e9313 100644
--- a/src/hydro/Default/hydro_part.h
+++ b/src/hydro/Default/hydro_part.h
@@ -19,6 +19,8 @@
#ifndef SWIFT_DEFAULT_HYDRO_PART_H
#define SWIFT_DEFAULT_HYDRO_PART_H
+#include "cooling_struct.h"
+
/* Extra particle data not needed during the SPH loops over neighbours. */
struct xpart {
@@ -28,12 +30,15 @@ struct xpart {
/* Velocity at the last full step. */
float v_full[3];
+ /* Additional data used to record cooling information */
+ struct cooling_xpart_data cooling_data;
+
float u_full;
/* Old density. */
float omega;
-} __attribute__((aligned(xpart_align)));
+} SWIFT_STRUCT_ALIGN;
/* Data of a single particle. */
struct part {
@@ -120,6 +125,6 @@ struct part {
/* Pointer to corresponding gravity part. */
struct gpart* gpart;
-} __attribute__((aligned(part_align)));
+} SWIFT_STRUCT_ALIGN;
#endif /* SWIFT_DEFAULT_HYDRO_PART_H */
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 09a8f50d5b2e9abd43c3b9bd43a12fad8a347258..2b17505c36f345779b269c90758931a13f9b4e0d 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -19,6 +19,18 @@
#ifndef SWIFT_GADGET2_HYDRO_H
#define SWIFT_GADGET2_HYDRO_H
+/**
+ * @file Gadget2/hydro.h
+ * @brief SPH interaction functions following the Gadget-2 version of SPH.
+ *
+ * The interactions computed here are the ones presented in the Gadget-2 paper
+ * Springel, V., MNRAS, Volume 364, Issue 4, pp. 1105-1134.
+ * We use the same numerical coefficients as the Gadget-2 code. When used with
+ * the Spline-3 kernel, the results should be equivalent to the ones obtained
+ * with Gadget-2 up to the rounding errors and interactions missed by the
+ * Gadget-2 tree-code neighbours search.
+ */
+
#include "adiabatic_index.h"
#include "approx_math.h"
#include "dimension.h"
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 8a4edfe62f59a3fae551fdb65f46987509f89251..5ee8cd0370a970aa83cef3d1c8909d923c12ba24 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -20,21 +20,20 @@
#ifndef SWIFT_GADGET2_HYDRO_IACT_H
#define SWIFT_GADGET2_HYDRO_IACT_H
-#include "minmax.h"
-
/**
+ * @file Gadget2/hydro_iact.h
* @brief SPH interaction functions following the Gadget-2 version of SPH.
*
* The interactions computed here are the ones presented in the Gadget-2 paper
- *and use the same
- * numerical coefficients as the Gadget-2 code. When used with the Spline-3
- *kernel, the results
- * should be equivalent to the ones obtained with Gadget-2 up to the rounding
- *errors and interactions
- * missed by the Gadget-2 tree-code neighbours search.
- *
+ * Springel, V., MNRAS, Volume 364, Issue 4, pp. 1105-1134.
+ * We use the same numerical coefficients as the Gadget-2 code. When used with
+ * the Spline-3 kernel, the results should be equivalent to the ones obtained
+ * with Gadget-2 up to the rounding errors and interactions missed by the
+ * Gadget-2 tree-code neighbours search.
*/
+#include "minmax.h"
+
/**
* @brief Density loop
*/
diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h
index 484792438d2717413c1ca8d4f429eac2e6d21b20..2c12a71ad6256d7372256b6e590790ee0ff4e22e 100644
--- a/src/hydro/Gadget2/hydro_part.h
+++ b/src/hydro/Gadget2/hydro_part.h
@@ -19,6 +19,20 @@
#ifndef SWIFT_GADGET2_HYDRO_PART_H
#define SWIFT_GADGET2_HYDRO_PART_H
+/**
+ * @file Gadget2/hydro_part.h
+ * @brief SPH interaction functions following the Gadget-2 version of SPH.
+ *
+ * The interactions computed here are the ones presented in the Gadget-2 paper
+ * Springel, V., MNRAS, Volume 364, Issue 4, pp. 1105-1134.
+ * We use the same numerical coefficients as the Gadget-2 code. When used with
+ * the Spline-3 kernel, the results should be equivalent to the ones obtained
+ * with Gadget-2 up to the rounding errors and interactions missed by the
+ * Gadget-2 tree-code neighbours search.
+ */
+
+#include "cooling_struct.h"
+
/* Extra particle data not needed during the SPH loops over neighbours. */
struct xpart {
@@ -28,7 +42,10 @@ struct xpart {
/* Velocity at the last full step. */
float v_full[3];
-} __attribute__((aligned(xpart_align)));
+ /* Additional data used to record cooling information */
+ struct cooling_xpart_data cooling_data;
+
+} SWIFT_STRUCT_ALIGN;
/* Data of a single particle. */
struct part {
@@ -110,6 +127,6 @@ struct part {
/* Pointer to corresponding gravity part. */
struct gpart* gpart;
-} __attribute__((aligned(part_align)));
+} SWIFT_STRUCT_ALIGN;
#endif /* SWIFT_GADGET2_HYDRO_PART_H */
diff --git a/src/hydro/Gizmo/hydro_part.h b/src/hydro/Gizmo/hydro_part.h
index d425294671d4bc172f45c928c2290f8cfa8e093c..c4919ff173c64a4a83a5d1bf61ab82697cc03096 100644
--- a/src/hydro/Gizmo/hydro_part.h
+++ b/src/hydro/Gizmo/hydro_part.h
@@ -16,6 +16,10 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
+#ifndef SWIFT_GIZMO_HYDRO_PART_H
+#define SWIFT_GIZMO_HYDRO_PART_H
+
+#include "cooling_struct.h"
/* Extra particle data not needed during the computation. */
struct xpart {
@@ -26,7 +30,10 @@ struct xpart {
/* Velocity at the last full step. */
float v_full[3];
-} __attribute__((aligned(xpart_align)));
+ /* Additional data used to record cooling information */
+ struct cooling_xpart_data cooling_data;
+
+} SWIFT_STRUCT_ALIGN;
/* Data of a single particle. */
struct part {
@@ -196,4 +203,6 @@ struct part {
/* Associated gravitas. */
struct gpart *gpart;
-} __attribute__((aligned(part_align)));
+} SWIFT_STRUCT_ALIGN;
+
+#endif /* SWIFT_GIZMO_HYDRO_PART_H */
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 0bbc77f4a2384c79f7d3329c20b92990598d5c63..eeb389537c56876126c60b2d29a728029c72844b 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -29,9 +29,8 @@
* term is implemented.
*
* This corresponds to equations (43), (44), (45), (101), (103) and (104) with
- * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
- * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
- * pp. 759-794.
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of Price, D., Journal of Computational
+ * Physics, 2012, Volume 231, Issue 3, pp. 759-794.
*/
#include "adiabatic_index.h"
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index 9e2028c978dc7cad03cfba17931f645bbfbfe1a0..edb060e4fd71fcc136e1bedf6e8a752d1d50d54f 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -28,9 +28,8 @@
* term is implemented.
*
* This corresponds to equations (43), (44), (45), (101), (103) and (104) with
- * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
- * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
- * pp. 759-794.
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of Price, D., Journal of Computational
+ * Physics, 2012, Volume 231, Issue 3, pp. 759-794.
*/
#include "adiabatic_index.h"
diff --git a/src/hydro/Minimal/hydro_part.h b/src/hydro/Minimal/hydro_part.h
index ad65f8b44fc67f4aae6470246cbab91bc3710007..8e23bddf5153043421319810683266b27d297f93 100644
--- a/src/hydro/Minimal/hydro_part.h
+++ b/src/hydro/Minimal/hydro_part.h
@@ -28,11 +28,12 @@
* term is implemented.
*
* This corresponds to equations (43), (44), (45), (101), (103) and (104) with
- * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of
- * Price, D., Journal of Computational Physics, 2012, Volume 231, Issue 3,
- * pp. 759-794.
+ * \f$\beta=3\f$ and \f$\alpha_u=0\f$ of Price, D., Journal of Computational
+ * Physics, 2012, Volume 231, Issue 3, pp. 759-794.
*/
+#include "cooling_struct.h"
+
/**
* @brief Particle fields not needed during the SPH loops over neighbours.
*
@@ -42,12 +43,16 @@
*/
struct xpart {
- float x_diff[3]; /*!< Offset between current position and position at last
- tree rebuild. */
+ /*! Offset between current position and position at last tree rebuild. */
+ float x_diff[3];
+
+ /*! Velocity at the last full step. */
+ float v_full[3];
- float v_full[3]; /*!< Velocity at the last full step. */
+ /*! Additional data used to record cooling information */
+ struct cooling_xpart_data cooling_data;
-} __attribute__((aligned(xpart_align)));
+} SWIFT_STRUCT_ALIGN;
/**
* @brief Particle fields for the SPH particles
@@ -58,27 +63,38 @@ struct xpart {
*/
struct part {
- double x[3]; /*!< Particle position. */
+ /*! Particle position. */
+ double x[3];
- float v[3]; /*!< Particle predicted velocity. */
+ /*! Particle predicted velocity. */
+ float v[3];
- float a_hydro[3]; /*!< Particle acceleration. */
+ /*! Particle acceleration. */
+ float a_hydro[3];
- float mass; /*!< Particle mass. */
+ /*! Particle mass. */
+ float mass;
- float h; /*!< Particle smoothing length. */
+ /*! Particle smoothing length. */
+ float h;
- int ti_begin; /*!< Time at the beginning of time-step. */
+ /*! Time at the beginning of time-step. */
+ int ti_begin;
- int ti_end; /*!< Time at the end of time-step. */
+ /*! Time at the end of time-step. */
+ int ti_end;
- float u; /*!< Particle internal energy. */
+ /*! Particle internal energy. */
+ float u;
- float u_dt; /*!< Time derivative of the internal energy. */
+ /*! Time derivative of the internal energy. */
+ float u_dt;
- float rho; /*!< Particle density. */
+ /*! Particle density. */
+ float rho;
- float rho_dh; /*!< Derivative of density with respect to h */
+ /*! Derivative of density with respect to h */
+ float rho_dh;
/* Store density/force specific stuff. */
union {
@@ -92,10 +108,12 @@ struct part {
*/
struct {
- float wcount; /*!< Neighbour number count. */
+ /*! Neighbour number count. */
+ float wcount;
+
+ /*! Derivative of the neighbour number with respect to h. */
+ float wcount_dh;
- float wcount_dh; /*!< Derivative of the neighbour number with respect to
- h. */
} density;
/**
@@ -107,19 +125,24 @@ struct part {
*/
struct {
- float pressure; /*!< Particle pressure. */
+ /*! Particle pressure. */
+ float pressure;
- float v_sig; /*!< Particle signal velocity */
+ /*! Particle signal velocity */
+ float v_sig;
- float h_dt; /*!< Time derivative of smoothing length */
+ /*! Time derivative of smoothing length */
+ float h_dt;
} force;
};
- long long id; /*!< Particle unique ID. */
+ /*! Particle unique ID. */
+ long long id;
- struct gpart* gpart; /*!< Pointer to corresponding gravity part. */
+ /*! Pointer to corresponding gravity part. */
+ struct gpart* gpart;
-} __attribute__((aligned(part_align)));
+} SWIFT_STRUCT_ALIGN;
#endif /* SWIFT_MINIMAL_HYDRO_PART_H */
diff --git a/src/part.h b/src/part.h
index ea895e6e0295d6a8b63309c7bd6855daa2cf7d64..af39d10fafc11d4435ddbcc087fbf08178b18959 100644
--- a/src/part.h
+++ b/src/part.h
@@ -31,12 +31,13 @@
#endif
/* Local headers. */
+#include "align.h"
#include "const.h"
/* Some constants. */
-#define part_align 64
-#define xpart_align 32
-#define gpart_align 32
+#define part_align 128
+#define xpart_align 128
+#define gpart_align 128
/* Import the right hydro particle definition */
#if defined(MINIMAL_SPH)
diff --git a/src/potential.c b/src/potential.c
new file mode 100644
index 0000000000000000000000000000000000000000..5433a05e3e7886ad88021d3916cae26adfe8b954
--- /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 0000000000000000000000000000000000000000..77bd41794a3a8cd244405493898d63b3f80ff3a6
--- /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 0000000000000000000000000000000000000000..21d168818e164ad3b3e18076ba824285e40956aa
--- /dev/null
+++ b/src/potential/disc_patch/potential.h
@@ -0,0 +1,199 @@
+/*******************************************************************************
+ * 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"
+
+/**
+ * @brief External Potential Properties - Disc patch case
+ *
+ * See Creasey, Theuns & Bower, 2013, MNRAS, Volume 429, Issue 3, p.1922-1948
+ */
+struct external_potential {
+
+ /*! Surface density of the disc */
+ double surface_density;
+
+ /*! Disc scale-height */
+ double scale_height;
+
+ /*! Position of the disc along the z-axis */
+ double z_disc;
+
+ /*! Dynamical time of the system */
+ double dynamical_time;
+
+ /*! Time over which to grow the disk in units of the dynamical time */
+ double growth_time;
+
+ /*! Time-step condition pre-factor */
+ 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.f * 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.f * 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 dz = g->x[2] - potential->z_disc;
+ const float t_dyn = potential->dynamical_time;
+
+ float reduction_factor = 1.f;
+ if (time < potential->growth_time * t_dyn)
+ reduction_factor = time / (potential->growth_time * t_dyn);
+
+ /* Accelerations. Note that they are multiplied by G later on */
+ const float z_accel = reduction_factor * 2.f * M_PI *
+ potential->surface_density *
+ tanh(fabs(dz) / potential->scale_height);
+
+ if (dz > 0) g->a_grav[2] -= z_accel;
+ if (dz < 0) g->a_grav[2] += z_accel;
+}
+
+/**
+ * @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 dynamical 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 0000000000000000000000000000000000000000..a993c09a978ca3692ec3359f7633df14760f263d
--- /dev/null
+++ b/src/potential/isothermal/potential.h
@@ -0,0 +1,162 @@
+/*******************************************************************************
+ * 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"
+
+/**
+ * @brief External Potential Properties - Isothermal sphere case
+ */
+struct external_potential {
+
+ /*! Position of the centre of potential */
+ double x, y, z;
+
+ /*! Rotation velocity */
+ double vrot;
+
+ /*! Square of vrot divided by G \f$ \frac{v_{rot}^2}{G} \f$ */
+ double vrot2_over_G;
+
+ /*! Time-step condition pre-factor */
+ double 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.f * drdv * dx * rinv2);
+ const float dota_y =
+ vrot * vrot * rinv2 * (g->v_full[1] - 2.f * drdv * dy * rinv2);
+ const float dota_z =
+ vrot * vrot * rinv2 * (g->v_full[2] - 2.f * 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 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 term = -potential->vrot2_over_G * rinv2;
+
+ 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");
+
+ potential->vrot2_over_G =
+ potential->vrot * potential->vrot / phys_const->const_newton_G;
+}
+
+/**
+ * @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 0000000000000000000000000000000000000000..8b1c3e841521f3fb42fbdf5c8922cead2ea7cbcb
--- /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 0000000000000000000000000000000000000000..f718af2e2c4ff91540e1834cb2072d321ce38705
--- /dev/null
+++ b/src/potential/point_mass/potential.h
@@ -0,0 +1,159 @@
+/*******************************************************************************
+ * 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 - Point mass case
+ */
+struct external_potential {
+
+ /*! Position of the point mass */
+ double x, y, z;
+
+ /*! Mass */
+ double mass;
+
+ /*! Time-step condition pre-factor */
+ 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 dd7aed8712c01921a01462bf9de713433c81f5c1..0000000000000000000000000000000000000000
--- 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 74f0fd28566f355962c83e5d743aeae9afe09c59..0000000000000000000000000000000000000000
--- 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 859af9199a4281f9671578948a43254534d87839..87d7ca8f5d16ffa9814c3eef55b7bc0794aab7cc 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);
}
}
@@ -170,9 +170,10 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
struct part *restrict parts = c->parts;
+ struct xpart *restrict xparts = c->xparts;
const int count = c->count;
const int ti_current = r->e->ti_current;
- const struct cooling_data *cooling = r->e->cooling_data;
+ const struct cooling_function_data *cooling_func = r->e->cooling_func;
const struct phys_const *constants = r->e->physical_constants;
const struct UnitSystem *us = r->e->internalUnits;
const double timeBase = r->e->timeBase;
@@ -195,13 +196,14 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
/* Get a direct pointer on the part. */
struct part *restrict p = &parts[i];
+ struct xpart *restrict xp = &xparts[i];
/* Kick has already updated ti_end, so need to check ti_begin */
if (p->ti_begin == ti_current) {
const double dt = (p->ti_end - p->ti_begin) * timeBase;
- cooling_cool_part(constants, us, cooling, p, dt);
+ cooling_cool_part(constants, us, cooling_func, p, xp, dt);
}
}
@@ -729,11 +731,15 @@ static void runner_do_drift(struct cell *c, struct engine *e) {
/* Do we need to drift ? */
if (!e->drift_all && !cell_is_drift_needed(c, ti_current)) return;
+ /* Check that we are actually going to move forward. */
+ if (ti_current == ti_old) return;
+
/* Drift from the last time the cell was drifted to the current time */
const double dt = (ti_current - ti_old) * timeBase;
float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
- double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, entropy = 0.0, mass = 0.0;
+ double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, e_rad = 0.0;
+ double entropy = 0.0, mass = 0.0;
double mom[3] = {0.0, 0.0, 0.0};
double ang_mom[3] = {0.0, 0.0, 0.0};
@@ -808,6 +814,7 @@ static void runner_do_drift(struct cell *c, struct engine *e) {
e_kin += 0.5 * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
e_pot += 0.;
e_int += m * hydro_get_internal_energy(p, half_dt);
+ e_rad += cooling_get_radiated_energy(xp);
/* Collect entropy */
entropy += m * hydro_get_entropy(p, half_dt);
@@ -836,6 +843,7 @@ static void runner_do_drift(struct cell *c, struct engine *e) {
e_kin += cp->e_kin;
e_int += cp->e_int;
e_pot += cp->e_pot;
+ e_rad += cp->e_rad;
entropy += cp->entropy;
mom[0] += cp->mom[0];
mom[1] += cp->mom[1];
@@ -853,6 +861,7 @@ static void runner_do_drift(struct cell *c, struct engine *e) {
c->e_kin = e_kin;
c->e_int = e_int;
c->e_pot = e_pot;
+ c->e_rad = e_rad;
c->entropy = entropy;
c->mom[0] = mom[0];
c->mom[1] = mom[1];
@@ -1316,7 +1325,6 @@ void *runner_main(void *data) {
break;
case task_type_sort:
runner_do_sort(r, ci, t->flags, 1);
- t->flags = 0;
break;
case task_type_sub_self:
if (t->subtype == task_subtype_density)
diff --git a/src/space.c b/src/space.c
index cdd5958cbc515003f4a86a41c9a7075fa3b4364f..a9958f6fbd7d85060db99a9682b0de10f507085d 100644
--- a/src/space.c
+++ b/src/space.c
@@ -42,6 +42,7 @@
/* Local headers. */
#include "atomic.h"
#include "const.h"
+#include "cooling.h"
#include "engine.h"
#include "error.h"
#include "gravity.h"
@@ -1471,6 +1472,23 @@ void space_init_parts(struct space *s) {
}
}
+/**
+ * @brief Initialises all the extra particle data
+ *
+ * Calls cooling_init_xpart() on all the particles
+ */
+void space_init_xparts(struct space *s) {
+
+ const size_t nr_parts = s->nr_parts;
+ struct part *restrict p = s->parts;
+ struct xpart *restrict xp = s->xparts;
+
+ for (size_t i = 0; i < nr_parts; ++i) {
+
+ cooling_init_part(&p[i], &xp[i]);
+ }
+}
+
/**
* @brief Initialises all the g-particles by setting them into a valid state
*
@@ -1630,6 +1648,7 @@ void space_init(struct space *s, const struct swift_params *params,
/* Set the particles in a state where they are ready for a run */
space_init_parts(s);
+ space_init_xparts(s);
space_init_gparts(s);
/* Init the space lock. */
diff --git a/src/swift.h b/src/swift.h
index 80a4b1a1f792ceab3184d671f9aace277a36f419..3119adfd4b95639a84137e243c6e971d80ef3825 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/task.c b/src/task.c
index 4bc50cd0b3daac48edea7ea3040c04e8c42dbb02..8da526f57886fa68121061c8eae2316912e73a30 100644
--- a/src/task.c
+++ b/src/task.c
@@ -63,7 +63,7 @@ const char *subtaskID_names[task_subtype_count] = {
* @param cj The second #cell.
*/
__attribute__((always_inline)) INLINE static size_t task_cell_overlap_part(
- const struct cell *ci, const struct cell *cj) {
+ const struct cell *restrict ci, const struct cell *restrict cj) {
if (ci == NULL || cj == NULL) return 0;
@@ -85,7 +85,7 @@ __attribute__((always_inline)) INLINE static size_t task_cell_overlap_part(
* @param cj The second #cell.
*/
__attribute__((always_inline)) INLINE static size_t task_cell_overlap_gpart(
- const struct cell *ci, const struct cell *cj) {
+ const struct cell *restrict ci, const struct cell *restrict cj) {
if (ci == NULL || cj == NULL) return 0;
@@ -168,6 +168,10 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
return task_action_none;
break;
}
+
+ /* Silence compile warnings */
+ error("Unknown task_action for task");
+ return task_action_none;
}
/**
@@ -177,7 +181,8 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
* @param ta The first #task.
* @param tb The second #task.
*/
-float task_overlap(const struct task *ta, const struct task *tb) {
+float task_overlap(const struct task *restrict ta,
+ const struct task *restrict tb) {
if (ta == NULL || tb == NULL) return 0.f;
diff --git a/src/task.h b/src/task.h
index f070451fe4e79e0c16dc3dcca1ce145c08841c47..4928dc00bcd7958efdd6987b5aec90ab4b3e92fa 100644
--- a/src/task.h
+++ b/src/task.h
@@ -26,6 +26,7 @@
#include "../config.h"
/* Includes. */
+#include "align.h"
#include "cell.h"
#include "cycle.h"
@@ -149,7 +150,7 @@ struct task {
/*! Is this task implicit (i.e. does not do anything) ? */
char implicit;
-} __attribute__((aligned(32)));
+} SWIFT_STRUCT_ALIGN;
/* Function prototypes. */
void task_unlock(struct task *t);
diff --git a/src/timestep.h b/src/timestep.h
index d92f88d06451892ce47db4b9468db9714bb52baa..599fb4762e11b08fc942fb02acbbf1970f477de4 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -68,18 +68,18 @@ __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
- float new_dt =
- (new_dt_external < new_dt_self) ? new_dt_external : new_dt_self;
+ float new_dt = min(new_dt_external, new_dt_self);
/* Limit timestep within the allowed range */
- new_dt = (new_dt < e->dt_max) ? new_dt : e->dt_max;
- new_dt = (new_dt > e->dt_min) ? new_dt : e->dt_min;
+ new_dt = min(new_dt, e->dt_max);
+ new_dt = max(new_dt, e->dt_min);
/* Convert to integer time */
const int new_dti =
@@ -102,22 +102,28 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
/* Compute the next timestep (hydro condition) */
const float new_dt_hydro = hydro_compute_timestep(p, xp, e->hydro_properties);
+ /* Compute the next timestep (cooling condition) */
+ float new_dt_cooling = FLT_MAX;
+ if (e->policy & engine_policy_cooling)
+ new_dt_cooling = cooling_timestep(e->cooling_func, e->physical_constants,
+ e->internalUnits, p);
+
/* Compute the next timestep (gravity condition) */
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
- new_dt_grav =
- (new_dt_external < new_dt_self) ? new_dt_external : new_dt_self;
+ new_dt_grav = min(new_dt_external, new_dt_self);
}
/* Final time-step is minimum of hydro and gravity */
- float new_dt = (new_dt_hydro < new_dt_grav) ? new_dt_hydro : new_dt_grav;
+ float new_dt = min(min(new_dt_hydro, new_dt_cooling), new_dt_grav);
/* Limit change in h */
const float dt_h_change =
@@ -125,11 +131,11 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
? fabsf(e->hydro_properties->log_max_h_change * p->h / p->force.h_dt)
: FLT_MAX;
- new_dt = (new_dt < dt_h_change) ? new_dt : dt_h_change;
+ new_dt = min(new_dt, dt_h_change);
/* Limit timestep within the allowed range */
- new_dt = (new_dt < e->dt_max) ? new_dt : e->dt_max;
- new_dt = (new_dt > e->dt_min) ? new_dt : e->dt_min;
+ new_dt = min(new_dt, e->dt_max);
+ new_dt = max(new_dt, e->dt_min);
/* Convert to integer time */
const int new_dti =
diff --git a/src/units.c b/src/units.c
index f598b5ddf0b1a4b165648d5378915cd6f10f0bba..2241d441ec9af9b6d5083191e8f61010ebaccb20 100644
--- a/src/units.c
+++ b/src/units.c
@@ -319,6 +319,9 @@ void units_get_base_unit_exponants_array(float baseUnitsExp[5],
break;
case UNIT_CONV_VOLUME:
+ baseUnitsExp[UNIT_LENGTH] = 3.f;
+
+ case UNIT_CONV_INV_VOLUME:
baseUnitsExp[UNIT_LENGTH] = -3.f;
}
}
diff --git a/src/units.h b/src/units.h
index 26fa15a66528dd39ea232cdf94da2ff0230300cd..78fdf1c23c3c276607d5353ee3437d8eb1e96537 100644
--- a/src/units.h
+++ b/src/units.h
@@ -90,7 +90,8 @@ enum UnitConversionFactor {
UNIT_CONV_MAGNETIC_FIELD,
UNIT_CONV_MAGNETIC_INDUCTANCE,
UNIT_CONV_TEMPERATURE,
- UNIT_CONV_VOLUME
+ UNIT_CONV_VOLUME,
+ UNIT_CONV_INV_VOLUME
};
void units_init_cgs(struct UnitSystem*);
diff --git a/tests/test125cells.c b/tests/test125cells.c
index a385a7c890fe27ed11d3c5d87d6903fa6d254516..e666658f43de135e3e72521b52f2a688c596a6f6 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -519,6 +519,7 @@ int main(int argc, char *argv[]) {
hp.CFL_condition = 0.1;
struct engine engine;
+ bzero(&engine, sizeof(struct engine));
engine.hydro_properties = &hp;
engine.physical_constants = &prog_const;
engine.s = &space;
diff --git a/tests/testSPHStep.c b/tests/testSPHStep.c
index fa49ed9d00c37393abd2f7e17ae628d79b4125f6..ff2ec841b27bd5ca6190517bc39f4da0c28fbc0c 100644
--- a/tests/testSPHStep.c
+++ b/tests/testSPHStep.c
@@ -102,10 +102,6 @@ int main() {
int i, j, k, offset[3];
struct part *p;
- struct hydro_props hp;
- hp.target_neighbours = 48.;
- hp.delta_neighbours = 1.;
- hp.max_smoothing_iterations = 30;
int N = 10;
float dim = 1.;
@@ -146,11 +142,24 @@ int main() {
/* Create the infrastructure */
struct space space;
space.periodic = 0;
- space.h_max = 1.;
+ space.cell_min = 1.;
+
+ struct phys_const prog_const;
+ prog_const.const_newton_G = 1.f;
+
+ struct hydro_props hp;
+ hp.target_neighbours = 48.f;
+ hp.delta_neighbours = 2.;
+ hp.max_smoothing_iterations = 1;
+ hp.CFL_condition = 0.1;
struct engine e;
- e.s = &space;
+ bzero(&e, sizeof(struct engine));
e.hydro_properties = &hp;
+ e.physical_constants = &prog_const;
+ e.s = &space;
+ e.time = 0.1f;
+ e.ti_current = 1;
struct runner r;
r.e = &e;