diff --git a/src/active.h b/src/active.h index edb276c5e9d7f420e8bcaa9cb9d15dae0ad4a880..19b0ed638df6bc76af7a548384760f1f24c9c76d 100644 --- a/src/active.h +++ b/src/active.h @@ -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 */ diff --git a/src/cell.c b/src/cell.c index f5ada2c6dfab4b8793df9f40b63deaf2b7e9009a..89ab7fb9add13db5c4bb613b456363809cc9add4 100644 --- a/src/cell.c +++ b/src/cell.c @@ -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 } diff --git a/src/drift.h b/src/drift.h index 5ad0a8ad26d431f03dc5332d79df2ff6c17dead8..687f8d8885a5fedca489f76d65ea8113101626c6 100644 --- a/src/drift.h +++ b/src/drift.h @@ -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 */ diff --git a/src/engine.c b/src/engine.c index 54e61a472250c99962d2028331ebe4010b301482..4092d95a0a06be582354a2384135f342b7101aed 100644 --- a/src/engine.c +++ b/src/engine.c @@ -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); } diff --git a/src/engine.h b/src/engine.h index 72d54f6d263eef9cc2f7d63a3b2aeddcf8b3c2e5..41d52cc6994f9e806878545a85dc33cbe3ef79bd 100644 --- a/src/engine.h +++ b/src/engine.h @@ -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; diff --git a/src/kick.h b/src/kick.h index b7dea7ffe9593d92a1eeef38b878c57328bad083..d6c85b5eab92a288f78f22fce2f03862bc34604f 100644 --- a/src/kick.h +++ b/src/kick.h @@ -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 */ diff --git a/src/runner.c b/src/runner.c index 5753973722f8ecb5c9299e5e2c9eb25cb8531926..d6f3d16276d3d8348c66794de890e1be86fdb575 100644 --- a/src/runner.c +++ b/src/runner.c @@ -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); } } } diff --git a/src/space.c b/src/space.c index d760222504b846c0e553e6433ce4b417f6dfa775..a1553bd18530fd1d43571fd7fc8c34e4d87762ec 100644 --- a/src/space.c +++ b/src/space.c @@ -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]); } diff --git a/src/space.h b/src/space.h index d4b2a753375d914f1bd08d884a8d41510f352792..53dc072f6b63c9c8658e61ac648d4eea0aaea848 100644 --- a/src/space.h +++ b/src/space.h @@ -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 */ diff --git a/src/stars/Default/star.h b/src/stars/Default/star.h index 5a33c45980e1bb99ca023f065c190b031cb2ab35..1f1b3940fe5da5468660262031001182ffc25cce 100644 --- a/src/stars/Default/star.h +++ b/src/stars/Default/star.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 */ diff --git a/src/timestep.h b/src/timestep.h index b3ef96493772e7e93853ad28a0536e4f8448dceb..432f0fd2c4eb713e11272546cfe84e8f6c342cbd 100644 --- a/src/timestep.h +++ b/src/timestep.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 */ + const integertime_t new_dti = make_integer_timestep( + new_dt, sp->time_bin, e->ti_current, e->timeBase_inv); + + return new_dti; +} + #endif /* SWIFT_TIMESTEP_H */