Time base in snap header

doc/RTD/source/Snapshots/index.rst

@@ -33,6 +33,15 @@ The ``RunName`` field contains the name of the simulation that was specified as
 the ``run_name`` in the :ref:`Parameters_meta_data` section of the YAML
 parameter file.
 
+The ``TimeBase_dloga`` field contains the change in logarithm of the
+scale-factor corresponding to a time-step of length 1 on the integer
+time-line. This is the smallest time-step size that the code can use. This field
+is zero in non-cosmological runs. Similarly, the field ``TimeBase_dt`` contains
+the smallest time-step size (in internal units) that the code can take. This
+would be the increase in time a particle in the time-bin one would have. Note
+that in cosmological runs this quantity evolves with redhsift as the (logarithm
+of the) scale-factor is used on the integer time-line.
+
 Meta-data about the code and run
 --------------------------------
 

src/cosmology.c

@@ -853,6 +853,54 @@ double cosmology_get_delta_time_from_scale_factors(const struct cosmology *c,
   return t2 - t1;
 }
 
+/**
+ * @brief Compute the time corresponding to the timebase interval at the current
+ * redshift
+ *
+ * This function is slow as it performs the actual integral.
+ *
+ * @param c The comology model.
+ * @param ti_current The current point on the time-line.
+ *
+ * @return The time corresponding to c->time_base in internal time units.
+ */
+double cosmology_get_timebase(struct cosmology *c,
+                              const integertime_t ti_current) {
+
+#ifdef HAVE_LIBGSL
+
+  /* We are going to average over a number of time-bins
+   * as in some cases the smallest bin is too small for double accuracy */
+  const int num_bins = 24;
+
+  /* Range in log(a) */
+  const double log_a_start = c->log_a_begin + (ti_current)*c->time_base;
+  const double log_a_end =
+      c->log_a_begin + (ti_current + (1LL << num_bins)) * c->time_base;
+
+  /* Range in (a) */
+  const double a_start = exp(log_a_start);
+  const double a_end = exp(log_a_end);
+
+  /* Initalise the GSL workspace and function */
+  gsl_integration_workspace *space =
+      gsl_integration_workspace_alloc(GSL_workspace_size);
+  gsl_function F = {&time_integrand, c};
+
+  /* Perform the integral */
+  double result, abserr;
+  gsl_integration_qag(&F, a_start, a_end, 0, 1.0e-10, GSL_workspace_size,
+                      GSL_INTEG_GAUSS61, space, &result, &abserr);
+
+  result /= (1LL << num_bins);
+
+  return result;
+
+#else
+  error("Code not compiled with GSL. Can't compute cosmology integrals.");
+#endif
+}
+
 /**
  * @brief Compute scale factor from time since big bang (in internal units).
  *

src/cosmology.h

@@ -213,6 +213,9 @@ double cosmology_get_delta_time_from_scale_factors(const struct cosmology *c,
                                                    const double a_start,
                                                    const double a_end);
 
+double cosmology_get_timebase(struct cosmology *c,
+                              const integertime_t ti_current);
+
 double cosmology_get_scale_factor(const struct cosmology *cosmo, double t);
 
 double cosmology_get_time_since_big_bang(const struct cosmology *c, double a);

src/distributed_io.c

@@ -394,6 +394,16 @@ void write_output_distributed(struct engine* e,
   io_write_attribute_s(h_grp, "Code", "SWIFT");
   io_write_attribute_s(h_grp, "RunName", e->run_name);
 
+  /* Write out the time-base */
+  if (with_cosmology) {
+    io_write_attribute_d(h_grp, "TimeBase_dloga", e->time_base);
+    const double delta_t = cosmology_get_timebase(e->cosmology, e->ti_current);
+    io_write_attribute_d(h_grp, "TimeBase_dt", delta_t);
+  } else {
+    io_write_attribute_d(h_grp, "TimeBase_dloga", 0);
+    io_write_attribute_d(h_grp, "TimeBase_dt", e->time_base);
+  }
+
   /* Store the time at which the snapshot was written */
   time_t tm = time(NULL);
   struct tm* timeinfo = localtime(&tm);

src/parallel_io.c

@@ -1175,6 +1175,16 @@ void prepare_file(struct engine* e, const char* fileName,
   io_write_attribute_s(h_grp, "Code", "SWIFT");
   io_write_attribute_s(h_grp, "RunName", e->run_name);
 
+  /* Write out the time-base */
+  if (with_cosmology) {
+    io_write_attribute_d(h_grp, "TimeBase_dloga", e->time_base);
+    const double delta_t = cosmology_get_timebase(e->cosmology, e->ti_current);
+    io_write_attribute_d(h_grp, "TimeBase_dt", delta_t);
+  } else {
+    io_write_attribute_d(h_grp, "TimeBase_dloga", 0);
+    io_write_attribute_d(h_grp, "TimeBase_dt", e->time_base);
+  }
+
   /* Store the time at which the snapshot was written */
   time_t tm = time(NULL);
   struct tm* timeinfo = localtime(&tm);

src/serial_io.c

@@ -1058,6 +1058,17 @@ void write_output_serial(struct engine* e,
     io_write_attribute_s(h_grp, "Code", "SWIFT");
     io_write_attribute_s(h_grp, "RunName", e->run_name);
 
+    /* Write out the time-base */
+    if (with_cosmology) {
+      io_write_attribute_d(h_grp, "TimeBase_dloga", e->time_base);
+      const double delta_t =
+          cosmology_get_timebase(e->cosmology, e->ti_current);
+      io_write_attribute_d(h_grp, "TimeBase_dt", delta_t);
+    } else {
+      io_write_attribute_d(h_grp, "TimeBase_dloga", 0);
+      io_write_attribute_d(h_grp, "TimeBase_dt", e->time_base);
+    }
+
     /* Store the time at which the snapshot was written */
     time_t tm = time(NULL);
     struct tm* timeinfo = localtime(&tm);

src/single_io.c

@@ -893,6 +893,16 @@ void write_output_single(struct engine* e,
   io_write_attribute_s(h_grp, "Code", "SWIFT");
   io_write_attribute_s(h_grp, "RunName", e->run_name);
 
+  /* Write out the time-base */
+  if (with_cosmology) {
+    io_write_attribute_d(h_grp, "TimeBase_dloga", e->time_base);
+    const double delta_t = cosmology_get_timebase(e->cosmology, e->ti_current);
+    io_write_attribute_d(h_grp, "TimeBase_dt", delta_t);
+  } else {
+    io_write_attribute_d(h_grp, "TimeBase_dloga", 0);
+    io_write_attribute_d(h_grp, "TimeBase_dt", e->time_base);
+  }
+
   /* Store the time at which the snapshot was written */
   time_t tm = time(NULL);
   struct tm* timeinfo = localtime(&tm);