diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
index d698c66e418a120e5e7ebbe1cba0a5a9af8cd1f1..ea188f64a60627c776604618329bb205cad37f7c 100644
--- a/src/gravity_derivatives.h
+++ b/src/gravity_derivatives.h
@@ -33,6 +33,7 @@
/* Local headers. */
#include "inline.h"
#include "kernel_gravity.h"
+#include "kernel_long_gravity.h"
/**
* @brief Structure containing all the derivatives of the potential field
@@ -40,7 +41,7 @@
*/
struct potential_derivatives_M2L {
- /* 0th order terms */
+ /* 0th order term */
float D_000;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
@@ -124,11 +125,14 @@ struct potential_derivatives_M2P {
* @param r_inv Inverse norm of distance vector
* @param eps Softening length.
* @param eps_inv Inverse of softening length.
+ * @param periodic Is the calculation periodic ?
+ * @param rs_inv
* @param pot (return) The structure containing all the derivatives.
*/
__attribute__((always_inline)) INLINE static void
compute_potential_derivatives_M2L(float r_x, float r_y, float r_z, float r2,
float r_inv, float eps, float eps_inv,
+ int periodic, float rs_inv,
struct potential_derivatives_M2L *pot) {
float Dt_1;
@@ -148,8 +152,8 @@ compute_potential_derivatives_M2L(float r_x, float r_y, float r_z, float r2,
float Dt_11;
#endif
- /* Un-softened case */
- if (r2 > eps * eps) {
+ /* Un-softened un-truncated case (Newtonian potential) */
+ if (!periodic && r2 > eps * eps) {
Dt_1 = r_inv;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
@@ -171,6 +175,52 @@ compute_potential_derivatives_M2L(float r_x, float r_y, float r_z, float r2,
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
#error "Missing implementation for order >5"
#endif
+ /* Un-softened truncated case */
+ } else if (periodic && r2 > eps * eps) {
+
+ /* Get the derivatives of the truncated potential */
+ const float r = r2 * r_inv;
+ struct truncated_derivatives derivs;
+ kernel_long_grav_derivatives(r, rs_inv, &derivs);
+
+ Dt_1 = derivs.chi_0 * r_inv;
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ const float r_inv2 = r_inv * r_inv;
+ const float r_inv3 = r_inv2 * r_inv;
+ Dt_3 = (r * derivs.chi_1 - derivs.chi_0) * r_inv3;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+ const float r_inv5 = r_inv2 * r_inv3;
+ Dt_5 =
+ (r * r * derivs.chi_2 - 3.f * r * derivs.chi_1 + 3.f * derivs.chi_0) *
+ r_inv5;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+ const float r_inv7 = r_inv2 * r_inv5;
+ Dt_7 = (r * r * r * derivs.chi_3 - 6.f * r * r * derivs.chi_2 +
+ 15.f * r * derivs.chi_1 - 15.f * derivs.chi_0) *
+ r_inv7;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+ const float r_inv9 = r_inv2 * r_inv7;
+ Dt_9 = (r * r * r * r * derivs.chi_4 - 10.f * r * r * r * derivs.chi_3 +
+ 45.f * r * r * derivs.chi_2 - 105.f * r * derivs.chi_1 +
+ 105.f * derivs.chi_0) *
+ r_inv9;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
+ const float r_inv11 = r_inv2 * r_inv9;
+ Dt_11 = (r * r * r * r * r * derivs.chi_5 -
+ 15.f * r * r * r * r * derivs.chi_4 +
+ 105.f * r * r * r * derivs.chi_3 - 420.f * r * r * derivs.chi_2 +
+ 945.f * r * derivs.chi_1 - 945.f * derivs.chi_0) *
+ r_inv11;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
+#error "Missing implementation for order >5"
+#endif
+
+ /* Softened case */
} else {
const float r = r2 * r_inv;
diff --git a/src/kernel_long_gravity.h b/src/kernel_long_gravity.h
index e2e3117c51523b9c2fa330e9591ce28805d3a42a..f0e239cb9bfcbdfbba17779dc262bf88e23d42fd 100644
--- a/src/kernel_long_gravity.h
+++ b/src/kernel_long_gravity.h
@@ -28,10 +28,55 @@
#include "inline.h"
/* Standard headers */
+#include <float.h>
#include <math.h>
//#define GADGET2_LONG_RANGE_CORRECTION
+struct truncated_derivatives {
+
+ float chi_0;
+ float chi_1;
+ float chi_2;
+ float chi_3;
+ float chi_4;
+ float chi_5;
+};
+
+__attribute__((always_inline)) INLINE static void kernel_long_grav_derivatives(
+ const float r, const float rs_inv,
+ struct truncated_derivatives *const derivs) {
+
+ const float constant1 = 2.f * rs_inv;
+ const float x = constant1 * r;
+
+ const float exp_x = expf(x); // good_approx_expf(x);
+ const float alpha_inv = 1.f + exp_x;
+
+ const float alpha1 = 1.f / alpha_inv;
+ const float alpha2 = alpha1 * alpha1;
+ const float alpha3 = alpha2 * alpha1;
+ const float alpha4 = alpha3 * alpha1;
+ const float alpha5 = alpha4 * alpha1;
+ const float alpha6 = alpha5 * alpha1;
+
+ const float constant2 = constant1 * constant1;
+ const float constant3 = constant2 * constant1;
+ const float constant4 = constant3 * constant1;
+ const float constant5 = constant4 * constant1;
+
+ derivs->chi_0 = 2.f * (1.f - exp_x * alpha1);
+ derivs->chi_1 = constant1 * (2.f * alpha2 - 2.f * alpha1);
+ derivs->chi_2 = constant2 * (4.f * alpha3 - 6.f * alpha2 + 2.f * alpha1);
+ derivs->chi_3 = constant3 * (12.f * alpha4 - 24.f * alpha3 + 14.f * alpha2 -
+ 2.f * alpha1);
+ derivs->chi_4 = constant4 * (48.f * alpha5 - 120.f * alpha4 + 100.f * alpha3 -
+ 30.f * alpha2 + 2.f * alpha1);
+ derivs->chi_5 =
+ constant5 * (240.f * alpha6 - 720.f * alpha5 + 780.f * alpha4 -
+ 360.f * alpha3 + 62.f * alpha2 - 2.f * alpha1);
+}
+
/**
* @brief Computes the long-range correction term for the potential calculation
* coming from FFT.
@@ -68,7 +113,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_pot_eval(
* @param W (return) The value of the kernel function.
*/
__attribute__((always_inline)) INLINE static void kernel_long_grav_force_eval(
- float u, float *const W) {
+ const float u, float *const W) {
#ifdef GADGET2_LONG_RANGE_CORRECTION
@@ -112,7 +157,7 @@ __attribute__((always_inline)) INLINE static void kernel_long_grav_force_eval(
* @param W (return) The value of the kernel function.
*/
__attribute__((always_inline)) INLINE static void fourier_kernel_long_grav_eval(
- double u2, double *const W) {
+ const double u2, double *const W) {
#ifdef GADGET2_LONG_RANGE_CORRECTION
*W = exp(-u2);
diff --git a/src/multipole.h b/src/multipole.h
index d5756c3ac060fc03faa9e841653dab5c9768f352..26857c3b877ffec1c3070e707d938d682f5c0069 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -1521,9 +1521,7 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
const struct multipole *m_a,
const double pos_b[3], const double pos_a[3],
const struct gravity_props *props, int periodic,
- const double dim[3]) {
-
- error("ooo");
+ const double dim[3], float rs_inv) {
/* Recover some constants */
const float eps = props->epsilon_cur;
@@ -1547,7 +1545,8 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
/* Compute all derivatives */
struct potential_derivatives_M2L pot;
- compute_potential_derivatives_M2L(dx, dy, dz, r2, r_inv, eps, eps_inv, &pot);
+ compute_potential_derivatives_M2L(dx, dy, dz, r2, r_inv, eps, eps_inv,
+ periodic, rs_inv, &pot);
#ifdef SWIFT_DEBUG_CHECKS
/* Count interactions */
@@ -2264,7 +2263,7 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
#endif
// MATTHIEU
- return;
+ // return;
/* Local accumulator */
double a_grav[3] = {0., 0., 0.};
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 0849a34b3f0647a5245306f3f08e741728444cfb..58b5b40cd50572c0cf399c2aab79c3b3e407f7e9 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -174,8 +174,8 @@ static INLINE void runner_dopair_grav_mm(struct runner *r,
const int periodic = s->periodic;
const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const struct gravity_props *props = e->gravity_properties;
- // const float a_smooth = e->gravity_properties->a_smooth;
- // const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
+ const float a_smooth = e->mesh->a_smooth;
+ const float rlr_inv = 1. / a_smooth;
TIMER_TIC;
@@ -206,11 +206,8 @@ static INLINE void runner_dopair_grav_mm(struct runner *r,
cj->ti_old_multipole, cj->nodeID, ci->nodeID, e->ti_current);
/* Let's interact at this level */
- if (0)
- gravity_M2L(&ci->multipole->pot, multi_j, ci->multipole->CoM,
- cj->multipole->CoM, props, periodic, dim);
-
- runner_dopair_grav_pp(r, ci, cj, 0);
+ gravity_M2L(&ci->multipole->pot, multi_j, ci->multipole->CoM,
+ cj->multipole->CoM, props, periodic, dim, rlr_inv);
TIMER_TOC(timer_dopair_grav_mm);
}