From 279470f549d74611a42556632f39cd69738f2f8e Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <schaller@strw.leidenuniv.nl>
Date: Tue, 22 Sep 2020 09:10:51 +0100
Subject: [PATCH] Time base in snap header
---
doc/RTD/source/Snapshots/index.rst | 9 ++++++
src/cosmology.c | 48 ++++++++++++++++++++++++++++++
src/cosmology.h | 3 ++
src/distributed_io.c | 10 +++++++
src/parallel_io.c | 10 +++++++
src/serial_io.c | 11 +++++++
src/single_io.c | 10 +++++++
7 files changed, 101 insertions(+)
diff --git a/doc/RTD/source/Snapshots/index.rst b/doc/RTD/source/Snapshots/index.rst
index da94cfaf8a..f3e81c6410 100644
--- a/doc/RTD/source/Snapshots/index.rst
+++ b/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
--------------------------------
diff --git a/src/cosmology.c b/src/cosmology.c
index 4c72bb4ed1..0946ff900b 100644
--- a/src/cosmology.c
+++ b/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).
*
diff --git a/src/cosmology.h b/src/cosmology.h
index 20036f4335..0fc602c71e 100644
--- a/src/cosmology.h
+++ b/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);
diff --git a/src/distributed_io.c b/src/distributed_io.c
index df1f22f9e4..33073426c3 100644
--- a/src/distributed_io.c
+++ b/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);
diff --git a/src/parallel_io.c b/src/parallel_io.c
index b34a06fb0f..5bac0214bb 100644
--- a/src/parallel_io.c
+++ b/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);
diff --git a/src/serial_io.c b/src/serial_io.c
index 5381e4c4fc..4764c6cf15 100644
--- a/src/serial_io.c
+++ b/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);
diff --git a/src/single_io.c b/src/single_io.c
index d978e4c472..d4c4f3ed17 100644
--- a/src/single_io.c
+++ b/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);
--
GitLab