diff --git a/.gitignore b/.gitignore
index ddeaaf8b235dc217fbfb6559e66bd665d1f31745..4c704f0665ed0b73d3abb0ca576bfd4709e080e2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -41,8 +41,13 @@ tests/input.hdf5
tests/testSingle
tests/testTimeIntegration
tests/testSPHStep
+tests/testKernel
+tests/testParser
theory/latex/swift.pdf
+theory/kernel/kernels.pdf
+theory/kernel/kernel_derivatives.pdf
+theory/kernel/kernel_definitions.pdf
m4/libtool.m4
m4/ltoptions.m4
diff --git a/doc/Doxyfile.in b/doc/Doxyfile.in
index 8f61a060b37b0e62189160d0a8c61e713cfd3b8f..802d8c31c251e006711934b6d30ace6c47eec4ac 100644
--- a/doc/Doxyfile.in
+++ b/doc/Doxyfile.in
@@ -759,7 +759,9 @@ WARN_LOGFILE =
# spaces.
# Note: If this tag is empty the current directory is searched.
-INPUT = @top_srcdir@ @top_srcdir@/src @top_srcdir@/src/hydro/Minimal @top_srcdir@/src/gravity/Default
+INPUT = @top_srcdir@ @top_srcdir@/src @top_srcdir@/tests @top_srcdir@/examples
+INPUT += @top_srcdir@/src/hydro/Minimal @top_srcdir@/src/gravity/Default
+INPUT += @top_srcdir@/src/riemann
# This tag can be used to specify the character encoding of the source files
# that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
diff --git a/examples/BigCosmoVolume/makeIC.py b/examples/BigCosmoVolume/makeIC.py
index 0994e1c95e053defe7766122c52bc405c7239776..3020feaf753f817f039d2fd09c4fa4f7fb69b896 100644
--- a/examples/BigCosmoVolume/makeIC.py
+++ b/examples/BigCosmoVolume/makeIC.py
@@ -77,7 +77,7 @@ indices = indices < numPart
coords = coords[indices,:]
v = v[indices,:]
m = m[indices]
-h = h[indices]
+h = h[indices] / 1.825742 # Correct from Gadget defintion of h to physical definition
u = u[indices]
ids = ids[indices]
diff --git a/examples/GreshoVortex/makeIC.py b/examples/GreshoVortex/makeIC.py
index 6aceeed559324f97b0b1e388ff0c3524498b52e4..12edcb6e8154ec6f865d28a6daeb02d385d14bbf 100644
--- a/examples/GreshoVortex/makeIC.py
+++ b/examples/GreshoVortex/makeIC.py
@@ -30,6 +30,7 @@ factor = 3
boxSize = [ 1.0 , 1.0, 1.0/factor ]
L = 120 # Number of particles along one axis
gamma = 5./3. # Gas adiabatic index
+eta = 1.2349 # 48 ngbs with cubic spline kernel
rho = 1 # Gas density
P0 = 0. # Constant additional pressure (should have no impact on the dynamics)
fileName = "greshoVortex.hdf5"
@@ -73,7 +74,7 @@ for i in range(L):
v[index,1] = v_phi * (x - boxSize[0] / 2) / r
v[index,2] = 0.
m[index] = mass
- h[index] = 2.251 * boxSize[0] / L
+ h[index] = eta * boxSize[0] / L
P = P0
if r < 0.2:
P = P + 5. + 12.5*r2
@@ -105,6 +106,14 @@ grp.attrs["Flag_Entropy_ICs"] = [0, 0, 0, 0, 0, 0]
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
#Particle group
grp = file.create_group("/PartType0")
ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
diff --git a/examples/MultiTypes/makeIC.py b/examples/MultiTypes/makeIC.py
index 3a41910c22c260086b5384b248a5c86ab6340a5e..cf889f9b6eab502f692cd6c8b4506c31664ecdcb 100644
--- a/examples/MultiTypes/makeIC.py
+++ b/examples/MultiTypes/makeIC.py
@@ -32,6 +32,7 @@ Lgas = int(sys.argv[1]) # Number of particles along one axis
rhoGas = 2. # Density
P = 1. # Pressure
gamma = 5./3. # Gas adiabatic index
+eta = 1.2349 # 48 ngbs with cubic spline kernel
rhoDM = 1.
Ldm = int(sys.argv[2]) # Number of particles along one axis
@@ -61,11 +62,18 @@ grp.attrs["NumFilesPerSnapshot"] = 1
grp.attrs["MassTable"] = [0.0, massDM, 0.0, 0.0, 0.0, 0.0]
grp.attrs["Flag_Entropy_ICs"] = 0
-
#Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
# Gas Particle group
grp = file.create_group("/PartType0")
@@ -80,7 +88,7 @@ ds = grp.create_dataset('Masses', (numGas,1), 'f')
ds[()] = m
m = zeros(1)
-h = full((numGas, 1), 1.1255 * boxSize / Lgas)
+h = full((numGas, 1), eta * boxSize / Lgas)
ds = grp.create_dataset('SmoothingLength', (numGas,1), 'f')
ds[()] = h
h = zeros(1)
diff --git a/examples/PerturbedBox/makeIC.py b/examples/PerturbedBox/makeIC.py
index 69c1a69199c9a5262f5ae6c4e95ca14699300fd4..ee1d845fc2149892909a54bf588046b0b1691b03 100644
--- a/examples/PerturbedBox/makeIC.py
+++ b/examples/PerturbedBox/makeIC.py
@@ -90,6 +90,14 @@ grp.attrs["NumPart_Total"] = numPart
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
#Particle group
grp = file.create_group("/PartType0")
ds = grp.create_dataset('Coordinates', (numPart, 3), 'd')
diff --git a/examples/SedovBlast/makeIC.py b/examples/SedovBlast/makeIC.py
index 75ff81165df51780848e3d8ac679a6dbeb17a039..e64942e8e92ee6fe67142f841f566019b1a668be 100644
--- a/examples/SedovBlast/makeIC.py
+++ b/examples/SedovBlast/makeIC.py
@@ -33,6 +33,7 @@ P = 1.e-5 # Pressure
E0= 1.e2 # Energy of the explosion
pert = 0.1
gamma = 5./3. # Gas adiabatic index
+eta = 1.2349 # 48 ngbs with cubic spline kernel
fileName = "sedov.hdf5"
@@ -67,7 +68,7 @@ for i in range(L):
v[index,1] = 0.
v[index,2] = 0.
m[index] = mass
- h[index] = 1.1255 * boxSize / L
+ h[index] = eta * boxSize / L
u[index] = internalEnergy
ids[index] = index + 1
if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 2.01 * boxSize/L:
@@ -98,6 +99,14 @@ grp.attrs["Flag_Entropy_ICs"] = 0
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
#Particle group
grp = file.create_group("/PartType0")
grp.create_dataset('Coordinates', data=coords, dtype='d')
diff --git a/examples/SedovBlast/makeIC_fcc.py b/examples/SedovBlast/makeIC_fcc.py
index 17f07440909cb5478d09a5b7a1444c72af2f3a47..0d3a017a9b7f3b30b61e723e3d1646d7797b40a4 100644
--- a/examples/SedovBlast/makeIC_fcc.py
+++ b/examples/SedovBlast/makeIC_fcc.py
@@ -33,6 +33,7 @@ P = 1.e-5 # Pressure
E0= 1.e2 # Energy of the explosion
pert = 0.025
gamma = 5./3. # Gas adiabatic index
+eta = 1.2349 # 48 ngbs with cubic spline kernel
fileName = "sedov.hdf5"
@@ -70,7 +71,7 @@ for i in range(L):
v[index,1] = 0.
v[index,2] = 0.
m[index] = mass
- h[index] = 1.1255 * hbox
+ h[index] = eta * hbox
u[index] = internalEnergy
ids[index] = index + 1
if sqrt((x - boxSize/2.)**2 + (y - boxSize/2.)**2 + (z - boxSize/2.)**2) < 1.2 * hbox:
@@ -101,6 +102,14 @@ grp.attrs["Flag_Entropy_ICs"] = 0
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
#Particle group
grp = file.create_group("/PartType0")
grp.create_dataset('Coordinates', data=coords, dtype='d')
diff --git a/examples/SodShock/makeIC.py b/examples/SodShock/makeIC.py
index 0ac0564116f8a6ceb57b4f41d23eb9907df0440d..8ae19050c11c0712579b44646c8870d7574d113b 100644
--- a/examples/SodShock/makeIC.py
+++ b/examples/SodShock/makeIC.py
@@ -43,14 +43,14 @@ vol = boxSize[0] * boxSize[1] * boxSize[2]
glass1 = h5py.File("glass_001.hdf5")
pos1 = glass1["/PartType0/Coordinates"][:,:]
pos1 = pos1 / factor # Particles are in [0:0.25, 0:0.25, 0:0.25]
-
+glass_h1 = glass1["/PartType0/SmoothingLength"][:] / factor
#Read in high density glass
# glass2 = h5py.File("../Glass/glass_50000.hdf5")
glass2 = h5py.File("glass_002.hdf5")
pos2 = glass2["/PartType0/Coordinates"][:,:]
pos2 = pos2 / factor # Particles are in [0:0.25, 0:0.25, 0:0.25]
-
+glass_h2 = glass2["/PartType0/SmoothingLength"][:] / factor
#Generate high density region
rho1 = 1.
@@ -61,9 +61,10 @@ coord1 = append(coord1, coord1 + [0.25, 0, 0], 0)
# coord1 = append(coord1, pos1 + [0, 0.5, 0.5], 0)
N1 = size(coord1)/3
v1 = zeros((N1, 3))
-h1 = ones(N1) * 2.251 * 0.5 * vol / (size(pos1)/3)**(1./3.)
u1 = ones(N1) * P1 / ((gamma - 1.) * rho1)
m1 = ones(N1) * vol * 0.5 * rho1 / N1
+h1 = append(glass_h1, glass_h1, 0)
+h1 = append(h1, h1, 0)
#Generate low density region
rho2 = 0.25
@@ -74,9 +75,10 @@ coord2 = append(coord2, coord2 + [0.25, 0, 0], 0)
# coord2 = append(coord2, pos2 + [0, 0.5, 0.5], 0)
N2 = size(coord2)/3
v2 = zeros((N2, 3))
-h2 = ones(N2) * 2.251 * 0.5 * vol / (size(pos2)/3)**(1./3.)
u2 = ones(N2) * P2 / ((gamma - 1.) * rho2)
m2 = ones(N2) * vol * 0.5 * rho2 / N2
+h2 = append(glass_h2, glass_h2, 0)
+h2 = append(h2, h2, 0)
#Merge arrays
numPart = N1 + N2
@@ -89,8 +91,8 @@ ids = zeros(numPart, dtype='L')
for i in range(1, numPart+1):
ids[i-1] = i
-#Final operations
-h /= 2
+#Final operation since we come from Gadget-2 cubic spline ICs
+h /= 1.825752
#File
file = h5py.File(fileName, 'w')
@@ -110,6 +112,14 @@ grp.attrs["Flag_Entropy_ICs"] = 0
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
#Particle group
grp = file.create_group("/PartType0")
grp.create_dataset('Coordinates', data=coords, dtype='d')
diff --git a/examples/UniformBox/makeIC.py b/examples/UniformBox/makeIC.py
index c175349e658799cbcb30dfe2619a1594bafc18b9..1484f60596e68734f0f98685ab2ab845f2e0b407 100644
--- a/examples/UniformBox/makeIC.py
+++ b/examples/UniformBox/makeIC.py
@@ -32,6 +32,7 @@ L = int(sys.argv[1]) # Number of particles along one axis
rho = 2. # Density
P = 1. # Pressure
gamma = 5./3. # Gas adiabatic index
+eta = 1.2349 # 48 ngbs with cubic spline kernel
fileName = "uniformBox.hdf5"
#---------------------------------------------------
@@ -55,11 +56,18 @@ grp.attrs["NumFilesPerSnapshot"] = 1
grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
grp.attrs["Flag_Entropy_ICs"] = 0
-
#Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
#Particle group
grp = file.create_group("/PartType0")
@@ -73,7 +81,7 @@ ds = grp.create_dataset('Masses', (numPart,1), 'f')
ds[()] = m
m = zeros(1)
-h = full((numPart, 1), 1.1255 * boxSize / L)
+h = full((numPart, 1), eta * boxSize / L)
ds = grp.create_dataset('SmoothingLength', (numPart,1), 'f')
ds[()] = h
h = zeros(1)
diff --git a/examples/UniformBox/makeICbig.py b/examples/UniformBox/makeICbig.py
index e475fdcbd9f3c4811e3dcfdf20bbd321be3d8b29..bd5cf627fb535595b3abb224cbc8de50589f12cf 100644
--- a/examples/UniformBox/makeICbig.py
+++ b/examples/UniformBox/makeICbig.py
@@ -32,6 +32,7 @@ N = int(sys.argv[2]) # Write N particles at a time to avoid requiring a lot of
rho = 2. # Density
P = 1. # Pressure
gamma = 5./3. # Gas adiabatic index
+eta = 1.2349 # 48 ngbs with cubic spline kernel
fileName = "uniformBox_%d.hdf5"%L
#---------------------------------------------------
@@ -62,11 +63,19 @@ grp.attrs["NumFilesPerSnapshot"] = 1
grp.attrs["MassTable"] = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
grp.attrs["Flag_Entropy_ICs"] = 0
-
#Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+
#Particle group
grp = file.create_group("/PartType0")
@@ -89,7 +98,7 @@ for n in range(n_iterations):
ds_m[offset:offset+N] = m
m = zeros(1)
- h = full((N, 1), 1.1255 * boxSize / L)
+ h = full((N, 1), eta * boxSize / L)
ds_h[offset:offset+N] = h
h = zeros(1)
@@ -122,7 +131,7 @@ m = full((remainder, 1), mass)
ds_m[offset:offset+remainder] = m
m = zeros(1)
-h = full((remainder, 1), 1.1255 * boxSize / L)
+h = full((remainder, 1), eta * boxSize / L)
ds_h[offset:offset+remainder] = h
h = zeros(1)
@@ -139,7 +148,7 @@ coords = zeros((remainder, 3))
coords[:,0] = z[:,0] * boxSize / L + boxSize / (2*L)
coords[:,1] = y[:,0] * boxSize / L + boxSize / (2*L)
coords[:,2] = x[:,0] * boxSize / L + boxSize / (2*L)
-ds_x[offset:offset+remainder,:] = coords
+ods_x[offset:offset+remainder,:] = coords
print "Done", offset+remainder,"/", numPart
diff --git a/examples/UniformDMBox/makeIC.py b/examples/UniformDMBox/makeIC.py
index 061b4d0ad1959d9e25356aff80e78adb9c1c4faa..449d780fb31bc23dd194f772be45d35e6b0bbe3f 100644
--- a/examples/UniformDMBox/makeIC.py
+++ b/examples/UniformDMBox/makeIC.py
@@ -52,11 +52,19 @@ grp.attrs["NumFilesPerSnapshot"] = 1
grp.attrs["MassTable"] = [0.0, mass, 0.0, 0.0, 0.0, 0.0]
grp.attrs["Flag_Entropy_ICs"] = 0
-
#Runtime parameters
grp = file.create_group("/RuntimePars")
grp.attrs["PeriodicBoundariesOn"] = periodic
+#Units
+grp = file.create_group("/Units")
+grp.attrs["Unit length in cgs (U_L)"] = 1.
+grp.attrs["Unit mass in cgs (U_M)"] = 1.
+grp.attrs["Unit time in cgs (U_t)"] = 1.
+grp.attrs["Unit current in cgs (U_I)"] = 1.
+grp.attrs["Unit temperature in cgs (U_T)"] = 1.
+
+
#Particle group
grp = file.create_group("/PartType1")
diff --git a/examples/main.c b/examples/main.c
index f770265c1b59dff93513a59ff7246258539325e3..9523af49ed30c54d256d287ea2846a854650bc05 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -1,4 +1,3 @@
-
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
@@ -56,7 +55,6 @@
* @brief Main routine that loads a few particles and generates some output.
*
*/
-
int main(int argc, char *argv[]) {
int c, icount, periodic = 1;
@@ -80,7 +78,8 @@ int main(int argc, char *argv[]) {
int nr_nodes = 1, myrank = 0;
FILE *file_thread;
int with_outputs = 1;
- int with_gravity = 0;
+ int with_external_gravity = 0;
+ int with_self_gravity = 0;
int engine_policies = 0;
int verbose = 0, talking = 0;
unsigned long long cpufreq = 0;
@@ -159,7 +158,7 @@ int main(int argc, char *argv[]) {
bzero(&s, sizeof(struct space));
/* Parse the options */
- while ((c = getopt(argc, argv, "a:c:d:e:f:gh:m:oP:q:R:s:t:v:w:y:z:")) != -1)
+ while ((c = getopt(argc, argv, "a:c:d:e:f:gGh:m:oP:q:R:s:t:v:w:y:z:")) != -1)
switch (c) {
case 'a':
if (sscanf(optarg, "%lf", &scaling) != 1)
@@ -189,7 +188,10 @@ int main(int argc, char *argv[]) {
if (!strcpy(ICfileName, optarg)) error("Error parsing IC file name.");
break;
case 'g':
- with_gravity = 1;
+ with_external_gravity = 1;
+ break;
+ case 'G':
+ with_self_gravity = 1;
break;
case 'h':
if (sscanf(optarg, "%llu", &cpufreq) != 1)
@@ -349,10 +351,6 @@ int main(int argc, char *argv[]) {
message("CPU frequency used for tick conversion: %llu Hz", cpufreq);
}
- /* Check we have sensible time step bounds */
- if (dt_min > dt_max)
- error("Minimal time step size must be large than maximal time step size ");
-
/* Check whether an IC file has been provided */
if (strcmp(ICfileName, "") == 0)
error("An IC file name must be provided via the option -f");
@@ -394,7 +392,7 @@ int main(int argc, char *argv[]) {
#endif
/* MATTHIEU: Temporary fix to preserve master */
- if (!with_gravity) {
+ if (!with_external_gravity && !with_self_gravity) {
free(gparts);
gparts = NULL;
for (size_t k = 0; k < Ngas; ++k) parts[k].gpart = NULL;
@@ -482,7 +480,8 @@ int main(int argc, char *argv[]) {
/* Construct the engine policy */
engine_policies = ENGINE_POLICY | engine_policy_steal | engine_policy_hydro;
- if (with_gravity) engine_policies |= engine_policy_external_gravity;
+ if (with_external_gravity) engine_policies |= engine_policy_external_gravity;
+ if (with_self_gravity) engine_policies |= engine_policy_self_gravity;
/* Initialize the engine with this space. */
if (myrank == 0) clocks_gettime(&tic);
@@ -545,8 +544,8 @@ int main(int argc, char *argv[]) {
/* Legend */
if (myrank == 0)
printf(
- "# Step Time time-step Number of updates CPU Wall-clock time "
- "[%s]\n",
+ "# Step Time time-step Number of updates Number of updates "
+ "CPU Wall-clock time [%s]\n",
clocks_getunit());
/* Let loose a runner on the space. */
diff --git a/examples/plot_tasks.py b/examples/plot_tasks.py
index eaff41ebae1bad0f1307d23a3204186ecbc63b2f..895c32ef9c3d1490e6d30b7dc79e40171a228ee9 100755
--- a/examples/plot_tasks.py
+++ b/examples/plot_tasks.py
@@ -60,7 +60,7 @@ pl.rcParams.update(PLOT_PARAMS)
# Tasks and subtypes. Indexed as in tasks.h.
TASKTYPES = ["none", "sort", "self", "pair", "sub", "init", "ghost", "drift", "kick",
"send", "recv", "grav_pp", "grav_mm", "grav_up", "grav_down",
- "psort", "split_cell", "rewait", "count"]
+ "part_sort", "gpart_sort", "split_cell", "rewait", "count"]
TASKCOLOURS = {"none": "black",
"sort": "lightblue",
@@ -77,7 +77,8 @@ TASKCOLOURS = {"none": "black",
"grav_mm": "mediumturquoise",
"grav_up": "mediumvioletred",
"grav_down": "mediumnightblue",
- "psort": "steelblue",
+ "part_sort": "steelblue",
+ "gpart_sort": "teal" ,
"split_cell": "seagreen",
"rewait": "olive",
"count": "powerblue"}
diff --git a/examples/plot_tasks_MPI.py b/examples/plot_tasks_MPI.py
index b7d1823ad746d6a10b5e67fc9f7315b13be4649f..d59fe6417b524b8cb3cf8f6117fca3b8b3f3c780 100755
--- a/examples/plot_tasks_MPI.py
+++ b/examples/plot_tasks_MPI.py
@@ -66,7 +66,7 @@ pl.rcParams.update(PLOT_PARAMS)
# Tasks and subtypes. Indexed as in tasks.h.
TASKTYPES = ["none", "sort", "self", "pair", "sub", "init", "ghost", "drift", "kick",
"send", "recv", "grav_pp", "grav_mm", "grav_up", "grav_down",
- "psort", "split_cell", "rewait", "count"]
+ "part_sort", "gpart_sort", "split_cell", "rewait", "count"]
TASKCOLOURS = {"none": "black",
"sort": "lightblue",
@@ -83,7 +83,8 @@ TASKCOLOURS = {"none": "black",
"grav_mm": "mediumturquoise",
"grav_up": "mediumvioletred",
"grav_down": "mediumnightblue",
- "psort": "steelblue",
+ "part_sort": "steelblue",
+ "gpart_sort": "teal",
"split_cell": "seagreen",
"rewait": "olive",
"count": "powerblue"}
diff --git a/src/Makefile.am b/src/Makefile.am
index f44d47819672d10445fd969fe2ff20dbcb49463b..a96f35b3cf0d8a23aec4f8c0f8d16bec8638cbcd 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -35,17 +35,17 @@ endif
# List required headers
include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
engine.h swift.h serial_io.h timers.h debug.h scheduler.h proxy.h parallel_io.h \
- common_io.h single_io.h multipole.h map.h tools.h partition.h clocks.h
+ common_io.h single_io.h multipole.h map.h tools.h partition.h clocks.h parser.h
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
serial_io.c timers.c debug.c scheduler.c proxy.c parallel_io.c \
units.c common_io.c single_io.c multipole.c version.c map.c \
- kernel.c tools.c part.c partition.c clocks.c
+ kernel_hydro.c kernel_gravity.c tools.c part.c partition.c clocks.c parser.c
# Include files for distribution, not installation.
-nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel.h vector.h \
- runner_doiact.h runner_doiact_grav.h units.h intrinsics.h minmax.h \
+nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
+ vector.h runner_doiact.h runner_doiact_grav.h units.h intrinsics.h minmax.h \
gravity.h gravity_io.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
diff --git a/src/cell.c b/src/cell.c
index b0e299999e9073cc14d92b769b3e157f50d64c8f..61acfaaea7a0af01a78ab773541564e9a2723f4e 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -45,6 +45,7 @@
/* Local headers. */
#include "atomic.h"
#include "error.h"
+#include "gravity.h"
#include "hydro.h"
#include "space.h"
#include "timers.h"
@@ -605,6 +606,27 @@ void cell_init_parts(struct cell *c, void *data) {
c->ti_end_max = 0;
}
+/**
+ * @brief Initialises all g-particles to a valid state even if the ICs were
+ *stupid
+ *
+ * @param c Cell to act upon
+ * @param data Unused parameter
+ */
+void cell_init_gparts(struct cell *c, void *data) {
+
+ struct gpart *gp = c->gparts;
+ const int gcount = c->gcount;
+
+ for (int i = 0; i < gcount; ++i) {
+ gp[i].ti_begin = 0;
+ gp[i].ti_end = 0;
+ gravity_first_init_gpart(&gp[i]);
+ }
+ c->ti_end_min = 0;
+ c->ti_end_max = 0;
+}
+
/**
* @brief Converts hydro quantities to a valid state after the initial density
*calculation
diff --git a/src/cell.h b/src/cell.h
index 05b870b144dd425ee9e72ee3b529122b6e721523..a471eac44bfd3533c4220ab8c5ff2ddec724e87f 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -141,7 +141,7 @@ struct cell {
double mass, e_pot, e_int, e_kin;
/* Number of particles updated in this cell. */
- int updated;
+ int updated, g_updated;
/* Linking pointer for "memory management". */
struct cell *next;
@@ -178,6 +178,7 @@ int cell_getsize(struct cell *c);
int cell_link_parts(struct cell *c, struct part *parts);
int cell_link_gparts(struct cell *c, struct gpart *gparts);
void cell_init_parts(struct cell *c, void *data);
+void cell_init_gparts(struct cell *c, void *data);
void cell_convert_hydro(struct cell *c, void *data);
void cell_clean_links(struct cell *c, void *data);
diff --git a/src/common_io.c b/src/common_io.c
index 6183effe9ce392ab930c581cbd118f025bbce773..bc981ecab31c56925ce08bfafb1a3a16aeee104b 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -42,7 +42,7 @@
/* Local includes. */
#include "const.h"
#include "error.h"
-#include "kernel.h"
+#include "kernel_hydro.h"
#include "version.h"
const char* particle_type_names[NUM_PARTICLE_TYPES] = {
diff --git a/src/common_io.h b/src/common_io.h
index 961f40e63d771e5e06ade525301caf59aae0bceb..b7f3a1a317d69937dde8692eead8f00c75649477 100644
--- a/src/common_io.h
+++ b/src/common_io.h
@@ -24,6 +24,7 @@
#include "../config.h"
/* Includes. */
+#include "kernel_hydro.h"
#include "part.h"
#include "units.h"
diff --git a/src/debug.c b/src/debug.c
index 4c1434118c98aab7def28d3a53493767d249d774..53a03d66aee2c169a555ed00a2efa2d5b984066a 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -60,7 +60,7 @@ void printParticle(struct part *parts, struct xpart *xparts, long long int id,
/* Look for the particle. */
for (size_t i = 0; i < N; i++)
if (parts[i].id == id) {
- printf("## Particle[%zd]:\n id=%lld", i, parts[i].id);
+ printf("## Particle[%zd]:\n id=%lld ", i, parts[i].id);
hydro_debug_particle(&parts[i], &xparts[i]);
found = 1;
break;
@@ -76,12 +76,12 @@ void printgParticle(struct gpart *gparts, long long int id, size_t N) {
/* Look for the particle. */
for (size_t i = 0; i < N; i++)
if (gparts[i].id == -id) {
- printf("## gParticle[%zd] (DM) :\n id=%lld", i, -gparts[i].id);
+ printf("## gParticle[%zd] (DM) :\n id=%lld ", i, -gparts[i].id);
gravity_debug_particle(&gparts[i]);
found = 1;
break;
} else if (gparts[i].id > 0 && gparts[i].part->id == id) {
- printf("## gParticle[%zd] (hydro) :\n id=%lld", i, gparts[i].id);
+ printf("## gParticle[%zd] (hydro) :\n id=%lld ", i, gparts[i].id);
gravity_debug_particle(&gparts[i]);
found = 1;
break;
diff --git a/src/engine.c b/src/engine.c
index d1bd044b791ca6fa9a9106dfdffb404839066b6f..51d4c617d0bf97beb8ebc305fec24fe406cf15d3 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -256,7 +256,7 @@ void engine_redistribute(struct engine *e) {
}
/* Sort the gparticles according to their cell index. */
- space_gparts_sort(gparts, g_dest, s->nr_gparts, 0, nr_nodes - 1);
+ space_gparts_sort(s, g_dest, s->nr_gparts, 0, nr_nodes - 1, e->verbose);
/* Get all the counts from all the nodes. */
if (MPI_Allreduce(MPI_IN_PLACE, counts, nr_nodes * nr_nodes, MPI_INT, MPI_SUM,
@@ -441,10 +441,10 @@ void engine_redistribute(struct engine *e) {
if (gparts_new[k].part->gpart != &gparts_new[k])
error("Linking problem !");
- /* if (gparts_new[k].x[0] != gparts_new[k].part->x[0] || */
- /* gparts_new[k].x[1] != gparts_new[k].part->x[1] || */
- /* gparts_new[k].x[2] != gparts_new[k].part->x[2]) */
- /* error("Linked particles are not at the same position !"); */
+ if (gparts_new[k].x[0] != gparts_new[k].part->x[0] ||
+ gparts_new[k].x[1] != gparts_new[k].part->x[1] ||
+ gparts_new[k].x[2] != gparts_new[k].part->x[2])
+ error("Linked particles are not at the same position !");
}
}
for (size_t k = 0; k < nr_parts; ++k) {
@@ -1365,9 +1365,12 @@ void engine_maketasks(struct engine *e) {
scheduler_reset(sched, s->tot_cells * engine_maxtaskspercell);
/* Add the space sorting tasks. */
- for (int i = 0; i < e->nr_threads; i++)
- scheduler_addtask(sched, task_type_psort, task_subtype_none, i, 0, NULL,
+ for (int i = 0; i < e->nr_threads; i++) {
+ scheduler_addtask(sched, task_type_part_sort, task_subtype_none, i, 0, NULL,
NULL, 0);
+ scheduler_addtask(sched, task_type_gpart_sort, task_subtype_none, i, 0,
+ NULL, NULL, 0);
+ }
/* Construct the firt hydro loop over neighbours */
engine_make_hydroloop_tasks(e);
@@ -1557,6 +1560,7 @@ int engine_marktasks(struct engine *e) {
else if (t->type == task_type_kick) {
t->skip = (t->ci->ti_end_min > ti_end);
t->ci->updated = 0;
+ t->ci->g_updated = 0;
}
/* Drift? */
@@ -1613,6 +1617,7 @@ void engine_print_task_counts(struct engine *e) {
printf(" skipped=%i ]\n", counts[task_type_count]);
fflush(stdout);
message("nr_parts = %zi.", e->s->nr_parts);
+ message("nr_gparts = %zi.", e->s->nr_gparts);
}
/**
@@ -1742,7 +1747,7 @@ void engine_collect_kick(struct cell *c) {
if (c->kick != NULL) return;
/* Counters for the different quantities. */
- int updated = 0;
+ int updated = 0, g_updated = 0;
double e_kin = 0.0, e_int = 0.0, e_pot = 0.0;
float mom[3] = {0.0f, 0.0f, 0.0f}, ang[3] = {0.0f, 0.0f, 0.0f};
int ti_end_min = max_nr_timesteps, ti_end_max = 0;
@@ -1765,6 +1770,7 @@ void engine_collect_kick(struct cell *c) {
ti_end_min = min(ti_end_min, cp->ti_end_min);
ti_end_max = max(ti_end_max, cp->ti_end_max);
updated += cp->updated;
+ g_updated += cp->g_updated;
e_kin += cp->e_kin;
e_int += cp->e_int;
e_pot += cp->e_pot;
@@ -1782,6 +1788,7 @@ void engine_collect_kick(struct cell *c) {
c->ti_end_min = ti_end_min;
c->ti_end_max = ti_end_max;
c->updated = updated;
+ c->g_updated = g_updated;
c->e_kin = e_kin;
c->e_int = e_int;
c->e_pot = e_pot;
@@ -1845,7 +1852,15 @@ void engine_init_particles(struct engine *e) {
/* Make sure all particles are ready to go */
/* i.e. clean-up any stupid state in the ICs */
- space_map_cells_pre(s, 1, cell_init_parts, NULL);
+ if ((e->policy & engine_policy_hydro) == engine_policy_hydro) {
+ space_map_cells_pre(s, 1, cell_init_parts, NULL);
+ }
+ if (((e->policy & engine_policy_self_gravity) ==
+ engine_policy_self_gravity) ||
+ ((e->policy & engine_policy_external_gravity) ==
+ engine_policy_external_gravity)) {
+ space_map_cells_pre(s, 1, cell_init_gparts, NULL);
+ }
engine_prepare(e);
@@ -1919,7 +1934,7 @@ void engine_init_particles(struct engine *e) {
*/
void engine_step(struct engine *e) {
- int updates = 0;
+ int updates = 0, g_updates = 0;
int ti_end_min = max_nr_timesteps, ti_end_max = 0;
double e_pot = 0.0, e_int = 0.0, e_kin = 0.0;
float mom[3] = {0.0, 0.0, 0.0};
@@ -1946,6 +1961,7 @@ void engine_step(struct engine *e) {
e_int += c->e_int;
e_pot += c->e_pot;
updates += c->updated;
+ g_updates += c->g_updated;
mom[0] += c->mom[0];
mom[1] += c->mom[1];
mom[2] += c->mom[2];
@@ -1971,18 +1987,20 @@ void engine_step(struct engine *e) {
ti_end_max = in_i[0];
}
{
- double in_d[4], out_d[4];
+ double in_d[5], out_d[5];
out_d[0] = updates;
- out_d[1] = e_kin;
- out_d[2] = e_int;
- out_d[3] = e_pot;
- if (MPI_Allreduce(out_d, in_d, 4, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) !=
+ out_d[1] = g_updates;
+ out_d[2] = e_kin;
+ out_d[3] = e_int;
+ out_d[4] = e_pot;
+ if (MPI_Allreduce(out_d, in_d, 5, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD) !=
MPI_SUCCESS)
error("Failed to aggregate energies.");
updates = in_d[0];
- e_kin = in_d[1];
- e_int = in_d[2];
- e_pot = in_d[3];
+ g_updates = in_d[1];
+ e_kin = in_d[2];
+ e_int = in_d[3];
+ e_pot = in_d[4];
}
#endif
@@ -2007,8 +2025,8 @@ void engine_step(struct engine *e) {
if (e->nodeID == 0) {
/* Print some information to the screen */
- printf("%d %e %e %d %.3f\n", e->step, e->time, e->timeStep, updates,
- e->wallclock_time);
+ printf("%d %e %e %d %d %.3f\n", e->step, e->time, e->timeStep, updates,
+ g_updates, e->wallclock_time);
fflush(stdout);
/* Write some energy statistics */
@@ -2250,6 +2268,28 @@ void engine_split(struct engine *e, struct partition *initial_partition) {
for (size_t k = 0; k < s->nr_gparts; k++)
if (s->gparts[k].id > 0) s->gparts[k].part->gpart = &s->gparts[k];
+ /* Verify that the links are correct */
+ /* MATTHIEU: To be commented out once we are happy */
+ for (size_t k = 0; k < s->nr_gparts; ++k) {
+
+ if (s->gparts[k].id > 0) {
+
+ if (s->gparts[k].part->gpart != &s->gparts[k]) error("Linking problem !");
+
+ if (s->gparts[k].x[0] != s->gparts[k].part->x[0] ||
+ s->gparts[k].x[1] != s->gparts[k].part->x[1] ||
+ s->gparts[k].x[2] != s->gparts[k].part->x[2])
+ error("Linked particles are not at the same position !");
+ }
+ }
+ for (size_t k = 0; k < s->nr_parts; ++k) {
+
+ if (s->parts[k].gpart != NULL) {
+
+ if (s->parts[k].gpart->part != &s->parts[k]) error("Linking problem !");
+ }
+ }
+
#else
error("SWIFT was not compiled with MPI support.");
#endif
@@ -2425,6 +2465,7 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
/* Print information about the hydro scheme */
if (e->nodeID == 0) message("Hydrodynamic scheme: %s", SPH_IMPLEMENTATION);
+ if (e->nodeID == 0) message("Hydrodynamic kernel: %s", kernel_name);
/* Check we have sensible time bounds */
if (timeBegin >= timeEnd)
@@ -2433,10 +2474,12 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
"(t_beg = %e)",
timeEnd, timeBegin);
- /* Check we have sensible time step bounds */
+ /* Check we have sensible time-step values */
if (e->dt_min > e->dt_max)
error(
- "Minimal time step size must be smaller than maximal time step size ");
+ "Minimal time-step size (%e) must be smaller than maximal time-step "
+ "size (%e)",
+ e->dt_min, e->dt_max);
/* Deal with timestep */
e->timeBase = (timeEnd - timeBegin) / max_nr_timesteps;
@@ -2503,9 +2546,13 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
s->nr_queues = nr_queues;
/* Create the sorting tasks. */
- for (int i = 0; i < e->nr_threads; i++)
- scheduler_addtask(&e->sched, task_type_psort, task_subtype_none, i, 0, NULL,
- NULL, 0);
+ for (int i = 0; i < e->nr_threads; i++) {
+ scheduler_addtask(&e->sched, task_type_part_sort, task_subtype_none, i, 0,
+ NULL, NULL, 0);
+
+ scheduler_addtask(&e->sched, task_type_gpart_sort, task_subtype_none, i, 0,
+ NULL, NULL, 0);
+ }
scheduler_ranktasks(&e->sched);
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 82bc52ad3e05794c8c05896075edc463a69197ff..92a9f64c1f84a9e949f4c0e9485f892b5c808cdc 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -22,14 +22,61 @@
/**
* @brief Computes the gravity time-step of a given particle
*
- * @param p Pointer to the particle data
- * @param xp Pointer to the extended particle data
+ * @param gp Pointer to the g-particle data
*
*/
-__attribute__((always_inline)) INLINE static float gravity_compute_timestep(
- struct part* p, struct xpart* xp) {
+__attribute__((always_inline))
+ INLINE static float gravity_compute_timestep(struct gpart* gp) {
/* Currently no limit is imposed */
return FLT_MAX;
}
+
+/**
+ * @brief Initialises the g-particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * @param gp The particle to act upon
+ */
+__attribute__((always_inline))
+ INLINE static void gravity_first_init_gpart(struct gpart* gp) {}
+
+/**
+ * @brief Prepares a g-particle for the gravity calculation
+ *
+ * Zeroes all the relevant arrays in preparation for the sums taking place in
+ * the variaous tasks
+ *
+ * @param gp The particle to act upon
+ */
+__attribute__((always_inline))
+ INLINE static void gravity_init_part(struct gpart* gp) {
+
+ /* Zero the acceleration */
+ gp->a_grav[0] = 0.f;
+ gp->a_grav[1] = 0.f;
+ gp->a_grav[2] = 0.f;
+}
+
+/**
+ * @brief Finishes the gravity calculation.
+ *
+ * Multiplies the forces and accelerations by the appropiate constants
+ *
+ * @param gp The particle to act upon
+ */
+__attribute__((always_inline))
+ INLINE static void gravity_end_force(struct gpart* gp) {}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * @param gp The particle to act upon
+ * @param dt The time-step for this kick
+ * @param half_dt The half time-step for this kick
+ */
+__attribute__((always_inline)) INLINE static void gravity_kick_extra(
+ struct gpart* gp, float dt, float half_dt) {}
diff --git a/src/gravity/Default/gravity_debug.h b/src/gravity/Default/gravity_debug.h
index 98e0c40a5700b4da70f27fb0955592bb5d2287c3..654745bfeb70dddba772af9e23797713376377a7 100644
--- a/src/gravity/Default/gravity_debug.h
+++ b/src/gravity/Default/gravity_debug.h
@@ -24,5 +24,5 @@ __attribute__((always_inline))
"v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
"mass=%.3e t_begin=%d, t_end=%d\n",
p->x[0], p->x[1], p->x[2], p->v_full[0], p->v_full[1], p->v_full[2],
- p->a[0], p->a[1], p->a[2], p->mass, p->ti_begin, p->ti_end);
+ p->a_grav[0], p->a_grav[1], p->a_grav[2], p->mass, p->ti_begin, p->ti_end);
}
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index 37c9e912fc38cddc1dac3bbd8c69cb11ae219f0d..62023345f174eb8cb9bae4d4438bdd50c9969494 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -22,7 +22,7 @@
/* Includes. */
#include "const.h"
-#include "kernel.h"
+#include "kernel_gravity.h"
#include "vector.h"
/**
@@ -49,8 +49,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav(
/* Aggregate the accelerations. */
for (k = 0; k < 3; k++) {
w = acc * dx[k];
- pi->a[k] -= w * mj;
- pj->a[k] += w * mi;
+ pi->a_grav[k] -= w * mj;
+ pj->a_grav[k] += w * mi;
}
}
@@ -100,8 +100,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_vec_grav(
ai.v = w.v * mj.v;
aj.v = w.v * mi.v;
for (j = 0; j < VEC_SIZE; j++) {
- pi[j]->a[k] -= ai.f[j];
- pj[j]->a[k] += aj.f[j];
+ pi[j]->a_grav[k] -= ai.f[j];
+ pj[j]->a_grav[k] += aj.f[j];
}
}
diff --git a/src/gravity/Default/gravity_part.h b/src/gravity/Default/gravity_part.h
index aaaab2055d8ce05d77a37d3f40e096bd91d87926..0dfdb82e4ec11c9153f77439027d7e4451ded7f4 100644
--- a/src/gravity/Default/gravity_part.h
+++ b/src/gravity/Default/gravity_part.h
@@ -29,7 +29,7 @@ struct gpart {
float v_full[3];
/* Particle acceleration. */
- float a[3];
+ float a_grav[3];
/* Particle mass. */
float mass;
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index fca4a346047d7dce0741924a69e95fdad5a5ce45..03953b07ad4e172d96b6e3382814e036a538e2bd 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -91,13 +91,16 @@ __attribute__((always_inline))
const float ih2 = ih * ih;
const float ih4 = ih2 * ih2;
- /* Final operation on the density. */
- p->rho = ih * ih2 * (p->rho + p->mass * kernel_root);
- p->rho_dh = (p->rho_dh - 3.0f * p->mass * kernel_root) * ih4;
- p->density.wcount =
- (p->density.wcount + kernel_root) * (4.0f / 3.0 * M_PI * kernel_gamma3);
- p->density.wcount_dh =
- p->density.wcount_dh * ih * (4.0f / 3.0 * M_PI * kernel_gamma3);
+ /* Final operation on the density (add self-contribution). */
+ p->rho += p->mass * kernel_root;
+ p->rho_dh -= 3.0f * p->mass * kernel_root * kernel_igamma;
+ p->density.wcount += kernel_root;
+
+ /* Finish the calculation by inserting the missing h-factors */
+ p->rho *= ih * ih2;
+ p->rho_dh *= ih4;
+ p->density.wcount *= (4.0f / 3.0f * M_PI * kernel_gamma3);
+ p->density.wcount_dh *= ih * (4.0f / 3.0f * M_PI * kernel_gamma4);
}
/**
diff --git a/src/hydro/Default/hydro_iact.h b/src/hydro/Default/hydro_iact.h
index b5b631501b2f9c398cf1f7e5ee32fd5c962ba86e..4f85299b9d61b3a66389bac3527a63068ab96db9 100644
--- a/src/hydro/Default/hydro_iact.h
+++ b/src/hydro/Default/hydro_iact.h
@@ -22,7 +22,7 @@
/* Includes. */
#include "const.h"
-#include "kernel.h"
+#include "kernel_hydro.h"
#include "part.h"
#include "vector.h"
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 8cc553363122099c748e3e3e1941611e986c8581..22c5734ed5762400285521b30f9aa60795c45325 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -101,7 +101,7 @@ __attribute__((always_inline))
p->rho *= ih * ih2;
p->rho_dh *= ih4;
p->density.wcount *= (4.0f / 3.0f * M_PI * kernel_gamma3);
- p->density.wcount_dh *= ih * (4.0f / 3.0f * M_PI * kernel_gamma3);
+ p->density.wcount_dh *= ih * (4.0f / 3.0f * M_PI * kernel_gamma4);
const float irho = 1.f / p->rho;
diff --git a/src/hydro/Gadget2/hydro_debug.h b/src/hydro/Gadget2/hydro_debug.h
index 46e156bb99015069f9958aeea05954e2be6db5e0..a4d1f7dd4397ebfc850b582e1ca81fc0d4edb76a 100644
--- a/src/hydro/Gadget2/hydro_debug.h
+++ b/src/hydro/Gadget2/hydro_debug.h
@@ -23,13 +23,13 @@ __attribute__((always_inline))
"x=[%.3e,%.3e,%.3e], "
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
"h=%.3e, "
- "wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, P=%.3e, S=%.3e, "
+ "wcount=%d, wcount_dh=%.3e, m=%.3e, dh_drho=%.3e, rho=%.3e, P=%.3e, S=%.3e, "
"dS/dt=%.3e, c=%.3e\n"
"divV=%.3e, curlV=%.3e, rotV=[%.3e,%.3e,%.3e] \n "
"v_sig=%e dh/dt=%.3e t_begin=%d, t_end=%d\n",
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
- p->h, (int)p->density.wcount, p->mass, p->rho_dh, p->rho,
+ p->h, (int)p->density.wcount, p->density.wcount_dh, p->mass, p->rho_dh, p->rho,
p->force.pressure, p->entropy, p->entropy_dt, p->force.soundspeed,
p->div_v, p->force.curl_v, p->density.rot_v[0], p->density.rot_v[1],
p->density.rot_v[2], p->force.v_sig, p->h_dt, p->ti_begin, p->ti_end);
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 09f796a8f37a9c015135f4aab3f821c2e862bdc9..d988c678affcf4ca722a965a7e52a7c120b4a924 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -22,7 +22,7 @@
/* Includes. */
#include "const.h"
-#include "kernel.h"
+#include "kernel_hydro.h"
#include "part.h"
#include "vector.h"
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index f4e3f1a70625430d9bd891c5f7596d71e7b8b231..7db3c275ce7e3389610e8297c287cbd5301c6c64 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -101,7 +101,12 @@ __attribute__((always_inline))
p->rho *= ih * ih2;
p->rho_dh *= ih4;
p->density.wcount *= (4.0f / 3.0f * M_PI * kernel_gamma3);
- p->density.wcount_dh *= ih * (4.0f / 3.0f * M_PI * kernel_gamma3);
+ p->density.wcount_dh *= ih * (4.0f / 3.0f * M_PI * kernel_gamma4);
+
+ const float irho = 1.f / p->rho;
+
+ /* Compute the derivative term */
+ p->rho_dh = 1.f / (1.f + 0.33333333f * p->h * p->rho_dh * irho);
}
/**
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index b3b81a9a0dfe41e7bfafe51050d6f7cf7157e31c..3427ec538613842f8fbcf0d8ba5f9ba5a0b8d540 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -21,7 +21,7 @@
/* Includes. */
#include "const.h"
-#include "kernel.h"
+#include "kernel_hydro.h"
#include "part.h"
#include "vector.h"
diff --git a/src/kernel.h b/src/kernel.h
deleted file mode 100644
index aead6a95adc35028834d671448223a31a57fc2b6..0000000000000000000000000000000000000000
--- a/src/kernel.h
+++ /dev/null
@@ -1,617 +0,0 @@
-/*******************************************************************************
- * This file is part of SWIFT.
- * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
- * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
- *
- * This program is free software: you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published
- * by the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program. If not, see <http://www.gnu.org/licenses/>.
- *
- ******************************************************************************/
-#ifndef SWIFT_KERNEL_H
-#define SWIFT_KERNEL_H
-
-/* Includes. */
-#include "const.h"
-#include "inline.h"
-#include "vector.h"
-
-/**
- * @file kernel.h
- * @brief SPH kernel functions. Compute W(x,h) and the gradient of W(x,h),
- * as well as the blending function used for gravity.
- */
-
-/* Gravity kernel stuff
- * -----------------------------------------------------------------------------------------------
- */
-
-/* The gravity kernel is defined as a degree 6 polynomial in the distance
- r. The resulting value should be post-multiplied with r^-3, resulting
- in a polynomial with terms ranging from r^-3 to r^3, which are
- sufficient to model both the direct potential as well as the splines
- near the origin. */
-
-/* Coefficients for the gravity kernel. */
-#define kernel_grav_degree 6
-#define kernel_grav_ivals 2
-#define kernel_grav_scale (2 * const_iepsilon)
-static float kernel_grav_coeffs
- [(kernel_grav_degree + 1) * (kernel_grav_ivals + 1)] = {
- 32.0f * const_iepsilon6, -192.0f / 5.0f * const_iepsilon5,
- 0.0f, 32.0f / 3.0f * const_iepsilon3,
- 0.0f, 0.0f,
- 0.0f, -32.0f / 3.0f * const_iepsilon6,
- 192.0f / 5.0f * const_iepsilon5, -48.0f * const_iepsilon4,
- 64.0f / 3.0f * const_iepsilon3, 0.0f,
- 0.0f, -1.0f / 15.0f,
- 0.0f, 0.0f,
- 0.0f, 0.0f,
- 0.0f, 0.0f,
- 1.0f};
-
-/**
- * @brief Computes the gravity cubic spline for a given distance x.
- */
-
-__attribute__((always_inline)) INLINE static void kernel_grav_eval(float x,
- float *W) {
- int ind = fmin(x * kernel_grav_scale, kernel_grav_ivals);
- float *coeffs = &kernel_grav_coeffs[ind * (kernel_grav_degree + 1)];
- float w = coeffs[0] * x + coeffs[1];
- for (int k = 2; k <= kernel_grav_degree; k++) w = x * w + coeffs[k];
- *W = w;
-}
-
-#ifdef VECTORIZE
-
-/**
- * @brief Computes the gravity cubic spline for a given distance x (Vectorized
- * version).
- */
-
-__attribute__((always_inline))
- INLINE static void kernel_grav_eval_vec(vector *x, vector *w) {
-
- vector ind, c[kernel_grav_degree + 1];
- int j, k;
-
- /* Load x and get the interval id. */
- ind.m = vec_ftoi(vec_fmin(x->v * vec_set1(kernel_grav_scale),
- vec_set1((float)kernel_grav_ivals)));
-
- /* load the coefficients. */
- for (k = 0; k < VEC_SIZE; k++)
- for (j = 0; j < kernel_grav_degree + 1; j++)
- c[j].f[k] = kernel_grav_coeffs[ind.i[k] * (kernel_grav_degree + 1) + j];
-
- /* Init the iteration for Horner's scheme. */
- w->v = (c[0].v * x->v) + c[1].v;
-
- /* And we're off! */
- for (int k = 2; k <= kernel_grav_degree; k++) w->v = (x->v * w->v) + c[k].v;
-}
-
-#endif
-
-/* Blending function stuff
- * --------------------------------------------------------------------------------------------
- */
-
-/* Coefficients for the blending function. */
-#define blender_degree 3
-#define blender_ivals 3
-#define blender_scale 4.0f
-static float blender_coeffs[(blender_degree + 1) * (blender_ivals + 1)] = {
- 0.0f, 0.0f, 0.0f, 1.0f, -32.0f, 24.0f, -6.0f, 1.5f,
- -32.0f, 72.0f, -54.0f, 13.5f, 0.0f, 0.0f, 0.0f, 0.0f};
-
-/**
- * @brief Computes the cubic spline blender for a given distance x.
- */
-
-__attribute__((always_inline)) INLINE static void blender_eval(float x,
- float *W) {
- int ind = fmin(x * blender_scale, blender_ivals);
- float *coeffs = &blender_coeffs[ind * (blender_degree + 1)];
- float w = coeffs[0] * x + coeffs[1];
- for (int k = 2; k <= blender_degree; k++) w = x * w + coeffs[k];
- *W = w;
-}
-
-/**
- * @brief Computes the cubic spline blender and its derivative for a given
- * distance x.
- */
-
-__attribute__((always_inline)) INLINE static void blender_deval(float x,
- float *W,
- float *dW_dx) {
- int ind = fminf(x * blender_scale, blender_ivals);
- float *coeffs = &blender_coeffs[ind * (blender_degree + 1)];
- float w = coeffs[0] * x + coeffs[1];
- float dw_dx = coeffs[0];
- for (int k = 2; k <= blender_degree; k++) {
- dw_dx = dw_dx * x + w;
- w = x * w + coeffs[k];
- }
- *W = w;
- *dW_dx = dw_dx;
-}
-
-#ifdef VECTORIZE
-
-/**
- * @brief Computes the cubic spline blender and its derivative for a given
- * distance x (Vectorized version). Gives a sensible answer only if x<2.
- */
-
-__attribute__((always_inline)) INLINE static void blender_eval_vec(vector *x,
- vector *w) {
-
- vector ind, c[blender_degree + 1];
- int j, k;
-
- /* Load x and get the interval id. */
- ind.m = vec_ftoi(
- vec_fmin(x->v * vec_set1(blender_scale), vec_set1((float)blender_ivals)));
-
- /* load the coefficients. */
- for (k = 0; k < VEC_SIZE; k++)
- for (j = 0; j < blender_degree + 1; j++)
- c[j].f[k] = blender_coeffs[ind.i[k] * (blender_degree + 1) + j];
-
- /* Init the iteration for Horner's scheme. */
- w->v = (c[0].v * x->v) + c[1].v;
-
- /* And we're off! */
- for (int k = 2; k <= blender_degree; k++) w->v = (x->v * w->v) + c[k].v;
-}
-
-/**
- * @brief Computes the cubic spline blender and its derivative for a given
- * distance x (Vectorized version). Gives a sensible answer only if x<2.
- */
-
-__attribute__((always_inline))
- INLINE static void blender_deval_vec(vector *x, vector *w, vector *dw_dx) {
-
- vector ind, c[blender_degree + 1];
- int j, k;
-
- /* Load x and get the interval id. */
- ind.m = vec_ftoi(
- vec_fmin(x->v * vec_set1(blender_scale), vec_set1((float)blender_ivals)));
-
- /* load the coefficients. */
- for (k = 0; k < VEC_SIZE; k++)
- for (j = 0; j < blender_degree + 1; j++)
- c[j].f[k] = blender_coeffs[ind.i[k] * (blender_degree + 1) + j];
-
- /* Init the iteration for Horner's scheme. */
- w->v = (c[0].v * x->v) + c[1].v;
- dw_dx->v = c[0].v;
-
- /* And we're off! */
- for (int k = 2; k <= blender_degree; k++) {
- dw_dx->v = (dw_dx->v * x->v) + w->v;
- w->v = (x->v * w->v) + c[k].v;
- }
-}
-
-#endif
-
-/* --------------------------------------------------------------------------------------------------------------------
- */
-
-#if defined(CUBIC_SPLINE_KERNEL)
-
-/* --------------------------------------------------------------------------------------------------------------------
- */
-
-/* Coefficients for the kernel. */
-#define kernel_name "Cubic spline"
-#define kernel_degree 3
-#define kernel_ivals 2
-#define kernel_gamma 2.0f
-#define kernel_gamma2 4.0f
-#define kernel_gamma3 8.0f
-#define kernel_igamma 0.5f
-#define kernel_nwneigh \
- (4.0 / 3.0 * M_PI *const_eta_kernel *const_eta_kernel *const_eta_kernel * \
- 6.0858f)
-static float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
- __attribute__((aligned(16))) = {
- 3.0 / 4.0 * M_1_PI, -3.0 / 2.0 * M_1_PI, 0.0, M_1_PI,
- -0.25 * M_1_PI, 3.0 / 2.0 * M_1_PI, -3.0 * M_1_PI, M_2_PI,
- 0.0, 0.0, 0.0, 0.0};
-#define kernel_root (kernel_coeffs[kernel_degree])
-#define kernel_wroot (4.0 / 3.0 * M_PI *kernel_coeffs[kernel_degree])
-
-/**
- * @brief Computes the cubic spline kernel and its derivative for a given
- * distance x. Gives a sensible answer only if x<2.
- */
-
-__attribute__((always_inline)) INLINE static void kernel_deval(float x,
- float *W,
- float *dW_dx) {
- int ind = fminf(x, kernel_ivals);
- float *coeffs = &kernel_coeffs[ind * (kernel_degree + 1)];
- float w = coeffs[0] * x + coeffs[1];
- float dw_dx = coeffs[0];
- for (int k = 2; k <= kernel_degree; k++) {
- dw_dx = dw_dx * x + w;
- w = x * w + coeffs[k];
- }
- *W = w;
- *dW_dx = dw_dx;
-}
-
-#ifdef VECTORIZE
-
-/**
- * @brief Computes the cubic spline kernel and its derivative for a given
- * distance x (Vectorized version). Gives a sensible answer only if x<2.
- */
-
-__attribute__((always_inline))
- INLINE static void kernel_deval_vec(vector *x, vector *w, vector *dw_dx) {
-
- vector ind, c[kernel_degree + 1];
- int j, k;
-
- /* Load x and get the interval id. */
- ind.m = vec_ftoi(vec_fmin(x->v, vec_set1((float)kernel_ivals)));
-
- /* load the coefficients. */
- for (k = 0; k < VEC_SIZE; k++)
- for (j = 0; j < kernel_degree + 1; j++)
- c[j].f[k] = kernel_coeffs[ind.i[k] * (kernel_degree + 1) + j];
-
- /* Init the iteration for Horner's scheme. */
- w->v = (c[0].v * x->v) + c[1].v;
- dw_dx->v = c[0].v;
-
- /* And we're off! */
- for (int k = 2; k <= kernel_degree; k++) {
- dw_dx->v = (dw_dx->v * x->v) + w->v;
- w->v = (x->v * w->v) + c[k].v;
- }
-}
-
-#endif
-
-/**
- * @brief Computes the cubic spline kernel for a given distance x. Gives a
- * sensible answer only if x<2.
- */
-
-__attribute__((always_inline)) INLINE static void kernel_eval(float x,
- float *W) {
- int ind = fmin(x, kernel_ivals);
- float *coeffs = &kernel_coeffs[ind * (kernel_degree + 1)];
- float w = coeffs[0] * x + coeffs[1];
- for (int k = 2; k <= kernel_degree; k++) w = x * w + coeffs[k];
- *W = w;
-}
-
-/* --------------------------------------------------------------------------------------------------------------------
- */
-
-#elif defined(QUARTIC_SPLINE_KERNEL)
-
-/* --------------------------------------------------------------------------------------------------------------------
- */
-
-/* Coefficients for the kernel. */
-#define kernel_name "Quartic spline"
-#define kernel_degree 4
-#define kernel_ivals 3
-#define kernel_gamma 2.5f
-#define kernel_gamma2 6.25f
-#define kernel_gamma3 15.625f
-#define kernel_igamma 0.4f
-#define kernel_nwneigh \
- (4.0 / 3.0 * M_PI *const_eta_kernel *const_eta_kernel *const_eta_kernel * \
- 8.2293f)
-static float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
- __attribute__((aligned(16))) = {
- 3.0 / 10.0 * M_1_PI, 0.0, -3.0 / 4.0 * M_1_PI,
- 0.0, 23.0 / 32.0 * M_1_PI, -1.0 / 5.0 * M_1_PI,
- M_1_PI, -3.0 / 2.0 * M_1_PI, 0.25 * M_1_PI,
- 11.0 / 16.0 * M_1_PI, 1.0 / 20.0 * M_1_PI, -0.5 * M_1_PI,
- 15.0 / 8.0 * M_1_PI, -25.0 / 8.0 * M_1_PI, 125.0 / 64.0 * M_1_PI,
- 0.0, 0.0, 0.0,
- 0.0, 0.0};
-#define kernel_root (kernel_coeffs[kernel_degree])
-#define kernel_wroot (4.0 / 3.0 * M_PI *kernel_coeffs[kernel_degree])
-
-/**
- * @brief Computes the quartic spline kernel and its derivative for a given
- * distance x. Gives a sensible answer only if x<2.5
- */
-
-__attribute__((always_inline)) INLINE static void kernel_deval(float x,
- float *W,
- float *dW_dx) {
- int ind = fminf(x + 0.5, kernel_ivals);
- float *coeffs = &kernel_coeffs[ind * (kernel_degree + 1)];
- float w = coeffs[0] * x + coeffs[1];
- float dw_dx = coeffs[0];
- for (int k = 2; k <= kernel_degree; k++) {
- dw_dx = dw_dx * x + w;
- w = x * w + coeffs[k];
- }
- *W = w;
- *dW_dx = dw_dx;
-}
-
-#ifdef VECTORIZE
-
-/**
- * @brief Computes the quartic spline kernel and its derivative for a given
- * distance x (Vectorized version). Gives a sensible answer only if x<2.5
- */
-
-__attribute__((always_inline))
- INLINE static void kernel_deval_vec(vector *x, vector *w, vector *dw_dx) {
-
- vector ind, c[kernel_degree + 1];
- int j, k;
-
- /* Load x and get the interval id. */
- ind.m = vec_ftoi(vec_fmin(x->v + 0.5f, vec_set1((float)kernel_ivals)));
-
- /* load the coefficients. */
- for (k = 0; k < VEC_SIZE; k++)
- for (j = 0; j < kernel_degree + 1; j++)
- c[j].f[k] = kernel_coeffs[ind.i[k] * (kernel_degree + 1) + j];
-
- /* Init the iteration for Horner's scheme. */
- w->v = (c[0].v * x->v) + c[1].v;
- dw_dx->v = c[0].v;
-
- /* And we're off! */
- for (int k = 2; k <= kernel_degree; k++) {
- dw_dx->v = (dw_dx->v * x->v) + w->v;
- w->v = (x->v * w->v) + c[k].v;
- }
-}
-
-#endif
-
-/**
- * @brief Computes the quartic spline kernel for a given distance x. Gives a
- * sensible answer only if x<2.5
- */
-
-__attribute__((always_inline)) INLINE static void kernel_eval(float x,
- float *W) {
- int ind = fmin(x + 0.5f, kernel_ivals);
- float *coeffs = &kernel_coeffs[ind * (kernel_degree + 1)];
- float w = coeffs[0] * x + coeffs[1];
- for (int k = 2; k <= kernel_degree; k++) w = x * w + coeffs[k];
- *W = w;
-}
-
-/* --------------------------------------------------------------------------------------------------------------------
- */
-
-#elif defined(QUINTIC_SPLINE_KERNEL)
-
-/* --------------------------------------------------------------------------------------------------------------------
- */
-
-/* Coefficients for the kernel. */
-#define kernel_name "Quintic spline"
-#define kernel_degree 5
-#define kernel_ivals 3
-#define kernel_gamma 3.f
-#define kernel_gamma2 9.f
-#define kernel_gamma3 27.f
-#define kernel_igamma 1.0f / 3.0f
-#define kernel_nwneigh \
- (4.0 / 3.0 * M_PI *const_eta_kernel *const_eta_kernel *const_eta_kernel * \
- 10.5868f)
-static float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
- __attribute__((aligned(16))) = {
- -1.0 / 12.0 * M_1_PI, 1.0 / 4.0 * M_1_PI, 0.0,
- -1.0 / 2.0 * M_1_PI, 0.0, 11.0 / 20.0 * M_1_PI,
- 1.0 / 24.0 * M_1_PI, -3.0 / 8.0 * M_1_PI, 5.0 / 4.0 * M_1_PI,
- -7.0 / 4.0 * M_1_PI, 5.0 / 8.0 * M_1_PI, 17.0 / 40.0 * M_1_PI,
- -1.0 / 120.0 * M_1_PI, 1.0 / 8.0 * M_1_PI, -3.0 / 4.0 * M_1_PI,
- 9.0 / 4.0 * M_1_PI, -27.0 / 8.0 * M_1_PI, 81.0 / 40.0 * M_1_PI,
- 0.0, 0.0, 0.0,
- 0.0, 0.0, 0.0};
-#define kernel_root (kernel_coeffs[kernel_degree])
-#define kernel_wroot (4.0 / 3.0 * M_PI *kernel_coeffs[kernel_degree])
-
-/**
- * @brief Computes the quintic spline kernel and its derivative for a given
- * distance x. Gives a sensible answer only if x<3.
- */
-
-__attribute__((always_inline)) INLINE static void kernel_deval(float x,
- float *W,
- float *dW_dx) {
- int ind = fminf(x, kernel_ivals);
- float *coeffs = &kernel_coeffs[ind * (kernel_degree + 1)];
- float w = coeffs[0] * x + coeffs[1];
- float dw_dx = coeffs[0];
- for (int k = 2; k <= kernel_degree; k++) {
- dw_dx = dw_dx * x + w;
- w = x * w + coeffs[k];
- }
- *W = w;
- *dW_dx = dw_dx;
-}
-
-#ifdef VECTORIZE
-
-/**
- * @brief Computes the quintic spline kernel and its derivative for a given
- * distance x (Vectorized version). Gives a sensible answer only if x<3.
- */
-
-__attribute__((always_inline))
- INLINE static void kernel_deval_vec(vector *x, vector *w, vector *dw_dx) {
-
- vector ind, c[kernel_degree + 1];
- int j, k;
-
- /* Load x and get the interval id. */
- ind.m = vec_ftoi(vec_fmin(x->v, vec_set1((float)kernel_ivals)));
-
- /* load the coefficients. */
- for (k = 0; k < VEC_SIZE; k++)
- for (j = 0; j < kernel_degree + 1; j++)
- c[j].f[k] = kernel_coeffs[ind.i[k] * (kernel_degree + 1) + j];
-
- /* Init the iteration for Horner's scheme. */
- w->v = (c[0].v * x->v) + c[1].v;
- dw_dx->v = c[0].v;
-
- /* And we're off! */
- for (int k = 2; k <= kernel_degree; k++) {
- dw_dx->v = (dw_dx->v * x->v) + w->v;
- w->v = (x->v * w->v) + c[k].v;
- }
-}
-
-#endif
-
-/**
- * @brief Computes the quintic spline kernel for a given distance x. Gives a
- * sensible answer only if x<3.
- */
-
-__attribute__((always_inline)) INLINE static void kernel_eval(float x,
- float *W) {
- int ind = fmin(x, kernel_ivals);
- float *coeffs = &kernel_coeffs[ind * (kernel_degree + 1)];
- float w = coeffs[0] * x + coeffs[1];
- for (int k = 2; k <= kernel_degree; k++) w = x * w + coeffs[k];
- *W = w;
-}
-
-/* --------------------------------------------------------------------------------------------------------------------
- */
-
-#elif defined(WENDLAND_C2_KERNEL)
-
-/* --------------------------------------------------------------------------------------------------------------------
- */
-
-/* Coefficients for the kernel. */
-#define kernel_name "Wendland C2"
-#define kernel_degree 5
-#define kernel_ivals 1
-#define kernel_gamma 2.f
-#define kernel_gamma2 4.f
-#define kernel_gamma3 8.f
-#define kernel_igamma 0.5f
-#define kernel_nwneigh \
- (4.0 / 3.0 * M_PI *const_eta_kernel *const_eta_kernel *const_eta_kernel * \
- 7.261825f)
-static float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
- __attribute__((aligned(16))) = {
- 0.05222272f, -0.39167037f, 1.04445431f, -1.04445431f, 0.f, 0.41778173f,
- 0.0f, 0.0f, 0.0f, 0.0f, 0.0f, 0.0f};
-#define kernel_root (kernel_coeffs[kernel_degree])
-#define kernel_wroot (4.0 / 3.0 * M_PI *kernel_coeffs[kernel_degree])
-
-/**
- * @brief Computes the quintic spline kernel and its derivative for a given
- * distance x. Gives a sensible answer only if x<1.
- */
-
-__attribute__((always_inline)) INLINE static void kernel_deval(float x,
- float *W,
- float *dW_dx) {
- int ind = fminf(0.5f * x, kernel_ivals);
- float *coeffs = &kernel_coeffs[ind * (kernel_degree + 1)];
- float w = coeffs[0] * x + coeffs[1];
- float dw_dx = coeffs[0];
- for (int k = 2; k <= kernel_degree; k++) {
- dw_dx = dw_dx * x + w;
- w = x * w + coeffs[k];
- }
- *W = w;
- *dW_dx = dw_dx;
-}
-
-#ifdef VECTORIZE
-
-/**
- * @brief Computes the Wendland C2 kernel and its derivative for a given
- * distance x (Vectorized version). Gives a sensible answer only if x<1.
- */
-
-__attribute__((always_inline))
- INLINE static void kernel_deval_vec(vector *x, vector *w, vector *dw_dx) {
-
- vector ind, c[kernel_degree + 1];
- int j, k;
-
- /* Load x and get the interval id. */
- ind.m = vec_ftoi(vec_fmin(0.5f * x->v, vec_set1((float)kernel_ivals)));
-
- /* load the coefficients. */
- for (k = 0; k < VEC_SIZE; k++)
- for (j = 0; j < kernel_degree + 1; j++)
- c[j].f[k] = kernel_coeffs[ind.i[k] * (kernel_degree + 1) + j];
-
- /* Init the iteration for Horner's scheme. */
- w->v = (c[0].v * x->v) + c[1].v;
- dw_dx->v = c[0].v;
-
- /* And we're off! */
- for (int k = 2; k <= kernel_degree; k++) {
- dw_dx->v = (dw_dx->v * x->v) + w->v;
- w->v = (x->v * w->v) + c[k].v;
- }
-}
-
-#endif
-
-/**
- * @brief Computes the Wendland C2 kernel for a given distance x. Gives a
- * sensible answer only if x<1.
- */
-
-__attribute__((always_inline)) INLINE static void kernel_eval(float x,
- float *W) {
- int ind = fmin(0.5f * x, kernel_ivals);
- float *coeffs = &kernel_coeffs[ind * (kernel_degree + 1)];
- float w = coeffs[0] * x + coeffs[1];
- for (int k = 2; k <= kernel_degree; k++) w = x * w + coeffs[k];
- *W = w;
-}
-
-/* --------------------------------------------------------------------------------------------------------------------
- */
-
-#else
-
-/* --------------------------------------------------------------------------------------------------------------------
- */
-
-#error "A kernel function must be chosen in const.h !!"
-
-#endif // Kernel choice
-
-/* Some cross-check functions */
-void SPH_kernel_dump(int N);
-void gravity_kernel_dump(float r_max, int N);
-
-#endif // SWIFT_KERNEL_H
diff --git a/src/kernel.c b/src/kernel_gravity.c
similarity index 78%
rename from src/kernel.c
rename to src/kernel_gravity.c
index 58f5b0c9fdaa62663c65d5af18afe0a15a813834..639a964c813ef7fd95008857ee17b7dd5ffafb27 100644
--- a/src/kernel.c
+++ b/src/kernel_gravity.c
@@ -21,32 +21,7 @@
#include <math.h>
#include <stdio.h>
-#include "kernel.h"
-
-/**
- * @brief Test the SPH kernel function by dumping it in the interval [0,1].
- *
- * @param N number of intervals in [0,1].
- */
-void SPH_kernel_dump(int N) {
-
- int k;
- float x, w, dw_dx;
- float x4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
- float w4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
- // float dw_dx4[4] __attribute__ ((aligned (16)));
-
- for (k = 0; k <= N; k++) {
- x = ((float)k) / N;
- x4[3] = x4[2];
- x4[2] = x4[1];
- x4[1] = x4[0];
- x4[0] = x;
- kernel_deval(x, &w, &dw_dx);
- // kernel_deval_vec( (vector *)x4 , (vector *)w4 , (vector *)dw_dx4 );
- printf(" %e %e %e %e %e %e %e\n", x, w, dw_dx, w4[0], w4[1], w4[2], w4[3]);
- }
-}
+#include "kernel_gravity.h"
/**
* @brief The Gadget-2 gravity kernel function
diff --git a/src/kernel_gravity.h b/src/kernel_gravity.h
new file mode 100644
index 0000000000000000000000000000000000000000..7fd4b061a7e94be01a11b06ad23d9113f579ebb8
--- /dev/null
+++ b/src/kernel_gravity.h
@@ -0,0 +1,209 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_KERNEL_GRAVITY_H
+#define SWIFT_KERNEL_GRAVITY_H
+
+/* Includes. */
+#include "const.h"
+#include "inline.h"
+#include "vector.h"
+
+/* Gravity kernel stuff
+ * -----------------------------------------------------------------------------------------------
+ */
+
+/* The gravity kernel is defined as a degree 6 polynomial in the distance
+ r. The resulting value should be post-multiplied with r^-3, resulting
+ in a polynomial with terms ranging from r^-3 to r^3, which are
+ sufficient to model both the direct potential as well as the splines
+ near the origin. */
+
+/* Coefficients for the gravity kernel. */
+#define kernel_grav_degree 6
+#define kernel_grav_ivals 2
+#define kernel_grav_scale (2 * const_iepsilon)
+static float kernel_grav_coeffs
+ [(kernel_grav_degree + 1) * (kernel_grav_ivals + 1)] = {
+ 32.0f * const_iepsilon6, -192.0f / 5.0f * const_iepsilon5,
+ 0.0f, 32.0f / 3.0f * const_iepsilon3,
+ 0.0f, 0.0f,
+ 0.0f, -32.0f / 3.0f * const_iepsilon6,
+ 192.0f / 5.0f * const_iepsilon5, -48.0f * const_iepsilon4,
+ 64.0f / 3.0f * const_iepsilon3, 0.0f,
+ 0.0f, -1.0f / 15.0f,
+ 0.0f, 0.0f,
+ 0.0f, 0.0f,
+ 0.0f, 0.0f,
+ 1.0f};
+
+/**
+ * @brief Computes the gravity cubic spline for a given distance x.
+ */
+
+__attribute__((always_inline)) INLINE static void kernel_grav_eval(float x,
+ float *W) {
+ int ind = fmin(x * kernel_grav_scale, kernel_grav_ivals);
+ float *coeffs = &kernel_grav_coeffs[ind * (kernel_grav_degree + 1)];
+ float w = coeffs[0] * x + coeffs[1];
+ for (int k = 2; k <= kernel_grav_degree; k++) w = x * w + coeffs[k];
+ *W = w;
+}
+
+#ifdef VECTORIZE
+
+/**
+ * @brief Computes the gravity cubic spline for a given distance x (Vectorized
+ * version).
+ */
+
+__attribute__((always_inline))
+ INLINE static void kernel_grav_eval_vec(vector *x, vector *w) {
+
+ vector ind, c[kernel_grav_degree + 1];
+ int j, k;
+
+ /* Load x and get the interval id. */
+ ind.m = vec_ftoi(vec_fmin(x->v * vec_set1(kernel_grav_scale),
+ vec_set1((float)kernel_grav_ivals)));
+
+ /* load the coefficients. */
+ for (k = 0; k < VEC_SIZE; k++)
+ for (j = 0; j < kernel_grav_degree + 1; j++)
+ c[j].f[k] = kernel_grav_coeffs[ind.i[k] * (kernel_grav_degree + 1) + j];
+
+ /* Init the iteration for Horner's scheme. */
+ w->v = (c[0].v * x->v) + c[1].v;
+
+ /* And we're off! */
+ for (int k = 2; k <= kernel_grav_degree; k++) w->v = (x->v * w->v) + c[k].v;
+}
+
+#endif
+
+/* Blending function stuff
+ * --------------------------------------------------------------------------------------------
+ */
+
+/* Coefficients for the blending function. */
+#define blender_degree 3
+#define blender_ivals 3
+#define blender_scale 4.0f
+static float blender_coeffs[(blender_degree + 1) * (blender_ivals + 1)] = {
+ 0.0f, 0.0f, 0.0f, 1.0f, -32.0f, 24.0f, -6.0f, 1.5f,
+ -32.0f, 72.0f, -54.0f, 13.5f, 0.0f, 0.0f, 0.0f, 0.0f};
+
+/**
+ * @brief Computes the cubic spline blender for a given distance x.
+ */
+
+__attribute__((always_inline)) INLINE static void blender_eval(float x,
+ float *W) {
+ int ind = fmin(x * blender_scale, blender_ivals);
+ float *coeffs = &blender_coeffs[ind * (blender_degree + 1)];
+ float w = coeffs[0] * x + coeffs[1];
+ for (int k = 2; k <= blender_degree; k++) w = x * w + coeffs[k];
+ *W = w;
+}
+
+/**
+ * @brief Computes the cubic spline blender and its derivative for a given
+ * distance x.
+ */
+
+__attribute__((always_inline)) INLINE static void blender_deval(float x,
+ float *W,
+ float *dW_dx) {
+ int ind = fminf(x * blender_scale, blender_ivals);
+ float *coeffs = &blender_coeffs[ind * (blender_degree + 1)];
+ float w = coeffs[0] * x + coeffs[1];
+ float dw_dx = coeffs[0];
+ for (int k = 2; k <= blender_degree; k++) {
+ dw_dx = dw_dx * x + w;
+ w = x * w + coeffs[k];
+ }
+ *W = w;
+ *dW_dx = dw_dx;
+}
+
+#ifdef VECTORIZE
+
+/**
+ * @brief Computes the cubic spline blender and its derivative for a given
+ * distance x (Vectorized version). Gives a sensible answer only if x<2.
+ */
+
+__attribute__((always_inline)) INLINE static void blender_eval_vec(vector *x,
+ vector *w) {
+
+ vector ind, c[blender_degree + 1];
+ int j, k;
+
+ /* Load x and get the interval id. */
+ ind.m = vec_ftoi(
+ vec_fmin(x->v * vec_set1(blender_scale), vec_set1((float)blender_ivals)));
+
+ /* load the coefficients. */
+ for (k = 0; k < VEC_SIZE; k++)
+ for (j = 0; j < blender_degree + 1; j++)
+ c[j].f[k] = blender_coeffs[ind.i[k] * (blender_degree + 1) + j];
+
+ /* Init the iteration for Horner's scheme. */
+ w->v = (c[0].v * x->v) + c[1].v;
+
+ /* And we're off! */
+ for (int k = 2; k <= blender_degree; k++) w->v = (x->v * w->v) + c[k].v;
+}
+
+/**
+ * @brief Computes the cubic spline blender and its derivative for a given
+ * distance x (Vectorized version). Gives a sensible answer only if x<2.
+ */
+
+__attribute__((always_inline))
+ INLINE static void blender_deval_vec(vector *x, vector *w, vector *dw_dx) {
+
+ vector ind, c[blender_degree + 1];
+ int j, k;
+
+ /* Load x and get the interval id. */
+ ind.m = vec_ftoi(
+ vec_fmin(x->v * vec_set1(blender_scale), vec_set1((float)blender_ivals)));
+
+ /* load the coefficients. */
+ for (k = 0; k < VEC_SIZE; k++)
+ for (j = 0; j < blender_degree + 1; j++)
+ c[j].f[k] = blender_coeffs[ind.i[k] * (blender_degree + 1) + j];
+
+ /* Init the iteration for Horner's scheme. */
+ w->v = (c[0].v * x->v) + c[1].v;
+ dw_dx->v = c[0].v;
+
+ /* And we're off! */
+ for (int k = 2; k <= blender_degree; k++) {
+ dw_dx->v = (dw_dx->v * x->v) + w->v;
+ w->v = (x->v * w->v) + c[k].v;
+ }
+}
+
+#endif
+
+void gravity_kernel_dump(float r_max, int N);
+
+#endif // SWIFT_KERNEL_GRAVITY_H
diff --git a/src/kernel_hydro.c b/src/kernel_hydro.c
new file mode 100644
index 0000000000000000000000000000000000000000..18a930d8ff7f792b2f9606787a6e4c547770629a
--- /dev/null
+++ b/src/kernel_hydro.c
@@ -0,0 +1,49 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2015 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#include <math.h>
+#include <stdio.h>
+
+#include "kernel_hydro.h"
+
+/**
+ * @brief Test the SPH kernel function by dumping it in the interval [0,1].
+ *
+ * @param N number of intervals in [0,1].
+ */
+void hydro_kernel_dump(int N) {
+
+ int k;
+ float x, w, dw_dx;
+ float x4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
+ float w4[4] = {0.0f, 0.0f, 0.0f, 0.0f};
+ // float dw_dx4[4] __attribute__ ((aligned (16)));
+
+ for (k = 0; k <= N; k++) {
+ x = ((float)k) / N;
+ x4[3] = x4[2];
+ x4[2] = x4[1];
+ x4[1] = x4[0];
+ x4[0] = x;
+ kernel_deval(x, &w, &dw_dx);
+ // kernel_deval_vec( (vector *)x4 , (vector *)w4 , (vector *)dw_dx4 );
+ printf(" %e %e %e %e %e %e %e\n", x, w, dw_dx, w4[0], w4[1], w4[2], w4[3]);
+ }
+}
diff --git a/src/kernel_hydro.h b/src/kernel_hydro.h
new file mode 100644
index 0000000000000000000000000000000000000000..66f51391fb9504ba30363b1980aaad1fcc9174b7
--- /dev/null
+++ b/src/kernel_hydro.h
@@ -0,0 +1,218 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_KERNEL_HYDRO_H
+#define SWIFT_KERNEL_HYDRO_H
+
+/* Includes. */
+#include "const.h"
+#include "error.h"
+#include "inline.h"
+#include "vector.h"
+
+/* ------------------------------------------------------------------------- */
+#if defined(CUBIC_SPLINE_KERNEL)
+
+/* Coefficients for the kernel. */
+#define kernel_name "Cubic spline (M4)"
+#define kernel_degree 3 /* Degree of the polynomial */
+#define kernel_ivals 2 /* Number of branches */
+#define kernel_gamma 1.825742
+#define kernel_constant 16. * M_1_PI
+static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
+ __attribute__((aligned(16))) = {3.f, -3.f, 0.f, 0.5f, /* 0 < u < 0.5 */
+ -1.f, 3.f, -3.f, 1.f, /* 0.5 < u < 1 */
+ 0.f, 0.f, 0.f, 0.f}; /* 1 < u */
+
+/* ------------------------------------------------------------------------- */
+#elif defined(QUARTIC_SPLINE_KERNEL)
+
+/* Coefficients for the kernel. */
+#define kernel_name "Quartic spline (M5)"
+#define kernel_degree 4
+#define kernel_ivals 5
+#define kernel_gamma 2.018932
+#define kernel_constant 15625. * M_1_PI / 512.
+static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
+ __attribute__((aligned(16))) = {
+ 6.f, 0.f, -2.4f, 0.f, 0.368f, /* 0 < u < 0.2 */
+ -4.f, 8.f, -4.8f, 0.32f, 0.352f, /* 0.2 < u < 0.4 */
+ -4.f, 8.f, -4.8f, 0.32f, 0.352f, /* 0.4 < u < 0.6 */
+ 1.f, -4.f, 6.f, -4.f, 1.f, /* 0.6 < u < 0.8 */
+ 1.f, -4.f, 6.f, -4.f, 1.f, /* 0.8 < u < 1 */
+ 0.f, 0.f, 0.f, 0.f, 0.f}; /* 1 < u */
+
+/* ------------------------------------------------------------------------- */
+#elif defined(QUINTIC_SPLINE_KERNEL)
+
+/* Coefficients for the kernel. */
+#define kernel_name "Quintic spline (M6)"
+#define kernel_degree 5
+#define kernel_ivals 3
+#define kernel_gamma 2.195775
+#define kernel_constant 2187. * M_1_PI / 40.
+static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
+ __attribute__((aligned(16))) = {
+ -10.f, 10.f, 0.f,
+ -2.2222222f, 0.f, 0.271604938f, /* 0 < u < 1/3 */
+ 5.f, -15.f, 16.666667f,
+ -7.77777777f, 0.925925f, 0.209876543f, /* 1/3 < u < 2/3 */
+ -1.f, 5.f, -10.f,
+ 10.f, -5.f, 1.f, /* 2/3 < u < 1. */
+ 0.f, 0.f, 0.f,
+ 0.f, 0.f, 0.f}; /* 1 < u */
+
+/* ------------------------------------------------------------------------- */
+#elif defined(WENDLAND_C2_KERNEL)
+
+/* Coefficients for the kernel. */
+#define kernel_name "Wendland C2"
+#define kernel_degree 5
+#define kernel_ivals 1
+#define kernel_gamma 1.936492
+#define kernel_constant 21. * M_1_PI / 2.
+static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
+ __attribute__((aligned(16))) = {
+ 4.f, -15.f, 20.f, -10.f, 0.f, 1.f, /* 0 < u < 1 */
+ 0.f, 0.f, 0.f, 0.f, 0.f, 0.f}; /* 1 < u */
+
+/* ------------------------------------------------------------------------- */
+#elif defined(WENDLAND_C4_KERNEL)
+
+/* Coefficients for the kernel. */
+#define kernel_name "Wendland C4"
+#define kernel_degree 8
+#define kernel_ivals 1
+#define kernel_gamma 2.207940
+#define kernel_constant 495. * M_1_PI / 32.
+static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
+ __attribute__((aligned(16))) = {
+ 11.666667f, -64.f, 140.f, -149.333333f, 70.f,
+ 0.f, -9.3333333f, 0.f, 1.f, /* 0 < u < 1 */
+ 0.f, 0.f, 0.f, 0.f, 0.f,
+ 0.f, 0.f, 0.f, 0.f}; /* 1 < u */
+
+/* ------------------------------------------------------------------------- */
+#elif defined(WENDLAND_C6_KERNEL)
+
+/* Coefficients for the kernel. */
+#define kernel_name "Wendland C6"
+#define kernel_degree 11
+#define kernel_ivals 1
+#define kernel_gamma 2.449490
+#define kernel_constant 1365. * M_1_PI / 64.
+static const float kernel_coeffs[(kernel_degree + 1) * (kernel_ivals + 1)]
+ __attribute__((aligned(16))) = {
+ 32.f, -231.f, 704.f, -1155.f, 1056.f, -462.f,
+ 0.f, 66.f, 0.f, -11.f, 0.f, 1.f, /* 0 < u < 1 */
+ 0.f, 0.f, 0.f, 0.f, 0.f, 0.f,
+ 0.f, 0.f, 0.f, 0.f, 0.f, 0.f}; /* 1 < u */
+
+/* ------------------------------------------------------------------------- */
+#else
+
+#error "A kernel function must be chosen in const.h !!"
+
+/* ------------------------------------------------------------------------- */
+#endif
+
+/* Ok, now comes the real deal. */
+
+/* First some powers of gamma = H/h */
+#define kernel_gamma2 kernel_gamma *kernel_gamma
+#define kernel_gamma3 kernel_gamma2 *kernel_gamma
+#define kernel_gamma4 kernel_gamma3 *kernel_gamma
+#define kernel_igamma 1. / kernel_gamma
+#define kernel_igamma2 kernel_igamma *kernel_igamma
+#define kernel_igamma3 kernel_igamma2 *kernel_igamma
+#define kernel_igamma4 kernel_igamma3 *kernel_igamma
+
+/* Some powers of eta */
+#define kernel_eta3 const_eta_kernel *const_eta_kernel *const_eta_kernel
+
+/* The number of neighbours (i.e. N_ngb) */
+#define kernel_nwneigh 4.0 * M_PI *kernel_gamma3 *kernel_eta3 / 3.0
+
+/* Kernel self contribution (i.e. W(0,h)) */
+#define kernel_root \
+ (kernel_coeffs[kernel_degree]) * kernel_constant *kernel_igamma3
+
+/**
+ * @brief Computes the kernel function and its derivative.
+ *
+ * Return 0 if $u > \\gamma = H/h$
+ *
+ * @param u The ratio of the distance to the smoothing length $u = x/h$.
+ * @param W (return) The value of the kernel function $W(x,h)$.
+ * @param dW_dx (return) The norm of the gradient of $|\\nabla W(x,h)|$.
+ */
+__attribute__((always_inline)) INLINE static void kernel_deval(
+ float u, float *const W, float *const dW_dx) {
+
+ /* Go to the range [0,1[ from [0,H[ */
+ const float x = u * (float)kernel_igamma;
+
+ /* Pick the correct branch of the kernel */
+ const int ind = (int)fminf(x * (float)kernel_ivals, kernel_ivals);
+ const float *const coeffs = &kernel_coeffs[ind * (kernel_degree + 1)];
+
+ /* First two terms of the polynomial ... */
+ float w = coeffs[0] * x + coeffs[1];
+ float dw_dx = coeffs[0];
+
+ /* ... and the rest of them */
+ for (int k = 2; k <= kernel_degree; k++) {
+ dw_dx = dw_dx * x + w;
+ w = x * w + coeffs[k];
+ }
+
+ /* Return everything */
+ *W = w * (float)kernel_constant * (float)kernel_igamma3;
+ *dW_dx = dw_dx * (float)kernel_constant * (float)kernel_igamma4;
+}
+
+/**
+ * @brief Computes the kernel function.
+ *
+ * @param u The ratio of the distance to the smoothing length $u = x/h$.
+ * @param W (return) The value of the kernel function $W(x,h)$.
+ */
+__attribute__((always_inline)) INLINE static void kernel_eval(float u,
+ float *const W) {
+ /* Go to the range [0,1[ from [0,H[ */
+ const float x = u * (float)kernel_igamma;
+
+ /* Pick the correct branch of the kernel */
+ const int ind = (int)fminf(x * (float)kernel_ivals, kernel_ivals);
+ const float *const coeffs = &kernel_coeffs[ind * (kernel_degree + 1)];
+
+ /* First two terms of the polynomial ... */
+ float w = coeffs[0] * x + coeffs[1];
+
+ /* ... and the rest of them */
+ for (int k = 2; k <= kernel_degree; k++) w = x * w + coeffs[k];
+
+ /* Return everything */
+ *W = w * (float)kernel_constant * (float)kernel_igamma3;
+}
+
+/* Some cross-check functions */
+void hydro_kernel_dump(int N);
+
+#endif // SWIFT_KERNEL_HYDRO_H
diff --git a/src/multipole.h b/src/multipole.h
index 91ba6df965ce9d3b088d538411b7f0a8555ba0e4..85ba44d3ce95d958b721d435ccd26b72e30a79c1 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -25,7 +25,7 @@
/* Includes. */
#include "const.h"
#include "inline.h"
-#include "kernel.h"
+#include "kernel_gravity.h"
#include "part.h"
/* Some constants. */
@@ -127,7 +127,7 @@ __attribute__((always_inline)) INLINE static void multipole_iact_mp(
/* Compute the forces on both multipoles. */
#if multipole_order == 1
- for (k = 0; k < 3; k++) p->a[k] += dx[k] * acc;
+ for (k = 0; k < 3; k++) p->a_grav[k] += dx[k] * acc;
#else
#error( "Multipoles of order %i not yet implemented." , multipole_order )
#endif
diff --git a/src/parser.c b/src/parser.c
new file mode 100644
index 0000000000000000000000000000000000000000..06dc819842d54d952704e4e0c40ebec5b561f691
--- /dev/null
+++ b/src/parser.c
@@ -0,0 +1,265 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 James Willis (james.s.willis@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"
+
+/* Some standard headers. */
+/* Needs to be included so that strtok returns char * instead of a int *. */
+#include <string.h>
+#include <stdlib.h>
+
+/* This object's header. */
+#include "parser.h"
+
+/* Local headers. */
+#include "error.h"
+
+/* Private functions. */
+static int count_char(char *str, char val);
+static void parse_line(FILE *fp, struct swift_params *params);
+
+/**
+ * @brief Reads an input file and stores each parameter in a structure.
+ *
+ * @param file_name Name of file to be read
+ * @param params Structure to be populated from file
+ */
+
+void parser_read_file(const char *file_name, struct swift_params *params) {
+
+ FILE *fp;
+
+ params->count = 0;
+
+ /* Open file for reading */
+ fp = fopen(file_name, "r");
+
+ if (fp == NULL) {
+ error("Error opening parameter file: %s", file_name);
+ }
+
+ /* Read until the end of the file is reached.*/
+ while (!feof(fp)) {
+ parse_line(fp, params);
+ }
+
+ fclose(fp);
+}
+
+/**
+ * @brief Counts the number of times a specific character appears in a string.
+ *
+ * @param str String to be checked
+ * @param val Character to be counted
+ */
+
+static int count_char(char *str, char val) {
+
+ int count = 0;
+
+ /* Check if the line contains the character */
+ while (*str) {
+ if (*str++ == val) ++count;
+ }
+
+ return count;
+}
+
+/**
+ * @brief Parses a line from a file and stores any parameters in a structure.
+ *
+ * @param fp File pointer to file to be read
+ * @param params Structure to be populated from file
+ *
+ */
+
+static void parse_line(FILE *fp, struct swift_params *params) {
+
+ char line[PARSER_MAX_LINE_SIZE];
+ char trim_line[PARSER_MAX_LINE_SIZE];
+
+ /* Read a line of the file */
+ if (fgets(line, PARSER_MAX_LINE_SIZE, fp) != NULL) {
+
+ char *token;
+ /* Remove comments */
+ token = strtok(line, PARSER_COMMENT_CHAR);
+ strcpy(trim_line, token);
+
+ /* Check if the line contains a value */
+ if (strchr(trim_line, PARSER_VALUE_CHAR)) {
+ /* Check for more than one parameter on the same line. */
+ if (count_char(trim_line, PARSER_VALUE_CHAR) > 1) {
+ error("Found more than one parameter in '%s', only one allowed.", line);
+ } else {
+ /* Take first token as the parameter name. */
+ token = strtok(trim_line, PARSER_VALUE_STRING);
+ strcpy(params->data[params->count].name, token);
+
+ /* Take second token as the parameter value. */
+ token = strtok(NULL, " #\n");
+ strcpy(params->data[params->count++].value, token);
+ }
+ }
+ }
+}
+
+/**
+ * @brief Retrieve integer parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @param retParam Value of the parameter found
+ *
+ */
+
+void parser_get_param_int(struct swift_params *params, char *name,
+ int *retParam) {
+
+ char str[128];
+
+ for (int i = 0; i < params->count; i++) {
+
+ /*strcmp returns 0 if both strings are the same.*/
+ if (!strcmp(name, params->data[i].name)) {
+
+ /* Check that exactly one number is parsed. */
+ if (sscanf(params->data[i].value, "%d%s", retParam, str) != 1) {
+ error(
+ "Tried parsing int '%s' but found '%s' with illegal integer "
+ "characters '%s'.",
+ params->data[i].name, params->data[i].value, str);
+ }
+
+ return;
+ }
+ }
+
+ message("Cannot find '%s' in the structure.", name);
+}
+
+/**
+ * @brief Retrieve float parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @param retParam Value of the parameter found
+ *
+ */
+
+void parser_get_param_float(struct swift_params *params, char *name,
+ float *retParam) {
+
+ char str[128];
+
+ for (int i = 0; i < params->count; i++) {
+
+ /*strcmp returns 0 if both strings are the same.*/
+ if (!strcmp(name, params->data[i].name)) {
+
+ /* Check that exactly one number is parsed. */
+ if (sscanf(params->data[i].value, "%f%s", retParam, str) != 1) {
+ error(
+ "Tried parsing float '%s' but found '%s' with illegal float "
+ "characters '%s'.",
+ params->data[i].name, params->data[i].value, str);
+ }
+
+ return;
+ }
+ }
+
+ message("Cannot find '%s' in the structure.", name);
+}
+
+/**
+ * @brief Retrieve double parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @param retParam Value of the parameter found
+ *
+ */
+
+void parser_get_param_double(struct swift_params *params, char *name,
+ double *retParam) {
+
+ char str[128];
+
+ for (int i = 0; i < params->count; i++) {
+
+ /*strcmp returns 0 if both strings are the same.*/
+ if (!strcmp(name, params->data[i].name)) {
+
+ /* Check that exactly one number is parsed. */
+ if (sscanf(params->data[i].value, "%lf", retParam) != 1) {
+ error(
+ "Tried parsing double '%s' but found '%s' with illegal double "
+ "characters '%s'.",
+ params->data[i].name, params->data[i].value, str);
+ }
+
+ return;
+ }
+ }
+
+ message("Cannot find '%s' in the structure.", name);
+}
+
+/**
+ * @brief Retrieve string parameter from structure.
+ *
+ * @param params Structure that holds the parameters
+ * @param name Name of the parameter to be found
+ * @param retParam Value of the parameter found
+ *
+ */
+
+void parser_get_param_string(struct swift_params *params, char *name,
+ char *retParam) {
+
+ for (int i = 0; i < params->count; i++) {
+
+ /*strcmp returns 0 if both strings are the same.*/
+ if (!strcmp(name, params->data[i].name)) {
+ strcpy(retParam, params->data[i].value);
+ return;
+ }
+ }
+}
+
+/**
+ * @brief Prints the contents of the parameter structure.
+ *
+ * @param params Structure that holds the parameters
+ *
+ */
+
+void parser_print_params(struct swift_params *params) {
+
+ printf("\n--------------------------\n");
+ printf("| SWIFT Parameter File |\n");
+ printf("--------------------------\n");
+
+ for (int i = 0; i < params->count; i++) {
+ printf("Parameter name: %s\n", params->data[i].name);
+ printf("Parameter value: %s\n", params->data[i].value);
+ }
+}
diff --git a/src/parser.h b/src/parser.h
new file mode 100644
index 0000000000000000000000000000000000000000..2fb4148944cd423da016341744cb6d58e222182e
--- /dev/null
+++ b/src/parser.h
@@ -0,0 +1,54 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 James Willis (james.s.willis@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_PARSER_H
+#define SWIFT_PARSER_H
+
+#include <stdio.h>
+
+#define PARSER_MAX_LINE_SIZE 128
+#define PARSER_MAX_NO_OF_PARAMS 512
+
+#define PARSER_COMMENT_CHAR "#"
+#define PARSER_VALUE_CHAR ':'
+#define PARSER_VALUE_STRING ":"
+#define PARSER_END_OF_FILE "..."
+
+struct parameter {
+ char name[PARSER_MAX_LINE_SIZE];
+ char value[PARSER_MAX_LINE_SIZE];
+};
+
+struct swift_params {
+ struct parameter data[PARSER_MAX_NO_OF_PARAMS];
+ int count;
+};
+
+/* Public API. */
+void parser_read_file(const char *file_name, struct swift_params *params);
+void parser_print_params(struct swift_params *params);
+void parser_get_param_int(struct swift_params *params, char *name,
+ int *retParam);
+void parser_get_param_float(struct swift_params *params, char *name,
+ float *retParam);
+void parser_get_param_double(struct swift_params *params, char *name,
+ double *retParam);
+void parser_get_param_string(struct swift_params *params, char *name,
+ char *retParam);
+
+#endif /* SWIFT_PARSER_H */
diff --git a/src/partition.c b/src/partition.c
index a55a5c75367887b6a0bc1875415c354a7e484371..e6422471841772eb36ecabaa4355fcbeb6882f04 100644
--- a/src/partition.c
+++ b/src/partition.c
@@ -35,7 +35,7 @@
#include <stdio.h>
#include <stdlib.h>
#include <strings.h>
-#include <values.h>
+#include <float.h>
/* MPI headers. */
#ifdef WITH_MPI
@@ -654,6 +654,7 @@ static void repart_edge_metis(int partweights, int bothweights, int nodeID,
split_metis(s, nr_nodes, celllist);
/* Clean up. */
+ free(inds);
if (bothweights) free(weights_v);
free(weights_e);
free(celllist);
diff --git a/src/riemann/riemann_exact.h b/src/riemann/riemann_exact.h
index 861dad9729794efb302638792fef6e3df43c700a..b768cde5f4f5dfd0463cc8a582a1af0a17607bbe 100644
--- a/src/riemann/riemann_exact.h
+++ b/src/riemann/riemann_exact.h
@@ -192,6 +192,8 @@ __attribute__((always_inline)) INLINE static GFLOAT riemann_guess_p(
*
* @param lower_limit Lower limit for the method (riemann_f(lower_limit) < 0)
* @param upper_limit Upper limit for the method (riemann_f(upper_limit) > 0)
+ * @param lowf ??? Bert?
+ * @param upf ??? Bert?
* @param error_tol Tolerance used to decide if the solution is converged
* @param WL Left state vector
* @param WR Right state vector
diff --git a/src/runner.c b/src/runner.c
index 09a68fa9c71b4e838aafb724fcc30b1aa66b0601..474acf624f499c20157be5840ac103f7dc721d1f 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -469,8 +469,10 @@ void runner_dogsort(struct runner *r, struct cell *c, int flags, int clock) {
void runner_doinit(struct runner *r, struct cell *c, int timer) {
- struct part *p, *parts = c->parts;
+ struct part *const parts = c->parts;
+ struct gpart *const gparts = c->gparts;
const int count = c->count;
+ const int gcount = c->gcount;
const int ti_current = r->e->ti_current;
TIMER_TIC;
@@ -486,7 +488,7 @@ void runner_doinit(struct runner *r, struct cell *c, int timer) {
for (int i = 0; i < count; i++) {
/* Get a direct pointer on the part. */
- p = &parts[i];
+ struct part *const p = &parts[i];
if (p->ti_end <= ti_current) {
@@ -494,6 +496,19 @@ void runner_doinit(struct runner *r, struct cell *c, int timer) {
hydro_init_part(p);
}
}
+
+ /* Loop over the gparts in this cell. */
+ for (int i = 0; i < gcount; i++) {
+
+ /* Get a direct pointer on the part. */
+ struct gpart *const gp = &gparts[i];
+
+ if (gp->ti_end <= ti_current) {
+
+ /* Get ready for a density calculation */
+ gravity_init_part(gp);
+ }
+ }
}
if (timer) TIMER_TOC(timer_init);
@@ -649,7 +664,7 @@ void runner_doghost(struct runner *r, struct cell *c) {
}
/**
- * @brief Drift particles forward in time
+ * @brief Drift particles and g-particles forward in time
*
* @param r The runner thread.
* @param c The cell.
@@ -657,12 +672,14 @@ void runner_doghost(struct runner *r, struct cell *c) {
*/
void runner_dodrift(struct runner *r, struct cell *c, int timer) {
+ const int nr_gparts = c->gcount;
const double timeBase = r->e->timeBase;
const double dt = (r->e->ti_current - r->e->ti_old) * timeBase;
- const float ti_old = r->e->ti_old;
- const float ti_current = r->e->ti_current;
- struct part *restrict parts = c->parts;
- struct xpart *restrict xparts = c->xparts;
+ const int ti_old = r->e->ti_old;
+ const int ti_current = r->e->ti_current;
+ struct part *const parts = c->parts;
+ struct xpart *const xparts = c->xparts;
+ struct gpart *const gparts = c->gparts;
float dx_max = 0.f, dx2_max = 0.f, h_max = 0.f;
TIMER_TIC
@@ -670,13 +687,23 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
/* No children? */
if (!c->split) {
- /* Loop over all the particles in the cell */
+ /* Loop over all the g-particles in the cell */
const size_t nr_parts = c->count;
for (size_t k = 0; k < nr_parts; k++) {
+ /* Get a handle on the gpart. */
+ struct gpart *const gp = &gparts[k];
+
+ /* Drift... */
+ gp->x[0] += gp->v_full[0] * dt;
+ gp->x[1] += gp->v_full[1] * dt;
+ gp->x[2] += gp->v_full[2] * dt;
+ }
+
+ /* Loop over all the particles in the cell (more work for these !) */
/* Get a handle on the part. */
- struct part *p = &parts[k];
- struct xpart *xp = &xparts[k];
+ struct part *const p = &parts[k];
+ struct xpart *const xp = &xparts[k];
/* Useful quantity */
const float h_inv = 1.0f / p->h;
@@ -692,18 +719,18 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
p->v[2] += p->a_hydro[2] * dt;
/* Predict smoothing length */
- float w = p->h_dt * h_inv * dt;
- if (fabsf(w) < 0.2f)
- p->h *= approx_expf(w); /* 4th order expansion of exp(w) */
+ const float w1 = p->h_dt * h_inv * dt;
+ if (fabsf(w1) < 0.2f)
+ p->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
else
- p->h *= expf(w);
+ p->h *= expf(w1);
/* Predict density */
- w = -3.0f * p->h_dt * h_inv * dt;
- if (fabsf(w) < 0.2f)
- p->rho *= approx_expf(w); /* 4th order expansion of exp(w) */
+ const float w2 = -3.0f * p->h_dt * h_inv * dt;
+ if (fabsf(w2) < 0.2f)
+ p->rho *= approx_expf(w2); /* 4th order expansion of exp(w) */
else
- p->rho *= expf(w);
+ p->rho *= expf(w2);
/* Predict the values of the extra fields */
hydro_predict_extra(p, xp, ti_old, ti_current, timeBase);
@@ -759,37 +786,105 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
const double timeBase = r->e->timeBase;
const double timeBase_inv = 1.0 / r->e->timeBase;
const int count = c->count;
+ const int gcount = c->gcount;
+ struct part *const parts = c->parts;
+ struct xpart *const xparts = c->xparts;
+ struct gpart *const gparts = c->gparts;
const int is_fixdt =
(r->e->policy & engine_policy_fixdt) == engine_policy_fixdt;
- int new_dti;
- int dti_timeline;
-
- int updated = 0;
+ int updated = 0, g_updated = 0;
int ti_end_min = max_nr_timesteps, ti_end_max = 0;
double e_kin = 0.0, e_int = 0.0, e_pot = 0.0, mass = 0.0;
float mom[3] = {0.0f, 0.0f, 0.0f};
float ang[3] = {0.0f, 0.0f, 0.0f};
- float x[3], v_full[3];
- struct part *restrict p, *restrict parts = c->parts;
- struct xpart *restrict xp, *restrict xparts = c->xparts;
TIMER_TIC
/* No children? */
if (!c->split) {
+ /* Loop over the g-particles and kick the active ones. */
+ for (int k = 0; k < gcount; k++) {
+
+ /* Get a handle on the part. */
+ struct gpart *const gp = &gparts[k];
+
+ /* If the g-particle has no counterpart and needs to be kicked */
+ if (gp->id < 0 && (is_fixdt || gp->ti_end <= ti_current)) {
+
+ /* First, finish the force calculation */
+ gravity_end_force(gp);
+
+ /* Now we are ready to compute the next time-step size */
+ int new_dti;
+
+ if (is_fixdt) {
+
+ /* Now we have a time step, proceed with the kick */
+ new_dti = global_dt_max * timeBase_inv;
+
+ } else {
+
+ /* Compute the next timestep (gravity condition) */
+ float new_dt = gravity_compute_timestep(gp);
+
+ /* Limit timestep within the allowed range */
+ new_dt = fminf(new_dt, global_dt_max);
+ new_dt = fmaxf(new_dt, global_dt_min);
+
+ /* Convert to integer time */
+ new_dti = new_dt * timeBase_inv;
+
+ /* Recover the current timestep */
+ const int current_dti = gp->ti_end - gp->ti_begin;
+
+ /* Limit timestep increase */
+ if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
+
+ /* Put this timestep on the time line */
+ int dti_timeline = max_nr_timesteps;
+ while (new_dti < dti_timeline) dti_timeline /= 2;
+
+ /* Now we have a time step, proceed with the kick */
+ new_dti = dti_timeline;
+ }
+
+ /* Compute the time step for this kick */
+ const int ti_start = (gp->ti_begin + gp->ti_end) / 2;
+ const int ti_end = gp->ti_end + new_dti / 2;
+ const double dt = (ti_end - ti_start) * timeBase;
+ const double half_dt = (ti_end - gp->ti_end) * timeBase;
+
+ /* Move particle forward in time */
+ gp->ti_begin = gp->ti_end;
+ gp->ti_end = gp->ti_begin + new_dti;
+
+ /* Kick particles in momentum space */
+ gp->v_full[0] += gp->a_grav[0] * dt;
+ gp->v_full[1] += gp->a_grav[1] * dt;
+ gp->v_full[2] += gp->a_grav[2] * dt;
+
+ /* Extra kick work */
+ gravity_kick_extra(gp, dt, half_dt);
+
+ /* Number of updated g-particles */
+ g_updated++;
+ }
+
+ /* Minimal time for next end of time-step */
+ ti_end_min = min(gp->ti_end, ti_end_min);
+ ti_end_max = max(gp->ti_end, ti_end_max);
+ }
+
+ /* Now do the hydro ones... */
+
/* Loop over the particles and kick the active ones. */
for (int k = 0; k < count; k++) {
/* Get a handle on the part. */
- p = &parts[k];
- xp = &xparts[k];
-
- const float m = p->mass;
- x[0] = p->x[0];
- x[1] = p->x[1];
- x[2] = p->x[2];
+ struct part *const p = &parts[k];
+ struct xpart *const xp = &xparts[k];
/* If particle needs to be kicked */
if (is_fixdt || p->ti_end <= ti_current) {
@@ -799,8 +894,10 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
/* And do the same of the extra variable */
hydro_end_force(p);
+ if (p->gpart != NULL) gravity_end_force(p->gpart);
/* Now we are ready to compute the next time-step size */
+ int new_dti;
if (is_fixdt) {
@@ -809,9 +906,13 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
} else {
- /* Compute the next timestep */
+ /* Compute the next timestep (hydro condition) */
const float new_dt_hydro = hydro_compute_timestep(p, xp);
- const float new_dt_grav = gravity_compute_timestep(p, xp);
+
+ /* Compute the next timestep (gravity condition) */
+ float new_dt_grav = FLT_MAX;
+ if (p->gpart != NULL)
+ new_dt_grav = gravity_compute_timestep(p->gpart);
float new_dt = fminf(new_dt_hydro, new_dt_grav);
@@ -836,7 +937,7 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
/* Put this timestep on the time line */
- dti_timeline = max_nr_timesteps;
+ int dti_timeline = max_nr_timesteps;
while (new_dti < dti_timeline) dti_timeline /= 2;
/* Now we have a time step, proceed with the kick */
@@ -846,34 +947,55 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
/* Compute the time step for this kick */
const int ti_start = (p->ti_begin + p->ti_end) / 2;
const int ti_end = p->ti_end + new_dti / 2;
- const float dt = (ti_end - ti_start) * timeBase;
- const float half_dt = (ti_end - p->ti_end) * timeBase;
+ const double dt = (ti_end - ti_start) * timeBase;
+ const double half_dt = (ti_end - p->ti_end) * timeBase;
/* Move particle forward in time */
p->ti_begin = p->ti_end;
p->ti_end = p->ti_begin + new_dti;
+ if (p->gpart != NULL) {
+ p->gpart->ti_begin = p->ti_begin;
+ p->gpart->ti_end = p->ti_end;
+ }
+
+ /* Get the acceleration */
+ float a_tot[3] = {p->a_hydro[0], p->a_hydro[1], p->a_hydro[2]};
+ if (p->gpart != NULL) {
+ a_tot[0] += p->gpart->a_grav[0];
+ a_tot[1] += p->gpart->a_grav[1];
+ a_tot[1] += p->gpart->a_grav[2];
+ }
/* Kick particles in momentum space */
- xp->v_full[0] += p->a_hydro[0] * dt;
- xp->v_full[1] += p->a_hydro[1] * dt;
- xp->v_full[2] += p->a_hydro[2] * dt;
+ xp->v_full[0] += a_tot[0] * dt;
+ xp->v_full[1] += a_tot[1] * dt;
+ xp->v_full[2] += a_tot[2] * dt;
+
+ if (p->gpart != NULL) {
+ 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];
+ }
- p->v[0] = xp->v_full[0] - half_dt * p->a_hydro[0];
- p->v[1] = xp->v_full[1] - half_dt * p->a_hydro[1];
- p->v[2] = xp->v_full[2] - half_dt * p->a_hydro[2];
+ /* Go back by half-step for the hydro velocity */
+ p->v[0] = xp->v_full[0] - half_dt * a_tot[0];
+ p->v[1] = xp->v_full[1] - half_dt * a_tot[1];
+ p->v[2] = xp->v_full[2] - half_dt * a_tot[2];
/* Extra kick work */
hydro_kick_extra(p, xp, dt, half_dt);
+ if (p->gpart != NULL) gravity_kick_extra(p->gpart, dt, half_dt);
/* Number of updated particles */
updated++;
+ if (p->gpart != NULL) g_updated++;
}
/* Now collect quantities for statistics */
- v_full[0] = xp->v_full[0];
- v_full[1] = xp->v_full[1];
- v_full[2] = xp->v_full[2];
+ const double x[3] = {p->x[0], p->x[1], p->x[2]};
+ const float v_full[3] = {xp->v_full[0], xp->v_full[1], xp->v_full[2]};
+ const float m = p->mass;
/* Collect mass */
mass += m;
@@ -907,13 +1029,14 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
/* Loop over the progeny. */
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) {
- struct cell *cp = c->progeny[k];
+ struct cell *const cp = c->progeny[k];
/* Recurse */
runner_dokick(r, cp, 0);
/* And aggregate */
updated += cp->updated;
+ g_updated += cp->g_updated;
e_kin += cp->e_kin;
e_int += cp->e_int;
e_pot += cp->e_pot;
@@ -931,6 +1054,7 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
/* Store the values. */
c->updated = updated;
+ c->g_updated = g_updated;
c->e_kin = e_kin;
c->e_int = e_int;
c->e_pot = e_pot;
@@ -1057,9 +1181,12 @@ void *runner_main(void *data) {
case task_type_grav_down:
runner_dograv_down(r, t->ci);
break;
- case task_type_psort:
+ case task_type_part_sort:
space_do_parts_sort();
break;
+ case task_type_gpart_sort:
+ space_do_gparts_sort();
+ break;
case task_type_split_cell:
space_do_split(e->s, t->ci);
break;
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index f374339da75e31b39a5295fcd8bbc23c34d8d67d..02626295a49f314fef840bc044a476f5c9cf332d 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -267,9 +267,9 @@ void runner_dograv_down(struct runner *r, struct cell *c) {
/* Apply the multipole acceleration to all gparts. */
for (int k = 0; k < c->gcount; k++) {
struct gpart *p = &c->gparts[k];
- p->a[0] += m->a[0];
- p->a[1] += m->a[1];
- p->a[2] += m->a[2];
+ p->a_grav[0] += m->a[0];
+ p->a_grav[1] += m->a[1];
+ p->a_grav[2] += m->a[2];
}
}
}
@@ -594,5 +594,4 @@ void runner_dosub_grav(struct runner *r, struct cell *ci, struct cell *cj,
if (gettimer) TIMER_TOC(timer_dosub_grav);
}
-
#endif /* SWIFT_RUNNER_DOIACT_GRAV_H */
diff --git a/src/scheduler.c b/src/scheduler.c
index 38a1cd8c663307e0c0378d8bec2e0cd3d8f37fa8..d1d343240b37f5afd5f41fecacf106b0e85f726f 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -43,7 +43,7 @@
#include "cycle.h"
#include "error.h"
#include "intrinsics.h"
-#include "kernel.h"
+#include "kernel_hydro.h"
#include "timers.h"
/**
@@ -124,7 +124,9 @@ void scheduler_splittasks(struct scheduler *s) {
}
/* Skip sorting tasks. */
- if (t->type == task_type_psort) continue;
+ if (t->type == task_type_part_sort) continue;
+
+ if (t->type == task_type_gpart_sort) continue;
/* Empty task? */
if (t->ci == NULL || (t->type == task_type_pair && t->cj == NULL)) {
diff --git a/src/space.c b/src/space.c
index 2bd68233e8f6e3d4dc214091b272b1770e3c0201..1be4e3171ae1c51543e08439bb9f1b40a9d50b0e 100644
--- a/src/space.c
+++ b/src/space.c
@@ -39,7 +39,7 @@
#include "atomic.h"
#include "engine.h"
#include "error.h"
-#include "kernel.h"
+#include "kernel_hydro.h"
#include "lock.h"
#include "minmax.h"
#include "runner.h"
@@ -490,7 +490,7 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
#endif
/* Sort the parts according to their cells. */
- space_gparts_sort(s->gparts, gind, nr_gparts, 0, s->nr_cells - 1);
+ space_gparts_sort(s, gind, nr_gparts, 0, s->nr_cells - 1, verbose);
/* Re-link the parts. */
for (int k = 0; k < nr_gparts; k++)
@@ -499,6 +499,28 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
/* We no longer need the indices as of here. */
free(gind);
+ /* Verify that the links are correct */
+ /* MATTHIEU: To be commented out once we are happy */
+ for (size_t k = 0; k < nr_gparts; ++k) {
+
+ if (s->gparts[k].id > 0) {
+
+ if (s->gparts[k].part->gpart != &s->gparts[k]) error("Linking problem !");
+
+ if (s->gparts[k].x[0] != s->gparts[k].part->x[0] ||
+ s->gparts[k].x[1] != s->gparts[k].part->x[1] ||
+ s->gparts[k].x[2] != s->gparts[k].part->x[2])
+ error("Linked particles are not at the same position !");
+ }
+ }
+ for (size_t k = 0; k < nr_parts; ++k) {
+
+ if (s->parts[k].gpart != NULL) {
+
+ if (s->parts[k].gpart->part != &s->parts[k]) error("Linking problem !");
+ }
+ }
+
/* Hook the cells up to the parts. */
// tic = getticks();
struct part *finger = s->parts;
@@ -561,7 +583,7 @@ void space_split(struct space *s, struct cell *cells, int verbose) {
void space_parts_sort(struct space *s, int *ind, size_t N, int min, int max,
int verbose) {
- ticks tic = getticks();
+ const ticks tic = getticks();
/*Populate the global parallel_sort structure with the input data */
space_sort_struct.parts = s->parts;
@@ -585,7 +607,7 @@ void space_parts_sort(struct space *s, int *ind, size_t N, int min, int max,
space_sort_struct.waiting = 1;
/* Launch the sorting tasks. */
- engine_launch(s->e, s->e->nr_threads, (1 << task_type_psort), 0);
+ engine_launch(s->e, s->e->nr_threads, (1 << task_type_part_sort), 0);
/* Verify space_sort_struct. */
/* for (int i = 1; i < N; i++)
@@ -728,103 +750,140 @@ void space_do_parts_sort() {
} /* main loop. */
}
-void space_gparts_sort(struct gpart *gparts, int *ind, size_t N, int min,
- int max) {
-
- struct qstack {
- volatile size_t i, j;
- volatile int min, max;
- volatile int ready;
- };
- struct qstack *qstack;
- int qstack_size = 2 * (max - min) + 10;
- volatile unsigned int first, last, waiting;
-
- int pivot;
- ptrdiff_t i, ii, j, jj, temp_i;
- int qid;
- struct gpart temp_p;
-
- /* for ( int k = 0 ; k < N ; k++ )
- if ( ind[k] > max || ind[k] < min )
- error( "ind[%i]=%i is not in [%i,%i]." , k , ind[k] , min , max ); */
-
- /* Allocate the stack. */
- if ((qstack = malloc(sizeof(struct qstack) * qstack_size)) == NULL)
- error("Failed to allocate qstack.");
-
- /* Init the interval stack. */
- qstack[0].i = 0;
- qstack[0].j = N - 1;
- qstack[0].min = min;
- qstack[0].max = max;
- qstack[0].ready = 1;
- for (i = 1; i < qstack_size; i++) qstack[i].ready = 0;
- first = 0;
- last = 1;
- waiting = 1;
+/**
+ * @brief Sort the g-particles and condensed particles according to the given
+ *indices.
+ *
+ * @param s The #space.
+ * @param ind The indices with respect to which the gparts are sorted.
+ * @param N The number of gparts
+ * @param min Lowest index.
+ * @param max highest index.
+ * @param verbose Are we talkative ?
+ */
+void space_gparts_sort(struct space *s, int *ind, size_t N, int min, int max,
+ int verbose) {
+
+ const ticks tic = getticks();
+
+ /*Populate the global parallel_sort structure with the input data */
+ space_sort_struct.gparts = s->gparts;
+ space_sort_struct.ind = ind;
+ space_sort_struct.stack_size = 2 * (max - min + 1) + 10 + s->e->nr_threads;
+ if ((space_sort_struct.stack = malloc(sizeof(struct qstack) *
+ space_sort_struct.stack_size)) == NULL)
+ error("Failed to allocate sorting stack.");
+ for (int i = 0; i < space_sort_struct.stack_size; i++)
+ space_sort_struct.stack[i].ready = 0;
+
+ /* Add the first interval. */
+ space_sort_struct.stack[0].i = 0;
+ space_sort_struct.stack[0].j = N - 1;
+ space_sort_struct.stack[0].min = min;
+ space_sort_struct.stack[0].max = max;
+ space_sort_struct.stack[0].ready = 1;
+ space_sort_struct.first = 0;
+ space_sort_struct.last = 1;
+ space_sort_struct.waiting = 1;
+
+ /* Launch the sorting tasks. */
+ engine_launch(s->e, s->e->nr_threads, (1 << task_type_gpart_sort), 0);
+
+ /* Verify space_sort_struct. */
+ /* for (int i = 1; i < N; i++)
+ if (ind[i - 1] > ind[i])
+ error("Sorting failed (ind[%i]=%i,ind[%i]=%i), min=%i, max=%i.", i - 1,
+ ind[i - 1], i,
+ ind[i], min, max);
+ message("Sorting succeeded."); */
+
+ /* Clean up. */
+ free(space_sort_struct.stack);
+
+ if (verbose)
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
+}
+
+void space_do_gparts_sort() {
+
+ /* Pointers to the sorting data. */
+ int *ind = space_sort_struct.ind;
+ struct gpart *gparts = space_sort_struct.gparts;
/* Main loop. */
- while (waiting > 0) {
+ while (space_sort_struct.waiting) {
/* Grab an interval off the queue. */
- qid = (first++) % qstack_size;
+ int qid =
+ atomic_inc(&space_sort_struct.first) % space_sort_struct.stack_size;
+
+ /* Wait for the entry to be ready, or for the sorting do be done. */
+ while (!space_sort_struct.stack[qid].ready)
+ if (!space_sort_struct.waiting) return;
/* Get the stack entry. */
- i = qstack[qid].i;
- j = qstack[qid].j;
- min = qstack[qid].min;
- max = qstack[qid].max;
- qstack[qid].ready = 0;
+ ptrdiff_t i = space_sort_struct.stack[qid].i;
+ ptrdiff_t j = space_sort_struct.stack[qid].j;
+ int min = space_sort_struct.stack[qid].min;
+ int max = space_sort_struct.stack[qid].max;
+ space_sort_struct.stack[qid].ready = 0;
/* Loop over sub-intervals. */
while (1) {
/* Bring beer. */
- pivot = (min + max) / 2;
+ const int pivot = (min + max) / 2;
+ /* message("Working on interval [%i,%i] with min=%i, max=%i, pivot=%i.",
+ i, j, min, max, pivot); */
/* One pass of QuickSort's partitioning. */
- ii = i;
- jj = j;
+ ptrdiff_t ii = i;
+ ptrdiff_t jj = j;
while (ii < jj) {
while (ii <= j && ind[ii] <= pivot) ii++;
while (jj >= i && ind[jj] > pivot) jj--;
if (ii < jj) {
- temp_i = ind[ii];
+ size_t temp_i = ind[ii];
ind[ii] = ind[jj];
ind[jj] = temp_i;
- temp_p = gparts[ii];
+ struct gpart temp_p = gparts[ii];
gparts[ii] = gparts[jj];
gparts[jj] = temp_p;
}
}
/* Verify space_sort_struct. */
- /* for ( int k = i ; k <= jj ; k++ )
- if ( ind[k] > pivot ) {
- message( "sorting failed at k=%i, ind[k]=%i, pivot=%i, i=%i, j=%i,
- N=%i." , k , ind[k] , pivot , i , j , N );
- error( "Partition failed (<=pivot)." );
- }
- for ( int k = jj+1 ; k <= j ; k++ )
- if ( ind[k] <= pivot ) {
- message( "sorting failed at k=%i, ind[k]=%i, pivot=%i, i=%i, j=%i,
- N=%i." , k , ind[k] , pivot , i , j , N );
- error( "Partition failed (>pivot)." );
- } */
+ /* for (int k = i; k <= jj; k++)
+ if (ind[k] > pivot) {
+ message("sorting failed at k=%i, ind[k]=%i, pivot=%i, i=%i, j=%i.", k,
+ ind[k], pivot, i, j);
+ error("Partition failed (<=pivot).");
+ }
+ for (int k = jj + 1; k <= j; k++)
+ if (ind[k] <= pivot) {
+ message("sorting failed at k=%i, ind[k]=%i, pivot=%i, i=%i, j=%i.", k,
+ ind[k], pivot, i, j);
+ error("Partition failed (>pivot).");
+ } */
/* Split-off largest interval. */
if (jj - i > j - jj + 1) {
/* Recurse on the left? */
if (jj > i && pivot > min) {
- qid = (last++) % qstack_size;
- qstack[qid].i = i;
- qstack[qid].j = jj;
- qstack[qid].min = min;
- qstack[qid].max = pivot;
- qstack[qid].ready = 1;
- if ((waiting++) >= qstack_size) error("Qstack overflow.");
+ qid = atomic_inc(&space_sort_struct.last) %
+ space_sort_struct.stack_size;
+ while (space_sort_struct.stack[qid].ready)
+ ;
+ space_sort_struct.stack[qid].i = i;
+ space_sort_struct.stack[qid].j = jj;
+ space_sort_struct.stack[qid].min = min;
+ space_sort_struct.stack[qid].max = pivot;
+ if (atomic_inc(&space_sort_struct.waiting) >=
+ space_sort_struct.stack_size)
+ error("Qstack overflow.");
+ space_sort_struct.stack[qid].ready = 1;
}
/* Recurse on the right? */
@@ -838,13 +897,18 @@ void space_gparts_sort(struct gpart *gparts, int *ind, size_t N, int min,
/* Recurse on the right? */
if (pivot + 1 < max) {
- qid = (last++) % qstack_size;
- qstack[qid].i = jj + 1;
- qstack[qid].j = j;
- qstack[qid].min = pivot + 1;
- qstack[qid].max = max;
- qstack[qid].ready = 1;
- if ((waiting++) >= qstack_size) error("Qstack overflow.");
+ qid = atomic_inc(&space_sort_struct.last) %
+ space_sort_struct.stack_size;
+ while (space_sort_struct.stack[qid].ready)
+ ;
+ space_sort_struct.stack[qid].i = jj + 1;
+ space_sort_struct.stack[qid].j = j;
+ space_sort_struct.stack[qid].min = pivot + 1;
+ space_sort_struct.stack[qid].max = max;
+ if (atomic_inc(&space_sort_struct.waiting) >=
+ space_sort_struct.stack_size)
+ error("Qstack overflow.");
+ space_sort_struct.stack[qid].ready = 1;
}
/* Recurse on the left? */
@@ -857,18 +921,9 @@ void space_gparts_sort(struct gpart *gparts, int *ind, size_t N, int min,
} /* loop over sub-intervals. */
- waiting--;
+ atomic_dec(&space_sort_struct.waiting);
} /* main loop. */
-
- /* Verify space_sort_struct. */
- /* for ( i = 1 ; i < N ; i++ )
- if ( ind[i-1] > ind[i] )
- error( "Sorting failed (ind[%i]=%i,ind[%i]=%i)." , i-1 , ind[i-1] , i
- , ind[i] ); */
-
- /* Clean up. */
- free(qstack);
}
/**
diff --git a/src/space.h b/src/space.h
index e761595838ae78b0d8a67cca676cfa59f3f700f6..8eb166d82fac2e1a551665680af95c2cf4a188ac 100644
--- a/src/space.h
+++ b/src/space.h
@@ -115,6 +115,7 @@ struct qstack {
};
struct parallel_sort {
struct part *parts;
+ struct gpart *gparts;
struct xpart *xparts;
int *ind;
struct qstack *stack;
@@ -126,8 +127,8 @@ extern struct parallel_sort space_sort_struct;
/* function prototypes. */
void space_parts_sort(struct space *s, int *ind, size_t N, int min, int max,
int verbose);
-void space_gparts_sort(struct gpart *gparts, int *ind, size_t N, int min,
- int max);
+void space_gparts_sort(struct space *s, int *ind, size_t N, int min, int max,
+ int verbose);
struct cell *space_getcell(struct space *s);
int space_getsid(struct space *s, struct cell **ci, struct cell **cj,
double *shift);
@@ -149,5 +150,6 @@ void space_recycle(struct space *s, struct cell *c);
void space_split(struct space *s, struct cell *cells, int verbose);
void space_do_split(struct space *s, struct cell *c);
void space_do_parts_sort();
+void space_do_gparts_sort();
void space_link_cleanup(struct space *s);
#endif /* SWIFT_SPACE_H */
diff --git a/src/swift.h b/src/swift.h
index 9ab090dccd195ff4927d3e614e446b36d273f824..e568a28c888295affc9ec45b6d059d34f5b4bf04 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -27,7 +27,6 @@
#include "cell.h"
#include "clocks.h"
#include "const.h"
-#include "const.h"
#include "cycle.h"
#include "debug.h"
#include "engine.h"
@@ -38,7 +37,9 @@
#include "map.h"
#include "multipole.h"
#include "parallel_io.h"
+#include "parser.h"
#include "part.h"
+#include "partition.h"
#include "queue.h"
#include "runner.h"
#include "scheduler.h"
@@ -47,9 +48,8 @@
#include "space.h"
#include "task.h"
#include "timers.h"
-#include "units.h"
#include "tools.h"
-#include "partition.h"
+#include "units.h"
#include "version.h"
#endif /* SWIFT_SWIFT_H */
diff --git a/src/task.c b/src/task.c
index 74f8451d5dfaec2454f2eeed8670765e1be5b658..5f1475a46e4626e1f51db673d73fd84f86e6edb6 100644
--- a/src/task.c
+++ b/src/task.c
@@ -43,9 +43,10 @@
/* Task type names. */
const char *taskID_names[task_type_count] = {
- "none", "sort", "self", "pair", "sub", "init",
- "ghost", "drift", "kick", "send", "recv", "grav_pp",
- "grav_mm", "grav_up", "grav_down", "psort", "split_cell", "rewait"};
+ "none", "sort", "self", "pair", "sub",
+ "init", "ghost", "drift", "kick", "send",
+ "recv", "grav_pp", "grav_mm", "grav_up", "grav_down",
+ "part_sort", "gpart_sort", "split_cell", "rewait"};
const char *subtaskID_names[task_type_count] = {"none", "density",
"force", "grav"};
@@ -78,9 +79,10 @@ float task_overlap(const struct task *ta, const struct task *tb) {
/* First check if any of the two tasks are of a type that don't
use cells. */
if (ta == NULL || tb == NULL || ta->type == task_type_none ||
- ta->type == task_type_psort || ta->type == task_type_split_cell ||
- ta->type == task_type_rewait || tb->type == task_type_none ||
- tb->type == task_type_psort || tb->type == task_type_split_cell ||
+ ta->type == task_type_part_sort || ta->type == task_type_gpart_sort ||
+ ta->type == task_type_split_cell || ta->type == task_type_rewait ||
+ tb->type == task_type_none || tb->type == task_type_part_sort ||
+ tb->type == task_type_gpart_sort || tb->type == task_type_split_cell ||
tb->type == task_type_rewait)
return 0.0f;
diff --git a/src/task.h b/src/task.h
index b86631cc49bfad102302e3bab380bfb5eb8ed1e0..9c0ba6087d772d7362a98bc40a838c6fa3713166 100644
--- a/src/task.h
+++ b/src/task.h
@@ -45,7 +45,8 @@ enum task_types {
task_type_grav_mm,
task_type_grav_up,
task_type_grav_down,
- task_type_psort,
+ task_type_part_sort,
+ task_type_gpart_sort,
task_type_split_cell,
task_type_rewait,
task_type_count
diff --git a/src/tools.c b/src/tools.c
index a384974fdc94452079838ae0eebf30e1815f644b..d25b7401a1e0515c650333b41193d54b5e155d39 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -300,9 +300,9 @@ void pairs_single_grav(double *dim, long long int pid,
break;
if (k == N) error("Part not found.");
pi = parts[k];
- pi.a[0] = 0.0f;
- pi.a[1] = 0.0f;
- pi.a[2] = 0.0f;
+ pi.a_grav[0] = 0.0f;
+ pi.a_grav[1] = 0.0f;
+ pi.a_grav[2] = 0.0f;
/* Loop over all particle pairs. */
for (k = 0; k < N; k++) {
@@ -320,15 +320,15 @@ void pairs_single_grav(double *dim, long long int pid,
}
r2 = fdx[0] * fdx[0] + fdx[1] * fdx[1] + fdx[2] * fdx[2];
runner_iact_grav(r2, fdx, &pi, &pj);
- a[0] += pi.a[0];
- a[1] += pi.a[1];
- a[2] += pi.a[2];
- aabs[0] += fabsf(pi.a[0]);
- aabs[1] += fabsf(pi.a[1]);
- aabs[2] += fabsf(pi.a[2]);
- pi.a[0] = 0.0f;
- pi.a[1] = 0.0f;
- pi.a[2] = 0.0f;
+ a[0] += pi.a_grav[0];
+ a[1] += pi.a_grav[1];
+ a[2] += pi.a_grav[2];
+ aabs[0] += fabsf(pi.a_grav[0]);
+ aabs[1] += fabsf(pi.a_grav[1]);
+ aabs[2] += fabsf(pi.a_grav[2]);
+ pi.a_grav[0] = 0.0f;
+ pi.a_grav[1] = 0.0f;
+ pi.a_grav[2] = 0.0f;
}
/* Dump the result. */
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 7edf6e0659c91e66c7eee3e3815e8f83e8551ec2..b53a08615c5a8c7c2c31475bf7207522f8b9a58c 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -22,11 +22,11 @@ AM_LDFLAGS = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
# List of programs and scripts to run in the test suite
TESTS = testGreetings testReading.sh testSingle testPair.sh testPairPerturbed.sh \
- test27cells.sh test27cellsPerturbed.sh
+ test27cells.sh test27cellsPerturbed.sh testParser.sh
# List of test programs to compile
check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
- testSPHStep testPair test27cells
+ testSPHStep testPair test27cells testParser testKernel
# Sources for the individual programs
testGreetings_SOURCES = testGreetings.c
@@ -43,6 +43,11 @@ testPair_SOURCES = testPair.c
test27cells_SOURCES = test27cells.c
+testParser_SOURCES = testParser.c
+
+testKernel_SOURCES = testKernel.c
+
# Files necessary for distribution
EXTRA_DIST = testReading.sh makeInput.py testPair.sh testPairPerturbed.sh \
- test27cells.sh test27cellsPerturbed.sh tolerance.dat
+ test27cells.sh test27cellsPerturbed.sh tolerance.dat testParser.sh \
+ testParserInput.yaml
diff --git a/tests/test27cells.c b/tests/test27cells.c
index 74c38996a81056b10633bf2bbf18cc7cff7e8f0d..7915511eed50a229a94eda6bb338607099303421 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -24,18 +24,37 @@
#include <unistd.h>
#include "swift.h"
+enum velocity_types {
+ velocity_zero,
+ velocity_random,
+ velocity_divergent,
+ velocity_rotating
+};
+
/**
- * Returns a random number (uniformly distributed) in [a,b[
+ * @brief Returns a random number (uniformly distributed) in [a,b[
*/
double random_uniform(double a, double b) {
return (rand() / (double)RAND_MAX) * (b - a) + a;
}
-/* n is both particles per axis and box size:
- * particles are generated on a mesh with unit spacing
+
+/**
+ * @brief Constructs a cell and all of its particle in a valid state prior to
+ * a DOPAIR or DOSELF calcuation.
+ *
+ * @param n The cube root of the number of particles.
+ * @param offset The position of the cell offset from (0,0,0).
+ * @param size The cell size.
+ * @param h The smoothing length of the particles in units of the inter-particle separation.
+ * @param density The density of the fluid.
+ * @param partId The running counter of IDs.
+ * @param pert The perturbation to apply to the particles in the cell in units of the inter-particle separation.
+ * @param vel The type of velocity field (0, random, divergent, rotating)
*/
struct cell *make_cell(size_t n, double *offset, double size, double h,
- double density, long long *partId, double pert) {
+ double density, long long *partId, double pert,
+ enum velocity_types vel) {
const size_t count = n * n * n;
const double volume = size * size * size;
struct cell *cell = malloc(sizeof(struct cell));
@@ -61,12 +80,28 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
part->x[2] =
offset[2] +
size * (z + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
- // part->v[0] = part->x[0] - 1.5;
- // part->v[1] = part->x[1] - 1.5;
- // part->v[2] = part->x[2] - 1.5;
- part->v[0] = random_uniform(-0.05, 0.05);
- part->v[1] = random_uniform(-0.05, 0.05);
- part->v[2] = random_uniform(-0.05, 0.05);
+ switch (vel) {
+ case velocity_zero:
+ part->v[0] = 0.f;
+ part->v[1] = 0.f;
+ part->v[2] = 0.f;
+ break;
+ case velocity_random:
+ part->v[0] = random_uniform(-0.05, 0.05);
+ part->v[1] = random_uniform(-0.05, 0.05);
+ part->v[2] = random_uniform(-0.05, 0.05);
+ break;
+ case velocity_divergent:
+ part->v[0] = part->x[0] - 1.5 * size;
+ part->v[1] = part->x[1] - 1.5 * size;
+ part->v[2] = part->x[2] - 1.5 * size;
+ break;
+ case velocity_rotating:
+ part->v[0] = part->x[1];
+ part->v[1] = -part->x[0];
+ part->v[2] = 0.f;
+ break;
+ }
part->h = size * h / (float)n;
part->id = ++(*partId);
part->mass = density * volume / count;
@@ -209,10 +244,11 @@ void runner_doself1_density(struct runner *r, struct cell *ci);
int main(int argc, char *argv[]) {
size_t runs = 0, particles = 0;
- double h = 1.12575, size = 1., rho = 1.;
+ double h = 1.2348, size = 1., rho = 1.;
double perturbation = 0.;
char outputFileNameExtension[200] = "";
char outputFileName[200] = "";
+ int vel = velocity_zero;
/* Initialize CPU frequency, this also starts time. */
unsigned long long cpufreq = 0;
@@ -222,7 +258,7 @@ int main(int argc, char *argv[]) {
srand(0);
char c;
- while ((c = getopt(argc, argv, "m:s:h:p:r:t:d:f:")) != -1) {
+ while ((c = getopt(argc, argv, "m:s:h:p:r:t:d:f:v:")) != -1) {
switch (c) {
case 'h':
sscanf(optarg, "%lf", &h);
@@ -245,6 +281,9 @@ int main(int argc, char *argv[]) {
case 'f':
strcpy(outputFileNameExtension, optarg);
break;
+ case 'v':
+ sscanf(optarg, "%d", &vel);
+ break;
case '?':
error("Unknown option.");
break;
@@ -257,18 +296,25 @@ int main(int argc, char *argv[]) {
"\nGenerates a cell pair, filled with particles on a Cartesian grid."
"\nThese are then interacted using runner_dopair1_density."
"\n\nOptions:"
- "\n-h DISTANCE=1.1255 - Smoothing length"
+ "\n-h DISTANCE=1.2348 - Smoothing length in units of <x>"
"\n-m rho - Physical density in the cell"
"\n-s size - Physical size of the cell"
"\n-d pert - Perturbation to apply to the particles [0,1["
+ "\n-v type (0,1,2,3) - Velocity field: (zero, random, divergent, "
+ "rotating)"
"\n-f fileName - Part of the file name used to save the dumps\n",
argv[0]);
exit(1);
}
/* Help users... */
- message("Smoothing length: h = %f", h);
- message("Neighbour target: N = %f", kernel_nwneigh);
+ message("Smoothing length: h = %f", h * size);
+ message("Kernel: %s", kernel_name);
+ message("Neighbour target: N = %f", h * h * h * kernel_nwneigh / 1.88273);
+ message("Density target: rho = %f", rho);
+ message("div_v target: div = %f", vel == 2 ? 3.f : 0.f);
+ message("curl_v target: curl = [0., 0., %f]", vel == 3 ? -2.f : 0.f);
+ printf("\n");
/* Build the infrastructure */
struct space space;
@@ -292,8 +338,8 @@ int main(int argc, char *argv[]) {
for (int k = 0; k < 3; ++k) {
double offset[3] = {i * size, j * size, k * size};
- cells[i * 9 + j * 3 + k] =
- make_cell(particles, offset, size, h, rho, &partId, perturbation);
+ cells[i * 9 + j * 3 + k] = make_cell(particles, offset, size, h, rho,
+ &partId, perturbation, vel);
}
}
}
diff --git a/tests/testKernel.c b/tests/testKernel.c
new file mode 100644
index 0000000000000000000000000000000000000000..5ad9cc81ea92e6ef9487489c5d560abf414e38df
--- /dev/null
+++ b/tests/testKernel.c
@@ -0,0 +1,37 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#include "swift.h"
+
+int main() {
+
+ const float h = const_eta_kernel;
+ const int numPoints = 30;
+
+ for (int i = 0; i < numPoints; ++i) {
+
+ const float x = i * 3.f / numPoints;
+ float W, dW;
+ kernel_deval(x / h, &W, &dW);
+
+ printf("h= %f H= %f x=%f W(x,h)=%f\n", h, h * kernel_gamma, x, W);
+ }
+
+ return 0;
+}
diff --git a/tests/testParser.c b/tests/testParser.c
new file mode 100644
index 0000000000000000000000000000000000000000..a4b8789fca056fef659bca78eae9d0effb2ceb66
--- /dev/null
+++ b/tests/testParser.c
@@ -0,0 +1,67 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2016 James Willis (james.s.willis@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/>.
+ *
+ ******************************************************************************/
+
+#include "parser.h"
+#include <assert.h>
+#include <string.h>
+#include <math.h>
+
+int main(int argc, char *argv[]) {
+
+ const char *input_file = argv[1];
+
+ /* Create a structure to read file into. */
+ struct swift_params param_file;
+
+ /* Create variables that will be set from the parameter file. */
+ int no_of_threads = 0;
+ int no_of_time_steps = 0;
+ float max_h = 0.0f;
+ double start_time = 0.0;
+ char ic_file[PARSER_MAX_LINE_SIZE];
+
+ /* Read the parameter file. */
+ parser_read_file(input_file, ¶m_file);
+
+ /* Print the contents of the structure. */
+ parser_print_params(¶m_file);
+
+ /* Retrieve parameters and store them in variables defined above.
+ * Have to specify the name of the parameter as it appears in the
+ * input file: testParserInput.yaml.*/
+ parser_get_param_int(¶m_file, "no_of_threads", &no_of_threads);
+ parser_get_param_int(¶m_file, "no_of_time_steps", &no_of_time_steps);
+ parser_get_param_float(¶m_file, "max_h", &max_h);
+ parser_get_param_double(¶m_file, "start_time", &start_time);
+ parser_get_param_string(¶m_file, "ic_file", ic_file);
+
+ /* Print the variables to check their values are correct. */
+ printf(
+ "no_of_threads: %d, no_of_time_steps: %d, max_h: %f, start_time: %lf, "
+ "ic_file: %s\n",
+ no_of_threads, no_of_time_steps, max_h, start_time, ic_file);
+
+ assert(no_of_threads == 16);
+ assert(no_of_time_steps == 10);
+ assert(fabs(max_h - 1.1255) < 0.00001);
+ assert(fabs(start_time - 1.23456789) < 0.00001);
+ assert(strcmp(ic_file, "ic_file.ini") == 0); /*strcmp returns 0 if correct.*/
+
+ return 0;
+}
diff --git a/tests/testParser.sh b/tests/testParser.sh
new file mode 100755
index 0000000000000000000000000000000000000000..3dad7f386f792ff2beb6e94eb093bad4085023a4
--- /dev/null
+++ b/tests/testParser.sh
@@ -0,0 +1,3 @@
+#!/bin/bash
+
+./testParser testParserInput.yaml
diff --git a/tests/testParserInput.yaml b/tests/testParserInput.yaml
new file mode 100644
index 0000000000000000000000000000000000000000..d695e6a8ddd327e31224f36a6e34767ea8d36408
--- /dev/null
+++ b/tests/testParserInput.yaml
@@ -0,0 +1,9 @@
+---
+no_of_threads: 16 # The number of threads that will be used.
+no_of_time_steps: 10
+max_h: 1.1255
+start_time: 1.23456789
+#Input file
+ic_file: ic_file.ini
+
+...
diff --git a/tests/tolerance.dat b/tests/tolerance.dat
index 48de4383eab6214812183be25d3036a324ccbc27..f5031c5f47dfa203300ebcc9a47fbac42f854d26 100644
--- a/tests/tolerance.dat
+++ b/tests/tolerance.dat
@@ -1,3 +1,3 @@
# ID pos_x pos_y pos_z v_x v_y v_z rho rho_dh wcount wcount_dh div_v curl_vx curl_vy curl_vz
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-5 1e-5 2e-5 3e-4 1e-5 1e-5 1e-5 1e-5
- 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-5 1.2e-5 1e-5 1e-5 1e-4 1e-4 1e-4 1e-4
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-5 1e-5 2e-5 3e-2 1e-5 1e-5 1e-5 1e-5
+ 0 1e-6 1e-6 1e-6 1e-6 1e-6 1e-6 1e-5 1.2e-5 1e-5 1e-2 1e-4 1e-4 1e-4 1e-4
diff --git a/theory/kernel/kernel.pdf b/theory/kernel/kernel.pdf
deleted file mode 100644
index b6e540dc61c36dd00e56f02e44ce33f9f91a7f01..0000000000000000000000000000000000000000
Binary files a/theory/kernel/kernel.pdf and /dev/null differ
diff --git a/theory/kernel/kernel.tex b/theory/kernel/kernel.tex
deleted file mode 100644
index 7087555d423afbe2745bb91010a17c52a32084f2..0000000000000000000000000000000000000000
--- a/theory/kernel/kernel.tex
+++ /dev/null
@@ -1,155 +0,0 @@
-\documentclass[a4paper,10pt]{article}
-\usepackage[utf8]{inputenc}
-\usepackage{amsmath}
-
-%opening
-\title{SPH kernels in SWIFT}
-\author{Matthieu Schaller}
-
-\begin{document}
-
-\maketitle
-
-\section{General Definitions}
-
-The smoothing kernels used in SPH are almost always isotropic and can hence be written in 3D as
-
-\begin{equation}
- W(\vec{x},h) = \frac{1}{h^3}f\left(\frac{|\vec{x}|}{h}\right),
-\end{equation}
-
-where $f(q)$ is a dimensionless function, usually a low-order polynomial, normalized to unity. For computational
-reasons, this kernel
-usually has a finite support of radius $H$. In other words,
-
-\begin{equation}
- W(\vec{x},h) = 0\quad \forall\quad |\vec{x}| > H.
-\end{equation}
- One can then define the weighted number of neighbours within $H$ as
-
-\begin{equation}
- N_{ngb} = \frac{4}{3}\pi H^3 \sum_j W(\vec{x}_i - \vec{x}_j,h).
-\end{equation}
-
-The value of $N_{ngb}$ is often used in the codes to find the smoothing length of each particle via Newton iterations
-or a bissection algorithm. $H$ is defined as \emph{the smoothing length} in the GADGET code. This definition is useful
-for implementation reasons but does not really correspond to a true physical quantity. \\
-The main question is the definition of the smoothing length. The function $W(\vec{x},h)$ is invariant under the
-rescaling $h\rightarrow \alpha h,~f(q)\rightarrow\alpha^{-3}f(\alpha q)$, which makes the definition of $h$ difficult.
-This ambiguity is present in the litterature with authors using different definition of the \emph{physical} smoothing
-length, $h=\frac{1}{2}H$ or $h=H$ for instance. \\
-A more physically motivated estimate is the standard deviation of the kernel:
-
-\begin{equation}
- \sigma^2 = \frac{1}{3} \int \vec{x}^2~W(\vec{x},h)~d^3\vec{x}
-\end{equation}
-
-which then allows us to set $h=2\sigma$. This definition of the smoothing length is more physical as one can
-demonstrate that the reconstruction of any smooth field $A(\vec{x})$ using interpolation of particles at the point
-$\vec{x}_i$ can be expanded as
-
-\begin{equation}
-A_i \approx A(\vec{x}_i) + \frac{1}{2}\sigma^2 \nabla^2A(\vec{x}_i) + \mathcal{O}\left(\sigma^4\right).
-\end{equation}
-
-The quantity $H/\sigma$ is independant of the choice of $h$ made and is purely a functional of $f(q)$. The number of
-neighbours (used in the code to construct the neighborhood of a given particle) can then be expressed as a function of
-this \emph{physical} $h$ (or $\sigma$). Or to relate it even more to the particle distribution, we can write
-$h=\eta\Delta x$, with $\Delta x$ the mean inter-particle separation:
-
-\begin{equation}
- N_{ngb} = \frac{4}{3}\pi \left(\frac{1}{2}\eta\frac{H}{\sigma}\right)^3 = \frac{4}{3}\pi
-\left(\eta\zeta\right)^3
-\end{equation}
-
-This definition of the number of neighbours only depends on $f(q)$ (via $\zeta$) and on the mean inter-particle
-separation. The problem is then fully specified by specifying a form for $f(q)$ and $\eta$. \\
-Experiments suggest that $\eta \approx 1.2 - 1.3$ is a reasonnable choice. The bigger $\eta$, the better the smoothing
-and hence the better the reconstruction of the field. This, however, comes at a higher computational cost as more
-interactions between neighbours will have to be computed. Also, spline kernels become instable when $\eta>1.5$.
-
-\section{Kernels available in SWIFT}
-
-The different kernels available are listed below.
-\paragraph{Cubic Spline Kernel}
-\begin{equation*}
- f(q) = \frac{1}{\pi}\left\lbrace \begin{array}{rcl}
- \frac{3}{4}q^3 - 15q^2 + 1 & \mbox{if} & 0 \leq q < 1 \\
- -\frac{1}{4}q^3 + \frac{3}{2}q^2-3q+2 & \mbox{if} & 1 \leq q < 2 \\
- 0 & \mbox{if} & q \geq 2 \\
- \end{array}\right.
-\end{equation*}
-with $\zeta = \frac{1}{2}\sqrt{\frac{40}{3}} \approx 1.825742$. Thus, for a resolution of $\eta = 1.235$, this kernel
-uses $N_{ngb} \approx 48$. The code uses $h = \frac{1}{2}H = \zeta \sigma$ internally.
-
-\paragraph{Quartic Spline Kernel}
-\begin{equation*}
- f(q) = \frac{1}{20\pi}\left\lbrace \begin{array}{rcl}
- 6q^4 - 15q^2 + \frac{115}{8} & \mbox{if} & 0 \leq q < \frac{1}{2} \\
- -4q^4 + 20q^3-30q^2 + 5q + \frac{55}{4} & \mbox{if} & \frac{1}{2} \leq q < \frac{3}{2} \\
- q^4-10q^3+\frac{75}{2}q^2-\frac{125}{2}q+\frac{625}{16} & \mbox{if} & \frac{3}{2} \leq q <
-\frac{5}{2} \\
- 0 & \mbox{if} & q \geq \frac{5}{2} \\
- \end{array}\right.
-\end{equation*}
-with $\zeta = \frac{1}{2}\sqrt{\frac{375}{23}} \approx 2.018932$. Thus, for a resolution of $\eta = 1.235$, this kernel
-uses $N_{ngb} \approx 64.9$. The code uses $h = \frac{2}{5}H =\frac{4}{5}\zeta \sigma$ internally.
-
-\paragraph{Quintic Spline Kernel}
-\begin{equation*}
- f(q) = \frac{1}{120\pi}\left\lbrace \begin{array}{rcl}
- -10q^5 + 30q^4 - 60q^2 + 66 & \mbox{if} & 0 \leq q < 1 \\
- 5q^5 - 45q^4 + 150q^3 - 210q^2 + 75q + 51 & \mbox{if} & 1 \leq q < 2 \\
- -q^5 + 15q^4 - 90q^3 + 270q^2 - 405q + 243 & \mbox{if} & 2 \leq q < 3 \\
- 0 & \mbox{if} & q \geq 3 \\
- \end{array}\right.
-\end{equation*}
-with $\zeta = \frac{1}{2}\sqrt{\frac{135}{7}} \approx 2.195775$. Thus, for a resolution of $\eta = 1.235$, this kernel
-uses $N_{ngb} \approx 83.5$. The code uses $h = \frac{1}{3}H = \frac{2}{3}\zeta \sigma$ internally.
-
-\paragraph{Wendland $C2$ Kernel}
-\begin{equation*}
- f(q) = \frac{21}{2\pi}\left\lbrace \begin{array}{rcl}
- 4 q^5-15 q^4+20 q^3-10 q^2+1 & \mbox{if} & 0 \leq q < 1 \\
- 0 & \mbox{if} & q \geq 1 \\
- \end{array}\right.
-\end{equation*}
- with $\zeta = \frac{1}{2}\sqrt{15} \approx 1.93649$. Thus, for a resolution of $\eta = 1.235$, this kernel
-uses $N_{ngb} \approx 57.3$. The code uses $h = H = 2\zeta \sigma$ internally.
-
-
-\paragraph{Wendland $C4$ Kernel}
-\begin{equation*}
- f(q) = \frac{495}{32\pi}\left\lbrace \begin{array}{rcl}
- \frac{35}{3} q^8-64 q^7+ 140 q^6-\frac{448}{3} q^5+70 q^4-\frac{28}{3} q^2+1 & \mbox{if} & 0
-\leq q < 1 \\
- 0 & \mbox{if} & q \geq 1 \\
- \end{array}\right.
-\end{equation*}
- with $\zeta = \frac{1}{2}\sqrt{\frac{39}{2}} \approx 2.20794$. Thus, for a resolution of $\eta = 1.235$, this kernel
-uses $N_{ngb} \approx 84.9$. The code uses $h = H = 2\zeta \sigma$ internally.
-
-\paragraph{Wendland $C6$ Kernel}
-\begin{equation*}
- f(q) = \frac{1365}{64\pi}\left\lbrace \begin{array}{rcl}
- 32 q^{11}-231 q^{10}+704 q^9-1155 q^8+1056 q^7-462 q^6+66 q^4-11 q^2+1 & \mbox{if} & 0
-\leq q < 1 \\
- 0 & \mbox{if} & q \geq 1 \\
- \end{array}\right.
-\end{equation*}
- with $\zeta = \frac{1}{2}\sqrt{24} \approx 2.44949$. Thus, for a resolution of $\eta = 1.235$, this kernel
-uses $N_{ngb} \approx 116$. The code uses $h = H = 2\zeta \sigma$ internally.
-
-\section{Kernel Derivatives}
-
-The derivatives of the kernel function have relatively simple expressions:
-
-\begin{eqnarray*}
- \vec\nabla_x W(\vec{x},h) &=& \frac{1}{h^4}f'\left(\frac{|\vec{x}|}{h}\right) \frac{\vec{x}}{|\vec{x}|} \\
- \frac{\partial W(\vec{x},h)}{\partial h} &=&- \frac{1}{h^4}\left[3f\left(\frac{|\vec{x}|}{h}\right) +
-\frac{|\vec{x}|}{h}f'\left(\frac{|\vec{x}|}{h}\right)\right]
-\end{eqnarray*}
-
-Note that for all the kernels listed above, $f'(0) = 0$.
-
-\end{document}
diff --git a/theory/kernel/kernel_definitions.tex b/theory/kernel/kernel_definitions.tex
new file mode 100644
index 0000000000000000000000000000000000000000..8999636109ffadcbf148ce3c1fbccdc44feafe65
--- /dev/null
+++ b/theory/kernel/kernel_definitions.tex
@@ -0,0 +1,242 @@
+\documentclass[a4paper]{mnras}
+\usepackage[utf8]{inputenc}
+\usepackage{amsmath}
+\usepackage{graphicx}
+\usepackage{xspace}
+
+\newcommand{\swift}{{\sc Swift}\xspace}
+
+
+
+%opening
+\title{SPH kernels in SWIFT}
+\author{Matthieu Schaller}
+
+\begin{document}
+
+\maketitle
+
+In here we follow the definitions of Dehnen \& Aly 2012.
+
+\section{General Definitions}
+
+The desirable properties of an SPH kernels $W(\vec{x},h)$ are:
+\begin{enumerate}
+\item $W(\vec{x},h)$ should be isotropic in $\vec{x}$.
+\item $W(\vec{x},h)$ should be positive and decrease monotonically.
+\item $W(\vec{x},h)$ should be twice differentiable.
+\item $W(\vec{x},h)$ should have a finite support and be cheap to
+ compute.
+\end{enumerate}
+
+As a consequence, the smoothing kernels used in SPH can
+hence be written (in 3D) as
+
+\begin{equation}
+ W(\vec{x},h) \equiv \frac{1}{H^3}f\left(\frac{|\vec{x}|}{H}\right),
+\end{equation}
+
+where $H=\gamma h$ is defined below and $f(u)$ is a dimensionless
+function, usually a low-order polynomial, such that $f(u \geq 1) = 0$
+and normalised such that
+
+\begin{equation}
+ \int f(|\vec{u}|){\rm d}^3u = 1.
+\end{equation}
+
+$H$ is the kernel's support radius and is used as the ``smoothing
+length'' in the Gadget code( {i.e.} $H=h$). This definition is,
+however, not very physical and makes comparison of kernels at a
+\emph{fixed resolution} difficult. A more sensible definition of the
+smoothing length, related to the Taylor expansion of the
+re-constructed density field is given in terms of the kernel's
+standard deviation
+
+\begin{equation}
+ \sigma^2 \equiv \frac{1}{3}\int \vec{u}^2 W(\vec{u},h) {\rm d}^3u.
+ \label{eq:sph:sigma}
+\end{equation}
+
+The smoothing length is then:
+\begin{equation}
+ h\equiv2\sigma.
+ \label{eq:sph:h}
+\end{equation}
+
+Each kernel, {\it i.e.} defintion of $f(u)$, will have a different
+ratio $\gamma = H/h$. So for a \emph{fixed resolution} $h$, one will
+have different kernel support sizes, $H$, and a different number of
+neighbours, $N_{\rm ngb}$ to interact with. One would typically choose
+$h$ for a simulation as a multiple $\eta$ of the mean-interparticle
+separation:
+
+\begin{equation}
+ h = \eta \langle x \rangle = \eta \left(\frac{m}{\rho}\right)^{1/3},
+\end{equation}
+
+where $\rho$ is the local density of the fluid and $m$ the SPH
+particle mass.
+
+The (weighted) number of neighbours within the kernel support is a
+useful quantity to use in implementations of SPH. It is defined as (in
+3D)
+
+\begin{equation}
+ N_{\rm ngb} \equiv \frac{4}{3}\pi \left(\frac{H}{h}\eta\right)^3.
+\end{equation}
+
+Once the fixed ratio $\gamma= H/h$ is known (via equations
+\ref{eq:sph:sigma} and \ref{eq:sph:h}) for a given kernel, the number
+of neighbours only depends on the resolution parameter $\eta$. For
+the usual cubic spline kernel (see below), setting the simulation
+resolution to $\eta=1.2348$ yields the commonly used value $N_{\rm
+ ngb} = 48$.
+
+\section{Kernels available in \swift}
+
+The \swift kernels are split into two categories, the B-splines
+($M_{4,5,6}$) and the Wendland kernels ($C2$, $C4$ and $C6$). In all
+cases we impose $f(u>1) = 0$.\\
+
+The spline kernels are defined as:
+
+\begin{align}
+ f(u) &= C M_n(u), \\
+ M_n(u) &\equiv \frac{1}{2\pi}
+ \int_{-\infty}^{\infty}
+ \left(\frac{\sin\left(k/n\right)}{k/n}\right)^n\cos\left(ku\right){\rm
+ d}k,
+\end{align}
+
+whilst the Wendland kernels read
+
+\begin{align}
+ f(u) &= C \Psi_{i,j}(u), \\
+ \Psi_{i,j}(u) &\equiv \mathcal{I}^k\left[\max\left(1-u,0\right)\right],\\
+ \mathcal{I}[f](u) &\equiv \int_u^\infty f\left(k\right)k{\rm d}k.
+\end{align}
+
+\subsubsection{Cubic spline ($M_4$) kernel}
+
+In 3D, we have $C=\frac{16}{\pi}$ and $\gamma=H/h = 1.825742$.\\
+The kernel function $f(u)$ reads:
+
+\begin{equation}
+ M_4(u) = \left\lbrace\begin{array}{rcl}
+ 3u^3 - 3u^2 + \frac{1}{2} & \mbox{if} & u<\frac{1}{2}\\
+ -u^3 + 3u^2 - 3u + 1 & \mbox{if} & u \geq \frac{1}{2}
+ \end{array}
+ \right.
+ \nonumber
+\end{equation}
+
+
+\subsubsection{Quartic spline ($M_5$) kernel}
+
+In 3D, we have $C=\frac{15625}{512\pi}$ and $\gamma=H/h = 2.018932$.\\
+The kernel function $f(u)$ reads:
+
+\begin{align}
+ M_5(u) &= \nonumber\\
+ &\left\lbrace\begin{array}{rcl}
+ 6u^4 - \frac{12}{5}u^2 + \frac{46}{125} & \mbox{if} & u < \frac{1}{5} \\
+ -4u^4 + 8u^3 - \frac{24}{5}u^2 + \frac{8}{25}u + \frac{44}{125} & \mbox{if} & \frac{1}{5} \leq u < \frac{3}{5}\\
+ u^4 - 4u^3 + 6u^2 - 4u + 1 & \mbox{if} & \frac{3}{5} \leq u \\
+ \end{array}
+ \right.
+ \nonumber
+\end{align}
+
+
+\subsubsection{Quintic spline ($M_6$) kernel}
+
+In 3D, we have $C=\frac{2187}{40\pi}$ and $\gamma=H/h = 2.195775$.\\
+The kernel function $f(u)$ reads:
+
+\begin{align}
+ M_6(u) &= \nonumber\\
+ &\left\lbrace\begin{array}{rcl}
+ -10u^5 + 10u^4 - \frac{20}{9}u^2 + \frac{22}{81} & \mbox{if} & u < \frac{1}{3} \\
+ 5u^5 - 15u^4 + \frac{50}{3}u^3 - \frac{70}{9}u^2 + \frac{25}{27}u + \frac{17}{81} & \mbox{if} & \frac{1}{3} \leq u < \frac{2}{3}\\
+ -1u^5 + 5u^4 - 10u^3 + 10u^2 - 5u + 1. & \mbox{if} & u \geq \frac{2}{3}
+ \end{array}
+ \right.
+ \nonumber
+\end{align}
+
+
+\subsubsection{Wendland C2 kernel}
+
+In 3D, we have $C=\frac{21}{2\pi}$ and $\gamma=H/h = 1.936492$.\\
+The kernel function $f(u)$ reads:
+
+\begin{align}
+ \Psi_{i,j}(u) &= 4u^5 - 15u^4 + 20u^3 - 10u^2 + 1.
+ \nonumber
+\end{align}
+
+
+\subsubsection{Wendland C4 kernel}
+
+In 3D, we have $C=\frac{495}{32\pi}$ and $\gamma=H/h = 2.207940$.\\
+The kernel function $f(u)$ reads:
+
+\begin{align}
+ \Psi_{i,j}(u) &= \frac{35}{3}u^8 - 64u^7 + 140u^6 \nonumber\\
+ & - \frac{448}{3}u^5 + 70u^4 - \frac{28}{3}u^2 + 1
+ \nonumber
+\end{align}
+
+
+\subsubsection{Wendland C6 kernel}
+
+In 3D, we have $C=\frac{1365}{64\pi}$ and $\gamma=H/h = 2.449490$.\\
+The kernel function $f(u)$ reads:
+
+\begin{align}
+ \Psi_{i,j}(u) &= 32u^{11} - 231u^{10} + 704u^9 - 1155u^8 \nonumber\\
+ & + 1056u^7 - 462u^6 + 66u^4 - 11u^2 + 1
+ \nonumber
+\end{align}
+
+
+\subsubsection{Summary}
+
+All kernels available in \swift are shown on Fig.~\ref{fig:sph:kernels}.
+
+\begin{figure}
+\includegraphics[width=\columnwidth]{kernels.pdf}
+\caption{The kernel functions available in \swift for a mean
+ inter-particle separation $\langle x\rangle=1.5$ and a resolution
+ $\eta=1.2348$. The corresponding kernel support radii $H$ (shown by
+ arrows) and number of neighours $N_{\rm ngb}$ are indicated on the
+ figure. A Gaussian kernel with the same smoothing length is shown
+ for comparison. Note that all these kernels have the \emph{same
+ resolution} despite having vastly different number of neighbours.}
+\label{fig:sph:kernels}
+\end{figure}
+
+\begin{figure}
+\includegraphics[width=\columnwidth]{kernel_derivatives.pdf}
+\caption{The first and secon derivatives of the kernel functions
+ available in \swift for a mean inter-particle separation $\langle
+ x\rangle=1.5$ and a resolution $\eta=1.2348$. A Gaussian kernel
+ with the same smoothing length is shown for comparison.}
+\label{fig:sph:kernel_derivatives}
+\end{figure}
+
+
+\section{Kernel Derivatives}
+
+The derivatives of the kernel function have relatively simple
+expressions and are shown on Fig.~\ref{fig:sph:kernel_derivatives}.
+
+\begin{eqnarray*}
+ \vec\nabla_x W(\vec{x},h) &=& \frac{1}{h^4}f'\left(\frac{|\vec{x}|}{h}\right) \frac{\vec{x}}{|\vec{x}|} \\
+ \frac{\partial W(\vec{x},h)}{\partial h} &=&- \frac{1}{h^4}\left[3f\left(\frac{|\vec{x}|}{h}\right) +
+\frac{|\vec{x}|}{h}f'\left(\frac{|\vec{x}|}{h}\right)\right]
+\end{eqnarray*}
+
+Note that for all the kernels listed above, $f'(0) = 0$.
+
+\end{document}
diff --git a/theory/kernel/kernels.py b/theory/kernel/kernels.py
index d7bdbe2bf9ba49a30f4c8a2ae136c4843ce5c2cf..184379e5eafbcd12a1a47560ee88e02066da3942 100644
--- a/theory/kernel/kernels.py
+++ b/theory/kernel/kernels.py
@@ -11,24 +11,24 @@ from matplotlib.font_manager import FontProperties
import numpy
params = {
- 'axes.labelsize': 8,
+ 'axes.labelsize': 10,
'axes.titlesize': 8,
- 'font.size': 8,
+ 'font.size': 10,
'legend.fontsize': 9,
- 'xtick.labelsize': 8,
- 'ytick.labelsize': 8,
+ 'xtick.labelsize': 10,
+ 'ytick.labelsize': 10,
'xtick.major.pad': 2.5,
'ytick.major.pad': 2.5,
'text.usetex': True,
-'figure.figsize' : (3.15,3.15),
-'figure.subplot.left' : 0.12,
+'figure.figsize' : (4.15,4.15),
+'figure.subplot.left' : 0.14,
'figure.subplot.right' : 0.99 ,
'figure.subplot.bottom' : 0.08 ,
'figure.subplot.top' : 0.99 ,
'figure.subplot.wspace' : 0. ,
'figure.subplot.hspace' : 0. ,
'lines.markersize' : 6,
-'lines.linewidth' : 2.,
+'lines.linewidth' : 1.5,
'text.latex.unicode': True
}
rcParams.update(params)
@@ -36,147 +36,277 @@ rc('font',**{'family':'sans-serif','sans-serif':['Times']})
#Parameters
-eta = 1.2349
-h = 2.1
+eta = 1.2348422195325 # Resolution (Gives 48 neighbours for a cubic spline kernel)
+dx = 1.5#4 #2.7162 # Mean inter-particle separation
#Constants
PI = math.pi
-#Cubic Spline
-cubic_kernel_degree = 3
-cubic_kernel_ivals = 2
-cubic_kernel_gamma = 2.
-cubic_kernel_ngb = 4.0 / 3.0 * PI * eta**3 * 6.0858
-cubic_kernel_coeffs = array([[3./(4.*PI) , -3./(2.*PI), 0., 1./PI],
- [-1./(4.*PI), 3./(2.*PI), -3./PI, 2./PI],
- [0., 0., 0., 0.]])
-def cubic_W(x):
- if size(x) == 1:
- x = array([x])
- ind = (minimum(x, cubic_kernel_ivals)).astype(int)
- coeffs = cubic_kernel_coeffs[ind,:]
- w = coeffs[:,0] * x + coeffs[:,1]
- for k in range(2, cubic_kernel_degree+1):
- w = x * w + coeffs[:,k]
- return w
-
-
-#Quartic Spline
-quartic_kernel_degree = 4
-quartic_kernel_ivals = 3
-quartic_kernel_gamma = 2.5
-quartic_kernel_ngb = 4.0 / 3.0 * PI * eta**3 * 8.2293
-quartic_kernel_coeffs = array([[3./(10.*PI) , 0., -3./(4.*PI) , 0. , 23./(32.*PI)],
- [-1./(5.*PI) , 1./PI , -3./(2.*PI) ,1./(4.*PI) , 11./(16.*PI)],
- [1./(20.*PI) , -1./(2.*PI) , 15./(8.*PI) , -25./(8.*PI), 125./(64.*PI)],
- [ 0. , 0., 0., 0., 0.]])
-def quartic_W(x):
- if size(x) == 1:
- x = array([x])
- ind = (minimum(x+0.5, quartic_kernel_ivals)).astype(int)
- coeffs = quartic_kernel_coeffs[ind,:]
- w = coeffs[:,0] * x + coeffs[:,1]
- for k in range(2, quartic_kernel_degree+1):
- w = x * w + coeffs[:,k]
- return w
-
-
-# Wendland kernel
-wendland2_kernel_degree = 5
-wendland2_kernel_ivals = 1
-wendland2_kernel_gamma = 2
-wendland2_kernel_ngb = 4.0 / 3.0 * PI * eta**3 * 7.261825
-wendland2_kernel_coeffs = 3.342253804929802 * array([[1./8, -30./32, 80./32, -80./32., 0., 1.],
- [ 0. , 0., 0., 0., 0., 0.]]) / 8.
-
-print wendland2_kernel_coeffs
-def wendland2_W(x):
- if size(x) == 1:
- x = array([x])
- ind = (minimum(0.5*x, wendland2_kernel_ivals)).astype(int)
- coeffs = wendland2_kernel_coeffs[ind,:]
- w = coeffs[:,0] * x + coeffs[:,1]
- for k in range(2, wendland2_kernel_degree+1):
- w = x * w + coeffs[:,k]
- return w
-
-#def wendland2_W(x):
-# if size(x) == 1:
-# x = array([x])
-# x /= 1.936492
-# x[x>1.] = 1.
-# oneminusu = 1.-x
-# oneminusu4 = oneminusu * oneminusu * oneminusu * oneminusu
-# return 3.342253804929802 * oneminusu4 * (1. + 4.*x) / h**3
-
-
-#Find H
-r = arange(0, 3.5*h, 1./1000.)
-xi = r/h
-cubic_Ws = cubic_W(xi)
-quartic_Ws = quartic_W(xi)
-wendland2_Ws = wendland2_W(xi)
-for j in range(size(r)):
- if cubic_Ws[j] == 0:
- cubic_H = r[j]
- break
-for j in range(size(r)):
- if quartic_Ws[j] == 0:
- quartic_H = r[j]
- break
-for j in range(size(r)):
- if wendland2_Ws[j] == 0:
- wendland2_H = r[j]
- break
-
-
-print "H=", cubic_H
-print "H=", quartic_H
-print "H=", wendland2_H
-
-
-# Compute sigma -----------------------------------------
-cubic_norm = 4.*PI*integrate.quad(lambda x: x**2*cubic_W(x), 0, cubic_H)[0]
-quartic_norm = 4.*PI*integrate.quad(lambda x: x**2*quartic_W(x), 0, quartic_H)[0]
-wendland2_norm = 4.*PI*integrate.quad(lambda x: x**2*wendland2_W(x), 0, wendland2_H)[0]
-
-print cubic_norm
-print quartic_norm
-print wendland2_norm
-
-
-# Plot kernels ------------------------------------------
-r = arange(0, 3.5*h, 1./100.)
-xi = r/h
-
-cubic_Ws = cubic_W(xi)
-quartic_Ws = quartic_W(xi)
-wendland2_Ws = wendland2_W(xi)
-
+# Compute expected moothing length
+h = eta * dx
+
+# Get kernel support (Dehnen & Aly 2012, table 1) for 3D kernels
+H_cubic = 1.825742 * h
+H_quartic = 2.018932 * h
+H_quintic = 2.195775 * h
+H_WendlandC2 = 1.936492 * h
+H_WendlandC4 = 2.207940 * h
+H_WendlandC6 = 2.449490 * h
+
+# Get the number of neighbours within kernel support:
+N_H_cubic = 4./3. * PI * H_cubic**3 / (dx)**3
+N_H_quartic = 4./3. * PI * H_quartic**3 / (dx)**3
+N_H_quintic = 4./3. * PI * H_quintic**3 / (dx)**3
+N_H_WendlandC2 = 4./3. * PI * H_WendlandC2**3 / (dx)**3
+N_H_WendlandC4 = 4./3. * PI * H_WendlandC4**3 / (dx)**3
+N_H_WendlandC6 = 4./3. * PI * H_WendlandC6**3 / (dx)**3
+
+
+print "Smoothing length: h =", h, "Cubic spline kernel support size: H =", H_cubic, "Number of neighbours N_H =", N_H_cubic
+print "Smoothing length: h =", h, "Quartic spline kernel support size: H =", H_quartic, "Number of neighbours N_H =", N_H_quartic
+print "Smoothing length: h =", h, "Quintic spline kernel support size: H =", H_quintic, "Number of neighbours N_H =", N_H_quintic
+print "Smoothing length: h =", h, "Wendland C2 kernel support size: H =", H_WendlandC2, "Number of neighbours N_H =", N_H_WendlandC2
+print "Smoothing length: h =", h, "Wendland C4 kernel support size: H =", H_WendlandC4, "Number of neighbours N_H =", N_H_WendlandC4
+print "Smoothing length: h =", h, "Wendland C6 kernel support size: H =", H_WendlandC6, "Number of neighbours N_H =", N_H_WendlandC6
+
+# Get kernel constants (Dehen & Aly 2012, table 1) for 3D kernel
+C_cubic = 16. / PI
+C_quartic = 5**6 / (512 * PI)
+C_quintic = 3**7 / (40 * PI)
+C_WendlandC2 = 21. / (2 * PI)
+C_WendlandC4 = 495. / (32 * PI)
+C_WendlandC6 = 1365. / (64 * PI)
+
+# Get the reduced kernel definitions (Dehen & Aly 2012, table 1) for 3D kernel
+#def plus(u) : return maximum(0., u)
+def cubic_spline(r): return where(r > 1., 0., where(r < 0.5,
+ 3.*r**3 - 3.*r**2 + 0.5,
+ -r**3 + 3.*r**2 - 3.*r + 1.) )
+
+#return plus(1. - r)**3 - 4.*plus(1./2. - r)**3
+def quartic_spline(r): return where(r > 1., 0., where(r < 0.2,
+ 6.*r**4 - 2.4*r**2 + 46./125.,
+ where(r < 0.6,
+ -4.*r**4 + 8.*r**3 - (24./5.)*r**2 + (8./25.)*r + 44./125.,
+ 1.*r**4 - 4.*r**3 + 6.*r**2 - 4.*r + 1.)))
+
+#return plus(1. - r)**4 - 5.*plus(3./5. - r)**4 + 10.*plus(1./5. - r)**4
+def quintic_spline(r): return where(r > 1., 0., where(r < 1./3.,
+ -10.*r**5 + 10.*r**4 - (20./9.)*r**2 + (22./81.),
+ where(r < 2./3.,
+ 5.*r**5 - 15.*r**4 + (50./3.)*r**3 - (70./9.)*r**2 + (25./27.)*r + (17./81.),
+ -1.*r**5 + 5.*r**4 - 10.*r**3 + 10.*r**2 - 5.*r + 1.)))
+
+#return plus(1. - r)**5 - 6.*plus(2./3. - r)**5 + 15.*plus(1./3. - r)**5
+def wendlandC2(r): return where(r > 1., 0., 4.*r**5 - 15.*r**4 + 20.*r**3 - 10*r**2 + 1.)
+def wendlandC4(r): return where(r > 1., 0., (35./3.)*r**8 - 64.*r**7 + 140.*r**6 - (448./3.)*r**5 + 70.*r**4 - (28. /3.)*r**2 + 1.)
+def wendlandC6(r): return where(r > 1., 0., 32.*r**11 - 231.*r**10 + 704.*r**9 - 1155.*r**8 + 1056.*r**7 - 462.*r**6 + 66.*r**4 - 11.*r**2 + 1.)
+def Gaussian(r,h): return (1./(0.5*pi*h**2)**(3./2.)) * exp(- 2.*r**2 / (h**2))
+
+
+# Kernel definitions (3D)
+def W_cubic_spline(r): return C_cubic * cubic_spline(r / H_cubic) / H_cubic**3
+def W_quartic_spline(r): return C_quartic * quartic_spline(r / H_quartic) / H_quartic**3
+def W_quintic_spline(r): return C_quintic * quintic_spline(r / H_quintic) / H_quintic**3
+def W_WendlandC2(r): return C_WendlandC2 * wendlandC2(r / H_WendlandC2) / H_WendlandC2**3
+def W_WendlandC4(r): return C_WendlandC4 * wendlandC4(r / H_WendlandC4) / H_WendlandC4**3
+def W_WendlandC6(r): return C_WendlandC6 * wendlandC6(r / H_WendlandC6) / H_WendlandC6**3
+
+# PLOT STUFF
+figure()
+subplot(211)
+xx = linspace(0., 5*h, 1000)
+maxY = 1.2*Gaussian(0, h)
+
+# Plot the kernels
+plot(xx, Gaussian(xx, h), 'k-', linewidth=0.7, label="${\\rm %14s\\quad H=\\infty}$"%("Gaussian~~~~~~"))
+plot(xx, W_cubic_spline(xx), 'b-', label="${\\rm %14s\\quad H=%4.3f}$"%("Cubic~spline~~", H_cubic))
+plot(xx, W_quartic_spline(xx), 'c-', label="${\\rm %14s\\quad H=%4.3f}$"%("Quartic~spline", H_quartic))
+plot(xx, W_quintic_spline(xx), 'g-', label="${\\rm %14s\\quad H=%4.3f}$"%("Quintic~spline", H_quintic))
+plot(xx, W_WendlandC2(xx), 'r-', label="${\\rm %14s\\quad H=%4.3f}$"%("Wendland~C2~", H_WendlandC2))
+plot(xx, W_WendlandC4(xx), 'm-', label="${\\rm %14s\\quad H=%4.3f}$"%("Wendland~C4~", H_WendlandC4))
+plot(xx, W_WendlandC6(xx), 'y-', label="${\\rm %14s\\quad H=%4.3f}$"%("Wendland~C6~", H_WendlandC6))
+
+# Indicate the position of H
+arrow(H_cubic, 0.12*maxY , 0., -0.12*maxY*0.9, fc='b', ec='b', length_includes_head=True, head_width=0.03, head_length=0.12*maxY*0.3)
+arrow(H_quartic, 0.12*maxY , 0., -0.12*maxY*0.9, fc='c', ec='c', length_includes_head=True, head_width=0.03, head_length=0.12*maxY*0.3)
+arrow(H_quintic, 0.12*maxY , 0., -0.12*maxY*0.9, fc='g', ec='g', length_includes_head=True, head_width=0.03, head_length=0.12*maxY*0.3)
+arrow(H_WendlandC2, 0.12*maxY , 0., -0.12*maxY*0.9, fc='r', ec='r', length_includes_head=True, head_width=0.03, head_length=0.12*maxY*0.3)
+arrow(H_WendlandC4, 0.12*maxY , 0., -0.12*maxY*0.9, fc='m', ec='m', length_includes_head=True, head_width=0.03, head_length=0.12*maxY*0.3)
+arrow(H_WendlandC6, 0.12*maxY , 0., -0.12*maxY*0.9, fc='y', ec='y', length_includes_head=True, head_width=0.03, head_length=0.12*maxY*0.3)
+
+# Show h
+plot([h, h], [0., maxY], 'k:', linewidth=0.5)
+text(h, maxY*0.35, "$h\\equiv\\eta\\langle x\\rangle = %.4f$"%h, rotation=90, backgroundcolor='w', ha='center', va='bottom')
+
+# Show <x>
+plot([dx, dx], [0., maxY], 'k:', linewidth=0.5)
+text(dx, maxY*0.35, "$\\langle x\\rangle = %.1f$"%dx, rotation=90, backgroundcolor='w', ha='center', va='bottom')
+
+xlim(0., 2.5*h)
+ylim(0., maxY)
+gca().xaxis.set_ticklabels([])
+ylabel("$W(r,h)$", labelpad=1.5)
+legend(loc="upper right", handlelength=1.2, handletextpad=0.2)
+
+
+# Same but now in log space
+subplot(212, yscale="log")
+plot(xx, Gaussian(xx, h), 'k-', linewidth=0.7, label="${\\rm Gaussian}$")
+plot(xx, W_cubic_spline(xx), 'b-', label="${\\rm Cubic~spline}$")
+plot(xx, W_quartic_spline(xx), 'c-', label="${\\rm Quartic~spline}$")
+plot(xx, W_quintic_spline(xx), 'g-', label="${\\rm Quintic~spline}$")
+plot(xx, W_WendlandC2(xx), 'r-', label="${\\rm Wendland~C2}$")
+plot(xx, W_WendlandC4(xx), 'm-', label="${\\rm Wendland~C4}$")
+plot(xx, W_WendlandC6(xx), 'y-', label="${\\rm Wendland~C6}$")
+
+# Show h
+plot([h, h], [0., 1.], 'k:', linewidth=0.5)
+
+# Show <x>
+plot([dx, dx], [0., 1.], 'k:', linewidth=0.5)
+
+
+# Show plot properties
+text(h/5., 1e-3, "$\\langle x \\rangle = %3.1f$"%(dx), va="top", backgroundcolor='w')
+text(h/5.+0.06, 3e-4, "$\\eta = %5.4f$"%(eta), va="top", backgroundcolor='w')
+
+# Show number of neighbours
+text(1.9*h, 2e-1/2.9**0, "$N_{\\rm ngb}=\\infty$", fontsize=10)
+text(1.9*h, 2e-1/2.9**1, "$N_{\\rm ngb}=%3.1f$"%(N_H_cubic), color='b', fontsize=9)
+text(1.9*h, 2e-1/2.9**2, "$N_{\\rm ngb}=%3.1f$"%(N_H_quartic), color='c', fontsize=9)
+text(1.9*h, 2e-1/2.9**3, "$N_{\\rm ngb}=%3.1f$"%(N_H_quintic), color='g', fontsize=9)
+text(1.9*h, 2e-1/2.9**4, "$N_{\\rm ngb}=%3.1f$"%(N_H_WendlandC2), color='r', fontsize=9)
+text(1.9*h, 2e-1/2.9**5, "$N_{\\rm ngb}=%3.1f$"%(N_H_WendlandC4), color='m', fontsize=9)
+text(1.9*h, 2e-1/2.9**6, "$N_{\\rm ngb}=%3.0f$"%(N_H_WendlandC6), color='y', fontsize=9)
+
+xlim(0., 2.5*h)
+ylim(1e-5, 0.7)
+xlabel("$r$", labelpad=0)
+ylabel("$W(r,h)$", labelpad=0.5)
+
+savefig("kernels.pdf")
+
+
+
+
+################################
+# Now, let's work on derivatives
+################################
+
+# Get the derivative of the reduced kernel definitions for 3D kernels
+def d_cubic_spline(r): return where(r > 1., 0., where(r < 0.5,
+ 9.*r**2 - 6.*r,
+ -3.*r**2 + 6.*r - 3.) )
+
+def d_quartic_spline(r): return where(r > 1., 0., where(r < 0.2,
+ 24.*r**3 - 4.8*r,
+ where(r < 0.6,
+ -16.*r**3 + 24.*r**2 - (48./5.)*r + (8./25.),
+ 4.*r**3 - 12.*r**2 + 12.*r - 4.)))
+
+def d_quintic_spline(r): return where(r > 1., 0., where(r < 1./3.,
+ -50.*r**4 + 40.*r**3 - (40./9.)*r,
+ where(r < 2./3.,
+ 25.*r**4 - 60.*r**3 + 50.*r**2 - (140./9.)*r + (25./27.),
+ -5.*r**4 + 20.*r**3 - 30.*r**2 + 20.*r - 5.)))
+
+def d_wendlandC2(r): return where(r > 1., 0., 20.*r**4 - 60.*r**3 + 60.*r**2 - 20.*r)
+def d_wendlandC4(r): return where(r > 1., 0., 93.3333*r**7 - 448.*r**6 + 840.*r**5 - 746.667*r**4 + 280.*r**3 - 18.6667*r)
+def d_wendlandC6(r): return where(r > 1., 0., 352.*r**10 - 2310.*r**9 + 6336.*r**8 - 9240.*r**7 + 7392.*r**6 - 2772.*r**5 + 264.*r**3 - 22.*r)
+def d_Gaussian(r,h): return (-8.*sqrt(2.)/(PI**(3./2.) * h**5)) * r * exp(- 2.*r**2 / (h**2))
+
+# Get the second derivative of the reduced kernel definitions for 3D kernels
+def d2_cubic_spline(r): return where(r > 1., 0., where(r < 0.5,
+ 18.*r - 6.,
+ -6.*r + 6.) )
+
+def d2_quartic_spline(r): return where(r > 1., 0., where(r < 0.2,
+ 72.*r**2 - 4.8,
+ where(r < 0.6,
+ -48.*r**2 + 48.*r - (48./5.),
+ 12.*r**2 - 24.*r + 12.)))
+
+def d2_quintic_spline(r): return where(r > 1., 0., where(r < 1./3.,
+ -200.*r**3 + 120.*r**2 - (40./9.),
+ where(r < 2./3.,
+ 100.*r**3 - 180.*r**2 + 100.*r - (140./9.),
+ -20.*r**3 + 60.*r**2 - 60.*r + 20)))
+def d2_wendlandC2(r): return where(r > 1., 0., 80.*r**3 - 180.*r**2 + 120.*r - 20.)
+def d2_wendlandC4(r): return where(r > 1., 0., 653.3333*r**6 - 2688.*r**5 + 4200.*r**4 - 2986.667*r**3 + 840.*r**2 - 18.6667)
+def d2_wendlandC6(r): return where(r > 1., 0., 3520.*r**9 - 20790.*r**8 + 50688.*r**7 - 64680.*r**6 + 44352.*r**5 - 13860.*r**4 + 792.*r**2 - 22)
+def d2_Gaussian(r,h): return (32*sqrt(2)/(PI**(3./2.)*h**7)) * r**2 * exp(-2.*r**2 / (h**2)) - 8.*sqrt(2.)/(PI**(3./2.) * h**5) * exp(- 2.*r**2 / (h**2))
+
+
+# Derivative of kernel definitions (3D)
+def dW_cubic_spline(r): return C_cubic * d_cubic_spline(r / H_cubic) / H_cubic**4
+def dW_quartic_spline(r): return C_quartic * d_quartic_spline(r / H_quartic) / H_quartic**4
+def dW_quintic_spline(r): return C_quintic * d_quintic_spline(r / H_quintic) / H_quintic**4
+def dW_WendlandC2(r): return C_WendlandC2 * d_wendlandC2(r / H_WendlandC2) / H_WendlandC2**4
+def dW_WendlandC4(r): return C_WendlandC4 * d_wendlandC4(r / H_WendlandC4) / H_WendlandC4**4
+def dW_WendlandC6(r): return C_WendlandC6 * d_wendlandC6(r / H_WendlandC6) / H_WendlandC6**4
+
+# Second derivative of kernel definitions (3D)
+def d2W_cubic_spline(r): return C_cubic * d2_cubic_spline(r / H_cubic) / H_cubic**5
+def d2W_quartic_spline(r): return C_quartic * d2_quartic_spline(r / H_quartic) / H_quartic**5
+def d2W_quintic_spline(r): return C_quintic * d2_quintic_spline(r / H_quintic) / H_quintic**5
+def d2W_WendlandC2(r): return C_WendlandC2 * d2_wendlandC2(r / H_WendlandC2) / H_WendlandC2**5
+def d2W_WendlandC4(r): return C_WendlandC4 * d2_wendlandC4(r / H_WendlandC4) / H_WendlandC4**5
+def d2W_WendlandC6(r): return C_WendlandC6 * d2_wendlandC6(r / H_WendlandC6) / H_WendlandC6**5
figure()
+subplot(211)
+
+plot([0, 2.5*h], [0., 0.], 'k--', linewidth=0.7)
+plot(xx, d_Gaussian(xx, h), 'k-', linewidth=0.7, label="${\\rm Gaussian}$")
+plot(xx, dW_cubic_spline(xx), 'b-', label="${\\rm Cubic~spline}$")
+plot(xx, dW_quartic_spline(xx), 'c-', label="${\\rm Quartic~spline}$")
+plot(xx, dW_quintic_spline(xx), 'g-', label="${\\rm Quintic~spline}$")
+plot(xx, dW_WendlandC2(xx), 'r-', label="${\\rm Wendland~C2}$")
+plot(xx, dW_WendlandC4(xx), 'm-', label="${\\rm Wendland~C4}$")
+plot(xx, dW_WendlandC6(xx), 'y-', label="${\\rm Wendland~C6}$")
+
+maxY = d_Gaussian(h/2, h)
+
+# Show h
+plot([h, h], [2*maxY, 0.1], 'k:', linewidth=0.5)
+
+# Show <x>
+plot([dx, dx], [2*maxY, 0.1], 'k:', linewidth=0.5)
+
+
+xlim(0., 2.5*h)
+gca().xaxis.set_ticklabels([])
+ylim(1.2*maxY, -0.1*maxY)
+xlabel("$r$", labelpad=0)
+ylabel("$\\partial W(r,h)/\\partial r$", labelpad=0.5)
+legend(loc="lower right")
-text(h-0.1, cubic_Ws[0]/20., "h", ha="right",va="center")
-arrow(h, cubic_Ws[0]/10., 0., -cubic_Ws[0]/10., fc='k', ec='k', length_includes_head=True, head_length=cubic_Ws[0]/30., head_width=0.1)
-plot(r,cubic_Ws, 'b-' ,label="Cubic")
-plot(r, quartic_Ws, 'r-', label="Quartic")
-plot(r, wendland2_Ws, 'g-', label="Wendland C2")
+subplot(212)
-text(cubic_H-0.1, cubic_Ws[0]/20., "H", ha="right",va="center", color='b')
-arrow(cubic_H, cubic_Ws[0]/10., 0., -cubic_Ws[0]/10., fc='b', ec='b', length_includes_head=True, head_length=cubic_Ws[0]/30., head_width=0.1)
+maxY = d2_Gaussian(h,h)
+plot([h, h], [-4*maxY, 1.4*maxY], 'k:', linewidth=0.5)
+text(h, -3.*maxY, "$h\\equiv\\eta\\langle x\\rangle = %.4f$"%h, rotation=90, backgroundcolor='w', ha='center', va='bottom')
-text(quartic_H-0.1, cubic_Ws[0]/20., "H", ha="right",va="center", color='r')
-arrow(quartic_H, cubic_Ws[0]/10., 0., -cubic_Ws[0]/10., fc='r', ec='r', length_includes_head=True, head_length=quartic_Ws[0]/30., head_width=0.1)
+plot([dx, dx], [-4*maxY, 1.4*maxY], 'k:', linewidth=0.5)
+text(dx, -3.*maxY, "$\\langle x\\rangle = %.1f$"%dx, rotation=90, backgroundcolor='w', ha='center', va='bottom')
-text(wendland2_H-0.1, cubic_Ws[0]/20., "H", ha="right",va="center", color='r')
-arrow(wendland2_H, cubic_Ws[0]/10., 0., -cubic_Ws[0]/10., fc='g', ec='g', length_includes_head=True, head_length=wendland2_Ws[0]/30., head_width=0.1)
+plot([0, 2.5*h], [0., 0.], 'k--', linewidth=0.7)
+plot(xx, d2_Gaussian(xx, h), 'k-', linewidth=0.7, label="${\\rm Gaussian}$")
+plot(xx, d2W_cubic_spline(xx), 'b-', label="${\\rm Cubic~spline}$")
+plot(xx, d2W_quartic_spline(xx), 'c-', label="${\\rm Quartic~spline}$")
+plot(xx, d2W_quintic_spline(xx), 'g-', label="${\\rm Quintic~spline}$")
+plot(xx, d2W_WendlandC2(xx), 'r-', label="${\\rm Wendland~C2}$")
+plot(xx, d2W_WendlandC4(xx), 'm-', label="${\\rm Wendland~C4}$")
+plot(xx, d2W_WendlandC6(xx), 'y-', label="${\\rm Wendland~C6}$")
+xlim(0., 2.5*h)
+ylim(-3.2*maxY, 1.4*maxY)
+xlabel("$r$", labelpad=0)
+ylabel("$\\partial^2 W(r,h)/\\partial r^2$", labelpad=0.5)
-xlabel("r", labelpad=0)
-ylabel("W(r,h)", labelpad=0)
-legend(loc="upper right")
-savefig("kernel.pdf")
+savefig("kernel_derivatives.pdf")
diff --git a/theory/kernel/spline_3.nb b/theory/kernel/spline_3.nb
deleted file mode 100644
index d59c7f43846fd6217c2b98e193e410f6b6268cc5..0000000000000000000000000000000000000000
--- a/theory/kernel/spline_3.nb
+++ /dev/null
@@ -1,871 +0,0 @@
-(* Content-type: application/vnd.wolfram.mathematica *)
-
-(*** Wolfram Notebook File ***)
-(* http://www.wolfram.com/nb *)
-
-(* CreatedBy='Mathematica 8.0' *)
-
-(*CacheID: 234*)
-(* Internal cache information:
-NotebookFileLineBreakTest
-NotebookFileLineBreakTest
-NotebookDataPosition[ 157, 7]
-NotebookDataLength[ 36970, 862]
-NotebookOptionsPosition[ 35595, 809]
-NotebookOutlinePosition[ 35934, 824]
-CellTagsIndexPosition[ 35891, 821]
-WindowFrame->Normal*)
-
-(* Beginning of Notebook Content *)
-Notebook[{
-
-Cell[CellGroupData[{
-Cell[BoxData[
- RowBox[{"\[IndentingNewLine]",
- RowBox[{
- RowBox[{
- RowBox[{"f", "[", "q_", "]"}], ":=",
- RowBox[{
- RowBox[{"1", "/", "Pi"}], "*",
- RowBox[{"If", "[",
- RowBox[{
- RowBox[{"q", ">", "2"}], ",", "0", ",",
- RowBox[{"If", "[",
- RowBox[{
- RowBox[{"q", ">", "1"}], ",",
- RowBox[{
- RowBox[{"1", "/", "4"}], "*",
- RowBox[{
- RowBox[{"(",
- RowBox[{"2", "-", "q"}], ")"}], "^", "3"}]}], ",",
- RowBox[{
- RowBox[{
- RowBox[{"1", "/", "4"}], "*",
- RowBox[{
- RowBox[{"(",
- RowBox[{"2", "-", "q"}], ")"}], "^", "3"}]}], "-",
- RowBox[{
- RowBox[{"(",
- RowBox[{"1", "-", "q"}], ")"}], "^", "3"}]}]}], "]"}]}],
- "]"}]}]}], "\[IndentingNewLine]",
- RowBox[{
- RowBox[{"W", "[",
- RowBox[{"r_", ",", "h_"}], "]"}], "=",
- RowBox[{
- RowBox[{"1", "/",
- RowBox[{"h", "^", "3"}]}], " ", "*",
- RowBox[{"f", "[",
- RowBox[{"r", "/", "h"}], "]"}]}]}]}]}]], "Input",
- CellChangeTimes->{{3.560154174311659*^9, 3.5601543108245993`*^9}}],
-
-Cell[BoxData[
- FractionBox[
- RowBox[{"If", "[",
- RowBox[{
- RowBox[{
- FractionBox["r", "h"], ">", "2"}], ",", "0", ",",
- RowBox[{"If", "[",
- RowBox[{
- RowBox[{
- FractionBox["r", "h"], ">", "1"}], ",",
- RowBox[{
- FractionBox["1", "4"], " ",
- SuperscriptBox[
- RowBox[{"(",
- RowBox[{"2", "-",
- FractionBox["r", "h"]}], ")"}], "3"]}], ",",
- RowBox[{
- RowBox[{
- FractionBox["1", "4"], " ",
- SuperscriptBox[
- RowBox[{"(",
- RowBox[{"2", "-",
- FractionBox["r", "h"]}], ")"}], "3"]}], "-",
- SuperscriptBox[
- RowBox[{"(",
- RowBox[{"1", "-",
- FractionBox["r", "h"]}], ")"}], "3"]}]}], "]"}]}], "]"}],
- RowBox[{
- SuperscriptBox["h", "3"], " ", "\[Pi]"}]]], "Output",
- CellChangeTimes->{{3.560154211258333*^9, 3.560154216293594*^9}, {
- 3.560154312540955*^9, 3.560154319804675*^9}}]
-}, Open ]],
-
-Cell[CellGroupData[{
-
-Cell[BoxData[
- RowBox[{"Plot", "[",
- RowBox[{
- RowBox[{"W", "[",
- RowBox[{"r", ",", "1"}], "]"}], ",",
- RowBox[{"{",
- RowBox[{"r", ",", "0", ",", "2.5"}], "}"}]}], "]"}]], "Input",
- CellChangeTimes->{{3.560154325145775*^9, 3.560154343883732*^9}, {
- 3.560154674704236*^9, 3.56015467532159*^9}}],
-
-Cell[BoxData[
- GraphicsBox[{{}, {},
- {Hue[0.67, 0.6, 0.6], LineBox[CompressedData["
-1:eJxF1nkwVW/8B3BLqWyFFspW2UopWZLqvLWoLGVLi2RfE0pJlvqKQvaQpULZ
-okUIlSi7EMW9IS3Knsi+XPe653d/M7/53fPPmdc8M+d5f+Z5PzNnvY27sT0X
-BweHPScHx/++k6/b8CjbuRCuH+pdt1tSUU2R1j8prYN1M1e0trFM1r4LlJY+
-jXKzZW8UWRaMKc7j7bBB/sUdHbIsJ5dHSC2Wdsb24vQUCZa3jNpFko7u8EsX
-FBBh+YdgCpPR7oGRQhkZHpZlspapz9M9UWWf0DtjQcX5vZ6us1LeOKwsatDH
-8kvqr4ypA354F3PKroXleRf97+OO/+HqYl7Ztyw/5Px+e1b7Bl4t+nknjeWD
-TXdkZ9oDUMkX9DSY5cHEwxVTzjfx5+5Lj3Msh9stmE/Sb8FVUnRQl2UK3Slu
-TOo2Ctbo/1rE8pU6ye2j+aFY9t91u59nqVgbS20cORCO7y3Pk4pYttmsxf3X
-MQoy8vViZ1nmovt8mvKLhoK0qMNWlqm/TCWmtO/g8/gmB7o5FVef8b2ZaI+B
-OL/5zWiWdWP6eSYexuK4Xlf6SZbFr1YcH3eOg7HhX29xlssPeI2N0u/ieaKt
-fuoZKpZ865YfkUpCVwqRctuMio7yMs/hwSQIhFy8RbD8JCux6m/+PWgIvVEY
-P03FMY9jlkMHHqBqTC/LgOX4ZSXxA46pcLAaih4/SUWhl67HJdGHiKxRPB/H
-8oOhnWKTfg9BZggsU2PZ7bOw44T2I3x+fzXU/QQVwg8+cI23p4Hvt+n7+uNU
-0AWLcy7sScfBuEVvzVjuvZFuOPYwHVwPuq79MaGi2PF6yqhzBrbE/zrDwbKZ
-iqrmP3omkO8VLmJExaOG1It/pXJgaX3ozGd9KqYSizZ3hedgHZ9lsTLLhxwb
-e1ppOVAxWRIQrUfFMPesaQnlCUKORbkd0qVCY+8xzdvBz1DSrVEUe5iKzy8Y
-XPKjL5CwtPVd4T5WXiyzEFDKw0vfvvARLSp2y8VpBtvl4UXTk68bWY6cypn0
-acmDdQplIoSgQvXOF3ubp/kI7jm9R203Ff81bNFVtnyJrTdFJ1arUpGd/0ru
-6d2XYH7QXaWgQkVr4n5u2Y+s9a1GNeo7WP11PFUquqsQ1B6bT/rbqWjgvqnE
-IVIEKbppxpktVKze+034U20xmKfM1iyWoSL3xe1v55VKQOHrMIwQomJIeChc
-17oE2985ljusoEL+ii6hEFeC6+srZYnlrH7v4XvUQyuB/9lX1H5+KmLqwx3M
-at7ixdnadqmlrL72RI0fNi9DqXC4iDKTglS+xkAh7XJcHrU592WQAmKvr7m/
-RTmkg+dm3QYo+O6mqDbqVY5vvWdiePopEG0N6//4pBxbTc+Pbeuh4E6ivk7I
-igq8eNtm7/KDgkDZZkHyWwW2qQ603WqhwBEtScMeVVh1Kb3z6GsKONJmSr6F
-VUF2Ii6zoJiCxEXi3xsyqpDuRslaVUTBh3oHyZy2KpAHR4615VOgcJyRZq9Z
-DfJJxJH9TykYdJJ/9pOzBlI+et41yRSci7n27nN0LTRW8Uddu0GBW9+mnsLc
-euyuULjeeYCCXxZfl1bU14PvnDBP1j4KTDpClJp661HoG5DlDgo0Gge8+9Y2
-wCrKRXVBkwLu/MwVq4Mb0Hkh6QHnDlY+v/WEl2UjMn9vtMqToqBCRCxp14om
-pGSOVzyltYJ3/U35Lo3P4FpuY6Wb3YrIFY2bH1lRIGfo5PeQ2QKbb8ffc39s
-A6+L1pTTthZcrLlp4y7fCd3FtWGvDD4j/wK/iajOD3zvaFEe72rGI1+Xoe0n
-f2GUNjNr/OsjAlIZRlPDv3F171fnqdEGcFB+aZX69KBgmmuF88IHyF12uHp8
-cx84pd7sdPpai97phBN8Zf2YLspPMPpaDd2+c9diLg5iv78ZbWlhJY71RitK
-Sg5hTajZO7PMcgxV2fh/eP4XFrF1a5/Hl0FJqE3igOUI+N4r83dVvwEvpq2/
-T/5DX/vFigD3YhweMZvdYDqGGef27YXSL9FRacdptmUcG3hEZDnU81AUUm0h
-PT2OGErPw+NhT+Gw2PNXc/MEGuqUHtq0ZcHWw7eI/94k5lYVBkhMp2Gbo2Kz
-vvUUBtLyt9xXT8Un4S18nrumodSuoqTSm4i3uQH3t3DPQKnzmPOy1ljw6Juu
-TeCZgWAkL29JWSyM/sgncfLOYOVKu6irObEY2NgU37ZiBs90BlWlbsRCJGlN
-jL/kDJhadHth5Vi4BD4Loe6agW2Uj3p4bAzWnWrz9Lswg92yp/Y/s7kDX87N
-hh9/zOA8b/EGeaNI9B/k65D7PQOevS3azdsiYRgybHmjdwZ+r2h9EYKRkFn+
-wl397wxGh1y9T3yMwEdx1aiHczPYI1E1xqcTAUkNovmyyCyWaEppcRwOR6Wb
-sb6Eziz+5sSc/nQ8FLzffY64Fc4iiH/ZzuXFQVD3yhGIejWL2Uu7ch0fBMFa
-uKP1Rcks6Ku9r1MCgvD6iNrZ8fJZ+K0d5u8zCIJD0ejFy02s79nyrRoduoWK
-SLv7Pv2zaLS94ZoocwtX9x37FyQ6Bxvln2ICWYHoy1ofl+o3h6eUSdW+Ln8M
-3/xN4/afg5aCuPbKen9M2DyydAqcw7P3d0dPFviDQ2q9onLoHMrC1P+J3/LH
-2gTpysqEOXC6HIh4oOgPw2Cpsb6COWx4an2COncdpY4S+lv+zIGRL3GSftoP
-Vdo/8qOH52Cxv2b0jLIfGjYmr5kenQNX1jOnnqV+aO8S7ymbmUNLdG6ew2tf
-TJwU9zFYRMOBxOyoPWK+UDiyLttDigaZyErNVd3eiFMQW/TGlIYGodlPs1Fe
-MLeQe9xzigbuoshIDzcvyMSp6Aqa07B+oldO4KgXCsmj0bY2NIybbsp+xOsF
-aluAuKA7DaLB8QcPBV+ByK1hVdsQGs7lKDy6EOSJmN/v7QXe0nC/1Lyh6/4l
-mK1pWqrxjoY2SXtrjeuXsOFo51ObChpWzvOP5FpdQsHrqfFXdTT470zInJa9
-hNaoTddtqDQY3Vk7OZDvASEiNv7VCA1vOrlKjT9dRPQ9hzpr6Xmsjw3tb5O6
-AK8rAZf8Ns7j/Podj0x4LsDCOEUqQW4ecl9mstq+umMLb5vXxy3zaOTNbGy7
-7I56b22FnbvmEXTb3Pb0czdwnZa5zW88jw6Df03HFFzhuaZb91XgPEq/GCbv
-0D4H88mFmZageVz/4ZJeIHoOBz+JpQ/fnsfmzvc89sPOEA42oq+PnoezTYuH
-VZwzcmfKn4Y/mIf7tZby3H4nDH5J5bcpmofeKYZadJwjzOLONvMPzOO3V15M
-ykp7iL/kXecyNI9gJwbf2JAdulpeOdaPzOOx8e77VyvsYLdcmDNoijVfW4Kz
-krsd3EJrVUhOOkzGJusUmmwR4L8taVycjlVCvg80Ym2Qc57Trs2YjswNnUo1
-h6xwPiw3T/UEHbTnBvkFG62g9OTMQsxpOnbq9r9p5bBCwUBRvIEVHeFjc/Em
-XpYotTlX/8GVjrmFP72Wjhb4fIqi9DaEtV9RQcXcaXPMHcykpb6jw695aeyW
-4FPwe5saWF1BR6i2cDjD8BQ4dtwT+FNNh+zJsUjxdaewRCpq/Y5GOv6GLKvR
-yzuJ1bSrOtXtdNj9bLhi9fMEVJ7rJw2O0XFUyFJom44p3FZNaShvZCAo3ndT
-5UljTIT9qzSVY8CwRfOGtKIxPLn+6PtsYiAy+YN3HdMIfqM/rKq2McC3aZvV
-vywjhNXXhZjuYWBH9f6gpHlDZF+73+FtykBb/Ms8GwsD9PTtv1oZzICq5pBa
-dIIu9qc2p20OY+Cqos7fYB1dPDxl1hQTyYDAWt6fogwdWDRe3GB7l4Hz1ACp
-RbY66Mx7+JE7nYFsHp6I0Z1H0OrLlNZ+x0Dow8WDdxe0oawWqve8goEEQzKz
-87U2ov+turKqhoG8uJTtLy5r45j11sa+RgZK3qUYLP53EA2HzD2DvjIg0X1m
-m+nAAVQKldTXTTFwJ9NXj5zch/WN2lPb5hj474YOz/yrffC/2SKZSGdAL4RD
-wtdvH4jZwUtOXAtQYsxzCi7Zh5LvaySXrVjAqpp2GYcNWih47Omhq7iAF2q1
-h19170HaXuV1TdYLyE7PdL7coY7KY7yBdLsFKE/Qv0bfUke3Zc/QJqcFxPLm
-l15RUcfGgLslQW4LeD5uxXn/jhoyammntXwXsGsp7eHmk6rINKhMfBm3ALWG
-VVxRPDtQa3WfozthAe3rlY9+eq+M/ouXnVbcX8C6x4ZP672VIRcnp+H6iJU3
-tubD4Oh2ZHWEtsvlLoB8b3jApX4bsq1NVt+rWwCXf4mCfZciPnhsuf6hYQH6
-r+Vv73NSxGDg4v6ZpgWIiAwaXprYDIWsV0Um1AWM/y50+MS7GTlD60wFfi/g
-fUXFDxU9BTy51BvrP7+AxR7vw8ilsmi4WUbPXVhAZs7yRr0sGQzdjbf9wcHE
-psaFguaDMlB8fURFcwkT96/xNF24tRHPGM9aJ0WYKLzGbz62cgOe3/IUctzK
-BEcG95h5mySG95wNuLWdieK8g+/Fd0liy9TByXQVJtQLFN/MJEvgic3KL792
-MTHy/FrekKs4svEy0ewQEwYa2w+5yq3F4My9Zd46TMSmpHInJYlBPjfAJ0Gf
-idYtXycqBMWQJW5sTjVmgvy6Z4jBXIMM2pjUMUsm0jfVfMmfWYXevI7o8zZM
-HPW/WWLtvQobnco5w+yZqFHe0TG+sBJpbVE9dS5MdJcN7cvgX4nUl0qP911l
-ImTrJ8dh1t/7PdfzSmp3mBCeKAotmxZAp8zxVJM4JoYH9x55ryIAse+7V3gk
-MJFipRIn4MmPRF2+idxkJtwXhf304OTDXfknRQpPmCgNGHK7q7UU1J935A4/
-Z+JG97eW7YlLsDLeO8E+jwlThcCWnxM8iFmk451WzAT/vTW5pvmLEf17YI94
-FRN+XqaCvUbc+Jz46blmLROCrpMhOhVcWG74SvJ0PRNluVG5wSpciHwXxBH/
-iYmYIwOcF6U5EX5ftnr5dyZ2eErPnFZnEh+NBdSUupig9+8mxfYtEHy805n6
-3UwcKO31WWPEIEKvVgffHmSielquI9p/ngg2tdXjnmadr6tT5MC6WSLN4Lu5
-/RwTkRGOfzv3zhClOqZudXQmsgN2f3SznyYm9h6+E8ZFYvHNHXJKNZOEuZxi
-u8gKEprCQncX14wRXtIZg54iJKJ5XM7NCo0RMWsl5ttXkzjDdaXkb90/ok5w
-ucQDCRLSB9cZJm8dJnbMTtjIKJLo8rU9RXUZII6Ou1wOUiLxa/2viCvG/YTT
-395bg8okGOphyrz7+4jkrrbsZxokBH4J/Li4u4d4/dWgRHAPidP8xeP0g91E
-K+VD4wWQ0Pu6NbLnxG9i6YeSf6qHWHl5z4y5yXURGytVOBJ0SNjR+DJf+vwg
-iNJnQjR9EvzhKQ/Cfn4jLuWlqJaZkFBlPgluonQQPxMDfLStSVjUKBcGnKcQ
-tBh62GM7Ek3+5f26aq3EyojLycucSCRnFFWe4GshdG44lDe5kTi/QvH4xtIm
-ws63q2WbB4ne8f0ppwsbieuep3rueJIod372qK66nih01uU54Uei7+ecJl2x
-lvhkW7Xm9X8kjkJyUdKNamLo7J5NawNJpM+LjUiOVBJSxkr6P2+TiDwsdiuW
-u5zQ1H98ViuCROCVd82FZmXE8UPS7mnRJHJSOB5e+lJChGoKxzgkkMjYma/W
-t7uYyFQNS/9wj4SJW4bVEdlColxpUdHmFBJq7ift9yoXEJ0K12rDH5E4mczv
-9cAij5jeMNP+L4OEfSF/jO+p58QKCfc/htms/EKCEpPKT4jNawbnC56SeCKa
-Nm1w+DFhxdcp4VVAYo/K3U2+4mmEz2KTbV+LSKhn0tOPv0gl7pKNWrvfkCg2
-F63mCLtP5NEOGieXkggVDfU/U5dANEyW2ZLvSXyxNcuqsYwj+kbUPa2rWPvV
-xnA9NLhDcAy+CKqqJbE0ybs6XiaCWNutkCjbwOoHWSVjdCiEUPv+KCe4iUQr
-ZU35dnogsULJae7FZxKyysJHZj39ieH/lA53UEj4JkuN51f6EB9apu5ytJNY
-tvvx4gwXTyJj49tehU4SDmv3bS9rv0D4e95QMfpBQki9P+vvqDNhXnc4wPsX
-iUqNoPRwH1tCQ0yw5VEPifDZb+oPbMyJlS5UqYZ+Ekv8jB/vGzYhxkrvuU38
-YZ0PTVC77aku8VHQumztCKtvozL2U7VaRLaVPP+BMRIR2YpB+2xViJsFI2Yu
-kyTWxe7/6N28kbBaVJgTO8O6X7836N7MXEHsOeEz95bG6nM5c/eB2Mm9otla
-h3sZJLiHXIV3SVbtnffj1uPh4NDi+L9nMNPB4x4n223N9UlK3GwXSEcPnuBh
-+1y1eNBjfra/8qlXHBFju1Y1afD7WrYLzzKWXxRnO/pFpUWSFNtHTAwYf2TZ
-fn3PaWfYDrYTN91/3qTH9tTkoczAo2wbvJt4sMuAbR4T3fBMY7Y9/Wjnrp1m
-2/jTSYWtjmzze65MDw9g+7/yiES9N2wfXWea2SDD+f/e51Z3Rug226niNjIf
-p9mWqXocf+wMF3verZWbW9+x7TC5qKV+M/f/2yZ/Qvp2FNvnA0wYunfY9jQp
-bOePZTtk2jMyOp7t3F00ekIy27QKsi3rKdvRrfwR1XVs30t3dQqqZzv9cvOB
-I41sF6+Onm9sZvu7mYgT5Qvb8j1iB7p72N5e6COZ0ce25q1vNPsBtvXlk/P+
-DLHt4SItOT7Ott+eG7SCSbZvCXRTL0+znfgiPWxuju00/0WOJfNsPzWy3+/H
-YLtwQ60EwWT73aQcjSTZ/h8gYzo2
- "]]}},
- AspectRatio->NCache[GoldenRatio^(-1), 0.6180339887498948],
- Axes->True,
- AxesOrigin->{0, 0},
- PlotRange->{{0, 2.5}, {0., 0.31830988618378947`}},
- PlotRangeClipping->True,
- PlotRangePadding->{
- Scaled[0.02],
- Scaled[0.02]}]], "Output",
- CellChangeTimes->{3.5601543446066847`*^9, 3.5601546760449047`*^9}]
-}, Open ]],
-
-Cell[CellGroupData[{
-
-Cell[BoxData[
- RowBox[{"N", "[",
- RowBox[{"Expand", "[",
- RowBox[{
- RowBox[{"(",
- RowBox[{
- RowBox[{
- RowBox[{"1", "/", "4"}], "*",
- RowBox[{
- RowBox[{"(",
- RowBox[{"2", "-", "q"}], ")"}], "^", "3"}]}], "-",
- RowBox[{
- RowBox[{"(",
- RowBox[{"1", "-", "q"}], ")"}], "^", "3"}]}], ")"}], "/", "Pi"}],
- "]"}], "]"}]], "Input",
- CellChangeTimes->{{3.560154431542004*^9, 3.560154500452031*^9}}],
-
-Cell[BoxData[
- RowBox[{"0.3183098861837907`", "\[VeryThinSpace]", "-",
- RowBox[{"0.477464829275686`", " ",
- SuperscriptBox["q", "2"]}], "+",
- RowBox[{"0.238732414637843`", " ",
- SuperscriptBox["q", "3"]}]}]], "Output",
- CellChangeTimes->{{3.560154427870244*^9, 3.560154500989884*^9}}]
-}, Open ]],
-
-Cell[CellGroupData[{
-
-Cell[BoxData[
- RowBox[{"N", "[",
- RowBox[{"Expand", "[",
- RowBox[{
- RowBox[{"(",
- RowBox[{
- RowBox[{"1", "/", "4"}], "*",
- RowBox[{
- RowBox[{"(",
- RowBox[{"2", "-", "q"}], ")"}], "^", "3"}]}], ")"}], "/", "Pi"}],
- "]"}], "]"}]], "Input",
- CellChangeTimes->{{3.560154530785256*^9, 3.56015454752137*^9}}],
-
-Cell[BoxData[
- RowBox[{"0.6366197723675814`", "\[VeryThinSpace]", "-",
- RowBox[{"0.954929658551372`", " ", "q"}], "+",
- RowBox[{"0.477464829275686`", " ",
- SuperscriptBox["q", "2"]}], "-",
- RowBox[{"0.07957747154594767`", " ",
- SuperscriptBox["q", "3"]}]}]], "Output",
- CellChangeTimes->{{3.560154539254085*^9, 3.560154548437131*^9}}]
-}, Open ]],
-
-Cell[BoxData[{
- RowBox[{
- RowBox[{"DWr", "[",
- RowBox[{"r_", ",", "h_"}], "]"}], ":=",
- RowBox[{
- RowBox[{
- RowBox[{"Derivative", "[",
- RowBox[{"1", ",", "0"}], "]"}], "[", "W", "]"}], "[",
- RowBox[{"r", ",", "h"}], "]"}]}], "\[IndentingNewLine]",
- RowBox[{
- RowBox[{"DWh", "[",
- RowBox[{"r_", ",", "h_"}], "]"}], ":=",
- RowBox[{
- RowBox[{
- RowBox[{"Derivative", "[",
- RowBox[{"0", ",", "1"}], "]"}], "[", "W", "]"}], "[",
- RowBox[{"r", ",", "h"}], "]"}]}]}], "Input",
- CellChangeTimes->{{3.5601545811631327`*^9, 3.5601545907204247`*^9}, {
- 3.5601546570572557`*^9, 3.56015471264272*^9}, {3.5601550735178423`*^9,
- 3.560155113042481*^9}, {3.560155146451144*^9, 3.560155154786213*^9}, {
- 3.5601552200011473`*^9, 3.56015522178111*^9}}],
-
-Cell[CellGroupData[{
-
-Cell[BoxData[
- RowBox[{"Plot", "[",
- RowBox[{
- RowBox[{"DWr", "[",
- RowBox[{"r", ",", "1"}], "]"}], ",",
- RowBox[{"{",
- RowBox[{"r", ",", "0", ",", "2.5"}], "}"}]}], "]"}]], "Input",
- CellChangeTimes->{{3.5601551158258877`*^9, 3.560155135295669*^9}}],
-
-Cell[BoxData[
- GraphicsBox[{{}, {},
- {Hue[0.67, 0.6, 0.6], LineBox[CompressedData["
-1:eJxF12k4VV/UAHBD+hdCpaRkLOlWQhMaFpEpSaZSSIYylFKUKSFkyFgiIRJp
-FFFR3WXKTNx7K5m5gykqlPne97zP8z7vPV/O8/twzll7n73WXlvO4aKZMx8P
-D08oLw/P/97TAx0Wqzq5H7j3xPWAopIbuYoqa3xc1hAWv5m6y5rVR041+aas
-rDUY7NFaQZo6jullMTICsq5wTFA3jPPrHJ7f731hSsYXtvOHPWntu4bUOZe7
-v2UiwfaxQlTk51v4X0f/plGZ+/AOLkrzxyRjVv1DzxGZp3A9K0Z5bt0TfJUf
-2XFeuRSuVhjOu616h+eg9f7Py5VwpemYZpFXGXowN9OLXtVBlLhJZPLXKhSU
-C93Uo94C0y3Tx9Ml6jBWrIGUZU+F2QnV26WpjejQYYH8jd+AfvXY7JHBLyj4
-oO+C549vkBV4wGfv5BcsPOkh1cP6Bjkdg55iPC3I337Lt5T3O7yy/Vs+JNSC
-j9tKd3ju+Q57bk1MuCu0IOurbF73o++w6VLow1fHWtCt5Wd8iW8bpNyfFt/2
-vAU9P4c6XNzUDreNSL4iZq0ozdBb/G9HOwifW9z00qoVG/iWPgvQagfTEyPp
-hqdacaNWzJ9I63ZYVfLl7VWnVvxRmhT0OLodkq/wfc692orar3Mzfvxqh8RG
-UXO71FZc/qC2Xed9B5iK1O/a0tOK5JKowPqqDqitKCOdpLeie5ux/LHWDuA3
-LksLG2jFz6taXe2GO6Co88dU669W9I3/Me0j1Qmrpq5GafFQsC9sePWroE6o
-FbWW85alYMElYfM1hl3EfH5ctc+Ogo2pK31fWnVBwGl1+8MOFGRVrX140KkL
-PMWfR1qdpeBayc0j5290Qekfb197DwqGlh+6WV7cBc99sxoPBlLQasWNN27y
-3bANj2prp1Nw9s2flZ/muqFYdqnTxjYKindPa5gt7YHM664xrzsoqLyEx35g
-dQ9weMqT9vRQ0NFG5MVytR6Q9LfU0WAR8SzaonvWpQcK94mvEJqkYIalk5fo
-1x4I/Tt4fVCEiln+7sMqx3vh6f5fsVHaVFQpWzp9y7YXSo3vqtvpUrFsUZ5A
-j2MvJLqwDZT1qdgbw5CNudQLYdVyGZXGVJTOtD0+FNkL0etyJCuPUzG16mhV
-1sdeECRlWK66QMXEZTszVsj3Ab3Num9xEhXlzSjPXZX6QC8o0/FcMhUL7l0q
-KVPug/M7XxhX3adii8xLmsfePmBXrDC/kkFFETVFoQaLPjjgK6Wa+4SKUVZr
-fG7e6oM9RvJqp0uoGPJw/tjkzz7gD5zSl26noqnKe88tE33gVrl5ireTiK/8
-SoLDTB9kP12k2t9FxdL+4ZYWgX4IO1Gkm9JHxfGNP46+lO4HCa25kfZBKjq8
-KD5y1rQfHHZv0iz9R8zHfs8L6Vb98C2CfufsNBXZTVtjaDb9IB1g4CM6S4z3
-V3bjQdd+kDok3Ge1QEWKWuJhmZB+YNr7q33kp6F2iYdh25t+eIRTtF5RGooY
-kVxFS/tBQebAD8vlNOxsZ0bolfWDr8w2l5oVNPSZt6ktauyHO3wCipmraPga
-DusnMPvhrlmanNI6GspWbzpkJEGHMknXNYkbaThmRXcKXk8Hs3CNmi+KNPw4
-kBH6XoEOxXcKeQSVaHhi6eoqRRU6jMRnrPQi0TDeeJEOvyEdNq6yK5HeTkMe
-aq/WRz86lOv6L2reQ8PQIXKWXzAdwofWhAyr03ApTzqfRgQdgtu9Rvk0abhi
-m3VVcRIdlmWdi9m0j4Ybb7Xq57+mQ1qZ+jho0fBZen7ehXd04K2p7dqrTUPl
-opilW8l0eCAp8W7HQRqq9xk25DXQwS17TE1Cl4ZGeytMHrHo4GcYefSlPg2b
-j2Xm24/S4UaB0NMIAxqauwSKyUzSocTZ+PdpQxraJGlSHvAyoKFwiQHfYRp6
-/Cq0TJJigG3PeanVJjQcF0h4a67AAIGnffQWwtekLkqsIDGgiyx+/9ZRGgYZ
-bvkRu4cBql+KGwdNaXgnO9smwowBh0LoJ/zMabimNPiTnjUDYl8axYpb0DCt
-5bS0gD0DhDZeKnhGOHdhXW/wBQZkTP8tq7ekIUl8VkvLiwHFFlYvLaxomE9q
-y2L7MYAz73WrnfD743cd/SMYsFttdFH7cRru87hcpRHHAGnUyzM/QcPyUNON
-00kM2Jtip1FHuL5AeMArmwGi2+s25FnT0KR2WH/HMwbU/Fx3Y/lJGlK7a/P+
-vGZA4b+tNVcJdwqFuXmQGbD5V6TsrlM0tJd3bNj6mQEFXV9VYwkz1bW3jjQw
-4NKaXhU64THnhdFzPxgQdnn7dJANDa8EdJgo9jKgM8GrrI7wdGJJPoPFAPeF
-a1dFbGnIW+bteWaSAQolJc+jCYd/M6fIzDHg8NJpUhVhwVHVHd28THiaM5k8
-TTiOX+xu2hImaMTmjyvZ0VB87djkSVEm8FQq7bUknKLSaCm5mgmzOvae1wmv
-13/29rsUE0jLrZKyCD+yjZC4p8CEKJLgk3LCil5nfSxITJCND8jtIrxNs8f8
-oQoTepJOrJ0kvJPnxPbh3UzA6Cq+/07TcG91i+Cu/Uwo8yw7vZrwwduGrBs6
-TGDoG6vKEzY0qyivNyS+J+JyjUTYdM3e9FWmTIivWqm2nfDx7jc+9lZMkHA/
-ekaFsN3jrRbPbZhQwbtOQJmws1vO9n8OTIiLuLpeifB5FWkhbVcm+HHOPJMm
-fOXfPVb0RSZEnO16vpyw30fRim/eTHj/aVCOl3BwSES6XAAThARChUeJ+CMM
-eH3PhzAhcF+++1fCcSJ+Fu8imCDm5AalhO/RxrfzxTGhzr/41gPC6anuQkeS
-mJB1M17Xl/BjewYr+QETUv3nvM0JP1e0rejPYkKx4791JMKFP7+mb8tjwvje
-4N0LxP8pKTTx9XnFBPNFaZ8bCZf51FhUFjGh9aNeXQrh5kUlQtblTNjNPqys
-SPhrvepAdg0TpG7lxg0Q66Mz/lnFWBMT5HgSnHMID0ul+Ya2MyGxYtXF1YT/
-9ItbfullAkfo9aM6Yr1N58WorB1gQuShdiNfwgK7bgzkTzBhaWT3dBOxfoXn
-pipmZ4j1kPiefZnwivJLGYd4WLAsRt5XnLDsEQfLDmEWJB/LbTAm8mGf86HK
-xYos0JK6MJBI5JPOFnLGsa0s2P+kOmEdYaM/u/3S1FhgJ//6dSaRfyeuK6mq
-AQt4h9SMM4l89UoSemh7ggXd76ZmvIh89z8V6pdnxwJWZYtBnxkNQ+QWLCec
-WLCqTEPYiHD8yzHhSE8WUCI+/BY7RsMXn1v9iqJY0HcqUMWLqC+Mv8lWQp9Y
-YLiENzqeqFePtHdc6apkQWhNx5dCop6djmmOy69nwW4IHG7Ro2H7BoE68zYW
-iJxYm8VziIat5pf3pk2wwGd7/kplol6SC4xlt5EGgOR9WKScqL8BCwP7OCoD
-IPAxKz+TqM8ahjetW/cMgM0HVAzYTcOi3pJE70MD4EU3u6e4k6i/opsEyPYD
-cObHVz9Dor4nX+AdNkkegPU2JtY1G2hoUZK2WC5jANJIssedFIj5FlBXmHg8
-AIkFJaQFORrGpnnYJBcOwNigjrm8DLEfNHQ09zQNwOmFoqeakjS8tPntm0uL
-BkHx09Lz2cLEfDPdrid6DoJuynU9y99U1Luu8ELXZxDeqz6Y9B6j4sFVne3/
-AgfB5YHCw8SfVNQ8dET91O1B8O/xGkJivyXlbJ9QeDII10iS5sO9VBRynnQp
-7hyEdw0ufoEtVGxkXLdo0x8C170bJB+/oGJtwO6bUSZD4Nz5wEL3GRWrxH8V
-7LMcgrQjgbd7if7go669aJbDEJRl3vwjnE3FF48P1rteH4KOpfHxu4j+Isbp
-P625wiGwGciXrw6logkjfou09DBEfRFKZpwgvreFQbPZQFgoYP1uKyoevrIn
-8AFpGI75blUMNaeiAV93y5o9w6AxfIS8yoQYn+yWqytNh8FV7s5SUR0q7rL5
-XLYkZBiOMhWFT26lohRtxnKCOQxKlVfKj81TMHXdEY7qz2HYoj4ct2OGgpKO
-mXmXxofhX3dm7/J/FFw9rjc3yh6GvvLmJVW/KCgmdjdzUGIEpLu2erH7Kchv
-rDzSZTgCkq7jbxtqKThc6RBU+3IEZL8/+lmTQEGZhDa32KIR0FwW0P0zhoIW
-diaWFh9GQDF5smhZFAXJ0xqk3toRcM+tFNUJoWDituW0KfoIlBT+CAi4TMG9
-ybhp09qfYBPC3Ekyo+BtN6kvYeE/oWnBX8pThOg3l39br3N6FHScH9e98GtF
-VgwrJMl5FJw61hpJE/13uuDUwID7KDBVMnbf9mxF4UVrCm/7jIKAX4O5rUsr
-Dv+z1vuWMAr1oWUZNMtWzOns8nCtGgW7APFrEiqtKJXHxLjNY5Az4ag61teC
-gvD3TOfEGFTzZTX2aLag/ujJKXnL3+C8X3DHlHETtlU48Z7c+gdoAvpp9/bX
-Y3FElZ3s3z+wdeOM+bvIGjwr4N3b3DwOt86mW4hUV6HjZf9i4dQJqCvbEqRp
-XYHbz21pNj4zCXVROtfm15fhlxVbhbw1/gLVzfuAwfsP+OFVyIOt/P9gp6bW
-gXMt79Cfl2Ta2PUPMrfBYvfhNyjY6WfgUTQFbbOsNrsLr5GZK3f3YcA0ZJjk
-lhVOPMO7SpKLSixnoG6l6D2+3TkYn3q25ozsLNwczU/+6JOJJ+/aNgsPzMIy
-Ez4X2ZYUlHojuM59eBY+Vz87Y5+fgj2t787Vjc7CvNSdxW9jU9BJdAVv+OQs
-CDsomD06koIeUdU7OLxzEG9y8HdgYzKGBG2//0dqDhQLfA0Mmu/h0/O8Tt/M
-5iAlrdllXf9dnNbNmXlInoOzwstdK40SMODDw5tV5XOgNnFB46BqAvKopS4b
-qpoDqnY1iyqRgP/JxMmpNcxB2mrXtu3MeFw942NY9X0O5DOuSdQHxeOOl8b3
-B3/Pga7RpaVi1+LQY9WkuqrCPORsjN1jHngbx6PHKiwV58F3UE5iielt9OYb
-MvbbPA9L/F+87pe7jQG/uuwrt8/DStHc5UurozG6ribCch/hCVL0WrFozLv+
-oM3Xch7mnvnbmD6PRDrzoE/FrXnY3K7u8x87HA8+bH5Eip4Hw1bR9H8t4Zh5
-4mRTYuw8qP4hS6x4HI52DZ7yjknzIP48R6TZKBzbX2c28mfPQ8p/SXwnU8OQ
-4s+WPUSeh2THOb4K7VCsWF5aVzM5D25PWM6UgmCUazg0uX16HtqKKw4ZRwVj
-UGirdMrcPPCMXvUfdwjGA1ODV1z4FiD467nBSvFgLO2UkF4qtgAQ517r5xeE
-hU+8LxttWQCbw3+1PdMC8dF+1XVNZxag/qOHy9A2P6wwEbw557QADSX7FS6N
-+WL/afrwZpcFkC54oqyb74sKIUml4R4L8CDWeA9N1RcfV89Ya/kvgJ1Au+hd
-TR/MOVqR8ubuAtT08XcnnryKeWfMV6fWLECL5Z2x27TLWHt5a2Bt/QK8SSv/
-+CnjMg7eFGD9a1oAVldUx17Xy6iU+67YnLYA6kfHBGU5nvh0eJ3lsr4FCBqd
-VolQ9sRnVxh3gmYXoFdVc2f+nov4Msx7+bltbFhldbCcae+KP/fZhoSpsIn+
-UELSUtAVt07qTmTvYEO7quNy6yIXfOYg/rVXgw3NTjZ0aUEXzIM3KSf12FD0
-5bvwo49n8fHMbxmT02ywKR33Dt/vhIzXbfHnHdiQ4rv5K/x2RAWXMt5oZza8
-MflGi8l2xEff4ug17myQ0kwsHhByxIdvlJ9o+7DhwOo784n0M5h64bzyrgQ2
-BFV7asfZnib2W4uH5nfZELuQJXLvqx1Kdu4Vu5zMBlN7slS1iR2mGAmNv0pn
-g27eh7ksXVtM2vSsWOkZG6oO1lYZHDiF8X0D+6Qq2WCYtfmmmetxbEn58lKz
-mk2cp6ROj05YoajpO2nrOjZovpjzL7hhhbHkcJ57X9hwzDnQpzjVEm8/2Fgl
-2smG4b6kgBfd5njL0vEw/182TCxaVmCVYYqPjnbaOE+z4bDO84c7d5jiR0NL
-j5o5NoRfTU3uTTyK4/v1E6L5OCCQr17CO3kEbRS3fF8pxoGq0OayZd+N8Jrs
-40HvlRzIu693NWe1ESauXT/7fTUH7shV506fMMQaEdH1aes58Caqw6eVqY9q
-U+MOG7ZwQLD8dJ+x7CE88sfdK1yZAwoeKy+WeeuiywgjbFCVA7qxzdd+Nelg
-es+3vBfqHKjfdGr2RvhBXFJbOrZTjwPj5KZKtzVa2J0S4nfoDIc476i0NAfu
-wZnEuegnTsT734cbbwzejeIxXulLXTiwwXqEPB6xCw2Dz5Y1eXCg5khf2Hje
-DixyNVpsFcABB86NtDCSCn5xrJR4f4OI9+v9xE0Fyjhsu2/z2pscKPc51qa6
-ZhvKmCkbd0dyIENhrZLYfySM0lyReDaZA1KCy39MnNuAOTujs2tTOcDQvn/n
-xjoFLFNeVEzK4EDDIdCc+SGHf+X/fR97zIHrTnXJ539Jo9j6i0OmeRwY3Fl7
-YeKlFJIkBmcLn3NgWwaflFXIWrQXal9/rZADwuok3vzzq9FPwHz7j2IOdLlK
-qOeEiGMSp0FrbwkHzHZo6V97tQLrJz45cpADH6QjO23jRJA5utv7TCUHWAe1
-hDe7CSHPYH54ZTUHgpM9fj6/tAR3dWY9vdXEgVqKCrmNhx/FlF2m81uI55HH
-UuUND/68oazfRuWA3tsyJQ2JBXJt62QSz3cOBB5+OK0pPkN+rPCBodTOARlI
-lzxr+5cc5B2841gXB0gJpdGl43/INjX6Ib69HNh3cyz9Ut4oWV1SpDWLzoHJ
-zNxqMv8gWdydJlPPIuLxT084PdVP/v0x1WN8iAPkjbkLOYu7yI0iZz6tHeWA
-dHOslNrRr+Q8+03COr85YNA7xKmObyKHFo6edJ/gQKpEwVEvRgXZflHR0zv/
-OHAhlGXaplZM3mflN/1hhgOeQeVKpeRH5DV5WvqMeQ68EL9l8u5KCHk2gP/w
-Yh4eLZ7/uwZzzl5O5eX6W3PdfWV+rgtl4wetFnPtViUV/kSY6x9Cu8sNJLmu
-3nl/sHMt10W286KeUlzH51fY3Zfh2sD86PzQRq7fp7rsiVbjOmXzg5dNh7me
-nNDLuXmE66Pk8TSNo1wvNje6nWPGtXfAjNt1a67NvhxX2naOa2Fv8ezbIVzf
-KItJOVzC9ZF1ljn1G3j/39oeNaeWR3L9UMphQ+NfrjdUPrlncoqPO95tFSQK
-meuzE4ta60j8/2+HgnHZyDiuz4eYzxslcO1tXvRd+A7XEX+9Y+Pvcf1KY2Yu
-OZ3rmXLOt9znXMdThGOqarhOzb7gEl7HdbZXs45BA9dvV8fPNjRz3XlypQv1
-K9eb6JI6/XSuVYr8pB8zudYM65hxHuDaeFP666Fhri+7y0r/+cN1wL7gmcIJ
-rsOW9dO8/nKdkp8dPT3N9aOgRedKZ7l+fsz5YMA810Xy1esPsLkmTyjOcDhc
-/w8yc+em
- "]]}},
- AspectRatio->NCache[GoldenRatio^(-1), 0.6180339887498948],
- Axes->True,
- AxesOrigin->{0, 0},
- PlotRange->{{0, 2.5}, {-0.3183098745627588, 0.}},
- PlotRangeClipping->True,
- PlotRangePadding->{
- Scaled[0.02],
- Scaled[0.02]}]], "Output",
- CellChangeTimes->{{3.560155132339073*^9, 3.560155136020277*^9}}]
-}, Open ]],
-
-Cell[CellGroupData[{
-
-Cell[BoxData[
- RowBox[{"Plot", "[",
- RowBox[{
- RowBox[{"DWh", "[",
- RowBox[{"r", ",", "1"}], "]"}], ",",
- RowBox[{"{",
- RowBox[{"r", ",", "0", ",", "2.5"}], "}"}], ",",
- RowBox[{"PlotRange", "\[Rule]", "Full"}]}], "]"}]], "Input",
- CellChangeTimes->{{3.560154655600813*^9, 3.560154655644041*^9}, {
- 3.560155161815674*^9, 3.5601551791128883`*^9}, {3.560158982837762*^9,
- 3.560159022468958*^9}}],
-
-Cell[BoxData[
- GraphicsBox[{{}, {},
- {Hue[0.67, 0.6, 0.6], LineBox[CompressedData["
-1:eJxF13k0FX/4B3BLqaRSStJCpUilTbK/lRZJpdAiUSJrtihfKSLJGiJLQiFb
-4V5UlowsIWtIksi+3mu71879ze+c3/nd+WfO65zPPDOf5/M8n5nZamh90ZiL
-g4Mji5OD43/Prx8Z8hwwslD2e58aV+ZPJ4rrRTUui55GqZlVRglp1rd8d1HR
-q3g5rNqeR3pl0Md03iZD9Kap7k8l/brAT2SxqBm8N7m9DSe9Z9jIn2VijTMJ
-AvIupP+ujFqY+2UHoUK+aX3SYu+WyczMOiC/9XG7LGlLJYc7kyL/YZ21J5OP
-dEbDvziGqjMYXQeO/PGjEzMWGi2jJi4ANSQxlnQMZ4vX5InHKFIrPXOb9PGq
-wB0Tv9xwRfjPtu2k+8JOfWWYPUFnWPe+Zl864Ws0rzc+64H9MrP2PqTrZ02D
-R0S88Cg4PK3Fh07cK92yf5jiDVm5M0nOpIVfNFTQVH2xrvhqpyBpQ0kV7kGT
-5/ig1SYq700nuGadahjOAViv6CJR5EUnGv7pbGacCIRAXZXjSdKO75dnj/0K
-Amd12azCMzqhHtTDMxbzAinrow6metKJTY5ftUfNgnGq8W6OMOkC1fsjw7Mh
-iP/PrqTDg04s+dMhThMJR07lFNc+dzrRVPDFYagvHE8PrJy760Ynkt+FFQ1S
-IjAXXnop4zGdOGd3zmBANRKWSo5yW13pxMtlOS97TaLRcXQFPduZTmTeV7e7
-KxSDgSMuUnUP6ETkwJEN484xuFMg297lRCesateYjJ14A3nNhvhpRzqxJrKM
-a/TXW7xs7e/+Z08nZld+TLJRjAVfS8NY5V060fU4VnMkJha5M8luGXZ04qPJ
-o6hhsziI3zI/ZGNDJ3QPScvTZ+MhoG0nkGhBJ958j7YdFEmCs5Z3tMBNOsEI
-y5Js803C3OTSNd4GdOKkSUVn3XQSzErsZGeu04kh7kmdnPpkXPB5k1KqSydk
-lc7Je3m+x66X0Q7rtelEbdocl/hwGnqU8vp9T5DPi2X6K6TSYbVppClalU4o
-7AyW9zQi3e3w7v1ROuHPSBp3+pEO6V7+3gwlOiEd+NPYMIUCB+fPav8dphMu
-3/eoHzDIQIa7ZVv7djqRSPm0MyUkAzJZX0dfbKUTdWHHuHdUZiAzqIZPRYSs
-X5MreUJymdjOlWbyVJhOfOd+IsUhkAUTj/SZCn46Iaj0Z03Nt48wz9DrVZyj
-EalpXn8spXJgf+2caVk1jRhYM+CrfjMHds5GaowKGiF+T11ZIjgHX5WoyhvL
-aUSM4vI3ndM52FEeck+3iEYElfve1i3JRbRw9/moTzTiXufz0VN6X2Bo1O2j
-FE0jopdXuK8+UYAuOzvVfnMaoaz0QM9VvwDnGT5CziY0osVq9+Hh+wWw0ulR
-4zWiEUJ1Pj2VyQUoa10oErpOIwLDNE4/4/+KMzvEbyycoxHuO6pXsv58hZ9S
-iAD9AI0wwY/wIbsitDnRXK0YQwTH24mcPz5FOLuRj9NjZIgIW7Sp5XtcEU41
-92uEDg0RZeW3tyQ1FkHk1iLNtK4hQkJ77q2xfDFyLd4ofmgYIvpMxd+3cpYg
-lfr7XkfGEGEe9DC/NuAbHPsOJ6tYDxFW3bs6M1PL0e90qEupeZD4p/976dfy
-crhvyPw51DBIaDU9k6rqKkdvjYBVaM0gIVvR+1+38Hfcj//xuq14kOCmxPML
-en5H3V/6UcW0QSLMeavyfYMKPKint0u4DxJfBTaEy/FX4a0WK7ldfJDg3fpE
-vE22FpxpMhRD4wHCn79C8s2NeqyukGhOrO8jDP9oE9yVjfhyWWJH9nAPYVvy
-xNBavBkWUjvPBdG6CIoNn5bQ6b8IWvGSbwujg3jzwGJg/+V/UJY/9C36yz/C
-LXruAmOoHdPdhy6vPP2X4Kj/p5Ln1Al78TUvrPf+Jnba33bUluxGV3R+SF9l
-A9HFDL20/EsP3jZaK6u21BLq3eYPg2z7MHshT1jgXTlxritg95YtA5B//O71
-gWuFxECRoWvZh0HoLalMePvkMyES2GTunzmIgNxHo1u0PhPa+ud0tHMHwWdm
-UJO99TORPyUn+a9sEHlTm1+dIT4RQXtXN0x2DsJx4GNO++xHQiGUEBcXHkIt
-f8C1Aacswtd8U43H0yHoyUTvZQVQCanVjZtVDWhQC9Uv9FdNJnr8etxCjGlw
-5awVezmRRLzmnezttaCBo4Xe8ScpieBbJET1daRhf+fg9LbVScTAxNWTjYE0
-jBzvUovtSCDiW/5amRXToJjQKOwdGE9sSuwmnu+i45ndolWtIm+IBrEJsY59
-dDDCy2LHmmMInzc83tIydGgrjf04/TKGmIkQ1/59jPT0HreCFTFEk59Z3zY9
-Olxj3GrjWa+JF3dp/B/96QgL/cO82B5O8IJ5s2WcDtnlBxwE+YIIx6MSSz/O
-0KFSk8Wj6xVI9KheS33OOYwbCR1TqksCiUK1rzPHVg0jMfrZFb7FAYTTRb8X
-SZLDKOCLGxrh9yUGjXeU3Ls5jJEPTpYC+k+ISl8didU1w1AZ1W0wGjQhTtF0
-J7fpjMC1ulN280N3RK11r3e5MoL9PkzLB3JPwFBISWu5NgKbXB/7g8wniPGe
-NXlpOIICf7GmFPOnmBGPbFpmM4L0qK09P8564b3h3+wR7xGI6gp+E2b6g7/J
-wDm/YAT/ZiWs/LuCcZvleWVj8QjCwlyfuMiGIG9nurRj6QiEvgtaivqGwMyB
-g3agmhx/j77R7uBLFK55cz3+zwhufPs1dflBKOzPdij7TpDxs6cenFgIR1Oh
-EafunlHwf09Qv10ehZGjzcIP941C6FrXg+cro7Hk63npmIOjUKMb/HbQioYM
-IW/SIzuKMPOntk4t0QjO5a+yOzGKIZumAqpDDM5n5oX56I/imbuyksWjN/gW
-v3b/l4BRjCxbUDe4H4tWMe/T/16Mgs9PQaAlJBbMWJYhd+go7Ctv1stmxkLs
-7UDI6dejKE54YxI0Egu3qILZn0nk+K5H60pN4qAcallKLxxF4ur9h1Q045H1
-rFhflDkKxiG75vxlCUhUrxDynBqFYlXmNZ+dCYjgq6ujzY7Ckfqh8ZBqAlwD
-2k7mco1hz2aRMSHnBGiEzUhd4h/D/lCv04eHEtCRsJ/ls3sMIyrfqgtLErGy
-NDJ68uYYIt3Ff77QTwanV+xVfeMxVGZsaZ65lwyGerJAiekY7Pmpr/A8Gb+r
-P3kGWY9BaDmvnyGRjLjGOuu9D8n4kkszVm1JgXzPUpVbYWOoPfj0hsDPFNxe
-7PCvunoMjh51ndskP0BYpPzX7A8yfkP/eR6lD6iW3Vwj8XMMZbQZ46rzHyBj
-WfLF7c8YRFcOfOFx+ACe+nWvZPrG0Me5QV88/wPiYz7qRHGOg8/9jankmVR0
-KUx9vyM9Dh6dVzRPzTSE6WgURhwZh+lxmR8TBmnQsI7JLpUfh+Oz0ip16zRk
-vlVL3Hp0HGtj3uOtXxqeLAv3+Hl2HJE8Uaf4y9Mg9ktORclkHIqEm+N/m9Nx
-y+5BFl/EOOS/BjkGeqRD5mBr1q7X4+BeZLHT0ycdvGMqH0/GjEOqpO+lTWA6
-0m15Prm+G8dqA9ufW16nY84m8DODSsaTZRn0Z6TjpXVCbkvFODw/nzdWaEuH
-2T7evOnqcbjINJdYd6VDcdgyT7BuHGFxzgjrT0en1cEvmk3jiKbFvygZT8d+
-qy/5xV3jYP1XdyloCQXllvVf38+Pwz3zJ9f0Hgoi98gUlnMwEOUl99D2AAU2
-Q2GFPdwMRBaIm7QepkDQUr9IlJeBmBFEhyhTYGjRXxwsyMBkdMIRp/MUzJmx
-Sp33MaBk/7JI1JqCml2GZeEHGbiTN5PfYUfB2/7iso+HGVBouuYafo+C02be
-5SMKDIiLU8y7H1Lw0nRdBfllAEWzEftJHwr2meyu1rjJQOVHbuFl7yhoNROx
-/mDEQFtKrKBeEgW+lgL8K00Z4DPcqPr2PQX9trMXaqwYODa069TSDApinSt+
-XnBmoNh2IE+BoOCiC3GP6sKAuZn1tEQhBRxuGesF3Bk49Gs2clkJBdc9X11t
-8GKgZrhQOv47BeuDLP5eCmWgbuUB05M/KfgWbPDoUwQD9ylCUTW/KLAP1RIR
-imJA2+WInUYzBT8iFW7+jmPg6H0q/9Y2CnwSlndfozLg96O4V6aPAvlklkde
-FgPJpb611wco6Hs/vnNzNgO8Ubn+jkMUnKD+MW0lGNj9rfud1wgFrLyUoRtV
-DATOBy3ZOEVBKhHt97WWgbDhz3rt0xToFb6Q2tbAAOOeQ/6rWQqySx/YdDYz
-cFBUvWhogQLT79arj7cycGWflL07BxWCVbeoce0M6KWqnF3JRcXdujMM4z4G
-HjU+e8GxmIptPxHybZCBLkvvlaY8VNT+OiQjPszA+7MXq4qWULH370bHPiYD
-w2WVB67wUtHbOzhrzs3EWPSLy8yVVLwcaIus4GFilzudwlpFxXFavdIeXiZ+
-GMNzgZ+KN2O5LjR+JjgWXcv/vYYKTWaa6Lm1TFS8H937UYCKhcnYr6nrmXi+
-T1Dacy0V1+Z9uG22MOH5OkZkiSAVvByucbVbmVhyr27NR9KfuexPHNjBBCe3
-5Oer66nkfmHaEyjBhKiCn+AY6XVL9TzHdjOxeF2HtIsQFcW8mhJa+5g457lS
-mHMDFXYrjpdnHGQi6tlU5X3Sovyy5mtlmIhY9VS7k3TNmj3LHeSYaF/z+pOq
-MBXKd+70v1Jkwubb7dgXpFNLU0sLwYRa1Oee36S3bBuJ7z/GxNOivZ6CG6nw
-dz7whP8kEycMbnqfJr3QaGd45DQTMRF0ph3pOwcyVfQ1mKBE2dQEkW7xYW7x
-OE/eP9htUxJpjR6Z+ZSL5P3iiY4s0nkqjn/qdJjQnikUzya951V29vQVJgzy
-xEeppCOZM6Giekykr/U6Hkt6uabivVMGTFzf4ybuTdop+aG2lSETXntSo01I
-DywiDoYYM6GhnEhVIK1rwLE6z5SJIu/lt3lIf88+OtxhQeZP4S5RSs5Pfq17
-1TJrJqQdnYpdSCdbFafst2PCwiLPSYq0cPli78sOTJQoNfTWk/nz2n7K9JEj
-E3pC11fbkJ56+Oxk/AMy/2Jjc1ykTZvKxSofMZHlLZXuS66Pmp9G+wYPJnSX
-8z55Qq7f514/QuUZE5X329LGyPUWP1bz2sSHCbHOsuzLpHkmL+hmBTLB0D3+
-iGcdFfcuvJBtCWbiXTbfcU2yfnpSGgS5w5hIzdVZCCDrq+TG5XrNKCb2nmu8
-PbmaCuncMMr9N0w07z8uupl03Lrm51FxZD6uRXfJkfX65LuexlAyE8O7veOu
-kPXNEIuSFEhl4pKYZoruCipuubQtlacwMb7Upl6Lj4pj0oYlnp+YkFOQjNpN
-9gfVPzY2NYcJ4T8zZ3iXkf3V3/X45xey3pYHHfhH9hNnlIny9mImuIpcc2zJ
-/itYcudTfi0TKrxXbuwh+3W/YWpIdz1ZX6k6ZiVkP8fkDd/l+8XEz93jjTrz
-FLjY2u3X/ctE32WnCb0ZCpT/3E9i9jMxOUbpXMSgIDfV7dUe7gls2+utX9BN
-AY+GjnAozwRm3r7xSuyk4EK/eDgn7wTO3LTm8mynoHd71ctG/gnESOZM7/xL
-gUD4+iDXLRM4ub7r23g9BRbu7581yE0gIvCYcCm5n34UdVkCpQlsXCP93/ov
-FHDmX3iapDKB84krwvRzKAidmnB/dGoCa70WH/yRSUHRnaMuEjpkPNW2KDVy
-/954pdHB2WYC5vXjEU8DKLjNTGT03J2A2MEWzZN+FFCCHty9cH8ClXdbZee9
-KDhVJWq789EEtgt8ydV0p+DuMQvLWp8JSPRnaz93oKBiD8ctsYQJJO62WGN7
-hYIHnJKalX8nYFXvQJ1aT0HP8eVNO9snECsUlFAgQIHmsyGDx10TEFTTmnNd
-RYHYqjRrmcEJGDWVtg3wUFC5Sfp5zNQEqO2BZ4yY6dgiq1xtLzCJuVTz6PHa
-dBRaXdTYfHoSvowjjzLcyPd5i5OaVeYkwq6sDDSKS4PM/aQVzz9N4m3d807d
-8DTcXNNUl5YziRkFrX8n/dPwWe3w9dGCSXBW3lvE4ZiG21nDtvZVk9hkvVpy
-Sj0NX/2NXjn1TCJgT/jjZcOpcDx6jv5UaAp2f9cJz+9LRfe7rcHRzlPIaWul
-+IW+R7DEhkXZOtOQbxbTXF6QiICI26U3RWeQOnSfK+x3HHSDr1fz9c5AeCvH
-qFFfDKaOx09H589Cv2JBuNEwHJ3dxxwLPeegPLzJKpIzCMeiq99K+sxh+4Yx
-0b2dgYi5olsV5D+Hc3ybraqKA8nrbLfdCpkDtn9RU3oWiOb0mEru2Dlc7dku
-9HdVIOoeLIieyJ/DIx2n3sPbAlC4Oqe8lDEHYuTRqh2JvnirdGBj1c15CKdG
-B2t1e6DwHK/7rNE8DJ31Z4yyPNBh0Dmwy3Qe211ODr728MB2t5Ccp1bzZL3y
-invu9EDct+mrKg/m0aCVslrR/AnizxeGZQTPo60iesWVaTck3tQSjCidR83v
-r3avjriizG7Po7Lv87Ds7A6d5HNFn/vinomqecjKDqznKnGBxLtPWVoN8xDI
-jDgot8kFSQMbdVa0z+P90+uJE9UPkXy364XrzDwea/Cos048wAcPh9UmexeQ
-tfCVt8r5HoYUr7t57F9AoNoukxnpe9jDOD4ee2gBOzcsGfGiOSDZcO3Pf3IL
-cIms7ag1cEAiMsJ0Ty7gjrpq9ZLT9oibHhE5Z7AAe+W+cmMpO3SlNwVYGi4g
-YMWHlNYhW2w3LeD0MV7A0ohACSLFlvxPft5ZarGA87JG5z9L2iI6QyrhqOMC
-ppN+1QXvtUHEHUupw4ELqPlk4+zhb4lmMe1oreAFZO9rEckTssSGFgV+u9AF
-HH3+u7Ey1gJh6svHUl8v4E90q4xLnjlCxJOzJJIX4HbSsFB+yhQB7b2Km4oW
-sC7qx9b4cPKlrnPrDDdzAZEoeL9qjQHenm/RM55agOaP49mRmvrIO61jVTq7
-gFMd9tL6gdcxpnQq0IeLhfixWGOPDXrQ27n7lwA/C6sGG+RDjl/FfdG4PgcB
-FpYHahyfCL6CIOHNM78EWXAMXxoR3nsZpStXbY7czIKrqvCG1BeXcHByzFBs
-NwuK6xVoyku0cXbUwv6pFAt+Q5b/Osy0YDrY5dF3gIX0Jc7zjdUX8bqtMfG9
-LAvqjxdiP0ZfwNKyHLr0SRbuyF3WXpp7Dq1hbk4nbrKQn5Cf9jnkJKaDZn0S
-jFjgoGaVx787gbV+9q+XmZLjf37/YJ53HKcf3y6osmLB+EwfxWviGDLN1Hku
-ObNwUVvRblYdqLlVtP6zCwt55qa31+YqYeC64i5hdxY4j0obrz2sCJGLUhqt
-XiyoMlqUuk/KwVt+TdDtUBaqO6r/8+04hHhpn9iyCBbOm+i0R789iAKpRVmS
-USzsk2sr/2F5AMxtE7/oceT1QtepVyEF/s3W/ZqJLDz6tyG3lH8PJNf3zVBT
-WAj/oaKYPLELN5Y3b75PZeH35/1bFnHshNNirX2/s1ig39stH7pdDCGsChWF
-bDKfA+PHG69tw/fxL7dYBAtGHLHyIuu2oJsm43CziAXluv0aZ7g2gqMv7WnR
-NxayJLdKKG8QwuGWN0meVSxo+Izs7i1eA34p06m0WhYkvCSsqx/zY8hF6lRT
-PQtWl2lrRr/xoewHI4TjFwuVOv58U/VLEbc9t0uimYXaB3wdrZKL4erw+NCF
-vyw0HVFJrqvghF7pKbf//rEQNMI9xWE1ryy7YeWPN50scF3WbpH8Pam81qJB
-5HsPC6Uu2/Y9/D2mPJIXYTXWz4Lac2XFuNYh5cqVN78I01j4JDj5YYlMt3Li
-DXE+1REWKqJiNlPv/lV+QqXpWoyz8L77sn34zzrlG4syk15MsHBiIHWv3s4S
-ZcVLTlO50yz0cfgMG36iKgslqpzqmiPrN6Y+dqtNkPKMM/cZHg4OFY7/O/ri
-b9tFcLLdWF0eLsXNNlU0oO8SD9vmxZueJvCx/Xu5zFe1DWx/kw7vaxFmO/P6
-3CrbTWwHpBXqh4uwraZ1fq5/B9ufI0yP+BxkO2zXqw9VZ9hmjJ+Mdz/L9vn8
-sUi582zzaKn7xl9k28F52vzhVbYv1lyW2GvCNp/D2lhfN7ZdCvzCzmSzfXaj
-Tvx3Mc7/91Gr0murvdiO3mQoVslkW6wo4eW5a1zs+e4tlKzLZ/v2+KIf5ZLc
-/29Dypio13O2Ld205tQD2XbQyvzF94LtZ0wH/4CXbKfKTc+GvmZ7+iur8V0K
-2wF1fH7FpWxHxN4xfVrOdqx9tapaBdsfBQNmKqrZbtEVMK3/ybZ45wbVjk62
-92c6bYnrZlve48+0cS/bGuKv0/sH2LazEN0yOsq2s+Ljaeo42x4rOhrsmWyH
-pcX6TE2x/dZ1kUnODNspF4yPOc+xnbnt22blBbbzx3dOs1hs/w+pgDiQ
- "]]}},
- AspectRatio->NCache[GoldenRatio^(-1), 0.6180339887498948],
- Axes->True,
- AxesOrigin->{0, 0},
- PlotRange->{{0, 2.5}, {-0.9549296585513659, 0.07073552342169737}},
- PlotRangeClipping->True,
- PlotRangePadding->{
- Scaled[0.02],
- Scaled[0.02]}]], "Output",
- CellChangeTimes->{{3.560155172170507*^9, 3.5601551796044407`*^9}, {
- 3.560158983450157*^9, 3.5601590230207157`*^9}}]
-}, Open ]],
-
-Cell[CellGroupData[{
-
-Cell[BoxData[
- RowBox[{"DWr", "[",
- RowBox[{"r", ",", "h"}], "]"}]], "Input",
- CellChangeTimes->{{3.5601552083271513`*^9, 3.560155253227319*^9}, {
- 3.560160694674526*^9, 3.560160694745482*^9}, {3.560161180197549*^9,
- 3.5601611810400257`*^9}}],
-
-Cell[BoxData[
- FractionBox[
- RowBox[{"If", "[",
- RowBox[{
- RowBox[{
- FractionBox["r", "h"], ">", "2"}], ",", "0", ",",
- RowBox[{"If", "[",
- RowBox[{
- RowBox[{
- FractionBox["r", "h"], ">", "1"}], ",",
- RowBox[{"-",
- FractionBox[
- RowBox[{"3", " ",
- SuperscriptBox[
- RowBox[{"(",
- RowBox[{"2", "-",
- FractionBox["r", "h"]}], ")"}], "2"]}],
- RowBox[{"4", " ", "h"}]]}], ",",
- RowBox[{
- FractionBox[
- RowBox[{"3", " ",
- SuperscriptBox[
- RowBox[{"(",
- RowBox[{"1", "-",
- FractionBox["r", "h"]}], ")"}], "2"]}], "h"], "-",
- FractionBox[
- RowBox[{"3", " ",
- SuperscriptBox[
- RowBox[{"(",
- RowBox[{"2", "-",
- FractionBox["r", "h"]}], ")"}], "2"]}],
- RowBox[{"4", " ", "h"}]]}]}], "]"}]}], "]"}],
- RowBox[{
- SuperscriptBox["h", "3"], " ", "\[Pi]"}]]], "Output",
- CellChangeTimes->{{3.560155230596974*^9, 3.560155253885023*^9},
- 3.5601606952114277`*^9, 3.56016118252979*^9}]
-}, Open ]],
-
-Cell[CellGroupData[{
-
-Cell[BoxData[
- RowBox[{"N", "[",
- RowBox[{"Expand", "[",
- RowBox[{
- RowBox[{"-",
- FractionBox["3", "4"]}], " ",
- RowBox[{
- SuperscriptBox[
- RowBox[{"(",
- RowBox[{"2", "-", "q"}], ")"}], "2"], "/", "Pi"}]}], "]"}],
- "]"}]], "Input",
- CellChangeTimes->{{3.560160709698295*^9, 3.560160723505558*^9}, {
- 3.560161185019305*^9, 3.560161189166279*^9}, 3.560237508328278*^9}],
-
-Cell[BoxData[
- RowBox[{
- RowBox[{"-", "0.954929658551372`"}], "+",
- RowBox[{"0.954929658551372`", " ", "q"}], "-",
- RowBox[{"0.238732414637843`", " ",
- SuperscriptBox["q", "2"]}]}]], "Output",
- CellChangeTimes->{{3.560160720336149*^9, 3.560160724559553*^9},
- 3.5601612038241377`*^9, 3.5602375092839746`*^9}]
-}, Open ]],
-
-Cell[CellGroupData[{
-
-Cell[BoxData[
- RowBox[{"N", "[",
- RowBox[{"Expand", "[",
- RowBox[{
- RowBox[{"(",
- RowBox[{
- RowBox[{"3", " ",
- SuperscriptBox[
- RowBox[{"(",
- RowBox[{"1", "-", "q"}], ")"}], "2"]}], "-",
- RowBox[{
- FractionBox["3", "4"], " ",
- SuperscriptBox[
- RowBox[{"(",
- RowBox[{"2", "-", "q"}], ")"}], "2"]}]}], ")"}], "/", "Pi"}], "]"}],
- "]"}]], "Input",
- CellChangeTimes->{{3.5601608470246563`*^9, 3.560160853545632*^9}, {
- 3.560161190598509*^9, 3.5601612011456413`*^9}}],
-
-Cell[BoxData[
- RowBox[{
- RowBox[{
- RowBox[{"-", "0.954929658551372`"}], " ", "q"}], "+",
- RowBox[{"0.716197243913529`", " ",
- SuperscriptBox["q", "2"]}]}]], "Output",
- CellChangeTimes->{3.560160854392119*^9, 3.560161202029501*^9}]
-}, Open ]],
-
-Cell[CellGroupData[{
-
-Cell[BoxData[
- RowBox[{"DWh", "[",
- RowBox[{"r", ",", "h"}], "]"}]], "Input"],
-
-Cell[BoxData[
- RowBox[{
- FractionBox[
- RowBox[{"If", "[",
- RowBox[{
- RowBox[{
- FractionBox["r", "h"], ">", "2"}], ",", "0", ",",
- RowBox[{"If", "[",
- RowBox[{
- RowBox[{
- FractionBox["r", "h"], ">", "1"}], ",",
- FractionBox[
- RowBox[{"3", " ", "r", " ",
- SuperscriptBox[
- RowBox[{"(",
- RowBox[{"2", "-",
- FractionBox["r", "h"]}], ")"}], "2"]}],
- RowBox[{"4", " ",
- SuperscriptBox["h", "2"]}]], ",",
- RowBox[{
- RowBox[{"-",
- FractionBox[
- RowBox[{"3", " ", "r", " ",
- SuperscriptBox[
- RowBox[{"(",
- RowBox[{"1", "-",
- FractionBox["r", "h"]}], ")"}], "2"]}],
- SuperscriptBox["h", "2"]]}], "+",
- FractionBox[
- RowBox[{"3", " ", "r", " ",
- SuperscriptBox[
- RowBox[{"(",
- RowBox[{"2", "-",
- FractionBox["r", "h"]}], ")"}], "2"]}],
- RowBox[{"4", " ",
- SuperscriptBox["h", "2"]}]]}]}], "]"}]}], "]"}],
- RowBox[{
- SuperscriptBox["h", "3"], " ", "\[Pi]"}]], "-",
- FractionBox[
- RowBox[{"3", " ",
- RowBox[{"If", "[",
- RowBox[{
- RowBox[{
- FractionBox["r", "h"], ">", "2"}], ",", "0", ",",
- RowBox[{"If", "[",
- RowBox[{
- RowBox[{
- FractionBox["r", "h"], ">", "1"}], ",",
- RowBox[{
- FractionBox["1", "4"], " ",
- SuperscriptBox[
- RowBox[{"(",
- RowBox[{"2", "-",
- FractionBox["r", "h"]}], ")"}], "3"]}], ",",
- RowBox[{
- RowBox[{
- FractionBox["1", "4"], " ",
- SuperscriptBox[
- RowBox[{"(",
- RowBox[{"2", "-",
- FractionBox["r", "h"]}], ")"}], "3"]}], "-",
- SuperscriptBox[
- RowBox[{"(",
- RowBox[{"1", "-",
- FractionBox["r", "h"]}], ")"}], "3"]}]}], "]"}]}], "]"}]}],
- RowBox[{
- SuperscriptBox["h", "4"], " ", "\[Pi]"}]]}]], "Output",
- CellChangeTimes->{3.560161212213023*^9}]
-}, Open ]]
-},
-WindowSize->{740, 867},
-WindowMargins->{{Automatic, -1324}, {Automatic, 61}},
-FrontEndVersion->"8.0 for Linux x86 (64-bit) (November 7, 2010)",
-StyleDefinitions->"Default.nb"
-]
-(* End of Notebook Content *)
-
-(* Internal cache information *)
-(*CellTagsOutline
-CellTagsIndex->{}
-*)
-(*CellTagsIndex
-CellTagsIndex->{}
-*)
-(*NotebookFileOutline
-Notebook[{
-Cell[CellGroupData[{
-Cell[579, 22, 1154, 36, 88, "Input"],
-Cell[1736, 60, 937, 30, 57, "Output"]
-}, Open ]],
-Cell[CellGroupData[{
-Cell[2710, 95, 309, 8, 30, "Input"],
-Cell[3022, 105, 7339, 126, 236, "Output"]
-}, Open ]],
-Cell[CellGroupData[{
-Cell[10398, 236, 459, 15, 30, "Input"],
-Cell[10860, 253, 294, 6, 30, "Output"]
-}, Open ]],
-Cell[CellGroupData[{
-Cell[11191, 264, 343, 11, 30, "Input"],
-Cell[11537, 277, 346, 7, 30, "Output"]
-}, Open ]],
-Cell[11898, 287, 774, 20, 50, "Input"],
-Cell[CellGroupData[{
-Cell[12697, 311, 265, 7, 30, "Input"],
-Cell[12965, 320, 7987, 137, 221, "Output"]
-}, Open ]],
-Cell[CellGroupData[{
-Cell[20989, 462, 414, 10, 30, "Input"],
-Cell[21406, 474, 9020, 153, 223, "Output"]
-}, Open ]],
-Cell[CellGroupData[{
-Cell[30463, 632, 247, 5, 30, "Input"],
-Cell[30713, 639, 1094, 35, 65, "Output"]
-}, Open ]],
-Cell[CellGroupData[{
-Cell[31844, 679, 404, 12, 54, "Input"],
-Cell[32251, 693, 318, 7, 30, "Output"]
-}, Open ]],
-Cell[CellGroupData[{
-Cell[32606, 705, 541, 17, 54, "Input"],
-Cell[33150, 724, 238, 6, 30, "Output"]
-}, Open ]],
-Cell[CellGroupData[{
-Cell[33425, 735, 79, 2, 30, "Input"],
-Cell[33507, 739, 2072, 67, 118, "Output"]
-}, Open ]]
-}
-]
-*)
-
-(* End of internal cache information *)