From 528af56a3bcfe1bd85c8160aef848b45d5c35e58 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Sat, 25 Apr 2020 14:58:42 +0200
Subject: [PATCH] Do not store the dipole terms since they are always zero.
---
src/cell.c | 5 -
src/multipole.h | 366 +++++++++++++++++++++--------------------
src/multipole_struct.h | 4 +-
src/space.c | 5 -
4 files changed, 187 insertions(+), 193 deletions(-)
diff --git a/src/cell.c b/src/cell.c
index 20631e5030..54bb5bd2c3 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -2283,11 +2283,6 @@ void cell_make_multipoles(struct cell *c, integertime_t ti_current,
/* Take minimum of both limits */
c->grav.multipole->r_max = min(r_max, sqrt(dx * dx + dy * dy + dz * dz));
- /* We know the first-order multipole (dipole) is 0. */
- c->grav.multipole->m_pole.M_100 = 0.f;
- c->grav.multipole->m_pole.M_010 = 0.f;
- c->grav.multipole->m_pole.M_001 = 0.f;
-
/* Compute the multipole power */
gravity_multipole_compute_power(&c->grav.multipole->m_pole);
diff --git a/src/multipole.h b/src/multipole.h
index 07825bfe15..950ec991da 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -305,8 +305,7 @@ __attribute__((nonnull)) INLINE static void gravity_multipole_print(
printf("M_000= %12.5e\n", m->M_000);
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
printf("-------------------------\n");
- printf("M_100= %12.5e M_010= %12.5e M_001= %12.5e\n", m->M_100, m->M_010,
- m->M_001);
+ printf("M_100= %12.5e M_010= %12.5e M_001= %12.5e\n", 0., 0., 0.);
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
printf("-------------------------\n");
@@ -364,10 +363,10 @@ __attribute__((nonnull)) INLINE static void gravity_multipole_add(
ma->M_000 += mb->M_000;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
- /* Add 1st order terms */
- ma->M_100 += mb->M_100;
- ma->M_010 += mb->M_010;
- ma->M_001 += mb->M_001;
+ /* Add 1st order terms (all 0 since we expand around CoM) */
+ /* ma->M_100 += mb->M_100; */
+ /* ma->M_010 += mb->M_010; */
+ /* ma->M_001 += mb->M_001; */
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* Add 2nd order terms */
@@ -527,27 +526,9 @@ __attribute__((nonnull)) INLINE static int gravity_multipole_equal(
return 0;
}
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
- /* Manhattan Norm of 1st order terms */
- const float order1_norm = fabsf(ma->M_001) + fabsf(mb->M_001) +
- fabsf(ma->M_010) + fabsf(mb->M_010) +
- fabsf(ma->M_100) + fabsf(mb->M_100);
-
- /* Compare 1st order terms above 1% of norm */
- if (fabsf(ma->M_001 + mb->M_001) > 0.01f * order1_norm &&
- fabsf(ma->M_001 - mb->M_001) / fabsf(ma->M_001 + mb->M_001) > tolerance) {
- message("M_001 term different");
- return 0;
- }
- if (fabsf(ma->M_010 + mb->M_010) > 0.01f * order1_norm &&
- fabsf(ma->M_010 - mb->M_010) / fabsf(ma->M_010 + mb->M_010) > tolerance) {
- message("M_010 term different");
- return 0;
- }
- if (fabsf(ma->M_100 + mb->M_100) > 0.01f * order1_norm &&
- fabsf(ma->M_100 - mb->M_100) / fabsf(ma->M_100 + mb->M_100) > tolerance) {
- message("M_100 term different");
- return 0;
- }
+ /* Manhattan Norm of 1st order terms */
+ /* Nothing to do here all the 1st order terms are 0 since we expand around
+ * CoM */
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* Manhattan Norm of 2nd order terms */
@@ -903,12 +884,13 @@ __attribute__((nonnull)) INLINE static void gravity_multipole_compute_power(
m->power[0] = m->M_000;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
- /* 1st order terms */
- power[1] += m->M_001 * m->M_001;
- power[1] += m->M_010 * m->M_010;
- power[1] += m->M_100 * m->M_100;
+ /* 1st order terms (all 0 since we expand around CoM) */
+ // power[1] += m->M_001 * m->M_001;
+ // power[1] += m->M_010 * m->M_010;
+ // power[1] += m->M_100 * m->M_100;
- m->power[1] = sqrt(power[1]);
+ // m->power[1] = sqrt(power[1]);
+ m->power[1] = 0.;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* 2nd order terms */
@@ -1177,12 +1159,6 @@ __attribute__((nonnull)) INLINE static void gravity_P2M(
#endif
}
-#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
-
- /* We know the first-order multipole (dipole) is 0. */
- M_100 = M_010 = M_001 = 0.f;
-#endif
-
/* Store the data on the multipole. */
multi->r_max = sqrt(r_max2);
multi->CoM[0] = com[0];
@@ -1203,10 +1179,10 @@ __attribute__((nonnull)) INLINE static void gravity_P2M(
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
- /* 1st order terms */
- multi->m_pole.M_100 = M_100;
- multi->m_pole.M_010 = M_010;
- multi->m_pole.M_001 = M_001;
+ /* 1st order terms (all 0 since we expand around CoM) */
+ // multi->m_pole.M_100 = M_100;
+ // multi->m_pole.M_010 = M_010;
+ // multi->m_pole.M_001 = M_001;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
@@ -1312,228 +1288,256 @@ __attribute__((nonnull)) INLINE static void gravity_M2M(
const double dx[3] = {pos_a[0] - pos_b[0], pos_a[1] - pos_b[1],
pos_a[2] - pos_b[2]};
- /* Shift 1st order term */
- m_a->M_100 = m_b->M_100 + X_100(dx) * m_b->M_000;
- m_a->M_010 = m_b->M_010 + X_010(dx) * m_b->M_000;
- m_a->M_001 = m_b->M_001 + X_001(dx) * m_b->M_000;
+ /* Shift 1st order term (all 0 (after add) since we expand around CoM) */
+ // m_a->M_100 = m_b->M_100 + X_100(dx) * m_b->M_000;
+ // m_a->M_010 = m_b->M_010 + X_010(dx) * m_b->M_000;
+ // m_a->M_001 = m_b->M_001 + X_001(dx) * m_b->M_000;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
- /* Shift 2nd order term */
- m_a->M_200 = m_b->M_200 + X_100(dx) * m_b->M_100 + X_200(dx) * m_b->M_000;
- m_a->M_020 = m_b->M_020 + X_010(dx) * m_b->M_010 + X_020(dx) * m_b->M_000;
- m_a->M_002 = m_b->M_002 + X_001(dx) * m_b->M_001 + X_002(dx) * m_b->M_000;
- m_a->M_110 = m_b->M_110 + X_100(dx) * m_b->M_010 + X_010(dx) * m_b->M_100 +
- X_110(dx) * m_b->M_000;
- m_a->M_101 = m_b->M_101 + X_100(dx) * m_b->M_001 + X_001(dx) * m_b->M_100 +
- X_101(dx) * m_b->M_000;
- m_a->M_011 = m_b->M_011 + X_010(dx) * m_b->M_001 + X_001(dx) * m_b->M_010 +
- X_011(dx) * m_b->M_000;
+ /* Shift 2nd order term (Simplify 1st order terms that are all 0) */
+ m_a->M_200 = m_b->M_200 + /*X_100(dx) * m_b->M_100 +*/ X_200(dx) * m_b->M_000;
+ m_a->M_020 = m_b->M_020 + /*X_010(dx) * m_b->M_010 +*/ X_020(dx) * m_b->M_000;
+ m_a->M_002 = m_b->M_002 + /*X_001(dx) * m_b->M_001 +*/ X_002(dx) * m_b->M_000;
+ m_a->M_110 = m_b->M_110 +
+ /*X_100(dx) * m_b->M_010 + X_010(dx) * m_b->M_100 +*/ X_110(dx) *
+ m_b->M_000;
+ m_a->M_101 = m_b->M_101 +
+ /*X_100(dx) * m_b->M_001 + X_001(dx) * m_b->M_100 +*/ X_101(dx) *
+ m_b->M_000;
+ m_a->M_011 = m_b->M_011 +
+ /*X_010(dx) * m_b->M_001 + X_001(dx) * m_b->M_010 +*/ X_011(dx) *
+ m_b->M_000;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
- /* Shift 3rd order term */
- m_a->M_300 = m_b->M_300 + X_100(dx) * m_b->M_200 + X_200(dx) * m_b->M_100 +
- X_300(dx) * m_b->M_000;
- m_a->M_030 = m_b->M_030 + X_010(dx) * m_b->M_020 + X_020(dx) * m_b->M_010 +
- X_030(dx) * m_b->M_000;
- m_a->M_003 = m_b->M_003 + X_001(dx) * m_b->M_002 + X_002(dx) * m_b->M_001 +
- X_003(dx) * m_b->M_000;
- m_a->M_210 = m_b->M_210 + X_100(dx) * m_b->M_110 + X_010(dx) * m_b->M_200 +
- X_200(dx) * m_b->M_010 + X_110(dx) * m_b->M_100 +
- X_210(dx) * m_b->M_000;
+ /* Shift 3rd order term (Simplify 1st order terms that are all 0) */
+ m_a->M_300 = m_b->M_300 + X_100(dx) * m_b->M_200 +
+ /*X_200(dx) * m_b->M_100 +*/ X_300(dx) * m_b->M_000;
+ m_a->M_030 = m_b->M_030 + X_010(dx) * m_b->M_020 +
+ /*X_020(dx) * m_b->M_010 +*/ X_030(dx) * m_b->M_000;
+ m_a->M_003 = m_b->M_003 + X_001(dx) * m_b->M_002 +
+ /*X_002(dx) * m_b->M_001 +*/ X_003(dx) * m_b->M_000;
+ m_a->M_210 =
+ m_b->M_210 + X_100(dx) * m_b->M_110 + X_010(dx) * m_b->M_200 +
+ /*X_200(dx) * m_b->M_010 + X_110(dx) * m_b->M_100 + */ X_210(dx) *
+ m_b->M_000;
m_a->M_201 = m_b->M_201 + X_100(dx) * m_b->M_101 + X_001(dx) * m_b->M_200 +
- X_200(dx) * m_b->M_001 + X_101(dx) * m_b->M_100 +
- X_201(dx) * m_b->M_000;
+ /*X_200(dx) * m_b->M_001 + X_101(dx) * m_b->M_100 +*/ X_201(dx) *
+ m_b->M_000;
m_a->M_120 = m_b->M_120 + X_010(dx) * m_b->M_110 + X_100(dx) * m_b->M_020 +
- X_020(dx) * m_b->M_100 + X_110(dx) * m_b->M_010 +
- X_120(dx) * m_b->M_000;
+ /*X_020(dx) * m_b->M_100 + X_110(dx) * m_b->M_010 +*/ X_120(dx) *
+ m_b->M_000;
m_a->M_021 = m_b->M_021 + X_010(dx) * m_b->M_011 + X_001(dx) * m_b->M_020 +
- X_020(dx) * m_b->M_001 + X_011(dx) * m_b->M_010 +
- X_021(dx) * m_b->M_000;
+ /*X_020(dx) * m_b->M_001 + X_011(dx) * m_b->M_010 +*/ X_021(dx) *
+ m_b->M_000;
m_a->M_102 = m_b->M_102 + X_001(dx) * m_b->M_101 + X_100(dx) * m_b->M_002 +
- X_002(dx) * m_b->M_100 + X_101(dx) * m_b->M_001 +
- X_102(dx) * m_b->M_000;
+ /*X_002(dx) * m_b->M_100 + X_101(dx) * m_b->M_001 +*/ X_102(dx) *
+ m_b->M_000;
m_a->M_012 = m_b->M_012 + X_001(dx) * m_b->M_011 + X_010(dx) * m_b->M_002 +
- X_002(dx) * m_b->M_010 + X_011(dx) * m_b->M_001 +
- X_012(dx) * m_b->M_000;
+ /*X_002(dx) * m_b->M_010 + X_011(dx) * m_b->M_001 +*/ X_012(dx) *
+ m_b->M_000;
m_a->M_111 = m_b->M_111 + X_100(dx) * m_b->M_011 + X_010(dx) * m_b->M_101 +
- X_001(dx) * m_b->M_110 + X_110(dx) * m_b->M_001 +
- X_101(dx) * m_b->M_010 + X_011(dx) * m_b->M_100 +
+ X_001(dx) * m_b->M_110 + /*X_110(dx) * m_b->M_001 + */
+ /*X_101(dx) * m_b->M_010 + X_011(dx) * m_b->M_100 +*/
X_111(dx) * m_b->M_000;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
- /* Shift 4th order terms */
+ /* Shift 4th order terms (Simplify 1st order terms that are all 0) */
m_a->M_004 = m_b->M_004 + X_001(dx) * m_b->M_003 + X_002(dx) * m_b->M_002 +
- X_003(dx) * m_b->M_001 + X_004(dx) * m_b->M_000;
+ /*X_003(dx) * m_b->M_001 +*/ X_004(dx) * m_b->M_000;
m_a->M_013 = m_b->M_013 + X_001(dx) * m_b->M_012 + X_002(dx) * m_b->M_011 +
- X_003(dx) * m_b->M_010 + X_010(dx) * m_b->M_003 +
- X_011(dx) * m_b->M_002 + X_012(dx) * m_b->M_001 +
+ /*X_003(dx) * m_b->M_010 +*/ X_010(dx) * m_b->M_003 +
+ X_011(dx) * m_b->M_002 + /*X_012(dx) * m_b->M_001 +*/
X_013(dx) * m_b->M_000;
m_a->M_022 = m_b->M_022 + X_001(dx) * m_b->M_021 + X_002(dx) * m_b->M_020 +
X_010(dx) * m_b->M_012 + X_011(dx) * m_b->M_011 +
- X_012(dx) * m_b->M_010 + X_020(dx) * m_b->M_002 +
- X_021(dx) * m_b->M_001 + X_022(dx) * m_b->M_000;
+ /*X_012(dx) * m_b->M_010*/ +X_020(dx) * m_b->M_002 +
+ /*X_021(dx) * m_b->M_001*/ +X_022(dx) * m_b->M_000;
m_a->M_031 = m_b->M_031 + X_001(dx) * m_b->M_030 + X_010(dx) * m_b->M_021 +
X_011(dx) * m_b->M_020 + X_020(dx) * m_b->M_011 +
- X_021(dx) * m_b->M_010 + X_030(dx) * m_b->M_001 +
+ /*X_021(dx) * m_b->M_010 + X_030(dx) * m_b->M_001 + */
X_031(dx) * m_b->M_000;
m_a->M_040 = m_b->M_040 + X_010(dx) * m_b->M_030 + X_020(dx) * m_b->M_020 +
- X_030(dx) * m_b->M_010 + X_040(dx) * m_b->M_000;
+ /*X_030(dx) * m_b->M_010 +*/ X_040(dx) * m_b->M_000;
m_a->M_103 = m_b->M_103 + X_001(dx) * m_b->M_102 + X_002(dx) * m_b->M_101 +
- X_003(dx) * m_b->M_100 + X_100(dx) * m_b->M_003 +
- X_101(dx) * m_b->M_002 + X_102(dx) * m_b->M_001 +
+ /*X_003(dx) * m_b->M_100 +*/ X_100(dx) * m_b->M_003 +
+ X_101(dx) * m_b->M_002 + /*X_102(dx) * m_b->M_001 +*/
X_103(dx) * m_b->M_000;
- m_a->M_112 =
- m_b->M_112 + X_001(dx) * m_b->M_111 + X_002(dx) * m_b->M_110 +
- X_010(dx) * m_b->M_102 + X_011(dx) * m_b->M_101 + X_012(dx) * m_b->M_100 +
- X_100(dx) * m_b->M_012 + X_101(dx) * m_b->M_011 + X_102(dx) * m_b->M_010 +
- X_110(dx) * m_b->M_002 + X_111(dx) * m_b->M_001 + X_112(dx) * m_b->M_000;
- m_a->M_121 =
- m_b->M_121 + X_001(dx) * m_b->M_120 + X_010(dx) * m_b->M_111 +
- X_011(dx) * m_b->M_110 + X_020(dx) * m_b->M_101 + X_021(dx) * m_b->M_100 +
- X_100(dx) * m_b->M_021 + X_101(dx) * m_b->M_020 + X_110(dx) * m_b->M_011 +
- X_111(dx) * m_b->M_010 + X_120(dx) * m_b->M_001 + X_121(dx) * m_b->M_000;
+ m_a->M_112 = m_b->M_112 + X_001(dx) * m_b->M_111 + X_002(dx) * m_b->M_110 +
+ X_010(dx) * m_b->M_102 +
+ X_011(dx) * m_b->M_101 + /*X_012(dx) * m_b->M_100 + */
+ X_100(dx) * m_b->M_012 +
+ X_101(dx) * m_b->M_011 + /*X_102(dx) * m_b->M_010 +*/
+ X_110(dx) * m_b->M_002 +
+ /*X_111(dx) * m_b->M_001 +*/ X_112(dx) * m_b->M_000;
+ m_a->M_121 = m_b->M_121 + X_001(dx) * m_b->M_120 + X_010(dx) * m_b->M_111 +
+ X_011(dx) * m_b->M_110 +
+ X_020(dx) * m_b->M_101 + /*X_021(dx) * m_b->M_100 +*/
+ X_100(dx) * m_b->M_021 + X_101(dx) * m_b->M_020 +
+ X_110(dx) * m_b->M_011 +
+ /*X_111(dx) * m_b->M_010 + X_120(dx) * m_b->M_001 +*/ X_121(dx) *
+ m_b->M_000;
m_a->M_130 = m_b->M_130 + X_010(dx) * m_b->M_120 + X_020(dx) * m_b->M_110 +
- X_030(dx) * m_b->M_100 + X_100(dx) * m_b->M_030 +
- X_110(dx) * m_b->M_020 + X_120(dx) * m_b->M_010 +
+ /*X_030(dx) * m_b->M_100 +*/ X_100(dx) * m_b->M_030 +
+ X_110(dx) * m_b->M_020 + /*X_120(dx) * m_b->M_010 +*/
X_130(dx) * m_b->M_000;
m_a->M_202 = m_b->M_202 + X_001(dx) * m_b->M_201 + X_002(dx) * m_b->M_200 +
X_100(dx) * m_b->M_102 + X_101(dx) * m_b->M_101 +
- X_102(dx) * m_b->M_100 + X_200(dx) * m_b->M_002 +
- X_201(dx) * m_b->M_001 + X_202(dx) * m_b->M_000;
- m_a->M_211 =
- m_b->M_211 + X_001(dx) * m_b->M_210 + X_010(dx) * m_b->M_201 +
- X_011(dx) * m_b->M_200 + X_100(dx) * m_b->M_111 + X_101(dx) * m_b->M_110 +
- X_110(dx) * m_b->M_101 + X_111(dx) * m_b->M_100 + X_200(dx) * m_b->M_011 +
- X_201(dx) * m_b->M_010 + X_210(dx) * m_b->M_001 + X_211(dx) * m_b->M_000;
+ /*X_102(dx) * m_b->M_100 +*/ X_200(dx) * m_b->M_002 +
+ /*X_201(dx) * m_b->M_001 +*/ X_202(dx) * m_b->M_000;
+ m_a->M_211 = m_b->M_211 + X_001(dx) * m_b->M_210 + X_010(dx) * m_b->M_201 +
+ X_011(dx) * m_b->M_200 + X_100(dx) * m_b->M_111 +
+ X_101(dx) * m_b->M_110 + X_110(dx) * m_b->M_101 +
+ /*X_111(dx) * m_b->M_100 +*/ X_200(dx) * m_b->M_011 +
+ /*X_201(dx) * m_b->M_010 + X_210(dx) * m_b->M_001 +*/ X_211(dx) *
+ m_b->M_000;
m_a->M_220 = m_b->M_220 + X_010(dx) * m_b->M_210 + X_020(dx) * m_b->M_200 +
X_100(dx) * m_b->M_120 + X_110(dx) * m_b->M_110 +
- X_120(dx) * m_b->M_100 + X_200(dx) * m_b->M_020 +
- X_210(dx) * m_b->M_010 + X_220(dx) * m_b->M_000;
+ /*X_120(dx) * m_b->M_100 +*/ X_200(dx) * m_b->M_020 +
+ /*X_210(dx) * m_b->M_010 +*/ X_220(dx) * m_b->M_000;
m_a->M_301 = m_b->M_301 + X_001(dx) * m_b->M_300 + X_100(dx) * m_b->M_201 +
X_101(dx) * m_b->M_200 + X_200(dx) * m_b->M_101 +
- X_201(dx) * m_b->M_100 + X_300(dx) * m_b->M_001 +
+ /*X_201(dx) * m_b->M_100 + X_300(dx) * m_b->M_001 + */
X_301(dx) * m_b->M_000;
m_a->M_310 = m_b->M_310 + X_010(dx) * m_b->M_300 + X_100(dx) * m_b->M_210 +
X_110(dx) * m_b->M_200 + X_200(dx) * m_b->M_110 +
- X_210(dx) * m_b->M_100 + X_300(dx) * m_b->M_010 +
+ /*X_210(dx) * m_b->M_100 + X_300(dx) * m_b->M_010 +*/
X_310(dx) * m_b->M_000;
m_a->M_400 = m_b->M_400 + X_100(dx) * m_b->M_300 + X_200(dx) * m_b->M_200 +
- X_300(dx) * m_b->M_100 + X_400(dx) * m_b->M_000;
+ /*X_300(dx) * m_b->M_100 +*/ X_400(dx) * m_b->M_000;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 4
/* Shift 5th order terms */
m_a->M_005 = m_b->M_005 + X_001(dx) * m_b->M_004 + X_002(dx) * m_b->M_003 +
- X_003(dx) * m_b->M_002 + X_004(dx) * m_b->M_001 +
+ X_003(dx) * m_b->M_002 + /*X_004(dx) * m_b->M_001 +*/
X_005(dx) * m_b->M_000;
m_a->M_014 = m_b->M_014 + X_001(dx) * m_b->M_013 + X_002(dx) * m_b->M_012 +
- X_003(dx) * m_b->M_011 + X_004(dx) * m_b->M_010 +
+ X_003(dx) * m_b->M_011 + /*X_004(dx) * m_b->M_010 +*/
X_010(dx) * m_b->M_004 + X_011(dx) * m_b->M_003 +
- X_012(dx) * m_b->M_002 + X_013(dx) * m_b->M_001 +
+ X_012(dx) * m_b->M_002 + /*X_013(dx) * m_b->M_001 +*/
X_014(dx) * m_b->M_000;
- m_a->M_023 =
- m_b->M_023 + X_001(dx) * m_b->M_022 + X_002(dx) * m_b->M_021 +
- X_003(dx) * m_b->M_020 + X_010(dx) * m_b->M_013 + X_011(dx) * m_b->M_012 +
- X_012(dx) * m_b->M_011 + X_013(dx) * m_b->M_010 + X_020(dx) * m_b->M_003 +
- X_021(dx) * m_b->M_002 + X_022(dx) * m_b->M_001 + X_023(dx) * m_b->M_000;
- m_a->M_032 =
- m_b->M_032 + X_001(dx) * m_b->M_031 + X_002(dx) * m_b->M_030 +
- X_010(dx) * m_b->M_022 + X_011(dx) * m_b->M_021 + X_012(dx) * m_b->M_020 +
- X_020(dx) * m_b->M_012 + X_021(dx) * m_b->M_011 + X_022(dx) * m_b->M_010 +
- X_030(dx) * m_b->M_002 + X_031(dx) * m_b->M_001 + X_032(dx) * m_b->M_000;
+ m_a->M_023 = m_b->M_023 + X_001(dx) * m_b->M_022 + X_002(dx) * m_b->M_021 +
+ X_003(dx) * m_b->M_020 + X_010(dx) * m_b->M_013 +
+ X_011(dx) * m_b->M_012 + X_012(dx) * m_b->M_011 +
+ /*X_013(dx) * m_b->M_010 +*/ X_020(dx) * m_b->M_003 +
+ X_021(dx) * m_b->M_002 +
+ /*X_022(dx) * m_b->M_001 +*/ X_023(dx) * m_b->M_000;
+ m_a->M_032 = m_b->M_032 + X_001(dx) * m_b->M_031 + X_002(dx) * m_b->M_030 +
+ X_010(dx) * m_b->M_022 + X_011(dx) * m_b->M_021 +
+ X_012(dx) * m_b->M_020 + X_020(dx) * m_b->M_012 +
+ X_021(dx) * m_b->M_011 + /*X_022(dx) * m_b->M_010 +*/
+ X_030(dx) * m_b->M_002 +
+ /*X_031(dx) * m_b->M_001 +*/ X_032(dx) * m_b->M_000;
m_a->M_041 = m_b->M_041 + X_001(dx) * m_b->M_040 + X_010(dx) * m_b->M_031 +
X_011(dx) * m_b->M_030 + X_020(dx) * m_b->M_021 +
X_021(dx) * m_b->M_020 + X_030(dx) * m_b->M_011 +
- X_031(dx) * m_b->M_010 + X_040(dx) * m_b->M_001 +
+ /*X_031(dx) * m_b->M_010 + X_040(dx) * m_b->M_001 +*/
X_041(dx) * m_b->M_000;
m_a->M_050 = m_b->M_050 + X_010(dx) * m_b->M_040 + X_020(dx) * m_b->M_030 +
- X_030(dx) * m_b->M_020 + X_040(dx) * m_b->M_010 +
+ X_030(dx) * m_b->M_020 + /* X_040(dx) * m_b->M_010 +*/
X_050(dx) * m_b->M_000;
m_a->M_104 = m_b->M_104 + X_001(dx) * m_b->M_103 + X_002(dx) * m_b->M_102 +
- X_003(dx) * m_b->M_101 + X_004(dx) * m_b->M_100 +
+ X_003(dx) * m_b->M_101 + /*X_004(dx) * m_b->M_100 +*/
X_100(dx) * m_b->M_004 + X_101(dx) * m_b->M_003 +
- X_102(dx) * m_b->M_002 + X_103(dx) * m_b->M_001 +
+ X_102(dx) * m_b->M_002 + /*X_103(dx) * m_b->M_001 +*/
X_104(dx) * m_b->M_000;
- m_a->M_113 =
- m_b->M_113 + X_001(dx) * m_b->M_112 + X_002(dx) * m_b->M_111 +
- X_003(dx) * m_b->M_110 + X_010(dx) * m_b->M_103 + X_011(dx) * m_b->M_102 +
- X_012(dx) * m_b->M_101 + X_013(dx) * m_b->M_100 + X_100(dx) * m_b->M_013 +
- X_101(dx) * m_b->M_012 + X_102(dx) * m_b->M_011 + X_103(dx) * m_b->M_010 +
- X_110(dx) * m_b->M_003 + X_111(dx) * m_b->M_002 + X_112(dx) * m_b->M_001 +
- X_113(dx) * m_b->M_000;
- m_a->M_122 =
- m_b->M_122 + X_001(dx) * m_b->M_121 + X_002(dx) * m_b->M_120 +
- X_010(dx) * m_b->M_112 + X_011(dx) * m_b->M_111 + X_012(dx) * m_b->M_110 +
- X_020(dx) * m_b->M_102 + X_021(dx) * m_b->M_101 + X_022(dx) * m_b->M_100 +
- X_100(dx) * m_b->M_022 + X_101(dx) * m_b->M_021 + X_102(dx) * m_b->M_020 +
- X_110(dx) * m_b->M_012 + X_111(dx) * m_b->M_011 + X_112(dx) * m_b->M_010 +
- X_120(dx) * m_b->M_002 + X_121(dx) * m_b->M_001 + X_122(dx) * m_b->M_000;
+ m_a->M_113 = m_b->M_113 + X_001(dx) * m_b->M_112 + X_002(dx) * m_b->M_111 +
+ X_003(dx) * m_b->M_110 + X_010(dx) * m_b->M_103 +
+ X_011(dx) * m_b->M_102 + X_012(dx) * m_b->M_101 +
+ /*X_013(dx) * m_b->M_100 +*/ X_100(dx) * m_b->M_013 +
+ X_101(dx) * m_b->M_012 +
+ X_102(dx) * m_b->M_011 + /*X_103(dx) * m_b->M_010 +*/
+ X_110(dx) * m_b->M_003 +
+ X_111(dx) * m_b->M_002 + /* X_112(dx) * m_b->M_001 +*/
+ X_113(dx) * m_b->M_000;
+ m_a->M_122 = m_b->M_122 + X_001(dx) * m_b->M_121 + X_002(dx) * m_b->M_120 +
+ X_010(dx) * m_b->M_112 + X_011(dx) * m_b->M_111 +
+ X_012(dx) * m_b->M_110 + X_020(dx) * m_b->M_102 +
+ X_021(dx) * m_b->M_101 + /*X_022(dx) * m_b->M_100 +*/
+ X_100(dx) * m_b->M_022 + X_101(dx) * m_b->M_021 +
+ X_102(dx) * m_b->M_020 + X_110(dx) * m_b->M_012 +
+ X_111(dx) * m_b->M_011 + /*X_112(dx) * m_b->M_010 +*/
+ X_120(dx) * m_b->M_002 +
+ /*X_121(dx) * m_b->M_001+*/ X_122(dx) * m_b->M_000;
m_a->M_131 =
m_b->M_131 + X_001(dx) * m_b->M_130 + X_010(dx) * m_b->M_121 +
X_011(dx) * m_b->M_120 + X_020(dx) * m_b->M_111 + X_021(dx) * m_b->M_110 +
- X_030(dx) * m_b->M_101 + X_031(dx) * m_b->M_100 + X_100(dx) * m_b->M_031 +
+ X_030(dx) * m_b->M_101 +
+ /*X_031(dx) * m_b->M_100 +*/ X_100(dx) * m_b->M_031 +
X_101(dx) * m_b->M_030 + X_110(dx) * m_b->M_021 + X_111(dx) * m_b->M_020 +
- X_120(dx) * m_b->M_011 + X_121(dx) * m_b->M_010 + X_130(dx) * m_b->M_001 +
+ X_120(dx) *
+ m_b->M_011 + /*X_121(dx) * m_b->M_010 + X_130(dx) * m_b->M_001 +*/
X_131(dx) * m_b->M_000;
m_a->M_140 = m_b->M_140 + X_010(dx) * m_b->M_130 + X_020(dx) * m_b->M_120 +
- X_030(dx) * m_b->M_110 + X_040(dx) * m_b->M_100 +
+ X_030(dx) * m_b->M_110 + /*X_040(dx) * m_b->M_100 +*/
X_100(dx) * m_b->M_040 + X_110(dx) * m_b->M_030 +
- X_120(dx) * m_b->M_020 + X_130(dx) * m_b->M_010 +
+ X_120(dx) * m_b->M_020 + /*X_130(dx) * m_b->M_010 +*/
X_140(dx) * m_b->M_000;
- m_a->M_203 =
- m_b->M_203 + X_001(dx) * m_b->M_202 + X_002(dx) * m_b->M_201 +
- X_003(dx) * m_b->M_200 + X_100(dx) * m_b->M_103 + X_101(dx) * m_b->M_102 +
- X_102(dx) * m_b->M_101 + X_103(dx) * m_b->M_100 + X_200(dx) * m_b->M_003 +
- X_201(dx) * m_b->M_002 + X_202(dx) * m_b->M_001 + X_203(dx) * m_b->M_000;
- m_a->M_212 =
- m_b->M_212 + X_001(dx) * m_b->M_211 + X_002(dx) * m_b->M_210 +
- X_010(dx) * m_b->M_202 + X_011(dx) * m_b->M_201 + X_012(dx) * m_b->M_200 +
- X_100(dx) * m_b->M_112 + X_101(dx) * m_b->M_111 + X_102(dx) * m_b->M_110 +
- X_110(dx) * m_b->M_102 + X_111(dx) * m_b->M_101 + X_112(dx) * m_b->M_100 +
- X_200(dx) * m_b->M_012 + X_201(dx) * m_b->M_011 + X_202(dx) * m_b->M_010 +
- X_210(dx) * m_b->M_002 + X_211(dx) * m_b->M_001 + X_212(dx) * m_b->M_000;
+ m_a->M_203 = m_b->M_203 + X_001(dx) * m_b->M_202 + X_002(dx) * m_b->M_201 +
+ X_003(dx) * m_b->M_200 + X_100(dx) * m_b->M_103 +
+ X_101(dx) * m_b->M_102 + X_102(dx) * m_b->M_101 +
+ /*X_103(dx) * m_b->M_100 +*/ X_200(dx) * m_b->M_003 +
+ X_201(dx) * m_b->M_002 +
+ /* X_202(dx) * m_b->M_001 +*/ X_203(dx) * m_b->M_000;
+ m_a->M_212 = m_b->M_212 + X_001(dx) * m_b->M_211 + X_002(dx) * m_b->M_210 +
+ X_010(dx) * m_b->M_202 + X_011(dx) * m_b->M_201 +
+ X_012(dx) * m_b->M_200 + X_100(dx) * m_b->M_112 +
+ X_101(dx) * m_b->M_111 + X_102(dx) * m_b->M_110 +
+ X_110(dx) * m_b->M_102 +
+ X_111(dx) * m_b->M_101 + /*X_112(dx) * m_b->M_100 +*/
+ X_200(dx) * m_b->M_012 +
+ X_201(dx) * m_b->M_011 + /*X_202(dx) * m_b->M_010 +*/
+ X_210(dx) * m_b->M_002 +
+ /*X_211(dx) * m_b->M_001 +*/ X_212(dx) * m_b->M_000;
m_a->M_221 =
m_b->M_221 + X_001(dx) * m_b->M_220 + X_010(dx) * m_b->M_211 +
X_011(dx) * m_b->M_210 + X_020(dx) * m_b->M_201 + X_021(dx) * m_b->M_200 +
X_100(dx) * m_b->M_121 + X_101(dx) * m_b->M_120 + X_110(dx) * m_b->M_111 +
- X_111(dx) * m_b->M_110 + X_120(dx) * m_b->M_101 + X_121(dx) * m_b->M_100 +
+ X_111(dx) * m_b->M_110 +
+ X_120(dx) * m_b->M_101 + /*X_121(dx) * m_b->M_100 +*/
X_200(dx) * m_b->M_021 + X_201(dx) * m_b->M_020 + X_210(dx) * m_b->M_011 +
- X_211(dx) * m_b->M_010 + X_220(dx) * m_b->M_001 + X_221(dx) * m_b->M_000;
- m_a->M_230 =
- m_b->M_230 + X_010(dx) * m_b->M_220 + X_020(dx) * m_b->M_210 +
- X_030(dx) * m_b->M_200 + X_100(dx) * m_b->M_130 + X_110(dx) * m_b->M_120 +
- X_120(dx) * m_b->M_110 + X_130(dx) * m_b->M_100 + X_200(dx) * m_b->M_030 +
- X_210(dx) * m_b->M_020 + X_220(dx) * m_b->M_010 + X_230(dx) * m_b->M_000;
- m_a->M_302 =
- m_b->M_302 + X_001(dx) * m_b->M_301 + X_002(dx) * m_b->M_300 +
- X_100(dx) * m_b->M_202 + X_101(dx) * m_b->M_201 + X_102(dx) * m_b->M_200 +
- X_200(dx) * m_b->M_102 + X_201(dx) * m_b->M_101 + X_202(dx) * m_b->M_100 +
- X_300(dx) * m_b->M_002 + X_301(dx) * m_b->M_001 + X_302(dx) * m_b->M_000;
+ /*X_211(dx) * m_b->M_010 + X_220(dx) * m_b->M_001 +*/ X_221(dx) *
+ m_b->M_000;
+ m_a->M_230 = m_b->M_230 + X_010(dx) * m_b->M_220 + X_020(dx) * m_b->M_210 +
+ X_030(dx) * m_b->M_200 + X_100(dx) * m_b->M_130 +
+ X_110(dx) * m_b->M_120 + X_120(dx) * m_b->M_110 +
+ /*X_130(dx) * m_b->M_100 +*/ X_200(dx) * m_b->M_030 +
+ X_210(dx) * m_b->M_020 +
+ /*X_220(dx) * m_b->M_010 +*/ X_230(dx) * m_b->M_000;
+ m_a->M_302 = m_b->M_302 + X_001(dx) * m_b->M_301 + X_002(dx) * m_b->M_300 +
+ X_100(dx) * m_b->M_202 + X_101(dx) * m_b->M_201 +
+ X_102(dx) * m_b->M_200 + X_200(dx) * m_b->M_102 +
+ X_201(dx) * m_b->M_101 + /*X_202(dx) * m_b->M_100 +*/
+ X_300(dx) * m_b->M_002 +
+ /*X_301(dx) * m_b->M_001 +*/ X_302(dx) * m_b->M_000;
m_a->M_311 =
m_b->M_311 + X_001(dx) * m_b->M_310 + X_010(dx) * m_b->M_301 +
X_011(dx) * m_b->M_300 + X_100(dx) * m_b->M_211 + X_101(dx) * m_b->M_210 +
X_110(dx) * m_b->M_201 + X_111(dx) * m_b->M_200 + X_200(dx) * m_b->M_111 +
- X_201(dx) * m_b->M_110 + X_210(dx) * m_b->M_101 + X_211(dx) * m_b->M_100 +
- X_300(dx) * m_b->M_011 + X_301(dx) * m_b->M_010 + X_310(dx) * m_b->M_001 +
+ X_201(dx) * m_b->M_110 +
+ X_210(dx) * m_b->M_101 + /*X_211(dx) * m_b->M_100 +*/
+ X_300(dx) *
+ m_b->M_011 + /*X_301(dx) * m_b->M_010 + X_310(dx) * m_b->M_001 +*/
X_311(dx) * m_b->M_000;
- m_a->M_320 =
- m_b->M_320 + X_010(dx) * m_b->M_310 + X_020(dx) * m_b->M_300 +
- X_100(dx) * m_b->M_220 + X_110(dx) * m_b->M_210 + X_120(dx) * m_b->M_200 +
- X_200(dx) * m_b->M_120 + X_210(dx) * m_b->M_110 + X_220(dx) * m_b->M_100 +
- X_300(dx) * m_b->M_020 + X_310(dx) * m_b->M_010 + X_320(dx) * m_b->M_000;
+ m_a->M_320 = m_b->M_320 + X_010(dx) * m_b->M_310 + X_020(dx) * m_b->M_300 +
+ X_100(dx) * m_b->M_220 + X_110(dx) * m_b->M_210 +
+ X_120(dx) * m_b->M_200 + X_200(dx) * m_b->M_120 +
+ X_210(dx) * m_b->M_110 + /*X_220(dx) * m_b->M_100 +*/
+ X_300(dx) * m_b->M_020 +
+ /*X_310(dx) * m_b->M_010 +*/ X_320(dx) * m_b->M_000;
m_a->M_401 = m_b->M_401 + X_001(dx) * m_b->M_400 + X_100(dx) * m_b->M_301 +
X_101(dx) * m_b->M_300 + X_200(dx) * m_b->M_201 +
X_201(dx) * m_b->M_200 + X_300(dx) * m_b->M_101 +
- X_301(dx) * m_b->M_100 + X_400(dx) * m_b->M_001 +
+ /*X_301(dx) * m_b->M_100 + X_400(dx) * m_b->M_001 +*/
X_401(dx) * m_b->M_000;
m_a->M_410 = m_b->M_410 + X_010(dx) * m_b->M_400 + X_100(dx) * m_b->M_310 +
X_110(dx) * m_b->M_300 + X_200(dx) * m_b->M_210 +
X_210(dx) * m_b->M_200 + X_300(dx) * m_b->M_110 +
- X_310(dx) * m_b->M_100 + X_400(dx) * m_b->M_010 +
+ /*X_310(dx) * m_b->M_100 + X_400(dx) * m_b->M_010 +*/
X_410(dx) * m_b->M_000;
m_a->M_500 = m_b->M_500 + X_100(dx) * m_b->M_400 + X_200(dx) * m_b->M_300 +
- X_300(dx) * m_b->M_200 + X_400(dx) * m_b->M_100 +
+ X_300(dx) * m_b->M_200 + /*X_400(dx) * m_b->M_100 +*/
X_500(dx) * m_b->M_000;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 5
diff --git a/src/multipole_struct.h b/src/multipole_struct.h
index f6f95576c2..ffd615cd79 100644
--- a/src/multipole_struct.h
+++ b/src/multipole_struct.h
@@ -132,8 +132,8 @@ struct multipole {
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
- /* 1st order terms */
- float M_100, M_010, M_001;
+ /* 1st order terms (all 0 since we expand around CoM) */
+ // float M_100, M_010, M_001;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
diff --git a/src/space.c b/src/space.c
index 5110873789..45dddc1b4b 100644
--- a/src/space.c
+++ b/src/space.c
@@ -3659,11 +3659,6 @@ void space_split_recursive(struct space *s, struct cell *c,
c->grav.multipole->CoM_rebuild[1] = c->grav.multipole->CoM[1];
c->grav.multipole->CoM_rebuild[2] = c->grav.multipole->CoM[2];
- /* We know the first-order multipole (dipole) is 0. */
- c->grav.multipole->m_pole.M_100 = 0.f;
- c->grav.multipole->m_pole.M_010 = 0.f;
- c->grav.multipole->m_pole.M_001 = 0.f;
-
/* Compute the multipole power */
gravity_multipole_compute_power(&c->grav.multipole->m_pole);
--
GitLab