Commit 60bca9f9 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Only use derivatives up to order 4 for P2M.

parent 976a13a8
......@@ -4193,8 +4193,8 @@ void engine_skip_force_and_kick(struct engine *e) {
t->type == task_type_timestep || t->subtype == task_subtype_force ||
t->subtype == task_subtype_grav || t->type == task_type_end_force ||
t->type == task_type_grav_long_range ||
t->type == task_type_grav_down ||
t->type == task_type_cooling || t->type == task_type_sourceterms)
t->type == task_type_grav_down || t->type == task_type_cooling ||
t->type == task_type_sourceterms)
t->skip = 1;
}
......
......@@ -168,8 +168,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
const float r_inv = 1.f / sqrtf(r2);
/* Compute the derivatives of the potential */
struct potential_derivatives_M2L d;
compute_potential_derivatives_M2L(r_x, r_y, r_z, r2, r_inv, h, h_inv, 0, 0.f,
struct potential_derivatives_M2P d;
compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 0, 0.f,
&d);
/* 0th order contributions */
......@@ -225,32 +225,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_full(
m->M_120 * d.D_120 + m->M_201 * d.D_201 + m->M_210 * d.D_210 +
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 */
......@@ -298,7 +272,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
*f_x = f_ij * r_x;
*f_y = f_ij * r_y;
*f_z = f_ij * r_z;
*pot =-pot_ij;
*pot = -pot_ij;
#else
......@@ -306,8 +280,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
const float r_inv = 1.f / sqrtf(r2);
/* Compute the derivatives of the potential */
struct potential_derivatives_M2L d;
compute_potential_derivatives_M2L(r_x, r_y, r_z, r2, r_inv, h, h_inv, 1,
struct potential_derivatives_M2P d;
compute_potential_derivatives_M2P(r_x, r_y, r_z, r2, r_inv, h, h_inv, 1,
r_s_inv, &d);
/* 0th order contributions */
......@@ -363,32 +337,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm_truncated(
m->M_120 * d.D_120 + m->M_201 * d.D_201 + m->M_210 * d.D_210 +
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 */
......
......@@ -113,11 +113,16 @@ struct potential_derivatives_M2P {
float D_102, D_012;
float D_111;
/* 4th order terms (terms used in the trace of the octupole interaction) */
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
/* 4th order terms */
float D_400, D_040, D_004;
float D_310, D_301;
float D_130, D_031;
float D_103, D_013;
float D_220, D_202, D_022;
float D_211, D_121, D_112;
#endif
};
/**
......@@ -498,8 +503,8 @@ compute_potential_derivatives_M2P(const float r_x, const float r_y,
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) */
#if SELF_GRAVITY_MULTIPOLE_ORDER > 3
/* 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;
pot->D_004 = r_z4 * Dt_9 + 6.f * r_z2 * Dt_7 + 3.f * Dt_5;
......@@ -509,6 +514,13 @@ compute_potential_derivatives_M2P(const float r_x, const float r_y,
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;
pot->D_220 = r_x2 * r_y2 * Dt_9 + r_x2 * Dt_7 + r_y2 * Dt_7 + Dt_5;
pot->D_202 = r_x2 * r_z2 * Dt_9 + r_x2 * Dt_7 + r_z2 * Dt_7 + Dt_5;
pot->D_022 = r_y2 * r_z2 * Dt_9 + r_y2 * Dt_7 + r_z2 * Dt_7 + Dt_5;
pot->D_211 = r_x2 * r_y * r_z * Dt_9 + r_y * r_z * Dt_7;
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
}
#endif /* SWIFT_GRAVITY_DERIVATIVE_H */
......@@ -1584,7 +1584,7 @@ static INLINE void runner_do_grav_long_range(struct runner *r, struct cell *ci,
/* Call the PM interaction fucntion on the active sub-cells of ci */
//runner_dopair_recursive_grav_pm(r, ci, cj);
// runner_dopair_recursive_grav_pm(r, ci, cj);
runner_dopair_grav_mm(r, ci, cj);
} /* We are in charge of this pair */
} /* Loop over top-level cells */
......
......@@ -48,15 +48,11 @@
/* Task type names. */
const char *taskID_names[task_type_count] = {
"none", "sort", "self",
"pair", "sub_self", "sub_pair",
"init_grav", "ghost_in", "ghost",
"ghost_out", "extra_ghost", "drift_part",
"drift_gpart", "end_force", "kick1",
"kick2", "timestep", "send",
"recv", "grav_long_range","grav_mm",
"grav_down", "grav_mesh", "cooling",
"sourceterms"};
"none", "sort", "self", "pair", "sub_self",
"sub_pair", "init_grav", "ghost_in", "ghost", "ghost_out",
"extra_ghost", "drift_part", "drift_gpart", "end_force", "kick1",
"kick2", "timestep", "send", "recv", "grav_long_range",
"grav_mm", "grav_down", "grav_mesh", "cooling", "sourceterms"};
/* Sub-task type names. */
const char *subtaskID_names[task_subtype_count] = {
......
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