diff --git a/examples/main.c b/examples/main.c
index ba0fb491d71a8c2d925d8113e1c95262fb722074..8c75c376dd7c607915681f8eef2d7a7791cd3eb9 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -62,6 +62,11 @@
#define ENGINE_HYDRO 0
#endif
+#if defined(EXTERNAL_POTENTIAL)
+#define ENGINE_EXTERNAL_GRAVITY engine_policy_external_gravity
+#else
+#define ENGINE_EXTERNAL_GRAVITY 0
+#endif
/**
* @brief Main routine that loads a few particles and generates some output.
*
@@ -415,25 +420,12 @@ int main(int argc, char *argv[]) {
// message( "cutoffs in [ %g %g ]." , s.h_min , s.h_max ); fflush(stdout);
}
- /* Verify that each particle is in it's proper cell. */
- if (myrank == 0) {
- icount = 0;
- space_map_cells_pre(&s, 0, &map_cellcheck, &icount);
- message("map_cellcheck picked up %i parts.", icount);
- }
-
- if (myrank == 0) {
- data[0] = s.maxdepth;
- data[1] = 0;
- space_map_cells_pre(&s, 0, &map_maxdepth, data);
- message("nr of cells at depth %i is %i.", data[0], data[1]);
- }
-
+
/* Initialize the engine with this space. */
tic = getticks();
if (myrank == 0) message("nr_nodes is %i.", nr_nodes);
engine_init(&e, &s, dt_max, nr_threads, nr_queues, nr_nodes, myrank,
- ENGINE_POLICY | engine_policy_steal | ENGINE_HYDRO, 0,
+ ENGINE_POLICY | engine_policy_steal | ENGINE_HYDRO | ENGINE_EXTERNAL_GRAVITY, 0,
time_end, dt_min, dt_max);
if (myrank == 0)
message("engine_init took %.3f ms.",
@@ -482,6 +474,21 @@ int main(int argc, char *argv[]) {
/* Initialise the particles */
engine_init_particles(&e);
+ /* Verify that each particle is in it's proper cell. */
+ if (myrank == 0) {
+ icount = 0;
+ space_map_cells_pre(&s, 0, &map_cellcheck, &icount);
+ message("map_cellcheck picked up %i parts.", icount);
+
+ }
+
+ if (myrank == 0) {
+ data[0] = s.maxdepth;
+ data[1] = 0;
+ space_map_cells_pre(&s, 0, &map_maxdepth, data);
+ message("nr of cells at depth %i is %i.", data[0], data[1]);
+ }
+
//exit(-99);
/* Legend */
if (myrank == 0)
@@ -504,7 +511,7 @@ int main(int argc, char *argv[]) {
/* Take a step. */
engine_step(&e);
-
+
if (with_outputs && j % 100 == 0) {
#if defined(WITH_MPI)
diff --git a/src/cell.c b/src/cell.c
index 0186846ecdbeec1cb62253f9a47df0f25fd1dd8f..35338dfa1d037779a50916a0f3cf636a1af45b32 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -558,7 +558,11 @@ void cell_init_parts(struct cell *c, void *data) {
struct part *p = c->parts;
struct xpart *xp = c->xparts;
- const int count = c->count;
+ int count = c->gcount;
+ if(!count) count = c->count;
+
+
+
for (int i = 0; i < count; ++i) {
p[i].ti_begin = 0;
diff --git a/src/engine.c b/src/engine.c
index e82f340dd8f9c97028954a43aa22b3406c8b584f..8048b8e68e54567064ba944acda030a74aa4479f 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1166,10 +1166,11 @@ void engine_print(struct engine *e) {
void engine_rebuild(struct engine *e) {
/* Clear the forcerebuild flag, whatever it was. */
e->forcerebuild = 0;
-
+
/* Re-build the space. */
// tic = getticks();
space_rebuild(e->s, 0.0, e->nodeID == 0);
+
// message( "space_rebuild took %.3f ms." , (double)(getticks() -
// tic)/CPU_TPS*1000 );
@@ -1216,7 +1217,7 @@ void engine_prepare(struct engine *e) {
rebuild = (e->forcerebuild || engine_marktasks(e));
// message( "space_marktasks took %.3f ms." , (double)(getticks() -
// tic)/CPU_TPS*1000 );
-
+
/* Collect the values of rebuild from all nodes. */
#ifdef WITH_MPI
// tic = getticks();
@@ -1237,7 +1238,7 @@ void engine_prepare(struct engine *e) {
// message( "engine_rebuild took %.3f ms." , (double)(getticks() -
// tic)/CPU_TPS*1000 );
}
-
+
/* Re-rank the tasks every now and then. */
if (e->tasks_age % engine_tasksreweight == 1) {
// tic = getticks();
@@ -1404,14 +1405,9 @@ void engine_init_particles(struct engine *e) {
space_map_cells_pre(s, 1, cell_init_parts, NULL);
engine_prepare(e);
-
+
engine_marktasks(e);
- // printParticle(e->s->parts, 1000, e->s->nr_parts);
- // printParticle(e->s->parts, 515050, e->s->nr_parts);
-
- // message("\n0th DENSITY CALC\n");
-
/* Build the masks corresponding to the policy */
unsigned int mask = 0;
unsigned int submask = 0;
@@ -1440,7 +1436,7 @@ void engine_init_particles(struct engine *e) {
/* Add the tasks corresponding to self-gravity to the masks */
if ((e->policy & engine_policy_external_gravity) ==
engine_policy_external_gravity) {
-
+ printf("%s: JR: Excellent lets add the external gravity tasks here.....\n", __FUNCTION__);
/* Nothing here for now */
}
@@ -1979,6 +1975,9 @@ void engine_init(struct engine *e, struct space *s, float dt, int nr_threads,
s->nr_queues = nr_queues;
/* Create the sorting tasks. */
+ /* At the moment this is only implemented as a task for sorting parts,
+ * a gparts implementation is in the works (March 2016).
+ */
for (i = 0; i < e->nr_threads; i++)
scheduler_addtask(&e->sched, task_type_psort, task_subtype_none, i, 0, NULL,
NULL, 0);
diff --git a/src/hydro.h b/src/hydro.h
index 48606793cd97b96390cdc4ee3b37110a7346e5be..4851523f23acb04acb4c0260be77fb148353ecbd 100644
--- a/src/hydro.h
+++ b/src/hydro.h
@@ -37,7 +37,7 @@
#elif defined(NO_SPH)
#include "./hydro/Gadget2/hydro_iact.h"
#include "./hydro/Gadget2/hydro.h"
-#define SPH_IMPLEMENTATION "No SPH"
+#define SPH_IMPLEMENTATION "No SPH i.e. use Gadget2 framework for now..."
#else
#error "Invalid choice of SPH variant"
#endif
diff --git a/src/map.c b/src/map.c
index dcaa53465767d414cc54fe05940069ae5ff06d77..fb4f4aaea15881992788d29c5aeb2fff91b066ae 100644
--- a/src/map.c
+++ b/src/map.c
@@ -96,14 +96,13 @@ void map_cells_plot(struct cell *c, void *data) {
void map_cellcheck(struct cell *c, void *data) {
- int k, *count = (int *)data;
- struct part *p;
+ int *count = (int *)data;
__sync_fetch_and_add(count, c->count);
/* Loop over all parts and check if they are in the cell. */
- for (k = 0; k < c->count; k++) {
- p = &c->parts[k];
+ for (int k = 0; k < c->count; k++) {
+ struct part *p = &c->parts[k];
if (p->x[0] < c->loc[0] || p->x[1] < c->loc[1] || p->x[2] < c->loc[2] ||
p->x[0] > c->loc[0] + c->h[0] || p->x[1] > c->loc[1] + c->h[1] ||
p->x[2] > c->loc[2] + c->h[2]) {
@@ -115,6 +114,22 @@ void map_cellcheck(struct cell *c, void *data) {
error("particle out of bounds!");
}
}
+ __sync_fetch_and_add(count, c->gcount);
+
+ /* Loop over all gparts and check if they are in the cell. */
+ for (int k = 0; k < c->gcount; k++) {
+ struct gpart *p = &c->gparts[k];
+ if (p->x[0] < c->loc[0] || p->x[1] < c->loc[1] || p->x[2] < c->loc[2] ||
+ p->x[0] > c->loc[0] + c->h[0] || p->x[1] > c->loc[1] + c->h[1] ||
+ p->x[2] > c->loc[2] + c->h[2]) {
+ printf(
+ "map_cellcheck: g-particle at [ %.16e %.16e %.16e ] outside of cell "
+ "[ %.16e %.16e %.16e ] - [ %.16e %.16e %.16e ].\n",
+ p->x[0], p->x[1], p->x[2], c->loc[0], c->loc[1], c->loc[2],
+ c->loc[0] + c->h[0], c->loc[1] + c->h[1], c->loc[2] + c->h[2]);
+ error("g-particle out of bounds!");
+ }
+ }
}
/**
diff --git a/src/single_io.c b/src/single_io.c
index b346663c904ab59c4a738b58d1d4b3fbeaaffb53..cde6c3726371d23ea7a19d0e22e4a94d5476ea4a 100644
--- a/src/single_io.c
+++ b/src/single_io.c
@@ -407,8 +407,12 @@ void read_ic_single(char* fileName, double dim[3], struct part** parts, struct g
/* Read particle fields into the particle structure */
if(GAS == ptype)
hydro_read_particles(h_grp, *Ngas, *Ngas, 0, *parts);
- else if(DARKMATTER == ptype)
+ else if(DARKMATTER == ptype) {
darkmatter_read_particles(h_grp, *Ndm, *Ndm, 0, *gparts);
+ for(int i=0; i<*Ndm; i++) {
+ (*gparts)[i].id = -abs( (*gparts)[i].id );
+ }
+ }
else
error("Particle Type %d not yet supported. Aborting", ptype);
diff --git a/src/space.c b/src/space.c
index 05f28ffadc50390f488df3cc61241315932414cf..8b3f09ab7cbf89b0b129a8e34fe40c27a62a8f02 100644
--- a/src/space.c
+++ b/src/space.c
@@ -162,7 +162,6 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
int i, j, k, cdim[3], nr_parts = s->nr_parts;
struct cell *restrict c;
// ticks tic;
-
/* Run through the parts and get the current h_max. */
// tic = getticks();
if(nr_parts) {
@@ -178,7 +177,8 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
}
}
else {
- h_max = s->dim[0]/16.0;
+ /* It would be nice to replace this with something more physical or meaningful */
+ h_max = s->dim[0]/16.0;
s->h_max = h_max;
}
@@ -211,7 +211,6 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
if (cdim[0] < s->cdim[0] || cdim[1] < s->cdim[1] || cdim[2] < s->cdim[2])
error("Root-level change of cell size not allowed.");
#endif
-
/* Do we need to re-build the upper-level cells? */
// tic = getticks();
if (s->cells == NULL || cdim[0] < s->cdim[0] || cdim[1] < s->cdim[1] ||
@@ -234,7 +233,7 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
s->ih[k] = 1.0 / s->h[k];
}
dmin = fminf(s->h[0], fminf(s->h[1], s->h[2]));
-
+
/* Allocate the highest level of cells. */
s->tot_cells = s->nr_cells = cdim[0] * cdim[1] * cdim[2];
if (posix_memalign((void *)&s->cells, 64,
@@ -268,7 +267,7 @@ void space_regrid(struct space *s, double cell_max, int verbose) {
message("set cell dimensions to [ %i %i %i ].", cdim[0], cdim[1],
cdim[2]);
fflush(stdout);
-
+
} /* re-build upper-level cells? */
// message( "rebuilding upper-level cells took %.3f ms." , (double)(getticks()
// - tic) / CPU_TPS * 1000 );
@@ -324,7 +323,6 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
/* Re-grid if necessary, or just re-set the cell data. */
space_regrid(s, cell_max, verbose);
cells = s->cells;
-
/* Run through the SPH particles and get their cell index. */
// tic = getticks();
@@ -428,7 +426,6 @@ void space_rebuild(struct space *s, double cell_max, int verbose) {
cell_getid(cdim, gp->x[0] * ih[0], gp->x[1] * ih[1], gp->x[2] * ih[2]);
cells[ind[k]].gcount++;
}
-
// message( "getting particle indices took %.3f ms." , (double)(getticks() -
// tic) / CPU_TPS * 1000 );