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

More time integration checks

parent 4b497bd7
......@@ -1289,14 +1289,6 @@ void cell_drift(struct cell *c, const struct engine *e) {
/* Drift... */
drift_gpart(gp, dt, timeBase, ti_old, ti_current);
#ifdef ICHECK
if (gp->id_or_neg_offset == ICHECK) {
message("--- ti_current=%lld time=%e dt=%e---", e->ti_current, e->time,
dt);
printgParticle_single(gp);
}
#endif
/* Compute (square of) motion since last cell construction */
const float dx2 = gp->x_diff[0] * gp->x_diff[0] +
gp->x_diff[1] * gp->x_diff[1] +
......
......@@ -81,7 +81,7 @@ __attribute__((always_inline)) INLINE static void drift_part(
#ifdef SWIFT_DEBUG_CHECKS
if (p->ti_drift != ti_old)
error(
"Particle has not been drifted to the current time p->ti_drift=%lld, "
"particle has not been drifted to the current time p->ti_drift=%lld, "
"c->ti_old=%lld, ti_current=%lld",
p->ti_drift, ti_old, ti_current);
......@@ -120,6 +120,17 @@ __attribute__((always_inline)) INLINE static void drift_spart(
struct spart *restrict sp, float dt, double timeBase, integertime_t ti_old,
integertime_t ti_current) {
#ifdef SWIFT_DEBUG_CHECKS
if (sp->ti_drift != ti_old)
error(
"s-particle has not been drifted to the current time "
"sp->ti_drift=%lld, "
"c->ti_old=%lld, ti_current=%lld",
sp->ti_drift, ti_old, ti_current);
sp->ti_drift = ti_current;
#endif
/* Drift... */
sp->x[0] += sp->v[0] * dt;
sp->x[1] += sp->v[1] * dt;
......
......@@ -85,6 +85,15 @@ __attribute__((always_inline)) INLINE static void gravity_end_force(
__attribute__((always_inline)) INLINE static void gravity_kick_extra(
struct gpart* gp, float dt) {}
/**
* @brief Sets the values to be predicted in the drifts to their values at a
* kick time
*
* @param gp The particle.
*/
__attribute__((always_inline)) INLINE static void gravity_reset_predicted_values(
struct gpart* gp) {}
/**
* @brief Initialises the g-particles for the first time
*
......
......@@ -46,7 +46,7 @@ __attribute__((always_inline)) INLINE static void kick_gpart(
#ifdef SWIFT_DEBUG_CHECKS
if (gp->ti_kick != ti_start)
error(
"Particle has not been kicked to the current time gp->ti_kick=%lld, "
"g-particle has not been kicked to the current time gp->ti_kick=%lld, "
"ti_start=%lld, ti_end=%lld",
gp->ti_kick, ti_start, ti_end);
......@@ -81,7 +81,7 @@ __attribute__((always_inline)) INLINE static void kick_part(
#ifdef SWIFT_DEBUG_CHECKS
if (p->ti_kick != ti_start)
error(
"Particle has not been kicked to the current time p->ti_kick=%lld, "
"particle has not been kicked to the current time p->ti_kick=%lld, "
"ti_start=%lld, ti_end=%lld",
p->ti_kick, ti_start, ti_end);
......@@ -126,6 +126,16 @@ __attribute__((always_inline)) INLINE static void kick_spart(
/* Time interval for this half-kick */
const float dt = (ti_end - ti_start) * timeBase;
#ifdef SWIFT_DEBUG_CHECKS
if (sp->ti_kick != ti_start)
error(
"s-particle has not been kicked to the current time sp->ti_kick=%lld, "
"ti_start=%lld, ti_end=%lld",
sp->ti_kick, ti_start, ti_end);
sp->ti_kick = ti_end;
#endif
/* Acceleration from gravity */
const float a[3] = {sp->gpart->a_grav[0], sp->gpart->a_grav[1],
sp->gpart->a_grav[2]};
......
......@@ -187,13 +187,6 @@ void runner_do_grav_external(struct runner *r, struct cell *c, int timer) {
/* Is this part within the time step? */
if (gpart_is_active(gp, e)) {
external_gravity_acceleration(time, potential, constants, gp);
#ifdef ICHECK
if (gp->id_or_neg_offset == ICHECK) {
message("--- ti_current=%lld time=%e ---", e->ti_current, e->time);
printgParticle_single(gp);
}
#endif
}
}
}
......@@ -531,13 +524,6 @@ void runner_do_init(struct runner *r, struct cell *c, int timer) {
/* Get ready for a density calculation */
gravity_init_gpart(gp);
#ifdef ICHECK
if (gp->id_or_neg_offset == ICHECK) {
message("--- ti_current=%lld time=%e ---", e->ti_current, e->time);
printgParticle_single(gp);
}
#endif
}
}
}
......@@ -935,13 +921,6 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
/* do the kick */
kick_gpart(gp, ti_begin, ti_begin + ti_step / 2, timeBase);
#ifdef ICHECK
if (gp->id_or_neg_offset == ICHECK) {
message("--- ti_current=%lld time=%e ---", e->ti_current, e->time);
printgParticle_single(gp);
}
#endif
}
}
......@@ -1036,6 +1015,12 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
/* Finish the time-step with a second half-kick */
kick_part(p, xp, ti_begin + ti_step / 2, ti_begin + ti_step, timeBase);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that kick and the drift are synchronized */
if(p->ti_drift != p->ti_kick)
error("Error integrating part in time.");
#endif
/* Prepare the values to be drifted */
hydro_reset_predicted_values(p, xp);
}
......@@ -1062,12 +1047,14 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
/* Finish the time-step with a second half-kick */
kick_gpart(gp, ti_begin + ti_step / 2, ti_begin + ti_step, timeBase);
#ifdef ICHECK
if (gp->id_or_neg_offset == ICHECK) {
message("--- ti_current=%lld time=%e ---", e->ti_current, e->time);
printgParticle_single(gp);
}
#ifdef SWIFT_DEBUG_CHECKS
/* Check that kick and the drift are synchronized */
if(gp->ti_drift != gp->ti_kick)
error("Error integrating g-part in time.");
#endif
/* Prepare the values to be drifted */
gravity_reset_predicted_values(gp);
}
}
......@@ -1092,6 +1079,12 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
/* Finish the time-step with a second half-kick */
kick_spart(sp, ti_begin + ti_step / 2, ti_begin + ti_step, timeBase);
#ifdef SWIFT_DEBUG_CHECKS
/* Check that kick and the drift are synchronized */
if(sp->ti_drift != sp->ti_kick)
error("Error integrating s-part in time.");
#endif
/* Prepare the values to be drifted */
star_reset_predicted_values(sp);
}
......@@ -1210,13 +1203,6 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
/* Update particle */
gp->time_bin = get_time_bin(ti_new_step);
#ifdef ICHECK
if (gp->id_or_neg_offset == ICHECK) {
message("--- ti_current=%lld time=%e ---", e->ti_current, e->time);
printgParticle_single(gp);
}
#endif
/* Number of updated g-particles */
g_updated++;
......@@ -1380,13 +1366,6 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
if (gp->type == swift_type_dark_matter) {
if (gpart_is_active(gp, e)) gravity_end_force(gp, const_G);
}
#ifdef ICHECK
if (gp->id_or_neg_offset == ICHECK) {
message("--- ti_current=%lld time=%e ---", e->ti_current, e->time);
printgParticle_single(gp);
}
#endif
}
/* Loop over the star particles in this cell. */
......
......@@ -2361,6 +2361,11 @@ void space_init_sparts(struct space *s) {
#endif
star_first_init_spart(&sp[i]);
#ifdef SWIFT_DEBUG_CHECKS
sp->ti_drift = 0;
sp->ti_kick = 0;
#endif
}
}
......
......@@ -47,6 +47,16 @@ struct spart {
/*! Particle time bin */
timebin_t time_bin;
#ifdef SWIFT_DEBUG_CHECKS
/* Time of the last drift */
integertime_t ti_drift;
/* Time of the last kick */
integertime_t ti_kick;
#endif
} SWIFT_STRUCT_ALIGN;
#endif /* SWIFT_DEFAULT_STAR_PART_H */
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