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/makeIC.py b/examples/UniformDMBox/makeIC.py
index 8e032500016eb6cc8e0decc54968bb5b841d7f93..8f3cd943b3cf19c4ae231d125c5ef97d076e0e8e 100644
--- a/examples/UniformDMBox/makeIC.py
+++ b/examples/UniformDMBox/makeIC.py
@@ -26,7 +26,7 @@ from numpy import *
# with a density of 1
# Parameters
-periodic= 0 # 1 For periodic box
+periodic= 1 # 1 For periodic box
boxSize = 1.
rho = 1.
L = int(sys.argv[1]) # Number of particles along one axis
diff --git a/examples/UniformDMBox/uniformBox.yml b/examples/UniformDMBox/uniformBox.yml
index 8d9ec300164a7bf8f3df257c34ee44d4f77fe94e..cffd442a9a5b16d8e042e41caf9991fcf0e1202e 100644
--- a/examples/UniformDMBox/uniformBox.yml
+++ b/examples/UniformDMBox/uniformBox.yml
@@ -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/Makefile.am b/src/Makefile.am
index 7bec5327f4759fcf7d3e1af9d041677ffbc7ab55..809f5c34642392d001f1987c1c8926b8b97eae1e 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -63,7 +63,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
kernel_long_gravity.h vector.h cache.h runner_doiact.h runner_doiact_vec.h runner_doiact_grav.h runner_doiact_fft.h \
runner_doiact_nosort.h units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
- dimension.h equation_of_state.h part_type.h \
+ dimension.h equation_of_state.h part_type.h periodic.h \
gravity.h gravity_io.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
diff --git a/src/active.h b/src/active.h
index 02e504f762735994e6c57f7e155071fede016713..58e88835b6f51ae15f9fd7270c0e1f89bbd6d61a 100644
--- a/src/active.h
+++ b/src/active.h
@@ -29,25 +29,48 @@
#include "timeline.h"
/**
- * @brief Check that a cell been drifted to the current time.
+ * @brief Check that the #part in a #cell have been drifted to the current time.
*
* @param c The #cell.
* @param e The #engine containing information about the current time.
* @return 1 if the #cell has been drifted to the current time, 0 otherwise.
*/
-__attribute__((always_inline)) INLINE static int cell_is_drifted(
+__attribute__((always_inline)) INLINE static int cell_are_part_drifted(
const struct cell *c, const struct engine *e) {
#ifdef SWIFT_DEBUG_CHECKS
- if (c->ti_old > e->ti_current)
+ if (c->ti_old_part > e->ti_current)
error(
"Cell has been drifted too far forward in time! c->ti_old=%lld (t=%e) "
"and e->ti_current=%lld (t=%e)",
- c->ti_old, c->ti_old * e->timeBase, e->ti_current,
+ c->ti_old_part, c->ti_old_part * e->timeBase, e->ti_current,
e->ti_current * e->timeBase);
#endif
- return (c->ti_old == e->ti_current);
+ return (c->ti_old_part == e->ti_current);
+}
+
+/**
+ * @brief Check that the #gpart in a #cell have been drifted to the current
+ * time.
+ *
+ * @param c The #cell.
+ * @param e The #engine containing information about the current time.
+ * @return 1 if the #cell has been drifted to the current time, 0 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int cell_are_gpart_drifted(
+ const struct cell *c, const struct engine *e) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->ti_old_gpart > e->ti_current)
+ error(
+ "Cell has been drifted too far forward in time! c->ti_old=%lld (t=%e) "
+ "and e->ti_current=%lld (t=%e)",
+ c->ti_old_gpart, c->ti_old_gpart * e->timeBase, e->ti_current,
+ e->ti_current * e->timeBase);
+#endif
+
+ return (c->ti_old_gpart == e->ti_current);
}
/* Are cells / particles active for regular tasks ? */
diff --git a/src/cell.c b/src/cell.c
index aa5f3687ae02f1b88ffc7d351e42d4dfc1ff1468..657cd85cea63635c532bd05e4425600e2f7fefb5 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -99,7 +99,8 @@ int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {
c->h_max = pc->h_max;
c->ti_end_min = pc->ti_end_min;
c->ti_end_max = pc->ti_end_max;
- c->ti_old = pc->ti_old;
+ c->ti_old_part = pc->ti_old_part;
+ c->ti_old_gpart = pc->ti_old_gpart;
c->count = pc->count;
c->gcount = pc->gcount;
c->scount = pc->scount;
@@ -128,7 +129,8 @@ int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {
if (k & 1) temp->loc[2] += temp->width[2];
temp->depth = c->depth + 1;
temp->split = 0;
- temp->dx_max = 0.f;
+ temp->dx_max_part = 0.f;
+ temp->dx_max_gpart = 0.f;
temp->dx_max_sort = 0.f;
temp->nodeID = c->nodeID;
temp->parent = c;
@@ -239,7 +241,8 @@ int cell_pack(struct cell *c, struct pcell *pc) {
pc->h_max = c->h_max;
pc->ti_end_min = c->ti_end_min;
pc->ti_end_max = c->ti_end_max;
- pc->ti_old = c->ti_old;
+ pc->ti_old_part = c->ti_old_part;
+ pc->ti_old_gpart = c->ti_old_gpart;
pc->count = c->count;
pc->gcount = c->gcount;
pc->scount = c->scount;
@@ -1018,7 +1021,7 @@ void cell_clean_links(struct cell *c, void *data) {
}
/**
- * @brief Checks that the particles in a cell are at the
+ * @brief Checks that the #part in a cell are at the
* current point in time
*
* Calls error() if the cell is not at the current time.
@@ -1026,7 +1029,7 @@ void cell_clean_links(struct cell *c, void *data) {
* @param c Cell to act upon
* @param data The current time on the integer time-line
*/
-void cell_check_particle_drift_point(struct cell *c, void *data) {
+void cell_check_part_drift_point(struct cell *c, void *data) {
#ifdef SWIFT_DEBUG_CHECKS
@@ -1035,14 +1038,40 @@ void cell_check_particle_drift_point(struct cell *c, void *data) {
/* Only check local cells */
if (c->nodeID != engine_rank) return;
- if (c->ti_old != ti_drift)
- error("Cell in an incorrect time-zone! c->ti_old=%lld ti_drift=%lld",
- c->ti_old, ti_drift);
+ if (c->ti_old_part != ti_drift)
+ error("Cell in an incorrect time-zone! c->ti_old_part=%lld ti_drift=%lld",
+ c->ti_old_part, ti_drift);
for (int i = 0; i < c->count; ++i)
if (c->parts[i].ti_drift != ti_drift)
error("part in an incorrect time-zone! p->ti_drift=%lld ti_drift=%lld",
c->parts[i].ti_drift, ti_drift);
+#else
+ error("Calling debugging code without debugging flag activated.");
+#endif
+}
+
+/**
+ * @brief Checks that the #gpart and #spart in a cell are at the
+ * current point in time
+ *
+ * Calls error() if the cell is not at the current time.
+ *
+ * @param c Cell to act upon
+ * @param data The current time on the integer time-line
+ */
+void cell_check_gpart_drift_point(struct cell *c, void *data) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+ const integertime_t ti_drift = *(integertime_t *)data;
+
+ /* Only check local cells */
+ if (c->nodeID != engine_rank) return;
+
+ if (c->ti_old_gpart != ti_drift)
+ error("Cell in an incorrect time-zone! c->ti_old_gpart=%lld ti_drift=%lld",
+ c->ti_old_gpart, ti_drift);
for (int i = 0; i < c->gcount; ++i)
if (c->gparts[i].ti_drift != ti_drift)
@@ -1327,7 +1356,7 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
error("bad flags in sort task.");
#endif
scheduler_activate(s, ci->sorts);
- if (ci->nodeID == engine_rank) scheduler_activate(s, ci->drift);
+ if (ci->nodeID == engine_rank) scheduler_activate(s, ci->drift_part);
}
if (!(cj->sorted & (1 << t->flags))) {
#ifdef SWIFT_DEBUG_CHECKS
@@ -1335,7 +1364,7 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
error("bad flags in sort task.");
#endif
scheduler_activate(s, cj->sorts);
- if (cj->nodeID == engine_rank) scheduler_activate(s, cj->drift);
+ if (cj->nodeID == engine_rank) scheduler_activate(s, cj->drift_part);
}
}
@@ -1344,7 +1373,8 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
/* Check whether there was too much particle motion, i.e. the
cell neighbour conditions were violated. */
- if (max(ci->h_max, cj->h_max) + ci->dx_max + cj->dx_max > cj->dmin)
+ if (max(ci->h_max, cj->h_max) + ci->dx_max_part + cj->dx_max_part >
+ cj->dmin)
rebuild = 1;
#ifdef WITH_MPI
@@ -1367,11 +1397,11 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
scheduler_activate(s, l->t);
/* Drift both cells, the foreign one at the level which it is sent. */
- if (l->t->ci->drift)
- scheduler_activate(s, l->t->ci->drift);
+ if (l->t->ci->drift_part)
+ scheduler_activate(s, l->t->ci->drift_part);
else
error("Drift task missing !");
- if (t->type == task_type_pair) scheduler_activate(s, cj->drift);
+ if (t->type == task_type_pair) scheduler_activate(s, cj->drift_part);
if (cell_is_active(cj, e)) {
for (l = cj->send_rho; l != NULL && l->t->cj->nodeID != ci->nodeID;
@@ -1405,11 +1435,11 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
scheduler_activate(s, l->t);
/* Drift both cells, the foreign one at the level which it is sent. */
- if (l->t->ci->drift)
- scheduler_activate(s, l->t->ci->drift);
+ if (l->t->ci->drift_part)
+ scheduler_activate(s, l->t->ci->drift_part);
else
error("Drift task missing !");
- if (t->type == task_type_pair) scheduler_activate(s, ci->drift);
+ if (t->type == task_type_pair) scheduler_activate(s, ci->drift_part);
if (cell_is_active(ci, e)) {
for (l = ci->send_rho; l != NULL && l->t->cj->nodeID != cj->nodeID;
@@ -1425,13 +1455,13 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
scheduler_activate(s, l->t);
}
} else if (t->type == task_type_pair) {
- scheduler_activate(s, ci->drift);
- scheduler_activate(s, cj->drift);
+ scheduler_activate(s, ci->drift_part);
+ scheduler_activate(s, cj->drift_part);
}
#else
if (t->type == task_type_pair) {
- scheduler_activate(s, ci->drift);
- scheduler_activate(s, cj->drift);
+ scheduler_activate(s, ci->drift_part);
+ scheduler_activate(s, cj->drift_part);
}
#endif
}
@@ -1447,13 +1477,15 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
if (c->extra_ghost != NULL) scheduler_activate(s, c->extra_ghost);
if (c->ghost != NULL) scheduler_activate(s, c->ghost);
if (c->init_grav != NULL) scheduler_activate(s, c->init_grav);
- if (c->drift != NULL) scheduler_activate(s, c->drift);
+ if (c->drift_part != NULL) scheduler_activate(s, c->drift_part);
+ if (c->drift_gpart != NULL) scheduler_activate(s, c->drift_gpart);
if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
if (c->kick2 != NULL) scheduler_activate(s, c->kick2);
if (c->timestep != NULL) scheduler_activate(s, c->timestep);
+ if (c->grav_ghost[0] != NULL) scheduler_activate(s, c->grav_ghost[0]);
+ if (c->grav_ghost[1] != NULL) scheduler_activate(s, c->grav_ghost[1]);
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);
@@ -1481,30 +1513,28 @@ void cell_set_super(struct cell *c, struct cell *super) {
}
/**
- * @brief Recursively drifts particles of all kinds in a cell hierarchy.
+ * @brief Recursively drifts the #part in a cell hierarchy.
*
* @param c The #cell.
* @param e The #engine (to get ti_current).
*/
-void cell_drift_particles(struct cell *c, const struct engine *e) {
+void cell_drift_part(struct cell *c, const struct engine *e) {
const float hydro_h_max = e->hydro_properties->h_max;
const double timeBase = e->timeBase;
- const integertime_t ti_old = c->ti_old;
+ const integertime_t ti_old_part = c->ti_old_part;
const integertime_t ti_current = e->ti_current;
struct part *const parts = c->parts;
struct xpart *const xparts = c->xparts;
- struct gpart *const gparts = c->gparts;
- struct spart *const sparts = c->sparts;
/* Drift from the last time the cell was drifted to the current time */
- const double dt = (ti_current - ti_old) * timeBase;
+ const double dt = (ti_current - ti_old_part) * timeBase;
float dx_max = 0.f, dx2_max = 0.f;
float dx_max_sort = 0.0f, dx2_max_sort = 0.f;
float cell_h_max = 0.f;
/* Check that we are actually going to move forward. */
- if (ti_current < ti_old) error("Attempt to drift to the past");
+ if (ti_current < ti_old_part) error("Attempt to drift to the past");
/* Are we not in a leaf ? */
if (c->split) {
@@ -1515,37 +1545,15 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
struct cell *cp = c->progeny[k];
/* Collect */
- cell_drift_particles(cp, e);
+ cell_drift_part(cp, e);
/* Update */
- dx_max = max(dx_max, cp->dx_max);
+ dx_max = max(dx_max, cp->dx_max_part);
dx_max_sort = max(dx_max_sort, cp->dx_max_sort);
cell_h_max = max(cell_h_max, cp->h_max);
}
- } else if (ti_current > ti_old) {
-
- /* Loop over all the g-particles in the cell */
- const size_t nr_gparts = c->gcount;
- for (size_t k = 0; k < nr_gparts; k++) {
-
- /* Get a handle on the gpart. */
- struct gpart *const gp = &gparts[k];
-
- /* Drift... */
- drift_gpart(gp, dt, timeBase, ti_old, ti_current);
-
- /* Compute (square of) motion since last cell construction */
- const float dx2 = gp->x_diff[0] * gp->x_diff[0] +
- gp->x_diff[1] * gp->x_diff[1] +
- gp->x_diff[2] * gp->x_diff[2];
- dx2_max = max(dx2_max, dx2);
-
- /* Init gravity force fields. */
- if (gpart_is_active(gp, e)) {
- gravity_init_gpart(gp);
- }
- }
+ } else if (ti_current > ti_old_part) {
/* Loop over all the gas particles in the cell */
const size_t nr_parts = c->count;
@@ -1556,7 +1564,7 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
struct xpart *const xp = &xparts[k];
/* Drift... */
- drift_part(p, xp, dt, timeBase, ti_old, ti_current);
+ drift_part(p, xp, dt, timeBase, ti_old_part, ti_current);
/* Limit h to within the allowed range */
p->h = min(p->h, hydro_h_max);
@@ -1580,6 +1588,86 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
}
}
+ /* Now, get the maximal particle motion from its square */
+ dx_max = sqrtf(dx2_max);
+ dx_max_sort = sqrtf(dx2_max_sort);
+
+ } else {
+
+ cell_h_max = c->h_max;
+ dx_max = c->dx_max_part;
+ dx_max_sort = c->dx_max_sort;
+ }
+
+ /* Store the values */
+ c->h_max = cell_h_max;
+ c->dx_max_part = dx_max;
+ c->dx_max_sort = dx_max_sort;
+
+ /* Update the time of the last drift */
+ c->ti_old_part = ti_current;
+}
+
+/**
+ * @brief Recursively drifts the #gpart in a cell hierarchy.
+ *
+ * @param c The #cell.
+ * @param e The #engine (to get ti_current).
+ */
+void cell_drift_gpart(struct cell *c, const struct engine *e) {
+
+ const double timeBase = e->timeBase;
+ const integertime_t ti_old_gpart = c->ti_old_gpart;
+ const integertime_t ti_current = e->ti_current;
+ struct gpart *const gparts = c->gparts;
+ struct spart *const sparts = c->sparts;
+
+ /* Drift from the last time the cell was drifted to the current time */
+ const double dt = (ti_current - ti_old_gpart) * timeBase;
+ float dx_max = 0.f, dx2_max = 0.f;
+
+ /* Check that we are actually going to move forward. */
+ if (ti_current < ti_old_gpart) error("Attempt to drift to the past");
+
+ /* Are we not in a leaf ? */
+ if (c->split) {
+
+ /* Loop over the progeny and collect their data. */
+ for (int k = 0; k < 8; k++)
+ if (c->progeny[k] != NULL) {
+ struct cell *cp = c->progeny[k];
+
+ /* Recurse */
+ cell_drift_gpart(cp, e);
+
+ /* Update */
+ dx_max = max(dx_max, cp->dx_max_gpart);
+ }
+
+ } else if (ti_current > ti_old_gpart) {
+
+ /* Loop over all the g-particles in the cell */
+ const size_t nr_gparts = c->gcount;
+ for (size_t k = 0; k < nr_gparts; k++) {
+
+ /* Get a handle on the gpart. */
+ struct gpart *const gp = &gparts[k];
+
+ /* Drift... */
+ drift_gpart(gp, dt, timeBase, ti_old_gpart, ti_current);
+
+ /* Compute (square of) motion since last cell construction */
+ const float dx2 = gp->x_diff[0] * gp->x_diff[0] +
+ gp->x_diff[1] * gp->x_diff[1] +
+ gp->x_diff[2] * gp->x_diff[2];
+ dx2_max = max(dx2_max, dx2);
+
+ /* Init gravity force fields. */
+ if (gpart_is_active(gp, e)) {
+ gravity_init_gpart(gp);
+ }
+ }
+
/* Loop over all the star particles in the cell */
const size_t nr_sparts = c->scount;
for (size_t k = 0; k < nr_sparts; k++) {
@@ -1588,29 +1676,24 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
struct spart *const sp = &sparts[k];
/* Drift... */
- drift_spart(sp, dt, timeBase, ti_old, ti_current);
+ drift_spart(sp, dt, timeBase, ti_old_gpart, ti_current);
/* Note: no need to compute dx_max as all spart have a gpart */
}
/* Now, get the maximal particle motion from its square */
dx_max = sqrtf(dx2_max);
- dx_max_sort = sqrtf(dx2_max_sort);
} else {
- cell_h_max = c->h_max;
- dx_max = c->dx_max;
- dx_max_sort = c->dx_max_sort;
+ dx_max = c->dx_max_gpart;
}
/* Store the values */
- c->h_max = cell_h_max;
- c->dx_max = dx_max;
- c->dx_max_sort = dx_max_sort;
+ c->dx_max_gpart = dx_max;
/* Update the time of the last drift */
- c->ti_old = ti_current;
+ c->ti_old_gpart = ti_current;
}
/**
diff --git a/src/cell.h b/src/cell.h
index c89e70960e84bb027f5e99a3cb362f2e0722b9bd..2e32533402110040310be88629d0fb33f0128c62 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -74,7 +74,7 @@ struct pcell {
/* Stats on this cell's particles. */
double h_max;
- integertime_t ti_end_min, ti_end_max, ti_beg_max, ti_old;
+ integertime_t ti_end_min, ti_end_max, ti_beg_max, ti_old_part, ti_old_gpart;
/* Number of particles in this cell. */
int count, gcount, scount;
@@ -157,8 +157,11 @@ struct cell {
/*! The extra ghost task for complex hydro schemes */
struct task *extra_ghost;
- /*! The drift task */
- struct task *drift;
+ /*! The drift task for parts */
+ struct task *drift_part;
+
+ /*! The drift task for gparts */
+ struct task *drift_gpart;
/*! The first kick task */
struct task *kick1;
@@ -169,10 +172,10 @@ struct cell {
/*! The task to compute time-steps */
struct task *timestep;
- /*! Task constructing the multipole from the particles */
- struct task *grav_top_level;
+ /*! Task linking the FFT mesh to the rest of gravity tasks */
+ struct task *grav_ghost[2];
- /*! 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 */
@@ -233,24 +236,30 @@ struct cell {
/*! Maximum beginning of (integer) time step in this cell. */
integertime_t ti_beg_max;
- /*! Last (integer) time the cell's particle was drifted forward in time. */
- integertime_t ti_old;
-
/*! Last (integer) time the cell's sort arrays were updated. */
integertime_t ti_sort;
+ /*! Last (integer) time the cell's part were drifted forward in time. */
+ integertime_t ti_old_part;
+
+ /*! Last (integer) time the cell's gpart were drifted forward in time. */
+ integertime_t ti_old_gpart;
+
/*! Last (integer) time the cell's multipole was drifted forward in time. */
integertime_t ti_old_multipole;
/*! Minimum dimension, i.e. smallest edge of this cell (min(width)). */
float dmin;
- /*! Maximum particle movement in this cell since last construction. */
- float dx_max;
-
/*! Maximum particle movement in this cell since the last sort. */
float dx_max_sort;
+ /*! Maximum part movement in this cell since last construction. */
+ float dx_max_part;
+
+ /*! Maximum gpart movement in this cell since last construction. */
+ float dx_max_gpart;
+
/*! Nr of #part in this cell. */
int count;
@@ -357,13 +366,15 @@ void cell_clean_links(struct cell *c, void *data);
void cell_make_multipoles(struct cell *c, integertime_t ti_current);
void cell_check_multipole(struct cell *c, void *data);
void cell_clean(struct cell *c);
-void cell_check_particle_drift_point(struct cell *c, void *data);
+void cell_check_part_drift_point(struct cell *c, void *data);
+void cell_check_gpart_drift_point(struct cell *c, void *data);
void cell_check_multipole_drift_point(struct cell *c, void *data);
void cell_reset_task_counters(struct cell *c);
int cell_is_drift_needed(struct cell *c, const struct engine *e);
int cell_unskip_tasks(struct cell *c, struct scheduler *s);
void cell_set_super(struct cell *c, struct cell *super);
-void cell_drift_particles(struct cell *c, const struct engine *e);
+void cell_drift_part(struct cell *c, const struct engine *e);
+void cell_drift_gpart(struct cell *c, const struct engine *e);
void cell_drift_multipole(struct cell *c, const struct engine *e);
void cell_drift_all_multipoles(struct cell *c, const struct engine *e);
void cell_check_timesteps(struct cell *c);
diff --git a/src/debug.c b/src/debug.c
index 3732ee5e769277deb393926ea2dc6f04fba93782..601f63d6e11bbbf95f62eaef1ec6ec7ec06d3ad9 100644
--- a/src/debug.c
+++ b/src/debug.c
@@ -259,8 +259,8 @@ int checkCellhdxmax(const struct cell *c, int *depth) {
message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
result = 0;
}
- if (c->dx_max != dx_max) {
- message("%d Inconsistent dx_max: %f != %f", *depth, c->dx_max, dx_max);
+ if (c->dx_max_part != dx_max) {
+ message("%d Inconsistent dx_max: %f != %f", *depth, c->dx_max_part, dx_max);
message("location: %f %f %f", c->loc[0], c->loc[1], c->loc[2]);
result = 0;
}
diff --git a/src/engine.c b/src/engine.c
index ee4838399924573a765466290b2fb7b1d37bcf87..6bdf3e6d5cfda053f151f080f0f9a0b975fa2e5d 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -133,6 +133,7 @@ void engine_addlink(struct engine *e, struct link **l, struct task *t) {
void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
struct scheduler *s = &e->sched;
+ const int periodic = e->s->periodic;
const int is_hydro = (e->policy & engine_policy_hydro);
const int is_self_gravity = (e->policy & engine_policy_self_gravity);
const int is_with_cooling = (e->policy & engine_policy_cooling);
@@ -168,18 +169,13 @@ 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);
+ if (periodic) scheduler_addunlock(s, c->init_grav, c->grav_ghost[0]);
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);
}
@@ -1049,10 +1045,10 @@ void engine_addtasks_send(struct engine *e, struct cell *ci, struct cell *cj,
#endif
/* Drift before you send */
- if (ci->drift == NULL)
- ci->drift = scheduler_addtask(s, task_type_drift, task_subtype_none, 0,
- 0, ci, NULL);
- scheduler_addunlock(s, ci->drift, t_xv);
+ if (ci->drift_part == NULL)
+ ci->drift_part = scheduler_addtask(s, task_type_drift_part,
+ task_subtype_none, 0, 0, ci, NULL);
+ scheduler_addunlock(s, ci->drift_part, t_xv);
/* The super-cell's timestep task should unlock the send_ti task. */
scheduler_addunlock(s, ci->super->timestep, t_ti);
@@ -1651,41 +1647,98 @@ void engine_make_self_gravity_tasks(struct engine *e) {
struct space *s = e->s;
struct scheduler *sched = &e->sched;
const int nodeID = e->nodeID;
+ const int periodic = s->periodic;
+ 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};
const double theta_crit_inv = e->gravity_properties->theta_crit_inv;
struct cell *cells = s->cells_top;
- const int nr_cells = s->nr_cells;
+ struct task **ghosts = NULL;
+ const int n_ghosts = cdim_ghost[0] * cdim_ghost[1] * cdim_ghost[2] * 2;
+
+ /* Create the top-level task if periodic */
+ if (periodic) {
+
+ /* Create the FFT task for this MPI rank */
+ s->grav_top_level = scheduler_addtask(sched, task_type_grav_top_level,
+ task_subtype_none, 0, 0, NULL, NULL);
+
+ /* Create a grid of ghosts to deal with the dependencies */
+ if ((ghosts = malloc(n_ghosts * sizeof(struct task *))) == 0)
+ error("Error allocating memory for gravity fft ghosts");
+
+ /* Make the ghosts implicit and add the dependencies */
+ for (int n = 0; n < n_ghosts / 2; ++n) {
+ ghosts[2 * n + 0] = scheduler_addtask(sched, task_type_grav_ghost,
+ task_subtype_none, 0, 0, NULL, NULL);
+ ghosts[2 * n + 1] = scheduler_addtask(sched, task_type_grav_ghost,
+ task_subtype_none, 0, 0, NULL, NULL);
+ ghosts[2 * n + 0]->implicit = 1;
+ ghosts[2 * n + 1]->implicit = 1;
+ scheduler_addunlock(sched, ghosts[2 * n + 0], s->grav_top_level);
+ scheduler_addunlock(sched, s->grav_top_level, ghosts[2 * n + 1]);
+ }
+ }
- for (int cid = 0; cid < nr_cells; ++cid) {
+ /* Run through the higher level cells */
+ for (int i = 0; i < cdim[0]; i++) {
+ for (int j = 0; j < cdim[1]; j++) {
+ for (int k = 0; k < cdim[2]; k++) {
- struct cell *ci = &cells[cid];
+ /* Get the cell */
+ const int cid = cell_getid(cdim, i, j, k);
+ struct cell *ci = &cells[cid];
- /* Skip cells without gravity particles */
- if (ci->gcount == 0) continue;
+ /* Skip cells without gravity particles */
+ if (ci->gcount == 0) continue;
- /* Is that cell local ? */
- if (ci->nodeID != nodeID) continue;
+ /* Is that cell local ? */
+ 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);
+
+ /* Deal with periodicity dependencies */
+ const int ghost_id = cell_getid(cdim_ghost, i / 4, j / 4, k / 4);
+ if (ghost_id > n_ghosts) error("Invalid ghost_id");
+ if (periodic) {
+ ci->grav_ghost[0] = ghosts[2 * ghost_id + 0];
+ ci->grav_ghost[1] = ghosts[2 * ghost_id + 1];
+ }
- /* If the cells is local build a self-interaction */
- scheduler_addtask(sched, task_type_self, task_subtype_grav, 0, 0, ci, NULL);
+ /* Loop over every other cell */
+ for (int ii = 0; ii < cdim[0]; ii++) {
+ for (int jj = 0; jj < cdim[1]; jj++) {
+ for (int kk = 0; kk < cdim[2]; kk++) {
- /* Loop over every other cell */
- for (int cjd = cid + 1; cjd < nr_cells; ++cjd) {
+ /* Get the cell */
+ const int cjd = cell_getid(cdim, ii, jj, kk);
+ struct cell *cj = &cells[cjd];
- struct cell *cj = &cells[cjd];
+ /* Avoid duplicates */
+ if (cid <= cjd) continue;
- /* Skip cells without gravity particles */
- if (cj->gcount == 0) continue;
+ /* Skip cells without gravity particles */
+ if (cj->gcount == 0) continue;
- /* Is that neighbour local ? */
- if (cj->nodeID != nodeID) continue; // MATTHIEU
+ /* Is that neighbour local ? */
+ 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))
- scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0, 0, ci,
- cj);
+ /* Are the cells to close for a MM interaction ? */
+ if (!gravity_multipole_accept(ci->multipole, cj->multipole,
+ theta_crit_inv, 1)) {
+
+ scheduler_addtask(sched, task_type_pair, task_subtype_grav, 0,
+ 0, ci, cj);
+ }
+ }
+ }
+ }
+ }
}
}
+ if (periodic) free(ghosts);
}
void engine_make_external_gravity_tasks(struct engine *e) {
@@ -1812,10 +1865,11 @@ void engine_count_and_link_tasks(struct engine *e) {
}
/* Link drift tasks to the next-higher drift task. */
- else if (t->type == task_type_drift) {
+ else if (t->type == task_type_drift_part) {
struct cell *finger = ci->parent;
- while (finger != NULL && finger->drift == NULL) finger = finger->parent;
- if (finger != NULL) scheduler_addunlock(sched, t, finger->drift);
+ while (finger != NULL && finger->drift_part == NULL)
+ finger = finger->parent;
+ if (finger != NULL) scheduler_addunlock(sched, t, finger->drift_part);
}
/* Link self tasks to cells. */
@@ -1906,7 +1960,7 @@ static inline void engine_make_external_gravity_dependencies(
struct scheduler *sched, struct task *gravity, struct cell *c) {
/* init --> external gravity --> kick */
- scheduler_addunlock(sched, c->drift, gravity);
+ scheduler_addunlock(sched, c->drift_gpart, gravity);
scheduler_addunlock(sched, gravity, c->super->kick2);
}
@@ -1920,6 +1974,7 @@ void engine_link_gravity_tasks(struct engine *e) {
struct scheduler *sched = &e->sched;
const int nodeID = e->nodeID;
const int nr_tasks = sched->nr_tasks;
+ const int periodic = e->s->periodic;
for (int k = 0; k < nr_tasks; k++) {
@@ -1930,6 +1985,7 @@ void engine_link_gravity_tasks(struct engine *e) {
if (t->type == task_type_self && t->subtype == task_subtype_grav) {
engine_make_self_gravity_dependencies(sched, t, t->ci);
+ if (periodic) scheduler_addunlock(sched, t->ci->super->grav_ghost[1], t);
}
/* Self-interaction for external gravity ? */
@@ -1945,11 +2001,15 @@ void engine_link_gravity_tasks(struct engine *e) {
if (t->ci->nodeID == nodeID) {
engine_make_self_gravity_dependencies(sched, t, t->ci);
+ if (periodic && t->ci->super < t->cj->super)
+ scheduler_addunlock(sched, t->ci->super->grav_ghost[1], t);
}
if (t->cj->nodeID == nodeID && t->ci->super != t->cj->super) {
engine_make_self_gravity_dependencies(sched, t, t->cj);
+ if (periodic && t->ci->super < t->cj->super)
+ scheduler_addunlock(sched, t->cj->super->grav_ghost[1], t);
}
}
@@ -2072,14 +2132,14 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
/* Sort tasks depend on the drift of the cell. */
if (t->type == task_type_sort && t->ci->nodeID == engine_rank) {
- scheduler_addunlock(sched, t->ci->drift, t);
+ scheduler_addunlock(sched, t->ci->drift_part, t);
}
/* Self-interaction? */
else if (t->type == task_type_self && t->subtype == task_subtype_density) {
/* Make all density tasks depend on the drift. */
- scheduler_addunlock(sched, t->ci->drift, t);
+ scheduler_addunlock(sched, t->ci->drift_part, t);
#ifdef EXTRA_HYDRO_LOOP
/* Start by constructing the task for the second and third hydro loop */
@@ -2115,9 +2175,9 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
/* Make all density tasks depend on the drift. */
if (t->ci->nodeID == engine_rank)
- scheduler_addunlock(sched, t->ci->drift, t);
+ scheduler_addunlock(sched, t->ci->drift_part, t);
if (t->cj->nodeID == engine_rank)
- scheduler_addunlock(sched, t->cj->drift, t);
+ scheduler_addunlock(sched, t->cj->drift_part, t);
#ifdef EXTRA_HYDRO_LOOP
/* Start by constructing the task for the second and third hydro loop */
@@ -2474,7 +2534,8 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
if (t->subtype != task_subtype_density) continue;
/* Too much particle movement? */
- if (max(ci->h_max, cj->h_max) + ci->dx_max + cj->dx_max > cj->dmin)
+ if (max(ci->h_max, cj->h_max) + ci->dx_max_part + cj->dx_max_part >
+ cj->dmin)
*rebuild_space = 1;
/* Set the correct sorting flags */
@@ -2495,7 +2556,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
error("bad flags in sort task.");
#endif
scheduler_activate(s, ci->sorts);
- if (ci->nodeID == engine_rank) scheduler_activate(s, ci->drift);
+ if (ci->nodeID == engine_rank) scheduler_activate(s, ci->drift_part);
}
if (!(cj->sorted & (1 << t->flags))) {
#ifdef SWIFT_DEBUG_CHECKS
@@ -2503,7 +2564,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
error("bad flags in sort task.");
#endif
scheduler_activate(s, cj->sorts);
- if (cj->nodeID == engine_rank) scheduler_activate(s, cj->drift);
+ if (cj->nodeID == engine_rank) scheduler_activate(s, cj->drift_part);
}
}
@@ -2527,11 +2588,11 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
scheduler_activate(s, l->t);
/* Drift both cells, the foreign one at the level which it is sent. */
- if (l->t->ci->drift)
- scheduler_activate(s, l->t->ci->drift);
+ if (l->t->ci->drift_part)
+ scheduler_activate(s, l->t->ci->drift_part);
else
error("Drift task missing !");
- if (t->type == task_type_pair) scheduler_activate(s, cj->drift);
+ if (t->type == task_type_pair) scheduler_activate(s, cj->drift_part);
if (cell_is_active(cj, e)) {
for (l = cj->send_rho; l != NULL && l->t->cj->nodeID != ci->nodeID;
@@ -2565,11 +2626,11 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
scheduler_activate(s, l->t);
/* Drift both cells, the foreign one at the level which it is sent. */
- if (l->t->ci->drift)
- scheduler_activate(s, l->t->ci->drift);
+ if (l->t->ci->drift_part)
+ scheduler_activate(s, l->t->ci->drift_part);
else
error("Drift task missing !");
- if (t->type == task_type_pair) scheduler_activate(s, ci->drift);
+ if (t->type == task_type_pair) scheduler_activate(s, ci->drift_part);
if (cell_is_active(ci, e)) {
for (l = ci->send_rho; l != NULL && l->t->cj->nodeID != cj->nodeID;
@@ -2586,30 +2647,37 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
}
} else if (t->type == task_type_pair) {
- scheduler_activate(s, ci->drift);
- scheduler_activate(s, cj->drift);
+ scheduler_activate(s, ci->drift_part);
+ scheduler_activate(s, cj->drift_part);
}
#else
if (t->type == task_type_pair) {
- scheduler_activate(s, ci->drift);
- scheduler_activate(s, cj->drift);
+ scheduler_activate(s, ci->drift_part);
+ scheduler_activate(s, cj->drift_part);
}
#endif
}
- /* Kick/Drift? */
+ /* Kick/Drift/init ? */
else if (t->type == task_type_kick1 || t->type == task_type_kick2 ||
- t->type == task_type_drift || t->type == task_type_init_grav) {
+ t->type == task_type_drift_part ||
+ t->type == task_type_drift_gpart ||
+ t->type == task_type_init_grav) {
if (cell_is_active(t->ci, e)) scheduler_activate(s, t);
}
/* 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 ||
+ t->type == task_type_grav_ghost) {
+ scheduler_activate(s, t);
+ }
+
/* Time-step? */
else if (t->type == task_type_timestep) {
t->ci->updated = 0;
@@ -3036,10 +3104,12 @@ void engine_skip_force_and_kick(struct engine *e) {
struct task *t = &tasks[i];
/* Skip everything that updates the particles */
- if (t->type == task_type_drift || t->type == task_type_kick1 ||
- t->type == task_type_kick2 || t->type == task_type_timestep ||
- t->subtype == task_subtype_force || t->subtype == task_subtype_grav ||
+ if (t->type == task_type_drift_part || t->type == task_type_drift_gpart ||
+ t->type == task_type_kick1 || t->type == task_type_kick2 ||
+ t->type == task_type_timestep || t->subtype == task_subtype_force ||
+ t->subtype == task_subtype_grav ||
t->type == task_type_grav_long_range ||
+ t->type == task_type_grav_ghost ||
t->type == task_type_grav_top_level || t->type == task_type_grav_down ||
t->type == task_type_cooling || t->type == task_type_sourceterms)
t->skip = 1;
@@ -3060,8 +3130,9 @@ void engine_skip_drift(struct engine *e) {
struct task *t = &tasks[i];
- /* Skip everything that updates the particles */
- if (t->type == task_type_drift) t->skip = 1;
+ /* Skip everything that moves the particles */
+ if (t->type == task_type_drift_part || t->type == task_type_drift_gpart)
+ t->skip = 1;
}
}
@@ -3290,8 +3361,8 @@ void engine_step(struct engine *e) {
if (e->policy & engine_policy_reconstruct_mpoles)
engine_reconstruct_multipoles(e);
- else
- engine_drift_top_multipoles(e);
+ // else
+ // engine_drift_top_multipoles(e);
}
/* Print the number of active tasks ? */
@@ -3403,9 +3474,15 @@ int engine_is_done(struct engine *e) {
void engine_unskip(struct engine *e) {
const ticks tic = getticks();
+
+ /* Activate all the regular tasks */
threadpool_map(&e->threadpool, runner_do_unskip_mapper, e->s->cells_top,
e->s->nr_cells, sizeof(struct cell), 1, e);
+ /* And the top level gravity FFT one */
+ if (e->s->periodic && (e->policy & engine_policy_self_gravity))
+ scheduler_activate(&e->sched, e->s->grav_top_level);
+
if (e->verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
clocks_getunit());
@@ -3429,7 +3506,10 @@ void engine_do_drift_all_mapper(void *map_data, int num_elements,
struct cell *c = &cells[ind];
if (c != NULL && c->nodeID == e->nodeID) {
/* Drift all the particles */
- cell_drift_particles(c, e);
+ cell_drift_part(c, e);
+
+ /* Drift all the g-particles */
+ cell_drift_gpart(c, e);
/* Drift the multipoles */
if (e->policy & engine_policy_self_gravity)
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/minmax.h b/src/minmax.h
index a53093663c79cf4280d136747663552e49c7f1b2..9d92cd71d849dba615fdb05bc342014e0593d989 100644
--- a/src/minmax.h
+++ b/src/minmax.h
@@ -43,18 +43,4 @@
_a > _b ? _a : _b; \
})
-/**
- * @brief Limits the value of x to be between a and b
- *
- * Only wraps once. If x > 2b, the returned value will be larger than b.
- * Similarly for x < -b.
- */
-#define box_wrap(x, a, b) \
- ({ \
- const __typeof__(x) _x = (x); \
- const __typeof__(a) _a = (a); \
- const __typeof__(b) _b = (b); \
- _x < _a ? (_x + _b) : ((_x > _b) ? (_x - _b) : _x); \
- })
-
#endif /* SWIFT_MINMAX_H */
diff --git a/src/multipole.h b/src/multipole.h
index b9d49dcf0fc3b605849f7b058aef14843b73517d..23f5194a30b7316aac15073cba36dc404efa21c1 100644
--- a/src/multipole.h
+++ b/src/multipole.h
@@ -28,7 +28,6 @@
#include <string.h>
/* Includes. */
-//#include "active.h"
#include "align.h"
#include "const.h"
#include "error.h"
@@ -37,8 +36,8 @@
#include "gravity_softened_derivatives.h"
#include "inline.h"
#include "kernel_gravity.h"
-#include "minmax.h"
#include "part.h"
+#include "periodic.h"
#include "vector_power.h"
#define multipole_align 128
diff --git a/src/partition.c b/src/partition.c
index 499efab263a9031b0116f073af8cebd5fef0c2eb..c57918745c11d2858b40eefc218e2551e635d6fb 100644
--- a/src/partition.c
+++ b/src/partition.c
@@ -524,7 +524,7 @@ static void repart_edge_metis(int partweights, int bothweights, int nodeID,
t->type != task_type_sub_self && t->type != task_type_sub_self &&
t->type != task_type_ghost && t->type != task_type_kick1 &&
t->type != task_type_kick2 && t->type != task_type_timestep &&
- t->type != task_type_drift)
+ t->type != task_type_drift_part && t->type != task_type_drift_gpart)
continue;
/* Get the task weight. */
@@ -557,7 +557,7 @@ static void repart_edge_metis(int partweights, int bothweights, int nodeID,
/* Different weights for different tasks. */
if (t->type == task_type_ghost || t->type == task_type_kick1 ||
t->type == task_type_kick2 || t->type == task_type_timestep ||
- t->type == task_type_drift) {
+ t->type == task_type_drift_part || t->type == task_type_drift_gpart) {
/* Particle updates add only to vertex weight. */
if (taskvweights) weights_v[cid] += w;
diff --git a/src/periodic.h b/src/periodic.h
new file mode 100644
index 0000000000000000000000000000000000000000..5874b8742e89c5c93727111adb5b289cff4cb6a6
--- /dev/null
+++ b/src/periodic.h
@@ -0,0 +1,75 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_PERIODIC_H
+#define SWIFT_PERIODIC_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Includes. */
+#include "inline.h"
+
+/**
+ * @brief Limits the value of x to be between a and b
+ *
+ * Only wraps once. If x > 2b, the returned value will be larger than b.
+ * Similarly for x < -b.
+ */
+#define box_wrap(x, a, b) \
+ ({ \
+ const __typeof__(x) _x = (x); \
+ const __typeof__(a) _a = (a); \
+ const __typeof__(b) _b = (b); \
+ _x < _a ? (_x + _b) : ((_x > _b) ? (_x - _b) : _x); \
+ })
+
+/**
+ * @brief Find the smallest distance dx along one axis within a box of size
+ * box_size
+ *
+ * This macro evaluates its arguments exactly once.
+ *
+ * Only wraps once. If dx > 2b, the returned value will be larger than b.
+ * Similarly for dx < -b.
+ *
+ */
+__attribute__((always_inline)) INLINE static double nearest(double dx,
+ double box_size) {
+ return dx > 0.5 * box_size ? (dx - box_size)
+ : ((dx < -0.5 * box_size) ? (dx + box_size) : dx);
+}
+
+/**
+ * @brief Find the smallest distance dx along one axis within a box of size
+ * box_size
+ *
+ * This macro evaluates its arguments exactly once.
+ *
+ * Only wraps once. If dx > 2b, the returned value will be larger than b.
+ * Similarly for dx < -b.
+ *
+ */
+__attribute__((always_inline)) INLINE static float nearestf(float dx,
+ float box_size) {
+ return dx > 0.5f * box_size
+ ? (dx - box_size)
+ : ((dx < -0.5f * box_size) ? (dx + box_size) : dx);
+}
+
+#endif /* SWIFT_PERIODIC_H */
diff --git a/src/runner.c b/src/runner.c
index b4dbe7b377255b27520b611e0ca9e24289b38ba7..dd8cd7f426d80b81bbe4522f5098120731053aad 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -53,6 +53,7 @@
#include "hydro_properties.h"
#include "kick.h"
#include "minmax.h"
+#include "runner_doiact_fft.h"
#include "runner_doiact_vec.h"
#include "scheduler.h"
#include "sort_part.h"
@@ -331,7 +332,7 @@ void runner_do_sort(struct runner *r, struct cell *c, int flags, int clock) {
TIMER_TIC;
/* Check that the particles have been moved to the current time */
- if (!cell_is_drifted(c, r->e)) error("Sorting un-drifted cell");
+ if (!cell_are_part_drifted(c, r->e)) error("Sorting un-drifted cell");
#ifdef SWIFT_DEBUG_CHECKS
/* Make sure the sort flags are consistent (downward). */
@@ -839,19 +840,35 @@ void runner_do_unskip_mapper(void *map_data, int num_elements,
}
}
/**
- * @brief Drift particles in real space.
+ * @brief Drift all part in a cell.
+ *
+ * @param r The runner thread.
+ * @param c The cell.
+ * @param timer Are we timing this ?
+ */
+void runner_do_drift_part(struct runner *r, struct cell *c, int timer) {
+
+ TIMER_TIC;
+
+ cell_drift_part(c, r->e);
+
+ if (timer) TIMER_TOC(timer_drift_part);
+}
+
+/**
+ * @brief Drift all gpart in a cell.
*
* @param r The runner thread.
* @param c The cell.
* @param timer Are we timing this ?
*/
-void runner_do_drift_particles(struct runner *r, struct cell *c, int timer) {
+void runner_do_drift_gpart(struct runner *r, struct cell *c, int timer) {
TIMER_TIC;
- cell_drift_particles(c, r->e);
+ cell_drift_gpart(c, r->e);
- if (timer) TIMER_TOC(timer_drift);
+ if (timer) TIMER_TOC(timer_drift_gpart);
}
/**
@@ -1522,7 +1539,7 @@ void runner_do_recv_part(struct runner *r, struct cell *c, int clear_sorts,
/* ... and store. */
c->ti_end_min = ti_end_min;
c->ti_end_max = ti_end_max;
- c->ti_old = ti_current;
+ c->ti_old_part = ti_current;
c->h_max = h_max;
if (timer) TIMER_TOC(timer_dorecv_part);
@@ -1596,7 +1613,7 @@ void runner_do_recv_gpart(struct runner *r, struct cell *c, int timer) {
/* ... and store. */
c->ti_end_min = ti_end_min;
c->ti_end_max = ti_end_max;
- c->ti_old = ti_current;
+ c->ti_old_gpart = ti_current;
if (timer) TIMER_TOC(timer_dorecv_gpart);
@@ -1669,7 +1686,7 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int timer) {
/* ... and store. */
c->ti_end_min = ti_end_min;
c->ti_end_max = ti_end_max;
- c->ti_old = ti_current;
+ c->ti_old_gpart = ti_current;
if (timer) TIMER_TOC(timer_dorecv_spart);
@@ -1727,15 +1744,18 @@ 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 &&
t->type != task_type_send && t->type != task_type_recv &&
- t->type != task_type_kick1 && t->type != task_type_drift)
+ t->type != task_type_kick1 && t->type != task_type_drift_part &&
+ t->type != task_type_drift_gpart)
error(
"Task (type='%s/%s') should have been skipped ti_current=%lld "
"c->ti_end_min=%lld",
@@ -1863,8 +1883,11 @@ void *runner_main(void *data) {
runner_do_extra_ghost(r, ci, 1);
break;
#endif
- case task_type_drift:
- runner_do_drift_particles(r, ci, 1);
+ case task_type_drift_part:
+ runner_do_drift_part(r, ci, 1);
+ break;
+ case task_type_drift_gpart:
+ runner_do_drift_gpart(r, ci, 1);
break;
case task_type_kick1:
runner_do_kick1(r, ci, 1);
@@ -1902,14 +1925,11 @@ void *runner_main(void *data) {
}
break;
#endif
- case task_type_grav_mm:
- // runner_do_grav_mm(r, t->ci, 1);
- break;
case task_type_grav_down:
runner_do_grav_down(r, t->ci, 1);
break;
case task_type_grav_top_level:
- // runner_do_grav_top_level(r);
+ runner_do_grav_fft(r, 1);
break;
case task_type_grav_long_range:
runner_do_grav_long_range(r, t->ci, 1);
diff --git a/src/runner.h b/src/runner.h
index 5cb08cc14222c074b229c9bdf96ac47a072cd34e..0c6edc3c0c1406855ac79c96617bbdaa310bb46d 100644
--- a/src/runner.h
+++ b/src/runner.h
@@ -62,13 +62,15 @@ struct runner {
void runner_do_ghost(struct runner *r, struct cell *c, int timer);
void runner_do_extra_ghost(struct runner *r, struct cell *c, int timer);
void runner_do_sort(struct runner *r, struct cell *c, int flag, int clock);
-void runner_do_drift_particles(struct runner *r, struct cell *c, int timer);
+void runner_do_drift_part(struct runner *r, struct cell *c, int timer);
+void runner_do_drift_gpart(struct runner *r, struct cell *c, int timer);
void runner_do_kick1(struct runner *r, struct cell *c, int timer);
void runner_do_kick2(struct runner *r, struct cell *c, int timer);
void runner_do_end_force(struct runner *r, struct cell *c, int timer);
void runner_do_init(struct runner *r, struct cell *c, int timer);
void runner_do_cooling(struct runner *r, struct cell *c, int timer);
void runner_do_grav_external(struct runner *r, struct cell *c, int timer);
+void runner_do_grav_fft(struct runner *r, int timer);
void *runner_main(void *data);
void runner_do_unskip_mapper(void *map_data, int num_elements,
void *extra_data);
diff --git a/src/runner_doiact.h b/src/runner_doiact.h
index 839e492956140bd2b633d1e7c99c2b2a43f89629..f4513ca75b994c51a74011f0da4ee6863a625968 100644
--- a/src/runner_doiact.h
+++ b/src/runner_doiact.h
@@ -903,7 +903,7 @@ void DOPAIR1(struct runner *r, struct cell *ci, struct cell *cj) {
/* Anything to do here? */
if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
- if (!cell_is_drifted(ci, e) || !cell_is_drifted(cj, e))
+ if (!cell_are_part_drifted(ci, e) || !cell_are_part_drifted(cj, e))
error("Interacting undrifted cells.");
/* Get the sort ID. */
@@ -1147,7 +1147,7 @@ void DOPAIR2(struct runner *r, struct cell *ci, struct cell *cj) {
/* Anything to do here? */
if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
- if (!cell_is_drifted(ci, e) || !cell_is_drifted(cj, e))
+ if (!cell_are_part_drifted(ci, e) || !cell_are_part_drifted(cj, e))
error("Interacting undrifted cells.");
/* Get the shift ID. */
@@ -1597,7 +1597,7 @@ void DOSELF1(struct runner *r, struct cell *restrict c) {
if (!cell_is_active(c, e)) return;
- if (!cell_is_drifted(c, e)) error("Interacting undrifted cell.");
+ if (!cell_are_part_drifted(c, e)) error("Interacting undrifted cell.");
struct part *restrict parts = c->parts;
const int count = c->count;
@@ -1846,7 +1846,7 @@ void DOSELF2(struct runner *r, struct cell *restrict c) {
if (!cell_is_active(c, e)) return;
- if (!cell_is_drifted(c, e)) error("Cell is not drifted");
+ if (!cell_are_part_drifted(c, e)) error("Cell is not drifted");
struct part *restrict parts = c->parts;
const int count = c->count;
@@ -2279,8 +2279,8 @@ void DOSUB_PAIR1(struct runner *r, struct cell *ci, struct cell *cj, int sid,
else if (cell_is_active(ci, e) || cell_is_active(cj, e)) {
/* Make sure both cells are drifted to the current timestep. */
- if (!cell_is_drifted(ci, e)) cell_drift_particles(ci, e);
- if (!cell_is_drifted(cj, e)) cell_drift_particles(cj, e);
+ if (!cell_are_part_drifted(ci, e)) cell_drift_part(ci, e);
+ if (!cell_are_part_drifted(cj, e)) cell_drift_part(cj, e);
/* Do any of the cells need to be sorted first? */
if (!(ci->sorted & (1 << sid)) ||
@@ -2333,7 +2333,7 @@ void DOSUB_SELF1(struct runner *r, struct cell *ci, int gettimer) {
else {
/* Drift the cell to the current timestep if needed. */
- if (!cell_is_drifted(ci, r->e)) cell_drift_particles(ci, r->e);
+ if (!cell_are_part_drifted(ci, r->e)) cell_drift_part(ci, r->e);
#if (DOSELF1 == runner_doself1_density) && defined(WITH_VECTORIZATION) && \
defined(GADGET2_SPH)
@@ -2586,8 +2586,8 @@ void DOSUB_PAIR2(struct runner *r, struct cell *ci, struct cell *cj, int sid,
else if (cell_is_active(ci, e) || cell_is_active(cj, e)) {
/* Make sure both cells are drifted to the current timestep. */
- if (!cell_is_drifted(ci, e)) cell_drift_particles(ci, e);
- if (!cell_is_drifted(cj, e)) cell_drift_particles(cj, e);
+ if (!cell_are_part_drifted(ci, e)) cell_drift_part(ci, e);
+ if (!cell_are_part_drifted(cj, e)) cell_drift_part(cj, e);
/* Do any of the cells need to be sorted first? */
if (!(ci->sorted & (1 << sid)) ||
@@ -3218,7 +3218,7 @@ void DOSUB_SUBSET(struct runner *r, struct cell *ci, struct part *parts,
new_sid = sortlistID[new_sid];
/* Do any of the cells need to be drifted first? */
- if (!cell_is_drifted(cj, e)) cell_drift_particles(cj, e);
+ if (!cell_are_part_drifted(cj, e)) cell_drift_part(cj, e);
DOPAIR_SUBSET(r, ci, parts, ind, count, cj);
}
diff --git a/src/runner_doiact_fft.c b/src/runner_doiact_fft.c
index 076ec2578361127266c637cc5ac224609b702c66..01ccb3446dda35135c01a13cd213c8b6368d79ec 100644
--- a/src/runner_doiact_fft.c
+++ b/src/runner_doiact_fft.c
@@ -31,6 +31,283 @@
#include "runner_doiact_fft.h"
/* Local includes. */
+#include "engine.h"
+#include "error.h"
#include "runner.h"
+#include "space.h"
+#include "timers.h"
-void runner_do_grav_fft(struct runner *r) {}
+/**
+ * @brief Returns 1D index of a 3D NxNxN array using row-major style.
+ *
+ * @param i Index along x.
+ * @param j Index along y.
+ * @param k Index along z.
+ * @param N Size of the array along one axis.
+ */
+__attribute__((always_inline)) INLINE static int row_major_id(int i, int j,
+ int k, int N) {
+ return ((i % N) * N * N + (j % N) * N + (k % N));
+}
+
+/**
+ * @brief Assigns a given multipole to a density mesh using the CIC method.
+ *
+ * @param m The #multipole.
+ * @param rho The density mesh.
+ * @param N the size of the mesh along one axis.
+ * @param fac The width of a mesh cell.
+ */
+__attribute__((always_inline)) INLINE static void multipole_to_mesh_CIC(
+ const struct gravity_tensors* m, double* rho, int N, double fac) {
+
+ int i = (int)(fac * m->CoM[0]);
+ if (i >= N) i = N - 1;
+ const double dx = fac * m->CoM[0] - i;
+ const double tx = 1. - dx;
+
+ int j = (int)(fac * m->CoM[1]);
+ if (j >= N) j = N - 1;
+ const double dy = fac * m->CoM[1] - j;
+ const double ty = 1. - dy;
+
+ int k = (int)(fac * m->CoM[2]);
+ if (k >= N) k = N - 1;
+ 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;
+ rho[row_major_id(i + 0, j + 1, k + 0, N)] += m->m_pole.M_000 * tx * dy * tz;
+ rho[row_major_id(i + 0, j + 1, k + 1, N)] += m->m_pole.M_000 * tx * dy * dz;
+ rho[row_major_id(i + 1, j + 0, k + 0, N)] += m->m_pole.M_000 * dx * ty * tz;
+ rho[row_major_id(i + 1, j + 0, k + 1, N)] += m->m_pole.M_000 * dx * ty * dz;
+ rho[row_major_id(i + 1, j + 1, k + 0, N)] += m->m_pole.M_000 * dx * dy * tz;
+ rho[row_major_id(i + 1, j + 1, k + 1, N)] += m->m_pole.M_000 * dx * dy * dz;
+}
+
+/**
+ * @brief Computes the potential on a multipole from a given mesh using the CIC
+ * method.
+ *
+ * @param m The #multipole.
+ * @param pot The potential mesh.
+ * @param N the size of the mesh along one axis.
+ * @param fac width of a mesh cell.
+ */
+__attribute__((always_inline)) INLINE static void mesh_to_multipole_CIC(
+ struct gravity_tensors* m, double* pot, int N, double fac) {
+
+ int i = (int)(fac * m->CoM[0]);
+ if (i >= N) i = N - 1;
+ const double dx = fac * m->CoM[0] - i;
+ const double tx = 1. - dx;
+
+ int j = (int)(fac * m->CoM[1]);
+ if (j >= N) j = N - 1;
+ const double dy = fac * m->CoM[1] - j;
+ const double ty = 1. - dy;
+
+ int k = (int)(fac * m->CoM[2]);
+ if (k >= N) k = N - 1;
+ 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;
+ m->pot.F_000 += pot[row_major_id(i + 0, j + 1, k + 0, N)] * tx * dy * tz;
+ m->pot.F_000 += pot[row_major_id(i + 0, j + 1, k + 1, N)] * tx * dy * dz;
+ m->pot.F_000 += pot[row_major_id(i + 1, j + 0, k + 0, N)] * dx * ty * tz;
+ m->pot.F_000 += pot[row_major_id(i + 1, j + 0, k + 1, N)] * dx * ty * dz;
+ m->pot.F_000 += pot[row_major_id(i + 1, j + 1, k + 0, N)] * dx * dy * tz;
+ m->pot.F_000 += pot[row_major_id(i + 1, j + 1, k + 1, N)] * dx * dy * dz;
+}
+
+/**
+ * @brief Computes the potential on the top multipoles using a Fourier transform
+ *
+ * @param r The #runner task
+ * @param timer Are we timing this ?
+ */
+void runner_do_grav_fft(struct runner* r, int timer) {
+
+#ifdef HAVE_FFTW
+
+ const struct engine* e = r->e;
+ const struct space* s = e->s;
+ const integertime_t ti_current = e->ti_current;
+ const double a_smooth = e->gravity_properties->a_smooth;
+ const double box_size = s->dim[0];
+ const int cdim[3] = {s->cdim[0], s->cdim[1], s->cdim[2]};
+
+ TIMER_TIC;
+
+ if (cdim[0] != cdim[1] || cdim[0] != cdim[2]) error("Non-square mesh");
+
+ /* Some useful constants */
+ const int N = cdim[0];
+ const int N_half = N / 2;
+ const double cell_fac = N / box_size;
+
+ /* Recover the list of top-level multipoles */
+ const int nr_cells = s->nr_cells;
+ struct gravity_tensors* restrict multipoles = s->multipoles_top;
+ struct cell* cells = s->cells_top;
+
+ /* Make sure everything has been drifted to the current point */
+ for (int i = 0; i < nr_cells; ++i)
+ if (cells[i].ti_old_multipole != ti_current)
+ cell_drift_multipole(&cells[i], e);
+ // error("Top-level multipole %d not drifted", i);
+
+ /* Allocates some memory for the density mesh */
+ 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 */
+ fftw_complex* restrict frho = fftw_alloc_complex(N * N * (N_half + 1));
+ if (frho == NULL)
+ error("Error allocating memory for transform of density mesh");
+
+ /* Prepare the FFT library */
+ fftw_plan forward_plan = fftw_plan_dft_r2c_3d(
+ N, N, N, rho, frho, FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
+ fftw_plan inverse_plan = fftw_plan_dft_c2r_3d(
+ N, N, N, frho, rho, FFTW_ESTIMATE | FFTW_DESTROY_INPUT);
+
+ /* Do a CIC mesh assignment of the multipoles */
+ 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);
+
+ /* Fourier transform to go to magic-land */
+ fftw_execute(forward_plan);
+
+ /* frho now contains the Fourier transform of the density field */
+ /* frho contains NxNx(N/2+1) complex numbers */
+
+ /* Some common factors */
+ const double green_fac = -1. / (M_PI * box_size);
+ const double a_smooth2 = 4. * M_PI * a_smooth * a_smooth / ((double)(N * N));
+ const double k_fac = M_PI / (double)N;
+
+ /* Now de-convolve the CIC kernel and apply the Green function */
+ for (int i = 0; i < N; ++i) {
+
+ /* kx component of vector in Fourier space and 1/sinc(kx) */
+ const int kx = (i > N_half ? i - N : i);
+ const double kx_d = (double)kx;
+ const double fx = k_fac * kx_d;
+ const double sinc_kx_inv = (kx != 0) ? fx / sin(fx) : 1.;
+
+ for (int j = 0; j < N; ++j) {
+
+ /* ky component of vector in Fourier space and 1/sinc(ky) */
+ const int ky = (j > N_half ? j - N : j);
+ const double ky_d = (double)ky;
+ const double fy = k_fac * ky_d;
+ const double sinc_ky_inv = (ky != 0) ? fy / sin(fy) : 1.;
+
+ for (int k = 0; k < N_half + 1; ++k) {
+
+ /* kz component of vector in Fourier space and 1/sinc(kz) */
+ const int kz = (k > N_half ? k - N : k);
+ const double kz_d = (double)kz;
+ const double fz = k_fac * kz_d;
+ const double sinc_kz_inv = (kz != 0) ? fz / sin(fz) : 1.;
+
+ /* Norm of vector in Fourier space */
+ const double k2 = (kx_d * kx_d + ky_d * ky_d + kz_d * kz_d);
+
+ /* Avoid FPEs... */
+ if (k2 == 0.) continue;
+
+ /* Green function */
+ const double green_cor = green_fac * exp(-k2 * a_smooth2) / k2;
+
+ /* Deconvolution of CIC */
+ const double CIC_cor = sinc_kx_inv * sinc_ky_inv * sinc_kz_inv;
+ const double CIC_cor2 = CIC_cor * CIC_cor;
+ const double CIC_cor4 = CIC_cor2 * CIC_cor2;
+
+ /* Combined correction */
+ const double total_cor = green_cor * CIC_cor4;
+
+ /* Apply to the mesh */
+ const int index = N * (N_half + 1) * i + (N_half + 1) * j + k;
+ frho[index][0] *= total_cor;
+ frho[index][1] *= total_cor;
+ }
+ }
+ }
+
+ /* Correct singularity at (0,0,0) */
+ frho[0][0] = 0.;
+ frho[0][1] = 0.;
+
+ /* Fourier transform to come back from magic-land */
+ fftw_execute(inverse_plan);
+
+ /* rho now contains the potential */
+ /* This array is now again NxNxN real numbers */
+
+ /* Get the potential from the mesh using CIC */
+ for (int i = 0; i < nr_cells; ++i)
+ mesh_to_multipole_CIC(&multipoles[i], rho, N, cell_fac);
+
+ /* Clean-up the mess */
+ fftw_destroy_plan(forward_plan);
+ fftw_destroy_plan(inverse_plan);
+ fftw_free(rho);
+ fftw_free(frho);
+
+ /* Time the whole thing */
+ if (timer) TIMER_TOC(timer_dograv_top_level);
+
+#else
+ error("No FFTW library found. Cannot compute periodic long-range forces.");
+#endif
+}
+
+#ifdef HAVE_FFTW
+void print_array(double* array, int N) {
+
+ for (int k = N - 1; k >= 0; --k) {
+ printf("--- z = %d ---------\n", k);
+ for (int j = N - 1; j >= 0; --j) {
+ for (int i = 0; i < N; ++i) {
+ printf("%f ", array[i * N * N + j * N + k]);
+ }
+ printf("\n");
+ }
+ }
+}
+
+void print_carray(fftw_complex* array, int N) {
+
+ for (int k = N - 1; k >= 0; --k) {
+ printf("--- z = %d ---------\n", k);
+ for (int j = N - 1; j >= 0; --j) {
+ for (int i = 0; i < N; ++i) {
+ printf("(%f %f) ", array[i * N * N + j * N + k][0],
+ array[i * N * N + j * N + k][1]);
+ }
+ printf("\n");
+ }
+ }
+}
+#endif /* HAVE_FFTW */
diff --git a/src/runner_doiact_fft.h b/src/runner_doiact_fft.h
index 263662383fb465dcf945e55494a569289b009ff9..e9836311e71803952969b9c9316e8c81676d2dd8 100644
--- a/src/runner_doiact_fft.h
+++ b/src/runner_doiact_fft.h
@@ -21,6 +21,6 @@
struct runner;
-void runner_do_grav_fft(struct runner *r);
+void runner_do_grav_fft(struct runner *r, int timer);
#endif /* SWIFT_RUNNER_DOIACT_FFT_H */
diff --git a/src/runner_doiact_grav.h b/src/runner_doiact_grav.h
index 13a55344d773e7fba000d680eae9866dffdd88e1..a66cc5e0c9ed241aba3bb1b4329016b8e505e280 100644
--- a/src/runner_doiact_grav.h
+++ b/src/runner_doiact_grav.h
@@ -182,8 +182,8 @@ void runner_dopair_grav_pp(struct runner *r, struct cell *ci, struct cell *cj) {
if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
/* Let's start by drifting things */
- if (!cell_is_drifted(ci, e)) cell_drift_particles(ci, e);
- if (!cell_is_drifted(cj, e)) cell_drift_particles(cj, e);
+ if (!cell_are_gpart_drifted(ci, e)) cell_drift_gpart(ci, e);
+ if (!cell_are_gpart_drifted(cj, e)) cell_drift_gpart(cj, e);
#if ICHECK > 0
for (int pid = 0; pid < gcount_i; pid++) {
@@ -318,7 +318,7 @@ void runner_doself_grav_pp(struct runner *r, struct cell *c) {
if (!cell_is_active(c, e)) return;
/* Do we need to start by drifting things ? */
- if (!cell_is_drifted(c, e)) cell_drift_particles(c, e);
+ if (!cell_are_gpart_drifted(c, e)) cell_drift_gpart(c, e);
#if ICHECK > 0
for (int pid = 0; pid < gcount; pid++) {
@@ -429,6 +429,11 @@ void runner_dopair_grav(struct runner *r, struct cell *ci, struct cell *cj,
/* Sanity check */
if (ci == cj) error("Pair interaction between a cell and itself.");
+
+ if (cell_is_active(ci, e) && ci->ti_old_multipole != e->ti_current)
+ error("ci->multipole not drifted.");
+ if (cell_is_active(cj, e) && cj->ti_old_multipole != e->ti_current)
+ error("cj->multipole not drifted.");
#endif
#if ICHECK > 0
diff --git a/src/runner_doiact_vec.c b/src/runner_doiact_vec.c
index 9f9db24d8e27e6ab64afc75fcd207b67aa34597c..a302817f7fec4089f1046eec5cf7ff09aadd25a6 100644
--- a/src/runner_doiact_vec.c
+++ b/src/runner_doiact_vec.c
@@ -381,7 +381,7 @@ __attribute__((always_inline)) INLINE void runner_doself1_density_vec(
if (!cell_is_active(c, e)) return;
- if (!cell_is_drifted(c, e)) error("Interacting undrifted cell.");
+ if (!cell_are_part_drifted(c, e)) error("Interacting undrifted cell.");
/* Get the particle cache from the runner and re-allocate
* the cache if it is not big enough for the cell. */
@@ -625,7 +625,7 @@ void runner_dopair1_density_vec(struct runner *r, struct cell *ci,
/* Anything to do here? */
if (!cell_is_active(ci, e) && !cell_is_active(cj, e)) return;
- if (!cell_is_drifted(ci, e) || !cell_is_drifted(cj, e))
+ if (!cell_are_part_drifted(ci, e) || !cell_are_part_drifted(cj, e))
error("Interacting undrifted cells.");
/* Get the sort ID. */
diff --git a/src/scheduler.c b/src/scheduler.c
index 9cf1598f6d3bfd31b1120058d85b18edf44c427f..7e42c3a214cff3fe30e7c885b06c48d25eac0e8d 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -112,22 +112,12 @@ void scheduler_addunlock(struct scheduler *s, struct task *ta,
}
/**
- * @brief Split a task if too large.
+ * @brief Split a hydrodynamic task if too large.
*
* @param t The #task
* @param s The #scheduler we are working in.
*/
-static void scheduler_splittask(struct task *t, struct scheduler *s) {
-
- /* Static constants. */
- static const int pts[7][8] = {
- {-1, 12, 10, 9, 4, 3, 1, 0}, {-1, -1, 11, 10, 5, 4, 2, 1},
- {-1, -1, -1, 12, 7, 6, 4, 3}, {-1, -1, -1, -1, 8, 7, 5, 4},
- {-1, -1, -1, -1, -1, 12, 10, 9}, {-1, -1, -1, -1, -1, -1, 11, 10},
- {-1, -1, -1, -1, -1, -1, -1, 12}};
- static const float sid_scale[13] = {
- 0.1897f, 0.4025f, 0.1897f, 0.4025f, 0.5788f, 0.4025f, 0.1897f,
- 0.4025f, 0.1897f, 0.4025f, 0.5788f, 0.4025f, 0.5788f};
+static void scheduler_splittask_hydro(struct task *t, struct scheduler *s) {
/* Iterate on this task until we're done with it. */
int redo = 1;
@@ -137,11 +127,7 @@ 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)) ||
- ((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) ||
- ((t->type == task_type_timestep) && t->ci->nodeID != s->nodeID)) {
+ if ((t->ci == NULL) || (t->type == task_type_pair && t->cj == NULL)) {
t->type = task_type_none;
t->subtype = task_subtype_none;
t->cj = NULL;
@@ -154,7 +140,7 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
/* Get a handle on the cell involved. */
struct cell *ci = t->ci;
- const double hi = ci->dmin;
+ const double width = ci->dmin;
/* Foreign task? */
if (ci->nodeID != s->nodeID) {
@@ -162,27 +148,16 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
break;
}
- /* Is this cell even split? */
- if (ci->split && 2.f * kernel_gamma * ci->h_max * space_stretch < hi) {
+ /* Is this cell even split and the task does not violate h ? */
+ if (ci->split && 2.f * kernel_gamma * ci->h_max * space_stretch < width) {
/* Make a sub? */
if (scheduler_dosub && /* Note division here to avoid overflow */
- ((ci->count > 0 && ci->count < space_subsize / ci->count) ||
- (ci->gcount > 0 && ci->gcount < space_subsize / ci->gcount))) {
+ (ci->count > 0 && ci->count < space_subsize / ci->count)) {
/* convert to a self-subtask. */
t->type = task_type_sub_self;
- /* Make sure we have a drift task (MATTHIEU temp. fix for gravity) */
- if (t->subtype == task_subtype_grav ||
- t->subtype == task_subtype_external_grav) {
- lock_lock(&ci->lock);
- if (ci->drift == NULL)
- ci->drift = scheduler_addtask(s, task_type_drift,
- task_subtype_none, 0, 0, ci, NULL);
- lock_unlock_blind(&ci->lock);
- }
-
/* Depend on local sorts on this cell. */
if (ci->sorts != NULL) scheduler_addunlock(s, ci->sorts, t);
@@ -198,44 +173,42 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
t->ci = ci->progeny[first_child];
for (int k = first_child + 1; k < 8; k++)
if (ci->progeny[k] != NULL)
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_self, t->subtype, 0, 0,
ci->progeny[k], NULL),
s);
- /* Make a task for each pair of progeny unless it's ext. gravity. */
- 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(
- scheduler_addtask(s, task_type_pair, t->subtype,
- pts[j][k], 0, ci->progeny[j],
- ci->progeny[k]),
- s);
- }
+ /* Make a task for each pair of progeny */
+ 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_hydro(
+ 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. */
+ /* Otherwise, make sure the self task has a drift task */
else {
+
lock_lock(&ci->lock);
- if (ci->drift == NULL)
- ci->drift = scheduler_addtask(s, task_type_drift, task_subtype_none,
- 0, 0, ci, NULL);
+
+ if (ci->drift_part == NULL)
+ ci->drift_part = scheduler_addtask(s, task_type_drift_part,
+ task_subtype_none, 0, 0, ci, NULL);
lock_unlock_blind(&ci->lock);
}
+ } /* Self interaction */
- /* Pair interaction? */
- } else if (t->type == task_type_pair && t->subtype != task_subtype_grav) {
+ /* Pair interaction? */
+ else if (t->type == task_type_pair) {
/* Get a handle on the cells involved. */
struct cell *ci = t->ci;
struct cell *cj = t->cj;
- const double hi = ci->dmin;
- const double hj = cj->dmin;
/* Foreign task? */
if (ci->nodeID != s->nodeID && cj->nodeID != s->nodeID) {
@@ -248,10 +221,13 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
double shift[3];
const int sid = space_getsid(s->space, &ci, &cj, shift);
+ const double width_i = ci->dmin;
+ const double width_j = cj->dmin;
+
/* Should this task be split-up? */
if (ci->split && cj->split &&
- 2.f * kernel_gamma * space_stretch * ci->h_max < hi &&
- 2.f * kernel_gamma * space_stretch * cj->h_max < hj) {
+ 2.f * kernel_gamma * space_stretch * ci->h_max < width_i &&
+ 2.f * kernel_gamma * space_stretch * cj->h_max < width_j) {
/* Replace by a single sub-task? */
if (scheduler_dosub &&
@@ -284,15 +260,15 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
t->ci = ci->progeny[6];
t->cj = cj->progeny[0];
t->flags = 1;
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 1, 0,
ci->progeny[7], cj->progeny[1]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 0, 0,
ci->progeny[6], cj->progeny[1]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 2, 0,
ci->progeny[7], cj->progeny[0]),
s);
@@ -308,15 +284,15 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
t->ci = ci->progeny[5];
t->cj = cj->progeny[0];
t->flags = 3;
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 3, 0,
ci->progeny[7], cj->progeny[2]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 0, 0,
ci->progeny[5], cj->progeny[2]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 6, 0,
ci->progeny[7], cj->progeny[0]),
s);
@@ -326,63 +302,63 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
t->ci = ci->progeny[4];
t->cj = cj->progeny[0];
t->flags = 4;
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 5, 0,
ci->progeny[5], cj->progeny[0]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 7, 0,
ci->progeny[6], cj->progeny[0]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 8, 0,
ci->progeny[7], cj->progeny[0]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 3, 0,
ci->progeny[4], cj->progeny[1]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 4, 0,
ci->progeny[5], cj->progeny[1]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 6, 0,
ci->progeny[6], cj->progeny[1]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 7, 0,
ci->progeny[7], cj->progeny[1]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 1, 0,
ci->progeny[4], cj->progeny[2]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 2, 0,
ci->progeny[5], cj->progeny[2]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 4, 0,
ci->progeny[6], cj->progeny[2]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 5, 0,
ci->progeny[7], cj->progeny[2]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 0, 0,
ci->progeny[4], cj->progeny[3]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 1, 0,
ci->progeny[5], cj->progeny[3]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 3, 0,
ci->progeny[6], cj->progeny[3]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 4, 0,
ci->progeny[7], cj->progeny[3]),
s);
@@ -392,15 +368,15 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
t->ci = ci->progeny[4];
t->cj = cj->progeny[1];
t->flags = 5;
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 5, 0,
ci->progeny[6], cj->progeny[3]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 2, 0,
ci->progeny[4], cj->progeny[3]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 8, 0,
ci->progeny[6], cj->progeny[1]),
s);
@@ -416,15 +392,15 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
t->ci = ci->progeny[4];
t->cj = cj->progeny[3];
t->flags = 6;
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 8, 0,
ci->progeny[5], cj->progeny[2]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 7, 0,
ci->progeny[4], cj->progeny[2]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 7, 0,
ci->progeny[5], cj->progeny[3]),
s);
@@ -440,15 +416,15 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
t->ci = ci->progeny[3];
t->cj = cj->progeny[0];
t->flags = 9;
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 9, 0,
ci->progeny[7], cj->progeny[4]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 0, 0,
ci->progeny[3], cj->progeny[4]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 8, 0,
ci->progeny[7], cj->progeny[0]),
s);
@@ -458,63 +434,63 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
t->ci = ci->progeny[2];
t->cj = cj->progeny[0];
t->flags = 10;
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 11, 0,
ci->progeny[3], cj->progeny[0]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 7, 0,
ci->progeny[6], cj->progeny[0]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 6, 0,
ci->progeny[7], cj->progeny[0]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 9, 0,
ci->progeny[2], cj->progeny[1]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 10, 0,
ci->progeny[3], cj->progeny[1]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 8, 0,
ci->progeny[6], cj->progeny[1]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 7, 0,
ci->progeny[7], cj->progeny[1]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 1, 0,
ci->progeny[2], cj->progeny[4]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 2, 0,
ci->progeny[3], cj->progeny[4]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 10, 0,
ci->progeny[6], cj->progeny[4]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 11, 0,
ci->progeny[7], cj->progeny[4]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 0, 0,
ci->progeny[2], cj->progeny[5]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 1, 0,
ci->progeny[3], cj->progeny[5]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 9, 0,
ci->progeny[6], cj->progeny[5]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 10, 0,
ci->progeny[7], cj->progeny[5]),
s);
@@ -524,15 +500,15 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
t->ci = ci->progeny[2];
t->cj = cj->progeny[1];
t->flags = 11;
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 11, 0,
ci->progeny[6], cj->progeny[5]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 2, 0,
ci->progeny[2], cj->progeny[5]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 6, 0,
ci->progeny[6], cj->progeny[1]),
s);
@@ -542,63 +518,63 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
t->ci = ci->progeny[1];
t->cj = cj->progeny[0];
t->flags = 12;
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 11, 0,
ci->progeny[3], cj->progeny[0]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 5, 0,
ci->progeny[5], cj->progeny[0]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 2, 0,
ci->progeny[7], cj->progeny[0]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 9, 0,
ci->progeny[1], cj->progeny[2]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 12, 0,
ci->progeny[3], cj->progeny[2]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 8, 0,
ci->progeny[5], cj->progeny[2]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 5, 0,
ci->progeny[7], cj->progeny[2]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 3, 0,
ci->progeny[1], cj->progeny[4]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 6, 0,
ci->progeny[3], cj->progeny[4]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 12, 0,
ci->progeny[5], cj->progeny[4]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 11, 0,
ci->progeny[7], cj->progeny[4]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 0, 0,
ci->progeny[1], cj->progeny[6]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 3, 0,
ci->progeny[3], cj->progeny[6]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 9, 0,
ci->progeny[5], cj->progeny[6]),
s);
- scheduler_splittask(
+ scheduler_splittask_hydro(
scheduler_addtask(s, task_type_pair, t->subtype, 12, 0,
ci->progeny[7], cj->progeny[6]),
s);
@@ -623,7 +599,7 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
struct task *tl =
scheduler_addtask(s, task_type_pair, t->subtype, 0, 0,
ci->progeny[j], cj->progeny[k]);
- scheduler_splittask(tl, s);
+ scheduler_splittask_hydro(tl, s);
tl->flags = space_getsid(s->space, &t->ci, &t->cj, shift);
}
@@ -632,9 +608,9 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
/* Create the drift and sort for ci. */
lock_lock(&ci->lock);
- if (ci->drift == NULL && ci->nodeID == engine_rank)
- ci->drift = scheduler_addtask(s, task_type_drift, task_subtype_none,
- 0, 0, ci, NULL);
+ if (ci->drift_part == NULL && ci->nodeID == engine_rank)
+ ci->drift_part = scheduler_addtask(s, task_type_drift_part,
+ task_subtype_none, 0, 0, ci, NULL);
if (ci->sorts == NULL)
ci->sorts = scheduler_addtask(s, task_type_sort, task_subtype_none,
1 << sid, 0, ci, NULL);
@@ -643,11 +619,11 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
lock_unlock_blind(&ci->lock);
scheduler_addunlock(s, ci->sorts, t);
- /* Create the sort for cj. */
+ /* Create the drift and sort for cj. */
lock_lock(&cj->lock);
- if (cj->drift == NULL && cj->nodeID == engine_rank)
- cj->drift = scheduler_addtask(s, task_type_drift, task_subtype_none,
- 0, 0, cj, NULL);
+ if (cj->drift_part == NULL && cj->nodeID == engine_rank)
+ cj->drift_part = scheduler_addtask(s, task_type_drift_part,
+ task_subtype_none, 0, 0, cj, NULL);
if (cj->sorts == NULL)
cj->sorts = scheduler_addtask(s, task_type_sort, task_subtype_none,
1 << sid, 0, cj, NULL);
@@ -656,19 +632,142 @@ static void scheduler_splittask(struct task *t, struct scheduler *s) {
lock_unlock_blind(&cj->lock);
scheduler_addunlock(s, cj->sorts, t);
}
-
} /* pair interaction? */
+ } /* iterate over the current task. */
+}
+
+/**
+ * @brief Split a gravity task if too large.
+ *
+ * @param t The #task
+ * @param s The #scheduler we are working in.
+ */
+static void scheduler_splittask_gravity(struct task *t, struct scheduler *s) {
+
+ /* Iterate on this task until we're done with it. */
+ int redo = 1;
+ while (redo) {
- /* Long-range gravity interaction ? */
- else if (t->type == task_type_grav_mm) {
+ /* Reset the redo flag. */
+ redo = 0;
+
+ /* Non-splittable task? */
+ if ((t->ci == NULL) || (t->type == task_type_pair && t->cj == NULL)) {
+ t->type = task_type_none;
+ t->subtype = task_subtype_none;
+ t->cj = NULL;
+ t->skip = 1;
+ break;
+ }
+
+ /* Self-interaction? */
+ if (t->type == task_type_self) {
+
+ /* Get a handle on the cell involved. */
+ struct cell *ci = t->ci;
+
+ /* Foreign task? */
+ if (ci->nodeID != s->nodeID) {
+ t->skip = 1;
+ break;
+ }
+
+ /* Is this cell even split? */
+ if (ci->split) {
+
+ /* Make a sub? */
+ if (scheduler_dosub && /* Note division here to avoid overflow */
+ (ci->gcount > 0 && ci->gcount < space_subsize / ci->gcount)) {
+
+ /* 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);
+ }
+ }
+ } /* Cell is split */
+
+ /* Otherwise, make sure the self task has a drift task */
+ else {
+
+ 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);
+ }
+ } /* Self interaction */
+
+ /* Pair interaction? */
+ else if (t->type == task_type_pair) {
/* Get a handle on the cells involved. */
struct cell *ci = t->ci;
+ struct cell *cj = t->cj;
+
+ /* Foreign task? */
+ if (ci->nodeID != s->nodeID && cj->nodeID != s->nodeID) {
+ t->skip = 1;
+ break;
+ }
+
+ /* Should this task be split-up? */
+ if (ci->split && cj->split) {
- /* Safety thing */
- if (ci->gcount == 0) t->type = task_type_none;
+ // MATTHIEU: nothing here for now
- } /* gravity interaction? */
+ } else {
+
+ /* Create the drift for ci. */
+ lock_lock(&ci->lock);
+ if (ci->drift_gpart == NULL && ci->nodeID == engine_rank)
+ ci->drift_gpart = scheduler_addtask(
+ s, task_type_drift_gpart, task_subtype_none, 0, 0, ci, NULL);
+ lock_unlock_blind(&ci->lock);
+
+ /* Create the drift for cj. */
+ lock_lock(&cj->lock);
+ if (cj->drift_gpart == NULL && cj->nodeID == engine_rank)
+ cj->drift_gpart = scheduler_addtask(
+ s, task_type_drift_gpart, task_subtype_none, 0, 0, cj, NULL);
+ lock_unlock_blind(&cj->lock);
+ }
+ } /* pair interaction? */
} /* iterate over the current task. */
}
@@ -688,7 +787,20 @@ void scheduler_splittasks_mapper(void *map_data, int num_elements,
for (int ind = 0; ind < num_elements; ind++) {
struct task *t = &tasks[ind];
- scheduler_splittask(t, s);
+
+ /* Invoke the correct splitting strategy */
+ if (t->subtype == task_subtype_density) {
+ scheduler_splittask_hydro(t, s);
+ } else if (t->subtype == task_subtype_external_grav) {
+ scheduler_splittask_gravity(t, s);
+ } else if (t->subtype == task_subtype_grav) {
+ scheduler_splittask_gravity(t, s);
+ } else if (t->type == task_type_grav_top_level ||
+ t->type == task_type_grav_ghost) {
+ // MATTHIEU: for the future
+ } else {
+ error("Unexpected task sub-type");
+ }
}
}
@@ -780,7 +892,8 @@ void scheduler_set_unlocks(struct scheduler *s) {
/* Check that we are not overflowing */
if (counts[s->unlock_ind[k]] < 0)
error("Task unlocking more than %d other tasks!",
- (1 << (sizeof(short int) - 1)) - 1);
+ (1 << (8 * sizeof(short int) - 1)) - 1);
+
#endif
}
@@ -963,9 +1076,6 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
int *tid = s->tasks_ind;
struct task *tasks = s->tasks;
const int nodeID = s->nodeID;
- const float sid_scale[13] = {0.1897, 0.4025, 0.1897, 0.4025, 0.5788,
- 0.4025, 0.1897, 0.4025, 0.1897, 0.4025,
- 0.5788, 0.4025, 0.5788};
const float wscale = 0.001;
const ticks tic = getticks();
@@ -1012,9 +1122,12 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
case task_type_ghost:
if (t->ci == t->ci->super) cost = wscale * t->ci->count;
break;
- case task_type_drift:
+ case task_type_drift_part:
cost = wscale * t->ci->count;
break;
+ case task_type_drift_gpart:
+ cost = wscale * t->ci->gcount;
+ break;
case task_type_kick1:
cost = wscale * t->ci->count;
break;
@@ -1142,6 +1255,11 @@ 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 ||
+ t->type == task_type_grav_ghost)
+ continue;
+
if (ci == NULL && cj == NULL) {
error("Task not associated with cells!");
@@ -1149,7 +1267,8 @@ void scheduler_start(struct scheduler *s) {
} else if (cj == NULL) { /* self */
if (ci->ti_end_min == ti_current && t->skip &&
- t->type != task_type_sort && t->type != task_type_drift && t->type)
+ t->type != task_type_sort && t->type != task_type_drift_part &&
+ t->type != task_type_drift_gpart)
error(
"Task (type='%s/%s') should not have been skipped "
"ti_current=%lld "
@@ -1240,7 +1359,8 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
case task_type_ghost:
case task_type_kick1:
case task_type_kick2:
- case task_type_drift:
+ case task_type_drift_part:
+ case task_type_drift_gpart:
case task_type_timestep:
qid = t->ci->super->owner;
break;
diff --git a/src/sort_part.h b/src/sort_part.h
index a243fcdfae8ec0aba606000e26bc18d35601215c..74116d7a8cada31c0663d5c5b70cfa978b11af8b 100644
--- a/src/sort_part.h
+++ b/src/sort_part.h
@@ -83,6 +83,18 @@ static const int sortlistID[27] = {
/* ( 1 , 1 , 0 ) */ 1,
/* ( 1 , 1 , 1 ) */ 0};
+/* Ratio of particles interacting assuming a uniform distribution */
+static const float sid_scale[13] = {0.1897f, 0.4025f, 0.1897f, 0.4025f, 0.5788f,
+ 0.4025f, 0.1897f, 0.4025f, 0.1897f, 0.4025f,
+ 0.5788f, 0.4025f, 0.5788f};
+
+/* Sid flags for every sub-pair of a self task. */
+static const int sub_sid_flag[7][8] = {
+ {-1, 12, 10, 9, 4, 3, 1, 0}, {-1, -1, 11, 10, 5, 4, 2, 1},
+ {-1, -1, -1, 12, 7, 6, 4, 3}, {-1, -1, -1, -1, 8, 7, 5, 4},
+ {-1, -1, -1, -1, -1, 12, 10, 9}, {-1, -1, -1, -1, -1, -1, 11, 10},
+ {-1, -1, -1, -1, -1, -1, -1, 12}};
+
/**
* @brief Determines whether a pair of cells are corner to corner.
*
diff --git a/src/space.c b/src/space.c
index 0c67fe27ac1c7b2cae3672e7f0b1ce82be5fb804..677219a0a201ca357d0e2e24fa9cc1b61baa3e31 100644
--- a/src/space.c
+++ b/src/space.c
@@ -205,7 +205,8 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->gradient = NULL;
c->force = NULL;
c->grav = NULL;
- c->dx_max = 0.0f;
+ c->dx_max_part = 0.0f;
+ c->dx_max_gpart = 0.0f;
c->dx_max_sort = 0.0f;
c->sorted = 0;
c->count = 0;
@@ -217,10 +218,12 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->kick1 = NULL;
c->kick2 = NULL;
c->timestep = NULL;
- c->drift = NULL;
+ c->drift_part = NULL;
+ c->drift_gpart = NULL;
c->cooling = NULL;
c->sourceterms = NULL;
- c->grav_top_level = NULL;
+ c->grav_ghost[0] = NULL;
+ c->grav_ghost[1] = NULL;
c->grav_long_range = NULL;
c->grav_down = NULL;
c->super = c;
@@ -420,7 +423,8 @@ void space_regrid(struct space *s, int verbose) {
c->gcount = 0;
c->scount = 0;
c->super = c;
- c->ti_old = ti_old;
+ c->ti_old_part = ti_old;
+ c->ti_old_gpart = ti_old;
c->ti_old_multipole = ti_old;
if (s->gravity) c->multipole = &s->multipoles_top[cid];
}
@@ -890,8 +894,9 @@ void space_rebuild(struct space *s, int verbose) {
#ifdef SWIFT_DEBUG_CHECKS
/* Verify that the links are correct */
- part_verify_links(s->parts, s->gparts, s->sparts, nr_parts, nr_gparts,
- nr_sparts, verbose);
+ if ((nr_gparts > 0 && nr_parts > 0) || (nr_gparts > 0 && nr_sparts > 0))
+ part_verify_links(s->parts, s->gparts, s->sparts, nr_parts, nr_gparts,
+ nr_sparts, verbose);
#endif
/* Hook the cells up to the parts. */
@@ -902,7 +907,8 @@ void space_rebuild(struct space *s, int verbose) {
struct spart *sfinger = s->sparts;
for (int k = 0; k < s->nr_cells; k++) {
struct cell *restrict c = &cells_top[k];
- c->ti_old = ti_old;
+ c->ti_old_part = ti_old;
+ c->ti_old_gpart = ti_old;
c->ti_old_multipole = ti_old;
c->parts = finger;
c->xparts = xfinger;
@@ -2011,7 +2017,8 @@ void space_split_recursive(struct space *s, struct cell *c,
cp->count = 0;
cp->gcount = 0;
cp->scount = 0;
- cp->ti_old = c->ti_old;
+ cp->ti_old_part = c->ti_old_part;
+ cp->ti_old_gpart = c->ti_old_gpart;
cp->ti_old_multipole = c->ti_old_multipole;
cp->loc[0] = c->loc[0];
cp->loc[1] = c->loc[1];
@@ -2025,8 +2032,9 @@ void space_split_recursive(struct space *s, struct cell *c,
if (k & 1) cp->loc[2] += cp->width[2];
cp->depth = c->depth + 1;
cp->split = 0;
- cp->h_max = 0.0;
- cp->dx_max = 0.f;
+ cp->h_max = 0.f;
+ cp->dx_max_part = 0.f;
+ cp->dx_max_gpart = 0.f;
cp->dx_max_sort = 0.f;
cp->nodeID = c->nodeID;
cp->parent = c;
@@ -2877,7 +2885,8 @@ void space_check_drift_point(struct space *s, integertime_t ti_drift,
int multipole) {
#ifdef SWIFT_DEBUG_CHECKS
/* Recursively check all cells */
- space_map_cells_pre(s, 1, cell_check_particle_drift_point, &ti_drift);
+ space_map_cells_pre(s, 1, cell_check_part_drift_point, &ti_drift);
+ space_map_cells_pre(s, 1, cell_check_gpart_drift_point, &ti_drift);
if (multipole)
space_map_cells_pre(s, 1, cell_check_multipole_drift_point, &ti_drift);
#else
diff --git a/src/space.h b/src/space.h
index c5f588563e5a9fb4b6cb73ac1446514f8149794f..69742eaa3319d6371006014ee2cab1ac5f10bfe7 100644
--- a/src/space.h
+++ b/src/space.h
@@ -130,6 +130,9 @@ struct space {
/*! The s-particle data (cells have pointers to this). */
struct spart *sparts;
+ /*! The top-level FFT task */
+ struct task *grav_top_level;
+
/*! General-purpose lock for this space. */
swift_lock_type lock;
diff --git a/src/swift.h b/src/swift.h
index 7f1b19b6066c2d55df1cb9101172ae94c9085583..20397eb24df478cba65a0e35d686b402f1d8ee70 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -45,6 +45,7 @@
#include "parser.h"
#include "part.h"
#include "partition.h"
+#include "periodic.h"
#include "physical_constants.h"
#include "potential.h"
#include "profiler.h"
diff --git a/src/task.c b/src/task.c
index e8c35e49a57595a6415c60ce7071ae1c0e3f09b7..43da1d35680783d977ea743dd4f43c52f0f291bc 100644
--- a/src/task.c
+++ b/src/task.c
@@ -47,27 +47,15 @@
#include "lock.h"
/* Task type names. */
-const char *taskID_names[task_type_count] = {"none",
- "sort",
- "self",
- "pair",
- "sub_self",
- "sub_pair",
- "init_grav",
- "ghost",
- "extra_ghost",
- "drift",
- "kick1",
- "kick2",
- "timestep",
- "send",
- "recv",
- "grav_top_level",
- "grav_long_range",
- "grav_mm",
- "grav_down",
- "cooling",
- "sourceterms"};
+const char *taskID_names[task_type_count] = {
+ "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_ghost", "grav_mm", "grav_down",
+ "cooling", "sourceterms"};
/* Sub-task type names. */
const char *subtaskID_names[task_subtype_count] = {
@@ -132,6 +120,7 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
return task_action_none;
break;
+ case task_type_drift_part:
case task_type_sort:
case task_type_ghost:
case task_type_extra_ghost:
@@ -169,7 +158,6 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
case task_type_timestep:
case task_type_send:
case task_type_recv:
- case task_type_drift:
if (t->ci->count > 0 && t->ci->gcount > 0)
return task_action_all;
else if (t->ci->count > 0)
@@ -187,8 +175,10 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
return task_action_multipole;
break;
+ case task_type_drift_gpart:
case task_type_grav_down:
return task_action_gpart;
+ break;
default:
error("Unknown task_action for task");
@@ -286,15 +276,19 @@ void task_unlock(struct task *t) {
case task_type_kick1:
case task_type_kick2:
case task_type_timestep:
- case task_type_drift:
cell_unlocktree(ci);
cell_gunlocktree(ci);
break;
+ case task_type_drift_part:
case task_type_sort:
cell_unlocktree(ci);
break;
+ case task_type_drift_gpart:
+ cell_gunlocktree(ci);
+ break;
+
case task_type_self:
case task_type_sub_self:
if (subtype == task_subtype_grav) {
@@ -323,7 +317,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);
@@ -372,7 +365,6 @@ int task_lock(struct task *t) {
case task_type_kick1:
case task_type_kick2:
case task_type_timestep:
- case task_type_drift:
if (ci->hold || ci->ghold) return 0;
if (cell_locktree(ci) != 0) return 0;
if (cell_glocktree(ci) != 0) {
@@ -381,10 +373,17 @@ int task_lock(struct task *t) {
}
break;
+ case task_type_drift_part:
case task_type_sort:
+ if (ci->hold) return 0;
if (cell_locktree(ci) != 0) return 0;
break;
+ case task_type_drift_gpart:
+ if (ci->ghold) return 0;
+ if (cell_glocktree(ci) != 0) return 0;
+ break;
+
case task_type_self:
case task_type_sub_self:
if (subtype == task_subtype_grav) {
@@ -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 */
diff --git a/src/task.h b/src/task.h
index 049f86bdd6b4baf0856745b2b53acda5cca8c9e1..052f3e8036381441e283d3f7847d09e98ec1dac2 100644
--- a/src/task.h
+++ b/src/task.h
@@ -47,7 +47,8 @@ enum task_types {
task_type_init_grav,
task_type_ghost,
task_type_extra_ghost,
- task_type_drift,
+ task_type_drift_part,
+ task_type_drift_gpart,
task_type_kick1,
task_type_kick2,
task_type_timestep,
@@ -55,6 +56,7 @@ enum task_types {
task_type_recv,
task_type_grav_top_level,
task_type_grav_long_range,
+ task_type_grav_ghost,
task_type_grav_mm,
task_type_grav_down,
task_type_cooling,
diff --git a/src/timers.c b/src/timers.c
index aa42eee14fc0df3edd5a18340c092b8eea2ffac1..62eac20596a082e411ced61a86f32bef9edcb636 100644
--- a/src/timers.c
+++ b/src/timers.c
@@ -40,7 +40,8 @@ const char* timers_names[timer_count] = {
"prepare",
"init",
"init_grav",
- "drift",
+ "drift_part",
+ "drift_gpart",
"kick1",
"kick2",
"timestep",
@@ -58,6 +59,7 @@ const char* timers_names[timer_count] = {
"dopair_grav_pp",
"dograv_external",
"dograv_down",
+ "dograv_top_level",
"dograv_long_range",
"dosource",
"dosub_self_density",
diff --git a/src/timers.h b/src/timers.h
index 08e983a947bc57d9dcc7a432df92c2a4b0a1f7d7..9248be4f3048e468deed476f822947eed3c4ce56 100644
--- a/src/timers.h
+++ b/src/timers.h
@@ -41,7 +41,8 @@ enum {
timer_prepare,
timer_init,
timer_init_grav,
- timer_drift,
+ timer_drift_part,
+ timer_drift_gpart,
timer_kick1,
timer_kick2,
timer_timestep,
@@ -59,6 +60,7 @@ enum {
timer_dopair_grav_pp,
timer_dograv_external,
timer_dograv_down,
+ timer_dograv_top_level,
timer_dograv_long_range,
timer_dosource,
timer_dosub_self_density,
diff --git a/tests/Makefile.am b/tests/Makefile.am
index a51b8eb82a17313818ff956ca3f15a232df8df65..09142a3744755b3d380235a0030d62dca487b7c9 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -94,4 +94,5 @@ EXTRA_DIST = testReading.sh makeInput.py testPair.sh testPairPerturbed.sh \
test27cells.sh test27cellsPerturbed.sh testParser.sh \
test125cells.sh testParserInput.yaml difffloat.py \
tolerance_125.dat tolerance_27_normal.dat tolerance_27_perturbed.dat \
- tolerance_pair_normal.dat tolerance_pair_perturbed.dat
+ tolerance_pair_normal.dat tolerance_pair_perturbed.dat \
+ fft_params.yml
diff --git a/tests/fft_params.yml b/tests/fft_params.yml
new file mode 100644
index 0000000000000000000000000000000000000000..05d6d8f0b0578d11645fc1d78c1a6322160ae87a
--- /dev/null
+++ b/tests/fft_params.yml
@@ -0,0 +1,10 @@
+Scheduler:
+ max_top_level_cells: 64
+
+# 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.00001 # Softening length (in internal units).
+ a_smooth: 0.
+ r_cut: 0.
diff --git a/tests/test125cells.c b/tests/test125cells.c
index 168b4838eab5b27f359ab927a7bae2240919e82f..91023ae76b70f8e997a162bed7402c3515673f74 100644
--- a/tests/test125cells.c
+++ b/tests/test125cells.c
@@ -1,3 +1,4 @@
+
/*******************************************************************************
* This file is part of SWIFT.
* Copyright (C) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
@@ -315,7 +316,7 @@ struct cell *make_cell(size_t n, const double offset[3], double size, double h,
cell->h_max = h;
cell->count = count;
cell->gcount = 0;
- cell->dx_max = 0.;
+ cell->dx_max_part = 0.;
cell->dx_max_sort = 0.;
cell->width[0] = size;
cell->width[1] = size;
@@ -324,7 +325,7 @@ struct cell *make_cell(size_t n, const double offset[3], double size, double h,
cell->loc[1] = offset[1];
cell->loc[2] = offset[2];
- cell->ti_old = 8;
+ cell->ti_old_part = 8;
cell->ti_end_min = 8;
cell->ti_end_max = 8;
cell->ti_sort = 0;
@@ -525,9 +526,9 @@ int main(int argc, char *argv[]) {
/* Build the infrastructure */
struct space space;
space.periodic = 1;
- space.dim[0] = 3.;
- space.dim[1] = 3.;
- space.dim[2] = 3.;
+ space.dim[0] = 5.;
+ space.dim[1] = 5.;
+ space.dim[2] = 5.;
hydro_space_init(&space.hs, &space);
struct phys_const prog_const;
@@ -592,8 +593,7 @@ int main(int argc, char *argv[]) {
const ticks tic = getticks();
/* Initialise the particles */
- for (int j = 0; j < 125; ++j)
- runner_do_drift_particles(&runner, cells[j], 0);
+ for (int j = 0; j < 125; ++j) runner_do_drift_part(&runner, cells[j], 0);
/* First, sort stuff */
for (int j = 0; j < 125; ++j) runner_do_sort(&runner, cells[j], 0x1FFF, 0);
@@ -670,6 +670,12 @@ int main(int argc, char *argv[]) {
outputFileNameExtension);
dump_particle_fields(outputFileName, main_cell, solution, 0);
}
+
+ /* Reset stuff */
+ for (int i = 0; i < 125; ++i) {
+ for (int n = 0; n < cells[i]->count; ++n)
+ hydro_init_part(&cells[i]->parts[n], &space.hs);
+ }
}
/* Output timing */
diff --git a/tests/test27cells.c b/tests/test27cells.c
index bd827b68e90ea5f4e9d5577612e6cecda2edf83a..a0f541d17100a13079580aabbef065fa5adbd5e1 100644
--- a/tests/test27cells.c
+++ b/tests/test27cells.c
@@ -158,7 +158,7 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
cell->split = 0;
cell->h_max = h;
cell->count = count;
- cell->dx_max = 0.;
+ cell->dx_max_part = 0.;
cell->dx_max_sort = 0.;
cell->width[0] = size;
cell->width[1] = size;
@@ -167,7 +167,7 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
cell->loc[1] = offset[1];
cell->loc[2] = offset[2];
- cell->ti_old = 8;
+ cell->ti_old_part = 8;
cell->ti_end_min = 8;
cell->ti_end_max = 8;
cell->ti_sort = 8;
@@ -438,7 +438,7 @@ int main(int argc, char *argv[]) {
cells[i * 9 + j * 3 + k] = make_cell(particles, offset, size, h, rho,
&partId, perturbation, vel);
- runner_do_drift_particles(&runner, cells[i * 9 + j * 3 + k], 0);
+ runner_do_drift_part(&runner, cells[i * 9 + j * 3 + k], 0);
runner_do_sort(&runner, cells[i * 9 + j * 3 + k], 0x1FFF, 0);
}
diff --git a/tests/testFFT.c b/tests/testFFT.c
index c4aeb2885c788bd769bda49bdd15ab121dd8e9d4..37b1b0ca00097593991dc551ae4dc5e90a11a3b4 100644
--- a/tests/testFFT.c
+++ b/tests/testFFT.c
@@ -18,8 +18,8 @@
******************************************************************************/
/* Some standard headers. */
-#include <stdlib.h>
-#include <string.h>
+
+#include "../config.h"
#ifndef HAVE_FFTW
@@ -27,169 +27,97 @@ int main() { return 0; }
#else
-#include <fftw3.h>
+/* Some standard headers. */
+#include <fenv.h>
+#include <stdlib.h>
+#include <string.h>
/* Includes. */
#include "swift.h"
-const double G = 1.;
-
-const size_t N = 16;
-const size_t PMGRID = 8;
-
-// const double asmth = 2. * M_PI * const_gravity_a_smooth / boxSize;
-// const double asmth2 = asmth * asmth;
-// const double fact = G / (M_PI * boxSize) * (1. / (2. * boxSize / PMGRID));
-
int main() {
/* Initialize CPU frequency, this also starts time. */
unsigned long long cpufreq = 0;
clocks_set_cpufreq(cpufreq);
- /* Simulation properties */
- const size_t count = N * N * N;
- const double boxSize = 1.;
-
- /* Create some particles */
- struct gpart* gparts = malloc(count * sizeof(struct gpart));
- bzero(gparts, count * sizeof(struct gpart));
- for (size_t i = 0; i < N; ++i) {
- for (size_t j = 0; j < N; ++j) {
- for (size_t k = 0; k < N; ++k) {
-
- struct gpart* gp = &gparts[i * N * N + j * N + k];
-
- gp->x[0] = i * boxSize / N + boxSize / (2 * N);
- gp->x[1] = j * boxSize / N + boxSize / (2 * N);
- gp->x[2] = k * boxSize / N + boxSize / (2 * N);
-
- gp->mass = 1. / count;
-
- gp->id_or_neg_offset = i * N * N + j * N + k;
- }
- }
- }
-
- /* Properties of the mesh */
- const size_t meshmin[3] = {0, 0, 0};
- const size_t meshmax[3] = {PMGRID - 1, PMGRID - 1, PMGRID - 1};
-
- const size_t dimx = meshmax[0] - meshmin[0] + 2;
- const size_t dimy = meshmax[1] - meshmin[1] + 2;
- const size_t dimz = meshmax[2] - meshmin[2] + 2;
-
- const double fac = PMGRID / boxSize;
- const size_t PMGRID2 = 2 * (PMGRID / 2 + 1);
-
- /* message("dimx=%zd dimy=%zd dimz=%zd", dimx, dimy, dimz); */
-
- /* Allocate and empty the workspace mesh */
- const size_t workspace_size = (dimx + 4) * (dimy + 4) * (dimz + 4);
- double* workspace = fftw_malloc(workspace_size * sizeof(double));
- bzero(workspace, workspace_size * sizeof(double));
-
- /* Do CIC with the particles */
- for (size_t pid = 0; pid < count; ++pid) {
-
- const struct gpart* const gp = &gparts[pid];
-
- const size_t slab_x =
- (fac * gp->x[0] >= PMGRID) ? PMGRID - 1 : fac * gp->x[0];
- const size_t slab_y =
- (fac * gp->x[1] >= PMGRID) ? PMGRID - 1 : fac * gp->x[1];
- const size_t slab_z =
- (fac * gp->x[2] >= PMGRID) ? PMGRID - 1 : fac * gp->x[2];
-
- const double dx = fac * gp->x[0] - (double)slab_x;
- const double dy = fac * gp->x[1] - (double)slab_y;
- const double dz = fac * gp->x[2] - (double)slab_z;
-
- const size_t slab_xx = slab_x + 1;
- const size_t slab_yy = slab_y + 1;
- const size_t slab_zz = slab_z + 1;
-
- workspace[(slab_x * dimy + slab_y) * dimz + slab_z] +=
- gp->mass * (1.0 - dx) * (1.0 - dy) * (1.0 - dz);
- workspace[(slab_x * dimy + slab_yy) * dimz + slab_z] +=
- gp->mass * (1.0 - dx) * dy * (1.0 - dz);
- workspace[(slab_x * dimy + slab_y) * dimz + slab_zz] +=
- gp->mass * (1.0 - dx) * (1.0 - dy) * dz;
- workspace[(slab_x * dimy + slab_yy) * dimz + slab_zz] +=
- gp->mass * (1.0 - dx) * dy * dz;
- workspace[(slab_xx * dimy + slab_y) * dimz + slab_z] +=
- gp->mass * (dx) * (1.0 - dy) * (1.0 - dz);
- workspace[(slab_xx * dimy + slab_yy) * dimz + slab_z] +=
- gp->mass * (dx)*dy * (1.0 - dz);
- workspace[(slab_xx * dimy + slab_y) * dimz + slab_zz] +=
- gp->mass * (dx) * (1.0 - dy) * dz;
- workspace[(slab_xx * dimy + slab_yy) * dimz + slab_zz] +=
- gp->mass * (dx)*dy * dz;
- }
-
- /* for(size_t i = 0 ; i < dimx*dimy*dimz; ++i) */
- /* message("workspace[%zd] = %f", i, workspace[i]); */
-
- /* Prepare the force grid */
- const size_t fft_size = workspace_size;
- double* forcegrid = fftw_malloc(fft_size * sizeof(double));
- bzero(forcegrid, fft_size * sizeof(double));
-
- const size_t sendmin = 0, recvmin = 0;
- const size_t sendmax = PMGRID, recvmax = PMGRID;
-
- memcpy(forcegrid, workspace + (sendmin - meshmin[0]) * dimy * dimz,
- (sendmax - sendmin + 1) * dimy * dimz * sizeof(double));
-
- /* for (size_t i = 0; i < fft_size; ++i) */
- /* if (forcegrid[i] != workspace[i]) error("wrong"); */
-
- /* Prepare the density grid */
- double* rhogrid = fftw_malloc(fft_size * sizeof(double));
- bzero(rhogrid, fft_size * sizeof(double));
-
- /* Now get the density */
- for (size_t slab_x = recvmin; slab_x <= recvmax; slab_x++) {
-
- const size_t slab_xx = slab_x % PMGRID;
-
- for (size_t slab_y = recvmin; slab_y <= recvmax; slab_y++) {
-
- const size_t slab_yy = slab_y % PMGRID;
-
- for (size_t slab_z = recvmin; slab_z <= recvmax; slab_z++) {
-
- const size_t slab_zz = slab_z % PMGRID;
-
- rhogrid[PMGRID * PMGRID2 * slab_xx + PMGRID2 * slab_yy + slab_zz] +=
- forcegrid[((slab_x - recvmin) * dimy + (slab_y - recvmin)) * dimz +
- (slab_z - recvmin)];
- }
- }
- }
-
- /* for (size_t i = 0; i < 640; i++) { */
- /* if (rhogrid[i] != workspace[i]) { */
- /* message("rhogrid[%zd]= %f workspace[%zd]= %f forcegrid[%zd]= %f", i, */
- /* rhogrid[i], i, workspace[i], i, forcegrid[i]); */
- /* } */
- /* } */
-
- /* FFT of the density field */
- fftw_complex* fftgrid = fftw_malloc(fft_size * sizeof(fftw_complex));
- fftw_plan plan_forward = fftw_plan_dft_r2c_3d(PMGRID, PMGRID, PMGRID, rhogrid,
- fftgrid, FFTW_ESTIMATE);
- fftw_execute(plan_forward);
-
- for (size_t i = 0; i < 640; i++) {
- message("workspace[%zd]= %f", i, fftgrid[i][0]);
+ /* Choke on FP-exceptions */
+ feenableexcept(FE_DIVBYZERO | FE_INVALID | FE_OVERFLOW);
+
+ /* Make one particle */
+ int nr_gparts = 1;
+ struct gpart *gparts = NULL;
+ if (posix_memalign((void **)&gparts, 64, nr_gparts * sizeof(struct gpart)) !=
+ 0)
+ error("Impossible to allocate memory for gparts.");
+ bzero(gparts, nr_gparts * sizeof(struct gpart));
+
+ gparts[0].x[0] = 0.3;
+ gparts[0].x[1] = 0.8;
+ gparts[0].x[2] = 0.2;
+ gparts[0].mass = 1.f;
+
+ /* Read the parameter file */
+ struct swift_params *params = malloc(sizeof(struct swift_params));
+ parser_read_file("fft_params.yml", params);
+
+ /* Initialise the gravity properties */
+ struct gravity_props gravity_properties;
+ gravity_props_init(&gravity_properties, params);
+
+ /* Build the infrastructure */
+ struct space space;
+ double dim[3] = {1., 1., 1.};
+ space_init(&space, params, dim, NULL, gparts, NULL, 0, nr_gparts, 0, 1, 1, 1,
+ 0, 0);
+
+ struct engine engine;
+ engine.s = &space;
+ space.e = &engine;
+ engine.time = 0.1f;
+ engine.ti_current = 0;
+ engine.ti_old = 0;
+ engine.max_active_bin = num_time_bins;
+ engine.gravity_properties = &gravity_properties;
+ engine.nr_threads = 1;
+
+ struct runner runner;
+ runner.e = &engine;
+
+ /* Initialize the threadpool. */
+ threadpool_init(&engine.threadpool, engine.nr_threads);
+
+ space_rebuild(&space, 0);
+
+ /* Run the FFT task */
+ runner_do_grav_fft(&runner, 1);
+
+ /* Now check that we got the right answer */
+ int nr_cells = space.nr_cells;
+ double *r = malloc(nr_cells * sizeof(double));
+ double *pot = malloc(nr_cells * sizeof(double));
+ double *pot_exact = malloc(nr_cells * sizeof(double));
+
+ // FILE *file = fopen("potential.dat", "w");
+ for (int i = 0; i < nr_cells; ++i) {
+ pot[i] = space.multipoles_top[i].pot.F_000;
+ double dx =
+ nearest(space.multipoles_top[i].CoM[0] - gparts[0].x[0], dim[0]);
+ double dy =
+ nearest(space.multipoles_top[i].CoM[1] - gparts[0].x[1], dim[1]);
+ double dz =
+ nearest(space.multipoles_top[i].CoM[2] - gparts[0].x[2], dim[2]);
+ r[i] = sqrt(dx * dx + dy * dy + dz * dz);
+ if (r[i] > 0) pot_exact[i] = -1. / r[i];
+ // fprintf(file, "%e %e %e\n", r[i], pot[i], pot_exact[i]);
}
+ // fclose(file);
- /* Clean-up */
- fftw_destroy_plan(plan_forward);
- fftw_free(forcegrid);
- fftw_free(rhogrid);
- fftw_free(workspace);
+ /* Clean up */
+ free(r);
+ free(pot);
+ free(pot_exact);
+ free(params);
free(gparts);
return 0;
}
diff --git a/tests/testPair.c b/tests/testPair.c
index c2533b63b902e3bdc7e7cae6fcbcf50c87dee4af..92987d2fdb625fec6e186a280837f145787f599b 100644
--- a/tests/testPair.c
+++ b/tests/testPair.c
@@ -84,7 +84,8 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
cell->split = 0;
cell->h_max = h;
cell->count = count;
- cell->dx_max = 0.;
+ cell->dx_max_part = 0.;
+ cell->dx_max_sort = 0.;
cell->width[0] = n;
cell->width[1] = n;
cell->width[2] = n;
@@ -92,7 +93,7 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
cell->loc[1] = offset[1];
cell->loc[2] = offset[2];
- cell->ti_old = 8;
+ cell->ti_old_part = 8;
cell->ti_end_min = 8;
cell->ti_end_max = 8;
diff --git a/tests/testSPHStep.c b/tests/testSPHStep.c
index 0c7ae1d0d8855371b8f8f9fbf51c7c63b3221aaa..014dacd1eb62040b03e6038b2c23183a24ec4850 100644
--- a/tests/testSPHStep.c
+++ b/tests/testSPHStep.c
@@ -71,7 +71,8 @@ struct cell *make_cell(size_t N, float cellSize, int offset[3], int id_offset) {
cell->h_max = h;
cell->count = count;
cell->gcount = 0;
- cell->dx_max = 0.;
+ cell->dx_max_part = 0.;
+ cell->dx_max_sort = 0.;
cell->width[0] = cellSize;
cell->width[1] = cellSize;
cell->width[2] = cellSize;