diff --git a/.gitignore b/.gitignore
index d69de3deae4cbe1dff23fd0f22704f6b4f5e40ed..14c768367ee9181c1459b57b93dd92d431ace981 100644
--- a/.gitignore
+++ b/.gitignore
@@ -29,6 +29,10 @@ tests/testGreetings
tests/testReading
tests/input.hdf5
tests/testSingle
+tests/testTimeIntegration
+tests/testSPHStep
+
+theory/latex/swift.pdf
m4/libtool.m4
m4/ltoptions.m4
diff --git a/src/engine.c b/src/engine.c
index a5b2a842f38232bdc75d55d51ad0b8a97c21083f..bdf430646a7632da5b1a847b6b22710eae115e4b 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -362,7 +362,7 @@ void engine_repartition(struct engine *e) {
if (t->type != task_type_self && t->type != task_type_pair &&
t->type != task_type_sub && t->type != task_type_ghost &&
t->type != task_type_drift && t->type != task_type_kick &&
- t->type != task_type_init)
+ t->type != task_type_init)
continue;
/* Get the task weight. */
@@ -1370,16 +1370,45 @@ int engine_marktasks(struct engine *e) {
}
/**
- * @brief Rebuild the space and tasks.
+ * @brief Prints the number of tasks in the engine
*
* @param e The #engine.
*/
-void engine_rebuild(struct engine *e) {
+void engine_print(struct engine *e) {
int k;
struct scheduler *sched = &e->sched;
+ /* Count and print the number of each task type. */
+ int counts[task_type_count + 1];
+ for (k = 0; k <= task_type_count; k++) counts[k] = 0;
+ for (k = 0; k < sched->nr_tasks; k++)
+ if (!sched->tasks[k].skip)
+ counts[(int)sched->tasks[k].type] += 1;
+ else
+ counts[task_type_count] += 1;
+#ifdef WITH_MPI
+ printf("[%03i] engine_print: task counts are [ %s=%i", e->nodeID,
+ taskID_names[0], counts[0]);
+#else
+ printf("engine_print: task counts are [ %s=%i", taskID_names[0], counts[0]);
+#endif
+ for (k = 1; k < task_type_count; k++)
+ printf(" %s=%i", taskID_names[k], counts[k]);
+ printf(" skipped=%i ]\n", counts[task_type_count]);
+ fflush(stdout);
+ message("nr_parts = %i.", e->s->nr_parts);
+}
+
+/**
+ * @brief Rebuild the space and tasks.
+ *
+ * @param e The #engine.
+ */
+
+void engine_rebuild(struct engine *e) {
+
/* Clear the forcerebuild flag, whatever it was. */
e->forcerebuild = 0;
@@ -1410,25 +1439,8 @@ void engine_rebuild(struct engine *e) {
// message( "engine_marktasks took %.3f ms." , (double)(getticks() -
// tic)/CPU_TPS*1000 );
- /* Count and print the number of each task type. */
- int counts[task_type_count + 1];
- for (k = 0; k <= task_type_count; k++) counts[k] = 0;
- for (k = 0; k < sched->nr_tasks; k++)
- if (!sched->tasks[k].skip)
- counts[(int)sched->tasks[k].type] += 1;
- else
- counts[task_type_count] += 1;
-#ifdef WITH_MPI
- printf("[%03i] engine_rebuild: task counts are [ %s=%i", e->nodeID,
- taskID_names[0], counts[0]);
-#else
- printf("engine_rebuild: task counts are [ %s=%i", taskID_names[0], counts[0]);
-#endif
- for (k = 1; k < task_type_count; k++)
- printf(" %s=%i", taskID_names[k], counts[k]);
- printf(" skipped=%i ]\n", counts[task_type_count]);
- fflush(stdout);
- message("nr_parts = %i.", e->s->nr_parts);
+ /* Print the status of the system */
+ engine_print(e);
}
/**
@@ -1746,9 +1758,9 @@ void engine_step(struct engine *e) {
struct cell *c;
struct space *s = e->s;
- TIMER_TIC2
-
+ TIMER_TIC2;
+ message("Begin step: %d", e->step);
/* Re-distribute the particles amongst the nodes? */
if (e->forcerepart) engine_repartition(e);
@@ -1756,33 +1768,20 @@ void engine_step(struct engine *e) {
/* Prepare the space. */
engine_prepare(e);
- // engine_single_density( e->s->dim , 3392063069037 , e->s->parts ,
- // e->s->nr_parts , e->s->periodic );
+ engine_print(e);
/* Send off the runners. */
- TIMER_TIC
- engine_launch(
- e, e->nr_threads,
- (1 << task_type_sort) | (1 << task_type_self) | (1 << task_type_pair) |
- (1 << task_type_sub) | (1 << task_type_ghost) |
- (1 << task_type_kick) | (1 << task_type_send) |
- (1 << task_type_recv) |
- /* (1 << task_type_grav_pp) | (1 << task_type_grav_mm) | */
- /* (1 << task_type_grav_up) | (1 << task_type_grav_down) | */
- (1 << task_type_link));
+ TIMER_TIC;
+ engine_launch(e, e->nr_threads,
+ (1 << task_type_sort) | (1 << task_type_self) |
+ (1 << task_type_pair) | (1 << task_type_sub) |
+ (1 << task_type_init) | (1 << task_type_ghost) |
+ (1 << task_type_kick) | (1 << task_type_send) |
+ (1 << task_type_recv) | (1 << task_type_link));
TIMER_TOC(timer_runners);
- // engine_single_force( e->s->dim , 8328423931905 , e->s->parts ,
- // e->s->nr_parts , e->s->periodic );
-
- // for(k=0; k<10; ++k)
- // printParticle(parts, k);
- // printParticle( parts , 432626 );
- // printParticle( e->s->parts , 3392063069037 , e->s->nr_parts );
- // printParticle( e->s->parts , 8328423931905 , e->s->nr_parts );
-
- /* Collect the cell data from the second kick. */
+ /* Collect the cell data from the kick. */
for (k = 0; k < s->nr_cells; k++)
if (s->cells[k].nodeID == e->nodeID) {
c = &s->cells[k];
@@ -1806,12 +1805,12 @@ void engine_step(struct engine *e) {
out[0] = t_end_min;
if (MPI_Allreduce(out, in, 1, MPI_DOUBLE, MPI_MIN, MPI_COMM_WORLD) !=
MPI_SUCCESS)
- error("Failed to aggregate dt_min.");
+ error("Failed to aggregate t_end_min.");
t_end_min = in[0];
out[0] = t_end_max;
if (MPI_Allreduce(out, in, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD) !=
MPI_SUCCESS)
- error("Failed to aggregate dt_max.");
+ error("Failed to aggregate t_end_max.");
t_end_max = in[0];
out[0] = count;
out[1] = ekin;
@@ -1830,71 +1829,16 @@ if ( e->nodeID == 0 )
message( "nr_parts=%i." , nr_parts ); */
#endif
- e->count_step = count;
- e->ekin = ekin;
- e->epot = epot;
- // printParticle( e->s->parts , 382557 , e->s->nr_parts );
- // message( "dt_min/dt_max is %e/%e." , dt_min , dt_max ); fflush(stdout);
- // message( "etot is %e (ekin=%e, epot=%e)." , ekin+epot , ekin , epot );
- // fflush(stdout);
- // message( "total momentum is [ %e , %e , %e ]." , mom[0] , mom[1] , mom[2]
- // ); fflush(stdout);
- // message( "total angular momentum is [ %e , %e , %e ]." , ang[0] , ang[1] ,
- // ang[2] ); fflush(stdout);
- // message( "updated %i parts (dt_step=%.3e)." , count , dt_step );
- // fflush(stdout);
-
- /* Increase the step. */
+ /* Move forward in time */
+ e->timeOld = e->time;
+ e->time = t_end_min;
e->step += 1;
- /* /\* Does the time step need adjusting? *\/ */
- /* if (e->policy & engine_policy_fixdt) { */
- /* dt = e->dt_orig; */
- /* } else { */
- /* if (dt == 0) { */
- /* e->nullstep += 1; */
- /* if (e->dt_orig > 0.0) { */
- /* dt = e->dt_orig; */
- /* while (dt_min < dt) dt *= 0.5; */
- /* while (dt_min > 2 * dt) dt *= 2.0; */
- /* } else */
- /* dt = dt_min; */
- /* for (k = 0; k < s->nr_parts; k++) { */
- /* /\* struct part *p = &s->parts[k]; */
- /* struct xpart *xp = &s->xparts[k]; */
- /* float dt_curr = dt; */
- /* for ( int j = (int)( p->dt / dt ) ; j > 1 ; j >>= 1 ) */
- /* dt_curr *= 2.0f; */
- /* xp->dt_curr = dt_curr; *\/ */
- /* s->parts[k].dt = dt; */
- /* s->xparts[k].dt_curr = dt; */
- /* } */
- /* // message( "dt_min=%.3e, adjusting time step to dt=%e." , dt_min ,
- * e->dt */
- /* // ); */
- /* } else { */
- /* while (dt_min < dt) { */
- /* dt *= 0.5; */
- /* e->step *= 2; */
- /* e->nullstep *= 2; */
- /* // message( "dt_min dropped below time step, adjusting to dt=%e." ,
- */
- /* // e->dt ); */
- /* } */
- /* while (dt_min > 2 * dt && (e->step & 1) == 0) { */
- /* dt *= 2.0; */
- /* e->step /= 2; */
- /* e->nullstep /= 2; */
- /* // message( "dt_min is larger than twice the time step, adjusting to
- */
- /* // dt=%e." , e->dt ); */
- /* } */
- /* } */
- /* } */
- /* e->dt = dt; */
+ message("e->time=%f", e->time);
+ message("Ready to drift");
- /* /\* Set the system time. *\/ */
- /* e->time = dt * (e->step - e->nullstep); */
+ /* Drift all particles */
+ engine_launch(e, e->nr_threads, (1 << task_type_drift));
TIMER_TOC2(timer_step);
}
@@ -2067,7 +2011,6 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
float timeBegin, float timeEnd) {
int k;
- float dt_min = dt;
#if defined(HAVE_SETAFFINITY)
int nr_cores = sysconf(_SC_NPROCESSORS_ONLN);
int i, j, cpuid[nr_cores];
diff --git a/src/engine.h b/src/engine.h
index 6ef59767a358ce950a2222201e1087352e995ec8..afc2f9ad157b5895c84811d6b61a9e7288131b7b 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -137,6 +137,7 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
int nr_queues, int nr_nodes, int nodeID, int policy,
float timeBegin, float timeEnd);
void engine_prepare(struct engine *e);
+void engine_print(struct engine *e);
void engine_step(struct engine *e);
void engine_maketasks(struct engine *e);
void engine_split(struct engine *e, int *grid);
diff --git a/src/map.h b/src/map.h
index 0fe04ad0158499626a16a0b18637ff190cf59c97..0753c2641af6deb050c1dcef6bcd3ae4621ae6aa 100644
--- a/src/map.h
+++ b/src/map.h
@@ -18,7 +18,6 @@
*
******************************************************************************/
-
/* Includes. */
#ifndef SWIFT_MAP_H
#define SWIFT_MAP_H
diff --git a/src/parallel_io.c b/src/parallel_io.c
index 3fff596718ef41caba9c81a71e912b9d19ac192a..f48aaa04e59cb084f6d1a4af0a904a7476b470ac 100644
--- a/src/parallel_io.c
+++ b/src/parallel_io.c
@@ -441,7 +441,8 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
* @param N_total Total number of particles across all cores
* @param mpi_rank The MPI task rank calling the function
* @param offset Offset in the array where this mpi task starts writing
- * @param part A (char*) pointer on the first occurrence of the field of interest
+ * @param part A (char*) pointer on the first occurrence of the field of
+ *interest
*in the parts array
* @param field The name (code name) of the field to read from.
* @param us The UnitSystem currently in use
@@ -582,7 +583,8 @@ void write_output_parallel(struct engine* e, struct UnitSystem* us,
N_total, mpi_rank, offset, u, us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
writeArray(h_grp, fileName, xmfFile, "ParticleIDs", ULONGLONG, N, 1, parts,
N_total, mpi_rank, offset, id, us, UNIT_CONV_NO_UNITS);
- /* writeArray(h_grp, fileName, xmfFile, "TimeStep", FLOAT, N, 1, parts, N_total, */
+ /* writeArray(h_grp, fileName, xmfFile, "TimeStep", FLOAT, N, 1, parts,
+ * N_total, */
/* mpi_rank, offset, dt, us, UNIT_CONV_TIME); */
writeArray(h_grp, fileName, xmfFile, "Acceleration", FLOAT, N, 3, parts,
N_total, mpi_rank, offset, a, us, UNIT_CONV_ACCELERATION);
diff --git a/src/part.h b/src/part.h
index 97977dadc988f361c7c6aae5b512ca857ad384da..2b5a28a897a59ec8984cef7941e6f9a8fb34e8f6 100644
--- a/src/part.h
+++ b/src/part.h
@@ -40,9 +40,6 @@ struct xpart {
/* Old density. */
float omega;
- /* particle's current time-step. */
- float dt_curr;
-
} __attribute__((aligned(32)));
/* Gravity particle. */
@@ -115,46 +112,46 @@ struct part {
#endif
/* Store density/force specific stuff. */
- union {
+ // union {
- struct {
+ struct {
- /* Particle velocity divergence. */
- float div_v;
+ /* Particle velocity divergence. */
+ float div_v;
- /* Derivative of particle number density. */
- float wcount_dh;
+ /* Derivative of particle number density. */
+ float wcount_dh;
- /* Particle velocity curl. */
- float curl_v[3];
+ /* Particle velocity curl. */
+ float curl_v[3];
- /* Particle number density. */
- float wcount;
+ /* Particle number density. */
+ float wcount;
- } density;
+ } density;
- struct {
+ struct {
- /* Balsara switch */
- float balsara;
+ /* Balsara switch */
+ float balsara;
- /* Aggregate quantities. */
- float POrho2;
+ /* Aggregate quantities. */
+ float POrho2;
- /* Change in particle energy over time. */
- float u_dt;
+ /* Change in particle energy over time. */
+ float u_dt;
- /* Change in smoothing length over time. */
- float h_dt;
+ /* Change in smoothing length over time. */
+ float h_dt;
- /* Signal velocity */
- float v_sig;
+ /* Signal velocity */
+ float v_sig;
- /* Sound speed */
- float c;
+ /* Sound speed */
+ float c;
- } force;
- };
+ } force;
+ //};
/* Particle pressure. */
// float P;
diff --git a/src/runner.c b/src/runner.c
index 2e8a13dbcb024ba0ffc6ab9a3359f64e63d836ff..9544450cd1de0446b125a3e01a880cf67bb8eb04 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -518,8 +518,6 @@ void runner_doinit(struct runner *r, struct cell *c) {
/* Get a direct pointer on the part. */
p = &parts[i];
- if (p->id == 0) message("init!");
-
if (p->t_end <= t_end) {
/* Get ready for a density calculation */
@@ -593,8 +591,8 @@ void runner_doghost(struct runner *r, struct cell *c) {
/* Final operation on the density. */
p->rho = rho = ih * ih2 * (p->rho + p->mass * kernel_root);
p->rho_dh = rho_dh = (p->rho_dh - 3.0f * p->mass * kernel_root) * ih4;
- wcount = (p->density.wcount + kernel_root) *
- (4.0f / 3.0 * M_PI * kernel_gamma3);
+ p->density.wcount = wcount = (p->density.wcount + kernel_root) *
+ (4.0f / 3.0 * M_PI * kernel_gamma3);
wcount_dh =
p->density.wcount_dh * ih * (4.0f / 3.0 * M_PI * kernel_gamma3);
@@ -610,20 +608,12 @@ void runner_doghost(struct runner *r, struct cell *c) {
h_corr = fmaxf(h_corr, -h / 2.f);
}
- if (p->id == 0)
- message("ghost! h=%f dh=%f wcount=%f rho=%f", p->h, h_corr, wcount,
- p->rho);
-
- /* Apply the correction to p->h and to the compact part. */
- h = p->h += h_corr;
-
/* Did we get the right number density? */
if (wcount > kernel_nwneigh + const_delta_nwneigh ||
wcount < kernel_nwneigh - const_delta_nwneigh) {
- // message( "particle %lli (h=%e,depth=%i) has bad wcount=%.3f." ,
- // p->id , p->h , c->depth , wcount ); fflush(stdout);
- // p->h += ( p->density.wcount + kernel_root - kernel_nwneigh ) /
- // p->density.wcount_dh;
+ h = p->h += h_corr;
+ p->h += (p->density.wcount + kernel_root - kernel_nwneigh) /
+ p->density.wcount_dh;
pid[redo] = pid[i];
redo += 1;
p->density.wcount = 0.0;
@@ -672,7 +662,7 @@ void runner_doghost(struct runner *r, struct cell *c) {
(const_viscosity_alpha_max - p->alpha) * S;
/* Update particle's viscosity paramter */
- p->alpha += alpha_dot * 1; // p->dt;
+ p->alpha += alpha_dot * (p->t_end - p->t_begin);
#endif
/* Reset the acceleration. */
@@ -768,8 +758,6 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
p = &parts[k];
xp = &xparts[k];
- if (p->id == 0) message("drift !");
-
/* Get local copies of particle data. */
h = p->h;
ih = 1.0f / h;
@@ -816,11 +804,9 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
/* Predict gradient term */
p->force.POrho2 = u * (const_hydro_gamma - 1.0f) / (rho * xp->omega);
-
}
}
-
if (timer) {
#ifdef TIMER_VERBOSE
message("runner %02i: %i parts at depth %i took %.3f ms.", r->id, c->count,
@@ -830,11 +816,8 @@ void runner_dodrift(struct runner *r, struct cell *c, int timer) {
TIMER_TOC(timer_drift);
#endif
}
-
-
}
-
/**
* @brief Combined second and first kick for fixed dt.
*
@@ -879,15 +862,13 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
m = p->mass;
x[0] = p->x[0], x[1] = p->x[1], x[2] = p->x[2];
- if (p->id == 0)
- message("Kick ! t_beg=%f t_end=%f t_cur=%f", p->t_begin, p->t_end,
- t_current);
+ /* if (p->id == 0) */
+ /* message("Kick ! t_beg=%f t_end=%f t_cur=%f", p->t_begin, p->t_end, */
+ /* t_current); */
/* If particle needs to be kicked */
if (p->t_end <= t_current) {
- if (p->id == 0) message("Kicking like a boss");
-
/* First, finish the force loop */
p->force.h_dt *= p->h * 0.333333333f;
@@ -919,16 +900,10 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
t_end = p->t_end + 0.5f * new_dt;
dt = t_end - t_start;
- if (p->id == 0)
- message("Kick ! dt=%f t_beg=%f t_end=%f", dt, p->t_begin, p->t_end);
-
/* Move particle forward in time */
p->t_begin = p->t_end;
p->t_end = p->t_begin + dt;
- if (p->id == 0)
- message("Kick ! dt=%f t_beg=%f t_end=%f", dt, p->t_begin, p->t_end);
-
/* Kick particles in momentum space */
xp->v_full[0] += p->a[0] * dt;
xp->v_full[1] += p->a[1] * dt;
diff --git a/src/serial_io.c b/src/serial_io.c
index c61774fed2e5991e5a1ed06bac572224b09ade9c..0feda5adfd1ac702014e00fb0c5adae89d2e6868 100644
--- a/src/serial_io.c
+++ b/src/serial_io.c
@@ -314,7 +314,8 @@ void read_ic_serial(char* fileName, double dim[3], struct part** parts, int* N,
u, COMPULSORY);
readArray(h_grp, "ParticleIDs", ULONGLONG, *N, 1, *parts, N_total, offset,
id, COMPULSORY);
- /* readArray(h_grp, "TimeStep", FLOAT, *N, 1, *parts, N_total, offset, dt, */
+ /* readArray(h_grp, "TimeStep", FLOAT, *N, 1, *parts, N_total, offset, dt,
+ */
/* OPTIONAL); */
readArray(h_grp, "Acceleration", FLOAT, *N, 3, *parts, N_total, offset, a,
OPTIONAL);
@@ -487,7 +488,8 @@ void writeArrayBackEnd(hid_t grp, char* name, enum DATA_TYPE type, int N,
* @param type The #DATA_TYPE of the array.
* @param N The number of particles to write.
* @param dim The dimension of the data (1 for scalar, 3 for vector)
- * @param part A (char*) pointer on the first occurrence of the field of interest
+ * @param part A (char*) pointer on the first occurrence of the field of
+ *interest
*in the parts array
* @param field The name (code name) of the field to read from.
* @param us The UnitSystem currently in use
@@ -627,7 +629,8 @@ void write_output_serial(struct engine* e, struct UnitSystem* us, int mpi_rank,
us, UNIT_CONV_ENERGY_PER_UNIT_MASS);
prepareArray(h_grp, fileName, xmfFile, "ParticleIDs", ULONGLONG, N_total, 1,
us, UNIT_CONV_NO_UNITS);
- /* prepareArray(h_grp, fileName, xmfFile, "TimeStep", FLOAT, N_total, 1, us, */
+ /* prepareArray(h_grp, fileName, xmfFile, "TimeStep", FLOAT, N_total, 1, us,
+ */
/* UNIT_CONV_TIME); */
prepareArray(h_grp, fileName, xmfFile, "Acceleration", FLOAT, N_total, 3,
us, UNIT_CONV_ACCELERATION);
diff --git a/src/single_io.c b/src/single_io.c
index 7e92949e24f6618861050597e2e9335c26c8e1b8..10b84afeeeef77e72031536876ff9ebb4f37ca70 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -353,7 +353,8 @@ void writeArrayBackEnd(hid_t grp, char* fileName, FILE* xmfFile, char* name,
* @param type The #DATA_TYPE of the array.
* @param N The number of particles to write.
* @param dim The dimension of the data (1 for scalar, 3 for vector)
- * @param part A (char*) pointer on the first occurrence of the field of interest
+ * @param part A (char*) pointer on the first occurrence of the field of
+ *interest
*in the parts array
* @param field The name (code name) of the field to read from.
* @param us The UnitSystem currently in use
diff --git a/src/task.c b/src/task.c
index 642e6833c48236395f20c20b8e83f93d74757ea0..f736a1b30e14955fff070f586ebea0c54df403dd 100644
--- a/src/task.c
+++ b/src/task.c
@@ -43,9 +43,9 @@
/* Task type names. */
const char *taskID_names[task_type_count] = {
- "none", "sort", "self", "pair", "sub",
- "ghost", "kick1", "kick2", "send", "recv",
- "link", "grav_pp", "grav_mm", "grav_up", "grav_down"};
+ "none", "sort", "self", "pair", "sub", "init",
+ "ghost", "drift", "kick", "send", "recv", "link",
+ "grav_pp", "grav_mm", "grav_up", "grav_down"};
/**
* @brief Unlock the cell held by this task.
diff --git a/src/tools.c b/src/tools.c
index 5f8abcbeb6d339bedc3df71e5ccf651535b8bc0c..3973696b98954e17c8e6c4ee6535f1be8016ee58 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -19,7 +19,6 @@
*
******************************************************************************/
-
#include <math.h>
#include <stdlib.h>
#include <stddef.h>
@@ -31,7 +30,6 @@
#include "tools.h"
#include "swift.h"
-
/**
* Factorize a given integer, attempts to keep larger pair of factors.
*/
@@ -49,8 +47,6 @@ void factor(int value, int *f1, int *f2) {
}
}
-
-
/**
* @brief Compute the average number of pairs per particle using
* a brute-force O(N^2) computation.
@@ -234,7 +230,6 @@ void pairs_single_grav(double *dim, long long int pid,
pi.part->id, a[0], a[1], a[2], aabs[0], aabs[1], aabs[2]);
}
-
/**
* @brief Test the density function by dumping it for two random parts.
*
diff --git a/src/tools.h b/src/tools.h
index 7354a51030fef3a1fc8445afce770f56d5866fd6..ea7138672a3d2a037f1da98368c185bc4633fe08 100644
--- a/src/tools.h
+++ b/src/tools.h
@@ -24,9 +24,7 @@ void density_dump(int N);
void pairs_single_grav(double *dim, long long int pid,
struct gpart *__restrict__ parts, int N, int periodic);
void pairs_single_density(double *dim, long long int pid,
- struct part *__restrict__ parts, int N,
- int periodic);
+ struct part *__restrict__ parts, int N, int periodic);
void pairs_n2(double *dim, struct part *__restrict__ parts, int N,
int periodic);
-
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 9ab592e6ccde6194ffcc85c510554d195e7da9c4..2fc5e37578414d83c9858bd7da7c20fb98ff11cf 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -20,10 +20,10 @@ AM_CFLAGS = -I../src -DCPU_TPS=2.67e9 $(HDF5_CPPFLAGS)
AM_LDFLAGS = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
# List of programs and scripts to run in the test suite
-TESTS = testGreetings testReading.sh testSingle testTimeIntegration
+TESTS = testGreetings testReading.sh testSingle testTimeIntegration testSPHStep
# List of test programs to compile
-check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration
+check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration testSPHStep
# Sources for the individual programs
testGreetings_SOURCES = testGreetings.c
@@ -32,7 +32,8 @@ testReading_SOURCES = testReading.c
testTimeIntegration_SOURCES = testTimeIntegration.c
-# Sources for test_single
+testSPHStep_SOURCES = testSPHStep.c
+
testSingle_SOURCES = testSingle.c
testSingle_CFLAGS = $(MYFLAGS) $(AM_CFLAGS)
testSingle_LDADD = ../src/.libs/libswiftsim.a $(HDF5_LDFLAGS) $(HDF5_LIBS)
diff --git a/tests/testSPHStep.c b/tests/testSPHStep.c
new file mode 100644
index 0000000000000000000000000000000000000000..0f33c9829bd9763d84bb528bca76c792e2720c86
--- /dev/null
+++ b/tests/testSPHStep.c
@@ -0,0 +1,153 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (C) 2015 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/>.
+ *
+ ******************************************************************************/
+
+#include "swift.h"
+
+#include <stdlib.h>
+#include <string.h>
+
+/**
+ * @brief Constructs a cell with N SPH particles
+ */
+struct cell *make_cell(size_t N, float cellSize, int offset[3], int id_offset) {
+ size_t count = N*N*N;
+ struct cell *cell = malloc(sizeof( struct cell) );
+ struct part *part;
+ struct xpart *xpart;
+ float h;
+ size_t x, y, z, size;
+
+ size = count*sizeof(struct part);
+ if (posix_memalign((void**)&cell->parts, 32, size) != 0) {
+ error("couldn't allocate particles");
+ }
+
+ size = count*sizeof(struct xpart);
+ if (posix_memalign((void**)&cell->xparts, 32, size) != 0) {
+ error("couldn't allocate extended particles");
+ }
+
+
+ h = 1.127 * cellSize / N;
+
+ part = cell->parts;
+ xpart = cell->xparts;
+ memset(part, 0, count*sizeof(struct part));
+ memset(xpart, 0, count*sizeof(struct xpart));
+ for (x = 0; x < N; ++x) {
+ for (y = 0; y < N; ++y) {
+ for (z = 0; z < N; ++z) {
+ part->x[0] = offset[0] * cellSize + x * cellSize / N + cellSize / (2*N);
+ part->x[1] = offset[1] * cellSize + y * cellSize / N + cellSize / (2*N);
+ part->x[2] = offset[2] * cellSize + z * cellSize / N + cellSize / (2*N);
+ part->h = h;
+ part->id = x*N*N + y*N + z + id_offset;
+ ++part;
+ }
+ }
+ }
+
+ cell->split = 0;
+ cell->h_max = h;
+ cell->count = count;
+ cell->h[0] = cellSize;
+ cell->h[1] = cellSize;
+ cell->h[2] = cellSize;
+
+ return cell;
+}
+
+
+
+/* Run a full time step integration for one cell */
+int main() {
+
+ int i,j,k, offset[3];
+ struct part *p;
+
+ int N = 10;
+ float dim = 1.;
+ float rho = 2.;
+ float P = 1.;
+
+ /* Create cells */
+ struct cell *cells[27];
+ for (i=0; i<3; i++)
+ for (j=0; j<3; j++)
+ for (k=0; k<3; k++) {
+ offset[0] = i;
+ offset[1] = j;
+ offset[2] = k;
+ cells[i*9 + j*3 + k] = make_cell(N, dim, offset, (i*9+j*3+k) * N*N*N);
+ }
+
+ /* Set particle properties */
+ for (j=0 ; j < 27; ++j)
+ for (i=0 ; i< cells[j]->count; ++i) {
+ cells[j]->parts[i].mass = dim * dim * dim * rho / (N * N * N);
+ cells[j]->parts[i].u = P / ((const_hydro_gamma - 1.)*rho);
+ }
+
+ message("m=%f", dim * dim * dim * rho / (N * N * N));
+
+ /* Pick the central cell */
+ struct cell *ci = cells[13];
+
+ /* Create the infrastructure */
+ struct engine e;
+ struct runner r;
+ r.e = &e;
+
+ /* Simulation properties */
+ e.timeBegin = 0.;
+ e.timeEnd = 1.;
+ e.timeOld = 0.;
+ e.time = 0.;
+ e.dt_min = 0.;
+ e.dt_max = 1e10;
+
+ /* The tracked particle */
+ p = &(ci->parts[N*N*N / 2 + N*N / 2 + N / 2]);
+
+ message("Studying particle p->id=%lld", p->id);
+
+
+ /* Initialise the particles */
+ for (j=0 ; j<27; ++j)
+ {
+ runner_doinit(&r, cells[j]);
+ }
+
+ /* Compute density */
+ runner_doself1_density(&r, ci);
+ runner_doghost(&r, ci);
+
+ message("h=%f rho=%f N_ngb=%f", p->h, p->rho, p->density.wcount);
+ message("c=%f", p->force.c);
+
+ runner_doself2_force(&r, ci);
+ runner_dokick(&r, ci, 1);
+
+ message("t_end=%f", p->t_end);
+
+ free(ci->parts);
+ free(ci->xparts);
+
+ return 0;
+}