diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index d4a95540de17631ad445075d672d03a1236e34e3..0d9d2b8c4af40556634e0cf67e9cab45e511636d 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -21,232 +21,94 @@
#define SWIFT_DEFAULT_GRAVITY_IACT_H
/* Includes. */
-#include "const.h"
#include "kernel_gravity.h"
#include "kernel_long_gravity.h"
-#include "multipole.h"
-#include "vector.h"
/**
- * @brief Gravity forces between particles truncated by the long-range kernel
+ * @brief Computes the intensity of the force at a point generated by a
+ * point-mass.
+ *
+ * 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 of the point-mass.
+ * @param f_ij (return) The force intensity.
*/
-__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
- float r2, const float *dx, struct gpart *gpi, struct gpart *gpj,
- float rlr_inv) {
-
- /* Apply the gravitational acceleration. */
- const float r = sqrtf(r2);
- const float ir = 1.f / r;
- const float mi = gpi->mass;
- const float mj = gpj->mass;
- const float hi = gpi->epsilon;
- const float hj = gpj->epsilon;
- const float u_lr = r * rlr_inv;
- float f_lr, fi, fj, W;
-
-#ifdef SWIFT_DEBUG_CHECKS
- if (r == 0.f) error("Interacting particles with 0 distance");
-#endif
-
- /* Get long-range correction */
- kernel_long_grav_eval(u_lr, &f_lr);
-
- if (r >= hi) {
-
- /* Get Newtonian gravity */
- fi = mj * ir * ir * ir * f_lr;
-
- } else {
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_full(
+ float r2, float h2, float h_inv, float h_inv3, float mass, float *f_ij) {
- const float hi_inv = 1.f / hi;
- const float hi_inv3 = hi_inv * hi_inv * hi_inv;
- const float ui = r * hi_inv;
+ /* Get the inverse distance */
+ const float r_inv = 1.f / sqrtf(r2);
- kernel_grav_eval(ui, &W);
-
- /* Get softened gravity */
- fi = mj * hi_inv3 * W * f_lr;
- }
-
- if (r >= hj) {
+ /* Should we soften ? */
+ if (r2 >= h2) {
/* Get Newtonian gravity */
- fj = mi * ir * ir * ir * f_lr;
+ *f_ij = mass * r_inv * r_inv * r_inv;
} else {
- const float hj_inv = 1.f / hj;
- const float hj_inv3 = hj_inv * hj_inv * hj_inv;
- const float uj = r * hj_inv;
+ const float r = r2 * r_inv;
+ const float ui = r * h_inv;
+ float W_ij;
- kernel_grav_eval(uj, &W);
+ kernel_grav_eval(ui, &W_ij);
/* Get softened gravity */
- fj = mi * hj_inv3 * W * f_lr;
+ *f_ij = mass * h_inv3 * W_ij;
}
-
- const float fidx[3] = {fi * dx[0], fi * dx[1], fi * dx[2]};
- gpi->a_grav[0] -= fidx[0];
- gpi->a_grav[1] -= fidx[1];
- gpi->a_grav[2] -= fidx[2];
-
- const float fjdx[3] = {fj * dx[0], fj * dx[1], fj * dx[2]};
- gpj->a_grav[0] += fjdx[0];
- gpj->a_grav[1] += fjdx[1];
- gpj->a_grav[2] += fjdx[2];
}
/**
- * @brief Gravity forces between particles
+ * @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 of the point-mass.
+ * @param rlr_inv Inverse of the mesh smoothing scale.
+ * @param f_ij (return) The force intensity.
*/
-__attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
- float r2, const float *dx, struct gpart *gpi, struct gpart *gpj) {
-
- /* Apply the gravitational acceleration. */
- const float r = sqrtf(r2);
- const float ir = 1.f / r;
- const float mi = gpi->mass;
- const float mj = gpj->mass;
- const float hi = gpi->epsilon;
- const float hj = gpj->epsilon;
- float fi, fj, W;
-
-#ifdef SWIFT_DEBUG_CHECKS
- if (r == 0.f) error("Interacting particles with 0 distance");
-#endif
-
- if (r >= hi) {
-
- /* Get Newtonian gravity */
- fi = mj * ir * ir * ir;
-
- } else {
-
- 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);
+__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_truncated(
+ float r2, float h2, float h_inv, float h_inv3, float mass, float rlr_inv,
+ float *f_ij) {
- /* Get softened gravity */
- fi = mj * hi_inv3 * W;
- }
+ /* Get the inverse distance */
+ const float r_inv = 1.f / sqrtf(r2);
+ const float r = r2 * r_inv;
- if (r >= hj) {
+ /* Should we soften ? */
+ if (r2 >= h2) {
/* Get Newtonian gravity */
- fj = mi * ir * ir * ir;
+ *f_ij = mass * r_inv * r_inv * r_inv;
} else {
- const float hj_inv = 1.f / hj;
- const float hj_inv3 = hj_inv * hj_inv * hj_inv;
- const float uj = r * hj_inv;
+ const float r = r2 * r_inv;
+ const float ui = r * h_inv;
+ float W_ij;
- kernel_grav_eval(uj, &W);
+ kernel_grav_eval(ui, &W_ij);
/* Get softened gravity */
- fj = mi * hj_inv3 * W;
+ *f_ij = mass * h_inv3 * W_ij;
}
- const float fidx[3] = {fi * dx[0], fi * dx[1], fi * dx[2]};
- gpi->a_grav[0] -= fidx[0];
- gpi->a_grav[1] -= fidx[1];
- gpi->a_grav[2] -= fidx[2];
-
- const float fjdx[3] = {fj * dx[0], fj * dx[1], fj * dx[2]};
- gpj->a_grav[0] += fjdx[0];
- gpj->a_grav[1] += fjdx[1];
- gpj->a_grav[2] += fjdx[2];
-}
-
-/**
- * @brief Gravity forces between particles truncated by the long-range kernel
- * (non-symmetric version)
- */
-__attribute__((always_inline)) INLINE static void
-runner_iact_grav_pp_truncated_nonsym(float r2, const float *dx,
- struct gpart *gpi, const struct gpart *gpj,
- float rlr_inv) {
-
- /* 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 u_lr = r * rlr_inv;
- float f_lr, f, W;
-
-#ifdef SWIFT_DEBUG_CHECKS
- if (r == 0.f) error("Interacting particles with 0 distance");
-#endif
-
/* Get long-range correction */
- kernel_long_grav_eval(u_lr, &f_lr);
-
- if (r >= hi) {
-
- /* Get Newtonian gravity */
- f = mj * ir * ir * ir * f_lr;
-
- } else {
-
- 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);
-
- /* Get softened gravity */
- f = mj * hi_inv3 * W * f_lr;
- }
-
- const float fdx[3] = {f * dx[0], f * dx[1], f * dx[2]};
-
- gpi->a_grav[0] -= fdx[0];
- gpi->a_grav[1] -= fdx[1];
- gpi->a_grav[2] -= fdx[2];
-}
-
-/**
- * @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) {
-
- /* 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;
- float f, W;
-
-#ifdef SWIFT_DEBUG_CHECKS
- if (r == 0.f) error("Interacting particles with 0 distance");
-#endif
-
- if (r >= hi) {
-
- /* Get Newtonian gravity */
- f = mj * ir * ir * ir;
-
- } else {
-
- 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);
-
- /* Get softened gravity */
- f = mj * hi_inv3 * W;
- }
-
- const float fdx[3] = {f * dx[0], f * dx[1], f * dx[2]};
-
- gpi->a_grav[0] -= fdx[0];
- gpi->a_grav[1] -= fdx[1];
- gpi->a_grav[2] -= fdx[2];
+ const float u_lr = r * rlr_inv;
+ float corr_lr;
+ kernel_long_grav_eval(u_lr, &corr_lr);
+ *f_ij *= corr_lr;
}
#endif /* SWIFT_DEFAULT_GRAVITY_IACT_H */
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 156d7c6facb09e9652b625758f3b5e4006e4f062..ab71ba20ee08c330ec5799717dabc8a6147b2553 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -247,26 +247,10 @@ void runner_dopair_grav_pp_full(struct runner *r, struct cell *ci,
/* Can we use the multipole in cj ? */
if (gravity_multipole_accept(0., rmax_j, theta_crit2, r2)) {
- /* Get the inverse distance */
- const float r_inv = 1.f / sqrtf(r2);
-
- float f_ij, W_ij;
-
- if (r2 >= h2_i) {
-
- /* Get Newtonian gravity */
- f_ij = multi_mass_j * r_inv * r_inv * r_inv;
-
- } else {
-
- const float r = r2 * r_inv;
- const float ui = r * h_inv_i;
-
- kernel_grav_eval(ui, &W_ij);
-
- /* Get softened gravity */
- f_ij = multi_mass_j * h_inv3_i * W_ij;
- }
+ /* Interact! */
+ float f_ij;
+ runner_iact_grav_pp_full(r2, h2_i, h_inv_i, h_inv3_i, multi_mass_j,
+ &f_ij);
/* Store it back */
ci_cache->a_x[pid] = -f_ij * dx;
@@ -318,26 +302,9 @@ void runner_dopair_grav_pp_full(struct runner *r, struct cell *ci,
error("gpj not drifted to current time");
#endif
- /* Get the inverse distance */
- const float r_inv = 1.f / sqrtf(r2);
-
- float f_ij, W_ij;
-
- if (r2 >= h2_i) {
-
- /* Get Newtonian gravity */
- f_ij = mass_j * r_inv * r_inv * r_inv;
-
- } else {
-
- const float r = r2 * r_inv;
- const float ui = r * h_inv_i;
-
- kernel_grav_eval(ui, &W_ij);
-
- /* Get softened gravity */
- f_ij = mass_j * h_inv3_i * W_ij;
- }
+ /* Interact! */
+ float f_ij;
+ runner_iact_grav_pp_full(r2, h2_i, h_inv_i, h_inv3_i, mass_j, &f_ij);
/* Store it back */
a_x -= f_ij * dx;
@@ -385,26 +352,10 @@ void runner_dopair_grav_pp_full(struct runner *r, struct cell *ci,
/* Can we use the multipole in cj ? */
if (gravity_multipole_accept(0., rmax_i, theta_crit2, r2)) {
- /* Get the inverse distance */
- const float r_inv = 1.f / sqrtf(r2);
-
- float f_ji, W_ji;
-
- if (r2 >= h2_j) {
-
- /* Get Newtonian gravity */
- f_ji = multi_mass_i * r_inv * r_inv * r_inv;
-
- } else {
-
- const float r = r2 * r_inv;
- const float uj = r * h_inv_j;
-
- kernel_grav_eval(uj, &W_ji);
-
- /* Get softened gravity */
- f_ji = multi_mass_i * h_inv3_j * W_ji;
- }
+ /* Interact! */
+ float f_ji;
+ runner_iact_grav_pp_full(r2, h2_j, h_inv_j, h_inv3_j, multi_mass_i,
+ &f_ji);
/* Store it back */
cj_cache->a_x[pjd] = -f_ji * dx;
@@ -456,26 +407,9 @@ void runner_dopair_grav_pp_full(struct runner *r, struct cell *ci,
error("gpi not drifted to current time");
#endif
- /* Get the inverse distance */
- const float r_inv = 1.f / sqrtf(r2);
-
- float f_ji, W_ji;
-
- if (r2 >= h2_j) {
-
- /* Get Newtonian gravity */
- f_ji = mass_i * r_inv * r_inv * r_inv;
-
- } else {
-
- const float r = r2 * r_inv;
- const float uj = r * h_inv_j;
-
- kernel_grav_eval(uj, &W_ji);
-
- /* Get softened gravity */
- f_ji = mass_i * h_inv3_j * W_ji;
- }
+ /* Interact! */
+ float f_ji;
+ runner_iact_grav_pp_full(r2, h2_j, h_inv_j, h_inv3_j, mass_i, &f_ji);
/* Store it back */
a_x -= f_ji * dx;
@@ -518,6 +452,7 @@ void runner_dopair_grav_pp_truncated(struct runner *r, struct cell *ci,
const struct space *s = e->s;
const double cell_width = s->width[0];
const double a_smooth = e->gravity_properties->a_smooth;
+ const float theta_crit2 = e->gravity_properties->theta_crit2;
const double rlr = cell_width * a_smooth;
const float rlr_inv = 1. / rlr;
@@ -556,6 +491,18 @@ void runner_dopair_grav_pp_truncated(struct runner *r, struct cell *ci,
gravity_cache_populate(cj_cache, gparts_j, gcount_j, gcount_padded_j,
loc_mean);
+ /* Recover the multipole info and shift the CoM locations */
+ const float rmax_i = ci->multipole->r_max;
+ const float rmax_j = cj->multipole->r_max;
+ const float multi_mass_i = ci->multipole->m_pole.M_000;
+ const float multi_mass_j = cj->multipole->m_pole.M_000;
+ const float CoM_i[3] = {ci->multipole->CoM[0] - loc_mean[0],
+ ci->multipole->CoM[1] - loc_mean[1],
+ ci->multipole->CoM[2] - loc_mean[2]};
+ const float CoM_j[3] = {cj->multipole->CoM[0] - loc_mean[0],
+ cj->multipole->CoM[1] - loc_mean[1],
+ cj->multipole->CoM[2] - loc_mean[2]};
+
/* Ok... Here we go ! */
if (ci_active) {
@@ -576,6 +523,35 @@ void runner_dopair_grav_pp_truncated(struct runner *r, struct cell *ci,
const float h_inv_i = 1.f / h_i;
const float h_inv3_i = h_inv_i * h_inv_i * h_inv_i;
+ /* Distance to the multipole in cj */
+ const float dx = x_i - CoM_j[0];
+ const float dy = y_i - CoM_j[1];
+ const float dz = z_i - CoM_j[2];
+ const float r2 = dx * dx + dy * dy + dz * dz;
+
+ /* Can we use the multipole in cj ? */
+ if (gravity_multipole_accept(0., rmax_j, theta_crit2, r2)) {
+
+ /* Interact! */
+ float f_ij;
+ runner_iact_grav_pp_truncated(r2, h2_i, h_inv_i, h_inv3_i, multi_mass_j,
+ rlr_inv, &f_ij);
+
+ /* Store it back */
+ ci_cache->a_x[pid] = -f_ij * dx;
+ ci_cache->a_y[pid] = -f_ij * dy;
+ ci_cache->a_z[pid] = -f_ij * dz;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Update the interaction counter */
+ gparts_i[pid].num_interacted += cj->multipole->m_pole.num_gpart;
+#endif
+ /* Done with this particle */
+ continue;
+ }
+
+ /* Ok, as of here we have no choice but directly interact everything */
+
/* Local accumulators for the acceleration */
float a_x = 0.f, a_y = 0.f, a_z = 0.f;
@@ -611,31 +587,10 @@ void runner_dopair_grav_pp_truncated(struct runner *r, struct cell *ci,
error("gpj not drifted to current time");
#endif
- /* Get the inverse distance */
- const float r_inv = 1.f / sqrtf(r2);
- const float r = r2 * r_inv;
-
- float f_ij, W_ij, corr_lr;
-
- if (r2 >= h2_i) {
-
- /* Get Newtonian gravity */
- f_ij = mass_j * r_inv * r_inv * r_inv;
-
- } else {
-
- const float ui = r * h_inv_i;
-
- kernel_grav_eval(ui, &W_ij);
-
- /* Get softened gravity */
- f_ij = mass_j * h_inv3_i * W_ij;
- }
-
- /* Get long-range correction */
- const float u_lr = r * rlr_inv;
- kernel_long_grav_eval(u_lr, &corr_lr);
- f_ij *= corr_lr;
+ /* Interact! */
+ float f_ij;
+ runner_iact_grav_pp_truncated(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
+ rlr_inv, &f_ij);
/* Store it back */
a_x -= f_ij * dx;
@@ -674,6 +629,35 @@ void runner_dopair_grav_pp_truncated(struct runner *r, struct cell *ci,
const float h_inv_j = 1.f / h_j;
const float h_inv3_j = h_inv_j * h_inv_j * h_inv_j;
+ /* Distance to the multipole in ci */
+ const float dx = x_j - CoM_i[0];
+ const float dy = y_j - CoM_i[1];
+ const float dz = z_j - CoM_i[2];
+ const float r2 = dx * dx + dy * dy + dz * dz;
+
+ /* Can we use the multipole in cj ? */
+ if (gravity_multipole_accept(0., rmax_i, theta_crit2, r2)) {
+
+ /* Interact! */
+ float f_ji;
+ runner_iact_grav_pp_truncated(r2, h2_j, h_inv_j, h_inv3_j, multi_mass_i,
+ rlr_inv, &f_ji);
+
+ /* Store it back */
+ cj_cache->a_x[pjd] = -f_ji * dx;
+ cj_cache->a_y[pjd] = -f_ji * dy;
+ cj_cache->a_z[pjd] = -f_ji * dz;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Update the interaction counter */
+ gparts_j[pjd].num_interacted += ci->multipole->m_pole.num_gpart;
+#endif
+ /* Done with this particle */
+ continue;
+ }
+
+ /* Ok, as of here we have no choice but directly interact everything */
+
/* Local accumulators for the acceleration */
float a_x = 0.f, a_y = 0.f, a_z = 0.f;
@@ -709,31 +693,10 @@ void runner_dopair_grav_pp_truncated(struct runner *r, struct cell *ci,
error("gpi not drifted to current time");
#endif
- /* Get the inverse distance */
- const float r_inv = 1.f / sqrtf(r2);
- const float r = r2 * r_inv;
-
- float f_ji, W_ji, corr_lr;
-
- if (r2 >= h2_j) {
-
- /* Get Newtonian gravity */
- f_ji = mass_i * r_inv * r_inv * r_inv;
-
- } else {
-
- const float uj = r * h_inv_j;
-
- kernel_grav_eval(uj, &W_ji);
-
- /* Get softened gravity */
- f_ji = mass_i * h_inv3_j * W_ji;
- }
-
- /* Get long-range correction */
- const float u_lr = r * rlr_inv;
- kernel_long_grav_eval(u_lr, &corr_lr);
- f_ji *= corr_lr;
+ /* Interact! */
+ float f_ji;
+ runner_iact_grav_pp_truncated(r2, h2_j, h_inv_j, h_inv3_j, mass_i,
+ rlr_inv, &f_ji);
/* Store it back */
a_x -= f_ji * dx;
@@ -907,26 +870,9 @@ void runner_doself_grav_pp_full(struct runner *r, struct cell *c) {
error("gpj not drifted to current time");
#endif
- /* Get the inverse distance */
- const float r_inv = 1.f / sqrtf(r2);
-
- float f_ij, W_ij;
-
- if (r2 >= h2_i) {
-
- /* Get Newtonian gravity */
- f_ij = mass_j * r_inv * r_inv * r_inv;
-
- } else {
-
- const float r = r2 * r_inv;
- const float ui = r * h_inv_i;
-
- kernel_grav_eval(ui, &W_ij);
-
- /* Get softened gravity */
- f_ij = mass_j * h_inv3_i * W_ij;
- }
+ /* Interact! */
+ float f_ij;
+ runner_iact_grav_pp_full(r2, h2_i, h_inv_i, h_inv3_i, mass_j, &f_ij);
/* Store it back */
a_x -= f_ij * dx;
@@ -1047,31 +993,10 @@ void runner_doself_grav_pp_truncated(struct runner *r, struct cell *c) {
error("gpj not drifted to current time");
#endif
- /* Get the inverse distance */
- const float r_inv = 1.f / sqrtf(r2);
- const float r = r2 * r_inv;
-
- float f_ij, W_ij, corr_lr;
-
- if (r2 >= h2_i) {
-
- /* Get Newtonian gravity */
- f_ij = mass_j * r_inv * r_inv * r_inv;
-
- } else {
-
- const float ui = r * h_inv_i;
-
- kernel_grav_eval(ui, &W_ij);
-
- /* Get softened gravity */
- f_ij = mass_j * h_inv3_i * W_ij;
- }
-
- /* Get long-range correction */
- const float u_lr = r * rlr_inv;
- kernel_long_grav_eval(u_lr, &corr_lr);
- f_ij *= corr_lr;
+ /* Interact! */
+ float f_ij;
+ runner_iact_grav_pp_truncated(r2, h2_i, h_inv_i, h_inv3_i, mass_j,
+ rlr_inv, &f_ij);
/* Store it back */
a_x -= f_ij * dx;
diff --git a/src/tools.c b/src/tools.c
index 68b23e1a815b510abafbd8f7c6285d3633ba0a63..3f866710c387a858c23ed052be1ad463222601d5 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -391,64 +391,6 @@ void self_all_force(struct runner *r, struct cell *ci) {
}
}
-void pairs_single_grav(double *dim, long long int pid,
- struct gpart *restrict gparts, const struct part *parts,
- int N, int periodic) {
-
- int i, k;
- // int mj, mk;
- // double maxratio = 1.0;
- double r2, dx[3];
- float fdx[3], a[3] = {0.0, 0.0, 0.0}, aabs[3] = {0.0, 0.0, 0.0};
- struct gpart pi, pj;
- // double ih = 12.0/6.25;
-
- /* Find "our" part. */
- for (k = 0; k < N; k++)
- if ((gparts[k].id_or_neg_offset < 0 &&
- parts[-gparts[k].id_or_neg_offset].id == pid) ||
- gparts[k].id_or_neg_offset == pid)
- break;
- if (k == N) error("Part not found.");
- pi = gparts[k];
- pi.a_grav[0] = 0.0f;
- pi.a_grav[1] = 0.0f;
- pi.a_grav[2] = 0.0f;
-
- /* Loop over all particle pairs. */
- for (k = 0; k < N; k++) {
- if (gparts[k].id_or_neg_offset == pi.id_or_neg_offset) continue;
- pj = gparts[k];
- for (i = 0; i < 3; i++) {
- dx[i] = pi.x[i] - pj.x[i];
- if (periodic) {
- if (dx[i] < -dim[i] / 2)
- dx[i] += dim[i];
- else if (dx[i] > dim[i] / 2)
- dx[i] -= dim[i];
- }
- 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);
- a[0] += pi.a_grav[0];
- a[1] += pi.a_grav[1];
- a[2] += pi.a_grav[2];
- aabs[0] += fabsf(pi.a_grav[0]);
- aabs[1] += fabsf(pi.a_grav[1]);
- aabs[2] += fabsf(pi.a_grav[2]);
- pi.a_grav[0] = 0.0f;
- pi.a_grav[1] = 0.0f;
- pi.a_grav[2] = 0.0f;
- }
-
- /* Dump the result. */
- message(
- "acceleration on gpart %lli is a=[ %e %e %e ], |a|=[ %.2e %.2e %.2e ].\n",
- parts[-pi.id_or_neg_offset].id, a[0], a[1], a[2], aabs[0], aabs[1],
- aabs[2]);
-}
-
/**
* @brief Compute the force on a single particle brute-force.
*/
@@ -738,69 +680,3 @@ int compare_particles(struct part a, struct part b, double threshold) {
#endif
}
-
-/**
- * @brief Computes the forces between all g-particles using the N^2 algorithm
- *
- * Overwrites the accelerations of the gparts with the values.
- * Do not use for actual runs.
- *
- * @brief gparts The array of particles.
- * @brief gcount The number of particles.
- * @brief constants Physical constants in internal units.
- * @brief gravity_properties Constants governing the gravity scheme.
- */
-void gravity_n2(struct gpart *gparts, const int gcount,
- const struct phys_const *constants,
- const struct gravity_props *gravity_properties, float rlr) {
-
- const float rlr_inv = 1. / rlr;
- const float r_cut = gravity_properties->r_cut_max;
- const float max_d = r_cut * rlr;
- const float max_d2 = max_d * max_d;
-
- message("rlr_inv= %f", rlr_inv);
- message("max_d: %f", max_d);
-
- /* Reset everything */
- for (int pid = 0; pid < gcount; pid++) {
- struct gpart *restrict gpi = &gparts[pid];
- gpi->a_grav[0] = 0.f;
- gpi->a_grav[1] = 0.f;
- gpi->a_grav[2] = 0.f;
- }
-
- /* Loop over all particles in ci... */
- for (int pid = 0; pid < gcount; pid++) {
-
- /* Get a hold of the ith part in ci. */
- struct gpart *restrict gpi = &gparts[pid];
-
- for (int pjd = pid + 1; pjd < gcount; pjd++) {
-
- /* Get a hold of the jth part in ci. */
- struct gpart *restrict gpj = &gparts[pjd];
-
- /* Compute the pairwise distance. */
- const float dx[3] = {gpi->x[0] - gpj->x[0], // x
- gpi->x[1] - gpj->x[1], // y
- gpi->x[2] - gpj->x[2]}; // z
- const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
-
- if (r2 < max_d2 || 1) {
-
- /* Apply the gravitational acceleration. */
- runner_iact_grav_pp(r2, dx, gpi, gpj);
- }
- }
- }
-
- /* Multiply by Newton's constant */
- const double const_G = constants->const_newton_G;
- for (int pid = 0; pid < gcount; pid++) {
- struct gpart *restrict gpi = &gparts[pid];
- gpi->a_grav[0] *= const_G;
- gpi->a_grav[1] *= const_G;
- gpi->a_grav[2] *= const_G;
- }
-}