diff --git a/src/const.h b/src/const.h
index 6432ef6a9e8107d94c93a32967e291d7e1e4d24d..8390da08980a47f45823fe7bcf3b784415305dd8 100644
--- a/src/const.h
+++ b/src/const.h
@@ -40,6 +40,7 @@
#define const_max_u_change 0.1f
/* Gravity stuff. */
+#define multipole_order 1
#define const_theta_max 0.57735f
#define const_G 6.672e-8f /* Gravitational constant. */
#define const_epsilon 0.0014f /* Gravity blending distance. */
diff --git a/src/engine.c b/src/engine.c
index d5fe48bf1f75d3afe17efef60391d9a9f099af48..e7c7e8aa363ca15203cc9c2f945b64de5d32c893 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1947,7 +1947,7 @@ void engine_init_particles(struct engine *e) {
/* Add the tasks corresponding to self-gravity to the masks */
if (e->policy & engine_policy_self_gravity) {
- /* Nothing here for now */
+ mask |= 1 << task_type_grav_up;
}
/* Add the tasks corresponding to external gravity to the masks */
@@ -2128,7 +2128,7 @@ void engine_step(struct engine *e) {
/* Add the tasks corresponding to self-gravity to the masks */
if (e->policy & engine_policy_self_gravity) {
- /* Nothing here for now */
+ mask |= 1 << task_type_grav_up;
}
/* Add the tasks corresponding to external gravity to the masks */
diff --git a/src/multipole.c b/src/multipole.c
index fc4701fead68e64e8742730325931282f80e8934..23d845c323efeab9f8bbf7206e925c6aef0d8d29 100644
--- a/src/multipole.c
+++ b/src/multipole.c
@@ -1,6 +1,7 @@
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
@@ -37,27 +38,28 @@
#include "multipole.h"
/**
- * @brief Merge two multipoles.
+ * @brief Add the second multipole to the first one.
*
- * @param ma The #multipole which will contain the merged result.
- * @param mb The other #multipole.
+ * @param m_sum The #multipole which will contain the sum.
+ * @param m_other The #multipole to add.
*/
-void multipole_merge(struct multipole *ma, struct multipole *mb) {
+void multipole_add(struct multipole *m_sum, const struct multipole *m_other) {
-#if multipole_order == 1
+#if multipole_order != 1
+#error "Multipoles of order >1 not yet implemented."
+#endif
/* Correct the position. */
- float mma = ma->coeffs[0], mmb = mb->coeffs[0];
- float w = 1.0f / (mma + mmb);
- for (int k = 0; k < 3; k++) ma->x[k] = (ma->x[k] * mma + mb->x[k] * mmb) * w;
+ const float M_tot = m_sum->mass + m_other->mass;
+ const float M_tot_inv = 1.f / M_tot;
+ for (int k = 0; k < 3; k++)
+ m_sum->CoM[k] =
+ (m_sum->CoM[k] * m_sum->mass + m_other->CoM[k] * m_other->mass) *
+ M_tot_inv;
/* Add the particle to the moments. */
- ma->coeffs[0] = mma + mmb;
-
-#else
-#error( "Multipoles of order %i not yet implemented." , multipole_order )
-#endif
+ m_sum->mass = M_tot;
}
/**
@@ -69,19 +71,21 @@ void multipole_merge(struct multipole *ma, struct multipole *mb) {
void multipole_addpart(struct multipole *m, struct gpart *p) {
-#if multipole_order == 1
+ /* #if multipole_order == 1 */
- /* Correct the position. */
- float mm = m->coeffs[0], mp = p->mass;
- float w = 1.0f / (mm + mp);
- for (int k = 0; k < 3; k++) m->x[k] = (m->x[k] * mm + p->x[k] * mp) * w;
+ /* /\* Correct the position. *\/ */
+ /* float mm = m->coeffs[0], mp = p->mass; */
+ /* float w = 1.0f / (mm + mp); */
+ /* for (int k = 0; k < 3; k++) m->x[k] = (m->x[k] * mm + p->x[k] * mp) * w;
+ */
- /* Add the particle to the moments. */
- m->coeffs[0] = mm + mp;
+ /* /\* Add the particle to the moments. *\/ */
+ /* m->coeffs[0] = mm + mp; */
-#else
-#error( "Multipoles of order %i not yet implemented." , multipole_order )
-#endif
+ /* #else */
+ /* #error( "Multipoles of order %i not yet implemented." , multipole_order )
+ */
+ /* #endif */
}
/**
@@ -94,76 +98,76 @@ void multipole_addpart(struct multipole *m, struct gpart *p) {
void multipole_addparts(struct multipole *m, struct gpart *p, int N) {
-#if multipole_order == 1
-
- /* Get the combined mass and positions. */
- double xp[3] = {0.0, 0.0, 0.0};
- float mp = 0.0f, w;
- for (int k = 0; k < N; k++) {
- w = p[k].mass;
- mp += w;
- xp[0] += p[k].x[0] * w;
- xp[1] += p[k].x[1] * w;
- xp[2] += p[k].x[2] * w;
- }
-
- /* Correct the position. */
- float mm = m->coeffs[0];
- w = 1.0f / (mm + mp);
- for (int k = 0; k < 3; k++) m->x[k] = (m->x[k] * mm + xp[k]) * w;
+ /* #if multipole_order == 1 */
+
+ /* /\* Get the combined mass and positions. *\/ */
+ /* double xp[3] = {0.0, 0.0, 0.0}; */
+ /* float mp = 0.0f, w; */
+ /* for (int k = 0; k < N; k++) { */
+ /* w = p[k].mass; */
+ /* mp += w; */
+ /* xp[0] += p[k].x[0] * w; */
+ /* xp[1] += p[k].x[1] * w; */
+ /* xp[2] += p[k].x[2] * w; */
+ /* } */
+
+ /* /\* Correct the position. *\/ */
+ /* float mm = m->coeffs[0]; */
+ /* w = 1.0f / (mm + mp); */
+ /* for (int k = 0; k < 3; k++) m->x[k] = (m->x[k] * mm + xp[k]) * w; */
+
+ /* /\* Add the particle to the moments. *\/ */
+ /* m->coeffs[0] = mm + mp; */
+
+ /* #else */
+ /* #error( "Multipoles of order %i not yet implemented." , multipole_order )
+ */
+ /* #endif */
+}
- /* Add the particle to the moments. */
- m->coeffs[0] = mm + mp;
+/**
+* @brief Reset the data of a #multipole.
+*
+* @param m The #multipole.
+*/
+void multipole_reset(struct multipole *m) {
-#else
-#error( "Multipoles of order %i not yet implemented." , multipole_order )
-#endif
+ /* Just bzero the struct. */
+ bzero(m, sizeof(struct multipole));
}
/**
- * @brief Init a multipole from a set of particles.
- *
- * @param m The #multipole.
- * @param parts The #gpart.
- * @param N The number of particles.
- */
-
-void multipole_init(struct multipole *m, struct gpart *parts, int N) {
+* @brief Init a multipole from a set of particles.
+*
+* @param m The #multipole.
+* @param parts The #gpart.
+* @param N The number of particles.
+*/
+void multipole_init(struct multipole *m, const struct gpart *gparts,
+ int gcount) {
+
+#if multipole_order != 1
+#error "Multipoles of order >1 not yet implemented."
+#endif
-#if multipole_order == 1
+ /* Zero everything */
+ multipole_reset(m);
- float mass = 0.0f, w;
+ float mass = 0.0f;
double x[3] = {0.0, 0.0, 0.0};
- int k;
/* Collect the particle data. */
- for (k = 0; k < N; k++) {
- w = parts[k].mass;
+ for (int k = 0; k < gcount; k++) {
+ const float w = gparts[k].mass;
mass += w;
- x[0] += parts[k].x[0] * w;
- x[1] += parts[k].x[1] * w;
- x[2] += parts[k].x[2] * w;
+ x[0] += gparts[k].x[0] * w;
+ x[1] += gparts[k].x[1] * w;
+ x[2] += gparts[k].x[2] * w;
}
/* Store the data on the multipole. */
- m->coeffs[0] = mass;
- m->x[0] = x[0] / mass;
- m->x[1] = x[1] / mass;
- m->x[2] = x[2] / mass;
-
-#else
-#error( "Multipoles of order %i not yet implemented." , multipole_order )
-#endif
-}
-
-/**
- * @brief Reset the data of a #multipole.
- *
- * @param m The #multipole.
- */
-
-void multipole_reset(struct multipole *m) {
-
- /* Just bzero the struct. */
- bzero(m, sizeof(struct multipole));
+ m->mass = mass;
+ m->CoM[0] = x[0] / mass;
+ m->CoM[1] = x[1] / mass;
+ m->CoM[2] = x[2] / mass;
}
diff --git a/src/multipole.h b/src/multipole.h
index 85ba44d3ce95d958b721d435ccd26b72e30a79c1..07e9ecbd48fec5844ef3df3555b48c5fa4b6c0e3 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -1,6 +1,7 @@
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
@@ -28,31 +29,33 @@
#include "kernel_gravity.h"
#include "part.h"
-/* Some constants. */
-#define multipole_order 1
-
/* Multipole struct. */
struct multipole {
/* Multipole location. */
- double x[3];
+ double CoM[3];
+
+ /* Multipole mass */
+ float mass;
- /* Acceleration on this multipole. */
- float a[3];
+ /* /\* Acceleration on this multipole. *\/ */
+ /* float a[3]; */
- /* Multipole coefficients. */
- float coeffs[multipole_order * multipole_order];
+ /* /\* Multipole coefficients. *\/ */
+ /* float coeffs[multipole_order * multipole_order]; */
};
/* Multipole function prototypes. */
-static void multipole_iact_mm(struct multipole *ma, struct multipole *mb,
- double *shift);
-void multipole_merge(struct multipole *ma, struct multipole *mb);
-void multipole_addpart(struct multipole *m, struct gpart *p);
-void multipole_addparts(struct multipole *m, struct gpart *p, int N);
-void multipole_init(struct multipole *m, struct gpart *parts, int N);
+void multipole_add(struct multipole *m_sum, const struct multipole *m_term);
+void multipole_init(struct multipole *m, const struct gpart *gparts,
+ int gcount);
void multipole_reset(struct multipole *m);
+/* static void multipole_iact_mm(struct multipole *ma, struct multipole *mb, */
+/* double *shift); */
+/* void multipole_addpart(struct multipole *m, struct gpart *p); */
+/* void multipole_addparts(struct multipole *m, struct gpart *p, int N); */
+
/**
* @brief Compute the pairwise interaction between two multipoles.
*
@@ -60,39 +63,38 @@ void multipole_reset(struct multipole *m);
* @param mb The second #multipole.
* @param shift The periodicity correction.
*/
-
__attribute__((always_inline)) INLINE static void multipole_iact_mm(
struct multipole *ma, struct multipole *mb, double *shift) {
-
- float dx[3], ir, r, r2 = 0.0f, acc;
- int k;
-
- /* Compute the multipole distance. */
- for (k = 0; k < 3; k++) {
- dx[k] = ma->x[k] - mb->x[k] - shift[k];
- r2 += dx[k] * dx[k];
- }
-
- /* Compute the normalized distance vector. */
- ir = 1.0f / sqrtf(r2);
- r = r2 * ir;
-
- /* Evaluate the gravity kernel. */
- kernel_grav_eval(r, &acc);
-
- /* Scale the acceleration. */
- acc *= const_G * ir * ir * ir;
-
-/* Compute the forces on both multipoles. */
-#if multipole_order == 1
- float mma = ma->coeffs[0], mmb = mb->coeffs[0];
- for (k = 0; k < 3; k++) {
- ma->a[k] -= dx[k] * acc * mmb;
- mb->a[k] += dx[k] * acc * mma;
- }
-#else
-#error( "Multipoles of order %i not yet implemented." , multipole_order )
-#endif
+ /* float dx[3], ir, r, r2 = 0.0f, acc; */
+ /* int k; */
+
+ /* /\* Compute the multipole distance. *\/ */
+ /* for (k = 0; k < 3; k++) { */
+ /* dx[k] = ma->x[k] - mb->x[k] - shift[k]; */
+ /* r2 += dx[k] * dx[k]; */
+ /* } */
+
+ /* /\* Compute the normalized distance vector. *\/ */
+ /* ir = 1.0f / sqrtf(r2); */
+ /* r = r2 * ir; */
+
+ /* /\* Evaluate the gravity kernel. *\/ */
+ /* kernel_grav_eval(r, &acc); */
+
+ /* /\* Scale the acceleration. *\/ */
+ /* acc *= const_G * ir * ir * ir; */
+
+ /* /\* Compute the forces on both multipoles. *\/ */
+ /* #if multipole_order == 1 */
+ /* float mma = ma->coeffs[0], mmb = mb->coeffs[0]; */
+ /* for (k = 0; k < 3; k++) { */
+ /* ma->a[k] -= dx[k] * acc * mmb; */
+ /* mb->a[k] += dx[k] * acc * mma; */
+ /* } */
+ /* #else */
+ /* #error( "Multipoles of order %i not yet implemented." , multipole_order )
+ */
+ /* #endif */
}
/**
@@ -102,35 +104,35 @@ __attribute__((always_inline)) INLINE static void multipole_iact_mm(
* @param p The #gpart.
* @param shift The periodicity correction.
*/
-
__attribute__((always_inline)) INLINE static void multipole_iact_mp(
struct multipole *m, struct gpart *p, double *shift) {
- float dx[3], ir, r, r2 = 0.0f, acc;
- int k;
+ /* float dx[3], ir, r, r2 = 0.0f, acc; */
+ /* int k; */
- /* Compute the multipole distance. */
- for (k = 0; k < 3; k++) {
- dx[k] = m->x[k] - p->x[k] - shift[k];
- r2 += dx[k] * dx[k];
- }
+ /* /\* Compute the multipole distance. *\/ */
+ /* for (k = 0; k < 3; k++) { */
+ /* dx[k] = m->x[k] - p->x[k] - shift[k]; */
+ /* r2 += dx[k] * dx[k]; */
+ /* } */
- /* Compute the normalized distance vector. */
- ir = 1.0f / sqrtf(r2);
- r = r2 * ir;
+ /* /\* Compute the normalized distance vector. *\/ */
+ /* ir = 1.0f / sqrtf(r2); */
+ /* r = r2 * ir; */
- /* Evaluate the gravity kernel. */
- kernel_grav_eval(r, &acc);
+ /* /\* Evaluate the gravity kernel. *\/ */
+ /* kernel_grav_eval(r, &acc); */
- /* Scale the acceleration. */
- acc *= const_G * ir * ir * ir * m->coeffs[0];
+ /* /\* Scale the acceleration. *\/ */
+ /* acc *= const_G * ir * ir * ir * m->coeffs[0]; */
-/* Compute the forces on both multipoles. */
-#if multipole_order == 1
- for (k = 0; k < 3; k++) p->a_grav[k] += dx[k] * acc;
-#else
-#error( "Multipoles of order %i not yet implemented." , multipole_order )
-#endif
+ /* /\* Compute the forces on both multipoles. *\/ */
+ /* #if multipole_order == 1 */
+ /* for (k = 0; k < 3; k++) p->a_grav[k] += dx[k] * acc; */
+ /* #else */
+ /* #error( "Multipoles of order %i not yet implemented." , multipole_order )
+ */
+ /* #endif */
}
#endif /* SWIFT_MULTIPOLE_H */
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 19e77879501c1e34a54e0d436e9a140a9b273d1a..31aac56827eb8f52af13a883ce8f5e002c6ded9e 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -1,6 +1,7 @@
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (c) 2013 Pedro Gonnet (pedro.gonnet@durham.ac.uk)
+ * 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
@@ -39,7 +40,8 @@
/* int pid, pjd, k, sid; */
/* double rshift, shift[3] = {0.0, 0.0, 0.0}, nshift[3]; */
/* struct entry *restrict sort_i, *restrict sort_j; */
-/* struct gpart *restrict pi, *restrict pj, *restrict parts_i, *restrict parts_j; */
+/* struct gpart *restrict pi, *restrict pj, *restrict parts_i, *restrict
+ * parts_j; */
/* double pix[3]; */
/* float dx[3], r2, h_max, di, dj; */
/* int count_i, count_j, cnj, cnj_new; */
@@ -77,7 +79,8 @@
/* /\* Get some other useful values. *\/ */
/* h_max = */
-/* sqrtf(ci->h[0] * ci->h[0] + ci->h[1] * ci->h[1] + ci->h[2] * ci->h[2]) * */
+/* sqrtf(ci->h[0] * ci->h[0] + ci->h[1] * ci->h[1] + ci->h[2] * ci->h[2])
+ * * */
/* const_theta_max; */
/* count_i = ci->gcount; */
/* count_j = cj->gcount; */
@@ -114,9 +117,12 @@
/* #ifndef VECTORIZE */
-/* // if ( pi->part->id == 3473472412525 || pj->part->id == 3473472412525 ) */
-/* // message( "interacting particles pi=%lli and pj=%lli with r=%.3e in */
-/* // cells %lli/%lli." , pi->part->id , pj->part->id , sqrtf(r2) , ((long */
+/* // if ( pi->part->id == 3473472412525 || pj->part->id == 3473472412525
+ * ) */
+/* // message( "interacting particles pi=%lli and pj=%lli with r=%.3e
+ * in */
+/* // cells %lli/%lli." , pi->part->id , pj->part->id , sqrtf(r2) , ((long
+ */
/* // long int)ci) / sizeof(struct cell) , ((long long int)cj) / */
/* // sizeof(struct cell) ); */
@@ -206,11 +212,7 @@
void runner_dograv_up(struct runner *r, struct cell *c) {
- /* Re-set this cell's multipole. */
- multipole_reset(&c->multipole);
-
- /* Split? */
- if (c->split) {
+ if (c->split) {/* Regular node */
/* Recurse. */
for (int k = 0; k < 8; k++)
@@ -218,17 +220,16 @@ void runner_dograv_up(struct runner *r, struct cell *c) {
/* Collect the multipoles from the progeny. */
multipole_reset(&c->multipole);
- for (int k = 0; k < 8; k++)
+ for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL)
- multipole_merge(&c->multipole, &c->progeny[k]->multipole);
+ multipole_add(&c->multipole, &c->progeny[k]->multipole);
+ }
- }
+ } else {/* Leaf node. */
- /* No, leaf node. */
- else
-
- /* Just collect the multipole. */
+ /* Just construct the multipole from the gparts. */
multipole_init(&c->multipole, c->gparts, c->gcount);
+ }
}
/* /\** */
@@ -303,7 +304,8 @@ void runner_dograv_up(struct runner *r, struct cell *c) {
/* } */
/* theta = */
/* sqrt((dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]) / */
-/* (ci->h[0] * ci->h[0] + ci->h[1] * ci->h[1] + ci->h[2] * ci->h[2])); */
+/* (ci->h[0] * ci->h[0] + ci->h[1] * ci->h[1] + ci->h[2] *
+ * ci->h[2])); */
/* /\* Do an MM or an MP/PM? *\/ */
/* if (theta > const_theta_max * 4) { */
@@ -335,7 +337,8 @@ void runner_dograv_up(struct runner *r, struct cell *c) {
/* struct engine *e = r->e; */
/* int pid, pjd, k, count_i = ci->gcount, count_j = cj->gcount; */
/* double shift[3] = {0.0, 0.0, 0.0}; */
-/* struct gpart *restrict parts_i = ci->gparts, *restrict parts_j = cj->gparts; */
+/* struct gpart *restrict parts_i = ci->gparts, *restrict parts_j =
+ * cj->gparts; */
/* struct gpart *restrict pi, *restrict pj; */
/* double pix[3]; */
/* float dx[3], r2; */
@@ -383,9 +386,12 @@ void runner_dograv_up(struct runner *r, struct cell *c) {
/* /\* Compute the interaction. *\/ */
/* #ifndef VECTORIZE */
-/* // if ( pi->part->id == 3473472412525 || pj->part->id == 3473472412525 ) */
-/* // message( "interacting particles pi=%lli and pj=%lli with r=%.3e in */
-/* // cells %lli/%lli." , pi->part->id , pj->part->id , sqrtf(r2) , ((long */
+/* // if ( pi->part->id == 3473472412525 || pj->part->id == 3473472412525
+ * ) */
+/* // message( "interacting particles pi=%lli and pj=%lli with r=%.3e
+ * in */
+/* // cells %lli/%lli." , pi->part->id , pj->part->id , sqrtf(r2) , ((long
+ */
/* // long int)ci) / sizeof(struct cell) , ((long long int)cj) / */
/* // sizeof(struct cell) ); */
@@ -473,8 +479,10 @@ void runner_dograv_up(struct runner *r, struct cell *c) {
/* /\* Compute the interaction. *\/ */
/* #ifndef VECTORIZE */
-/* // if ( pi->part->id == 3473472412525 || pj->part->id == 3473472412525 ) */
-/* // message( "interacting particles pi=%lli and pj=%lli with r=%.3e." , */
+/* // if ( pi->part->id == 3473472412525 || pj->part->id == 3473472412525
+ * ) */
+/* // message( "interacting particles pi=%lli and pj=%lli with
+ * r=%.3e." , */
/* // pi->part->id , pj->part->id , sqrtf(r2) ); */
/* runner_iact_grav(r2, dx, pi, pj); */
@@ -563,7 +571,8 @@ void runner_dograv_up(struct runner *r, struct cell *c) {
/* if (dx[k] > 0.0f) dx[k] -= ci->h[k]; */
/* } */
/* theta = (dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]) / */
-/* (ci->h[0] * ci->h[0] + ci->h[1] * ci->h[1] + ci->h[2] * ci->h[2]); */
+/* (ci->h[0] * ci->h[0] + ci->h[1] * ci->h[1] + ci->h[2] *
+ * ci->h[2]); */
/* /\* Split the interaction? *\/ */
/* if (theta < const_theta_max * const_theta_max) { */