Commit 2cb4ba76 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'updated_vectorisation_tests' into 'master'

Updated vectorisation tests



See merge request !133
parents 7ff9b72d a28fec06
......@@ -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
......
......@@ -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;
......
......@@ -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;
......
......@@ -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) {
......
......@@ -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);
......
......@@ -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
###############################################################################
# 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)
/*******************************************************************************
* 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;
}
#!/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 $?
#!/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 $?
/*******************************************************************************
* 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;