From 60bca9f9efee4f7cb0acd860b117b38e4a34be9e Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Fri, 22 Jun 2018 14:42:37 +0100
Subject: [PATCH] Only use derivatives up to order 4 for P2M.
---
src/engine.c | 4 +-
src/gravity/Default/gravity_iact.h | 62 +++---------------------------
src/gravity_derivatives.h | 18 +++++++--
src/runner_doiact_grav.h | 2 +-
src/task.c | 14 +++----
5 files changed, 28 insertions(+), 72 deletions(-)
diff --git a/src/engine.c b/src/engine.c
index d9c7212d4a..6fdc5ee88e 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -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;
}
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index 3926f95f7d..4619ce0b12 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -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 */
diff --git a/src/gravity_derivatives.h b/src/gravity_derivatives.h
index 6e6c1505a0..55ea7e6a71 100644
--- a/src/gravity_derivatives.h
+++ b/src/gravity_derivatives.h
@@ -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 */
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index b8a668bf98..240635922b 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.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 */
diff --git a/src/task.c b/src/task.c
index 25e77da89f..ef7c5f0eef 100644
--- a/src/task.c
+++ b/src/task.c
@@ -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] = {
--
GitLab