diff --git a/examples/SedovBlast_3D/makeIC.py b/examples/SedovBlast_3D/makeIC.py index 7d1a78188c8bc7034d5cd859c213ab838f52b3e8..e1b743c6cdcd8dcc2f8da14d1d5589fb9ed111f0 100644 --- a/examples/SedovBlast_3D/makeIC.py +++ b/examples/SedovBlast_3D/makeIC.py @@ -54,7 +54,6 @@ u[:] = P0 / (rho0 * (gamma - 1)) # Make the central particles detonate index = argsort(r) u[index[0:N_inject]] = E0 / (N_inject * m[0]) -print ids[index[0:N_inject]] #-------------------------------------------------- diff --git a/src/statistics.c b/src/statistics.c index 85b130f4da6aa95f8c2bf3a12e374431cb9be8d4..55b76eba974df19cdcd1e8ef24e5b6f8d0679c31 100644 --- a/src/statistics.c +++ b/src/statistics.c @@ -104,8 +104,8 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) { const struct part *restrict parts = (struct part *)map_data; const struct xpart *restrict xparts = s->xparts + (ptrdiff_t)(parts - s->parts); - const int ti_current = s->e->ti_current; - const double timeBase = s->e->timeBase; + // const integertime_t ti_current = s->e->ti_current; + // const double timeBase = s->e->timeBase; const double time = s->e->time; struct statistics *const global_stats = data->stats; @@ -126,10 +126,13 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) { const struct gpart *gp = (p->gpart != NULL) ? gp = p->gpart : NULL; /* Get useful variables */ - const integertime_t ti_begin = - get_integer_time_begin(ti_current, p->time_bin); - const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin); - const float dt = (ti_current - (ti_begin + ti_end) / 2) * timeBase; + // const integertime_t ti_begin = + // get_integer_time_begin(ti_current, p->time_bin); + // const integertime_t ti_end = get_integer_time_end(ti_current, + // p->time_bin); + // const integertime_t dti = get_integer_timestep(p->time_bin); + const float dt = + 0.f; //(ti_current - (ti_begin + ti_end) / 2) * timeBase; //MATTHIEU const double x[3] = {p->x[0], p->x[1], p->x[2]}; float a_tot[3] = {p->a_hydro[0], p->a_hydro[1], p->a_hydro[2]}; if (gp != NULL) { @@ -143,6 +146,12 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) { const float m = hydro_get_mass(p); + /* if (p->id == ICHECK) */ + /* message("bin=%d dti=%lld ti_begin=%lld ti_end=%lld dt=%e v=[%e %e %e]", + */ + /* p->time_bin, dti, ti_begin, ti_end, */ + /* dt, v[0], v[1], v[2]); */ + /* Collect mass */ stats.mass += m; @@ -188,8 +197,8 @@ void stats_collect_gpart_mapper(void *map_data, int nr_gparts, const struct index_data *data = (struct index_data *)extra_data; const struct space *s = data->s; const struct gpart *restrict gparts = (struct gpart *)map_data; - const int ti_current = s->e->ti_current; - const double timeBase = s->e->timeBase; + // const integertime_t ti_current = s->e->ti_current; + // const double timeBase = s->e->timeBase; const double time = s->e->time; struct statistics *const global_stats = data->stats; @@ -211,10 +220,12 @@ void stats_collect_gpart_mapper(void *map_data, int nr_gparts, if (gp->id_or_neg_offset < 0) continue; /* Get useful variables */ - const integertime_t ti_begin = - get_integer_time_begin(ti_current, gp->time_bin); - const integertime_t ti_end = get_integer_time_end(ti_current, gp->time_bin); - const float dt = (ti_current - (ti_begin + ti_end) / 2) * timeBase; + // const integertime_t ti_begin = + // get_integer_time_begin(ti_current, gp->time_bin); + // const integertime_t ti_end = get_integer_time_end(ti_current, + // gp->time_bin); + const float dt = + 0.f; //(ti_current - (ti_begin + ti_end) / 2) * timeBase; // MATTHIEU const double x[3] = {gp->x[0], gp->x[1], gp->x[2]}; const float v[3] = {gp->v_full[0] + gp->a_grav[0] * dt, gp->v_full[1] + gp->a_grav[1] * dt, diff --git a/src/timeline.h b/src/timeline.h index 88f0c6f48b781674bbf22798dce2df181e0c1378..b6c411411c845cd2dacce7b336a3c6471577163e 100644 --- a/src/timeline.h +++ b/src/timeline.h @@ -44,7 +44,7 @@ typedef char timebin_t; */ static INLINE integertime_t get_integer_timestep(timebin_t bin) { - if (bin == 0) return 0; + if (bin <= 0) return 0; return 1LL << (bin + 1); } @@ -101,7 +101,7 @@ static INLINE integertime_t get_integer_time_end(integertime_t ti_current, if (dti == 0) return 0; else - return dti * ceill((double)ti_current / (double)dti); + return dti * ceil((double)ti_current / (double)dti); } /** diff --git a/src/timestep.h b/src/timestep.h index c7529bd042eaa20814c14c5f570532026ebdb27b..10e450987e7b02f84243a9001e903e6930d07b2c 100644 --- a/src/timestep.h +++ b/src/timestep.h @@ -42,15 +42,16 @@ make_integer_timestep(float new_dt, timebin_t old_bin, integertime_t ti_current, /* Convert to integer time */ integertime_t new_dti = (integertime_t)(new_dt * timeBase_inv); - if (verbose) message("new_dti=%lld", new_dti); + /* if (verbose) message("new_dti=%lld", new_dti); */ /* Current time-step */ integertime_t current_dti = get_integer_timestep(old_bin); integertime_t ti_end = get_integer_time_end(ti_current, old_bin); - if (verbose) - message("current_dti=%lld old_bin=%d ti_end=%lld", current_dti, old_bin, - ti_end); + /* if (verbose) */ + /* message("current_dti=%lld old_bin=%d ti_end=%lld", current_dti, old_bin, + */ + /* ti_end); */ /* Limit timestep increase */ if (old_bin > 0) new_dti = min(new_dti, 2 * current_dti); @@ -60,14 +61,14 @@ make_integer_timestep(float new_dt, timebin_t old_bin, integertime_t ti_current, while (new_dti < dti_timeline) dti_timeline /= 2LL; new_dti = dti_timeline; - if (verbose) message("new_dti=%lld", new_dti); + /* if (verbose) message("new_dti=%lld", new_dti); */ /* Make sure we are allowed to increase the timestep size */ if (new_dti > current_dti) { if ((max_nr_timesteps - ti_end) % new_dti > 0) new_dti = current_dti; } - if (verbose) message("new_dti=%lld", new_dti); + /* if (verbose) message("new_dti=%lld", new_dti); */ return new_dti; } @@ -150,13 +151,13 @@ __attribute__((always_inline)) INLINE static integertime_t get_part_timestep( new_dt = min(new_dt, e->dt_max); new_dt = max(new_dt, e->dt_min); - if (p->id == ICHECK) message("new_dt=%e", new_dt); + /* if (p->id == ICHECK) message("new_dt=%e", new_dt); */ /* Convert to integer time */ const integertime_t new_dti = make_integer_timestep( new_dt, p->time_bin, e->ti_current, e->timeBase_inv, p->id == ICHECK); - if (p->id == ICHECK) message("new_dti=%lld", new_dti); + /* if (p->id == ICHECK) message("new_dti=%lld", new_dti); */ return new_dti; }