diff --git a/tests/fft_params.yml b/tests/fft_params.yml index 05d6d8f0b0578d11645fc1d78c1a6322160ae87a..6938e3658148874e58f50ab768a5e1fbc41d9573 100644 --- a/tests/fft_params.yml +++ b/tests/fft_params.yml @@ -3,8 +3,9 @@ Scheduler: # Parameters for the self-gravity scheme Gravity: - eta: 0.025 # Constant dimensionless multiplier for time integration. - theta: 0.7 # Opening angle (Multipole acceptance criterion) - epsilon: 0.00001 # Softening length (in internal units). - a_smooth: 0. - r_cut: 0. + eta: 0.025 # Constant dimensionless multiplier for time integration. + theta: 0.7 # Opening angle (Multipole acceptance criterion) + comoving_softening: 0.00001 # Comoving softening length (in internal units). + max_physical_softening: 0.00001 # Physical softening length (in internal units). + a_smooth: 0. + r_cut: 0. diff --git a/tests/testFFT.c b/tests/testFFT.c index ba737935afc4568a1509d2c43a3534b60a6ef867..7b67181ebd3e29bffbf564d00f702e6c15669fab 100644 --- a/tests/testFFT.c +++ b/tests/testFFT.c @@ -57,15 +57,18 @@ int main() { struct swift_params *params = malloc(sizeof(struct swift_params)); parser_read_file("fft_params.yml", params); + struct cosmology cosmo; + cosmology_init_no_cosmo(&cosmo); + /* Initialise the gravity properties */ struct gravity_props gravity_properties; - gravity_props_init(&gravity_properties, params); + gravity_props_init(&gravity_properties, params, &cosmo); /* Build the infrastructure */ struct space space; double dim[3] = {1., 1., 1.}; - space_init(&space, params, dim, NULL, gparts, NULL, 0, nr_gparts, 0, 1, 1, 1, - 0, 0); + space_init(&space, params, &cosmo, dim, NULL, gparts, NULL, 0, nr_gparts, 0, + 1, 1, 0, 1, 0, 0); struct engine engine; engine.s = &space; diff --git a/tests/testPotentialPair.c b/tests/testPotentialPair.c index 5fbfb04ecd0b4a8c05c07d622dd77bd38af970ae..53fc54ccdd63a9a9150b6701c1a76ac20af91d4c 100644 --- a/tests/testPotentialPair.c +++ b/tests/testPotentialPair.c @@ -106,6 +106,7 @@ int main() { props.a_smooth = 1.25; props.r_cut_min = 0.; props.theta_crit2 = 0.; + props.epsilon_cur = eps; e.gravity_properties = &props; struct runner r; @@ -178,7 +179,6 @@ int main() { gp->x[1] = 0.5; gp->x[2] = 0.5; gp->mass = 0.; - gp->epsilon = eps; gp->time_bin = 1; gp->type = swift_type_dark_matter; gp->id_or_neg_offset = n + 1; @@ -196,7 +196,6 @@ int main() { ci.gparts[0].x[1] = 0.5; ci.gparts[0].x[2] = 0.5; ci.gparts[0].mass = 1.; - ci.gparts[0].epsilon = eps; ci.gparts[0].time_bin = 1; ci.gparts[0].type = swift_type_dark_matter; ci.gparts[0].id_or_neg_offset = 1; @@ -211,11 +210,12 @@ int main() { for (int n = 0; n < num_tests; ++n) { const struct gpart *gp = &cj.gparts[n]; const struct gpart *gp2 = &ci.gparts[0]; + const double epsilon = gravity_get_softening(gp, &props); double pot_true = - potential(ci.gparts[0].mass, gp->x[0] - gp2->x[0], gp->epsilon, rlr); + potential(ci.gparts[0].mass, gp->x[0] - gp2->x[0], epsilon, rlr); double acc_true = - acceleration(ci.gparts[0].mass, gp->x[0] - gp2->x[0], gp->epsilon, rlr); + acceleration(ci.gparts[0].mass, gp->x[0] - gp2->x[0], epsilon, rlr); message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] - gp2->x[0], gp->a_grav[0], acc_true, gp->potential, pot_true); @@ -252,12 +252,12 @@ int main() { for (int n = 0; n < num_tests; ++n) { const struct gpart *gp = &cj.gparts[n]; const struct gravity_tensors *mpole = ci.multipole; + const double epsilon = gravity_get_softening(gp, &props); double pot_true = potential(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], - gp->epsilon, rlr * FLT_MAX); - double acc_true = - acceleration(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], gp->epsilon, - rlr * FLT_MAX); + epsilon, rlr * FLT_MAX); + double acc_true = acceleration( + mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], epsilon, rlr * FLT_MAX); message("x=%e f=%e f_true=%e pot=%e pot_true=%e", gp->x[0] - mpole->CoM[0], gp->a_grav[0], acc_true, gp->potential, pot_true); @@ -304,7 +304,6 @@ int main() { ci.gparts[n].mass = 1. / 8.; - ci.gparts[n].epsilon = eps; ci.gparts[n].time_bin = 1; ci.gparts[n].type = swift_type_dark_matter; ci.gparts[n].id_or_neg_offset = 1; @@ -333,18 +332,19 @@ int main() { for (int i = 0; i < 8; ++i) { const struct gpart *gp2 = &ci.gparts[i]; + const double epsilon = gravity_get_softening(gp, &props); const double dx[3] = {gp2->x[0] - gp->x[0], gp2->x[1] - gp->x[1], gp2->x[2] - gp->x[2]}; const double d = sqrt(dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]); - pot_true += potential(gp2->mass, d, gp->epsilon, rlr * FLT_MAX); + pot_true += potential(gp2->mass, d, epsilon, rlr * FLT_MAX); acc_true[0] -= - acceleration(gp2->mass, d, gp->epsilon, rlr * FLT_MAX) * dx[0] / d; + acceleration(gp2->mass, d, epsilon, rlr * FLT_MAX) * dx[0] / d; acc_true[1] -= - acceleration(gp2->mass, d, gp->epsilon, rlr * FLT_MAX) * dx[1] / d; + acceleration(gp2->mass, d, epsilon, rlr * FLT_MAX) * dx[1] / d; acc_true[2] -= - acceleration(gp2->mass, d, gp->epsilon, rlr * FLT_MAX) * dx[2] / d; + acceleration(gp2->mass, d, epsilon, rlr * FLT_MAX) * dx[2] / d; } message("x=%e f=%e f_true=%e pot=%e pot_true=%e %e %e", diff --git a/tests/testPotentialSelf.c b/tests/testPotentialSelf.c index 33fe574e34bcc49ff8bac6f19f22b55fda19f186..246139ba243a37dbc5f6cafbc488cb1bf6cdd8eb 100644 --- a/tests/testPotentialSelf.c +++ b/tests/testPotentialSelf.c @@ -103,6 +103,7 @@ int main() { struct gravity_props props; props.a_smooth = 1.25; + props.epsilon_cur = eps; e.gravity_properties = &props; struct runner r; @@ -139,7 +140,6 @@ int main() { c.gparts[0].x[1] = 0.5; c.gparts[0].x[2] = 0.5; c.gparts[0].mass = 1.; - c.gparts[0].epsilon = eps; c.gparts[0].time_bin = 1; c.gparts[0].type = swift_type_dark_matter; c.gparts[0].id_or_neg_offset = 1; @@ -156,7 +156,6 @@ int main() { gp->x[1] = 0.5; gp->x[2] = 0.5; gp->mass = 0.; - gp->epsilon = eps; gp->time_bin = 1; gp->type = swift_type_dark_matter; gp->id_or_neg_offset = n + 1; @@ -172,9 +171,10 @@ int main() { for (int n = 1; n < num_tests + 1; ++n) { const struct gpart *gp = &c.gparts[n]; - double pot_true = potential(c.gparts[0].mass, gp->x[0], gp->epsilon, rlr); - double acc_true = - acceleration(c.gparts[0].mass, gp->x[0], gp->epsilon, rlr); + const double epsilon = gravity_get_softening(gp, &props); + + double pot_true = potential(c.gparts[0].mass, gp->x[0], epsilon, rlr); + double acc_true = acceleration(c.gparts[0].mass, gp->x[0], epsilon, rlr); // message("x=%e f=%e f_true=%e", gp->x[0], gp->a_grav[0], acc_true); diff --git a/tests/testReading.c b/tests/testReading.c index f5e2757c2fe3bd507e877f134e8c768aba0ae645..ca1e0ef69078c5e384a9cd4eab1098923ce9f279 100644 --- a/tests/testReading.c +++ b/tests/testReading.c @@ -48,7 +48,8 @@ int main() { /* Read data */ read_ic_single("input.hdf5", &us, dim, &parts, &gparts, &sparts, &Ngas, - &Ngpart, &Nspart, &periodic, &flag_entropy_ICs, 1, 1, 0, 1, 0); + &Ngpart, &Nspart, &periodic, &flag_entropy_ICs, 1, 1, 0, 0, 1., + 1, 0); /* Check global properties read are correct */ assert(dim[0] == boxSize);