Commit 8a357b65 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'updated_vectorisation_tests' into 'master'

Updated vectorisation tests



See merge request !134
parents 2cb4ba76 728d7eec
......@@ -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);
......
......@@ -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 */
......@@ -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
......@@ -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)
......
......@@ -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;
......
......@@ -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 $?
......@@ -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 $?
......@@ -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) {
......
......@@ -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 $?
......@@ -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 $?
......@@ -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