From 303c62a2873d6d0fd9da083cc9b42ee0abfb96db Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Mon, 18 Jun 2018 16:57:03 +0200
Subject: [PATCH] Add a parameter to the gravity pair P-P interaction function
to forbid the use of P2M interactions
---
.../ZeldovichPancake_3D/zeldovichPancake.yml | 2 +-
src/approx_math.h | 8 ++-
src/gravity_cache.h | 12 ++--
src/gravity_properties.c | 3 +-
src/kernel_long_gravity.h | 2 +-
src/multipole.h | 4 +-
src/runner_doiact_grav.h | 60 +++++++++++--------
7 files changed, 53 insertions(+), 38 deletions(-)
diff --git a/examples/ZeldovichPancake_3D/zeldovichPancake.yml b/examples/ZeldovichPancake_3D/zeldovichPancake.yml
index ca415706ee..3f82af2c1f 100644
--- a/examples/ZeldovichPancake_3D/zeldovichPancake.yml
+++ b/examples/ZeldovichPancake_3D/zeldovichPancake.yml
@@ -46,7 +46,7 @@ Cosmology:
Gravity:
mesh_side_length: 16
eta: 0.025
- theta: 0.5
+ theta: 0.7
r_cut_max: 6.
comoving_softening: 0.0001
max_physical_softening: 0.0001
diff --git a/src/approx_math.h b/src/approx_math.h
index aa9ed2b4ef..1445a84ab0 100644
--- a/src/approx_math.h
+++ b/src/approx_math.h
@@ -32,7 +32,7 @@
*
* @param x The number to take the exponential of.
*/
-__attribute__((always_inline)) INLINE static float approx_expf(float x) {
+__attribute__((always_inline, const)) 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)));
}
@@ -40,7 +40,8 @@ __attribute__((always_inline)) INLINE static float approx_expf(float x) {
* @brief Approximate version of expf(x) using a 6th order Taylor expansion
*
*/
-__attribute__((always_inline)) INLINE static float good_approx_expf(float x) {
+__attribute__((always_inline, const)) INLINE static float good_approx_expf(
+ float x) {
return 1.f +
x * (1.f +
x * (0.5f +
@@ -52,7 +53,8 @@ __attribute__((always_inline)) INLINE static float good_approx_expf(float x) {
/**
* @brief Approximate version of exp(x) using a 6th order Taylor expansion
*/
-__attribute__((always_inline)) INLINE static double good_approx_exp(double x) {
+__attribute__((always_inline, const)) INLINE static double good_approx_exp(
+ double x) {
return 1. +
x * (1. +
x * (0.5 +
diff --git a/src/gravity_cache.h b/src/gravity_cache.h
index 7eebc396b1..5dc65dabaf 100644
--- a/src/gravity_cache.h
+++ b/src/gravity_cache.h
@@ -171,6 +171,7 @@ __attribute__((always_inline)) INLINE static void gravity_cache_zero_output(
* more expensive P2P.
*
* @param max_active_bin The largest active bin in the current time-step.
+ * @param allow_mpole Are we allowing the use of multipoles?
* @param periodic Are we using periodic BCs ?
* @param dim The size of the simulation volume along each dimension.
* @param c The #gravity_cache to fill.
@@ -185,10 +186,11 @@ __attribute__((always_inline)) INLINE static void gravity_cache_zero_output(
* @param grav_props The global gravity properties.
*/
__attribute__((always_inline)) INLINE static void gravity_cache_populate(
- const timebin_t max_active_bin, const int periodic, const float dim[3],
- struct gravity_cache *c, const struct gpart *restrict gparts,
- const int gcount, const int gcount_padded, const double shift[3],
- const float CoM[3], const float r_max2, const struct cell *cell,
+ const timebin_t max_active_bin, const int allow_mpole, const int periodic,
+ const float dim[3], struct gravity_cache *c,
+ const struct gpart *restrict gparts, const int gcount,
+ const int gcount_padded, const double shift[3], const float CoM[3],
+ const float r_max2, const struct cell *cell,
const struct gravity_props *grav_props) {
const float theta_crit2 = grav_props->theta_crit2;
@@ -227,7 +229,7 @@ __attribute__((always_inline)) INLINE static void gravity_cache_populate(
const float r2 = dx * dx + dy * dy + dz * dz;
/* Check whether we can use the multipole instead of P-P */
- use_mpole[i] = gravity_M2P_accept(r_max2, theta_crit2, r2);
+ use_mpole[i] = allow_mpole && gravity_M2P_accept(r_max2, theta_crit2, r2);
}
#ifdef SWIFT_DEBUG_CHECKS
diff --git a/src/gravity_properties.c b/src/gravity_properties.c
index f3dc6c373f..5e0fa238a8 100644
--- a/src/gravity_properties.c
+++ b/src/gravity_properties.c
@@ -113,7 +113,8 @@ void gravity_props_print(const struct gravity_props *p) {
message("Self-gravity opening angle: theta=%.4f", p->theta_crit);
- message("Self-gravity softening: %s", kernel_gravity_softening_name);
+ message("Self-gravity softening functional form: %s",
+ kernel_gravity_softening_name);
message(
"Self-gravity comoving softening: epsilon=%.4f (Plummer equivalent: "
diff --git a/src/kernel_long_gravity.h b/src/kernel_long_gravity.h
index 2c14c0d888..1744f2cd04 100644
--- a/src/kernel_long_gravity.h
+++ b/src/kernel_long_gravity.h
@@ -31,7 +31,7 @@
#include <float.h>
#include <math.h>
-//#define GADGET2_LONG_RANGE_CORRECTION
+#define GADGET2_LONG_RANGE_CORRECTION
#ifdef GADGET2_LONG_RANGE_CORRECTION
#define kernel_long_gravity_truncation_name "Gadget-like (using erfc())"
diff --git a/src/multipole.h b/src/multipole.h
index f85a2fe63d..ac746728ee 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -2396,7 +2396,7 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
* @param r2 Square of the distance (periodically wrapped) between the
* multipoles.
*/
-__attribute__((always_inline)) INLINE static int gravity_M2L_accept(
+__attribute__((always_inline, const)) INLINE static int gravity_M2L_accept(
const double r_crit_a, const double r_crit_b, const double theta_crit2,
const double r2) {
@@ -2421,7 +2421,7 @@ __attribute__((always_inline)) INLINE static int gravity_M2L_accept(
* @param r2 Square of the distance (periodically wrapped) between the
* particle and the multipole.
*/
-__attribute__((always_inline)) INLINE static int gravity_M2P_accept(
+__attribute__((always_inline, const)) INLINE static int gravity_M2P_accept(
const float r_max2, const float theta_crit2, const float r2) {
// MATTHIEU: Make this mass-dependent ?
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 2c12d81223..f843394b0e 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -621,9 +621,11 @@ static INLINE void runner_dopair_grav_pm_truncated(
* @param ci The first #cell.
* @param cj The other #cell.
* @param symmetric Are we updating both cells (1) or just ci (0) ?
+ * @param Are we allowing the use of P2M interactions ?
*/
static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
- struct cell *cj, const int symmetric) {
+ struct cell *cj, const int symmetric,
+ const int allow_mpole) {
/* Recover some useful constants */
const struct engine *e = r->e;
@@ -691,12 +693,12 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
#endif
/* Fill the caches */
- gravity_cache_populate(e->max_active_bin, periodic, dim, ci_cache, ci->gparts,
- gcount_i, gcount_padded_i, shift_i, CoM_j, rmax2_j, ci,
- e->gravity_properties);
- gravity_cache_populate(e->max_active_bin, periodic, dim, cj_cache, cj->gparts,
- gcount_j, gcount_padded_j, shift_j, CoM_i, rmax2_i, cj,
- e->gravity_properties);
+ gravity_cache_populate(e->max_active_bin, allow_mpole, periodic, dim,
+ ci_cache, ci->gparts, gcount_i, gcount_padded_i,
+ shift_i, CoM_j, rmax2_j, ci, e->gravity_properties);
+ gravity_cache_populate(e->max_active_bin, allow_mpole, periodic, dim,
+ cj_cache, cj->gparts, gcount_j, gcount_padded_j,
+ shift_j, CoM_i, rmax2_i, cj, e->gravity_properties);
/* Can we use the Newtonian version or do we need the truncated one ? */
if (!periodic) {
@@ -712,8 +714,9 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
cj->gparts);
/* Then the M2P */
- runner_dopair_grav_pm_full(ci_cache, gcount_padded_i, CoM_j, multi_j,
- periodic, dim, e, ci->gparts, gcount_i, cj);
+ if (allow_mpole)
+ runner_dopair_grav_pm_full(ci_cache, gcount_padded_i, CoM_j, multi_j,
+ periodic, dim, e, ci->gparts, gcount_i, cj);
}
if (cj_active && symmetric) {
@@ -723,8 +726,9 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
ci->gparts);
/* Then the M2P */
- runner_dopair_grav_pm_full(cj_cache, gcount_padded_j, CoM_i, multi_i,
- periodic, dim, e, cj->gparts, gcount_j, ci);
+ if (allow_mpole)
+ runner_dopair_grav_pm_full(cj_cache, gcount_padded_j, CoM_i, multi_i,
+ periodic, dim, e, cj->gparts, gcount_j, ci);
}
} else { /* Periodic BC */
@@ -751,9 +755,10 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
ci->gparts, cj->gparts);
/* Then the M2P */
- runner_dopair_grav_pm_truncated(ci_cache, gcount_padded_i, CoM_j,
- multi_j, dim, r_s_inv, e, ci->gparts,
- gcount_i, cj);
+ if (allow_mpole)
+ runner_dopair_grav_pm_truncated(ci_cache, gcount_padded_i, CoM_j,
+ multi_j, dim, r_s_inv, e, ci->gparts,
+ gcount_i, cj);
}
if (cj_active && symmetric) {
@@ -763,9 +768,10 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
cj->gparts, ci->gparts);
/* Then the M2P */
- runner_dopair_grav_pm_truncated(cj_cache, gcount_padded_j, CoM_i,
- multi_i, dim, r_s_inv, e, cj->gparts,
- gcount_j, ci);
+ if (allow_mpole)
+ runner_dopair_grav_pm_truncated(cj_cache, gcount_padded_j, CoM_i,
+ multi_i, dim, r_s_inv, e, cj->gparts,
+ gcount_j, ci);
}
} else {
@@ -781,8 +787,10 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
ci->gparts, cj->gparts);
/* Then the M2P */
- runner_dopair_grav_pm_full(ci_cache, gcount_padded_i, CoM_j, multi_j,
- periodic, dim, e, ci->gparts, gcount_i, cj);
+ if (allow_mpole)
+ runner_dopair_grav_pm_full(ci_cache, gcount_padded_i, CoM_j, multi_j,
+ periodic, dim, e, ci->gparts, gcount_i,
+ cj);
}
if (cj_active && symmetric) {
@@ -792,8 +800,10 @@ static INLINE void runner_dopair_grav_pp(struct runner *r, struct cell *ci,
cj->gparts, ci->gparts);
/* Then the M2P */
- runner_dopair_grav_pm_full(cj_cache, gcount_padded_j, CoM_i, multi_i,
- periodic, dim, e, cj->gparts, gcount_j, ci);
+ if (allow_mpole)
+ runner_dopair_grav_pm_full(cj_cache, gcount_padded_j, CoM_i, multi_i,
+ periodic, dim, e, cj->gparts, gcount_j,
+ ci);
}
}
}
@@ -1337,7 +1347,7 @@ static INLINE void runner_dopair_recursive_grav(struct runner *r,
} else if (!ci->split && !cj->split) {
/* We have two leaves. Go P-P. */
- runner_dopair_grav_pp(r, ci, cj, 1);
+ runner_dopair_grav_pp(r, ci, cj, /*symmetric*/ 1, /*allow_mpoles*/ 0);
} else {
@@ -1546,13 +1556,13 @@ static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
/* Call the PM interaction fucntion on the active sub-cells of ci */
// runner_dopair_recursive_grav_pm(r, ci, cj);
- runner_dopair_grav_mm(r, ci, cj);
- // runner_dopair_grav_pp(r, ci, cj, 0);
+ // runner_dopair_grav_mm(r, ci, cj);
+ runner_dopair_grav_pp(r, ci, cj, 0, 1);
} /* We are in charge of this pair */
} /* Loop over top-level cells */
- message("cc=%d", cc);
+ // message("cc=%d", cc);
if (timer) TIMER_TOC(timer_dograv_long_range);
}
--
GitLab