From 50c3e8250ef2143e6f175b6cde1dc347449c4aa3 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Wed, 29 Jun 2016 21:42:32 +0100
Subject: [PATCH] Long range forces OK
---
src/engine.c | 5 +++--
src/gravity/Default/gravity.h | 2 +-
src/gravity/Default/gravity_iact.h | 17 +++++++++++------
src/runner_doiact_grav.h | 6 ++++--
src/tools.c | 2 +-
5 files changed, 20 insertions(+), 12 deletions(-)
diff --git a/src/engine.c b/src/engine.c
index 5d8dd1ac16..727f9fa8f8 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1343,10 +1343,11 @@ static inline void engine_make_gravity_dependencies(struct scheduler *sched,
struct cell *c) {
/* init --> gravity --> kick */
- /* grav_up --> gravity ( --> kick) */
scheduler_addunlock(sched, c->super->init, gravity);
- scheduler_addunlock(sched, c->super->grav_up, gravity);
scheduler_addunlock(sched, gravity, c->super->kick);
+
+ /* grav_up --> gravity ( --> kick) */
+ scheduler_addunlock(sched, c->super->grav_up, gravity);
}
/**
diff --git a/src/gravity/Default/gravity.h b/src/gravity/Default/gravity.h
index 7349ad8e22..e3a181cebe 100644
--- a/src/gravity/Default/gravity.h
+++ b/src/gravity/Default/gravity.h
@@ -77,7 +77,7 @@ gravity_compute_timestep_self(const struct phys_const* const phys_const,
*/
__attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
struct gpart* gp) {
- gp->epsilon = 0.1; // MATTHIEU
+ gp->epsilon = 0.; // MATTHIEU
}
/**
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index 2b3dfc8555..048d47e737 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -147,14 +147,19 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
/* Apply the gravitational acceleration. */
const float r = sqrtf(r2);
const float ir = 1.f / r;
+ const float u = r * rlr_inv;
const float mrinv3 = multi->mass * ir * ir * ir;
+ /* Get long-range correction */
+ float f_lr;
+ kernel_long_grav_eval(u, &f_lr);
+
#if const_gravity_multipole_order < 2
/* 0th and 1st order terms */
- gp->a_grav[0] += mrinv3 * dx[0];
- gp->a_grav[1] += mrinv3 * dx[1];
- gp->a_grav[2] += mrinv3 * dx[2];
+ gp->a_grav[0] += mrinv3 * f_lr * dx[0];
+ gp->a_grav[1] += mrinv3 * f_lr * dx[1];
+ gp->a_grav[2] += mrinv3 * f_lr * dx[2];
gp->mass_interacted += multi->mass;
#elif const_gravity_multipole_order == 2
@@ -178,9 +183,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
const float qRR = qRx * dx[0] + qRy * dx[1] + qRz * dx[2];
const float C = D1 + 0.5f * D2 * q + 0.5f * D3 * qRR;
- gp->a_grav[0] -= C * dx[0] + D2 * qRx;
- gp->a_grav[1] -= C * dx[1] + D2 * qRy;
- gp->a_grav[2] -= C * dx[2] + D2 * qRz;
+ gp->a_grav[0] -= f_lr * (C * dx[0] + D2 * qRx);
+ gp->a_grav[1] -= f_lr * (C * dx[1] + D2 * qRy);
+ gp->a_grav[2] -= f_lr * (C * dx[2] + D2 * qRz);
gp->mass_interacted += multi->mass;
#else
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index cf41c735e6..f6bced124c 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -76,6 +76,8 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
TIMER_TIC;
+// message("rlr_inv= %f", rlr_inv);
+
#ifdef SANITY_CHECKS
if (gcount == 0) error("Empty cell!"); // MATTHIEU sanity check
@@ -487,9 +489,9 @@ static void runner_do_grav_mm(struct runner *r, struct cell *ci, int timer) {
cj->loc[2] - pos_i[2]}; // z
const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
- if (!cell_are_neighbours(ci, cj)) runner_dopair_grav_pm(r, ci, cj);
-
if (r2 > max_d2) continue;
+
+ if (!cell_are_neighbours(ci, cj)) runner_dopair_grav_pm(r, ci, cj);
}
}
diff --git a/src/tools.c b/src/tools.c
index f655e5e14c..a60c785104 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -496,7 +496,7 @@ void shuffle_particles(struct part *parts, const int count) {
void gravity_n2(struct gpart *gparts, const int gcount,
const struct phys_const *constants, float rlr) {
- const float rlr_inv = 0.; // 1. / rlr;
+ const float rlr_inv = 1. / rlr;
const float max_d = const_gravity_r_cut * rlr;
const float max_d2 = max_d * max_d;
--
GitLab