Commit 533e4a00 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Restored non-periodic BH-like gravity calculation.

parent fb757174
......@@ -167,6 +167,7 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
scheduler_addunlock(s, c->drift, c->init);
if (is_self_gravity) {
/* Gravity non-neighbouring pm calculations */
c->grav_long_range = scheduler_addtask(
s, task_type_grav_long_range, task_subtype_none, 0, 0, c, NULL, 0);
......@@ -175,8 +176,15 @@ void engine_make_hierarchical_tasks(struct engine *e, struct cell *c) {
c->grav_top_level = scheduler_addtask(
s, task_type_grav_top_level, task_subtype_none, 0, 0, c, NULL, 0);
/* Gravity recursive down-pass */
c->grav_down = scheduler_addtask(s, task_type_grav_down,
task_subtype_none, 0, 0, c, NULL, 0);
scheduler_addunlock(s, c->init, c->grav_long_range);
scheduler_addunlock(s, c->init, 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);
}
/* Generate the ghost task. */
......@@ -1799,8 +1807,7 @@ static inline void engine_make_self_gravity_dependencies(
/* init --> gravity --> grav_down --> kick */
scheduler_addunlock(sched, c->super->init, gravity);
// scheduler_addunlock(sched, gravity, c->super->grav_down);
scheduler_addunlock(sched, gravity, c->super->kick2);
scheduler_addunlock(sched, gravity, c->super->grav_down);
}
/**
......@@ -2169,27 +2176,28 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
*/
void engine_make_gravityrecursive_tasks(struct engine *e) {
struct space *s = e->s;
struct scheduler *sched = &e->sched;
const int nodeID = e->nodeID;
const int nr_cells = s->nr_cells;
struct cell *cells = s->cells_top;
/* struct space *s = e->s; */
/* struct scheduler *sched = &e->sched; */
/* const int nodeID = e->nodeID; */
/* const int nr_cells = s->nr_cells; */
/* struct cell *cells = s->cells_top; */
for (int k = 0; k < nr_cells; k++) {
/* for (int k = 0; k < nr_cells; k++) { */
/* Only do this for local cells containing gravity particles */
if (cells[k].nodeID == nodeID && cells[k].gcount > 0) {
/* /\* Only do this for local cells containing gravity particles *\/ */
/* if (cells[k].nodeID == nodeID && cells[k].gcount > 0) { */
/* Create tasks at top level. */
struct task *up = NULL;
struct task *down =
scheduler_addtask(sched, task_type_grav_down, task_subtype_none, 0, 0,
&cells[k], NULL, 0);
/* /\* Create tasks at top level. *\/ */
/* struct task *up = NULL; */
/* struct task *down = NULL; */
/* /\* scheduler_addtask(sched, task_type_grav_down,
* task_subtype_none, 0, 0, *\/ */
/* /\* &cells[k], NULL, 0); *\/ */
/* Push tasks down the cell hierarchy. */
engine_addtasks_grav(e, &cells[k], up, down);
}
}
/* /\* Push tasks down the cell hierarchy. *\/ */
/* engine_addtasks_grav(e, &cells[k], up, down); */
/* } */
/* } */
}
/**
......@@ -2842,6 +2850,8 @@ void engine_skip_force_and_kick(struct engine *e) {
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 ||
t->type == task_type_grav_long_range ||
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;
}
......
......@@ -113,7 +113,7 @@ const char runner_flip[27] = {1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
/* Import the gravity loop functions. */
#include "runner_doiact_fft.h"
//#include "runner_doiact_grav.h"
#include "runner_doiact_grav.h"
/**
* @brief Perform source terms
......@@ -1350,11 +1350,11 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
#ifdef SWIFT_DEBUG_CHECKS
if (e->policy & engine_policy_self_gravity) {
gp->mass_interacted += gp->mass;
if (gp->mass_interacted != e->s->total_mass)
if (fabs(gp->mass_interacted - e->s->total_mass) > gp->mass)
error(
"g-particle did not interact gravitationally with all other "
"particles gp->mass_interacted=%e, total_mass=%e",
gp->mass_interacted, e->s->total_mass);
"particles gp->mass_interacted=%e, total_mass=%e, gp->mass=%e",
gp->mass_interacted, e->s->total_mass, gp->mass);
}
#endif
}
......@@ -1686,8 +1686,7 @@ void *runner_main(void *data) {
else if (t->subtype == task_subtype_force)
runner_doself2_force(r, ci);
else if (t->subtype == task_subtype_grav)
// runner_doself_grav(r, ci, 1);
;
runner_doself_grav(r, ci, 1);
else if (t->subtype == task_subtype_external_grav)
runner_do_grav_external(r, ci, 1);
else
......@@ -1704,8 +1703,7 @@ void *runner_main(void *data) {
else if (t->subtype == task_subtype_force)
runner_dopair2_force(r, ci, cj);
else if (t->subtype == task_subtype_grav)
// runner_dopair_grav(r, ci, cj, 1);#
;
runner_dopair_grav(r, ci, cj, 1);
else
error("Unknown/invalid task subtype (%d).", t->subtype);
break;
......@@ -1720,8 +1718,7 @@ void *runner_main(void *data) {
else if (t->subtype == task_subtype_force)
runner_dosub_self2_force(r, ci, 1);
else if (t->subtype == task_subtype_grav)
// runner_dosub_grav(r, ci, cj, 1);#
;
runner_dosub_grav(r, ci, cj, 1);
else if (t->subtype == task_subtype_external_grav)
runner_do_grav_external(r, ci, 1);
else
......@@ -1738,8 +1735,7 @@ void *runner_main(void *data) {
else if (t->subtype == task_subtype_force)
runner_dosub_pair2_force(r, ci, cj, t->flags, 1);
else if (t->subtype == task_subtype_grav)
// runner_dosub_grav(r, ci, cj, 1);
;
runner_dosub_grav(r, ci, cj, 1);
else
error("Unknown/invalid task subtype (%d).", t->subtype);
break;
......@@ -1806,7 +1802,7 @@ void *runner_main(void *data) {
// runner_do_grav_top_level(r);
break;
case task_type_grav_long_range:
// runner_do_grav_fft(r);
runner_do_grav_long_range(r, t->ci, 1);
break;
case task_type_cooling:
if (e->policy & engine_policy_cooling) runner_do_end_force(r, ci, 1);
......
......@@ -34,23 +34,23 @@
*/
void runner_do_grav_up(struct runner *r, struct cell *c) {
if (c->split) { /* Regular node */
/* Recurse. */
for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) runner_do_grav_up(r, c->progeny[k]);
/* Collect the multipoles from the progeny. */
multipole_reset(&c->multipole);
for (int k = 0; k < 8; k++) {
if (c->progeny[k] != NULL)
multipole_add(&c->multipole, &c->progeny[k]->multipole);
}
} else { /* Leaf node. */
/* Just construct the multipole from the gparts. */
multipole_init(&c->multipole, c->gparts, c->gcount);
}
/* if (c->split) { /\* Regular node *\/ */
/* /\* Recurse. *\/ */
/* for (int k = 0; k < 8; k++) */
/* if (c->progeny[k] != NULL) runner_do_grav_up(r, c->progeny[k]); */
/* /\* Collect the multipoles from the progeny. *\/ */
/* multipole_reset(&c->multipole); */
/* for (int k = 0; k < 8; k++) { */
/* if (c->progeny[k] != NULL) */
/* multipole_add(&c->multipole, &c->progeny[k]->multipole); */
/* } */
/* } else { /\* Leaf node. *\/ */
/* /\* Just construct the multipole from the gparts. *\/ */
/* multipole_init(&c->multipole, c->gparts, c->gcount); */
/* } */
}
/**
......@@ -68,16 +68,15 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
const struct engine *e = r->e;
const int gcount = ci->gcount;
struct gpart *restrict gparts = ci->gparts;
const struct multipole multi = cj->multipole;
const struct multipole *multi = cj->multipole;
const float rlr_inv = 1. / (const_gravity_a_smooth * ci->super->width[0]);
TIMER_TIC;
#ifdef SWIFT_DEBUG_CHECKS
if (gcount == 0) error("Empty cell!"); // MATTHIEU sanity check
if (gcount == 0) error("Empty cell!");
if (multi.mass == 0.0) // MATTHIEU sanity check
error("Multipole does not seem to have been set.");
if (multi->mass == 0.0) error("Multipole does not seem to have been set.");
#endif
/* Anything to do here? */
......@@ -105,13 +104,17 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pm(
if (!gpart_is_active(gp, e)) continue;
/* Compute the pairwise distance. */
const float dx[3] = {multi.CoM[0] - gp->x[0], // x
multi.CoM[1] - gp->x[1], // y
multi.CoM[2] - gp->x[2]}; // z
const float dx[3] = {multi->CoM[0] - gp->x[0], // x
multi->CoM[1] - gp->x[1], // y
multi->CoM[2] - gp->x[2]}; // z
const float r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* Interact !*/
runner_iact_grav_pm(rlr_inv, r2, dx, gp, &multi);
runner_iact_grav_pm(rlr_inv, r2, dx, gp, multi);
#ifdef SWIFT_DEBUG_CHECKS
gp->mass_interacted += multi->mass;
#endif
}
TIMER_TOC(timer_dopair_grav_pm);
......@@ -195,18 +198,30 @@ __attribute__((always_inline)) INLINE static void runner_dopair_grav_pp(
runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
#ifdef SWIFT_DEBUG_CHECKS
gpi->mass_interacted += gpj->mass;
gpj->mass_interacted += gpi->mass;
#endif
} else {
if (gpart_is_active(gpi, e)) {
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
#ifdef SWIFT_DEBUG_CHECKS
gpi->mass_interacted += gpj->mass;
#endif
} else if (gpart_is_active(gpj, e)) {
dx[0] = -dx[0];
dx[1] = -dx[1];
dx[2] = -dx[2];
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
#ifdef SWIFT_DEBUG_CHECKS
gpj->mass_interacted += gpi->mass;
#endif
}
}
}
......@@ -277,18 +292,31 @@ __attribute__((always_inline)) INLINE static void runner_doself_grav_pp(
runner_iact_grav_pp(rlr_inv, r2, dx, gpi, gpj);
#ifdef SWIFT_DEBUG_CHECKS
gpi->mass_interacted += gpj->mass;
gpj->mass_interacted += gpi->mass;
#endif
} else {
if (gpart_is_active(gpi, e)) {
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpi, gpj);
#ifdef SWIFT_DEBUG_CHECKS
gpi->mass_interacted += gpj->mass;
#endif
} else if (gpart_is_active(gpj, e)) {
dx[0] = -dx[0];
dx[1] = -dx[1];
dx[2] = -dx[2];
runner_iact_grav_pp_nonsym(rlr_inv, r2, dx, gpj, gpi);
#ifdef SWIFT_DEBUG_CHECKS
gpj->mass_interacted += gpi->mass;
#endif
}
}
}
......@@ -466,7 +494,8 @@ static void runner_dosub_grav(struct runner *r, struct cell *ci,
}
}
static void runner_do_grav_mm(struct runner *r, struct cell *ci, int timer) {
static void runner_do_grav_long_range(struct runner *r, struct cell *ci,
int timer) {
#if ICHECK > 0
for (int pid = 0; pid < ci->gcount; pid++) {
......@@ -485,10 +514,10 @@ static void runner_do_grav_mm(struct runner *r, struct cell *ci, int timer) {
const struct engine *e = r->e;
struct cell *cells = e->s->cells_top;
const int nr_cells = e->s->nr_cells;
const double max_d =
const_gravity_a_smooth * const_gravity_r_cut * ci->width[0];
const double max_d2 = max_d * max_d;
const double pos_i[3] = {ci->loc[0], ci->loc[1], ci->loc[2]};
/* const double max_d = */
/* const_gravity_a_smooth * const_gravity_r_cut * ci->width[0]; */
/* const double max_d2 = max_d * max_d; */
// const double pos_i[3] = {ci->loc[0], ci->loc[1], ci->loc[2]};
/* Anything to do here? */
if (!cell_is_active(ci, e)) return;
......@@ -501,12 +530,12 @@ static void runner_do_grav_mm(struct runner *r, struct cell *ci, int timer) {
if (ci == cj) continue;
const double dx[3] = {cj->loc[0] - pos_i[0], // x
cj->loc[1] - pos_i[1], // y
cj->loc[2] - pos_i[2]}; // z
const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2];
/* const double dx[3] = {cj->loc[0] - pos_i[0], // x */
/* cj->loc[1] - pos_i[1], // y */
/* cj->loc[2] - pos_i[2]}; // z */
/* const double r2 = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; */
if (r2 > max_d2) continue;
// if (r2 > max_d2) continue;
if (!cell_are_neighbours(ci, cj)) runner_dopair_grav_pm(r, ci, cj);
}
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment