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/EAGLE_25/README b/examples/EAGLE_25/README
index 077fd9cf06ce64be98fa0d8a736125c474fb7a76..88fc1ea3eede1e254907dd5ba1dbf2eaa81fb694 100644
--- a/examples/EAGLE_25/README
+++ b/examples/EAGLE_25/README
@@ -13,4 +13,4 @@ The particle load of the main EAGLE simulation can be reproduced by
running these ICs on 64 cores.
MD5 checksum of the ICs:
-ada2c728db2bd2d77a20c4eef52dfaf1 EAGLE_ICs_25.hdf5
+02cd1c353b86230af047b5d4ab22afcf EAGLE_ICs_25.hdf5
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 28d02fefa8168e35af696975a7c73a1bf767155e..51a6d2b478681e2e1c61e199f758e35c507ec195 100644
--- a/examples/MultiTypes/multiTypes.yml
+++ b/examples/MultiTypes/multiTypes.yml
@@ -35,8 +35,9 @@ 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.
mass: 1e10 # mass of external point mass in internal units
+ timestep_mult: 1e-2
diff --git a/examples/SedovBlast_1D/sedov.yml b/examples/SedovBlast_1D/sedov.yml
index 6f519835d26ff5aa851ffb8999e650815c522cd3..1ecfeb32452d05f299b98124c4fdfc79126f7504 100644
--- a/examples/SedovBlast_1D/sedov.yml
+++ b/examples/SedovBlast_1D/sedov.yml
@@ -21,7 +21,7 @@ Snapshots:
# Parameters governing the conserved quantities statistics
Statistics:
- delta_time: 1e-3 # Time between statistics output
+ delta_time: 1e-5 # Time between statistics output
# Parameters for the hydrodynamics scheme
SPH:
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 8591c0c70010075e62be73d29ac3644bc8b1f315..15d0fdb78fede79c24efc49f19082e25d58694ae 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -541,7 +541,7 @@ int main(int argc, char *argv[]) {
/* Make sure output file is empty, only on one rank. */
char dumpfile[30];
- snprintf(dumpfile, 30, "thread_info_MPI-step%d.dat", j);
+ snprintf(dumpfile, 30, "thread_info_MPI-step%d.dat", j + 1);
FILE *file_thread;
if (myrank == 0) {
file_thread = fopen(dumpfile, "w");
@@ -593,7 +593,7 @@ int main(int argc, char *argv[]) {
#else
char dumpfile[30];
- snprintf(dumpfile, 30, "thread_info-step%d.dat", j);
+ snprintf(dumpfile, 30, "thread_info-step%d.dat", j + 1);
FILE *file_thread;
file_thread = fopen(dumpfile, "w");
/* Add some information to help with the plots */
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 10bd09a18f7e8099e4db559e76957d19d8f90164..fb193a81baf410360172793763875891f654211a 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -12,6 +12,7 @@ Scheduler:
cell_max_size: 8000000 # (Optional) Maximal number of interactions per task if we force the split (this is the default value).
cell_sub_size: 64000000 # (Optional) Maximal number of interactions per sub-task (this is the default value).
cell_split_size: 400 # (Optional) Maximal number of particles per cell (this is the default value).
+ cell_max_count: 10000 # (Optional) Maximal number of particles per cell allowed before triggering a sanitizing (this is the default value).
# Parameters governing the time integration
TimeIntegration:
@@ -66,7 +67,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 +83,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/examples/plot_scaling_results.py b/examples/plot_scaling_results.py
index 5a76e9870bd3ec55807c7b79c475c62b14119e5c..938d204bea8a9da81a93cf2d6d7334ac99d10e3a 100755
--- a/examples/plot_scaling_results.py
+++ b/examples/plot_scaling_results.py
@@ -131,8 +131,10 @@ def parse_files():
file_list = [ file_list[j] for j in sorted_indices]
parse_header(file_list[0])
+
+ branch[i] = branch[i].replace("_", "\\_")
- version.append(branch[i] + " " + revision[i] + "\n" + hydro_scheme[i] +
+ version.append("$\\textrm{%s}$"%str(branch[i]) + " " + revision[i] + "\n" + hydro_scheme[i] +
"\n" + hydro_kernel[i] + r", $N_{ngb}=%d$"%float(hydro_neighbours[i]) +
r", $\eta=%.3f$"%float(hydro_eta[i]))
times.append([])
@@ -215,14 +217,21 @@ def plot_results(times,totalTime,speedUp,parallelEff):
pts = [1, 10**np.floor(np.log10(threadList[i][-1])+1)]
totalTimePlot.loglog(pts,totalTime[i][0]/pts, 'k--', lw=1., color='0.2')
totalTimePlot.loglog(threadList[i],totalTime[i],linestyle[i],label=version[i])
-
+
+ y_min = 10**np.floor(np.log10(np.min(totalTime[:][-1])*0.6))
+ y_max = 1.2*10**np.floor(np.log10(np.max(totalTime[:][0]) * 1.5)+1)
totalTimePlot.set_xscale('log')
totalTimePlot.set_xlabel("${\\rm Threads}$", labelpad=0.)
totalTimePlot.set_ylabel("${\\rm Time~to~solution}~[{\\rm ms}]$", labelpad=0.)
totalTimePlot.set_xlim([0.9, 10**(np.floor(np.log10(threadList[i][-1]))+0.5)])
- totalTimePlot.set_ylim([10**np.floor(np.log10(np.min(totalTime)*0.6)), 1.2*10**np.floor(np.log10(np.max(totalTime) * 1.5)+1)])
+ totalTimePlot.set_ylim(y_min, y_max)
+
+ ax2 = totalTimePlot.twinx()
+ ax2.set_yscale('log')
+ ax2.set_ylabel("${\\rm Time~to~solution}~[{\\rm hr}]$", labelpad=0.)
+ ax2.set_ylim(y_min / 3.6e6, y_max / 3.6e6)
- totalTimePlot.legend(bbox_to_anchor=(1.14, 0.97), loc=2, borderaxespad=0.,prop={'size':12}, frameon=False)
+ totalTimePlot.legend(bbox_to_anchor=(1.21, 0.97), loc=2, borderaxespad=0.,prop={'size':12}, frameon=False)
emptyPlot.axis('off')
for i, txt in enumerate(threadList[0]):
diff --git a/examples/plot_tasks.py b/examples/plot_tasks.py
index fb8b2ce57a47b4d397284bba9960b098c1e3ce62..7b4683422725f206a3f582e00d82712c7e3c3f59 100755
--- a/examples/plot_tasks.py
+++ b/examples/plot_tasks.py
@@ -56,8 +56,9 @@ pl.rcParams.update(PLOT_PARAMS)
# Tasks and subtypes. Indexed as in tasks.h.
TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair", "init", "ghost",
- "kick", "kick_fixdt", "send", "recv", "grav_gather_m", "grav_fft",
- "grav_mm", "grav_up", "grav_external", "count"]
+ "extra_ghost", "kick", "kick_fixdt", "send", "recv",
+ "grav_gather_m", "grav_fft", "grav_mm", "grav_up",
+ "grav_external", "cooling", "count"]
TASKCOLOURS = {"none": "black",
"sort": "lightblue",
@@ -67,6 +68,7 @@ TASKCOLOURS = {"none": "black",
"sub_pair": "navy",
"init": "indigo",
"ghost": "cyan",
+ "extra_ghost": "cyan",
"kick": "green",
"kick_fixdt": "green",
"send": "yellow",
@@ -76,6 +78,7 @@ TASKCOLOURS = {"none": "black",
"grav_mm": "mediumturquoise",
"grav_up": "mediumvioletred",
"grav_external": "darkred",
+ "cooling": "darkblue",
"count": "powerblue"}
SUBTYPES = ["none", "density", "gradient", "force", "grav", "tend", "count"]
@@ -118,7 +121,7 @@ toc_step = int(full_step[5])
CPU_CLOCK = float(full_step[-1])
data = data[1:,:]
-print "CPU frequency:", CPU_CLOCK / 1.e9
+print "CPU frequency:", CPU_CLOCK
# Avoid start and end times of zero.
data = data[data[:,4] != 0]
@@ -130,6 +133,7 @@ if delta_t == 0:
dt = max(data[:,5]) - min(data[:,4])
if dt > delta_t:
delta_t = dt
+ print "Data range: ", delta_t / CPU_CLOCK * 1000, "ms"
# Once more doing the real gather and plots this time.
start_t = tic_step
diff --git a/examples/plot_tasks_MPI.py b/examples/plot_tasks_MPI.py
index 398324cdc773d1dc4b7f26c58866c9df6469cc0b..02bc4a03510afce396429316fb4489926bc41b12 100755
--- a/examples/plot_tasks_MPI.py
+++ b/examples/plot_tasks_MPI.py
@@ -2,7 +2,7 @@
"""
Usage:
plot_tasks_MPI.py input.dat png-output-prefix [time-range-ms]
-
+
where input.dat is a thread info file for a step of an MPI run. Use the '-y
interval' flag of the swift MPI commands to create these. The output plots
will be called 'png-output-prefix<mpi-rank>.png', i.e. one each for all the
@@ -40,6 +40,8 @@ matplotlib.use("Agg")
import pylab as pl
import numpy as np
import sys
+#import warnings
+#warnings.simplefilter("error")
# Basic plot configuration.
PLOT_PARAMS = {"axes.labelsize": 10,
@@ -61,9 +63,10 @@ PLOT_PARAMS = {"axes.labelsize": 10,
pl.rcParams.update(PLOT_PARAMS)
# Tasks and subtypes. Indexed as in tasks.h.
-TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair", "init", "ghost",
- "kick", "kick_fixdt", "send", "recv", "grav_gather_m", "grav_fft",
- "grav_mm", "grav_up", "grav_external", "count"]
+TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair", "init",
+ "ghost", "extra_ghost", "kick", "kick_fixdt", "send", "recv",
+ "grav_gather_m", "grav_fft", "grav_mm", "grav_up",
+ "grav_external", "cooling", "count"]
TASKCOLOURS = {"none": "black",
"sort": "lightblue",
@@ -73,6 +76,7 @@ TASKCOLOURS = {"none": "black",
"sub_pair": "navy",
"init": "indigo",
"ghost": "cyan",
+ "extra_ghost": "cyan",
"kick": "green",
"kick_fixdt": "green",
"send": "yellow",
@@ -82,6 +86,7 @@ TASKCOLOURS = {"none": "black",
"grav_mm": "mediumturquoise",
"grav_up": "mediumvioletred",
"grav_external": "darkred",
+ "cooling": "darkblue",
"count": "powerblue"}
SUBTYPES = ["none", "density", "gradient", "force", "grav", "tend", "count"]
@@ -111,7 +116,7 @@ outbase = sys.argv[2]
delta_t = 0
if len( sys.argv ) == 4:
delta_t = int(sys.argv[3])
-
+
# Read input.
data = pl.loadtxt( infile )
@@ -121,7 +126,7 @@ tic_step = int(full_step[5])
toc_step = int(full_step[6])
CPU_CLOCK = float(full_step[-1])
-print "CPU frequency:", CPU_CLOCK / 1.e9
+print "CPU frequency:", CPU_CLOCK
nranks = int(max(data[:,0])) + 1
@@ -143,6 +148,8 @@ if delta_t == 0:
dt = max(data[:,6]) - min(data[:,5])
if dt > delta_t:
delta_t = dt
+ print "Data range: ", delta_t / CPU_CLOCK * 1000, "ms"
+
# Once more doing the real gather and plots this time.
for rank in range(nranks):
@@ -152,93 +159,107 @@ for rank in range(nranks):
tic_step = int(full_step[5])
toc_step = int(full_step[6])
data = data[1:,:]
+ typesseen = []
- start_t = tic_step
- data[:,5] -= start_t
- data[:,6] -= start_t
- end_t = (toc_step - start_t) / CPU_CLOCK * 1000
-
- tasks = {}
- tasks[-1] = []
- for i in range(nthread):
- tasks[i] = []
-
- num_lines = pl.shape(data)[0]
- for line in range(num_lines):
- thread = int(data[line,1])
- tasks[thread].append({})
- tasks[thread][-1]["type"] = TASKTYPES[int(data[line,2])]
- tasks[thread][-1]["subtype"] = SUBTYPES[int(data[line,3])]
- tic = int(data[line,5]) / CPU_CLOCK * 1000
- toc = int(data[line,6]) / CPU_CLOCK * 1000
- tasks[thread][-1]["tic"] = tic
- tasks[thread][-1]["toc"] = toc
- tasks[thread][-1]["t"] = (toc + tic)/ 2
-
- combtasks = {}
- combtasks[-1] = []
- for i in range(nthread):
- combtasks[i] = []
-
- for thread in range(nthread):
- tasks[thread] = sorted(tasks[thread], key=lambda l: l["t"])
- lasttype = ""
- types = []
- for task in tasks[thread]:
- if task["type"] not in types:
- types.append(task["type"])
- if lasttype == "" or not lasttype == task["type"]:
- combtasks[thread].append({})
- combtasks[thread][-1]["type"] = task["type"]
- combtasks[thread][-1]["subtype"] = task["subtype"]
- combtasks[thread][-1]["tic"] = task["tic"]
- combtasks[thread][-1]["toc"] = task["toc"]
- if task["type"] == "self" or task["type"] == "pair" or task["type"] == "sub":
- combtasks[thread][-1]["colour"] = SUBCOLOURS[task["subtype"]]
+ # Dummy image for ranks that have no tasks.
+ if data.size == 0:
+ print "rank ", rank, " has no tasks"
+ fig = pl.figure()
+ ax = fig.add_subplot(1,1,1)
+ ax.set_xlim(-delta_t * 0.03 * 1000 / CPU_CLOCK, delta_t * 1.03 * 1000 / CPU_CLOCK)
+ ax.set_ylim(0, nthread)
+ start_t = tic_step
+ end_t = (toc_step - start_t) / CPU_CLOCK * 1000
+ else:
+
+ start_t = tic_step
+ data[:,5] -= start_t
+ data[:,6] -= start_t
+ end_t = (toc_step - start_t) / CPU_CLOCK * 1000
+
+ tasks = {}
+ tasks[-1] = []
+ for i in range(nthread):
+ tasks[i] = []
+
+ num_lines = pl.shape(data)[0]
+ for line in range(num_lines):
+ thread = int(data[line,1])
+ tasks[thread].append({})
+ tasks[thread][-1]["type"] = TASKTYPES[int(data[line,2])]
+ tasks[thread][-1]["subtype"] = SUBTYPES[int(data[line,3])]
+ tic = int(data[line,5]) / CPU_CLOCK * 1000
+ toc = int(data[line,6]) / CPU_CLOCK * 1000
+ tasks[thread][-1]["tic"] = tic
+ tasks[thread][-1]["toc"] = toc
+ tasks[thread][-1]["t"] = (toc + tic)/ 2
+
+ combtasks = {}
+ combtasks[-1] = []
+ for i in range(nthread):
+ combtasks[i] = []
+
+ for thread in range(nthread):
+ tasks[thread] = sorted(tasks[thread], key=lambda l: l["t"])
+ lasttype = ""
+ types = []
+ for task in tasks[thread]:
+ if task["type"] not in types:
+ types.append(task["type"])
+ if lasttype == "" or not lasttype == task["type"]:
+ combtasks[thread].append({})
+ combtasks[thread][-1]["type"] = task["type"]
+ combtasks[thread][-1]["subtype"] = task["subtype"]
+ combtasks[thread][-1]["tic"] = task["tic"]
+ combtasks[thread][-1]["toc"] = task["toc"]
+ if task["type"] == "self" or task["type"] == "pair" or task["type"] == "sub":
+ combtasks[thread][-1]["colour"] = SUBCOLOURS[task["subtype"]]
+ else:
+ combtasks[thread][-1]["colour"] = TASKCOLOURS[task["type"]]
+ lasttype = task["type"]
else:
- combtasks[thread][-1]["colour"] = TASKCOLOURS[task["type"]]
- lasttype = task["type"]
- else:
- combtasks[thread][-1]["toc"] = task["toc"]
+ combtasks[thread][-1]["toc"] = task["toc"]
+
+ fig = pl.figure()
+ ax = fig.add_subplot(1,1,1)
+ ax.set_xlim(-delta_t * 0.03 * 1000 / CPU_CLOCK, delta_t * 1.03 * 1000 / CPU_CLOCK)
+ ax.set_ylim(0, nthread)
+ tictoc = np.zeros(2)
+ for i in range(nthread):
+
+ # Collect ranges and colours into arrays.
+ tictocs = np.zeros(len(combtasks[i])*2)
+ colours = np.empty(len(combtasks[i])*2, dtype='object')
+ coloursseen = []
+ j = 0
+ for task in combtasks[i]:
+ tictocs[j] = task["tic"]
+ tictocs[j+1] = task["toc"]
+ colours[j] = task["colour"]
+ colours[j+1] = task["colour"]
+ j = j + 2
+ if task["colour"] not in coloursseen:
+ coloursseen.append(task["colour"])
+
+ # Legend support, collections don't add to this.
+ if task["subtype"] != "none":
+ qtask = task["type"] + "/" + task["subtype"]
+ else:
+ qtask = task["type"]
- typesseen = []
- fig = pl.figure()
- ax = fig.add_subplot(1,1,1)
- ax.set_xlim(-delta_t * 0.03 * 1000 / CPU_CLOCK, delta_t * 1.03 * 1000 / CPU_CLOCK)
- ax.set_ylim(0, nthread)
- tictoc = np.zeros(2)
- for i in range(nthread):
-
- # Collect ranges and colours into arrays.
- tictocs = np.zeros(len(combtasks[i])*2)
- colours = np.empty(len(combtasks[i])*2, dtype='object')
- coloursseen = []
- j = 0
- for task in combtasks[i]:
- tictocs[j] = task["tic"]
- tictocs[j+1] = task["toc"]
- colours[j] = task["colour"]
- colours[j+1] = task["colour"]
- j = j + 2
- if task["colour"] not in coloursseen:
- coloursseen.append(task["colour"])
-
- # Legend support, collections don't add to this.
- if task["subtype"] != "none":
- qtask = task["type"] + "/" + task["subtype"]
- else:
- qtask = task["type"]
- if qtask not in typesseen:
- pl.plot([], [], color=task["colour"], label=qtask)
- typesseen.append(qtask)
-
- # Now plot each colour, faster to use a mask to select colour ranges.
- for colour in coloursseen:
- collection = collections.BrokenBarHCollection.span_where(tictocs, ymin=i+0.05, ymax=i+0.95,
- where=colours == colour,
- facecolor=colour,
- linewidths=0)
- ax.add_collection(collection)
+ if qtask not in typesseen:
+ pl.plot([], [], color=task["colour"], label=qtask)
+ typesseen.append(qtask)
+
+ # Now plot each colour, faster to use a mask to select colour ranges.
+ for colour in coloursseen:
+ collection = collections.BrokenBarHCollection.span_where(tictocs,
+ ymin=i+0.05,
+ ymax=i+0.95,
+ where=colours == colour,
+ facecolor=colour,
+ linewidths=0)
+ ax.add_collection(collection)
# Legend and room for it.
@@ -247,7 +268,8 @@ for rank in range(nranks):
nrow = nrow + 1
ax.fill_between([0, 0], nthread+0.5, nthread + nrow + 0.5, facecolor="white")
ax.set_ylim(0, nthread + nrow + 1)
- ax.legend(loc=1, shadow=True, mode="expand", ncol=5)
+ if data.size > 0:
+ ax.legend(loc=1, shadow=True, mode="expand", ncol=5)
# Start and end of time-step
ax.plot([0, 0], [0, nthread + nrow + 1], 'k--', linewidth=1)
diff --git a/src/Makefile.am b/src/Makefile.am
index 97db6eca672974a62561388cda3a80ba659475b6..fb606d7bd36cc9593958cb385caef3318960967e 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -42,19 +42,20 @@ 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 sourceterms.h
+ physical_constants.h physical_constants_cgs.h potential.h version.h \
+ hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.h
+
# Common source files
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 \
- runner_doiact_fft.c threadpool.c cooling.c sourceterms.c
+ physical_constants.c potential.c hydro_properties.c \
+ runner_doiact_fft.c threadpool.c cooling.c sourceterms.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 \
@@ -69,11 +70,17 @@ nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_h
hydro/Default/hydro_debug.h hydro/Default/hydro_part.h \
hydro/Gadget2/hydro.h hydro/Gadget2/hydro_iact.h hydro/Gadget2/hydro_io.h \
hydro/Gadget2/hydro_debug.h hydro/Gadget2/hydro_part.h \
+ hydro/PressureEntropy/hydro.h hydro/PressureEntropy/hydro_iact.h hydro/PressureEntropy/hydro_io.h \
+ hydro/PressureEntropy/hydro_debug.h hydro/PressureEntropy/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/adiabatic_index.h b/src/adiabatic_index.h
index a0c9ce09e3e004af07e8b208ef9f1af5f46c9e81..59937db3c8d09fc7e6e969de0aee60f01a2e5a2c 100644
--- a/src/adiabatic_index.h
+++ b/src/adiabatic_index.h
@@ -410,7 +410,7 @@ __attribute__((always_inline)) INLINE static float pow_one_over_gamma(float x) {
#if defined(HYDRO_GAMMA_5_3)
- return powf(x, 0.6f); /* x^(3/5) */
+ return powf(x, hydro_one_over_gamma); /* x^(3/5) */
#elif defined(HYDRO_GAMMA_7_5)
@@ -418,7 +418,7 @@ __attribute__((always_inline)) INLINE static float pow_one_over_gamma(float x) {
#elif defined(HYDRO_GAMMA_4_3)
- return powf(x, 0.75f); /* x^(3/4) */
+ return powf(x, hydro_one_over_gamma); /* x^(3/4) */
#elif defined(HYDRO_GAMMA_2_1)
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.c b/src/cell.c
index 7ce6fb81a8fa6875884d3f5c840c36e5177cdf6b..8728f32cb04f2f1de389386b632371a27d523bb9 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -679,6 +679,59 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset) {
part_relink_parts(gparts, gcount, parts - parts_offset);
}
+/**
+ * @brief Sanitizes the smoothing length values of cells by setting large
+ * outliers to more sensible values.
+ *
+ * We compute the mean and standard deviation of the smoothing lengths in
+ * logarithmic space and limit values to mean + 4 sigma.
+ *
+ * @param c The cell.
+ */
+void cell_sanitize(struct cell *c) {
+
+ const int count = c->count;
+ struct part *parts = c->parts;
+
+ /* First collect some statistics */
+ float h_mean = 0.f, h_mean2 = 0.f;
+ float h_min = FLT_MAX, h_max = 0.f;
+ for (int i = 0; i < count; ++i) {
+
+ const float h = logf(parts[i].h);
+ h_mean += h;
+ h_mean2 += h * h;
+ h_max = max(h_max, h);
+ h_min = min(h_min, h);
+ }
+ h_mean /= count;
+ h_mean2 /= count;
+ const float h_var = h_mean2 - h_mean * h_mean;
+ const float h_std = (h_var > 0.f) ? sqrtf(h_var) : 0.1f * h_mean;
+
+ /* Choose a cut */
+ const float h_limit = expf(h_mean + 4.f * h_std);
+
+ /* Be verbose this is not innocuous */
+ message("Cell properties: h_min= %f h_max= %f geometric mean= %f.",
+ expf(h_min), expf(h_max), expf(h_mean));
+
+ if (c->h_max > h_limit) {
+
+ message("Smoothing lengths will be limited to (mean + 4sigma)= %f.",
+ h_limit);
+
+ /* Apply the cut */
+ for (int i = 0; i < count; ++i) parts->h = min(parts[i].h, h_limit);
+
+ c->h_max = h_limit;
+
+ } else {
+
+ message("Smoothing lengths will not be limited.");
+ }
+}
+
/**
* @brief Converts hydro quantities to a valid state after the initial density
* calculation
diff --git a/src/cell.h b/src/cell.h
index b1abd0e2134687c123ba0a2e7b397df5572349ff..f87420c86c5765f42fe04e8dee95558edff90cf9 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 {
@@ -185,7 +187,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;
@@ -211,7 +213,7 @@ struct cell {
#endif
-} __attribute__((aligned(cell_align)));
+} SWIFT_STRUCT_ALIGN;
/* Convert cell location to ID. */
#define cell_getid(cdim, i, j, k) \
@@ -219,6 +221,7 @@ struct cell {
/* Function prototypes. */
void cell_split(struct cell *c, ptrdiff_t parts_offset);
+void cell_sanitize(struct cell *c);
int cell_locktree(struct cell *c);
void cell_unlocktree(struct cell *c);
int cell_glocktree(struct cell *c);
diff --git a/src/const.h b/src/const.h
index 817c645ba2354614caa6d70bf42d54bca41a284f..bffdaa7f773a160178017f9ba8a77c796028bf55 100644
--- a/src/const.h
+++ b/src/const.h
@@ -65,6 +65,7 @@
/* SPH variant to use */
//#define MINIMAL_SPH
#define GADGET2_SPH
+//#define HOPKINS_PE_SPH
//#define DEFAULT_SPH
//#define GIZMO_SPH
@@ -90,16 +91,18 @@
#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
/* Source terms */
#define SN_FEEDBACK
/* 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/debug.c b/src/debug.c
index d73bc86a92cf5ca28c202e7b567cf7c40ba6eccb..be42485a38ea8d560797e8f1ccc5936456febcd8 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -44,6 +44,8 @@
#include "./hydro/Minimal/hydro_debug.h"
#elif defined(GADGET2_SPH)
#include "./hydro/Gadget2/hydro_debug.h"
+#elif defined(HOPKINS_PE_SPH)
+#include "./hydro/PressureEntropy/hydro_debug.h"
#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_debug.h"
#elif defined(GIZMO_SPH)
@@ -65,7 +67,7 @@
*
* (Should be used for debugging only as it runs in O(N).)
*/
-void printParticle(const struct part *parts, struct xpart *xparts,
+void printParticle(const struct part *parts, const struct xpart *xparts,
long long int id, size_t N) {
int found = 0;
@@ -125,7 +127,7 @@ void printgParticle(const struct gpart *gparts, const struct part *parts,
*/
void printParticle_single(const struct part *p, const struct xpart *xp) {
- printf("## Particle: id=%lld", p->id);
+ printf("## Particle: id=%lld ", p->id);
hydro_debug_particle(p, xp);
printf("\n");
}
diff --git a/src/debug.h b/src/debug.h
index 22b63820745ca7282b7699f0be09e493238d83c2..2142a22eca91338580d8f50197a57de0cf248bee 100644
--- a/src/debug.h
+++ b/src/debug.h
@@ -24,7 +24,7 @@
#include "part.h"
#include "space.h"
-void printParticle(const struct part *parts, struct xpart *xparts,
+void printParticle(const struct part *parts, const struct xpart *xparts,
long long int id, size_t N);
void printgParticle(const struct gpart *gparts, const struct part *parts,
long long int id, size_t N);
diff --git a/src/engine.c b/src/engine.c
index f877bf5f28335b3c937fc653951d64286333bb68..175eaec6de5fe47665f0f6276143ab840ed2f75a 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1946,14 +1946,14 @@ void engine_maketasks(struct engine *e) {
scheduler_ranktasks(sched);
/* Weight the tasks. */
- scheduler_reweight(sched);
+ scheduler_reweight(sched, e->verbose);
/* Set the tasks age. */
e->tasks_age = 0;
if (e->verbose)
- message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
- clocks_getunit());
+ message("took %.3f %s (including reweight).",
+ clocks_from_ticks(getticks() - tic), clocks_getunit());
}
/**
@@ -2254,6 +2254,8 @@ void engine_rebuild(struct engine *e) {
/* Re-build the space. */
space_rebuild(e->s, 0.0, e->verbose);
+ if (e->ti_current == 0) space_sanitize(e->s);
+
/* If in parallel, exchange the cell structure. */
#ifdef WITH_MPI
engine_exchange_cells(e);
@@ -2278,8 +2280,11 @@ void engine_rebuild(struct engine *e) {
* @brief Prepare the #engine by re-building the cells and tasks.
*
* @param e The #engine to prepare.
+ * @param nodrift Whether to drift particles before rebuilding or not. Will
+ * not be necessary if all particles have already been
+ * drifted (before repartitioning for instance).
*/
-void engine_prepare(struct engine *e) {
+void engine_prepare(struct engine *e, int nodrift) {
TIMER_TIC;
@@ -2295,32 +2300,32 @@ void engine_prepare(struct engine *e) {
rebuild = buff;
#endif
- /* Did this not go through? */
+ /* And rebuild if necessary. */
if (rebuild) {
- /* First drift all particles to the current time */
- e->drift_all = 1;
- threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells_top,
- e->s->nr_cells, sizeof(struct cell), 1, e);
+ /* Drift all particles to the current time if needed. */
+ if (!nodrift) {
+ e->drift_all = 1;
+ engine_drift(e);
- /* Restore the default drifting policy */
- e->drift_all = (e->policy & engine_policy_drift_all);
+ /* Restore the default drifting policy */
+ e->drift_all = (e->policy & engine_policy_drift_all);
+ }
- /* And now rebuild */
engine_rebuild(e);
}
/* Re-rank the tasks every now and then. */
if (e->tasks_age % engine_tasksreweight == 1) {
- scheduler_reweight(&e->sched);
+ scheduler_reweight(&e->sched, e->verbose);
}
e->tasks_age += 1;
TIMER_TOC(timer_prepare);
if (e->verbose)
- message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
- clocks_getunit());
+ message("took %.3f %s (including marktask, rebuild and reweight).",
+ clocks_from_ticks(getticks() - tic), clocks_getunit());
}
/**
@@ -2412,6 +2417,7 @@ void engine_collect_kick(struct cell *c) {
*/
void engine_collect_timestep(struct engine *e) {
+ const ticks tic = getticks();
int updates = 0, g_updates = 0;
int ti_end_min = max_nr_timesteps;
const struct space *s = e->s;
@@ -2456,6 +2462,10 @@ void engine_collect_timestep(struct engine *e) {
e->ti_end_min = ti_end_min;
e->updates = updates;
e->g_updates = g_updates;
+
+ if (e->verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
}
/**
@@ -2465,9 +2475,11 @@ void engine_collect_timestep(struct engine *e) {
*/
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. */
@@ -2478,6 +2490,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];
@@ -2490,33 +2503,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
@@ -2524,13 +2539,17 @@ 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);
}
+
+ if (e->verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
}
/**
@@ -2544,6 +2563,8 @@ void engine_print_stats(struct engine *e) {
void engine_launch(struct engine *e, int nr_runners, unsigned int mask,
unsigned int submask) {
+ const ticks tic = getticks();
+
/* Prepare the scheduler. */
atomic_inc(&e->sched.waiting);
@@ -2568,6 +2589,10 @@ void engine_launch(struct engine *e, int nr_runners, unsigned int mask,
while (e->barrier_launch || e->barrier_running)
if (pthread_cond_wait(&e->barrier_cond, &e->barrier_mutex) != 0)
error("Error while waiting for barrier.");
+
+ if (e->verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
}
/**
@@ -2587,7 +2612,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
if (e->nodeID == 0) message("Running initialisation fake time-step.");
- engine_prepare(e);
+ engine_prepare(e, 1);
engine_marktasks(e);
@@ -2646,7 +2671,15 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
TIMER_TOC(timer_runners);
/* Apply some conversions (e.g. internal energy -> entropy) */
- if (!flag_entropy_ICs) space_map_cells_pre(s, 0, cell_convert_hydro, NULL);
+ if (!flag_entropy_ICs) {
+
+ /* Apply the conversion */
+ space_map_cells_pre(s, 0, cell_convert_hydro, NULL);
+
+ /* Correct what we did (e.g. in PE-SPH, need to recompute rho_bar) */
+ if (hydro_need_extra_init_loop)
+ engine_launch(e, e->nr_threads, mask, submask);
+ }
clocks_gettime(&time2);
@@ -2688,8 +2721,7 @@ void engine_step(struct engine *e) {
/* Drift everybody to the snapshot position */
e->drift_all = 1;
- threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells_top,
- e->s->nr_cells, sizeof(struct cell), 1, e);
+ engine_drift(e);
/* Restore the default drifting policy */
e->drift_all = (e->policy & engine_policy_drift_all);
@@ -2727,15 +2759,20 @@ void engine_step(struct engine *e) {
e->timeLastStatistics += e->deltaTimeStatistics;
}
- /* Drift only the necessary particles */
- threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells_top,
- e->s->nr_cells, sizeof(struct cell), 1, e);
+ /* Drift only the necessary particles, that all means all particles
+ * if we are about to repartition. */
+ int repart = (e->forcerepart != REPART_NONE);
+ e->drift_all = repart || e->drift_all;
+ engine_drift(e);
/* Re-distribute the particles amongst the nodes? */
- if (e->forcerepart != REPART_NONE) engine_repartition(e);
+ if (repart) engine_repartition(e);
/* Prepare the space. */
- engine_prepare(e);
+ engine_prepare(e, e->drift_all);
+
+ /* Restore the default drifting policy */
+ e->drift_all = (e->policy & engine_policy_drift_all);
/* Build the masks corresponding to the policy */
unsigned int mask = 0, submask = 0;
@@ -2807,6 +2844,8 @@ void engine_step(struct engine *e) {
submask |= 1 << task_subtype_tend;
}
+ if (e->verbose) engine_print_task_counts(e);
+
/* Send off the runners. */
TIMER_TIC;
engine_launch(e, e->nr_threads, mask, submask);
@@ -2827,6 +2866,21 @@ int engine_is_done(struct engine *e) {
return !(e->ti_current < max_nr_timesteps);
}
+/**
+ * @brief Drift particles using the current engine drift policy.
+ *
+ * @param e The #engine.
+ */
+void engine_drift(struct engine *e) {
+
+ const ticks tic = getticks();
+ threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells_top,
+ e->s->nr_cells, sizeof(struct cell), 1, e);
+ if (e->verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
+}
+
/**
* @brief Create and fill the proxies.
*
@@ -3339,11 +3393,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 1005dbb46aeefd6843e2fab278d0c17115cf869a..914025c125b9646be93b463e76606bf5813366f4 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 "sourceterms.h"
@@ -219,6 +219,7 @@ struct engine {
/* Function prototypes. */
void engine_barrier(struct engine *e, int tid);
void engine_compute_next_snapshot_time(struct engine *e);
+void engine_drift(struct engine *e);
void engine_dump_snapshot(struct engine *e);
void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, int nr_nodes, int nodeID,
@@ -231,7 +232,7 @@ void engine_init(struct engine *e, struct space *s,
struct sourceterms *sourceterms);
void engine_launch(struct engine *e, int nr_runners, unsigned int mask,
unsigned int submask);
-void engine_prepare(struct engine *e);
+void engine_prepare(struct engine *e, int nodrift);
void engine_print(struct engine *e);
void engine_init_particles(struct engine *e, int flag_entropy_ICs);
void engine_step(struct engine *e);
diff --git a/src/equation_of_state.h b/src/equation_of_state.h
index af59d8a2cad1632c67b6d377b5ed9dfe9484b4aa..5b19f99ab85d8a2b4e5e6f0b4eeb542925b4ee50 100644
--- a/src/equation_of_state.h
+++ b/src/equation_of_state.h
@@ -126,6 +126,21 @@ gas_soundspeed_from_internal_energy(float density, float u) {
return sqrtf(u * hydro_gamma * hydro_gamma_minus_one);
}
+/**
+ * @brief Returns the sound speed given density and pressure
+ *
+ * Computes \f$c = \sqrt{\frac{\gamma P}{\rho} }\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param P The pressure \f$P\f$
+ */
+__attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
+ float density, float P) {
+
+ const float density_inv = 1.f / density;
+ return sqrtf(hydro_gamma * P * density_inv);
+}
+
/* ------------------------------------------------------------------------- */
#elif defined(EOS_ISOTHERMAL_GAS)
@@ -221,6 +236,22 @@ gas_soundspeed_from_internal_energy(float density, float u) {
hydro_gamma_minus_one);
}
+/**
+ * @brief Returns the sound speed given density and pressure
+ *
+ * Since we are using an isothermal EoS, the pressure value is ignored
+ * Computes \f$c = \sqrt{u_{cst} \gamma \rho^{\gamma-1}}\f$.
+ *
+ * @param density The density \f$\rho\f$
+ * @param P The pressure \f$P\f$
+ */
+__attribute__((always_inline)) INLINE static float gas_soundspeed_from_pressure(
+ float density, float P) {
+
+ return sqrtf(const_isothermal_internal_energy * hydro_gamma *
+ hydro_gamma_minus_one);
+}
+
/* ------------------------------------------------------------------------- */
#else
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.h b/src/hydro.h
index 4a2b0bd029d494e2091b9081d22b7949cec5648c..9e02c2009e307f0623ffb535ff9068603c2d4147 100644
--- a/src/hydro.h
+++ b/src/hydro.h
@@ -34,6 +34,10 @@
#include "./hydro/Gadget2/hydro.h"
#include "./hydro/Gadget2/hydro_iact.h"
#define SPH_IMPLEMENTATION "Gadget-2 version of SPH (Springel 2005)"
+#elif defined(HOPKINS_PE_SPH)
+#include "./hydro/PressureEntropy/hydro.h"
+#include "./hydro/PressureEntropy/hydro_iact.h"
+#define SPH_IMPLEMENTATION "Pressure-Entropy SPH (Hopkins 2013)"
#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro.h"
#include "./hydro/Default/hydro_iact.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..0bdf56fe844309998db1fad4cf9581edb6b88305 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -19,12 +19,25 @@
#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"
#include "equation_of_state.h"
#include "hydro_properties.h"
#include "kernel_hydro.h"
+#include "minmax.h"
/**
* @brief Returns the internal energy of a particle
@@ -104,8 +117,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
* @brief Modifies the thermal state of a particle to the imposed internal
* energy
*
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * This overwrites the current state of the particle but does *not* change its
+ * time-derivatives. Entropy, pressure, sound-speed and signal velocity will be
+ * updated.
*
* @param p The particle
* @param u The new internal energy
@@ -114,13 +128,31 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
struct part *restrict p, float u) {
p->entropy = gas_entropy_from_internal_energy(p->rho, u);
+
+ /* Compute the new pressure */
+ const float pressure = gas_pressure_from_internal_energy(p->rho, u);
+
+ /* Compute the new sound speed */
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+
+ /* Update the signal velocity */
+ const float v_sig_old = p->force.v_sig;
+ const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
+ const float v_sig = max(v_sig_old, v_sig_new);
+
+ const float rho_inv = 1.f / p->rho;
+
+ p->force.soundspeed = soundspeed;
+ p->force.P_over_rho2 = pressure * rho_inv * rho_inv;
+ p->force.v_sig = v_sig;
}
/**
* @brief Modifies the thermal state of a particle to the imposed entropy
*
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * This overwrites the current state of the particle but does *not* change its
+ * time-derivatives. Entropy, pressure, sound-speed and signal velocity will be
+ * updated.
*
* @param p The particle
* @param S The new entropy
@@ -129,6 +161,23 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
struct part *restrict p, float S) {
p->entropy = S;
+
+ /* Compute the pressure */
+ const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+
+ /* Compute the new sound speed */
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+
+ /* Update the signal velocity */
+ const float v_sig_old = p->force.v_sig;
+ const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
+ const float v_sig = max(v_sig_old, v_sig_new);
+
+ const float rho_inv = 1.f / p->rho;
+
+ p->force.soundspeed = soundspeed;
+ p->force.P_over_rho2 = pressure * rho_inv * rho_inv;
+ p->force.v_sig = v_sig;
}
/**
@@ -183,7 +232,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
- p->rho_dh = 0.f;
+ p->density.rho_dh = 0.f;
p->density.div_v = 0.f;
p->density.rot_v[0] = 0.f;
p->density.rot_v[1] = 0.f;
@@ -210,27 +259,24 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* Final operation on the density (add self-contribution). */
p->rho += p->mass * kernel_root;
- p->rho_dh -= hydro_dimension * p->mass * kernel_root;
+ p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
p->density.wcount += kernel_root;
/* Finish the calculation by inserting the missing h-factors */
p->rho *= h_inv_dim;
- p->rho_dh *= h_inv_dim_plus_one;
+ p->density.rho_dh *= h_inv_dim_plus_one;
p->density.wcount *= kernel_norm;
p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
- const float irho = 1.f / p->rho;
-
- /* Compute the derivative term */
- p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * irho);
+ const float rho_inv = 1.f / p->rho;
/* Finish calculation of the velocity curl components */
- p->density.rot_v[0] *= h_inv_dim_plus_one * irho;
- p->density.rot_v[1] *= h_inv_dim_plus_one * irho;
- p->density.rot_v[2] *= h_inv_dim_plus_one * irho;
+ p->density.rot_v[0] *= h_inv_dim_plus_one * rho_inv;
+ p->density.rot_v[1] *= h_inv_dim_plus_one * rho_inv;
+ p->density.rot_v[2] *= h_inv_dim_plus_one * rho_inv;
/* Finish calculation of the velocity divergence */
- p->density.div_v *= h_inv_dim_plus_one * irho;
+ p->density.div_v *= h_inv_dim_plus_one * rho_inv;
}
/**
@@ -261,19 +307,23 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
const float pressure = hydro_get_pressure(p, half_dt);
- const float irho = 1.f / p->rho;
-
- /* Divide the pressure by the density and density gradient */
- const float P_over_rho2 = pressure * irho * irho * p->rho_dh;
-
/* Compute the sound speed */
- const float soundspeed = sqrtf(hydro_gamma * pressure * irho);
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+
+ /* Divide the pressure by the density squared to get the SPH term */
+ const float rho_inv = 1.f / p->rho;
+ const float P_over_rho2 = pressure * rho_inv * rho_inv;
/* Compute the Balsara switch */
const float balsara =
abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed / fac_mu / p->h);
+ /* Compute the "grad h" term */
+ const float grad_h_term =
+ 1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
+
/* Update variables. */
+ p->force.f = grad_h_term;
p->force.P_over_rho2 = P_over_rho2;
p->force.soundspeed = soundspeed;
p->force.balsara = balsara;
@@ -337,17 +387,16 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
const float pressure = hydro_get_pressure(p, dt_entr);
- const float irho = 1.f / p->rho;
-
- /* Divide the pressure by the density and density gradient */
- const float P_over_rho2 = pressure * irho * irho * p->rho_dh;
-
/* Compute the new sound speed */
- const float soundspeed = sqrtf(hydro_gamma * pressure * irho);
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+
+ /* Divide the pressure by the density squared to get the SPH term */
+ const float rho_inv = 1.f / p->rho;
+ const float P_over_rho2 = pressure * rho_inv * rho_inv;
/* Update variables */
- p->force.P_over_rho2 = P_over_rho2;
p->force.soundspeed = soundspeed;
+ p->force.P_over_rho2 = P_over_rho2;
}
/**
@@ -388,6 +437,19 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
/* Do not 'overcool' when timestep increases */
if (p->entropy + p->entropy_dt * half_dt < 0.5f * p->entropy)
p->entropy_dt = -0.5f * p->entropy / half_dt;
+
+ /* Compute the pressure */
+ const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+
+ /* Compute the new sound speed */
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+
+ /* Divide the pressure by the density squared to get the SPH term */
+ const float rho_inv = 1.f / p->rho;
+ const float P_over_rho2 = pressure * rho_inv * rho_inv;
+
+ p->force.soundspeed = soundspeed;
+ p->force.P_over_rho2 = P_over_rho2;
}
/**
@@ -402,6 +464,19 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
/* We read u in the entropy field. We now get S from u */
p->entropy = gas_entropy_from_internal_energy(p->rho, p->entropy);
+
+ /* Compute the pressure */
+ const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
+
+ /* Compute the sound speed */
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+
+ /* Divide the pressure by the density squared to get the SPH term */
+ const float rho_inv = 1.f / p->rho;
+ const float P_over_rho2 = pressure * rho_inv * rho_inv;
+
+ p->force.soundspeed = soundspeed;
+ p->force.P_over_rho2 = P_over_rho2;
}
#endif /* SWIFT_GADGET2_HYDRO_H */
diff --git a/src/hydro/Gadget2/hydro_debug.h b/src/hydro/Gadget2/hydro_debug.h
index 7c8a6eaba96929b01f1901393d7d0498d58badf4..656299b38374f68824ec20d85ece169d5f1fd599 100644
--- a/src/hydro/Gadget2/hydro_debug.h
+++ b/src/hydro/Gadget2/hydro_debug.h
@@ -30,8 +30,8 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
"v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%d\n",
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
- p->h, p->density.wcount, p->density.wcount_dh, p->mass, p->rho_dh, p->rho,
- hydro_get_pressure(p, 0.), p->force.P_over_rho2, p->entropy,
+ p->h, p->density.wcount, p->density.wcount_dh, p->mass, p->density.rho_dh,
+ p->rho, hydro_get_pressure(p, 0.), p->force.P_over_rho2, p->entropy,
p->entropy_dt, p->force.soundspeed, p->density.div_v, p->density.rot_v[0],
p->density.rot_v[1], p->density.rot_v[2], p->force.balsara,
p->force.v_sig, p->force.h_dt, p->ti_begin, p->ti_end);
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 8a4edfe62f59a3fae551fdb65f46987509f89251..fca1bcff91a71fb2a26adc117382630a576f9090 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
*/
@@ -60,7 +59,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* Compute contribution to the density */
pi->rho += mj * wi;
- pi->rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+ pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
/* Compute contribution to the number of neighbours */
pi->density.wcount += wi;
@@ -73,7 +72,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
/* Compute contribution to the density */
pj->rho += mi * wj;
- pj->rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
+ pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
/* Compute contribution to the number of neighbours */
pj->density.wcount += wj;
@@ -210,13 +209,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_density(
/* Update particles. */
for (k = 0; k < VEC_SIZE; k++) {
pi[k]->rho += rhoi.f[k];
- pi[k]->rho_dh -= rhoi_dh.f[k];
+ pi[k]->density.rho_dh -= rhoi_dh.f[k];
pi[k]->density.wcount += wcounti.f[k];
pi[k]->density.wcount_dh -= wcounti_dh.f[k];
pi[k]->density.div_v -= div_vi.f[k];
for (j = 0; j < 3; j++) pi[k]->density.rot_v[j] += curl_vi[j].f[k];
pj[k]->rho += rhoj.f[k];
- pj[k]->rho_dh -= rhoj_dh.f[k];
+ pj[k]->density.rho_dh -= rhoj_dh.f[k];
pj[k]->density.wcount += wcountj.f[k];
pj[k]->density.wcount_dh -= wcountj_dh.f[k];
pj[k]->density.div_v -= div_vj.f[k];
@@ -249,17 +248,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
const float ri = 1.0f / r;
/* Compute the kernel function */
- const float h_inv = 1.0f / hi;
- const float u = r * h_inv;
- kernel_deval(u, &wi, &wi_dx);
+ const float hi_inv = 1.0f / hi;
+ const float ui = r * hi_inv;
+ kernel_deval(ui, &wi, &wi_dx);
/* Compute contribution to the density */
pi->rho += mj * wi;
- pi->rho_dh -= mj * (hydro_dimension * wi + u * wi_dx);
+ pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
/* Compute contribution to the number of neighbours */
pi->density.wcount += wi;
- pi->density.wcount_dh -= u * wi_dx;
+ pi->density.wcount_dh -= ui * wi_dx;
const float fac = mj * wi_dx * ri;
@@ -367,7 +366,7 @@ runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
/* Update particles. */
for (k = 0; k < VEC_SIZE; k++) {
pi[k]->rho += rhoi.f[k];
- pi[k]->rho_dh -= rhoi_dh.f[k];
+ pi[k]->density.rho_dh -= rhoi_dh.f[k];
pi[k]->density.wcount += wcounti.f[k];
pi[k]->density.wcount_dh -= wcounti_dh.f[k];
pi[k]->density.div_v -= div_vi.f[k];
@@ -416,7 +415,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
kernel_deval(xj, &wj, &wj_dx);
const float wj_dr = hjd_inv * wj_dx;
- /* Compute gradient terms */
+ /* Compute h-gradient terms */
+ const float f_i = pi->force.f;
+ const float f_j = pj->force.f;
+
+ /* Compute pressure terms */
const float P_over_rho2_i = pi->force.P_over_rho2;
const float P_over_rho2_j = pj->force.P_over_rho2;
@@ -448,7 +451,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Now, convolve with the kernel */
const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
const float sph_term =
- (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
+ (f_i * P_over_rho2_i * wi_dr + f_j * P_over_rho2_j * wj_dr) * r_inv;
/* Eventually got the acceleration */
const float acc = visc_term + sph_term;
@@ -691,7 +694,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
kernel_deval(xj, &wj, &wj_dx);
const float wj_dr = hjd_inv * wj_dx;
- /* Compute gradient terms */
+ /* Compute h-gradient terms */
+ const float f_i = pi->force.f;
+ const float f_j = pj->force.f;
+
+ /* Compute pressure terms */
const float P_over_rho2_i = pi->force.P_over_rho2;
const float P_over_rho2_j = pj->force.P_over_rho2;
@@ -723,7 +730,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Now, convolve with the kernel */
const float visc_term = 0.5f * visc * (wi_dr + wj_dr) * r_inv;
const float sph_term =
- (P_over_rho2_i * wi_dr + P_over_rho2_j * wj_dr) * r_inv;
+ (f_i * P_over_rho2_i * wi_dr + f_j * P_over_rho2_j * wj_dr) * r_inv;
/* Eventually got the acceleration */
const float acc = visc_term + sph_term;
diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h
index 484792438d2717413c1ca8d4f429eac2e6d21b20..0c14da957463720dcee5349e88be91986cbb3f54 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 {
@@ -63,9 +80,6 @@ struct part {
/* Entropy time derivative */
float entropy_dt;
- /* Derivative of the density with respect to smoothing length. */
- float rho_dh;
-
union {
struct {
@@ -76,6 +90,9 @@ struct part {
/* Number of neighbours spatial derivative. */
float wcount_dh;
+ /* Derivative of the density with respect to h. */
+ float rho_dh;
+
/* Particle velocity curl. */
float rot_v[3];
@@ -89,6 +106,9 @@ struct part {
/* Balsara switch */
float balsara;
+ /*! "Grad h" term */
+ float f;
+
/* Pressure over density squared (including drho/dh term) */
float P_over_rho2;
@@ -110,6 +130,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.h b/src/hydro/Gizmo/hydro.h
index 37ce34e5910a033fab82dbc69839888a28c5ab12..a383a2fdce9e64a5673d0136184368220dc08458 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -150,18 +150,16 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* compute primitive variables */
/* eqns (3)-(5) */
const float m = p->conserved.mass;
- if (m > 0.f) {
- float momentum[3];
- momentum[0] = p->conserved.momentum[0];
- momentum[1] = p->conserved.momentum[1];
- momentum[2] = p->conserved.momentum[2];
- p->primitives.rho = m / volume;
- p->primitives.v[0] = momentum[0] / m;
- p->primitives.v[1] = momentum[1] / m;
- p->primitives.v[2] = momentum[2] / m;
- const float energy = p->conserved.energy;
- p->primitives.P = hydro_gamma_minus_one * energy / volume;
- }
+ float momentum[3];
+ momentum[0] = p->conserved.momentum[0];
+ momentum[1] = p->conserved.momentum[1];
+ momentum[2] = p->conserved.momentum[2];
+ p->primitives.rho = m / volume;
+ p->primitives.v[0] = momentum[0] / m;
+ p->primitives.v[1] = momentum[1] / m;
+ p->primitives.v[2] = momentum[2] / m;
+ const float energy = p->conserved.energy;
+ p->primitives.P = hydro_gamma_minus_one * energy / volume;
}
/**
@@ -259,6 +257,12 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
const float m = p->conserved.mass;
p->primitives.rho = m / volume;
+ /* first get the initial velocities, as they were overwritten in end_density
+ */
+ p->primitives.v[0] = p->v[0];
+ p->primitives.v[1] = p->v[1];
+ p->primitives.v[2] = p->v[2];
+
p->conserved.momentum[0] = m * p->primitives.v[0];
p->conserved.momentum[1] = m * p->primitives.v[1];
p->conserved.momentum[2] = m * p->primitives.v[2];
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..3015f26c6bad7f006e8cda16695c750ff40d74df 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"
@@ -40,6 +39,7 @@
#include "equation_of_state.h"
#include "hydro_properties.h"
#include "kernel_hydro.h"
+#include "minmax.h"
/**
* @brief Returns the internal energy of a particle
@@ -66,9 +66,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
__attribute__((always_inline)) INLINE static float hydro_get_pressure(
const struct part *restrict p, float dt) {
- const float u = p->u + p->u_dt * dt;
-
- return gas_pressure_from_internal_energy(p->rho, u);
+ return p->force.pressure;
}
/**
@@ -98,9 +96,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
const struct part *restrict p, float dt) {
- const float u = p->u + p->u_dt * dt;
-
- return gas_soundspeed_from_internal_energy(p->rho, u);
+ return p->force.soundspeed;
}
/**
@@ -129,8 +125,9 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
* @brief Modifies the thermal state of a particle to the imposed internal
* energy
*
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * This overwrites the current state of the particle but does *not* change its
+ * time-derivatives. Internal energy, pressure, sound-speed and signal velocity
+ * will be updated.
*
* @param p The particle
* @param u The new internal energy
@@ -139,13 +136,29 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
struct part *restrict p, float u) {
p->u = u;
+
+ /* Compute the new pressure */
+ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+
+ /* Compute the new sound speed */
+ const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);
+
+ /* Update the signal velocity */
+ const float v_sig_old = p->force.v_sig;
+ const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
+ const float v_sig = max(v_sig_old, v_sig_new);
+
+ p->force.soundspeed = soundspeed;
+ p->force.pressure = pressure;
+ p->force.v_sig = v_sig;
}
/**
* @brief Modifies the thermal state of a particle to the imposed entropy
*
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
+ * This overwrites the current state of the particle but does *not* change its
+ * time-derivatives. Internal energy, pressure, sound-speed and signal velocity
+ * will be updated.
*
* @param p The particle
* @param S The new entropy
@@ -154,6 +167,21 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
struct part *restrict p, float S) {
p->u = gas_internal_energy_from_entropy(p->rho, S);
+
+ /* Compute the pressure */
+ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+
+ /* Compute the new sound speed */
+ const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);
+
+ /* Update the signal velocity */
+ const float v_sig_old = p->force.v_sig;
+ const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
+ const float v_sig = max(v_sig_old, v_sig_new);
+
+ p->force.soundspeed = soundspeed;
+ p->force.pressure = pressure;
+ p->force.v_sig = v_sig;
}
/**
@@ -215,7 +243,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->density.wcount = 0.f;
p->density.wcount_dh = 0.f;
p->rho = 0.f;
- p->rho_dh = 0.f;
+ p->density.rho_dh = 0.f;
}
/**
@@ -241,19 +269,14 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* Final operation on the density (add self-contribution). */
p->rho += p->mass * kernel_root;
- p->rho_dh -= hydro_dimension * p->mass * kernel_root;
+ p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
p->density.wcount += kernel_root;
/* Finish the calculation by inserting the missing h-factors */
p->rho *= h_inv_dim;
- p->rho_dh *= h_inv_dim_plus_one;
+ p->density.rho_dh *= h_inv_dim_plus_one;
p->density.wcount *= kernel_norm;
p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
-
- const float irho = 1.f / p->rho;
-
- /* Compute the derivative term */
- p->rho_dh = 1.f / (1.f + hydro_dimension_inv * p->h * p->rho_dh * irho);
}
/**
@@ -275,10 +298,22 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part *restrict p, struct xpart *restrict xp, int ti_current,
double timeBase) {
+ /* Compute the pressure */
const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
const float pressure = hydro_get_pressure(p, half_dt);
+ /* Compute the sound speed */
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+
+ /* Compute the "grad h" term */
+ const float rho_inv = 1.f / p->rho;
+ const float grad_h_term =
+ 1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
+
+ /* Update variables. */
+ p->force.f = grad_h_term;
p->force.pressure = pressure;
+ p->force.soundspeed = soundspeed;
}
/**
@@ -338,7 +373,13 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Drift the pressure */
const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
- p->force.pressure = hydro_get_pressure(p, dt_entr);
+ const float pressure = hydro_get_pressure(p, dt_entr);
+
+ /* Compute the new sound speed */
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
+
+ p->force.pressure = pressure;
+ p->force.soundspeed = soundspeed;
}
/**
@@ -380,6 +421,15 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
/* Do not 'overcool' when timestep increases */
if (p->u + p->u_dt * half_dt < 0.5f * p->u) p->u_dt = -0.5f * p->u / half_dt;
+
+ /* Compute the pressure */
+ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+
+ /* Compute the sound speed */
+ const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);
+
+ p->force.pressure = pressure;
+ p->force.soundspeed = soundspeed;
}
/**
@@ -393,6 +443,16 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
* @param p The particle to act upon
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
- struct part *restrict p) {}
+ struct part *restrict p) {
+
+ /* Compute the pressure */
+ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+
+ /* Compute the sound speed */
+ const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u);
+
+ p->force.pressure = pressure;
+ p->force.soundspeed = soundspeed;
+}
#endif /* SWIFT_MINIMAL_HYDRO_H */
diff --git a/src/hydro/Minimal/hydro_debug.h b/src/hydro/Minimal/hydro_debug.h
index 0a6e50da6051fd961f4bf35f32caee311888a13e..16ae62413a0d76b7bf871e615fe5684219752fee 100644
--- a/src/hydro/Minimal/hydro_debug.h
+++ b/src/hydro/Minimal/hydro_debug.h
@@ -44,7 +44,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
p->u, p->u_dt, p->force.v_sig, p->force.pressure, p->h, p->force.h_dt,
- (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->ti_begin,
+ (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho, p->ti_begin,
p->ti_end);
}
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index 9e2028c978dc7cad03cfba17931f645bbfbfe1a0..169947b99e92d9bd1b0870d502a49e311820ff81 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"
@@ -56,7 +55,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
kernel_deval(xi, &wi, &wi_dx);
pi->rho += mj * wi;
- pi->rho_dh -= mj * (hydro_dimension * wi + xi * wi_dx);
+ pi->density.rho_dh -= mj * (hydro_dimension * wi + xi * wi_dx);
pi->density.wcount += wi;
pi->density.wcount_dh -= xi * wi_dx;
@@ -66,7 +65,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
kernel_deval(xj, &wj, &wj_dx);
pj->rho += mi * wj;
- pj->rho_dh -= mi * (hydro_dimension * wj + xj * wj_dx);
+ pj->density.rho_dh -= mi * (hydro_dimension * wj + xj * wj_dx);
pj->density.wcount += wj;
pj->density.wcount_dh -= xj * wj_dx;
}
@@ -101,7 +100,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
kernel_deval(xi, &wi, &wi_dx);
pi->rho += mj * wi;
- pi->rho_dh -= mj * (hydro_dimension * wi + xi * wi_dx);
+ pi->density.rho_dh -= mj * (hydro_dimension * wi + xi * wi_dx);
pi->density.wcount += wi;
pi->density.wcount_dh -= xi * wi_dx;
}
@@ -153,8 +152,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */
- const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
- const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
+ const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f;
+ const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f;
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
@@ -166,8 +165,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Compute sound speeds and signal velocity */
- const float ci = sqrtf(hydro_gamma * pressurei / rhoi);
- const float cj = sqrtf(hydro_gamma * pressurej / rhoj);
+ const float ci = pi->force.soundspeed;
+ const float cj = pj->force.soundspeed;
const float v_sig = ci + cj - 3.f * mu_ij;
/* Construct the full viscosity term */
@@ -264,8 +263,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float wj_dr = hjd_inv * wj_dx;
/* Compute gradient terms */
- const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->rho_dh;
- const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
+ const float P_over_rho2_i = pressurei / (rhoi * rhoi) * pi->force.f;
+ const float P_over_rho2_j = pressurej / (rhoj * rhoj) * pj->force.f;
/* Compute dv dot r. */
const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
@@ -277,8 +276,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Compute sound speeds and signal velocity */
- const float ci = sqrtf(hydro_gamma * pressurei / rhoi);
- const float cj = sqrtf(hydro_gamma * pressurej / rhoj);
+ const float ci = pi->force.soundspeed;
+ const float cj = pj->force.soundspeed;
const float v_sig = ci + cj - 3.f * mu_ij;
/* Construct the full viscosity term */
diff --git a/src/hydro/Minimal/hydro_part.h b/src/hydro/Minimal/hydro_part.h
index ad65f8b44fc67f4aae6470246cbab91bc3710007..8542177278998d5e0b830dc164988611549ef24d 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,35 @@ 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. */
-
- float rho_dh; /*!< Derivative of density with respect to h */
+ /*! Particle density. */
+ float rho;
/* Store density/force specific stuff. */
union {
@@ -92,10 +105,15 @@ 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;
+
+ /*! Derivative of density with respect to h */
+ float rho_dh;
- float wcount_dh; /*!< Derivative of the neighbour number with respect to
- h. */
} density;
/**
@@ -107,19 +125,30 @@ struct part {
*/
struct {
- float pressure; /*!< Particle pressure. */
+ /*! "Grad h" term */
+ float f;
+
+ /*! Particle pressure. */
+ float pressure;
+
+ /*! Particle soundspeed. */
+ float soundspeed;
- 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/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
new file mode 100644
index 0000000000000000000000000000000000000000..ede16fe5181a546921db7acf8ffcbc130c5e91c5
--- /dev/null
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -0,0 +1,515 @@
+/*******************************************************************************
+ * 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_PRESSURE_ENTROPY_HYDRO_H
+#define SWIFT_PRESSURE_ENTROPY_HYDRO_H
+
+/**
+ * @file PressureEntropy/hydro.h
+ * @brief Pressure-Entropy implementation of SPH (Non-neighbour loop
+ * equations)
+ *
+ * The thermal variable is the entropy (S) and the entropy is smoothed over
+ * contact discontinuities to prevent spurious surface tension.
+ *
+ * Follows eqautions (19), (21) and (22) of Hopkins, P., MNRAS, 2013,
+ * Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
+ */
+
+#include "adiabatic_index.h"
+#include "approx_math.h"
+#include "dimension.h"
+#include "equation_of_state.h"
+#include "hydro_properties.h"
+#include "kernel_hydro.h"
+#include "minmax.h"
+
+/**
+ * @brief Returns the internal energy of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
+ const struct part *restrict p, float dt) {
+
+ const float entropy = p->entropy + p->entropy_dt * dt;
+
+ return gas_internal_energy_from_entropy(p->rho_bar, entropy);
+}
+
+/**
+ * @brief Returns the pressure of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_pressure(
+ const struct part *restrict p, float dt) {
+
+ const float entropy = p->entropy + p->entropy_dt * dt;
+
+ return gas_pressure_from_entropy(p->rho_bar, entropy);
+}
+
+/**
+ * @brief Returns the entropy of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_entropy(
+ const struct part *restrict p, float dt) {
+
+ return p->entropy + p->entropy_dt * dt;
+}
+
+/**
+ * @brief Returns the sound speed of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
+ const struct part *restrict p, float dt) {
+
+ return p->force.soundspeed;
+}
+
+/**
+ * @brief Returns the physical density of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_density(
+ const struct part *restrict p) {
+
+ return p->rho;
+}
+
+/**
+ * @brief Returns the mass of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_mass(
+ const struct part *restrict p) {
+
+ return p->mass;
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed internal
+ * energy
+ *
+ * This overwrites the current state of the particle but does *not* change its
+ * time-derivatives. Entropy, pressure, sound-speed and signal velocity will be
+ * updated.
+ *
+ * @param p The particle
+ * @param u The new internal energy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
+ struct part *restrict p, float u) {
+
+ p->entropy = gas_entropy_from_internal_energy(p->rho_bar, u);
+ p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
+
+ /* Compute the pressure */
+ const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
+
+ /* Compute the sound speed from the pressure*/
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
+
+ /* Update the signal velocity */
+ const float v_sig_old = p->force.v_sig;
+ const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
+ const float v_sig = max(v_sig_old, v_sig_new);
+
+ const float rho_bar_inv = 1.f / p->rho_bar;
+
+ p->force.soundspeed = soundspeed;
+ p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
+ p->force.v_sig = v_sig;
+}
+
+/**
+ * @brief Modifies the thermal state of a particle to the imposed entropy
+ *
+ * This overwrites the current state of the particle but does *not* change its
+ * time-derivatives. Entropy, pressure, sound-speed and signal velocity will be
+ * updated.
+ *
+ * @param p The particle
+ * @param S The new entropy
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_entropy(
+ struct part *restrict p, float S) {
+
+ p->entropy = S;
+ p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
+
+ /* Compute the pressure */
+ const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
+
+ /* Compute the sound speed from the pressure*/
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
+
+ /* Update the signal velocity */
+ const float v_sig_old = p->force.v_sig;
+ const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
+ const float v_sig = max(v_sig_old, v_sig_new);
+
+ const float rho_bar_inv = 1.f / p->rho_bar;
+
+ p->force.soundspeed = soundspeed;
+ p->force.P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
+ p->force.v_sig = v_sig;
+}
+
+/**
+ * @brief Computes the hydro time-step of a given particle
+ *
+ * @param p Pointer to the particle data
+ * @param xp Pointer to the extended particle data
+ *
+ */
+__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
+ const struct part *restrict p, const struct xpart *restrict xp,
+ const struct hydro_props *restrict hydro_properties) {
+
+ const float CFL_condition = hydro_properties->CFL_condition;
+
+ /* CFL condition */
+ const float dt_cfl =
+ 2.f * kernel_gamma * CFL_condition * p->h / p->force.v_sig;
+
+ return dt_cfl;
+}
+
+/**
+ * @brief Initialises the particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_first_init_part(
+ struct part *restrict p, struct xpart *restrict xp) {
+
+ p->ti_begin = 0;
+ p->ti_end = 0;
+ p->rho_bar = 0.f;
+ p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
+ xp->v_full[0] = p->v[0];
+ xp->v_full[1] = p->v[1];
+ xp->v_full[2] = p->v[2];
+}
+
+/**
+ * @brief Prepares a particle for the density calculation.
+ *
+ * Zeroes all the relevant arrays in preparation for the sums taking place in
+ * the variaous density tasks
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_init_part(
+ struct part *restrict p) {
+
+ p->rho = 0.f;
+ p->rho_bar = 0.f;
+ p->density.wcount = 0.f;
+ p->density.wcount_dh = 0.f;
+ p->density.rho_dh = 0.f;
+ p->density.pressure_dh = 0.f;
+
+ p->density.div_v = 0.f;
+ p->density.rot_v[0] = 0.f;
+ p->density.rot_v[1] = 0.f;
+ p->density.rot_v[2] = 0.f;
+}
+
+/**
+ * @brief Finishes the density calculation.
+ *
+ * Multiplies the density and number of neighbours by the appropiate constants
+ * and add the self-contribution term.
+ *
+ * @param p The particle to act upon
+ * @param ti_current The current time (on the integer timeline)
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_density(
+ struct part *restrict p, int ti_current) {
+
+ /* Some smoothing length multiples. */
+ const float h = p->h;
+ const float h_inv = 1.0f / h; /* 1/h */
+ const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+ const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
+
+ /* Final operation on the density (add self-contribution). */
+ p->rho += p->mass * kernel_root;
+ p->rho_bar += p->mass * p->entropy_one_over_gamma * kernel_root;
+ p->density.rho_dh -= hydro_dimension * p->mass * kernel_root;
+ p->density.pressure_dh -=
+ hydro_dimension * p->mass * p->entropy_one_over_gamma * kernel_root;
+ p->density.wcount += kernel_root;
+
+ /* Finish the calculation by inserting the missing h-factors */
+ p->rho *= h_inv_dim;
+ p->rho_bar *= h_inv_dim;
+ p->density.rho_dh *= h_inv_dim_plus_one;
+ p->density.pressure_dh *= h_inv_dim_plus_one;
+ p->density.wcount *= kernel_norm;
+ p->density.wcount_dh *= h_inv * kernel_gamma * kernel_norm;
+
+ const float rho_inv = 1.f / p->rho;
+ const float entropy_minus_one_over_gamma = 1.f / p->entropy_one_over_gamma;
+
+ /* Final operation on the weighted density */
+ p->rho_bar *= entropy_minus_one_over_gamma;
+
+ /* Finish calculation of the velocity curl components */
+ p->density.rot_v[0] *= h_inv_dim_plus_one * rho_inv;
+ p->density.rot_v[1] *= h_inv_dim_plus_one * rho_inv;
+ p->density.rot_v[2] *= h_inv_dim_plus_one * rho_inv;
+
+ /* Finish calculation of the velocity divergence */
+ p->density.div_v *= h_inv_dim_plus_one * rho_inv;
+}
+
+/**
+ * @brief Prepare a particle for the force calculation.
+ *
+ * Computes viscosity term, conduction term and smoothing length gradient terms.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param ti_current The current time (on the timeline)
+ * @param timeBase The minimal time-step size
+ */
+__attribute__((always_inline)) INLINE static void hydro_prepare_force(
+ struct part *restrict p, struct xpart *restrict xp, int ti_current,
+ double timeBase) {
+
+ const float fac_mu = 1.f; /* Will change with cosmological integration */
+
+ /* Compute the norm of the curl */
+ const float curl_v = sqrtf(p->density.rot_v[0] * p->density.rot_v[0] +
+ p->density.rot_v[1] * p->density.rot_v[1] +
+ p->density.rot_v[2] * p->density.rot_v[2]);
+
+ /* Compute the norm of div v */
+ const float abs_div_v = fabsf(p->density.div_v);
+
+ /* Compute the pressure */
+ const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
+ const float entropy = hydro_get_entropy(p, half_dt);
+ const float pressure = gas_pressure_from_entropy(p->rho_bar, entropy);
+
+ /* Compute the sound speed from the pressure*/
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
+
+ /* Compute the Balsara switch */
+ const float balsara =
+ abs_div_v / (abs_div_v + curl_v + 0.0001f * soundspeed / fac_mu / p->h);
+
+ /* Divide the pressure by the density squared to get the SPH term */
+ const float rho_bar_inv = 1.f / p->rho_bar;
+ const float P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
+
+ /* Compute "grad h" term (note we use rho here and not rho_bar !)*/
+ const float rho_inv = 1.f / p->rho;
+ const float entropy_minus_one_over_gamma = 1.f / p->entropy_one_over_gamma;
+ const float rho_dh =
+ 1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
+ const float pressure_dh = p->density.pressure_dh * rho_inv * p->h *
+ hydro_dimension_inv * entropy_minus_one_over_gamma;
+
+ const float grad_h_term = rho_dh * pressure_dh;
+
+ /* Update variables. */
+ p->force.soundspeed = soundspeed;
+ p->force.P_over_rho2 = P_over_rho2;
+ p->force.balsara = balsara;
+ p->force.f = grad_h_term;
+}
+
+/**
+ * @brief Reset acceleration fields of a particle
+ *
+ * Resets all hydro acceleration and time derivative fields in preparation
+ * for the sums taking place in the variaous force tasks
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
+ struct part *restrict p) {
+
+ /* Reset the acceleration. */
+ p->a_hydro[0] = 0.0f;
+ p->a_hydro[1] = 0.0f;
+ p->a_hydro[2] = 0.0f;
+
+ /* Reset the time derivatives. */
+ p->entropy_dt = 0.0f;
+ p->force.h_dt = 0.0f;
+
+ /* Reset maximal signal velocity */
+ p->force.v_sig = 0.0f;
+}
+
+/**
+ * @brief Predict additional particle fields forward in time when drifting
+ *
+ * @param p The particle
+ * @param xp The extended data of the particle
+ * @param dt The drift time-step.
+ * @param t0 The time at the start of the drift (on the timeline).
+ * @param t1 The time at the end of the drift (on the timeline).
+ * @param timeBase The minimal time-step size
+ */
+__attribute__((always_inline)) INLINE static void hydro_predict_extra(
+ struct part *restrict p, const struct xpart *restrict xp, float dt, int t0,
+ int t1, double timeBase) {
+
+ const float h_inv = 1.f / p->h;
+
+ /* Predict smoothing length */
+ const float w1 = p->force.h_dt * h_inv * dt;
+ if (fabsf(w1) < 0.2f)
+ p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
+ else
+ p->h *= expf(w1);
+
+ /* Predict density */
+ const float w2 = -hydro_dimension * w1;
+ if (fabsf(w2) < 0.2f) {
+ p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
+ p->rho_bar *= approx_expf(w2);
+ } else {
+ p->rho *= expf(w2);
+ p->rho_bar *= expf(w2);
+ }
+
+ /* Drift the entropy */
+ const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
+ const float entropy = hydro_get_entropy(p, dt_entr);
+
+ /* Compute the pressure */
+ const float pressure = gas_pressure_from_entropy(p->rho_bar, entropy);
+
+ /* Compute the new sound speed */
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
+
+ /* Divide the pressure by the density squared to get the SPH term */
+ const float rho_bar_inv = 1.f / p->rho_bar;
+ const float P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
+
+ /* Update the variables */
+ p->entropy_one_over_gamma = pow_one_over_gamma(entropy);
+ p->force.soundspeed = soundspeed;
+ p->force.P_over_rho2 = P_over_rho2;
+}
+
+/**
+ * @brief Finishes the force calculation.
+ *
+ * Multiplies the forces and accelerationsby the appropiate constants
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_force(
+ struct part *restrict p) {
+
+ p->force.h_dt *= p->h * hydro_dimension_inv;
+
+ p->entropy_dt *=
+ 0.5f * hydro_gamma_minus_one * pow_minus_gamma_minus_one(p->rho_bar);
+}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * @param p The particle to act upon
+ * @param xp The particle extended data to act upon
+ * @param dt The time-step for this kick
+ * @param half_dt The half time-step for this kick
+ */
+__attribute__((always_inline)) INLINE static void hydro_kick_extra(
+ struct part *restrict p, struct xpart *restrict xp, float dt,
+ float half_dt) {
+
+ /* Do not decrease the entropy (temperature) by more than a factor of 2*/
+ const float entropy_change = p->entropy_dt * dt;
+ if (entropy_change > -0.5f * p->entropy)
+ p->entropy += entropy_change;
+ else
+ p->entropy *= 0.5f;
+
+ /* Do not 'overcool' when timestep increases */
+ if (p->entropy + p->entropy_dt * half_dt < 0.5f * p->entropy)
+ p->entropy_dt = -0.5f * p->entropy / half_dt;
+
+ /* Compute the pressure */
+ const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
+
+ /* Compute the new sound speed */
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
+
+ /* Divide the pressure by the density squared to get the SPH term */
+ const float rho_bar_inv = 1.f / p->rho_bar;
+ const float P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
+
+ p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
+ p->force.soundspeed = soundspeed;
+ p->force.P_over_rho2 = P_over_rho2;
+}
+
+/**
+ * @brief Converts hydro quantity of a particle at the start of a run
+ *
+ * Requires the density to be known
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
+ struct part *restrict p) {
+
+ /* We read u in the entropy field. We now get S from u */
+ p->entropy = gas_entropy_from_internal_energy(p->rho_bar, p->entropy);
+ p->entropy_one_over_gamma = pow_one_over_gamma(p->entropy);
+
+ /* Compute the pressure */
+ const float pressure = gas_pressure_from_entropy(p->rho_bar, p->entropy);
+
+ /* Compute the sound speed */
+ const float soundspeed = gas_soundspeed_from_pressure(p->rho_bar, pressure);
+
+ /* Divide the pressure by the density squared to get the SPH term */
+ const float rho_bar_inv = 1.f / p->rho_bar;
+ const float P_over_rho2 = pressure * rho_bar_inv * rho_bar_inv;
+
+ p->force.soundspeed = soundspeed;
+ p->force.P_over_rho2 = P_over_rho2;
+}
+
+#endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_H */
diff --git a/src/hydro/PressureEntropy/hydro_debug.h b/src/hydro/PressureEntropy/hydro_debug.h
new file mode 100644
index 0000000000000000000000000000000000000000..486543793515795092e7cc97fe7b567b8230be3b
--- /dev/null
+++ b/src/hydro/PressureEntropy/hydro_debug.h
@@ -0,0 +1,50 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (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_PRESSURE_ENTROPY_HYDRO_DEBUG_H
+#define SWIFT_PRESSURE_ENTROPY_HYDRO_DEBUG_H
+
+/**
+ * @file PressureEntropy/hydro_debug.h
+ * @brief Pressure-Entropy implementation of SPH (Debugging routines)
+ *
+ * The thermal variable is the entropy (S) and the entropy is smoothed over
+ * contact discontinuities to prevent spurious surface tension.
+ *
+ * Follows eqautions (19), (21) and (22) of Hopkins, P., MNRAS, 2013,
+ * Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
+ */
+
+__attribute__((always_inline)) INLINE static void hydro_debug_particle(
+ const struct part* p, const struct xpart* xp) {
+ printf(
+ "x=[%.3e,%.3e,%.3e], "
+ "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
+ "h=%.3e, wcount=%.3f, wcount_dh=%.3e, m=%.3e, dh_drho=%.3e, rho=%.3e, "
+ "rho_bar=%.3e, P=%.3e, dP_dh=%.3e, P_over_rho2=%.3e, S=%.3e, S^1/g=%.3e, "
+ "dS/dt=%.3e,\nc=%.3e v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%d\n",
+ p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
+ xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
+ p->h, p->density.wcount, p->density.wcount_dh, p->mass, p->density.rho_dh,
+ p->rho, p->rho_bar, hydro_get_pressure(p, 0.), p->density.pressure_dh,
+ p->force.P_over_rho2, p->entropy, p->entropy_one_over_gamma,
+ p->entropy_dt, p->force.soundspeed, p->force.v_sig, p->force.h_dt,
+ p->ti_begin, p->ti_end);
+}
+
+#endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_DEBUG_H */
diff --git a/src/hydro/PressureEntropy/hydro_iact.h b/src/hydro/PressureEntropy/hydro_iact.h
new file mode 100644
index 0000000000000000000000000000000000000000..18e22233f6f53e488119a58a988ba4c9faf3db36
--- /dev/null
+++ b/src/hydro/PressureEntropy/hydro_iact.h
@@ -0,0 +1,399 @@
+/*******************************************************************************
+ * 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_PRESSURE_ENTROPY_HYDRO_IACT_H
+#define SWIFT_PRESSURE_ENTROPY_HYDRO_IACT_H
+
+/**
+ * @file PressureEntropy/hydro_iact.h
+ * @brief Pressure-Entropy implementation of SPH (Neighbour loop equations)
+ *
+ * The thermal variable is the entropy (S) and the entropy is smoothed over
+ * contact discontinuities to prevent spurious surface tension.
+ *
+ * Follows eqautions (19), (21) and (22) of Hopkins, P., MNRAS, 2013,
+ * Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
+ */
+
+/**
+ * @brief Density loop
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_density(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+ float wi, wi_dx;
+ float wj, wj_dx;
+ float dv[3], curlvr[3];
+
+ /* Get the masses. */
+ const float mi = pi->mass;
+ const float mj = pj->mass;
+
+ /* Get r and r inverse. */
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
+
+ /* Compute the kernel function for pi */
+ const float hi_inv = 1.f / hi;
+ const float ui = r * hi_inv;
+ kernel_deval(ui, &wi, &wi_dx);
+
+ /* Compute contribution to the density */
+ pi->rho += mj * wi;
+ pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+
+ /* Compute contribution to the number of neighbours */
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= ui * wi_dx;
+
+ /* Compute contribution to the weighted density */
+ pi->rho_bar += mj * pj->entropy_one_over_gamma * wi;
+ pi->density.pressure_dh -=
+ mj * pj->entropy_one_over_gamma * (hydro_dimension * wi + ui * wi_dx);
+
+ /* Compute the kernel function for pj */
+ const float hj_inv = 1.f / hj;
+ const float uj = r * hj_inv;
+ kernel_deval(uj, &wj, &wj_dx);
+
+ /* Compute contribution to the density */
+ pj->rho += mi * wj;
+ pj->density.rho_dh -= mi * (hydro_dimension * wj + uj * wj_dx);
+
+ /* Compute contribution to the number of neighbours */
+ pj->density.wcount += wj;
+ pj->density.wcount_dh -= uj * wj_dx;
+
+ /* Compute contribution to the weighted density */
+ pj->rho_bar += mi * pi->entropy_one_over_gamma * wj;
+ pj->density.pressure_dh -=
+ mi * pi->entropy_one_over_gamma * (hydro_dimension * wj + uj * wj_dx);
+
+ const float faci = mj * wi_dx * r_inv;
+ const float facj = mi * wj_dx * r_inv;
+
+ /* Compute dv dot r */
+ dv[0] = pi->v[0] - pj->v[0];
+ dv[1] = pi->v[1] - pj->v[1];
+ dv[2] = pi->v[2] - pj->v[2];
+ const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
+
+ pi->density.div_v -= faci * dvdr;
+ pj->density.div_v -= facj * dvdr;
+
+ /* Compute dv cross r */
+ curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
+ curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
+ curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
+
+ pi->density.rot_v[0] += faci * curlvr[0];
+ pi->density.rot_v[1] += faci * curlvr[1];
+ pi->density.rot_v[2] += faci * curlvr[2];
+
+ pj->density.rot_v[0] += facj * curlvr[0];
+ pj->density.rot_v[1] += facj * curlvr[1];
+ pj->density.rot_v[2] += facj * curlvr[2];
+}
+
+/**
+ * @brief Density loop (Vectorized version)
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_vec_density(
+ float *R2, float *Dx, float *Hi, float *Hj, struct part **pi,
+ struct part **pj) {
+
+ error("Vectorized version of Pressure-Entropy SPH routine not existant yet.");
+}
+
+/**
+ * @brief Density loop (non-symmetric version)
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+ float wi, wi_dx;
+ float dv[3], curlvr[3];
+
+ /* Get the masses. */
+ const float mj = pj->mass;
+
+ /* Get r and r inverse. */
+ const float r = sqrtf(r2);
+ const float ri = 1.0f / r;
+
+ /* Compute the kernel function */
+ const float h_inv = 1.0f / hi;
+ const float ui = r * h_inv;
+ kernel_deval(ui, &wi, &wi_dx);
+
+ /* Compute contribution to the density */
+ pi->rho += mj * wi;
+ pi->density.rho_dh -= mj * (hydro_dimension * wi + ui * wi_dx);
+
+ /* Compute contribution to the number of neighbours */
+ pi->density.wcount += wi;
+ pi->density.wcount_dh -= ui * wi_dx;
+
+ /* Compute contribution to the weighted density */
+ pi->rho_bar += mj * pj->entropy_one_over_gamma * wi;
+ pi->density.pressure_dh -=
+ mj * pj->entropy_one_over_gamma * (hydro_dimension * wi + ui * wi_dx);
+
+ const float fac = mj * wi_dx * ri;
+
+ /* Compute dv dot r */
+ dv[0] = pi->v[0] - pj->v[0];
+ dv[1] = pi->v[1] - pj->v[1];
+ dv[2] = pi->v[2] - pj->v[2];
+ const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
+ pi->density.div_v -= fac * dvdr;
+
+ /* Compute dv cross r */
+ curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
+ curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];
+ curlvr[2] = dv[0] * dx[1] - dv[1] * dx[0];
+
+ pi->density.rot_v[0] += fac * curlvr[0];
+ pi->density.rot_v[1] += fac * curlvr[1];
+ pi->density.rot_v[2] += fac * curlvr[2];
+}
+
+/**
+ * @brief Density loop (non-symmetric vectorized version)
+ */
+__attribute__((always_inline)) INLINE static void
+runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
+ struct part **pi, struct part **pj) {
+
+ error("Vectorized version of Pressure-Entropy SPH routine not existant yet.");
+}
+
+/**
+ * @brief Force loop
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_force(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+ float wi, wj, wi_dx, wj_dx;
+
+ const float fac_mu = 1.f; /* Will change with cosmological integration */
+
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
+
+ /* Get some values in local variables. */
+ const float mi = pi->mass;
+ const float mj = pj->mass;
+ const float rhoi = pi->rho;
+ const float rhoj = pj->rho;
+
+ /* Get the kernel for hi. */
+ const float hi_inv = 1.0f / hi;
+ const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
+ const float ui = r * hi_inv;
+ kernel_deval(ui, &wi, &wi_dx);
+ const float wi_dr = hid_inv * wi_dx;
+
+ /* Get the kernel for hj. */
+ const float hj_inv = 1.0f / hj;
+ const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
+ const float xj = r * hj_inv;
+ kernel_deval(xj, &wj, &wj_dx);
+ const float wj_dr = hjd_inv * wj_dx;
+
+ /* Compute gradient terms */
+ const float f_i = pi->force.f;
+ const float f_j = pj->force.f;
+
+ /* Compute Pressure terms */
+ const float P_over_rho2_i = pi->force.P_over_rho2;
+ const float P_over_rho2_j = pj->force.P_over_rho2;
+
+ /* Compute entropy terms */
+ const float S_gamma_i = pi->entropy_one_over_gamma;
+ const float S_gamma_j = pj->entropy_one_over_gamma;
+
+ /* Compute sound speeds */
+ const float ci = pi->force.soundspeed;
+ const float cj = pj->force.soundspeed;
+
+ /* Compute dv dot r. */
+ const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
+ (pi->v[1] - pj->v[1]) * dx[1] +
+ (pi->v[2] - pj->v[2]) * dx[2];
+
+ /* Balsara term */
+ const float balsara_i = pi->force.balsara;
+ const float balsara_j = pj->force.balsara;
+
+ /* Are the particles moving towards each others ? */
+ const float omega_ij = fminf(dvdr, 0.f);
+ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+ /* Signal velocity */
+ const float v_sig = ci + cj - 3.f * mu_ij;
+
+ /* Now construct the full viscosity term */
+ const float rho_ij = 0.5f * (rhoi + rhoj);
+ const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
+ (balsara_i + balsara_j) / rho_ij;
+
+ /* Now, convolve with the kernel */
+ const float visc_term = 0.5f * visc * (wi_dr + wj_dr);
+ const float sph_term = (S_gamma_j / S_gamma_i - f_i) * P_over_rho2_i * wi_dr +
+ (S_gamma_i / S_gamma_j - f_j) * P_over_rho2_j * wj_dr;
+
+ /* Eventually got the acceleration */
+ const float acc = (visc_term + sph_term) * r_inv;
+
+ /* Use the force Luke ! */
+ pi->a_hydro[0] -= mj * acc * dx[0];
+ pi->a_hydro[1] -= mj * acc * dx[1];
+ pi->a_hydro[2] -= mj * acc * dx[2];
+
+ pj->a_hydro[0] += mi * acc * dx[0];
+ pj->a_hydro[1] += mi * acc * dx[1];
+ pj->a_hydro[2] += mi * acc * dx[2];
+
+ /* Get the time derivative for h. */
+ pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
+ pj->force.h_dt -= mi * dvdr * r_inv / rhoi * wj_dr;
+
+ /* Update the signal velocity. */
+ pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
+ pj->force.v_sig = fmaxf(pj->force.v_sig, v_sig);
+
+ /* Change in entropy */
+ pi->entropy_dt += mj * visc_term * r_inv * dvdr;
+ pj->entropy_dt += mi * visc_term * r_inv * dvdr;
+}
+
+/**
+ * @brief Force loop (Vectorized version)
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_vec_force(
+ float *R2, float *Dx, float *Hi, float *Hj, struct part **pi,
+ struct part **pj) {
+
+ error("Vectorized version of Pressure-Entropy SPH routine not existant yet.");
+}
+
+/**
+ * @brief Force loop (non-symmetric version)
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
+ float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+ float wi, wj, wi_dx, wj_dx;
+
+ const float fac_mu = 1.f; /* Will change with cosmological integration */
+
+ const float r = sqrtf(r2);
+ const float r_inv = 1.0f / r;
+
+ /* Get some values in local variables. */
+ // const float mi = pi->mass;
+ const float mj = pj->mass;
+ const float rhoi = pi->rho;
+ const float rhoj = pj->rho;
+
+ /* Get the kernel for hi. */
+ const float hi_inv = 1.0f / hi;
+ const float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
+ const float ui = r * hi_inv;
+ kernel_deval(ui, &wi, &wi_dx);
+ const float wi_dr = hid_inv * wi_dx;
+
+ /* Get the kernel for hj. */
+ const float hj_inv = 1.0f / hj;
+ const float hjd_inv = pow_dimension_plus_one(hj_inv); /* 1/h^(d+1) */
+ const float xj = r * hj_inv;
+ kernel_deval(xj, &wj, &wj_dx);
+ const float wj_dr = hjd_inv * wj_dx;
+
+ /* Compute gradient terms */
+ const float f_i = pi->force.f;
+ const float f_j = pj->force.f;
+
+ /* Compute Pressure terms */
+ const float P_over_rho2_i = pi->force.P_over_rho2;
+ const float P_over_rho2_j = pj->force.P_over_rho2;
+
+ /* Compute entropy terms */
+ const float S_gamma_i = pi->entropy_one_over_gamma;
+ const float S_gamma_j = pj->entropy_one_over_gamma;
+
+ /* Compute sound speeds */
+ const float ci = pi->force.soundspeed;
+ const float cj = pj->force.soundspeed;
+
+ /* Compute dv dot r. */
+ const float dvdr = (pi->v[0] - pj->v[0]) * dx[0] +
+ (pi->v[1] - pj->v[1]) * dx[1] +
+ (pi->v[2] - pj->v[2]) * dx[2];
+
+ /* Balsara term */
+ const float balsara_i = pi->force.balsara;
+ const float balsara_j = pj->force.balsara;
+
+ /* Are the particles moving towards each others ? */
+ const float omega_ij = fminf(dvdr, 0.f);
+ const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
+
+ /* Signal velocity */
+ const float v_sig = ci + cj - 3.f * mu_ij;
+
+ /* Now construct the full viscosity term */
+ const float rho_ij = 0.5f * (rhoi + rhoj);
+ const float visc = -0.25f * const_viscosity_alpha * v_sig * mu_ij *
+ (balsara_i + balsara_j) / rho_ij;
+
+ /* Now, convolve with the kernel */
+ const float visc_term = 0.5f * visc * (wi_dr + wj_dr);
+ const float sph_term = (S_gamma_j / S_gamma_i - f_i) * P_over_rho2_i * wi_dr +
+ (S_gamma_i / S_gamma_j - f_j) * P_over_rho2_j * wj_dr;
+
+ /* Eventually got the acceleration */
+ const float acc = (visc_term + sph_term) * r_inv;
+
+ /* Use the force Luke ! */
+ pi->a_hydro[0] -= mj * acc * dx[0];
+ pi->a_hydro[1] -= mj * acc * dx[1];
+ pi->a_hydro[2] -= mj * acc * dx[2];
+
+ /* Get the time derivative for h. */
+ pi->force.h_dt -= mj * dvdr * r_inv / rhoj * wi_dr;
+
+ /* Update the signal velocity. */
+ pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
+
+ /* Change in entropy */
+ pi->entropy_dt += mj * visc_term * r_inv * dvdr;
+}
+
+/**
+ * @brief Force loop (Vectorized non-symmetric version)
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
+ float *R2, float *Dx, float *Hi, float *Hj, struct part **pi,
+ struct part **pj) {
+
+ error("Vectorized version of Pressure-Entropy SPH routine not existant yet.");
+}
+
+#endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_IACT_H */
diff --git a/src/hydro/PressureEntropy/hydro_io.h b/src/hydro/PressureEntropy/hydro_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..9914a656466f3f0d0a5eeb79b511706d7068ffc6
--- /dev/null
+++ b/src/hydro/PressureEntropy/hydro_io.h
@@ -0,0 +1,145 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (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_PRESSURE_ENTROPY_HYDRO_IO_H
+#define SWIFT_PRESSURE_ENTROPY_HYDRO_IO_H
+
+/**
+ * @file PressureEntropy/hydro_io.h
+ * @brief Pressure-Entropy implementation of SPH (i/o routines)
+ *
+ * The thermal variable is the entropy (S) and the entropy is smoothed over
+ * contact discontinuities to prevent spurious surface tension.
+ *
+ * Follows eqautions (19), (21) and (22) of Hopkins, P., MNRAS, 2013,
+ * Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
+ */
+
+#include "adiabatic_index.h"
+#include "hydro.h"
+#include "io_properties.h"
+#include "kernel_hydro.h"
+
+/**
+ * @brief Specifies which particle fields to read from a dataset
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to read.
+ * @param num_fields The number of i/o fields to read.
+ */
+void hydro_read_particles(struct part* parts, struct io_props* list,
+ int* num_fields) {
+
+ *num_fields = 8;
+
+ /* List what we want to read */
+ list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
+ UNIT_CONV_LENGTH, parts, x);
+ list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
+ UNIT_CONV_SPEED, parts, v);
+ list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
+ parts, mass);
+ list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
+ UNIT_CONV_LENGTH, parts, h);
+ list[4] =
+ io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
+ UNIT_CONV_ENTROPY_PER_UNIT_MASS, parts, entropy);
+ list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
+ UNIT_CONV_NO_UNITS, parts, id);
+ list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL,
+ UNIT_CONV_ACCELERATION, parts, a_hydro);
+ list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL,
+ UNIT_CONV_DENSITY, parts, rho);
+}
+
+float convert_u(struct engine* e, struct part* p) {
+
+ return hydro_get_internal_energy(p, 0);
+}
+
+float convert_P(struct engine* e, struct part* p) {
+
+ return hydro_get_pressure(p, 0);
+}
+
+/**
+ * @brief Specifies which particle fields to write to a dataset
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to write.
+ * @param num_fields The number of i/o fields to write.
+ */
+void hydro_write_particles(struct part* parts, struct io_props* list,
+ int* num_fields) {
+
+ *num_fields = 11;
+
+ /* List what we want to write */
+ list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
+ parts, x);
+ list[1] =
+ io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, v);
+ list[2] =
+ io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
+ list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
+ parts, h);
+ list[4] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
+ UNIT_CONV_ENERGY_PER_UNIT_MASS,
+ parts, entropy, convert_u);
+ list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
+ UNIT_CONV_NO_UNITS, parts, id);
+ list[6] = io_make_output_field("Acceleration", FLOAT, 3,
+ UNIT_CONV_ACCELERATION, parts, a_hydro);
+ list[7] =
+ io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
+ list[8] = io_make_output_field(
+ "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, parts, entropy);
+ list[9] = io_make_output_field_convert_part(
+ "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, rho, convert_P);
+ list[10] = io_make_output_field("WeightedDensity", FLOAT, 1,
+ UNIT_CONV_DENSITY, parts, rho_bar);
+}
+
+/**
+ * @brief Writes the current model of SPH to the file
+ * @param h_grpsph The HDF5 group in which to write
+ */
+void writeSPHflavour(hid_t h_grpsph) {
+
+ /* Viscosity and thermal conduction */
+ /* Nothing in this minimal model... */
+ writeAttribute_s(h_grpsph, "Thermal Conductivity Model", "No treatment");
+ writeAttribute_s(
+ h_grpsph, "Viscosity Model",
+ "as in Springel (2005), i.e. Monaghan (1992) with Balsara (1995) switch");
+ writeAttribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
+ writeAttribute_f(h_grpsph, "Viscosity beta", 3.f);
+
+ /* Time integration properties */
+ writeAttribute_f(h_grpsph, "Maximal Delta u change over dt",
+ const_max_u_change);
+}
+
+/**
+ * @brief Are we writing entropy in the internal energy field ?
+ *
+ * @return 1 if entropy is in 'internal energy', 0 otherwise.
+ */
+int writeEntropyFlag() { return 0; }
+
+#endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_IO_H */
diff --git a/src/hydro/PressureEntropy/hydro_part.h b/src/hydro/PressureEntropy/hydro_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..cac585ff79bae737f0e5c09860a38536cbf3a38c
--- /dev/null
+++ b/src/hydro/PressureEntropy/hydro_part.h
@@ -0,0 +1,143 @@
+/*******************************************************************************
+ * 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_PRESSURE_ENTROPY_HYDRO_PART_H
+#define SWIFT_PRESSURE_ENTROPY_HYDRO_PART_H
+
+/**
+ * @file PressureEntropy/hydro_part.h
+ * @brief Pressure-Entropy implementation of SPH (Particle definition)
+ *
+ * The thermal variable is the entropy (S) and the entropy is smoothed over
+ * contact discontinuities to prevent spurious surface tension.
+ *
+ * Follows eqautions (19), (21) and (22) of Hopkins, P., MNRAS, 2013,
+ * Volume 428, Issue 4, pp. 2840-2856 with a simple Balsara viscosity term.
+ */
+
+#include "cooling_struct.h"
+
+/* Extra particle data not needed during the SPH loops over neighbours. */
+struct xpart {
+
+ /*! Offset between current position and position at last tree rebuild. */
+ float x_diff[3];
+
+ /*! Velocity at the last full step. */
+ float v_full[3];
+
+ /*! Additional data used to record cooling information */
+ struct cooling_xpart_data cooling_data;
+
+} SWIFT_STRUCT_ALIGN;
+
+/* Data of a single particle. */
+struct part {
+
+ /*! Particle position. */
+ double x[3];
+
+ /*! Particle predicted velocity. */
+ float v[3];
+
+ /*! Particle acceleration. */
+ float a_hydro[3];
+
+ /*! Particle cutoff radius. */
+ float h;
+
+ /*! Particle mass. */
+ float mass;
+
+ /*! Particle time of beginning of time-step. */
+ int ti_begin;
+
+ /*! Particle time of end of time-step. */
+ int ti_end;
+
+ /*! Particle density. */
+ float rho;
+
+ /*! Particle weighted density */
+ float rho_bar;
+
+ /*! Particle entropy. */
+ float entropy;
+
+ /*! Entropy time derivative */
+ float entropy_dt;
+
+ /*! Particle entropy to the power 1/gamma. */
+ float entropy_one_over_gamma;
+
+ union {
+
+ struct {
+
+ /*! Number of neighbours. */
+ float wcount;
+
+ /*! Number of neighbours spatial derivative. */
+ float wcount_dh;
+
+ /*! Derivative of density with respect to h */
+ float rho_dh;
+
+ /*! Derivative of pressure with respect to h */
+ float pressure_dh;
+
+ /*! Particle velocity curl. */
+ float rot_v[3];
+
+ /*! Particle velocity divergence. */
+ float div_v;
+
+ } density;
+
+ struct {
+
+ /*! Balsara switch */
+ float balsara;
+
+ /*! "Grad h" term */
+ float f;
+
+ /*! Pressure term */
+ float P_over_rho2;
+
+ /*! Particle sound speed. */
+ float soundspeed;
+
+ /*! Signal velocity. */
+ float v_sig;
+
+ /*! Time derivative of the smoothing length */
+ float h_dt;
+
+ } force;
+ };
+
+ /*! Particle ID. */
+ long long id;
+
+ /*! Pointer to corresponding gravity part. */
+ struct gpart* gpart;
+
+} SWIFT_STRUCT_ALIGN;
+
+#endif /* SWIFT_PRESSURE_ENTROPY_HYDRO_PART_H */
diff --git a/src/hydro_io.h b/src/hydro_io.h
index f0a619b90b774574c434007b1c01a0e55e75e464..05ae94ade7b103ff1b584dc2447cbab40479d1fc 100644
--- a/src/hydro_io.h
+++ b/src/hydro_io.h
@@ -26,6 +26,8 @@
#include "./hydro/Minimal/hydro_io.h"
#elif defined(GADGET2_SPH)
#include "./hydro/Gadget2/hydro_io.h"
+#elif defined(HOPKINS_PE_SPH)
+#include "./hydro/PressureEntropy/hydro_io.h"
#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_io.h"
#elif defined(GIZMO_SPH)
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 0da34d4dad114db0920c8e8f3bb617295ff3da96..66c9203e39e56d520eeace8858b0c618b45e6a22 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -486,6 +486,11 @@ void read_ic_parallel(char* fileName, const struct UnitSystem* internal_units,
}
}
+ /* Convert the dimensions of the box */
+ for (int j = 0; j < 3; j++)
+ dim[j] *=
+ units_conversion_factor(ic_units, internal_units, UNIT_CONV_LENGTH);
+
/* Allocate memory to store SPH particles */
*Ngas = N[0];
if (posix_memalign((void*)parts, part_align, (*Ngas) * sizeof(struct part)) !=
diff --git a/src/part.h b/src/part.h
index ea895e6e0295d6a8b63309c7bd6855daa2cf7d64..6832e0c6c2e0f2c324d90629447cf6a5e809d6fb 100644
--- a/src/part.h
+++ b/src/part.h
@@ -31,22 +31,30 @@
#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)
#include "./hydro/Minimal/hydro_part.h"
+#define hydro_need_extra_init_loop 0
#elif defined(GADGET2_SPH)
#include "./hydro/Gadget2/hydro_part.h"
+#define hydro_need_extra_init_loop 0
+#elif defined(HOPKINS_PE_SPH)
+#include "./hydro/PressureEntropy/hydro_part.h"
+#define hydro_need_extra_init_loop 1
#elif defined(DEFAULT_SPH)
#include "./hydro/Default/hydro_part.h"
+#define hydro_need_extra_init_loop 0
#elif defined(GIZMO_SPH)
#include "./hydro/Gizmo/hydro_part.h"
+#define hydro_need_extra_init_loop 0
#define EXTRA_HYDRO_LOOP
#else
#error "Invalid choice of SPH variant"
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/riemann/riemann_hllc.h b/src/riemann/riemann_hllc.h
index b8b1239d7799221c98522c06631aba5cabe69183..1a4c3f5f8338f6504a8c5f1d9eab8451eb20246c 100644
--- a/src/riemann/riemann_hllc.h
+++ b/src/riemann/riemann_hllc.h
@@ -145,11 +145,13 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
we add the extra velocity flux due to the absolute motion of the fluid
similarly, we need to add the energy fluxes due to the absolute motion */
v2 = vij[0] * vij[0] + vij[1] * vij[1] + vij[2] * vij[2];
+ // order is important: we first use the momentum fluxes to update the energy
+ // flux and then de-boost the momentum fluxes!
+ totflux[4] += vij[0] * totflux[1] + vij[1] * totflux[2] +
+ vij[2] * totflux[3] + 0.5 * v2 * totflux[0];
totflux[1] += vij[0] * totflux[0];
totflux[2] += vij[1] * totflux[0];
totflux[3] += vij[2] * totflux[0];
- totflux[4] += vij[0] * totflux[1] + vij[1] * totflux[2] +
- vij[2] * totflux[3] + 0.5 * v2 * totflux[0];
}
#endif /* SWIFT_RIEMANN_HLLC_H */
diff --git a/src/runner.c b/src/runner.c
index ea8f87a6b35ac2367493c2f07a7c80a89beebaa9..b77e6978f0b4acbaa2cd6fe7ff3f4a2fc81a3400 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -279,7 +279,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);
}
}
@@ -297,9 +297,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;
@@ -322,13 +323,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);
}
}
@@ -859,7 +861,8 @@ static void runner_do_drift(struct cell *c, struct engine *e) {
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};
@@ -931,6 +934,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);
@@ -957,6 +961,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];
@@ -974,6 +979,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];
diff --git a/src/scheduler.c b/src/scheduler.c
index fa065879f7daeb60cfeeb67e52c64eb2036cf3cb..f98f8824c00186956e8abe5e5c2aa3aae73abf4e 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -141,6 +141,7 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
/* Get a handle on the cell involved. */
struct cell *ci = t->ci;
+ const double hi = ci->dmin;
/* Foreign task? */
if (ci->nodeID != s->nodeID) {
@@ -149,7 +150,7 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
}
/* Is this cell even split? */
- if (ci->split) {
+ if (ci->split && ci->h_max * kernel_gamma * space_stretch < hi / 2) {
/* Make a sub? */
if (scheduler_dosub &&
@@ -890,10 +891,11 @@ void scheduler_reset(struct scheduler *s, int size) {
* @brief Compute the task weights
*
* @param s The #scheduler.
+ * @param verbose Are we talkative ?
*/
+void scheduler_reweight(struct scheduler *s, int verbose) {
-void scheduler_reweight(struct scheduler *s) {
-
+ const ticks tic = getticks();
const int nr_tasks = s->nr_tasks;
int *tid = s->tasks_ind;
struct task *tasks = s->tasks;
@@ -902,11 +904,8 @@ void scheduler_reweight(struct scheduler *s) {
0.4025, 0.1897, 0.4025, 0.1897, 0.4025,
0.5788, 0.4025, 0.5788};
const float wscale = 0.001;
- // ticks tic;
- /* Run through the tasks backwards and set their waits and
- weights. */
- // tic = getticks();
+ /* Run through the tasks backwards and set their weights. */
for (int k = nr_tasks - 1; k >= 0; k--) {
struct task *t = &tasks[tid[k]];
t->weight = 0;
@@ -963,8 +962,10 @@ void scheduler_reweight(struct scheduler *s) {
break;
}
}
- // message( "weighting tasks took %.3f %s." ,
- // clocks_from_ticks( getticks() - tic ), clocks_getunit());
+
+ if (verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
/* int min = tasks[0].weight, max = tasks[0].weight;
for ( int k = 1 ; k < nr_tasks ; k++ )
diff --git a/src/scheduler.h b/src/scheduler.h
index 3e51197de148ee8f0bf3090551568715d3000e98..c4eb5e99447d623e5fb8e442efc1c254c00bfadd 100644
--- a/src/scheduler.h
+++ b/src/scheduler.h
@@ -116,7 +116,7 @@ void scheduler_start(struct scheduler *s, unsigned int mask,
unsigned int submask);
void scheduler_reset(struct scheduler *s, int nr_tasks);
void scheduler_ranktasks(struct scheduler *s);
-void scheduler_reweight(struct scheduler *s);
+void scheduler_reweight(struct scheduler *s, int verbose);
struct task *scheduler_addtask(struct scheduler *s, enum task_types type,
enum task_subtypes subtype, int flags, int wait,
struct cell *ci, struct cell *cj, int tight);
diff --git a/src/serial_io.c b/src/serial_io.c
index 6e26be1a33fbc2c74ae1b8f7af2b83db285c962e..b9ad0fbaa856a889d3f84bb42013282f3640fd5e 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -528,6 +528,11 @@ void read_ic_serial(char* fileName, const struct UnitSystem* internal_units,
H5Fclose(h_file);
}
+ /* Convert the dimensions of the box */
+ for (int j = 0; j < 3; j++)
+ dim[j] *=
+ units_conversion_factor(ic_units, internal_units, UNIT_CONV_LENGTH);
+
/* Now need to broadcast that information to all ranks. */
MPI_Bcast(flag_entropy, 1, MPI_INT, 0, comm);
MPI_Bcast(periodic, 1, MPI_INT, 0, comm);
diff --git a/src/single_io.c b/src/single_io.c
index 6cb7e830209b0d58919fe6f529f675b4c611a51d..ceeba4eb80a47c3feed7e898deb5e1fe7e427c0c 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -433,6 +433,11 @@ void read_ic_single(char* fileName, const struct UnitSystem* internal_units,
internal_units->UnitTemperature_in_cgs);
}
+ /* Convert the dimensions of the box */
+ for (int j = 0; j < 3; j++)
+ dim[j] *=
+ units_conversion_factor(ic_units, internal_units, UNIT_CONV_LENGTH);
+
/* Allocate memory to store SPH particles */
*Ngas = N[0];
if (posix_memalign((void*)parts, part_align, *Ngas * sizeof(struct part)) !=
diff --git a/src/space.c b/src/space.c
index cdd5958cbc515003f4a86a41c9a7075fa3b4364f..1d2733ab2d4356c19cc0e65a0b9441edcf0acb80 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"
@@ -57,6 +58,7 @@
int space_splitsize = space_splitsize_default;
int space_subsize = space_subsize_default;
int space_maxsize = space_maxsize_default;
+int space_maxcount = space_maxcount_default;
/* Map shift vector to sortlist. */
const int sortlistID[27] = {
@@ -227,7 +229,14 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
if (s->periodic && (cdim[0] < 3 || cdim[1] < 3 || cdim[2] < 3))
error(
"Must have at least 3 cells in each spatial dimension when periodicity "
- "is switched on.");
+ "is switched on.\nThis error is often caused by any of the "
+ "followings:\n"
+ " - too few particles to generate a sensible grid,\n"
+ " - the initial value of 'SPH:max_smoothing_length' is too large,\n"
+ " - the (minimal) time-step is too large leading to particles with "
+ "predicted smoothing lengths too large for the box size,\n"
+ " - particle with velocities so large that they move by more than two "
+ "box sizes per time-step.\n");
/* Check if we have enough cells for gravity. */
if (s->gravity && (cdim[0] < 8 || cdim[1] < 8 || cdim[2] < 8))
@@ -690,7 +699,7 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
/* At this point, we have the upper-level cells, old or new. Now make
sure that the parts in each cell are ok. */
- space_split(s, cells_top, verbose);
+ space_split(s, cells_top, s->nr_cells, verbose);
if (verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
@@ -703,14 +712,16 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
* This is done in parallel using threads in the #threadpool.
*
* @param s The #space.
- * @param cells The cell hierarchy
+ * @param cells The cell hierarchy.
+ * @param nr_cells The number of cells.
* @param verbose Are we talkative ?
*/
-void space_split(struct space *s, struct cell *cells, int verbose) {
+void space_split(struct space *s, struct cell *cells, int nr_cells,
+ int verbose) {
const ticks tic = getticks();
- threadpool_map(&s->e->threadpool, space_split_mapper, cells, s->nr_cells,
+ threadpool_map(&s->e->threadpool, space_split_mapper, cells, nr_cells,
sizeof(struct cell), 1, s);
if (verbose)
@@ -718,6 +729,29 @@ void space_split(struct space *s, struct cell *cells, int verbose) {
clocks_getunit());
}
+/**
+ * @brief Runs through the top-level cells and checks whether tasks associated
+ * with them can be split. If not, try to sanitize the cells.
+ *
+ * @param s The #space to act upon.
+ */
+void space_sanitize(struct space *s) {
+
+ for (int k = 0; k < s->nr_cells; k++) {
+
+ struct cell *c = &s->cells_top[k];
+ const double min_width = c->dmin;
+
+ /* Do we have a problem ? */
+ if (c->h_max * kernel_gamma * space_stretch > min_width * 0.5 &&
+ c->count > space_maxcount) {
+
+ /* Ok, clean-up the mess */
+ cell_sanitize(c);
+ }
+ }
+}
+
/**
* @brief Sort the particles and condensed particles according to the given
* indices.
@@ -1239,17 +1273,17 @@ void space_map_cells_pre(struct space *s, int full,
* too many particles.
*
* @param map_data Pointer towards the top-cells.
- * @param num_elements The number of cells to treat.
+ * @param num_cells The number of cells to treat.
* @param extra_data Pointers to the #space.
*/
-void space_split_mapper(void *map_data, int num_elements, void *extra_data) {
+void space_split_mapper(void *map_data, int num_cells, void *extra_data) {
/* Unpack the inputs. */
struct space *s = (struct space *)extra_data;
struct cell *restrict cells_top = (struct cell *)map_data;
struct engine *e = s->e;
- for (int ind = 0; ind < num_elements; ind++) {
+ for (int ind = 0; ind < num_cells; ind++) {
struct cell *c = &cells_top[ind];
@@ -1268,6 +1302,12 @@ void space_split_mapper(void *map_data, int num_elements, void *extra_data) {
atomic_cas(&s->maxdepth, maxdepth, c->depth);
}
+ /* If the depth is too large, we have a problem and should stop. */
+ if (maxdepth > space_cell_maxdepth) {
+ error("Exceeded maximum depth (%d) when splitting cells, aborting",
+ space_cell_maxdepth);
+ }
+
/* Split or let it be? */
if (count > space_splitsize || gcount > space_splitsize) {
@@ -1471,6 +1511,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
*
@@ -1547,6 +1604,8 @@ void space_init(struct space *s, const struct swift_params *params,
space_subsize_default);
space_splitsize = parser_get_opt_param_int(
params, "Scheduler:cell_split_size", space_splitsize_default);
+ space_maxcount = parser_get_opt_param_int(params, "Scheduler:cell_max_count",
+ space_maxcount_default);
if (verbose)
message("max_size set to %d, sub_size set to %d, split_size set to %d",
space_maxsize, space_subsize, space_splitsize);
@@ -1630,6 +1689,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/space.h b/src/space.h
index 72b17405f13766ad2ccc9d53712068f28172067b..85f32bd7580945de6d9713316b557c92667987c0 100644
--- a/src/space.h
+++ b/src/space.h
@@ -41,13 +41,18 @@
#define space_splitsize_default 400
#define space_maxsize_default 8000000
#define space_subsize_default 64000000
+#define space_maxcount_default 10000
#define space_stretch 1.10f
#define space_maxreldx 0.25f
+/* Maximum allowed depth of cell splits. */
+#define space_cell_maxdepth 52
+
/* Split size. */
extern int space_splitsize;
extern int space_maxsize;
extern int space_subsize;
+extern int space_maxcount;
/* Map shift vector to sortlist. */
extern const int sortlistID[27];
@@ -145,6 +150,7 @@ void space_init(struct space *s, const struct swift_params *params,
double dim[3], struct part *parts, struct gpart *gparts,
size_t Npart, size_t Ngpart, int periodic, int gravity,
int verbose, int dry_run);
+void space_sanitize(struct space *s);
void space_map_cells_pre(struct space *s, int full,
void (*fun)(struct cell *c, void *data), void *data);
void space_map_parts(struct space *s,
@@ -161,7 +167,8 @@ void space_gparts_sort_mapper(void *map_data, int num_elements,
void *extra_data);
void space_rebuild(struct space *s, double h_max, int verbose);
void space_recycle(struct space *s, struct cell *c);
-void space_split(struct space *s, struct cell *cells, int verbose);
+void space_split(struct space *s, struct cell *cells, int nr_cells,
+ int verbose);
void space_split_mapper(void *map_data, int num_elements, void *extra_data);
void space_do_parts_sort();
void space_do_gparts_sort();
diff --git a/src/swift.h b/src/swift.h
index 9a425dbac34e6e6ae6ecc874d5994eaab4c9836b..58745f3111fc3f29490d4591ee549907392cb87e 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 18197f9c5f0b239562e015c6b6b88de307444c1c..00068f45769b4ca606cc729bd5e89c13ae729eef 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;
@@ -169,6 +169,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;
}
/**
@@ -178,7 +182,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 d7d8b162f985b475e2087b90b4e2777513e670f9..bc4df3dc2a4cfee3c382e9f2059cba84f29299f7 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"
@@ -150,7 +151,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..d423efba3fda764310721586e14192afdee18a85 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -95,6 +95,8 @@ void set_energy_state(struct part *part, enum pressure_field press, float size,
#if defined(GADGET2_SPH)
part->entropy = pressure / pow_gamma(density);
+#elif defined(HOPKINS_PE_SPH)
+ part->entropy = pressure / pow_gamma(density);
#elif defined(DEFAULT_SPH)
part->u = pressure / (hydro_gamma_minus_one * density);
#elif defined(MINIMAL_SPH)
@@ -519,6 +521,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/test27cells.c b/tests/test27cells.c
index 1a1ab88748d922b3e7fbb30a73a10809dca10863..22619af53c39218da34d771fab6ed2d10993689c 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -104,11 +104,18 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
}
part->h = size * h / (float)n;
part->id = ++(*partId);
+
#ifdef GIZMO_SPH
part->conserved.mass = density * volume / count;
#else
part->mass = density * volume / count;
#endif
+
+#if defined(HOPKINS_PE_SPH)
+ part->entropy = 1.f;
+ part->entropy_one_over_gamma = 1.f;
+#endif
+
part->ti_begin = 0;
part->ti_end = 1;
++part;
@@ -193,16 +200,11 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
#if defined(GIZMO_SPH)
0.f,
#else
- main_cell->parts[pid].rho_dh,
+ main_cell->parts[pid].density.rho_dh,
#endif
main_cell->parts[pid].density.wcount,
main_cell->parts[pid].density.wcount_dh,
-#if defined(GADGET2_SPH)
- main_cell->parts[pid].density.div_v,
- main_cell->parts[pid].density.rot_v[0],
- main_cell->parts[pid].density.rot_v[1],
- main_cell->parts[pid].density.rot_v[2]
-#elif defined(DEFAULT_SPH)
+#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH)
main_cell->parts[pid].density.div_v,
main_cell->parts[pid].density.rot_v[0],
main_cell->parts[pid].density.rot_v[1],
@@ -235,13 +237,10 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
#if defined(GIZMO_SPH)
0.f,
#else
- main_cell->parts[pjd].rho_dh,
+ main_cell->parts[pjd].density.rho_dh,
#endif
cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh,
-#if defined(GADGET2_SPH)
- cj->parts[pjd].density.div_v, cj->parts[pjd].density.rot_v[0],
- cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2]
-#elif defined(DEFAULT_SPH)
+#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH)
cj->parts[pjd].density.div_v, cj->parts[pjd].density.rot_v[0],
cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2]
#else
diff --git a/tests/testInteractions.c b/tests/testInteractions.c
index 52ad0c54258848883a9025bbcd9d68133eddc4b9..d14c840ec77819bbef5750b897c72139f4d7b2b4 100644
--- a/tests/testInteractions.c
+++ b/tests/testInteractions.c
@@ -110,9 +110,9 @@ void dump_indv_particle_fields(char *fileName, struct part *p) {
"%13e %13e %13e %13e "
"%13e %13e %13e %10f\n",
p->id, p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2],
- p->a_hydro[0], p->a_hydro[1], p->a_hydro[2], p->rho, p->rho_dh,
- p->density.wcount, p->density.wcount_dh, p->force.h_dt,
- p->force.v_sig,
+ p->a_hydro[0], p->a_hydro[1], p->a_hydro[2], p->rho,
+ p->density.rho_dh, p->density.wcount, p->density.wcount_dh,
+ p->force.h_dt, p->force.v_sig,
#if defined(GADGET2_SPH)
p->density.div_v, p->density.rot_v[0], p->density.rot_v[1],
p->density.rot_v[2], p->entropy_dt
diff --git a/tests/testPair.c b/tests/testPair.c
index efa1e628c2d57bf7922be8affdd5436ebca2f9cf..8b272b866431db3bfe36239222cd87d669961ae7 100644
--- a/tests/testPair.c
+++ b/tests/testPair.c
@@ -139,10 +139,10 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
#if defined(GIZMO_SPH)
0.f,
#else
- cj->parts[pid].rho_dh,
+ ci->parts[pid].density.rho_dh,
#endif
ci->parts[pid].density.wcount, ci->parts[pid].density.wcount_dh,
-#if defined(GADGET2_SPH) || defined(DEFAULT_SPH)
+#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH)
ci->parts[pid].density.div_v, ci->parts[pid].density.rot_v[0],
ci->parts[pid].density.rot_v[1], ci->parts[pid].density.rot_v[2]
#else
@@ -163,10 +163,10 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
#if defined(GIZMO_SPH)
0.f,
#else
- cj->parts[pjd].rho_dh,
+ cj->parts[pjd].density.rho_dh,
#endif
cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh,
-#if defined(GADGET2_SPH) || defined(DEFAULT_SPH)
+#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH)
cj->parts[pjd].density.div_v, cj->parts[pjd].density.rot_v[0],
cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2]
#else
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;
diff --git a/tests/tolerance_27_normal.dat b/tests/tolerance_27_normal.dat
index 141ed3baeedd5dad7a2fda0730c2e9c828ae4b2e..71acaa89be231d02fc33e47c96a7bacf623bbf48 100644
--- a/tests/tolerance_27_normal.dat
+++ b/tests/tolerance_27_normal.dat
@@ -1,3 +1,3 @@
# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 2e-5 2e-3 2e-6 2e-6 2e-6 2e-6
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-5 1e-4 2e-5 2e-5 2e-5 2e-5
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 2e-6 1e-5 1e-4 2e-5 2e-5 2e-5 2e-5
diff --git a/tests/tolerance_27_perturbed.dat b/tests/tolerance_27_perturbed.dat
index f0417d845f83d171dacb6b66024cf9a5dc41c6f1..45293cbaa223b5887f3b0ce05cd9430d0db7440b 100644
--- a/tests/tolerance_27_perturbed.dat
+++ b/tests/tolerance_27_perturbed.dat
@@ -1,3 +1,3 @@
# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1.2e-6 2e-6 2.1e-5 2e-3 2.1e-6 2e-6 2e-6 2e-6
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-5 1e-4 2e-5 4e-4 4e-4 4e-4
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1.2e-6 1e-5 2.1e-5 2e-3 2.1e-6 2e-6 2e-6 2e-6
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 3e-3 1e-5 1e-4 2e-5 4e-4 4e-4 4e-4