Start working on logger

examples/main.c

@@ -85,6 +85,8 @@ void print_help_message(void) {
   printf("  %2s %14s %s\n", "-g", "",
          "Run with an external gravitational potential.");
   printf("  %2s %14s %s\n", "-G", "", "Run with self-gravity.");
+  printf("  %2s %14s %s\n", "-l", "", "Use logger for the output "
+	 "(Optimized binary snapshot).");
   printf("  %2s %14s %s\n", "-M", "",
          "Reconstruct the multipoles every time-step.");
   printf("  %2s %14s %s\n", "-n", "{int}",
@@ -188,6 +190,7 @@ int main(int argc, char *argv[]) {
   int nsteps = -2;
   int restart = 0;
   int with_cosmology = 0;
+  int with_logger = 0;
   int with_external_gravity = 0;
   int with_sourceterms = 0;
   int with_cooling = 0;

examples/parameter_example.yml

@@ -90,6 +90,7 @@ Snapshots:
   UnitTemp_in_cgs:     1  # (Optional) Unit system for the outputs (Kelvin)
   output_list_on:      0  # (Optional) Enable the output list
   output_list:         snaplist.txt # (Optional) File containing the output times (see documentation in "Parameter File" section)
+  logger_max_steps:    10 # (Optional) Number of particle steps between two chunk writing if using logger
 
 # Parameters governing the conserved quantities statistics
 Statistics:

src/active.h

@@ -443,4 +443,44 @@ __attribute__((always_inline)) INLINE static int spart_is_starting(
 
   return (spart_bin <= max_active_bin);
 }
+
+/**
+ * @brief Should this particle write its data now ?
+ *
+ * @param p The #part.
+ * @param e The #engine containing information about the current time.
+ * @return 1 if the #part should write, 0 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int part_should_write(
+    const struct part *p, const struct engine *e) {
+
+  return (p->last_output > e->logger_max_steps);  
+}
+
+/**
+ * @brief Should this particle write its data now ?
+ *
+ * @param p The #gpart.
+ * @param e The #engine containing information about the current time.
+ * @return 1 if the #gpart should write, 0 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int gpart_should_write(
+    const struct gpart *p, const struct engine *e) {
+
+  return (p->last_output > e->logger_max_steps);  
+}
+
+/**
+ * @brief Should this particle write its data now ?
+ *
+ * @param p The #spart.
+ * @param e The #engine containing information about the current time.
+ * @return 1 if the #spart should write, 0 otherwise.
+ */
+__attribute__((always_inline)) INLINE static int spart_should_write(
+    const struct spart *p, const struct engine *e) {
+
+  return (p->last_output > e->logger_max_steps);  
+}
+
 #endif /* SWIFT_ACTIVE_H */

src/cell.c

@@ -2904,6 +2904,7 @@ int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s) {
     if (c->grav.down_in != NULL) scheduler_activate(s, c->grav.down_in);
     if (c->grav.mesh != NULL) scheduler_activate(s, c->grav.mesh);
     if (c->grav.long_range != NULL) scheduler_activate(s, c->grav.long_range);
+    if (c->logger != NULL) scheduler_activate(s, c->logger);
   }
 
   return rebuild;

src/cell.h

@@ -279,6 +279,12 @@ struct cell {
     /*! Values of h_max before the drifts, used for sub-cell tasks. */
     float h_max_old;
 
+    /*! The logger task */
+    struct task *logger;
+
+    /*! The task to compute time-steps */
+    struct task *timestep;
+
     /*! Values of dx_max before the drifts, used for sub-cell tasks. */
     float dx_max_part_old;
 

src/engine.c

@@ -112,7 +112,9 @@ const char *engine_policy_names[] = {"none",
                                      "stars",
                                      "structure finding",
                                      "star formation",
-                                     "feedback"};
+                                     "feedback",
+				     "logger"
+};
 
 /** The rank of the engine as a global variable (for messages). */
 int engine_rank;
@@ -217,6 +219,7 @@ void engine_make_hierarchical_tasks_common(struct engine *e, struct cell *c) {
   struct scheduler *s = &e->sched;
   const int is_with_cooling = (e->policy & engine_policy_cooling);
   const int is_with_star_formation = (e->policy & engine_policy_star_formation);
+  const int is_logger = (e->policy & engine_policy_logger);
 
   /* Are we in a super-cell ? */
   if (c->super == c) {
@@ -228,6 +231,11 @@ void engine_make_hierarchical_tasks_common(struct engine *e, struct cell *c) {
       c->kick1 = scheduler_addtask(s, task_type_kick1, task_subtype_none, 0, 0,
                                    c, NULL);
 
+      if (is_logger)
+	c->logger = scheduler_addtask(s, task_type_logger, task_subtype_none, 0, 0,
+				      c, NULL);
+	
+
       c->kick2 = scheduler_addtask(s, task_type_kick2, task_subtype_none, 0, 0,
                                    c, NULL);
 
@@ -5589,7 +5597,10 @@ void engine_check_for_dumps(struct engine *e) {
 
         /* Dump everything */
         engine_print_stats(e);
-        engine_dump_snapshot(e);
+	if (e->policy & engine_policy_logger)
+	  engine_dump_index(e);
+	else
+	  engine_dump_snapshot(e);
 
       } else if (e->ti_next_stats < e->ti_next_snapshot) {
 
@@ -5619,7 +5630,10 @@ void engine_check_for_dumps(struct engine *e) {
         engine_drift_all(e);
 
         /* Dump snapshot */
-        engine_dump_snapshot(e);
+	if (e->policy & engine_policy_logger)
+	  engine_dump_index(e);
+	else
+	  engine_dump_snapshot(e);
 
       } else if (e->ti_next_stats > e->ti_next_snapshot) {
 
@@ -5637,7 +5651,10 @@ void engine_check_for_dumps(struct engine *e) {
         engine_drift_all(e);
 
         /* Dump snapshot */
-        engine_dump_snapshot(e);
+	if (e->policy & engine_policy_logger)
+	  engine_dump_index(e);
+	else
+	  engine_dump_snapshot(e);
 
         /* Let's fake that we are at the stats dump time */
         e->ti_current = e->ti_next_stats;
@@ -5672,7 +5689,10 @@ void engine_check_for_dumps(struct engine *e) {
       engine_drift_all(e);
 
       /* Dump... */
-      engine_dump_snapshot(e);
+      if (e->policy & engine_policy_logger)
+	engine_dump_index(e);
+      else
+	engine_dump_snapshot(e);
 
       /* ... and find the next output time */
       engine_compute_next_snapshot_time(e);
@@ -6447,6 +6467,30 @@ void engine_dump_snapshot(struct engine *e) {
             (float)clocks_diff(&time1, &time2), clocks_getunit());
 }
 
+/**
+ * @brief Writes an index file with the current state of the engine
+ *
+ * @param e The #engine.
+ */
+void engine_dump_index(struct engine *e) {
+
+  struct clocks_time time1, time2;
+  clocks_gettime(&time1);
+
+  if (e->verbose) message("writing index at t=%e.", e->time);
+
+  /* Dump... */
+  /* Should use snapshotBaseName for the index name */
+  message("I am dumping an index file...");
+
+  e->dump_snapshot = 0;
+
+  clocks_gettime(&time2);
+  if (e->verbose)
+    message("writing particle properties took %.3f %s.",
+            (float)clocks_diff(&time1, &time2), clocks_getunit());
+}
+
 #ifdef HAVE_SETAFFINITY
 /**
  * @brief Returns the initial affinity the main thread is using.
@@ -6517,6 +6561,7 @@ void engine_unpin(void) {
  * @param Nstars total number of star particles in the simulation.
  * @param policy The queuing policy to use.
  * @param verbose Is this #engine talkative ?
+ * @param logger_max_steps Max number of particle steps before writing with logger
  * @param reparttype What type of repartition algorithm are we using ?
  * @param internal_units The system of units used internally.
  * @param physical_constants The #phys_const used for this run.
@@ -6612,6 +6657,14 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
   e->last_repartition = 0;
 #endif
 
+
+  /* Logger params */
+  char logger_name_file[PARSER_MAX_LINE_SIZE];
+  e->logger_max_steps = parser_get_opt_param_int(params, "Snapshots:logger_max_steps", 10);
+  parser_get_opt_param_string(params, "Snapshots:dump_file", logger_name_file, "dump.smew");
+  e->logger_dump = malloc(sizeof(struct dump));
+  dump_init(e->logger_dump, logger_name_file, 1024 * 1024 * 10);
+
   /* Make the space link back to the engine. */
   s->e = e;
 
@@ -7069,6 +7122,11 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
 #endif
 
   /* Find the time of the first snapshot  output */
+  if (e->policy & engine_policy_logger)
+    if (e->nodeID == 0)
+      message("Expected output of over 9000\n Should write a real message...");
+
+  /* Find the time of the first output */
   engine_compute_next_snapshot_time(e);
 
   /* Find the time of the first statistics output */
@@ -7631,6 +7689,7 @@ void engine_clean(struct engine *e) {
 
   free(e->links);
   free(e->cell_loc);
+  free(e->logger_dump);
   scheduler_clean(&e->sched);
   space_clean(e->s);
   threadpool_clean(&e->threadpool);

src/engine.h

@@ -38,6 +38,7 @@
 #include "clocks.h"
 #include "collectgroup.h"
 #include "cooling_struct.h"
+#include "dump.h"
 #include "gravity_properties.h"
 #include "mesh_gravity.h"
 #include "parser.h"
@@ -74,9 +75,10 @@ enum engine_policy {
   engine_policy_stars = (1 << 15),
   engine_policy_structure_finding = (1 << 16),
   engine_policy_star_formation = (1 << 17),
-  engine_policy_feedback = (1 << 18)
+  engine_policy_feedback = (1 << 18),
+  engine_policy_logger = (1 << 19)
 };
-#define engine_maxpolicy 18
+#define engine_maxpolicy 19
 extern const char *engine_policy_names[];
 
 /**
@@ -312,6 +314,12 @@ struct engine {
   int forcerepart;
   struct repartition *reparttype;
 
+  /* Number of particle steps between dumping a chunk of data */
+  int logger_max_steps;
+
+  /* File name of the dump file */
+  struct dump *logger_dump;
+
   /* How many steps have we done with the same set of tasks? */
   int tasks_age;
 
@@ -415,6 +423,7 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
 void engine_config(int restart, struct engine *e, struct swift_params *params,
                    int nr_nodes, int nodeID, int nr_threads, int with_aff,
                    int verbose, const char *restart_file);
+void engine_dump_index(struct engine *e);
 void engine_launch(struct engine *e);
 void engine_prepare(struct engine *e);
 void engine_init_particles(struct engine *e, int flag_entropy_ICs,

src/gravity/Default/gravity.h

@@ -227,6 +227,11 @@ __attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
   gp->time_bin = 0;
 
   gravity_init_gpart(gp);
+
+#ifdef WITH_LOGGER
+  gp->last_output = 0;
+  gp->last_offset = 0;
+#endif
 }
 
 #endif /* SWIFT_DEFAULT_GRAVITY_H */

src/gravity/Default/gravity_part.h

@@ -75,6 +75,14 @@ struct gpart {
   double potential_exact;
 #endif
 
+#ifdef WITH_LOGGER
+  /* Number of time step since last output */
+  short int last_output;
+  
+  /* offset at last writing */
+  size_t last_offset;
+#endif
+
 } SWIFT_STRUCT_ALIGN;
 
 #endif /* SWIFT_DEFAULT_GRAVITY_PART_H */

src/hydro/Gadget2/hydro.h

@@ -743,6 +743,8 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
     struct part *restrict p, struct xpart *restrict xp) {
 
   p->time_bin = 0;
+  p->last_output = 0;
+  p->last_offset = 0;
   xp->v_full[0] = p->v[0];
   xp->v_full[1] = p->v[1];
   xp->v_full[2] = p->v[2];

src/hydro/Gadget2/hydro_part.h

@@ -90,6 +90,12 @@ struct part {
   /* Entropy time derivative */
   float entropy_dt;
 
+  /* Number of time step since last output */
+  short int last_output;
+
+  /* offset at last writing */
+  size_t last_offset;
+
   union {
 
     struct {

src/runner.c

@@ -53,6 +53,7 @@
 #include "hydro.h"
 #include "hydro_properties.h"
 #include "kick.h"
+#include "logger.h"
 #include "minmax.h"
 #include "runner_doiact_vec.h"
 #include "scheduler.h"
@@ -2610,6 +2611,9 @@ void *runner_main(void *data) {
         case task_type_end_force:
           runner_do_end_force(r, ci, 1);
           break;
+        case task_type_logger:
+	  runner_do_logger(r, ci, 1);
+	  break;
         case task_type_timestep:
           runner_do_timestep(r, ci, 1);
           break;
@@ -2688,3 +2692,112 @@ void *runner_main(void *data) {
   /* Be kind, rewind. */
   return NULL;
 }
+
+
+/**
+ * @brief Write the required particles through the logger.
+ *
+ * @param r The runner thread.
+ * @param c The cell.
+ * @param timer Are we timing this ?
+ */
+void runner_do_logger(struct runner *r, struct cell *c, int timer) {
+
+  const struct engine *e = r->e;
+  struct part *restrict parts = c->parts;
+  //struct xpart *restrict xparts = c->xparts;
+  struct gpart *restrict gparts = c->gparts;
+  struct spart *restrict sparts = c->sparts;
+  const int count = c->count;
+  const int gcount = c->gcount;
+  const int scount = c->scount;
+  //const integertime_t ti_current = e->ti_current;
+  //const double timeBase = e->timeBase;
+
+  TIMER_TIC;
+
+  /* Anything to do here? */
+  if (!cell_is_starting(c, e)) return;
+
+  /* Recurse? */
+  if (c->split) {
+    for (int k = 0; k < 8; k++)
+      if (c->progeny[k] != NULL) runner_do_logger(r, c->progeny[k], 0);
+  } else {
+
+    /* Loop over the parts in this cell. */
+    for (int k = 0; k < count; k++) {
+
+      /* Get a handle on the part. */
+      struct part *restrict p = &parts[k];
+      //struct xpart *restrict xp = &xparts[k];
+
+      /* If particle needs to be kicked */
+      if (part_is_starting(p, e)) {
+
+	if (part_should_write(p, e))
+	  {
+	    /* Write particle */
+	    logger_log_part(p, logger_mask_x | logger_mask_v | logger_mask_a |
+			    logger_mask_u | logger_mask_h | logger_mask_rho |
+			    logger_mask_consts,
+			    &p->last_offset, e->logger_dump);
+	    message("Offset: %lu", p->last_offset);
+	    /* Set counter back to zero */
+	    p->last_output = 0;
+	  }
+	else
+	  /* Update counter */
+	  p->last_output += 1;
+      }
+    }
+
+    /* Loop over the gparts in this cell. */
+    for (int k = 0; k < gcount; k++) {
+
+      /* Get a handle on the part. */
+      struct gpart *restrict gp = &gparts[k];
+
+      /* If particle needs to be kicked */
+      if (gpart_is_starting(gp, e)) {
+
+	if (gpart_should_write(gp, e))
+	  {
+	    /* Write particle */
+	    logger_log_gpart(gp, logger_mask_x | logger_mask_v | logger_mask_a |
+			     logger_mask_h | logger_mask_consts,
+			     &gp->last_offset, e->logger_dump);
+	    /* Set counter back to zero */
+	    gp->last_output = 0;
+	  }
+	else
+	  {
+	    /* Update counter */
+	    gp->last_output += 1;
+	  }
+      }
+    }
+
+    /* Loop over the star particles in this cell. */
+    for (int k = 0; k < scount; k++) {
+
+      /* Get a handle on the s-part. */
+      struct spart *restrict sp = &sparts[k];
+
+      /* If particle needs to be kicked */
+      if (spart_is_starting(sp, e)) {
+
+	if (spart_should_write(sp, e))
+	  {
+	    message("I am writing sparticle %lli", sp->id);
+	    error("Not implemented");
+	    sp->last_output = 0;
+	  }
+	else
+	  sp->last_output += 1;
+      }
+    }
+  }
+
+  if (timer) TIMER_TOC(timer_logger);
+}

src/runner.h

@@ -80,6 +80,7 @@ void runner_do_init(struct runner *r, struct cell *c, int timer);
 void runner_do_cooling(struct runner *r, struct cell *c, int timer);
 void runner_do_grav_external(struct runner *r, struct cell *c, int timer);
 void runner_do_grav_fft(struct runner *r, int timer);
+void runner_do_logger(struct runner *r, struct cell *c, int timer);
 void *runner_main(void *data);
 void runner_do_unskip_mapper(void *map_data, int num_elements,
                              void *extra_data);

src/scheduler.c

@@ -2123,6 +2123,7 @@ void scheduler_enqueue(struct scheduler *s, struct task *t) {
       case task_type_kick1:
       case task_type_kick2:
       case task_type_stars_ghost:
+      case task_type_logger:
       case task_type_timestep:
         qid = t->ci->super->owner;
         break;

src/stars/Default/stars_part.h

@@ -62,6 +62,9 @@ struct spart {
 
   } density;
 
+  /* Number of time step since last output */
+  short int last_output;
+    
 #ifdef SWIFT_DEBUG_CHECKS
 
   /* Time of the last drift */

src/task.c

@@ -58,7 +58,8 @@ const char *taskID_names[task_type_count] = {
     "grav_mm",        "grav_down_in",   "grav_down",
     "grav_mesh",      "cooling",        "star_formation",
     "sourceterms",    "stars_ghost_in", "stars_ghost",
-    "stars_ghost_out"};
+    "stars_ghost_out","logger"
+};
 
 /* Sub-task type names. */
 const char *subtaskID_names[task_subtype_count] = {
@@ -161,6 +162,7 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
     case task_type_end_force:
     case task_type_kick1:
     case task_type_kick2:
+    case task_type_logger:
     case task_type_timestep:
     case task_type_send:
     case task_type_recv:
@@ -304,6 +306,7 @@ void task_unlock(struct task *t) {
 
     case task_type_end_force:
     case task_type_kick1:
+    case task_type_logger:
     case task_type_kick2:
     case task_type_timestep:
       cell_unlocktree(ci);
@@ -400,6 +403,7 @@ int task_lock(struct task *t) {
     case task_type_end_force:
     case task_type_kick1:
     case task_type_kick2:
+    case task_type_logger:
     case task_type_timestep:
       if (ci->hydro.hold || ci->grav.phold) return 0;
       if (cell_locktree(ci) != 0) return 0;

src/task.h

@@ -70,6 +70,7 @@ enum task_types {
   task_type_stars_ghost_in,
   task_type_stars_ghost,
   task_type_stars_ghost_out,
+  task_type_logger,
   task_type_count
 } __attribute__((packed));
 

src/timers.c

@@ -94,6 +94,7 @@ const char* timers_names[timer_count] = {
     "dopair_subset_stars_density",
     "dosubpair_stars_density",
     "dosub_self_stars_density",
+    "logger"
 };
 
 /* File to store the timers */

src/timers.h

@@ -95,6 +95,7 @@ enum {
   timer_dopair_subset_stars_density,
   timer_dosubpair_stars_density,
   timer_dosub_self_stars_density,
+  timer_logger,
   timer_count,
 };