diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
index 337ae67f4620677e7a3dfdd1b0c527eb016504e2..d34fc2884645f6176f8f689bf1d276755f1b659b 100644
--- a/src/gravity_derivatives.h
+++ b/src/gravity_derivatives.h
@@ -21,7 +21,7 @@
/**
* @file gravity_derivatives.h
- * @brief Derivatives (up to 3rd order) of the gravitational potential.
+ * @brief Derivatives (up to 5th order) of the gravitational potential.
*
* We use the notation of Dehnen, Computational Astrophysics and Cosmology,
* 1, 1, pp. 24 (2014), arXiv:1405.2255
diff --git a/src/gravity_softened_derivatives.h b/src/gravity_softened_derivatives.h
new file mode 100644
index 0000000000000000000000000000000000000000..2ae927fa9008dbec69df0df413bf8770c16bd68f
--- /dev/null
+++ b/src/gravity_softened_derivatives.h
@@ -0,0 +1,256 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GRAVITY_SOFTENED_DERIVATIVE_H
+#define SWIFT_GRAVITY_SOFTENED_DERIVATIVE_H
+
+/**
+ * @file gravity_softened_derivatives.h
+ * @brief Derivatives of the softened gravitational potential.
+ *
+ * We use the notation of Dehnen, Computational Astrophysics and Cosmology,
+ * 1, 1, pp. 24 (2014), arXiv:1405.2255
+ */
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local headers. */
+#include "inline.h"
+
+/*************************/
+/* 0th order derivatives */
+/*************************/
+
+__attribute__((always_inline)) INLINE static double D_soft_0(double u) {
+
+ /* phi(u) = 3u^7 - 15u^6 + 28u^5 - 21u^4 + 7u^2 - 3 */
+ double phi = 3. * u - 15.;
+ phi = phi * u + 28.;
+ phi = phi * u - 21.;
+ phi = phi * u;
+ phi = phi * u + 7.;
+ phi = phi * u;
+ phi = phi * u - 3.;
+
+ return phi;
+}
+
+/**
+ * @brief \f$ \phi(r_x, r_y, r_z, h) \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r Norm of the distance vector (\f$ |r| \f$).
+ * @param eps_inv Inverse of the softening length (\f$ 1/h \f$).
+ */
+__attribute__((always_inline)) INLINE static double D_soft_000(
+ double r_x, double r_y, double r_z, double r, double eps_inv) {
+
+ const double u = r * eps_inv;
+ return eps_inv * D_soft_0(u);
+}
+
+/*************************/
+/* 1st order derivatives */
+/*************************/
+
+__attribute__((always_inline)) INLINE static double D_soft_1(double u) {
+
+ /* phi(u) = 21u^6 - 90u^5 + 140u^4 - 84u^3 + 14u */
+ double phi = 21. * u - 90.;
+ phi = phi * u + 140.;
+ phi = phi * u - 84.;
+ phi = phi * u;
+ phi = phi * u + 14.;
+ phi = phi * u;
+
+ return phi;
+}
+
+/**
+ * @brief \f$ \frac{\partial\phi(r_x, r_y, r_z, h)}{\partial r_x} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r Norm of the distance vector (\f$ |r| \f$).
+ * @param eps_inv Inverse of the softening length (\f$ 1/h \f$).
+ */
+__attribute__((always_inline)) INLINE static double D_soft_100(
+ double r_x, double r_y, double r_z, double r, double eps_inv) {
+
+ const double u = r * eps_inv;
+ return -r_x * eps_inv * eps_inv * eps_inv * D_soft_1(u);
+}
+
+/**
+ * @brief \f$ \frac{\partial\phi(r_x, r_y, r_z, h)}{\partial r_x} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r Norm of the distance vector (\f$ |r| \f$).
+ * @param eps_inv Inverse of the softening length (\f$ 1/h \f$).
+ */
+__attribute__((always_inline)) INLINE static double D_soft_010(
+ double r_x, double r_y, double r_z, double r, double eps_inv) {
+
+ const double u = r * eps_inv;
+ return -r_y * eps_inv * eps_inv * eps_inv * D_soft_1(u);
+}
+
+/**
+ * @brief \f$ \frac{\partial\phi(r_x, r_y, r_z, h)}{\partial r_x} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r Norm of the distance vector (\f$ |r| \f$).
+ * @param eps_inv Inverse of the softening length (\f$ 1/h \f$).
+ */
+__attribute__((always_inline)) INLINE static double D_soft_001(
+ double r_x, double r_y, double r_z, double r, double eps_inv) {
+
+ const double u = r * eps_inv;
+ return -r_z * eps_inv * eps_inv * eps_inv * D_soft_1(u);
+}
+
+/*************************/
+/* 2nd order derivatives */
+/*************************/
+
+__attribute__((always_inline)) INLINE static double D_soft_2(double u) {
+
+ /* phi(u) = 126u^5 - 450u^4 + 560u^3 - 252u^2 + 14 */
+ double phi = 126. * u - 450.;
+ phi = phi * u + 560.;
+ phi = phi * u - 252.;
+ phi = phi * u;
+ phi = phi * u + 14.;
+
+ return phi;
+}
+
+/**
+ * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z, h)}{\partial r_x^2} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r Norm of the distance vector (\f$ |r| \f$).
+ * @param eps_inv Inverse of the softening length (\f$ 1/h \f$).
+ */
+__attribute__((always_inline)) INLINE static double D_soft_200(
+ double r_x, double r_y, double r_z, double r, double eps_inv) {
+
+ const double u = r * eps_inv;
+ return r_x * r_x * eps_inv * eps_inv * D_soft_2(u) +
+ (r_y * r_y + r_z * r_z) * eps_inv * eps_inv * eps_inv * D_soft_1(u);
+}
+
+/**
+ * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z, h)}{\partial r_y^2} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r Norm of the distance vector (\f$ |r| \f$).
+ * @param eps_inv Inverse of the softening length (\f$ 1/h \f$).
+ */
+__attribute__((always_inline)) INLINE static double D_soft_020(
+ double r_x, double r_y, double r_z, double r, double eps_inv) {
+
+ const double u = r * eps_inv;
+ return r_y * r_y * eps_inv * eps_inv * D_soft_2(u) +
+ (r_x * r_x + r_z * r_z) * eps_inv * eps_inv * eps_inv * D_soft_1(u);
+}
+
+/**
+ * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z, h)}{\partial r_z^2} \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r Norm of the distance vector (\f$ |r| \f$).
+ * @param eps_inv Inverse of the softening length (\f$ 1/h \f$).
+ */
+__attribute__((always_inline)) INLINE static double D_soft_002(
+ double r_x, double r_y, double r_z, double r, double eps_inv) {
+
+ const double u = r * eps_inv;
+ return r_z * r_z * eps_inv * eps_inv * D_soft_2(u) +
+ (r_x * r_x + r_y * r_y) * eps_inv * eps_inv * eps_inv * D_soft_1(u);
+}
+
+/**
+ * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z, h)}{\partial r_x\partial r_y}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r Norm of the distance vector (\f$ |r| \f$).
+ * @param eps_inv Inverse of the softening length (\f$ 1/h \f$).
+ */
+__attribute__((always_inline)) INLINE static double D_soft_110(
+ double r_x, double r_y, double r_z, double r, double eps_inv) {
+
+ const double u = r * eps_inv;
+ return r_x * r_y * eps_inv * eps_inv * D_soft_2(u) -
+ r_x * r_y * eps_inv * eps_inv * eps_inv * D_soft_1(u);
+}
+
+/**
+ * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z, h)}{\partial r_x\partial r_z}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r Norm of the distance vector (\f$ |r| \f$).
+ * @param eps_inv Inverse of the softening length (\f$ 1/h \f$).
+ */
+__attribute__((always_inline)) INLINE static double D_soft_101(
+ double r_x, double r_y, double r_z, double r, double eps_inv) {
+
+ const double u = r * eps_inv;
+ return r_x * r_z * eps_inv * eps_inv * D_soft_2(u) -
+ r_x * r_z * eps_inv * eps_inv * eps_inv * D_soft_1(u);
+}
+
+/**
+ * @brief \f$ \frac{\partial^2\phi(r_x, r_y, r_z, h)}{\partial r_y\partial r_z}
+ * \f$.
+ *
+ * @param r_x x-coordinate of the distance vector (\f$ r_x \f$).
+ * @param r_y y-coordinate of the distance vector (\f$ r_y \f$).
+ * @param r_z z-coordinate of the distance vector (\f$ r_z \f$).
+ * @param r Norm of the distance vector (\f$ |r| \f$).
+ * @param eps_inv Inverse of the softening length (\f$ 1/h \f$).
+ */
+__attribute__((always_inline)) INLINE static double D_soft_011(
+ double r_x, double r_y, double r_z, double r, double eps_inv) {
+
+ const double u = r * eps_inv;
+ return r_y * r_z * eps_inv * eps_inv * D_soft_2(u) -
+ r_y * r_z * eps_inv * eps_inv * eps_inv * D_soft_1(u);
+}
+
+#endif /* SWIFT_GRAVITY_SOFTENED_DERIVATIVE_H */
diff --git a/src/multipole.h b/src/multipole.h
index 24a526f5758db8cf11fffa0b8362c5e24e045c1b..ea4105d37a54b30d82c722e89a416509f2032f83 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -34,6 +34,7 @@
#include "error.h"
#include "gravity_derivatives.h"
#include "gravity_properties.h"
+#include "gravity_softened_derivatives.h"
#include "inline.h"
#include "kernel_gravity.h"
#include "minmax.h"
@@ -1505,6 +1506,9 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
const struct gravity_props *props,
int periodic) {
+ /* Recover some constants */
+ const double eps2 = props->epsilon2;
+
/* Compute distance vector */
const double dx =
periodic ? box_wrap(pos_b[0] - pos_a[0], 0., 1.) : pos_b[0] - pos_a[0];
@@ -1523,7 +1527,7 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
#endif
/* Un-softened case */
- if (r2 > props->epsilon2) {
+ if (r2 > eps2) {
/* 0th order term */
l_b->F_000 += m_a->M_000 * D_000(dx, dy, dz, r_inv);
@@ -2045,7 +2049,54 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
/* Softened case */
} else {
- message("Un-supported softened M2L interaction.");
+
+ const double eps_inv = props->epsilon_inv;
+ const double r = r2 * r_inv;
+
+ /* 0th order term */
+ l_b->F_000 += m_a->M_000 * D_soft_000(dx, dy, dz, r, eps_inv);
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+
+ /* 1st order multipole term (addition to rank 0)*/
+ l_b->F_000 += m_a->M_100 * D_soft_100(dx, dy, dz, r, eps_inv) +
+ m_a->M_010 * D_soft_010(dx, dy, dz, r, eps_inv) +
+ m_a->M_001 * D_soft_001(dx, dy, dz, r, eps_inv);
+
+ /* 1st order multipole term (addition to rank 1)*/
+ l_b->F_100 += m_a->M_000 * D_soft_100(dx, dy, dz, r, eps_inv);
+ l_b->F_010 += m_a->M_000 * D_soft_010(dx, dy, dz, r, eps_inv);
+ l_b->F_001 += m_a->M_000 * D_soft_001(dx, dy, dz, r, eps_inv);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+
+ /* 2nd order multipole term (addition to rank 0)*/
+ l_b->F_000 += m_a->M_200 * D_soft_200(dx, dy, dz, r, eps_inv) +
+ m_a->M_020 * D_soft_020(dx, dy, dz, r, eps_inv) +
+ m_a->M_002 * D_soft_002(dx, dy, dz, r, eps_inv);
+ l_b->F_000 += m_a->M_110 * D_soft_110(dx, dy, dz, r, eps_inv) +
+ m_a->M_101 * D_soft_101(dx, dy, dz, r, eps_inv) +
+ m_a->M_011 * D_soft_011(dx, dy, dz, r, eps_inv);
+
+ /* 2nd order multipole term (addition to rank 1)*/
+ l_b->F_100 += m_a->M_100 * D_soft_200(dx, dy, dz, r, eps_inv) +
+ m_a->M_010 * D_soft_110(dx, dy, dz, r, eps_inv) +
+ m_a->M_001 * D_soft_101(dx, dy, dz, r, eps_inv);
+ l_b->F_010 += m_a->M_100 * D_soft_110(dx, dy, dz, r, eps_inv) +
+ m_a->M_010 * D_soft_020(dx, dy, dz, r, eps_inv) +
+ m_a->M_001 * D_soft_011(dx, dy, dz, r, eps_inv);
+ l_b->F_001 += m_a->M_100 * D_soft_101(dx, dy, dz, r, eps_inv) +
+ m_a->M_010 * D_soft_011(dx, dy, dz, r, eps_inv) +
+ m_a->M_001 * D_soft_002(dx, dy, dz, r, eps_inv);
+
+ /* 2nd order multipole term (addition to rank 2)*/
+ l_b->F_200 += m_a->M_000 * D_soft_200(dx, dy, dz, r, eps_inv);
+ l_b->F_020 += m_a->M_000 * D_soft_020(dx, dy, dz, r, eps_inv);
+ l_b->F_002 += m_a->M_000 * D_soft_002(dx, dy, dz, r, eps_inv);
+ l_b->F_110 += m_a->M_000 * D_soft_110(dx, dy, dz, r, eps_inv);
+ l_b->F_101 += m_a->M_000 * D_soft_101(dx, dy, dz, r, eps_inv);
+ l_b->F_011 += m_a->M_000 * D_soft_011(dx, dy, dz, r, eps_inv);
+#endif
}
}
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 8a75182323859cb5662b0f7d0e83bd0f25f78b4e..10343c05db2d6d1133ccdc847c6774c169414d83 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -23,7 +23,7 @@ AM_LDFLAGS = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS) $(FFTW_LIBS
# List of programs and scripts to run in the test suite
TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetry \
testPair.sh testPairPerturbed.sh test27cells.sh test27cellsPerturbed.sh \
- testParser.sh testSPHStep test125cells.sh testKernelGrav testFFT \
+ testParser.sh testSPHStep test125cells.sh testFFT \
testAdiabaticIndex testRiemannExact testRiemannTRRS testRiemannHLLC \
testMatrixInversion testThreadpool testDump testLogger \
testVoronoi1D testVoronoi2D testVoronoi3D
@@ -31,7 +31,7 @@ TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetr
# List of test programs to compile
check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
testSPHStep testPair test27cells test125cells testParser \
- testKernel testKernelGrav testFFT testInteractions testMaths \
+ testKernel testFFT testInteractions testMaths \
testSymmetry testThreadpool benchmarkInteractions \
testAdiabaticIndex testRiemannExact testRiemannTRRS \
testRiemannHLLC testMatrixInversion testDump testLogger \
@@ -62,8 +62,6 @@ testParser_SOURCES = testParser.c
testKernel_SOURCES = testKernel.c
-testKernelGrav_SOURCES = testKernelGrav.c
-
testFFT_SOURCES = testFFT.c
testInteractions_SOURCES = testInteractions.c