Commit 4340e825 authored by Stuart McAlpine's avatar Stuart McAlpine Committed by Matthieu Schaller
Browse files

Accelerations from M2P now use derivatives from one order higher than is computed for the multipole

parent abe15cbd
......@@ -150,7 +150,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
/* In the case where the order is < 3, 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,
......@@ -185,8 +185,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 +
......@@ -197,8 +196,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 +
......@@ -214,6 +212,31 @@ __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;
*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 */
......@@ -255,7 +278,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
/* In the case where the order is < 3, 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,
......@@ -290,8 +313,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 +
......@@ -302,8 +324,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 +
......@@ -319,7 +340,31 @@ __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;
*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;
......
......@@ -152,7 +152,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
/* In the case where the order is < 3, 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,
......@@ -178,7 +178,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 */
......@@ -190,8 +190,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 +
......@@ -204,8 +203,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 +
......@@ -226,7 +224,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;
......@@ -264,7 +286,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
/* In the case where the order is < 3, 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,
......@@ -290,7 +312,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 */
......@@ -302,8 +324,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 +
......@@ -316,8 +337,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 +
......@@ -338,7 +358,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;
......
......@@ -415,6 +415,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
......@@ -439,69 +442,104 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
float Dt_1 = 0.f;
float Dt_3 = 0.f;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
float Dt_5 = 0.f;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
float Dt_7 = 0.f;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
float Dt_9 = 0.f;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
float Dt_11 = 0.f;
#endif
/* Un-truncated case (Newtonian potential) */
if (!periodic) {
const float r_inv2 = r_inv * r_inv;
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
/* 0th order derivative */
/* Get the 0th order term */
pot->D_000 = Dt_1;
/* 1st order derivatives */
......@@ -509,6 +547,7 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
pot->D_010 = r_y * Dt_3;
pot->D_001 = r_z * Dt_3;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 0
/* 2nd order derivatives */
pot->D_200 = r_x2 * Dt_5 + Dt_3;
pot->D_020 = r_y2 * Dt_5 + Dt_3;
......@@ -516,7 +555,8 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
pot->D_110 = r_x * r_y * Dt_5;
pot->D_101 = r_x * r_z * Dt_5;
pot->D_011 = r_y * r_z * Dt_5;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 1
/* 3rd order derivatives */
pot->D_300 = r_x3 * Dt_7 + 3.f * r_x * Dt_5;
pot->D_030 = r_y3 * Dt_7 + 3.f * r_y * Dt_5;
......@@ -528,8 +568,8 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
pot->D_102 = r_z2 * r_x * Dt_7 + r_x * Dt_5;
pot->D_012 = r_z2 * r_y * Dt_7 + r_y * Dt_5;
pot->D_111 = r_x * r_y * r_z * Dt_7;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 2
/* 4th order derivatives */
pot->D_400 = r_x4 * Dt_9 + 6.f * r_x2 * Dt_7 + 3.f * Dt_5;
pot->D_040 = r_y4 * Dt_9 + 6.f * r_y2 * Dt_7 + 3.f * Dt_5;
......@@ -547,6 +587,39 @@ potential_derivatives_compute_M2P(const float r_x, const float r_y,
pot->D_121 = r_y2 * r_x * r_z * Dt_9 + r_x * r_z * Dt_7;
pot->D_112 = r_z2 * r_x * r_y * Dt_9 + r_x * r_y * Dt_7;
#endif
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
/* 5th order derivatives */
pot->D_500 = r_x5 * Dt_11 + 10.f * r_x3 * Dt_9 + 15.f * r_x * Dt_7;
pot->D_050 = r_y5 * Dt_11 + 10.f * r_y3 * Dt_9 + 15.f * r_y * Dt_7;
pot->D_005 = r_z5 * Dt_11 + 10.f * r_z3 * Dt_9 + 15.f * r_z * Dt_7;
pot->D_410 = r_x4 * r_y * Dt_11 + 6.f * r_x2 * r_y * Dt_9 + 3.f * r_y * Dt_7;
pot->D_401 = r_x4 * r_z * Dt_11 + 6.f * r_x2 * r_z * Dt_9 + 3.f * r_z * Dt_7;
pot->D_140 = r_y4 * r_x * Dt_11 + 6.f * r_y2 * r_x * Dt_9 + 3.f * r_x * Dt_7;
pot->D_041 = r_y4 * r_z * Dt_11 + 6.f * r_y2 * r_z * Dt_9 + 3.f * r_z * Dt_7;
pot->D_104 = r_z4 * r_x * Dt_11 + 6.f * r_z2 * r_x * Dt_9 + 3.f * r_x * Dt_7;
pot->D_014 = r_z4 * r_y * Dt_11 + 6.f * r_z2 * r_y * Dt_9 + 3.f * r_y * Dt_7;
pot->D_320 = r_x3 * r_y2 * Dt_11 + r_x3 * Dt_9 + 3.f * r_x * r_y2 * Dt_9 +
3.f * r_x * Dt_7;
pot->D_302 = r_x3 * r_z2 * Dt_11 + r_x3 * Dt_9 + 3.f * r_x * r_z2 * Dt_9 +
3.f * r_x * Dt_7;
pot->D_230 = r_y3 * r_x2 * Dt_11 + r_y3 * Dt_9 + 3.f * r_y * r_x2 * Dt_9 +
3.f * r_y * Dt_7;
pot->D_032 = r_y3 * r_z2 * Dt_11 + r_y3 * Dt_9 + 3.f * r_y * r_z2 * Dt_9 +
3.f * r_y * Dt_7;
pot->D_203 = r_z3 * r_x2 * Dt_11 + r_z3 * Dt_9 + 3.f * r_z * r_x2 * Dt_9 +
3.f * r_z * Dt_7;
pot->D_023 = r_z3 * r_y2 * Dt_11 + r_z3 * Dt_9 + 3.f * r_z * r_y2 * Dt_9 +
3.f * r_z * Dt_7;
pot->D_311 = r_x3 * r_y * r_z * Dt_11 + 3.f * r_x * r_y * r_z * Dt_9;
pot->D_131 = r_y3 * r_x * r_z * Dt_11 + 3.f * r_x * r_y * r_z * Dt_9;
pot->D_113 = r_z3 * r_x * r_y * Dt_11 + 3.f * r_x * r_y * r_z * Dt_9;
pot->D_122 = r_x * r_y2 * r_z2 * Dt_11 + r_x * r_y2 * Dt_9 +
r_x * r_z2 * Dt_9 + r_x * Dt_7;
pot->D_212 = r_y * r_x2 * r_z2 * Dt_11 + r_y * r_x2 * Dt_9 +
r_y * r_z2 * Dt_9 + r_y * Dt_7;
pot->D_221 = r_z * r_x2 * r_y2 * Dt_11 + r_z * r_x2 * Dt_9 +
r_z * r_y2 * Dt_9 + r_z * Dt_7;
#endif
}
#endif /* SWIFT_GRAVITY_DERIVATIVE_H */
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment