diff --git a/examples/EAGLE_12/eagle_12.yml b/examples/EAGLE_12/eagle_12.yml
index e4eb43672a2e80da2f7ee558cebbebe92a28e269..6afffed0f9d39b34588b89569a85ab56223fc548 100644
--- a/examples/EAGLE_12/eagle_12.yml
+++ b/examples/EAGLE_12/eagle_12.yml
@@ -31,8 +31,6 @@ 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).
- a_smooth: 1000.
- r_cut: 4.
# Parameters for the hydrodynamics scheme
SPH:
diff --git a/examples/UniformDMBox/uniformBox.yml b/examples/UniformDMBox/uniformBox.yml
index 8d9ec300164a7bf8f3df257c34ee44d4f77fe94e..1572e55662d2ba6b113ab5f21ba357db1bc811d4 100644
--- a/examples/UniformDMBox/uniformBox.yml
+++ b/examples/UniformDMBox/uniformBox.yml
@@ -14,7 +14,7 @@ TimeIntegration:
dt_max: 1e-3 # The maximal time-step size of the simulation (in internal units).
Scheduler:
- max_top_level_cells: 8
+ max_top_level_cells: 16
cell_split_size: 50
# Parameters governing the snapshots
@@ -35,4 +35,4 @@ Statistics:
# Parameters related to the initial conditions
InitialConditions:
- file_name: ./uniformDMBox_100.hdf5 # The file to read
+ file_name: ./uniformDMBox_50.hdf5 # The file to read
diff --git a/src/cell.c b/src/cell.c
index aa5f3687ae02f1b88ffc7d351e42d4dfc1ff1468..1c20add31d8bb07c2cbc015e4dff3cfa98b2310d 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -1453,7 +1453,6 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
if (c->timestep != NULL) scheduler_activate(s, c->timestep);
if (c->grav_down != NULL) scheduler_activate(s, c->grav_down);
if (c->grav_long_range != NULL) scheduler_activate(s, c->grav_long_range);
- if (c->grav_top_level != NULL) scheduler_activate(s, c->grav_top_level);
if (c->cooling != NULL) scheduler_activate(s, c->cooling);
if (c->sourceterms != NULL) scheduler_activate(s, c->sourceterms);
diff --git a/src/cell.h b/src/cell.h
index c89e70960e84bb027f5e99a3cb362f2e0722b9bd..d04e20ca33cd1f0ffff5759977f9fc52eab630fc 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -169,10 +169,7 @@ struct cell {
/*! The task to compute time-steps */
struct task *timestep;
- /*! Task constructing the multipole from the particles */
- struct task *grav_top_level;
-
- /*! Task constructing the multipole from the particles */
+ /*! Task computing long range non-periodic gravity interactions */
struct task *grav_long_range;
/*! Task propagating the multipole to the particles */
diff --git a/src/engine.c b/src/engine.c
index ab14af5c2e3fee19fcf2ec5963788aa8338d0568..03d5b4be461bbe9667c2963559df62a81ca17dc5 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -168,18 +168,12 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
c->grav_long_range = scheduler_addtask(
s, task_type_grav_long_range, task_subtype_none, 0, 0, c, NULL);
- /* Gravity top-level periodic calculation */
- c->grav_top_level = scheduler_addtask(s, task_type_grav_top_level,
- task_subtype_none, 0, 0, c, NULL);
-
/* Gravity recursive down-pass */
c->grav_down = scheduler_addtask(s, task_type_grav_down,
task_subtype_none, 0, 0, c, NULL);
scheduler_addunlock(s, c->init_grav, c->grav_long_range);
- scheduler_addunlock(s, c->init_grav, c->grav_top_level);
scheduler_addunlock(s, c->grav_long_range, c->grav_down);
- scheduler_addunlock(s, c->grav_top_level, c->grav_down);
scheduler_addunlock(s, c->grav_down, c->kick2);
}
@@ -1655,6 +1649,10 @@ void engine_make_self_gravity_tasks(struct engine *e) {
struct cell *cells = s->cells_top;
const int nr_cells = s->nr_cells;
+ /* Create the FFT task for this MPI rank */
+ struct task *fft = scheduler_addtask(sched, task_type_grav_top_level,
+ task_type_none, 0, 0, NULL, NULL);
+
for (int cid = 0; cid < nr_cells; ++cid) {
struct cell *ci = &cells[cid];
@@ -1666,7 +1664,10 @@ void engine_make_self_gravity_tasks(struct engine *e) {
if (ci->nodeID != nodeID) continue;
/* If the cells is local build a self-interaction */
- scheduler_addtask(sched, task_type_self, task_subtype_grav, 0, 0, ci, NULL);
+ struct task *self = scheduler_addtask(sched, task_type_self,
+ task_subtype_grav, 0, 0, ci, NULL);
+
+ scheduler_addunlock(sched, fft, self);
/* Loop over every other cell */
for (int cjd = cid + 1; cjd < nr_cells; ++cjd) {
@@ -1681,9 +1682,12 @@ void engine_make_self_gravity_tasks(struct engine *e) {
/* Are the cells to close for a MM interaction ? */
if (!gravity_multipole_accept(ci->multipole, cj->multipole,
- theta_crit_inv, 1))
+ theta_crit_inv, 1)) {
scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0, ci,
cj);
+
+ // scheduler_addunlock(sched, fft, pair);
+ }
}
}
}
@@ -2605,11 +2609,15 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
/* Gravity ? */
else if (t->type == task_type_grav_down ||
- t->type == task_type_grav_long_range ||
- t->type == task_type_grav_top_level) {
+ t->type == task_type_grav_long_range) {
if (cell_is_active(t->ci, e)) scheduler_activate(s, t);
}
+ /* Periodic gravity ? */
+ else if (t->type == task_type_grav_top_level) {
+ scheduler_activate(s, t);
+ }
+
/* Time-step? */
else if (t->type == task_type_timestep) {
t->ci->updated = 0;
diff --git a/src/gravity_properties.c b/src/gravity_properties.c
index 7b9b8cd7c35f8fa9b21ff34ce2589b5d45ce8393..b1098888b96cdef2205ed513e60a3799c63e8b9f 100644
--- a/src/gravity_properties.c
+++ b/src/gravity_properties.c
@@ -69,11 +69,9 @@ void gravity_props_print(const struct gravity_props *p) {
message("Self-gravity softening: epsilon=%.4f (Plummer equivalent: %.4f)",
p->epsilon, p->epsilon / 3.);
- if (p->a_smooth != gravity_props_default_a_smooth)
- message("Self-gravity MM smoothing-scale: a_smooth=%f", p->a_smooth);
+ message("Self-gravity MM smoothing-scale: a_smooth=%f", p->a_smooth);
- if (p->r_cut != gravity_props_default_r_cut)
- message("Self-gravity MM cut-off: r_cut=%f", p->r_cut);
+ message("Self-gravity MM cut-off: r_cut=%f", p->r_cut);
}
#if defined(HAVE_HDF5)
diff --git a/src/runner.c b/src/runner.c
index 8d54dc07c89754f14caf39583f1123f620b03547..ca94d410bd739f439091326e336ff1f9820e489b 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -1728,10 +1728,12 @@ void *runner_main(void *data) {
#ifdef SWIFT_DEBUG_CHECKS
t->ti_run = e->ti_current;
#ifndef WITH_MPI
- if (ci == NULL && cj == NULL) {
+ if (t->type == task_type_grav_top_level) {
+ if (ci != NULL || cj != NULL)
+ error("Top-level gravity task associated with a cell");
+ } else if (ci == NULL && cj == NULL) {
error("Task not associated with cells!");
-
} else if (cj == NULL) { /* self */
if (!cell_is_active(ci, e) && t->type != task_type_sort &&
diff --git a/src/runner_doiact_fft.c b/src/runner_doiact_fft.c
index 1e9a4d3be684a6fdec921a4b1ef345895a41d5a9..f0d49a6d3a4c643eb319c08bd2c21af1e2031465 100644
--- a/src/runner_doiact_fft.c
+++ b/src/runner_doiact_fft.c
@@ -47,7 +47,7 @@
*/
__attribute__((always_inline)) INLINE static int row_major_id(int i, int j,
int k, int N) {
- return (i * N * N + j * N + k);
+ return ((i % N) * N * N + (j % N) * N + (k % N));
}
/**
@@ -76,6 +76,12 @@ __attribute__((always_inline)) INLINE static void multipole_to_mesh_CIC(
const double dz = fac * m->CoM[2] - k;
const double tz = 1. - dz;
+#ifdef SWIFT_DEBUG_CHECKS
+ if (i < 0 || i >= N) error("Invalid multipole position in x");
+ if (j < 0 || j >= N) error("Invalid multipole position in y");
+ if (k < 0 || k >= N) error("Invalid multipole position in z");
+#endif
+
/* CIC ! */
rho[row_major_id(i + 0, j + 0, k + 0, N)] += m->m_pole.M_000 * tx * ty * tz;
rho[row_major_id(i + 0, j + 0, k + 1, N)] += m->m_pole.M_000 * tx * ty * dz;
@@ -114,6 +120,12 @@ __attribute__((always_inline)) INLINE static void mesh_to_multipole_CIC(
const double dz = fac * m->CoM[2] - k;
const double tz = 1. - dz;
+#ifdef SWIFT_DEBUG_CHECKS
+ if (i < 0 || i >= N) error("Invalid multipole position in x");
+ if (j < 0 || j >= N) error("Invalid multipole position in y");
+ if (k < 0 || k >= N) error("Invalid multipole position in z");
+#endif
+
/* CIC ! */
m->pot.F_000 += pot[row_major_id(i + 0, j + 0, k + 0, N)] * tx * ty * tz;
m->pot.F_000 += pot[row_major_id(i + 0, j + 0, k + 1, N)] * tx * ty * dz;
@@ -149,7 +161,6 @@ void runner_do_grav_fft(struct runner* r, int timer) {
/* Some useful constants */
const int N = cdim[0];
const int N_half = N / 2;
- const int Ncells = N * N * N;
const double cell_fac = N / box_size;
/* Recover the list of top-level multipoles */
@@ -164,7 +175,7 @@ void runner_do_grav_fft(struct runner* r, int timer) {
}
/* Allocates some memory for the density mesh */
- double* restrict rho = fftw_alloc_real(Ncells);
+ double* restrict rho = fftw_alloc_real(N * N * N);
if (rho == NULL) error("Error allocating memory for density mesh");
/* Allocates some memory for the mesh in Fourier space */
@@ -179,7 +190,7 @@ void runner_do_grav_fft(struct runner* r, int timer) {
N, N, N, frho, rho, FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
/* Do a CIC mesh assignment of the multipoles */
- bzero(rho, Ncells * sizeof(double));
+ bzero(rho, N * N * N * sizeof(double));
for (int i = 0; i < nr_cells; ++i)
multipole_to_mesh_CIC(&multipoles[i], rho, N, cell_fac);
diff --git a/src/scheduler.c b/src/scheduler.c
index e1074a32e9459d0b4c6439055fcc92ac41e6e514..2d1e0093cc2b61fe2485c2295f7a3ff41de30cf6 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -137,7 +137,8 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
redo = 0;
/* Non-splittable task? */
- if ((t->ci == NULL || (t->type == task_type_pair && t->cj == NULL)) ||
+ if ((t->ci == NULL && t->type != task_type_grav_top_level) ||
+ ((t->type == task_type_pair && t->cj == NULL)) ||
((t->type == task_type_kick1) && t->ci->nodeID != s->nodeID) ||
((t->type == task_type_kick2) && t->ci->nodeID != s->nodeID) ||
((t->type == task_type_drift) && t->ci->nodeID != s->nodeID) ||
@@ -1133,6 +1134,9 @@ void scheduler_start(struct scheduler *s) {
/* Don't check MPI stuff */
if (t->type == task_type_send || t->type == task_type_recv) continue;
+ /* Don't check the FFT task */
+ if (t->type == task_type_grav_top_level) continue;
+
if (ci == NULL && cj == NULL) {
error("Task not associated with cells!");
diff --git a/src/space.c b/src/space.c
index 0c67fe27ac1c7b2cae3672e7f0b1ce82be5fb804..1d957e7e72461185b0ddbbe60da3cb2b5b8832b1 100644
--- a/src/space.c
+++ b/src/space.c
@@ -26,6 +26,7 @@
/* Some standard headers. */
#include <float.h>
+#include <hdf5.h>
#include <limits.h>
#include <math.h>
#include <stdlib.h>
@@ -220,7 +221,6 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->drift = NULL;
c->cooling = NULL;
c->sourceterms = NULL;
- c->grav_top_level = NULL;
c->grav_long_range = NULL;
c->grav_down = NULL;
c->super = c;
@@ -2943,3 +2943,47 @@ void space_clean(struct space *s) {
free(s->gparts);
free(s->sparts);
}
+
+void space_dump_multipoles(struct space *s) {
+
+ hid_t h_file =
+ H5Fcreate("multipoles.hdf5", H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
+ hid_t h_grp =
+ H5Gcreate(h_file, "/Multipoles", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
+
+ const hid_t h_space = H5Screate(H5S_SIMPLE);
+ const hid_t h_space2 = H5Screate(H5S_SIMPLE);
+ hsize_t shape[2] = {s->nr_cells, 3};
+ hsize_t shape2[1] = {s->nr_cells};
+ H5Sset_extent_simple(h_space, 2, shape, NULL);
+ H5Sset_extent_simple(h_space2, 1, shape2, NULL);
+
+ const hid_t h_data =
+ H5Dcreate(h_grp, "Coordinates", H5T_NATIVE_DOUBLE, h_space, H5P_DEFAULT,
+ H5P_DEFAULT, H5P_DEFAULT);
+
+ const hid_t h_data2 =
+ H5Dcreate(h_grp, "Potential", H5T_NATIVE_DOUBLE, h_space2, H5P_DEFAULT,
+ H5P_DEFAULT, H5P_DEFAULT);
+
+ double *temp = malloc(s->nr_cells * 3 * sizeof(double));
+ double *temp2 = malloc(s->nr_cells * 1 * sizeof(double));
+ for (int i = 0; i < s->nr_cells; ++i) {
+ temp[i * 3 + 0] = s->multipoles_top[i].CoM[0];
+ temp[i * 3 + 1] = s->multipoles_top[i].CoM[1];
+ temp[i * 3 + 2] = s->multipoles_top[i].CoM[2];
+ temp2[i] = s->multipoles_top[i].pot.F_000;
+ }
+
+ H5Dwrite(h_data, H5T_NATIVE_DOUBLE, h_space, H5S_ALL, H5P_DEFAULT, temp);
+ H5Dwrite(h_data2, H5T_NATIVE_DOUBLE, h_space2, H5S_ALL, H5P_DEFAULT, temp2);
+
+ free(temp);
+ free(temp2);
+ H5Dclose(h_data);
+ H5Dclose(h_data2);
+ H5Sclose(h_space);
+ H5Sclose(h_space2);
+ H5Gclose(h_grp);
+ H5Fclose(h_file);
+}
diff --git a/src/space.h b/src/space.h
index c5f588563e5a9fb4b6cb73ac1446514f8149794f..86dc6ed5d56a5328e9e49555ac7d4b2bf5990c5e 100644
--- a/src/space.h
+++ b/src/space.h
@@ -221,5 +221,6 @@ void space_check_timesteps(struct space *s);
void space_replicate(struct space *s, int replicate, int verbose);
void space_reset_task_counters(struct space *s);
void space_clean(struct space *s);
+void space_dump_multipoles(struct space *s);
#endif /* SWIFT_SPACE_H */
diff --git a/src/task.c b/src/task.c
index e8c35e49a57595a6415c60ce7071ae1c0e3f09b7..1bb3df7ea9817488f6391ef2bff7686835deed46 100644
--- a/src/task.c
+++ b/src/task.c
@@ -323,7 +323,6 @@ void task_unlock(struct task *t) {
cell_munlocktree(ci);
break;
- case task_type_grav_top_level:
case task_type_grav_long_range:
case task_type_grav_mm:
cell_munlocktree(ci);
@@ -442,7 +441,6 @@ int task_lock(struct task *t) {
}
break;
- case task_type_grav_top_level:
case task_type_grav_long_range:
case task_type_grav_mm:
/* Lock the m-poles */