Skip to content
Snippets Groups Projects
Commit e7fe7d56 authored by James Willis's avatar James Willis
Browse files

Only tests cells with particles that are all active.

parent b68d8ffc
No related branches found
No related tags found
1 merge request!396Avx512 fixes
......@@ -28,14 +28,27 @@
/* Local headers. */
#include "swift.h"
/* n is both particles per axis and box size:
* particles are generated on a mesh with unit spacing
/**
* @brief Constructs a cell and all of its particle in a valid state prior to
* a DOPAIR or DOSELF calcuation.
*
* @param n The cube root of the number of particles.
* @param offset The position of the cell offset from (0,0,0).
* @param size The cell size.
* @param h The smoothing length of the particles in units of the inter-particle
* separation.
* @param density The density of the fluid.
* @param partId The running counter of IDs.
* @param pert The perturbation to apply to the particles in the cell in units
* of the inter-particle separation.
* @param h_pert The perturbation to apply to the smoothing length.
*/
struct cell *make_cell(size_t n, double *offset, double size, double h,
double density, unsigned long long *partId,
double pert) {
double density, long long *partId, double pert,
double h_pert) {
const size_t count = n * n * n;
const double volume = size * size * size;
float h_max = 0.f;
struct cell *cell = malloc(sizeof(struct cell));
bzero(cell, sizeof(struct cell));
......@@ -51,26 +64,44 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
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;
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;
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;
offset[2] +
size * (z + 0.5 + random_uniform(-0.5, 0.5) * pert) / (float)n;
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;
if (h_pert)
part->h = size * h * random_uniform(1.f, h_pert) / (float)n;
else
part->h = size * h / (float)n;
h_max = fmaxf(h_max, part->h);
part->id = ++(*partId);
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
part->conserved.mass = density * volume / count;
#ifdef SHADOWFAX_SPH
double anchor[3] = {0., 0., 0.};
double side[3] = {1., 1., 1.};
voronoi_cell_init(&part->cell, part->x, anchor, side);
#endif
#else
part->mass = density * volume / count;
#endif
#if defined(HOPKINS_PE_SPH)
part->entropy = 1.f;
part->entropy_one_over_gamma = 1.f;
#endif
part->time_bin = 1;
//part->time_bin = num_time_bins + 1;
#ifdef SWIFT_DEBUG_CHECKS
part->ti_drift = 8;
......@@ -84,13 +115,13 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
/* Cell properties */
cell->split = 0;
cell->h_max = h;
cell->h_max = h_max;
cell->count = count;
cell->dx_max_part = 0.;
cell->dx_max_sort = 0.;
cell->width[0] = n;
cell->width[1] = n;
cell->width[2] = n;
cell->width[0] = size;
cell->width[1] = size;
cell->width[2] = size;
cell->loc[0] = offset[0];
cell->loc[1] = offset[1];
cell->loc[2] = offset[2];
......@@ -98,6 +129,7 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
cell->ti_old_part = 8;
cell->ti_end_min = 8;
cell->ti_end_max = 8;
cell->ti_sort = 8;
shuffle_particles(cell->parts, cell->count);
......@@ -122,6 +154,19 @@ void zero_particle_fields(struct cell *c) {
}
}
/**
* @brief Ends the loop by adding the appropriate coefficients
*/
void end_calculation(struct cell *c) {
for (int pid = 0; pid < c->count; pid++) {
hydro_end_density(&c->parts[pid]);
/* Recover the common "Neighbour number" definition */
c->parts[pid].density.wcount *= pow_dimension(c->parts[pid].h);
c->parts[pid].density.wcount *= kernel_norm;
}
}
/**
* @brief Dump all the particles to a file
*/
......@@ -130,57 +175,24 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
/* Write header */
fprintf(file,
"# %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");
"# %4s %13s\n", "ID", "wcount");
fprintf(file, "# ci --------------------------------------------\n");
for (int pid = 0; pid < ci->count; pid++) {
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], hydro_get_density(&ci->parts[pid]),
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
0.f,
#else
ci->parts[pid].density.rho_dh,
#endif
ci->parts[pid].density.wcount, ci->parts[pid].density.wcount_dh,
#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH)
ci->parts[pid].density.div_v, ci->parts[pid].density.rot_v[0],
ci->parts[pid].density.rot_v[1], ci->parts[pid].density.rot_v[2]
#else
0., 0., 0., 0.
#endif
);
"%6llu %13e\n",
ci->parts[pid].id,
ci->parts[pid].density.wcount);
}
fprintf(file, "# cj --------------------------------------------\n");
for (int 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].v[0], cj->parts[pjd].v[1],
cj->parts[pjd].v[2], hydro_get_density(&cj->parts[pjd]),
#if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH)
0.f,
#else
cj->parts[pjd].density.rho_dh,
#endif
cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh,
#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH)
cj->parts[pjd].density.div_v, cj->parts[pjd].density.rot_v[0],
cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2]
#else
0., 0., 0., 0.
#endif
);
"%6llu %13e\n",
cj->parts[pjd].id,
cj->parts[pjd].density.wcount);
}
fclose(file);
......@@ -192,19 +204,64 @@ void runner_doself1_density_vec(struct runner *r, struct cell *ci);
void runner_dopair1_branch_density(struct runner *r, struct cell *ci,
struct cell *cj);
void test_pair_interactions(struct runner runner,
struct cell **ci, struct cell **cj, char *swiftOutputFileName,
char *bruteForceOutputFileName) {
runner_do_sort(&runner, *ci, 0x1FFF, 0, 0);
runner_do_sort(&runner, *cj, 0x1FFF, 0, 0);
/* Zero the fields */
zero_particle_fields(*ci);
zero_particle_fields(*cj);
#ifdef WITH_VECTORIZATION
runner.ci_cache.count = 0;
cache_init(&runner.ci_cache, 512);
runner.cj_cache.count = 0;
cache_init(&runner.cj_cache, 512);
#endif
/* Run the test */
runner_dopair1_branch_density(&runner, *ci, *cj);
/* Let's get physical ! */
end_calculation(*ci);
end_calculation(*cj);
/* Dump if necessary */
dump_particle_fields(swiftOutputFileName, *ci, *cj);
/* Now perform a brute-force version for accuracy tests */
/* Zero the fields */
zero_particle_fields(*ci);
zero_particle_fields(*cj);
/* Run the brute-force test */
pairs_all_density(&runner, *ci, *cj);
/* Let's get physical ! */
end_calculation(*ci);
end_calculation(*cj);
dump_particle_fields(bruteForceOutputFileName, *ci, *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, size = 1., rho = 1.;
double perturbation = 0.1;
size_t particles = 0, runs = 0, type = 0;
double offset[3] = {0, 0, 0}, h = 1.23485, size = 1., rho = 1.;
double perturbation = 0.1, h_pert = 1.1;
struct cell *ci, *cj;
struct space space;
struct engine engine;
struct runner runner;
char c;
static unsigned long long partId = 0;
static long long partId = 0;
char outputFileNameExtension[200] = "";
char outputFileName[200] = "";
ticks tic, toc, time;
char swiftOutputFileName[200] = "";
char bruteForceOutputFileName[200] = "";
/* Initialize CPU frequency, this also starts time. */
unsigned long long cpufreq = 0;
......@@ -215,12 +272,15 @@ int main(int argc, char *argv[]) {
srand(0);
while ((c = getopt(argc, argv, "h:p:r:t:d:f:")) != -1) {
while ((c = getopt(argc, argv, "h:p:n:r:t:d:f:")) != -1) {
switch (c) {
case 'h':
sscanf(optarg, "%lf", &h);
break;
case 'p':
sscanf(optarg, "%lf", &h_pert);
break;
case 'n':
sscanf(optarg, "%zu", &particles);
break;
case 'r':
......@@ -243,12 +303,13 @@ int main(int argc, char *argv[]) {
if (h < 0 || particles == 0 || runs == 0 || type > 2) {
printf(
"\nUsage: %s -p PARTICLES_PER_AXIS -r NUMBER_OF_RUNS [OPTIONS...]\n"
"\nUsage: %s -n 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-h DISTANCE=1.2348 - smoothing length"
"\n-p - Random fractional change in h, h=h*random(1,p)"
"\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]);
......@@ -263,62 +324,17 @@ int main(int argc, char *argv[]) {
engine.max_active_bin = num_time_bins;
runner.e = &engine;
volume = particles * particles * particles;
message("particles: %zu B\npositions: 0 B", 2 * volume * sizeof(struct part));
ci = make_cell(particles, offset, size, h, rho, &partId, perturbation);
/* Create output file names. */
sprintf(swiftOutputFileName, "swift_dopair_%s.dat",
outputFileNameExtension);
sprintf(bruteForceOutputFileName, "brute_force_%s.dat",
outputFileNameExtension);
ci = make_cell(particles, offset, size, h, rho, &partId, perturbation, h_pert);
for (size_t i = 0; i < type + 1; ++i) offset[i] = 1.;
cj = make_cell(particles, offset, size, h, rho, &partId, perturbation);
runner_do_sort(&runner, ci, 0x1FFF, 0, 0);
runner_do_sort(&runner, cj, 0x1FFF, 0, 0);
time = 0;
/* Zero the fields */
zero_particle_fields(ci);
zero_particle_fields(cj);
#ifdef WITH_VECTORIZATION
runner.ci_cache.count = 0;
cache_init(&runner.ci_cache, 512);
runner.cj_cache.count = 0;
cache_init(&runner.cj_cache, 512);
#endif
tic = getticks();
/* Run the test */
runner_dopair1_branch_density(&runner, ci, cj);
toc = getticks();
time += toc - tic;
/* Dump if necessary */
sprintf(outputFileName, "swift_dopair_%s.dat", outputFileNameExtension);
dump_particle_fields(outputFileName, 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 brute-force test */
pairs_all_density(&runner, ci, cj);
toc = getticks();
/* Dump */
sprintf(outputFileName, "brute_force_%s.dat", outputFileNameExtension);
dump_particle_fields(outputFileName, ci, cj);
cj = make_cell(particles, offset, size, h, rho, &partId, perturbation, h_pert);
/* Output timing */
message("Brute force calculation took %lli ticks.", toc - tic);
test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName);
/* Clean things to make the sanitizer happy ... */
clean_up(ci);
......
......@@ -4,7 +4,7 @@ echo ""
rm -f brute_force_pair_active.dat swift_dopair_active.dat
./testActivePair -p 6 -r 1 -d 0 -f active
./testActivePair -n 6 -r 1 -d 0 -f active
python @srcdir@/difffloat.py brute_force_active.dat swift_dopair_active.dat @srcdir@/tolerance_pair_active.dat
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment