Commit d40f8ebb authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Merge branch 'master' of gitlab.cosma.dur.ac.uk:swift/swiftsim

parents c5cdfc60 0aaf4486
......@@ -148,10 +148,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
float *restrict f_x, float *restrict f_y, float *restrict f_z,
float *restrict pot) {
/* In the case where the order is < 3, then there is only a monopole term left.
/* In the case where the order is < 2, then there is only a monopole term left.
* We can default to the normal P-P interaction with the mass of the multipole
* and its CoM as the "particle" property */
#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
#if SELF_GRAVITY_MULTIPOLE_ORDER < 2
float f_ij, pot_ij;
runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000,
......@@ -186,8 +186,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
/* *f_z = m->M_001 * d.D_002 + m->M_010 * d.D_011 + m->M_100 * d.D_101 ; */
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* 2nd order contributions */
*f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
......@@ -198,8 +197,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
m->M_101 * d.D_102 + m->M_110 * d.D_111 + m->M_200 * d.D_201;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
/* 3rd order contributions */
*f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
......@@ -215,6 +213,25 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
m->M_120 * d.D_121 + m->M_201 * d.D_202 + m->M_210 * d.D_211 +
m->M_300 * d.D_301;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
/* 4th order contributions */
*f_x += m->M_004 * d.D_104 + m->M_013 * d.D_113 + m->M_022 * d.D_122 +
m->M_031 * d.D_131 + m->M_040 * d.D_140 + m->M_103 * d.D_203 +
m->M_112 * d.D_212 + m->M_121 * d.D_221 + m->M_130 * d.D_230 +
m->M_202 * d.D_302 + m->M_211 * d.D_311 + m->M_220 * d.D_320 +
m->M_301 * d.D_401 + m->M_310 * d.D_410 + m->M_400 * d.D_500;
*f_y += m->M_004 * d.D_014 + m->M_013 * d.D_023 + m->M_022 * d.D_032 +
m->M_031 * d.D_041 + m->M_040 * d.D_050 + m->M_103 * d.D_113 +
m->M_112 * d.D_122 + m->M_121 * d.D_131 + m->M_130 * d.D_140 +
m->M_202 * d.D_212 + m->M_211 * d.D_221 + m->M_220 * d.D_230 +
m->M_301 * d.D_311 + m->M_310 * d.D_320 + m->M_400 * d.D_410;
*f_z += m->M_004 * d.D_005 + m->M_013 * d.D_014 + m->M_022 * d.D_023 +
m->M_031 * d.D_032 + m->M_040 * d.D_041 + m->M_103 * d.D_104 +
m->M_112 * d.D_113 + m->M_121 * d.D_122 + m->M_130 * d.D_131 +
m->M_202 * d.D_203 + m->M_211 * d.D_212 + m->M_220 * d.D_221 +
m->M_301 * d.D_302 + m->M_310 * d.D_311 + m->M_400 * d.D_401;
#endif
/* Take care of the the sign convention */
......@@ -254,10 +271,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
const struct multipole *m, float *restrict f_x, float *restrict f_y,
float *restrict f_z, float *restrict pot) {
/* In the case where the order is < 3, then there is only a monopole term left.
/* In the case where the order is < 2, then there is only a monopole term left.
* We can default to the normal P-P interaction with the mass of the multipole
* and its CoM as the "particle" property */
#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
#if SELF_GRAVITY_MULTIPOLE_ORDER < 2
float f_ij, pot_ij;
runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv,
......@@ -292,8 +309,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
/* *f_z = m->M_001 * d.D_002 + m->M_010 * d.D_011 + m->M_100 * d.D_101 ; */
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* 2nd order contributions */
*f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
......@@ -304,8 +320,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
m->M_101 * d.D_102 + m->M_110 * d.D_111 + m->M_200 * d.D_201;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
/* 3rd order contributions */
*f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
......@@ -321,7 +336,25 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
m->M_120 * d.D_121 + m->M_201 * d.D_202 + m->M_210 * d.D_211 +
m->M_300 * d.D_301;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
/* 4th order contributions */
*f_x += m->M_004 * d.D_104 + m->M_013 * d.D_113 + m->M_022 * d.D_122 +
m->M_031 * d.D_131 + m->M_040 * d.D_140 + m->M_103 * d.D_203 +
m->M_112 * d.D_212 + m->M_121 * d.D_221 + m->M_130 * d.D_230 +
m->M_202 * d.D_302 + m->M_211 * d.D_311 + m->M_220 * d.D_320 +
m->M_301 * d.D_401 + m->M_310 * d.D_410 + m->M_400 * d.D_500;
*f_y += m->M_004 * d.D_014 + m->M_013 * d.D_023 + m->M_022 * d.D_032 +
m->M_031 * d.D_041 + m->M_040 * d.D_050 + m->M_103 * d.D_113 +
m->M_112 * d.D_122 + m->M_121 * d.D_131 + m->M_130 * d.D_140 +
m->M_202 * d.D_212 + m->M_211 * d.D_221 + m->M_220 * d.D_230 +
m->M_301 * d.D_311 + m->M_310 * d.D_320 + m->M_400 * d.D_410;
*f_z += m->M_004 * d.D_005 + m->M_013 * d.D_014 + m->M_022 * d.D_023 +
m->M_031 * d.D_032 + m->M_040 * d.D_041 + m->M_103 * d.D_104 +
m->M_112 * d.D_113 + m->M_121 * d.D_122 + m->M_130 * d.D_131 +
m->M_202 * d.D_203 + m->M_211 * d.D_212 + m->M_220 * d.D_221 +
m->M_301 * d.D_302 + m->M_310 * d.D_311 + m->M_400 * d.D_401;
#endif
/* Take care of the the sign convention */
*f_x *= -1.f;
*f_y *= -1.f;
......
......@@ -150,10 +150,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
float *restrict f_x, float *restrict f_y, float *restrict f_z,
float *restrict pot) {
/* In the case where the order is < 3, then there is only a monopole term left.
/* In the case where the order is < 2, then there is only a monopole term left.
* We can default to the normal P-P interaction with the mass of the multipole
* and its CoM as the "particle" property */
#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
#if SELF_GRAVITY_MULTIPOLE_ORDER < 2
float f_ij, pot_ij;
runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000,
......@@ -179,7 +179,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
*f_z = m->M_000 * d.D_001;
*pot = m->M_000 * d.D_000;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
/* 1st order contributions */
......@@ -191,8 +191,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
/* *pot = m->M_001 * d.D_001 + m->M_010 * d.D_010 + m->M_100 * d.D_100 ; */
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* 2nd order contributions */
*f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
......@@ -205,8 +204,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
m->M_101 * d.D_101 + m->M_110 * d.D_110 + m->M_200 * d.D_200;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
/* 3rd order contributions */
*f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
......@@ -227,7 +225,31 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
m->M_300 * d.D_300;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
/* 4th order contributions */
*f_x += m->M_004 * d.D_104 + m->M_013 * d.D_113 + m->M_022 * d.D_122 +
m->M_031 * d.D_131 + m->M_040 * d.D_140 + m->M_103 * d.D_203 +
m->M_112 * d.D_212 + m->M_121 * d.D_221 + m->M_130 * d.D_230 +
m->M_202 * d.D_302 + m->M_211 * d.D_311 + m->M_220 * d.D_320 +
m->M_301 * d.D_401 + m->M_310 * d.D_410 + m->M_400 * d.D_500;
*f_y += m->M_004 * d.D_014 + m->M_013 * d.D_023 + m->M_022 * d.D_032 +
m->M_031 * d.D_041 + m->M_040 * d.D_050 + m->M_103 * d.D_113 +
m->M_112 * d.D_122 + m->M_121 * d.D_131 + m->M_130 * d.D_140 +
m->M_202 * d.D_212 + m->M_211 * d.D_221 + m->M_220 * d.D_230 +
m->M_301 * d.D_311 + m->M_310 * d.D_320 + m->M_400 * d.D_410;
*f_z += m->M_004 * d.D_005 + m->M_013 * d.D_014 + m->M_022 * d.D_023 +
m->M_031 * d.D_032 + m->M_040 * d.D_041 + m->M_103 * d.D_104 +
m->M_112 * d.D_113 + m->M_121 * d.D_122 + m->M_130 * d.D_131 +
m->M_202 * d.D_203 + m->M_211 * d.D_212 + m->M_220 * d.D_221 +
m->M_301 * d.D_302 + m->M_310 * d.D_311 + m->M_400 * d.D_401;
*pot += m->M_004 * d.D_004 + m->M_013 * d.D_013 + m->M_022 * d.D_022 +
m->M_031 * d.D_031 + m->M_040 * d.D_040 + m->M_103 * d.D_103 +
m->M_112 * d.D_112 + m->M_121 * d.D_121 + m->M_130 * d.D_130 +
m->M_202 * d.D_202 + m->M_211 * d.D_211 + m->M_220 * d.D_220 +
m->M_301 * d.D_301 + m->M_310 * d.D_310 + m->M_400 * d.D_400;
#endif
/* Take care of the the sign convention */
*f_x *= -1.f;
*f_y *= -1.f;
......@@ -263,10 +285,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
const struct multipole *m, float *restrict f_x, float *restrict f_y,
float *restrict f_z, float *restrict pot) {
/* In the case where the order is < 3, then there is only a monopole term left.
/* In the case where the order is < 2, then there is only a monopole term left.
* We can default to the normal P-P interaction with the mass of the multipole
* and its CoM as the "particle" property */
#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
#if SELF_GRAVITY_MULTIPOLE_ORDER < 2
float f_ij, pot_ij;
runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv,
......@@ -292,7 +314,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
*f_z = m->M_000 * d.D_001;
*pot = m->M_000 * d.D_000;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
/* 1st order contributions */
......@@ -304,8 +326,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
/* *pot = m->M_001 * d.D_001 + m->M_010 * d.D_010 + m->M_100 * d.D_100 ; */
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* 2nd order contributions */
*f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
......@@ -318,8 +339,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
m->M_101 * d.D_101 + m->M_110 * d.D_110 + m->M_200 * d.D_200;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
/* 3rd order contributions */
*f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
......@@ -340,7 +360,31 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
m->M_300 * d.D_300;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
/* 4th order contributions */
*f_x += m->M_004 * d.D_104 + m->M_013 * d.D_113 + m->M_022 * d.D_122 +
m->M_031 * d.D_131 + m->M_040 * d.D_140 + m->M_103 * d.D_203 +
m->M_112 * d.D_212 + m->M_121 * d.D_221 + m->M_130 * d.D_230 +
m->M_202 * d.D_302 + m->M_211 * d.D_311 + m->M_220 * d.D_320 +
m->M_301 * d.D_401 + m->M_310 * d.D_410 + m->M_400 * d.D_500;
*f_y += m->M_004 * d.D_014 + m->M_013 * d.D_023 + m->M_022 * d.D_032 +
m->M_031 * d.D_041 + m->M_040 * d.D_050 + m->M_103 * d.D_113 +
m->M_112 * d.D_122 + m->M_121 * d.D_131 + m->M_130 * d.D_140 +
m->M_202 * d.D_212 + m->M_211 * d.D_221 + m->M_220 * d.D_230 +
m->M_301 * d.D_311 + m->M_310 * d.D_320 + m->M_400 * d.D_410;
*f_z += m->M_004 * d.D_005 + m->M_013 * d.D_014 + m->M_022 * d.D_023 +
m->M_031 * d.D_032 + m->M_040 * d.D_041 + m->M_103 * d.D_104 +
m->M_112 * d.D_113 + m->M_121 * d.D_122 + m->M_130 * d.D_131 +
m->M_202 * d.D_203 + m->M_211 * d.D_212 + m->M_220 * d.D_221 +
m->M_301 * d.D_302 + m->M_310 * d.D_311 + m->M_400 * d.D_401;
*pot += m->M_004 * d.D_004 + m->M_013 * d.D_013 + m->M_022 * d.D_022 +
m->M_031 * d.D_031 + m->M_040 * d.D_040 + m->M_103 * d.D_103 +
m->M_112 * d.D_112 + m->M_121 * d.D_121 + m->M_130 * d.D_130 +
m->M_202 * d.D_202 + m->M_211 * d.D_211 + m->M_220 * d.D_220 +
m->M_301 * d.D_301 + m->M_310 * d.D_310 + m->M_400 * d.D_400;
#endif
/* Take care of the the sign convention */
*f_x *= -1.f;
*f_y *= -1.f;
......
......@@ -150,10 +150,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
float *restrict f_x, float *restrict f_y, float *restrict f_z,
float *restrict pot) {
/* In the case where the order is < 3, then there is only a monopole term left.
/* In the case where the order is < 2, then there is only a monopole term left.
* We can default to the normal P-P interaction with the mass of the multipole
* and its CoM as the "particle" property */
#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
#if SELF_GRAVITY_MULTIPOLE_ORDER < 2
float f_ij, pot_ij;
runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000,
......@@ -179,7 +179,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
*f_z = m->M_000 * d.D_001;
*pot = m->M_000 * d.D_000;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
/* 1st order contributions */
......@@ -191,8 +191,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
/* *pot = m->M_001 * d.D_001 + m->M_010 * d.D_010 + m->M_100 * d.D_100 ; */
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* 2nd order contributions */
*f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
......@@ -205,8 +204,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
m->M_101 * d.D_101 + m->M_110 * d.D_110 + m->M_200 * d.D_200;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
/* 3rd order contributions */
*f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
......@@ -227,7 +225,31 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
m->M_300 * d.D_300;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
/* 4th order contributions */
*f_x += m->M_004 * d.D_104 + m->M_013 * d.D_113 + m->M_022 * d.D_122 +
m->M_031 * d.D_131 + m->M_040 * d.D_140 + m->M_103 * d.D_203 +
m->M_112 * d.D_212 + m->M_121 * d.D_221 + m->M_130 * d.D_230 +
m->M_202 * d.D_302 + m->M_211 * d.D_311 + m->M_220 * d.D_320 +
m->M_301 * d.D_401 + m->M_310 * d.D_410 + m->M_400 * d.D_500;
*f_y += m->M_004 * d.D_014 + m->M_013 * d.D_023 + m->M_022 * d.D_032 +
m->M_031 * d.D_041 + m->M_040 * d.D_050 + m->M_103 * d.D_113 +
m->M_112 * d.D_122 + m->M_121 * d.D_131 + m->M_130 * d.D_140 +
m->M_202 * d.D_212 + m->M_211 * d.D_221 + m->M_220 * d.D_230 +
m->M_301 * d.D_311 + m->M_310 * d.D_320 + m->M_400 * d.D_410;
*f_z += m->M_004 * d.D_005 + m->M_013 * d.D_014 + m->M_022 * d.D_023 +
m->M_031 * d.D_032 + m->M_040 * d.D_041 + m->M_103 * d.D_104 +
m->M_112 * d.D_113 + m->M_121 * d.D_122 + m->M_130 * d.D_131 +
m->M_202 * d.D_203 + m->M_211 * d.D_212 + m->M_220 * d.D_221 +
m->M_301 * d.D_302 + m->M_310 * d.D_311 + m->M_400 * d.D_401;
*pot += m->M_004 * d.D_004 + m->M_013 * d.D_013 + m->M_022 * d.D_022 +
m->M_031 * d.D_031 + m->M_040 * d.D_040 + m->M_103 * d.D_103 +
m->M_112 * d.D_112 + m->M_121 * d.D_121 + m->M_130 * d.D_130 +
m->M_202 * d.D_202 + m->M_211 * d.D_211 + m->M_220 * d.D_220 +
m->M_301 * d.D_301 + m->M_310 * d.D_310 + m->M_400 * d.D_400;
#endif
/* Take care of the the sign convention */
*f_x *= -1.f;
*f_y *= -1.f;
......@@ -263,10 +285,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
const struct multipole *m, float *restrict f_x, float *restrict f_y,
float *restrict f_z, float *restrict pot) {
/* In the case where the order is < 3, then there is only a monopole term left.
/* In the case where the order is < 2, then there is only a monopole term left.
* We can default to the normal P-P interaction with the mass of the multipole
* and its CoM as the "particle" property */
#if SELF_GRAVITY_MULTIPOLE_ORDER < 3
#if SELF_GRAVITY_MULTIPOLE_ORDER < 2
float f_ij, pot_ij;
runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv,
......@@ -292,7 +314,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
*f_z = m->M_000 * d.D_001;
*pot = m->M_000 * d.D_000;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
/* 1st order contributions */
......@@ -304,8 +326,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
/* *pot = m->M_001 * d.D_001 + m->M_010 * d.D_010 + m->M_100 * d.D_100 ; */
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* 2nd order contributions */
*f_x += m->M_002 * d.D_102 + m->M_011 * d.D_111 + m->M_020 * d.D_120 +
......@@ -318,8 +339,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
m->M_101 * d.D_101 + m->M_110 * d.D_110 + m->M_200 * d.D_200;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
/* 3rd order contributions */
*f_x += m->M_003 * d.D_103 + m->M_012 * d.D_112 + m->M_021 * d.D_121 +
......@@ -340,7 +360,31 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
m->M_300 * d.D_300;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
/* 4th order contributions */
*f_x += m->M_004 * d.D_104 + m->M_013 * d.D_113 + m->M_022 * d.D_122 +
m->M_031 * d.D_131 + m->M_040 * d.D_140 + m->M_103 * d.D_203 +
m->M_112 * d.D_212 + m->M_121 * d.D_221 + m->M_130 * d.D_230 +
m->M_202 * d.D_302 + m->M_211 * d.D_311 + m->M_220 * d.D_320 +
m->M_301 * d.D_401 + m->M_310 * d.D_410 + m->M_400 * d.D_500;
*f_y += m->M_004 * d.D_014 + m->M_013 * d.D_023 + m->M_022 * d.D_032 +
m->M_031 * d.D_041 + m->M_040 * d.D_050 + m->M_103 * d.D_113 +
m->M_112 * d.D_122 + m->M_121 * d.D_131 + m->M_130 * d.D_140 +
m->M_202 * d.D_212 + m->M_211 * d.D_221 + m->M_220 * d.D_230 +
m->M_301 * d.D_311 + m->M_310 * d.D_320 + m->M_400 * d.D_410;
*f_z += m->M_004 * d.D_005 + m->M_013 * d.D_014 + m->M_022 * d.D_023 +
m->M_031 * d.D_032 + m->M_040 * d.D_041 + m->M_103 * d.D_104 +
m->M_112 * d.D_113 + m->M_121 * d.D_122 + m->M_130 * d.D_131 +
m->M_202 * d.D_203 + m->M_211 * d.D_212 + m->M_220 * d.D_221 +
m->M_301 * d.D_302 + m->M_310 * d.D_311 + m->M_400 * d.D_401;
*pot += m->M_004 * d.D_004 + m->M_013 * d.D_013 + m->M_022 * d.D_022 +
m->M_031 * d.D_031 + m->M_040 * d.D_040 + m->M_103 * d.D_103 +
m->M_112 * d.D_112 + m->M_121 * d.D_121 + m->M_130 * d.D_130 +
m->M_202 * d.D_202 + m->M_211 * d.D_211 + m->M_220 * d.D_220 +
m->M_301 * d.D_301 + m->M_310 * d.D_310 + m->M_400 * d.D_400;
#endif
/* Take care of the the sign convention */
*f_x *= -1.f;
*f_y *= -1.f;
......
......@@ -102,9 +102,13 @@ struct potential_derivatives_M2P {
/* 1st order terms */
float D_100, D_010, D_001;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
/* 2nd order terms */
float D_200, D_020, D_002;
float D_110, D_101, D_011;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* 3rd order terms */
float D_300, D_030, D_003;
......@@ -112,8 +116,8 @@ struct potential_derivatives_M2P {
float D_120, D_021;
float D_102, D_012;
float D_111;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
/* 4th order terms */
float D_400, D_040, D_004;
......@@ -123,6 +127,17 @@ struct potential_derivatives_M2P {
float D_220, D_202, D_022;
float D_211, D_121, D_112;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
/* 5th order terms */
float D_005, D_014, D_023;
float D_032, D_041, D_050;
float D_104, D_113, D_122;
float D_131, D_140, D_203;
float D_212, D_221, D_230;
float D_302, D_311, D_320;
float D_401, D_410, D_500;
#endif
};
/**
......@@ -415,6 +430,9 @@ potential_derivatives_compute_M2L(const float r_x, const float r_y,
* @brief Compute all the relevent derivatives of the softened and truncated
* gravitational potential for the M2P kernel.
*
* For M2P, we compute the derivatives to one order higher than
* SELF_GRAVITY_MULTIPOLE_ORDER, as these are needed for the accelerations.
*
* @param r_x x-component of distance vector
* @param r_y y-component of distance vector
* @param r_z z-component of distance vector
......@@ -437,71 +455,106 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
error("Computing M2P derivatives below softening length");
#endif
float Dt_1 = 0.f;
float Dt_3 = 0.f;
float Dt_5 = 0.f;
float Dt_7 = 0.f;
float Dt_1;
float Dt_3;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
float Dt_5;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
float Dt_7;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
float Dt_9;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
float Dt_9 = 0.f;
float Dt_11;
#endif
/* Un-truncated case (Newtonian potential) */
if (!periodic) {
const float r_inv2 = r_inv * r_inv;
Dt_1 = r_inv;
Dt_3 = -1.f * Dt_1 * r_inv2; /* -1 / r^3 */
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
Dt_5 = -3.f * Dt_3 * r_inv2; /* 3 / r^5 */
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
Dt_7 = -5.f * Dt_5 * r_inv2; /* -15 / r^7 */
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
Dt_9 = -7.f * Dt_7 * r_inv2; /* 105 / r^9 */
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
Dt_9 = -7.f * Dt_7 * r_inv2; /* -105 / r^9 */
Dt_11 = -9.f * Dt_9 * r_inv2; /* -945 / r^11 */
#endif
/* Truncated case */
} else if (periodic) {
} else {
/* Get the derivatives of the truncated potential */
const float r = r2 * r_inv;
struct chi_derivatives d;
kernel_long_grav_derivatives(r, r_s_inv, &d);
struct chi_derivatives derivs;
kernel_long_grav_derivatives(r, r_s_inv, &derivs);
const float r_inv2 = r_inv * r_inv;
Dt_1 = d.chi_0 * r_inv;
const float r_inv3 = r_inv2 * r_inv;
Dt_3 = (r * d.chi_1 - d.chi_0) * r_inv3;
Dt_1 = derivs.chi_0 * r_inv;
Dt_3 = (r * derivs.chi_1 - derivs.chi_0) * r_inv3;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
const float r_inv5 = r_inv2 * r_inv3;
Dt_5 = (r * r * d.chi_2 - 3.f * r * d.chi_1 + 3.f * d.chi_0) * r_inv5;
Dt_5 =
(r * r * derivs.chi_2 - 3.f * r * derivs.chi_1 + 3.f * derivs.chi_0) *
r_inv5;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
const float r_inv7 = r_inv2 * r_inv5;
Dt_7 = (r * r * r * d.chi_3 - 6.f * r * r * d.chi_2 + 15.f * r * d.chi_1 -
15.f * d.chi_0) *
Dt_7 = (r * r * r * derivs.chi_3 - 6.f * r * r * derivs.chi_2 +
15.f * r * derivs.chi_1 - 15.f * derivs.chi_0) *
r_inv7;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
const float r_inv9 = r_inv2 * r_inv7;
Dt_9 = (r * r * r * r * d.chi_4 - 10.f * r * r * r * d.chi_3 +
45.f * r * r * d.chi_2 - 105.f * r * d.chi_1 + 105.f * d.chi_0) *
Dt_9 = (r * r * r * r * derivs.chi_4 - 10.f * r * r * r * derivs.chi_3 +
45.f * r * r * derivs.chi_2 - 105.f * r * derivs.chi_1 +
105.f * derivs.chi_0) *
r_inv9;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
const float r_inv11 = r_inv2 * r_inv9;
Dt_11 = (r * r * r * r * r * derivs.chi_5 -
15.f * r * r * r * r * derivs.chi_4 +
105.f * r * r * r * derivs.chi_3 - 420.f * r * r * derivs.chi_2 +
945.f * r * derivs.chi_1 - 945.f * derivs.chi_0) *
r_inv11;
#endif
}
/* Compute some powers of r_x, r_y and r_z */
/* Alright, let's get the full terms */
/* Compute some powers of r_x, r_y and r_z */
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
const float r_x2 = r_x * r_x;
const float r_y2 = r_y * r_y;
const float r_z2 = r_z * r_z;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
const float r_x3 = r_x2 * r_x;
const float r_y3 = r_y2 * r_y;
const float r_z3 = r_z2 * r_z;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
const float r_x4 = r_x3 * r_x;
const float r_y4 = r_y3 * r_y;
const float r_z4 = r_z3 * r_z;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
const float r_x5 = r_x4 * r_x;
const float r_y5 = r_y4 * r_y;
const float r_z5 = r_z4 * r_z;
#endif