Particles only carry their time bin.

src/Makefile.am

@@ -60,7 +60,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
 nobase_noinst_HEADERS = align.h approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
 		 kernel_long_gravity.h vector.h runner_doiact.h runner_doiact_grav.h runner_doiact_fft.h \
                  units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
-		 dimension.h equation_of_state.h active.h \
+		 dimension.h equation_of_state.h active.h timeline.h \
 		 gravity.h gravity_io.h \
 		 gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
 		 gravity/Default/gravity_debug.h gravity/Default/gravity_part.h  \

src/active.h

@@ -27,6 +27,7 @@
 #include "const.h"
 #include "engine.h"
 #include "part.h"
+#include "timeline.h"
 
 /**
  * @brief Check that a cell been drifted to the current time.
@@ -107,7 +108,8 @@ __attribute__((always_inline)) INLINE static int part_is_active(
           p->ti_end, e->ti_current);
 #endif
 
-  return (p->ti_end == e->ti_current);
+  // return (p->ti_end == e->ti_current);
+  return (p->time_bin < get_max_active_bin(e->ti_current));
 }
 
 /**
@@ -126,7 +128,8 @@ __attribute__((always_inline)) INLINE static int gpart_is_active(
         gp->ti_end, e->ti_current);
 #endif
 
-  return (gp->ti_end == e->ti_current);
+  // return (gp->ti_end == e->ti_current);
+  return (gp->time_bin < get_max_active_bin(e->ti_current));
 }
 
 #endif /* SWIFT_ACTIVE_H */

src/engine.h

@@ -82,9 +82,6 @@ extern const char *engine_policy_names[];
 /* The rank of the engine as a global variable (for messages). */
 extern int engine_rank;
 
-/* The maximal number of timesteps in a simulation */
-#define max_nr_timesteps (1 << 28)
-
 /* Data structure for the engine. */
 struct engine {
 

src/gravity/Default/gravity.h

@@ -53,8 +53,7 @@ gravity_compute_timestep_self(const struct gpart* const gp) {
 __attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
     struct gpart* gp) {
 
-  gp->ti_begin = 0;
-  gp->ti_end = 0;
+  gp->time_bin = -1;
   gp->epsilon = 0.;  // MATTHIEU
 }
 

src/gravity/Default/gravity_part.h

@@ -19,12 +19,13 @@
 #ifndef SWIFT_DEFAULT_GRAVITY_PART_H
 #define SWIFT_DEFAULT_GRAVITY_PART_H
 
-/* Some standard headers. */
-#include <stdlib.h>
-
 /* Gravity particle. */
 struct gpart {
 
+  /* Particle ID. If negative, it is the negative offset of the #part with
+     which this gpart is linked. */
+  long long id_or_neg_offset;
+
   /* Particle position. */
   double x[3];
 
@@ -44,14 +45,13 @@ struct gpart {
   float epsilon;
 
   /* Particle time of beginning of time-step. */
-  int ti_begin;
+  // int ti_begin;
 
   /* Particle time of end of time-step. */
-  int ti_end;
+  // int ti_end;
 
-  /* Particle ID. If negative, it is the negative offset of the #part with
-     which this gpart is linked. */
-  long long id_or_neg_offset;
+  /* Time-step length */
+  timebin_t time_bin;
 
 } SWIFT_STRUCT_ALIGN;
 

src/hydro/Gadget2/hydro.h

@@ -212,8 +212,7 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
 __attribute__((always_inline)) INLINE static void hydro_first_init_part(
     struct part *restrict p, struct xpart *restrict xp) {
 
-  p->ti_begin = 0;
-  p->ti_end = 0;
+  p->time_bin = -1;
   xp->v_full[0] = p->v[0];
   xp->v_full[1] = p->v[1];
   xp->v_full[2] = p->v[2];
@@ -303,7 +302,10 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   const float abs_div_v = fabsf(p->density.div_v);
 
   /* Compute the pressure */
-  const float half_dt = (ti_current - (p->ti_begin + p->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 float half_dt = (ti_current - (ti_begin + ti_end) / 2) * timeBase;
   const float pressure = hydro_get_pressure(p, half_dt);
 
   /* Compute the sound speed */
@@ -383,7 +385,9 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
     p->rho *= expf(w2);
 
   /* Drift the pressure */
-  const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
+  const integertime_t ti_begin = get_integer_time_begin(t1, p->time_bin);
+  const integertime_t ti_end = get_integer_time_end(t1, p->time_bin);
+  const float dt_entr = (t1 - (ti_begin + ti_end) / 2) * timeBase;
   const float pressure = hydro_get_pressure(p, dt_entr);
 
   /* Compute the new sound speed */

src/hydro/Gadget2/hydro_part.h

@@ -50,6 +50,12 @@ struct xpart {
 /* Data of a single particle. */
 struct part {
 
+  /* Particle ID. */
+  long long id;
+
+  /* Pointer to corresponding gravity part. */
+  struct gpart* gpart;
+
   /* Particle position. */
   double x[3];
 
@@ -66,10 +72,10 @@ struct part {
   float mass;
 
   /* Particle time of beginning of time-step. */
-  int ti_begin;
+  // int ti_begin;
 
   /* Particle time of end of time-step. */
-  int ti_end;
+  // int ti_end;
 
   /* Particle density. */
   float rho;
@@ -124,11 +130,8 @@ struct part {
     } force;
   };
 
-  /* Particle ID. */
-  long long id;
-
-  /* Pointer to corresponding gravity part. */
-  struct gpart* gpart;
+  /* Time-step length */
+  timebin_t time_bin;
 
 } SWIFT_STRUCT_ALIGN;
 

src/kick.h

@@ -25,6 +25,7 @@
 /* Local headers. */
 #include "const.h"
 #include "debug.h"
+#include "timeline.h"
 
 /**
  * @brief Perform the 'kick' operation on a #gpart
@@ -68,8 +69,10 @@ __attribute__((always_inline)) INLINE static void kick_part(
     double timeBase) {
 
   /* Compute the time step for this kick */
-  const int ti_start = (p->ti_begin + p->ti_end) / 2;
-  const int ti_end = p->ti_end + new_dti / 2;
+  const integertime_t ti_begin = get_integer_time_begin(t1, p->time_bin);
+  const integertime_t ti_end = get_integer_time_end(t1, p->time_bin);
+  const int ti_start = (ti_begin + ti_end) / 2;
+  const int ti_end = ti_end + new_dti / 2;
   const float dt = (ti_end - ti_start) * timeBase;
   const float half_dt = (ti_end - p->ti_end) * timeBase;
 

src/part.h

@@ -33,6 +33,7 @@
 /* Local headers. */
 #include "align.h"
 #include "const.h"
+#include "timeline.h"
 
 /* Some constants. */
 #define part_align 128

src/runner.c

@@ -1029,7 +1029,7 @@ void runner_do_kick(struct runner *r, struct cell *c, int timer) {
 void runner_do_recv_cell(struct runner *r, struct cell *c, int timer) {
 
   const struct part *restrict parts = c->parts;
-  const struct gpart *restrict gparts = c->gparts;
+  // const struct gpart *restrict gparts = c->gparts;
   const size_t nr_parts = c->count;
   const size_t nr_gparts = c->gcount;
   // const int ti_current = r->e->ti_current;
@@ -1045,14 +1045,14 @@ void runner_do_recv_cell(struct runner *r, struct cell *c, int timer) {
 
     /* Collect everything... */
     for (size_t k = 0; k < nr_parts; k++) {
-      const int ti_end = parts[k].ti_end;
+      const int ti_end = 0;  // parts[k].ti_end; //MATTHIEU
       // if(ti_end < ti_current) error("Received invalid particle !");
       ti_end_min = min(ti_end_min, ti_end);
       ti_end_max = max(ti_end_max, ti_end);
       h_max = max(h_max, parts[k].h);
     }
     for (size_t k = 0; k < nr_gparts; k++) {
-      const int ti_end = gparts[k].ti_end;
+      const int ti_end = 0;  // gparts[k].ti_end; //MATTHIEU
       // if(ti_end < ti_current) error("Received invalid particle !");
       ti_end_min = min(ti_end_min, ti_end);
       ti_end_max = max(ti_end_max, ti_end);

src/space.c

@@ -1529,7 +1529,7 @@ void space_split_mapper(void *map_data, int num_cells, void *extra_data) {
         struct part *p = &parts[k];
         struct xpart *xp = &xparts[k];
         const float h = p->h;
-        const int ti_end = p->ti_end;
+        const int ti_end = get_integer_time_end(e->ti_current, p->time_bin);
         xp->x_diff[0] = 0.f;
         xp->x_diff[1] = 0.f;
         xp->x_diff[2] = 0.f;
@@ -1539,7 +1539,7 @@ void space_split_mapper(void *map_data, int num_cells, void *extra_data) {
       }
       for (int k = 0; k < gcount; k++) {
         struct gpart *gp = &gparts[k];
-        const int ti_end = gp->ti_end;
+        const int ti_end = get_integer_time_end(e->ti_current, gp->time_bin);
         gp->x_diff[0] = 0.f;
         gp->x_diff[1] = 0.f;
         gp->x_diff[2] = 0.f;

src/timeline.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_TIMELINE_H
+#define SWIFT_TIMELINE_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local headers. */
+#include "inline.h"
+
+typedef long long integertime_t;
+typedef char timebin_t;
+
+/*! The number of time bins */
+#define num_time_bins 20
+
+/*! The maximal number of timesteps in a simulation */
+#define max_nr_timesteps (1LL << (num_time_bins + 1))
+
+/**
+ * @brief Returns the integer time interval corresponding to a time bin
+ *
+ * @param bin The time bin of interest.
+ */
+static INLINE integertime_t get_integer_timestep(timebin_t bin) {
+
+  if (bin == 0) return 0;
+  return 1LL << (bin + 1);
+}
+
+/**
+ * @brief Returns the physical time interval corresponding to a time bin.
+ *
+ * @param bin The time bin of interest.
+ * @param timeBase the minimal time-step size of the simulation.
+ */
+static INLINE double get_timestep(timebin_t bin, double timeBase) {
+
+  return get_integer_timestep(bin) * timeBase;
+}
+
+/**
+ * @brief Returns the integer time corresponding to the start of the time-step
+ * given by a time-bin.
+ *
+ * @param ti_current The current time on the integer time line.
+ * @param bin The time bin of interest.
+ */
+static INLINE integertime_t get_integer_time_begin(integertime_t ti_current,
+                                                   timebin_t bin) {
+
+  const integertime_t dti = get_integer_timestep(bin);
+  return ti_current % dti;
+}
+
+/**
+ * @brief Returns the integer time corresponding to the start of the time-step
+ * given by a time-bin.
+ *
+ * @param ti_current The current time on the integer time line.
+ * @param bin The time bin of interest.
+ */
+static INLINE integertime_t get_integer_time_end(integertime_t ti_current,
+                                                 timebin_t bin) {
+
+  const integertime_t dti = get_integer_timestep(bin);
+  return ti_current % dti + dti;
+}
+
+/**
+ * @brief Returns the highest active time bin at a given point on the time line.
+ *
+ * @param time The current point on the time line.
+ */
+static INLINE timebin_t get_max_active_bin(integertime_t time) {
+
+  if (time == 0) return num_time_bins;
+
+  timebin_t bin = 1;
+  while (!((1LL << (bin + 1)) & time)) ++bin;
+
+  return bin;
+}
+
+#endif /* SWIFT_TIMELINE_H */