diff --git a/tests/testActivePair.c b/tests/testActivePair.c index c647d3ad5f27e07640fd2058d375187e6744aa6c..038d40a44afe49a12034cd9c408eba30cd026055 100644 --- a/tests/testActivePair.c +++ b/tests/testActivePair.c @@ -33,6 +33,8 @@ /* Typdef function pointer for interaction function. */ typedef void (*interaction_func)(struct runner *, struct cell *, struct cell *); +typedef void (*init_func)(struct cell *, const struct cosmology *); +typedef void (*finalise_func)(struct cell *, const struct cosmology *); /** * @brief Constructs a cell and all of its particle in a valid state prior to @@ -91,6 +93,7 @@ struct cell *make_cell(size_t n, double *offset, double size, double h, h_max = fmaxf(h_max, part->h); part->id = ++(*partId); +/* Set the mass */ #if defined(GIZMO_SPH) || defined(SHADOWFAX_SPH) part->conserved.mass = density * volume / count; @@ -98,27 +101,28 @@ struct cell *make_cell(size_t n, double *offset, double size, double h, double anchor[3] = {0., 0., 0.}; double side[3] = {1., 1., 1.}; voronoi_cell_init(&part->cell, part->x, anchor, side); -#endif +#endif /* SHADOWFAX_SPH */ #else part->mass = density * volume / count; -#endif +#endif /* GIZMO_SPH */ -#if defined(HOPKINS_PE_SPH) +/* Set the thermodynamic variable */ +#if defined(GADGET2_SPH) + part->entropy = 1.f; +#elif defined(MINIMAL_SPH) + part->u = 1.f; +#elif defined(HOPKINS_PE_SPH) part->entropy = 1.f; part->entropy_one_over_gamma = 1.f; #endif + + /* Set the time-bin */ if (random_uniform(0, 1.f) < fraction_active) part->time_bin = 1; else part->time_bin = num_time_bins + 1; - part->rho = 1.f; - part->force.f = 1.f; - part->force.P_over_rho2 = 1.f; - part->force.balsara = 1.f; - part->force.soundspeed = 1.f; - #ifdef SWIFT_DEBUG_CHECKS part->ti_drift = 8; part->ti_kick = 8; @@ -165,21 +169,57 @@ void clean_up(struct cell *ci) { /** * @brief Initializes all particles field to be ready for a density calculation */ -void zero_particle_fields(struct cell *c) { +void zero_particle_fields_density(struct cell *c, + const struct cosmology *cosmo) { for (int pid = 0; pid < c->count; pid++) { hydro_init_part(&c->parts[pid], NULL); - c->parts[pid].rho = 1.f; - c->parts[pid].force.f = 1.f; - c->parts[pid].force.P_over_rho2 = 1.f; - c->parts[pid].force.balsara = 1.f; - c->parts[pid].force.soundspeed = 1.f; } } /** - * @brief Ends the loop by adding the appropriate coefficients + * @brief Initializes all particles field to be ready for a force calculation */ -void end_calculation(struct cell *c, const struct cosmology *cosmo) { +void zero_particle_fields_force(struct cell *c, const struct cosmology *cosmo) { + for (int pid = 0; pid < c->count; pid++) { + struct part *p = &c->parts[pid]; + struct xpart *xp = &c->xparts[pid]; + +/* Mimic the result of a density calculation */ +#ifdef GADGET2_SPH + p->rho = 1.f; + p->density.rho_dh = 0.f; + p->density.wcount = 48.f / (kernel_norm * pow_dimension(p->h)); + p->density.wcount_dh = 0.f; + p->density.rot_v[0] = 0.f; + p->density.rot_v[1] = 0.f; + p->density.rot_v[2] = 0.f; + p->density.div_v = 0.f; +#endif +#ifdef MINIMAL_SPH + p->rho = 1.f; + p->density.rho_dh = 0.f; + p->density.wcount = 48.f / (kernel_norm * pow_dimension(p->h)); + p->density.wcount_dh = 0.f; +#endif +#ifdef HOPKINS_PE_SPH + p->rho = 1.f; + p->rho_bar = 1.f; + p->density.rho_dh = 0.f; + p->density.pressure_dh = 0.f; + p->density.wcount = 48.f / (kernel_norm * pow_dimension(p->h)); + p->density.wcount_dh = 0.f; +#endif + + /* And prepare for a round of force tasks. */ + hydro_prepare_force(p, xp, cosmo); + hydro_reset_acceleration(p); + } +} + +/** + * @brief Ends the density loop by adding the appropriate coefficients + */ +void end_calculation_density(struct cell *c, const struct cosmology *cosmo) { for (int pid = 0; pid < c->count; pid++) { hydro_end_density(&c->parts[pid], cosmo); @@ -189,6 +229,15 @@ void end_calculation(struct cell *c, const struct cosmology *cosmo) { } } +/** + * @brief Ends the force loop by adding the appropriate coefficients + */ +void end_calculation_force(struct cell *c, const struct cosmology *cosmo) { + for (int pid = 0; pid < c->count; pid++) { + hydro_end_force(&c->parts[pid], cosmo); + } +} + /** * @brief Dump all the particles to a file */ @@ -233,21 +282,22 @@ void test_pair_interactions(struct runner *runner, struct cell **ci, struct cell **cj, char *swiftOutputFileName, char *bruteForceOutputFileName, interaction_func serial_interaction, - interaction_func vec_interaction) { + interaction_func vec_interaction, init_func init, + finalise_func finalise) { 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); + init(*ci, runner->e->cosmology); + init(*cj, runner->e->cosmology); /* Run the test */ vec_interaction(runner, *ci, *cj); /* Let's get physical ! */ - end_calculation(*ci, runner->e->cosmology); - end_calculation(*cj, runner->e->cosmology); + finalise(*ci, runner->e->cosmology); + finalise(*cj, runner->e->cosmology); /* Dump if necessary */ dump_particle_fields(swiftOutputFileName, *ci, *cj); @@ -255,15 +305,15 @@ void test_pair_interactions(struct runner *runner, struct cell **ci, /* Now perform a brute-force version for accuracy tests */ /* Zero the fields */ - zero_particle_fields(*ci); - zero_particle_fields(*cj); + init(*ci, runner->e->cosmology); + init(*cj, runner->e->cosmology); /* Run the brute-force test */ serial_interaction(runner, *ci, *cj); /* Let's get physical ! */ - end_calculation(*ci, runner->e->cosmology); - end_calculation(*cj, runner->e->cosmology); + finalise(*ci, runner->e->cosmology); + finalise(*cj, runner->e->cosmology); dump_particle_fields(bruteForceOutputFileName, *ci, *cj); } @@ -275,7 +325,8 @@ 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, - interaction_func serial_interaction, interaction_func vec_interaction) { + interaction_func serial_interaction, interaction_func vec_interaction, + init_func init, finalise_func finalise) { double offset1[3] = {0, 0, 0}; struct cell *ci, *cj; @@ -286,7 +337,7 @@ void test_all_pair_interactions( test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName, serial_interaction, - vec_interaction); + vec_interaction, init, finalise); clean_up(ci); clean_up(cj); @@ -299,7 +350,7 @@ void test_all_pair_interactions( test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName, serial_interaction, - vec_interaction); + vec_interaction, init, finalise); clean_up(ci); clean_up(cj); @@ -312,7 +363,7 @@ void test_all_pair_interactions( test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName, serial_interaction, - vec_interaction); + vec_interaction, init, finalise); clean_up(ci); clean_up(cj); @@ -325,7 +376,7 @@ void test_all_pair_interactions( test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName, serial_interaction, - vec_interaction); + vec_interaction, init, finalise); clean_up(ci); clean_up(cj); @@ -338,7 +389,7 @@ void test_all_pair_interactions( test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName, serial_interaction, - vec_interaction); + vec_interaction, init, finalise); clean_up(ci); clean_up(cj); @@ -351,7 +402,7 @@ void test_all_pair_interactions( test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName, serial_interaction, - vec_interaction); + vec_interaction, init, finalise); clean_up(ci); clean_up(cj); @@ -364,7 +415,7 @@ void test_all_pair_interactions( test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName, serial_interaction, - vec_interaction); + vec_interaction, init, finalise); clean_up(ci); clean_up(cj); @@ -375,7 +426,7 @@ void test_all_pair_interactions( test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName, serial_interaction, - vec_interaction); + vec_interaction, init, finalise); clean_up(ci); clean_up(cj); @@ -386,7 +437,7 @@ void test_all_pair_interactions( test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName, serial_interaction, - vec_interaction); + vec_interaction, init, finalise); clean_up(ci); clean_up(cj); @@ -399,7 +450,7 @@ void test_all_pair_interactions( test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName, serial_interaction, - vec_interaction); + vec_interaction, init, finalise); clean_up(ci); clean_up(cj); @@ -412,7 +463,7 @@ void test_all_pair_interactions( test_pair_interactions(runner, &ci, &cj, swiftOutputFileName, bruteForceOutputFileName, serial_interaction, - vec_interaction); + vec_interaction, init, finalise); /* Clean things to make the sanitizer happy ... */ clean_up(ci); @@ -539,12 +590,14 @@ int main(int argc, char *argv[]) { /* Define which interactions to call */ interaction_func serial_inter_func = &pairs_all_density; interaction_func vec_inter_func = &runner_dopair1_branch_density; + init_func init = &zero_particle_fields_density; + finalise_func finalise = &end_calculation_density; /* Test a pair of cells face-on. */ test_all_pair_interactions(runner, offset, particles, size, h, rho, &partId, perturbation, h_pert, swiftOutputFileName, bruteForceOutputFileName, serial_inter_func, - vec_inter_func); + vec_inter_func, init, finalise); /* Test a pair of cells edge-on. */ offset[0] = 1.; @@ -553,7 +606,7 @@ int main(int argc, char *argv[]) { test_all_pair_interactions(runner, offset, particles, size, h, rho, &partId, perturbation, h_pert, swiftOutputFileName, bruteForceOutputFileName, serial_inter_func, - vec_inter_func); + vec_inter_func, init, finalise); /* Test a pair of cells corner-on. */ offset[0] = 1.; @@ -562,11 +615,13 @@ int main(int argc, char *argv[]) { test_all_pair_interactions(runner, offset, particles, size, h, rho, &partId, perturbation, h_pert, swiftOutputFileName, bruteForceOutputFileName, serial_inter_func, - vec_inter_func); + vec_inter_func, init, finalise); /* Re-assign function pointers. */ serial_inter_func = &pairs_all_force; vec_inter_func = &runner_dopair2_branch_force; + init = &zero_particle_fields_force; + finalise = &end_calculation_force; /* Create new output file names. */ sprintf(swiftOutputFileName, "swift_dopair2_force_%s.dat", @@ -585,7 +640,7 @@ int main(int argc, char *argv[]) { test_all_pair_interactions(runner, offset, particles, size, h, rho, &partId, perturbation, h_pert, swiftOutputFileName, bruteForceOutputFileName, serial_inter_func, - vec_inter_func); + vec_inter_func, init, finalise); /* Test a pair of cells edge-on. */ offset[0] = 1.; @@ -594,7 +649,7 @@ int main(int argc, char *argv[]) { test_all_pair_interactions(runner, offset, particles, size, h, rho, &partId, perturbation, h_pert, swiftOutputFileName, bruteForceOutputFileName, serial_inter_func, - vec_inter_func); + vec_inter_func, init, finalise); /* Test a pair of cells corner-on. */ offset[0] = 1.; @@ -603,6 +658,6 @@ int main(int argc, char *argv[]) { test_all_pair_interactions(runner, offset, particles, size, h, rho, &partId, perturbation, h_pert, swiftOutputFileName, bruteForceOutputFileName, serial_inter_func, - vec_inter_func); + vec_inter_func, init, finalise); return 0; }