diff --git a/examples/main.c b/examples/main.c
index b3e85987d753d1f54c01dc68e64efce12d34d4f5..5f7ac3ff9e7442348185f6ba81036613eae1c81b 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -1076,7 +1076,7 @@ int main(int argc, char *argv[]) {
gravity_props_init(&gravity_properties, params, &prog_const, &cosmo,
with_cosmology, with_external_gravity,
with_baryon_particles, with_DM_particles,
- with_DM_background_particles, periodic);
+ with_DM_background_particles, periodic, s.dim);
/* Initialise the external potential properties */
bzero(&potential, sizeof(struct external_potential));
diff --git a/examples/main_fof.c b/examples/main_fof.c
index aa5de58a5616282399d82716a4733e9289539b1e..6204d2d9caa9284b7c5e762f305deb7d989bcb7e 100644
--- a/examples/main_fof.c
+++ b/examples/main_fof.c
@@ -526,7 +526,7 @@ int main(int argc, char *argv[]) {
gravity_props_init(&gravity_properties, params, &prog_const, &cosmo,
/*with_cosmology=*/1, /*with_external_gravity=*/0,
with_baryon_particles, with_DM_particles,
- with_DM_background_particles, periodic);
+ with_DM_background_particles, periodic, s.dim);
/* Initialise the long-range gravity mesh */
if (periodic) {
diff --git a/src/gravity_properties.c b/src/gravity_properties.c
index 1c8d2bbfaee6b8525932058006ad9fadc270ba26..b1c1bbac0696df04cbc94e36fc9233d02365e575 100644
--- a/src/gravity_properties.c
+++ b/src/gravity_properties.c
@@ -45,7 +45,8 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
const struct cosmology *cosmo, const int with_cosmology,
const int with_external_potential,
const int has_baryons, const int has_DM,
- const int is_zoom_simulation, const int periodic) {
+ const int is_zoom_simulation, const int periodic,
+ const double dim[3]) {
/* Tree updates */
p->rebuild_frequency =
@@ -65,6 +66,9 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
p->r_cut_min_ratio = parser_get_opt_param_float(
params, "Gravity:r_cut_min", gravity_props_default_r_cut_min);
+ const double r_s = p->a_smooth * dim[0] / p->mesh_size;
+ p->r_s_inv = 1. / r_s;
+
/* Some basic checks of what we read */
if (p->mesh_size % 2 != 0)
error("The mesh side-length must be an even number.");
@@ -101,7 +105,7 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
}
/* We always start with the geometric MAC */
- p->use_advanced_mac = 0;
+ p->use_advanced_MAC = 0;
/* Geometric opening angle */
p->theta_crit = parser_get_param_double(params, "Gravity:theta_cr");
@@ -112,6 +116,11 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
p->adaptive_tolerance =
parser_get_param_float(params, "Gravity:epsilon_fmm");
+ /* Consider truncated forces in the MAC? */
+ if (p->use_adaptive_tolerance)
+ p->consider_truncation_in_MAC =
+ parser_get_param_int(params, "Gravity:allow_truncation_in_MAC");
+
/* Are we allowing tree use below softening? */
p->use_tree_below_softening =
parser_get_opt_param_int(params, "Gravity:use_tree_below_softening", 1);
@@ -206,7 +215,7 @@ void gravity_props_update_MAC_choice(struct gravity_props *p) {
/* Now that we have run initial accelerations,
* switch to the better MAC */
- if (p->use_adaptive_tolerance) p->use_advanced_mac = 1;
+ if (p->use_adaptive_tolerance) p->use_advanced_MAC = 1;
}
void gravity_props_update(struct gravity_props *p,
diff --git a/src/gravity_properties.h b/src/gravity_properties.h
index 2941e224921c1c7557d32533491875fd5608bb37..161c5008468a9a9819982401e281bb8151695945 100644
--- a/src/gravity_properties.h
+++ b/src/gravity_properties.h
@@ -55,7 +55,7 @@ struct gravity_props {
/* -------------- Properties of the FFM gravity ---------------------- */
/*! What MAC are we currently using? */
- int use_advanced_mac;
+ int use_advanced_MAC;
/*! Are we using the adaptive opening angle? (as read from param file) */
int use_adaptive_tolerance;
@@ -69,6 +69,9 @@ struct gravity_props {
/*! Are we allowing tree gravity below softening? */
int use_tree_below_softening;
+ /*! Are we applying long-range truncation to the forces in the MAC? */
+ int consider_truncation_in_MAC;
+
/* ------------- Properties of the softened gravity ------------------ */
/*! Co-moving softening length for for high-res. DM particles */
@@ -113,12 +116,17 @@ struct gravity_props {
* a_smooth */
float r_cut_max_ratio;
+ /*! Inverse of the long-range gravity mesh scale. */
+ float r_s_inv;
+
/*! Are we dithering the particles at every rebuild? */
int with_dithering;
/*! Fraction of the top-level cell size used to normalize the dithering */
double dithering_ratio;
+ /* ------------- Physical constants ---------------------------------- */
+
/*! Gravitational constant (in internal units, copied from the physical
* constants) */
float G_Newton;
@@ -130,7 +138,8 @@ void gravity_props_init(struct gravity_props *p, struct swift_params *params,
const struct cosmology *cosmo, const int with_cosmology,
const int with_external_potential,
const int has_baryons, const int has_DM,
- const int is_zoom_simulation, const int periodic);
+ const int is_zoom_simulation, const int periodic,
+ const double dim[3]);
void gravity_props_update(struct gravity_props *p,
const struct cosmology *cosmo);
void gravity_props_update_MAC_choice(struct gravity_props *p);
diff --git a/src/mesh_gravity.c b/src/mesh_gravity.c
index bbd4496112114277f650582432799b5743422a14..bd5eb6a385198782cc39ec52a996a435837005d7 100644
--- a/src/mesh_gravity.c
+++ b/src/mesh_gravity.c
@@ -726,7 +726,7 @@ void pm_mesh_free(struct pm_mesh* mesh) {
* @param nr_threads The number of threads on this MPI rank.
*/
void pm_mesh_init(struct pm_mesh* mesh, const struct gravity_props* props,
- double dim[3], int nr_threads) {
+ const double dim[3], int nr_threads) {
#ifdef HAVE_FFTW
diff --git a/src/mesh_gravity.h b/src/mesh_gravity.h
index 3c05eb04d99b828b07af9e524897d82f30c3bb82..79c4d1b619cb3f73bc8aa39e97f1c5b2f6386386 100644
--- a/src/mesh_gravity.h
+++ b/src/mesh_gravity.h
@@ -68,7 +68,7 @@ struct pm_mesh {
};
void pm_mesh_init(struct pm_mesh *mesh, const struct gravity_props *props,
- double dim[3], int nr_threads);
+ const double dim[3], int nr_threads);
void pm_mesh_init_no_mesh(struct pm_mesh *mesh, double dim[3]);
void pm_mesh_compute_potential(struct pm_mesh *mesh, const struct space *s,
struct threadpool *tp, int verbose);
diff --git a/src/multipole_accept.h b/src/multipole_accept.h
index 4cdee6131870af0e736971e0f5412391e9fd5fa6..9a3f482763fa370c4036cd87c54d76da1d9bd708 100644
--- a/src/multipole_accept.h
+++ b/src/multipole_accept.h
@@ -25,6 +25,7 @@
/* Local includes */
#include "binomial.h"
#include "integer_power.h"
+#include "kernel_long_gravity.h"
#include "minmax.h"
#include "multipole_struct.h"
@@ -74,12 +75,24 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept(
E_BA_term /= (rho_A + rho_B);
}
- /* Compute r^(p+2) */
+ float r_to_p_plus2, W;
+ if (periodic && props->consider_truncation_in_MAC) {
+
+ /* Compute r^(p+2) and the long-range correction */
+ const float r = sqrtf(r2);
+ r_to_p_plus2 = integer_powf(r, (p + 2));
+ W = kernel_long_grav_force_eval(r * props->r_s_inv);
+
+ } else {
+
+ /* Compute r^(p+2) */
#if SELF_GRAVITY_MULTIPOLE_ORDER % 2 == 1
- const float r_to_p_plus2 = integer_powf(sqrtf(r2), (p + 2));
+ r_to_p_plus2 = integer_powf(sqrtf(r2), (p + 2));
#else
- const float r_to_p_plus2 = integer_powf(r2, ((p / 2) + 1));
+ r_to_p_plus2 = integer_powf(r2, ((p / 2) + 1));
#endif
+ W = 1.f;
+ }
/* Get the mimimal acceleration in A */
const float min_a_grav = A->m_pole.min_old_a_grav_norm;
@@ -98,7 +111,7 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept(
/* Get the sum of the multipole sizes */
const float rho_sum = rho_A + rho_B;
- if (props->use_advanced_mac) {
+ if (props->use_advanced_MAC) {
#ifdef SWIFT_DEBUG_CHECKS
if (min_a_grav == 0.) error("Acceleration is 0");
@@ -114,8 +127,8 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2L_accept(
props->use_tree_below_softening || max_softening * max_softening < r2;
/* Condition 3: The contribution is accurate enough
- * (E_BA / r^(p+2) < eps * a_min) */
- const int cond_3 = E_BA_term < eps * min_a_grav * r_to_p_plus2;
+ * (E_BA * (1 / r^(p)) * ((1 / r^2) * W) < eps * a_min) */
+ const int cond_3 = E_BA_term * W < eps * min_a_grav * r_to_p_plus2;
return cond_1 && cond_2 && cond_3;
@@ -185,12 +198,24 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept(
/* Compute the error estimator (without the 1/M_B term that cancels out) */
const float E_BA_term = 8.f * B->m_pole.power[p];
+ float r_to_p_plus2, W;
+ if (periodic && props->consider_truncation_in_MAC) {
+
+ /* Compute r^(p+2) and the long-range correction */
+ const float r = sqrtf(r2);
+ r_to_p_plus2 = integer_powf(r, (p + 2));
+ W = kernel_long_grav_force_eval(r * props->r_s_inv);
+
+ } else {
+
/* Compute r^(p+2) */
#if SELF_GRAVITY_MULTIPOLE_ORDER % 2 == 1
- const float r_to_p_plus2 = integer_powf(sqrtf(r2), (p + 2));
+ r_to_p_plus2 = integer_powf(sqrtf(r2), (p + 2));
#else
- const float r_to_p_plus2 = integer_powf(r2, ((p / 2) + 1));
+ r_to_p_plus2 = integer_powf(r2, ((p / 2) + 1));
#endif
+ W = 1.f;
+ }
/* Get the estimate of the acceleration */
const float old_a_grav = pa->old_a_grav_norm;
@@ -206,7 +231,7 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept(
const float theta_crit = props->theta_crit;
const float theta_crit2 = theta_crit * theta_crit;
- if (props->use_advanced_mac) {
+ if (props->use_advanced_MAC) {
#ifdef SWIFT_DEBUG_CHECKS
if (old_a_grav == 0.) error("Acceleration is 0");
@@ -222,8 +247,8 @@ __attribute__((nonnull, pure)) INLINE static int gravity_M2P_accept(
props->use_tree_below_softening || max_softening * max_softening < r2;
/* Condition 3: The contribution is accurate enough
- * (E_BA / r^(p+2) < eps * a) */
- const int cond_3 = E_BA_term < eps * old_a_grav * r_to_p_plus2;
+ * (E_BA * (1 / r^(p)) * ((1 / r^2) * W) < eps * a_min) */
+ const int cond_3 = E_BA_term * W < eps * old_a_grav * r_to_p_plus2;
return cond_1 && cond_2 && cond_3;