diff --git a/.gitignore b/.gitignore
index 8137ea759b24b3f4ec9909a460da4bcb47b0a1ac..ddeaaf8b235dc217fbfb6559e66bd665d1f31745 100644
--- a/.gitignore
+++ b/.gitignore
@@ -25,9 +25,16 @@ examples/swift_mindt
examples/swift_mindt_mpi
examples/swift_mpi
-tests/testVectorize
-tests/brute_force.dat
-tests/swift_dopair.dat
+tests/testPair
+tests/brute_force_standard.dat
+tests/swift_dopair_standard.dat
+tests/brute_force_perturbed.dat
+tests/swift_dopair_perturbed.dat
+tests/test27cells
+tests/brute_force_27_standard.dat
+tests/swift_dopair_27_standard.dat
+tests/brute_force_27_perturbed.dat
+tests/swift_dopair_27_perturbed.dat
tests/testGreetings
tests/testReading
tests/input.hdf5
diff --git a/src/runner_doiact.h b/src/runner_doiact.h
index cf5d56e94169b44e6cd2974a3422a0bc5e4610ac..de339db6133fcc829bdc6ee0ce9e537b68982422 100644
--- a/src/runner_doiact.h
+++ b/src/runner_doiact.h
@@ -1235,7 +1235,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
#else
/* Does pi need to be updated too? */
- if (pi->dt <= dt_step) {
+ if (pi->ti_end <= ti_current) {
/* Add this interaction to the symmetric queue. */
r2q2[icount2] = r2;
diff --git a/src/space.h b/src/space.h
index 2f54b5d74ce961c1486929698eb0e8f50b72c748..e761595838ae78b0d8a67cca676cfa59f3f700f6 100644
--- a/src/space.h
+++ b/src/space.h
@@ -64,9 +64,6 @@ struct space {
/* The minimum and maximum cutoff radii. */
double h_max, cell_min;
- /* Current time step for particles. */
- float dt_step;
-
/* Current maximum displacement for particles. */
float dx_max;
diff --git a/src/tools.c b/src/tools.c
index 5feba7759f730faea1f38ceb9835f2076bc37a56..a384974fdc94452079838ae0eebf30e1815f644b 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -236,6 +236,53 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
}
}
+void self_all_density(struct runner *r, struct cell *ci) {
+ float r2, hi, hj, hig2, hjg2, dxi[3]; //, dxj[3];
+ struct part *pi, *pj;
+
+ /* Implements a double-for loop and checks every interaction */
+ for (int i = 0; i < ci->count; ++i) {
+
+ pi = &ci->parts[i];
+ hi = pi->h;
+ hig2 = hi * hi * kernel_gamma2;
+
+ for (int j = i + 1; j < ci->count; ++j) {
+
+ pj = &ci->parts[j];
+ hj = pj->h;
+ hjg2 = hj * hj * kernel_gamma2;
+
+ if (pi == pj) continue;
+
+ /* Pairwise distance */
+ r2 = 0.0f;
+ for (int k = 0; k < 3; k++) {
+ dxi[k] = ci->parts[i].x[k] - ci->parts[j].x[k];
+ r2 += dxi[k] * dxi[k];
+ }
+
+ /* Hit or miss? */
+ if (r2 < hig2) {
+
+ /* Interact */
+ runner_iact_nonsym_density(r2, dxi, hi, hj, pi, pj);
+ }
+
+ /* Hit or miss? */
+ if (r2 < hjg2) {
+
+ dxi[0] = -dxi[0];
+ dxi[1] = -dxi[1];
+ dxi[2] = -dxi[2];
+
+ /* Interact */
+ runner_iact_nonsym_density(r2, dxi, hj, hi, pj, pi);
+ }
+ }
+ }
+}
+
void pairs_single_grav(double *dim, long long int pid,
struct gpart *__restrict__ parts, int N, int periodic) {
diff --git a/src/tools.h b/src/tools.h
index 59646291bda46a7dd0f5a34e158e3e0a6f21d3ca..ccffc77ceb8a967fd40c3737651ba75d529eee0f 100644
--- a/src/tools.h
+++ b/src/tools.h
@@ -33,6 +33,7 @@ void pairs_single_density(double *dim, long long int pid,
struct part *__restrict__ parts, int N, int periodic);
void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj);
+void self_all_density(struct runner *r, struct cell *ci);
void pairs_n2(double *dim, struct part *__restrict__ parts, int N,
int periodic);
diff --git a/tests/Makefile.am b/tests/Makefile.am
index f0bfbefd3c7f4591134d1707c4ac9bf63278e855..a35344790967d3cc3f5d299c6f9af3cb4c2fab43 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -21,10 +21,12 @@ AM_CFLAGS = -I../src $(HDF5_CPPFLAGS) -DTIMER
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 testTimeIntegration
+TESTS = testGreetings testReading.sh testSingle testPair.sh testPairPerturbed.sh \
+ test27cells.sh test27cellsPerturbed.sh
# List of test programs to compile
-check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration testSPHStep testVectorize
+check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
+ testSPHStep testPair test27cells
# Sources for the individual programs
testGreetings_SOURCES = testGreetings.c
@@ -37,7 +39,10 @@ testSPHStep_SOURCES = testSPHStep.c
testSingle_SOURCES = testSingle.c
-testVectorize_SOURCES = testVectorize.c
+testPair_SOURCES = testPair.c
+
+test27cells_SOURCES = test27cells.c
# Files necessary for distribution
-EXTRA_DIST = testReading.sh makeInput.py
+EXTRA_DIST = testReading.sh makeInput.py testPair.sh testPairPerturbed.sh \
+ test27cells.sh test27cellsPerturbed.sh
diff --git a/tests/difffloat.py b/tests/difffloat.py
new file mode 100644
index 0000000000000000000000000000000000000000..fe1118dfa050d3507137febf790eb5a698cfdb57
--- /dev/null
+++ b/tests/difffloat.py
@@ -0,0 +1,83 @@
+###############################################################################
+ # 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/>.
+ #
+ ##############################################################################
+
+from numpy import *
+import sys
+
+abs_tol = 1e-7
+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
+
+file1 = sys.argv[1]
+file2 = sys.argv[2]
+
+if len(sys.argv) >= 5:
+ abs_tol = float(sys.argv[3])
+ rel_tol = float(sys.argv[4])
+
+print "Absolute difference tolerance:", abs_tol
+print "Relative difference tolerance:", rel_tol
+
+data1 = loadtxt(file1)
+data2 = loadtxt(file2)
+
+if shape(data1) != shape(data2):
+ print "Non-matching array sizes in the files", file1, "and", file2, "."
+ sys.exit(1)
+
+n_lines = shape(data1)[0]
+n_columns = shape(data1)[1]
+
+error = False
+for i in range(n_lines):
+ for j in range(n_columns):
+
+ abs_diff = abs(data1[i,j] - data2[i,j])
+
+ sum = abs(data1[i,j] + data2[i,j])
+ if sum > 0:
+ 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)
+ 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)
+ 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)
+ print ""
+ error = True
+
+
+if error:
+ exit(1)
+else:
+ print "No differences found"
+ exit(0)
diff --git a/tests/test27cells.c b/tests/test27cells.c
new file mode 100644
index 0000000000000000000000000000000000000000..ad0391c981abca599f506590e75bc58b4c433e59
--- /dev/null
+++ b/tests/test27cells.c
@@ -0,0 +1,334 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+#include <fenv.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdio.h>
+#include <unistd.h>
+#include "swift.h"
+
+/**
+ * Returns a random number (uniformly distributed) in [a,b[
+ */
+double random_uniform(double a, double b) {
+ return (rand() / (double)RAND_MAX) * (a - b) + 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 size, double h,
+ double density, 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));
+
+ 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));
+
+ /* 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) {
+ // 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;
+ 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] = 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->h = size * h / (float)n;
+ part->id = ++(*partId);
+ part->mass = density * volume / count;
+ part->ti_begin = 0;
+ part->ti_end = 1;
+ ++part;
+ }
+ }
+ }
+
+ /* Cell properties */
+ cell->split = 0;
+ cell->h_max = h;
+ cell->count = count;
+ cell->dx_max = 0.;
+ cell->h[0] = size;
+ cell->h[1] = size;
+ cell->h[2] = size;
+ cell->loc[0] = offset[0];
+ cell->loc[1] = offset[1];
+ cell->loc[2] = offset[2];
+
+ cell->ti_end_min = 1;
+ cell->ti_end_max = 1;
+
+ cell->sorted = 0;
+ cell->sort = NULL;
+ cell->sortsize = 0;
+ runner_dosort(NULL, cell, 0x1FFF, 0);
+
+ return cell;
+}
+
+void clean_up(struct cell *ci) {
+ free(ci->parts);
+ free(ci->sort);
+ free(ci);
+}
+
+/**
+ * @brief Initializes all particles field to be ready for a density calculation
+ */
+void zero_particle_fields(struct cell *c) {
+
+ for (size_t pid = 0; pid < c->count; pid++) {
+ c->parts[pid].rho = 0.f;
+ c->parts[pid].rho_dh = 0.f;
+ hydro_init_part(&c->parts[pid]);
+ }
+}
+
+/**
+ * @brief Ends the loop by adding the appropriate coefficients
+ */
+void end_calculation(struct cell *c) {
+
+ for (size_t pid = 0; pid < c->count; pid++) {
+ hydro_end_density(&c->parts[pid], 1);
+ }
+}
+
+/**
+ * @brief Dump all the particles to a file
+ */
+void dump_particle_fields(char *fileName, struct cell *main_cell,
+ struct cell **cells) {
+
+ FILE *file = fopen(fileName, "w");
+
+ fprintf(file,
+ "# ID pos:[x y z] rho rho_dh wcount wcount_dh div_v curl_v:[x "
+ "y z]\n");
+
+ fprintf(file, "# -----------------------------------\n");
+
+ 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",
+ 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].density.wcount_dh,
+ 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]);
+ }
+
+ for (int j = 0; j < 27; ++j) {
+
+ struct cell *cj = cells[j];
+ if (cj == main_cell) continue;
+
+ fprintf(file, "# -----------------------------------\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,
+ cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh,
+ 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]);
+ }
+ }
+
+ fclose(file);
+}
+
+/* Just a forward declaration... */
+void runner_dopair1_density(struct runner *r, struct cell *ci, struct cell *cj);
+void runner_doself1_density(struct runner *r, struct cell *ci);
+
+/* And go... */
+int main(int argc, char *argv[]) {
+
+ size_t runs = 0, particles = 0;
+ double h = 1.1255, size = 1., rho = 1.;
+ double perturbation = 0.;
+ char outputFileNameExtension[200] = "";
+ char outputFileName[200] = "";
+
+ /* Initialize CPU frequency, this also starts time. */
+ unsigned long long cpufreq = 0;
+ clocks_set_cpufreq(cpufreq);
+
+ /* Get some randomness going */
+ srand(0);
+
+ char c;
+ while ((c = getopt(argc, argv, "m:s:h:p:r:t:d:f:")) != -1) {
+ switch (c) {
+ case 'h':
+ sscanf(optarg, "%lf", &h);
+ break;
+ case 's':
+ sscanf(optarg, "%lf", &size);
+ break;
+ case 'p':
+ sscanf(optarg, "%zu", &particles);
+ break;
+ case 'r':
+ sscanf(optarg, "%zu", &runs);
+ break;
+ case 'd':
+ sscanf(optarg, "%lf", &perturbation);
+ break;
+ case 'm':
+ sscanf(optarg, "%lf", &rho);
+ break;
+ case 'f':
+ strcpy(outputFileNameExtension, optarg);
+ break;
+ case '?':
+ error("Unknown option.");
+ break;
+ }
+ }
+
+ if (h < 0 || particles == 0 || runs == 0) {
+ printf(
+ "\nUsage: %s -p PARTICLES_PER_AXIS -r NUMBER_OF_RUNS [OPTIONS...]\n"
+ "\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-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);
+ }
+
+ /* Build the infrastructure */
+ struct space space;
+ space.periodic = 0;
+ space.h_max = h;
+
+ struct engine engine;
+ engine.s = &space;
+ engine.time = 0.1f;
+ engine.ti_current = 1;
+
+ struct runner runner;
+ runner.e = &engine;
+
+ /* Construct some cells */
+ struct cell *cells[27];
+ struct cell *main_cell;
+ static long long partId = 0;
+ for (int i = 0; i < 3; ++i) {
+ for (int j = 0; j < 3; ++j) {
+ 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);
+ }
+ }
+ }
+
+ main_cell = cells[13];
+
+ ticks time = 0;
+ for (size_t i = 0; i < runs; ++i) {
+
+ /* Zero the fields */
+ for (int j = 0; j < 27; ++j) zero_particle_fields(cells[j]);
+
+ const ticks tic = getticks();
+
+ /* Run all the pairs */
+ for (int j = 0; j < 27; ++j)
+ if (cells[j] != main_cell)
+ runner_dopair1_density(&runner, main_cell, cells[j]);
+
+ /* And now the self-interaction */
+ runner_doself1_density(&runner, main_cell);
+
+ const ticks toc = getticks();
+ time += toc - tic;
+
+ /* Let's get physical ! */
+ end_calculation(main_cell);
+
+ /* Dump if necessary */
+ if (i % 50 == 0) {
+ sprintf(outputFileName, "swift_dopair_27_%s.dat",
+ outputFileNameExtension);
+ dump_particle_fields(outputFileName, main_cell, cells);
+ }
+ }
+
+ /* Output timing */
+ message("SWIFT calculation took : %15lli ticks.", time / runs);
+
+ /* Now perform a brute-force version for accuracy tests */
+
+ /* Zero the fields */
+ for (int i = 0; i < 27; ++i) zero_particle_fields(cells[i]);
+
+ const ticks tic = getticks();
+
+ /* Run all the brute-force pairs */
+ for (int j = 0; j < 27; ++j)
+ if (cells[j] != main_cell) pairs_all_density(&runner, main_cell, cells[j]);
+
+ /* And now the self-interaction */
+ self_all_density(&runner, main_cell);
+
+ const ticks toc = getticks();
+
+ /* Let's get physical ! */
+ end_calculation(main_cell);
+
+ /* Dump */
+ sprintf(outputFileName, "brute_force_27_%s.dat", outputFileNameExtension);
+ dump_particle_fields(outputFileName, main_cell, cells);
+
+ /* Output timing */
+ message("Brute force calculation took : %15lli ticks.", toc - tic);
+
+ /* Clean things to make the sanitizer happy ... */
+ for (int i = 0; i < 27; ++i) clean_up(cells[i]);
+
+ return 0;
+}
diff --git a/tests/test27cells.sh b/tests/test27cells.sh
new file mode 100755
index 0000000000000000000000000000000000000000..03fb5e093525c22495d61997a08f038af99c4675
--- /dev/null
+++ b/tests/test27cells.sh
@@ -0,0 +1,8 @@
+#!/bin/bash
+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
+
+exit $?
diff --git a/tests/test27cellsPerturbed.sh b/tests/test27cellsPerturbed.sh
new file mode 100755
index 0000000000000000000000000000000000000000..0d2f6d4762386aa53d2321fdaff396454fb2ff1f
--- /dev/null
+++ b/tests/test27cellsPerturbed.sh
@@ -0,0 +1,8 @@
+#!/bin/bash
+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
+
+exit $?
diff --git a/tests/testVectorize.c b/tests/testPair.c
similarity index 56%
rename from tests/testVectorize.c
rename to tests/testPair.c
index a18b6e8af5ac3f7b94bd7be3bdf8fd21e49681ff..38b4f37260c41956236915a85244050ee281699b 100644
--- a/tests/testVectorize.c
+++ b/tests/testPair.c
@@ -1,3 +1,22 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
#include <fenv.h>
#include <stdlib.h>
#include <string.h>
@@ -5,13 +24,21 @@
#include <unistd.h>
#include "swift.h"
+/**
+ * Returns a random number (uniformly distributed) in [a,b[
+ */
+double random_uniform(double a, double b) {
+ return (rand() / (double)RAND_MAX) * (a - b) + 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) {
+ unsigned long long *partId, double pert) {
size_t count = n * n * n;
- struct cell *cell = malloc(sizeof *cell);
+ struct cell *cell = malloc(sizeof(struct cell));
+ bzero(cell, sizeof(struct cell));
struct part *part;
size_t x, y, z, size;
@@ -19,18 +46,19 @@ struct cell *make_cell(size_t n, double *offset, double h,
if (posix_memalign((void **)&cell->parts, part_align, size) != 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;
- part->x[1] = y + offset[1] + 0.5;
- part->x[2] = z + offset[2] + 0.5;
- part->v[0] = 1.0f;
- part->v[1] = 1.0f;
- part->v[2] = 1.0f;
+ 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;
part->id = ++(*partId);
part->mass = 1.0f;
@@ -44,12 +72,20 @@ struct cell *make_cell(size_t n, double *offset, double h,
cell->split = 0;
cell->h_max = h;
cell->count = count;
- cell->dx_max = 1.;
+ cell->dx_max = 0.;
cell->h[0] = n;
cell->h[1] = n;
cell->h[2] = n;
+ cell->loc[0] = offset[0];
+ cell->loc[1] = offset[1];
+ cell->loc[2] = offset[2];
+
+ cell->ti_end_min = 1;
+ cell->ti_end_max = 1;
- cell->sort = malloc(13 * count * sizeof *cell->sort);
+ cell->sorted = 0;
+ cell->sort = NULL;
+ cell->sortsize = 0;
runner_dosort(NULL, cell, 0x1FFF, 0);
return cell;
@@ -81,10 +117,12 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
FILE *file = fopen(fileName, "w");
fprintf(file,
- "# ID rho rho_dh wcount wcount_dh div_v curl_v:[x y z]\n");
+ "# ID pos:[x y z] rho rho_dh wcount wcount_dh div_v curl_v:[x "
+ "y z]\n");
for (size_t pid = 0; pid < ci->count; pid++) {
- fprintf(file, "%6llu %f %f %f %f %f %f %f %f\n", ci->parts[pid].id,
+ 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,
ci->parts[pid].density.wcount, ci->parts[pid].density.wcount_dh,
ci->parts[pid].div_v, ci->parts[pid].density.rot_v[0],
@@ -94,7 +132,8 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
fprintf(file, "# -----------------------------------\n");
for (size_t pjd = 0; pjd < cj->count; pjd++) {
- fprintf(file, "%6llu %f %f %f %f %f %f %f %f\n", cj->parts[pjd].id,
+ 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,
cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh,
cj->parts[pjd].div_v, cj->parts[pjd].density.rot_v[0],
@@ -109,16 +148,25 @@ 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; // * DIM/PARTS_PER_AXIS == * 1
+ double offset[3] = {0, 0, 0}, h = 1.1255;
+ double perturbation = 0.1;
struct cell *ci, *cj;
struct space space;
struct engine engine;
struct runner runner;
char c;
static unsigned long long partId = 0;
+ char outputFileNameExtension[200] = "";
+ char outputFileName[200] = "";
ticks tic, toc, time;
- while ((c = getopt(argc, argv, "h:p:r:t:")) != -1) {
+ /* Initialize CPU frequency, this also starts time. */
+ unsigned long long cpufreq = 0;
+ clocks_set_cpufreq(cpufreq);
+
+ srand(0);
+
+ while ((c = getopt(argc, argv, "h:p:r:t:d:f:")) != -1) {
switch (c) {
case 'h':
sscanf(optarg, "%lf", &h);
@@ -132,6 +180,15 @@ int main(int argc, char *argv[]) {
case 't':
sscanf(optarg, "%zu", &type);
break;
+ case 'd':
+ sscanf(optarg, "%lf", &perturbation);
+ break;
+ case 'f':
+ strcpy(outputFileNameExtension, optarg);
+ break;
+ case '?':
+ error("Unknown option.");
+ break;
}
}
@@ -142,27 +199,28 @@ int main(int argc, char *argv[]) {
"\nThese are then interacted using runner_dopair1_density."
"\n\nOptions:"
"\n-t TYPE=0 - cells share face (0), edge (1) or corner (2)"
- "\n-h DISTANCE=1.1255 - smoothing length\n",
+ "\n-h DISTANCE=1.1255 - smoothing length"
+ "\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);
}
- volume = particles * particles * particles;
- message("particles: %zu B\npositions: 0 B", 2 * volume * sizeof(struct part));
-
- ci = make_cell(particles, offset, h, &partId);
- for (size_t i = 0; i < type + 1; ++i) offset[i] = particles;
- cj = make_cell(particles, offset, h, &partId);
-
- for (int i = 0; i < 3; ++i) {
- space.h_max = h;
- space.dt_step = 0.1;
- }
+ space.periodic = 0;
+ space.h_max = h;
engine.s = &space;
engine.time = 0.1f;
+ engine.ti_current = 1;
runner.e = &engine;
+ 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);
+
time = 0;
for (size_t i = 0; i < runs; ++i) {
@@ -179,7 +237,10 @@ int main(int argc, char *argv[]) {
time += toc - tic;
/* Dump if necessary */
- if (i % 50 == 0) dump_particle_fields("swift_dopair.dat", ci, cj);
+ if (i % 50 == 0) {
+ sprintf(outputFileName, "swift_dopair_%s.dat", outputFileNameExtension);
+ dump_particle_fields(outputFileName, ci, cj);
+ }
}
/* Output timing */
@@ -193,13 +254,14 @@ int main(int argc, char *argv[]) {
tic = getticks();
- /* Run the test */
+ /* Run the brute-force test */
pairs_all_density(&runner, ci, cj);
toc = getticks();
/* Dump */
- dump_particle_fields("brute_force.dat", ci, cj);
+ sprintf(outputFileName, "brute_force_%s.dat", outputFileNameExtension);
+ dump_particle_fields(outputFileName, ci, cj);
/* Output timing */
message("Brute force calculation took %lli ticks.", toc - tic);
diff --git a/tests/testPair.sh b/tests/testPair.sh
new file mode 100755
index 0000000000000000000000000000000000000000..d39ad74e3f6b306dd12026fc095841d085c583c0
--- /dev/null
+++ b/tests/testPair.sh
@@ -0,0 +1,8 @@
+#!/bin/bash
+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
+
+exit $?
diff --git a/tests/testPairPerturbed.sh b/tests/testPairPerturbed.sh
new file mode 100755
index 0000000000000000000000000000000000000000..c3c6fc82eb482821ff6dbb3ea4e0d640e0e86f49
--- /dev/null
+++ b/tests/testPairPerturbed.sh
@@ -0,0 +1,8 @@
+#!/bin/bash
+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
+
+exit $?
diff --git a/tests/testReading.c b/tests/testReading.c
index d2a2a766171a85ace486914f0f39a987d9d8c3d3..9dda4c7bad75d35a8a93e0c2acb0619409a91afd 100644
--- a/tests/testReading.c
+++ b/tests/testReading.c
@@ -22,7 +22,7 @@
int main() {
- int Ngas = -1, Ngpart = -1;
+ size_t Ngas = 0, Ngpart = 0;
int periodic = -1;
int i, j, k, n;
double dim[3];