diff --git a/src/runner_doiact_vec.c b/src/runner_doiact_vec.c index 26bd2840efb418db1f166569d2734823abcafd81..1ab6cc3a50423be6d9e3f72c26be77ba7432a70e 100644 --- a/src/runner_doiact_vec.c +++ b/src/runner_doiact_vec.c @@ -265,10 +265,11 @@ __attribute__((always_inline)) INLINE static void populate_max_d_no_cache( int first_pi = 0, last_pj = cj->count - 1; int temp; - /* Find the leftmost active particle in cell i that interacts with any particle in cell j. */ + /* Find the leftmost active particle in cell i that interacts with any + * particle in cell j. */ first_pi = ci->count; int active_id = first_pi; - while(first_pi > 0 && sort_i[first_pi - 1].d + dx_max + hi_max > dj_min) { + while (first_pi > 0 && sort_i[first_pi - 1].d + dx_max + hi_max > dj_min) { first_pi--; /* Store the index of the particle if it is active. */ if (part_is_active(&parts_i[sort_i[first_pi].i], e)) active_id = first_pi; @@ -278,7 +279,7 @@ __attribute__((always_inline)) INLINE static void populate_max_d_no_cache( first_pi = active_id; /* Find the maximum index into cell j for each particle in range in cell i. */ - if(first_pi < ci->count) { + if (first_pi < ci->count) { /* Start from the first particle in cell j. */ temp = 0; @@ -286,31 +287,35 @@ __attribute__((always_inline)) INLINE static void populate_max_d_no_cache( const struct part *pi = &parts_i[sort_i[first_pi].i]; /* Loop through particles in cell j until they are not in range of pi. */ - while(temp <= cj->count && (sort_i[first_pi].d + (pi->h * kernel_gamma + dx_max - rshift) > sort_j[temp].d)) + while (temp <= cj->count && + (sort_i[first_pi].d + (pi->h * kernel_gamma + dx_max - rshift) > + sort_j[temp].d)) temp++; max_index_i[first_pi] = temp; /* Populate max_index_i for remaining particles that are within range. */ - for(int i = first_pi + 1; i<ci->count; i++) { + for (int i = first_pi + 1; i < ci->count; i++) { temp = max_index_i[i - 1]; - while(temp <= cj->count && (sort_i[i].d + (pi->h * kernel_gamma + dx_max - rshift) > sort_j[temp].d)) + while (temp <= cj->count && + (sort_i[i].d + (pi->h * kernel_gamma + dx_max - rshift) > + sort_j[temp].d)) temp++; max_index_i[i] = temp; - } - } - else { + } else { /* Make sure that max index is set to first particle in cj.*/ max_index_i[ci->count - 1] = 0; } - /* Find the rightmost active particle in cell j that interacts with any particle in cell i. */ + /* Find the rightmost active particle in cell j that interacts with any + * particle in cell i. */ last_pj = -1; active_id = last_pj; - while(last_pj < cj->count && sort_j[last_pj + 1].d - hj_max - dx_max < di_max) { + while (last_pj < cj->count && + sort_j[last_pj + 1].d - hj_max - dx_max < di_max) { last_pj++; /* Store the index of the particle if it is active. */ if (part_is_active(&parts_j[sort_j[last_pj].i], e)) active_id = last_pj; @@ -320,7 +325,7 @@ __attribute__((always_inline)) INLINE static void populate_max_d_no_cache( last_pj = active_id; /* Find the maximum index into cell i for each particle in range in cell j. */ - if(last_pj > 0 ) { + if (last_pj > 0) { /* Start from the last particle in cell i. */ temp = ci->count - 1; @@ -328,25 +333,27 @@ __attribute__((always_inline)) INLINE static void populate_max_d_no_cache( const struct part *pj = &parts_j[sort_j[last_pj].i]; /* Loop through particles in cell i until they are not in range of pj. */ - while(temp > 0 && sort_j[last_pj].d - dx_max - (pj->h * kernel_gamma) < sort_i[temp].d - rshift) + while (temp > 0 && + sort_j[last_pj].d - dx_max - (pj->h * kernel_gamma) < + sort_i[temp].d - rshift) temp--; max_index_j[last_pj] = temp; /* Populate max_index_j for remaining particles that are within range. */ - for(int i = last_pj - 1; i>=0; i--) { + for (int i = last_pj - 1; i >= 0; i--) { temp = max_index_j[i + 1]; - while(temp > 0 && sort_j[i].d - dx_max - (pj->h * kernel_gamma) < sort_i[temp].d - rshift) + while (temp > 0 && + sort_j[i].d - dx_max - (pj->h * kernel_gamma) < + sort_i[temp].d - rshift) temp--; max_index_j[i] = temp; - } - } - else { + } else { /* Make sure that max index is set to last particle in ci.*/ - max_index_j[0] = ci->count - 1; + max_index_j[0] = ci->count - 1; } *init_pi = first_pi; @@ -695,8 +702,8 @@ void runner_dopair1_density_vec(struct runner *r, struct cell *ci, /* Also find the first pi that interacts with any particle in cj and the last * pj that interacts with any particle in ci. */ populate_max_d_no_cache(ci, cj, sort_i, sort_j, dx_max, rshift, hi_max, - hj_max, di_max, dj_min, max_index_i, max_index_j, &first_pi, - &last_pj, e); + hj_max, di_max, dj_min, max_index_i, max_index_j, + &first_pi, &last_pj, e); /* Limits of the outer loops. */ int first_pi_loop = first_pi; @@ -720,7 +727,8 @@ void runner_dopair1_density_vec(struct runner *r, struct cell *ci, if (cell_is_active(ci, e)) { /* Loop over the parts in ci until nothing is within range in cj. */ - //for (int pid = count_i - 1; pid >= first_pi_loop && max_index_i[pid] >= 0; pid--) { + // for (int pid = count_i - 1; pid >= first_pi_loop && max_index_i[pid] >= + // 0; pid--) { for (int pid = count_i - 1; pid >= first_pi_loop; pid--) { /* Get a hold of the ith part in ci. */ @@ -849,7 +857,8 @@ void runner_dopair1_density_vec(struct runner *r, struct cell *ci, if (cell_is_active(cj, e)) { /* Loop over the parts in cj until nothing is within range in ci. */ - //for (int pjd = 0; pjd <= last_pj_loop && max_index_j[pjd] < count_i; pjd++) { + // for (int pjd = 0; pjd <= last_pj_loop && max_index_j[pjd] < count_i; + // pjd++) { for (int pjd = 0; pjd <= last_pj_loop; pjd++) { /* Get a hold of the jth part in cj. */ diff --git a/src/tools.c b/src/tools.c index 6aa73276cbcbdd23bc6b47132751fc5214d3ed19..7d69ebc6c476312081d8a8c34c76c6592da5cab0 100644 --- a/src/tools.c +++ b/src/tools.c @@ -32,6 +32,7 @@ #include "tools.h" /* Local includes. */ +#include "active.h" #include "cell.h" #include "error.h" #include "gravity.h" @@ -39,7 +40,6 @@ #include "part.h" #include "periodic.h" #include "runner.h" -#include "active.h" /** * Factorize a given integer, attempts to keep larger pair of factors. @@ -194,7 +194,7 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) { hig2 = hi * hi * kernel_gamma2; /* Skip inactive particles. */ - if(!part_is_active(pi,e)) continue; + if (!part_is_active(pi, e)) continue; for (int j = 0; j < cj->count; ++j) { @@ -225,8 +225,8 @@ void pairs_all_density(struct runner *r, struct cell *ci, struct cell *cj) { hjg2 = hj * hj * kernel_gamma2; /* Skip inactive particles. */ - if(!part_is_active(pj,e)) continue; - + if (!part_is_active(pj, e)) continue; + for (int i = 0; i < ci->count; ++i) { pi = &ci->parts[i]; @@ -344,14 +344,14 @@ void self_all_density(struct runner *r, struct cell *ci) { } /* Hit or miss? */ - if (r2 < hig2 && part_is_active(pi,e)) { + if (r2 < hig2 && part_is_active(pi, e)) { /* Interact */ runner_iact_nonsym_density(r2, dxi, hi, hj, pi, pj); } /* Hit or miss? */ - if (r2 < hjg2 && part_is_active(pj,e)) { + if (r2 < hjg2 && part_is_active(pj, e)) { dxi[0] = -dxi[0]; dxi[1] = -dxi[1]; diff --git a/tests/testActivePair.c b/tests/testActivePair.c index 7740d24b594905ab2ff23b27c1a6447db9192dec..0f8ecdea875557d972ff37dd9bc9a51d81911050 100644 --- a/tests/testActivePair.c +++ b/tests/testActivePair.c @@ -23,8 +23,8 @@ #include <stdio.h> #include <stdlib.h> #include <string.h> -#include <unistd.h> #include <time.h> +#include <unistd.h> /* Local headers. */ #include "swift.h" @@ -65,14 +65,14 @@ 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); @@ -101,8 +101,8 @@ struct cell *make_cell(size_t n, double *offset, double size, double h, part->entropy = 1.f; part->entropy_one_over_gamma = 1.f; #endif - if (random_uniform(0,1.f) < fraction_active) - part->time_bin = 1; + if (random_uniform(0, 1.f) < fraction_active) + part->time_bin = 1; else part->time_bin = num_time_bins + 1; @@ -137,7 +137,7 @@ struct cell *make_cell(size_t n, double *offset, double size, double h, cell->sorted = 0; for (int k = 0; k < 13; k++) cell->sort[k] = NULL; - + return cell; } @@ -177,24 +177,19 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) { FILE *file = fopen(fileName, "a"); /* Write header */ - fprintf(file, - "# %4s %13s\n", "ID", "wcount"); + fprintf(file, "# %4s %13s\n", "ID", "wcount"); fprintf(file, "# ci --------------------------------------------\n"); for (int pid = 0; pid < ci->count; pid++) { - fprintf(file, - "%6llu %13e\n", - ci->parts[pid].id, + fprintf(file, "%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 %13e\n", - cj->parts[pjd].id, + fprintf(file, "%6llu %13e\n", cj->parts[pjd].id, cj->parts[pjd].density.wcount); } @@ -208,10 +203,11 @@ void runner_dopair1_branch_density(struct runner *r, struct cell *ci, struct cell *cj); /** - * @brief Computes the pair interactions of two cells using SWIFT and a brute force implementation. + * @brief Computes the pair interactions of two cells using SWIFT and a brute + * force implementation. */ -void test_pair_interactions(struct runner *runner, - struct cell **ci, struct cell **cj, char *swiftOutputFileName, +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); @@ -227,7 +223,7 @@ void test_pair_interactions(struct runner *runner, /* Let's get physical ! */ end_calculation(*ci); end_calculation(*cj); - + /* Dump if necessary */ dump_particle_fields(swiftOutputFileName, *ci, *cj); @@ -245,70 +241,90 @@ void test_pair_interactions(struct runner *runner, end_calculation(*cj); dump_particle_fields(bruteForceOutputFileName, *ci, *cj); - } /** * @brief Computes the pair interactions of two cells in various configurations. */ -void test_all_pair_interactions(struct runner *runner, - double *offset2, size_t particles, double size, double h, double rho, long long *partId, double perturbation, double h_pert, char *swiftOutputFileName, - char *bruteForceOutputFileName) { +void test_all_pair_interactions(struct runner *runner, double *offset2, + size_t particles, double size, double h, + double rho, long long *partId, + double perturbation, double h_pert, + char *swiftOutputFileName, + char *bruteForceOutputFileName) { double offset1[3] = {0, 0, 0}; struct cell *ci, *cj; /* All active particles. */ - ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert, 1.); - cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert, 1.); + ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert, + 1.); + cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert, + 1.); - test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName); + test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, + bruteForceOutputFileName); clean_up(ci); clean_up(cj); /* Half particles are active. */ - ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert, 0.5); - cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert, 0.5); + ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert, + 0.5); + cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert, + 0.5); + + test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, + bruteForceOutputFileName); - test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName); - clean_up(ci); clean_up(cj); /* All particles inactive. */ - ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert, 0.); - cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert, 0.); + ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert, + 0.); + cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert, + 0.); - test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName); + test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, + bruteForceOutputFileName); clean_up(ci); clean_up(cj); /* 10% of particles active. */ - ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert, 0.1); - cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert, 0.1); + ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert, + 0.1); + cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert, + 0.1); + + test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, + bruteForceOutputFileName); - test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName); - clean_up(ci); clean_up(cj); /* One active cell one inactive cell. */ - ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert, 1.0); - cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert, 0.); + ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert, + 1.0); + cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert, + 0.); + + test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, + bruteForceOutputFileName); - test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName); - clean_up(ci); clean_up(cj); /* One active cell one inactive cell. */ - ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert, 0.); - cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert, 1.0); + ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert, + 0.); + cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert, + 1.0); + + test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, + bruteForceOutputFileName); - test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName); - clean_up(ci); clean_up(cj); @@ -316,8 +332,9 @@ void test_all_pair_interactions(struct runner *runner, ci = make_cell(2, offset1, size, h, rho, partId, perturbation, h_pert, 1.0); cj = make_cell(2, offset2, size, h, rho, partId, perturbation, h_pert, 1.0); - test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName); - + test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, + bruteForceOutputFileName); + clean_up(ci); clean_up(cj); @@ -325,26 +342,33 @@ void test_all_pair_interactions(struct runner *runner, ci = make_cell(10, offset1, size, h, rho, partId, perturbation, h_pert, 0.5); cj = make_cell(3, offset2, size, h, rho, partId, perturbation, h_pert, 0.75); - test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName); - + test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, + bruteForceOutputFileName); + clean_up(ci); clean_up(cj); /* One cell inactive and the other only half active. */ - ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert, 0.5); - cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert, 0.); + ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert, + 0.5); + cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert, + 0.); + + test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, + bruteForceOutputFileName); - test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName); - clean_up(ci); clean_up(cj); - + /* One cell inactive and the other only half active. */ - ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert, 0.); - cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert, 0.5); + ci = make_cell(particles, offset1, size, h, rho, partId, perturbation, h_pert, + 0.); + cj = make_cell(particles, offset2, size, h, rho, partId, perturbation, h_pert, + 0.5); + + test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, + bruteForceOutputFileName); - test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName); - /* Clean things to make the sanitizer happy ... */ clean_up(ci); clean_up(cj); @@ -424,27 +448,26 @@ int main(int argc, char *argv[]) { /* Seed RNG. */ message("Seed used for RNG: %d", seed); srand(seed); - + space.periodic = 0; engine.s = &space; engine.time = 0.1f; engine.ti_current = 8; engine.max_active_bin = num_time_bins; - - if (posix_memalign((void **)&runner, part_align, - sizeof(struct runner)) != 0) { + + if (posix_memalign((void **)&runner, part_align, sizeof(struct runner)) != + 0) { error("couldn't allocate runner"); } - + runner->e = &engine; /* Create output file names. */ - sprintf(swiftOutputFileName, "swift_dopair_%s.dat", - outputFileNameExtension); + sprintf(swiftOutputFileName, "swift_dopair_%s.dat", outputFileNameExtension); sprintf(bruteForceOutputFileName, "brute_force_%s.dat", outputFileNameExtension); - + /* Delete files if they already exist. */ remove(swiftOutputFileName); remove(bruteForceOutputFileName); @@ -456,17 +479,27 @@ int main(int argc, char *argv[]) { cache_init(&runner->cj_cache, 512); #endif - double offset[3] = {1.,0.,0.}; + double offset[3] = {1., 0., 0.}; /* Test a pair of cells face-on. */ - test_all_pair_interactions(runner, offset, particles, size, h, rho, &partId, perturbation, h_pert, swiftOutputFileName, bruteForceOutputFileName); - + test_all_pair_interactions(runner, offset, particles, size, h, rho, &partId, + perturbation, h_pert, swiftOutputFileName, + bruteForceOutputFileName); + /* Test a pair of cells edge-on. */ - offset[0] = 1.; offset[1] = 1.; offset[2] = 0.; - test_all_pair_interactions(runner, offset, particles, size, h, rho, &partId, perturbation, h_pert, swiftOutputFileName, bruteForceOutputFileName); + offset[0] = 1.; + offset[1] = 1.; + offset[2] = 0.; + test_all_pair_interactions(runner, offset, particles, size, h, rho, &partId, + perturbation, h_pert, swiftOutputFileName, + bruteForceOutputFileName); /* Test a pair of cells corner-on. */ - offset[0] = 1.; offset[1] = 1.; offset[2] = 1.; - test_all_pair_interactions(runner, offset, particles, size, h, rho, &partId, perturbation, h_pert, swiftOutputFileName, bruteForceOutputFileName); + offset[0] = 1.; + offset[1] = 1.; + offset[2] = 1.; + test_all_pair_interactions(runner, offset, particles, size, h, rho, &partId, + perturbation, h_pert, swiftOutputFileName, + bruteForceOutputFileName); return 0; } diff --git a/tests/testInteractions.c b/tests/testInteractions.c index 781bdda74e8c750a899e84d7eacbd733caa3a2dd..fc7aef63329c5665da2bc899c765d8ed7ad43807 100644 --- a/tests/testInteractions.c +++ b/tests/testInteractions.c @@ -310,7 +310,7 @@ void test_interactions(struct part test_part, struct part *parts, size_t count, vjzq[i] = pj_vec[i].v[2]; } -/* Perform vector interaction. */ + /* Perform vector interaction. */ vector hi_vec, hi_inv_vec, vix_vec, viy_vec, viz_vec; vector rhoSum, rho_dhSum, wcountSum, wcount_dhSum, div_vSum, curlvxSum, curlvySum, curlvzSum; @@ -377,7 +377,7 @@ void test_interactions(struct part test_part, struct part *parts, size_t count, dump_indv_particle_fields(vec_filename, piq[0]); for (size_t i = 0; i < count; i++) dump_indv_particle_fields(vec_filename, pjq[i]); - + /* Check serial results against the vectorised results. */ if (check_results(pi_serial, pj_serial, pi_vec, pj_vec, count)) message("Differences found..."); @@ -454,6 +454,6 @@ int main(int argc, char *argv[]) { #else -int main() {return 1;} +int main() { return 1; } #endif diff --git a/tests/testPeriodicBC.c b/tests/testPeriodicBC.c index 0f350f607b54e790c8db9209b8fb18335f5009d1..6fa2dc607b996b9e8508338a9806633c5a4d1a89 100644 --- a/tests/testPeriodicBC.c +++ b/tests/testPeriodicBC.c @@ -517,7 +517,8 @@ int main(int argc, char *argv[]) { const int half_dim = (dim - 1) / 2; - /* Test the periodic boundary conditions for each of the 8 corners. Interact each corner with all of its 26 neighbours.*/ + /* Test the periodic boundary conditions for each of the 8 corners. Interact + * each corner with all of its 26 neighbours.*/ test_boundary_conditions(cells, runner, 0, 0, 0, dim, swiftOutputFileName, bruteForceOutputFileName); test_boundary_conditions(cells, runner, dim - 1, 0, 0, dim,