diff --git a/configure.ac b/configure.ac
index 82382447fde7c411f61dbd62f7db388a6a8d9cf9..977a86a03d5e4c82139dfd56d3f1ee7afe3ea6e8 100644
--- a/configure.ac
+++ b/configure.ac
@@ -155,6 +155,21 @@ LT_INIT
AC_PROG_CC_C99
AC_C_INLINE
+# If debugging try to show inlined functions.
+if test "x$enable_debug" = "xyes"; then
+ # Show inlined functions.
+ if test "$ax_cv_c_compiler_vendor" = "gnu"; then
+ # Would like to use -gdwarf and let the compiler pick a good version
+ # but that doesn't always work.
+ AX_CHECK_COMPILE_FLAG([-gdwarf -fvar-tracking-assignments],
+ [inline_EXTRA_FLAGS="-gdwarf -fvar-tracking-assignments"],
+ [inline_EXTRA_FLAGS="-gdwarf-2 -fvar-tracking-assignments"])
+ CFLAGS="$CFLAGS $inline_EXTRA_FLAGS"
+ elif test "$ax_cv_c_compiler_vendor" = "intel"; then
+ CFLAGS="$CFLAGS -debug inline-debug-info"
+ fi
+fi
+
# Define HAVE_POSIX_MEMALIGN if it works.
AX_FUNC_POSIX_MEMALIGN
@@ -451,7 +466,7 @@ AC_LINK_IFELSE([AC_LANG_PROGRAM(
AC_MSG_RESULT($rtc_ok)
# Add warning flags by default, if these can be used. Option =error adds
-# -Werror to GCC, clang and Intel. Note do this last as compiler tests may
+# -Werror to GCC, clang and Intel. Note do this last as compiler tests may
# become errors, if that's an issue don't use CFLAGS for these, use an AC_SUBST().
AC_ARG_ENABLE([compiler-warnings],
[AS_HELP_STRING([--enable-compiler-warnings],
@@ -461,7 +476,7 @@ AC_ARG_ENABLE([compiler-warnings],
[enable_warn="error"]
)
if test "$enable_warn" != "no"; then
-
+
# AX_CFLAGS_WARN_ALL does not give good warning flags for the Intel compiler
# We will do this by hand instead and only default to the macro for unknown compilers
case "$ax_cv_c_compiler_vendor" in
@@ -475,7 +490,7 @@ if test "$enable_warn" != "no"; then
AX_CFLAGS_WARN_ALL
;;
esac
-
+
# Add a "choke on warning" flag if it exists
if test "$enable_warn" = "error"; then
case "$ax_cv_c_compiler_vendor" in
diff --git a/examples/DiscPatch/HydroStatic/dynamic.pro b/examples/DiscPatch/HydroStatic/dynamic.pro
index c02c65fe418e84cdd62978dbddcf5a641fa4c156..00ee3f7a8d2dc435be2093af959efd2c49903637 100644
--- a/examples/DiscPatch/HydroStatic/dynamic.pro
+++ b/examples/DiscPatch/HydroStatic/dynamic.pro
@@ -8,7 +8,8 @@ iplot = 1 ; if iplot = 1, make plot of E/Lz conservation, else, simply compare f
@physunits
indir = './'
-basefile = 'Disk-Patch-dynamic_'
+;basefile = 'Disc-Patch-dynamic_'
+basefile = 'Disc-Patch_'
; set properties of potential
uL = phys.pc ; unit of length
@@ -16,18 +17,27 @@ uM = phys.msun ; unit of mass
uV = 1d5 ; unit of velocity
; properties of patch
-surface_density = 10.
+surface_density = 100. ; surface density of all mass, which generates the gravitational potential
scale_height = 100.
-z_disk = 200.;
+z_disk = 200. ;
+fgas = 0.1 ; gas fraction
gamma = 5./3.
; derived units
constG = 10.^(alog10(phys.g)+alog10(uM)-2d0*alog10(uV)-alog10(uL)) ;
pcentre = [0.,0.,z_disk] * pc / uL
utherm = !pi * constG * surface_density * scale_height / (gamma-1.)
+temp = (utherm*uV^2)*phys.m_h/phys.kb
soundspeed = sqrt(gamma * (gamma-1.) * utherm)
t_dyn = sqrt(scale_height / (constG * surface_density))
-
+rho0 = fgas*(surface_density)/(2.*scale_height)
+print,' dynamical time = ',t_dyn,' = ',t_dyn*UL/uV/(1d6*phys.yr),' Myr'
+print,' thermal energy per unit mass = ',utherm
+print,' central density = ',rho0,' = ',rho0*uM/uL^3/m_h,' particles/cm^3'
+print,' central temperature = ',temp
+lambda = 2 * !pi * phys.G^1.5 * (scale_height*uL)^1.5 * (surface_density * uM/uL^2)^0.5 * phys.m_h^2 / (gamma-1) / fgas
+print,' lambda = ',lambda
+stop
;
infile = indir + basefile + '*'
spawn,'ls -1 '+infile,res
diff --git a/examples/DiscPatch/HydroStatic/makeIC.py b/examples/DiscPatch/HydroStatic/makeIC.py
index 48cc578658a9520477d40bf504e3eb7c3c8e5164..e2846d82a8cfa8bf08de83632b19ae2e7818f3c1 100644
--- a/examples/DiscPatch/HydroStatic/makeIC.py
+++ b/examples/DiscPatch/HydroStatic/makeIC.py
@@ -56,9 +56,10 @@ print "UnitVelocity_in_cgs: ", const_unit_velocity_in_cgs
# parameters of potential
-surface_density = 10.
+surface_density = 100. # surface density of all mass, which generates the gravitational potential
scale_height = 100.
gamma = 5./3.
+fgas = 0.1 # gas fraction
# derived units
const_unit_time_in_cgs = (const_unit_length_in_cgs / const_unit_velocity_in_cgs)
@@ -131,7 +132,7 @@ h = glass_h[0:numGas]
numGas = numpy.shape(pos)[0]
# compute furthe properties of ICs
-column_density = surface_density * numpy.tanh(boxSize/2./scale_height)
+column_density = fgas * surface_density * numpy.tanh(boxSize/2./scale_height)
enclosed_mass = column_density * boxSize * boxSize
pmass = enclosed_mass / numGas
meanrho = enclosed_mass / boxSize**3
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/Feedback/feedback.pro b/examples/Feedback/feedback.pro
new file mode 100644
index 0000000000000000000000000000000000000000..02d616fc82f0aeb7011d022d13db9d1d1030e89c
--- /dev/null
+++ b/examples/Feedback/feedback.pro
@@ -0,0 +1,24 @@
+base = 'Feedback'
+inf = 'Feedback_005.hdf5'
+
+blast = [5.650488e-01, 5.004371e-01, 5.010494e-01] ; location of blast
+pos = h5rd(inf,'PartType0/Coordinates')
+vel = h5rd(inf,'PartType0/Velocities')
+rho = h5rd(inf,'PartType0/Density')
+utherm = h5rd(inf,'PartType0/InternalEnergy')
+
+; shift to centre
+for ic=0,2 do pos[ic,*] = pos[ic,*] - blast[ic]
+
+;; distance from centre
+dist = fltarr(n_elements(rho))
+for ic=0,2 do dist = dist + pos[ic,*]^2
+dist = sqrt(dist)
+
+; radial velocity
+vr = fltarr(n_elements(rho))
+for ic=0,2 do vr = vr + pos[ic,*]*vel[ic,*]
+vr = vr / dist
+
+;
+end
diff --git a/examples/Feedback/feedback.yml b/examples/Feedback/feedback.yml
new file mode 100644
index 0000000000000000000000000000000000000000..9ff7eea325b43fe3c670d16bdbd4ccc66f33333e
--- /dev/null
+++ b/examples/Feedback/feedback.yml
@@ -0,0 +1,44 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1 # Grams
+ UnitLength_in_cgs: 1 # Centimeters
+ UnitVelocity_in_cgs: 1 # Centimeters per second
+ UnitCurrent_in_cgs: 1 # Amperes
+ UnitTemp_in_cgs: 1 # Kelvin
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 5e-2 # The end time of the simulation (in internal units).
+ dt_min: 1e-7 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-4 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: Feedback # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 1e-2 # Time difference between consecutive outputs (in internal units)
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-3 # Time between statistics output
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
+ max_smoothing_length: 0.1 # Maximal smoothing length allowed (in internal units).
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: ./Feedback.hdf5 # The file to read
+
+# Parameters for feedback
+
+SN:
+ time: 0.001 # time the SN explodes (internal units)
+ energy: 1.0 # energy of the explosion (internal units)
+ x: 0.5 # x-position of explostion (internal units)
+ y: 0.5 # y-position of explostion (internal units)
+ z: 0.5 # z-position of explostion (internal units)
diff --git a/examples/Feedback/makeIC.py b/examples/Feedback/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..bd1081a9c275616038f5fa4e3eb943c36cb4c3eb
--- /dev/null
+++ b/examples/Feedback/makeIC.py
@@ -0,0 +1,109 @@
+###############################################################################
+ # This file is part of SWIFT.
+ # Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ # Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ # 2016 Tom Theuns (tom.theuns@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/>.
+ #
+ ##############################################################################
+
+import h5py
+import sys
+from numpy import *
+
+# Generates a swift IC file containing a cartesian distribution of particles
+# at a constant density and pressure in a cubic box
+
+# Parameters
+periodic= 1 # 1 For periodic box
+boxSize = 1.
+L = int(sys.argv[1]) # Number of particles along one axis
+rho = 1. # Density
+P = 1.e-6 # Pressure
+gamma = 5./3. # Gas adiabatic index
+eta = 1.2349 # 48 ngbs with cubic spline kernel
+fileName = "Feedback.hdf5"
+
+#---------------------------------------------------
+numPart = L**3
+mass = boxSize**3 * rho / numPart
+internalEnergy = P / ((gamma - 1.)*rho)
+
+#--------------------------------------------------
+
+#File
+file = h5py.File(fileName, 'w')
+
+# Header
+grp = file.create_group("/Header")
+grp.attrs["BoxSize"] = boxSize
+grp.attrs["NumPart_Total"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_Total_HighWord"] = [0, 0, 0, 0, 0, 0]
+grp.attrs["NumPart_ThisFile"] = [numPart, 0, 0, 0, 0, 0]
+grp.attrs["Time"] = 0.0
+grp.attrs["NumFilesPerSnapshot"] = 1
+grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+grp.attrs["Flag_Entropy_ICs"] = 0
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = periodic
+
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+#Particle group
+grp = file.create_group("/PartType0")
+
+v = zeros((numPart, 3))
+ds = grp.create_dataset('Velocities', (numPart, 3), 'f')
+ds[()] = v
+v = zeros(1)
+
+m = full((numPart, 1), mass)
+ds = grp.create_dataset('Masses', (numPart,1), 'f')
+ds[()] = m
+m = zeros(1)
+
+h = full((numPart, 1), eta * boxSize / L)
+ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
+ds[()] = h
+h = zeros(1)
+
+u = full((numPart, 1), internalEnergy)
+ds = grp.create_dataset('InternalEnergy', (numPart,1), 'f')
+ds[()] = u
+u = zeros(1)
+
+
+ids = linspace(0, numPart, numPart, endpoint=False).reshape((numPart,1))
+ds = grp.create_dataset('ParticleIDs', (numPart, 1), 'L')
+ds[()] = ids + 1
+x = ids % L;
+y = ((ids - x) / L) % L;
+z = (ids - x - L * y) / L**2;
+coords = zeros((numPart, 3))
+coords[:,0] = z[:,0] * boxSize / L + boxSize / (2*L)
+coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L)
+coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L)
+ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
+ds[()] = coords
+
+file.close()
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/main.c b/examples/main.c
index cf252c16f7efc08ee5bc6bf36f58458583631e87..2833ec3169ea93a99eefe8f126356d5848f9fd13 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -76,6 +76,7 @@ void print_help_message() {
"Overwrite the CPU frequency (Hz) to be used for time measurements");
printf(" %2s %8s %s\n", "-g", "",
"Run with an external gravitational potential");
+ printf(" %2s %8s %s\n", "-F", "", "Run with feedback ");
printf(" %2s %8s %s\n", "-G", "", "Run with self-gravity");
printf(" %2s %8s %s\n", "-n", "{int}",
"Execute a fixed number of time steps. When unset use the time_end "
@@ -147,6 +148,7 @@ int main(int argc, char *argv[]) {
int nsteps = -2;
int with_cosmology = 0;
int with_external_gravity = 0;
+ int with_sourceterms = 0;
int with_cooling = 0;
int with_self_gravity = 0;
int with_hydro = 0;
@@ -159,7 +161,7 @@ int main(int argc, char *argv[]) {
/* Parse the parameters */
int c;
- while ((c = getopt(argc, argv, "acCdDef:gGhn:st:v:y:")) != -1) switch (c) {
+ while ((c = getopt(argc, argv, "acCdDef:FgGhn:st:v:y:")) != -1) switch (c) {
case 'a':
with_aff = 1;
break;
@@ -185,6 +187,9 @@ int main(int argc, char *argv[]) {
return 1;
}
break;
+ case 'F':
+ with_sourceterms = 1;
+ break;
case 'g':
with_external_gravity = 1;
break;
@@ -444,6 +449,11 @@ int main(int argc, char *argv[]) {
if (with_cooling) cooling_init(params, &us, &prog_const, &cooling_func);
if (with_cooling && myrank == 0) cooling_print(&cooling_func);
+ /* Initialise the feedback properties */
+ struct sourceterms sourceterms;
+ if (with_sourceterms) sourceterms_init(params, &us, &sourceterms);
+ if (with_sourceterms && myrank == 0) sourceterms_print(&sourceterms);
+
/* Construct the engine policy */
int engine_policies = ENGINE_POLICY | engine_policy_steal;
if (with_drift_all) engine_policies |= engine_policy_drift_all;
@@ -452,13 +462,14 @@ int main(int argc, char *argv[]) {
if (with_external_gravity) engine_policies |= engine_policy_external_gravity;
if (with_cosmology) engine_policies |= engine_policy_cosmology;
if (with_cooling) engine_policies |= engine_policy_cooling;
+ if (with_sourceterms) engine_policies |= engine_policy_sourceterms;
/* Initialize the engine with the space and policies. */
if (myrank == 0) clocks_gettime(&tic);
struct engine e;
engine_init(&e, &s, params, nr_nodes, myrank, nr_threads, with_aff,
engine_policies, talking, &us, &prog_const, &hydro_properties,
- &potential, &cooling_func);
+ &potential, &cooling_func, &sourceterms);
if (myrank == 0) {
clocks_gettime(&toc);
message("engine_init took %.3f %s.", clocks_diff(&tic, &toc),
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index be2f512613edffd72eaf03689bccee0dd755726b..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:
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 f7cb52ba40a34269173ab8c5019d3011b0c34b61..5859c59d9ddbc7e701145c40f8f0a35b92a075c3 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -43,7 +43,8 @@ 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 potential.h version.h \
- hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h
+ hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.h \
+ sourceterms_struct.h
# Common source files
@@ -52,7 +53,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.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 potential.c hydro_properties.c \
- runner_doiact_fft.c threadpool.c cooling.c
+ runner_doiact_fft.c threadpool.c cooling.c sourceterms.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
@@ -62,6 +63,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h
gravity.h gravity_io.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
+ sourceterms.h \
hydro.h hydro_io.h \
hydro/Minimal/hydro.h hydro/Minimal/hydro_iact.h hydro/Minimal/hydro_io.h \
hydro/Minimal/hydro_debug.h hydro/Minimal/hydro_part.h \
@@ -69,6 +71,8 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.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/riemann_hllc.h riemann/riemann_trrs.h \
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/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 b78cc0a8f842770f60777e3986616a175d2f33ca..f87420c86c5765f42fe04e8dee95558edff90cf9 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -168,6 +168,9 @@ struct cell {
/* Task for cooling */
struct task *cooling;
+ /* Task for source terms */
+ struct task *sourceterms;
+
/* Number of tasks that are associated with this cell. */
int nr_tasks;
@@ -218,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 316a13c8ebba4785e8e1a66b8307f8681697c938..ef26525b7865bac4c0164f7678398ca307fc1d67 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
@@ -95,6 +96,10 @@
//#define EXTERNAL_POTENTIAL_ISOTHERMALPOTENTIAL
//#define EXTERNAL_POTENTIAL_DISC_PATCH
+/* Source terms */
+#define SOURCETERMS_NONE
+//#define SOURCETERMS_SN_FEEDBACK
+
/* Cooling properties */
#define COOLING_NONE
//#define COOLING_CONST_DU
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 507155e3edb468d268d3a0fd7170414b9cebfb84..0319dfd81fa5c3ebf68b00a04946499326eda261 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -68,7 +68,7 @@
#include "units.h"
#include "version.h"
-const char *engine_policy_names[15] = {"none",
+const char *engine_policy_names[16] = {"none",
"rand",
"steal",
"keep",
@@ -82,7 +82,8 @@ const char *engine_policy_names[15] = {"none",
"external_gravity",
"cosmology_integration",
"drift_all",
- "cooling"};
+ "cooling",
+ "sourceterms"};
/** The rank of the engine as a global variable (for messages). */
int engine_rank;
@@ -187,6 +188,8 @@ void engine_make_hydro_hierarchical_tasks(struct engine *e, struct cell *c,
const int is_fixdt = (e->policy & engine_policy_fixdt) == engine_policy_fixdt;
const int is_with_cooling =
(e->policy & engine_policy_cooling) == engine_policy_cooling;
+ const int is_with_sourceterms =
+ (e->policy & engine_policy_sourceterms) == engine_policy_sourceterms;
/* Is this the super-cell? */
if (super == NULL && (c->density != NULL || (c->count > 0 && !c->split))) {
@@ -216,16 +219,18 @@ void engine_make_hydro_hierarchical_tasks(struct engine *e, struct cell *c,
/* Generate the ghost task. */
c->ghost = scheduler_addtask(s, task_type_ghost, task_subtype_none, 0, 0,
c, NULL, 0);
-
#ifdef EXTRA_HYDRO_LOOP
/* Generate the extra ghost task. */
c->extra_ghost = scheduler_addtask(s, task_type_extra_ghost,
task_subtype_none, 0, 0, c, NULL, 0);
#endif
-
if (is_with_cooling)
c->cooling = scheduler_addtask(s, task_type_cooling, task_subtype_none,
0, 0, c, NULL, 0);
+ /* add source terms */
+ if (is_with_sourceterms)
+ c->sourceterms = scheduler_addtask(s, task_type_sourceterms,
+ task_subtype_none, 0, 0, c, NULL, 0);
}
}
@@ -1782,12 +1787,18 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
scheduler_addunlock(sched, t->ci->init, t);
scheduler_addunlock(sched, t, t->ci->kick);
}
-
- /* Cooling tasks should depend on kick and does not unlock anything since
- it is the last task*/
+ /* Cooling tasks should depend on kick and unlock sourceterms */
else if (t->type == task_type_cooling) {
scheduler_addunlock(sched, t->ci->kick, t);
}
+ /* source terms depend on cooling if performed, else on kick. It is the last
+ task */
+ else if (t->type == task_type_sourceterms) {
+ if (e->policy == engine_policy_cooling)
+ scheduler_addunlock(sched, t->ci->cooling, t);
+ else
+ scheduler_addunlock(sched, t->ci->kick, t);
+ }
}
}
@@ -2196,6 +2207,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);
@@ -2611,7 +2624,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);
@@ -2763,6 +2784,11 @@ void engine_step(struct engine *e) {
mask |= 1 << task_type_cooling;
}
+ /* Add the tasks corresponding to sourceterms to the masks */
+ if (e->policy & engine_policy_sourceterms) {
+ mask |= 1 << task_type_sourceterms;
+ }
+
/* Add MPI tasks if need be */
if (e->policy & engine_policy_mpi) {
@@ -2773,6 +2799,8 @@ void engine_step(struct engine *e) {
engine_print_task_counts(e);
+ if (e->verbose) engine_print_task_counts(e);
+
/* Send off the runners. */
TIMER_TIC;
engine_launch(e, e->nr_threads, mask, submask);
@@ -3106,6 +3134,7 @@ void engine_unpin() {
* @param hydro The #hydro_props used for this run.
* @param potential The properties of the external potential.
* @param cooling_func The properties of the cooling function.
+ * @param sourceterms The properties of the source terms function.
*/
void engine_init(struct engine *e, struct space *s,
const struct swift_params *params, int nr_nodes, int nodeID,
@@ -3114,7 +3143,8 @@ void engine_init(struct engine *e, struct space *s,
const struct phys_const *physical_constants,
const struct hydro_props *hydro,
const struct external_potential *potential,
- const struct cooling_function_data *cooling_func) {
+ const struct cooling_function_data *cooling_func,
+ struct sourceterms *sourceterms) {
/* Clean-up everything */
bzero(e, sizeof(struct engine));
@@ -3167,6 +3197,7 @@ void engine_init(struct engine *e, struct space *s,
e->hydro_properties = hydro;
e->external_potential = potential;
e->cooling_func = cooling_func;
+ e->sourceterms = sourceterms;
e->parameter_file = params;
engine_rank = nodeID;
diff --git a/src/engine.h b/src/engine.h
index d36914af611723bc4d496d5c2dc68050eea6ffe6..ec9e16553c3172f22e9bfe60971163488947a47e 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -44,6 +44,7 @@
#include "potential.h"
#include "runner.h"
#include "scheduler.h"
+#include "sourceterms_struct.h"
#include "space.h"
#include "task.h"
#include "units.h"
@@ -65,6 +66,7 @@ enum engine_policy {
engine_policy_cosmology = (1 << 11),
engine_policy_drift_all = (1 << 12),
engine_policy_cooling = (1 << 13),
+ engine_policy_sourceterms = (1 << 14)
};
extern const char *engine_policy_names[];
@@ -207,6 +209,9 @@ struct engine {
/* Properties of the cooling scheme */
const struct cooling_function_data *cooling_func;
+ /* Properties of source terms */
+ struct sourceterms *sourceterms;
+
/* The (parsed) parameter file */
const struct swift_params *parameter_file;
};
@@ -223,7 +228,8 @@ void engine_init(struct engine *e, struct space *s,
const struct phys_const *physical_constants,
const struct hydro_props *hydro,
const struct external_potential *potential,
- const struct cooling_function_data *cooling);
+ const struct cooling_function_data *cooling,
+ struct sourceterms *sourceterms);
void engine_launch(struct engine *e, int nr_runners, unsigned int mask,
unsigned int submask);
void engine_prepare(struct engine *e, int nodrift);
diff --git a/src/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/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/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 2b17505c36f345779b269c90758931a13f9b4e0d..0bdf56fe844309998db1fad4cf9581edb6b88305 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -37,6 +37,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
@@ -116,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
@@ -126,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
@@ -141,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;
}
/**
@@ -195,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;
@@ -222,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;
}
/**
@@ -273,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;
@@ -349,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;
}
/**
@@ -400,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;
}
/**
@@ -414,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 5ee8cd0370a970aa83cef3d1c8909d923c12ba24..fca1bcff91a71fb2a26adc117382630a576f9090 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -59,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;
@@ -72,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;
@@ -209,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];
@@ -248,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;
@@ -366,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];
@@ -415,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;
@@ -447,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;
@@ -690,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;
@@ -722,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 2c12a71ad6256d7372256b6e590790ee0ff4e22e..0c14da957463720dcee5349e88be91986cbb3f54 100644
--- a/src/hydro/Gadget2/hydro_part.h
+++ b/src/hydro/Gadget2/hydro_part.h
@@ -80,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 {
@@ -93,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];
@@ -106,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;
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/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index eeb389537c56876126c60b2d29a728029c72844b..3015f26c6bad7f006e8cda16695c750ff40d74df 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -39,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
@@ -65,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;
}
/**
@@ -97,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;
}
/**
@@ -128,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
@@ -138,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
@@ -153,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;
}
/**
@@ -214,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;
}
/**
@@ -240,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);
}
/**
@@ -274,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;
}
/**
@@ -337,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;
}
/**
@@ -379,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;
}
/**
@@ -392,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 edb060e4fd71fcc136e1bedf6e8a752d1d50d54f..169947b99e92d9bd1b0870d502a49e311820ff81 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -55,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;
@@ -65,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;
}
@@ -100,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;
}
@@ -152,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] +
@@ -165,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 */
@@ -263,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] +
@@ -276,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 8e23bddf5153043421319810683266b27d297f93..8542177278998d5e0b830dc164988611549ef24d 100644
--- a/src/hydro/Minimal/hydro_part.h
+++ b/src/hydro/Minimal/hydro_part.h
@@ -93,9 +93,6 @@ struct part {
/*! Particle density. */
float rho;
- /*! Derivative of density with respect to h */
- float rho_dh;
-
/* Store density/force specific stuff. */
union {
@@ -114,6 +111,9 @@ struct part {
/*! Derivative of the neighbour number with respect to h. */
float wcount_dh;
+ /*! Derivative of density with respect to h */
+ float rho_dh;
+
} density;
/**
@@ -125,9 +125,15 @@ struct part {
*/
struct {
+ /*! "Grad h" term */
+ float f;
+
/*! Particle pressure. */
float pressure;
+ /*! Particle soundspeed. */
+ float soundspeed;
+
/*! Particle signal velocity */
float v_sig;
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 af39d10fafc11d4435ddbcc087fbf08178b18959..6832e0c6c2e0f2c324d90629447cf6a5e809d6fb 100644
--- a/src/part.h
+++ b/src/part.h
@@ -42,12 +42,19 @@
/* 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/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 274976077fe285c749816800fda166795da094b3..a554d059607aad10fea8010906ad4d2617759729 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -53,6 +53,7 @@
#include "kick.h"
#include "minmax.h"
#include "scheduler.h"
+#include "sourceterms.h"
#include "space.h"
#include "task.h"
#include "timers.h"
@@ -111,6 +112,42 @@ const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
#include "runner_doiact_fft.h"
#include "runner_doiact_grav.h"
+/**
+ * @brief Perform source terms
+ *
+ * @param r runner task
+ * @param c cell
+ * @param timer 1 if the time is to be recorded.
+ */
+void runner_do_sourceterms(struct runner *r, struct cell *c, int timer) {
+ const int count = c->count;
+ const double cell_min[3] = {c->loc[0], c->loc[1], c->loc[2]};
+ const double cell_width[3] = {c->width[0], c->width[1], c->width[2]};
+ struct sourceterms *sourceterms = r->e->sourceterms;
+ const int dimen = 3;
+
+ TIMER_TIC;
+
+ /* Recurse? */
+ if (c->split) {
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) runner_do_sourceterms(r, c->progeny[k], 0);
+ return;
+ }
+
+ if (count > 0) {
+
+ /* do sourceterms in this cell? */
+ const int incell =
+ sourceterms_test_cell(cell_min, cell_width, sourceterms, dimen);
+ if (incell == 1) {
+ sourceterms_apply(r, sourceterms, c);
+ }
+ }
+
+ if (timer) TIMER_TOC(timer_dosource);
+}
+
/**
* @brief Calculate gravity acceleration from external potential
*
@@ -1397,6 +1434,9 @@ void *runner_main(void *data) {
case task_type_cooling:
runner_do_cooling(r, t->ci, 1);
break;
+ case task_type_sourceterms:
+ runner_do_sourceterms(r, t->ci, 1);
+ break;
default:
error("Unknown task type.");
}
diff --git a/src/scheduler.c b/src/scheduler.c
index a45aa6570d647d0aa0ad95e102d717b9895d5340..29cdbe86ddc9f6eff31766433148b36bf64d38f6 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 &&
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/sourceterms.c b/src/sourceterms.c
new file mode 100644
index 0000000000000000000000000000000000000000..2a5f326688b52b0e4abfde8dca823d765e0e8a5e
--- /dev/null
+++ b/src/sourceterms.c
@@ -0,0 +1,60 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@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"
+
+/* Local includes. */
+#include "const.h"
+#include "hydro.h"
+#include "parser.h"
+#include "units.h"
+
+/* This object's header. */
+#include "sourceterms.h"
+
+/**
+ * @brief Initialises the sourceterms
+ *
+ * @param parameter_file The parsed parameter file
+ * @param us The current internal system of units
+ * @param source the structure that has all the source term properties
+ */
+void sourceterms_init(const struct swift_params* parameter_file,
+ struct UnitSystem* us, struct sourceterms* source) {
+#ifdef SOURCETERMS_SN_FEEDBACK
+ supernova_init(parameter_file, us, source);
+#endif /* SOURCETERMS_SN_FEEDBACK */
+};
+
+/**
+ * @brief Prints the properties of the source terms to stdout
+ * @param source the structure that has all the source term properties
+ */
+void sourceterms_print(struct sourceterms* source) {
+#ifdef SOURCETERMS_NONE
+ error(" no sourceterms defined yet you ran with -F");
+#ifdef SOURCETERMS_SN_FEEDBACK
+#error "can't have sourceterms when defined SOURCETERMS_NONE"
+#endif
+#endif
+#ifdef SOURCETERMS_SN_FEEDBACK
+ supernova_print(source);
+#endif /* SOURCETERMS_SN_FEEDBACK */
+};
diff --git a/src/sourceterms.h b/src/sourceterms.h
new file mode 100644
index 0000000000000000000000000000000000000000..95c37e9bba9881702a6cf6caea1004b8cfb52636
--- /dev/null
+++ b/src/sourceterms.h
@@ -0,0 +1,74 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2015 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_SOURCETERMS_H
+#define SWIFT_SOURCETERMS_H
+
+/**
+ * @file src/sourceterms.h
+ * @brief Branches between the different sourceterms functions.
+ */
+
+#include "./const.h"
+#include "runner.h"
+
+#ifdef SOURCETERMS_SN_FEEDBACK
+#include "sourceterms/sn_feedback/sn_feedback_struct.h"
+#endif
+
+/* So far only one model here */
+struct sourceterms {
+#ifdef SOURCETERMS_SN_FEEDBACK
+ struct supernova_struct supernova;
+#endif
+};
+#ifdef SOURCETERMS_SN_FEEDBACK
+#include "sourceterms/sn_feedback/sn_feedback.h"
+#endif
+
+void sourceterms_init(const struct swift_params* parameter_file,
+ struct UnitSystem* us, struct sourceterms* source);
+void sourceterms_print(struct sourceterms* source);
+
+/**
+ * @brief Routines related to source terms
+ * @param cell_min: corner of cell to test
+ * @param cell_width: width of cell to test
+ * @param sourceterms: properties of source terms to test
+ * @param dimen: dimensionality of the problem
+ *
+ * This routine tests whether a source term should be applied to this cell
+ * return: 1 if yes, return: 0 if no
+ */
+
+__attribute__((always_inline)) INLINE static int sourceterms_test_cell(
+ const double cell_min[], const double cell_width[],
+ struct sourceterms* sourceterms, const int dimen) {
+#ifdef SOURCETERMS_SN_FEEDBACK
+ return supernova_feedback_test_cell(cell_min, cell_width, sourceterms, dimen);
+#endif
+ return 0;
+};
+
+__attribute__((always_inline)) INLINE static void sourceterms_apply(
+ struct runner* r, struct sourceterms* sourceterms, struct cell* c) {
+#ifdef SOURCETERMS_SN_FEEDBACK
+ supernova_feedback_apply(r, sourceterms, c);
+#endif
+};
+#endif /* SWIFT_SOURCETERMS_H */
diff --git a/src/sourceterms/sn_feedback/sn_feedback.h b/src/sourceterms/sn_feedback/sn_feedback.h
new file mode 100644
index 0000000000000000000000000000000000000000..667c79cb1eef42f7fde79c43059d1baa5c61268e
--- /dev/null
+++ b/src/sourceterms/sn_feedback/sn_feedback.h
@@ -0,0 +1,192 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@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_SN_FEEDBACK_H
+#define SWIFT_SN_FEEDBACK_H
+#include <float.h>
+/* Config parameters. */
+#include "../config.h"
+
+#include "engine.h"
+#include "equation_of_state.h"
+#include "hydro.h"
+#include "runner.h"
+#include "timestep.h"
+
+/**
+ * @file src/sourceterms/sn_feedback.h
+ *
+ * @brief Routines related to sourceterms (supernova feedback): determine if
+ * feedback occurs in this cell
+ *
+ * @param cell_min: corner of cell to test
+ * @param cell_width: width of cell to test
+ * @param sourceterms: properties of source terms to test
+ * @param dimen: dimensionality of the problem
+ *
+ * This routine tests whether a source term should be applied to this cell
+ * return: 1 if yes, return: 0 if no
+ */
+__attribute__((always_inline)) INLINE static int supernova_feedback_test_cell(
+ const double cell_min[], const double cell_width[],
+ struct sourceterms* sourceterms, const int dimen) {
+ if (sourceterms->supernova.status == supernova_is_done) return 0;
+
+ const double location[3] = {sourceterms->supernova.x,
+ sourceterms->supernova.y,
+ sourceterms->supernova.z};
+ for (int i = 0; i < dimen; i++) {
+ if (cell_min[i] > location[i]) return 0;
+ if ((cell_min[i] + cell_width[i]) <= location[i]) return 0;
+ };
+ return 1;
+};
+
+/**
+ * @file src/sourceterms/sn_feedback.h
+ *
+ * @brief Routines related to source terms (supernova feedback): perform
+ * feedback in this cell
+ * @param r: the runner
+ * @param sourceterms the structure describing the source terms properties
+ * @param c the cell to apply feedback to
+ *
+ * This routine heats an individual particle (p), increasing its thermal energy
+ * per unit mass
+ * by supernova energy / particle mass.
+ */
+__attribute__((always_inline)) INLINE static void supernova_feedback_apply(
+ struct runner* restrict r, struct sourceterms* restrict sourceterms,
+ struct cell* restrict c) {
+
+ const int count = c->count;
+ struct part* restrict parts = c->parts;
+ struct xpart* restrict xparts = c->xparts;
+ const double timeBase = r->e->timeBase;
+ const int ti_current = r->e->ti_current;
+
+ /* inject SN energy into the particle with highest id in this cell if it is
+ * active */
+ int imax = 0;
+ struct part* restrict p_sn = NULL;
+ struct xpart* restrict xp_sn = NULL;
+
+ for (int i = 0; i < count; i++) {
+
+ /* Get a direct pointer on the part. */
+ struct part* restrict p = &parts[i];
+ if (p->id > imax) {
+ imax = p->id;
+ p_sn = p;
+ xp_sn = &xparts[i];
+ }
+ }
+
+ /* Is this part within the time step? */
+ if (p_sn->ti_begin == ti_current) {
+
+ /* Does this time step straddle the feedback injection time? */
+ const float t_begin = p_sn->ti_begin * timeBase;
+ const float t_end = p_sn->ti_end * timeBase;
+ if (t_begin <= sourceterms->supernova.time &&
+ t_end > sourceterms->supernova.time) {
+
+ /* store old time step */
+ const int dti_old = p_sn->ti_end - p_sn->ti_begin;
+
+ /* add supernova feedback */
+ const float u_old = hydro_get_internal_energy(p_sn, 0);
+ const float ent_old = hydro_get_entropy(p_sn, 0.0);
+ const float u_new =
+ u_old + sourceterms->supernova.energy / hydro_get_mass(p_sn);
+ hydro_set_internal_energy(p_sn, u_new);
+ const float u_set = hydro_get_internal_energy(p_sn, 0.0);
+ const float ent_set = hydro_get_entropy(p_sn, 0.0);
+ message(
+ " applied super nova, time = %e, location= %e %e %e velocity= %e %e "
+ "%e",
+ ti_current * timeBase, p_sn->x[0], p_sn->x[1], p_sn->x[2], p_sn->v[0],
+ p_sn->v[1], p_sn->v[2]);
+ message(
+ " injected SN energy in particle = %lld, increased energy from %e to "
+ "%e and is notw %e, entropy from %e to %e",
+ p_sn->id, u_old, u_new, u_set, ent_old, ent_set);
+
+ /* label supernova as done */
+ sourceterms->supernova.status = supernova_is_done;
+
+ /* update timestep if new time step shorter than old time step */
+ const int dti = get_part_timestep(p_sn, xp_sn, r->e);
+ if (dti < dti_old) {
+ p_sn->ti_end = p_sn->ti_begin + dti;
+ message(" changed timestep from %d to %d", dti_old, dti);
+
+ /* apply simple time-step limiter on all particles in same cell:
+ */
+ int i_limit = 0;
+ for (int i = 0; i < count; i++) {
+ struct part* restrict p = &parts[i];
+ const int dti_old = p->ti_end - p->ti_begin;
+ if (dti_old > 2 * dti) {
+ i_limit++;
+ const int dti_new = 2 * dti;
+ p->ti_end = p->ti_begin + dti_new;
+ message(" old step = %d new step = %d", dti_old, dti_new);
+ } else
+ message(" old step = %d", dti_old);
+ }
+ message(" count= %d limited timestep of %d particles ", count, i_limit);
+ } /* end of limiter */
+ error("end");
+ }
+ }
+};
+
+/**
+ * @file src/sourceterms/sn_feedback.h
+ *
+ * @brief Routine to initialise supernova feedback
+ * @param parameterfile: the parse parmeter file
+ * @param us: the unit system in use
+ * @param sourceterms the structure describing the source terms properties
+ *
+ * This routine heats an individual particle (p), increasing its thermal energy
+ * per unit mass
+ * by supernova energy / particle mass.
+ */
+
+__attribute__((always_inline)) INLINE static void supernova_init(
+ const struct swift_params* parameter_file, struct UnitSystem* us,
+ struct sourceterms* source) {
+ source->supernova.time = parser_get_param_double(parameter_file, "SN:time");
+ source->supernova.energy =
+ parser_get_param_double(parameter_file, "SN:energy");
+ source->supernova.x = parser_get_param_double(parameter_file, "SN:x");
+ source->supernova.y = parser_get_param_double(parameter_file, "SN:y");
+ source->supernova.z = parser_get_param_double(parameter_file, "SN:z");
+ source->supernova.status = supernova_is_not_done;
+}
+__attribute__((always_inline)) INLINE static void supernova_print(
+ struct sourceterms* source) {
+ message(
+ " Single SNe of energy= %e will explode at time= %e at location "
+ "(%e,%e,%e)",
+ source->supernova.energy, source->supernova.time, source->supernova.x,
+ source->supernova.y, source->supernova.z);
+}
+#endif /* SWIFT_SN_FEEDBACK_H */
diff --git a/src/sourceterms/sn_feedback/sn_feedback_struct.h b/src/sourceterms/sn_feedback/sn_feedback_struct.h
new file mode 100644
index 0000000000000000000000000000000000000000..2b4dbe97b0a687cc48d0ba7e12b105f306cb7e96
--- /dev/null
+++ b/src/sourceterms/sn_feedback/sn_feedback_struct.h
@@ -0,0 +1,45 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Tom Theuns (tom.theuns@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/>.
+ *
+ ******************************************************************************/
+/**
+ * @file src/sourceterms/sn_feedback_struct.h
+ * @brief Routines related to source terms (feedback)
+ *
+ * enumeration type that sets if supernova explosion is done (is_done) or still
+ * needs doing (is_not_done)
+ */
+#ifndef SWIFT_SN_FEEDBACK_STRUCT_H
+#define SWIFT_SN_FEEDBACK_STRUCT_H
+enum supernova_status { supernova_is_done, supernova_is_not_done };
+
+/**
+ * @file src/sourceterms/sn_feedback_struct.h
+ * @brief Routines related to source terms (feedback)
+ *
+ * The structure that describes the source term (supernova feedback)
+ * It specifies the time, energy and location of the desired supernova
+ * explosion, and a status (supernova_is_done/supernova_is_not_done)
+ * that records the status of the supernova
+ */
+struct supernova_struct {
+ double time;
+ double energy;
+ double x, y, z;
+ enum supernova_status status;
+};
+#endif SWIFT_SN_FEEDBACK_STRUCT_H
diff --git a/src/sourceterms_struct.h b/src/sourceterms_struct.h
new file mode 100644
index 0000000000000000000000000000000000000000..b3c38986db52d72df825fda97b36c985dff922b6
--- /dev/null
+++ b/src/sourceterms_struct.h
@@ -0,0 +1,26 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2015 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_SOURCETERMS_STRUCT_H
+#define SWIFT_SOURCETERMS_STRUCT_H
+#include "./const.h"
+#ifdef SOURCETERMS_SN_FEEDBACK
+#include "sourceterms/sn_feedback/sn_feedback_struct.h"
+#endif
+
+#endif /* SWIFT_SOURCETERMS_STRUCT_H */
diff --git a/src/space.c b/src/space.c
index a9958f6fbd7d85060db99a9682b0de10f507085d..1d2733ab2d4356c19cc0e65a0b9441edcf0acb80 100644
--- a/src/space.c
+++ b/src/space.c
@@ -58,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] = {
@@ -228,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))
@@ -691,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),
@@ -704,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)
@@ -719,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.
@@ -1240,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];
@@ -1269,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) {
@@ -1565,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);
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 3119adfd4b95639a84137e243c6e971d80ef3825..58745f3111fc3f29490d4591ee549907392cb87e 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -49,6 +49,7 @@
#include "scheduler.h"
#include "serial_io.h"
#include "single_io.h"
+#include "sourceterms.h"
#include "space.h"
#include "task.h"
#include "timers.h"
diff --git a/src/task.c b/src/task.c
index 8da526f57886fa68121061c8eae2316912e73a30..00068f45769b4ca606cc729bd5e89c13ae729eef 100644
--- a/src/task.c
+++ b/src/task.c
@@ -51,7 +51,7 @@ const char *taskID_names[task_type_count] = {
"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"};
+ "grav_mm", "grav_up", "grav_external", "cooling", "sourceterms"};
const char *subtaskID_names[task_subtype_count] = {
"none", "density", "gradient", "force", "grav", "tend"};
@@ -118,6 +118,7 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
case task_type_ghost:
case task_type_extra_ghost:
case task_type_cooling:
+ case task_type_sourceterms:
return task_action_part;
break;
diff --git a/src/task.h b/src/task.h
index 4928dc00bcd7958efdd6987b5aec90ab4b3e92fa..bc4df3dc2a4cfee3c382e9f2059cba84f29299f7 100644
--- a/src/task.h
+++ b/src/task.h
@@ -55,6 +55,7 @@ enum task_types {
task_type_grav_up,
task_type_grav_external,
task_type_cooling,
+ task_type_sourceterms,
task_type_count
} __attribute__((packed));
diff --git a/src/timers.h b/src/timers.h
index b93a34df4d90251d4686afaef0be51b36dc9a25e..d75508afd88cf2fd57d06f9530bff8607d79d127 100644
--- a/src/timers.h
+++ b/src/timers.h
@@ -45,6 +45,7 @@ enum {
timer_dopair_grav_pm,
timer_dopair_grav_pp,
timer_dograv_external,
+ timer_dosource,
timer_dosub_self_density,
timer_dosub_self_gradient,
timer_dosub_self_force,
diff --git a/src/timestep.h b/src/timestep.h
index 599fb4762e11b08fc942fb02acbbf1970f477de4..db52911ec1e8fbf31f35e8877e0a7ae7ba5ee478 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -24,8 +24,8 @@
/* Local headers. */
#include "const.h"
+#include "cooling.h"
#include "debug.h"
-
/**
* @brief Compute a valid integer time-step form a given time-step
*
diff --git a/tests/test125cells.c b/tests/test125cells.c
index e666658f43de135e3e72521b52f2a688c596a6f6..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)
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/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