Statistics are now correct as well.

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]]
 
 #--------------------------------------------------
 

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,

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);
 }
 
 /**

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;
 }