From 5ed1998c3f4c179ebdda49272ea39d0de16d6f86 Mon Sep 17 00:00:00 2001 From: Matthieu Schaller <schaller@strw.leidenuniv.nl> Date: Sat, 25 Apr 2020 14:16:32 +0200 Subject: [PATCH] Check for multipoles of size 0 in M2L --- src/multipole_accept.h | 24 +++++++++++++++++------- 1 file changed, 17 insertions(+), 7 deletions(-) diff --git a/src/multipole_accept.h b/src/multipole_accept.h index cf78bb3bb0..17027339f5 100644 --- a/src/multipole_accept.h +++ b/src/multipole_accept.h @@ -56,6 +56,11 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept( const float rho_A = use_rebuild_sizes ? A->r_max_rebuild : A->r_max; const float rho_B = use_rebuild_sizes ? B->r_max_rebuild : B->r_max; +#ifdef SWIFT_DEBUG_CHECKS + if (rho_A == 0.) error("Size of multipole A is 0!"); + if (rho_B == 0.) error("Size of multipole B is 0!"); +#endif + /* Compute the error estimator (without the 1/M_B term that cancels out) */ float E_BA_term = 0.f; for (int n = 0; n <= p; ++n) { @@ -63,12 +68,10 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept( binomial(p, n) * B->m_pole.power[n] * integer_powf(rho_A, p - n); } E_BA_term *= 8.f; - if (rho_A != 0. && rho_B != 0.) { - E_BA_term *= max(rho_A, rho_B); - E_BA_term /= (rho_A + rho_B); - } + E_BA_term *= max(rho_A, rho_B); + E_BA_term /= (rho_A + rho_B); - /* Compute r^(p+2) */ + /* Compute r^(p+2) */ #if SELF_GRAVITY_MULTIPOLE_ORDER % 2 == 1 const float r_to_p_plus2 = integer_powf(sqrtf(r2), (p + 2)); #else @@ -164,6 +167,13 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept( /* Order of the expansion */ const int p = SELF_GRAVITY_MULTIPOLE_ORDER; + /* Sizes of the multipoles */ + const float rho_B = B->r_max; + +#ifdef SWIFT_DEBUG_CHECKS + if (rho_B == 0.) error("Size of multipole B is 0!"); +#endif + /* Compute the error estimator (without the 1/M_B term that cancels out) */ const float E_BA_term = 8.f * B->m_pole.power[p]; @@ -196,7 +206,7 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept( /* Test the different conditions */ /* Condition 1: We are in the converging part of the Taylor expansion */ - const int cond_1 = (B->r_max) * (B->r_max) < r2; + const int cond_1 = rho_B * rho_B < r2; /* Condition 2: We are not below softening */ const int cond_2 = max_softening * max_softening < r2; @@ -210,7 +220,7 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept( } else { /* Condition 1: We are obeying the purely geometric criterion */ - const int cond_1 = (B->r_max) * (B->r_max) < theta_crit2 * r2; + const int cond_1 = rho_B * rho_B < theta_crit2 * r2; /* Condition 2: We are not below softening */ const int cond_2 = max_softening * max_softening < r2; -- GitLab