From ebb5d65aa00411ceb2c0d44c5f935739b10dedea Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Fri, 8 May 2020 17:02:41 +0200
Subject: [PATCH] Impose a mininum r/H in the softened gravity derivatives for
M2P and M2L
---
src/gravity_derivatives.h | 18 ++++--------------
1 file changed, 4 insertions(+), 14 deletions(-)
diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
index 40caa4a436..45dcd24051 100644
--- a/src/gravity_derivatives.h
+++ b/src/gravity_derivatives.h
@@ -220,11 +220,6 @@ potential_derivatives_compute_M2L(const float r_x, const float r_y,
const int periodic, const float r_s_inv,
struct potential_derivatives_M2L *pot) {
-#ifdef SWIFT_DEBUG_CHECKS
-// if (r2 < 0.99f * eps * eps)
-// error("Computing M2L derivatives below softening length");
-#endif
-
float Dt_1;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
float Dt_3;
@@ -247,8 +242,8 @@ potential_derivatives_compute_M2L(const float r_x, const float r_y,
const float eps_inv = 1.f / eps;
const float r = r2 * r_inv;
- const float u = r * eps_inv;
- const float u_inv = r_inv * eps;
+ const float u = max(r * eps_inv, 0.01f);
+ const float u_inv = min(r_inv * eps, 100.f);
Dt_1 = -eps_inv * D_soft_1(u, u_inv);
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
@@ -484,11 +479,6 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
const int periodic, const float r_s_inv,
struct potential_derivatives_M2P *pot) {
-#ifdef SWIFT_DEBUG_CHECKS
-// if (r2 < 0.99f * eps * eps)
-// error("Computing M2P derivatives below softening length");
-#endif
-
float Dt_1;
float Dt_3;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
@@ -509,8 +499,8 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
const float eps_inv = 1.f / eps;
const float r = r2 * r_inv;
- const float u = r * eps_inv;
- const float u_inv = r_inv * eps;
+ const float u = max(r * eps_inv, 0.01f);
+ const float u_inv = min(r_inv * eps, 100.f);
Dt_1 = -eps_inv * D_soft_1(u, u_inv);
const float eps_inv2 = eps_inv * eps_inv;
--
GitLab