Commit eef5614e authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'periodic_gravity' into 'master'

Periodic gravity speed and accuracy improvements



See merge request !393
parents 59abd8af 9571b2a3
......@@ -23,6 +23,12 @@ Snapshots:
Statistics:
delta_time: 1e-2 # Time between statistics output
# Parameters for the self-gravity scheme
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
epsilon: 0.0001 # Softening length (in internal units).
theta: 0.7 # Opening angle (Multipole acceptance criterion)
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
......
......@@ -26,8 +26,8 @@ Statistics:
# Parameters for the self-gravity scheme
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
epsilon: 0.0001 # Softening length (in internal units).
theta: 0.7 # Opening angle (Multipole acceptance criterion)
# Parameters for the hydrodynamics scheme
SPH:
......
......@@ -27,8 +27,7 @@ Statistics:
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
epsilon: 0.0001 # Softening length (in internal units).
a_smooth: 1000.
r_cut: 4.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
# Parameters for the hydrodynamics scheme
SPH:
......
......@@ -23,6 +23,12 @@ Snapshots:
Statistics:
delta_time: 1e-2 # Time between statistics output
# Parameters for the self-gravity scheme
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
epsilon: 0.0001 # Softening length (in internal units).
theta: 0.7 # Opening angle (Multipole acceptance criterion)
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
......
......@@ -28,7 +28,7 @@ Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
epsilon: 0.00001 # Softening length (in internal units).
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-2 # Time between statistics output
......
......@@ -50,8 +50,8 @@ infile = args.input
TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair",
"init_grav", "ghost", "extra_ghost", "drift_part",
"drift_gpart", "kick1", "kick2", "timestep", "send", "recv",
"grav_top_level", "grav_long_range", "grav_mm", "grav_down",
"cooling", "sourceterms", "count"]
"grav_top_level", "grav_long_range", "grav_ghost", "grav_mm",
"grav_down", "cooling", "sourceterms", "count"]
SUBTYPES = ["none", "density", "gradient", "force", "grav", "external_grav",
"tend", "xv", "rho", "gpart", "multipole", "spart", "count"]
......
......@@ -51,11 +51,12 @@ SPH:
# Parameters for the self-gravity scheme
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
epsilon: 0.1 # Softening length (in internal units).
a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
r_cut: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
epsilon: 0.1 # Softening length (in internal units).
a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
r_cut_max: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
r_cut_min: 0.1 # (Optional) Cut-off in number of top-level cells below which no truncation of FMM forces are performed (this is the default value).
# Parameters related to the initial conditions
InitialConditions:
......
......@@ -91,17 +91,18 @@ pl.rcParams.update(PLOT_PARAMS)
TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair",
"init_grav", "ghost", "extra_ghost", "drift_part",
"drift_gpart", "kick1", "kick2", "timestep", "send", "recv",
"grav_top_level", "grav_long_range", "grav_mm", "grav_down",
"cooling", "sourceterms", "count"]
"grav_top_level", "grav_long_range", "grav_ghost", "grav_mm",
"grav_down", "cooling", "sourceterms", "count"]
SUBTYPES = ["none", "density", "gradient", "force", "grav", "external_grav",
"tend", "xv", "rho", "gpart", "multipole", "spart", "count"]
# Task/subtypes of interest.
FULLTYPES = ["self/force", "self/density", "self/grav", "sub_self/force",
"sub_self/density", "pair/force", "pair/density", "pair/grav",
"sub_pair/force",
"sub_pair/density", "recv/xv", "send/xv", "recv/rho", "send/rho",
"sub_self/density", "sub_self/grav", "pair/force", "pair/density",
"pair/grav", "sub_pair/force",
"sub_pair/density", "sub_pair/grav", "recv/xv", "send/xv",
"recv/rho", "send/rho",
"recv/tend", "send/tend"]
# A number of colours for the various types. Recycled when there are
......@@ -109,7 +110,7 @@ FULLTYPES = ["self/force", "self/density", "self/grav", "sub_self/force",
colours = ["cyan", "lightgray", "darkblue", "yellow", "tan", "dodgerblue",
"sienna", "aquamarine", "bisque", "blue", "green", "lightgreen",
"brown", "purple", "moccasin", "olivedrab", "chartreuse",
"darksage", "darkgreen", "green", "mediumseagreen",
"steelblue", "darkgreen", "green", "mediumseagreen",
"mediumaquamarine", "darkslategrey", "mediumturquoise",
"black", "cadetblue", "skyblue", "red", "slategray", "gold",
"slateblue", "blueviolet", "mediumorchid", "firebrick",
......
......@@ -36,4 +36,17 @@ __attribute__((always_inline)) INLINE static float approx_expf(float x) {
return 1.f + x * (1.f + x * (0.5f + x * (1.f / 6.f + 1.f / 24.f * x)));
}
/**
* @brief Approximate version of expf(x) using a 6th order Taylor expansion
*
*/
__attribute__((always_inline)) INLINE static float good_approx_expf(float x) {
return 1.f +
x * (1.f +
x * (0.5f +
x * ((1.f / 6.f) +
x * ((1.f / 24.f) +
x * ((1.f / 120.f) + (1.f / 720.f) * x)))));
}
#endif /* SWIFT_APPROX_MATH_H */
......@@ -1682,6 +1682,16 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
#endif
}
/**
* @brief Constructs the top-level tasks for the short-range gravity
* and long-range gravity interactions.
*
* - One FTT task per MPI rank.
* - Multiple gravity ghosts for dependencies.
* - All top-cells get a self task.
* - All pairs within range according to the multipole acceptance
* criterion get a pair task.
*/
void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
void *extra_data) {
......@@ -1692,6 +1702,7 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
struct scheduler *sched = &e->sched;
const int nodeID = e->nodeID;
const int periodic = s->periodic;
const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
const int cdim_ghost[3] = {s->cdim[0] / 4 + 1, s->cdim[1] / 4 + 1,
s->cdim[2] / 4 + 1};
......@@ -1747,8 +1758,9 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
if (cj->nodeID != nodeID) continue; // MATTHIEU
/* Are the cells to close for a MM interaction ? */
if (!gravity_multipole_accept(ci->multipole, cj->multipole,
theta_crit_inv, 1)) {
if (!gravity_multipole_accept_rebuild(ci->multipole, cj->multipole,
theta_crit_inv, periodic,
dim)) {
scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0,
ci, cj);
......@@ -1812,6 +1824,11 @@ void engine_make_self_gravity_tasks(struct engine *e) {
if (periodic) free(ghosts);
}
/**
* @brief Constructs the top-level tasks for the external gravity.
*
* @param e The #engine.
*/
void engine_make_external_gravity_tasks(struct engine *e) {
struct space *s = e->s;
......
......@@ -63,7 +63,7 @@ int gravity_exact_force_file_exits(const struct engine *e) {
char line[100];
char dummy1[10], dummy2[10];
double epsilon, newton_G;
int N;
int N, periodic;
/* Reads file header */
if (fgets(line, 100, file) != line) error("Problem reading title");
if (fgets(line, 100, file) != line) error("Problem reading G");
......@@ -72,10 +72,12 @@ int gravity_exact_force_file_exits(const struct engine *e) {
sscanf(line, "%s %s %d", dummy1, dummy2, &N);
if (fgets(line, 100, file) != line) error("Problem reading epsilon");
sscanf(line, "%s %s %le", dummy1, dummy2, &epsilon);
if (fgets(line, 100, file) != line) error("Problem reading BC");
sscanf(line, "%s %s %d", dummy1, dummy2, &periodic);
fclose(file);
/* Check whether it matches the current parameters */
if (N == SWIFT_GRAVITY_FORCE_CHECKS &&
if (N == SWIFT_GRAVITY_FORCE_CHECKS && periodic == e->s->periodic &&
(fabs(epsilon - e->gravity_properties->epsilon) / epsilon < 1e-5) &&
(fabs(newton_G - e->physical_constants->const_newton_G) / newton_G <
1e-5)) {
......@@ -101,6 +103,8 @@ void gravity_exact_force_compute_mapper(void *map_data, int nr_gparts,
struct exact_force_data *data = (struct exact_force_data *)extra_data;
const struct space *s = data->s;
const struct engine *e = data->e;
const int periodic = s->periodic;
const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const double const_G = data->const_G;
int counter = 0;
......@@ -124,11 +128,18 @@ void gravity_exact_force_compute_mapper(void *map_data, int nr_gparts,
if (gpi == gpj) continue;
/* Compute the pairwise distance. */
const double dx[3] = {gpi->x[0] - gpj->x[0], // x
gpi->x[1] - gpj->x[1], // y
gpi->x[2] - gpj->x[2]}; // z
const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
double dx = gpi->x[0] - gpj->x[0];
double dy = gpi->x[1] - gpj->x[1];
double dz = gpi->x[2] - gpj->x[2];
/* Now apply periodic BC */
if (periodic) {
dx = nearest(dx, dim[0]);
dy = nearest(dy, dim[1]);
dz = nearest(dz, dim[2]);
}
const double r2 = dx * dx + dy * dy + dz * dz;
const double r = sqrt(r2);
const double ir = 1. / r;
const double mj = gpj->mass;
......@@ -152,15 +163,11 @@ void gravity_exact_force_compute_mapper(void *map_data, int nr_gparts,
/* Get softened gravity */
f = mj * hi_inv3 * W * f_lr;
// printf("r=%e hi=%e W=%e fac=%e\n", r, hi, W, f);
}
const double fdx[3] = {f * dx[0], f * dx[1], f * dx[2]};
a_grav[0] -= fdx[0];
a_grav[1] -= fdx[1];
a_grav[2] -= fdx[2];
a_grav[0] -= f * dx;
a_grav[1] -= f * dy;
a_grav[2] -= f * dz;
}
/* Store the exact answer */
......@@ -245,8 +252,9 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
fprintf(file_swift, "# Gravity accuracy test - SWIFT FORCES\n");
fprintf(file_swift, "# G= %16.8e\n", e->physical_constants->const_newton_G);
fprintf(file_swift, "# N= %d\n", SWIFT_GRAVITY_FORCE_CHECKS);
fprintf(file_swift, "# epsilon=%16.8e\n", e->gravity_properties->epsilon);
fprintf(file_swift, "# theta=%16.8e\n", e->gravity_properties->theta_crit);
fprintf(file_swift, "# epsilon= %16.8e\n", e->gravity_properties->epsilon);
fprintf(file_swift, "# periodic= %d\n", s->periodic);
fprintf(file_swift, "# theta= %16.8e\n", e->gravity_properties->theta_crit);
fprintf(file_swift, "# Git Branch: %s\n", git_branch());
fprintf(file_swift, "# Git Revision: %s\n", git_revision());
fprintf(file_swift, "# %16s %16s %16s %16s %16s %16s %16s\n", "id", "pos[0]",
......
......@@ -28,11 +28,11 @@
#include "vector.h"
/**
* @brief Gravity forces between particles
* @brief Gravity forces between particles truncated by the long-range kernel
*/
__attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
float rlr_inv, float r2, const float *dx, struct gpart *gpi,
struct gpart *gpj) {
__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
float r2, const float *dx, struct gpart *gpi, struct gpart *gpj,
float rlr_inv) {
/* Apply the gravitational acceleration. */
const float r = sqrtf(r2);
......@@ -41,7 +41,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
const float mj = gpj->mass;
const float hi = gpi->epsilon;
const float hj = gpj->epsilon;
const float u = r * rlr_inv;
const float u_lr = r * rlr_inv;
float f_lr, fi, fj, W;
#ifdef SWIFT_DEBUG_CHECKS
......@@ -49,7 +49,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
#endif
/* Get long-range correction */
kernel_long_grav_eval(u, &f_lr);
kernel_long_grav_eval(u_lr, &f_lr);
if (r >= hi) {
......@@ -97,18 +97,84 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
}
/**
* @brief Gravity forces between particles (non-symmetric version)
* @brief Gravity forces between particles
*/
__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_nonsym(
float rlr_inv, float r2, const float *dx, struct gpart *gpi,
const struct gpart *gpj) {
__attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
float r2, const float *dx, struct gpart *gpi, struct gpart *gpj) {
/* Apply the gravitational acceleration. */
const float r = sqrtf(r2);
const float ir = 1.f / r;
const float mi = gpi->mass;
const float mj = gpj->mass;
const float hi = gpi->epsilon;
const float hj = gpj->epsilon;
float fi, fj, W;
#ifdef SWIFT_DEBUG_CHECKS
if (r == 0.f) error("Interacting particles with 0 distance");
#endif
if (r >= hi) {
/* Get Newtonian gravity */
fi = mj * ir * ir * ir;
} else {
const float hi_inv = 1.f / hi;
const float hi_inv3 = hi_inv * hi_inv * hi_inv;
const float ui = r * hi_inv;
kernel_grav_eval(ui, &W);
/* Get softened gravity */
fi = mj * hi_inv3 * W;
}
if (r >= hj) {
/* Get Newtonian gravity */
fj = mi * ir * ir * ir;
} else {
const float hj_inv = 1.f / hj;
const float hj_inv3 = hj_inv * hj_inv * hj_inv;
const float uj = r * hj_inv;
kernel_grav_eval(uj, &W);
/* Get softened gravity */
fj = mi * hj_inv3 * W;
}
const float fidx[3] = {fi * dx[0], fi * dx[1], fi * dx[2]};
gpi->a_grav[0] -= fidx[0];
gpi->a_grav[1] -= fidx[1];
gpi->a_grav[2] -= fidx[2];
const float fjdx[3] = {fj * dx[0], fj * dx[1], fj * dx[2]};
gpj->a_grav[0] += fjdx[0];
gpj->a_grav[1] += fjdx[1];
gpj->a_grav[2] += fjdx[2];
}
/**
* @brief Gravity forces between particles truncated by the long-range kernel
* (non-symmetric version)
*/
__attribute__((always_inline)) INLINE static void
runner_iact_grav_pp_truncated_nonsym(float r2, const float *dx,
struct gpart *gpi, const struct gpart *gpj,
float rlr_inv) {
/* Apply the gravitational acceleration. */
const float r = sqrtf(r2);
const float ir = 1.f / r;
const float mj = gpj->mass;
const float hi = gpi->epsilon;
const float u = r * rlr_inv;
const float u_lr = r * rlr_inv;
float f_lr, f, W;
#ifdef SWIFT_DEBUG_CHECKS
......@@ -116,7 +182,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_nonsym(
#endif
/* Get long-range correction */
kernel_long_grav_eval(u, &f_lr);
kernel_long_grav_eval(u_lr, &f_lr);
if (r >= hi) {
......@@ -143,13 +209,44 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_nonsym(
}
/**
* @brief Gravity forces between particle and multipole
* @brief Gravity forces between particles (non-symmetric version)
*/
__attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
float rlr_inv, float r2, const float *dx, struct gpart *gp,
const struct multipole *multi) {
__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_nonsym(
float r2, const float *dx, struct gpart *gpi, const struct gpart *gpj) {
/* Apply the gravitational acceleration. */
const float r = sqrtf(r2);
const float ir = 1.f / r;
const float mj = gpj->mass;
const float hi = gpi->epsilon;
float f, W;
error("Dead function");
#ifdef SWIFT_DEBUG_CHECKS
if (r == 0.f) error("Interacting particles with 0 distance");
#endif
if (r >= hi) {
/* Get Newtonian gravity */
f = mj * ir * ir * ir;
} else {
const float hi_inv = 1.f / hi;
const float hi_inv3 = hi_inv * hi_inv * hi_inv;
const float ui = r * hi_inv;
kernel_grav_eval(ui, &W);
/* Get softened gravity */
f = mj * hi_inv3 * W;
}
const float fdx[3] = {f * dx[0], f * dx[1], f * dx[2]};
gpi->a_grav[0] -= fdx[0];
gpi->a_grav[1] -= fdx[1];
gpi->a_grav[2] -= fdx[2];
}
#endif /* SWIFT_DEFAULT_GRAVITY_IACT_H */
......@@ -33,7 +33,8 @@
#include "kernel_gravity.h"
#define gravity_props_default_a_smooth 1.25f
#define gravity_props_default_r_cut 4.5f
#define gravity_props_default_r_cut_max 4.5f
#define gravity_props_default_r_cut_min 0.1f
void gravity_props_init(struct gravity_props *p,
const struct swift_params *params) {
......@@ -41,8 +42,10 @@ void gravity_props_init(struct gravity_props *p,
/* Tree-PM parameters */
p->a_smooth = parser_get_opt_param_float(params, "Gravity:a_smooth",
gravity_props_default_a_smooth);
p->r_cut = parser_get_opt_param_float(params, "Gravity:r_cut",
gravity_props_default_r_cut);
p->r_cut_max = parser_get_opt_param_float(params, "Gravity:r_cut_max",
gravity_props_default_r_cut_max);
p->r_cut_min = parser_get_opt_param_float(params, "Gravity:r_cut_min",
gravity_props_default_r_cut_min);
/* Time integration */
p->eta = parser_get_param_float(params, "Gravity:eta");
......@@ -69,9 +72,10 @@ void gravity_props_print(const struct gravity_props *p) {
message("Self-gravity softening: epsilon=%.4f (Plummer equivalent: %.4f)",
p->epsilon, p->epsilon / 3.);
message("Self-gravity MM smoothing-scale: a_smooth=%f", p->a_smooth);
message("Self-gravity mesh smoothing-scale: a_smooth=%f", p->a_smooth);
message("Self-gravity MM cut-off: r_cut=%f", p->r_cut);
message("Self-gravity tree cut-off: r_cut_max=%f", p->r_cut_max);
message("Self-gravity truncation cut-off: r_cut_min=%f", p->r_cut_min);
}
#if defined(HAVE_HDF5)
......@@ -84,7 +88,8 @@ void gravity_props_print_snapshot(hid_t h_grpgrav,
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);
io_write_attribute_f(h_grpgrav, "MM r_cut", p->r_cut);
io_write_attribute_f(h_grpgrav, "Mesh a_smooth", p->a_smooth);
io_write_attribute_f(h_grpgrav, "Mesh r_cut_max", p->r_cut_max);
io_write_attribute_f(h_grpgrav, "Mesh r_cut_min", p->r_cut_min);
}
#endif
......@@ -34,9 +34,16 @@
*/
struct gravity_props {
/* Tree-PM parameters */
/*! Mesh smoothing scale in units of top-level cell size */
float a_smooth;
float r_cut;
/*! Distance below which the truncated mesh force is Newtonian in units of
* a_smooth */
float r_cut_min;
/*! Distance above which the truncated mesh force is negligible in units of
* a_smooth */
float r_cut_max;
/*! Time integration dimensionless multiplier */
float eta;
......
......@@ -19,33 +19,67 @@
#ifndef SWIFT_KERNEL_LONG_GRAVITY_H
#define SWIFT_KERNEL_LONG_GRAVITY_H
#include <math.h>
/* Config parameters. */
#include "../config.h"
/* Includes. */
/* Local headers. */
#include "approx_math.h"
#include "const.h"
#include "inline.h"
#include "vector.h"
#define one_over_sqrt_pi ((float)(M_2_SQRTPI * 0.5))
/* Standard headers */
#include <math.h>
/**
* @brief Computes the long-range correction term for the FFT calculation.
*
* @param u The ratio of the distance to the FFT cell scale $u = x/A$.
* @param u The ratio of the distance to the FFT cell scale \f$u = r/r_s\f$.
* @param W (return) The value of the kernel function.
*/
__attribute__((always_inline)) INLINE static void kernel_long_grav_eval(
float u, float *const W) {
/* const float arg1 = u * 0.5f; */
/* const float arg2 = u * one_over_sqrt_pi; */
/* const float arg3 = -arg1 * arg1; */
#ifdef GADGET2_LONG_RANGE_CORRECTION
const float one_over_sqrt_pi = ((float)(M_2_SQRTPI * 0.5));
const float arg1 = u * 0.5f;
const float arg2 = u * one_over_sqrt_pi;
const float arg3 = -arg1 * arg1;
const float term1 = erfcf(arg1);
const float term2 = arg2 * expf(arg3);
*W = term1 + term2;
#else
const float arg = 2.f * u;
const float exp_arg = good_approx_expf(arg);
const float term = 1.f / (1.f + exp_arg);
/* const float term1 = erfcf(arg1); */
/* const float term2 = arg2 * expf(arg3); */
*W = arg * exp_arg * term * term - exp_arg * term + 1.f;
*W *= 2.f;
#endif
}
/**
* @brief Returns the long-range truncation of the Poisson potential in Fourier
* space.
*
* @param u2 The square of the Fourier mode times the cell scale
* \f$u^2 = k^2r_s^2\f$.
* @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) {
/* *W = term1 + term2; */
*W = 1.f;
#ifdef GADGET2_LONG_RANGE_CORRECTION
*W = exp(-u2);
#else
const double u = sqrt(u2);
const double arg = M_PI_2 * u;
*W = arg / sinh(arg);
#endif
}
#endif // SWIFT_KERNEL_LONG_GRAVITY_H
......@@ -1498,23 +1498,28 @@ INLINE static void gravity_M2M(struct multipole *m_a,
* @param pos_a The position of the multipole.
* @param props The #gravity_props of this calculation.
* @param periodic Is the calculation periodic ?