diff --git a/doc/RTD/source/Snapshots/index.rst b/doc/RTD/source/Snapshots/index.rst index da94cfaf8ae377160916854cf1264db66aeb00ed..f3e81c6410da14f859eb8843e91f7f45889ddfd1 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 4c72bb4ed13bdf412d531f16f1be8c4967810a03..0946ff900baa4c8b3685c7f02c2351b4afd18818 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 20036f4335ef25d029a60143ac7e4376e45032a9..0fc602c71e7dfb9819d8f8189de80b12e0a63821 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 df1f22f9e4fe1436e134bae735d126a6df92fc4f..33073426c382f732b1cf0d035c616a7db96243ad 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 b34a06fb0fa648f7cbc951b8624feb1c756576e1..5bac0214bb4791db035a2ed8ea3744ef32122c0a 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 5381e4c4fcc9805b7cd57cc9c72877ad6c99eed8..4764c6cf1528cdc294c04123175283b07fb79930 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 d978e4c4726cc6a1d7bdd102cdfa70d52f675c38..d4c4f3ed1742274745d6b5dc2f74d7809063fa9b 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);