Commit 6777c37c authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'gravity_multi_dt' into 'master'

Added softened gravity to the M2L kernels



See merge request !336
parents f973c37e 58a41f8e
......@@ -24,6 +24,39 @@
#include "gravity_properties.h"
#include "minmax.h"
/**
* @brief Returns the mass of a particle
*
* @param gp The particle of interest
*/
__attribute__((always_inline)) INLINE static float gravity_get_mass(
const struct gpart* restrict gp) {
return gp->mass;
}
/**
* @brief Returns the softening of a particle
*
* @param gp The particle of interest
*/
__attribute__((always_inline)) INLINE static float gravity_get_softening(
const struct gpart* restrict gp) {
return gp->epsilon;
}
/**
* @brief Returns the potential of a particle
*
* @param gp The particle of interest
*/
__attribute__((always_inline)) INLINE static float gravity_get_potential(
const struct gpart* restrict gp) {
return gp->potential;
}
/**
* @brief Computes the gravity time-step of a given particle due to self-gravity
*
......@@ -42,8 +75,10 @@ gravity_compute_timestep_self(const struct gpart* const gp,
const float ac_inv = (ac2 > 0.f) ? 1.f / sqrtf(ac2) : FLT_MAX;
/* Note that 0.714285714 = 2. (from Gadget) / 2.8 (Plummer softening) */
const float dt = sqrtf(0.714285714f * grav_props->eta * gp->epsilon * ac_inv);
const float epsilon = gravity_get_softening(gp);
/* Note that 0.66666667 = 2. (from Gadget) / 3. (Plummer softening) */
const float dt = sqrtf(0.66666667f * grav_props->eta * epsilon * ac_inv);
return dt;
}
......@@ -63,6 +98,7 @@ __attribute__((always_inline)) INLINE static void gravity_init_gpart(
gp->a_grav[0] = 0.f;
gp->a_grav[1] = 0.f;
gp->a_grav[2] = 0.f;
gp->potential = 0.f;
#ifdef SWIFT_DEBUG_CHECKS
gp->num_interacted = 0;
......@@ -130,8 +166,8 @@ __attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
__attribute__((always_inline)) INLINE static void gravity_init_softening(
struct gpart* gp, const struct gravity_props* grav_props) {
/* Note 2.8 is the Plummer-equivalent correction */
gp->epsilon = 2.8f * grav_props->epsilon;
/* Note 3 is the Plummer-equivalent correction */
gp->epsilon = grav_props->epsilon;
}
#endif /* SWIFT_DEFAULT_GRAVITY_H */
......@@ -41,6 +41,9 @@ struct gpart {
/* Particle mass. */
float mass;
/* Gravitational potential */
float potential;
/* Softening length */
float epsilon;
......
......@@ -27,8 +27,8 @@
* 1, 1, pp. 24 (2014), arXiv:1405.2255
*/
/* Some standard headers. */
#include <math.h>
/* Config parameters. */
#include "../config.h"
/* Local headers. */
#include "inline.h"
......
......@@ -52,7 +52,7 @@ void gravity_props_init(struct gravity_props *p,
p->theta_crit_inv = 1. / p->theta_crit;
/* Softening lengths */
p->epsilon = parser_get_param_double(params, "Gravity:epsilon");
p->epsilon = 3. * parser_get_param_double(params, "Gravity:epsilon");
p->epsilon2 = p->epsilon * p->epsilon;
p->epsilon_inv = 1. / p->epsilon;
}
......@@ -66,7 +66,8 @@ void gravity_props_print(const struct gravity_props *p) {
message("Self-gravity opening angle: theta=%.4f", p->theta_crit);
message("Self-gravity softening: epsilon=%.4f", p->epsilon);
message("Self-gravity softening: epsilon=%.4f (Plummer equivalent: %.4f)",
p->epsilon, p->epsilon / 3.);
if (p->a_smooth != gravity_props_default_a_smooth)
message("Self-gravity MM smoothing-scale: a_smooth=%f", p->a_smooth);
......@@ -81,6 +82,8 @@ void gravity_props_print_snapshot(hid_t h_grpgrav,
io_write_attribute_f(h_grpgrav, "Time integration eta", p->eta);
io_write_attribute_f(h_grpgrav, "Softening length", p->epsilon);
io_write_attribute_f(h_grpgrav, "Softening length (Plummer equivalent)",
p->epsilon / 3.);
io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit);
io_write_attribute_d(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
io_write_attribute_f(h_grpgrav, "MM a_smooth", p->a_smooth);
......
......@@ -27,30 +27,17 @@
* 1, 1, pp. 24 (2014), arXiv:1405.2255
*/
/* Some standard headers. */
#include <math.h>
/* Config parameters. */
#include "../config.h"
/* Local headers. */
#include "inline.h"
#include "kernel_gravity.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$.
*
......@@ -71,19 +58,6 @@ __attribute__((always_inline)) INLINE static double D_soft_000(
/* 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$.
*
......@@ -136,18 +110,6 @@ __attribute__((always_inline)) INLINE static double D_soft_001(
/* 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$.
*
......@@ -161,8 +123,10 @@ __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);
const double eps_inv2 = eps_inv * eps_inv;
const double eps_inv3 = eps_inv * eps_inv2;
const double eps_inv5 = eps_inv3 * eps_inv2;
return r_x * r_x * eps_inv5 * D_soft_2(u) - eps_inv3 * D_soft_1(u);
}
/**
......@@ -178,8 +142,10 @@ __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);
const double eps_inv2 = eps_inv * eps_inv;
const double eps_inv3 = eps_inv * eps_inv2;
const double eps_inv5 = eps_inv3 * eps_inv2;
return r_y * r_y * eps_inv5 * D_soft_2(u) - eps_inv3 * D_soft_1(u);
}
/**
......@@ -195,8 +161,10 @@ __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);
const double eps_inv2 = eps_inv * eps_inv;
const double eps_inv3 = eps_inv * eps_inv2;
const double eps_inv5 = eps_inv3 * eps_inv2;
return r_z * r_z * eps_inv5 * D_soft_2(u) - eps_inv3 * D_soft_1(u);
}
/**
......@@ -213,8 +181,9 @@ __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);
const double eps_inv2 = eps_inv * eps_inv;
const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
return r_x * r_y * eps_inv5 * D_soft_2(u);
}
/**
......@@ -231,8 +200,9 @@ __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);
const double eps_inv2 = eps_inv * eps_inv;
const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
return r_x * r_z * eps_inv5 * D_soft_2(u);
}
/**
......@@ -249,8 +219,225 @@ __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);
const double eps_inv2 = eps_inv * eps_inv;
const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
return r_y * r_z * eps_inv5 * D_soft_2(u);
}
/*************************/
/* 3rd order derivatives */
/*************************/
/**
* @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z, h)}{\partial r_x^3} \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_300(
double r_x, double r_y, double r_z, double r, double eps_inv) {
const double u = r * eps_inv;
const double eps_inv2 = eps_inv * eps_inv;
const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
const double eps_inv7 = eps_inv5 * eps_inv2;
return -r_x * r_x * r_x * eps_inv7 * D_soft_3(u) +
3. * r_x * eps_inv5 * D_soft_2(u);
}
/**
* @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z, h)}{\partial r_y^3} \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_030(
double r_x, double r_y, double r_z, double r, double eps_inv) {
const double u = r * eps_inv;
const double eps_inv2 = eps_inv * eps_inv;
const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
const double eps_inv7 = eps_inv5 * eps_inv2;
return -r_y * r_y * r_y * eps_inv7 * D_soft_3(u) +
3. * r_y * eps_inv5 * D_soft_2(u);
}
/**
* @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z, h)}{\partial r_z^3} \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_003(
double r_x, double r_y, double r_z, double r, double eps_inv) {
const double u = r * eps_inv;
const double eps_inv2 = eps_inv * eps_inv;
const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
const double eps_inv7 = eps_inv5 * eps_inv2;
return -r_z * r_z * r_z * eps_inv7 * D_soft_3(u) +
3. * r_z * eps_inv5 * D_soft_2(u);
}
/**
* @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z, h)}{\partial r_x^2\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_210(
double r_x, double r_y, double r_z, double r, double eps_inv) {
const double u = r * eps_inv;
const double eps_inv2 = eps_inv * eps_inv;
const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
const double eps_inv7 = eps_inv5 * eps_inv2;
return -r_x * r_x * r_y * eps_inv7 * D_soft_3(u) +
r_y * eps_inv5 * D_soft_2(u);
}
/**
* @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z, h)}{\partial r_x^2\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_201(
double r_x, double r_y, double r_z, double r, double eps_inv) {
const double u = r * eps_inv;
const double eps_inv2 = eps_inv * eps_inv;
const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
const double eps_inv7 = eps_inv5 * eps_inv2;
return -r_x * r_x * r_z * eps_inv7 * D_soft_3(u) +
r_z * eps_inv5 * D_soft_2(u);
}
/**
* @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z, h)}{\partial r_x\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_120(
double r_x, double r_y, double r_z, double r, double eps_inv) {
const double u = r * eps_inv;
const double eps_inv2 = eps_inv * eps_inv;
const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
const double eps_inv7 = eps_inv5 * eps_inv2;
return -r_x * r_y * r_y * eps_inv7 * D_soft_3(u) +
r_x * eps_inv5 * D_soft_2(u);
}
/**
* @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z, h)}{\partial r_y^2\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_021(
double r_x, double r_y, double r_z, double r, double eps_inv) {
const double u = r * eps_inv;
const double eps_inv2 = eps_inv * eps_inv;
const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
const double eps_inv7 = eps_inv5 * eps_inv2;
return -r_y * r_y * r_z * eps_inv7 * D_soft_3(u) +
r_z * eps_inv5 * D_soft_2(u);
}
/**
* @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z, h)}{\partial r_x\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_102(
double r_x, double r_y, double r_z, double r, double eps_inv) {
const double u = r * eps_inv;
const double eps_inv2 = eps_inv * eps_inv;
const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
const double eps_inv7 = eps_inv5 * eps_inv2;
return -r_x * r_z * r_z * eps_inv7 * D_soft_3(u) +
r_x * eps_inv5 * D_soft_2(u);
}
/**
* @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z, h)}{\partial r_y\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_012(
double r_x, double r_y, double r_z, double r, double eps_inv) {
const double u = r * eps_inv;
const double eps_inv2 = eps_inv * eps_inv;
const double eps_inv5 = eps_inv2 * eps_inv2 * eps_inv;
const double eps_inv7 = eps_inv5 * eps_inv2;
return -r_y * r_z * r_z * eps_inv7 * D_soft_3(u) +
r_y * eps_inv5 * D_soft_2(u);
}
/**
* @brief \f$ \frac{\partial^3\phi(r_x, r_y, r_z, h)}{\partial r_z\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_111(
double r_x, double r_y, double r_z, double r, double eps_inv) {
const double u = r * eps_inv;
const double eps_inv2 = eps_inv * eps_inv;
const double eps_inv4 = eps_inv2 * eps_inv2;
const double eps_inv7 = eps_inv4 * eps_inv2 * eps_inv;
return -r_x * r_y * r_z * eps_inv7 * D_soft_3(u);
}
#endif /* SWIFT_GRAVITY_SOFTENED_DERIVATIVE_H */
......@@ -23,10 +23,8 @@
#include <math.h>
/* Includes. */
#include "const.h"
#include "inline.h"
#include "minmax.h"
#include "vector.h"
/**
* @brief Computes the gravity softening function.
......@@ -39,13 +37,12 @@
__attribute__((always_inline)) INLINE static void kernel_grav_eval(
float u, float *const W) {
/* W(u) = 21u^6 - 90u^5 + 140u^4 - 84u^3 + 14u */
/* W(u) = 21u^5 - 90u^4 + 140u^3 - 84u^2 + 14 */
*W = 21.f * u - 90.f;
*W = *W * u + 140.f;
*W = *W * u - 84.f;
*W = *W * u;
*W = *W * u + 14.f;
*W = *W * u;
}
#ifdef SWIFT_GRAVITY_FORCE_CHECKS
......@@ -61,14 +58,59 @@ __attribute__((always_inline)) INLINE static void kernel_grav_eval(
__attribute__((always_inline)) INLINE static void kernel_grav_eval_double(
double u, double *const W) {
/* W(u) = 21u^6 - 90u^5 + 140u^4 - 84u^3 + 14u */
/* W(u) = 21u^5 - 90u^4 + 140u^3 - 84u^2 + 14 */
*W = 21. * u - 90.;
*W = *W * u + 140.;
*W = *W * u - 84.;
*W = *W * u;
*W = *W * u + 14.;
*W = *W * u;
}
#endif /* SWIFT_GRAVITY_FORCE_CHECKS */
/************************************************/
/* Derivatives of softening kernel used for FMM */
/************************************************/
__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;
}
__attribute__((always_inline)) INLINE static double D_soft_1(double u) {
/* phi'(u)/u = 21u^5 - 90u^4 + 140u^3 - 84u^2 + 14 */
double phi = 21. * u - 90.;
phi = phi * u + 140.;
phi = phi * u - 84.;
phi = phi * u;
phi = phi * u + 14.;
return phi;
}
__attribute__((always_inline)) INLINE static double D_soft_2(double u) {
/* (phi'(u)/u)'/u = -105u^3 + 360u^2 - 420u + 168 */
double phi = -105. * u + 360.;
phi = phi * u - 420.;
phi = phi * u + 168.;
return phi;
}
__attribute__((always_inline)) INLINE static double D_soft_3(double u) {
/* ((phi'(u)/u)'/u)'/u = 315u - 720 + 420/u */
return 315. * u - 720. + 420. / u;
}
#endif /* SWIFT_KERNEL_GRAVITY_H */
......@@ -1626,7 +1626,7 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
m_a->M_010 * D_021(dx, dy, dz, r_inv) +
m_a->M_001 * D_012(dx, dy, dz, r_inv);
/* 3rd order multipole term (addition to rank 2)*/
/* 3rd order multipole term (addition to rank 3)*/
l_b->F_300 += m_a->M_000 * D_300(dx, dy, dz, r_inv);
l_b->F_030 += m_a->M_000 * D_030(dx, dy, dz, r_inv);
l_b->F_003 += m_a->M_000 * D_003(dx, dy, dz, r_inv);
......@@ -2096,6 +2096,72 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
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
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
/* 3rd order multipole term (addition to rank 0)*/
l_b->F_000 += m_a->M_300 * D_soft_300(dx, dy, dz, r, eps_inv) +
m_a->M_030 * D_soft_030(dx, dy, dz, r, eps_inv) +
m_a->M_003 * D_soft_003(dx, dy, dz, r, eps_inv);
l_b->F_000 += m_a->M_210 * D_soft_210(dx, dy, dz, r, eps_inv) +
m_a->M_201 * D_soft_201(dx, dy, dz, r, eps_inv) +
m_a->M_120 * D_soft_120(dx, dy, dz, r, eps_inv);
l_b->F_000 += m_a->M_021 * D_soft_021(dx, dy, dz, r, eps_inv) +
m_a->M_102 * D_soft_102(dx, dy, dz, r, eps_inv) +
m_a->M_012 * D_soft_012(dx, dy, dz, r, eps_inv);
l_b->F_000 += m_a->M_111 * D_soft_111(dx, dy, dz, r, eps_inv);
/* 3rd order multipole term (addition to rank 1)*/
l_b->F_100 += m_a->M_200 * D_soft_300(dx, dy, dz, r, eps_inv) +
m_a->M_020 * D_soft_120(dx, dy, dz, r, eps_inv) +
m_a->M_002 * D_soft_102(dx, dy, dz, r, eps_inv);
l_b->F_100 += m_a->M_110 * D_soft_210(dx, dy, dz, r, eps_inv) +
m_a->M_101 * D_soft_201(dx, dy, dz, r, eps_inv) +
m_a->M_011 * D_soft_111(dx, dy, dz, r, eps_inv);
l_b->F_010 += m_a->M_200 * D_soft_210(dx, dy, dz, r, eps_inv) +
m_a->M_020 * D_soft_030(dx, dy, dz, r, eps_inv) +
m_a->M_002 * D_soft_012(dx, dy, dz, r, eps_inv);