From a49b8bd5e11e65ea85b5cc03101167fca1b97e7a Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Sun, 18 Jun 2017 12:36:58 +0200
Subject: [PATCH] Made the M2L kernel use periodic boundary conditions.
---
src/multipole.h | 21 +++++++++++++--------
src/runner_doiact_grav.h | 11 ++++++++---
2 files changed, 21 insertions(+), 11 deletions(-)
diff --git a/src/multipole.h b/src/multipole.h
index 8e8ecae54e..a0728e47b8 100644
--- a/src/multipole.h
+++ b/src/multipole.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 ?
+ * @param dim The size of the simulation box.
*/
INLINE static void gravity_M2L(struct grav_tensor *l_b,
const struct multipole *m_a,
const double pos_b[3], const double pos_a[3],
- const struct gravity_props *props,
- int periodic) {
+ const struct gravity_props *props, int periodic,
+ const double dim[3]) {
/* 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];
- const double dy =
- periodic ? box_wrap(pos_b[1] - pos_a[1], 0., 1.) : pos_b[1] - pos_a[1];
- const double dz =
- periodic ? box_wrap(pos_b[2] - pos_a[2], 0., 1.) : pos_b[2] - pos_a[2];
+ double dx = pos_b[0] - pos_a[0];
+ double dy = pos_b[1] - pos_a[1];
+ double dz = pos_b[2] - pos_a[2];
+
+ /* Apply BC */
+ if (periodic) {
+ dx = nearest(dx, dim[0]);
+ dy = nearest(dy, dim[1]);
+ dz = nearest(dz, dim[2]);
+ }
/* Compute distance */
const double r2 = dx * dx + dy * dy + dz * dz;
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 0a5915cd10..daf0914dac 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -110,10 +110,12 @@ void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
struct cell *restrict cj) {
+ /* Some constants */
const struct engine *e = r->e;
+ const struct space *s = e->s;
+ const int periodic = s->periodic;
+ const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const struct gravity_props *props = e->gravity_properties;
- const int periodic = e->s->periodic;
- const struct multipole *multi_j = &cj->multipole->m_pole;
// const float a_smooth = e->gravity_properties->a_smooth;
// const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
@@ -122,6 +124,9 @@ void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
/* Anything to do here? */
if (!cell_is_active(ci, e)) return;
+ /* Short-cut to the multipole */
+ const struct multipole *multi_j = &cj->multipole->m_pole;
+
#ifdef SWIFT_DEBUG_CHECKS
if (ci == cj) error("Interacting a cell with itself using M2L");
@@ -136,7 +141,7 @@ void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
/* Let's interact at this level */
gravity_M2L(&ci->multipole->pot, multi_j, ci->multipole->CoM,
- cj->multipole->CoM, props, periodic * 0);
+ cj->multipole->CoM, props, periodic, dim);
TIMER_TOC(timer_dopair_grav_mm);
}
--
GitLab