diff --git a/tests/testActivePair.c b/tests/testActivePair.c index 990370d44d124a8065b95870b17e939b3618913c..bd4d027f3771d63e5ae450be4aa0aa56604641f0 100644 --- a/tests/testActivePair.c +++ b/tests/testActivePair.c @@ -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); diff --git a/tests/testActivePair.sh.in b/tests/testActivePair.sh.in index eb3cf5da3e9f456f2bb7b12f2e884e778fa86bca..ff8d027a469bd9bc78286b843cf2dffd3ef27ad3 100755 --- a/tests/testActivePair.sh.in +++ b/tests/testActivePair.sh.in @@ -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