diff --git a/examples/EAGLE_100/eagle_100.yml b/examples/EAGLE_100/eagle_100.yml
index a9b83b81f085e66b36d115c5265b66d6093ffdfb..1ea1518825debe56cb8462c4a1b398c03c257bfe 100644
--- a/examples/EAGLE_100/eagle_100.yml
+++ b/examples/EAGLE_100/eagle_100.yml
@@ -23,6 +23,12 @@ Snapshots:
Statistics:
delta_time: 1e-2 # Time between statistics output
+# Parameters for the self-gravity scheme
+Gravity:
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ epsilon: 0.0001 # Softening length (in internal units).
+ theta: 0.7 # Opening angle (Multipole acceptance criterion)
+
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
diff --git a/examples/EAGLE_12/eagle_12.yml b/examples/EAGLE_12/eagle_12.yml
index b4697f96d4eed6c3768d492f0b1c52eaed878d5d..fd9569c50f4ad0d87eccd54e34636180c0337596 100644
--- a/examples/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_12/eagle_12.yml
@@ -26,8 +26,8 @@ Statistics:
# Parameters for the self-gravity scheme
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
- theta: 0.7 # Opening angle (Multipole acceptance criterion)
epsilon: 0.0001 # Softening length (in internal units).
+ theta: 0.7 # Opening angle (Multipole acceptance criterion)
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/examples/EAGLE_25/eagle_25.yml b/examples/EAGLE_25/eagle_25.yml
index c755768bcfafebf3efe6307080e9e85d3a0a4bf5..5dee9dad0b5d7f694c61fa4c983ead0f1cd6e5e2 100644
--- a/examples/EAGLE_25/eagle_25.yml
+++ b/examples/EAGLE_25/eagle_25.yml
@@ -27,8 +27,7 @@ Statistics:
Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
epsilon: 0.0001 # Softening length (in internal units).
- a_smooth: 1000.
- r_cut: 4.
+ theta: 0.7 # Opening angle (Multipole acceptance criterion)
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/examples/EAGLE_50/eagle_50.yml b/examples/EAGLE_50/eagle_50.yml
index b84b1eb7c362f85d8cd6a08ff2a15f72d1337396..898c28935abd02ec115ce107bdcfa4006c41dc48 100644
--- a/examples/EAGLE_50/eagle_50.yml
+++ b/examples/EAGLE_50/eagle_50.yml
@@ -23,6 +23,12 @@ Snapshots:
Statistics:
delta_time: 1e-2 # Time between statistics output
+# Parameters for the self-gravity scheme
+Gravity:
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ epsilon: 0.0001 # Softening length (in internal units).
+ theta: 0.7 # Opening angle (Multipole acceptance criterion)
+
# Parameters for the hydrodynamics scheme
SPH:
resolution_eta: 1.2348 # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
diff --git a/examples/UniformDMBox/uniformBox.yml b/examples/UniformDMBox/uniformBox.yml
index cffd442a9a5b16d8e042e41caf9991fcf0e1202e..e59d677b308ca70f212f74c7e4d8b79f015c77a9 100644
--- a/examples/UniformDMBox/uniformBox.yml
+++ b/examples/UniformDMBox/uniformBox.yml
@@ -28,7 +28,7 @@ Gravity:
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion)
epsilon: 0.00001 # Softening length (in internal units).
-
+
# Parameters governing the conserved quantities statistics
Statistics:
delta_time: 1e-2 # Time between statistics output
diff --git a/examples/analyse_tasks.py b/examples/analyse_tasks.py
index d28553173ebb636a99c96082c1eb35ffedca7d6d..970c4a91042b8c61185727f27ef898f93af81fdc 100755
--- a/examples/analyse_tasks.py
+++ b/examples/analyse_tasks.py
@@ -50,8 +50,8 @@ infile = args.input
TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair",
"init_grav", "ghost", "extra_ghost", "drift_part",
"drift_gpart", "kick1", "kick2", "timestep", "send", "recv",
- "grav_top_level", "grav_long_range", "grav_mm", "grav_down",
- "cooling", "sourceterms", "count"]
+ "grav_top_level", "grav_long_range", "grav_ghost", "grav_mm",
+ "grav_down", "cooling", "sourceterms", "count"]
SUBTYPES = ["none", "density", "gradient", "force", "grav", "external_grav",
"tend", "xv", "rho", "gpart", "multipole", "spart", "count"]
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 7ca15576bc9da20a858db6213bb9f79d5a880bf7..9c3cee7630edf1be1e161a3e70547f06e6108ebd 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -51,11 +51,12 @@ SPH:
# Parameters for the self-gravity scheme
Gravity:
- eta: 0.025 # Constant dimensionless multiplier for time integration.
- theta: 0.7 # Opening angle (Multipole acceptance criterion)
- epsilon: 0.1 # Softening length (in internal units).
- a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
- r_cut: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
+ eta: 0.025 # Constant dimensionless multiplier for time integration.
+ theta: 0.7 # Opening angle (Multipole acceptance criterion)
+ epsilon: 0.1 # Softening length (in internal units).
+ a_smooth: 1.25 # (Optional) Smoothing scale in top-level cell sizes to smooth the long-range forces over (this is the default value).
+ r_cut_max: 4.5 # (Optional) Cut-off in number of top-level cells beyond which no FMM forces are computed (this is the default value).
+ r_cut_min: 0.1 # (Optional) Cut-off in number of top-level cells below which no truncation of FMM forces are performed (this is the default value).
# Parameters related to the initial conditions
InitialConditions:
diff --git a/examples/plot_tasks.py b/examples/plot_tasks.py
index e80c3635cf02d409c82960395261179e27cff853..c49020939cca8f744db352631b2ec47267d7bd20 100755
--- a/examples/plot_tasks.py
+++ b/examples/plot_tasks.py
@@ -91,17 +91,18 @@ pl.rcParams.update(PLOT_PARAMS)
TASKTYPES = ["none", "sort", "self", "pair", "sub_self", "sub_pair",
"init_grav", "ghost", "extra_ghost", "drift_part",
"drift_gpart", "kick1", "kick2", "timestep", "send", "recv",
- "grav_top_level", "grav_long_range", "grav_mm", "grav_down",
- "cooling", "sourceterms", "count"]
+ "grav_top_level", "grav_long_range", "grav_ghost", "grav_mm",
+ "grav_down", "cooling", "sourceterms", "count"]
SUBTYPES = ["none", "density", "gradient", "force", "grav", "external_grav",
"tend", "xv", "rho", "gpart", "multipole", "spart", "count"]
# Task/subtypes of interest.
FULLTYPES = ["self/force", "self/density", "self/grav", "sub_self/force",
- "sub_self/density", "pair/force", "pair/density", "pair/grav",
- "sub_pair/force",
- "sub_pair/density", "recv/xv", "send/xv", "recv/rho", "send/rho",
+ "sub_self/density", "sub_self/grav", "pair/force", "pair/density",
+ "pair/grav", "sub_pair/force",
+ "sub_pair/density", "sub_pair/grav", "recv/xv", "send/xv",
+ "recv/rho", "send/rho",
"recv/tend", "send/tend"]
# A number of colours for the various types. Recycled when there are
@@ -109,7 +110,7 @@ FULLTYPES = ["self/force", "self/density", "self/grav", "sub_self/force",
colours = ["cyan", "lightgray", "darkblue", "yellow", "tan", "dodgerblue",
"sienna", "aquamarine", "bisque", "blue", "green", "lightgreen",
"brown", "purple", "moccasin", "olivedrab", "chartreuse",
- "darksage", "darkgreen", "green", "mediumseagreen",
+ "steelblue", "darkgreen", "green", "mediumseagreen",
"mediumaquamarine", "darkslategrey", "mediumturquoise",
"black", "cadetblue", "skyblue", "red", "slategray", "gold",
"slateblue", "blueviolet", "mediumorchid", "firebrick",
diff --git a/src/approx_math.h b/src/approx_math.h
index ad07adeb4f3b1b54ca5f33d80eabb6a004d2a3aa..48319ddfd7a86c132a1cd18b4a08fa849a36a15a 100644
--- a/src/approx_math.h
+++ b/src/approx_math.h
@@ -36,4 +36,17 @@ __attribute__((always_inline)) INLINE static float approx_expf(float x) {
return 1.f + x * (1.f + x * (0.5f + x * (1.f / 6.f + 1.f / 24.f * x)));
}
+/**
+ * @brief Approximate version of expf(x) using a 6th order Taylor expansion
+ *
+ */
+__attribute__((always_inline)) INLINE static float good_approx_expf(float x) {
+ return 1.f +
+ x * (1.f +
+ x * (0.5f +
+ x * ((1.f / 6.f) +
+ x * ((1.f / 24.f) +
+ x * ((1.f / 120.f) + (1.f / 720.f) * x)))));
+}
+
#endif /* SWIFT_APPROX_MATH_H */
diff --git a/src/engine.c b/src/engine.c
index 8c26099204a69a68af121d44776cf5d1ad60a813..b5b2a96c0b7d1fe74c900f25efa3d61d23cefcd6 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1682,6 +1682,16 @@ void engine_exchange_strays(struct engine *e, size_t offset_parts,
#endif
}
+/**
+ * @brief Constructs the top-level tasks for the short-range gravity
+ * and long-range gravity interactions.
+ *
+ * - One FTT task per MPI rank.
+ * - Multiple gravity ghosts for dependencies.
+ * - All top-cells get a self task.
+ * - All pairs within range according to the multipole acceptance
+ * criterion get a pair task.
+ */
void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
void *extra_data) {
@@ -1692,6 +1702,7 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
struct scheduler *sched = &e->sched;
const int nodeID = e->nodeID;
const int periodic = s->periodic;
+ const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
const int cdim_ghost[3] = {s->cdim[0] / 4 + 1, s->cdim[1] / 4 + 1,
s->cdim[2] / 4 + 1};
@@ -1747,8 +1758,9 @@ void engine_make_self_gravity_tasks_mapper(void *map_data, int num_elements,
if (cj->nodeID != nodeID) continue; // MATTHIEU
/* Are the cells to close for a MM interaction ? */
- if (!gravity_multipole_accept(ci->multipole, cj->multipole,
- theta_crit_inv, 1)) {
+ if (!gravity_multipole_accept_rebuild(ci->multipole, cj->multipole,
+ theta_crit_inv, periodic,
+ dim)) {
scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0,
ci, cj);
@@ -1812,6 +1824,11 @@ void engine_make_self_gravity_tasks(struct engine *e) {
if (periodic) free(ghosts);
}
+/**
+ * @brief Constructs the top-level tasks for the external gravity.
+ *
+ * @param e The #engine.
+ */
void engine_make_external_gravity_tasks(struct engine *e) {
struct space *s = e->s;
diff --git a/src/gravity.c b/src/gravity.c
index 49bbaca39b5278009543204a0d9f5e72d69806c4..579546182445ce3fe58c2bbd43c328c0df2d377b 100644
--- a/src/gravity.c
+++ b/src/gravity.c
@@ -63,7 +63,7 @@ int gravity_exact_force_file_exits(const struct engine *e) {
char line[100];
char dummy1[10], dummy2[10];
double epsilon, newton_G;
- int N;
+ int N, periodic;
/* Reads file header */
if (fgets(line, 100, file) != line) error("Problem reading title");
if (fgets(line, 100, file) != line) error("Problem reading G");
@@ -72,10 +72,12 @@ int gravity_exact_force_file_exits(const struct engine *e) {
sscanf(line, "%s %s %d", dummy1, dummy2, &N);
if (fgets(line, 100, file) != line) error("Problem reading epsilon");
sscanf(line, "%s %s %le", dummy1, dummy2, &epsilon);
+ if (fgets(line, 100, file) != line) error("Problem reading BC");
+ sscanf(line, "%s %s %d", dummy1, dummy2, &periodic);
fclose(file);
/* Check whether it matches the current parameters */
- if (N == SWIFT_GRAVITY_FORCE_CHECKS &&
+ if (N == SWIFT_GRAVITY_FORCE_CHECKS && periodic == e->s->periodic &&
(fabs(epsilon - e->gravity_properties->epsilon) / epsilon < 1e-5) &&
(fabs(newton_G - e->physical_constants->const_newton_G) / newton_G <
1e-5)) {
@@ -101,6 +103,8 @@ void gravity_exact_force_compute_mapper(void *map_data, int nr_gparts,
struct exact_force_data *data = (struct exact_force_data *)extra_data;
const struct space *s = data->s;
const struct engine *e = data->e;
+ const int periodic = s->periodic;
+ const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const double const_G = data->const_G;
int counter = 0;
@@ -124,11 +128,18 @@ void gravity_exact_force_compute_mapper(void *map_data, int nr_gparts,
if (gpi == gpj) continue;
/* Compute the pairwise distance. */
- const double dx[3] = {gpi->x[0] - gpj->x[0], // x
- gpi->x[1] - gpj->x[1], // y
- gpi->x[2] - gpj->x[2]}; // z
- const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+ double dx = gpi->x[0] - gpj->x[0];
+ double dy = gpi->x[1] - gpj->x[1];
+ double dz = gpi->x[2] - gpj->x[2];
+
+ /* Now apply periodic BC */
+ if (periodic) {
+ dx = nearest(dx, dim[0]);
+ dy = nearest(dy, dim[1]);
+ dz = nearest(dz, dim[2]);
+ }
+ const double r2 = dx * dx + dy * dy + dz * dz;
const double r = sqrt(r2);
const double ir = 1. / r;
const double mj = gpj->mass;
@@ -152,15 +163,11 @@ void gravity_exact_force_compute_mapper(void *map_data, int nr_gparts,
/* Get softened gravity */
f = mj * hi_inv3 * W * f_lr;
-
- // printf("r=%e hi=%e W=%e fac=%e\n", r, hi, W, f);
}
- const double fdx[3] = {f * dx[0], f * dx[1], f * dx[2]};
-
- a_grav[0] -= fdx[0];
- a_grav[1] -= fdx[1];
- a_grav[2] -= fdx[2];
+ a_grav[0] -= f * dx;
+ a_grav[1] -= f * dy;
+ a_grav[2] -= f * dz;
}
/* Store the exact answer */
@@ -245,8 +252,9 @@ void gravity_exact_force_check(struct space *s, const struct engine *e,
fprintf(file_swift, "# Gravity accuracy test - SWIFT FORCES\n");
fprintf(file_swift, "# G= %16.8e\n", e->physical_constants->const_newton_G);
fprintf(file_swift, "# N= %d\n", SWIFT_GRAVITY_FORCE_CHECKS);
- fprintf(file_swift, "# epsilon=%16.8e\n", e->gravity_properties->epsilon);
- fprintf(file_swift, "# theta=%16.8e\n", e->gravity_properties->theta_crit);
+ fprintf(file_swift, "# epsilon= %16.8e\n", e->gravity_properties->epsilon);
+ fprintf(file_swift, "# periodic= %d\n", s->periodic);
+ fprintf(file_swift, "# theta= %16.8e\n", e->gravity_properties->theta_crit);
fprintf(file_swift, "# Git Branch: %s\n", git_branch());
fprintf(file_swift, "# Git Revision: %s\n", git_revision());
fprintf(file_swift, "# %16s %16s %16s %16s %16s %16s %16s\n", "id", "pos[0]",
diff --git a/src/gravity/Default/gravity_iact.h b/src/gravity/Default/gravity_iact.h
index eca5c2491cbdcf5f0eca01417c8e6b29efc53459..d4a95540de17631ad445075d672d03a1236e34e3 100644
--- a/src/gravity/Default/gravity_iact.h
+++ b/src/gravity/Default/gravity_iact.h
@@ -28,11 +28,11 @@
#include "vector.h"
/**
- * @brief Gravity forces between particles
+ * @brief Gravity forces between particles truncated by the long-range kernel
*/
-__attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
- float rlr_inv, float r2, const float *dx, struct gpart *gpi,
- struct gpart *gpj) {
+__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);
@@ -41,7 +41,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
const float mj = gpj->mass;
const float hi = gpi->epsilon;
const float hj = gpj->epsilon;
- const float u = r * rlr_inv;
+ const float u_lr = r * rlr_inv;
float f_lr, fi, fj, W;
#ifdef SWIFT_DEBUG_CHECKS
@@ -49,7 +49,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
#endif
/* Get long-range correction */
- kernel_long_grav_eval(u, &f_lr);
+ kernel_long_grav_eval(u_lr, &f_lr);
if (r >= hi) {
@@ -97,18 +97,84 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp(
}
/**
- * @brief Gravity forces between particles (non-symmetric version)
+ * @brief Gravity forces between particles
*/
-__attribute__((always_inline)) INLINE static void runner_iact_grav_pp_nonsym(
- float rlr_inv, float r2, const float *dx, struct gpart *gpi,
- const struct gpart *gpj) {
+__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);
+
+ /* Get softened gravity */
+ fi = mj * hi_inv3 * W;
+ }
+
+ if (r >= hj) {
+
+ /* Get Newtonian gravity */
+ fj = mi * ir * ir * ir;
+
+ } else {
+
+ 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);
+
+ /* Get softened gravity */
+ fj = mi * hj_inv3 * W;
+ }
+
+ 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 = r * rlr_inv;
+ const float u_lr = r * rlr_inv;
float f_lr, f, W;
#ifdef SWIFT_DEBUG_CHECKS
@@ -116,7 +182,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_nonsym(
#endif
/* Get long-range correction */
- kernel_long_grav_eval(u, &f_lr);
+ kernel_long_grav_eval(u_lr, &f_lr);
if (r >= hi) {
@@ -143,13 +209,44 @@ __attribute__((always_inline)) INLINE static void runner_iact_grav_pp_nonsym(
}
/**
- * @brief Gravity forces between particle and multipole
+ * @brief Gravity forces between particles (non-symmetric version)
*/
-__attribute__((always_inline)) INLINE static void runner_iact_grav_pm(
- float rlr_inv, float r2, const float *dx, struct gpart *gp,
- const struct multipole *multi) {
+__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;
- error("Dead function");
+#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];
}
#endif /* SWIFT_DEFAULT_GRAVITY_IACT_H */
diff --git a/src/gravity_properties.c b/src/gravity_properties.c
index b1098888b96cdef2205ed513e60a3799c63e8b9f..18cf044434f7840a5a76f483540bb924a2365e26 100644
--- a/src/gravity_properties.c
+++ b/src/gravity_properties.c
@@ -33,7 +33,8 @@
#include "kernel_gravity.h"
#define gravity_props_default_a_smooth 1.25f
-#define gravity_props_default_r_cut 4.5f
+#define gravity_props_default_r_cut_max 4.5f
+#define gravity_props_default_r_cut_min 0.1f
void gravity_props_init(struct gravity_props *p,
const struct swift_params *params) {
@@ -41,8 +42,10 @@ void gravity_props_init(struct gravity_props *p,
/* Tree-PM parameters */
p->a_smooth = parser_get_opt_param_float(params, "Gravity:a_smooth",
gravity_props_default_a_smooth);
- p->r_cut = parser_get_opt_param_float(params, "Gravity:r_cut",
- gravity_props_default_r_cut);
+ p->r_cut_max = parser_get_opt_param_float(params, "Gravity:r_cut_max",
+ gravity_props_default_r_cut_max);
+ p->r_cut_min = parser_get_opt_param_float(params, "Gravity:r_cut_min",
+ gravity_props_default_r_cut_min);
/* Time integration */
p->eta = parser_get_param_float(params, "Gravity:eta");
@@ -69,9 +72,10 @@ void gravity_props_print(const struct gravity_props *p) {
message("Self-gravity softening: epsilon=%.4f (Plummer equivalent: %.4f)",
p->epsilon, p->epsilon / 3.);
- message("Self-gravity MM smoothing-scale: a_smooth=%f", p->a_smooth);
+ message("Self-gravity mesh smoothing-scale: a_smooth=%f", p->a_smooth);
- message("Self-gravity MM cut-off: r_cut=%f", p->r_cut);
+ message("Self-gravity tree cut-off: r_cut_max=%f", p->r_cut_max);
+ message("Self-gravity truncation cut-off: r_cut_min=%f", p->r_cut_min);
}
#if defined(HAVE_HDF5)
@@ -84,7 +88,8 @@ void gravity_props_print_snapshot(hid_t h_grpgrav,
p->epsilon / 3.);
io_write_attribute_f(h_grpgrav, "Opening angle", p->theta_crit);
io_write_attribute_d(h_grpgrav, "MM order", SELF_GRAVITY_MULTIPOLE_ORDER);
- io_write_attribute_f(h_grpgrav, "MM a_smooth", p->a_smooth);
- io_write_attribute_f(h_grpgrav, "MM r_cut", p->r_cut);
+ io_write_attribute_f(h_grpgrav, "Mesh a_smooth", p->a_smooth);
+ io_write_attribute_f(h_grpgrav, "Mesh r_cut_max", p->r_cut_max);
+ io_write_attribute_f(h_grpgrav, "Mesh r_cut_min", p->r_cut_min);
}
#endif
diff --git a/src/gravity_properties.h b/src/gravity_properties.h
index be26f0d1d23b8cec71fa3cbbeedac9f61f337b2c..2a5e4cb1e07ea591e2e3821704ec55abe7980360 100644
--- a/src/gravity_properties.h
+++ b/src/gravity_properties.h
@@ -34,9 +34,16 @@
*/
struct gravity_props {
- /* Tree-PM parameters */
+ /*! Mesh smoothing scale in units of top-level cell size */
float a_smooth;
- float r_cut;
+
+ /*! Distance below which the truncated mesh force is Newtonian in units of
+ * a_smooth */
+ float r_cut_min;
+
+ /*! Distance above which the truncated mesh force is negligible in units of
+ * a_smooth */
+ float r_cut_max;
/*! Time integration dimensionless multiplier */
float eta;
diff --git a/src/kernel_long_gravity.h b/src/kernel_long_gravity.h
index 7b1c5984647c3be232770dc32fc1b112ad8bee94..ec31c2743079da22d1f3dd0c8683adf674aca1e3 100644
--- a/src/kernel_long_gravity.h
+++ b/src/kernel_long_gravity.h
@@ -19,33 +19,67 @@
#ifndef SWIFT_KERNEL_LONG_GRAVITY_H
#define SWIFT_KERNEL_LONG_GRAVITY_H
-#include <math.h>
+/* Config parameters. */
+#include "../config.h"
-/* Includes. */
+/* Local headers. */
+#include "approx_math.h"
#include "const.h"
#include "inline.h"
-#include "vector.h"
-#define one_over_sqrt_pi ((float)(M_2_SQRTPI * 0.5))
+/* Standard headers */
+#include <math.h>
/**
* @brief Computes the long-range correction term for the FFT calculation.
*
- * @param u The ratio of the distance to the FFT cell scale $u = x/A$.
+ * @param u The ratio of the distance to the FFT cell scale \f$u = r/r_s\f$.
* @param W (return) The value of the kernel function.
*/
__attribute__((always_inline)) INLINE static void kernel_long_grav_eval(
float u, float *const W) {
- /* const float arg1 = u * 0.5f; */
- /* const float arg2 = u * one_over_sqrt_pi; */
- /* const float arg3 = -arg1 * arg1; */
+#ifdef GADGET2_LONG_RANGE_CORRECTION
+
+ const float one_over_sqrt_pi = ((float)(M_2_SQRTPI * 0.5));
+
+ const float arg1 = u * 0.5f;
+ const float arg2 = u * one_over_sqrt_pi;
+ const float arg3 = -arg1 * arg1;
+
+ const float term1 = erfcf(arg1);
+ const float term2 = arg2 * expf(arg3);
+
+ *W = term1 + term2;
+#else
+
+ const float arg = 2.f * u;
+ const float exp_arg = good_approx_expf(arg);
+ const float term = 1.f / (1.f + exp_arg);
- /* const float term1 = erfcf(arg1); */
- /* const float term2 = arg2 * expf(arg3); */
+ *W = arg * exp_arg * term * term - exp_arg * term + 1.f;
+ *W *= 2.f;
+#endif
+}
+
+/**
+ * @brief Returns the long-range truncation of the Poisson potential in Fourier
+ * space.
+ *
+ * @param u2 The square of the Fourier mode times the cell scale
+ * \f$u^2 = k^2r_s^2\f$.
+ * @param W (return) The value of the kernel function.
+ */
+__attribute__((always_inline)) INLINE static void fourier_kernel_long_grav_eval(
+ double u2, double *const W) {
- /* *W = term1 + term2; */
- *W = 1.f;
+#ifdef GADGET2_LONG_RANGE_CORRECTION
+ *W = exp(-u2);
+#else
+ const double u = sqrt(u2);
+ const double arg = M_PI_2 * u;
+ *W = arg / sinh(arg);
+#endif
}
#endif // SWIFT_KERNEL_LONG_GRAVITY_H
diff --git a/src/multipole.h b/src/multipole.h
index 23f5194a30b7316aac15073cba36dc404efa21c1..3932b265de2dc7416a96d0a02695fe5973b6fb4e 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -1498,23 +1498,28 @@ INLINE static void gravity_M2M(struct multipole *m_a,
* @param pos_a The position of the multipole.
* @param props The #gravity_props of this calculation.
* @param periodic Is the calculation periodic ?
+ * @param dim The size of the simulation box.
*/
INLINE static void gravity_M2L(struct grav_tensor *l_b,
const struct multipole *m_a,
const double pos_b[3], const double pos_a[3],
- const struct gravity_props *props,
- int periodic) {
+ const struct gravity_props *props, int periodic,
+ const double dim[3]) {
/* Recover some constants */
const double eps2 = props->epsilon2;
/* Compute distance vector */
- const double dx =
- periodic ? box_wrap(pos_b[0] - pos_a[0], 0., 1.) : pos_b[0] - pos_a[0];
- const double dy =
- periodic ? box_wrap(pos_b[1] - pos_a[1], 0., 1.) : pos_b[1] - pos_a[1];
- const double dz =
- periodic ? box_wrap(pos_b[2] - pos_a[2], 0., 1.) : pos_b[2] - pos_a[2];
+ double dx = pos_b[0] - pos_a[0];
+ double dy = pos_b[1] - pos_a[1];
+ double dz = pos_b[2] - pos_a[2];
+
+ /* Apply BC */
+ if (periodic) {
+ dx = nearest(dx, dim[0]);
+ dy = nearest(dy, dim[1]);
+ dz = nearest(dz, dim[2]);
+ }
/* Compute distance */
const double r2 = dx * dx + dy * dy + dz * dz;
@@ -2174,12 +2179,10 @@ INLINE static void gravity_M2L(struct grav_tensor *l_b,
* @param lb The #grav_tensor to shift.
* @param pos_a The position to which m_b will be shifted.
* @param pos_b The current postion of the multipole to shift.
- * @param periodic Is the calculation periodic ?
*/
INLINE static void gravity_L2L(struct grav_tensor *la,
const struct grav_tensor *lb,
- const double pos_a[3], const double pos_b[3],
- int periodic) {
+ const double pos_a[3], const double pos_b[3]) {
/* Initialise everything to zero */
gravity_field_tensors_init(la);
@@ -2636,29 +2639,76 @@ INLINE static void gravity_L2P(const struct grav_tensor *lb,
/**
* @brief Checks whether a cell-cell interaction can be appromixated by a M-M
- * interaction.
+ * interaction using the CoM and cell radius at rebuild.
+ *
+ * We use the multipole acceptance criterion of Dehnen, 2002, JCoPh, Volume 179,
+ * Issue 1, pp.27-42, equation 10.
+ *
+ * @param ma The #multipole of the first #cell.
+ * @param mb The #multipole of the second #cell.
+ * @param theta_crit_inv The inverse of the critical opening angle.
+ * @param periodic Are we using periodic boundary conditions ?
+ * @param dim The dimensions of the box.
+ */
+__attribute__((always_inline)) INLINE static int
+gravity_multipole_accept_rebuild(const struct gravity_tensors *const ma,
+ const struct gravity_tensors *const mb,
+ double theta_crit_inv, int periodic,
+ const double dim[3]) {
+
+ const double r_crit_a = ma->r_max_rebuild * theta_crit_inv;
+ const double r_crit_b = mb->r_max_rebuild * theta_crit_inv;
+
+ double dx = ma->CoM_rebuild[0] - mb->CoM_rebuild[0];
+ double dy = ma->CoM_rebuild[1] - mb->CoM_rebuild[1];
+ double dz = ma->CoM_rebuild[2] - mb->CoM_rebuild[2];
+
+ /* Apply BC */
+ if (periodic) {
+ dx = nearest(dx, dim[0]);
+ dy = nearest(dy, dim[1]);
+ dz = nearest(dz, dim[2]);
+ }
+
+ const double r2 = dx * dx + dy * dy + dz * dz;
+
+ // MATTHIEU: Make this mass-dependent ?
+
+ /* Multipole acceptance criterion (Dehnen 2002, eq.10) */
+ return (r2 > (r_crit_a + r_crit_b) * (r_crit_a + r_crit_b));
+}
+
+/**
+ * @brief Checks whether a cell-cell interaction can be appromixated by a M-M
+ * interaction using the CoM and cell radius at the current time.
+ *
+ * We use the multipole acceptance criterion of Dehnen, 2002, JCoPh, Volume 179,
+ * Issue 1, pp.27-42, equation 10.
*
* @param ma The #multipole of the first #cell.
* @param mb The #multipole of the second #cell.
* @param theta_crit_inv The inverse of the critical opening angle.
- * @param rebuild Are we using the current value of CoM or the ones from
- * the last rebuild ?
+ * @param periodic Are we using periodic boundary conditions ?
+ * @param dim The dimensions of the box.
*/
__attribute__((always_inline)) INLINE static int gravity_multipole_accept(
- const struct gravity_tensors *ma, const struct gravity_tensors *mb,
- double theta_crit_inv, int rebuild) {
-
- const double r_crit_a =
- (rebuild ? ma->r_max_rebuild : ma->r_max) * theta_crit_inv;
- const double r_crit_b =
- (rebuild ? mb->r_max_rebuild : mb->r_max) * theta_crit_inv;
-
- const double dx = rebuild ? ma->CoM_rebuild[0] - mb->CoM_rebuild[0]
- : ma->CoM[0] - mb->CoM[0];
- const double dy = rebuild ? ma->CoM_rebuild[1] - mb->CoM_rebuild[1]
- : ma->CoM[1] - mb->CoM[1];
- const double dz = rebuild ? ma->CoM_rebuild[2] - mb->CoM_rebuild[2]
- : ma->CoM[2] - mb->CoM[2];
+ const struct gravity_tensors *const ma,
+ const struct gravity_tensors *const mb, double theta_crit_inv, int periodic,
+ const double dim[3]) {
+
+ const double r_crit_a = ma->r_max * theta_crit_inv;
+ const double r_crit_b = mb->r_max * theta_crit_inv;
+
+ double dx = ma->CoM[0] - mb->CoM[0];
+ double dy = ma->CoM[1] - mb->CoM[1];
+ double dz = ma->CoM[2] - mb->CoM[2];
+
+ /* Apply BC */
+ if (periodic) {
+ dx = nearest(dx, dim[0]);
+ dy = nearest(dy, dim[1]);
+ dz = nearest(dz, dim[2]);
+ }
const double r2 = dx * dx + dy * dy + dz * dz;
diff --git a/src/runner_doiact_fft.c b/src/runner_doiact_fft.c
index a3e3f38fba920c0c58d600bb25feda88d4a3cf84..26b59f9f6b864445df9190c6041ee684c456ba22 100644
--- a/src/runner_doiact_fft.c
+++ b/src/runner_doiact_fft.c
@@ -20,9 +20,6 @@
/* Config parameters. */
#include "../config.h"
-/* Some standard headers. */
-#include <pthread.h>
-
#ifdef HAVE_FFTW
#include <fftw3.h>
#endif
@@ -33,6 +30,7 @@
/* Local includes. */
#include "engine.h"
#include "error.h"
+#include "kernel_long_gravity.h"
#include "runner.h"
#include "space.h"
#include "timers.h"
@@ -179,11 +177,12 @@ void runner_do_grav_fft(struct runner* r, int timer) {
// error("Top-level multipole %d not drifted", i);
/* Allocates some memory for the density mesh */
- double* restrict rho = fftw_alloc_real(N * N * N);
+ double* restrict rho = fftw_malloc(sizeof(double) * N * N * N);
if (rho == NULL) error("Error allocating memory for density mesh");
/* Allocates some memory for the mesh in Fourier space */
- fftw_complex* restrict frho = fftw_alloc_complex(N * N * (N_half + 1));
+ fftw_complex* restrict frho =
+ fftw_malloc(sizeof(fftw_complex) * N * N * (N_half + 1));
if (frho == NULL)
error("Error allocating memory for transform of density mesh");
@@ -241,7 +240,9 @@ void runner_do_grav_fft(struct runner* r, int timer) {
if (k2 == 0.) continue;
/* Green function */
- const double green_cor = green_fac * exp(-k2 * a_smooth2) / k2;
+ double W;
+ fourier_kernel_long_grav_eval(k2 * a_smooth2, &W);
+ const double green_cor = green_fac * W / k2;
/* Deconvolution of CIC */
const double CIC_cor = sinc_kx_inv * sinc_ky_inv * sinc_kz_inv;
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index a66cc5e0c9ed241aba3bb1b4329016b8e505e280..a5b803ae62ad199dce7168aa9b9fb270d1e7042c 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -36,8 +36,10 @@
*/
void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
+ /* Some constants */
const struct engine *e = r->e;
- const int periodic = e->s->periodic;
+
+ /* Cell properties */
struct gpart *gparts = c->gparts;
const int gcount = c->gcount;
@@ -52,7 +54,6 @@ void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
/* Add the field-tensor to all the 8 progenitors */
for (int k = 0; k < 8; ++k) {
struct cell *cp = c->progeny[k];
- struct grav_tensor temp;
/* Do we have a progenitor with any active g-particles ? */
if (cp != NULL && cell_is_active(cp, e)) {
@@ -61,13 +62,14 @@ void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
if (cp->ti_old_multipole != e->ti_current)
error("cp->multipole not drifted.");
#endif
+ struct grav_tensor shifted_tensor;
/* Shift the field tensor */
- gravity_L2L(&temp, &c->multipole->pot, cp->multipole->CoM,
- c->multipole->CoM, 0 * periodic);
+ gravity_L2L(&shifted_tensor, &c->multipole->pot, cp->multipole->CoM,
+ c->multipole->CoM);
/* Add it to this level's tensor */
- gravity_field_tensors_add(&cp->multipole->pot, &temp);
+ gravity_field_tensors_add(&cp->multipole->pot, &shifted_tensor);
/* Recurse */
runner_do_grav_down(r, cp, 0);
@@ -91,6 +93,7 @@ void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
error("gpart not drifted to current time");
#endif
+ /* Apply the kernel */
gravity_L2P(&c->multipole->pot, c->multipole->CoM, gp);
}
}
@@ -110,10 +113,12 @@ void runner_do_grav_down(struct runner *r, struct cell *c, int timer) {
void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
struct cell *restrict cj) {
+ /* Some constants */
const struct engine *e = r->e;
+ const struct space *s = e->s;
+ const int periodic = s->periodic;
+ const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const struct gravity_props *props = e->gravity_properties;
- const int periodic = e->s->periodic;
- const struct multipole *multi_j = &cj->multipole->m_pole;
// const float a_smooth = e->gravity_properties->a_smooth;
// const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
@@ -122,6 +127,9 @@ void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
/* Anything to do here? */
if (!cell_is_active(ci, e)) return;
+ /* Short-cut to the multipole */
+ const struct multipole *multi_j = &cj->multipole->m_pole;
+
#ifdef SWIFT_DEBUG_CHECKS
if (ci == cj) error("Interacting a cell with itself using M2L");
@@ -136,79 +144,149 @@ void runner_dopair_grav_mm(const struct runner *r, struct cell *restrict ci,
/* Let's interact at this level */
gravity_M2L(&ci->multipole->pot, multi_j, ci->multipole->CoM,
- cj->multipole->CoM, props, periodic * 0);
+ cj->multipole->CoM, props, periodic, dim);
TIMER_TOC(timer_dopair_grav_mm);
}
-/**
- * @brief Computes the interaction of all the particles in a cell with the
- * multipole of another cell.
- *
- * @param r The #runner.
- * @param ci The #cell with particles to interct.
- * @param cj The #cell with the multipole.
- */
-void runner_dopair_grav_pm(const struct runner *r,
- const struct cell *restrict ci,
- const struct cell *restrict cj) {
-
- error("Function should not be called");
-}
-
/**
* @brief Computes the interaction of all the particles in a cell with all the
- * particles of another cell.
+ * particles of another cell using the full Newtonian potential
*
* @param r The #runner.
* @param ci The first #cell.
* @param cj The other #cell.
- *
- * @todo Use a local cache for the particles.
+ * @param shift The distance vector (periodically wrapped) between the cell
+ * centres.
*/
-void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
+void runner_dopair_grav_pp_full(struct runner *r, struct cell *ci,
+ struct cell *cj, double shift[3]) {
- const struct engine *e = r->e;
+ /* Some constants */
+ const struct engine *const e = r->e;
+
+ /* Cell properties */
const int gcount_i = ci->gcount;
const int gcount_j = cj->gcount;
struct gpart *restrict gparts_i = ci->gparts;
struct gpart *restrict gparts_j = cj->gparts;
- const float a_smooth = e->gravity_properties->a_smooth;
- const float rlr_inv = 1. / (a_smooth * ci->super->width[0]);
- TIMER_TIC;
+ /* MATTHIEU: Should we use local DP accumulators ? */
- /* Anything to do here? */
- if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
+ /* Loop over all particles in ci... */
+ if (cell_is_active(ci, e)) {
+ for (int pid = 0; pid < gcount_i; pid++) {
- /* Let's start by drifting things */
- if (!cell_are_gpart_drifted(ci, e)) cell_drift_gpart(ci, e);
- if (!cell_are_gpart_drifted(cj, e)) cell_drift_gpart(cj, e);
+ /* Get a hold of the ith part in ci. */
+ struct gpart *restrict gpi = &gparts_i[pid];
-#if ICHECK > 0
- for (int pid = 0; pid < gcount_i; pid++) {
+ if (!gpart_is_active(gpi, e)) continue;
- /* Get a hold of the ith part in ci. */
- struct gpart *restrict gp = &gparts_i[pid];
+ /* Apply boundary condition */
+ const double pix[3] = {gpi->x[0] - shift[0], gpi->x[1] - shift[1],
+ gpi->x[2] - shift[2]};
- if (gp->id_or_neg_offset == ICHECK)
- message("id=%lld loc=[ %f %f %f ] size= %f count= %d",
- gp->id_or_neg_offset, cj->loc[0], cj->loc[1], cj->loc[2],
- cj->width[0], cj->gcount);
- }
+ /* Loop over every particle in the other cell. */
+ for (int pjd = 0; pjd < gcount_j; pjd++) {
- for (int pid = 0; pid < gcount_j; pid++) {
+ /* Get a hold of the jth part in cj. */
+ const struct gpart *restrict gpj = &gparts_j[pjd];
- /* Get a hold of the ith part in ci. */
- struct gpart *restrict gp = &gparts_j[pid];
+ /* Compute the pairwise distance. */
+ const float dx[3] = {pix[0] - gpj->x[0], // x
+ pix[1] - gpj->x[1], // y
+ pix[2] - gpj->x[2]}; // z
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
- if (gp->id_or_neg_offset == ICHECK)
- message("id=%lld loc=[ %f %f %f ] size= %f count=%d",
- gp->id_or_neg_offset, ci->loc[0], ci->loc[1], ci->loc[2],
- ci->width[0], ci->gcount);
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that particles have been drifted to the current time */
+ if (gpi->ti_drift != e->ti_current)
+ error("gpi not drifted to current time");
+ if (gpj->ti_drift != e->ti_current)
+ error("gpj not drifted to current time");
+#endif
+
+ /* Interact ! */
+ runner_iact_grav_pp_nonsym(r2, dx, gpi, gpj);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gpi->num_interacted++;
+#endif
+ }
+ }
}
+
+ /* Loop over all particles in cj... */
+ if (cell_is_active(cj, e)) {
+ for (int pjd = 0; pjd < gcount_j; pjd++) {
+
+ /* Get a hold of the ith part in ci. */
+ struct gpart *restrict gpj = &gparts_j[pjd];
+
+ if (!gpart_is_active(gpj, e)) continue;
+
+ /* Apply boundary condition */
+ const double pjx[3] = {gpj->x[0] + shift[0], gpj->x[1] + shift[1],
+ gpj->x[2] + shift[2]};
+
+ /* Loop over every particle in the other cell. */
+ for (int pid = 0; pid < gcount_i; pid++) {
+
+ /* Get a hold of the ith part in ci. */
+ const struct gpart *restrict gpi = &gparts_i[pid];
+
+ /* Compute the pairwise distance. */
+ const float dx[3] = {pjx[0] - gpi->x[0], // x
+ pjx[1] - gpi->x[1], // y
+ pjx[2] - gpi->x[2]}; // z
+ const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that particles have been drifted to the current time */
+ if (gpi->ti_drift != e->ti_current)
+ error("gpi not drifted to current time");
+ if (gpj->ti_drift != e->ti_current)
+ error("gpj not drifted to current time");
#endif
+ /* Interact ! */
+ runner_iact_grav_pp_nonsym(r2, dx, gpj, gpi);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gpj->num_interacted++;
+#endif
+ }
+ }
+ }
+}
+
+/**
+ * @brief Computes the interaction of all the particles in a cell with all the
+ * particles of another cell using the truncated Newtonian potential
+ *
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param cj The other #cell.
+ * @param shift The distance vector (periodically wrapped) between the cell
+ * centres.
+ */
+void runner_dopair_grav_pp_truncated(struct runner *r, struct cell *ci,
+ struct cell *cj, double shift[3]) {
+
+ /* Some constants */
+ const struct engine *const e = r->e;
+ const struct space *s = e->s;
+ const double cell_width = s->width[0];
+ const double a_smooth = e->gravity_properties->a_smooth;
+ const double rlr = cell_width * a_smooth;
+ const float rlr_inv = 1. / rlr;
+
+ /* Cell properties */
+ const int gcount_i = ci->gcount;
+ const int gcount_j = cj->gcount;
+ struct gpart *restrict gparts_i = ci->gparts;
+ struct gpart *restrict gparts_j = cj->gparts;
+
/* MATTHIEU: Should we use local DP accumulators ? */
/* Loop over all particles in ci... */
@@ -220,6 +298,10 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
if (!gpart_is_active(gpi, e)) continue;
+ /* Apply boundary condition */
+ const double pix[3] = {gpi->x[0] - shift[0], gpi->x[1] - shift[1],
+ gpi->x[2] - shift[2]};
+
/* Loop over every particle in the other cell. */
for (int pjd = 0; pjd < gcount_j; pjd++) {
@@ -227,9 +309,9 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
const struct gpart *restrict gpj = &gparts_j[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 dx[3] = {pix[0] - gpj->x[0], // x
+ pix[1] - gpj->x[1], // y
+ pix[2] - gpj->x[2]}; // z
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
#ifdef SWIFT_DEBUG_CHECKS
@@ -241,7 +323,7 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
#endif
/* Interact ! */
- runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
+ runner_iact_grav_pp_truncated_nonsym(r2, dx, gpi, gpj, rlr_inv);
#ifdef SWIFT_DEBUG_CHECKS
gpi->num_interacted++;
@@ -259,6 +341,10 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
if (!gpart_is_active(gpj, e)) continue;
+ /* Apply boundary condition */
+ const double pjx[3] = {gpj->x[0] + shift[0], gpj->x[1] + shift[1],
+ gpj->x[2] + shift[2]};
+
/* Loop over every particle in the other cell. */
for (int pid = 0; pid < gcount_i; pid++) {
@@ -266,9 +352,9 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
const struct gpart *restrict gpi = &gparts_i[pid];
/* Compute the pairwise distance. */
- const float dx[3] = {gpj->x[0] - gpi->x[0], // x
- gpj->x[1] - gpi->x[1], // y
- gpj->x[2] - gpi->x[2]}; // z
+ const float dx[3] = {pjx[0] - gpi->x[0], // x
+ pjx[1] - gpi->x[1], // y
+ pjx[2] - gpi->x[2]}; // z
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
#ifdef SWIFT_DEBUG_CHECKS
@@ -280,7 +366,7 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
#endif
/* Interact ! */
- runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
+ runner_iact_grav_pp_truncated_nonsym(r2, dx, gpj, gpi, rlr_inv);
#ifdef SWIFT_DEBUG_CHECKS
gpj->num_interacted++;
@@ -288,51 +374,168 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
}
}
}
+}
+
+/**
+ * @brief Computes the interaction of all the particles in a cell with all the
+ * particles of another cell (switching function between full and truncated).
+ *
+ * @param r The #runner.
+ * @param ci The first #cell.
+ * @param cj The other #cell.
+ */
+void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
+
+ /* Some properties of the space */
+ const struct engine *e = r->e;
+ const struct space *s = e->s;
+ const int periodic = s->periodic;
+ const double cell_width = s->width[0];
+ const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
+ const double a_smooth = e->gravity_properties->a_smooth;
+ const double r_cut_min = e->gravity_properties->r_cut_min;
+ const double min_trunc = cell_width * r_cut_min * a_smooth;
+ double shift[3] = {0.0, 0.0, 0.0};
+
+ TIMER_TIC;
+
+ /* Anything to do here? */
+ if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
+
+ /* Let's start by drifting things */
+ if (!cell_are_gpart_drifted(ci, e)) cell_drift_gpart(ci, e);
+ if (!cell_are_gpart_drifted(cj, e)) cell_drift_gpart(cj, e);
+
+ /* Can we use the Newtonian version or do we need the truncated one ? */
+ if (!periodic) {
+ runner_dopair_grav_pp_full(r, ci, cj, shift);
+ } else {
+
+ /* Get the relative distance between the pairs, wrapping. */
+ shift[0] = nearest(cj->loc[0] - ci->loc[0], dim[0]);
+ shift[1] = nearest(cj->loc[1] - ci->loc[1], dim[1]);
+ shift[2] = nearest(cj->loc[2] - ci->loc[2], dim[2]);
+ const double r2 =
+ shift[0] * shift[0] + shift[1] * shift[1] + shift[2] * shift[2];
+
+ /* Get the maximal distance between any two particles */
+ const double max_r = sqrt(r2) + ci->multipole->r_max + cj->multipole->r_max;
+
+ /* Do we need to use the truncated interactions ? */
+ if (max_r > min_trunc)
+ runner_dopair_grav_pp_truncated(r, ci, cj, shift);
+ else
+ runner_dopair_grav_pp_full(r, ci, cj, shift);
+ }
TIMER_TOC(timer_dopair_grav_pp);
}
/**
- * @brief Computes the interaction of all the particles in a cell directly
+ * @brief Computes the interaction of all the particles in a cell using the
+ * full Newtonian potential.
*
* @param r The #runner.
* @param c The #cell.
*
* @todo Use a local cache for the particles.
*/
-void runner_doself_grav_pp(struct runner *r, struct cell *c) {
+void runner_doself_grav_pp_full(struct runner *r, struct cell *c) {
- const struct engine *e = r->e;
+ /* Some constants */
+ const struct engine *const e = r->e;
+
+ /* Cell properties */
const int gcount = c->gcount;
struct gpart *restrict gparts = c->gparts;
- const float a_smooth = e->gravity_properties->a_smooth;
- const float rlr_inv = 1. / (a_smooth * c->super->width[0]);
- TIMER_TIC;
+ /* MATTHIEU: Should we use local DP accumulators ? */
+
+ /* 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];
+
+ /* Loop over every particle in the other cell. */
+ 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. */
+ 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];
#ifdef SWIFT_DEBUG_CHECKS
- if (c->gcount == 0) error("Doing self gravity on an empty cell !");
+ /* Check that particles have been drifted to the current time */
+ if (gpi->ti_drift != e->ti_current)
+ error("gpi not drifted to current time");
+ if (gpj->ti_drift != e->ti_current)
+ error("gpj not drifted to current time");
#endif
- /* Anything to do here? */
- if (!cell_is_active(c, e)) return;
+ /* Interact ! */
+ if (gpart_is_active(gpi, e) && gpart_is_active(gpj, e)) {
- /* Do we need to start by drifting things ? */
- if (!cell_are_gpart_drifted(c, e)) cell_drift_gpart(c, e);
+ runner_iact_grav_pp(r2, dx, gpi, gpj);
-#if ICHECK > 0
- for (int pid = 0; pid < gcount; pid++) {
+#ifdef SWIFT_DEBUG_CHECKS
+ gpi->num_interacted++;
+ gpj->num_interacted++;
+#endif
- /* Get a hold of the ith part in ci. */
- struct gpart *restrict gp = &gparts[pid];
+ } else {
- if (gp->id_or_neg_offset == ICHECK)
- message("id=%lld loc=[ %f %f %f ] size= %f count= %d",
- gp->id_or_neg_offset, c->loc[0], c->loc[1], c->loc[2],
- c->width[0], c->gcount);
- }
+ if (gpart_is_active(gpi, e)) {
+
+ runner_iact_grav_pp_nonsym(r2, dx, gpi, gpj);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gpi->num_interacted++;
#endif
+ } else if (gpart_is_active(gpj, e)) {
+
+ dx[0] = -dx[0];
+ dx[1] = -dx[1];
+ dx[2] = -dx[2];
+ runner_iact_grav_pp_nonsym(r2, dx, gpj, gpi);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ gpj->num_interacted++;
+#endif
+ }
+ }
+ }
+ }
+}
+
+/**
+ * @brief Computes the interaction of all the particles in a cell using the
+ * truncated Newtonian potential.
+ *
+ * @param r The #runner.
+ * @param c The #cell.
+ *
+ * @todo Use a local cache for the particles.
+ */
+void runner_doself_grav_pp_truncated(struct runner *r, struct cell *c) {
+
+ /* Some constants */
+ const struct engine *const e = r->e;
+ const struct space *s = e->s;
+ const double cell_width = s->width[0];
+ const double a_smooth = e->gravity_properties->a_smooth;
+ const double rlr = cell_width * a_smooth;
+ const float rlr_inv = 1. / rlr;
+
+ /* Cell properties */
+ const int gcount = c->gcount;
+ struct gpart *restrict gparts = c->gparts;
+
/* MATTHIEU: Should we use local DP accumulators ? */
/* Loop over all particles in ci... */
@@ -364,7 +567,7 @@ void runner_doself_grav_pp(struct runner *r, struct cell *c) {
/* Interact ! */
if (gpart_is_active(gpi, e) && gpart_is_active(gpj, e)) {
- runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
+ runner_iact_grav_pp_truncated(r2, dx, gpi, gpj, rlr_inv);
#ifdef SWIFT_DEBUG_CHECKS
gpi->num_interacted++;
@@ -375,7 +578,7 @@ void runner_doself_grav_pp(struct runner *r, struct cell *c) {
if (gpart_is_active(gpi, e)) {
- runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
+ runner_iact_grav_pp_truncated_nonsym(r2, dx, gpi, gpj, rlr_inv);
#ifdef SWIFT_DEBUG_CHECKS
gpi->num_interacted++;
@@ -386,7 +589,7 @@ void runner_doself_grav_pp(struct runner *r, struct cell *c) {
dx[0] = -dx[0];
dx[1] = -dx[1];
dx[2] = -dx[2];
- runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
+ runner_iact_grav_pp_truncated_nonsym(r2, dx, gpj, gpi, rlr_inv);
#ifdef SWIFT_DEBUG_CHECKS
gpj->num_interacted++;
@@ -395,6 +598,52 @@ void runner_doself_grav_pp(struct runner *r, struct cell *c) {
}
}
}
+}
+
+/**
+ * @brief Computes the interaction of all the particles in a cell directly
+ * (Switching function between truncated and full)
+ *
+ * @param r The #runner.
+ * @param c The #cell.
+ */
+void runner_doself_grav_pp(struct runner *r, struct cell *c) {
+
+ /* Some properties of the space */
+ const struct engine *e = r->e;
+ const struct space *s = e->s;
+ const int periodic = s->periodic;
+ const double cell_width = s->width[0];
+ const double a_smooth = e->gravity_properties->a_smooth;
+ const double r_cut_min = e->gravity_properties->r_cut_min;
+ const double min_trunc = cell_width * r_cut_min * a_smooth;
+
+ TIMER_TIC;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->gcount == 0) error("Doing self gravity on an empty cell !");
+#endif
+
+ /* Anything to do here? */
+ if (!cell_is_active(c, e)) return;
+
+ /* Do we need to start by drifting things ? */
+ if (!cell_are_gpart_drifted(c, e)) cell_drift_gpart(c, e);
+
+ /* Can we use the Newtonian version or do we need the truncated one ? */
+ if (!periodic) {
+ runner_doself_grav_pp_full(r, c);
+ } else {
+
+ /* Get the maximal distance between any two particles */
+ const double max_r = 2 * c->multipole->r_max;
+
+ /* Do we need to use the truncated interactions ? */
+ if (max_r > min_trunc)
+ runner_doself_grav_pp_truncated(r, c);
+ else
+ runner_doself_grav_pp_full(r, c);
+ }
TIMER_TOC(timer_doself_grav_pp);
}
@@ -415,6 +664,9 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
/* Some constants */
const struct engine *e = r->e;
+ const struct space *s = e->s;
+ const int periodic = s->periodic;
+ const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const struct gravity_props *props = e->gravity_properties;
const double theta_crit_inv = props->theta_crit_inv;
@@ -460,11 +712,15 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
}
#endif
+ /* Anything to do here? */
+ if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
+
TIMER_TIC;
/* Can we use M-M interactions ? */
if (gravity_multipole_accept(ci->multipole, cj->multipole, theta_crit_inv,
- 0)) {
+ periodic, dim)) {
+
/* MATTHIEU: make a symmetric M-M interaction function ! */
runner_dopair_grav_mm(r, ci, cj);
runner_dopair_grav_mm(r, cj, ci);
@@ -543,6 +799,9 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
*/
void runner_doself_grav(struct runner *r, struct cell *c, int gettimer) {
+ /* Some constants */
+ const struct engine *e = r->e;
+
#ifdef SWIFT_DEBUG_CHECKS
/* Early abort? */
if (c->gcount == 0) error("Doing self gravity on an empty cell !");
@@ -550,6 +809,9 @@ void runner_doself_grav(struct runner *r, struct cell *c, int gettimer) {
TIMER_TIC;
+ /* Anything to do here? */
+ if (!cell_is_active(c, e)) return;
+
/* If the cell is split, interact each progeny with itself, and with
each of its siblings. */
if (c->split) {
@@ -617,8 +879,16 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
/* Some constants */
const struct engine *e = r->e;
+ const struct space *s = e->s;
const struct gravity_props *props = e->gravity_properties;
+ const int periodic = s->periodic;
+ const double cell_width = s->width[0];
+ const double dim[3] = {s->dim[0], s->dim[1], s->dim[2]};
const double theta_crit_inv = props->theta_crit_inv;
+ const double max_distance = props->a_smooth * props->r_cut_max * cell_width;
+ const double max_distance2 = max_distance * max_distance;
+ struct gravity_tensors *const mi = ci->multipole;
+ const double CoM[3] = {mi->CoM[0], mi->CoM[1], mi->CoM[2]};
TIMER_TIC;
@@ -627,7 +897,7 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
const int nr_cells = e->s->nr_cells;
/* Anything to do here? */
- if (!cell_is_active(ci, e)) return; // MATTHIEU (should never happen)
+ if (!cell_is_active(ci, e)) return;
/* Check multipole has been drifted */
if (ci->ti_old_multipole != e->ti_current)
@@ -637,22 +907,41 @@ void runner_do_grav_long_range(struct runner *r, struct cell *ci, int timer) {
* well-separated */
for (int i = 0; i < nr_cells; ++i) {
- /* Handle on the top-level cell */
+ /* Handle on the top-level cell and it's gravity business*/
struct cell *cj = &cells[i];
+ const struct gravity_tensors *const mj = cj->multipole;
/* Avoid stupid cases */
if (ci == cj || cj->gcount == 0) continue;
+ /* Is this interaction part of the periodic mesh calculation ?*/
+ if (periodic) {
+
+ const double dx = nearest(CoM[0] - mj->CoM[0], dim[0]);
+ const double dy = nearest(CoM[1] - mj->CoM[1], dim[1]);
+ const double dz = nearest(CoM[2] - mj->CoM[2], dim[2]);
+ const double r2 = dx * dx + dy * dy + dz * dz;
+
+ /* Are we beyond the distance where the truncated forces are 0 ?*/
+ if (r2 > max_distance2) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Need to account for the interactions we missed */
+ mi->pot.num_interacted += mj->m_pole.num_gpart;
+#endif
+ continue;
+ }
+ }
+
/* Check the multipole acceptance criterion */
- if (gravity_multipole_accept(ci->multipole, cj->multipole, theta_crit_inv,
- 0)) {
+ if (gravity_multipole_accept(mi, mj, theta_crit_inv, periodic, dim)) {
/* Go for a (non-symmetric) M-M calculation */
runner_dopair_grav_mm(r, ci, cj);
}
/* Is the criterion violated now but was OK at the last rebuild ? */
- else if (gravity_multipole_accept(ci->multipole, cj->multipole,
- theta_crit_inv, 1)) {
+ else if (gravity_multipole_accept_rebuild(mi, mj, theta_crit_inv, periodic,
+ dim)) {
/* Alright, we have to take charge of that pair in a different way. */
// MATTHIEU: We should actually open the tree-node here and recurse.
diff --git a/src/scheduler.c b/src/scheduler.c
index 0ae09b44bfd9eb9c5e432dd9d8cc35be6ae154c7..1a1d36d9814c8b1ace73d668311d1986fc304a13 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -620,53 +620,36 @@ static void scheduler_splittask_gravity(struct task *t, struct scheduler *s) {
break;
}
- /* Is this cell even split? */
- if (ci->split) {
-
- /* Make a sub? */
- if (scheduler_dosub && ci->gcount < space_subsize_self) {
-
- /* convert to a self-subtask. */
- t->type = task_type_sub_self;
-
- /* Make sure we have a drift task (MATTHIEU temp. fix) */
- lock_lock(&ci->lock);
- if (ci->drift_gpart == NULL)
- ci->drift_gpart = scheduler_addtask(
- s, task_type_drift_gpart, task_subtype_none, 0, 0, ci, NULL);
- lock_unlock_blind(&ci->lock);
-
- /* Otherwise, make tasks explicitly. */
- } else {
-
- /* Take a step back (we're going to recycle the current task)... */
- redo = 1;
-
- /* Add the self tasks. */
- int first_child = 0;
- while (ci->progeny[first_child] == NULL) first_child++;
- t->ci = ci->progeny[first_child];
- for (int k = first_child + 1; k < 8; k++)
- if (ci->progeny[k] != NULL)
- scheduler_splittask_gravity(
- scheduler_addtask(s, task_type_self, t->subtype, 0, 0,
- ci->progeny[k], NULL),
- s);
-
- /* Make a task for each pair of progeny */
- if (t->subtype != task_subtype_external_grav) {
- for (int j = 0; j < 8; j++)
- if (ci->progeny[j] != NULL)
- for (int k = j + 1; k < 8; k++)
- if (ci->progeny[k] != NULL)
- scheduler_splittask_gravity(
- scheduler_addtask(s, task_type_pair, t->subtype,
- sub_sid_flag[j][k], 0, ci->progeny[j],
- ci->progeny[k]),
- s);
- }
+ /* Should we split this task? */
+ if (ci->split && ci->gcount > space_subsize_self_grav) {
+
+ /* Take a step back (we're going to recycle the current task)... */
+ redo = 1;
+
+ /* Add the self tasks. */
+ int first_child = 0;
+ while (ci->progeny[first_child] == NULL) first_child++;
+ t->ci = ci->progeny[first_child];
+ for (int k = first_child + 1; k < 8; k++)
+ if (ci->progeny[k] != NULL)
+ scheduler_splittask_gravity(
+ scheduler_addtask(s, task_type_self, t->subtype, 0, 0,
+ ci->progeny[k], NULL),
+ s);
+
+ /* Make a task for each pair of progeny */
+ if (t->subtype != task_subtype_external_grav) {
+ for (int j = 0; j < 8; j++)
+ if (ci->progeny[j] != NULL)
+ for (int k = j + 1; k < 8; k++)
+ if (ci->progeny[k] != NULL)
+ scheduler_splittask_gravity(
+ scheduler_addtask(s, task_type_pair, t->subtype,
+ sub_sid_flag[j][k], 0, ci->progeny[j],
+ ci->progeny[k]),
+ s);
}
- } /* Cell is split */
+ }
/* Otherwise, make sure the self task has a drift task */
else {
@@ -694,7 +677,7 @@ static void scheduler_splittask_gravity(struct task *t, struct scheduler *s) {
}
/* Should this task be split-up? */
- if (ci->split && cj->split) {
+ if (0 && ci->split && cj->split) {
// MATTHIEU: nothing here for now
@@ -1364,6 +1347,11 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
} else if (t->subtype == task_subtype_xv ||
t->subtype == task_subtype_rho ||
t->subtype == task_subtype_gradient) {
+#ifdef SWIFT_DEBUG_CHECKS
+ for (int k = 0; k < t->ci->count; k++)
+ if (t->ci->parts[k].ti_drift != s->space->e->ti_current)
+ error("Sending un-drifted particle !");
+#endif
err = MPI_Isend(t->ci->parts, t->ci->count, part_mpi_type,
t->cj->nodeID, t->flags, MPI_COMM_WORLD, &t->req);
// message( "sending %i parts with tag=%i from %i to %i." ,
diff --git a/src/space.c b/src/space.c
index 7fc5efdb14ebaf6e8b00aa90d31730a2ac2d97c5..fc7732f45350b37600198b057583baf5b761c28c 100644
--- a/src/space.c
+++ b/src/space.c
@@ -62,6 +62,7 @@
int space_splitsize = space_splitsize_default;
int space_subsize_pair = space_subsize_pair_default;
int space_subsize_self = space_subsize_self_default;
+int space_subsize_self_grav = space_subsize_self_grav_default;
int space_maxsize = space_maxsize_default;
/**
@@ -2705,6 +2706,9 @@ void space_init(struct space *s, const struct swift_params *params,
params, "Scheduler:cell_sub_size_pair", space_subsize_pair_default);
space_subsize_self = parser_get_opt_param_int(
params, "Scheduler:cell_sub_size_self", space_subsize_self_default);
+ space_subsize_self_grav =
+ parser_get_opt_param_int(params, "Scheduler:cell_sub_size_self_grav",
+ space_subsize_self_grav_default);
space_splitsize = parser_get_opt_param_int(
params, "Scheduler:cell_split_size", space_splitsize_default);
diff --git a/src/space.h b/src/space.h
index fb22d83acd9e416c5a2fc38e4856ba27905d038e..7fe8f27ecaf671738e03881d547d6cb29a61503a 100644
--- a/src/space.h
+++ b/src/space.h
@@ -45,6 +45,7 @@ struct cell;
#define space_maxsize_default 8000000
#define space_subsize_pair_default 256000000
#define space_subsize_self_default 32000
+#define space_subsize_self_grav_default 32000
#define space_max_top_level_cells_default 12
#define space_stretch 1.10f
#define space_maxreldx 0.1f
@@ -57,6 +58,7 @@ extern int space_splitsize;
extern int space_maxsize;
extern int space_subsize_pair;
extern int space_subsize_self;
+extern int space_subsize_self_grav;
/**
* @brief The space in which the cells and particles reside.
diff --git a/src/task.h b/src/task.h
index cd15e09cc5a1b597f96dad7dfeaa479f721f0c8b..dee888c9f16dd69785a31371da15078af4e0af0c 100644
--- a/src/task.h
+++ b/src/task.h
@@ -36,6 +36,7 @@
* @brief The different task types.
*
* Be sure to update the taskID_names array in tasks.c if you modify this list!
+ * Also update the python task plotting scripts!
*/
enum task_types {
task_type_none = 0,
diff --git a/src/tools.c b/src/tools.c
index b2932f3eb09bfe2acaf435a18e25b0421984542a..68b23e1a815b510abafbd8f7c6285d3633ba0a63 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -430,7 +430,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(0.f, r2, fdx, &pi, &pj);
+ 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];
@@ -755,7 +755,7 @@ void gravity_n2(struct gpart *gparts, const int gcount,
const struct gravity_props *gravity_properties, float rlr) {
const float rlr_inv = 1. / rlr;
- const float r_cut = gravity_properties->r_cut;
+ const float r_cut = gravity_properties->r_cut_max;
const float max_d = r_cut * rlr;
const float max_d2 = max_d * max_d;
@@ -790,7 +790,7 @@ void gravity_n2(struct gpart *gparts, const int gcount,
if (r2 < max_d2 || 1) {
/* Apply the gravitational acceleration. */
- runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
+ runner_iact_grav_pp(r2, dx, gpi, gpj);
}
}
}