Commit cdad747a authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Add the s-particles to the time integration tasks.

parent c6af1fa8
......@@ -139,4 +139,28 @@ __attribute__((always_inline)) INLINE static int gpart_is_active(
return (ti_end == ti_current);
}
/**
* @brief Is this s-particle active ?
*
* @param sp The #spart.
* @param e The #engine containing information about the current time.
* @return 1 if the #spart is active, 0 otherwise.
*/
__attribute__((always_inline)) INLINE static int spart_is_active(
const struct spart *sp, const struct engine *e) {
const integertime_t ti_current = e->ti_current;
const integertime_t ti_end = get_integer_time_end(ti_current, sp->time_bin);
#ifdef SWIFT_DEBUG_CHECKS
if (ti_end < ti_current)
error(
"s-particle in an impossible time-zone! gp->ti_end=%lld "
"e->ti_current=%lld",
ti_end, ti_current);
#endif
return (ti_end == ti_current);
}
#endif /* SWIFT_ACTIVE_H */
......@@ -1217,6 +1217,7 @@ void cell_drift(struct cell *c, const struct engine *e) {
struct part *const parts = c->parts;
struct xpart *const xparts = c->xparts;
struct gpart *const gparts = c->gparts;
struct spart *const sparts = c->sparts;
/* Drift from the last time the cell was drifted to the current time */
const double dt = (ti_current - ti_old) * timeBase;
......@@ -1256,7 +1257,7 @@ void cell_drift(struct cell *c, const struct engine *e) {
dx2_max = (dx2_max > dx2) ? dx2_max : dx2;
}
/* Loop over all the particles in the cell */
/* Loop over all the gas particles in the cell */
const size_t nr_parts = c->count;
for (size_t k = 0; k < nr_parts; k++) {
......@@ -1277,6 +1278,19 @@ void cell_drift(struct cell *c, const struct engine *e) {
h_max = (h_max > p->h) ? h_max : p->h;
}
/* Loop over all the star particles in the cell */
const size_t nr_sparts = c->scount;
for (size_t k = 0; k < nr_sparts; k++) {
/* Get a handle on the spart. */
struct spart *const sp = &sparts[k];
/* Drift... */
drift_spart(sp, dt, timeBase, ti_old, ti_current);
/* Note: no need to compute dx_max as all spart have a gpart */
}
/* Now, get the maximal particle motion from its square */
dx_max = sqrtf(dx2_max);
......@@ -1300,19 +1314,18 @@ void cell_drift(struct cell *c, const struct engine *e) {
void cell_check_timesteps(struct cell *c) {
#ifdef SWIFT_DEBUG_CHECKS
if(c->nodeID != engine_rank) return;
if (c->nodeID != engine_rank) return;
if(c->ti_end_min == 0) error("Cell without assigned time-step");
if(c->split) {
for(int k=0; k<8; ++k)
if(c->progeny[k] != NULL) cell_check_timesteps(c->progeny[k]);
} else {
if (c->ti_end_min == 0) error("Cell without assigned time-step");
for(int i=0; i<c->count; ++i)
if(c->parts[i].time_bin == 0)
error("Particle without assigned time-bin");
if (c->split) {
for (int k = 0; k < 8; ++k)
if (c->progeny[k] != NULL) cell_check_timesteps(c->progeny[k]);
} else {
for (int i = 0; i < c->count; ++i)
if (c->parts[i].time_bin == 0)
error("Particle without assigned time-bin");
}
#endif
}
......@@ -95,4 +95,23 @@ __attribute__((always_inline)) INLINE static void drift_part(
xp->x_diff[2] -= xp->v_full[2] * dt;
}
/**
* @brief Perform the 'drift' operation on a #spart
*
* @param sp The #spart to drift.
* @param dt The drift time-step
* @param timeBase The minimal allowed time-step size.
* @param ti_old Integer start of time-step
* @param ti_current Integer end of time-step
*/
__attribute__((always_inline)) INLINE static void drift_spart(
struct spart *restrict sp, float dt, double timeBase, integertime_t ti_old,
integertime_t ti_current) {
/* Drift... */
sp->x[0] += sp->v[0] * dt;
sp->x[1] += sp->v[1] * dt;
sp->x[2] += sp->v[2] * dt;
}
#endif /* SWIFT_DRIFT_H */
......@@ -2144,6 +2144,7 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
else if (t->type == task_type_timestep) {
t->ci->updated = 0;
t->ci->g_updated = 0;
t->ci->s_updated = 0;
if (t->ci->ti_end_min <= ti_end) scheduler_activate(s, t);
}
......@@ -2376,7 +2377,7 @@ void engine_collect_kick(struct cell *c) {
if (c->timestep != NULL) return;
/* Counters for the different quantities. */
int updated = 0, g_updated = 0;
int updated = 0, g_updated = 0, s_updated = 0;
integertime_t ti_end_min = max_nr_timesteps;
/* Only do something is the cell is non-empty */
......@@ -2397,10 +2398,12 @@ void engine_collect_kick(struct cell *c) {
ti_end_min = min(ti_end_min, cp->ti_end_min);
updated += cp->updated;
g_updated += cp->g_updated;
s_updated += cp->s_updated;
/* Collected, so clear for next time. */
cp->updated = 0;
cp->g_updated = 0;
cp->s_updated = 0;
}
}
}
......@@ -2409,6 +2412,7 @@ void engine_collect_kick(struct cell *c) {
c->ti_end_min = ti_end_min;
c->updated = updated;
c->g_updated = g_updated;
c->s_updated = s_updated;
}
/**
......@@ -2420,7 +2424,7 @@ void engine_collect_kick(struct cell *c) {
void engine_collect_timestep(struct engine *e) {
const ticks tic = getticks();
int updates = 0, g_updates = 0;
int updates = 0, g_updates = 0, s_updates = 0;
integertime_t ti_end_min = max_nr_timesteps;
const struct space *s = e->s;
......@@ -2436,10 +2440,12 @@ void engine_collect_timestep(struct engine *e) {
ti_end_min = min(ti_end_min, c->ti_end_min);
updates += c->updated;
g_updates += c->g_updated;
s_updates += c->s_updated;
/* Collected, so clear for next time. */
c->updated = 0;
c->g_updated = 0;
c->s_updated = 0;
}
/* Aggregate the data from the different nodes. */
......@@ -2454,20 +2460,23 @@ void engine_collect_timestep(struct engine *e) {
ti_end_min = in_i[0];
}
{
long long in_ll[2], out_ll[2];
long long in_ll[3], out_ll[3];
out_ll[0] = updates;
out_ll[1] = g_updates;
if (MPI_Allreduce(out_ll, in_ll, 2, MPI_LONG_LONG_INT, MPI_SUM,
out_ll[2] = s_updates;
if (MPI_Allreduce(out_ll, in_ll, 3, MPI_LONG_LONG_INT, MPI_SUM,
MPI_COMM_WORLD) != MPI_SUCCESS)
error("Failed to aggregate energies.");
updates = in_ll[0];
g_updates = in_ll[1];
s_updates = in_ll[2];
}
#endif
e->ti_end_min = ti_end_min;
e->updates = updates;
e->g_updates = g_updates;
e->s_updates = s_updates;
if (e->verbose)
message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
......@@ -2650,7 +2659,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs) {
#ifdef SWIFT_DEBUG_CHECKS
space_check_timesteps(e->s);
#endif
#endif
/* Ready to go */
e->step = -1;
......@@ -2710,12 +2719,14 @@ void engine_step(struct engine *e) {
if (e->nodeID == 0) {
/* Print some information to the screen */
printf(" %6d %14e %14e %10zu %10zu %21.3f\n", e->step, e->time,
e->timeStep, e->updates, e->g_updates, e->wallclock_time);
printf(" %6d %14e %14e %10zu %10zu %10zu %21.3f\n", e->step, e->time,
e->timeStep, e->updates, e->g_updates, e->s_updates,
e->wallclock_time);
fflush(stdout);
fprintf(e->file_timesteps, " %6d %14e %14e %10zu %10zu %21.3f\n", e->step,
e->time, e->timeStep, e->updates, e->g_updates, e->wallclock_time);
fprintf(e->file_timesteps, " %6d %14e %14e %10zu %10zu %10zu %21.3f\n",
e->step, e->time, e->timeStep, e->updates, e->g_updates,
e->s_updates, e->wallclock_time);
fflush(e->file_timesteps);
}
......
......@@ -131,7 +131,7 @@ struct engine {
integertime_t ti_end_min;
/* Number of particles updated */
size_t updates, g_updates;
size_t updates, g_updates, s_updates;
/* The internal system of units */
const struct UnitSystem *internalUnits;
......
......@@ -25,6 +25,7 @@
/* Local headers. */
#include "const.h"
#include "debug.h"
#include "stars.h"
#include "timeline.h"
/**
......@@ -100,4 +101,35 @@ __attribute__((always_inline)) INLINE static void kick_part(
if (p->gpart != NULL) gravity_kick_extra(p->gpart, dt);
}
/**
* @brief Perform the 'kick' operation on a #spart
*
* @param sp The #spart to kick.
* @param ti_start The starting (integer) time of the kick
* @param ti_end The ending (integer) time of the kick
* @param timeBase The minimal allowed time-step size.
*/
__attribute__((always_inline)) INLINE static void kick_spart(
struct spart *restrict sp, integertime_t ti_start, integertime_t ti_end,
double timeBase) {
/* Time interval for this half-kick */
const float dt = (ti_end - ti_start) * timeBase;
/* Acceleration from gravity */
const float a[3] = {sp->gpart->a_grav[0], sp->gpart->a_grav[1],
sp->gpart->a_grav[2]};
/* Kick particles in momentum space */
sp->v[0] += a[0] * dt;
sp->v[1] += a[1] * dt;
sp->v[2] += a[2] * dt;
sp->gpart->v_full[0] = sp->v[0];
sp->gpart->v_full[1] = sp->v[1];
sp->gpart->v_full[2] = sp->v[2];
/* Kick extra variables */
star_kick_extra(sp, dt);
}
#endif /* SWIFT_KICK_H */
......@@ -57,6 +57,7 @@
#include "scheduler.h"
#include "sourceterms.h"
#include "space.h"
#include "stars.h"
#include "task.h"
#include "timers.h"
#include "timestep.h"
......@@ -843,8 +844,10 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
struct part *restrict parts = c->parts;
struct xpart *restrict xparts = c->xparts;
struct gpart *restrict gparts = c->gparts;
struct spart *restrict sparts = c->sparts;
const int count = c->count;
const int gcount = c->gcount;
const int scount = c->scount;
const integertime_t ti_current = e->ti_current;
const double timeBase = e->timeBase;
......@@ -892,7 +895,7 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
struct gpart *restrict gp = &gparts[k];
/* If the g-particle has no counterpart and needs to be kicked */
if (gp->id_or_neg_offset > 0 && gpart_is_active(gp, e)) {
if (gp->type == swift_type_dark_matter && gpart_is_active(gp, e)) {
const integertime_t ti_step = get_integer_timestep(gp->time_bin);
const integertime_t ti_begin =
......@@ -909,7 +912,33 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
kick_gpart(gp, ti_begin, ti_begin + ti_step / 2, timeBase);
}
}
/* Loop over the star particles in this cell. */
for (int k = 0; k < scount; k++) {
/* Get a handle on the s-part. */
struct spart *restrict sp = &sparts[k];
/* If particle needs to be kicked */
if (spart_is_active(sp, e)) {
const integertime_t ti_step = get_integer_timestep(sp->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(ti_current, sp->time_bin);
#ifdef SWIFT_DEBUG_CHECKS
const integertime_t ti_end =
get_integer_time_end(ti_current, sp->time_bin);
if (ti_end - ti_begin != ti_step) error("Particle in wrong time-bin");
#endif
/* do the kick */
kick_spart(sp, ti_begin, ti_begin + ti_step / 2, timeBase);
}
}
}
if (timer) TIMER_TOC(timer_kick1);
}
......@@ -929,9 +958,11 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
const double timeBase = e->timeBase;
const int count = c->count;
const int gcount = c->gcount;
const int scount = c->scount;
struct part *restrict parts = c->parts;
struct xpart *restrict xparts = c->xparts;
struct gpart *restrict gparts = c->gparts;
struct spart *restrict sparts = c->sparts;
TIMER_TIC;
......@@ -978,7 +1009,7 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
struct gpart *restrict gp = &gparts[k];
/* If the g-particle has no counterpart and needs to be kicked */
if (gp->id_or_neg_offset > 0 && gpart_is_active(gp, e)) {
if (gp->type == swift_type_dark_matter && gpart_is_active(gp, e)) {
const integertime_t ti_step = get_integer_timestep(gp->time_bin);
const integertime_t ti_begin =
......@@ -993,6 +1024,32 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
kick_gpart(gp, ti_begin + ti_step / 2, ti_begin + ti_step, timeBase);
}
}
/* Loop over the particles in this cell. */
for (int k = 0; k < scount; k++) {
/* Get a handle on the part. */
struct spart *restrict sp = &sparts[k];
/* If particle needs to be kicked */
if (spart_is_active(sp, e)) {
const integertime_t ti_step = get_integer_timestep(sp->time_bin);
const integertime_t ti_begin =
get_integer_time_begin(ti_current, sp->time_bin);
#ifdef SWIFT_DEBUG_CHECKS
if (ti_begin + ti_step != ti_current)
error("Particle in wrong time-bin");
#endif
/* Finish the time-step with a second half-kick */
kick_spart(sp, ti_begin + ti_step / 2, ti_begin + ti_step, timeBase);
/* Prepare the values to be drifted */
star_reset_predicted_values(sp);
}
}
}
if (timer) TIMER_TOC(timer_kick2);
}
......@@ -1011,13 +1068,15 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
const integertime_t ti_current = e->ti_current;
const int count = c->count;
const int gcount = c->gcount;
const int scount = c->scount;
struct part *restrict parts = c->parts;
struct xpart *restrict xparts = c->xparts;
struct gpart *restrict gparts = c->gparts;
struct spart *restrict sparts = c->sparts;
TIMER_TIC;
int updated = 0, g_updated = 0;
int updated = 0, g_updated = 0, s_updated = 0;
integertime_t ti_end_min = max_nr_timesteps, ti_end_max = 0;
/* No children? */
......@@ -1076,7 +1135,7 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
struct gpart *restrict gp = &gparts[k];
/* If the g-particle has no counterpart */
if (gp->id_or_neg_offset > 0) {
if (gp->type == swift_type_dark_matter) {
/* need to be updated ? */
if (gpart_is_active(gp, e)) {
......@@ -1089,7 +1148,6 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
if (ti_end != ti_current)
error("Computing time-step of rogue particle.");
#endif
/* Get new time-step */
const integertime_t ti_new_step = get_gpart_timestep(gp, e);
......@@ -1102,6 +1160,7 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
/* What is the next sync-point ? */
ti_end_min = min(ti_current + ti_new_step, ti_end_min);
ti_end_max = max(ti_current + ti_new_step, ti_end_max);
} else { /* gpart is inactive */
const integertime_t ti_end =
......@@ -1113,6 +1172,48 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
}
}
}
/* Loop over the star particles in this cell. */
for (int k = 0; k < scount; k++) {
/* Get a handle on the part. */
struct spart *restrict sp = &sparts[k];
/* need to be updated ? */
if (spart_is_active(sp, e)) {
#ifdef SWIFT_DEBUG_CHECKS
/* Current end of time-step */
const integertime_t ti_end =
get_integer_time_end(ti_current, sp->time_bin);
if (ti_end != ti_current)
error("Computing time-step of rogue particle.");
#endif
/* Get new time-step */
const integertime_t ti_new_step = get_spart_timestep(sp, e);
/* Update particle */
sp->time_bin = get_time_bin(ti_new_step);
/* Number of updated s-particles */
s_updated++;
g_updated++;
/* What is the next sync-point ? */
ti_end_min = min(ti_current + ti_new_step, ti_end_min);
ti_end_max = max(ti_current + ti_new_step, ti_end_max);
} else { /* star particle is inactive */
const integertime_t ti_end =
get_integer_time_end(ti_current, sp->time_bin);
/* What is the next sync-point ? */
ti_end_min = min(ti_end, ti_end_min);
ti_end_max = max(ti_end, ti_end_max);
}
}
} else {
/* Loop over the progeny. */
......@@ -1126,6 +1227,7 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
/* And aggregate */
updated += cp->updated;
g_updated += cp->g_updated;
s_updated += cp->s_updated;
ti_end_min = min(cp->ti_end_min, ti_end_min);
ti_end_max = max(cp->ti_end_max, ti_end_max);
}
......@@ -1134,6 +1236,7 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
/* Store the values. */
c->updated = updated;
c->g_updated = g_updated;
c->s_updated = s_updated;
c->ti_end_min = ti_end_min;
c->ti_end_max = ti_end_max;
......@@ -1153,8 +1256,10 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
const struct engine *e = r->e;
const int count = c->count;
const int gcount = c->gcount;
const int scount = c->scount;
struct part *restrict parts = c->parts;
struct gpart *restrict gparts = c->gparts;
struct spart *restrict sparts = c->sparts;
const double const_G = e->physical_constants->const_newton_G;
TIMER_TIC;
......@@ -1168,7 +1273,7 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
if (c->progeny[k] != NULL) runner_do_end_force(r, c->progeny[k], 0);
} else {
/* Loop over the particles in this cell. */
/* Loop over the gas particles in this cell. */
for (int k = 0; k < count; k++) {
/* Get a handle on the part. */
......@@ -1188,8 +1293,22 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
/* Get a handle on the gpart. */
struct gpart *restrict gp = &gparts[k];
if (gp->id_or_neg_offset > 0 && gpart_is_active(gp, e)) {
gravity_end_force(gp, const_G);
if (gp->type == swift_type_dark_matter) {
if (gpart_is_active(gp, e)) gravity_end_force(gp, const_G);
}
}
/* Loop over the star particles in this cell. */
for (int k = 0; k < scount; k++) {
/* Get a handle on the part. */
struct spart *restrict sp = &sparts[k];
if (spart_is_active(sp, e)) {
/* First, finish the force loop */
star_end_force(sp);
gravity_end_force(sp->gpart, const_G);
}
}
}
......
......@@ -2394,11 +2394,11 @@ void space_check_drift_point(struct space *s, integertime_t ti_current) {
/**
* @brief Checks that all particles and local cells have a non-zero time-step.
*/
*/
void space_check_timesteps(struct space *s) {
#ifdef SWIFT_DEBUG_CHECKS
for (int i=0; i<s->nr_cells; ++i) {
for (int i = 0; i < s->nr_cells; ++i) {
cell_check_timesteps(&s->cells_top[i]);
}
......
......@@ -204,5 +204,4 @@ void space_check_drift_point(struct space *s, integertime_t ti_current);
void space_check_timesteps(struct space *s);
void space_clean(struct space *s);
#endif /* SWIFT_SPACE_H */
......@@ -27,7 +27,7 @@
*
* @param sp Pointer to the s-particle data.
*/
__attribute__((always_inline)) INLINE static float star_compute_timestep_self(
__attribute__((always_inline)) INLINE static float star_compute_timestep(
const struct spart* const sp) {
return FLT_MAX;
......@@ -55,6 +55,15 @@ __attribute__((always_inline)) INLINE static void star_first_init_spart(
__attribute__((always_inline)) INLINE static void star_init_spart(
struct spart* sp) {}
/**
* @brief Sets the values to be predicted in the drifts to their values at a
* kick time
*
* @param s The particle.
*/
__attribute__((always_inline)) INLINE static void star_reset_predicted_values(
struct spart* restrict sp) {}
/**
* @brief Finishes the calculation of (non-gravity) forces acting on stars
*
......@@ -73,6 +82,6 @@ __attribute__((always_inline)) INLINE static void star_end_force(
* @param half_dt The half time-step for this kick
*/
__attribute__((always_inline)) INLINE static void star_kick_extra(
struct spart* sp, float dt, float half_dt) {}
struct spart* sp, float dt) {}
#endif /* SWIFT_DEFAULT_STAR_H */
......@@ -146,4 +146,34 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep(
return new_dti;
}
/**
* @brief Compute the new (integer) time-step of a given #spart
*
* @param sp The #spart.
* @param e The #engine (used to get some constants).
*/
__attribute__((always_inline)) INLINE static integertime_t get_spart_timestep(
const struct spart *restrict sp, const struct engine *restrict e) {
float new_dt = star_compute_timestep(sp);
if (e->policy & engine_policy_external_gravity)
new_dt = min(new_dt,
external_gravity_timestep(e->time, e->external_potential,
e->physical_constants, sp->gpart));
if (e->policy & engine_policy_self_gravity)
new_dt = min(new_dt, gravity_compute_timestep_self(sp->gpart));
/* Limit timestep within the allowed range */
new_dt = min(new_dt, e->dt_max);
new_dt = max(new_dt, e->dt_min);
/* Convert to integer time */