Create interpolation tables for drift factor, kick factor and time since Big Bang.

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) {

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);