Allow for more general expansion histories with evolving euqation of state for dark-energy.

examples/main.c

@@ -474,9 +474,26 @@ int main(int argc, char *argv[]) {
 
   /* Initialise the cosmology */
   struct cosmology cosmo;
-  if(with_cosmology) cosmology_init(params, &us, &prog_const, &cosmo);
-  if(with_cosmology) cosmology_print(&cosmo);
-  
+  if (with_cosmology) cosmology_init(params, &us, &prog_const, &cosmo);
+  if (with_cosmology) cosmology_print(&cosmo);
+
+  struct engine ee;
+  ee.ti_current = 0;
+  ee.timeBase = (log(cosmo.a_end) - log(cosmo.a_begin)) / max_nr_timesteps;
+  cosmology_update(&cosmo, &ee);
+
+  for (int i = 0; i <= 16; ++i) {
+
+    ee.ti_current = (max_nr_timesteps / 16) * i;
+
+    cosmology_update(&cosmo, &ee);
+
+    message("z=%e H(z)=%e w=%f t=%e [yrs]", cosmo.z, cosmo.H, cosmo.w,
+            cosmo.time / prog_const.const_year);
+  }
+
+  return 0;
+
   /* Initialise the hydro properties */
   struct hydro_props hydro_properties;
   if (with_hydro) hydro_props_init(&hydro_properties, params);
@@ -573,13 +590,13 @@ int main(int argc, char *argv[]) {
     fflush(stdout);
   }
 
-  /* Also update the total counts (in case of changes due to replication) */
+/* Also update the total counts (in case of changes due to replication) */
 #if defined(WITH_MPI)
   N_long[0] = s.nr_parts;
   N_long[1] = s.nr_gparts;
   N_long[2] = s.nr_sparts;
   MPI_Allreduce(&N_long, &N_total, 3, MPI_LONG_LONG_INT, MPI_SUM,
-             MPI_COMM_WORLD);
+                MPI_COMM_WORLD);
 #else
   N_total[0] = s.nr_parts;
   N_total[1] = s.nr_gparts;
@@ -614,7 +631,7 @@ int main(int argc, char *argv[]) {
     message("nr of cells at depth %i is %i.", data[0], data[1]);
   }
 
-  /* Initialise the table of Ewald corrections for the gravity checks */
+/* Initialise the table of Ewald corrections for the gravity checks */
 #ifdef SWIFT_GRAVITY_FORCE_CHECKS
   if (periodic) gravity_exact_force_ewald_init(dim[0]);
 #endif

src/cosmology.c

@@ -37,6 +37,8 @@
 
 /*! Number of values stored in the cosmological interpolation tables */
 const int cosmology_table_length = 10000;
+
+/*! Size of the GSL workspace */
 const size_t GSL_workspace_size = 10000;
 
 /**
@@ -63,14 +65,29 @@ static INLINE double interp_table(double *table, double x, double x_min,
     return table[ii - 1] + (table[ii] - table[ii - 1]) * (xx - ii);
 }
 
+/**
+ * @brief Computes the dark-energy equation of state at a given scale-factor a.
+ *
+ * We follow the convention of Linder & Jenkins, MNRAS, 346, 573, 2003
+ *
+ * @param a The current scale-factor
+ * @param w_0 The equation of state parameter at z=0
+ * @param w_a The equation of state evolution parameter
+ */
+static INLINE double cosmology_dark_energy_EoS(double a, double w_0,
+                                               double w_a) {
+
+  return w_0 + w_a * (1. - a);
+}
+
 /**
  * @brief Compute \f$ E(z) \f$.
  */
-static INLINE double E(double Or, double Om, double Ok, double Ol,
+static INLINE double E(double Or, double Om, double Ok, double Ol, double w,
                        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);
+              Ok * a_inv * a_inv + Ol * pow(a_inv, 3. * (1. + w)));
 }
 
 /**
@@ -106,12 +123,16 @@ void cosmology_update(struct cosmology *c, const struct engine *e) {
   /* Redshift */
   c->z = a_inv - 1.;
 
+  /* Dark-energy equation of state */
+  c->w = cosmology_dark_energy_EoS(a, c->w_0, c->w_a);
+
   /* 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);
+  const double w = c->w;
+  const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, w, a_inv);
 
   /* H(z) */
   c->H = c->H0 * E_z;
@@ -132,6 +153,8 @@ void cosmology_init_tables(struct cosmology *c) {
   const double Omega_m = c->Omega_m;
   const double Omega_k = c->Omega_k;
   const double Omega_l = c->Omega_lambda;
+  const double w_0 = c->w_0;
+  const double w_a = c->w_a;
   const double H0 = c->H0;
 
   /* Allocate memory for the interpolation tables */
@@ -146,7 +169,8 @@ void cosmology_init_tables(struct cosmology *c) {
   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 w = cosmology_dark_energy_EoS(a, w_0, w_a);
+    const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, w, a_inv);
     const double H = H0 * E_z;
 
     return (1. / H) * a_inv * a_inv * a_inv;
@@ -156,7 +180,8 @@ void cosmology_init_tables(struct cosmology *c) {
   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 w = cosmology_dark_energy_EoS(a, w_0, w_a);
+    const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, w, a_inv);
     const double H = H0 * E_z;
 
     return (1. / H) * a_inv * a_inv;
@@ -166,7 +191,8 @@ void cosmology_init_tables(struct cosmology *c) {
   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 w = cosmology_dark_energy_EoS(a, w_0, w_a);
+    const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, w, a_inv);
     const double H = H0 * E_z;
 
     return (1. / H) * a_inv;
@@ -240,9 +266,11 @@ void cosmology_init(const struct swift_params *params,
 
   /* Read in the cosmological parameters */
   c->Omega_m = parser_get_param_double(params, "Cosmology:Omega_m");
-  c->Omega_r = parser_get_param_double(params, "Cosmology:Omega_r");
+  c->Omega_r = parser_get_param_double(params, "Cosmology:Omega_r", 0.);
   c->Omega_lambda = parser_get_param_double(params, "Cosmology:Omega_lambda");
   c->Omega_b = parser_get_param_double(params, "Cosmology:Omega_b");
+  c->w_0 = parser_get_opt_param_double(params, "Cosmology:w_0", -1.);
+  c->w_a = parser_get_opt_param_double(params, "Cosmology:w_a", 0.);
   c->h = parser_get_param_double(params, "Cosmology:h");
 
   /* Read the start and end of the simulation */
@@ -256,6 +284,9 @@ void cosmology_init(const struct swift_params *params,
   /* Curvature density (for closure) */
   c->Omega_k = 1. - (c->Omega_m + c->Omega_r + c->Omega_lambda);
 
+  /* Dark-energy equation of state */
+  c->w = c->w_0 + c->w_a * (1. - c->a_begin);
+
   /* Hubble constant in internal units */
   const double km = 1.e5 / units_cgs_conversion_factor(us, UNIT_CONV_LENGTH);
   const double H0_cgs =
@@ -282,6 +313,7 @@ void cosmology_print(const struct cosmology *c) {
   message(
       "Density parameters: [O_m, O_l, O_b, O_k, O_r] = [%f, %f, %f, %f, %f]",
       c->Omega_m, c->Omega_lambda, c->Omega_b, c->Omega_k, c->Omega_r);
+  message("Dark energy equation of state: w_0=%f w_a=%f", c->w_0, c->w_a);
   message("Hubble constant: h = %f, H_0 = %e U_t^(-1) (internal units)", c->h,
           c->H0);
 }

src/cosmology.h

@@ -45,6 +45,9 @@ struct cosmology {
   /*! Time (in internal units) since the Big Bang */
   double time;
 
+  /*! Dark-energy equation of state at the current time */
+  double w;
+
   /*! Starting expansion factor */
   double a_begin;
 
@@ -72,6 +75,12 @@ struct cosmology {
   /*! Curvature density parameter */
   double Omega_k;
 
+  /*! Dark-energy equation of state at z=0 */
+  double w_0;
+
+  /*! Dark-energy evolution parameter */
+  double w_a;
+
   /*! Log of starting expansion factor */
   double log_a_begin;