diff --git a/src/cosmology.c b/src/cosmology.c
index f33a41e1b3370ad8f6c84bbacfdcaa87c31558aa..3302bdcb519c2bfcb34d03b97508edacb09928dd 100644
--- a/src/cosmology.c
+++ b/src/cosmology.c
@@ -28,6 +28,44 @@
/* Some standard headers */
#include <math.h>
+/* Local headers */
+#include "inline.h"
+
+#ifdef HAVE_LIBGSL
+#include <gsl/gsl_integration.h>
+#endif
+
+const int cosmology_table_length = 10000;
+const size_t GSL_workspace_size = 10000;
+
+/**
+ * @brief Compute \f$ E(z) \f$.
+ */
+static INLINE double E(double Or, double Om, double Ok, double Ol,
+ double a_inv) {
+
+ return sqrt(Or * a_inv * a_inv * a_inv * a_inv + Om * a_inv * a_inv * a_inv +
+ Ok * a_inv * a_inv + Ol);
+}
+
+double cosmology_get_time(const struct cosmology *c, double a) {
+
+ const double log_a = log(a);
+
+ /* Position in the table */
+ const double x =
+ ((log_a - c->log_a_begin) / (c->log_a_end - c->log_a_begin)) *
+ cosmology_table_length;
+ const int i =
+ ((int)x >= cosmology_table_length) ? cosmology_table_length - 1 : (int)x;
+
+ if (i <= 1)
+ return c->time_interp_table_offset + x * c->time_interp_table[0];
+ else
+ return c->time_interp_table_offset + c->time_interp_table[i - 1] +
+ (c->time_interp_table[i] - c->time_interp_table[i - 1]) * (x - i);
+}
+
/**
* @brief Update the cosmological parameters to the current simulation time.
*
@@ -48,17 +86,129 @@ void cosmology_update(struct cosmology *c, const struct engine *e) {
/* E(z) */
const double Omega_r = c->Omega_r;
const double Omega_m = c->Omega_m;
+ const double Omega_k = c->Omega_k;
const double Omega_l = c->Omega_lambda;
+ const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, a_inv);
+
+ /* H(z) */
+ c->H = c->H0 * E_z;
+
+ /* Time since Big Bang */
+ c->time = cosmology_get_time(c, a);
+}
+
+void cosmology_init_tables(struct cosmology *c) {
+
+ const ticks tic = getticks();
+
+#ifdef HAVE_LIBGSL
+
+ /* Retrieve some constants */
+ const double a_begin = c->a_begin;
+ const double Omega_r = c->Omega_r;
+ const double Omega_m = c->Omega_m;
const double Omega_k = c->Omega_k;
+ const double Omega_l = c->Omega_lambda;
+ const double H0 = c->H0;
- const double E_z2 = Omega_r * a_inv * a_inv * a_inv * a_inv +
- Omega_m * a_inv * a_inv * a_inv +
- Omega_k * a_inv * a_inv + Omega_l;
+ /* Allocate memory for the interpolation tables */
+ c->drift_fac_interp_table = malloc(cosmology_table_length * sizeof(double));
+ c->grav_kick_fac_interp_table =
+ malloc(cosmology_table_length * sizeof(double));
+ c->time_interp_table = malloc(cosmology_table_length * sizeof(double));
- const double E_z = sqrt(E_z2);
+ /* Define the integrands we need */
- /* H(z) */
- c->H = c->H0 * E_z;
+ /* dt / a^2 */
+ double drift_integrand(double a, void *param) {
+
+ const double a_inv = 1. / a;
+ const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, a_inv);
+ const double H = H0 * E_z;
+
+ return (1. / H) * a_inv * a_inv * a_inv;
+ }
+
+ /* dt / a */
+ double gravity_kick_integrand(double a, void *param) {
+
+ const double a_inv = 1. / a;
+ const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, a_inv);
+ const double H = H0 * E_z;
+
+ return (1. / H) * a_inv * a_inv;
+ }
+
+ /* dt */
+ double time_integrand(double a, void *param) {
+
+ const double a_inv = 1. / a;
+ const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, a_inv);
+ const double H = H0 * E_z;
+
+ return (1. / H) * a_inv;
+ }
+
+ /* Prepare a table of scale factors for the integral bounds */
+ const double delta_a =
+ (c->log_a_end - c->log_a_begin) / cosmology_table_length;
+ double *a_table = malloc(cosmology_table_length * sizeof(double));
+ for (int i = 0; i < cosmology_table_length; i++)
+ a_table[i] = exp(c->log_a_begin + delta_a * (i + 1));
+
+ /* Initalise the GSL workspace */
+ gsl_integration_workspace *space =
+ gsl_integration_workspace_alloc(GSL_workspace_size);
+
+ double result, abserr;
+
+ /* Integrate the drift factor \int_{a_begin}^{a_table[i]} dt/a^2 */
+ gsl_function F = {&drift_integrand, NULL};
+ for (int i = 0; i < cosmology_table_length; i++) {
+ gsl_integration_qag(&F, a_begin, a_table[i], 0, 1.0e-10, GSL_workspace_size,
+ GSL_INTEG_GAUSS61, space, &result, &abserr);
+
+ /* Store result */
+ c->drift_fac_interp_table[i] = result;
+ }
+
+ /* Integrate the kick factor \int_{a_begin}^{a_table[i]} dt/a */
+ F.function = &gravity_kick_integrand;
+ for (int i = 0; i < cosmology_table_length; i++) {
+ gsl_integration_qag(&F, a_begin, a_table[i], 0, 1.0e-10, GSL_workspace_size,
+ GSL_INTEG_GAUSS61, space, &result, &abserr);
+
+ /* Store result */
+ c->grav_kick_fac_interp_table[i] = result;
+ }
+
+ /* Integrate the time \int_{a_begin}^{a_table[i]} dt */
+ F.function = &time_integrand;
+ for (int i = 0; i < cosmology_table_length; i++) {
+ gsl_integration_qag(&F, a_begin, a_table[i], 0, 1.0e-10, GSL_workspace_size,
+ GSL_INTEG_GAUSS61, space, &result, &abserr);
+
+ /* Store result */
+ c->time_interp_table[i] = result;
+ }
+
+ /* Integrate the time \int_{0}^{a_begin} dt */
+ gsl_integration_qag(&F, 0., a_begin, 0, 1.0e-10, GSL_workspace_size,
+ GSL_INTEG_GAUSS61, space, &result, &abserr);
+ c->time_interp_table_offset = result;
+
+ /* Free the workspace and temp array */
+ gsl_integration_workspace_free(space);
+ free(a_table);
+
+#else
+
+ error("Code not compiled with GSL. Can't compute cosmology integrals.");
+
+#endif
+
+ message("took %.3f %s.", clocks_from_ticks(getticks() - tic),
+ clocks_getunit());
}
void cosmology_init(const struct swift_params *params,
@@ -75,6 +225,8 @@ void cosmology_init(const struct swift_params *params,
/* Read the start and end of the simulation */
c->a_begin = parser_get_param_double(params, "Cosmology:a_begin");
c->a_end = parser_get_param_double(params, "Cosmology:a_end");
+ c->log_a_begin = log(c->a_begin);
+ c->log_a_end = log(c->a_end);
/* Construct derived quantities */
@@ -82,8 +234,9 @@ void cosmology_init(const struct swift_params *params,
c->Omega_k = 1. - (c->Omega_m + c->Omega_r + c->Omega_lambda);
/* Hubble constant in internal units */
- const double km = 1.e5 / units_cgs_conversion_factor(us, UNIT_CONV_LENGTH);
- const double H0_cgs = 100. * c->h * (km / (1.e6 * phys_const->const_parsec)); /* s^-1 */
+ const double km = 1.e5 / units_cgs_conversion_factor(us, UNIT_CONV_LENGTH);
+ const double H0_cgs =
+ 100. * c->h * (km / (1.e6 * phys_const->const_parsec)); /* s^-1 */
c->H0 = H0_cgs * units_cgs_conversion_factor(us, UNIT_CONV_TIME);
/* Set remaining variables to invalid values */
@@ -91,6 +244,14 @@ void cosmology_init(const struct swift_params *params,
c->a = -1.;
c->a_inv = -1;
c->z = -1.;
+ c->time = -1.;
+
+ /* Initialise the interpolation tables */
+ c->drift_fac_interp_table = NULL;
+ c->grav_kick_fac_interp_table = NULL;
+ c->time_interp_table = NULL;
+ c->time_interp_table_offset = 0.;
+ cosmology_init_tables(c);
}
void cosmology_print(const struct cosmology *c) {
diff --git a/src/cosmology.h b/src/cosmology.h
index df8cea4ec91687d2ca4f49ed6a71db89a6eb6672..6029d2be42b3d8c27a613a20008d2a556d7f4b7e 100644
--- a/src/cosmology.h
+++ b/src/cosmology.h
@@ -42,10 +42,13 @@ struct cosmology {
/*! Hubble constant at the current redshift (in internal units) */
double H;
+ /*! Time (in internal units) since the Big Bang */
+ double time;
+
/*! Starting expansion factor */
double a_begin;
- /*! Ending expansion factor */
+ /*! Final expansion factor */
double a_end;
/*! Reduced Hubble constant (H0 / (100km/s/Mpc)) */
@@ -68,6 +71,24 @@ struct cosmology {
/*! Curvature density parameter */
double Omega_k;
+
+ /*! Log of starting expansion factor */
+ double log_a_begin;
+
+ /*! Log of final expansion factor */
+ double log_a_end;
+
+ /*! Drift factor interpolation table */
+ double *drift_fac_interp_table;
+
+ /*! Kick factor interpolation table */
+ double *grav_kick_fac_interp_table;
+
+ /*! Time interpolation table */
+ double *time_interp_table;
+
+ /*! Time between Big Bang and first entry in the table */
+ double time_interp_table_offset;
};
void cosmology_update(struct cosmology *c, const struct engine *e);