Skip to content
Snippets Groups Projects
Commit 002f4be0 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Use the correct += and not = whend adding the 4th-order contribution to M2P.

parent 801b082f
No related branches found
No related tags found
1 merge request!914Improvements to the particle-multipole interactions
Hide whitespace changes
Inline Side-by-side
src/gravity/MultiSoftening/gravity_iact.h
+ 42
− 42
View file @ 002f4be0
......@@ -149,7 +149,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
const float h, const float h_inv, const struct multipole *m, float *f_x,
float *f_y, float *f_z, float *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 < 2
......@@ -227,26 +227,26 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
#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;
*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 */
......@@ -283,7 +283,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
const float h, const float h_inv, const float r_s_inv,
const struct multipole *m, float *f_x, float *f_y, float *f_z, float *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 < 2
......@@ -361,26 +361,26 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
#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;
*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 */
......
src/gravity/Potential/gravity_iact.h
+ 40
− 40
View file @ 002f4be0
......@@ -227,26 +227,26 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
#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;
*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 */
......@@ -361,26 +361,26 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
#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;
*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 */
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment