Commit 1d1e86c2 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Misplaced if-statement leading to 0-sized time-step for gparts with counterparts

parent f3c535b4
...@@ -929,75 +929,78 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) { ...@@ -929,75 +929,78 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
struct gpart *const gp = &gparts[k]; struct gpart *const gp = &gparts[k];
/* If the g-particle has no counterpart and needs to be kicked */ /* If the g-particle has no counterpart and needs to be kicked */
if (gp->id < 0 && (is_fixdt || gp->ti_end <= ti_current)) { if (gp->id < 0) {
/* First, finish the force calculation */ if (is_fixdt || gp->ti_end <= ti_current) {
gravity_end_force(gp);
/* First, finish the force calculation */
/* Now we are ready to compute the next time-step size */ gravity_end_force(gp);
int new_dti;
/* Now we are ready to compute the next time-step size */
if (is_fixdt) { int new_dti;
/* Now we have a time step, proceed with the kick */ if (is_fixdt) {
new_dti = global_dt_max * timeBase_inv;
/* Now we have a time step, proceed with the kick */
} else { new_dti = global_dt_max * timeBase_inv;
/* Compute the next timestep (gravity condition) */ } else {
const float new_dt_external =
/* Compute the next timestep (gravity condition) */
const float new_dt_external =
gravity_compute_timestep_external(potential, constants, gp); gravity_compute_timestep_external(potential, constants, gp);
const float new_dt_self = const float new_dt_self =
gravity_compute_timestep_self(constants, gp); gravity_compute_timestep_self(constants, gp);
float new_dt = fminf(new_dt_external, new_dt_self); float new_dt = fminf(new_dt_external, new_dt_self);
/* Limit timestep within the allowed range */ /* Limit timestep within the allowed range */
new_dt = fminf(new_dt, global_dt_max); new_dt = fminf(new_dt, global_dt_max);
new_dt = fmaxf(new_dt, global_dt_min); new_dt = fmaxf(new_dt, global_dt_min);
/* Convert to integer time */ /* Convert to integer time */
new_dti = new_dt * timeBase_inv; new_dti = new_dt * timeBase_inv;
/* Recover the current timestep */ /* Recover the current timestep */
const int current_dti = gp->ti_end - gp->ti_begin; const int current_dti = gp->ti_end - gp->ti_begin;
/* Limit timestep increase */ /* Limit timestep increase */
if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti); if (current_dti > 0) new_dti = min(new_dti, 2 * current_dti);
/* Put this timestep on the time line */ /* Put this timestep on the time line */
int dti_timeline = max_nr_timesteps; int dti_timeline = max_nr_timesteps;
while (new_dti < dti_timeline) dti_timeline /= 2; while (new_dti < dti_timeline) dti_timeline /= 2;
/* Now we have a time step, proceed with the kick */ /* Now we have a time step, proceed with the kick */
new_dti = dti_timeline; new_dti = dti_timeline;
} }
/* Compute the time step for this kick */ /* Compute the time step for this kick */
const int ti_start = (gp->ti_begin + gp->ti_end) / 2; const int ti_start = (gp->ti_begin + gp->ti_end) / 2;
const int ti_end = gp->ti_end + new_dti / 2; const int ti_end = gp->ti_end + new_dti / 2;
const double dt = (ti_end - ti_start) * timeBase; const double dt = (ti_end - ti_start) * timeBase;
const double half_dt = (ti_end - gp->ti_end) * timeBase; const double half_dt = (ti_end - gp->ti_end) * timeBase;
/* Move particle forward in time */ /* Move particle forward in time */
gp->ti_begin = gp->ti_end; gp->ti_begin = gp->ti_end;
gp->ti_end = gp->ti_begin + new_dti; gp->ti_end = gp->ti_begin + new_dti;
/* Kick particles in momentum space */ /* Kick particles in momentum space */
gp->v_full[0] += gp->a_grav[0] * dt; gp->v_full[0] += gp->a_grav[0] * dt;
gp->v_full[1] += gp->a_grav[1] * dt; gp->v_full[1] += gp->a_grav[1] * dt;
gp->v_full[2] += gp->a_grav[2] * dt; gp->v_full[2] += gp->a_grav[2] * dt;
/* Extra kick work */ /* Extra kick work */
gravity_kick_extra(gp, dt, half_dt); gravity_kick_extra(gp, dt, half_dt);
/* Number of updated g-particles */ /* Number of updated g-particles */
g_updated++; g_updated++;
}
/* Minimal time for next end of time-step */
ti_end_min = min(gp->ti_end, ti_end_min);
ti_end_max = max(gp->ti_end, ti_end_max);
} }
/* Minimal time for next end of time-step */
ti_end_min = min(gp->ti_end, ti_end_min);
ti_end_max = max(gp->ti_end, ti_end_max);
} }
/* Now do the hydro ones... */ /* Now do the hydro ones... */
......
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