diff --git a/src/engine.c b/src/engine.c
index 7dd11d7bd20e96a18bd9d839056e1de2c435f586..e21e476cc30a3436bb7f26c932f87c5834de9e33 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -628,6 +628,7 @@ void engine_repartition(struct engine *e) {
void engine_addtasks_send(struct engine *e, struct cell *ci, struct cell *cj) {
+#ifdef WITH_MPI
int k;
struct link *l = NULL;
struct scheduler *s = &e->sched;
@@ -664,6 +665,11 @@ void engine_addtasks_send(struct engine *e, struct cell *ci, struct cell *cj) {
else if (ci->split)
for (k = 0; k < 8; k++)
if (ci->progeny[k] != NULL) engine_addtasks_send(e, ci->progeny[k], cj);
+
+#else
+ error("SWIFT was not compiled with MPI support.");
+#endif
+
}
/**
@@ -678,6 +684,7 @@ void engine_addtasks_send(struct engine *e, struct cell *ci, struct cell *cj) {
void engine_addtasks_recv(struct engine *e, struct cell *c, struct task *t_xv,
struct task *t_rho) {
+#ifdef WITH_MPI
int k;
struct scheduler *s = &e->sched;
@@ -707,6 +714,11 @@ void engine_addtasks_recv(struct engine *e, struct cell *c, struct task *t_xv,
for (k = 0; k < 8; k++)
if (c->progeny[k] != NULL)
engine_addtasks_recv(e, c->progeny[k], t_xv, t_rho);
+
+#else
+ error("SWIFT was not compiled with MPI support.");
+#endif
+
}
/**
@@ -1418,6 +1430,8 @@ void engine_print(struct engine *e) {
void engine_rebuild(struct engine *e) {
+ message("REBUILD !!!");
+
/* Clear the forcerebuild flag, whatever it was. */
e->forcerebuild = 0;
diff --git a/src/tools.c b/src/tools.c
index 3973696b98954e17c8e6c4ee6535f1be8016ee58..13a1df4187698337d7959918b12ac1407d924177 100644
--- a/src/tools.c
+++ b/src/tools.c
@@ -37,7 +37,7 @@ void factor(int value, int *f1, int *f2) {
int j;
int i;
- j = (int)sqrt(value);
+ j = (int) sqrt(value);
for (i = j; i > 0; i--) {
if ((value % i) == 0) {
*f1 = i;
@@ -110,7 +110,7 @@ void pairs_n2(double *dim, struct part *__restrict__ parts, int N,
/* Dump the result. */
printf("pairs_n2: avg. density per part is %.3f (nr. pairs %.3f).\n",
- rho / N + 32.0 / 3, ((double)count) / N);
+ rho / N + 32.0 / 3, ((double) count) / N);
printf("pairs_n2: densities are in [ %e , %e ].\n", rho_min / N + 32.0 / 3,
rho_max / N + 32.0 / 3);
/* printf( "pairs_n2: maximum ratio between parts %i [%e,%e,%e] and %i
@@ -183,7 +183,7 @@ void pairs_single_grav(double *dim, long long int pid,
// int mj, mk;
// double maxratio = 1.0;
double r2, dx[3];
- float fdx[3], a[3] = {0.0, 0.0, 0.0}, aabs[3] = {0.0, 0.0, 0.0};
+ float fdx[3], a[3] = { 0.0, 0.0, 0.0 }, aabs[3] = { 0.0, 0.0, 0.0 };
struct gpart pi, pj;
// double ih = 12.0/6.25;
@@ -239,7 +239,7 @@ void pairs_single_grav(double *dim, long long int pid,
void density_dump(int N) {
int k;
- float r2[4] = {0.0f, 0.0f, 0.0f, 0.0f}, hi[4], hj[4];
+ float r2[4] = { 0.0f, 0.0f, 0.0f, 0.0f }, hi[4], hj[4];
struct part *pi[4], *pj[4], Pi[4], Pj[4];
/* Init the interaction parameters. */
@@ -260,7 +260,7 @@ void density_dump(int N) {
r2[3] = r2[2];
r2[2] = r2[1];
r2[1] = r2[0];
- r2[0] = ((float)k) / N;
+ r2[0] = ((float) k) / N;
Pi[0].density.wcount = 0;
Pj[0].density.wcount = 0;
runner_iact_density(r2[0], NULL, hi[0], hj[0], &Pi[0], &Pj[0]);
@@ -278,3 +278,118 @@ void density_dump(int N) {
Pi[2].density.wcount, Pi[3].density.wcount);
}
}
+
+/**
+ * @brief Compute the force on a single particle brute-force.
+ */
+
+void engine_single_density(double *dim, long long int pid,
+ struct part *__restrict__ parts, int N,
+ int periodic) {
+
+ int i, k;
+ double r2, dx[3];
+ float fdx[3], ih;
+ struct part p;
+
+ /* Find "our" part. */
+ for (k = 0; k < N && parts[k].id != pid; k++)
+ ;
+ if (k == N) error("Part not found.");
+ p = parts[k];
+
+ /* Clear accumulators. */
+ ih = 1.0f / p.h;
+ p.rho = 0.0f;
+ p.rho_dh = 0.0f;
+ p.density.wcount = 0.0f;
+ p.density.wcount_dh = 0.0f;
+ p.density.div_v = 0.0;
+ for (k = 0; k < 3; k++)
+ p.density.curl_v[k] = 0.0;
+
+ /* Loop over all particle pairs (force). */
+ for (k = 0; k < N; k++) {
+ if (parts[k].id == p.id) continue;
+ for (i = 0; i < 3; i++) {
+ dx[i] = p.x[i] - parts[k].x[i];
+ if (periodic) {
+ if (dx[i] < -dim[i] / 2)
+ dx[i] += dim[i];
+ else if (dx[i] > dim[i] / 2)
+ dx[i] -= dim[i];
+ }
+ fdx[i] = dx[i];
+ }
+ r2 = fdx[0] * fdx[0] + fdx[1] * fdx[1] + fdx[2] * fdx[2];
+ if (r2 < p.h * p.h * kernel_gamma2) {
+ runner_iact_nonsym_density(r2, fdx, p.h, parts[k].h, &p, &parts[k]);
+ }
+ }
+
+ /* Dump the result. */
+ p.rho = ih * ih * ih * (p.rho + p.mass * kernel_root);
+ p.rho_dh = p.rho_dh * ih * ih * ih * ih;
+ p.density.wcount =
+ (p.density.wcount + kernel_root) * (4.0f / 3.0 * M_PI * kernel_gamma3);
+ message("part %lli (h=%e) has wcount=%e, rho=%e, rho_dh=%e.", p.id, p.h,
+ p.density.wcount, p.rho, p.rho_dh);
+ fflush(stdout);
+
+}
+
+void engine_single_force(double *dim, long long int pid,
+ struct part *__restrict__ parts, int N, int periodic) {
+
+ int i, k;
+ double r2, dx[3];
+ float fdx[3];
+ struct part p;
+
+ /* Find "our" part. */
+ for (k = 0; k < N && parts[k].id != pid; k++)
+ ;
+ if (k == N) error("Part not found.");
+ p = parts[k];
+
+ /* Clear accumulators. */
+ p.a[0] = 0.0f;
+ p.a[1] = 0.0f;
+ p.a[2] = 0.0f;
+ p.force.u_dt = 0.0f;
+ p.force.h_dt = 0.0f;
+ p.force.v_sig = 0.0f;
+
+ /* Loop over all particle pairs (force). */
+ for (k = 0; k < N; k++) {
+ // for ( k = N-1 ; k >= 0 ; k-- ) {
+ if (parts[k].id == p.id) continue;
+ for (i = 0; i < 3; i++) {
+ dx[i] = p.x[i] - parts[k].x[i];
+ if (periodic) {
+ if (dx[i] < -dim[i] / 2)
+ dx[i] += dim[i];
+ else if (dx[i] > dim[i] / 2)
+ dx[i] -= dim[i];
+ }
+ fdx[i] = dx[i];
+ }
+ r2 = fdx[0] * fdx[0] + fdx[1] * fdx[1] + fdx[2] * fdx[2];
+ if (r2 < p.h * p.h * kernel_gamma2 ||
+ r2 < parts[k].h * parts[k].h * kernel_gamma2) {
+ p.a[0] = 0.0f;
+ p.a[1] = 0.0f;
+ p.a[2] = 0.0f;
+ p.force.u_dt = 0.0f;
+ p.force.h_dt = 0.0f;
+ p.force.v_sig = 0.0f;
+ runner_iact_nonsym_force(r2, fdx, p.h, parts[k].h, &p, &parts[k]);
+ }
+ }
+
+ /* Dump the result. */
+ message( "part %lli (h=%e) has a=[%.3e,%.3e,%.3e], udt=%e." , p.id ,
+ p.h , p.a[0] , p.a[1] , p.a[2] , p.force.u_dt );
+ fflush(stdout);
+
+}