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

Add quadrupole term in the M2P kernel.

parent 3606f7e1
No related branches found
No related tags found
2 merge requests!566Periodic gravity calculation,!565Mesh force task
...@@ -218,6 +218,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full( ...@@ -218,6 +218,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
float r_x, float r_y, float r_z, float r2, float h, float h_inv, float r_x, float r_y, float r_z, float r2, float h, float h_inv,
const struct multipole *m, float *f_x, float *f_y, float *f_z, float *pot) { 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. */
/* 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 < 3
float f_ij; float f_ij;
runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000, runner_iact_grav_pp_full(r2, h * h, h_inv, h_inv * h_inv * h_inv, m->M_000,
...@@ -230,8 +233,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full( ...@@ -230,8 +233,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
/* Get the inverse distance */ /* Get the inverse distance */
const float r_inv = 1.f / sqrtf(r2); const float r_inv = 1.f / sqrtf(r2);
/* Compute the derivatives of the potential */
struct potential_derivatives_M2P d; struct potential_derivatives_M2P d;
compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, &d); compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 0, FLT_MAX, &d);
/* 1st order terms (monopole) */ /* 1st order terms (monopole) */
*f_x = m->M_000 * d.D_100; *f_x = m->M_000 * d.D_100;
...@@ -250,6 +254,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full( ...@@ -250,6 +254,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
*f_z += m->M_110 * d.D_111 + m->M_101 * d.D_102 + m->M_011 * d.D_012; *f_z += m->M_110 * d.D_111 + m->M_101 * d.D_102 + m->M_011 * d.D_012;
*pot -= m->M_110 * d.D_110 + m->M_101 * d.D_101 + m->M_011 * d.D_011; *pot -= m->M_110 * d.D_110 + m->M_101 * d.D_101 + m->M_011 * d.D_011;
/* Take care of the the sign convention */
*f_x *= -1.f;
*f_y *= -1.f;
*f_z *= -1.f;
*pot *= -1.f;
#endif #endif
} }
...@@ -258,12 +267,49 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated( ...@@ -258,12 +267,49 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
float rlr_inv, const struct multipole *m, float *f_x, float *f_y, float rlr_inv, const struct multipole *m, float *f_x, float *f_y,
float *f_z, float *pot) { float *f_z, float *pot) {
/* 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
float f_ij; float f_ij;
runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv, runner_iact_grav_pp_truncated(r2, h * h, h_inv, h_inv * h_inv * h_inv,
m->M_000, rlr_inv, &f_ij, pot); m->M_000, rlr_inv, &f_ij, pot);
*f_x = f_ij * r_x; *f_x = f_ij * r_x;
*f_y = f_ij * r_y; *f_y = f_ij * r_y;
*f_z = f_ij * r_z; *f_z = f_ij * r_z;
#else
/* Get the inverse distance */
const float r_inv = 1.f / sqrtf(r2);
/* Compute the derivatives of the potential */
struct potential_derivatives_M2P d;
compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 1, rlr_inv, &d);
/* 1st order terms (monopole) */
*f_x = m->M_000 * d.D_100;
*f_y = m->M_000 * d.D_010;
*f_z = m->M_000 * d.D_001;
*pot = m->M_000 * d.D_000;
/* 3rd order terms (quadrupole) */
*f_x += m->M_200 * d.D_300 + m->M_020 * d.D_120 + m->M_002 * d.D_102;
*f_y += m->M_200 * d.D_210 + m->M_020 * d.D_030 + m->M_002 * d.D_012;
*f_z += m->M_200 * d.D_201 + m->M_020 * d.D_021 + m->M_002 * d.D_003;
*pot -= m->M_200 * d.D_100 + m->M_020 * d.D_020 + m->M_002 * d.D_002;
*f_x += m->M_110 * d.D_210 + m->M_101 * d.D_201 + m->M_011 * d.D_111;
*f_y += m->M_110 * d.D_120 + m->M_101 * d.D_111 + m->M_011 * d.D_021;
*f_z += m->M_110 * d.D_111 + m->M_101 * d.D_102 + m->M_011 * d.D_012;
*pot -= m->M_110 * d.D_110 + m->M_101 * d.D_101 + m->M_011 * d.D_011;
/* Take care of the the sign convention */
*f_x *= -1.f;
*f_y *= -1.f;
*f_z *= -1.f;
*pot *= -1.f;
#endif
} }
#endif /* SWIFT_DEFAULT_GRAVITY_IACT_H */ #endif /* SWIFT_DEFAULT_GRAVITY_IACT_H */
...@@ -96,7 +96,7 @@ struct potential_derivatives_M2L { ...@@ -96,7 +96,7 @@ struct potential_derivatives_M2L {
*/ */
struct potential_derivatives_M2P { struct potential_derivatives_M2P {
/* 0th order terms */ /* 0th order term */
float D_000; float D_000;
/* 1st order terms */ /* 1st order terms */
...@@ -130,9 +130,9 @@ struct potential_derivatives_M2P { ...@@ -130,9 +130,9 @@ struct potential_derivatives_M2P {
* @param pot (return) The structure containing all the derivatives. * @param pot (return) The structure containing all the derivatives.
*/ */
__attribute__((always_inline)) INLINE static void __attribute__((always_inline)) INLINE static void
compute_potential_derivatives_M2L(float r_x, float r_y, float r_z, float r2, compute_potential_derivatives_M2L(const float r_x, const float r_y, const float r_z, const float r2,
float r_inv, float eps, float eps_inv, const float r_inv, const float eps, const float eps_inv,
int periodic, float rs_inv, const int periodic, const float rs_inv,
struct potential_derivatives_M2L *pot) { struct potential_derivatives_M2L *pot) {
float Dt_1; float Dt_1;
...@@ -379,11 +379,14 @@ compute_potential_derivatives_M2L(float r_x, float r_y, float r_z, float r2, ...@@ -379,11 +379,14 @@ compute_potential_derivatives_M2L(float r_x, float r_y, float r_z, float r2,
* @param r_inv Inverse norm of distance vector * @param r_inv Inverse norm of distance vector
* @param eps Softening length. * @param eps Softening length.
* @param eps_inv Inverse of softening length. * @param eps_inv Inverse of softening length.
* @param periodic Is the calculation using periodic BCs?
* @param rs_inv The inverse of the gravity mesh-smoothing scale.
* @param pot (return) The structure containing all the derivatives. * @param pot (return) The structure containing all the derivatives.
*/ */
__attribute__((always_inline)) INLINE static void __attribute__((always_inline)) INLINE static void
compute_potential_derivatives_M2P(float r_x, float r_y, float r_z, float r2, compute_potential_derivatives_M2P(const float r_x, const float r_y, const float r_z, const float r2,
float r_inv, float eps, float eps_inv, const float r_inv, const float eps, const float eps_inv,
const int periodic, const float rs_inv,
struct potential_derivatives_M2P *pot) { struct potential_derivatives_M2P *pot) {
float Dt_1; float Dt_1;
...@@ -391,8 +394,8 @@ compute_potential_derivatives_M2P(float r_x, float r_y, float r_z, float r2, ...@@ -391,8 +394,8 @@ compute_potential_derivatives_M2P(float r_x, float r_y, float r_z, float r2,
float Dt_5; float Dt_5;
float Dt_7; float Dt_7;
/* Un-softened case */ /* Un-softened un-truncated case (Newtonian potential) */
if (r2 > eps * eps) { if (!periodic && r2 > eps * eps) {
const float r_inv2 = r_inv * r_inv; const float r_inv2 = r_inv * r_inv;
...@@ -401,6 +404,25 @@ compute_potential_derivatives_M2P(float r_x, float r_y, float r_z, float r2, ...@@ -401,6 +404,25 @@ compute_potential_derivatives_M2P(float r_x, float r_y, float r_z, float r2,
Dt_5 = -3.f * Dt_3 * r_inv2; /* 3 / r^5 */ Dt_5 = -3.f * Dt_3 * r_inv2; /* 3 / r^5 */
Dt_7 = -5.f * Dt_5 * r_inv2; /* -15 / r^7 */ Dt_7 = -5.f * Dt_5 * r_inv2; /* -15 / r^7 */
/* Un-softened truncated case */
} else if (periodic && r2 > eps * eps) {
/* Get the derivatives of the truncated potential */
const float r = r2 * r_inv;
struct truncated_derivatives d;
kernel_long_grav_derivatives(r, rs_inv, &d);
const float r_inv2 = r_inv * r_inv;
const float r_inv3 = r_inv2 * r_inv;
const float r_inv5 = r_inv2 * r_inv3;
const float r_inv7 = r_inv2 * r_inv5;
Dt_1 = d.chi_0 * r_inv;
Dt_3 = (r * d.chi_1 - d.chi_0) * r_inv3;
Dt_5 = (r * r * d.chi_2 - 3.f * r * d.chi_1 + 3.f * d.chi_0) * 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) * r_inv7;
/* Softened case */
} else { } else {
const float r = r2 * r_inv; const float r = r2 * r_inv;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment