From 2e7b35b16f206beeaa10bbac3c0c93f27b2f4a3c Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Fri, 3 Mar 2017 08:56:31 +0000
Subject: [PATCH] Don't initialise multipoles if they don't exist.
---
src/multipole.h | 49 ++++++++++++++++++++--------------------
src/runner.c | 3 ++-
src/runner_doiact_grav.h | 38 ++++++++++++++++++++++++++++++-
src/timers.h | 1 +
4 files changed, 64 insertions(+), 27 deletions(-)
diff --git a/src/multipole.h b/src/multipole.h
index 37b5defc34..6065a85633 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -289,31 +289,30 @@ INLINE static void multipole_M2M(struct multipole *m_a,
* @param pos_b The position of the multipole.
* @param periodic Is the calculation periodic ?
*/
-/* INLINE static void multipole_M2L(struct field_tensors *l_a, */
-/* const struct multipole m_b, */
-/* const double pos_a[3], const double
- * pos_b[3], */
-/* int periodic) { */
-
-/* /\* double dx, dy, dz; *\/ */
-/* /\* if (periodic) { *\/ */
-/* /\* dx = box_wrap(pos_a[0] - pos_b[0], 0., 1.); *\/ */
-/* /\* dy = box_wrap(pos_a[1] - pos_b[1], 0., 1.); *\/ */
-/* /\* dz = box_wrap(pos_a[2] - pos_b[2], 0., 1.); *\/ */
-/* /\* } else { *\/ */
-/* /\* dx = pos_a[0] - pos_b[0]; *\/ */
-/* /\* dy = pos_a[1] - pos_b[1]; *\/ */
-/* /\* dz = pos_a[2] - pos_b[2]; *\/ */
-/* /\* } *\/ */
-/* /\* const double r2 = dx * dx + dy * dy + dz * dz; *\/ */
-
-/* /\* const double r_inv = 1. / sqrt(r2); *\/ */
-
-/* /\* /\\* 1st order multipole term *\\/ *\/ */
-/* /\* l_a->x.F_000 = D_100(dx, dy, dz, r_inv); *\/ */
-/* /\* l_a->y.F_000 = D_010(dx, dy, dz, r_inv); *\/ */
-/* /\* l_a->z.F_000 = D_001(dx, dy, dz, r_inv); *\/ */
-/* } */
+INLINE static void multipole_M2L(struct gravity_tensors *l_a,
+ const struct multipole *m_b,
+ const double pos_a[3], const double pos_b[3],
+ int periodic) {
+
+ double dx, dy, dz;
+ if (periodic) {
+ dx = box_wrap(pos_a[0] - pos_b[0], 0., 1.);
+ dy = box_wrap(pos_a[1] - pos_b[1], 0., 1.);
+ dz = box_wrap(pos_a[2] - pos_b[2], 0., 1.);
+ } else {
+ dx = pos_a[0] - pos_b[0];
+ dy = pos_a[1] - pos_b[1];
+ dz = pos_a[2] - pos_b[2];
+ }
+ const double r2 = dx * dx + dy * dy + dz * dz;
+
+ const double r_inv = 1. / sqrt(r2);
+
+ /* 1st order multipole term */
+ l_a->a_x.F_000 = D_100(dx, dy, dz, r_inv) * m_b->mass;
+ l_a->a_y.F_000 = D_010(dx, dy, dz, r_inv) * m_b->mass;
+ l_a->a_z.F_000 = D_001(dx, dy, dz, r_inv) * m_b->mass;
+}
#if 0
diff --git a/src/runner.c b/src/runner.c
index 9241e630c3..2e581e029f 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -496,7 +496,8 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
if (!cell_is_active(c, e)) return;
/* Reset the gravity acceleration tensors */
- if (c->multipole) multipole_init(c->multipole);
+ if (e->policy & engine_policy_self_gravity)
+ multipole_init(c->multipole);
/* Recurse? */
if (c->split) {
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 6704631716..c9e95dc47e 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -53,6 +53,42 @@ void runner_do_grav_up(struct runner *r, struct cell *c) {
/* } */
}
+/**
+ * @brief Computes the interaction of the field tensor in a cell with the
+ * multipole of another cell.
+ *
+ * @param r The #runner.
+ * @param ci The #cell with field tensor to interact.
+ * @param cj The #cell with the multipole.
+ */
+__attribute__((always_inline)) INLINE static void runner_dopair_grav_mm(
+ const struct runner *r, const struct cell *restrict ci,
+ const struct cell *restrict cj) {
+
+ const struct engine *e = r->e;
+ 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]);
+
+ TIMER_TIC;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (multi_j->mass == 0.0)
+ error("Multipole does not seem to have been set.");
+#endif
+
+ /* Anything to do here? */
+ if (!cell_is_active(ci, e)) return;
+
+ multipole_M2L(ci->multipole, multi_j, ci->multipole->CoM, cj->multipole->CoM,
+ periodic);
+
+
+ TIMER_TOC(timer_dopair_grav_mm);
+}
+
+
/**
* @brief Computes the interaction of all the particles in a cell with the
* multipole of another cell.
@@ -541,7 +577,7 @@ static void runner_do_grav_long_range(struct runner *r, struct cell *ci,
// if (r2 > max_d2) continue;
- if (!cell_are_neighbours(ci, cj)) runner_dopair_grav_pm(r, ci, cj);
+ if (!cell_are_neighbours(ci, cj)) runner_dopair_grav_mm(r, ci, cj);
}
}
diff --git a/src/timers.h b/src/timers.h
index 50f630e7fc..4cb4d7e0ba 100644
--- a/src/timers.h
+++ b/src/timers.h
@@ -46,6 +46,7 @@ enum {
timer_dopair_gradient,
timer_dopair_force,
timer_dopair_grav_pm,
+ timer_dopair_grav_mm,
timer_dopair_grav_pp,
timer_dograv_external,
timer_dosource,
--
GitLab