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

Added interactions due to the trace of the octupole in M2P

parent 5e01d104
......@@ -133,69 +133,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
*pot_ij *= corr_pot_lr;
}
/**
* @brief Computes the intensity of the force at a point generated by a
* point-mass truncated for long-distance periodicity.
*
* The returned quantity needs to be multiplied by the distance vector to obtain
* the force vector.
*
* @param r2 Square of the distance to the point-mass.
* @param h2 Square of the softening length.
* @param h_inv Inverse of the softening length.
* @param h_inv3 Cube of the inverse of the softening length.
* @param mass Mass of the point-mass.
* @param rlr_inv Inverse of the mesh smoothing scale.
* @param f_ij (return) The force intensity.
* @param pot_ij (return) The potential.
*/
__attribute__((always_inline)) INLINE static void
runner_iact_grav_pp_truncated_debug(float r2, float h2, float h_inv,
float h_inv3, float mass, float rlr_inv,
float *f_ij, float *f1_ij, float *corr,
float *pot_ij) {
/* *f_ij = 0.f; */
/* *pot_ij = 0.f; */
/* return; */
/* Get the inverse distance */
const float r_inv = 1.f / sqrtf(r2 + FLT_MIN);
const float r = r2 * r_inv;
/* Should we soften ? */
if (r2 >= h2) {
/* Get Newtonian gravity */
*f_ij = mass * r_inv * r_inv * r_inv;
*pot_ij = -mass * r_inv;
} else {
const float ui = r * h_inv;
float W_f_ij, W_pot_ij;
kernel_grav_force_eval(ui, &W_f_ij);
kernel_grav_pot_eval(ui, &W_pot_ij);
/* Get softened gravity */
*f_ij = mass * h_inv3 * W_f_ij;
*pot_ij = mass * h_inv * W_pot_ij;
}
*f1_ij = *pot_ij;
/* Get long-range correction */
const float u_lr = r * rlr_inv;
float corr_f_lr, corr_pot_lr;
kernel_long_grav_force_eval(u_lr, &corr_f_lr);
kernel_long_grav_pot_eval(u_lr, &corr_pot_lr);
*f_ij *= corr_f_lr;
*pot_ij *= corr_pot_lr;
*corr = corr_pot_lr;
}
/**
* @brief Computes the force at a point generated by a multipole.
*
......@@ -254,6 +191,15 @@ __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;
*pot -= m->M_110 * d.D_110 + m->M_101 * d.D_101 + m->M_011 * d.D_011;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
/* Some 4th order terms (trace of octupole) */
*f_x += m->M_300 * d.D_400 + m->M_030 * d.D_130 + m->M_003 * d.D_103;
*f_y += m->M_300 * d.D_310 + m->M_030 * d.D_040 + m->M_003 * d.D_013;
*f_z += m->M_300 * d.D_301 + m->M_030 * d.D_031 + m->M_003 * d.D_004;
*pot -= m->M_300 * d.D_100 + m->M_030 * d.D_030 + m->M_003 * d.D_003;
#endif
/* Take care of the the sign convention */
*f_x *= -1.f;
*f_y *= -1.f;
......@@ -304,6 +250,15 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
*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;
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
/* Some 4th order terms (trace of octupole) */
*f_x += m->M_300 * d.D_400 + m->M_030 * d.D_130 + m->M_003 * d.D_103;
*f_y += m->M_300 * d.D_310 + m->M_030 * d.D_040 + m->M_003 * d.D_013;
*f_z += m->M_300 * d.D_301 + m->M_030 * d.D_031 + m->M_003 * d.D_004;
*pot -= m->M_300 * d.D_100 + m->M_030 * d.D_030 + m->M_003 * d.D_003;
#endif
/* Take care of the the sign convention */
*f_x *= -1.f;
*f_y *= -1.f;
......
......@@ -112,6 +112,12 @@ struct potential_derivatives_M2P {
float D_120, D_021;
float D_102, D_012;
float D_111;
/* 4th order terms (terms used in the trace of the octupole interaction) */
float D_400, D_040, D_004;
float D_310, D_301;
float D_130, D_031;
float D_103, D_013;
};
/**
......@@ -393,7 +399,8 @@ compute_potential_derivatives_M2P(const float r_x, const float r_y, const float
float Dt_3;
float Dt_5;
float Dt_7;
float Dt_9;
/* Un-softened un-truncated case (Newtonian potential) */
if (!periodic && r2 > eps * eps) {
......@@ -403,6 +410,7 @@ compute_potential_derivatives_M2P(const float r_x, const float r_y, const float
Dt_3 = -1.f * Dt_1 * r_inv2; /* -1 / r^3 */
Dt_5 = -3.f * Dt_3 * r_inv2; /* 3 / r^5 */
Dt_7 = -5.f * Dt_5 * r_inv2; /* -15 / r^7 */
Dt_9 = -7.f * Dt_7 * r_inv2; /* -105 / r^9 */
/* Un-softened truncated case */
} else if (periodic && r2 > eps * eps) {
......@@ -416,11 +424,14 @@ compute_potential_derivatives_M2P(const float r_x, const float r_y, const float
const float r_inv3 = r_inv2 * r_inv;
const float r_inv5 = r_inv2 * r_inv3;
const float r_inv7 = r_inv2 * r_inv5;
const float r_inv9 = r_inv2 * r_inv7;
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;
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) *
r_inv9;
/* Softened case */
} else {
......@@ -432,11 +443,13 @@ compute_potential_derivatives_M2P(const float r_x, const float r_y, const float
const float eps_inv3 = eps_inv * eps_inv2;
const float eps_inv5 = eps_inv3 * eps_inv2;
const float eps_inv7 = eps_inv5 * eps_inv2;
const float eps_inv9 = eps_inv7 * eps_inv2;
Dt_1 = eps_inv * D_soft_1(u, u_inv);
Dt_3 = -eps_inv3 * D_soft_3(u, u_inv);
Dt_5 = eps_inv5 * D_soft_5(u, u_inv);
Dt_7 = -eps_inv7 * D_soft_7(u, u_inv);
Dt_9 = eps_inv9 * D_soft_9(u, u_inv);
}
/* Compute some powers of r_x, r_y and r_z */
......@@ -446,6 +459,9 @@ compute_potential_derivatives_M2P(const float r_x, const float r_y, const float
const float r_x3 = r_x2 * r_x;
const float r_y3 = r_y2 * r_y;
const float r_z3 = r_z2 * r_z;
const float r_x4 = r_x3 * r_x;
const float r_y4 = r_y3 * r_y;
const float r_z4 = r_z3 * r_z;
/* 0th order derivative */
pot->D_000 = Dt_1;
......@@ -474,6 +490,18 @@ compute_potential_derivatives_M2P(const float r_x, const float r_y, const float
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;
/* Selection of 4th order derivatives (the ones used by terms */
/* based on the trace of the octupole) */
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;
pot->D_004 = r_z4 * Dt_9 + 6.f * r_z2 * Dt_7 + 3.f * Dt_5;
pot->D_310 = r_x3 * r_y * Dt_9 + 3.f * r_x * r_y * Dt_7;
pot->D_301 = r_x3 * r_z * Dt_9 + 3.f * r_x * r_z * Dt_7;
pot->D_130 = r_y3 * r_x * Dt_9 + 3.f * r_y * r_x * Dt_7;
pot->D_031 = r_y3 * r_z * Dt_9 + 3.f * r_y * r_z * Dt_7;
pot->D_103 = r_z3 * r_x * Dt_9 + 3.f * r_z * r_x * Dt_7;
pot->D_013 = r_z3 * r_y * Dt_9 + 3.f * r_z * r_y * Dt_7;
}
#endif /* SWIFT_GRAVITY_DERIVATIVE_H */
......@@ -401,10 +401,9 @@ static INLINE void runner_dopair_grav_pp_truncated(
/* long long id_j = e->s->parts[-gparts_j[pjd].id_or_neg_offset].id; */
/* Interact! */
float f_ij, pot_ij, f1_ij, corr;
runner_iact_grav_pp_truncated_debug(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
rlr_inv, &f_ij, &f1_ij, &corr,
&pot_ij);
float f_ij, pot_ij;
runner_iact_grav_pp_truncated(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
rlr_inv, &f_ij, &pot_ij);
/* if (id_i == 1000 && id_j == 900 && pjd < gcount_j) { */
/* message("--- Interacting part"); */
......@@ -1060,10 +1059,9 @@ static INLINE void runner_doself_grav_pp_truncated(struct runner *r,
/* long long id_j = e->s->parts[-gparts[pjd].id_or_neg_offset].id; */
/* Interact! */
float f_ij, pot_ij, f1_ij, corr;
runner_iact_grav_pp_truncated_debug(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
rlr_inv, &f_ij, &f1_ij, &corr,
&pot_ij);
float f_ij, pot_ij;
runner_iact_grav_pp_truncated(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
rlr_inv, &f_ij, &pot_ij);
/* if (id_i == 1000 && id_j == 901) { */
/* message("--- Interacting part"); */
......
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