From d44155a8bc3cb9c48dcafe22c6a3f0d9aa127e9b Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Mon, 13 Jun 2016 00:04:57 +0100
Subject: [PATCH] Use the long-range smoothing kernel
---
src/const.h | 4 +-
src/gravity/Default/gravity_iact.h | 66 +++++++++++++++++++-----------
src/multipole.c | 24 ++++++-----
src/multipole.h | 12 +++---
src/runner_doiact_grav.h | 17 ++++----
src/tools.c | 4 +-
6 files changed, 76 insertions(+), 51 deletions(-)
diff --git a/src/const.h b/src/const.h
index f7cf611336..ba71e3aeb1 100644
--- a/src/const.h
+++ b/src/const.h
@@ -53,7 +53,9 @@
//#define DEFAULT_SPH
/* Self gravity stuff. */
-#define multipole_order 2
+#define const_gravity_multipole_order 2
+#define const_gravity_a_smooth 1.25f
+#define const_gravity_r_cut 4.5f
#define const_gravity_eta 0.025f
/* External gravity properties */
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index 997bff0fb1..2b3dfc8555 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -31,7 +31,8 @@
* @brief Gravity forces between particles
*/
__attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
- float r2, const float *dx, struct gpart *gpi, struct gpart *gpj) {
+ float rlr_inv, float r2, const float *dx, struct gpart *gpi,
+ struct gpart *gpj) {
/* Apply the gravitational acceleration. */
const float r = sqrtf(r2);
@@ -39,37 +40,45 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
const float mi = gpi->mass;
const float mj = gpj->mass;
const float hi = gpi->epsilon;
- const float hi_inv = 1.f / hi;
- const float hi_inv3 = hi_inv * hi_inv * hi_inv;
const float hj = gpj->epsilon;
- const float hj_inv = 1.f / hj;
- const float hj_inv3 = hj_inv * hj_inv * hj_inv;
- const float ui = r * hi_inv;
- const float uj = r * hj_inv;
- float fi, fj, W;
+ const float u = r * rlr_inv;
+ float f_lr, fi, fj, W;
+
+ /* Get long-range correction */
+ kernel_long_grav_eval(u, &f_lr);
if (r >= hi) {
/* Get Newtonian gravity */
- fi = mj * ir * ir * ir;
+ fi = mj * ir * ir * ir * f_lr;
} else {
- /* Get softened gravity */
+ const float hi_inv = 1.f / hi;
+ const float hi_inv3 = hi_inv * hi_inv * hi_inv;
+ const float ui = r * hi_inv;
+
kernel_grav_eval(ui, &W);
- fi = mj * hi_inv3 * W;
+
+ /* Get softened gravity */
+ fi = mj * hi_inv3 * W * f_lr;
}
if (r >= hj) {
/* Get Newtonian gravity */
- fj = mi * ir * ir * ir;
+ fj = mi * ir * ir * ir * f_lr;
} else {
- /* Get softened gravity */
+ const float hj_inv = 1.f / hj;
+ const float hj_inv3 = hj_inv * hj_inv * hj_inv;
+ const float uj = r * hj_inv;
+
kernel_grav_eval(uj, &W);
- fj = mi * hj_inv3 * W;
+
+ /* Get softened gravity */
+ fj = mi * hj_inv3 * W * f_lr;
}
const float fidx[3] = {fi * dx[0], fi * dx[1], fi * dx[2]};
@@ -89,28 +98,35 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
* @brief Gravity forces between particles (non-symmetric version)
*/
__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_nonsym(
- float r2, const float *dx, struct gpart *gpi, const struct gpart *gpj) {
+ float rlr_inv, float r2, const float *dx, struct gpart *gpi,
+ const struct gpart *gpj) {
/* Apply the gravitational acceleration. */
const float r = sqrtf(r2);
const float ir = 1.f / r;
const float mj = gpj->mass;
const float hi = gpi->epsilon;
- const float hi_inv = 1.f / hi;
- const float hi_inv3 = hi_inv * hi_inv * hi_inv;
- const float ui = r * hi_inv;
- float f, W;
+ const float u = r * rlr_inv;
+ float f_lr, f, W;
+
+ /* Get long-range correction */
+ kernel_long_grav_eval(u, &f_lr);
if (r >= hi) {
/* Get Newtonian gravity */
- f = mj * ir * ir * ir;
+ f = mj * ir * ir * ir * f_lr;
} else {
- /* Get softened gravity */
+ const float hi_inv = 1.f / hi;
+ const float hi_inv3 = hi_inv * hi_inv * hi_inv;
+ const float ui = r * hi_inv;
+
kernel_grav_eval(ui, &W);
- f = mj * hi_inv3 * W;
+
+ /* Get softened gravity */
+ f = mj * hi_inv3 * W * f_lr;
}
const float fdx[3] = {f * dx[0], f * dx[1], f * dx[2]};
@@ -125,7 +141,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_nonsym(
* @brief Gravity forces between particle and multipole
*/
__attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
- float r2, const float *dx, struct gpart *gp,
+ float rlr_inv, float r2, const float *dx, struct gpart *gp,
const struct multipole *multi) {
/* Apply the gravitational acceleration. */
@@ -133,7 +149,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
const float ir = 1.f / r;
const float mrinv3 = multi->mass * ir * ir * ir;
-#if multipole_order < 2
+#if const_gravity_multipole_order < 2
/* 0th and 1st order terms */
gp->a_grav[0] += mrinv3 * dx[0];
@@ -141,7 +157,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
gp->a_grav[2] += mrinv3 * dx[2];
gp->mass_interacted += multi->mass;
-#elif multipole_order == 2
+#elif const_gravity_multipole_order == 2
/* Terms up to 2nd order (quadrupole) */
/* Follows the notation in Bonsai */
diff --git a/src/multipole.c b/src/multipole.c
index dcd4f120d1..f0d0c7084f 100644
--- a/src/multipole.c
+++ b/src/multipole.c
@@ -48,7 +48,7 @@ void multipole_reset(struct multipole *m) {
void multipole_init(struct multipole *m, const struct gpart *gparts,
int gcount) {
-#if multipole_order > 2
+#if const_gravity_multipole_order > 2
#error "Multipoles of order >2 not yet implemented."
#endif
@@ -59,7 +59,7 @@ void multipole_init(struct multipole *m, const struct gpart *gparts,
double mass = 0.0;
double com[3] = {0.0, 0.0, 0.0};
-#if multipole_order >= 2
+#if const_gravity_multipole_order >= 2
double I_xx = 0.0, I_yy = 0.0, I_zz = 0.0;
double I_xy = 0.0, I_xz = 0.0, I_yz = 0.0;
#endif
@@ -73,7 +73,7 @@ void multipole_init(struct multipole *m, const struct gpart *gparts,
com[1] += gparts[k].x[1] * w;
com[2] += gparts[k].x[2] * w;
-#if multipole_order >= 2
+#if const_gravity_multipole_order >= 2
I_xx += gparts[k].x[0] * gparts[k].x[0] * w;
I_yy += gparts[k].x[1] * gparts[k].x[1] * w;
I_zz += gparts[k].x[2] * gparts[k].x[2] * w;
@@ -91,7 +91,7 @@ void multipole_init(struct multipole *m, const struct gpart *gparts,
m->CoM[1] = com[1] * imass;
m->CoM[2] = com[2] * imass;
-#if multipole_order >= 2
+#if const_gravity_multipole_order >= 2
m->I_xx = I_xx - imass * com[0] * com[0];
m->I_yy = I_yy - imass * com[1] * com[1];
m->I_zz = I_zz - imass * com[2] * com[2];
@@ -110,7 +110,7 @@ void multipole_init(struct multipole *m, const struct gpart *gparts,
void multipole_add(struct multipole *ma, const struct multipole *mb) {
-#if multipole_order > 2
+#if const_gravity_multipole_order > 2
#error "Multipoles of order >2 not yet implemented."
#endif
@@ -123,7 +123,7 @@ void multipole_add(struct multipole *ma, const struct multipole *mb) {
const double ma_CoM[3] = {ma->CoM[0], ma->CoM[1], ma->CoM[2]};
const double mb_CoM[3] = {mb->CoM[0], mb->CoM[1], mb->CoM[2]};
-#if multipole_order >= 2
+#if const_gravity_multipole_order >= 2
const double ma_I_xx = (double)ma->I_xx + ma_mass * ma_CoM[0] * ma_CoM[0];
const double ma_I_yy = (double)ma->I_yy + ma_mass * ma_CoM[1] * ma_CoM[1];
const double ma_I_zz = (double)ma->I_zz + ma_mass * ma_CoM[2] * ma_CoM[2];
@@ -148,7 +148,7 @@ void multipole_add(struct multipole *ma, const struct multipole *mb) {
ma->CoM[2] = (ma_CoM[2] * ma_mass + mb_CoM[2] * mb_mass) * M_tot_inv;
/* New quadrupole */
-#if multipole_order >= 2
+#if const_gravity_multipole_order >= 2
ma->I_xx = (ma_I_xx + mb_I_xx) - M_tot * ma->CoM[0] * ma->CoM[0];
ma->I_yy = (ma_I_yy + mb_I_yy) - M_tot * ma->CoM[1] * ma->CoM[1];
ma->I_zz = (ma_I_zz + mb_I_zz) - M_tot * ma->CoM[2] * ma->CoM[2];
@@ -167,7 +167,7 @@ void multipole_add(struct multipole *ma, const struct multipole *mb) {
void multipole_addpart(struct multipole *m, struct gpart *p) {
- /* #if multipole_order == 1 */
+ /* #if const_gravity_multipole_order == 1 */
/* /\* Correct the position. *\/ */
/* float mm = m->coeffs[0], mp = p->mass; */
@@ -179,7 +179,8 @@ void multipole_addpart(struct multipole *m, struct gpart *p) {
/* m->coeffs[0] = mm + mp; */
/* #else */
- /* #error( "Multipoles of order %i not yet implemented." , multipole_order )
+ /* #error( "Multipoles of order %i not yet implemented." ,
+ * const_gravity_multipole_order )
*/
/* #endif */
}
@@ -194,7 +195,7 @@ void multipole_addpart(struct multipole *m, struct gpart *p) {
void multipole_addparts(struct multipole *m, struct gpart *p, int N) {
- /* #if multipole_order == 1 */
+ /* #if const_gravity_multipole_order == 1 */
/* /\* Get the combined mass and positions. *\/ */
/* double xp[3] = {0.0, 0.0, 0.0}; */
@@ -216,7 +217,8 @@ void multipole_addparts(struct multipole *m, struct gpart *p, int N) {
/* m->coeffs[0] = mm + mp; */
/* #else */
- /* #error( "Multipoles of order %i not yet implemented." , multipole_order )
+ /* #error( "Multipoles of order %i not yet implemented." ,
+ * const_gravity_multipole_order )
*/
/* #endif */
}
diff --git a/src/multipole.h b/src/multipole.h
index e9fa348b35..dc1f914dd7 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -38,7 +38,7 @@ struct multipole {
/* Multipole mass */
float mass;
-#if multipole_order >= 2
+#if const_gravity_multipole_order >= 2
/* Quadrupole terms */
float I_xx, I_yy, I_zz;
float I_xy, I_xz, I_yz;
@@ -85,14 +85,15 @@ __attribute__((always_inline)) INLINE static void multipole_iact_mm(
/* acc *= const_G * ir * ir * ir; */
/* /\* Compute the forces on both multipoles. *\/ */
- /* #if multipole_order == 1 */
+ /* #if const_gravity_multipole_order == 1 */
/* float mma = ma->coeffs[0], mmb = mb->coeffs[0]; */
/* for (k = 0; k < 3; k++) { */
/* ma->a[k] -= dx[k] * acc * mmb; */
/* mb->a[k] += dx[k] * acc * mma; */
/* } */
/* #else */
- /* #error( "Multipoles of order %i not yet implemented." , multipole_order )
+ /* #error( "Multipoles of order %i not yet implemented." ,
+ * const_gravity_multipole_order )
*/
/* #endif */
}
@@ -127,10 +128,11 @@ __attribute__((always_inline)) INLINE static void multipole_iact_mp(
/* acc *= const_G * ir * ir * ir * m->coeffs[0]; */
/* /\* Compute the forces on both multipoles. *\/ */
- /* #if multipole_order == 1 */
+ /* #if const_gravity_multipole_order == 1 */
/* for (k = 0; k < 3; k++) p->a_grav[k] += dx[k] * acc; */
/* #else */
- /* #error( "Multipoles of order %i not yet implemented." , multipole_order )
+ /* #error( "Multipoles of order %i not yet implemented." ,
+ * const_gravity_multipole_order )
*/
/* #endif */
}
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 86bb47ae66..b8c820ea9f 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -97,6 +97,7 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
struct gpart *restrict gparts = ci->gparts;
const struct multipole multi = cj->multipole;
const int ti_current = e->ti_current;
+ const float rlr_inv = 1.f;
TIMER_TIC;
@@ -137,7 +138,7 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Interact !*/
- runner_iact_grav_pm(r2, dx, gp, &multi);
+ runner_iact_grav_pm(rlr_inv, r2, dx, gp, &multi);
}
TIMER_TOC(TIMER_DOPAIR); // MATTHIEU
@@ -162,6 +163,7 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
struct gpart *restrict gparts_i = ci->gparts;
struct gpart *restrict gparts_j = cj->gparts;
const int ti_current = e->ti_current;
+ const float rlr_inv = 1.f;
TIMER_TIC;
@@ -216,20 +218,20 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
/* Interact ! */
if (gpi->ti_end <= ti_current && gpj->ti_end <= ti_current) {
- runner_iact_grav_pp(r2, dx, gpi, gpj);
+ runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
} else {
if (gpi->ti_end <= ti_current) {
- runner_iact_grav_pp_nonsym(r2, dx, gpi, gpj);
+ runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
} else if (gpj->ti_end <= ti_current) {
dx[0] = -dx[0];
dx[1] = -dx[1];
dx[2] = -dx[2];
- runner_iact_grav_pp_nonsym(r2, dx, gpj, gpi);
+ runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
}
}
}
@@ -253,6 +255,7 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
const int gcount = c->gcount;
struct gpart *restrict gparts = c->gparts;
const int ti_current = e->ti_current;
+ const float rlr_inv = 1.f;
TIMER_TIC;
@@ -297,20 +300,20 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
/* Interact ! */
if (gpi->ti_end <= ti_current && gpj->ti_end <= ti_current) {
- runner_iact_grav_pp(r2, dx, gpi, gpj);
+ runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
} else {
if (gpi->ti_end <= ti_current) {
- runner_iact_grav_pp_nonsym(r2, dx, gpi, gpj);
+ runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
} else if (gpj->ti_end <= ti_current) {
dx[0] = -dx[0];
dx[1] = -dx[1];
dx[2] = -dx[2];
- runner_iact_grav_pp_nonsym(r2, dx, gpj, gpi);
+ runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
}
}
}
diff --git a/src/tools.c b/src/tools.c
index 3af35051cd..743b282288 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -305,7 +305,7 @@ void pairs_single_grav(double *dim, long long int pid,
fdx[i] = dx[i];
}
r2 = fdx[0] * fdx[0] + fdx[1] * fdx[1] + fdx[2] * fdx[2];
- runner_iact_grav_pp(r2, fdx, &pi, &pj);
+ runner_iact_grav_pp(0.f, r2, fdx, &pi, &pj);
a[0] += pi.a_grav[0];
a[1] += pi.a_grav[1];
a[2] += pi.a_grav[2];
@@ -519,7 +519,7 @@ void gravity_n2(struct gpart *gparts, const int gcount,
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Apply the gravitational acceleration. */
- runner_iact_grav_pp(r2, dx, gpi, gpj);
+ runner_iact_grav_pp(0.f, r2, dx, gpi, gpj);
}
}
--
GitLab