diff --git a/INSTALL.swift b/INSTALL.swift
index 8e1635b0715426512503fd9dcde32f59a7ad1b62..90096bc1b88d34ab86c04f623f1d3f04ca5a2997 100644
--- a/INSTALL.swift
+++ b/INSTALL.swift
@@ -143,6 +143,8 @@ before you can build it.
- DOXYGEN: the doxygen library is required to create the SWIFT API
documentation.
+ - python: Examples and solution script use python and rely on the
+ numpy library version 1.8.2 or higher.
SWIFT Coding style
diff --git a/README b/README
index 0d658c333f1328b423851031c5b5d202f43df3c2..bec538bfc35239945d60487d66ed2e02b5d363a2 100644
--- a/README
+++ b/README
@@ -15,25 +15,26 @@ Usage: swift [OPTION]... PARAMFILE
swift_mpi [OPTION]... PARAMFILE
Valid options are:
- -a Pin runners using processor affinity
- -c Run with cosmological time integration
- -C Run with cooling
+ -a Pin runners using processor affinity.
+ -c Run with cosmological time integration.
+ -C Run with cooling.
-d Dry run. Read the parameter file, allocate memory but does not read
the particles from ICs and exit before the start of time integration.
Allows user to check validy of parameter and IC files as well as memory limits.
- -D Always drift all particles even the ones far from active particles.
- -e Enable floating-point exceptions (debugging mode)
- -f {int} Overwrite the CPU frequency (Hz) to be used for time measurements
- -g Run with an external gravitational potential
- -G Run with self-gravity
+ -D Always drift all particles even the ones far from active particles. This emulates
+ Gadget-[23] and GIZMO's default behaviours.
+ -e Enable floating-point exceptions (debugging mode).
+ -f {int} Overwrite the CPU frequency (Hz) to be used for time measurements.
+ -g Run with an external gravitational potential.
+ -G Run with self-gravity.
-n {int} Execute a fixed number of time steps. When unset use the time_end parameter to stop.
- -s Run with SPH
- -S Run with stars
+ -s Run with hydrodynamics.
+ -S Run with stars.
-t {int} The number of threads to use on each MPI rank. Defaults to 1 if not specified.
-v [12] Increase the level of verbosity
1: MPI-rank 0 writes
2: All MPI-ranks write
- -y {int} Time-step frequency at which task graphs are dumped
- -h Print this help message and exit
+ -y {int} Time-step frequency at which task graphs are dumped.
+ -h Print this help message and exit.
See the file examples/parameter_example.yml for an example of parameter file.
diff --git a/configure.ac b/configure.ac
index 4b6308e96b81bbfd0a9256bb4914f1356fbfa6f8..c94f5c23b9cc5868a199ad4d177053c6f94f639d 100644
--- a/configure.ac
+++ b/configure.ac
@@ -16,7 +16,7 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# Init the project.
-AC_INIT([SWIFT],[0.4.0],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
+AC_INIT([SWIFT],[0.5.0],[https://gitlab.cosma.dur.ac.uk/swift/swiftsim])
swift_config_flags="$*"
# Need to define this, instead of using fifth argument of AC_INIT, until 2.64.
@@ -200,6 +200,20 @@ if test "$enable_debugging_checks" = "yes"; then
AC_DEFINE([SWIFT_DEBUG_CHECKS],1,[Enable expensive debugging])
fi
+# Check if gravity force checks are on for some particles.
+AC_ARG_ENABLE([gravity-force-checks],
+ [AS_HELP_STRING([--enable-gravity-force-checks],
+ [Activate expensive brute-force gravity checks for a fraction 1/N of all particles @<:@N@:>@]
+ )],
+ [gravity_force_checks="$enableval"],
+ [gravity_force_checks="no"]
+)
+if test "$gravity_force_checks" == "yes"; then
+ AC_MSG_ERROR(Need to specify the fraction of particles to check when using --enable-gravity-force-checks!)
+elif test "$gravity_force_checks" != "no"; then
+ AC_DEFINE_UNQUOTED([SWIFT_GRAVITY_FORCE_CHECKS], [$enableval] ,[Enable gravity brute-force checks])
+fi
+
# Define HAVE_POSIX_MEMALIGN if it works.
AX_FUNC_POSIX_MEMALIGN
@@ -637,13 +651,13 @@ AC_ARG_WITH([hydro-dimension],
)
case "$with_dimension" in
1)
- AC_DEFINE([HYDRO_DIMENSION_1D], [1], [1D analysis])
+ AC_DEFINE([HYDRO_DIMENSION_1D], [1], [1D solver])
;;
2)
- AC_DEFINE([HYDRO_DIMENSION_2D], [2], [2D analysis])
+ AC_DEFINE([HYDRO_DIMENSION_2D], [2], [2D solver])
;;
3)
- AC_DEFINE([HYDRO_DIMENSION_3D], [3], [3D analysis])
+ AC_DEFINE([HYDRO_DIMENSION_3D], [3], [3D solver])
;;
*)
AC_MSG_ERROR([Dimensionality must be 1, 2 or 3])
@@ -751,7 +765,7 @@ esac
# External potential
AC_ARG_WITH([ext-potential],
[AS_HELP_STRING([--with-ext-potential=<pot>],
- [external potential @<:@none, point-mass, isothermal, softened-isothermal, disc-patch default: none@:>@]
+ [external potential @<:@none, point-mass, isothermal, softened-isothermal, disc-patch, sine-wave default: none@:>@]
)],
[with_potential="$withval"],
[with_potential="none"]
@@ -769,6 +783,9 @@ case "$with_potential" in
disc-patch)
AC_DEFINE([EXTERNAL_POTENTIAL_DISC_PATCH], [1], [Disc-patch external potential])
;;
+ sine-wave)
+ AC_DEFINE([EXTERNAL_POTENTIAL_SINE_WAVE], [1], [Sine wave external potential in 1D])
+ ;;
*)
AC_MSG_ERROR([Unknown external potential: $with_potential])
;;
@@ -821,8 +838,10 @@ AC_MSG_RESULT([
Riemann solver : $with_riemann
Cooling function : $with_cooling
External potential : $with_potential
+
Task debugging : $enable_task_debugging
Debugging checks : $enable_debugging_checks
+ Gravity checks : $gravity_force_checks
])
# Make sure the latest git revision string gets included
diff --git a/examples/BigCosmoVolume/makeIC.py b/examples/BigCosmoVolume/makeIC.py
index c141337c06fb28aa4049e2823fcc7cd3e9d5513c..8f83b564a3f747d164fd03b7dddf3cd9c609d9eb 100644
--- a/examples/BigCosmoVolume/makeIC.py
+++ b/examples/BigCosmoVolume/makeIC.py
@@ -107,7 +107,7 @@ if n_copy > 1:
for i in range(n_copy):
for j in range(n_copy):
for k in range(n_copy):
- coords = np.append(coords, coords_tile + np.array([ i * boxSize, j * boxSize, k * boxSize ]), axis=0)
+ coords = np.append(coords, coords_tile + np.array([ i * boxSize[0], j * boxSize[1], k * boxSize[2] ]), axis=0)
v = np.append(v, v_tile, axis=0)
m = np.append(m, m_tile)
h = np.append(h, h_tile)
diff --git a/examples/CoolingHalo/makeIC.py b/examples/CoolingHalo/makeIC.py
index 0b542e200da709e2cc7f668ab8b62b94e0bf95ee..3ec1be6f7b5e568ebe8e0fefe508ef8287edb29c 100644
--- a/examples/CoolingHalo/makeIC.py
+++ b/examples/CoolingHalo/makeIC.py
@@ -227,7 +227,7 @@ ds[()] = u
u = np.zeros(1)
# Particle IDs
-ids = 1 + np.linspace(0, N, N, endpoint=False, dtype='L')
+ids = 1 + np.linspace(0, N, N, endpoint=False)
ds = grp.create_dataset('ParticleIDs', (N, ), 'L')
ds[()] = ids
diff --git a/examples/CoolingHaloWithSpin/makeIC.py b/examples/CoolingHaloWithSpin/makeIC.py
index a6d57868ad7542498b27007a5c3ef9234b9feb84..2cf3127c743f61756b3ff6c4a7738c83d185f9cd 100644
--- a/examples/CoolingHaloWithSpin/makeIC.py
+++ b/examples/CoolingHaloWithSpin/makeIC.py
@@ -233,7 +233,7 @@ ds[()] = u
u = np.zeros(1)
# Particle IDs
-ids = 1 + np.linspace(0, N, N, endpoint=False, dtype='L')
+ids = 1 + np.linspace(0, N, N, endpoint=False)
ds = grp.create_dataset('ParticleIDs', (N, ), 'L')
ds[()] = ids
diff --git a/examples/DiscPatch/GravityOnly/makeIC.py b/examples/DiscPatch/GravityOnly/makeIC.py
index 42cd26e235deb17a899a65050ef5caa9c832c59c..6e4e148392eb7ca2fbf8c29c3f737d029916c59b 100644
--- a/examples/DiscPatch/GravityOnly/makeIC.py
+++ b/examples/DiscPatch/GravityOnly/makeIC.py
@@ -150,7 +150,7 @@ ds[()] = m
m = numpy.zeros(1)
-ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False, dtype='L')
+ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False)
ds = grp1.create_dataset('ParticleIDs', (numPart, ), 'L')
ds[()] = ids
diff --git a/examples/DiscPatch/HydroStatic/makeIC.py b/examples/DiscPatch/HydroStatic/makeIC.py
index e2846d82a8cfa8bf08de83632b19ae2e7818f3c1..6ba1ccd06fed84ca728aadaa5922dbba536b6881 100644
--- a/examples/DiscPatch/HydroStatic/makeIC.py
+++ b/examples/DiscPatch/HydroStatic/makeIC.py
@@ -205,7 +205,7 @@ if (entropy_flag == 1):
else:
ds[()] = u
-ids = 1 + numpy.linspace(0, numGas, numGas, endpoint=False, dtype='L')
+ids = 1 + numpy.linspace(0, numGas, numGas, endpoint=False)
ds = grp0.create_dataset('ParticleIDs', (numGas, ), 'L')
ds[()] = ids
diff --git a/examples/HydrostaticHalo/density_profile.py b/examples/HydrostaticHalo/density_profile.py
index 52bebb9ffefa77dae66af155fb31fed539dcde13..d0afd399f951cf3b727e869ca8571a3a802c2e8d 100644
--- a/examples/HydrostaticHalo/density_profile.py
+++ b/examples/HydrostaticHalo/density_profile.py
@@ -28,7 +28,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
unit_time_cgs = unit_length_cgs / unit_velocity_cgs
-v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
+v_c = float(params.attrs["IsothermalPotential:vrot"])
v_c_cgs = v_c * unit_velocity_cgs
#lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
#X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])
diff --git a/examples/HydrostaticHalo/internal_energy_profile.py b/examples/HydrostaticHalo/internal_energy_profile.py
index a1b2bda314a66eb965974d34519f66c544ee8aed..ea52cf8fc5fd098a46f05eaa58494529a868000c 100644
--- a/examples/HydrostaticHalo/internal_energy_profile.py
+++ b/examples/HydrostaticHalo/internal_energy_profile.py
@@ -46,7 +46,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
unit_time_cgs = unit_length_cgs / unit_velocity_cgs
-v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
+v_c = float(params.attrs["IsothermalPotential:vrot"])
v_c_cgs = v_c * unit_velocity_cgs
#lambda_cgs = float(params.attrs["LambdaCooling:lambda_cgs"])
#X_H = float(params.attrs["LambdaCooling:hydrogen_mass_abundance"])
diff --git a/examples/HydrostaticHalo/makeIC.py b/examples/HydrostaticHalo/makeIC.py
index f33387e18dd0ab523684227f5c745b5c8b807b7f..d5081ac84473edc87857c6872278b4d0ca6389b1 100644
--- a/examples/HydrostaticHalo/makeIC.py
+++ b/examples/HydrostaticHalo/makeIC.py
@@ -227,7 +227,7 @@ ds[()] = u
u = np.zeros(1)
# Particle IDs
-ids = 1 + np.linspace(0, N, N, endpoint=False, dtype='L')
+ids = 1 + np.linspace(0, N, N, endpoint=False)
ds = grp.create_dataset('ParticleIDs', (N, ), 'L')
ds[()] = ids
diff --git a/examples/HydrostaticHalo/velocity_profile.py b/examples/HydrostaticHalo/velocity_profile.py
index f6a7350b9731d660b2092266d4d6ad3730bab48c..9133195d942233514148aa419003ee0ab7923494 100644
--- a/examples/HydrostaticHalo/velocity_profile.py
+++ b/examples/HydrostaticHalo/velocity_profile.py
@@ -46,7 +46,7 @@ unit_mass_cgs = float(params.attrs["InternalUnitSystem:UnitMass_in_cgs"])
unit_length_cgs = float(params.attrs["InternalUnitSystem:UnitLength_in_cgs"])
unit_velocity_cgs = float(params.attrs["InternalUnitSystem:UnitVelocity_in_cgs"])
unit_time_cgs = unit_length_cgs / unit_velocity_cgs
-v_c = float(params.attrs["SoftenedIsothermalPotential:vrot"])
+v_c = float(params.attrs["IsothermalPotential:vrot"])
v_c_cgs = v_c * unit_velocity_cgs
header = f["Header"]
N = header.attrs["NumPart_Total"][0]
diff --git a/examples/IsothermalPotential/makeIC.py b/examples/IsothermalPotential/makeIC.py
index 7d1c5361f9a255365517226e49c55a8a50c4d6ce..27ddf15fe63d69d4172b927cfca8b8662d11ea4e 100644
--- a/examples/IsothermalPotential/makeIC.py
+++ b/examples/IsothermalPotential/makeIC.py
@@ -138,7 +138,7 @@ ds = grp1.create_dataset('Masses', (numPart,), 'f')
ds[()] = m
m = numpy.zeros(1)
-ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False, dtype='L')
+ids = 1 + numpy.linspace(0, numPart, numPart, endpoint=False)
ds = grp1.create_dataset('ParticleIDs', (numPart, ), 'L')
ds[()] = ids
diff --git a/examples/IsothermalPotential/run.sh b/examples/IsothermalPotential/run.sh
index 976fbddc01cf7a3dcbb114d437ddb8f03b4d54bd..a5f03f32f82e27660d0a950335d731cf0ff7401d 100755
--- a/examples/IsothermalPotential/run.sh
+++ b/examples/IsothermalPotential/run.sh
@@ -4,7 +4,7 @@
if [ ! -e Isothermal.hdf5 ]
then
echo "Generating initial conditions for the isothermal potential box example..."
- python makeIC.py 1000 1
+ python makeIC.py 1000 0
fi
rm -rf Isothermal_*.hdf5
diff --git a/examples/Noh_3D/makeIC.py b/examples/Noh_3D/makeIC.py
index ec8d46639eecdefd55b917a955f13efc8ffe4126..0c25a5c8b3e967185cf16bae4b1f21c215266def 100644
--- a/examples/Noh_3D/makeIC.py
+++ b/examples/Noh_3D/makeIC.py
@@ -35,8 +35,8 @@ glass = h5py.File("glassCube_64.hdf5", "r")
vol = 8.
-pos = glass["/PartType0/Coordinates"][:,:] * cbrt(vol)
-h = glass["/PartType0/SmoothingLength"][:] * cbrt(vol)
+pos = glass["/PartType0/Coordinates"][:,:] * vol**(1./3.)
+h = glass["/PartType0/SmoothingLength"][:] * vol**(1./3.)
numPart = size(h)
# Generate extra arrays
@@ -65,7 +65,7 @@ file = h5py.File(fileName, 'w')
# Header
grp = file.create_group("/Header")
-grp.attrs["BoxSize"] = [cbrt(vol), cbrt(vol), cbrt(vol)]
+grp.attrs["BoxSize"] = [vol**(1./3.), vol**(1./3.), vol**(1./3.)]
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]
diff --git a/examples/SineWavePotential_1D/makeIC.py b/examples/SineWavePotential_1D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..321af0714219dfa0c7cbb3d80845881dcbb8416d
--- /dev/null
+++ b/examples/SineWavePotential_1D/makeIC.py
@@ -0,0 +1,99 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# 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/>.
+#
+##############################################################################
+
+# Generates a distorted 1D grid with a density profile that balances out the
+# external sine wave potential if run with an isothermal equation of state.
+
+import numpy as np
+import h5py
+
+# constant thermal energy
+# the definition below assumes the same thermal energy is defined in const.h,
+# and also that the code was configured with an adiabatic index of 5./3.
+uconst = 20.2615290634
+cs2 = 2.*uconst/3.
+A = 10.
+
+fileName = "sineWavePotential.hdf5"
+numPart = 100
+boxSize = 1.
+
+coords = np.zeros((numPart, 3))
+v = np.zeros((numPart, 3))
+m = np.zeros(numPart) + 1.
+h = np.zeros(numPart) + 2./numPart
+u = np.zeros(numPart) + uconst
+ids = np.arange(numPart, dtype = 'L')
+rho = np.zeros(numPart)
+
+# first set the positions, as we try to do a reasonable volume estimate to
+# set the masses
+for i in range(numPart):
+ coords[i,0] = (i+np.random.random())/numPart
+
+for i in range(numPart):
+ # reasonable mass estimate (actually not really good, but better than assuming
+ # a constant volume)
+ if i == 0:
+ V = 0.5*(coords[1,0]-coords[-1,0]+1.)
+ elif i == numPart-1:
+ V = 0.5*(coords[0,0]+1.-coords[-2,0])
+ else:
+ V = 0.5*(coords[i+1,0] - coords[i-1,0])
+ rho[i] = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*coords[i,0]))
+ m[i] = rho[i]*V
+
+#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
+grp.attrs["Dimension"] = 1
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
+
+#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")
+grp.create_dataset('Coordinates', data=coords, dtype='d')
+grp.create_dataset('Velocities', data=v, dtype='f')
+grp.create_dataset('Masses', data=m, dtype='f')
+grp.create_dataset('SmoothingLength', data=h, dtype='f')
+grp.create_dataset('InternalEnergy', data=u, dtype='f')
+grp.create_dataset('ParticleIDs', data=ids, dtype='L')
+grp.create_dataset('Density', data=rho, dtype='f')
+
+file.close()
diff --git a/examples/SineWavePotential_1D/plotSolution.py b/examples/SineWavePotential_1D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..65e981e4648fe0fe5d1da6cf3e753fb8a34f0fb4
--- /dev/null
+++ b/examples/SineWavePotential_1D/plotSolution.py
@@ -0,0 +1,77 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# 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/>.
+#
+##############################################################################
+
+# Plots some quantities for the snapshot file which is passed on as a command
+# line argument (full name)
+
+import numpy as np
+import h5py
+import sys
+import pylab as pl
+
+# these should be the same as in makeIC.py
+uconst = 20.2615290634
+cs2 = 2.*uconst/3.
+A = 10.
+
+if len(sys.argv) < 2:
+ print "Need to provide a filename argument!"
+ exit()
+
+fileName = sys.argv[1]
+
+file = h5py.File(fileName, 'r')
+coords = np.array(file["/PartType0/Coordinates"])
+rho = np.array(file["/PartType0/Density"])
+u = np.array(file["/PartType0/InternalEnergy"])
+agrav = np.array(file["/PartType0/GravAcceleration"])
+m = np.array(file["/PartType0/Masses"])
+ids = np.array(file["/PartType0/ParticleIDs"])
+
+x = np.linspace(0., 1., 1000)
+rho_x = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*x))
+
+P = cs2*rho
+
+ids_reverse = np.argsort(ids)
+
+gradP = np.zeros(P.shape)
+for i in range(len(P)):
+ iself = int(ids[i])
+ corr = 0.
+ im1 = iself-1
+ if im1 < 0:
+ im1 = len(P)-1
+ corr = 1.
+ ip1 = iself+1
+ if ip1 == len(P):
+ ip1 = 0
+ corr = 1.
+ idxp1 = ids_reverse[ip1]
+ idxm1 = ids_reverse[im1]
+ gradP[i] = (P[idxp1]-P[idxm1])/(coords[idxp1,0]-coords[idxm1,0])
+
+fig, ax = pl.subplots(2, 2)
+
+ax[0][0].plot(coords[:,0], rho, "r.", markersize = 0.5)
+ax[0][0].plot(x, rho_x, "g-")
+ax[0][1].plot(coords[:,0], gradP/rho, "b.")
+ax[1][0].plot(coords[:,0], agrav[:,0], "g.", markersize = 0.5)
+ax[1][1].plot(coords[:,0], m, "y.")
+pl.savefig("{fileName}.png".format(fileName = fileName[:-5]))
diff --git a/examples/SineWavePotential_1D/run.sh b/examples/SineWavePotential_1D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..077cf1c0cc64ef7a85cfd0e67f8f490b0de4ba37
--- /dev/null
+++ b/examples/SineWavePotential_1D/run.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+
+if [ ! -e sineWavePotential.hdf5 ]
+then
+ echo "Generating initial conditions for the 1D SineWavePotential example..."
+ python makeIC.py
+fi
+
+../swift -g -s -t 2 sineWavePotential.yml 2>&1 | tee output.log
+
+for f in sineWavePotential_*.hdf5
+do
+ python plotSolution.py $f
+done
diff --git a/examples/SineWavePotential_1D/sineWavePotential.yml b/examples/SineWavePotential_1D/sineWavePotential.yml
new file mode 100644
index 0000000000000000000000000000000000000000..1f98b47a0e3c72e77a2024aebf4dcd863ed02e9b
--- /dev/null
+++ b/examples/SineWavePotential_1D/sineWavePotential.yml
@@ -0,0 +1,39 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.
+ UnitLength_in_cgs: 1.
+ UnitVelocity_in_cgs: 1.
+ UnitCurrent_in_cgs: 1.
+ UnitTemp_in_cgs: 1.
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 10. # The end time of the simulation (in internal units).
+ dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1e-2 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-2 # Time between statistics output
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: sineWavePotential # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 1. # Time difference between consecutive outputs (in internal units)
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2349 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: sineWavePotential.hdf5 # The file to read
+
+# External potential parameters
+SineWavePotential:
+ amplitude: 10.
+ timestep_limit: 1.
diff --git a/examples/SineWavePotential_2D/makeIC.py b/examples/SineWavePotential_2D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..62ae89f8f52bff9c0db37cd537f286ab817da3fe
--- /dev/null
+++ b/examples/SineWavePotential_2D/makeIC.py
@@ -0,0 +1,95 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# 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/>.
+#
+##############################################################################
+
+# Generates a distorted 1D grid with a density profile that balances out the
+# external sine wave potential if run with an isothermal equation of state.
+
+import numpy as np
+import h5py
+
+# constant thermal energy
+# the definition below assumes the same thermal energy is defined in const.h,
+# and also that the code was configured with an adiabatic index of 5./3.
+uconst = 20.2615290634
+cs2 = 2.*uconst/3.
+A = 10.
+
+fileName = "sineWavePotential.hdf5"
+numPart_1D = 50
+boxSize = [1., 1.]
+numPart = numPart_1D**2
+
+coords = np.zeros((numPart, 3))
+v = np.zeros((numPart, 3))
+m = np.zeros(numPart) + 1.
+h = np.zeros(numPart) + 2./numPart
+u = np.zeros(numPart) + uconst
+ids = np.arange(numPart, dtype = 'L')
+rho = np.zeros(numPart)
+
+# first set the positions, as we try to do a reasonable volume estimate to
+# set the masses
+for i in range(numPart_1D):
+ for j in range(numPart_1D):
+ coords[numPart_1D*i+j,0] = (i+0.5)/numPart_1D
+ coords[numPart_1D*i+j,1] = (j+0.5)/numPart_1D
+
+V = 1./numPart
+for i in range(numPart):
+ rho[i] = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*coords[i,0]))
+ m[i] = rho[i]*V
+
+#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
+grp.attrs["Dimension"] = 2
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
+
+#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")
+grp.create_dataset('Coordinates', data=coords, dtype='d')
+grp.create_dataset('Velocities', data=v, dtype='f')
+grp.create_dataset('Masses', data=m, dtype='f')
+grp.create_dataset('SmoothingLength', data=h, dtype='f')
+grp.create_dataset('InternalEnergy', data=u, dtype='f')
+grp.create_dataset('ParticleIDs', data=ids, dtype='L')
+grp.create_dataset('Density', data=rho, dtype='f')
+
+file.close()
diff --git a/examples/SineWavePotential_2D/plotSolution.py b/examples/SineWavePotential_2D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..ee02f59c404db87a790465d2786e6296525e36b0
--- /dev/null
+++ b/examples/SineWavePotential_2D/plotSolution.py
@@ -0,0 +1,83 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# 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/>.
+#
+##############################################################################
+
+# Plots some quantities for the snapshot file which is passed on as a command
+# line argument (full name)
+
+import numpy as np
+import h5py
+import sys
+import pylab as pl
+
+# these should be the same as in makeIC.py
+uconst = 20.2615290634
+cs2 = 2.*uconst/3.
+A = 10.
+
+if len(sys.argv) < 2:
+ print "Need to provide a filename argument!"
+ exit()
+
+fileName = sys.argv[1]
+
+file = h5py.File(fileName, 'r')
+coords = np.array(file["/PartType0/Coordinates"])
+rho = np.array(file["/PartType0/Density"])
+u = np.array(file["/PartType0/InternalEnergy"])
+agrav = np.array(file["/PartType0/GravAcceleration"])
+m = np.array(file["/PartType0/Masses"])
+ids = np.array(file["/PartType0/ParticleIDs"])
+
+# ids_reverse gives the index original particle 0 now has in the particle arrays
+# and so on
+# note that this will only work if the particles do not move away too much from
+# there initial positions
+ids_reverse = np.argsort(ids)
+
+x = np.linspace(0., 1., 1000)
+rho_x = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*x))
+
+P = cs2*rho
+
+n1D = int(np.sqrt(len(P)))
+gradP = np.zeros(P.shape)
+for i in range(len(P)):
+ iself = int(ids[i]/n1D)
+ jself = int(ids[i]-n1D*iself)
+ corr = 0.
+ im1 = iself-1
+ if im1 < 0:
+ im1 = n1D-1
+ corr = 1.
+ ip1 = iself+1
+ if ip1 == n1D:
+ ip1 = 0
+ corr = 1.
+ idxp1 = ids_reverse[ip1*n1D+jself]
+ idxm1 = ids_reverse[im1*n1D+jself]
+ gradP[i] = (P[idxp1]-P[idxm1])/(coords[idxp1,0]-coords[idxm1,0]+corr)
+
+fig, ax = pl.subplots(2, 2)
+
+ax[0][0].plot(coords[:,0], rho, "r.", markersize = 0.5)
+ax[0][0].plot(x, rho_x, "g-")
+ax[0][1].plot(coords[:,0], gradP/rho, "b.")
+ax[1][0].plot(coords[:,0], agrav[:,0], "g.", markersize = 0.5)
+ax[1][1].plot(coords[:,0], m, "y.")
+pl.savefig("{fileName}.png".format(fileName = fileName[:-5]))
diff --git a/examples/SineWavePotential_2D/run.sh b/examples/SineWavePotential_2D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..077cf1c0cc64ef7a85cfd0e67f8f490b0de4ba37
--- /dev/null
+++ b/examples/SineWavePotential_2D/run.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+
+if [ ! -e sineWavePotential.hdf5 ]
+then
+ echo "Generating initial conditions for the 1D SineWavePotential example..."
+ python makeIC.py
+fi
+
+../swift -g -s -t 2 sineWavePotential.yml 2>&1 | tee output.log
+
+for f in sineWavePotential_*.hdf5
+do
+ python plotSolution.py $f
+done
diff --git a/examples/SineWavePotential_2D/sineWavePotential.yml b/examples/SineWavePotential_2D/sineWavePotential.yml
new file mode 100644
index 0000000000000000000000000000000000000000..f5a565ce4db46676c1f1f340f00025decea19644
--- /dev/null
+++ b/examples/SineWavePotential_2D/sineWavePotential.yml
@@ -0,0 +1,39 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.
+ UnitLength_in_cgs: 1.
+ UnitVelocity_in_cgs: 1.
+ UnitCurrent_in_cgs: 1.
+ UnitTemp_in_cgs: 1.
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 1. # The end time of the simulation (in internal units).
+ dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1.e-2 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-2 # Time between statistics output
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: sineWavePotential # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 0.1 # Time difference between consecutive outputs (in internal units)
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2349 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: sineWavePotential.hdf5 # The file to read
+
+# External potential parameters
+SineWavePotential:
+ amplitude: 10.
+ timestep_limit: 1.
diff --git a/examples/SineWavePotential_3D/makeIC.py b/examples/SineWavePotential_3D/makeIC.py
new file mode 100644
index 0000000000000000000000000000000000000000..4833ec1b055e27b63751136f0491e972fb9e492a
--- /dev/null
+++ b/examples/SineWavePotential_3D/makeIC.py
@@ -0,0 +1,106 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# 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/>.
+#
+##############################################################################
+
+# Generates a distorted 1D grid with a density profile that balances out the
+# external sine wave potential if run with an isothermal equation of state.
+
+import numpy as np
+import h5py
+
+# constant thermal energy
+# the definition below assumes the same thermal energy is defined in const.h,
+# and also that the code was configured with an adiabatic index of 5./3.
+uconst = 20.2615290634
+cs2 = 2.*uconst/3.
+A = 10.
+
+fileName = "sineWavePotential.hdf5"
+numPart_1D = 20
+boxSize = [1., 1.]
+numPart = numPart_1D**3
+
+coords = np.zeros((numPart, 3))
+v = np.zeros((numPart, 3))
+m = np.zeros(numPart) + 1.
+h = np.zeros(numPart) + 2./numPart
+u = np.zeros(numPart) + uconst
+ids = np.arange(numPart, dtype = 'L')
+rho = np.zeros(numPart)
+
+# first set the positions, as we try to do a reasonable volume estimate to
+# set the masses
+for i in range(numPart_1D):
+# coords[i,0] = (i+np.random.random())/numPart
+ for j in range(numPart_1D):
+ for k in range(numPart_1D):
+ coords[numPart_1D**2*i+numPart_1D*j+k,0] = (i+0.5)/numPart_1D
+ coords[numPart_1D**2*i+numPart_1D*j+k,1] = (j+0.5)/numPart_1D
+ coords[numPart_1D**2*i+numPart_1D*j+k,2] = (k+0.5)/numPart_1D
+
+V = 1./numPart
+for i in range(numPart):
+ # reasonable mass estimate (actually not really good, but better than assuming
+ # a constant volume)
+# if i == 0:
+# V = 0.5*(coords[1,0]-coords[-1,0]+1.)
+# elif i == numPart-1:
+# V = 0.5*(coords[0,0]+1.-coords[-2,0])
+# else:
+# V = 0.5*(coords[i+1,0] - coords[i-1,0])
+ rho[i] = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*coords[i,0]))
+ m[i] = rho[i]*V
+
+#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
+grp.attrs["Dimension"] = 3
+
+#Runtime parameters
+grp = file.create_group("/RuntimePars")
+grp.attrs["PeriodicBoundariesOn"] = 1
+
+#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")
+grp.create_dataset('Coordinates', data=coords, dtype='d')
+grp.create_dataset('Velocities', data=v, dtype='f')
+grp.create_dataset('Masses', data=m, dtype='f')
+grp.create_dataset('SmoothingLength', data=h, dtype='f')
+grp.create_dataset('InternalEnergy', data=u, dtype='f')
+grp.create_dataset('ParticleIDs', data=ids, dtype='L')
+grp.create_dataset('Density', data=rho, dtype='f')
+
+file.close()
diff --git a/examples/SineWavePotential_3D/plotSolution.py b/examples/SineWavePotential_3D/plotSolution.py
new file mode 100644
index 0000000000000000000000000000000000000000..13cae037b64eff4ad4fec0003bf0f5ad3ce94896
--- /dev/null
+++ b/examples/SineWavePotential_3D/plotSolution.py
@@ -0,0 +1,84 @@
+###############################################################################
+# This file is part of SWIFT.
+# Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+#
+# 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/>.
+#
+##############################################################################
+
+# Plots some quantities for the snapshot file which is passed on as a command
+# line argument (full name)
+
+import numpy as np
+import h5py
+import sys
+import pylab as pl
+
+# these should be the same as in makeIC.py
+uconst = 20.2615290634
+cs2 = 2.*uconst/3.
+A = 10.
+
+if len(sys.argv) < 2:
+ print "Need to provide a filename argument!"
+ exit()
+
+fileName = sys.argv[1]
+
+file = h5py.File(fileName, 'r')
+coords = np.array(file["/PartType0/Coordinates"])
+rho = np.array(file["/PartType0/Density"])
+u = np.array(file["/PartType0/InternalEnergy"])
+agrav = np.array(file["/PartType0/GravAcceleration"])
+m = np.array(file["/PartType0/Masses"])
+ids = np.array(file["/PartType0/ParticleIDs"])
+
+# ids_reverse gives the index original particle 0 now has in the particle arrays
+# and so on
+# note that this will only work if the particles do not move away too much from
+# there initial positions
+ids_reverse = np.argsort(ids)
+
+x = np.linspace(0., 1., 1000)
+rho_x = 1000.*np.exp(-0.5*A/np.pi/cs2*np.cos(2.*np.pi*x))
+
+P = cs2*rho
+
+n1D = np.ceil(len(P)**(1./3.))
+gradP = np.zeros(P.shape)
+for i in range(len(P)):
+ iself = int(ids[i]/n1D/n1D)
+ jself = int(int(ids[i]-n1D*iself)/n1D)
+ kself = int(ids[i]-n1D**2*iself-n1D*jself)
+ corr = 0.
+ im1 = iself-1
+ if im1 < 0:
+ im1 = n1D-1
+ corr = 1.
+ ip1 = iself+1
+ if ip1 == n1D:
+ ip1 = 0
+ corr = 1.
+ idxp1 = ids_reverse[ip1*n1D**2+jself*n1D+kself]
+ idxm1 = ids_reverse[im1*n1D**2+jself*n1D+kself]
+ gradP[i] = (P[idxp1]-P[idxm1])/(coords[idxp1,0]-coords[idxm1,0]+corr)
+
+fig, ax = pl.subplots(2, 2)
+
+ax[0][0].plot(coords[:,0], rho, "r.", markersize = 0.5)
+ax[0][0].plot(x, rho_x, "g-")
+ax[0][1].plot(coords[:,0], gradP/rho, "b.")
+ax[1][0].plot(coords[:,0], agrav[:,0], "g.", markersize = 0.5)
+ax[1][1].plot(coords[:,0], m, "y.")
+pl.savefig("{fileName}.png".format(fileName = fileName[:-5]))
diff --git a/examples/SineWavePotential_3D/run.sh b/examples/SineWavePotential_3D/run.sh
new file mode 100755
index 0000000000000000000000000000000000000000..077cf1c0cc64ef7a85cfd0e67f8f490b0de4ba37
--- /dev/null
+++ b/examples/SineWavePotential_3D/run.sh
@@ -0,0 +1,14 @@
+#!/bin/bash
+
+if [ ! -e sineWavePotential.hdf5 ]
+then
+ echo "Generating initial conditions for the 1D SineWavePotential example..."
+ python makeIC.py
+fi
+
+../swift -g -s -t 2 sineWavePotential.yml 2>&1 | tee output.log
+
+for f in sineWavePotential_*.hdf5
+do
+ python plotSolution.py $f
+done
diff --git a/examples/SineWavePotential_3D/sineWavePotential.yml b/examples/SineWavePotential_3D/sineWavePotential.yml
new file mode 100644
index 0000000000000000000000000000000000000000..391383568a5a94bd492cb228da4a7d4d24db413f
--- /dev/null
+++ b/examples/SineWavePotential_3D/sineWavePotential.yml
@@ -0,0 +1,39 @@
+# Define the system of units to use internally.
+InternalUnitSystem:
+ UnitMass_in_cgs: 1.
+ UnitLength_in_cgs: 1.
+ UnitVelocity_in_cgs: 1.
+ UnitCurrent_in_cgs: 1.
+ UnitTemp_in_cgs: 1.
+
+# Parameters governing the time integration
+TimeIntegration:
+ time_begin: 0. # The starting time of the simulation (in internal units).
+ time_end: 0.1 # The end time of the simulation (in internal units).
+ dt_min: 1e-6 # The minimal time-step size of the simulation (in internal units).
+ dt_max: 1.e-2 # The maximal time-step size of the simulation (in internal units).
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+ delta_time: 1e-2 # Time between statistics output
+
+# Parameters governing the snapshots
+Snapshots:
+ basename: sineWavePotential # Common part of the name of output files
+ time_first: 0. # Time of the first output (in internal units)
+ delta_time: 0.01 # Time difference between consecutive outputs (in internal units)
+
+# Parameters for the hydrodynamics scheme
+SPH:
+ resolution_eta: 1.2349 # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).
+ delta_neighbours: 0.1 # The tolerance for the targetted number of neighbours.
+ CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
+
+# Parameters related to the initial conditions
+InitialConditions:
+ file_name: sineWavePotential.hdf5 # The file to read
+
+# External potential parameters
+SineWavePotential:
+ amplitude: 10.
+ timestep_limit: 1.
diff --git a/examples/main.c b/examples/main.c
index 7f2331821271e203743f69154499fcc67db0384c..8b367cba14c96414d6eb5eb758a2e46f01a24774 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -57,9 +57,9 @@ void print_help_message() {
printf(" swift_mpi [OPTION]... PARAMFILE\n\n");
printf("Valid options are:\n");
- printf(" %2s %8s %s\n", "-a", "", "Pin runners using processor affinity");
- printf(" %2s %8s %s\n", "-c", "", "Run with cosmological time integration");
- printf(" %2s %8s %s\n", "-C", "", "Run with cooling");
+ printf(" %2s %8s %s\n", "-a", "", "Pin runners using processor affinity.");
+ printf(" %2s %8s %s\n", "-c", "", "Run with cosmological time integration.");
+ printf(" %2s %8s %s\n", "-C", "", "Run with cooling.");
printf(
" %2s %8s %s\n", "-d", "",
"Dry run. Read the parameter file, allocate memory but does not read ");
@@ -70,29 +70,32 @@ void print_help_message() {
"Allows user to check validy of parameter and IC files as well as "
"memory limits.");
printf(" %2s %8s %s\n", "-D", "",
- "Always drift all particles even the ones far from active particles.");
+ "Always drift all particles even the ones far from active particles. "
+ "This emulates");
+ printf(" %2s %8s %s\n", "", "",
+ "Gadget-[23] and GIZMO's default behaviours.");
printf(" %2s %8s %s\n", "-e", "",
- "Enable floating-point exceptions (debugging mode)");
+ "Enable floating-point exceptions (debugging mode).");
printf(" %2s %8s %s\n", "-f", "{int}",
- "Overwrite the CPU frequency (Hz) to be used for time measurements");
+ "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");
+ "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 "
"parameter to stop.");
- printf(" %2s %8s %s\n", "-s", "", "Run with SPH");
- printf(" %2s %8s %s\n", "-S", "", "Run with stars");
+ printf(" %2s %8s %s\n", "-s", "", "Run with hydrodynamics.");
+ printf(" %2s %8s %s\n", "-S", "", "Run with stars.");
printf(" %2s %8s %s\n", "-t", "{int}",
"The number of threads to use on each MPI rank. Defaults to 1 if not "
"specified.");
- printf(" %2s %8s %s\n", "-v", "[12]", "Increase the level of verbosity");
+ printf(" %2s %8s %s\n", "-v", "[12]", "Increase the level of verbosity.");
printf(" %2s %8s %s\n", "", "", "1: MPI-rank 0 writes ");
printf(" %2s %8s %s\n", "", "", "2: All MPI-ranks write");
printf(" %2s %8s %s\n", "-y", "{int}",
- "Time-step frequency at which task graphs are dumped");
- printf(" %2s %8s %s\n", "-h", "", "Print this help message and exit");
+ "Time-step frequency at which task graphs are dumped.");
+ printf(" %2s %8s %s\n", "-h", "", "Print this help message and exit.");
printf(
"\nSee the file parameter_example.yml for an example of "
"parameter file.\n");
@@ -302,6 +305,15 @@ int main(int argc, char *argv[]) {
message("WARNING: Debugging checks activated. Code will be slower !");
#endif
+/* Do we have gravity accuracy checks ? */
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+ if (myrank == 0)
+ message(
+ "WARNING: Checking 1/%d of all gpart for gravity accuracy. Code will "
+ "be slower !",
+ SWIFT_GRAVITY_FORCE_CHECKS);
+#endif
+
/* Do we choke on FP-exceptions ? */
if (with_fp_exceptions) {
feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
@@ -311,12 +323,14 @@ int main(int argc, char *argv[]) {
/* How large are the parts? */
if (myrank == 0) {
- message("sizeof(struct part) is %4zi bytes.", sizeof(struct part));
- message("sizeof(struct xpart) is %4zi bytes.", sizeof(struct xpart));
- message("sizeof(struct spart) is %4zi bytes.", sizeof(struct spart));
- message("sizeof(struct gpart) is %4zi bytes.", sizeof(struct gpart));
- message("sizeof(struct task) is %4zi bytes.", sizeof(struct task));
- message("sizeof(struct cell) is %4zi bytes.", sizeof(struct cell));
+ message("sizeof(part) is %4zi bytes.", sizeof(struct part));
+ message("sizeof(xpart) is %4zi bytes.", sizeof(struct xpart));
+ message("sizeof(spart) is %4zi bytes.", sizeof(struct spart));
+ message("sizeof(gpart) is %4zi bytes.", sizeof(struct gpart));
+ message("sizeof(multipole) is %4zi bytes.", sizeof(struct multipole));
+ message("sizeof(acc_tensor) is %4zi bytes.", sizeof(struct acc_tensor));
+ message("sizeof(task) is %4zi bytes.", sizeof(struct task));
+ message("sizeof(cell) is %4zi bytes.", sizeof(struct cell));
}
/* Read the parameter file */
@@ -588,7 +602,7 @@ int main(int argc, char *argv[]) {
clocks_getunit());
/* Main simulation loop */
- for (int j = 0; !engine_is_done(&e) && e.step != nsteps; j++) {
+ for (int j = 0; !engine_is_done(&e) && e.step - 1 != nsteps; j++) {
/* Reset timers */
timers_reset(timers_mask_all);
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index c6836713574da33056ac63558bf2d8518b07e418..bc893c6d27b67a2687e7766977c61f22df83d290 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -47,6 +47,7 @@ SPH:
CFL_condition: 0.1 # Courant-Friedrich-Levy condition for time integration.
max_ghost_iterations: 30 # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
max_volume_change: 2. # (Optional) Maximal allowed change of kernel volume over one time-step
+ h_max: 10. # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
# Parameters for the self-gravity scheme
Gravity:
@@ -55,7 +56,6 @@ Gravity:
a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
r_cut: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
-
# Parameters related to the initial conditions
InitialConditions:
file_name: SedovBlast/sedov.hdf5 # The file to read
diff --git a/src/Makefile.am b/src/Makefile.am
index 3e759d8fd15cfad6b38fe4ed469fb3b745aff786..29c110218416cc4bece0d766516eeb0a97fe4810 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -45,7 +45,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
physical_constants.h physical_constants_cgs.h potential.h version.h \
hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.h \
sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h profiler.h \
- dump.h logger.h active.h timeline.h xmf.h gravity_properties.h
+ dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
@@ -55,7 +55,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
physical_constants.c potential.c hydro_properties.c \
runner_doiact_fft.c threadpool.c cooling.c sourceterms.c \
statistics.c runner_doiact_vec.c profiler.c dump.c logger.c \
- part_type.c xmf.c gravity_properties.c
+ part_type.c xmf.c gravity_properties.c gravity.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 \
@@ -84,6 +84,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h
stars/Default/star_debug.h stars/Default/star_part.h \
potential/none/potential.h potential/point_mass/potential.h \
potential/isothermal/potential.h potential/disc_patch/potential.h \
+ potential/sine_wave/potential.h \
cooling/none/cooling.h cooling/none/cooling_struct.h \
cooling/const_du/cooling.h cooling/const_du/cooling_struct.h \
cooling/const_lambda/cooling.h cooling/const_lambda/cooling_struct.h \
diff --git a/src/active.h b/src/active.h
index 38020fde8ed53a638231097643476b7ef72d6d49..02e504f762735994e6c57f7e155071fede016713 100644
--- a/src/active.h
+++ b/src/active.h
@@ -105,10 +105,12 @@ __attribute__((always_inline)) INLINE static int cell_is_all_active(
__attribute__((always_inline)) INLINE static int part_is_active(
const struct part *p, const struct engine *e) {
- const integertime_t ti_current = e->ti_current;
- const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+ const timebin_t max_active_bin = e->max_active_bin;
+ const timebin_t part_bin = p->time_bin;
#ifdef SWIFT_DEBUG_CHECKS
+ const integertime_t ti_current = e->ti_current;
+ const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
if (ti_end < ti_current)
error(
"particle in an impossible time-zone! p->ti_end=%lld "
@@ -116,7 +118,7 @@ __attribute__((always_inline)) INLINE static int part_is_active(
ti_end, ti_current);
#endif
- return (ti_end == ti_current);
+ return (part_bin <= max_active_bin);
}
/**
@@ -129,10 +131,13 @@ __attribute__((always_inline)) INLINE static int part_is_active(
__attribute__((always_inline)) INLINE static int gpart_is_active(
const struct gpart *gp, const struct engine *e) {
+ const timebin_t max_active_bin = e->max_active_bin;
+ const timebin_t gpart_bin = gp->time_bin;
+
+#ifdef SWIFT_DEBUG_CHECKS
const integertime_t ti_current = e->ti_current;
const integertime_t ti_end = get_integer_time_end(ti_current, gp->time_bin);
-#ifdef SWIFT_DEBUG_CHECKS
if (ti_end < ti_current)
error(
"g-particle in an impossible time-zone! gp->ti_end=%lld "
@@ -140,7 +145,7 @@ __attribute__((always_inline)) INLINE static int gpart_is_active(
ti_end, ti_current);
#endif
- return (ti_end == ti_current);
+ return (gpart_bin <= max_active_bin);
}
/**
@@ -153,10 +158,13 @@ __attribute__((always_inline)) INLINE static int gpart_is_active(
__attribute__((always_inline)) INLINE static int spart_is_active(
const struct spart *sp, const struct engine *e) {
+ const timebin_t max_active_bin = e->max_active_bin;
+ const timebin_t spart_bin = sp->time_bin;
+
+#ifdef SWIFT_DEBUG_CHECKS
const integertime_t ti_current = e->ti_current;
const integertime_t ti_end = get_integer_time_end(ti_current, sp->time_bin);
-#ifdef SWIFT_DEBUG_CHECKS
if (ti_end < ti_current)
error(
"s-particle in an impossible time-zone! sp->ti_end=%lld "
@@ -164,7 +172,7 @@ __attribute__((always_inline)) INLINE static int spart_is_active(
ti_end, ti_current);
#endif
- return (ti_end == ti_current);
+ return (spart_bin <= max_active_bin);
}
/* Are cells / particles active for kick1 tasks ? */
@@ -201,11 +209,14 @@ __attribute__((always_inline)) INLINE static int cell_is_starting(
__attribute__((always_inline)) INLINE static int part_is_starting(
const struct part *p, const struct engine *e) {
+ const timebin_t max_active_bin = e->max_active_bin;
+ const timebin_t part_bin = p->time_bin;
+
+#ifdef SWIFT_DEBUG_CHECKS
const integertime_t ti_current = e->ti_current;
const integertime_t ti_beg =
get_integer_time_begin(ti_current + 1, p->time_bin);
-#ifdef SWIFT_DEBUG_CHECKS
if (ti_beg > ti_current)
error(
"particle in an impossible time-zone! p->ti_beg=%lld "
@@ -213,7 +224,7 @@ __attribute__((always_inline)) INLINE static int part_is_starting(
ti_beg, ti_current);
#endif
- return (ti_beg == ti_current);
+ return (part_bin <= max_active_bin);
}
/**
@@ -226,11 +237,14 @@ __attribute__((always_inline)) INLINE static int part_is_starting(
__attribute__((always_inline)) INLINE static int gpart_is_starting(
const struct gpart *gp, const struct engine *e) {
+ const timebin_t max_active_bin = e->max_active_bin;
+ const timebin_t gpart_bin = gp->time_bin;
+
+#ifdef SWIFT_DEBUG_CHECKS
const integertime_t ti_current = e->ti_current;
const integertime_t ti_beg =
get_integer_time_begin(ti_current + 1, gp->time_bin);
-#ifdef SWIFT_DEBUG_CHECKS
if (ti_beg > ti_current)
error(
"g-particle in an impossible time-zone! gp->ti_beg=%lld "
@@ -238,7 +252,7 @@ __attribute__((always_inline)) INLINE static int gpart_is_starting(
ti_beg, ti_current);
#endif
- return (ti_beg == ti_current);
+ return (gpart_bin <= max_active_bin);
}
/**
@@ -251,11 +265,14 @@ __attribute__((always_inline)) INLINE static int gpart_is_starting(
__attribute__((always_inline)) INLINE static int spart_is_starting(
const struct spart *sp, const struct engine *e) {
+ const timebin_t max_active_bin = e->max_active_bin;
+ const timebin_t spart_bin = sp->time_bin;
+
+#ifdef SWIFT_DEBUG_CHECKS
const integertime_t ti_current = e->ti_current;
const integertime_t ti_beg =
get_integer_time_begin(ti_current + 1, sp->time_bin);
-#ifdef SWIFT_DEBUG_CHECKS
if (ti_beg > ti_current)
error(
"s-particle in an impossible time-zone! sp->ti_beg=%lld "
@@ -263,6 +280,6 @@ __attribute__((always_inline)) INLINE static int spart_is_starting(
ti_beg, ti_current);
#endif
- return (ti_beg == ti_current);
+ return (spart_bin <= max_active_bin);
}
#endif /* SWIFT_ACTIVE_H */
diff --git a/src/cell.c b/src/cell.c
index 078310206b50ae65242a3a6b3dac738e2713e8d4..753bdd55061ebd2b2cdd2067691bd8319ca5a623 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -453,6 +453,70 @@ int cell_glocktree(struct cell *c) {
}
}
+/**
+ * @brief Lock a cell for access to its #multipole and hold its parents.
+ *
+ * @param c The #cell.
+ * @return 0 on success, 1 on failure
+ */
+int cell_mlocktree(struct cell *c) {
+
+ TIMER_TIC
+
+ /* First of all, try to lock this cell. */
+ if (c->mhold || lock_trylock(&c->mlock) != 0) {
+ TIMER_TOC(timer_locktree);
+ return 1;
+ }
+
+ /* Did somebody hold this cell in the meantime? */
+ if (c->mhold) {
+
+ /* Unlock this cell. */
+ if (lock_unlock(&c->mlock) != 0) error("Failed to unlock cell.");
+
+ /* Admit defeat. */
+ TIMER_TOC(timer_locktree);
+ return 1;
+ }
+
+ /* Climb up the tree and lock/hold/unlock. */
+ struct cell *finger;
+ for (finger = c->parent; finger != NULL; finger = finger->parent) {
+
+ /* Lock this cell. */
+ if (lock_trylock(&finger->mlock) != 0) break;
+
+ /* Increment the hold. */
+ atomic_inc(&finger->mhold);
+
+ /* Unlock the cell. */
+ if (lock_unlock(&finger->mlock) != 0) error("Failed to unlock cell.");
+ }
+
+ /* If we reached the top of the tree, we're done. */
+ if (finger == NULL) {
+ TIMER_TOC(timer_locktree);
+ return 0;
+ }
+
+ /* Otherwise, we hit a snag. */
+ else {
+
+ /* Undo the holds up to finger. */
+ for (struct cell *finger2 = c->parent; finger2 != finger;
+ finger2 = finger2->parent)
+ atomic_dec(&finger2->mhold);
+
+ /* Unlock this cell. */
+ if (lock_unlock(&c->mlock) != 0) error("Failed to unlock cell.");
+
+ /* Admit defeat. */
+ TIMER_TOC(timer_locktree);
+ return 1;
+ }
+}
+
/**
* @brief Lock a cell for access to its array of #spart and hold its parents.
*
@@ -555,6 +619,25 @@ void cell_gunlocktree(struct cell *c) {
TIMER_TOC(timer_locktree);
}
+/**
+ * @brief Unlock a cell's parents for access to its #multipole.
+ *
+ * @param c The #cell.
+ */
+void cell_munlocktree(struct cell *c) {
+
+ TIMER_TIC
+
+ /* First of all, try to unlock this cell. */
+ if (lock_unlock(&c->mlock) != 0) error("Failed to unlock cell.");
+
+ /* Climb up the tree and unhold the parents. */
+ for (struct cell *finger = c->parent; finger != NULL; finger = finger->parent)
+ atomic_dec(&finger->mhold);
+
+ TIMER_TOC(timer_locktree);
+}
+
/**
* @brief Unlock a cell's parents for access to #spart array.
*
@@ -973,6 +1056,23 @@ void cell_check_drift_point(struct cell *c, void *data) {
#endif
}
+/**
+ * @brief Resets all the individual cell task counters to 0.
+ *
+ * Should only be used for debugging purposes.
+ *
+ * @param c The #cell to reset.
+ */
+void cell_reset_task_counters(struct cell *c) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ for (int t = 0; t < task_type_count; ++t) c->tasks_executed[t] = 0;
+ for (int t = 0; t < task_subtype_count; ++t) c->subtasks_executed[t] = 0;
+#else
+ error("Calling debugging code without debugging flag activated.");
+#endif
+}
+
/**
* @brief Checks whether the cells are direct neighbours ot not. Both cells have
* to be of the same size
@@ -1012,52 +1112,34 @@ int cell_are_neighbours(const struct cell *restrict ci,
*/
void cell_check_multipole(struct cell *c, void *data) {
- struct multipole ma;
+#ifdef SWIFT_DEBUG_CHECKS
+ struct gravity_tensors ma;
+ const double tolerance = 1e-5; /* Relative */
+
+ /* First recurse */
+ if (c->split)
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) cell_check_multipole(c->progeny[k], NULL);
if (c->gcount > 0) {
/* Brute-force calculation */
- multipole_init(&ma, c->gparts, c->gcount);
-
- /* Compare with recursive one */
- struct multipole mb = c->multipole;
-
- if (fabsf(ma.mass - mb.mass) / fabsf(ma.mass + mb.mass) > 1e-5)
- error("Multipole masses are different (%12.15e vs. %12.15e)", ma.mass,
- mb.mass);
-
- for (int k = 0; k < 3; ++k)
- if (fabs(ma.CoM[k] - mb.CoM[k]) / fabs(ma.CoM[k] + mb.CoM[k]) > 1e-5)
- error("Multipole CoM are different (%12.15e vs. %12.15e", ma.CoM[k],
- mb.CoM[k]);
-
-#if const_gravity_multipole_order >= 2
- if (fabsf(ma.I_xx - mb.I_xx) / fabsf(ma.I_xx + mb.I_xx) > 1e-5 &&
- ma.I_xx > 1e-9)
- error("Multipole I_xx are different (%12.15e vs. %12.15e)", ma.I_xx,
- mb.I_xx);
- if (fabsf(ma.I_yy - mb.I_yy) / fabsf(ma.I_yy + mb.I_yy) > 1e-5 &&
- ma.I_yy > 1e-9)
- error("Multipole I_yy are different (%12.15e vs. %12.15e)", ma.I_yy,
- mb.I_yy);
- if (fabsf(ma.I_zz - mb.I_zz) / fabsf(ma.I_zz + mb.I_zz) > 1e-5 &&
- ma.I_zz > 1e-9)
- error("Multipole I_zz are different (%12.15e vs. %12.15e)", ma.I_zz,
- mb.I_zz);
- if (fabsf(ma.I_xy - mb.I_xy) / fabsf(ma.I_xy + mb.I_xy) > 1e-5 &&
- ma.I_xy > 1e-9)
- error("Multipole I_xy are different (%12.15e vs. %12.15e)", ma.I_xy,
- mb.I_xy);
- if (fabsf(ma.I_xz - mb.I_xz) / fabsf(ma.I_xz + mb.I_xz) > 1e-5 &&
- ma.I_xz > 1e-9)
- error("Multipole I_xz are different (%12.15e vs. %12.15e)", ma.I_xz,
- mb.I_xz);
- if (fabsf(ma.I_yz - mb.I_yz) / fabsf(ma.I_yz + mb.I_yz) > 1e-5 &&
- ma.I_yz > 1e-9)
- error("Multipole I_yz are different (%12.15e vs. %12.15e)", ma.I_yz,
- mb.I_yz);
-#endif
+ gravity_P2M(&ma, c->gparts, c->gcount);
+
+ /* Now compare the multipole expansion */
+ if (!gravity_multipole_equal(&ma.m_pole, &c->multipole->m_pole,
+ tolerance)) {
+ message("Multipoles are not equal at depth=%d!", c->depth);
+ message("Correct answer:");
+ gravity_multipole_print(&ma.m_pole);
+ message("Recursive multipole:");
+ gravity_multipole_print(&c->multipole->m_pole);
+ error("Aborting");
+ }
}
+#else
+ error("Calling debugging code without debugging flag activated.");
+#endif
}
/**
@@ -1239,6 +1321,9 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
if (c->kick2 != NULL) scheduler_activate(s, c->kick2);
if (c->timestep != NULL) scheduler_activate(s, c->timestep);
+ if (c->grav_down != NULL) scheduler_activate(s, c->grav_down);
+ if (c->grav_long_range != NULL) scheduler_activate(s, c->grav_long_range);
+ if (c->grav_top_level != NULL) scheduler_activate(s, c->grav_top_level);
if (c->cooling != NULL) scheduler_activate(s, c->cooling);
if (c->sourceterms != NULL) scheduler_activate(s, c->sourceterms);
@@ -1266,13 +1351,14 @@ void cell_set_super(struct cell *c, struct cell *super) {
}
/**
- * @brief Recursively drifts all particles and g-particles in a cell hierarchy.
+ * @brief Recursively drifts particles of all kinds in a cell hierarchy.
*
* @param c The #cell.
* @param e The #engine (to get ti_current).
*/
void cell_drift_particles(struct cell *c, const struct engine *e) {
+ const float hydro_h_max = e->hydro_properties->h_max;
const double timeBase = e->timeBase;
const integertime_t ti_old = c->ti_old;
const integertime_t ti_current = e->ti_current;
@@ -1283,7 +1369,7 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
/* Drift from the last time the cell was drifted to the current time */
const double dt = (ti_current - ti_old) * timeBase;
- float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
+ float dx_max = 0.f, dx2_max = 0.f, cell_h_max = 0.f;
/* Check that we are actually going to move forward. */
if (ti_current < ti_old) error("Attempt to drift to the past");
@@ -1297,7 +1383,7 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
struct cell *cp = c->progeny[k];
cell_drift_particles(cp, e);
dx_max = max(dx_max, cp->dx_max);
- h_max = max(h_max, cp->h_max);
+ cell_h_max = max(cell_h_max, cp->h_max);
}
} else if (ti_current > ti_old) {
@@ -1316,7 +1402,7 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
const float dx2 = gp->x_diff[0] * gp->x_diff[0] +
gp->x_diff[1] * gp->x_diff[1] +
gp->x_diff[2] * gp->x_diff[2];
- dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
+ dx2_max = max(dx2_max, dx2);
}
/* Loop over all the gas particles in the cell */
@@ -1330,14 +1416,17 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
/* Drift... */
drift_part(p, xp, dt, timeBase, ti_old, ti_current);
+ /* Limit h to within the allowed range */
+ p->h = min(p->h, hydro_h_max);
+
/* Compute (square of) motion since last cell construction */
const float dx2 = xp->x_diff[0] * xp->x_diff[0] +
xp->x_diff[1] * xp->x_diff[1] +
xp->x_diff[2] * xp->x_diff[2];
- dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
+ dx2_max = max(dx2_max, dx2);
/* Maximal smoothing length */
- h_max = (h_max > p->h) ? h_max : p->h;
+ cell_h_max = max(cell_h_max, p->h);
}
/* Loop over all the star particles in the cell */
@@ -1358,18 +1447,66 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
} else {
- h_max = c->h_max;
+ cell_h_max = c->h_max;
dx_max = c->dx_max;
}
/* Store the values */
- c->h_max = h_max;
+ c->h_max = cell_h_max;
c->dx_max = dx_max;
/* Update the time of the last drift */
c->ti_old = ti_current;
}
+/**
+ * @brief Recursively drifts all multipoles in a cell hierarchy.
+ *
+ * @param c The #cell.
+ * @param e The #engine (to get ti_current).
+ */
+void cell_drift_all_multipoles(struct cell *c, const struct engine *e) {
+
+ const double timeBase = e->timeBase;
+ const integertime_t ti_old_multipole = c->ti_old_multipole;
+ const integertime_t ti_current = e->ti_current;
+
+ /* Drift from the last time the cell was drifted to the current time */
+ const double dt = (ti_current - ti_old_multipole) * timeBase;
+
+ /* Check that we are actually going to move forward. */
+ if (ti_current < ti_old_multipole) error("Attempt to drift to the past");
+
+ /* Are we not in a leaf ? */
+ if (c->split) {
+
+ /* Loop over the progeny and drift the multipoles. */
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) cell_drift_particles(c->progeny[k], e);
+
+ } else if (ti_current > ti_old_multipole) {
+
+ /* Drift the multipole */
+ gravity_multipole_drift(c->multipole, dt);
+ }
+
+ /* Update the time of the last drift */
+ c->ti_old_multipole = ti_current;
+}
+
+/**
+ * @brief Drifts the multipole of a cell to the current time.
+ *
+ * Only drifts the multipole at this level. Multipoles deeper in the
+ * tree are not updated.
+ *
+ * @param c The #cell.
+ * @param e The #engine (to get ti_current).
+ */
+void cell_drift_multipole(struct cell *c, const struct engine *e) {
+ error("To be implemented");
+}
+
/**
* @brief Recursively checks that all particles in a cell have a time-step
*/
diff --git a/src/cell.h b/src/cell.h
index c97449f9151b8954ee54cbd6bff5d8ab73cbf53b..b9fbb6abffe427b81daec033447fbf35804c242d 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -94,9 +94,6 @@ struct pcell {
*/
struct cell {
- /*! This cell's multipole. */
- struct multipole multipole;
-
/*! The cell location on the grid. */
double loc[3];
@@ -106,6 +103,9 @@ struct cell {
/*! Max smoothing length in this cell. */
double h_max;
+ /*! This cell's multipole. */
+ struct gravity_tensors *multipole;
+
/*! Linking pointer for "memory management". */
struct cell *next;
@@ -170,7 +170,10 @@ struct cell {
struct task *timestep;
/*! Task constructing the multipole from the particles */
- struct task *grav_up;
+ struct task *grav_top_level;
+
+ /*! Task constructing the multipole from the particles */
+ struct task *grav_long_range;
/*! Task propagating the multipole to the particles */
struct task *grav_down;
@@ -230,9 +233,12 @@ struct cell {
/*! Maximum beginning of (integer) time step in this cell. */
integertime_t ti_beg_max;
- /*! Last (integer) time the cell's content was drifted forward in time. */
+ /*! Last (integer) time the cell's particle was drifted forward in time. */
integertime_t ti_old;
+ /*! Last (integer) time the cell's multipole was drifted forward in time. */
+ integertime_t ti_old_multipole;
+
/*! Minimum dimension, i.e. smallest edge of this cell (min(width)). */
float dmin;
@@ -260,6 +266,9 @@ struct cell {
/*! Spin lock for various uses (#gpart case). */
swift_lock_type glock;
+ /*! Spin lock for various uses (#multipole case). */
+ swift_lock_type mlock;
+
/*! Spin lock for various uses (#spart case). */
swift_lock_type slock;
@@ -284,6 +293,9 @@ struct cell {
/*! Is the #gpart data of this cell being used in a sub-cell? */
int ghold;
+ /*! Is the #multipole data of this cell being used in a sub-cell? */
+ int mhold;
+
/*! Is the #spart data of this cell being used in a sub-cell? */
int shold;
@@ -299,6 +311,14 @@ struct cell {
/*! The maximal depth of this cell and its progenies */
char maxdepth;
+#ifdef SWIFT_DEBUG_CHECKS
+ /*! The list of tasks that have been executed on this cell */
+ char tasks_executed[64];
+
+ /*! The list of sub-tasks that have been executed on this cell */
+ char subtasks_executed[64];
+#endif
+
} SWIFT_STRUCT_ALIGN;
/* Convert cell location to ID. */
@@ -314,6 +334,8 @@ int cell_locktree(struct cell *c);
void cell_unlocktree(struct cell *c);
int cell_glocktree(struct cell *c);
void cell_gunlocktree(struct cell *c);
+int cell_mlocktree(struct cell *c);
+void cell_munlocktree(struct cell *c);
int cell_slocktree(struct cell *c);
void cell_sunlocktree(struct cell *c);
int cell_pack(struct cell *c, struct pcell *pc);
@@ -331,10 +353,13 @@ int cell_are_neighbours(const struct cell *restrict ci,
void cell_check_multipole(struct cell *c, void *data);
void cell_clean(struct cell *c);
void cell_check_drift_point(struct cell *c, void *data);
+void cell_reset_task_counters(struct cell *c);
int cell_is_drift_needed(struct cell *c, const struct engine *e);
int cell_unskip_tasks(struct cell *c, struct scheduler *s);
void cell_set_super(struct cell *c, struct cell *super);
void cell_drift_particles(struct cell *c, const struct engine *e);
+void cell_drift_multipole(struct cell *c, const struct engine *e);
+void cell_drift_all_multipoles(struct cell *c, const struct engine *e);
void cell_check_timesteps(struct cell *c);
#endif /* SWIFT_CELL_H */
diff --git a/src/const.h b/src/const.h
index c290a3e73a524e9003cadb63f8bd45e8b3c51dac..88c1a1af1cfc36cc401fdfea0b077f79fcd13bc0 100644
--- a/src/const.h
+++ b/src/const.h
@@ -41,9 +41,6 @@
/* Self gravity stuff. */
#define const_gravity_multipole_order 1
-#define const_gravity_a_smooth 1.25f
-#define const_gravity_r_cut 4.5f
-#define const_gravity_eta 0.025f
/* Type of gradients to use (GIZMO_SPH only) */
/* If no option is chosen, no gradients are used (first order scheme) */
@@ -55,6 +52,11 @@
#define SLOPE_LIMITER_PER_FACE
#define SLOPE_LIMITER_CELL_WIDE
+/* Options to control the movement of particles for GIZMO_SPH. */
+/* This option disables particle movement */
+//#define GIZMO_FIX_PARTICLES
+//#define GIZMO_TOTAL_ENERGY
+
/* Source terms */
#define SOURCETERMS_NONE
//#define SOURCETERMS_SN_FEEDBACK
diff --git a/src/engine.c b/src/engine.c
index 4fb988e41c66413470b4b20f3c0f2d71589ff529..9b5f771f80ff269dfd39d93ad59e5d81d613c89c 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -55,6 +55,7 @@
#include "cycle.h"
#include "debug.h"
#include "error.h"
+#include "gravity.h"
#include "hydro.h"
#include "minmax.h"
#include "parallel_io.h"
@@ -132,6 +133,7 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
struct scheduler *s = &e->sched;
const int is_hydro = (e->policy & engine_policy_hydro);
+ const int is_self_gravity = (e->policy & engine_policy_self_gravity);
const int is_with_cooling = (e->policy & engine_policy_cooling);
const int is_with_sourceterms = (e->policy & engine_policy_sourceterms);
@@ -165,6 +167,27 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
scheduler_addunlock(s, c->drift, c->init);
+ if (is_self_gravity) {
+
+ /* Gravity non-neighbouring pm calculations */
+ c->grav_long_range = scheduler_addtask(
+ s, task_type_grav_long_range, task_subtype_none, 0, 0, c, NULL, 0);
+
+ /* Gravity top-level periodic calculation */
+ c->grav_top_level = scheduler_addtask(
+ s, task_type_grav_top_level, task_subtype_none, 0, 0, c, NULL, 0);
+
+ /* Gravity recursive down-pass */
+ c->grav_down = scheduler_addtask(s, task_type_grav_down,
+ task_subtype_none, 0, 0, c, NULL, 0);
+
+ scheduler_addunlock(s, c->init, c->grav_long_range);
+ scheduler_addunlock(s, c->init, c->grav_top_level);
+ scheduler_addunlock(s, c->grav_long_range, c->grav_down);
+ scheduler_addunlock(s, c->grav_top_level, c->grav_down);
+ scheduler_addunlock(s, c->grav_down, c->kick2);
+ }
+
/* Generate the ghost task. */
if (is_hydro)
c->ghost = scheduler_addtask(s, task_type_ghost, task_subtype_none, 0,
@@ -871,7 +894,7 @@ void engine_addtasks_grav(struct engine *e, struct cell *c, struct task *up,
struct task *down) {
/* Link the tasks to this cell. */
- c->grav_up = up;
+ // c->grav_up = up;
c->grav_down = down;
/* Recurse? */
@@ -1550,7 +1573,7 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
*
* @param e The #engine.
*/
-void engine_make_gravity_tasks(struct engine *e) {
+void engine_make_self_gravity_tasks(struct engine *e) {
struct space *s = e->s;
struct scheduler *sched = &e->sched;
@@ -1572,10 +1595,6 @@ void engine_make_gravity_tasks(struct engine *e) {
scheduler_addtask(sched, task_type_self, task_subtype_grav, 0, 0, ci, NULL,
0);
- /* Let's also build a task for all the non-neighbouring pm calculations */
- scheduler_addtask(sched, task_type_grav_mm, task_subtype_none, 0, 0, ci,
- NULL, 0);
-
for (int cjd = cid + 1; cjd < nr_cells; ++cjd) {
struct cell *cj = &cells[cjd];
@@ -1784,16 +1803,12 @@ void engine_count_and_link_tasks(struct engine *e) {
* @param gravity The gravity task to link.
* @param c The cell.
*/
-static inline void engine_make_gravity_dependencies(struct scheduler *sched,
- struct task *gravity,
- struct cell *c) {
+static inline void engine_make_self_gravity_dependencies(
+ struct scheduler *sched, struct task *gravity, struct cell *c) {
- /* init --> gravity --> kick */
+ /* init --> gravity --> grav_down --> kick */
scheduler_addunlock(sched, c->super->init, gravity);
- scheduler_addunlock(sched, gravity, c->super->kick2);
-
- /* grav_up --> gravity ( --> kick) */
- scheduler_addunlock(sched, c->super->grav_up, gravity);
+ scheduler_addunlock(sched, gravity, c->super->grav_down);
}
/**
@@ -1823,47 +1838,19 @@ void engine_link_gravity_tasks(struct engine *e) {
const int nodeID = e->nodeID;
const int nr_tasks = sched->nr_tasks;
- /* Add one task gathering all the multipoles */
- struct task *gather = scheduler_addtask(
- sched, task_type_grav_gather_m, task_subtype_none, 0, 0, NULL, NULL, 0);
-
- /* And one task performing the FFT */
- struct task *fft = scheduler_addtask(sched, task_type_grav_fft,
- task_subtype_none, 0, 0, NULL, NULL, 0);
-
- scheduler_addunlock(sched, gather, fft);
-
for (int k = 0; k < nr_tasks; k++) {
/* Get a pointer to the task. */
struct task *t = &sched->tasks[k];
- /* Multipole construction */
- if (t->type == task_type_grav_up) {
- scheduler_addunlock(sched, t, gather);
- }
-
- /* Long-range interaction */
- if (t->type == task_type_grav_mm) {
-
- /* Gather the multipoles --> mm interaction --> kick */
- scheduler_addunlock(sched, gather, t);
- scheduler_addunlock(sched, t, t->ci->super->kick2);
-
- /* init --> mm interaction */
- scheduler_addunlock(sched, t->ci->super->init, t);
- }
-
/* Self-interaction for self-gravity? */
if (t->type == task_type_self && t->subtype == task_subtype_grav) {
- engine_make_gravity_dependencies(sched, t, t->ci);
-
+ engine_make_self_gravity_dependencies(sched, t, t->ci);
}
/* Self-interaction for external gravity ? */
- else if (t->type == task_type_self &&
- t->subtype == task_subtype_external_grav) {
+ if (t->type == task_type_self && t->subtype == task_subtype_external_grav) {
engine_make_external_gravity_dependencies(sched, t, t->ci);
@@ -1874,12 +1861,12 @@ void engine_link_gravity_tasks(struct engine *e) {
if (t->ci->nodeID == nodeID) {
- engine_make_gravity_dependencies(sched, t, t->ci);
+ engine_make_self_gravity_dependencies(sched, t, t->ci);
}
if (t->cj->nodeID == nodeID && t->ci->super != t->cj->super) {
- engine_make_gravity_dependencies(sched, t, t->cj);
+ engine_make_self_gravity_dependencies(sched, t, t->cj);
}
}
@@ -1888,7 +1875,7 @@ void engine_link_gravity_tasks(struct engine *e) {
else if (t->type == task_type_sub_self && t->subtype == task_subtype_grav) {
if (t->ci->nodeID == nodeID) {
- engine_make_gravity_dependencies(sched, t, t->ci);
+ engine_make_self_gravity_dependencies(sched, t, t->ci);
}
}
@@ -1906,11 +1893,11 @@ void engine_link_gravity_tasks(struct engine *e) {
if (t->ci->nodeID == nodeID) {
- engine_make_gravity_dependencies(sched, t, t->ci);
+ engine_make_self_gravity_dependencies(sched, t, t->ci);
}
if (t->cj->nodeID == nodeID && t->ci->super != t->cj->super) {
- engine_make_gravity_dependencies(sched, t, t->cj);
+ engine_make_self_gravity_dependencies(sched, t, t->cj);
}
}
}
@@ -1926,10 +1913,12 @@ void engine_link_gravity_tasks(struct engine *e) {
* @param gradient The gradient task to link.
* @param force The force task to link.
* @param c The cell.
+ * @param with_cooling Do we have a cooling task ?
*/
static inline void engine_make_hydro_loops_dependencies(
struct scheduler *sched, struct task *density, struct task *gradient,
struct task *force, struct cell *c, int with_cooling) {
+
/* init --> density loop --> ghost --> gradient loop --> extra_ghost */
/* extra_ghost --> force loop */
scheduler_addunlock(sched, c->super->init, density);
@@ -2188,32 +2177,28 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
*/
void engine_make_gravityrecursive_tasks(struct engine *e) {
- struct space *s = e->s;
- struct scheduler *sched = &e->sched;
- const int nodeID = e->nodeID;
- const int nr_cells = s->nr_cells;
- struct cell *cells = s->cells_top;
-
- for (int k = 0; k < nr_cells; k++) {
+ /* struct space *s = e->s; */
+ /* struct scheduler *sched = &e->sched; */
+ /* const int nodeID = e->nodeID; */
+ /* const int nr_cells = s->nr_cells; */
+ /* struct cell *cells = s->cells_top; */
- /* Only do this for local cells containing gravity particles */
- if (cells[k].nodeID == nodeID && cells[k].gcount > 0) {
+ /* for (int k = 0; k < nr_cells; k++) { */
- /* Create tasks at top level. */
- struct task *up =
- scheduler_addtask(sched, task_type_grav_up, task_subtype_none, 0, 0,
- &cells[k], NULL, 0);
+ /* /\* Only do this for local cells containing gravity particles *\/ */
+ /* if (cells[k].nodeID == nodeID && cells[k].gcount > 0) { */
- struct task *down = NULL;
- /* struct task *down = */
- /* scheduler_addtask(sched, task_type_grav_down, task_subtype_none, 0,
- * 0, */
- /* &cells[k], NULL, 0); */
+ /* /\* Create tasks at top level. *\/ */
+ /* struct task *up = NULL; */
+ /* struct task *down = NULL; */
+ /* /\* scheduler_addtask(sched, task_type_grav_down,
+ * task_subtype_none, 0, 0, *\/ */
+ /* /\* &cells[k], NULL, 0); *\/ */
- /* Push tasks down the cell hierarchy. */
- engine_addtasks_grav(e, &cells[k], up, down);
- }
- }
+ /* /\* Push tasks down the cell hierarchy. *\/ */
+ /* engine_addtasks_grav(e, &cells[k], up, down); */
+ /* } */
+ /* } */
}
/**
@@ -2235,8 +2220,8 @@ void engine_maketasks(struct engine *e) {
/* Construct the firt hydro loop over neighbours */
if (e->policy & engine_policy_hydro) engine_make_hydroloop_tasks(e);
- /* Add the gravity mm tasks. */
- if (e->policy & engine_policy_self_gravity) engine_make_gravity_tasks(e);
+ /* Add the self gravity tasks. */
+ if (e->policy & engine_policy_self_gravity) engine_make_self_gravity_tasks(e);
/* Add the external gravity tasks. */
if (e->policy & engine_policy_external_gravity)
@@ -2479,6 +2464,13 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
if (cell_is_active(t->ci, e)) scheduler_activate(s, t);
}
+ /* Gravity ? */
+ else if (t->type == task_type_grav_down ||
+ t->type == task_type_grav_long_range ||
+ t->type == task_type_grav_top_level) {
+ if (cell_is_active(t->ci, e)) scheduler_activate(s, t);
+ }
+
/* Time-step? */
else if (t->type == task_type_timestep) {
t->ci->updated = 0;
@@ -2486,12 +2478,6 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
t->ci->s_updated = 0;
if (cell_is_active(t->ci, e)) scheduler_activate(s, t);
}
-
- /* Tasks with no cells should not be skipped? */
- else if (t->type == task_type_grav_gather_m ||
- t->type == task_type_grav_fft) {
- scheduler_activate(s, t);
- }
}
}
@@ -2555,6 +2541,7 @@ void engine_print_task_counts(struct engine *e) {
fflush(stdout);
message("nr_parts = %zu.", e->s->nr_parts);
message("nr_gparts = %zu.", e->s->nr_gparts);
+ message("nr_sparts = %zu.", e->s->nr_sparts);
if (e->verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
@@ -2817,6 +2804,20 @@ void engine_print_stats(struct engine *e) {
const ticks tic = getticks();
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that all cells have been drifted to the current time.
+ * That can include cells that have not
+ * previously been active on this rank. */
+ space_check_drift_point(e->s, e->ti_current);
+
+ /* Be verbose about this */
+ message("Saving statistics at t=%e.", e->time);
+#else
+ if (e->verbose) message("Saving statistics at t=%e.", e->time);
+#endif
+
+ e->save_stats = 0;
+
struct statistics stats;
stats_init(&stats);
@@ -2861,8 +2862,10 @@ void engine_skip_force_and_kick(struct engine *e) {
/* Skip everything that updates the particles */
if (t->type == task_type_drift || t->type == task_type_kick1 ||
t->type == task_type_kick2 || t->type == task_type_timestep ||
- t->subtype == task_subtype_force || t->type == task_type_cooling ||
- t->type == task_type_sourceterms)
+ t->subtype == task_subtype_force || t->subtype == task_subtype_grav ||
+ t->type == task_type_grav_long_range ||
+ t->type == task_type_grav_top_level || t->type == task_type_grav_down ||
+ t->type == task_type_cooling || t->type == task_type_sourceterms)
t->skip = 1;
}
}
@@ -2896,6 +2899,11 @@ void engine_launch(struct engine *e, int nr_runners) {
const ticks tic = getticks();
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Re-set all the cell task counters to 0 */
+ space_reset_task_counters(e->s);
+#endif
+
/* Prepare the scheduler. */
atomic_inc(&e->sched.waiting);
@@ -2974,6 +2982,19 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
}
}
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Let's store the total mass in the system for future checks */
+ e->s->total_mass = 0.;
+ for (size_t i = 0; i < s->nr_gparts; ++i)
+ e->s->total_mass += s->gparts[i].mass;
+#ifdef WITH_MPI
+ if (MPI_Allreduce(MPI_IN_PLACE, &e->s->total_mass, 1, MPI_DOUBLE, MPI_SUM,
+ MPI_COMM_WORLD) != MPI_SUCCESS)
+ error("Failed to all-reduce total mass in the system.");
+#endif
+ message("Total mass in the system: %e", e->s->total_mass);
+#endif
+
/* Now time to get ready for the first time-step */
if (e->nodeID == 0) message("Running initial fake time-step.");
@@ -3001,7 +3022,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
#endif
/* Ready to go */
- e->step = -1;
+ e->step = 0;
e->forcerebuild = 1;
e->wallclock_time = (float)clocks_diff(&time1, &time2);
@@ -3024,14 +3045,6 @@ void engine_step(struct engine *e) {
e->tic_step = getticks();
#endif
- /* Move forward in time */
- e->ti_old = e->ti_current;
- e->ti_current = e->ti_end_min;
- e->step += 1;
- e->time = e->ti_current * e->timeBase + e->timeBegin;
- e->timeOld = e->ti_old * e->timeBase + e->timeBegin;
- e->timeStep = (e->ti_current - e->ti_old) * e->timeBase;
-
if (e->nodeID == 0) {
/* Print some information to the screen */
@@ -3048,6 +3061,18 @@ void engine_step(struct engine *e) {
fflush(e->file_timesteps);
}
+ /* Move forward in time */
+ e->ti_old = e->ti_current;
+ e->ti_current = e->ti_end_min;
+ e->max_active_bin = get_max_active_bin(e->ti_end_min);
+ e->step += 1;
+ e->time = e->ti_current * e->timeBase + e->timeBegin;
+ e->timeOld = e->ti_old * e->timeBase + e->timeBegin;
+ e->timeStep = (e->ti_current - e->ti_old) * e->timeBase;
+
+ /* Prepare the tasks to be launched, rebuild or repartition if needed. */
+ engine_prepare(e);
+
/* Repartition the space amongst the nodes? */
#ifdef WITH_MPI
@@ -3138,24 +3163,27 @@ void engine_step(struct engine *e) {
e->lastrepart = e->forcerepart;
#endif
- /* Prepare the tasks to be launched, rebuild or repartition if needed. */
- e->lastrebuild = e->forcerebuild;
- engine_prepare(e);
+ /* Are we drifting everything (a la Gadget/GIZMO) ? */
+ if (e->policy & engine_policy_drift_all) engine_drift_all(e);
/* Print the number of active tasks ? */
if (e->verbose) engine_print_task_counts(e);
- /* Save some statistics ? */
- if (e->time - e->timeLastStatistics >= e->deltaTimeStatistics) {
- engine_print_stats(e);
- e->timeLastStatistics += e->deltaTimeStatistics;
- }
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+ /* Run the brute-force gravity calculation for some gparts */
+ gravity_exact_force_compute(e->s, e);
+#endif
/* Start all the tasks. */
TIMER_TIC;
engine_launch(e, e->nr_threads);
TIMER_TOC(timer_runners);
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+ /* Check the accuracy of the gravity calculation */
+ gravity_exact_force_check(e->s, e, 1e-1);
+#endif
+
/* Collect the values of rebuild from all nodes. */
#ifdef WITH_MPI
int buff = 0;
@@ -3165,13 +3193,18 @@ void engine_step(struct engine *e) {
e->forcerebuild = buff;
#endif
+ /* Save some statistics ? */
+ if (e->time - e->timeLastStatistics >= e->deltaTimeStatistics) {
+ e->save_stats = 1;
+ }
+
/* Do we want a snapshot? */
if (e->ti_end_min >= e->ti_nextSnapshot && e->ti_nextSnapshot > 0)
e->dump_snapshot = 1;
/* Drift everybody (i.e. what has not yet been drifted) */
/* to the current time */
- if (e->dump_snapshot || e->forcerebuild || e->forcerepart)
+ if (e->dump_snapshot || e->forcerebuild || e->forcerepart || e->save_stats)
engine_drift_all(e);
/* Write a snapshot ? */
@@ -3184,6 +3217,16 @@ void engine_step(struct engine *e) {
engine_compute_next_snapshot_time(e);
}
+ /* Save some statistics */
+ if (e->save_stats) {
+
+ /* Dump */
+ engine_print_stats(e);
+
+ /* and move on */
+ e->timeLastStatistics += e->deltaTimeStatistics;
+ }
+
/* Recover the (integer) end of the next time-step */
engine_collect_timestep(e);
@@ -3226,6 +3269,33 @@ void engine_unskip(struct engine *e) {
clocks_getunit());
}
+/**
+ * @brief Mapper function to drift ALL particle types and multipoles forward in
+ * time.
+ *
+ * @param map_data An array of #cell%s.
+ * @param num_elements Chunk size.
+ * @param extra_data Pointer to an #engine.
+ */
+void engine_do_drift_all_mapper(void *map_data, int num_elements,
+ void *extra_data) {
+
+ struct engine *e = (struct engine *)extra_data;
+ struct cell *cells = (struct cell *)map_data;
+
+ for (int ind = 0; ind < num_elements; ind++) {
+ struct cell *c = &cells[ind];
+ if (c != NULL && c->nodeID == e->nodeID) {
+ /* Drift all the particles */
+ cell_drift_particles(c, e);
+
+ /* Drift the multipole */
+ if (e->policy & engine_policy_self_gravity)
+ cell_drift_all_multipoles(c, e);
+ }
+ }
+}
+
/**
* @brief Drift *all* particles forward to the current time.
*
@@ -3234,7 +3304,12 @@ void engine_unskip(struct engine *e) {
void engine_drift_all(struct engine *e) {
const ticks tic = getticks();
- threadpool_map(&e->threadpool, runner_do_drift_mapper, e->s->cells_top,
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (e->nodeID == 0) message("Drifting all");
+#endif
+
+ threadpool_map(&e->threadpool, engine_do_drift_all_mapper, e->s->cells_top,
e->s->nr_cells, sizeof(struct cell), 1, e);
#ifdef SWIFT_DEBUG_CHECKS
@@ -3456,9 +3531,12 @@ void engine_dump_snapshot(struct engine *e) {
* That can include cells that have not
* previously been active on this rank. */
space_check_drift_point(e->s, e->ti_current);
-#endif
+ /* Be verbose about this */
+ message("writing snapshot at t=%e.", e->time);
+#else
if (e->verbose) message("writing snapshot at t=%e.", e->time);
+#endif
/* Dump... */
#if defined(WITH_MPI)
@@ -3595,6 +3673,7 @@ void engine_init(struct engine *e, struct space *s,
e->lastrepart = 0;
e->reparttype = reparttype;
e->dump_snapshot = 0;
+ e->save_stats = 0;
e->links = NULL;
e->nr_links = 0;
e->timeBegin = parser_get_param_double(params, "TimeIntegration:time_begin");
@@ -3603,6 +3682,7 @@ void engine_init(struct engine *e, struct space *s,
e->time = e->timeBegin;
e->ti_old = 0;
e->ti_current = 0;
+ e->max_active_bin = num_time_bins;
e->timeStep = 0.;
e->timeBase = 0.;
e->timeBase_inv = 0.;
diff --git a/src/engine.h b/src/engine.h
index bbb3ae2aa96a2a61c4b52fe41b084ad21819c5a4..47940bb7cce49f5d27afa5b380a63b3cebbad7b8 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -123,6 +123,9 @@ struct engine {
double time;
integertime_t ti_current;
+ /* The highest active bin at this time */
+ timebin_t max_active_bin;
+
/* Time step */
double timeStep;
@@ -207,6 +210,9 @@ struct engine {
int lastrepart;
struct repartition *reparttype;
+ /* Need to dump some statistics ? */
+ int save_stats;
+
/* Need to dump a snapshot ? */
int dump_snapshot;
diff --git a/src/gravity.c b/src/gravity.c
new file mode 100644
index 0000000000000000000000000000000000000000..86f9fa82e3eb693cc3c051420fb9c7bff277eb9f
--- /dev/null
+++ b/src/gravity.c
@@ -0,0 +1,187 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* This object's header. */
+#include "gravity.h"
+
+/* Local headers. */
+#include "active.h"
+#include "error.h"
+
+/**
+ * @brief Run a brute-force gravity calculation for a subset of particles.
+ *
+ * All gpart with ID modulo SWIFT_GRAVITY_FORCE_CHECKS will get their forces
+ * computed.
+ *
+ * @param s The #space to use.
+ * @param e The #engine (to access the current time).
+ */
+void gravity_exact_force_compute(struct space *s, const struct engine *e) {
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+
+ const ticks tic = getticks();
+ const double const_G = e->physical_constants->const_newton_G;
+ int counter = 0;
+
+ for (size_t i = 0; i < s->nr_gparts; ++i) {
+
+ struct gpart *gpi = &s->gparts[i];
+
+ /* Is the particle active and part of the subset to be tested ? */
+ if (gpi->id_or_neg_offset % SWIFT_GRAVITY_FORCE_CHECKS == 0 &&
+ gpart_is_active(gpi, e)) {
+
+ /* Be ready for the calculation */
+ gpi->a_grav[0] = 0.f;
+ gpi->a_grav[1] = 0.f;
+ gpi->a_grav[2] = 0.f;
+
+ /* Interact it with all other particles in the space.*/
+ for (size_t j = 0; j < s->nr_gparts; ++j) {
+
+ /* No self interaction */
+ if (i == j) continue;
+
+ struct gpart *gpj = &s->gparts[j];
+
+ /* Compute the pairwise distance. */
+ float dx[3] = {gpi->x[0] - gpj->x[0], // x
+ gpi->x[1] - gpj->x[1], // y
+ gpi->x[2] - gpj->x[2]}; // z
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+
+ runner_iact_grav_pp_nonsym(0.f, r2, dx, gpi, gpj);
+ }
+
+ /* Finish the calculation */
+ gravity_end_force(gpi, const_G);
+
+ /* Store the exact answer */
+ gpi->a_grav_exact[0] = gpi->a_grav[0];
+ gpi->a_grav_exact[1] = gpi->a_grav[1];
+ gpi->a_grav_exact[2] = gpi->a_grav[2];
+
+ /* Restore everything */
+ gpi->a_grav[0] = 0.f;
+ gpi->a_grav[1] = 0.f;
+ gpi->a_grav[2] = 0.f;
+
+ counter++;
+ }
+ }
+
+ message("Computed exact gravity for %d gparts.", counter);
+
+ if (e->verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
+#else
+ error("Gravity checking function called without the corresponding flag.");
+#endif
+}
+
+/**
+ * @brief Check the accuracy of the gravity calculation by comparing the
+ * accelerations
+ * to the brute-force computed ones.
+ *
+ * All gpart with ID modulo SWIFT_GRAVITY_FORCE_CHECKS will be checked.
+ *
+ * @param s The #space to use.
+ * @param e The #engine (to access the current time).
+ * @param rel_tol The maximal relative error. Will call error() if one particle
+ * has a larger error.
+ */
+void gravity_exact_force_check(struct space *s, const struct engine *e,
+ float rel_tol) {
+
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+
+ const double const_G = e->physical_constants->const_newton_G;
+
+ int counter = 0;
+
+ /* Some accumulators */
+ float err_rel[3];
+ float err_rel_max[3] = {0.f, 0.f, 0.f};
+ float err_rel_min[3] = {FLT_MAX, FLT_MAX, FLT_MAX};
+ float err_rel_mean[3] = {0.f, 0.f, 0.f};
+ float err_rel_mean2[3] = {0.f, 0.f, 0.f};
+ float err_rel_std[3] = {0.f, 0.f, 0.f};
+
+ for (size_t i = 0; i < s->nr_gparts; ++i) {
+
+ struct gpart *gpi = &s->gparts[i];
+
+ /* Is the particle was active and part of the subset to be tested ? */
+ if (gpi->id_or_neg_offset % SWIFT_GRAVITY_FORCE_CHECKS == 0 &&
+ gpart_is_starting(gpi, e)) {
+
+ /* Compute relative error */
+ for (int k = 0; k < 3; ++k)
+ if (fabsf(gpi->a_grav_exact[k]) > FLT_EPSILON * const_G)
+ err_rel[k] = (gpi->a_grav[k] - gpi->a_grav_exact[k]) /
+ fabsf(gpi->a_grav_exact[k]);
+ else
+ err_rel[k] = 0.f;
+
+ /* Check that we are not above tolerance */
+ if (fabsf(err_rel[0]) > rel_tol || fabsf(err_rel[1]) > rel_tol ||
+ fabsf(err_rel[2]) > rel_tol)
+ error("Error too large ! gp->a_grav=[%e %e %e] gp->a_exact=[%e %e %e]",
+ gpi->a_grav[0], gpi->a_grav[1], gpi->a_grav[2],
+ gpi->a_grav_exact[0], gpi->a_grav_exact[1], gpi->a_grav_exact[2]);
+
+ /* Construct some statistics */
+ for (int k = 0; k < 3; ++k) {
+ err_rel_max[k] = max(err_rel_max[k], fabsf(err_rel[k]));
+ err_rel_min[k] = min(err_rel_min[k], fabsf(err_rel[k]));
+ err_rel_mean[k] += err_rel[k];
+ err_rel_mean2[k] += err_rel[k] * err_rel[k];
+ }
+
+ counter++;
+ }
+ }
+
+ /* Final operation on the stats */
+ if (counter > 0) {
+ for (int k = 0; k < 3; ++k) {
+ err_rel_mean[k] /= counter;
+ err_rel_mean2[k] /= counter;
+ err_rel_std[k] =
+ sqrtf(err_rel_mean2[k] - err_rel_mean[k] * err_rel_mean[k]);
+ }
+ }
+
+ /* Report on the findings */
+ message("Checked gravity for %d gparts.", counter);
+ for (int k = 0; k < 3; ++k)
+ message("Error on a_grav[%d]: min=%e max=%e mean=%e std=%e", k,
+ err_rel_min[k], err_rel_max[k], err_rel_mean[k], err_rel_std[k]);
+
+#else
+ error("Gravity checking function called without the corresponding flag.");
+#endif
+}
diff --git a/src/gravity.h b/src/gravity.h
index 6edcd59d9955029217364d4fc67874216b0e24dc..00b930c00fb2558f274feb2991b78e96dc8b990b 100644
--- a/src/gravity.h
+++ b/src/gravity.h
@@ -22,9 +22,20 @@
/* Config parameters. */
#include "../config.h"
+/* Local headers. */
+#include "const.h"
+#include "engine.h"
+#include "inline.h"
+#include "part.h"
+#include "space.h"
+
/* So far only one model here */
/* Straight-forward import */
#include "./gravity/Default/gravity.h"
#include "./gravity/Default/gravity_iact.h"
+void gravity_exact_force_compute(struct space *s, const struct engine *e);
+void gravity_exact_force_check(struct space *s, const struct engine *e,
+ float rel_tol);
+
#endif
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 3f97070fbb80c8920a929e6df8d93def6ac6041a..93a65e2f5a70e09ad14280cf9c334753359fb8b5 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -21,15 +21,18 @@
#define SWIFT_DEFAULT_GRAVITY_H
#include <float.h>
+#include "gravity_properties.h"
#include "minmax.h"
/**
* @brief Computes the gravity time-step of a given particle due to self-gravity
*
* @param gp Pointer to the g-particle data.
+ * @param grav_props Constants used in the gravity scheme.
*/
__attribute__((always_inline)) INLINE static float
-gravity_compute_timestep_self(const struct gpart* const gp) {
+gravity_compute_timestep_self(const struct gpart* const gp,
+ const struct gravity_props* restrict grav_props) {
const float ac2 = gp->a_grav[0] * gp->a_grav[0] +
gp->a_grav[1] * gp->a_grav[1] +
@@ -37,7 +40,7 @@ gravity_compute_timestep_self(const struct gpart* const gp) {
const float ac = (ac2 > 0.f) ? sqrtf(ac2) : FLT_MIN;
- const float dt = sqrtf(2.f * const_gravity_eta * gp->epsilon / ac);
+ const float dt = sqrtf(2.f * grav_props->eta * gp->epsilon / ac);
return dt;
}
@@ -57,6 +60,10 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(
gp->a_grav[0] = 0.f;
gp->a_grav[1] = 0.f;
gp->a_grav[2] = 0.f;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gp->mass_interacted = 0.;
+#endif
}
/**
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index f484b13663059fa5f4f822aa78748fe4ef9d5926..00ae0f5b05cd95750c34b60e2353a9fc1d0a5c32 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -52,6 +52,9 @@ struct gpart {
#ifdef SWIFT_DEBUG_CHECKS
+ /* Total mass this gpart interacted with */
+ double mass_interacted;
+
/* Time of the last drift */
integertime_t ti_drift;
@@ -60,6 +63,13 @@ struct gpart {
#endif
+#ifdef SWIFT_GRAVITY_FORCE_CHECKS
+
+ /* Brute-force particle acceleration. */
+ float a_grav_exact[3];
+
+#endif
+
} SWIFT_STRUCT_ALIGN;
#endif /* SWIFT_DEFAULT_GRAVITY_PART_H */
diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
new file mode 100644
index 0000000000000000000000000000000000000000..4730d9df5dc573de74fde422e1b7dafc0ee0994a
--- /dev/null
+++ b/src/gravity_derivatives.h
@@ -0,0 +1,72 @@
+/*******************************************************************************
+ * 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_GRAVITY_DERIVATIVE_H
+#define SWIFT_GRAVITY_DERIVATIVE_H
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local headers. */
+#include "inline.h"
+
+/**
+ * @brief \f$ \phi(r_x, r_y, r_z) \f$.
+ */
+__attribute__((always_inline)) INLINE static double D_000(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+
+ return r_inv;
+}
+
+/**
+ * @brief \f$ \frac{\partial\phi(r_x, r_y, r_z)}{\partial r_x} \f$.
+ */
+__attribute__((always_inline)) INLINE static double D_100(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+
+ return -r_x * r_inv * r_inv * r_inv;
+}
+
+/**
+ * @brief \f$ \frac{\partial\phi(r_x, r_y, r_z)}{\partial r_x} \f$.
+ */
+__attribute__((always_inline)) INLINE static double D_010(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+
+ return -r_y * r_inv * r_inv * r_inv;
+}
+
+/**
+ * @brief \f$ \frac{\partial\phi(r_x, r_y, r_z)}{\partial r_x} \f$.
+ */
+__attribute__((always_inline)) INLINE static double D_001(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+
+ return -r_z * r_inv * r_inv * r_inv;
+}
+
+#endif /* SWIFT_GRAVITY_DERIVATIVE_H */
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index bfb5cd1ce39a9908573c66406f41b56561a870d6..a614d08c30b21f9e7d422bf6b6a09d10d2e89799 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -148,6 +148,16 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
return min(dt_cfl, dt_u_change);
}
+/**
+ * @brief Does some extra hydro operations once the actual physical time step
+ * for the particle is known.
+ *
+ * @param p The particle to act upon.
+ * @param dt Physical time step of the particle during the next step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
+ struct part *p, float dt) {}
+
/**
* @brief Prepares a particle for the density calculation.
*
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 160a2d8b5d25a97cefb2afd5e22d8e6bcea0006e..cc7b422ccbe7c678969df5779a4d4a054c65528e 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -152,6 +152,16 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
return dt_cfl;
}
+/**
+ * @brief Does some extra hydro operations once the actual physical time step
+ * for the particle is known.
+ *
+ * @param p The particle to act upon.
+ * @param dt Physical time step of the particle during the next step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
+ struct part *p, float dt) {}
+
/**
* @brief Prepares a particle for the density calculation.
*
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index c59af05460157a756c15d8ca84af8a7834fde2d3..643489912a6c6b1db921e73b508910cc670d49ae 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -1,3 +1,4 @@
+
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
@@ -23,6 +24,7 @@
#include "equation_of_state.h"
#include "hydro_gradients.h"
#include "minmax.h"
+#include "riemann.h"
/**
* @brief Computes the hydro time-step of a given particle
@@ -37,7 +39,46 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
const float CFL_condition = hydro_properties->CFL_condition;
- return CFL_condition * p->h / fabsf(p->timestepvars.vmax);
+ if (p->timestepvars.vmax == 0.) {
+ /* vmax can be zero in vacuum cells that only have vacuum neighbours */
+ /* in this case, the time step should be limited by the maximally
+ allowed time step. Since we do not know what that value is here, we set
+ the time step to a very large value */
+ return FLT_MAX;
+ } else {
+ return CFL_condition * p->h / fabsf(p->timestepvars.vmax);
+ }
+}
+
+/**
+ * @brief Does some extra hydro operations once the actual physical time step
+ * for the particle is known.
+ *
+ * We use this to store the physical time step, since it is used for the flux
+ * exchange during the force loop.
+ *
+ * We also set the active flag of the particle to inactive. It will be set to
+ * active in hydro_init_part, which is called the next time the particle becomes
+ * active.
+ *
+ * @param p The particle to act upon.
+ * @param dt Physical time step of the particle during the next step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
+ struct part* p, float dt) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (dt == 0.) {
+ error("Zero time step assigned to particle!");
+ }
+
+ if (dt != dt) {
+ error("NaN time step assigned to particle!");
+ }
+#endif
+
+ p->force.dt = dt;
+ p->force.active = 0;
}
/**
@@ -58,13 +99,44 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
__attribute__((always_inline)) INLINE static void hydro_first_init_part(
struct part* p, struct xpart* xp) {
- xp->v_full[0] = p->v[0];
- xp->v_full[1] = p->v[1];
- xp->v_full[2] = p->v[2];
+ const float mass = p->conserved.mass;
p->primitives.v[0] = p->v[0];
p->primitives.v[1] = p->v[1];
p->primitives.v[2] = p->v[2];
+
+ /* we can already initialize the momentum */
+ p->conserved.momentum[0] = mass * p->primitives.v[0];
+ p->conserved.momentum[1] = mass * p->primitives.v[1];
+ p->conserved.momentum[2] = mass * p->primitives.v[2];
+
+/* and the thermal energy */
+/* remember that we store the total thermal energy, not the specific thermal
+ energy (as in Gadget) */
+#if defined(EOS_ISOTHERMAL_GAS)
+ /* this overwrites the internal energy from the initial condition file */
+ p->conserved.energy = mass * const_isothermal_internal_energy;
+#else
+ p->conserved.energy *= mass;
+#endif
+
+#ifdef GIZMO_TOTAL_ENERGY
+ /* add the total kinetic energy */
+ p->conserved.energy += 0.5f * (p->conserved.momentum[0] * p->primitives.v[0] +
+ p->conserved.momentum[1] * p->primitives.v[1] +
+ p->conserved.momentum[2] * p->primitives.v[2]);
+#endif
+
+#if defined(GIZMO_FIX_PARTICLES)
+ /* make sure the particles are initially at rest */
+ p->v[0] = 0.;
+ p->v[1] = 0.;
+ p->v[2] = 0.;
+#endif
+
+ xp->v_full[0] = p->v[0];
+ xp->v_full[1] = p->v[1];
+ xp->v_full[2] = p->v[2];
}
/**
@@ -89,6 +161,9 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
p->geometry.matrix_E[2][0] = 0.0f;
p->geometry.matrix_E[2][1] = 0.0f;
p->geometry.matrix_E[2][2] = 0.0f;
+
+ /* Set the active flag to active. */
+ p->force.active = 1;
}
/**
@@ -149,6 +224,17 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
/* compute primitive variables */
/* eqns (3)-(5) */
const float m = p->conserved.mass;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (m == 0.) {
+ error("Mass is 0!");
+ }
+
+ if (volume == 0.) {
+ error("Volume is 0!");
+ }
+#endif
+
float momentum[3];
momentum[0] = p->conserved.momentum[0];
momentum[1] = p->conserved.momentum[1];
@@ -157,8 +243,35 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
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;
+
+#ifdef EOS_ISOTHERMAL_GAS
+ /* although the pressure is not formally used anywhere if an isothermal eos
+ has been selected, we still make sure it is set to the correct value */
+ p->primitives.P = const_isothermal_soundspeed * const_isothermal_soundspeed *
+ p->primitives.rho;
+#else
+
+ float energy = p->conserved.energy;
+
+#ifdef GIZMO_TOTAL_ENERGY
+ /* subtract the kinetic energy; we want the thermal energy */
+ energy -= 0.5f * (momentum[0] * p->primitives.v[0] +
+ momentum[1] * p->primitives.v[1] +
+ momentum[2] * p->primitives.v[2]);
+#endif
+
+ /* energy contains the total thermal energy, we want the specific energy.
+ this is why we divide by the volume, and not by the density */
p->primitives.P = hydro_gamma_minus_one * energy / volume;
+#endif
+
+ /* sanity checks */
+ /* it would probably be safer to throw a warning if netive densities or
+ pressures occur */
+ if (p->primitives.rho < 0.0f || p->primitives.P < 0.0f) {
+ p->primitives.rho = 0.0f;
+ p->primitives.P = 0.0f;
+ }
}
/**
@@ -180,13 +293,11 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
__attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part* restrict p, struct xpart* restrict xp) {
- /* Set the physical time step */
- p->force.dt = get_timestep(p->time_bin, 0.); // MATTHIEU 0
-
/* Initialize time step criterion variables */
p->timestepvars.vmax = 0.0f;
/* Set the actual velocity of the particle */
+ /* if GIZMO_FIX_PARTICLES has been selected, v_full will always be zero */
p->force.v_full[0] = xp->v_full[0];
p->force.v_full[1] = xp->v_full[1];
p->force.v_full[2] = xp->v_full[2];
@@ -246,40 +357,14 @@ __attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
* @brief Converts the hydrodynamic variables from the initial condition file to
* conserved variables that can be used during the integration
*
- * Requires the volume to be known.
- *
- * The initial condition file contains a mixture of primitive and conserved
- * variables. Mass is a conserved variable, and we just copy the particle
- * mass into the corresponding conserved quantity. We need the volume to
- * also derive a density, which is then used to convert the internal energy
- * to a pressure. However, we do not actually use these variables anymore.
- * We do need to initialize the linear momentum, based on the mass and the
- * velocity of the particle.
+ * We no longer do this, as the mass needs to be provided in the initial
+ * condition file, and the mass alone is enough to initialize all conserved
+ * variables. This is now done in hydro_first_init_part.
*
* @param p The particle to act upon.
*/
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
- struct part* p, struct xpart* xp) {
-
- const float volume = p->geometry.volume;
- 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];
-
- p->primitives.P =
- hydro_gamma_minus_one * p->conserved.energy * p->primitives.rho;
-
- p->conserved.energy *= m;
-}
+ struct part* p, struct xpart* xp) {}
/**
* @brief Extra operations to be done during the drift
@@ -303,19 +388,33 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
else
p->h *= expf(w1);
- const float w2 = -hydro_dimension * w1;
- if (fabsf(w2) < 0.2f) {
- p->primitives.rho *= approx_expf(w2);
- } else {
- p->primitives.rho *= expf(w2);
+/* we temporarily disabled the primitive variable drift.
+ This should be reenabled once gravity works, and we have time to check that
+ drifting works properly. */
+// const float w2 = -hydro_dimension * w1;
+// if (fabsf(w2) < 0.2f) {
+// p->primitives.rho *= approx_expf(w2);
+// } else {
+// p->primitives.rho *= expf(w2);
+// }
+
+// p->primitives.v[0] += (p->a_hydro[0] + p->gravity.old_a[0]) * dt;
+// p->primitives.v[1] += (p->a_hydro[1] + p->gravity.old_a[1]) * dt;
+// p->primitives.v[2] += (p->a_hydro[2] + p->gravity.old_a[2]) * dt;
+
+//#if !defined(EOS_ISOTHERMAL_GAS)
+// if (p->conserved.mass > 0.) {
+// const float u = p->conserved.energy + p->du_dt * dt;
+// p->primitives.P =
+// hydro_gamma_minus_one * u * p->primitives.rho / p->conserved.mass;
+// }
+//#endif
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (p->h <= 0.) {
+ error("Zero or negative smoothing length (%g)!", p->h);
}
-
- p->primitives.v[0] += (p->a_hydro[0] + p->gravity.old_a[0]) * dt;
- p->primitives.v[1] += (p->a_hydro[1] + p->gravity.old_a[1]) * dt;
- p->primitives.v[2] += (p->a_hydro[2] + p->gravity.old_a[2]) * dt;
- const float u = p->conserved.energy + p->du_dt * dt;
- p->primitives.P =
- hydro_gamma_minus_one * u * p->primitives.rho / p->conserved.mass;
+#endif
}
/**
@@ -333,35 +432,24 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
__attribute__((always_inline)) INLINE static void hydro_end_force(
struct part* p) {
+ /* set the variables that are used to drift the primitive variables */
+
/* Add normalization to h_dt. */
p->force.h_dt *= p->h * hydro_dimension_inv;
- /* Set the hydro acceleration, based on the new momentum and mass */
- /* NOTE: the momentum and mass are only correct for active particles, since
- only active particles have received flux contributions from all their
- neighbours. Since this method is only called for active particles,
- this is indeed the case. */
if (p->force.dt) {
- float mnew;
- float vnew[3];
-
- mnew = p->conserved.mass + p->conserved.flux.mass;
- vnew[0] = (p->conserved.momentum[0] + p->conserved.flux.momentum[0]) / mnew;
- vnew[1] = (p->conserved.momentum[1] + p->conserved.flux.momentum[1]) / mnew;
- vnew[2] = (p->conserved.momentum[2] + p->conserved.flux.momentum[2]) / mnew;
-
- p->a_hydro[0] = (vnew[0] - p->force.v_full[0]) / p->force.dt;
- p->a_hydro[1] = (vnew[1] - p->force.v_full[1]) / p->force.dt;
- p->a_hydro[2] = (vnew[2] - p->force.v_full[2]) / p->force.dt;
-
p->du_dt = p->conserved.flux.energy / p->force.dt;
} else {
- p->a_hydro[0] = 0.0f;
- p->a_hydro[1] = 0.0f;
- p->a_hydro[2] = 0.0f;
-
p->du_dt = 0.0f;
}
+
+#if defined(GIZMO_FIX_PARTICLES)
+ p->du_dt = 0.0f;
+
+ /* disable the smoothing length update, since the smoothing lengths should
+ stay the same for all steps (particles don't move) */
+ p->force.h_dt = 0.0f;
+#endif
}
/**
@@ -377,63 +465,54 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
__attribute__((always_inline)) INLINE static void hydro_kick_extra(
struct part* p, struct xpart* xp, float dt) {
- float oldm, oldp[3], anew[3];
- const float half_dt = 0.5f * dt; // MATTHIEU
-
- /* Retrieve the current value of the gravitational acceleration from the
- gpart. We are only allowed to do this because this is the kick. We still
- need to check whether gpart exists though.*/
- if (p->gpart) {
- anew[0] = p->gpart->a_grav[0];
- anew[1] = p->gpart->a_grav[1];
- anew[2] = p->gpart->a_grav[2];
-
- /* Copy the old mass and momentum before updating the conserved variables */
- oldm = p->conserved.mass;
- oldp[0] = p->conserved.momentum[0];
- oldp[1] = p->conserved.momentum[1];
- oldp[2] = p->conserved.momentum[2];
- }
+ float a_grav[3];
/* Update conserved variables. */
p->conserved.mass += p->conserved.flux.mass;
p->conserved.momentum[0] += p->conserved.flux.momentum[0];
p->conserved.momentum[1] += p->conserved.flux.momentum[1];
p->conserved.momentum[2] += p->conserved.flux.momentum[2];
+#if defined(EOS_ISOTHERMAL_GAS)
+ p->conserved.energy = p->conserved.mass * const_isothermal_internal_energy;
+#else
p->conserved.energy += p->conserved.flux.energy;
+#endif
/* Add gravity. We only do this if we have gravity activated. */
if (p->gpart) {
- p->conserved.momentum[0] +=
- half_dt * (oldm * p->gravity.old_a[0] + p->conserved.mass * anew[0]);
- p->conserved.momentum[1] +=
- half_dt * (oldm * p->gravity.old_a[1] + p->conserved.mass * anew[1]);
- p->conserved.momentum[2] +=
- half_dt * (oldm * p->gravity.old_a[2] + p->conserved.mass * anew[2]);
-
- float paold, panew;
- paold = oldp[0] * p->gravity.old_a[0] + oldp[1] * p->gravity.old_a[1] +
- oldp[2] * p->gravity.old_a[2];
- panew = p->conserved.momentum[0] * anew[0] +
- p->conserved.momentum[1] * anew[1] +
- p->conserved.momentum[2] * anew[2];
- p->conserved.energy += half_dt * (paold + panew);
-
- float fluxaold, fluxanew;
- fluxaold = p->gravity.old_a[0] * p->gravity.old_mflux[0] +
- p->gravity.old_a[1] * p->gravity.old_mflux[1] +
- p->gravity.old_a[2] * p->gravity.old_mflux[2];
- fluxanew = anew[0] * p->gravity.mflux[0] + anew[1] * p->gravity.mflux[1] +
- anew[2] * p->gravity.mflux[2];
- p->conserved.energy += half_dt * (fluxaold + fluxanew);
-
- /* Store gravitational acceleration and mass flux for next step */
- p->gravity.old_a[0] = anew[0];
- p->gravity.old_a[1] = anew[1];
- p->gravity.old_a[2] = anew[2];
- p->gravity.old_mflux[0] = p->gravity.mflux[0];
- p->gravity.old_mflux[1] = p->gravity.mflux[1];
- p->gravity.old_mflux[2] = p->gravity.mflux[2];
+ /* Retrieve the current value of the gravitational acceleration from the
+ gpart. We are only allowed to do this because this is the kick. We still
+ need to check whether gpart exists though.*/
+ a_grav[0] = p->gpart->a_grav[0];
+ a_grav[1] = p->gpart->a_grav[1];
+ a_grav[2] = p->gpart->a_grav[2];
+
+ /* Store the gravitational acceleration for later use. */
+ /* This is currently only used for output purposes. */
+ p->gravity.old_a[0] = a_grav[0];
+ p->gravity.old_a[1] = a_grav[1];
+ p->gravity.old_a[2] = a_grav[2];
+
+ /* Make sure the gpart knows the mass has changed. */
+ p->gpart->mass = p->conserved.mass;
+
+ /* Kick the momentum for half a time step */
+ /* Note that this also affects the particle movement, as the velocity for
+ the particles is set after this. */
+ p->conserved.momentum[0] += dt * p->conserved.mass * a_grav[0];
+ p->conserved.momentum[1] += dt * p->conserved.mass * a_grav[1];
+ p->conserved.momentum[2] += dt * p->conserved.mass * a_grav[2];
+
+#if !defined(EOS_ISOTHERMAL_GAS) && defined(GIZMO_TOTAL_ENERGY)
+ /* This part still needs to be tested! */
+ p->conserved.energy += dt * (p->conserved.momentum[0] * a_grav[0] +
+ p->conserved.momentum[1] * a_grav[1] +
+ p->conserved.momentum[2] * a_grav[2]);
+
+ p->conserved.energy += dt * (a_grav[0] * p->gravity.mflux[0] +
+ a_grav[1] * p->gravity.mflux[1] +
+ a_grav[2] * p->gravity.mflux[2]);
+#endif
}
/* reset fluxes */
@@ -444,6 +523,47 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->conserved.flux.momentum[1] = 0.0f;
p->conserved.flux.momentum[2] = 0.0f;
p->conserved.flux.energy = 0.0f;
+
+#if defined(GIZMO_FIX_PARTICLES)
+ xp->v_full[0] = 0.;
+ xp->v_full[1] = 0.;
+ xp->v_full[2] = 0.;
+
+ p->v[0] = 0.;
+ p->v[1] = 0.;
+ p->v[2] = 0.;
+
+ if (p->gpart) {
+ p->gpart->v_full[0] = 0.;
+ p->gpart->v_full[1] = 0.;
+ p->gpart->v_full[2] = 0.;
+ }
+#else
+ /* Set particle movement */
+ if (p->conserved.mass > 0.) {
+ xp->v_full[0] = p->conserved.momentum[0] / p->conserved.mass;
+ xp->v_full[1] = p->conserved.momentum[1] / p->conserved.mass;
+ xp->v_full[2] = p->conserved.momentum[2] / p->conserved.mass;
+ } else {
+ /* vacuum particles don't move */
+ xp->v_full[0] = 0.;
+ xp->v_full[1] = 0.;
+ xp->v_full[2] = 0.;
+ }
+ p->v[0] = xp->v_full[0];
+ p->v[1] = xp->v_full[1];
+ p->v[2] = xp->v_full[2];
+
+ /* Update gpart! */
+ /* This is essential, as the gpart drift is done independently from the part
+ drift, and we don't want the gpart and the part to have different
+ positions! */
+ if (p->gpart) {
+ p->gpart->v_full[0] = xp->v_full[0];
+ p->gpart->v_full[1] = xp->v_full[1];
+ p->gpart->v_full[2] = xp->v_full[2];
+ }
+#endif
}
/**
@@ -454,7 +574,15 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
const struct part* restrict p) {
- return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
+ if (p->primitives.rho > 0.) {
+#ifdef EOS_ISOTHERMAL_GAS
+ return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
+#else
+ return const_isothermal_internal_energy;
+#endif
+ } else {
+ return 0.;
+ }
}
/**
@@ -465,7 +593,11 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
__attribute__((always_inline)) INLINE static float hydro_get_entropy(
const struct part* restrict p) {
- return p->primitives.P / pow_gamma(p->primitives.rho);
+ if (p->primitives.rho > 0.) {
+ return p->primitives.P / pow_gamma(p->primitives.rho);
+ } else {
+ return 0.;
+ }
}
/**
@@ -476,7 +608,16 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
const struct part* restrict p) {
- return sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
+ if (p->primitives.rho > 0.) {
+#ifdef EOS_ISOTHERMAL_GAS
+ return sqrtf(const_isothermal_internal_energy * hydro_gamma *
+ hydro_gamma_minus_one);
+#else
+ return sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
+#endif
+ } else {
+ return 0.;
+ }
}
/**
@@ -528,6 +669,13 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
/* conserved.energy is NOT the specific energy (u), but the total thermal
energy (u*m) */
p->conserved.energy = u * p->conserved.mass;
+#ifdef GIZMO_TOTAL_ENERGY
+ /* add the kinetic energy */
+ p->conserved.energy += 0.5f * p->conserved.mass *
+ (p->conserved.momentum[0] * p->primitives.v[0] +
+ p->conserved.momentum[1] * p->primitives.v[1] +
+ p->conserved.momentum[2] * p->primitives.v[2]);
+#endif
p->primitives.P = hydro_gamma_minus_one * p->primitives.rho * u;
}
@@ -543,7 +691,14 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
__attribute__((always_inline)) INLINE static void hydro_set_entropy(
struct part* restrict p, float S) {
- p->conserved.energy = gas_internal_energy_from_entropy(p->primitives.rho, S) *
- p->conserved.mass;
- p->primitives.P = gas_pressure_from_entropy(p->primitives.rho, S);
+ p->conserved.energy = S * pow_gamma_minus_one(p->primitives.rho) *
+ hydro_one_over_gamma_minus_one * p->conserved.mass;
+#ifdef GIZMO_TOTAL_ENERGY
+ /* add the kinetic energy */
+ p->conserved.energy += 0.5f * p->conserved.mass *
+ (p->conserved.momentum[0] * p->primitives.v[0] +
+ p->conserved.momentum[1] * p->primitives.v[1] +
+ p->conserved.momentum[2] * p->primitives.v[2]);
+#endif
+ p->primitives.P = S * pow_gamma(p->primitives.rho);
}
diff --git a/src/hydro/Gizmo/hydro_gradients.h b/src/hydro/Gizmo/hydro_gradients.h
index 90448efc7adb8ccecaaa98c7388f89eaa8d16bcd..a5c1e9038d0d3de6896afe773e3193a2304a6b6b 100644
--- a/src/hydro/Gizmo/hydro_gradients.h
+++ b/src/hydro/Gizmo/hydro_gradients.h
@@ -1,3 +1,4 @@
+
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
@@ -21,6 +22,7 @@
#define SWIFT_HYDRO_GRADIENTS_H
#include "hydro_slope_limiters.h"
+#include "riemann.h"
#if defined(GRADIENTS_SPH)
@@ -140,69 +142,144 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
hydro_slope_limit_face(Wi, Wj, dWi, dWj, xij_i, xij_j, r);
/* time */
- dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
- Wi[2] * pi->primitives.gradients.rho[1] +
- Wi[3] * pi->primitives.gradients.rho[2] +
- Wi[0] * (pi->primitives.gradients.v[0][0] +
- pi->primitives.gradients.v[1][1] +
- pi->primitives.gradients.v[2][2]));
- dWi[1] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[0][0] +
- Wi[2] * pi->primitives.gradients.v[0][1] +
- Wi[3] * pi->primitives.gradients.v[0][2] +
- pi->primitives.gradients.P[0] / Wi[0]);
- dWi[2] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[1][0] +
- Wi[2] * pi->primitives.gradients.v[1][1] +
- Wi[3] * pi->primitives.gradients.v[1][2] +
- pi->primitives.gradients.P[1] / Wi[0]);
- dWi[3] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[2][0] +
- Wi[2] * pi->primitives.gradients.v[2][1] +
- Wi[3] * pi->primitives.gradients.v[2][2] +
- pi->primitives.gradients.P[2] / Wi[0]);
- dWi[4] -=
- 0.5 * mindt * (Wi[1] * pi->primitives.gradients.P[0] +
- Wi[2] * pi->primitives.gradients.P[1] +
- Wi[3] * pi->primitives.gradients.P[2] +
- hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
- pi->primitives.gradients.v[1][1] +
- pi->primitives.gradients.v[2][2]));
-
- dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
- Wj[2] * pj->primitives.gradients.rho[1] +
- Wj[3] * pj->primitives.gradients.rho[2] +
- Wj[0] * (pj->primitives.gradients.v[0][0] +
- pj->primitives.gradients.v[1][1] +
- pj->primitives.gradients.v[2][2]));
- dWj[1] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[0][0] +
- Wj[2] * pj->primitives.gradients.v[0][1] +
- Wj[3] * pj->primitives.gradients.v[0][2] +
- pj->primitives.gradients.P[0] / Wj[0]);
- dWj[2] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[1][0] +
- Wj[2] * pj->primitives.gradients.v[1][1] +
- Wj[3] * pj->primitives.gradients.v[1][2] +
- pj->primitives.gradients.P[1] / Wj[0]);
- dWj[3] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[2][0] +
- Wj[2] * pj->primitives.gradients.v[2][1] +
- Wj[3] * pj->primitives.gradients.v[2][2] +
- pj->primitives.gradients.P[2] / Wj[0]);
- dWj[4] -=
- 0.5 * mindt * (Wj[1] * pj->primitives.gradients.P[0] +
- Wj[2] * pj->primitives.gradients.P[1] +
- Wj[3] * pj->primitives.gradients.P[2] +
- hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
- pj->primitives.gradients.v[1][1] +
- pj->primitives.gradients.v[2][2]));
-
- Wi[0] += dWi[0];
+ if (Wi[0] > 0.0f) {
+#ifdef EOS_ISOTHERMAL_GAS
+ dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
+ Wi[2] * pi->primitives.gradients.rho[1] +
+ Wi[3] * pi->primitives.gradients.rho[2] +
+ Wi[0] * (pi->primitives.gradients.v[0][0] +
+ pi->primitives.gradients.v[1][1] +
+ pi->primitives.gradients.v[2][2]));
+ dWi[1] -= 0.5 * mindt *
+ (Wi[1] * pi->primitives.gradients.v[0][0] +
+ Wi[2] * pi->primitives.gradients.v[0][1] +
+ Wi[3] * pi->primitives.gradients.v[0][2] +
+ const_isothermal_soundspeed * const_isothermal_soundspeed *
+ pi->primitives.gradients.rho[0] / Wi[0]);
+ dWi[2] -= 0.5 * mindt *
+ (Wi[1] * pi->primitives.gradients.v[1][0] +
+ Wi[2] * pi->primitives.gradients.v[1][1] +
+ Wi[3] * pi->primitives.gradients.v[1][2] +
+ const_isothermal_soundspeed * const_isothermal_soundspeed *
+ pi->primitives.gradients.rho[1] / Wi[0]);
+ dWi[3] -= 0.5 * mindt *
+ (Wi[1] * pi->primitives.gradients.v[2][0] +
+ Wi[2] * pi->primitives.gradients.v[2][1] +
+ Wi[3] * pi->primitives.gradients.v[2][2] +
+ const_isothermal_soundspeed * const_isothermal_soundspeed *
+ pi->primitives.gradients.rho[2] / Wi[0]);
+/* we don't care about P in this case */
+#else
+ dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
+ Wi[2] * pi->primitives.gradients.rho[1] +
+ Wi[3] * pi->primitives.gradients.rho[2] +
+ Wi[0] * (pi->primitives.gradients.v[0][0] +
+ pi->primitives.gradients.v[1][1] +
+ pi->primitives.gradients.v[2][2]));
+ dWi[1] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[0][0] +
+ Wi[2] * pi->primitives.gradients.v[0][1] +
+ Wi[3] * pi->primitives.gradients.v[0][2] +
+ pi->primitives.gradients.P[0] / Wi[0]);
+ dWi[2] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[1][0] +
+ Wi[2] * pi->primitives.gradients.v[1][1] +
+ Wi[3] * pi->primitives.gradients.v[1][2] +
+ pi->primitives.gradients.P[1] / Wi[0]);
+ dWi[3] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[2][0] +
+ Wi[2] * pi->primitives.gradients.v[2][1] +
+ Wi[3] * pi->primitives.gradients.v[2][2] +
+ pi->primitives.gradients.P[2] / Wi[0]);
+ dWi[4] -= 0.5 * mindt *
+ (Wi[1] * pi->primitives.gradients.P[0] +
+ Wi[2] * pi->primitives.gradients.P[1] +
+ Wi[3] * pi->primitives.gradients.P[2] +
+ hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
+ pi->primitives.gradients.v[1][1] +
+ pi->primitives.gradients.v[2][2]));
+#endif
+ }
+
+ if (Wj[0] > 0.0f) {
+#ifdef EOS_ISOTHERMAL_GAS
+ dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
+ Wj[2] * pj->primitives.gradients.rho[1] +
+ Wj[3] * pj->primitives.gradients.rho[2] +
+ Wj[0] * (pj->primitives.gradients.v[0][0] +
+ pj->primitives.gradients.v[1][1] +
+ pj->primitives.gradients.v[2][2]));
+ dWj[1] -= 0.5 * mindt *
+ (Wj[1] * pj->primitives.gradients.v[0][0] +
+ Wj[2] * pj->primitives.gradients.v[0][1] +
+ Wj[3] * pj->primitives.gradients.v[0][2] +
+ const_isothermal_soundspeed * const_isothermal_soundspeed *
+ pj->primitives.gradients.rho[0] / Wj[0]);
+ dWj[2] -= 0.5 * mindt *
+ (Wj[1] * pj->primitives.gradients.v[1][0] +
+ Wj[2] * pj->primitives.gradients.v[1][1] +
+ Wj[3] * pj->primitives.gradients.v[1][2] +
+ const_isothermal_soundspeed * const_isothermal_soundspeed *
+ pj->primitives.gradients.rho[1] / Wj[0]);
+ dWj[3] -= 0.5 * mindt *
+ (Wj[1] * pj->primitives.gradients.v[2][0] +
+ Wj[2] * pj->primitives.gradients.v[2][1] +
+ Wj[3] * pj->primitives.gradients.v[2][2] +
+ const_isothermal_soundspeed * const_isothermal_soundspeed *
+ pj->primitives.gradients.rho[2] / Wj[0]);
+#else
+ dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
+ Wj[2] * pj->primitives.gradients.rho[1] +
+ Wj[3] * pj->primitives.gradients.rho[2] +
+ Wj[0] * (pj->primitives.gradients.v[0][0] +
+ pj->primitives.gradients.v[1][1] +
+ pj->primitives.gradients.v[2][2]));
+ dWj[1] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[0][0] +
+ Wj[2] * pj->primitives.gradients.v[0][1] +
+ Wj[3] * pj->primitives.gradients.v[0][2] +
+ pj->primitives.gradients.P[0] / Wj[0]);
+ dWj[2] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[1][0] +
+ Wj[2] * pj->primitives.gradients.v[1][1] +
+ Wj[3] * pj->primitives.gradients.v[1][2] +
+ pj->primitives.gradients.P[1] / Wj[0]);
+ dWj[3] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[2][0] +
+ Wj[2] * pj->primitives.gradients.v[2][1] +
+ Wj[3] * pj->primitives.gradients.v[2][2] +
+ pj->primitives.gradients.P[2] / Wj[0]);
+ dWj[4] -= 0.5 * mindt *
+ (Wj[1] * pj->primitives.gradients.P[0] +
+ Wj[2] * pj->primitives.gradients.P[1] +
+ Wj[3] * pj->primitives.gradients.P[2] +
+ hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
+ pj->primitives.gradients.v[1][1] +
+ pj->primitives.gradients.v[2][2]));
+#endif
+ }
+
+ if (-dWi[0] > Wi[0]) {
+ Wi[0] = 0.0f;
+ } else {
+ Wi[0] += dWi[0];
+ }
Wi[1] += dWi[1];
Wi[2] += dWi[2];
Wi[3] += dWi[3];
- Wi[4] += dWi[4];
+ if (-dWi[4] > Wi[4]) {
+ Wi[4] = 0.0f;
+ } else {
+ Wi[4] += dWi[4];
+ }
- Wj[0] += dWj[0];
+ if (-dWj[0] > Wj[0]) {
+ Wj[0] = 0.0f;
+ } else {
+ Wj[0] += dWj[0];
+ }
Wj[1] += dWj[1];
Wj[2] += dWj[2];
Wj[3] += dWj[3];
- Wj[4] += dWj[4];
+ if (-dWj[4] > Wj[4]) {
+ Wj[4] = 0.0f;
+ } else {
+ Wj[4] += dWj[4];
+ }
}
#endif // SWIFT_HYDRO_GRADIENTS_H
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index aba6bd53c1c9557929426c11a0986e5f02888874..bdcea49e4fdc2a9ce17c43b99dab55d61f4095c8 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/src/hydro/Gizmo/hydro_iact.h
@@ -231,9 +231,15 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
dtj = pj->force.dt;
/* calculate the maximal signal velocity */
- if (Wi[0] && Wj[0]) {
+ if (Wi[0] > 0.0f && Wj[0] > 0.0f) {
+#ifdef EOS_ISOTHERMAL_GAS
+ /* we use a value that is slightly higher than necessary, since the correct
+ value does not always work */
+ vmax = 2.5 * const_isothermal_soundspeed;
+#else
vmax =
sqrtf(hydro_gamma * Wi[4] / Wi[0]) + sqrtf(hydro_gamma * Wj[4] / Wj[0]);
+#endif
} else {
vmax = 0.0f;
}
@@ -274,8 +280,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
float wi_dr = hidp1 * wi_dx;
float wj_dr = hjdp1 * wj_dx;
dvdr *= ri;
- pi->force.h_dt -= pj->conserved.mass * dvdr / pj->primitives.rho * wi_dr;
- if (mode == 1) {
+ if (pj->primitives.rho > 0.) {
+ pi->force.h_dt -= pj->conserved.mass * dvdr / pj->primitives.rho * wi_dr;
+ }
+ if (mode == 1 && pi->primitives.rho > 0.) {
pj->force.h_dt -= pi->conserved.mass * dvdr / pi->primitives.rho * wj_dr;
}
@@ -391,6 +399,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
pi->conserved.flux.momentum[2] -= dti * Anorm * totflux[3];
pi->conserved.flux.energy -= dti * Anorm * totflux[4];
+#ifndef GIZMO_TOTAL_ENERGY
float ekin = 0.5f * (pi->primitives.v[0] * pi->primitives.v[0] +
pi->primitives.v[1] * pi->primitives.v[1] +
pi->primitives.v[2] * pi->primitives.v[2]);
@@ -398,6 +407,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
pi->conserved.flux.energy += dti * Anorm * totflux[2] * pi->primitives.v[1];
pi->conserved.flux.energy += dti * Anorm * totflux[3] * pi->primitives.v[2];
pi->conserved.flux.energy -= dti * Anorm * totflux[0] * ekin;
+#endif
/* here is how it works:
Mode will only be 1 if both particles are ACTIVE and they are in the same
@@ -412,11 +422,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
==> we update particle j if (MODE IS 1) OR (j IS INACTIVE)
*/
- // MATTHIEU
- const integertime_t pj_ti_end = 0; // get_integer_time_end(pj->time_bin);
- const integertime_t pi_ti_end = 0; // get_integer_time_end(pi->time_bin);
-
- if (mode == 1 || pj_ti_end > pi_ti_end) {
+ if (mode == 1 || pj->force.active == 0) {
/* Store mass flux */
mflux = dtj * Anorm * totflux[0];
pj->gravity.mflux[0] -= mflux * dx[0];
@@ -429,6 +435,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
pj->conserved.flux.momentum[2] += dtj * Anorm * totflux[3];
pj->conserved.flux.energy += dtj * Anorm * totflux[4];
+#ifndef GIZMO_TOTAL_ENERGY
ekin = 0.5f * (pj->primitives.v[0] * pj->primitives.v[0] +
pj->primitives.v[1] * pj->primitives.v[1] +
pj->primitives.v[2] * pj->primitives.v[2]);
@@ -436,6 +443,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
pj->conserved.flux.energy -= dtj * Anorm * totflux[2] * pj->primitives.v[1];
pj->conserved.flux.energy -= dtj * Anorm * totflux[3] * pj->primitives.v[2];
pj->conserved.flux.energy += dtj * Anorm * totflux[0] * ekin;
+#endif
}
}
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index dcf8cd8fd6cf37c0d7f4ba718746ec7dc3e79b01..236106a1fb04cc2e5b84f997a2389d583ce17cff 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/src/hydro/Gizmo/hydro_io.h
@@ -23,6 +23,13 @@
#include "io_properties.h"
#include "riemann.h"
+/* Set the description of the particle movement. */
+#if defined(GIZMO_FIX_PARTICLES)
+#define GIZMO_PARTICLE_MOVEMENT "Fixed particles."
+#else
+#define GIZMO_PARTICLE_MOVEMENT "Particles move with flow velocity."
+#endif
+
/**
* @brief Specifies which particle fields to read from a dataset
*
@@ -63,7 +70,15 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
* @return Internal energy of the particle
*/
float convert_u(struct engine* e, struct part* p) {
- return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
+ if (p->primitives.rho > 0.) {
+#ifdef EOS_ISOTHERMAL_GAS
+ return const_isothermal_internal_energy;
+#else
+ return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
+#endif
+ } else {
+ return 0.;
+ }
}
/**
@@ -74,7 +89,11 @@ float convert_u(struct engine* e, struct part* p) {
* @return Entropic function of the particle
*/
float convert_A(struct engine* e, struct part* p) {
- return p->primitives.P / pow_gamma(p->primitives.rho);
+ if (p->primitives.rho > 0.) {
+ return p->primitives.P / pow_gamma(p->primitives.rho);
+ } else {
+ return 0.;
+ }
}
/**
@@ -85,6 +104,9 @@ float convert_A(struct engine* e, struct part* p) {
* @return Total energy of the particle
*/
float convert_Etot(struct engine* e, struct part* p) {
+#ifdef GIZMO_TOTAL_ENERGY
+ return p->conserved.energy;
+#else
float momentum2;
momentum2 = p->conserved.momentum[0] * p->conserved.momentum[0] +
@@ -92,6 +114,7 @@ float convert_Etot(struct engine* e, struct part* p) {
p->conserved.momentum[2] * p->conserved.momentum[2];
return p->conserved.energy + 0.5f * momentum2 / p->conserved.mass;
+#endif
}
/**
@@ -104,7 +127,7 @@ float convert_Etot(struct engine* e, struct part* p) {
void hydro_write_particles(struct part* parts, struct io_props* list,
int* num_fields) {
- *num_fields = 13;
+ *num_fields = 14;
/* List what we want to write */
list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -135,6 +158,8 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
list[12] =
io_make_output_field_convert_part("TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY,
parts, conserved.energy, convert_Etot);
+ list[13] = io_make_output_field("GravAcceleration", FLOAT, 3,
+ UNIT_CONV_ACCELERATION, parts, gravity.old_a);
}
/**
@@ -155,6 +180,9 @@ void writeSPHflavour(hid_t h_grpsph) {
/* Riemann solver information */
io_write_attribute_s(h_grpsph, "Riemann solver type",
RIEMANN_SOLVER_IMPLEMENTATION);
+
+ /* Particle movement information */
+ io_write_attribute_s(h_grpsph, "Particle movement", GIZMO_PARTICLE_MOVEMENT);
}
/**
diff --git a/src/hydro/Gizmo/hydro_part.h b/src/hydro/Gizmo/hydro_part.h
index f6592ca107d8d2c6970f34ebd3929e226b53a355..928011d8201f3f355b00d5fe064267d379278e64 100644
--- a/src/hydro/Gizmo/hydro_part.h
+++ b/src/hydro/Gizmo/hydro_part.h
@@ -178,6 +178,9 @@ struct part {
/* Physical time step of the particle. */
float dt;
+ /* Flag keeping track of whether this is an active or inactive particle. */
+ char active;
+
/* Actual velocity of the particle. */
float v_full[3];
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 20856b7e038855e22aa3776a74ba9f495ff6c93f..56078a82569fb0bc30347d5c01831e9eecfd48f4 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -165,6 +165,16 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
return dt_cfl;
}
+/**
+ * @brief Does some extra hydro operations once the actual physical time step
+ * for the particle is known.
+ *
+ * @param p The particle to act upon.
+ * @param dt Physical time step of the particle during the next step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
+ struct part *p, float dt) {}
+
/**
* @brief Prepares a particle for the density calculation.
*
diff --git a/src/hydro/PressureEntropy/hydro.h b/src/hydro/PressureEntropy/hydro.h
index f22bb8a13a8ba4d896a77bd4c4f5e86bed5a5960..20238896f1458d0abebacca4865968a3a671c886 100644
--- a/src/hydro/PressureEntropy/hydro.h
+++ b/src/hydro/PressureEntropy/hydro.h
@@ -152,6 +152,16 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
return dt_cfl;
}
+/**
+ * @brief Does some extra hydro operations once the actual physical time step
+ * for the particle is known.
+ *
+ * @param p The particle to act upon.
+ * @param dt Physical time step of the particle during the next step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
+ struct part *p, float dt) {}
+
/**
* @brief Prepares a particle for the density calculation.
*
diff --git a/src/hydro_properties.c b/src/hydro_properties.c
index 8db23a3de0123400fe691b0d60f076929e1b0b48..46785b4b2d5b958f6db3bd9813d139575217d6fe 100644
--- a/src/hydro_properties.c
+++ b/src/hydro_properties.c
@@ -34,6 +34,7 @@
#define hydro_props_default_max_iterations 30
#define hydro_props_default_volume_change 2.0f
+#define hydro_props_default_h_max FLT_MAX
void hydro_props_init(struct hydro_props *p,
const struct swift_params *params) {
@@ -43,7 +44,11 @@ void hydro_props_init(struct hydro_props *p,
p->target_neighbours = pow_dimension(p->eta_neighbours) * kernel_norm;
p->delta_neighbours = parser_get_param_float(params, "SPH:delta_neighbours");
- /* Ghost stuff */
+ /* Maximal smoothing length */
+ p->h_max = parser_get_opt_param_float(params, "SPH:h_max",
+ hydro_props_default_h_max);
+
+ /* Number of iterations to converge h */
p->max_smoothing_iterations = parser_get_opt_param_int(
params, "SPH:max_ghost_iterations", hydro_props_default_max_iterations);
@@ -81,6 +86,9 @@ void hydro_props_print(const struct hydro_props *p) {
"(max|dlog(h)/dt|=%f).",
pow_dimension(expf(p->log_max_h_change)), p->log_max_h_change);
+ if (p->h_max != hydro_props_default_h_max)
+ message("Maximal smoothing length allowed: %.4f", p->h_max);
+
if (p->max_smoothing_iterations != hydro_props_default_max_iterations)
message("Maximal iterations in ghost task set to %d (default is %d)",
p->max_smoothing_iterations, hydro_props_default_max_iterations);
@@ -96,6 +104,7 @@ void hydro_props_print_snapshot(hid_t h_grpsph, const struct hydro_props *p) {
io_write_attribute_f(h_grpsph, "Kernel target N_ngb", p->target_neighbours);
io_write_attribute_f(h_grpsph, "Kernel delta N_ngb", p->delta_neighbours);
io_write_attribute_f(h_grpsph, "Kernel eta", p->eta_neighbours);
+ io_write_attribute_f(h_grpsph, "Maximal smoothing length", p->h_max);
io_write_attribute_f(h_grpsph, "CFL parameter", p->CFL_condition);
io_write_attribute_f(h_grpsph, "Volume log(max(delta h))",
p->log_max_h_change);
diff --git a/src/hydro_properties.h b/src/hydro_properties.h
index 859ae10f18a71ba920f44b8c59e63abd57cd7c92..716c4c060c21eb95d05f9d50e13d4681a958a6fd 100644
--- a/src/hydro_properties.h
+++ b/src/hydro_properties.h
@@ -40,7 +40,10 @@ struct hydro_props {
float target_neighbours;
float delta_neighbours;
- /* Kernel properties */
+ /* Maximal smoothing length */
+ float h_max;
+
+ /* Number of iterations to converge h */
int max_smoothing_iterations;
/* Time integration properties */
diff --git a/src/multipole.c b/src/multipole.c
index f0d0c7084fd9bec366f2185cb24887da40591b17..bd5c6d6546fa0546108dcd53d7fe4060293c37a7 100644
--- a/src/multipole.c
+++ b/src/multipole.c
@@ -20,205 +20,3 @@
/* Config parameters. */
#include "../config.h"
-
-/* Some standard headers. */
-#include <strings.h>
-
-/* This object's header. */
-#include "multipole.h"
-
-/**
-* @brief Reset the data of a #multipole.
-*
-* @param m The #multipole.
-*/
-void multipole_reset(struct multipole *m) {
-
- /* Just bzero the struct. */
- bzero(m, sizeof(struct multipole));
-}
-
-/**
-* @brief Init a multipole from a set of particles.
-*
-* @param m The #multipole.
-* @param gparts The #gpart.
-* @param gcount The number of particles.
-*/
-void multipole_init(struct multipole *m, const struct gpart *gparts,
- int gcount) {
-
-#if const_gravity_multipole_order > 2
-#error "Multipoles of order >2 not yet implemented."
-#endif
-
- /* Zero everything */
- multipole_reset(m);
-
- /* Temporary variables */
- double mass = 0.0;
- double com[3] = {0.0, 0.0, 0.0};
-
-#if const_gravity_multipole_order >= 2
- double I_xx = 0.0, I_yy = 0.0, I_zz = 0.0;
- double I_xy = 0.0, I_xz = 0.0, I_yz = 0.0;
-#endif
-
- /* Collect the particle data. */
- for (int k = 0; k < gcount; k++) {
- const float w = gparts[k].mass;
-
- mass += w;
- com[0] += gparts[k].x[0] * w;
- com[1] += gparts[k].x[1] * w;
- com[2] += gparts[k].x[2] * w;
-
-#if const_gravity_multipole_order >= 2
- I_xx += gparts[k].x[0] * gparts[k].x[0] * w;
- I_yy += gparts[k].x[1] * gparts[k].x[1] * w;
- I_zz += gparts[k].x[2] * gparts[k].x[2] * w;
- I_xy += gparts[k].x[0] * gparts[k].x[1] * w;
- I_xz += gparts[k].x[0] * gparts[k].x[2] * w;
- I_yz += gparts[k].x[1] * gparts[k].x[2] * w;
-#endif
- }
-
- const double imass = 1.0 / mass;
-
- /* Store the data on the multipole. */
- m->mass = mass;
- m->CoM[0] = com[0] * imass;
- m->CoM[1] = com[1] * imass;
- m->CoM[2] = com[2] * imass;
-
-#if const_gravity_multipole_order >= 2
- m->I_xx = I_xx - imass * com[0] * com[0];
- m->I_yy = I_yy - imass * com[1] * com[1];
- m->I_zz = I_zz - imass * com[2] * com[2];
- m->I_xy = I_xy - imass * com[0] * com[1];
- m->I_xz = I_xz - imass * com[0] * com[2];
- m->I_yz = I_yz - imass * com[1] * com[2];
-#endif
-}
-
-/**
- * @brief Add the second multipole to the first one.
- *
- * @param ma The #multipole which will contain the sum.
- * @param mb The #multipole to add.
- */
-
-void multipole_add(struct multipole *ma, const struct multipole *mb) {
-
-#if const_gravity_multipole_order > 2
-#error "Multipoles of order >2 not yet implemented."
-#endif
-
- /* Correct the position. */
- const double ma_mass = ma->mass;
- const double mb_mass = mb->mass;
- const double M_tot = ma_mass + mb_mass;
- const double M_tot_inv = 1.0 / M_tot;
-
- const double ma_CoM[3] = {ma->CoM[0], ma->CoM[1], ma->CoM[2]};
- const double mb_CoM[3] = {mb->CoM[0], mb->CoM[1], mb->CoM[2]};
-
-#if const_gravity_multipole_order >= 2
- const double ma_I_xx = (double)ma->I_xx + ma_mass * ma_CoM[0] * ma_CoM[0];
- const double ma_I_yy = (double)ma->I_yy + ma_mass * ma_CoM[1] * ma_CoM[1];
- const double ma_I_zz = (double)ma->I_zz + ma_mass * ma_CoM[2] * ma_CoM[2];
- const double ma_I_xy = (double)ma->I_xy + ma_mass * ma_CoM[0] * ma_CoM[1];
- const double ma_I_xz = (double)ma->I_xz + ma_mass * ma_CoM[0] * ma_CoM[2];
- const double ma_I_yz = (double)ma->I_yz + ma_mass * ma_CoM[1] * ma_CoM[2];
-
- const double mb_I_xx = (double)mb->I_xx + mb_mass * mb_CoM[0] * mb_CoM[0];
- const double mb_I_yy = (double)mb->I_yy + mb_mass * mb_CoM[1] * mb_CoM[1];
- const double mb_I_zz = (double)mb->I_zz + mb_mass * mb_CoM[2] * mb_CoM[2];
- const double mb_I_xy = (double)mb->I_xy + mb_mass * mb_CoM[0] * mb_CoM[1];
- const double mb_I_xz = (double)mb->I_xz + mb_mass * mb_CoM[0] * mb_CoM[2];
- const double mb_I_yz = (double)mb->I_yz + mb_mass * mb_CoM[1] * mb_CoM[2];
-#endif
-
- /* New mass */
- ma->mass = M_tot;
-
- /* New CoM */
- ma->CoM[0] = (ma_CoM[0] * ma_mass + mb_CoM[0] * mb_mass) * M_tot_inv;
- ma->CoM[1] = (ma_CoM[1] * ma_mass + mb_CoM[1] * mb_mass) * M_tot_inv;
- ma->CoM[2] = (ma_CoM[2] * ma_mass + mb_CoM[2] * mb_mass) * M_tot_inv;
-
-/* New quadrupole */
-#if const_gravity_multipole_order >= 2
- ma->I_xx = (ma_I_xx + mb_I_xx) - M_tot * ma->CoM[0] * ma->CoM[0];
- ma->I_yy = (ma_I_yy + mb_I_yy) - M_tot * ma->CoM[1] * ma->CoM[1];
- ma->I_zz = (ma_I_zz + mb_I_zz) - M_tot * ma->CoM[2] * ma->CoM[2];
- ma->I_xy = (ma_I_xy + mb_I_xy) - M_tot * ma->CoM[0] * ma->CoM[1];
- ma->I_xz = (ma_I_xz + mb_I_xz) - M_tot * ma->CoM[0] * ma->CoM[2];
- ma->I_yz = (ma_I_yz + mb_I_yz) - M_tot * ma->CoM[1] * ma->CoM[2];
-#endif
-}
-
-/**
- * @brief Add a particle to the given multipole.
- *
- * @param m The #multipole.
- * @param p The #gpart.
- */
-
-void multipole_addpart(struct multipole *m, struct gpart *p) {
-
- /* #if const_gravity_multipole_order == 1 */
-
- /* /\* Correct the position. *\/ */
- /* float mm = m->coeffs[0], mp = p->mass; */
- /* float w = 1.0f / (mm + mp); */
- /* for (int k = 0; k < 3; k++) m->x[k] = (m->x[k] * mm + p->x[k] * mp) * w;
- */
-
- /* /\* Add the particle to the moments. *\/ */
- /* m->coeffs[0] = mm + mp; */
-
- /* #else */
- /* #error( "Multipoles of order %i not yet implemented." ,
- * const_gravity_multipole_order )
- */
- /* #endif */
-}
-
-/**
- * @brief Add a group of particles to the given multipole.
- *
- * @param m The #multipole.
- * @param p The #gpart array.
- * @param N Number of parts to add.
- */
-
-void multipole_addparts(struct multipole *m, struct gpart *p, int N) {
-
- /* #if const_gravity_multipole_order == 1 */
-
- /* /\* Get the combined mass and positions. *\/ */
- /* double xp[3] = {0.0, 0.0, 0.0}; */
- /* float mp = 0.0f, w; */
- /* for (int k = 0; k < N; k++) { */
- /* w = p[k].mass; */
- /* mp += w; */
- /* xp[0] += p[k].x[0] * w; */
- /* xp[1] += p[k].x[1] * w; */
- /* xp[2] += p[k].x[2] * w; */
- /* } */
-
- /* /\* Correct the position. *\/ */
- /* float mm = m->coeffs[0]; */
- /* w = 1.0f / (mm + mp); */
- /* for (int k = 0; k < 3; k++) m->x[k] = (m->x[k] * mm + xp[k]) * w; */
-
- /* /\* Add the particle to the moments. *\/ */
- /* m->coeffs[0] = mm + mp; */
-
- /* #else */
- /* #error( "Multipoles of order %i not yet implemented." ,
- * const_gravity_multipole_order )
- */
- /* #endif */
-}
diff --git a/src/multipole.h b/src/multipole.h
index dc1f914dd73969bc63f3beffca7f34d0f889edb6..b5c9335ee8fabf46740cefe310fcfecbea3fd77e 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -22,34 +22,359 @@
/* Some standard headers. */
#include <math.h>
+#include <string.h>
/* Includes. */
+//#include "active.h"
+#include "align.h"
#include "const.h"
+#include "error.h"
+#include "gravity_derivatives.h"
#include "inline.h"
#include "kernel_gravity.h"
+#include "minmax.h"
#include "part.h"
-/* Multipole struct. */
+#define multipole_align 128
+
+struct acc_tensor {
+
+ double F_000;
+};
+
+struct pot_tensor {
+
+ double F_000;
+};
+
struct multipole {
- /* Multipole location. */
+ float mass;
+
+ /*! Bulk velocity */
+ float vel[3];
+};
+
+/**
+ * @brief The multipole expansion of a mass distribution.
+ */
+struct gravity_tensors {
+
+ /*! Linking pointer for "memory management". */
+ struct gravity_tensors *next;
+
+ /*! Centre of mass of the matter dsitribution */
double CoM[3];
- /* Multipole mass */
- float mass;
+ /*! The actual content */
+ struct {
+
+ /*! Multipole mass */
+ struct multipole m_pole;
+
+ /*! Field tensor for acceleration along x */
+ struct acc_tensor a_x;
+
+ /*! Field tensor for acceleration along y */
+ struct acc_tensor a_y;
+
+ /*! Field tensor for acceleration along z */
+ struct acc_tensor a_z;
+
+ /*! Field tensor for the potential */
+ struct pot_tensor pot;
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+ /* Total mass this gpart interacted with */
+ double mass_interacted;
+
+#endif
+ };
+} SWIFT_STRUCT_ALIGN;
+
+/**
+ * @brief Reset the data of a #multipole.
+ *
+ * @param m The #multipole.
+ */
+INLINE static void gravity_reset(struct gravity_tensors *m) {
+
+ /* Just bzero the struct. */
+ bzero(m, sizeof(struct gravity_tensors));
+}
+
+INLINE static void gravity_field_tensor_init(struct gravity_tensors *m) {
+
+ bzero(&m->a_x, sizeof(struct acc_tensor));
+ bzero(&m->a_y, sizeof(struct acc_tensor));
+ bzero(&m->a_z, sizeof(struct acc_tensor));
+ bzero(&m->pot, sizeof(struct pot_tensor));
+
+#ifdef SWIFT_DEBUG_CHECKS
+ m->mass_interacted = 0.;
+#endif
+}
+
+/**
+ * @brief Prints the content of a #multipole to stdout.
+ *
+ * Note: Uses directly printf(), not a call to message().
+ *
+ * @param m The #multipole to print.
+ */
+INLINE static void gravity_multipole_print(const struct multipole *m) {
+
+ // printf("CoM= [%12.5e %12.5e %12.5e\n", m->CoM[0], m->CoM[1], m->CoM[2]);
+ printf("Mass= %12.5e\n", m->mass);
+ printf("Vel= [%12.5e %12.5e %12.5e\n", m->vel[0], m->vel[1], m->vel[2]);
+}
+
+/**
+ * @brief Adds a #multipole to another one (i.e. does ma += mb).
+ *
+ * @param ma The multipole to add to.
+ * @param mb The multipole to add.
+ */
+INLINE static void gravity_multipole_add(struct multipole *ma,
+ const struct multipole *mb) {
+
+ const float mass = ma->mass + mb->mass;
+ const float imass = 1.f / mass;
+
+ /* Add the bulk velocities */
+ ma->vel[0] = (ma->vel[0] * ma->mass + mb->vel[0] * mb->mass) * imass;
+ ma->vel[1] = (ma->vel[1] * ma->mass + mb->vel[1] * mb->mass) * imass;
+ ma->vel[2] = (ma->vel[2] * ma->mass + mb->vel[2] * mb->mass) * imass;
+
+ /* Add the masses */
+ ma->mass = mass;
+}
+
+/**
+ * @brief Verifies whether two #multipole's are equal or not.
+ *
+ * @param ma The first #multipole.
+ * @param mb The second #multipole.
+ * @param tolerance The maximal allowed difference for the fields.
+ * @return 1 if the multipoles are equal 0 otherwise.
+ */
+INLINE static int gravity_multipole_equal(const struct multipole *ma,
+ const struct multipole *mb,
+ double tolerance) {
+
+ /* Check CoM */
+ /* if (fabs(ma->CoM[0] - mb->CoM[0]) / fabs(ma->CoM[0] + mb->CoM[0]) >
+ * tolerance) */
+ /* return 0; */
+ /* if (fabs(ma->CoM[1] - mb->CoM[1]) / fabs(ma->CoM[1] + mb->CoM[1]) >
+ * tolerance) */
+ /* return 0; */
+ /* if (fabs(ma->CoM[2] - mb->CoM[2]) / fabs(ma->CoM[2] + mb->CoM[2]) >
+ * tolerance) */
+ /* return 0; */
+
+ /* Check bulk velocity (if non-zero)*/
+ if (fabsf(ma->vel[0] + mb->vel[0]) > 0.f &&
+ fabsf(ma->vel[0] - mb->vel[0]) / fabsf(ma->vel[0] + mb->vel[0]) >
+ tolerance)
+ return 0;
+ if (fabsf(ma->vel[1] + mb->vel[1]) > 0.f &&
+ fabsf(ma->vel[1] - mb->vel[1]) / fabsf(ma->vel[1] + mb->vel[1]) >
+ tolerance)
+ return 0;
+ if (fabsf(ma->vel[2] + mb->vel[2]) > 0.f &&
+ fabsf(ma->vel[2] - mb->vel[2]) / fabsf(ma->vel[2] + mb->vel[2]) >
+ tolerance)
+ return 0;
+
+ /* Check mass */
+ if (fabsf(ma->mass - mb->mass) / fabsf(ma->mass + mb->mass) > tolerance)
+ return 0;
+
+ /* All is good */
+ return 1;
+}
+
+/**
+ * @brief Drifts a #multipole forward in time.
+ *
+ * @param m The #multipole.
+ * @param dt The drift time-step.
+ */
+INLINE static void gravity_multipole_drift(struct gravity_tensors *m,
+ double dt) {
+
+ /* Move the whole thing according to bulk motion */
+ m->CoM[0] += m->m_pole.vel[0];
+ m->CoM[1] += m->m_pole.vel[1];
+ m->CoM[2] += m->m_pole.vel[2];
+}
+
+/**
+ * @brief Applies the forces due to particles j onto particles i directly.
+ *
+ * @param gparts_i The #gpart to update.
+ * @param gcount_i The number of particles to update.
+ * @param gparts_j The #gpart that source the gravity field.
+ * @param gcount_j The number of sources.
+ */
+INLINE static void gravity_P2P(struct gpart *gparts_i, int gcount_i,
+ const struct gpart *gparts_j, int gcount_j) {}
+
+/**
+ * @brief Constructs the #multipole of a bunch of particles around their
+ * centre of mass.
+ *
+ * Corresponds to equation (28c).
+ *
+ * @param m The #multipole (content will be overwritten).
+ * @param gparts The #gpart.
+ * @param gcount The number of particles.
+ */
+INLINE static void gravity_P2M(struct gravity_tensors *m,
+ const struct gpart *gparts, int gcount) {
#if const_gravity_multipole_order >= 2
- /* Quadrupole terms */
- float I_xx, I_yy, I_zz;
- float I_xy, I_xz, I_yz;
+#error "Implementation of P2M kernel missing for this order."
#endif
-};
+
+ /* Temporary variables */
+ double mass = 0.0;
+ double com[3] = {0.0, 0.0, 0.0};
+ float vel[3] = {0.f, 0.f, 0.f};
+
+ /* Collect the particle data. */
+ for (int k = 0; k < gcount; k++) {
+ const float m = gparts[k].mass;
+
+ mass += m;
+ com[0] += gparts[k].x[0] * m;
+ com[1] += gparts[k].x[1] * m;
+ com[2] += gparts[k].x[2] * m;
+ vel[0] += gparts[k].v_full[0] * m;
+ vel[1] += gparts[k].v_full[1] * m;
+ vel[2] += gparts[k].v_full[2] * m;
+ }
+
+ const double imass = 1.0 / mass;
+
+ /* Store the data on the multipole. */
+ m->m_pole.mass = mass;
+ m->CoM[0] = com[0] * imass;
+ m->CoM[1] = com[1] * imass;
+ m->CoM[2] = com[2] * imass;
+ m->m_pole.vel[0] = vel[0] * imass;
+ m->m_pole.vel[1] = vel[1] * imass;
+ m->m_pole.vel[2] = vel[2] * imass;
+}
+
+/**
+ * @brief Creates a copy of #multipole shifted to a new location.
+ *
+ * Corresponds to equation (28d).
+ *
+ * @param m_a The #multipole copy (content will be overwritten).
+ * @param m_b The #multipole to shift.
+ * @param pos_a The position to which m_b will be shifted.
+ * @param pos_b The current postion of the multipole to shift.
+ * @param periodic Is the calculation periodic ?
+ */
+INLINE static void gravity_M2M(struct multipole *m_a,
+ const struct multipole *m_b,
+ const double pos_a[3], const double pos_b[3],
+ int periodic) {
+
+ m_a->mass = m_b->mass;
+
+ m_a->vel[0] = m_b->vel[0];
+ m_a->vel[1] = m_b->vel[1];
+ m_a->vel[2] = m_b->vel[2];
+}
+/**
+ * @brief Compute the field tensors due to a multipole.
+ *
+ * Corresponds to equation (28b).
+ *
+ * @param l_a The field tensor to compute.
+ * @param m_b The multipole creating the field.
+ * @param pos_a The position of the field tensor.
+ * @param pos_b The position of the multipole.
+ * @param periodic Is the calculation periodic ?
+ */
+INLINE static void gravity_M2L(struct gravity_tensors *l_a,
+ const struct multipole *m_b,
+ const double pos_a[3], const double pos_b[3],
+ int periodic) {
+
+ double dx, dy, dz;
+ if (periodic) {
+ dx = box_wrap(pos_a[0] - pos_b[0], 0., 1.);
+ dy = box_wrap(pos_a[1] - pos_b[1], 0., 1.);
+ dz = box_wrap(pos_a[2] - pos_b[2], 0., 1.);
+ } else {
+ dx = pos_a[0] - pos_b[0];
+ dy = pos_a[1] - pos_b[1];
+ dz = pos_a[2] - pos_b[2];
+ }
+ const double r2 = dx * dx + dy * dy + dz * dz;
+
+ const double r_inv = 1. / sqrt(r2);
+
+ /* 1st order multipole term */
+ l_a->a_x.F_000 = D_100(dx, dy, dz, r_inv) * m_b->mass;
+ l_a->a_y.F_000 = D_010(dx, dy, dz, r_inv) * m_b->mass;
+ l_a->a_z.F_000 = D_001(dx, dy, dz, r_inv) * m_b->mass;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ l_a->mass_interacted += m_b->mass;
+#endif
+}
+
+/**
+ * @brief Creates a copy of #acc_tensor shifted to a new location.
+ *
+ * Corresponds to equation (28e).
+ *
+ * @param l_a The #acc_tensor copy (content will be overwritten).
+ * @param l_b The #acc_tensor to shift.
+ * @param pos_a The position to which m_b will be shifted.
+ * @param pos_b The current postion of the multipole to shift.
+ * @param periodic Is the calculation periodic ?
+ */
+INLINE static void gravity_L2L(struct gravity_tensors *l_a,
+ const struct gravity_tensors *l_b,
+ const double pos_a[3], const double pos_b[3],
+ int periodic) {}
+
+/**
+ * @brief Applies the #acc_tensor to a set of #gpart.
+ *
+ * Corresponds to equation (28a).
+ */
+INLINE static void gravity_L2P(const struct gravity_tensors *l,
+ struct gpart *gparts, int gcount) {
+
+ for (int i = 0; i < gcount; ++i) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ struct gpart *gp = &gparts[i];
+
+ // if(gpart_is_active(gp, e)){
+
+ gp->mass_interacted += l->mass_interacted;
+#endif
+ //}
+ }
+}
+
+#if 0
/* Multipole function prototypes. */
-void multipole_add(struct multipole *m_sum, const struct multipole *m_term);
-void multipole_init(struct multipole *m, const struct gpart *gparts,
+void multipole_add(struct gravity_tensors *m_sum, const struct gravity_tensors *m_term);
+void multipole_init(struct gravity_tensors *m, const struct gpart *gparts,
int gcount);
-void multipole_reset(struct multipole *m);
+void multipole_reset(struct gravity_tensors *m);
/* static void multipole_iact_mm(struct multipole *ma, struct multipole *mb, */
/* double *shift); */
@@ -64,7 +389,7 @@ void multipole_reset(struct multipole *m);
* @param shift The periodicity correction.
*/
__attribute__((always_inline)) INLINE static void multipole_iact_mm(
- struct multipole *ma, struct multipole *mb, double *shift) {
+ struct gravity_tensors *ma, struct gravity_tensors *mb, double *shift) {
/* float dx[3], ir, r, r2 = 0.0f, acc; */
/* int k; */
@@ -106,7 +431,7 @@ __attribute__((always_inline)) INLINE static void multipole_iact_mm(
* @param shift The periodicity correction.
*/
__attribute__((always_inline)) INLINE static void multipole_iact_mp(
- struct multipole *m, struct gpart *p, double *shift) {
+ struct gravity_tensors *m, struct gpart *p, double *shift) {
/* float dx[3], ir, r, r2 = 0.0f, acc; */
/* int k; */
@@ -137,4 +462,6 @@ __attribute__((always_inline)) INLINE static void multipole_iact_mp(
/* #endif */
}
+#endif
+
#endif /* SWIFT_MULTIPOLE_H */
diff --git a/src/part.c b/src/part.c
index ecc5ca977ae0716b395cdd61af97382e6603186e..712a13d42a629355290bcf89f31b94fb9670a215 100644
--- a/src/part.c
+++ b/src/part.c
@@ -27,6 +27,7 @@
/* This object's header. */
#include "error.h"
+#include "multipole.h"
#include "part.h"
/**
@@ -233,6 +234,7 @@ MPI_Datatype part_mpi_type;
MPI_Datatype xpart_mpi_type;
MPI_Datatype gpart_mpi_type;
MPI_Datatype spart_mpi_type;
+MPI_Datatype multipole_mpi_type;
/**
* @brief Registers MPI particle types.
@@ -265,5 +267,11 @@ void part_create_mpi_types() {
MPI_Type_commit(&spart_mpi_type) != MPI_SUCCESS) {
error("Failed to create MPI type for sparts.");
}
+ if (MPI_Type_contiguous(
+ sizeof(struct gravity_tensors) / sizeof(unsigned char), MPI_BYTE,
+ &multipole_mpi_type) != MPI_SUCCESS ||
+ MPI_Type_commit(&multipole_mpi_type) != MPI_SUCCESS) {
+ error("Failed to create MPI type for multipole.");
+ }
}
#endif
diff --git a/src/part.h b/src/part.h
index 4ed4b490964b59239faf170218cd099d225f5edd..e9a151f5b6dc6c2670ced942d195e108555c5a4b 100644
--- a/src/part.h
+++ b/src/part.h
@@ -86,6 +86,7 @@ extern MPI_Datatype part_mpi_type;
extern MPI_Datatype xpart_mpi_type;
extern MPI_Datatype gpart_mpi_type;
extern MPI_Datatype spart_mpi_type;
+extern MPI_Datatype multipole_mpi_type;
void part_create_mpi_types();
#endif
diff --git a/src/potential.h b/src/potential.h
index 4c20f4c7b4eaa7ca134b3d90ae790a4710a1f217..e6ab9a5bd6bd91801eae0b3f1e3d8f65778f5065 100644
--- a/src/potential.h
+++ b/src/potential.h
@@ -36,6 +36,8 @@
#include "./potential/isothermal/potential.h"
#elif defined(EXTERNAL_POTENTIAL_DISC_PATCH)
#include "./potential/disc_patch/potential.h"
+#elif defined(EXTERNAL_POTENTIAL_SINE_WAVE)
+#include "./potential/sine_wave/potential.h"
#else
#error "Invalid choice of external potential"
#endif
diff --git a/src/potential/point_mass/potential.h b/src/potential/point_mass/potential.h
index 160b74cae82f8dba584e15eac743e6ddc2002c0a..adea9d912056fd07e134fb98a7603030e897ec7a 100644
--- a/src/potential/point_mass/potential.h
+++ b/src/potential/point_mass/potential.h
@@ -24,6 +24,7 @@
#include "../config.h"
/* Some standard headers. */
+#include <float.h>
#include <math.h>
/* Local includes. */
@@ -84,7 +85,10 @@ __attribute__((always_inline)) INLINE static float external_gravity_timestep(
const float a_2 = g->a_grav[0] * g->a_grav[0] + g->a_grav[1] * g->a_grav[1] +
g->a_grav[2] * g->a_grav[2];
- return potential->timestep_mult * sqrtf(a_2 / dota_2);
+ if (fabsf(dota_2) > 0.f)
+ return potential->timestep_mult * sqrtf(a_2 / dota_2);
+ else
+ return FLT_MAX;
}
/**
diff --git a/src/potential/sine_wave/potential.h b/src/potential/sine_wave/potential.h
new file mode 100644
index 0000000000000000000000000000000000000000..e2e2b8ffcc170c28a5facc8373a81746811a9991
--- /dev/null
+++ b/src/potential/sine_wave/potential.h
@@ -0,0 +1,133 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2017 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
+ *
+ * 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_SINE_WAVE_H
+#define SWIFT_SINE_WAVE_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <float.h>
+#include <math.h>
+
+/* Local includes. */
+#include "const.h"
+#include "error.h"
+#include "parser.h"
+#include "part.h"
+#include "physical_constants.h"
+#include "space.h"
+#include "units.h"
+
+/**
+ * @brief External Potential Properties - Sine wave case
+ */
+struct external_potential {
+
+ /*! Amplitude of the sine wave. */
+ double amplitude;
+
+ /*! Time-step limiting factor. */
+ double timestep_limit;
+};
+
+/**
+ * @brief Computes the time-step from the acceleration due to a sine wave.
+ *
+ * @param time The current time.
+ * @param potential The properties of the potential.
+ * @param phys_const The physical constants in internal units.
+ * @param g Pointer to the g-particle data.
+ */
+__attribute__((always_inline)) INLINE static float external_gravity_timestep(
+ double time, const struct external_potential* restrict potential,
+ const struct phys_const* restrict phys_const,
+ const struct gpart* restrict g) {
+
+ return potential->timestep_limit;
+}
+
+/**
+ * @brief Computes the gravitational acceleration along x given by the sine
+ * wave.
+ *
+ * @param time The current time in internal units.
+ * @param potential The properties of the potential.
+ * @param phys_const The physical constants in internal units.
+ * @param g Pointer to the g-particle data.
+ */
+__attribute__((always_inline)) INLINE static void external_gravity_acceleration(
+ double time, const struct external_potential* restrict potential,
+ const struct phys_const* restrict phys_const, struct gpart* restrict g) {
+
+ g->a_grav[0] = potential->amplitude * sin(2. * M_PI * g->x[0]) /
+ phys_const->const_newton_G;
+}
+
+/**
+ * @brief Computes the gravitational potential energy of a particle in the
+ * sine wave.
+ *
+ * @param time The current time.
+ * @param potential The #external_potential used in the run.
+ * @param phys_const Physical constants in internal units.
+ * @param gp Pointer to the particle data.
+ */
+__attribute__((always_inline)) INLINE static float
+external_gravity_get_potential_energy(
+ double time, const struct external_potential* potential,
+ const struct phys_const* const phys_const, const struct gpart* gp) {
+
+ /* this potential does not really have a potential energy */
+ return 0.;
+}
+
+/**
+ * @brief Initialises the external potential properties in the internal system
+ * of units.
+ *
+ * @param parameter_file The parsed parameter file
+ * @param phys_const Physical constants in internal units
+ * @param us The current internal system of units
+ * @param potential The external potential properties to initialize
+ */
+static INLINE void potential_init_backend(
+ const struct swift_params* parameter_file,
+ const struct phys_const* phys_const, const struct unit_system* us,
+ const struct space* s, struct external_potential* potential) {
+
+ potential->amplitude =
+ parser_get_param_double(parameter_file, "SineWavePotential:amplitude");
+ potential->timestep_limit = parser_get_param_double(
+ parameter_file, "SineWavePotential:timestep_limit");
+}
+
+/**
+ * @brief Prints the properties of the external potential to stdout.
+ *
+ * @param potential The external potential properties.
+ */
+static INLINE void potential_print_backend(
+ const struct external_potential* potential) {
+
+ message("External potential is a sine wave with amplitude %g",
+ potential->amplitude);
+}
+
+#endif /* SWIFT_SINE_WAVE_H */
diff --git a/src/riemann/riemann_exact_isothermal.h b/src/riemann/riemann_exact_isothermal.h
index 8e8a6a4dcce456b55f85148bed7342ad4e651c1b..26fb649f02ddd62c102e31c9191c20f9003c0133 100644
--- a/src/riemann/riemann_exact_isothermal.h
+++ b/src/riemann/riemann_exact_isothermal.h
@@ -112,7 +112,7 @@ __attribute__((always_inline)) INLINE static float riemann_guess_rho(float* WL,
/* Currently three possibilities and not really an algorithm to decide which
one to choose: */
/* just the average */
- return 0.5f * (WL[0] + WR[0]);
+ // return 0.5f * (WL[0] + WR[0]);
/* two rarefaction approximation */
return sqrtf(WL[0] * WR[0] * expf((vL - vR) / const_isothermal_soundspeed));
@@ -276,12 +276,24 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
vL = WL[1] * n_unit[0] + WL[2] * n_unit[1] + WL[3] * n_unit[2];
vR = WR[1] * n_unit[0] + WR[2] * n_unit[1] + WR[3] * n_unit[2];
- /* VACUUM... */
+ if (WL[0] == 0. || WR[0] == 0.) {
+ error(
+ "One of the states is vacuum, the isothermal solver cannot solve "
+ "this!");
+ }
rho = 0.;
/* obtain a first guess for p */
rhoguess = riemann_guess_rho(WL, WR, vL, vR);
- frho = riemann_f(rho, WL, WR, vL, vR);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (rhoguess <= 0.) {
+ error("Zero or negative initial density guess.");
+ }
+#endif
+
+ /* we know that the value of riemann_f for rho=0 is negative (it is -inf). */
+ frho = -1.;
frhoguess = riemann_f(rhoguess, WL, WR, vL, vR);
/* ok, rhostar is close to 0, better use Brent's method... */
/* we use Newton-Raphson until we find a suitable interval */
@@ -289,14 +301,18 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
/* Newton-Raphson until convergence or until suitable interval is found
to use Brent's method */
unsigned int counter = 0;
- while (fabs(rho - rhoguess) > 1.e-6f * 0.5f * (rho + rhoguess) &&
+ while (fabs(rho - rhoguess) > 5.e-7f * (rho + rhoguess) &&
frhoguess < 0.0f) {
rho = rhoguess;
rhoguess = rhoguess - frhoguess / riemann_fprime(rhoguess, WL, WR);
frhoguess = riemann_f(rhoguess, WL, WR, vL, vR);
counter++;
if (counter > 1000) {
- error("Stuck in Newton-Raphson!\n");
+ error(
+ "Stuck in Newton-Raphson (rho: %g, rhoguess: %g, frhoguess: %g, "
+ "fprime: %g, rho-rhoguess: %g, WL: %g %g %g, WR: %g %g %g)!\n",
+ rho, rhoguess, frhoguess, riemann_fprime(rhoguess, WL, WR),
+ (rho - rhoguess), WL[0], vL, WL[4], WR[0], vR, WR[4]);
}
}
}
@@ -304,7 +320,7 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
if (1.e6 * fabs(rho - rhoguess) > 0.5f * (rho + rhoguess) &&
frhoguess > 0.0f) {
rho = 0.0f;
- frho = riemann_f(rho, WL, WR, vL, vR);
+ frho = -1.;
/* use Brent's method to find the zeropoint */
rho = riemann_solve_brent(rho, rhoguess, frho, frhoguess, 1.e-6, WL, WR, vL,
vR);
diff --git a/src/runner.c b/src/runner.c
index 7cea20de87d8f1602aa1a058cfb8128c8186f680..792f03f9b454d2afe0116f9bda194bb00ce79543 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -495,6 +495,10 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
/* Anything to do here? */
if (!cell_is_active(c, e)) return;
+ /* Reset the gravity acceleration tensors */
+ if (e->policy & engine_policy_self_gravity)
+ gravity_field_tensor_init(c->multipole);
+
/* Recurse? */
if (c->split) {
for (int k = 0; k < 8; k++)
@@ -592,6 +596,7 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
struct part *restrict parts = c->parts;
struct xpart *restrict xparts = c->xparts;
const struct engine *e = r->e;
+ const float hydro_h_max = e->hydro_properties->h_max;
const float target_wcount = e->hydro_properties->target_neighbours;
const float max_wcount =
target_wcount + e->hydro_properties->delta_neighbours;
@@ -663,15 +668,23 @@ void runner_do_ghost(struct runner *r, struct cell *c, int timer) {
/* Ok, correct then */
p->h += h_corr;
- /* Flag for another round of fun */
- pid[redo] = pid[i];
- redo += 1;
+ /* If below the absolute maximum, try again */
+ if (p->h < hydro_h_max) {
+
+ /* Flag for another round of fun */
+ pid[redo] = pid[i];
+ redo += 1;
+
+ /* Re-initialise everything */
+ hydro_init_part(p);
- /* Re-initialise everything */
- hydro_init_part(p);
+ /* Off we go ! */
+ continue;
+ } else {
- /* Off we go ! */
- continue;
+ /* Ok, this particle is a lost cause... */
+ p->h = hydro_h_max;
+ }
}
/* We now have a particle whose smoothing length has converged */
@@ -819,25 +832,6 @@ void runner_do_drift_particles(struct runner *r, struct cell *c, int timer) {
if (timer) TIMER_TOC(timer_drift);
}
-/**
- * @brief Mapper function to drift ALL particle and g-particles forward in time.
- *
- * @param map_data An array of #cell%s.
- * @param num_elements Chunk size.
- * @param extra_data Pointer to an #engine.
- */
-void runner_do_drift_mapper(void *map_data, int num_elements,
- void *extra_data) {
-
- struct engine *e = (struct engine *)extra_data;
- struct cell *cells = (struct cell *)map_data;
-
- for (int ind = 0; ind < num_elements; ind++) {
- struct cell *c = &cells[ind];
- if (c != NULL && c->nodeID == e->nodeID) cell_drift_particles(c, e);
- }
-}
-
/**
* @brief Perform the first half-kick on all the active particles in a cell.
*
@@ -1115,6 +1109,7 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
struct xpart *restrict xparts = c->xparts;
struct gpart *restrict gparts = c->gparts;
struct spart *restrict sparts = c->sparts;
+ const double timeBase = e->timeBase;
TIMER_TIC;
@@ -1150,6 +1145,10 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
p->time_bin = get_time_bin(ti_new_step);
if (p->gpart != NULL) p->gpart->time_bin = get_time_bin(ti_new_step);
+ /* Tell the particle what the new physical time step is */
+ float dt = get_timestep(p->time_bin, timeBase);
+ hydro_timestep_extra(p, dt);
+
/* Number of updated particles */
updated++;
if (p->gpart != NULL) g_updated++;
@@ -1367,16 +1366,27 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
struct gpart *restrict gp = &gparts[k];
if (gp->type == swift_type_dark_matter) {
+
if (gpart_is_active(gp, e)) gravity_end_force(gp, const_G);
}
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (e->policy & engine_policy_self_gravity) {
+ gp->mass_interacted += gp->mass;
+ if (fabs(gp->mass_interacted - e->s->total_mass) > gp->mass)
+ error(
+ "g-particle did not interact gravitationally with all other "
+ "particles gp->mass_interacted=%e, total_mass=%e, gp->mass=%e",
+ gp->mass_interacted, e->s->total_mass, gp->mass);
+ }
+#endif
}
/* Loop over the star particles in this cell. */
for (int k = 0; k < scount; k++) {
- /* Get a handle on the part. */
+ /* Get a handle on the spart. */
struct spart *restrict sp = &sparts[k];
-
if (spart_is_active(sp, e)) {
/* First, finish the force loop */
@@ -1408,6 +1418,8 @@ void runner_do_recv_part(struct runner *r, struct cell *c, int timer) {
integertime_t ti_end_min = max_nr_timesteps;
integertime_t ti_end_max = 0;
+ timebin_t time_bin_min = num_time_bins;
+ timebin_t time_bin_max = 0;
float h_max = 0.f;
/* If this cell is a leaf, collect the particle data. */
@@ -1415,10 +1427,9 @@ void runner_do_recv_part(struct runner *r, struct cell *c, int timer) {
/* Collect everything... */
for (size_t k = 0; k < nr_parts; k++) {
- const integertime_t ti_end =
- get_integer_time_end(ti_current, parts[k].time_bin);
- ti_end_min = min(ti_end_min, ti_end);
- ti_end_max = max(ti_end_max, ti_end);
+ if (parts[k].time_bin == time_bin_inhibited) continue;
+ time_bin_min = min(time_bin_min, parts[k].time_bin);
+ time_bin_max = max(time_bin_max, parts[k].time_bin);
h_max = max(h_max, parts[k].h);
#ifdef SWIFT_DEBUG_CHECKS
@@ -1426,6 +1437,10 @@ void runner_do_recv_part(struct runner *r, struct cell *c, int timer) {
error("Received un-drifted particle !");
#endif
}
+
+ /* Convert into a time */
+ ti_end_min = get_integer_time_end(ti_current, time_bin_min);
+ ti_end_max = get_integer_time_end(ti_current, time_bin_max);
}
/* Otherwise, recurse and collect. */
@@ -1480,22 +1495,27 @@ void runner_do_recv_gpart(struct runner *r, struct cell *c, int timer) {
integertime_t ti_end_min = max_nr_timesteps;
integertime_t ti_end_max = 0;
+ timebin_t time_bin_min = num_time_bins;
+ timebin_t time_bin_max = 0;
/* If this cell is a leaf, collect the particle data. */
if (!c->split) {
/* Collect everything... */
for (size_t k = 0; k < nr_gparts; k++) {
- const integertime_t ti_end =
- get_integer_time_end(ti_current, gparts[k].time_bin);
- ti_end_min = min(ti_end_min, ti_end);
- ti_end_max = max(ti_end_max, ti_end);
+ if (gparts[k].time_bin == time_bin_inhibited) continue;
+ time_bin_min = min(time_bin_min, gparts[k].time_bin);
+ time_bin_max = max(time_bin_max, gparts[k].time_bin);
#ifdef SWIFT_DEBUG_CHECKS
if (gparts[k].ti_drift != ti_current)
error("Received un-drifted g-particle !");
#endif
}
+
+ /* Convert into a time */
+ ti_end_min = get_integer_time_end(ti_current, time_bin_min);
+ ti_end_max = get_integer_time_end(ti_current, time_bin_max);
}
/* Otherwise, recurse and collect. */
@@ -1548,17 +1568,27 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int timer) {
integertime_t ti_end_min = max_nr_timesteps;
integertime_t ti_end_max = 0;
+ timebin_t time_bin_min = num_time_bins;
+ timebin_t time_bin_max = 0;
/* If this cell is a leaf, collect the particle data. */
if (!c->split) {
/* Collect everything... */
for (size_t k = 0; k < nr_sparts; k++) {
- const integertime_t ti_end =
- get_integer_time_end(ti_current, sparts[k].time_bin);
- ti_end_min = min(ti_end_min, ti_end);
- ti_end_max = max(ti_end_max, ti_end);
+ if (sparts[k].time_bin == time_bin_inhibited) continue;
+ time_bin_min = min(time_bin_min, sparts[k].time_bin);
+ time_bin_max = max(time_bin_max, sparts[k].time_bin);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (sparts[k].ti_drift != ti_current)
+ error("Received un-drifted s-particle !");
+#endif
}
+
+ /* Convert into a time */
+ ti_end_min = get_integer_time_end(ti_current, time_bin_min);
+ ti_end_max = get_integer_time_end(ti_current, time_bin_max);
}
/* Otherwise, recurse and collect. */
@@ -1631,6 +1661,8 @@ void *runner_main(void *data) {
/* Get the cells. */
struct cell *ci = t->ci;
struct cell *cj = t->cj;
+
+/* Mark the thread we run on */
#ifdef SWIFT_DEBUG_TASKS
t->rid = r->cpuid;
#endif
@@ -1641,8 +1673,7 @@ void *runner_main(void *data) {
#ifndef WITH_MPI
if (ci == NULL && cj == NULL) {
- if (t->type != task_type_grav_gather_m && t->type != task_type_grav_fft)
- error("Task not associated with cells!");
+ error("Task not associated with cells!");
} else if (cj == NULL) { /* self */
@@ -1802,21 +1833,24 @@ void *runner_main(void *data) {
runner_do_recv_gpart(r, ci, 1);
} else if (t->subtype == task_subtype_spart) {
runner_do_recv_spart(r, ci, 1);
+ } else if (t->subtype == task_subtype_multipole) {
+ ci->ti_old_multipole = e->ti_current;
} else {
error("Unknown/invalid task subtype (%d).", t->subtype);
}
break;
#endif
case task_type_grav_mm:
- runner_do_grav_mm(r, t->ci, 1);
+ // runner_do_grav_mm(r, t->ci, 1);
break;
- case task_type_grav_up:
- runner_do_grav_up(r, t->ci);
+ case task_type_grav_down:
+ runner_do_grav_down(r, t->ci);
break;
- case task_type_grav_gather_m:
+ case task_type_grav_top_level:
+ // runner_do_grav_top_level(r);
break;
- case task_type_grav_fft:
- runner_do_grav_fft(r);
+ case task_type_grav_long_range:
+ runner_do_grav_long_range(r, t->ci, 1);
break;
case task_type_cooling:
if (e->policy & engine_policy_cooling) runner_do_end_force(r, ci, 1);
@@ -1829,6 +1863,18 @@ void *runner_main(void *data) {
error("Unknown/invalid task type (%d).", t->type);
}
+/* Mark that we have run this task on these cells */
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ci != NULL) {
+ ci->tasks_executed[t->type]++;
+ ci->subtasks_executed[t->subtype]++;
+ }
+ if (cj != NULL) {
+ cj->tasks_executed[t->type]++;
+ cj->subtasks_executed[t->subtype]++;
+ }
+#endif
+
/* We're done with this task, see if we get a next one. */
prev = t;
t = scheduler_done(sched, t);
diff --git a/src/runner.h b/src/runner.h
index 49ad5c98354dac157d30ba292f4e48ce05301586..ec63eb3ec98547859f7da75809555f231f277f62 100644
--- a/src/runner.h
+++ b/src/runner.h
@@ -66,6 +66,7 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer);
void *runner_main(void *data);
void runner_do_unskip_mapper(void *map_data, int num_elements,
void *extra_data);
-void runner_do_drift_mapper(void *map_data, int num_elements, void *extra_data);
+void runner_do_drift_all_mapper(void *map_data, int num_elements,
+ void *extra_data);
#endif /* SWIFT_RUNNER_H */
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 9d2606ceb06fd6d32592010376e867a6ae582bf0..9988b7d553a962ee541f20cab64cc534c991957a 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -34,25 +34,78 @@
*/
void runner_do_grav_up(struct runner *r, struct cell *c) {
- if (c->split) { /* Regular node */
+ /* if (c->split) { /\* Regular node *\/ */
+
+ /* /\* Recurse. *\/ */
+ /* for (int k = 0; k < 8; k++) */
+ /* if (c->progeny[k] != NULL) runner_do_grav_up(r, c->progeny[k]); */
+
+ /* /\* Collect the multipoles from the progeny. *\/ */
+ /* multipole_reset(&c->multipole); */
+ /* for (int k = 0; k < 8; k++) { */
+ /* if (c->progeny[k] != NULL) */
+ /* multipole_add(&c->multipole, &c->progeny[k]->multipole); */
+ /* } */
+
+ /* } else { /\* Leaf node. *\/ */
+ /* /\* Just construct the multipole from the gparts. *\/ */
+ /* multipole_init(&c->multipole, c->gparts, c->gcount); */
+ /* } */
+}
+
+void runner_do_grav_down(struct runner *r, struct cell *c) {
+
+ if (c->split) {
- /* Recurse. */
- for (int k = 0; k < 8; k++)
- if (c->progeny[k] != NULL) runner_do_grav_up(r, c->progeny[k]);
+ for (int k = 0; k < 8; ++k) {
+ struct cell *cp = c->progeny[k];
+ struct gravity_tensors temp;
- /* Collect the multipoles from the progeny. */
- multipole_reset(&c->multipole);
- for (int k = 0; k < 8; k++) {
- if (c->progeny[k] != NULL)
- multipole_add(&c->multipole, &c->progeny[k]->multipole);
+ if (cp != NULL) {
+ gravity_L2L(&temp, c->multipole, cp->multipole->CoM, c->multipole->CoM,
+ 1);
+ }
}
- } else { /* Leaf node. */
- /* Just construct the multipole from the gparts. */
- multipole_init(&c->multipole, c->gparts, c->gcount);
+ } else {
+
+ gravity_L2P(c->multipole, c->gparts, c->gcount);
}
}
+/**
+ * @brief Computes the interaction of the field tensor in a cell with the
+ * multipole of another cell.
+ *
+ * @param r The #runner.
+ * @param ci The #cell with field tensor to interact.
+ * @param cj The #cell with the multipole.
+ */
+__attribute__((always_inline)) INLINE static void runner_dopair_grav_mm(
+ const struct runner *r, const struct cell *restrict ci,
+ const struct cell *restrict cj) {
+
+ const struct engine *e = r->e;
+ const int periodic = e->s->periodic;
+ const struct multipole *multi_j = &cj->multipole->m_pole;
+ // const float a_smooth = e->gravity_properties->a_smooth;
+ // const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
+
+ TIMER_TIC;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (multi_j->mass == 0.0) error("Multipole does not seem to have been set.");
+#endif
+
+ /* Anything to do here? */
+ if (!cell_is_active(ci, e)) return;
+
+ gravity_M2L(ci->multipole, multi_j, ci->multipole->CoM, cj->multipole->CoM,
+ periodic);
+
+ TIMER_TOC(timer_dopair_grav_mm);
+}
+
/**
* @brief Computes the interaction of all the particles in a cell with the
* multipole of another cell.
@@ -68,15 +121,16 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
const struct engine *e = r->e;
const int gcount = ci->gcount;
struct gpart *restrict gparts = ci->gparts;
- const struct multipole multi = cj->multipole;
- const float rlr_inv = 1. / (const_gravity_a_smooth * ci->super->width[0]);
+ const struct gravity_tensors *multi = cj->multipole;
+ const float a_smooth = e->gravity_properties->a_smooth;
+ const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
TIMER_TIC;
#ifdef SWIFT_DEBUG_CHECKS
- if (gcount == 0) error("Empty cell!"); // MATTHIEU sanity check
+ if (gcount == 0) error("Empty cell!");
- if (multi.mass == 0.0) // MATTHIEU sanity check
+ if (multi->m_pole.mass == 0.0)
error("Multipole does not seem to have been set.");
#endif
@@ -105,13 +159,17 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
if (!gpart_is_active(gp, e)) continue;
/* Compute the pairwise distance. */
- const float dx[3] = {multi.CoM[0] - gp->x[0], // x
- multi.CoM[1] - gp->x[1], // y
- multi.CoM[2] - gp->x[2]}; // z
+ const float dx[3] = {multi->CoM[0] - gp->x[0], // x
+ multi->CoM[1] - gp->x[1], // y
+ multi->CoM[2] - gp->x[2]}; // z
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Interact !*/
- runner_iact_grav_pm(rlr_inv, r2, dx, gp, &multi);
+ runner_iact_grav_pm(rlr_inv, r2, dx, gp, &multi->m_pole);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gp->mass_interacted += multi->m_pole.mass;
+#endif
}
TIMER_TOC(timer_dopair_grav_pm);
@@ -135,7 +193,8 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
const int gcount_j = cj->gcount;
struct gpart *restrict gparts_i = ci->gparts;
struct gpart *restrict gparts_j = cj->gparts;
- const float rlr_inv = 1. / (const_gravity_a_smooth * ci->super->width[0]);
+ const float a_smooth = e->gravity_properties->a_smooth;
+ const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
TIMER_TIC;
@@ -195,18 +254,30 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
+#ifdef SWIFT_DEBUG_CHECKS
+ gpi->mass_interacted += gpj->mass;
+ gpj->mass_interacted += gpi->mass;
+#endif
+
} else {
if (gpart_is_active(gpi, e)) {
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
+#ifdef SWIFT_DEBUG_CHECKS
+ gpi->mass_interacted += gpj->mass;
+#endif
} else if (gpart_is_active(gpj, e)) {
dx[0] = -dx[0];
dx[1] = -dx[1];
dx[2] = -dx[2];
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gpj->mass_interacted += gpi->mass;
+#endif
}
}
}
@@ -229,7 +300,8 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
const struct engine *e = r->e;
const int gcount = c->gcount;
struct gpart *restrict gparts = c->gparts;
- const float rlr_inv = 1. / (const_gravity_a_smooth * c->super->width[0]);
+ const float a_smooth = e->gravity_properties->a_smooth;
+ const float rlr_inv = 1. / (a_smooth * c->super->width[0]);
TIMER_TIC;
@@ -277,18 +349,31 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
+#ifdef SWIFT_DEBUG_CHECKS
+ gpi->mass_interacted += gpj->mass;
+ gpj->mass_interacted += gpi->mass;
+#endif
+
} else {
if (gpart_is_active(gpi, e)) {
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
+#ifdef SWIFT_DEBUG_CHECKS
+ gpi->mass_interacted += gpj->mass;
+#endif
+
} else if (gpart_is_active(gpj, e)) {
dx[0] = -dx[0];
dx[1] = -dx[1];
dx[2] = -dx[2];
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gpj->mass_interacted += gpi->mass;
+#endif
}
}
}
@@ -466,7 +551,8 @@ static void runner_dosub_grav(struct runner *r, struct cell *ci,
}
}
-static void runner_do_grav_mm(struct runner *r, struct cell *ci, int timer) {
+static void runner_do_grav_long_range(struct runner *r, struct cell *ci,
+ int timer) {
#if ICHECK > 0
for (int pid = 0; pid < ci->gcount; pid++) {
@@ -485,10 +571,10 @@ static void runner_do_grav_mm(struct runner *r, struct cell *ci, int timer) {
const struct engine *e = r->e;
struct cell *cells = e->s->cells_top;
const int nr_cells = e->s->nr_cells;
- const double max_d =
- const_gravity_a_smooth * const_gravity_r_cut * ci->width[0];
- const double max_d2 = max_d * max_d;
- const double pos_i[3] = {ci->loc[0], ci->loc[1], ci->loc[2]};
+ /* const double max_d = */
+ /* const_gravity_a_smooth * const_gravity_r_cut * ci->width[0]; */
+ /* const double max_d2 = max_d * max_d; */
+ // const double pos_i[3] = {ci->loc[0], ci->loc[1], ci->loc[2]};
/* Anything to do here? */
if (!cell_is_active(ci, e)) return;
@@ -501,14 +587,14 @@ static void runner_do_grav_mm(struct runner *r, struct cell *ci, int timer) {
if (ci == cj) continue;
- const double dx[3] = {cj->loc[0] - pos_i[0], // x
- cj->loc[1] - pos_i[1], // y
- cj->loc[2] - pos_i[2]}; // z
- const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+ /* const double dx[3] = {cj->loc[0] - pos_i[0], // x */
+ /* cj->loc[1] - pos_i[1], // y */
+ /* cj->loc[2] - pos_i[2]}; // z */
+ /* const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; */
- if (r2 > max_d2) continue;
+ // if (r2 > max_d2) continue;
- if (!cell_are_neighbours(ci, cj)) runner_dopair_grav_pm(r, ci, cj);
+ if (!cell_are_neighbours(ci, cj)) runner_dopair_grav_mm(r, ci, cj);
}
}
diff --git a/src/scheduler.c b/src/scheduler.c
index f98c1082afbf7ec029a7556e36eb9d18ed37bd0a..9f0f0fd944ebeb1399591214489272d38e150a72 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -1092,8 +1092,7 @@ void scheduler_start(struct scheduler *s) {
if (ci == NULL && cj == NULL) {
- if (t->type != task_type_grav_gather_m && t->type != task_type_grav_fft)
- error("Task not associated with cells!");
+ error("Task not associated with cells!");
} else if (cj == NULL) { /* self */
@@ -1221,6 +1220,9 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
} else if (t->subtype == task_subtype_spart) {
err = MPI_Irecv(t->ci->sparts, t->ci->scount, spart_mpi_type,
t->ci->nodeID, t->flags, MPI_COMM_WORLD, &t->req);
+ } else if (t->subtype == task_subtype_multipole) {
+ err = MPI_Irecv(t->ci->multipole, 1, multipole_mpi_type,
+ t->ci->nodeID, t->flags, MPI_COMM_WORLD, &t->req);
} else {
error("Unknown communication sub-type");
}
@@ -1258,6 +1260,9 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
} else if (t->subtype == task_subtype_spart) {
err = MPI_Isend(t->ci->sparts, t->ci->scount, spart_mpi_type,
t->cj->nodeID, t->flags, MPI_COMM_WORLD, &t->req);
+ } else if (t->subtype == task_subtype_multipole) {
+ err = MPI_Isend(t->ci->multipole, 1, multipole_mpi_type,
+ t->cj->nodeID, t->flags, MPI_COMM_WORLD, &t->req);
} else {
error("Unknown communication sub-type");
}
diff --git a/src/space.c b/src/space.c
index 5deaa308af842436b370c78ad773062550ff242e..00e5b87ad6d3587b02d14d514c9788f8639c5dd5 100644
--- a/src/space.c
+++ b/src/space.c
@@ -51,6 +51,7 @@
#include "lock.h"
#include "memswap.h"
#include "minmax.h"
+#include "multipole.h"
#include "runner.h"
#include "stars.h"
#include "threadpool.h"
@@ -175,18 +176,38 @@ int space_getsid(struct space *s, struct cell **ci, struct cell **cj,
*
* @param s The #space.
* @param c The #cell to recycle.
- * @param rec_begin Pointer to the start of the list of cells to recycle.
- * @param rec_end Pointer to the end of the list of cells to recycle.
+ * @param cell_rec_begin Pointer to the start of the list of cells to recycle.
+ * @param cell_rec_end Pointer to the end of the list of cells to recycle.
+ * @param multipole_rec_begin Pointer to the start of the list of multipoles to
+ * recycle.
+ * @param multipole_rec_end Pointer to the end of the list of multipoles to
+ * recycle.
*/
void space_rebuild_recycle_rec(struct space *s, struct cell *c,
- struct cell **rec_begin, struct cell **rec_end) {
+ struct cell **cell_rec_begin,
+ struct cell **cell_rec_end,
+ struct gravity_tensors **multipole_rec_begin,
+ struct gravity_tensors **multipole_rec_end) {
if (c->split)
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) {
- space_rebuild_recycle_rec(s, c->progeny[k], rec_begin, rec_end);
- c->progeny[k]->next = *rec_begin;
- *rec_begin = c->progeny[k];
- if (*rec_end == NULL) *rec_end = *rec_begin;
+ space_rebuild_recycle_rec(s, c->progeny[k], cell_rec_begin,
+ cell_rec_end, multipole_rec_begin,
+ multipole_rec_end);
+
+ c->progeny[k]->next = *cell_rec_begin;
+ *cell_rec_begin = c->progeny[k];
+
+ if (s->gravity) {
+ c->progeny[k]->multipole->next = *multipole_rec_begin;
+ *multipole_rec_begin = c->progeny[k]->multipole;
+ }
+
+ if (*cell_rec_end == NULL) *cell_rec_end = *cell_rec_begin;
+ if (s->gravity && *multipole_rec_end == NULL)
+ *multipole_rec_end = *multipole_rec_begin;
+
+ c->progeny[k]->multipole = NULL;
c->progeny[k] = NULL;
}
}
@@ -199,9 +220,14 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
for (int k = 0; k < num_elements; k++) {
struct cell *c = &cells[k];
- struct cell *rec_begin = NULL, *rec_end = NULL;
- space_rebuild_recycle_rec(s, c, &rec_begin, &rec_end);
- if (rec_begin != NULL) space_recycle_list(s, rec_begin, rec_end);
+ struct cell *cell_rec_begin = NULL, *cell_rec_end = NULL;
+ struct gravity_tensors *multipole_rec_begin = NULL,
+ *multipole_rec_end = NULL;
+ space_rebuild_recycle_rec(s, c, &cell_rec_begin, &cell_rec_end,
+ &multipole_rec_begin, &multipole_rec_end);
+ if (cell_rec_begin != NULL)
+ space_recycle_list(s, cell_rec_begin, cell_rec_end, multipole_rec_begin,
+ multipole_rec_end);
c->sorts = NULL;
c->nr_tasks = 0;
c->density = NULL;
@@ -222,6 +248,9 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->drift = NULL;
c->cooling = NULL;
c->sourceterms = NULL;
+ c->grav_top_level = NULL;
+ c->grav_long_range = NULL;
+ c->grav_down = NULL;
c->super = c;
if (c->sort != NULL) {
free(c->sort);
@@ -362,6 +391,7 @@ void space_regrid(struct space *s, int verbose) {
threadpool_map(&s->e->threadpool, space_rebuild_recycle_mapper,
s->cells_top, s->nr_cells, sizeof(struct cell), 100, s);
free(s->cells_top);
+ free(s->multipoles_top);
s->maxdepth = 0;
}
@@ -377,19 +407,35 @@ void space_regrid(struct space *s, int verbose) {
s->tot_cells = s->nr_cells = cdim[0] * cdim[1] * cdim[2];
if (posix_memalign((void *)&s->cells_top, cell_align,
s->nr_cells * sizeof(struct cell)) != 0)
- error("Failed to allocate cells.");
+ error("Failed to allocate top-level cells.");
bzero(s->cells_top, s->nr_cells * sizeof(struct cell));
+ /* Allocate the multipoles for the top-level cells. */
+ if (s->gravity) {
+ if (posix_memalign((void *)&s->multipoles_top, multipole_align,
+ s->nr_cells * sizeof(struct gravity_tensors)) != 0)
+ error("Failed to allocate top-level multipoles.");
+ bzero(s->multipoles_top, s->nr_cells * sizeof(struct gravity_tensors));
+ }
+
/* Set the cells' locks */
- for (int k = 0; k < s->nr_cells; k++)
+ for (int k = 0; k < s->nr_cells; k++) {
if (lock_init(&s->cells_top[k].lock) != 0)
- error("Failed to init spinlock.");
+ error("Failed to init spinlock for hydro.");
+ if (lock_init(&s->cells_top[k].glock) != 0)
+ error("Failed to init spinlock for gravity.");
+ if (lock_init(&s->cells_top[k].mlock) != 0)
+ error("Failed to init spinlock for multipoles.");
+ if (lock_init(&s->cells_top[k].slock) != 0)
+ error("Failed to init spinlock for stars.");
+ }
/* Set the cell location and sizes. */
for (int i = 0; i < cdim[0]; i++)
for (int j = 0; j < cdim[1]; j++)
for (int k = 0; k < cdim[2]; k++) {
- struct cell *restrict c = &s->cells_top[cell_getid(cdim, i, j, k)];
+ const size_t cid = cell_getid(cdim, i, j, k);
+ struct cell *restrict c = &s->cells_top[cid];
c->loc[0] = i * s->width[0];
c->loc[1] = j * s->width[1];
c->loc[2] = k * s->width[2];
@@ -403,7 +449,8 @@ void space_regrid(struct space *s, int verbose) {
c->scount = 0;
c->super = c;
c->ti_old = ti_old;
- lock_init(&c->lock);
+ c->ti_old_multipole = ti_old;
+ if (s->gravity) c->multipole = &s->multipoles_top[cid];
}
/* Be verbose about the change. */
@@ -884,6 +931,7 @@ void space_rebuild(struct space *s, int verbose) {
for (int k = 0; k < s->nr_cells; k++) {
struct cell *restrict c = &cells_top[k];
c->ti_old = ti_old;
+ c->ti_old_multipole = ti_old;
c->parts = finger;
c->xparts = xfinger;
c->gparts = gfinger;
@@ -900,6 +948,13 @@ void space_rebuild(struct space *s, int verbose) {
sure that the parts in each cell are ok. */
space_split(s, cells_top, s->nr_cells, verbose);
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that the multipole construction went OK */
+ if (s->gravity)
+ for (int k = 0; k < s->nr_cells; k++)
+ cell_check_multipole(&s->cells_top[k], NULL);
+#endif
+
if (verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
clocks_getunit());
@@ -1985,6 +2040,7 @@ void space_split_recursive(struct space *s, struct cell *c,
cp->gcount = 0;
cp->scount = 0;
cp->ti_old = c->ti_old;
+ cp->ti_old_multipole = c->ti_old_multipole;
cp->loc[0] = c->loc[0];
cp->loc[1] = c->loc[1];
cp->loc[2] = c->loc[2];
@@ -2011,7 +2067,7 @@ void space_split_recursive(struct space *s, struct cell *c,
/* Remove any progeny with zero parts. */
struct cell_buff *progeny_buff = buff, *progeny_gbuff = gbuff,
*progeny_sbuff = sbuff;
- for (int k = 0; k < 8; k++)
+ for (int k = 0; k < 8; k++) {
if (c->progeny[k]->count == 0 && c->progeny[k]->gcount == 0 &&
c->progeny[k]->scount == 0) {
space_recycle(s, c->progeny[k]);
@@ -2029,7 +2085,43 @@ void space_split_recursive(struct space *s, struct cell *c,
if (c->progeny[k]->maxdepth > maxdepth)
maxdepth = c->progeny[k]->maxdepth;
}
+ }
+
+ /* Deal with multipole */
+ if (s->gravity) {
+ /* Reset everything */
+ gravity_reset(c->multipole);
+
+ /* Compute CoM of all progenies */
+ double CoM[3] = {0., 0., 0.};
+ double mass = 0.;
+
+ for (int k = 0; k < 8; ++k) {
+ if (c->progeny[k] != NULL) {
+ const struct gravity_tensors *m = c->progeny[k]->multipole;
+ CoM[0] += m->CoM[0] * m->m_pole.mass;
+ CoM[1] += m->CoM[1] * m->m_pole.mass;
+ CoM[2] += m->CoM[2] * m->m_pole.mass;
+ mass += m->m_pole.mass;
+ }
+ }
+ c->multipole->CoM[0] = CoM[0] / mass;
+ c->multipole->CoM[1] = CoM[1] / mass;
+ c->multipole->CoM[2] = CoM[2] / mass;
+
+ /* Now shift progeny multipoles and add them up */
+ struct multipole temp;
+ for (int k = 0; k < 8; ++k) {
+ if (c->progeny[k] != NULL) {
+ const struct cell *cp = c->progeny[k];
+ const struct multipole *m = &cp->multipole->m_pole;
+ gravity_M2M(&temp, m, c->multipole->CoM, cp->multipole->CoM,
+ s->periodic);
+ gravity_multipole_add(&c->multipole->m_pole, &temp);
+ }
+ }
+ }
}
/* Otherwise, collect the data for this cell. */
@@ -2040,43 +2132,54 @@ void space_split_recursive(struct space *s, struct cell *c,
c->split = 0;
maxdepth = c->depth;
- /* Get dt_min/dt_max. */
+ timebin_t time_bin_min = num_time_bins, time_bin_max = 0;
+
+ /* parts: Get dt_min/dt_max and h_max. */
for (int k = 0; k < count; k++) {
- struct part *p = &parts[k];
- struct xpart *xp = &xparts[k];
- const float h = p->h;
- const integertime_t ti_end =
- get_integer_time_end(e->ti_current, p->time_bin);
- const integertime_t ti_beg =
- get_integer_time_begin(e->ti_current + 1, p->time_bin);
- xp->x_diff[0] = 0.f;
- xp->x_diff[1] = 0.f;
- xp->x_diff[2] = 0.f;
- if (h > h_max) h_max = h;
- if (ti_end < ti_end_min) ti_end_min = ti_end;
- if (ti_end > ti_end_max) ti_end_max = ti_end;
- if (ti_beg > ti_beg_max) ti_beg_max = ti_beg;
+#ifdef SWIFT_DEBUG_CHECKS
+ if (parts[k].time_bin == time_bin_inhibited)
+ error("Inhibited particle present in space_split()");
+#endif
+ time_bin_min = min(time_bin_min, parts[k].time_bin);
+ time_bin_max = max(time_bin_max, parts[k].time_bin);
+ h_max = max(h_max, parts[k].h);
+ }
+ /* parts: Reset x_diff */
+ for (int k = 0; k < count; k++) {
+ xparts[k].x_diff[0] = 0.f;
+ xparts[k].x_diff[1] = 0.f;
+ xparts[k].x_diff[2] = 0.f;
}
+ /* gparts: Get dt_min/dt_max, reset x_diff. */
for (int k = 0; k < gcount; k++) {
- struct gpart *gp = &gparts[k];
- const integertime_t ti_end =
- get_integer_time_end(e->ti_current, gp->time_bin);
- const integertime_t ti_beg =
- get_integer_time_begin(e->ti_current + 1, gp->time_bin);
- gp->x_diff[0] = 0.f;
- gp->x_diff[1] = 0.f;
- gp->x_diff[2] = 0.f;
- if (ti_end < ti_end_min) ti_end_min = ti_end;
- if (ti_end > ti_end_max) ti_end_max = ti_end;
- if (ti_beg > ti_beg_max) ti_beg_max = ti_beg;
+#ifdef SWIFT_DEBUG_CHECKS
+ if (gparts[k].time_bin == time_bin_inhibited)
+ error("Inhibited s-particle present in space_split()");
+#endif
+ time_bin_min = min(time_bin_min, gparts[k].time_bin);
+ time_bin_max = max(time_bin_max, gparts[k].time_bin);
+
+ gparts[k].x_diff[0] = 0.f;
+ gparts[k].x_diff[1] = 0.f;
+ gparts[k].x_diff[2] = 0.f;
}
+ /* sparts: Get dt_min/dt_max */
for (int k = 0; k < scount; k++) {
- struct spart *sp = &sparts[k];
- const integertime_t ti_end =
- get_integer_time_end(e->ti_current, sp->time_bin);
- if (ti_end < ti_end_min) ti_end_min = ti_end;
- if (ti_end > ti_end_max) ti_end_max = ti_end;
+#ifdef SWIFT_DEBUG_CHECKS
+ if (sparts[k].time_bin == time_bin_inhibited)
+ error("Inhibited g-particle present in space_split()");
+#endif
+ time_bin_min = min(time_bin_min, sparts[k].time_bin);
+ time_bin_max = max(time_bin_max, sparts[k].time_bin);
}
+
+ /* Convert into integer times */
+ ti_end_min = get_integer_time_end(e->ti_current, time_bin_min);
+ ti_end_max = get_integer_time_end(e->ti_current, time_bin_max);
+ ti_beg_max = get_integer_time_begin(e->ti_current + 1, time_bin_max);
+
+ /* Construct the multipole and the centre of mass*/
+ if (s->gravity) gravity_P2M(c->multipole, c->gparts, c->gcount);
}
/* Set the values for this cell. */
@@ -2145,8 +2248,9 @@ void space_split_mapper(void *map_data, int num_cells, void *extra_data) {
void space_recycle(struct space *s, struct cell *c) {
/* Clear the cell. */
- if (lock_destroy(&c->lock) != 0 || lock_destroy(&c->glock) != 0)
- error("Failed to destroy spinlock.");
+ if (lock_destroy(&c->lock) != 0 || lock_destroy(&c->glock) != 0 ||
+ lock_destroy(&c->mlock) != 0 || lock_destroy(&c->slock) != 0)
+ error("Failed to destroy spinlocks.");
/* Clear this cell's sort arrays. */
if (c->sort != NULL) free(c->sort);
@@ -2154,6 +2258,12 @@ void space_recycle(struct space *s, struct cell *c) {
/* Lock the space. */
lock_lock(&s->lock);
+ /* Hook the multipole back in the buffer */
+ if (s->gravity) {
+ c->multipole->next = s->multipoles_sub;
+ s->multipoles_sub = c->multipole;
+ }
+
/* Hook this cell into the buffer. */
c->next = s->cells_sub;
s->cells_sub = c;
@@ -2167,22 +2277,32 @@ void space_recycle(struct space *s, struct cell *c) {
* @brief Return a list of used cells to the buffer of unused sub-cells.
*
* @param s The #space.
- * @param list_begin Pointer to the first #cell in the linked list of
+ * @param cell_list_begin Pointer to the first #cell in the linked list of
* cells joined by their @c next pointers.
- * @param list_end Pointer to the last #cell in the linked list of
+ * @param cell_list_end Pointer to the last #cell in the linked list of
* cells joined by their @c next pointers. It is assumed that this
* cell's @c next pointer is @c NULL.
+ * @param multipole_list_begin Pointer to the first #multipole in the linked
+ * list of
+ * multipoles joined by their @c next pointers.
+ * @param multipole_list_end Pointer to the last #multipole in the linked list
+ * of
+ * multipoles joined by their @c next pointers. It is assumed that this
+ * multipole's @c next pointer is @c NULL.
*/
-void space_recycle_list(struct space *s, struct cell *list_begin,
- struct cell *list_end) {
+void space_recycle_list(struct space *s, struct cell *cell_list_begin,
+ struct cell *cell_list_end,
+ struct gravity_tensors *multipole_list_begin,
+ struct gravity_tensors *multipole_list_end) {
int count = 0;
/* Clean up the list of cells. */
- for (struct cell *c = list_begin; c != NULL; c = c->next) {
+ for (struct cell *c = cell_list_begin; c != NULL; c = c->next) {
/* Clear the cell. */
- if (lock_destroy(&c->lock) != 0 || lock_destroy(&c->glock) != 0)
- error("Failed to destroy spinlock.");
+ if (lock_destroy(&c->lock) != 0 || lock_destroy(&c->glock) != 0 ||
+ lock_destroy(&c->mlock) != 0 || lock_destroy(&c->slock) != 0)
+ error("Failed to destroy spinlocks.");
/* Clear this cell's sort arrays. */
if (c->sort != NULL) free(c->sort);
@@ -2194,11 +2314,17 @@ void space_recycle_list(struct space *s, struct cell *list_begin,
/* Lock the space. */
lock_lock(&s->lock);
- /* Hook this cell into the buffer. */
- list_end->next = s->cells_sub;
- s->cells_sub = list_begin;
+ /* Hook the cells into the buffer. */
+ cell_list_end->next = s->cells_sub;
+ s->cells_sub = cell_list_begin;
s->tot_cells -= count;
+ /* Hook the multipoles into the buffer. */
+ if (s->gravity) {
+ multipole_list_end->next = s->multipoles_sub;
+ s->multipoles_sub = multipole_list_begin;
+ }
+
/* Unlock the space. */
lock_unlock_blind(&s->lock);
}
@@ -2222,7 +2348,7 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells) {
/* For each requested cell... */
for (int j = 0; j < nr_cells; j++) {
- /* Is the buffer empty? */
+ /* Is the cell buffer empty? */
if (s->cells_sub == NULL) {
if (posix_memalign((void *)&s->cells_sub, cell_align,
space_cellallocchunk * sizeof(struct cell)) != 0)
@@ -2234,10 +2360,29 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells) {
s->cells_sub[space_cellallocchunk - 1].next = NULL;
}
+ /* Is the multipole buffer empty? */
+ if (s->gravity && s->multipoles_sub == NULL) {
+ if (posix_memalign(
+ (void *)&s->multipoles_sub, multipole_align,
+ space_cellallocchunk * sizeof(struct gravity_tensors)) != 0)
+ error("Failed to allocate more multipoles.");
+
+ /* Constructed a linked list */
+ for (int k = 0; k < space_cellallocchunk - 1; k++)
+ s->multipoles_sub[k].next = &s->multipoles_sub[k + 1];
+ s->multipoles_sub[space_cellallocchunk - 1].next = NULL;
+ }
+
/* Pick off the next cell. */
cells[j] = s->cells_sub;
s->cells_sub = cells[j]->next;
s->tot_cells += 1;
+
+ /* Hook the multipole */
+ if (s->gravity) {
+ cells[j]->multipole = s->multipoles_sub;
+ s->multipoles_sub = cells[j]->multipole->next;
+ }
}
/* Unlock the space. */
@@ -2245,9 +2390,12 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells) {
/* Init some things in the cell we just got. */
for (int j = 0; j < nr_cells; j++) {
+ struct gravity_tensors *temp = cells[j]->multipole;
bzero(cells[j], sizeof(struct cell));
+ cells[j]->multipole = temp;
cells[j]->nodeID = -1;
- if (lock_init(&cells[j]->lock) != 0 || lock_init(&cells[j]->glock) != 0)
+ if (lock_init(&cells[j]->lock) != 0 || lock_init(&cells[j]->glock) != 0 ||
+ lock_init(&cells[j]->mlock) != 0 || lock_init(&cells[j]->slock) != 0)
error("Failed to initialize cell spinlocks.");
}
}
@@ -2553,7 +2701,7 @@ void space_init(struct space *s, const struct swift_params *params,
/* Init the space lock. */
if (lock_init(&s->lock) != 0) error("Failed to create space spin-lock.");
- /* Build the cells and the tasks. */
+ /* Build the cells recursively. */
if (!dry_run) space_regrid(s, verbose);
}
@@ -2694,27 +2842,51 @@ void space_link_cleanup(struct space *s) {
/**
* @brief Checks that all cells have been drifted to a given point in time
*
- * Expensive function. Should only be used for debugging purposes.
+ * Should only be used for debugging purposes.
*
* @param s The #space to check.
* @param ti_drift The (integer) time.
*/
void space_check_drift_point(struct space *s, integertime_t ti_drift) {
-
+#ifdef SWIFT_DEBUG_CHECKS
/* Recursively check all cells */
space_map_cells_pre(s, 1, cell_check_drift_point, &ti_drift);
+#else
+ error("Calling debugging code without debugging flag activated.");
+#endif
}
/**
* @brief Checks that all particles and local cells have a non-zero time-step.
+ *
+ * Should only be used for debugging purposes.
+ *
+ * @param s The #space to check.
*/
void space_check_timesteps(struct space *s) {
#ifdef SWIFT_DEBUG_CHECKS
-
for (int i = 0; i < s->nr_cells; ++i) {
cell_check_timesteps(&s->cells_top[i]);
}
+#else
+ error("Calling debugging code without debugging flag activated.");
+#endif
+}
+/**
+ * @brief Resets all the individual cell task counters to 0.
+ *
+ * Should only be used for debugging purposes.
+ *
+ * @param s The #space to reset.
+ */
+void space_reset_task_counters(struct space *s) {
+#ifdef SWIFT_DEBUG_CHECKS
+ for (int i = 0; i < s->nr_cells; ++i) {
+ cell_reset_task_counters(&s->cells_top[i]);
+ }
+#else
+ error("Calling debugging code without debugging flag activated.");
#endif
}
@@ -2725,6 +2897,7 @@ void space_clean(struct space *s) {
for (int i = 0; i < s->nr_cells; ++i) cell_clean(&s->cells_top[i]);
free(s->cells_top);
+ free(s->multipoles_top);
free(s->parts);
free(s->xparts);
free(s->gparts);
diff --git a/src/space.h b/src/space.h
index 7f53fbca232e72dd2d5a76fb4f23b0f815d5e421..e3886364bfa7b2f85e724993e5fe3c2d30663bd2 100644
--- a/src/space.h
+++ b/src/space.h
@@ -72,6 +72,9 @@ struct space {
/*! Are we doing gravity? */
int gravity;
+ /*! Total mass in the system */
+ double total_mass;
+
/*! Width of the top-level cells. */
double width[3];
@@ -102,6 +105,12 @@ struct space {
/*! Buffer of unused cells for the sub-cells. */
struct cell *cells_sub;
+ /*! The multipoles associated with the top-level (level 0) cells */
+ struct gravity_tensors *multipoles_top;
+
+ /*! Buffer of unused multipoles for the sub-cells. */
+ struct gravity_tensors *multipoles_sub;
+
/*! The total number of parts in the space. */
size_t nr_parts, size_parts;
@@ -186,8 +195,10 @@ void space_sparts_sort_mapper(void *map_data, int num_elements,
void *extra_data);
void space_rebuild(struct space *s, int verbose);
void space_recycle(struct space *s, struct cell *c);
-void space_recycle_list(struct space *s, struct cell *list_begin,
- struct cell *list_end);
+void space_recycle_list(struct space *s, struct cell *cell_list_begin,
+ struct cell *cell_list_end,
+ struct gravity_tensors *multipole_list_begin,
+ struct gravity_tensors *multipole_list_end);
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);
@@ -207,6 +218,7 @@ void space_link_cleanup(struct space *s);
void space_check_drift_point(struct space *s, integertime_t ti_drift);
void space_check_timesteps(struct space *s);
void space_replicate(struct space *s, int replicate, int verbose);
+void space_reset_task_counters(struct space *s);
void space_clean(struct space *s);
#endif /* SWIFT_SPACE_H */
diff --git a/src/task.c b/src/task.c
index b05d782af305b25bf95b25279c6abc2e1f4037c2..dfc3dc538c75297cbe7e79b7d64c1b7e13a015dc 100644
--- a/src/task.c
+++ b/src/task.c
@@ -47,16 +47,32 @@
#include "lock.h"
/* Task type names. */
-const char *taskID_names[task_type_count] = {
- "none", "sort", "self", "pair", "sub_self",
- "sub_pair", "init", "ghost", "extra_ghost", "drift",
- "kick1", "kick2", "timestep", "send", "recv",
- "grav_gather_m", "grav_fft", "grav_mm", "grav_up", "cooling",
- "sourceterms"};
-
+const char *taskID_names[task_type_count] = {"none",
+ "sort",
+ "self",
+ "pair",
+ "sub_self",
+ "sub_pair",
+ "init",
+ "ghost",
+ "extra_ghost",
+ "drift",
+ "kick1",
+ "kick2",
+ "timestep",
+ "send",
+ "recv",
+ "grav_top_level",
+ "grav_long_range",
+ "grav_mm",
+ "grav_down",
+ "cooling",
+ "sourceterms"};
+
+/* Sub-task type names. */
const char *subtaskID_names[task_subtype_count] = {
- "none", "density", "gradient", "force", "grav", "external_grav",
- "tend", "xv", "rho", "gpart", "spart"};
+ "none", "density", "gradient", "force", "grav", "external_grav",
+ "tend", "xv", "rho", "gpart", "multipole", "spart"};
/**
* @brief Computes the overlap between the parts array of two given cells.
@@ -165,20 +181,22 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
error("Task without particles");
break;
- case task_type_grav_gather_m:
- case task_type_grav_fft:
+ case task_type_grav_top_level:
+ case task_type_grav_long_range:
case task_type_grav_mm:
- case task_type_grav_up:
return task_action_multipole;
break;
+ case task_type_grav_down:
+ return task_action_gpart;
+
default:
error("Unknown task_action for task");
return task_action_none;
break;
}
- /* Silence compile warnings */
+ /* Silence compiler warnings */
error("Unknown task_action for task");
return task_action_none;
}
@@ -265,6 +283,10 @@ void task_unlock(struct task *t) {
/* Act based on task type. */
switch (type) {
+ case task_type_init:
+ case task_type_kick1:
+ case task_type_kick2:
+ case task_type_timestep:
case task_type_drift:
cell_unlocktree(ci);
cell_gunlocktree(ci);
@@ -337,6 +359,10 @@ int task_lock(struct task *t) {
#endif
break;
+ case task_type_init:
+ case task_type_kick1:
+ case task_type_kick2:
+ case task_type_timestep:
case task_type_drift:
if (ci->hold || ci->ghold) return 0;
if (cell_locktree(ci) != 0) return 0;
diff --git a/src/task.h b/src/task.h
index f2733318a34421fa39f3130f9e76f1ed09246d55..5aa1f75e622bf174623980199cb0cf115d2e2dea 100644
--- a/src/task.h
+++ b/src/task.h
@@ -51,10 +51,10 @@ enum task_types {
task_type_timestep,
task_type_send,
task_type_recv,
- task_type_grav_gather_m,
- task_type_grav_fft,
+ task_type_grav_top_level,
+ task_type_grav_long_range,
task_type_grav_mm,
- task_type_grav_up,
+ task_type_grav_down,
task_type_cooling,
task_type_sourceterms,
task_type_count
@@ -74,6 +74,7 @@ enum task_subtypes {
task_subtype_xv,
task_subtype_rho,
task_subtype_gpart,
+ task_subtype_multipole,
task_subtype_spart,
task_subtype_count
} __attribute__((packed));
diff --git a/src/timeline.h b/src/timeline.h
index c73b2432b219a8ab0254d21c59102841557a57b9..0f38ff3e9421c122b61caf74290c170b62b2dd92 100644
--- a/src/timeline.h
+++ b/src/timeline.h
@@ -37,6 +37,9 @@ typedef char timebin_t;
/*! The maximal number of timesteps in a simulation */
#define max_nr_timesteps (1LL << (num_time_bins + 1))
+/*! Fictious time-bin to hold inhibited particles */
+#define time_bin_inhibited (num_time_bins + 2)
+
/**
* @brief Returns the integer time interval corresponding to a time bin
*
diff --git a/src/timers.h b/src/timers.h
index 50f630e7fc355808596967d8d7d887583674d24a..4cb4d7e0ba60003ba4caefffe257c929c59a8d9e 100644
--- a/src/timers.h
+++ b/src/timers.h
@@ -46,6 +46,7 @@ enum {
timer_dopair_gradient,
timer_dopair_force,
timer_dopair_grav_pm,
+ timer_dopair_grav_mm,
timer_dopair_grav_pp,
timer_dograv_external,
timer_dosource,
diff --git a/src/timestep.h b/src/timestep.h
index 6497e25984ef36d759afd12616b45c8166816dfb..99ca1b4f90cc7b8894d4dfe0e0d2670da8f01d65 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -78,7 +78,8 @@ __attribute__((always_inline)) INLINE static integertime_t get_gpart_timestep(
e->physical_constants, gp));
if (e->policy & engine_policy_self_gravity)
- new_dt = min(new_dt, gravity_compute_timestep_self(gp));
+ new_dt =
+ min(new_dt, gravity_compute_timestep_self(gp, e->gravity_properties));
/* Limit timestep within the allowed range */
new_dt = min(new_dt, e->dt_max);
@@ -121,7 +122,8 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
e->physical_constants, p->gpart));
if (e->policy & engine_policy_self_gravity)
- new_dt_grav = min(new_dt_grav, gravity_compute_timestep_self(p->gpart));
+ new_dt_grav = min(new_dt_grav, gravity_compute_timestep_self(
+ p->gpart, e->gravity_properties));
}
/* Final time-step is minimum of hydro and gravity */
@@ -163,7 +165,8 @@ __attribute__((always_inline)) INLINE static integertime_t get_spart_timestep(
e->physical_constants, sp->gpart));
if (e->policy & engine_policy_self_gravity)
- new_dt = min(new_dt, gravity_compute_timestep_self(sp->gpart));
+ new_dt = min(new_dt, gravity_compute_timestep_self(sp->gpart,
+ e->gravity_properties));
/* Limit timestep within the allowed range */
new_dt = min(new_dt, e->dt_max);
diff --git a/src/tools.c b/src/tools.c
index ab11d1f5930cf5319aaf6424f1559f144718e154..89ac286fb435c01b361bdea66e62dd2d7f41ee24 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -732,19 +732,24 @@ int compare_particles(struct part a, struct part b, double threshold) {
#endif
}
-/** @brief Computes the forces between all g-particles using the N^2 algorithm
+/**
+ * @brief Computes the forces between all g-particles using the N^2 algorithm
*
* Overwrites the accelerations of the gparts with the values.
* Do not use for actual runs.
*
* @brief gparts The array of particles.
* @brief gcount The number of particles.
+ * @brief constants Physical constants in internal units.
+ * @brief gravity_properties Constants governing the gravity scheme.
*/
void gravity_n2(struct gpart *gparts, const int gcount,
- const struct phys_const *constants, float rlr) {
+ const struct phys_const *constants,
+ const struct gravity_props *gravity_properties, float rlr) {
const float rlr_inv = 1. / rlr;
- const float max_d = const_gravity_r_cut * rlr;
+ const float r_cut = gravity_properties->r_cut;
+ const float max_d = r_cut * rlr;
const float max_d2 = max_d * max_d;
message("rlr_inv= %f", rlr_inv);
diff --git a/src/tools.h b/src/tools.h
index ece3078dce7cc8ab4b15538a1e5d9a990d81b36d..4d9e8d3ef86f9ad2661118acf008797893ea5bd7 100644
--- a/src/tools.h
+++ b/src/tools.h
@@ -23,6 +23,7 @@
#define SWIFT_TOOL_H
#include "cell.h"
+#include "gravity_properties.h"
#include "part.h"
#include "physical_constants.h"
#include "runner.h"
@@ -45,8 +46,8 @@ void pairs_n2(double *dim, struct part *restrict parts, int N, int periodic);
double random_uniform(double a, double b);
void shuffle_particles(struct part *parts, const int count);
void gravity_n2(struct gpart *gparts, const int gcount,
- const struct phys_const *constants, float rlr);
-
+ const struct phys_const *constants,
+ const struct gravity_props *gravity_properties, float rlr);
int compare_values(double a, double b, double threshold, double *absDiff,
double *absSum, double *relDiff);
int compare_particles(struct part a, struct part b, double threshold);
diff --git a/tests/benchmarkInteractions.c b/tests/benchmarkInteractions.c
index 6d6d345bee743d28fb4bdda911bd4bcc4c78205f..e3f558f88dffbab252bf7c06f9e943ff568b6fff 100644
--- a/tests/benchmarkInteractions.c
+++ b/tests/benchmarkInteractions.c
@@ -91,7 +91,10 @@ struct part *make_particles(size_t count, double *offset, double spacing,
p->h = h;
p->id = ++(*partId);
+
+#if !defined(GIZMO_SPH)
p->mass = 1.0f;
+#endif
/* Place rest of particles around the test particle
* with random position within a unit sphere. */
@@ -110,7 +113,9 @@ struct part *make_particles(size_t count, double *offset, double spacing,
p->h = h;
p->id = ++(*partId);
+#if !defined(GIZMO_SPH)
p->mass = 1.0f;
+#endif
}
return particles;
}
@@ -120,6 +125,7 @@ struct part *make_particles(size_t count, double *offset, double spacing,
*/
void prepare_force(struct part *parts, size_t count) {
+#if !defined(GIZMO_SPH)
struct part *p;
for (size_t i = 0; i < count; ++i) {
p = &parts[i];
@@ -130,6 +136,7 @@ void prepare_force(struct part *parts, size_t count) {
p->force.v_sig = 0.0f;
p->force.h_dt = 0.0f;
}
+#endif
}
/**
@@ -144,10 +151,19 @@ void dump_indv_particle_fields(char *fileName, struct part *p) {
"%13e %13e %13e %13e "
"%13e %13e %13e\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->density.rho_dh, p->density.wcount, p->density.wcount_dh,
- p->force.h_dt, p->force.v_sig,
-#if defined(MINIMAL_SPH)
+ p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
+#if defined(GIZMO_SPH)
+ 0., 0.,
+#else
+ p->rho, p->density.rho_dh,
+#endif
+ p->density.wcount, p->density.wcount_dh, p->force.h_dt,
+#if defined(GIZMO_SPH)
+ 0.,
+#else
+ p->force.v_sig,
+#endif
+#if defined(MINIMAL_SPH) || defined(GIZMO_SPH)
0., 0., 0., 0.
#else
p->density.div_v, p->density.rot_v[0], p->density.rot_v[1],
diff --git a/tests/test125cells.c b/tests/test125cells.c
index aedf0e0b77898cadf97383dfc956a3d6b5ebf39f..b3895ffa9fc0e4c147bfbf58dc1b5a6301b02763 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -532,6 +532,7 @@ int main(int argc, char *argv[]) {
engine.s = &space;
engine.time = 0.1f;
engine.ti_current = 8;
+ engine.max_active_bin = num_time_bins;
struct runner runner;
runner.e = &engine;
diff --git a/tests/test27cells.c b/tests/test27cells.c
index 929a148d1f5730b63de79e9a1ab7e25f1ca7311e..599c8713a4b2077444b10f85ed37ee76f5219c5a 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -396,6 +396,7 @@ int main(int argc, char *argv[]) {
engine.s = &space;
engine.time = 0.1f;
engine.ti_current = 8;
+ engine.max_active_bin = num_time_bins;
struct runner runner;
runner.e = &engine;
diff --git a/tests/testKernelGrav.c b/tests/testKernelGrav.c
index 2733a4a4d041149499d354ef217ae85b1cd35b7f..41e085945693efaf658c219d0e5f992ddf023d74 100644
--- a/tests/testKernelGrav.c
+++ b/tests/testKernelGrav.c
@@ -87,7 +87,7 @@ int main() {
printf("\nAll values are consistent\n");
/* Now test the long range function */
- const float a_smooth = const_gravity_a_smooth;
+ const float a_smooth = 4.5f;
for (int k = 1; k < numPoints; ++k) {
diff --git a/tests/testPair.c b/tests/testPair.c
index 8b23cc419a661f4d50ea53948302729784a129f9..c734424bca58b7a8ce05bba2c3392ca5c600fa9e 100644
--- a/tests/testPair.c
+++ b/tests/testPair.c
@@ -252,6 +252,7 @@ int main(int argc, char *argv[]) {
engine.s = &space;
engine.time = 0.1f;
engine.ti_current = 8;
+ engine.max_active_bin = num_time_bins;
runner.e = &engine;
volume = particles * particles * particles;