From d9d8742a1b58dbd811c28b0607dce0a15e07bc32 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Mon, 18 Jun 2018 15:04:55 +0200 Subject: [PATCH] Correct P-M forces for multipole order < 3 --- src/gravity/Default/gravity_iact.h | 8 +-- src/kernel_long_gravity.h | 8 +-- src/runner_doiact_grav.h | 26 ++++----- tests/testPotentialPair.c | 90 +++++++++++++++++++++--------- tests/testPotentialSelf.c | 2 +- 5 files changed, 84 insertions(+), 50 deletions(-) diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h index 8f90db42ae..e73e1b55fe 100644 --- a/src/gravity/Default/gravity_iact.h +++ b/src/gravity/Default/gravity_iact.h @@ -152,7 +152,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full( /* In the case where the order is < 3, then there is only a monopole term left. * We can default to the normal P-P interaction with the mass of the multipole * and its CoM as the "particle" property */ -#if 1 //SELF_GRAVITY_MULTIPOLE_ORDER < 3 +#if 1 // SELF_GRAVITY_MULTIPOLE_ORDER < 3 float f_ij, pot_ij; runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000, @@ -240,9 +240,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated( float f_ij, pot_ij; runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000, r_s_inv, &f_ij, &pot_ij); - *f_x = -f_ij * r_x; - *f_y = -f_ij * r_y; - *f_z = -f_ij * r_z; + *f_x = f_ij * r_x; + *f_y = f_ij * r_y; + *f_z = f_ij * r_z; *pot = pot_ij; #else diff --git a/src/kernel_long_gravity.h b/src/kernel_long_gravity.h index f11f760022..2c14c0d888 100644 --- a/src/kernel_long_gravity.h +++ b/src/kernel_long_gravity.h @@ -31,7 +31,7 @@ #include <float.h> #include <math.h> -#define GADGET2_LONG_RANGE_CORRECTION +//#define GADGET2_LONG_RANGE_CORRECTION #ifdef GADGET2_LONG_RANGE_CORRECTION #define kernel_long_gravity_truncation_name "Gadget-like (using erfc())" @@ -121,7 +121,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_derivatives( const float x = c1 * r; /* e^(2r / r_s) */ - const float exp_x = good_approx_expf(x); + const float exp_x = expf(x); // good_approx_expf(x); /* 1 / alpha(w) */ const float a_inv = 1.f + exp_x; @@ -164,7 +164,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_pot_eval( #else const float x = 2.f * u; - const float exp_x = good_approx_expf(x); + const float exp_x = expf(x); // good_approx_expf(x); const float alpha = 1.f / (1.f + exp_x); /* We want 2 - 2 exp(x) * alpha */ @@ -197,7 +197,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_force_eval( #else const float x = 2.f * u; - const float exp_x = good_approx_expf(x); + const float exp_x = expf(x); // good_approx_expf(x); const float alpha = 1.f / (1.f + exp_x); /* We want 2*(x*alpha - x*alpha^2 - exp(x)*alpha + 1) */ diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h index 3d5aa78bf1..2c12d81223 100644 --- a/src/runner_doiact_grav.h +++ b/src/runner_doiact_grav.h @@ -751,10 +751,9 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci, ci->gparts, cj->gparts); /* Then the M2P */ - if (0) - runner_dopair_grav_pm_truncated(ci_cache, gcount_padded_i, CoM_j, - multi_j, dim, r_s_inv, e, ci->gparts, - gcount_i, cj); + runner_dopair_grav_pm_truncated(ci_cache, gcount_padded_i, CoM_j, + multi_j, dim, r_s_inv, e, ci->gparts, + gcount_i, cj); } if (cj_active && symmetric) { @@ -764,10 +763,9 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci, cj->gparts, ci->gparts); /* Then the M2P */ - if (0) - runner_dopair_grav_pm_truncated(cj_cache, gcount_padded_j, CoM_i, - multi_i, dim, r_s_inv, e, cj->gparts, - gcount_j, ci); + runner_dopair_grav_pm_truncated(cj_cache, gcount_padded_j, CoM_i, + multi_i, dim, r_s_inv, e, cj->gparts, + gcount_j, ci); } } else { @@ -783,10 +781,8 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci, ci->gparts, cj->gparts); /* Then the M2P */ - if (0) - runner_dopair_grav_pm_full(ci_cache, gcount_padded_i, CoM_j, multi_j, - periodic, dim, e, ci->gparts, gcount_i, - cj); + runner_dopair_grav_pm_full(ci_cache, gcount_padded_i, CoM_j, multi_j, + periodic, dim, e, ci->gparts, gcount_i, cj); } if (cj_active && symmetric) { @@ -796,10 +792,8 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci, cj->gparts, ci->gparts); /* Then the M2P */ - if (0) - runner_dopair_grav_pm_full(cj_cache, gcount_padded_j, CoM_i, multi_i, - periodic, dim, e, cj->gparts, gcount_j, - ci); + runner_dopair_grav_pm_full(cj_cache, gcount_padded_j, CoM_i, multi_i, + periodic, dim, e, cj->gparts, gcount_j, ci); } } } diff --git a/tests/testPotentialPair.c b/tests/testPotentialPair.c index b282606ffe..b180efe159 100644 --- a/tests/testPotentialPair.c +++ b/tests/testPotentialPair.c @@ -42,7 +42,7 @@ const double eps = 0.1; */ void check_value(double a, double b, const char *s) { if (fabs(a - b) / fabs(a + b) > 2e-6 && fabs(a - b) > 1.e-6) - error("Values are inconsistent: %12.15e %12.15e (%s)!", a, b, s); + error("Values are inconsistent: SWIFT:%12.15e true:%12.15e (%s)!", a, b, s); } /* Definitions of the potential and force that match @@ -92,25 +92,34 @@ int main(int argc, char *argv[]) { #ifdef HAVE_FE_ENABLE_EXCEPT feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); #endif - + /* Initialise a few things to get us going */ + + /* Non-truncated forces first */ + double rlr = FLT_MAX; + struct engine e; e.max_active_bin = num_time_bins; e.time = 0.1f; e.ti_current = 8; e.time_base = 1e-10; + struct space s; + s.periodic = 0; + e.s = &s; + struct pm_mesh mesh; mesh.periodic = 0; mesh.dim[0] = 10.; mesh.dim[1] = 10.; mesh.dim[2] = 10.; - mesh.r_s_inv = 0.; + mesh.r_s = rlr; + mesh.r_s_inv = 1. / rlr; mesh.r_cut_min = 0.; + mesh.r_cut_max = FLT_MAX; e.mesh = &mesh; struct gravity_props props; - props.a_smooth = 1.25; props.theta_crit2 = 0.; props.epsilon_cur = eps; e.gravity_properties = &props; @@ -119,8 +128,6 @@ int main(int argc, char *argv[]) { bzero(&r, sizeof(struct runner)); r.e = &e; - const double rlr = FLT_MAX; - /* Init the cache for gravity interaction */ gravity_cache_init(&r.ci_gravity_cache, num_tests); gravity_cache_init(&r.cj_gravity_cache, num_tests); @@ -225,8 +232,8 @@ int main(int argc, char *argv[]) { double acc_true = 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); + /* 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); */ check_value(gp->potential, pot_true, "potential"); check_value(gp->a_grav[0], acc_true, "acceleration"); @@ -262,13 +269,13 @@ int main(int argc, char *argv[]) { 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], - epsilon, rlr * FLT_MAX); - double acc_true = acceleration( - mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], epsilon, rlr * FLT_MAX); + double pot_true = + potential(mpole->m_pole.M_000, gp->x[0] - mpole->CoM[0], epsilon, rlr); + double acc_true = acceleration(mpole->m_pole.M_000, + gp->x[0] - mpole->CoM[0], epsilon, rlr); - 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); + /* 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); */ check_value(gp->potential, pot_true, "potential"); check_value(gp->a_grav[0], acc_true, "acceleration"); @@ -279,11 +286,47 @@ int main(int argc, char *argv[]) { /* Reset the accelerations */ for (int n = 0; n < num_tests; ++n) gravity_init_gpart(&cj.gparts[n]); -/***************************************/ -/* Test the high-order PM interactions */ -/***************************************/ + /***************************************/ + /* Test the truncated PM interactions */ + /***************************************/ + rlr = 2.; + mesh.r_s = rlr; + mesh.r_s_inv = 1. / rlr; + mesh.periodic = 1; + s.periodic = 1; + props.epsilon_cur = FLT_MIN; /* No softening */ + + /* Now compute the forces */ + runner_dopair_grav_pp(&r, &ci, &cj, 1); + + /* Verify everything */ + 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], epsilon, rlr); + double acc_true = acceleration(mpole->m_pole.M_000, + gp->x[0] - mpole->CoM[0], epsilon, rlr); + + /* 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); */ + + check_value(gp->potential, pot_true, "potential"); + check_value(gp->a_grav[0], acc_true, "acceleration"); + } + + message("\n\t\t truncated P-M interactions all good\n"); + + /***************************************/ + /* Test the high-order PM interactions */ + /***************************************/ + + /* Reset the accelerations */ + for (int n = 0; n < num_tests; ++n) gravity_init_gpart(&cj.gparts[n]); -#if SELF_GRAVITY_MULTIPOLE_ORDER >= 3 +#if 0 // SELF_GRAVITY_MULTIPOLE_ORDER >= 3 /* Let's make ci more interesting */ free(ci.gparts); @@ -346,13 +389,10 @@ int main(int argc, char *argv[]) { 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, epsilon, rlr * FLT_MAX); - acc_true[0] -= - acceleration(gp2->mass, d, epsilon, rlr * FLT_MAX) * dx[0] / d; - acc_true[1] -= - acceleration(gp2->mass, d, epsilon, rlr * FLT_MAX) * dx[1] / d; - acc_true[2] -= - acceleration(gp2->mass, d, epsilon, rlr * FLT_MAX) * dx[2] / d; + pot_true += potential(gp2->mass, d, epsilon, rlr); + acc_true[0] -= acceleration(gp2->mass, d, epsilon, rlr) * dx[0] / d; + acc_true[1] -= acceleration(gp2->mass, d, epsilon, rlr) * dx[1] / d; + acc_true[2] -= acceleration(gp2->mass, d, epsilon, rlr) * 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 9d57770c86..6c68825689 100644 --- a/tests/testPotentialSelf.c +++ b/tests/testPotentialSelf.c @@ -95,7 +95,7 @@ int main(int argc, char *argv[]) { #ifdef HAVE_FE_ENABLE_EXCEPT feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW); #endif - + /* Initialise a few things to get us going */ struct engine e; e.max_active_bin = num_time_bins; -- GitLab