Commit f01a4e0c authored by James Willis's avatar James Willis
Browse files

Merge branch 'master' into dopair-vectorisation-merge

parents ed2260b9 c2c61436
...@@ -91,6 +91,7 @@ theory/SPH/*.pdf ...@@ -91,6 +91,7 @@ theory/SPH/*.pdf
theory/paper_pasc/pasc_paper.pdf theory/paper_pasc/pasc_paper.pdf
theory/Multipoles/fmm.pdf theory/Multipoles/fmm.pdf
theory/Multipoles/fmm_standalone.pdf theory/Multipoles/fmm_standalone.pdf
theory/Multipoles/potential.pdf
m4/libtool.m4 m4/libtool.m4
m4/ltoptions.m4 m4/ltoptions.m4
......
...@@ -27,6 +27,7 @@ Valid options are: ...@@ -27,6 +27,7 @@ Valid options are:
-f {int} Overwrite the CPU frequency (Hz) to be used for time measurements. -f {int} Overwrite the CPU frequency (Hz) to be used for time measurements.
-g Run with an external gravitational potential. -g Run with an external gravitational potential.
-G Run with self-gravity. -G Run with self-gravity.
-M Reconstruct the multipoles every time-step.
-n {int} Execute a fixed number of time steps. When unset use the time_end parameter to stop. -n {int} Execute a fixed number of time steps. When unset use the time_end parameter to stop.
-s Run with hydrodynamics. -s Run with hydrodynamics.
-S Run with stars. -S Run with stars.
......
...@@ -215,6 +215,21 @@ elif test "$gravity_force_checks" != "no"; then ...@@ -215,6 +215,21 @@ elif test "$gravity_force_checks" != "no"; then
AC_DEFINE_UNQUOTED([SWIFT_GRAVITY_FORCE_CHECKS], [$enableval] ,[Enable gravity brute-force checks]) AC_DEFINE_UNQUOTED([SWIFT_GRAVITY_FORCE_CHECKS], [$enableval] ,[Enable gravity brute-force checks])
fi fi
# Check if we want to zero the gravity forces for all particles below some ID.
AC_ARG_ENABLE([no-gravity-below-id],
[AS_HELP_STRING([--enable-no-gravity-below-id],
[Zeros the gravitational acceleration of all particles with an ID smaller than @<:@N@:>@]
)],
[no_gravity_below_id="$enableval"],
[no_gravity_below_id="no"]
)
if test "$no_gravity_below_id" == "yes"; then
AC_MSG_ERROR(Need to specify the ID below which particles get zero forces when using --enable-no-gravity-below-id!)
elif test "$no_gravity_below_id" != "no"; then
AC_DEFINE_UNQUOTED([SWIFT_NO_GRAVITY_BELOW_ID], [$enableval] ,[Particles with smaller ID than this will have zero gravity forces])
fi
# Define HAVE_POSIX_MEMALIGN if it works. # Define HAVE_POSIX_MEMALIGN if it works.
AX_FUNC_POSIX_MEMALIGN AX_FUNC_POSIX_MEMALIGN
...@@ -854,12 +869,14 @@ AC_MSG_RESULT([ ...@@ -854,12 +869,14 @@ AC_MSG_RESULT([
Adiabatic index : $with_gamma Adiabatic index : $with_gamma
Riemann solver : $with_riemann Riemann solver : $with_riemann
Cooling function : $with_cooling Cooling function : $with_cooling
External potential : $with_potential
Multipole order : $with_multipole_order External potential : $with_potential
Multipole order : $with_multipole_order
Task debugging : $enable_task_debugging No gravity below ID : $no_gravity_below_id
Debugging checks : $enable_debugging_checks
Gravity checks : $gravity_force_checks Task debugging : $enable_task_debugging
Debugging checks : $enable_debugging_checks
Gravity checks : $gravity_force_checks
]) ])
# Make sure the latest git revision string gets included # Make sure the latest git revision string gets included
......
...@@ -82,6 +82,8 @@ void print_help_message() { ...@@ -82,6 +82,8 @@ void print_help_message() {
"Run with an external gravitational potential."); "Run with an external gravitational potential.");
printf(" %2s %8s %s\n", "-F", "", "Run with feedback."); printf(" %2s %8s %s\n", "-F", "", "Run with feedback.");
printf(" %2s %8s %s\n", "-G", "", "Run with self-gravity."); printf(" %2s %8s %s\n", "-G", "", "Run with self-gravity.");
printf(" %2s %8s %s\n", "-M", "",
"Reconstruct the multipoles every time-step.");
printf(" %2s %8s %s\n", "-n", "{int}", printf(" %2s %8s %s\n", "-n", "{int}",
"Execute a fixed number of time steps. When unset use the time_end " "Execute a fixed number of time steps. When unset use the time_end "
"parameter to stop."); "parameter to stop.");
...@@ -164,6 +166,7 @@ int main(int argc, char *argv[]) { ...@@ -164,6 +166,7 @@ int main(int argc, char *argv[]) {
int with_stars = 0; int with_stars = 0;
int with_fp_exceptions = 0; int with_fp_exceptions = 0;
int with_drift_all = 0; int with_drift_all = 0;
int with_mpole_reconstruction = 0;
int verbose = 0; int verbose = 0;
int nr_threads = 1; int nr_threads = 1;
int with_verbose_timers = 0; int with_verbose_timers = 0;
...@@ -172,7 +175,8 @@ int main(int argc, char *argv[]) { ...@@ -172,7 +175,8 @@ int main(int argc, char *argv[]) {
/* Parse the parameters */ /* Parse the parameters */
int c; int c;
while ((c = getopt(argc, argv, "acCdDef:FgGhn:sSt:Tv:y:")) != -1) switch (c) { while ((c = getopt(argc, argv, "acCdDef:FgGhMn:sSt:Tv:y:")) != -1)
switch (c) {
case 'a': case 'a':
with_aff = 1; with_aff = 1;
break; break;
...@@ -210,6 +214,9 @@ int main(int argc, char *argv[]) { ...@@ -210,6 +214,9 @@ int main(int argc, char *argv[]) {
case 'h': case 'h':
if (myrank == 0) print_help_message(); if (myrank == 0) print_help_message();
return 0; return 0;
case 'M':
with_mpole_reconstruction = 1;
break;
case 'n': case 'n':
if (sscanf(optarg, "%d", &nsteps) != 1) { if (sscanf(optarg, "%d", &nsteps) != 1) {
if (myrank == 0) printf("Error parsing fixed number of steps.\n"); if (myrank == 0) printf("Error parsing fixed number of steps.\n");
...@@ -521,6 +528,8 @@ int main(int argc, char *argv[]) { ...@@ -521,6 +528,8 @@ int main(int argc, char *argv[]) {
/* Construct the engine policy */ /* Construct the engine policy */
int engine_policies = ENGINE_POLICY | engine_policy_steal; int engine_policies = ENGINE_POLICY | engine_policy_steal;
if (with_drift_all) engine_policies |= engine_policy_drift_all; if (with_drift_all) engine_policies |= engine_policy_drift_all;
if (with_mpole_reconstruction)
engine_policies |= engine_policy_reconstruct_mpoles;
if (with_hydro) engine_policies |= engine_policy_hydro; if (with_hydro) engine_policies |= engine_policy_hydro;
if (with_self_gravity) engine_policies |= engine_policy_self_gravity; if (with_self_gravity) engine_policies |= engine_policy_self_gravity;
if (with_external_gravity) engine_policies |= engine_policy_external_gravity; if (with_external_gravity) engine_policies |= engine_policy_external_gravity;
...@@ -628,7 +637,7 @@ int main(int argc, char *argv[]) { ...@@ -628,7 +637,7 @@ int main(int argc, char *argv[]) {
for (int k = 0; k < timer_count; k++) for (int k = 0; k < timer_count; k++)
printf("%.3f\t", clocks_from_ticks(timers[k])); printf("%.3f\t", clocks_from_ticks(timers[k]));
printf("\n"); printf("\n");
timers_reset(0xFFFFFFFFllu); timers_reset(timers_mask_all);
} }
#ifdef SWIFT_DEBUG_TASKS #ifdef SWIFT_DEBUG_TASKS
......
...@@ -46,7 +46,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \ ...@@ -46,7 +46,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.h \ hydro_properties.h riemann.h threadpool.h cooling.h cooling_struct.h sourceterms.h \
sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h profiler.h \ sourceterms_struct.h statistics.h memswap.h cache.h runner_doiact_vec.h profiler.h \
dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \ dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \
vector_power.h collectgroup.h hydro_space.h sort_part.h gravity_softened_derivatives.h vector_power.h collectgroup.h hydro_space.h sort_part.h
# Common source files # Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
......
...@@ -129,6 +129,7 @@ int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) { ...@@ -129,6 +129,7 @@ int cell_unpack(struct pcell *pc, struct cell *c, struct space *s) {
temp->depth = c->depth + 1; temp->depth = c->depth + 1;
temp->split = 0; temp->split = 0;
temp->dx_max = 0.f; temp->dx_max = 0.f;
temp->dx_max_sort = 0.f;
temp->nodeID = c->nodeID; temp->nodeID = c->nodeID;
temp->parent = c; temp->parent = c;
c->progeny[k] = temp; c->progeny[k] = temp;
...@@ -1103,33 +1104,93 @@ void cell_reset_task_counters(struct cell *c) { ...@@ -1103,33 +1104,93 @@ void cell_reset_task_counters(struct cell *c) {
} }
/** /**
* @brief Checks whether the cells are direct neighbours ot not. Both cells have * @brief Recursively construct all the multipoles in a cell hierarchy.
* to be of the same size
* *
* @param ci First #cell. * @param c The #cell.
* @param cj Second #cell. * @param ti_current The current integer time.
*
* @todo Deal with periodicity.
*/ */
int cell_are_neighbours(const struct cell *restrict ci, void cell_make_multipoles(struct cell *c, integertime_t ti_current) {
const struct cell *restrict cj) {
#ifdef SWIFT_DEBUG_CHECKS /* Reset everything */
if (ci->width[0] != cj->width[0]) error("Cells of different size !"); gravity_reset(c->multipole);
#endif
/* Maximum allowed distance */ if (c->split) {
const double min_dist =
1.2 * ci->width[0]; /* 1.2 accounts for rounding errors */
/* (Manhattan) Distance between the cells */ /* Compute CoM of all progenies */
for (int k = 0; k < 3; k++) { double CoM[3] = {0., 0., 0.};
const double center_i = ci->loc[k]; double mass = 0.;
const double center_j = cj->loc[k];
if (fabs(center_i - center_j) > min_dist) return 0; for (int k = 0; k < 8; ++k) {
if (c->progeny[k] != NULL) {
const struct gravity_tensors *m = c->progeny[k]->multipole;
CoM[0] += m->CoM[0] * m->m_pole.M_000;
CoM[1] += m->CoM[1] * m->m_pole.M_000;
CoM[2] += m->CoM[2] * m->m_pole.M_000;
mass += m->m_pole.M_000;
}
}
c->multipole->CoM[0] = CoM[0] / mass;
c->multipole->CoM[1] = CoM[1] / mass;
c->multipole->CoM[2] = CoM[2] / mass;
/* Now shift progeny multipoles and add them up */
struct multipole temp;
double r_max = 0.;
for (int k = 0; k < 8; ++k) {
if (c->progeny[k] != NULL) {
const struct cell *cp = c->progeny[k];
const struct multipole *m = &cp->multipole->m_pole;
/* Contribution to multipole */
gravity_M2M(&temp, m, c->multipole->CoM, cp->multipole->CoM);
gravity_multipole_add(&c->multipole->m_pole, &temp);
/* Upper limit of max CoM<->gpart distance */
const double dx = c->multipole->CoM[0] - cp->multipole->CoM[0];
const double dy = c->multipole->CoM[1] - cp->multipole->CoM[1];
const double dz = c->multipole->CoM[2] - cp->multipole->CoM[2];
const double r2 = dx * dx + dy * dy + dz * dz;
r_max = max(r_max, cp->multipole->r_max + sqrt(r2));
}
}
/* Alternative upper limit of max CoM<->gpart distance */
const double dx = c->multipole->CoM[0] > c->loc[0] + c->width[0] / 2.
? c->multipole->CoM[0] - c->loc[0]
: c->loc[0] + c->width[0] - c->multipole->CoM[0];
const double dy = c->multipole->CoM[1] > c->loc[1] + c->width[1] / 2.
? c->multipole->CoM[1] - c->loc[1]
: c->loc[1] + c->width[1] - c->multipole->CoM[1];
const double dz = c->multipole->CoM[2] > c->loc[2] + c->width[2] / 2.
? c->multipole->CoM[2] - c->loc[2]
: c->loc[2] + c->width[2] - c->multipole->CoM[2];
/* Take minimum of both limits */
c->multipole->r_max = min(r_max, sqrt(dx * dx + dy * dy + dz * dz));
} else {
if (c->gcount > 0) {
gravity_P2M(c->multipole, c->gparts, c->gcount);
const double dx = c->multipole->CoM[0] > c->loc[0] + c->width[0] / 2.
? c->multipole->CoM[0] - c->loc[0]
: c->loc[0] + c->width[0] - c->multipole->CoM[0];
const double dy = c->multipole->CoM[1] > c->loc[1] + c->width[1] / 2.
? c->multipole->CoM[1] - c->loc[1]
: c->loc[1] + c->width[1] - c->multipole->CoM[1];
const double dz = c->multipole->CoM[2] > c->loc[2] + c->width[2] / 2.
? c->multipole->CoM[2] - c->loc[2]
: c->loc[2] + c->width[2] - c->multipole->CoM[2];
c->multipole->r_max = sqrt(dx * dx + dy * dy + dz * dz);
} else {
gravity_multipole_init(&c->multipole->m_pole);
c->multipole->CoM[0] = c->loc[0] + c->width[0] / 2.;
c->multipole->CoM[1] = c->loc[1] + c->width[1] / 2.;
c->multipole->CoM[2] = c->loc[2] + c->width[2] / 2.;
c->multipole->r_max = 0.;
}
} }
return 1; c->ti_old_multipole = ti_current;
} }
/** /**
...@@ -1145,6 +1206,8 @@ void cell_check_multipole(struct cell *c, void *data) { ...@@ -1145,6 +1206,8 @@ void cell_check_multipole(struct cell *c, void *data) {
struct gravity_tensors ma; struct gravity_tensors ma;
const double tolerance = 1e-3; /* Relative */ const double tolerance = 1e-3; /* Relative */
return;
/* First recurse */ /* First recurse */
if (c->split) if (c->split)
for (int k = 0; k < 8; k++) for (int k = 0; k < 8; k++)
...@@ -1244,28 +1307,45 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) { ...@@ -1244,28 +1307,45 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
/* Un-skip the density tasks involved with this cell. */ /* Un-skip the density tasks involved with this cell. */
for (struct link *l = c->density; l != NULL; l = l->next) { for (struct link *l = c->density; l != NULL; l = l->next) {
struct task *t = l->t; struct task *t = l->t;
const struct cell *ci = t->ci; struct cell *ci = t->ci;
const struct cell *cj = t->cj; struct cell *cj = t->cj;
scheduler_activate(s, t); scheduler_activate(s, t);
/* Set the correct sorting flags */ /* Set the correct sorting flags */
if (t->type == task_type_pair) { if (t->type == task_type_pair) {
if (ci->dx_max_sort > space_maxreldx * ci->dmin) {
for (struct cell *finger = ci; finger != NULL; finger = finger->parent)
finger->sorted = 0;
}
if (cj->dx_max_sort > space_maxreldx * cj->dmin) {
for (struct cell *finger = cj; finger != NULL; finger = finger->parent)
finger->sorted = 0;
}
if (!(ci->sorted & (1 << t->flags))) { if (!(ci->sorted & (1 << t->flags))) {
atomic_or(&ci->sorts->flags, (1 << t->flags)); #ifdef SWIFT_DEBUG_CHECKS
if (!(ci->sorts->flags & (1 << t->flags)))
error("bad flags in sort task.");
#endif
scheduler_activate(s, ci->sorts); scheduler_activate(s, ci->sorts);
if (ci->nodeID == engine_rank) scheduler_activate(s, ci->drift);
} }
if (!(cj->sorted & (1 << t->flags))) { if (!(cj->sorted & (1 << t->flags))) {
atomic_or(&cj->sorts->flags, (1 << t->flags)); #ifdef SWIFT_DEBUG_CHECKS
if (!(cj->sorts->flags & (1 << t->flags)))
error("bad flags in sort task.");
#endif
scheduler_activate(s, cj->sorts); scheduler_activate(s, cj->sorts);
if (cj->nodeID == engine_rank) scheduler_activate(s, cj->drift);
} }
} }
/* Check whether there was too much particle motion */ /* Only interested in pair interactions as of here. */
if (t->type == task_type_pair || t->type == task_type_sub_pair) { if (t->type == task_type_pair || t->type == task_type_sub_pair) {
/* Check whether there was too much particle motion, i.e. the
cell neighbour conditions were violated. */
if (t->tight && if (t->tight &&
(max(ci->h_max, cj->h_max) + ci->dx_max + cj->dx_max > cj->dmin || max(ci->h_max, cj->h_max) + ci->dx_max + cj->dx_max > cj->dmin)
ci->dx_max > space_maxreldx * ci->h_max ||
cj->dx_max > space_maxreldx * cj->h_max))
rebuild = 1; rebuild = 1;
#ifdef WITH_MPI #ifdef WITH_MPI
...@@ -1287,10 +1367,12 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) { ...@@ -1287,10 +1367,12 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
if (l == NULL) error("Missing link to send_xv task."); if (l == NULL) error("Missing link to send_xv task.");
scheduler_activate(s, l->t); scheduler_activate(s, l->t);
if (cj->super->drift) /* Drift both cells, the foreign one at the level which it is sent. */
scheduler_activate(s, cj->super->drift); if (l->t->ci->drift)
scheduler_activate(s, l->t->ci->drift);
else else
error("Drift task missing !"); error("Drift task missing !");
if (t->type == task_type_pair) scheduler_activate(s, cj->drift);
if (cell_is_active(cj, e)) { if (cell_is_active(cj, e)) {
for (l = cj->send_rho; l != NULL && l->t->cj->nodeID != ci->nodeID; for (l = cj->send_rho; l != NULL && l->t->cj->nodeID != ci->nodeID;
...@@ -1323,10 +1405,12 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) { ...@@ -1323,10 +1405,12 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
if (l == NULL) error("Missing link to send_xv task."); if (l == NULL) error("Missing link to send_xv task.");
scheduler_activate(s, l->t); scheduler_activate(s, l->t);
if (ci->super->drift) /* Drift both cells, the foreign one at the level which it is sent. */
scheduler_activate(s, ci->super->drift); if (l->t->ci->drift)
scheduler_activate(s, l->t->ci->drift);
else else
error("Drift task missing !"); error("Drift task missing !");
if (t->type == task_type_pair) scheduler_activate(s, ci->drift);
if (cell_is_active(ci, e)) { if (cell_is_active(ci, e)) {
for (l = ci->send_rho; l != NULL && l->t->cj->nodeID != cj->nodeID; for (l = ci->send_rho; l != NULL && l->t->cj->nodeID != cj->nodeID;
...@@ -1341,6 +1425,14 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) { ...@@ -1341,6 +1425,14 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
if (l == NULL) error("Missing link to send_ti task."); if (l == NULL) error("Missing link to send_ti task.");
scheduler_activate(s, l->t); scheduler_activate(s, l->t);
} }
} else if (t->type == task_type_pair) {
scheduler_activate(s, ci->drift);
scheduler_activate(s, cj->drift);
}
#else
if (t->type == task_type_pair) {
scheduler_activate(s, ci->drift);
scheduler_activate(s, cj->drift);
} }
#endif #endif
} }
...@@ -1355,7 +1447,6 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) { ...@@ -1355,7 +1447,6 @@ int cell_unskip_tasks(struct cell *c, struct scheduler *s) {
scheduler_activate(s, l->t); scheduler_activate(s, l->t);
if (c->extra_ghost != NULL) scheduler_activate(s, c->extra_ghost); if (c->extra_ghost != NULL) scheduler_activate(s, c->extra_ghost);
if (c->ghost != NULL) scheduler_activate(s, c->ghost); if (c->ghost != NULL) scheduler_activate(s, c->ghost);
if (c->init != NULL) scheduler_activate(s, c->init);
if (c->init_grav != NULL) scheduler_activate(s, c->init_grav); if (c->init_grav != NULL) scheduler_activate(s, c->init_grav);
if (c->drift != NULL) scheduler_activate(s, c->drift); if (c->drift != NULL) scheduler_activate(s, c->drift);
if (c->kick1 != NULL) scheduler_activate(s, c->kick1); if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
...@@ -1409,7 +1500,9 @@ void cell_drift_particles(struct cell *c, const struct engine *e) { ...@@ -1409,7 +1500,9 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
/* Drift from the last time the cell was drifted to the current time */ /* 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) * timeBase;
float dx_max = 0.f, dx2_max = 0.f, cell_h_max = 0.f; 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. */ /* 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) error("Attempt to drift to the past");
...@@ -1421,8 +1514,13 @@ void cell_drift_particles(struct cell *c, const struct engine *e) { ...@@ -1421,8 +1514,13 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
for (int k = 0; k < 8; k++) for (int k = 0; k < 8; k++)
if (c->progeny[k] != NULL) { if (c->progeny[k] != NULL) {
struct cell *cp = c->progeny[k]; struct cell *cp = c->progeny[k];
/* Collect */
cell_drift_particles(cp, e); cell_drift_particles(cp, e);
/* Update */
dx_max = max(dx_max, cp->dx_max); dx_max = max(dx_max, cp->dx_max);
dx_max_sort = max(dx_max_sort, cp->dx_max_sort);
cell_h_max = max(cell_h_max, cp->h_max); cell_h_max = max(cell_h_max, cp->h_max);
} }
...@@ -1443,6 +1541,11 @@ void cell_drift_particles(struct cell *c, const struct engine *e) { ...@@ -1443,6 +1541,11 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
gp->x_diff[1] * gp->x_diff[1] + gp->x_diff[1] * gp->x_diff[1] +
gp->x_diff[2] * gp->x_diff[2]; gp->x_diff[2] * gp->x_diff[2];
dx2_max = max(dx2_max, dx2); dx2_max = max(dx2_max, dx2);
/* Init gravity force fields. */
if (gpart_is_active(gp, e)) {
gravity_init_gpart(gp);
}
} }
/* Loop over all the gas particles in the cell */ /* Loop over all the gas particles in the cell */
...@@ -1464,9 +1567,18 @@ void cell_drift_particles(struct cell *c, const struct engine *e) { ...@@ -1464,9 +1567,18 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
xp->x_diff[1] * xp->x_diff[1] + xp->x_diff[1] * xp->x_diff[1] +
xp->x_diff[2] * xp->x_diff[2]; xp->x_diff[2] * xp->x_diff[2];
dx2_max = max(dx2_max, dx2); dx2_max = max(dx2_max, dx2);
const float dx2_sort = xp->x_diff_sort[0] * xp->x_diff_sort[0] +
xp->x_diff_sort[1] * xp->x_diff_sort[1] +
xp->x_diff_sort[2] * xp->x_diff_sort[2];
dx2_max_sort = max(dx2_max_sort, dx2_sort);
/* Maximal smoothing length */ /* Maximal smoothing length */
cell_h_max = max(cell_h_max, p->h); cell_h_max = max(cell_h_max, p->h);
/* Get ready for a density calculation */
if (part_is_active(p, e)) {
hydro_init_part(p, &e->s->hs);
}
} }
/* Loop over all the star particles in the cell */ /* Loop over all the star particles in the cell */
...@@ -1484,16 +1596,19 @@ void cell_drift_particles(struct cell *c, const struct engine *e) { ...@@ -1484,16 +1596,19 @@ void cell_drift_particles(struct cell *c, const struct engine *e) {
/* Now, get the maximal particle motion from its square */ /* Now, get the maximal particle motion from its square */
dx_max = sqrtf(dx2_max); dx_max = sqrtf(dx2_max);
dx_max_sort = sqrtf(dx2_max_sort);
} else { } else {
cell_h_max = c->h_max; cell_h_max = c->h_max;
dx_max = c->dx_max; dx_max = c->dx_max;
dx_max_sort = c->dx_max_sort;
} }
/* Store the values */ /* Store the values */
c->h_max = cell_h_max; c->h_max = cell_h_max;
c->dx_max = dx_max; c->dx_max = dx_max;
c->dx_max_sort = dx_max_sort;
/* Update the time of the last drift */ /* Update the time of the last drift */
c->ti_old = ti_current; c->ti_old = ti_current;
......
...@@ -148,9 +148,6 @@ struct cell { ...@@ -148,9 +148,6 @@ struct cell {
/*! Linked list of the tasks computing this cell's gravity forces. */ /*! Linked list of the tasks computing this cell's gravity forces. */
struct link *grav; struct link *grav;
/*! The particle initialistation task */
struct task *init;
/*! The multipole initialistation task */ /*! The multipole initialistation task */
struct task *init_grav; struct task *init_grav;
...@@ -239,6 +236,9 @@ struct cell { ...@@ -239,6 +236,9 @@ struct cell {
/*! Last (integer) time the cell's particle was drifted forward in time. */ /*! Last (integer) time the cell's particle was drifted forward in time. */
integertime_t ti_old; integertime_t ti_old;
/*! Last (integer) time the cell's sort arrays were updated. */
integertime_t ti_sort;
/*! Last (integer) time the cell's multipole was drifted forward in time. */ /*! Last (integer) time the cell's multipole was drifted forward in time. */
integertime_t ti_old_multipole; integertime_t ti_old_multipole;
...@@ -248,6 +248,9 @@ struct cell { ...@@ -248,6 +248,9 @@ struct cell {
/*! Maximum particle movement in this cell since last construction. */ /*! Maximum particle movement in this cell since last construction. */
float dx_max; float dx_max;