Read snapshot units from parameter file. Drift towards output time.

examples/UniformBox/uniformBox.yml

 
 # Define the system of units to use internally. 
-UnitSystem:
+InternalUnitSystem:
   UnitMass_in_cgs:     1   # Grams
   UnitLength_in_cgs:   1   # Centimeters
   UnitVelocity_in_cgs: 1   # Centimeters per second
@@ -23,10 +23,16 @@ TimeIntegration:
 
 # Parameters governing the snapshots
 Snapshots:
-  basename:   uniformBox # Common part of the name of output files
-  time_first: 0.         # Time of the frist output (in internal units)
-  delta_time: 0.01       # Time difference between consecutive outputs (in internal units)
+  basename:            uniformBox # Common part of the name of output files
+  time_first:          0.         # Time of the frist output (in internal units)
+  delta_time:          0.01       # Time difference between consecutive outputs (in internal units)
+  UnitMass_in_cgs:     1   # Grams
+  UnitLength_in_cgs:   1   # Centimeters
+  UnitVelocity_in_cgs: 1   # Centimeters per second
+  UnitCurrent_in_cgs:  1   # Amperes
+  UnitTemp_in_cgs:     1   # Kelvin
 
+  
 # Parameters for the hydrodynamics scheme
 SPH:
   resolution_eta:        1.2349   # Target smoothing length in units of the mean inter-particle separation (1.2349 == 48Ngbs with the cubic spline kernel).

examples/main.c

@@ -291,7 +291,7 @@ int main(int argc, char *argv[]) {
   /* Initialize unit system and constants */
   struct UnitSystem us;
   struct phys_const prog_const;
-  units_init(&us, params);
+  units_init(&us, params, "InternalUnitSystem");
   phys_const_init(&us, &prog_const);
   if (myrank == 0 && verbose > 0) {
     message("Unit system: U_M = %e g.", us.UnitMass_in_cgs);
@@ -423,9 +423,7 @@ int main(int argc, char *argv[]) {
   }
 
   /* Write the state of the system before starting time integration. */
-  char baseName[200];
-  parser_get_param_string(params, "Snapshots:basename", baseName);
-  if (!dry_run) engine_dump_snapshot(&e, &us, baseName);
+  if (!dry_run) engine_dump_snapshot(&e);
 
   /* Now that everything is ready, no need for the parameters any more */
   free(params);
@@ -487,9 +485,6 @@ int main(int argc, char *argv[]) {
     /* Take a step. */
     engine_step(&e);
 
-    /* Snapshot if need be */
-    if (j % 10 == 0) engine_dump_snapshot(&e, &us, baseName);
-
     /* Dump the task data using the given frequency. */
     if (dump_tasks && (dump_tasks == 1 || j % dump_tasks == 1)) {
 #ifdef WITH_MPI
@@ -572,7 +567,7 @@ int main(int argc, char *argv[]) {
 #endif
 
   /* Write final output. */
-  engine_dump_snapshot(&e, &us, baseName);
+  engine_dump_snapshot(&e);
 
 #ifdef WITH_MPI
   if ((res = MPI_Finalize()) != MPI_SUCCESS)

src/engine.c

@@ -2049,6 +2049,25 @@ void engine_step(struct engine *e) {
   }
 #endif
 
+  /* Check for output */
+  while (ti_end_min >= e->ti_nextSnapshot && e->ti_nextSnapshot > 0) {
+
+    e->ti_old = e->ti_current;
+    e->ti_current = e->ti_nextSnapshot;
+    e->time = e->ti_current * e->timeBase + e->timeBegin;
+    e->timeOld = e->ti_old * e->timeBase + e->timeBegin;
+    e->timeStep = (e->ti_current - e->ti_old) * e->timeBase;
+
+    /* Drift everybody to the snapshot position */
+    engine_launch(e, e->nr_threads, 1 << task_type_drift, 0);
+
+    /* Dump... */
+    engine_dump_snapshot(e);
+
+    /* ... and find the next output time */
+    engine_compute_next_snapshot_time(e);
+  }
+
   /* Move forward in time */
   e->ti_old = e->ti_current;
   e->ti_current = ti_end_min;
@@ -2331,35 +2350,32 @@ void engine_split(struct engine *e, struct partition *initial_partition) {
  * @brief Writes a snapshot with the current state of the engine
  *
  * @param e The #engine.
- * @param us The unit system to use for the snapshots.
- * @param baseName The common part of the snapshot file names.
  */
-void engine_dump_snapshot(struct engine *e, struct UnitSystem *us,
-                          const char *baseName) {
+void engine_dump_snapshot(struct engine *e) {
 
   struct clocks_time time1, time2;
   clocks_gettime(&time1);
 
+  if (e->verbose)
+    message("writing snapshot at t=%f", e->time);
+  
 /* Dump... */
 #if defined(WITH_MPI)
 #if defined(HAVE_PARALLEL_HDF5)
-  write_output_parallel(e, baseName, us, myrank, nr_nodes, MPI_COMM_WORLD,
-                        MPI_INFO_NULL);
+  write_output_parallel(e, e->snapshotBaseName, e->snapshotUnits, e->nodeID,
+                        e->nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
 #else
-  write_output_serial(e, baseName, us, myrank, nr_nodes, MPI_COMM_WORLD,
-                      MPI_INFO_NULL);
+  write_output_serial(e, e->snapshotBaseName, e->snapshotUnits, e->nodeID,
+                      e->nr_nodes, MPI_COMM_WORLD, MPI_INFO_NULL);
 #endif
 #else
-  write_output_single(e, baseName, us);
+  write_output_single(e, e->snapshotBaseName, e->snapshotUnits);
 #endif
 
   clocks_gettime(&time2);
   if (e->verbose)
     message("writing particle properties took %.3f %s.",
             (float)clocks_diff(&time1, &time2), clocks_getunit());
-
-  /* ... and find the next output time */
-  engine_compute_next_snapshot_time(e);
 }
 
 #if defined(HAVE_LIBNUMA) && defined(_GNU_SOURCE)
@@ -2431,6 +2447,9 @@ void engine_init(struct engine *e, struct space *s,
   e->deltaTimeSnapshot =
       parser_get_param_double(params, "Snapshots:delta_time");
   e->ti_nextSnapshot = 0;
+  parser_get_param_string(params, "Snapshots:basename", e->snapshotBaseName);
+  e->snapshotUnits = malloc(sizeof(struct UnitSystem));
+  units_init(e->snapshotUnits, params, "Snapshots");
   e->dt_min = parser_get_param_double(params, "TimeIntegration:dt_min");
   e->dt_max = parser_get_param_double(params, "TimeIntegration:dt_max");
   e->file_stats = NULL;
@@ -2617,6 +2636,9 @@ void engine_init(struct engine *e, struct space *s,
         "Time of first snapshot (%e) must be after the simulation start t=%e.",
         e->timeFirstSnapshot, e->timeBegin);
 
+  /* Find the time of the first output */
+  engine_compute_next_snapshot_time(e);
+
 /* Construct types for MPI communications */
 #ifdef WITH_MPI
   part_create_mpi_types();

src/engine.h

@@ -47,6 +47,7 @@
 #include "partition.h"
 #include "physical_constants.h"
 #include "potentials.h"
+#include "units.h"
 
 /* Some constants. */
 enum engine_policy {
@@ -131,10 +132,12 @@ struct engine {
   /* Time base */
   double timeBase;
 
-  /* Snapshot times */
+  /* Snapshot information */
   double timeFirstSnapshot;
   double deltaTimeSnapshot;
   int ti_nextSnapshot;
+  char snapshotBaseName[200];
+  struct UnitSystem *snapshotUnits;
 
   /* File for statistics */
   FILE *file_stats;
@@ -187,8 +190,7 @@ struct engine {
 /* Function prototypes. */
 void engine_barrier(struct engine *e, int tid);
 void engine_compute_next_snapshot_time(struct engine *e);
-void engine_dump_snapshot(struct engine *e, struct UnitSystem *us,
-                          const char *baseName);
+void engine_dump_snapshot(struct engine *e);
 void engine_init(struct engine *e, struct space *s,
                  const struct swift_params *params, int nr_nodes, int nodeID,
                  int nr_threads, int policy, int verbose,

src/units.c

@@ -46,20 +46,23 @@
  *
  * @param us The UnitSystem to initialize.
  * @param params The parsed parameter file.
+ * @param category The section of the parameter file to read from.
  */
-void units_init(struct UnitSystem* us, const struct swift_params* params) {
-
-  us->UnitMass_in_cgs =
-      parser_get_param_double(params, "UnitSystem:UnitMass_in_cgs");
-  us->UnitLength_in_cgs =
-      parser_get_param_double(params, "UnitSystem:UnitLength_in_cgs");
-  const double unitVelocity =
-      parser_get_param_double(params, "UnitSystem:UnitVelocity_in_cgs");
+void units_init(struct UnitSystem* us, const struct swift_params* params,
+                const char* category) {
+
+  char buffer[200];
+  sprintf(buffer, "%s:UnitMass_in_cgs", category);
+  us->UnitMass_in_cgs = parser_get_param_double(params, buffer);
+  sprintf(buffer, "%s:UnitLength_in_cgs", category);
+  us->UnitLength_in_cgs = parser_get_param_double(params, buffer);
+  sprintf(buffer, "%s:UnitVelocity_in_cgs", category);
+  const double unitVelocity = parser_get_param_double(params, buffer);
   us->UnitTime_in_cgs = us->UnitLength_in_cgs / unitVelocity;
-  us->UnitCurrent_in_cgs =
-      parser_get_param_double(params, "UnitSystem:UnitCurrent_in_cgs");
-  us->UnitTemperature_in_cgs =
-      parser_get_param_double(params, "UnitSystem:UnitTemp_in_cgs");
+  sprintf(buffer, "%s:UnitCurrent_in_cgs", category);
+  us->UnitCurrent_in_cgs = parser_get_param_double(params, buffer);
+  sprintf(buffer, "%s:UnitTemp_in_cgs", category);
+  us->UnitTemperature_in_cgs = parser_get_param_double(params, buffer);
 }
 
 /**

src/units.h

@@ -92,7 +92,8 @@ enum UnitConversionFactor {
   UNIT_CONV_TEMPERATURE
 };
 
-void units_init(struct UnitSystem*, const struct swift_params*);
+void units_init(struct UnitSystem*, const struct swift_params*,
+                const char* category);
 double units_get_base_unit(const struct UnitSystem*, enum BaseUnits);
 const char* units_get_base_unit_symbol(enum BaseUnits);
 const char* units_get_base_unit_CGS_symbol(enum BaseUnits);