Commit 357f87d5 authored by Peter W. Draper's avatar Peter W. Draper
Browse files

Merge branch 'Test-Vectorise' into 'master'

Agreed upon version of the density vectorization test

Wrote a unified version of the vectorization test to be used as a basis for the tests of @jwillis, @rgb and @alepper 

The test also contains a brute-force calculation to check that any approximation made even in the stock DOPAIR1() function is correct.

See merge request !64
parents 0689a0cd 9ef8c274
......@@ -25,6 +25,9 @@ examples/swift_mindt
examples/swift_mindt_mpi
examples/swift_mpi
tests/testVectorize
tests/brute_force.dat
tests/swift_dopair.dat
tests/testGreetings
tests/testReading
tests/input.hdf5
......
......@@ -185,8 +185,8 @@ void engine_redistribute(struct engine *e) {
/* Get the new number of parts for this node, be generous in allocating. */
int nr_parts = 0;
for (int k = 0; k < nr_nodes; k++) nr_parts += counts[k * nr_nodes + nodeID];
struct part *parts_new;
struct xpart *xparts_new, *xparts = s->xparts;
struct part *parts_new = NULL;
struct xpart *xparts_new = NULL, *xparts = s->xparts;
if (posix_memalign((void **)&parts_new, part_align,
sizeof(struct part) * nr_parts * 1.2) != 0 ||
posix_memalign((void **)&xparts_new, part_align,
......@@ -951,8 +951,8 @@ int engine_exchange_strays(struct engine *e, int offset, int *ind, int N) {
message("sent out %i particles, got %i back.", N, count_in);
if (offset + count_in > s->size_parts) {
s->size_parts = (offset + count_in) * 1.05;
struct part *parts_new;
struct xpart *xparts_new;
struct part *parts_new = NULL;
struct xpart *xparts_new = NULL;
if (posix_memalign((void **)&parts_new, part_align,
sizeof(struct part) * s->size_parts) != 0 ||
posix_memalign((void **)&xparts_new, part_align,
......@@ -2098,8 +2098,8 @@ void engine_split(struct engine *e, int *grid) {
message("Re-allocating parts array from %i to %i.", s->size_parts,
(int)(s->nr_parts * 1.2));
s->size_parts = s->nr_parts * 1.2;
struct part *parts_new;
struct xpart *xparts_new;
struct part *parts_new = NULL;
struct xpart *xparts_new = NULL;
if (posix_memalign((void **)&parts_new, part_align,
sizeof(struct part) * s->size_parts) != 0 ||
posix_memalign((void **)&xparts_new, part_align,
......
......@@ -19,6 +19,8 @@
#ifndef SWIFT_PART_H
#define SWIFT_PART_H
#include "const.h"
/* Some constants. */
#define part_maxwait 3
#define part_maxunlock 39
......
......@@ -19,7 +19,6 @@
*
******************************************************************************/
#include <math.h>
#include <stdlib.h>
#include <stddef.h>
......@@ -31,7 +30,6 @@
#include "tools.h"
#include "swift.h"
/**
* Factorize a given integer, attempts to keep larger pair of factors.
*/
......@@ -49,8 +47,6 @@ void factor(int value, int *f1, int *f2) {
}
}
/**
* @brief Compute the average number of pairs per particle using
* a brute-force O(N^2) computation.
......@@ -180,6 +176,66 @@ void pairs_single_density(double *dim, long long int pid,
fflush(stdout);
}
void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) {
float r2, hi, hj, hig2, hjg2, dx[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 = 0; j < cj->count; ++j) {
pj = &cj->parts[j];
/* Pairwise distance */
r2 = 0.0f;
for (int k = 0; k < 3; k++) {
dx[k] = ci->parts[i].x[k] - cj->parts[j].x[k];
r2 += dx[k] * dx[k];
}
/* Hit or miss? */
if (r2 < hig2) {
/* Interact */
runner_iact_nonsym_density(r2, dx, hi, pj->h, pi, pj);
}
}
}
/* Reverse double-for loop and checks every interaction */
for (int j = 0; j < cj->count; ++j) {
pj = &cj->parts[j];
hj = pj->h;
hjg2 = hj * hj * kernel_gamma2;
for (int i = 0; i < ci->count; ++i) {
pi = &ci->parts[i];
/* Pairwise distance */
r2 = 0.0f;
for (int k = 0; k < 3; k++) {
dx[k] = cj->parts[j].x[k] - ci->parts[i].x[k];
r2 += dx[k] * dx[k];
}
/* Hit or miss? */
if (r2 < hjg2) {
/* Interact */
runner_iact_nonsym_density(r2, dx, hj, pi->h, pj, pi);
}
}
}
}
void pairs_single_grav(double *dim, long long int pid,
struct gpart *__restrict__ parts, int N, int periodic) {
......@@ -234,7 +290,6 @@ void pairs_single_grav(double *dim, long long int pid,
pi.part->id, a[0], a[1], a[2], aabs[0], aabs[1], aabs[2]);
}
/**
* @brief Test the density function by dumping it for two random parts.
*
......
......@@ -19,14 +19,22 @@
*
******************************************************************************/
#ifndef SWIFT_TOOL_H
#define SWIFT_TOOL_H
#include "runner.h"
#include "cell.h"
void factor(int value, int *f1, int *f2);
void density_dump(int N);
void pairs_single_grav(double *dim, long long int pid,
struct gpart *__restrict__ parts, int N, int periodic);
void pairs_single_density(double *dim, long long int pid,
struct part *__restrict__ parts, int N,
int periodic);
struct part *__restrict__ parts, int N, int periodic);
void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj);
void pairs_n2(double *dim, struct part *__restrict__ parts, int N,
int periodic);
#endif /* SWIFT_TOOL_H */
......@@ -15,22 +15,23 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
# Add the source directory and debug to CFLAGS
AM_CFLAGS = -I../src -DCPU_TPS=2.67e9 $(HDF5_CPPFLAGS)
AM_CFLAGS = -I../src -DCPU_TPS=2.67e9 $(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
TESTS = testGreetings testReading.sh testSingle
# List of test programs to compile
check_PROGRAMS = testGreetings testReading testSingle
check_PROGRAMS = testGreetings testReading testSingle testVectorize
# Sources for the individual programs
testGreetings_SOURCES = testGreetings.c
testReading_SOURCES = testReading.c
# Sources for test_single
testSingle_SOURCES = testSingle.c
testSingle_CFLAGS = $(MYFLAGS) $(AM_CFLAGS)
testSingle_LDADD = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
testVectorize_SOURCES = testVectorize.c
#include <fenv.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <unistd.h>
#include "swift.h"
/* 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) {
size_t count = n * n * n;
struct cell *cell = malloc(sizeof *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) {
error("couldn't allocate particles, no. of particles: %d", (int)count);
}
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->dt = 0.01;
part->h = h;
part->id = ++(*partId);
part->mass = 1.0f;
++part;
}
}
}
cell->split = 0;
cell->h_max = h;
cell->count = count;
cell->dx_max = 1.;
cell->h[0] = n;
cell->h[1] = n;
cell->h[2] = n;
cell->sort = malloc(13 * count * sizeof *cell->sort);
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;
c->parts[pid].density.wcount = 0.f;
c->parts[pid].density.wcount_dh = 0.f;
c->parts[pid].density.div_v = 0.f;
c->parts[pid].density.curl_v[0] = 0.f;
c->parts[pid].density.curl_v[1] = 0.f;
c->parts[pid].density.curl_v[2] = 0.f;
}
}
/**
* @brief Dump all the particles to a file
*/
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");
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,
ci->parts[pid].rho, ci->parts[pid].rho_dh,
ci->parts[pid].density.wcount, ci->parts[pid].density.wcount_dh,
ci->parts[pid].density.div_v, ci->parts[pid].density.curl_v[0],
ci->parts[pid].density.curl_v[1], ci->parts[pid].density.curl_v[2]);
}
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,
cj->parts[pjd].rho, cj->parts[pjd].rho_dh,
cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh,
cj->parts[pjd].density.div_v, cj->parts[pjd].density.curl_v[0],
cj->parts[pjd].density.curl_v[1], cj->parts[pjd].density.curl_v[2]);
}
fclose(file);
}
/* Just a forward declaration... */
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
struct cell *ci, *cj;
struct space space;
struct engine engine;
struct runner runner;
char c;
static unsigned long long partId = 0;
ticks tic, toc, time;
while ((c = getopt(argc, argv, "h:p:r:t:")) != -1) {
switch (c) {
case 'h':
sscanf(optarg, "%lf", &h);
break;
case 'p':
sscanf(optarg, "%zu", &particles);
break;
case 'r':
sscanf(optarg, "%zu", &runs);
break;
case 't':
sscanf(optarg, "%zu", &type);
break;
}
}
if (h < 0 || particles == 0 || runs == 0 || type > 2) {
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-t TYPE=0 - cells share face (0), edge (1) or corner (2)"
"\n-h DISTANCE=1.1255 - smoothing length\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;
}
engine.s = &space;
engine.dt_step = 0.1;
runner.e = &engine;
time = 0;
for (size_t i = 0; i < runs; ++i) {
/* Zero the fields */
zero_particle_fields(ci);
zero_particle_fields(cj);
tic = getticks();
/* Run the test */
runner_dopair1_density(&runner, ci, cj);
toc = getticks();
time += toc - tic;
/* Dump if necessary */
if (i % 50 == 0) dump_particle_fields("swift_dopair.dat", ci, cj);
}
/* Output timing */
message("SWIFT calculation took %lli ticks." , time / runs);
/* Now perform a brute-force version for accuracy tests */
/* Zero the fields */
zero_particle_fields(ci);
zero_particle_fields(cj);
tic = getticks();
/* Run the test */
pairs_all_density(&runner, ci, cj);
toc = getticks();
/* Dump */
dump_particle_fields("brute_force.dat", ci, cj);
/* Output timing */
message("Brute force calculation took %lli ticks." , toc - tic);
/* Clean things to make the sanitizer happy ... */
clean_up(ci);
clean_up(cj);
return 0;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment