diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index d31b6be383b80a2698b63d27308f6fee9b23518f..09f796a8f37a9c015135f4aab3f821c2e862bdc9 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -93,8 +93,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
dv[2] = pi->v[2] - pj->v[2];
const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2];
- pi->div_v += faci * dvdr;
- pj->div_v += facj * dvdr;
+ pi->div_v -= faci * dvdr;
+ pj->div_v -= facj * dvdr;
/* Compute dv cross r */
curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1];
@@ -211,10 +211,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Balsara term */
const float balsara_i =
fabsf(pi->div_v) /
- (fabsf(pi->div_v) + pi->force.curl_v + 0.0001 * ci / fac_mu / hi);
+ (fabsf(pi->div_v) + pi->force.curl_v + 0.0001f * ci / fac_mu / hi);
const float balsara_j =
fabsf(pj->div_v) /
- (fabsf(pj->div_v) + pj->force.curl_v + 0.0001 * cj / fac_mu / hj);
+ (fabsf(pj->div_v) + pj->force.curl_v + 0.0001f * cj / fac_mu / hj);
/* Are the particles moving towards each others ? */
const float omega_ij = fminf(dvdr, 0.f);
@@ -309,10 +309,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Balsara term */
const float balsara_i =
fabsf(pi->div_v) /
- (fabsf(pi->div_v) + pi->force.curl_v + 0.0001 * ci / fac_mu / hi);
+ (fabsf(pi->div_v) + pi->force.curl_v + 0.0001f * ci / fac_mu / hi);
const float balsara_j =
fabsf(pj->div_v) /
- (fabsf(pj->div_v) + pj->force.curl_v + 0.0001 * cj / fac_mu / hj);
+ (fabsf(pj->div_v) + pj->force.curl_v + 0.0001f * cj / fac_mu / hj);
/* Are the particles moving towards each others ? */
const float omega_ij = fminf(dvdr, 0.f);
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index 6afb9d8d38a4fc7f1d38b7286720ddb7f3c51ab4..b3b81a9a0dfe41e7bfafe51050d6f7cf7157e31c 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -16,8 +16,8 @@
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
-#ifndef SWIFT_RUNNER_IACT_H
-#define SWIFT_RUNNER_IACT_H
+#ifndef SWIFT_RUNNER_IACT_MINIMAL_H
+#define SWIFT_RUNNER_IACT_MINIMAL_H
/* Includes. */
#include "const.h"
@@ -38,33 +38,31 @@
__attribute__((always_inline)) INLINE static void runner_iact_density(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
- float r = sqrtf(r2);
- float xi, xj;
- float h_inv;
float wi, wj, wi_dx, wj_dx;
- float mi, mj;
+
+ const float r = sqrtf(r2);
/* Get the masses. */
- mi = pi->mass;
- mj = pj->mass;
+ const float mi = pi->mass;
+ const float mj = pj->mass;
/* Compute density of pi. */
- h_inv = 1.0 / hi;
- xi = r * h_inv;
+ const float hi_inv = 1.f / hi;
+ const float xi = r * hi_inv;
kernel_deval(xi, &wi, &wi_dx);
pi->rho += mj * wi;
- pi->rho_dh -= mj * (3.0 * wi + xi * wi_dx);
+ pi->rho_dh -= mj * (3.f * wi + xi * wi_dx);
pi->density.wcount += wi;
pi->density.wcount_dh -= xi * wi_dx;
/* Compute density of pj. */
- h_inv = 1.f / hj;
- xj = r * h_inv;
+ const float hj_inv = 1.f / hj;
+ const float xj = r * hj_inv;
kernel_deval(xj, &wj, &wj_dx);
pj->rho += mi * wj;
- pj->rho_dh -= mi * (3.0 * wj + xj * wj_dx);
+ pj->rho_dh -= mi * (3.f * wj + xj * wj_dx);
pj->density.wcount += wj;
pj->density.wcount_dh -= xj * wj_dx;
}
@@ -76,24 +74,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_density(
__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
- float r;
- float xi;
- float h_inv;
float wi, wi_dx;
- float mj;
/* Get the masses. */
- mj = pj->mass;
+ const float mj = pj->mass;
/* Get r and r inverse. */
- r = sqrtf(r2);
+ const float r = sqrtf(r2);
- h_inv = 1.f / hi;
- xi = r * h_inv;
+ const float h_inv = 1.f / hi;
+ const float xi = r * h_inv;
kernel_deval(xi, &wi, &wi_dx);
pi->rho += mj * wi;
- pi->rho_dh -= mj * (3.0 * wi + xi * wi_dx);
+ pi->rho_dh -= mj * (3.f * wi + xi * wi_dx);
pi->density.wcount += wi;
pi->density.wcount_dh -= xi * wi_dx;
}
@@ -148,7 +142,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
/* Compute sound speeds */
const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
- float v_sig = ci + cj + 3.f * omega_ij;
+ const float v_sig = ci + cj + 3.f * omega_ij;
/* SPH acceleration term */
const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
@@ -225,7 +219,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
/* Compute sound speeds */
const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
- float v_sig = ci + cj + 3.f * omega_ij;
+ const float v_sig = ci + cj + 3.f * omega_ij;
/* SPH acceleration term */
const float sph_term = (P_over_rho_i * wi_dr + P_over_rho_j * wj_dr) * r_inv;
@@ -245,4 +239,4 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->force.v_sig = fmaxf(pi->force.v_sig, v_sig);
}
-#endif /* SWIFT_RUNNER_IACT_H */
+#endif /* SWIFT_RUNNER_IACT_MINIMAL_H */
diff --git a/tests/Makefile.am b/tests/Makefile.am
index a35344790967d3cc3f5d299c6f9af3cb4c2fab43..7edf6e0659c91e66c7eee3e3815e8f83e8551ec2 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -45,4 +45,4 @@ test27cells_SOURCES = test27cells.c
# Files necessary for distribution
EXTRA_DIST = testReading.sh makeInput.py testPair.sh testPairPerturbed.sh \
- test27cells.sh test27cellsPerturbed.sh
+ test27cells.sh test27cellsPerturbed.sh tolerance.dat
diff --git a/tests/difffloat.py b/tests/difffloat.py
index fe1118dfa050d3507137febf790eb5a698cfdb57..bbb7c95a1e77e04bbe21bec6dc6c5d529cd77c70 100644
--- a/tests/difffloat.py
+++ b/tests/difffloat.py
@@ -25,20 +25,24 @@ rel_tol = 1e-7
# Compares the content of two ASCII tables of floats line by line and
# reports all differences beyond the given tolerances
-# Comparisons are done both in absolute and relative values
+# Comparisons are done both in absolute and relative terms
+
+# Individual tolerances for each column can be provided in a file
file1 = sys.argv[1]
file2 = sys.argv[2]
+fileTol = ""
-if len(sys.argv) >= 5:
- abs_tol = float(sys.argv[3])
- rel_tol = float(sys.argv[4])
+if len(sys.argv) == 4:
+ fileTol = sys.argv[3]
-print "Absolute difference tolerance:", abs_tol
-print "Relative difference tolerance:", rel_tol
-
data1 = loadtxt(file1)
data2 = loadtxt(file2)
+if fileTol != "":
+ dataTol = loadtxt(fileTol)
+ n_linesTol = shape(dataTol)[0]
+ n_columnsTol = shape(dataTol)[1]
+
if shape(data1) != shape(data2):
print "Non-matching array sizes in the files", file1, "and", file2, "."
@@ -47,6 +51,22 @@ if shape(data1) != shape(data2):
n_lines = shape(data1)[0]
n_columns = shape(data1)[1]
+if fileTol != "":
+ if n_linesTol != 2:
+ print "Incorrect number of lines in tolerance file '%s'."%fileTol
+ if n_columnsTol != n_columns:
+ print "Incorrect number of columns in tolerance file '%s'."%fileTol
+
+if fileTol == "":
+ print "Absolute difference tolerance:", abs_tol
+ print "Relative difference tolerance:", rel_tol
+ absTol = ones(n_columns) * abs_tol
+ relTol = ones(n_columns) * rel_tol
+else:
+ print "Tolerances read from file"
+ absTol = dataTol[0,:]
+ relTol = dataTol[1,:]
+
error = False
for i in range(n_lines):
for j in range(n_columns):
@@ -58,17 +78,17 @@ for i in range(n_lines):
rel_diff = abs(data1[i,j] - data2[i,j]) / sum
else:
rel_diff = 0.
-
- if( abs_diff > abs_tol):
- print "Absolute difference larger than tolerance (%e) on line %d, column %d:"%(abs_tol, i,j)
+
+ if( abs_diff > absTol[j]):
+ print "Absolute difference larger than tolerance (%e) on line %d, column %d:"%(absTol[j], i,j)
print "%10s: a = %e"%("File 1", data1[i,j])
print "%10s: b = %e"%("File 2", data2[i,j])
print "%10s: |a-b| = %e"%("Difference", abs_diff)
print ""
error = True
- if( rel_diff > rel_tol):
- print "Relative difference larger than tolerance (%e) on line %d, column %d:"%(rel_tol, i,j)
+ if( rel_diff > relTol[j]):
+ print "Relative difference larger than tolerance (%e) on line %d, column %d:"%(relTol[j], i,j)
print "%10s: a = %e"%("File 1", data1[i,j])
print "%10s: b = %e"%("File 2", data2[i,j])
print "%10s: |a-b|/|a+b| = %e"%("Difference", rel_diff)
diff --git a/tests/test27cells.c b/tests/test27cells.c
index ad0391c981abca599f506590e75bc58b4c433e59..74c38996a81056b10633bf2bbf18cc7cff7e8f0d 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -28,7 +28,7 @@
* Returns a random number (uniformly distributed) in [a,b[
*/
double random_uniform(double a, double b) {
- return (rand() / (double)RAND_MAX) * (a - b) + a;
+ return (rand() / (double)RAND_MAX) * (b - a) + a;
}
/* n is both particles per axis and box size:
@@ -52,7 +52,6 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
for (size_t x = 0; x < n; ++x) {
for (size_t y = 0; y < n; ++y) {
for (size_t z = 0; z < n; ++z) {
- // Add .5 for symmetry: 0.5, 1.5, 2.5 vs. 0, 1, 2
part->x[0] =
offset[0] +
size * (x + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
@@ -62,9 +61,12 @@ 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] = 1. * random_uniform(-0.1, 0.1);
- part->v[1] = 1. * random_uniform(-0.1, 0.1);
- part->v[2] = 1. * random_uniform(-0.1, 0.1);
+ // 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);
part->h = size * h / (float)n;
part->id = ++(*partId);
part->mass = density * volume / count;
@@ -134,41 +136,68 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
FILE *file = fopen(fileName, "w");
+ /* Write header */
fprintf(file,
- "# ID pos:[x y z] rho rho_dh wcount wcount_dh div_v curl_v:[x "
- "y z]\n");
+ "# %4s %10s %10s %10s %10s %10s %10s %13s %13s %13s %13s %13s "
+ "%13s %13s %13s\n",
+ "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");
- fprintf(file, "# -----------------------------------\n");
+ fprintf(file, "# Main cell --------------------------------------------\n");
+ /* Write main cell */
for (size_t pid = 0; pid < main_cell->count; pid++) {
- fprintf(file, "%6llu %f %f %f %f %f %f %f %f %f %f %f\n",
+ fprintf(file,
+ "%6llu %10f %10f %10f %10f %10f %10f %13e %13e %13e %13e %13e "
+ "%13e %13e %13e\n",
main_cell->parts[pid].id, main_cell->parts[pid].x[0],
main_cell->parts[pid].x[1], main_cell->parts[pid].x[2],
- main_cell->parts[pid].rho, main_cell->parts[pid].rho_dh,
- main_cell->parts[pid].density.wcount,
+ main_cell->parts[pid].v[0], main_cell->parts[pid].v[1],
+ main_cell->parts[pid].v[2], main_cell->parts[pid].rho,
+ main_cell->parts[pid].rho_dh, main_cell->parts[pid].density.wcount,
main_cell->parts[pid].density.wcount_dh,
+#ifdef GADGET2_SPH
main_cell->parts[pid].div_v, main_cell->parts[pid].density.rot_v[0],
main_cell->parts[pid].density.rot_v[1],
- main_cell->parts[pid].density.rot_v[2]);
+ main_cell->parts[pid].density.rot_v[2]
+#else
+ 0., 0., 0., 0.
+#endif
+ );
}
- for (int j = 0; j < 27; ++j) {
+ /* Write all other cells */
+ for (int i = 0; i < 3; ++i) {
+ for (int j = 0; j < 3; ++j) {
+ for (int k = 0; k < 3; ++k) {
- struct cell *cj = cells[j];
- if (cj == main_cell) continue;
+ struct cell *cj = cells[i * 9 + j * 3 + k];
+ if (cj == main_cell) continue;
- fprintf(file, "# -----------------------------------\n");
+ fprintf(file,
+ "# Offset: [%2d %2d %2d] -----------------------------------\n",
+ i - 1, j - 1, k - 1);
- for (size_t pjd = 0; pjd < cj->count; pjd++) {
- fprintf(file, "%6llu %f %f %f %f %f %f %f %f %f %f %f\n",
+ for (size_t pjd = 0; pjd < cj->count; pjd++) {
+ fprintf(
+ file,
+ "%6llu %10f %10f %10f %10f %10f %10f %13e %13e %13e %13e %13e "
+ "%13e %13e %13e\n",
cj->parts[pjd].id, cj->parts[pjd].x[0], cj->parts[pjd].x[1],
- cj->parts[pjd].x[2], cj->parts[pjd].rho, cj->parts[pjd].rho_dh,
+ cj->parts[pjd].x[2], cj->parts[pjd].v[0], cj->parts[pjd].v[1],
+ cj->parts[pjd].v[2], cj->parts[pjd].rho, cj->parts[pjd].rho_dh,
cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh,
+#ifdef GADGET2_SPH
cj->parts[pjd].div_v, cj->parts[pjd].density.rot_v[0],
- cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2]);
+ cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2]
+#else
+ 0., 0., 0., 0.
+#endif
+ );
+ }
+ }
}
}
-
fclose(file);
}
@@ -180,7 +209,7 @@ 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.1255, size = 1., rho = 1.;
+ double h = 1.12575, size = 1., rho = 1.;
double perturbation = 0.;
char outputFileNameExtension[200] = "";
char outputFileName[200] = "";
@@ -229,14 +258,18 @@ int main(int argc, char *argv[]) {
"\nThese are then interacted using runner_dopair1_density."
"\n\nOptions:"
"\n-h DISTANCE=1.1255 - Smoothing length"
- "\n-m rho - Physical density in the cell"
- "\n-s size - Physical size of the cell"
+ "\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-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);
+
/* Build the infrastructure */
struct space space;
space.periodic = 0;
@@ -259,13 +292,13 @@ 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);
}
}
}
+ /* Store the main cell for future use */
main_cell = cells[13];
ticks time = 0;
diff --git a/tests/test27cells.sh b/tests/test27cells.sh
index 03fb5e093525c22495d61997a08f038af99c4675..09d2513bd3ef404c7bf434948af7f10306c98ede 100755
--- a/tests/test27cells.sh
+++ b/tests/test27cells.sh
@@ -3,6 +3,6 @@ rm brute_force_27_standard.dat swift_dopair_27_standard.dat
./test27cells -p 6 -r 1 -d 0 -f standard
-python difffloat.py brute_force_27_standard.dat swift_dopair_27_standard.dat 1e-5 5e-6
+python difffloat.py brute_force_27_standard.dat swift_dopair_27_standard.dat tolerance.dat
exit $?
diff --git a/tests/test27cellsPerturbed.sh b/tests/test27cellsPerturbed.sh
index 0d2f6d4762386aa53d2321fdaff396454fb2ff1f..73d2933984d38f7dcc992f07ec2e016f3544b636 100755
--- a/tests/test27cellsPerturbed.sh
+++ b/tests/test27cellsPerturbed.sh
@@ -3,6 +3,6 @@ rm brute_force_27_perturbed.dat swift_dopair_27_perturbed.dat
./test27cells -p 6 -r 1 -d 0.1 -f perturbed
-python difffloat.py brute_force_27_perturbed.dat swift_dopair_27_perturbed.dat 1e-5 5e-6
+python difffloat.py brute_force_27_perturbed.dat swift_dopair_27_perturbed.dat tolerance.dat
exit $?
diff --git a/tests/testPair.c b/tests/testPair.c
index 38b4f37260c41956236915a85244050ee281699b..23ce4eb3de460f4e17b7b6f81cb39a628f3d100f 100644
--- a/tests/testPair.c
+++ b/tests/testPair.c
@@ -28,40 +28,49 @@
* Returns a random number (uniformly distributed) in [a,b[
*/
double random_uniform(double a, double b) {
- return (rand() / (double)RAND_MAX) * (a - b) + a;
+ 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
*/
-struct cell *make_cell(size_t n, double *offset, double h,
- unsigned long long *partId, double pert) {
- size_t count = n * n * n;
+struct cell *make_cell(size_t n, double *offset, double size, double h,
+ double density, unsigned long long *partId,
+ double pert) {
+ const size_t count = n * n * n;
+ const double volume = size * size * size;
struct cell *cell = malloc(sizeof(struct cell));
bzero(cell, sizeof(struct cell));
- struct part *part;
- size_t x, y, z, size;
- size = count * sizeof(struct part);
- if (posix_memalign((void **)&cell->parts, part_align, size) != 0) {
+ if (posix_memalign((void **)&cell->parts, part_align,
+ count * sizeof(struct part)) != 0) {
error("couldn't allocate particles, no. of particles: %d", (int)count);
}
bzero(cell->parts, count * sizeof(struct part));
- part = cell->parts;
- for (x = 0; x < n; ++x) {
- for (y = 0; y < n; ++y) {
- for (z = 0; z < n; ++z) {
- // Add .5 for symmetry: 0.5, 1.5, 2.5 vs. 0, 1, 2
- part->x[0] = x + offset[0] + 0.5 + random_uniform(-0.5, 0.5) * pert;
- part->x[1] = y + offset[1] + 0.5 + random_uniform(-0.5, 0.5) * pert;
- part->x[2] = z + offset[2] + 0.5 + random_uniform(-0.5, 0.5) * pert;
- part->v[0] = 0.0f;
- part->v[1] = 0.0f;
- part->v[2] = 0.0f;
- part->h = h;
+ /* Construct the parts */
+ struct part *part = cell->parts;
+ for (size_t x = 0; x < n; ++x) {
+ for (size_t y = 0; y < n; ++y) {
+ for (size_t z = 0; z < n; ++z) {
+ part->x[0] =
+ offset[0] +
+ size * (x + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
+ part->x[1] =
+ offset[1] +
+ size * (y + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
+ 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);
+ part->h = size * h / (float)n;
part->id = ++(*partId);
- part->mass = 1.0f;
+ part->mass = density * volume / count;
part->ti_begin = 0;
part->ti_end = 1;
++part;
@@ -69,6 +78,7 @@ struct cell *make_cell(size_t n, double *offset, double h,
}
}
+ /* Cell properties */
cell->split = 0;
cell->h_max = h;
cell->count = count;
@@ -116,28 +126,49 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
FILE *file = fopen(fileName, "w");
+ /* Write header */
fprintf(file,
- "# ID pos:[x y z] rho rho_dh wcount wcount_dh div_v curl_v:[x "
- "y z]\n");
+ "# %4s %10s %10s %10s %10s %10s %10s %13s %13s %13s %13s %13s "
+ "%13s %13s %13s\n",
+ "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");
+
+ fprintf(file, "# ci --------------------------------------------\n");
for (size_t pid = 0; pid < ci->count; pid++) {
- fprintf(file, "%6llu %f %f %f %f %f %f %f %f %f %f %f\n", ci->parts[pid].id,
- ci->parts[pid].x[0], ci->parts[pid].x[1], ci->parts[pid].x[2],
- ci->parts[pid].rho, ci->parts[pid].rho_dh,
+ fprintf(file,
+ "%6llu %10f %10f %10f %10f %10f %10f %13e %13e %13e %13e %13e "
+ "%13e %13e %13e\n",
+ ci->parts[pid].id, ci->parts[pid].x[0], ci->parts[pid].x[1],
+ ci->parts[pid].x[2], ci->parts[pid].v[0], ci->parts[pid].v[1],
+ ci->parts[pid].v[2], ci->parts[pid].rho, ci->parts[pid].rho_dh,
ci->parts[pid].density.wcount, ci->parts[pid].density.wcount_dh,
+#ifdef GADGET2_SPH
ci->parts[pid].div_v, ci->parts[pid].density.rot_v[0],
- ci->parts[pid].density.rot_v[1], ci->parts[pid].density.rot_v[2]);
+ ci->parts[pid].density.rot_v[1], ci->parts[pid].density.rot_v[2]
+#else
+ 0., 0., 0., 0.
+#endif
+ );
}
- fprintf(file, "# -----------------------------------\n");
+ fprintf(file, "# cj --------------------------------------------\n");
for (size_t pjd = 0; pjd < cj->count; pjd++) {
- fprintf(file, "%6llu %f %f %f %f %f %f %f %f %f %f %f\n", cj->parts[pjd].id,
- cj->parts[pjd].x[0], cj->parts[pjd].x[1], cj->parts[pjd].x[2],
- cj->parts[pjd].rho, cj->parts[pjd].rho_dh,
+ fprintf(file,
+ "%6llu %10f %10f %10f %10f %10f %10f %13e %13e %13e %13e %13e "
+ "%13e %13e %13e\n",
+ cj->parts[pjd].id, cj->parts[pjd].x[0], cj->parts[pjd].x[1],
+ cj->parts[pjd].x[2], cj->parts[pjd].v[0], cj->parts[pjd].v[1],
+ cj->parts[pjd].v[2], cj->parts[pjd].rho, cj->parts[pjd].rho_dh,
cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh,
+#ifdef GADGET2_SPH
cj->parts[pjd].div_v, cj->parts[pjd].density.rot_v[0],
- cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2]);
+ cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2]
+#else
+ 0., 0., 0., 0.
+#endif
+ );
}
fclose(file);
@@ -148,7 +179,7 @@ void runner_dopair1_density(struct runner *r, struct cell *ci, struct cell *cj);
int main(int argc, char *argv[]) {
size_t particles = 0, runs = 0, volume, type = 0;
- double offset[3] = {0, 0, 0}, h = 1.1255;
+ double offset[3] = {0, 0, 0}, h = 1.1255, size = 1., rho = 1.;
double perturbation = 0.1;
struct cell *ci, *cj;
struct space space;
@@ -217,9 +248,9 @@ int main(int argc, char *argv[]) {
volume = particles * particles * particles;
message("particles: %zu B\npositions: 0 B", 2 * volume * sizeof(struct part));
- ci = make_cell(particles, offset, h, &partId, perturbation);
- for (size_t i = 0; i < type + 1; ++i) offset[i] = particles;
- cj = make_cell(particles, offset, h, &partId, perturbation);
+ ci = make_cell(particles, offset, size, h, rho, &partId, perturbation);
+ for (size_t i = 0; i < type + 1; ++i) offset[i] = 1.;
+ cj = make_cell(particles, offset, size, h, rho, &partId, perturbation);
time = 0;
for (size_t i = 0; i < runs; ++i) {
diff --git a/tests/testPair.sh b/tests/testPair.sh
index d39ad74e3f6b306dd12026fc095841d085c583c0..f6f505e56a2c7a5c3cff0ec04bd871278634193c 100755
--- a/tests/testPair.sh
+++ b/tests/testPair.sh
@@ -3,6 +3,6 @@ rm brute_force_standard.dat swift_dopair_standard.dat
./testPair -p 6 -r 1 -d 0 -f standard
-python difffloat.py brute_force_standard.dat swift_dopair_standard.dat 1e-5 5e-6
+python difffloat.py brute_force_standard.dat swift_dopair_standard.dat tolerance.dat
exit $?
diff --git a/tests/testPairPerturbed.sh b/tests/testPairPerturbed.sh
index c3c6fc82eb482821ff6dbb3ea4e0d640e0e86f49..544ba1b032da8426c065dcfb2ce3ee554c5e76a1 100755
--- a/tests/testPairPerturbed.sh
+++ b/tests/testPairPerturbed.sh
@@ -3,6 +3,6 @@ rm brute_force_perturbed.dat swift_dopair_perturbed.dat
./testPair -p 6 -r 1 -d 0.1 -f perturbed
-python difffloat.py brute_force_perturbed.dat swift_dopair_perturbed.dat 1e-5 5e-6
+python difffloat.py brute_force_perturbed.dat swift_dopair_perturbed.dat tolerance.dat
exit $?
diff --git a/tests/testSPHStep.c b/tests/testSPHStep.c
index 984b8ea867250d0bda1bc14d2600279a27321b2c..223078ecb637e64d94e37cdf8c0f60a86bdd5ff7 100644
--- a/tests/testSPHStep.c
+++ b/tests/testSPHStep.c
@@ -77,6 +77,10 @@ struct cell *make_cell(size_t N, float cellSize, int offset[3], int id_offset) {
#ifdef DEFAULT_SPH
+/* Just a forward declaration... */
+void runner_doself1_density(struct runner *r, struct cell *ci);
+void runner_doself2_force(struct runner *r, struct cell *ci);
+
/* Run a full time step integration for one cell */
int main() {
@@ -132,7 +136,7 @@ int main() {
/* Initialise the particles */
for (j = 0; j < 27; ++j) {
- runner_doinit(&r, cells[j]);
+ runner_doinit(&r, cells[j], 0);
}
/* Compute density */
@@ -145,7 +149,7 @@ int main() {
runner_doself2_force(&r, ci);
runner_dokick(&r, ci, 1);
- message("t_end=%f", p->t_end);
+ message("ti_end=%d", p->ti_end);
free(ci->parts);
free(ci->xparts);
diff --git a/tests/testSingle.c b/tests/testSingle.c
index c85b77ff1c5b2285c33fa7787bbd53deab463039..8771fba0c1912905d3936562fa9dad0223d89220 100644
--- a/tests/testSingle.c
+++ b/tests/testSingle.c
@@ -91,8 +91,8 @@ int main(int argc, char *argv[]) {
p2.force.POrho2 = p2.u * (const_hydro_gamma - 1.0f) / p2.rho;
/* Dump a header. */
- printParticle_single(&p1);
- printParticle_single(&p2);
+ //printParticle_single(&p1, NULL);
+ //printParticle_single(&p2, NULL);
printf("# r a_1 udt_1 a_2 udt_2\n");
/* Loop over the different radii. */
@@ -103,9 +103,9 @@ int main(int argc, char *argv[]) {
r2 = dx[0] * dx[0];
/* Clear the particle fields. */
- p1.a[0] = 0.0f;
+ p1.a_hydro[0] = 0.0f;
p1.force.u_dt = 0.0f;
- p2.a[0] = 0.0f;
+ p2.a_hydro[0] = 0.0f;
p2.force.u_dt = 0.0f;
/* Interact the particles. */
@@ -130,8 +130,8 @@ int main(int argc, char *argv[]) {
/* Output the results. */
printf(
- "%.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e\n", -dx[0], p1.a[0],
- p1.a[1], p1.a[2], p1.force.u_dt,
+ "%.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e %.3e\n", -dx[0],
+ p1.a_hydro[0], p1.a_hydro[1], p1.a_hydro[2], p1.force.u_dt,
/// -dx[0] , p1.rho , p1.density.wcount , p2.rho , p2.density.wcount ,
w, dwdx, gradw[0], gradw[1], gradw[2]);
diff --git a/tests/tolerance.dat b/tests/tolerance.dat
new file mode 100644
index 0000000000000000000000000000000000000000..48de4383eab6214812183be25d3036a324ccbc27
--- /dev/null
+++ b/tests/tolerance.dat
@@ -0,0 +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