diff --git a/src/engine.c b/src/engine.c
index 1f51380c7c6a2a4bd3cea2470967647b29be1f0c..04124ec5aa8435b2551a19de999f37886020385e 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -3005,7 +3005,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
if (e->policy & engine_policy_self_gravity) {
double mass = 0.;
for (int i = 0; i < e->s->nr_cells; ++i)
- mass += e->s->cells_top[i].multipole->m_pole.mass;
+ mass += e->s->cells_top[i].multipole->m_pole.M_000;
if (fabs(mass - e->s->total_mass) > e->s->total_mass / e->s->nr_gparts)
error(
"Total mass in multipoles does not match particle content. part=%e "
@@ -3104,7 +3104,7 @@ void engine_step(struct engine *e) {
if (e->policy & engine_policy_self_gravity) {
double mass = 0.;
for (int i = 0; i < e->s->nr_cells; ++i)
- mass += e->s->cells_top[i].multipole->m_pole.mass;
+ mass += e->s->cells_top[i].multipole->m_pole.M_000;
if (fabs(mass - e->s->total_mass) > e->s->total_mass / e->s->nr_gparts)
error(
"Total mass in multipoles does not match particle content. part=%e "
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index 7c7c5b9d244534ef3f6b5f509062019bfcd5f9fb..3f1beb9562b054a331c6461484f940a3f9a0afb2 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -141,51 +141,53 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
float rlr_inv, float r2, const float *dx, struct gpart *gp,
const struct multipole *multi) {
- /* Apply the gravitational acceleration. */
- const float r = sqrtf(r2);
- const float ir = 1.f / r;
- const float u = r * rlr_inv;
- const float mrinv3 = multi->mass * ir * ir * ir;
-
- /* Get long-range correction */
- float f_lr;
- kernel_long_grav_eval(u, &f_lr);
-
-#if const_gravity_multipole_order < 2
-
- /* 0th and 1st order terms */
- gp->a_grav[0] += mrinv3 * f_lr * dx[0];
- gp->a_grav[1] += mrinv3 * f_lr * dx[1];
- gp->a_grav[2] += mrinv3 * f_lr * dx[2];
-
-#elif const_gravity_multipole_order == 2
- /* Terms up to 2nd order (quadrupole) */
-
- /* Follows the notation in Bonsai */
- const float mrinv5 = mrinv3 * ir * ir;
- const float mrinv7 = mrinv5 * ir * ir;
-
- const float D1 = -mrinv3;
- const float D2 = 3.f * mrinv5;
- const float D3 = -15.f * mrinv7;
-
- const float q = multi->I_xx + multi->I_yy + multi->I_zz;
- const float qRx =
- multi->I_xx * dx[0] + multi->I_xy * dx[1] + multi->I_xz * dx[2];
- const float qRy =
- multi->I_xy * dx[0] + multi->I_yy * dx[1] + multi->I_yz * dx[2];
- const float qRz =
- multi->I_xz * dx[0] + multi->I_yz * dx[1] + multi->I_zz * dx[2];
- const float qRR = qRx * dx[0] + qRy * dx[1] + qRz * dx[2];
- const float C = D1 + 0.5f * D2 * q + 0.5f * D3 * qRR;
-
- gp->a_grav[0] -= f_lr * (C * dx[0] + D2 * qRx);
- gp->a_grav[1] -= f_lr * (C * dx[1] + D2 * qRy);
- gp->a_grav[2] -= f_lr * (C * dx[2] + D2 * qRz);
-
-#else
-#error "Multipoles of order >2 not yet implemented."
-#endif
+ error("Dead function");
+
+ /* /\* Apply the gravitational acceleration. *\/ */
+ /* const float r = sqrtf(r2); */
+ /* const float ir = 1.f / r; */
+ /* const float u = r * rlr_inv; */
+ /* const float mrinv3 = multi->mass * ir * ir * ir; */
+
+ /* /\* Get long-range correction *\/ */
+ /* float f_lr; */
+ /* kernel_long_grav_eval(u, &f_lr); */
+
+ /* #if const_gravity_multipole_order < 2 */
+
+ /* /\* 0th and 1st order terms *\/ */
+ /* gp->a_grav[0] += mrinv3 * f_lr * dx[0]; */
+ /* gp->a_grav[1] += mrinv3 * f_lr * dx[1]; */
+ /* gp->a_grav[2] += mrinv3 * f_lr * dx[2]; */
+
+ /* #elif const_gravity_multipole_order == 2 */
+ /* /\* Terms up to 2nd order (quadrupole) *\/ */
+
+ /* /\* Follows the notation in Bonsai *\/ */
+ /* const float mrinv5 = mrinv3 * ir * ir; */
+ /* const float mrinv7 = mrinv5 * ir * ir; */
+
+ /* const float D1 = -mrinv3; */
+ /* const float D2 = 3.f * mrinv5; */
+ /* const float D3 = -15.f * mrinv7; */
+
+ /* const float q = multi->I_xx + multi->I_yy + multi->I_zz; */
+ /* const float qRx = */
+ /* multi->I_xx * dx[0] + multi->I_xy * dx[1] + multi->I_xz * dx[2]; */
+ /* const float qRy = */
+ /* multi->I_xy * dx[0] + multi->I_yy * dx[1] + multi->I_yz * dx[2]; */
+ /* const float qRz = */
+ /* multi->I_xz * dx[0] + multi->I_yz * dx[1] + multi->I_zz * dx[2]; */
+ /* const float qRR = qRx * dx[0] + qRy * dx[1] + qRz * dx[2]; */
+ /* const float C = D1 + 0.5f * D2 * q + 0.5f * D3 * qRR; */
+
+ /* gp->a_grav[0] -= f_lr * (C * dx[0] + D2 * qRx); */
+ /* gp->a_grav[1] -= f_lr * (C * dx[1] + D2 * qRy); */
+ /* gp->a_grav[2] -= f_lr * (C * dx[2] + D2 * qRz); */
+
+ /* #else */
+ /* #error "Multipoles of order >2 not yet implemented." */
+ /* #endif */
}
#endif /* SWIFT_DEFAULT_GRAVITY_IACT_H */
diff --git a/src/multipole.h b/src/multipole.h
index 3cbb55301f3175d25de9c5523d02564a19dcd702..613bbd16f6a2402dc21352fb3f61a7281b6ae77b 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -49,10 +49,36 @@ struct pot_tensor {
struct multipole {
- float mass;
-
- /*! Bulk velocity */
+ /* Bulk velocity */
float vel[3];
+
+ /* 0th order terms */
+ float M_000;
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+
+ /* 1st order terms */
+ float M_100, M_010, M_001;
+#endif
+
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+
+ /* 2nd order terms */
+ float M_200, M_020, M_002;
+ float M_110, M_101, M_011;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+
+ /* 3rd order terms */
+ float M_300, M_030, M_003;
+ float M_210, M_201;
+ float M_120, M_021;
+ float M_102, M_012;
+ float M_111;
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
+#error "Missing implementation for order >3"
+#endif
};
/**
@@ -157,9 +183,27 @@ INLINE static void gravity_field_tensors_add(struct gravity_tensors *la,
*/
INLINE static void gravity_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]);
+ printf("M_000= %12.5e\n", m->M_000);
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
+ printf("M_100= %8.5e M_010= %8.5e M_001= %8.5e\n", m->M_100, m->M_010,
+ m->M_001);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
+ printf("M_200= %8.5e M_020= %8.5e M_002= %8.5e\n", m->M_200, m->M_020,
+ m->M_002);
+ printf("M_110= %8.5e M_101= %8.5e M_011= %8.5e\n", m->M_110, m->M_101,
+ m->M_011);
+#endif
+#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
+ printf("M_300= %8.5e M_030= %8.5e M_003= %8.5e\n", m->M_300, m->M_030,
+ m->M_003);
+ printf("M_210= %8.5e M_201= %8.5e M_120= %8.5e\n", m->M_210, m->M_201,
+ m->M_120);
+ printf("M_021= %8.5e M_102= %8.5e M_012= %8.5e\n", m->M_021, m->M_102,
+ m->M_012);
+ printf("M_111= %8.5e\n", m->M_111);
+#endif
}
/**
@@ -171,16 +215,16 @@ INLINE static void gravity_multipole_print(const struct multipole *m) {
INLINE static void gravity_multipole_add(struct multipole *ma,
const struct multipole *mb) {
- const float mass = ma->mass + mb->mass;
- const float imass = 1.f / mass;
+ const float M_000 = ma->M_000 + mb->M_000;
+ const float inv_M_000 = 1.f / M_000;
/* 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;
+ ma->vel[0] = (ma->vel[0] * ma->M_000 + mb->vel[0] * mb->M_000) * inv_M_000;
+ ma->vel[1] = (ma->vel[1] * ma->M_000 + mb->vel[1] * mb->M_000) * inv_M_000;
+ ma->vel[2] = (ma->vel[2] * ma->M_000 + mb->vel[2] * mb->M_000) * inv_M_000;
- /* Add the masses */
- ma->mass = mass;
+ /* Add 0th order terms */
+ ma->M_000 = M_000;
}
/**
@@ -220,8 +264,8 @@ INLINE static int gravity_multipole_equal(const struct multipole *ma,
tolerance)
return 0;
- /* Check mass */
- if (fabsf(ma->mass - mb->mass) / fabsf(ma->mass + mb->mass) > tolerance)
+ /* Check 0th order terms */
+ if (fabsf(ma->M_000 - mb->M_000) / fabsf(ma->M_000 + mb->M_000) > tolerance)
return 0;
/* All is good */
@@ -266,7 +310,7 @@ INLINE static void gravity_P2M(struct gravity_tensors *m,
const double imass = 1.0 / mass;
/* Store the data on the multipole. */
- m->m_pole.mass = mass;
+ m->m_pole.M_000 = mass;
m->CoM[0] = com[0] * imass;
m->CoM[1] = com[1] * imass;
m->CoM[2] = com[2] * imass;
@@ -290,12 +334,13 @@ INLINE static void gravity_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;
-
+ /* Shift bulk velocity */
m_a->vel[0] = m_b->vel[0];
m_a->vel[1] = m_b->vel[1];
m_a->vel[2] = m_b->vel[2];
+
+ /* Shift 0th order term */
+ m_a->M_000 = m_b->M_000;
}
/**
* @brief Compute the field tensors due to a multipole.
@@ -327,13 +372,13 @@ INLINE static void gravity_M2L(struct gravity_tensors *l_a,
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;
+ /* 0st order multipole term */
+ l_a->a_x.F_000 += D_100(dx, dy, dz, r_inv) * m_b->M_000;
+ l_a->a_y.F_000 += D_010(dx, dy, dz, r_inv) * m_b->M_000;
+ l_a->a_z.F_000 += D_001(dx, dy, dz, r_inv) * m_b->M_000;
#ifdef SWIFT_DEBUG_CHECKS
- l_a->mass_interacted += m_b->mass;
+ l_a->mass_interacted += m_b->M_000;
#endif
}
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index d03fdaea71bff848682605ec8aaf14c80a412c5f..bdfa16ec7e5774bb8eee6568eed12d0b22bed10a 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -82,7 +82,7 @@ void runner_dopair_grav_mm(const struct runner *r,
TIMER_TIC;
#ifdef SWIFT_DEBUG_CHECKS
- if (multi_j->mass == 0.0) error("Multipole does not seem to have been set.");
+ if (multi_j->M_000 == 0.f) error("Multipole does not seem to have been set.");
#endif
/* Anything to do here? */
diff --git a/src/space.c b/src/space.c
index 625fe944c488c871d02b14c72d15051b3e53b3c4..b34571368bc7b5fe6644f3af931dfab3812bbf90 100644
--- a/src/space.c
+++ b/src/space.c
@@ -2100,10 +2100,10 @@ void space_split_recursive(struct space *s, struct cell *c,
for (int k = 0; k < 8; ++k) {
if (c->progeny[k] != NULL) {
const struct gravity_tensors *m = c->progeny[k]->multipole;
- CoM[0] += m->CoM[0] * m->m_pole.mass;
- CoM[1] += m->CoM[1] * m->m_pole.mass;
- CoM[2] += m->CoM[2] * m->m_pole.mass;
- mass += m->m_pole.mass;
+ CoM[0] += m->CoM[0] * m->m_pole.M_000;
+ CoM[1] += m->CoM[1] * m->m_pole.M_000;
+ CoM[2] += m->CoM[2] * m->m_pole.M_000;
+ mass += m->m_pole.M_000;
}
}
c->multipole->CoM[0] = CoM[0] / mass;