diff --git a/src/Makefile.am b/src/Makefile.am
index 3e759d8fd15cfad6b38fe4ed469fb3b745aff786..2fc7ac5134c793f9476136759676a3ffaaf9816d 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -45,7 +45,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
physical_constants.h physical_constants_cgs.h potential.h version.h \
hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.h \
sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h profiler.h \
- dump.h logger.h active.h timeline.h xmf.h gravity_properties.h
+ dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
diff --git a/src/cell.c b/src/cell.c
index 126e6c7e3d1fd9f73a49055c565b0d4fde937321..941ee7a7623b0006f0291878da58f7c8fa966565 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1073,50 +1073,27 @@ int cell_are_neighbours(const struct cell *restrict ci,
void cell_check_multipole(struct cell *c, void *data) {
struct multipole ma;
+ const double tolerance = 1e-5; /* Relative */
+
+ /* First recurse */
+ if (c->split)
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) cell_check_multipole(c->progeny[k], NULL);
if (c->gcount > 0) {
/* Brute-force calculation */
- multipole_init(&ma, c->gparts, c->gcount);
-
- /* Compare with recursive one */
- struct multipole mb = c->multipole;
-
- if (fabsf(ma.mass - mb.mass) / fabsf(ma.mass + mb.mass) > 1e-5)
- error("Multipole masses are different (%12.15e vs. %12.15e)", ma.mass,
- mb.mass);
-
- for (int k = 0; k < 3; ++k)
- if (fabs(ma.CoM[k] - mb.CoM[k]) / fabs(ma.CoM[k] + mb.CoM[k]) > 1e-5)
- error("Multipole CoM are different (%12.15e vs. %12.15e", ma.CoM[k],
- mb.CoM[k]);
-
-#if const_gravity_multipole_order >= 2
- if (fabsf(ma.I_xx - mb.I_xx) / fabsf(ma.I_xx + mb.I_xx) > 1e-5 &&
- ma.I_xx > 1e-9)
- error("Multipole I_xx are different (%12.15e vs. %12.15e)", ma.I_xx,
- mb.I_xx);
- if (fabsf(ma.I_yy - mb.I_yy) / fabsf(ma.I_yy + mb.I_yy) > 1e-5 &&
- ma.I_yy > 1e-9)
- error("Multipole I_yy are different (%12.15e vs. %12.15e)", ma.I_yy,
- mb.I_yy);
- if (fabsf(ma.I_zz - mb.I_zz) / fabsf(ma.I_zz + mb.I_zz) > 1e-5 &&
- ma.I_zz > 1e-9)
- error("Multipole I_zz are different (%12.15e vs. %12.15e)", ma.I_zz,
- mb.I_zz);
- if (fabsf(ma.I_xy - mb.I_xy) / fabsf(ma.I_xy + mb.I_xy) > 1e-5 &&
- ma.I_xy > 1e-9)
- error("Multipole I_xy are different (%12.15e vs. %12.15e)", ma.I_xy,
- mb.I_xy);
- if (fabsf(ma.I_xz - mb.I_xz) / fabsf(ma.I_xz + mb.I_xz) > 1e-5 &&
- ma.I_xz > 1e-9)
- error("Multipole I_xz are different (%12.15e vs. %12.15e)", ma.I_xz,
- mb.I_xz);
- if (fabsf(ma.I_yz - mb.I_yz) / fabsf(ma.I_yz + mb.I_yz) > 1e-5 &&
- ma.I_yz > 1e-9)
- error("Multipole I_yz are different (%12.15e vs. %12.15e)", ma.I_yz,
- mb.I_yz);
-#endif
+ multipole_P2M(&ma, c->gparts, c->gcount);
+
+ /* Now compare the multipole expansion */
+ if (!multipole_equal(&ma, c->multipole, tolerance)) {
+ message("Multipoles are not equal at depth=%d!", c->depth);
+ message("Correct answer:");
+ multipole_print(&ma);
+ message("Recursive multipole:");
+ multipole_print(c->multipole);
+ error("Aborting");
+ }
}
}
diff --git a/src/cell.h b/src/cell.h
index 86748161c79b39e038adc455e58376229f49a841..0d5a46006bd34cce4e65e9c2f9761292c31f40be 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -94,9 +94,6 @@ struct pcell {
*/
struct cell {
- /*! This cell's multipole. */
- struct multipole multipole;
-
/*! The cell location on the grid. */
double loc[3];
@@ -106,6 +103,9 @@ struct cell {
/*! Max smoothing length in this cell. */
double h_max;
+ /*! This cell's multipole. */
+ struct multipole *multipole;
+
/*! Linking pointer for "memory management". */
struct cell *next;
diff --git a/src/engine.c b/src/engine.c
index 284b9cfe83bee77948b3f0c741bedc4ef3697136..531b84c471514341736b0e846efec4e6c734e43e 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -3448,6 +3448,7 @@ void engine_unpin() {
* @param internal_units The system of units used internally.
* @param physical_constants The #phys_const used for this run.
* @param hydro The #hydro_props used for this run.
+ * @param gravity The #gravity_props used for this run.
* @param potential The properties of the external potential.
* @param cooling_func The properties of the cooling function.
* @param sourceterms The properties of the source terms function.
diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
new file mode 100644
index 0000000000000000000000000000000000000000..4730d9df5dc573de74fde422e1b7dafc0ee0994a
--- /dev/null
+++ b/src/gravity_derivatives.h
@@ -0,0 +1,72 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 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
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_GRAVITY_DERIVATIVE_H
+#define SWIFT_GRAVITY_DERIVATIVE_H
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local headers. */
+#include "inline.h"
+
+/**
+ * @brief \f$ \phi(r_x, r_y, r_z) \f$.
+ */
+__attribute__((always_inline)) INLINE static double D_000(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+
+ return r_inv;
+}
+
+/**
+ * @brief \f$ \frac{\partial\phi(r_x, r_y, r_z)}{\partial r_x} \f$.
+ */
+__attribute__((always_inline)) INLINE static double D_100(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+
+ return -r_x * r_inv * r_inv * r_inv;
+}
+
+/**
+ * @brief \f$ \frac{\partial\phi(r_x, r_y, r_z)}{\partial r_x} \f$.
+ */
+__attribute__((always_inline)) INLINE static double D_010(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+
+ return -r_y * r_inv * r_inv * r_inv;
+}
+
+/**
+ * @brief \f$ \frac{\partial\phi(r_x, r_y, r_z)}{\partial r_x} \f$.
+ */
+__attribute__((always_inline)) INLINE static double D_001(double r_x,
+ double r_y,
+ double r_z,
+ double r_inv) {
+
+ return -r_z * r_inv * r_inv * r_inv;
+}
+
+#endif /* SWIFT_GRAVITY_DERIVATIVE_H */
diff --git a/src/multipole.c b/src/multipole.c
index f0d0c7084fd9bec366f2185cb24887da40591b17..bd5c6d6546fa0546108dcd53d7fe4060293c37a7 100644
--- a/src/multipole.c
+++ b/src/multipole.c
@@ -20,205 +20,3 @@
/* Config parameters. */
#include "../config.h"
-
-/* Some standard headers. */
-#include <strings.h>
-
-/* This object's header. */
-#include "multipole.h"
-
-/**
-* @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));
-}
-
-/**
-* @brief Init a multipole from a set of particles.
-*
-* @param m The #multipole.
-* @param gparts The #gpart.
-* @param gcount The number of particles.
-*/
-void multipole_init(struct multipole *m, const struct gpart *gparts,
- int gcount) {
-
-#if const_gravity_multipole_order > 2
-#error "Multipoles of order >2 not yet implemented."
-#endif
-
- /* Zero everything */
- multipole_reset(m);
-
- /* Temporary variables */
- double mass = 0.0;
- double com[3] = {0.0, 0.0, 0.0};
-
-#if const_gravity_multipole_order >= 2
- double I_xx = 0.0, I_yy = 0.0, I_zz = 0.0;
- double I_xy = 0.0, I_xz = 0.0, I_yz = 0.0;
-#endif
-
- /* Collect the particle data. */
- for (int k = 0; k < gcount; k++) {
- const float w = gparts[k].mass;
-
- mass += w;
- com[0] += gparts[k].x[0] * w;
- com[1] += gparts[k].x[1] * w;
- com[2] += gparts[k].x[2] * w;
-
-#if const_gravity_multipole_order >= 2
- I_xx += gparts[k].x[0] * gparts[k].x[0] * w;
- I_yy += gparts[k].x[1] * gparts[k].x[1] * w;
- I_zz += gparts[k].x[2] * gparts[k].x[2] * w;
- I_xy += gparts[k].x[0] * gparts[k].x[1] * w;
- I_xz += gparts[k].x[0] * gparts[k].x[2] * w;
- I_yz += gparts[k].x[1] * gparts[k].x[2] * w;
-#endif
- }
-
- const double imass = 1.0 / mass;
-
- /* Store the data on the multipole. */
- m->mass = mass;
- m->CoM[0] = com[0] * imass;
- m->CoM[1] = com[1] * imass;
- m->CoM[2] = com[2] * imass;
-
-#if const_gravity_multipole_order >= 2
- m->I_xx = I_xx - imass * com[0] * com[0];
- m->I_yy = I_yy - imass * com[1] * com[1];
- m->I_zz = I_zz - imass * com[2] * com[2];
- m->I_xy = I_xy - imass * com[0] * com[1];
- m->I_xz = I_xz - imass * com[0] * com[2];
- m->I_yz = I_yz - imass * com[1] * com[2];
-#endif
-}
-
-/**
- * @brief Add the second multipole to the first one.
- *
- * @param ma The #multipole which will contain the sum.
- * @param mb The #multipole to add.
- */
-
-void multipole_add(struct multipole *ma, const struct multipole *mb) {
-
-#if const_gravity_multipole_order > 2
-#error "Multipoles of order >2 not yet implemented."
-#endif
-
- /* Correct the position. */
- const double ma_mass = ma->mass;
- const double mb_mass = mb->mass;
- const double M_tot = ma_mass + mb_mass;
- const double M_tot_inv = 1.0 / M_tot;
-
- const double ma_CoM[3] = {ma->CoM[0], ma->CoM[1], ma->CoM[2]};
- const double mb_CoM[3] = {mb->CoM[0], mb->CoM[1], mb->CoM[2]};
-
-#if const_gravity_multipole_order >= 2
- const double ma_I_xx = (double)ma->I_xx + ma_mass * ma_CoM[0] * ma_CoM[0];
- const double ma_I_yy = (double)ma->I_yy + ma_mass * ma_CoM[1] * ma_CoM[1];
- const double ma_I_zz = (double)ma->I_zz + ma_mass * ma_CoM[2] * ma_CoM[2];
- const double ma_I_xy = (double)ma->I_xy + ma_mass * ma_CoM[0] * ma_CoM[1];
- const double ma_I_xz = (double)ma->I_xz + ma_mass * ma_CoM[0] * ma_CoM[2];
- const double ma_I_yz = (double)ma->I_yz + ma_mass * ma_CoM[1] * ma_CoM[2];
-
- const double mb_I_xx = (double)mb->I_xx + mb_mass * mb_CoM[0] * mb_CoM[0];
- const double mb_I_yy = (double)mb->I_yy + mb_mass * mb_CoM[1] * mb_CoM[1];
- const double mb_I_zz = (double)mb->I_zz + mb_mass * mb_CoM[2] * mb_CoM[2];
- const double mb_I_xy = (double)mb->I_xy + mb_mass * mb_CoM[0] * mb_CoM[1];
- const double mb_I_xz = (double)mb->I_xz + mb_mass * mb_CoM[0] * mb_CoM[2];
- const double mb_I_yz = (double)mb->I_yz + mb_mass * mb_CoM[1] * mb_CoM[2];
-#endif
-
- /* New mass */
- ma->mass = M_tot;
-
- /* New CoM */
- ma->CoM[0] = (ma_CoM[0] * ma_mass + mb_CoM[0] * mb_mass) * M_tot_inv;
- ma->CoM[1] = (ma_CoM[1] * ma_mass + mb_CoM[1] * mb_mass) * M_tot_inv;
- ma->CoM[2] = (ma_CoM[2] * ma_mass + mb_CoM[2] * mb_mass) * M_tot_inv;
-
-/* New quadrupole */
-#if const_gravity_multipole_order >= 2
- ma->I_xx = (ma_I_xx + mb_I_xx) - M_tot * ma->CoM[0] * ma->CoM[0];
- ma->I_yy = (ma_I_yy + mb_I_yy) - M_tot * ma->CoM[1] * ma->CoM[1];
- ma->I_zz = (ma_I_zz + mb_I_zz) - M_tot * ma->CoM[2] * ma->CoM[2];
- ma->I_xy = (ma_I_xy + mb_I_xy) - M_tot * ma->CoM[0] * ma->CoM[1];
- ma->I_xz = (ma_I_xz + mb_I_xz) - M_tot * ma->CoM[0] * ma->CoM[2];
- ma->I_yz = (ma_I_yz + mb_I_yz) - M_tot * ma->CoM[1] * ma->CoM[2];
-#endif
-}
-
-/**
- * @brief Add a particle to the given multipole.
- *
- * @param m The #multipole.
- * @param p The #gpart.
- */
-
-void multipole_addpart(struct multipole *m, struct gpart *p) {
-
- /* #if const_gravity_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;
- */
-
- /* /\* Add the particle to the moments. *\/ */
- /* m->coeffs[0] = mm + mp; */
-
- /* #else */
- /* #error( "Multipoles of order %i not yet implemented." ,
- * const_gravity_multipole_order )
- */
- /* #endif */
-}
-
-/**
- * @brief Add a group of particles to the given multipole.
- *
- * @param m The #multipole.
- * @param p The #gpart array.
- * @param N Number of parts to add.
- */
-
-void multipole_addparts(struct multipole *m, struct gpart *p, int N) {
-
- /* #if const_gravity_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." ,
- * const_gravity_multipole_order )
- */
- /* #endif */
-}
diff --git a/src/multipole.h b/src/multipole.h
index dc1f914dd73969bc63f3beffca7f34d0f889edb6..69195f8ec72cb18c1f661f0e94802ad4a2014739 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -22,28 +22,286 @@
/* Some standard headers. */
#include <math.h>
+#include <string.h>
/* Includes. */
+#include "align.h"
#include "const.h"
+#include "error.h"
+#include "gravity_derivatives.h"
#include "inline.h"
#include "kernel_gravity.h"
+#include "minmax.h"
#include "part.h"
-/* Multipole struct. */
+#define multipole_align 128
+
+/**
+ * @brief The multipole expansion of a mass distribution.
+ */
struct multipole {
- /* Multipole location. */
- double CoM[3];
+ union {
+
+ /*! Linking pointer for "memory management". */
+ struct multipole *next;
+
+ /*! The actual content */
+ struct {
+
+ /*! Multipole mass */
+ float mass;
+
+ /*! Centre of mass of the matter dsitribution */
+ double CoM[3];
+
+ /*! Bulk velocity */
+ float vel[3];
+ };
+ };
+} SWIFT_STRUCT_ALIGN;
+
+struct acc_tensor {
+
+ double F_000;
+};
+
+struct pot_tensor {
+
+ double F_000;
+};
+
+struct field_tensors {
+
+ union {
+
+ /*! Linking pointer for "memory management". */
+ struct field_tensors *next;
+
+ /*! The actual content */
+ struct {
+
+ /*! Field tensor for acceleration along x */
+ struct acc_tensor a_x;
+
+ /*! Field tensor for acceleration along y */
+ struct acc_tensor a_y;
+
+ /*! Field tensor for acceleration along z */
+ struct acc_tensor a_z;
+
+ /*! Field tensor for the potential */
+ struct pot_tensor pot;
+ };
+ };
+
+} SWIFT_STRUCT_ALIGN;
+
+/**
+ * @brief Reset the data of a #multipole.
+ *
+ * @param m The #multipole.
+ */
+INLINE static void multipole_init(struct multipole *m) {
+
+ /* Just bzero the struct. */
+ bzero(m, sizeof(struct multipole));
+}
+
+/**
+ * @brief Prints the content of a #multipole to stdout.
+ *
+ * Note: Uses directly printf(), not a call to message().
+ *
+ * @param m The #multipole to print.
+ */
+INLINE static void multipole_print(const struct multipole *m) {
+
+ printf("CoM= [%12.5e %12.5e %12.5e\n", m->CoM[0], m->CoM[1], m->CoM[2]);
+ printf("Mass= %12.5e\n", m->mass);
+ printf("Vel= [%12.5e %12.5e %12.5e\n", m->vel[0], m->vel[1], m->vel[2]);
+}
+
+/**
+ * @brief Adds a #multipole to another one (i.e. does ma += mb).
+ *
+ * @param ma The multipole to add to.
+ * @param mb The multipole to add.
+ */
+INLINE static void multipole_add(struct multipole *ma,
+ const struct multipole *mb) {
+
+ const float mass = ma->mass + mb->mass;
+ const float imass = 1.f / mass;
+
+ /* Add the bulk velocities */
+ ma->vel[0] = (ma->vel[0] * ma->mass + mb->vel[0] * mb->mass) * imass;
+ ma->vel[1] = (ma->vel[1] * ma->mass + mb->vel[1] * mb->mass) * imass;
+ ma->vel[2] = (ma->vel[2] * ma->mass + mb->vel[2] * mb->mass) * imass;
+
+ /* Add the masses */
+ ma->mass = mass;
+}
+
+/**
+ * @brief Verifies whether two #multipole's are equal or not.
+ *
+ * @param ma The first #multipole.
+ * @param mb The second #multipole.
+ * @param tolerance The maximal allowed difference for the fields.
+ * @return 1 if the multipoles are equal 0 otherwise.
+ */
+INLINE static int multipole_equal(const struct multipole *ma,
+ const struct multipole *mb,
+ double tolerance) {
+
+ /* Check CoM */
+ if (fabs(ma->CoM[0] - mb->CoM[0]) / fabs(ma->CoM[0] + mb->CoM[0]) > tolerance)
+ return 0;
+ if (fabs(ma->CoM[1] - mb->CoM[1]) / fabs(ma->CoM[1] + mb->CoM[1]) > tolerance)
+ return 0;
+ if (fabs(ma->CoM[2] - mb->CoM[2]) / fabs(ma->CoM[2] + mb->CoM[2]) > tolerance)
+ return 0;
+
+ /* Check bulk velocity (if non-zero)*/
+ if (fabsf(ma->vel[0] + mb->vel[0]) > 0.f &&
+ fabsf(ma->vel[0] - mb->vel[0]) / fabsf(ma->vel[0] + mb->vel[0]) >
+ tolerance)
+ return 0;
+ if (fabsf(ma->vel[1] + mb->vel[1]) > 0.f &&
+ fabsf(ma->vel[1] - mb->vel[1]) / fabsf(ma->vel[1] + mb->vel[1]) >
+ tolerance)
+ return 0;
+ if (fabsf(ma->vel[2] + mb->vel[2]) > 0.f &&
+ fabsf(ma->vel[2] - mb->vel[2]) / fabsf(ma->vel[2] + mb->vel[2]) >
+ tolerance)
+ return 0;
+
+ /* Check mass */
+ if (fabsf(ma->mass - mb->mass) / fabsf(ma->mass + mb->mass) > tolerance)
+ return 0;
- /* Multipole mass */
- float mass;
+ /* All is good */
+ return 1;
+}
+
+/**
+ * @brief Applies the forces due to particles j onto particles i directly.
+ *
+ * @param gparts_i The #gpart to update.
+ * @param gcount_i The number of particles to update.
+ * @param gparts_j The #gpart that source the gravity field.
+ * @param gcount_j The number of sources.
+ */
+INLINE static void multipole_P2P(struct gpart *gparts_i, int gcount_i,
+ const struct gpart *gparts_j, int gcount_j) {}
+
+/**
+ * @brief Constructs the #multipole of a bunch of particles around their
+ * centre of mass.
+ *
+ * Corresponds to equation (28c).
+ *
+ * @param m The #multipole (content will be overwritten).
+ * @param gparts The #gpart.
+ * @param gcount The number of particles.
+ */
+INLINE static void multipole_P2M(struct multipole *m,
+ const struct gpart *gparts, int gcount) {
#if const_gravity_multipole_order >= 2
- /* Quadrupole terms */
- float I_xx, I_yy, I_zz;
- float I_xy, I_xz, I_yz;
+#error "Implementation of P2M kernel missing for this order."
#endif
-};
+
+ /* Temporary variables */
+ double mass = 0.0;
+ double com[3] = {0.0, 0.0, 0.0};
+ float vel[3] = {0.f, 0.f, 0.f};
+
+ /* Collect the particle data. */
+ for (int k = 0; k < gcount; k++) {
+ const float m = gparts[k].mass;
+
+ mass += m;
+ com[0] += gparts[k].x[0] * m;
+ com[1] += gparts[k].x[1] * m;
+ com[2] += gparts[k].x[2] * m;
+ vel[0] += gparts[k].v_full[0] * m;
+ vel[1] += gparts[k].v_full[1] * m;
+ vel[2] += gparts[k].v_full[2] * m;
+ }
+
+ const double imass = 1.0 / mass;
+
+ /* Store the data on the multipole. */
+ m->mass = mass;
+ m->CoM[0] = com[0] * imass;
+ m->CoM[1] = com[1] * imass;
+ m->CoM[2] = com[2] * imass;
+ m->vel[0] = vel[0] * imass;
+ m->vel[1] = vel[1] * imass;
+ m->vel[2] = vel[2] * imass;
+}
+
+/**
+ * @brief Creates a copy of #multipole shifted to a new location.
+ *
+ * Corresponds to equation (28d).
+ *
+ * @param m_a The #multipole copy (content will be overwritten).
+ * @param m_b The #multipole to shift.
+ * @param pos_a The position to which m_b will be shifted.
+ * @param pos_b The current postion of the multipole to shift.
+ * @param periodic Is the calculation periodic ?
+ */
+INLINE static void multipole_M2M(struct multipole *m_a,
+ const struct multipole *m_b,
+ const double pos_a[3], const double pos_b[3],
+ int periodic) {
+
+ m_a->mass = m_b->mass;
+
+ m_a->vel[0] = m_b->vel[0];
+ m_a->vel[1] = m_b->vel[1];
+ m_a->vel[2] = m_b->vel[2];
+}
+/**
+ * @brief Compute the field tensors due to a multipole.
+ *
+ * Corresponds to equation (28b).
+ *
+ * @param l_a The field tensor to compute.
+ * @param m_b The multipole creating the field.
+ * @param pos_a The position of the field tensor.
+ * @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); */
+}
+
+#if 0
/* Multipole function prototypes. */
void multipole_add(struct multipole *m_sum, const struct multipole *m_term);
@@ -137,4 +395,6 @@ __attribute__((always_inline)) INLINE static void multipole_iact_mp(
/* #endif */
}
+#endif
+
#endif /* SWIFT_MULTIPOLE_H */
diff --git a/src/space.c b/src/space.c
index 503727cfef81bffabd1d1cf568be6954bad175cf..68d96397a3c85ed11d5e506519af699dfdd2cecf 100644
--- a/src/space.c
+++ b/src/space.c
@@ -51,6 +51,7 @@
#include "lock.h"
#include "memswap.h"
#include "minmax.h"
+#include "multipole.h"
#include "runner.h"
#include "stars.h"
#include "threadpool.h"
@@ -179,14 +180,27 @@ int space_getsid(struct space *s, struct cell **ci, struct cell **cj,
* @param rec_end Pointer to the end of the list of cells to recycle.
*/
void space_rebuild_recycle_rec(struct space *s, struct cell *c,
- struct cell **rec_begin, struct cell **rec_end) {
+ struct cell **cell_rec_begin,
+ struct cell **cell_rec_end,
+ struct multipole **multipole_rec_begin,
+ struct multipole **multipole_rec_end) {
if (c->split)
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) {
- space_rebuild_recycle_rec(s, c->progeny[k], rec_begin, rec_end);
- c->progeny[k]->next = *rec_begin;
- *rec_begin = c->progeny[k];
- if (*rec_end == NULL) *rec_end = *rec_begin;
+ space_rebuild_recycle_rec(s, c->progeny[k], cell_rec_begin,
+ cell_rec_end, multipole_rec_begin,
+ multipole_rec_end);
+
+ c->progeny[k]->next = *cell_rec_begin;
+ *cell_rec_begin = c->progeny[k];
+ c->progeny[k]->multipole->next = *multipole_rec_begin;
+ *multipole_rec_begin = c->progeny[k]->multipole;
+
+ if (*cell_rec_end == NULL) *cell_rec_end = *cell_rec_begin;
+ if (*multipole_rec_end == NULL)
+ *multipole_rec_end = *multipole_rec_begin;
+
+ c->progeny[k]->multipole = NULL;
c->progeny[k] = NULL;
}
}
@@ -199,9 +213,13 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
for (int k = 0; k < num_elements; k++) {
struct cell *c = &cells[k];
- struct cell *rec_begin = NULL, *rec_end = NULL;
- space_rebuild_recycle_rec(s, c, &rec_begin, &rec_end);
- if (rec_begin != NULL) space_recycle_list(s, rec_begin, rec_end);
+ struct cell *cell_rec_begin = NULL, *cell_rec_end = NULL;
+ struct multipole *multipole_rec_begin = NULL, *multipole_rec_end = NULL;
+ space_rebuild_recycle_rec(s, c, &cell_rec_begin, &cell_rec_end,
+ &multipole_rec_begin, &multipole_rec_end);
+ if (cell_rec_begin != NULL)
+ space_recycle_list(s, cell_rec_begin, cell_rec_end, multipole_rec_begin,
+ multipole_rec_end);
c->sorts = NULL;
c->nr_tasks = 0;
c->density = NULL;
@@ -362,6 +380,7 @@ void space_regrid(struct space *s, int verbose) {
threadpool_map(&s->e->threadpool, space_rebuild_recycle_mapper,
s->cells_top, s->nr_cells, sizeof(struct cell), 100, s);
free(s->cells_top);
+ free(s->multipoles_top);
s->maxdepth = 0;
}
@@ -377,19 +396,35 @@ void space_regrid(struct space *s, int verbose) {
s->tot_cells = s->nr_cells = cdim[0] * cdim[1] * cdim[2];
if (posix_memalign((void *)&s->cells_top, cell_align,
s->nr_cells * sizeof(struct cell)) != 0)
- error("Failed to allocate cells.");
+ error("Failed to allocate top-level cells.");
bzero(s->cells_top, s->nr_cells * sizeof(struct cell));
+ /* Allocate the multipoles for the top-level cells. */
+ if (s->gravity) {
+ if (posix_memalign((void *)&s->multipoles_top, multipole_align,
+ s->nr_cells * sizeof(struct multipole)) != 0)
+ error("Failed to allocate top-level multipoles.");
+ bzero(s->multipoles_top, s->nr_cells * sizeof(struct multipole));
+ }
+
/* Set the cells' locks */
- for (int k = 0; k < s->nr_cells; k++)
+ for (int k = 0; k < s->nr_cells; k++) {
if (lock_init(&s->cells_top[k].lock) != 0)
- error("Failed to init spinlock.");
+ error("Failed to init spinlock for hydro.");
+ if (lock_init(&s->cells_top[k].glock) != 0)
+ error("Failed to init spinlock for gravity.");
+ if (lock_init(&s->cells_top[k].mlock) != 0)
+ error("Failed to init spinlock for multipoles.");
+ if (lock_init(&s->cells_top[k].slock) != 0)
+ error("Failed to init spinlock for stars.");
+ }
/* Set the cell location and sizes. */
for (int i = 0; i < cdim[0]; i++)
for (int j = 0; j < cdim[1]; j++)
for (int k = 0; k < cdim[2]; k++) {
- struct cell *restrict c = &s->cells_top[cell_getid(cdim, i, j, k)];
+ const size_t cid = cell_getid(cdim, i, j, k);
+ struct cell *restrict c = &s->cells_top[cid];
c->loc[0] = i * s->width[0];
c->loc[1] = j * s->width[1];
c->loc[2] = k * s->width[2];
@@ -403,7 +438,7 @@ void space_regrid(struct space *s, int verbose) {
c->scount = 0;
c->super = c;
c->ti_old = ti_current;
- lock_init(&c->lock);
+ if (s->gravity) c->multipole = &s->multipoles_top[cid];
}
/* Be verbose about the change. */
@@ -900,6 +935,12 @@ void space_rebuild(struct space *s, int verbose) {
sure that the parts in each cell are ok. */
space_split(s, cells_top, s->nr_cells, verbose);
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that the multipole construction went OK */
+ for (int k = 0; k < s->nr_cells; k++)
+ cell_check_multipole(&s->cells_top[k], NULL);
+#endif
+
if (verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
clocks_getunit());
@@ -2011,7 +2052,7 @@ void space_split_recursive(struct space *s, struct cell *c,
/* Remove any progeny with zero parts. */
struct cell_buff *progeny_buff = buff, *progeny_gbuff = gbuff,
*progeny_sbuff = sbuff;
- for (int k = 0; k < 8; k++)
+ for (int k = 0; k < 8; k++) {
if (c->progeny[k]->count == 0 && c->progeny[k]->gcount == 0 &&
c->progeny[k]->scount == 0) {
space_recycle(s, c->progeny[k]);
@@ -2028,7 +2069,41 @@ void space_split_recursive(struct space *s, struct cell *c,
if (c->progeny[k]->maxdepth > maxdepth)
maxdepth = c->progeny[k]->maxdepth;
}
+ }
+ /* Deal with multipole */
+ if (s->gravity) {
+
+ /* Reset everything */
+ multipole_init(c->multipole);
+
+ /* Compute CoM of all progenies */
+ double CoM[3] = {0., 0., 0.};
+ double mass = 0.;
+
+ for (int k = 0; k < 8; ++k) {
+ if (c->progeny[k] != NULL) {
+ const struct multipole *m = c->progeny[k]->multipole;
+ CoM[0] += m->CoM[0] * m->mass;
+ CoM[1] += m->CoM[1] * m->mass;
+ CoM[2] += m->CoM[2] * m->mass;
+ mass += m->mass;
+ }
+ }
+ c->multipole->CoM[0] = CoM[0] / mass;
+ c->multipole->CoM[1] = CoM[1] / mass;
+ c->multipole->CoM[2] = CoM[2] / mass;
+
+ /* Now shift progeny multipoles and add them up */
+ struct multipole temp;
+ for (int k = 0; k < 8; ++k) {
+ if (c->progeny[k] != NULL) {
+ const struct multipole *m = c->progeny[k]->multipole;
+ multipole_M2M(&temp, m, c->multipole->CoM, m->CoM, s->periodic);
+ multipole_add(c->multipole, &temp);
+ }
+ }
+ }
}
/* Otherwise, collect the data for this cell. */
@@ -2070,6 +2145,9 @@ void space_split_recursive(struct space *s, struct cell *c,
if (ti_end < ti_end_min) ti_end_min = ti_end;
if (ti_end > ti_end_max) ti_end_max = ti_end;
}
+
+ /* Construct the multipole and the centre of mass*/
+ multipole_P2M(c->multipole, c->gparts, c->gcount);
}
/* Set the values for this cell. */
@@ -2137,8 +2215,9 @@ void space_split_mapper(void *map_data, int num_cells, void *extra_data) {
void space_recycle(struct space *s, struct cell *c) {
/* Clear the cell. */
- if (lock_destroy(&c->lock) != 0 || lock_destroy(&c->glock) != 0)
- error("Failed to destroy spinlock.");
+ if (lock_destroy(&c->lock) != 0 || lock_destroy(&c->glock) != 0 ||
+ lock_destroy(&c->mlock) != 0 || lock_destroy(&c->slock) != 0)
+ error("Failed to destroy spinlocks.");
/* Clear this cell's sort arrays. */
if (c->sort != NULL) free(c->sort);
@@ -2146,6 +2225,12 @@ void space_recycle(struct space *s, struct cell *c) {
/* Lock the space. */
lock_lock(&s->lock);
+ /* Hook the multipole back in the buffer */
+ if (s->gravity) {
+ c->multipole->next = s->multipoles_sub;
+ s->multipoles_sub = c->multipole;
+ }
+
/* Hook this cell into the buffer. */
c->next = s->cells_sub;
s->cells_sub = c;
@@ -2159,22 +2244,32 @@ void space_recycle(struct space *s, struct cell *c) {
* @brief Return a list of used cells to the buffer of unused sub-cells.
*
* @param s The #space.
- * @param list_begin Pointer to the first #cell in the linked list of
+ * @param cell_list_begin Pointer to the first #cell in the linked list of
* cells joined by their @c next pointers.
- * @param list_end Pointer to the last #cell in the linked list of
+ * @param cell_list_end Pointer to the last #cell in the linked list of
* cells joined by their @c next pointers. It is assumed that this
* cell's @c next pointer is @c NULL.
+ * @param multipole_list_begin Pointer to the first #multipole in the linked
+ * list of
+ * multipoles joined by their @c next pointers.
+ * @param multipole_list_end Pointer to the last #multipole in the linked list
+ * of
+ * multipoles joined by their @c next pointers. It is assumed that this
+ * multipole's @c next pointer is @c NULL.
*/
-void space_recycle_list(struct space *s, struct cell *list_begin,
- struct cell *list_end) {
+void space_recycle_list(struct space *s, struct cell *cell_list_begin,
+ struct cell *cell_list_end,
+ struct multipole *multipole_list_begin,
+ struct multipole *multipole_list_end) {
int count = 0;
/* Clean up the list of cells. */
- for (struct cell *c = list_begin; c != NULL; c = c->next) {
+ for (struct cell *c = cell_list_begin; c != NULL; c = c->next) {
/* Clear the cell. */
- if (lock_destroy(&c->lock) != 0 || lock_destroy(&c->glock) != 0)
- error("Failed to destroy spinlock.");
+ if (lock_destroy(&c->lock) != 0 || lock_destroy(&c->glock) != 0 ||
+ lock_destroy(&c->mlock) != 0 || lock_destroy(&c->slock) != 0)
+ error("Failed to destroy spinlocks.");
/* Clear this cell's sort arrays. */
if (c->sort != NULL) free(c->sort);
@@ -2186,11 +2281,15 @@ void space_recycle_list(struct space *s, struct cell *list_begin,
/* Lock the space. */
lock_lock(&s->lock);
- /* Hook this cell into the buffer. */
- list_end->next = s->cells_sub;
- s->cells_sub = list_begin;
+ /* Hook the cells into the buffer. */
+ cell_list_end->next = s->cells_sub;
+ s->cells_sub = cell_list_begin;
s->tot_cells -= count;
+ /* Hook the multipoles into the buffer. */
+ multipole_list_end->next = s->multipoles_sub;
+ s->multipoles_sub = multipole_list_begin;
+
/* Unlock the space. */
lock_unlock_blind(&s->lock);
}
@@ -2214,7 +2313,7 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells) {
/* For each requested cell... */
for (int j = 0; j < nr_cells; j++) {
- /* Is the buffer empty? */
+ /* Is the cell buffer empty? */
if (s->cells_sub == NULL) {
if (posix_memalign((void *)&s->cells_sub, cell_align,
space_cellallocchunk * sizeof(struct cell)) != 0)
@@ -2226,10 +2325,28 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells) {
s->cells_sub[space_cellallocchunk - 1].next = NULL;
}
+ /* Is the multipole buffer empty? */
+ if (s->gravity && s->multipoles_sub == NULL) {
+ if (posix_memalign((void *)&s->multipoles_sub, multipole_align,
+ space_cellallocchunk * sizeof(struct multipole)) != 0)
+ error("Failed to allocate more multipoles.");
+
+ /* Constructed a linked list */
+ for (int k = 0; k < space_cellallocchunk - 1; k++)
+ s->multipoles_sub[k].next = &s->multipoles_sub[k + 1];
+ s->multipoles_sub[space_cellallocchunk - 1].next = NULL;
+ }
+
/* Pick off the next cell. */
cells[j] = s->cells_sub;
s->cells_sub = cells[j]->next;
s->tot_cells += 1;
+
+ /* Hook the multipole */
+ if (s->gravity) {
+ cells[j]->multipole = s->multipoles_sub;
+ s->multipoles_sub = cells[j]->multipole->next;
+ }
}
/* Unlock the space. */
@@ -2237,9 +2354,12 @@ void space_getcells(struct space *s, int nr_cells, struct cell **cells) {
/* Init some things in the cell we just got. */
for (int j = 0; j < nr_cells; j++) {
+ struct multipole *temp = cells[j]->multipole;
bzero(cells[j], sizeof(struct cell));
+ cells[j]->multipole = temp;
cells[j]->nodeID = -1;
- if (lock_init(&cells[j]->lock) != 0 || lock_init(&cells[j]->glock) != 0)
+ if (lock_init(&cells[j]->lock) != 0 || lock_init(&cells[j]->glock) != 0 ||
+ lock_init(&cells[j]->mlock) != 0 || lock_init(&cells[j]->slock) != 0)
error("Failed to initialize cell spinlocks.");
}
}
@@ -2535,7 +2655,7 @@ void space_init(struct space *s, const struct swift_params *params,
/* Init the space lock. */
if (lock_init(&s->lock) != 0) error("Failed to create space spin-lock.");
- /* Build the cells and the tasks. */
+ /* Build the cells recursively. */
if (!dry_run) space_regrid(s, verbose);
}
@@ -2707,6 +2827,7 @@ void space_clean(struct space *s) {
for (int i = 0; i < s->nr_cells; ++i) cell_clean(&s->cells_top[i]);
free(s->cells_top);
+ free(s->multipoles_top);
free(s->parts);
free(s->xparts);
free(s->gparts);
diff --git a/src/space.h b/src/space.h
index 327f0225085663d1859bcdbd82a6ac2583fc2b43..cb16d7df56d4fb4c9f0206de7abb8216aacc1de9 100644
--- a/src/space.h
+++ b/src/space.h
@@ -102,6 +102,12 @@ struct space {
/*! Buffer of unused cells for the sub-cells. */
struct cell *cells_sub;
+ /*! The multipoles associated with the top-level (level 0) cells */
+ struct multipole *multipoles_top;
+
+ /*! Buffer of unused multipoles for the sub-cells. */
+ struct multipole *multipoles_sub;
+
/*! The total number of parts in the space. */
size_t nr_parts, size_parts;
@@ -186,8 +192,10 @@ void space_sparts_sort_mapper(void *map_data, int num_elements,
void *extra_data);
void space_rebuild(struct space *s, int verbose);
void space_recycle(struct space *s, struct cell *c);
-void space_recycle_list(struct space *s, struct cell *list_begin,
- struct cell *list_end);
+void space_recycle_list(struct space *s, struct cell *cell_list_begin,
+ struct cell *cell_list_end,
+ struct multipole *multipole_list_begin,
+ struct multipole *multipole_list_end);
void space_split(struct space *s, struct cell *cells, int nr_cells,
int verbose);
void space_split_mapper(void *map_data, int num_elements, void *extra_data);