cosmology.c 28 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
/*******************************************************************************
 * This file is part of SWIFT.
 * Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
 *
 * This program is free software: you can redistribute it and/or modify
 * it under the terms of the GNU Lesser General Public License as published
 * by the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * This program is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU Lesser General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 ******************************************************************************/

/**
 *  @file cosmology.c
 *  @brief Functions relating cosmological parameters
 */

/* This object's header. */
#include "cosmology.h"

/* Some standard headers */
#include <math.h>

31
/* Local headers */
32
#include "adiabatic_index.h"
33
#include "common_io.h"
34
#include "inline.h"
35
#include "restart.h"
36
37
38
39
40

#ifdef HAVE_LIBGSL
#include <gsl/gsl_integration.h>
#endif

41
/*! Number of values stored in the cosmological interpolation tables */
42
const int cosmology_table_length = 10000;
43

44
#ifdef HAVE_LIBGSL
45
/*! Size of the GSL workspace */
46
const size_t GSL_workspace_size = 100000;
47
#endif
48

49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
/**
 * @brief Returns the interpolated value from a table.
 *
 * Uses linear interpolation.
 *
 * @brief table The table of value to interpolate from (should be of length
 * cosmology_table_length).
 * @brief x The value to interpolate at.
 * @brief x_min The mininum of the range of x.
 * @brief x_max The maximum of the range of x.
 */
static INLINE double interp_table(double *table, double x, double x_min,
                                  double x_max) {

  const double xx = ((x - x_min) / (x_max - x_min)) * cosmology_table_length;
  const int i = (int)xx;
  const int ii = (i >= cosmology_table_length) ? cosmology_table_length - 1 : i;

  if (ii <= 1)
    return table[0] * xx;
  else
    return table[ii - 1] + (table[ii] - table[ii - 1]) * (xx - ii);
}

73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
/**
 * @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);
}

88
/**
89
 * @brief Computes the integral of the dark-energy equation of state
90
91
92
93
94
95
 * up to a scale-factor a.
 *
 * We follow the convention of Linder & Jenkins, MNRAS, 346, 573, 2003
 * and compute \f$ \tilde{w}(a) = \int_0^a\frac{1 + w(z)}{1+z}dz \f$.
 *
 * @param a The current scale-factor.
Matthieu Schaller's avatar
Matthieu Schaller committed
96
97
 * @param w0 The equation of state parameter at z=0
 * @param wa The equation of state evolution parameter
98
 */
99
100
101
static INLINE double w_tilde(double a, double w0, double wa) {
  return (a - 1.) * wa - (1. + w0 + wa) * log(a);
}
102

103
104
105
106
/**
 * @brief Compute \f$ E(z) \f$.
 */
static INLINE double E(double Or, double Om, double Ok, double Ol, double w0,
107
                       double wa, double a) {
108
  const double a_inv = 1. / a;
109
  return sqrt(Or * a_inv * a_inv * a_inv * a_inv + Om * a_inv * a_inv * a_inv +
110
              Ok * a_inv * a_inv + Ol * exp(3. * w_tilde(a, w0, wa)));
111
112
}

113
114
115
116
117
118
119
120
/**
 * @brief Returns the time (in internal units) since Big Bang at a given
 * scale-factor.
 *
 * @param c The current #cosmology.
 * @param a Scale-factor of interest.
 */
double cosmology_get_time_since_big_bang(const struct cosmology *c, double a) {
121

122
#ifdef SWIFT_DEBUG_CHECKS
123
  if (a < c->a_begin) error("Error a can't be smaller than a_begin");
124
125
#endif

126
127
128
  /* Time between a_begin and a */
  const double delta_t =
      interp_table(c->time_interp_table, log(a), c->log_a_begin, c->log_a_end);
129

130
  return c->time_interp_table_offset + delta_t;
131
132
}

133
134
135
136
/**
 * @brief Update the cosmological parameters to the current simulation time.
 *
 * @param c The #cosmology struct.
137
 * @param phys_const The physical constants in the internal units.
138
 * @param ti_current The current (integer) time.
139
 */
140
141
void cosmology_update(struct cosmology *c, const struct phys_const *phys_const,
                      integertime_t ti_current) {
142

143
  /* Get scale factor and powers of it */
144
  const double a = c->a_begin * exp(ti_current * c->time_base);
145
146
147
  const double a_inv = 1. / a;
  c->a = a;
  c->a_inv = a_inv;
148
  c->a2_inv = a_inv * a_inv;
149
  c->a3_inv = a_inv * a_inv * a_inv;
150
151
152
  c->a_factor_internal_energy =
      pow(a, -3. * hydro_gamma_minus_one);          /* a^{3*(1-gamma)} */
  c->a_factor_pressure = pow(a, -3. * hydro_gamma); /* a^{-3*gamma} */
153
  c->a_factor_sound_speed =
154
      pow(a, -1.5 * hydro_gamma_minus_one); /* a^{3*(1-gamma)/2} */
155
156
  c->a_factor_grav_accel = a_inv * a_inv;   /* 1 / a^2 */
  c->a_factor_hydro_accel =
157
158
159
      pow(a, -3. * hydro_gamma + 2.); /* 1 / a^(3*gamma - 2) */
  c->a_factor_mu =
      pow(a, 0.5 * (3. * hydro_gamma - 5.)); /* a^{(3*gamma - 5) / 2} */
160
161

  /* Redshift */
Matthieu Schaller's avatar
Matthieu Schaller committed
162
  c->z = a_inv - 1.;
163

164
165
166
  /* Dark-energy equation of state */
  c->w = cosmology_dark_energy_EoS(a, c->w_0, c->w_a);

167
168
169
  /* E(z) */
  const double Omega_r = c->Omega_r;
  const double Omega_m = c->Omega_m;
170
  const double Omega_k = c->Omega_k;
171
  const double Omega_l = c->Omega_lambda;
172
173
174
  const double w0 = c->w_0;
  const double wa = c->w_a;
  const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, w0, wa, a);
175
176
177
178

  /* H(z) */
  c->H = c->H0 * E_z;

179
180
181
  /* Expansion rate */
  c->a_dot = c->H * c->a;

182
183
184
185
  /* Critical density */
  c->critical_density =
      3. * c->H * c->H / (8. * M_PI * phys_const->const_newton_G);

186
187
188
  /* Time-step conversion factor */
  c->time_step_factor = c->H;

189
  /* Time */
190
  c->time = cosmology_get_time_since_big_bang(c, a);
191
  c->lookback_time = c->universe_age_at_present_day - c->time;
192
193
}

194
195
196
197
/**
 * @brief Computes \f$ dt / a^2 \f$ for the current cosmology.
 *
 * @param a The scale-factor of interest.
198
 * @param param The current #cosmology.
199
200
 */
double drift_integrand(double a, void *param) {
201

Matthieu Schaller's avatar
Matthieu Schaller committed
202
  const struct cosmology *c = (const struct cosmology *)param;
203
204
  const double Omega_r = c->Omega_r;
  const double Omega_m = c->Omega_m;
205
  const double Omega_k = c->Omega_k;
206
  const double Omega_l = c->Omega_lambda;
207
208
  const double w_0 = c->w_0;
  const double w_a = c->w_a;
209
  const double H0 = c->H0;
Matthieu Schaller's avatar
Matthieu Schaller committed
210

211
  const double a_inv = 1. / a;
212
  const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, w_0, w_a, a);
213
  const double H = H0 * E_z;
Matthieu Schaller's avatar
Matthieu Schaller committed
214

215
216
  return (1. / H) * a_inv * a_inv * a_inv;
}
217

218
219
220
221
/**
 * @brief Computes \f$ dt / a \f$ for the current cosmology.
 *
 * @param a The scale-factor of interest.
222
 * @param param The current #cosmology.
223
224
 */
double gravity_kick_integrand(double a, void *param) {
Matthieu Schaller's avatar
Matthieu Schaller committed
225
226

  const struct cosmology *c = (const struct cosmology *)param;
227
228
229
230
231
232
233
  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 w_0 = c->w_0;
  const double w_a = c->w_a;
  const double H0 = c->H0;
Matthieu Schaller's avatar
Matthieu Schaller committed
234

235
  const double a_inv = 1. / a;
236
  const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, w_0, w_a, a);
237
  const double H = H0 * E_z;
Matthieu Schaller's avatar
Matthieu Schaller committed
238

239
240
  return (1. / H) * a_inv * a_inv;
}
241

242
243
244
245
/**
 * @brief Computes \f$ dt / a^{3(\gamma - 1) + 1} \f$ for the current cosmology.
 *
 * @param a The scale-factor of interest.
246
 * @param param The current #cosmology.
247
248
 */
double hydro_kick_integrand(double a, void *param) {
Matthieu Schaller's avatar
Matthieu Schaller committed
249
250

  const struct cosmology *c = (const struct cosmology *)param;
251
252
253
254
255
256
257
  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 w_0 = c->w_0;
  const double w_a = c->w_a;
  const double H0 = c->H0;
Matthieu Schaller's avatar
Matthieu Schaller committed
258

259
  const double a_inv = 1. / a;
260
  const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, w_0, w_a, a);
261
  const double H = H0 * E_z;
Matthieu Schaller's avatar
Matthieu Schaller committed
262

263
264
265
266
267
  /* Note: we can't use the pre-defined pow_gamma_xxx() function as
     as we need double precision accuracy for the GSL routine. */
  return (1. / H) * pow(a_inv, 3. * hydro_gamma_minus_one) * a_inv;
}

Bert Vandenbroucke's avatar
Bert Vandenbroucke committed
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
/**
 * @brief Computes \f$a dt\f$ for the current cosmology.
 *
 * @param a The scale-factor of interest.
 * @param param The current #cosmology.
 */
double hydro_kick_corr_integrand(double a, void *param) {

  const struct cosmology *c = (const struct cosmology *)param;
  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 w_0 = c->w_0;
  const double w_a = c->w_a;
  const double H0 = c->H0;

  const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, w_0, w_a, a);
  const double H = H0 * E_z;

  return 1. / H;
}

291
292
293
294
/**
 * @brief Computes \f$ dt \f$ for the current cosmology.
 *
 * @param a The scale-factor of interest.
295
 * @param param The current #cosmology.
296
297
 */
double time_integrand(double a, void *param) {
298

Matthieu Schaller's avatar
Matthieu Schaller committed
299
  const struct cosmology *c = (const struct cosmology *)param;
300
301
302
303
304
305
306
  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 w_0 = c->w_0;
  const double w_a = c->w_a;
  const double H0 = c->H0;
Matthieu Schaller's avatar
Matthieu Schaller committed
307

308
  const double a_inv = 1. / a;
309
  const double E_z = E(Omega_r, Omega_m, Omega_k, Omega_l, w_0, w_a, a);
310
  const double H = H0 * E_z;
Matthieu Schaller's avatar
Matthieu Schaller committed
311

312
313
  return (1. / H) * a_inv;
}
314

315
316
317
/**
 * @brief Initialise the interpolation tables for the integrals.
 */
318
void cosmology_init_tables(struct cosmology *c) {
319

320
#ifdef HAVE_LIBGSL
321

322
323
  /* Retrieve some constants */
  const double a_begin = c->a_begin;
324

325
  /* Allocate memory for the interpolation tables */
326
327
  c->drift_fac_interp_table =
      (double *)malloc(cosmology_table_length * sizeof(double));
328
  c->grav_kick_fac_interp_table =
329
      (double *)malloc(cosmology_table_length * sizeof(double));
330
  c->hydro_kick_fac_interp_table =
331
      (double *)malloc(cosmology_table_length * sizeof(double));
Bert Vandenbroucke's avatar
Bert Vandenbroucke committed
332
333
  c->hydro_kick_corr_interp_table =
      (double *)malloc(cosmology_table_length * sizeof(double));
334
335
  c->time_interp_table =
      (double *)malloc(cosmology_table_length * sizeof(double));
336
  c->scale_factor_interp_table =
Matthieu Schaller's avatar
Matthieu Schaller committed
337
      (double *)malloc(cosmology_table_length * sizeof(double));
338
339
340
341

  /* 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;
342
  double *a_table = (double *)malloc(cosmology_table_length * sizeof(double));
343
344
345
346
347
348
349
350
351
352
  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 */
353
  gsl_function F = {&drift_integrand, c};
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
  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;
  }

372
373
374
375
376
377
378
379
380
381
  /* Integrate the kick factor \int_{a_begin}^{a_table[i]} dt/a^(3(g-1)+1) */
  F.function = &hydro_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->hydro_kick_fac_interp_table[i] = result;
  }

Bert Vandenbroucke's avatar
Bert Vandenbroucke committed
382
383
384
385
386
387
388
389
390
391
  /* Integrate the kick correction factor \int_{a_begin}^{a_table[i]} a dt */
  F.function = &hydro_kick_corr_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->hydro_kick_corr_interp_table[i] = result;
  }

392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
  /* 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;

407
408
409
410
411
  /* Integrate the time \int_{0}^{1} dt */
  gsl_integration_qag(&F, 0., 1, 0, 1.0e-13, GSL_workspace_size,
                      GSL_INTEG_GAUSS61, space, &result, &abserr);
  c->universe_age_at_present_day = result;

Loic Hausammann's avatar
Loic Hausammann committed
412
413
414
  /* Update the times */
  c->time_begin = cosmology_get_time_since_big_bang(c, c->a_begin);
  c->time_end = cosmology_get_time_since_big_bang(c, c->a_end);
415

Loic Hausammann's avatar
Loic Hausammann committed
416
417
418
419
  /*
   * Inverse t(a)
   */

Bert Vandenbroucke's avatar
Bert Vandenbroucke committed
420
  const double delta_t = (c->time_end - c->time_begin) / cosmology_table_length;
421

Loic Hausammann's avatar
Loic Hausammann committed
422
423
424
425
426
427
428
  /* index in the time_interp_table */
  int i_a = 0;

  for (int i_time = 0; i_time < cosmology_table_length; i_time++) {
    /* Current time
     * time_interp_table = \int_a_begin^a => no need of time_begin */
    double time_interp = delta_t * (i_time + 1);
429
430

    /* Find next time in time_interp_table */
Loic Hausammann's avatar
Loic Hausammann committed
431
432
433
    while (i_a < cosmology_table_length &&
           c->time_interp_table[i_a] <= time_interp) {
      i_a++;
434
435
436
    }

    /* Find linear interpolation scaling */
Loic Hausammann's avatar
Loic Hausammann committed
437
438
439
440
441
442
443
    double scale = 0;
    if (i_a != cosmology_table_length) {
      scale = time_interp - c->time_interp_table[i_a - 1];
      scale /= c->time_interp_table[i_a] - c->time_interp_table[i_a - 1];
    }

    scale += i_a;
444
445

    /* Compute interpolated scale factor */
Matthieu Schaller's avatar
Matthieu Schaller committed
446
447
    double log_a = c->log_a_begin + scale * (c->log_a_end - c->log_a_begin) /
                                        cosmology_table_length;
Loic Hausammann's avatar
Loic Hausammann committed
448
    c->scale_factor_interp_table[i_time] = exp(log_a) - c->a_begin;
449
450
  }

451
452
453
454
455
456
457
458
459
  /* 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
460
461
}

462
463
464
465
/**
 * @brief Initialises the #cosmology from the values read in the parameter file.
 *
 * @param params The parsed values.
466
 * @param us The current internal system of units.
467
468
469
 * @param phys_const The physical constants in the current system of units.
 * @param c The #cosmology to initialise.
 */
Matthieu Schaller's avatar
Matthieu Schaller committed
470
void cosmology_init(struct swift_params *params, const struct unit_system *us,
Matthieu Schaller's avatar
Matthieu Schaller committed
471
                    const struct phys_const *phys_const, struct cosmology *c) {
472
473
474

  /* Read in the cosmological parameters */
  c->Omega_m = parser_get_param_double(params, "Cosmology:Omega_m");
475
  c->Omega_r = parser_get_opt_param_double(params, "Cosmology:Omega_r", 0.);
476
477
  c->Omega_lambda = parser_get_param_double(params, "Cosmology:Omega_lambda");
  c->Omega_b = parser_get_param_double(params, "Cosmology:Omega_b");
478
479
  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.);
480
  c->h = parser_get_param_double(params, "Cosmology:h");
Matthieu Schaller's avatar
Matthieu Schaller committed
481
482

  /* Read the start and end of the simulation */
483
484
  c->a_begin = parser_get_param_double(params, "Cosmology:a_begin");
  c->a_end = parser_get_param_double(params, "Cosmology:a_end");
485
486
  c->log_a_begin = log(c->a_begin);
  c->log_a_end = log(c->a_end);
487
488
  c->time_base = (c->log_a_end - c->log_a_begin) / max_nr_timesteps;
  c->time_base_inv = 1. / c->time_base;
489

490
491
492
493
494
  /* If a_begin == a_end we hang */

  if (c->a_begin >= c->a_end)
    error("a_begin must be strictly before (and not equal to) a_end");

495
  /* Construct derived quantities */
Matthieu Schaller's avatar
Matthieu Schaller committed
496
497

  /* Curvature density (for closure) */
498
  c->Omega_k = 1. - (c->Omega_m + c->Omega_r + c->Omega_lambda);
Matthieu Schaller's avatar
Matthieu Schaller committed
499

500
  /* Dark-energy equation of state */
501
  c->w = cosmology_dark_energy_EoS(c->a_begin, c->w_0, c->w_a);
502

Matthieu Schaller's avatar
Matthieu Schaller committed
503
  /* Hubble constant in internal units */
504
505
506
  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 */
Matthieu Schaller's avatar
Matthieu Schaller committed
507
  c->H0 = H0_cgs * units_cgs_conversion_factor(us, UNIT_CONV_TIME);
508
  c->Hubble_time = 1. / c->H0;
509

510
511
512
  /* Initialise the interpolation tables */
  c->drift_fac_interp_table = NULL;
  c->grav_kick_fac_interp_table = NULL;
513
  c->hydro_kick_fac_interp_table = NULL;
514
515
516
  c->time_interp_table = NULL;
  c->time_interp_table_offset = 0.;
  cosmology_init_tables(c);
517

518
  /* Set remaining variables to alid values */
519
  cosmology_update(c, phys_const, 0);
520

521
522
523
  /* Update the times */
  c->time_begin = cosmology_get_time_since_big_bang(c, c->a_begin);
  c->time_end = cosmology_get_time_since_big_bang(c, c->a_end);
524
}
525

526
527
528
/**
 * @brief Initialise the #cosmology for non-cosmological time-integration
 *
529
 * Essentially sets all constants to 1 or 0.
530
531
532
 *
 * @param c The #cosmology to initialise.
 */
533
void cosmology_init_no_cosmo(struct cosmology *c) {
534
535
536
537
538
539
540
541

  c->Omega_m = 0.;
  c->Omega_r = 0.;
  c->Omega_k = 0.;
  c->Omega_lambda = 0.;
  c->Omega_b = 0.;
  c->w_0 = 0.;
  c->w_a = 0.;
542
  c->h = 1.;
543
  c->w = -1.;
544
545
546
547
548
549

  c->a_begin = 1.;
  c->a_end = 1.;
  c->log_a_begin = 0.;
  c->log_a_end = 0.;

550
  c->H = 0.;
551
  c->a = 1.;
552
  c->z = 0.;
553
  c->a_inv = 1.;
554
  c->a2_inv = 1.;
555
  c->a3_inv = 1.;
556
557
  c->a_factor_internal_energy = 1.;
  c->a_factor_pressure = 1.;
558
  c->a_factor_sound_speed = 1.;
559
  c->a_factor_mu = 1.;
560
561
  c->a_factor_hydro_accel = 1.;
  c->a_factor_grav_accel = 1.;
562

563
564
  c->critical_density = 0.;

565
566
  c->time_step_factor = 1.;

567
568
569
570
571
572
573
574
  c->a_dot = 0.;
  c->time = 0.;
  c->universe_age_at_present_day = 0.;

  /* Initialise the interpolation tables */
  c->drift_fac_interp_table = NULL;
  c->grav_kick_fac_interp_table = NULL;
  c->hydro_kick_fac_interp_table = NULL;
575
  c->hydro_kick_corr_interp_table = NULL;
576
577
  c->time_interp_table = NULL;
  c->time_interp_table_offset = 0.;
578
  c->scale_factor_interp_table = NULL;
579
580
581

  c->time_begin = 0.;
  c->time_end = 0.;
582
583
}

584
585
586
587
588
/**
 * @brief Computes the cosmology factor that enters the drift operator.
 *
 * Computes \f$ \int_{a_start}^{a_end} dt/a^2 \f$ using the interpolation table.
 *
589
 * @param c The current #cosmology.
590
591
592
 * @param ti_start the (integer) time of the start of the drift.
 * @param ti_end the (integer) time of the end of the drift.
 */
593
double cosmology_get_drift_factor(const struct cosmology *c,
Matthieu Schaller's avatar
Matthieu Schaller committed
594
595
                                  integertime_t ti_start,
                                  integertime_t ti_end) {
596

597
598
599
600
#ifdef SWIFT_DEBUG_CHECKS
  if (ti_end < ti_start) error("ti_end must be >= ti_start");
#endif

601
602
  const double a_start = c->log_a_begin + ti_start * c->time_base;
  const double a_end = c->log_a_begin + ti_end * c->time_base;
603

Matthieu Schaller's avatar
Matthieu Schaller committed
604
605
606
607
  const double int_start = interp_table(c->drift_fac_interp_table, a_start,
                                        c->log_a_begin, c->log_a_end);
  const double int_end = interp_table(c->drift_fac_interp_table, a_end,
                                      c->log_a_begin, c->log_a_end);
608
609
610
611
612
613
614
615
616

  return int_end - int_start;
}

/**
 * @brief Computes the cosmology factor that enters the gravity kick operator.
 *
 * Computes \f$ \int_{a_start}^{a_end} dt/a \f$ using the interpolation table.
 *
617
 * @param c The current #cosmology.
618
619
620
 * @param ti_start the (integer) time of the start of the drift.
 * @param ti_end the (integer) time of the end of the drift.
 */
621
double cosmology_get_grav_kick_factor(const struct cosmology *c,
Matthieu Schaller's avatar
Matthieu Schaller committed
622
623
                                      integertime_t ti_start,
                                      integertime_t ti_end) {
624

625
626
627
628
#ifdef SWIFT_DEBUG_CHECKS
  if (ti_end < ti_start) error("ti_end must be >= ti_start");
#endif

629
630
  const double a_start = c->log_a_begin + ti_start * c->time_base;
  const double a_end = c->log_a_begin + ti_end * c->time_base;
631

Matthieu Schaller's avatar
Matthieu Schaller committed
632
633
634
635
  const double int_start = interp_table(c->grav_kick_fac_interp_table, a_start,
                                        c->log_a_begin, c->log_a_end);
  const double int_end = interp_table(c->grav_kick_fac_interp_table, a_end,
                                      c->log_a_begin, c->log_a_end);
636
637
638
639
640
641
642

  return int_end - int_start;
}

/**
 * @brief Computes the cosmology factor that enters the hydro kick operator.
 *
643
644
 * Computes \f$ \int_{a_start}^{a_end} dt/a^{3(gamma - 1)} \f$ using the
 * interpolation table.
645
 *
646
 * @param c The current #cosmology.
647
648
649
 * @param ti_start the (integer) time of the start of the drift.
 * @param ti_end the (integer) time of the end of the drift.
 */
650
double cosmology_get_hydro_kick_factor(const struct cosmology *c,
Matthieu Schaller's avatar
Matthieu Schaller committed
651
652
                                       integertime_t ti_start,
                                       integertime_t ti_end) {
653

654
655
656
657
#ifdef SWIFT_DEBUG_CHECKS
  if (ti_end < ti_start) error("ti_end must be >= ti_start");
#endif

658
659
  const double a_start = c->log_a_begin + ti_start * c->time_base;
  const double a_end = c->log_a_begin + ti_end * c->time_base;
660

661
  const double int_start = interp_table(c->hydro_kick_fac_interp_table, a_start,
662
                                        c->log_a_begin, c->log_a_end);
663
  const double int_end = interp_table(c->hydro_kick_fac_interp_table, a_end,
664
665
666
667
668
                                      c->log_a_begin, c->log_a_end);

  return int_end - int_start;
}

Bert Vandenbroucke's avatar
Bert Vandenbroucke committed
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
/**
 * @brief Computes the cosmology factor that enters the hydro kick correction
 * operator for the meshless schemes (GIZMO-MFV).
 *
 * Computes \f$ \int_{a_start}^{a_end} a dt \f$ using the interpolation table.
 *
 * @param c The current #cosmology.
 * @param ti_start the (integer) time of the start of the drift.
 * @param ti_end the (integer) time of the end of the drift.
 */
double cosmology_get_corr_kick_factor(const struct cosmology *c,
                                      integertime_t ti_start,
                                      integertime_t ti_end) {

#ifdef SWIFT_DEBUG_CHECKS
  if (ti_end < ti_start) error("ti_end must be >= ti_start");
#endif

  const double a_start = c->log_a_begin + ti_start * c->time_base;
  const double a_end = c->log_a_begin + ti_end * c->time_base;

  const double int_start = interp_table(c->hydro_kick_corr_interp_table,
                                        a_start, c->log_a_begin, c->log_a_end);
  const double int_end = interp_table(c->hydro_kick_corr_interp_table, a_end,
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
                                      c->log_a_begin, c->log_a_end);

  return int_end - int_start;
}

/**
 * @brief Computes the cosmology factor that enters the thermal variable kick
 * operator.
 *
 * Computes \f$ \int_{a_start}^{a_end} dt/a^2 \f$ using the interpolation table.
 *
 * @param c The current #cosmology.
 * @param ti_start the (integer) time of the start of the drift.
 * @param ti_end the (integer) time of the end of the drift.
 */
708
double cosmology_get_therm_kick_factor(const struct cosmology *c,
709
710
711
                                       integertime_t ti_start,
                                       integertime_t ti_end) {

712
713
714
715
#ifdef SWIFT_DEBUG_CHECKS
  if (ti_end < ti_start) error("ti_end must be >= ti_start");
#endif

716
717
718
  const double a_start = c->log_a_begin + ti_start * c->time_base;
  const double a_end = c->log_a_begin + ti_end * c->time_base;

719
  const double int_start = interp_table(c->drift_fac_interp_table, a_start,
Matthieu Schaller's avatar
Matthieu Schaller committed
720
                                        c->log_a_begin, c->log_a_end);
721
  const double int_end = interp_table(c->drift_fac_interp_table, a_end,
Matthieu Schaller's avatar
Matthieu Schaller committed
722
                                      c->log_a_begin, c->log_a_end);
723
724
725

  return int_end - int_start;
}
726
727

/**
728
729
 * @brief Compute the cosmic time (in internal units) between two points
 * on the integer time line.
730
 *
731
 * @param c The current #cosmology.
732
733
 * @param ti_start the (integer) time of the start.
 * @param ti_end the (integer) time of the end.
734
 */
735
736
double cosmology_get_delta_time(const struct cosmology *c,
                                integertime_t ti_start, integertime_t ti_end) {
737

738
#ifdef SWIFT_DEBUG_CHECKS
739
  if (ti_end < ti_start) error("ti_end must be >= ti_start");
740
741
#endif

742
743
744
745
746
747
  const double log_a_start = c->log_a_begin + ti_start * c->time_base;
  const double log_a_end = c->log_a_begin + ti_end * c->time_base;

  /* Time between a_begin and a_start */
  const double t1 = interp_table(c->time_interp_table, log_a_start,
                                 c->log_a_begin, c->log_a_end);
748

749
750
751
  /* Time between a_begin and a_end */
  const double t2 = interp_table(c->time_interp_table, log_a_end,
                                 c->log_a_begin, c->log_a_end);
752
753
754

  return t2 - t1;
}
755

756
/**
Loic Hausammann's avatar
Loic Hausammann committed
757
 * @brief Compute scale factor from time since big bang (in internal units).
758
 *
759
760
 * @param c The current #cosmology.
 * @param t time since the big bang
761
762
 * @return The scale factor.
 */
763
764
765
766
double cosmology_get_scale_factor(const struct cosmology *c, double t) {
  /* scale factor between time_begin and t */
  const double a =
      interp_table(c->scale_factor_interp_table, t, c->time_interp_table_offset,
Matthieu Schaller's avatar
Matthieu Schaller committed
767
                   c->universe_age_at_present_day);
768
  return a + c->a_begin;
769
770
}

771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
/**
 * @brief Prints the #cosmology model to stdout.
 */
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)", c->h, c->H0);
  message("Hubble time: 1/H0 = %e U_t", c->Hubble_time);
  message("Universe age at present day: %e U_t",
          c->universe_age_at_present_day);
}

786
787
788
789
790
void cosmology_clean(struct cosmology *c) {

  free(c->drift_fac_interp_table);
  free(c->grav_kick_fac_interp_table);
  free(c->hydro_kick_fac_interp_table);
Bert Vandenbroucke's avatar
Bert Vandenbroucke committed
791
  free(c->hydro_kick_corr_interp_table);
792
  free(c->time_interp_table);
793
  free(c->scale_factor_interp_table);
794
795
}

796
797
798
799
800
801
802
#ifdef HAVE_HDF5
void cosmology_write_model(hid_t h_grp, const struct cosmology *c) {

  io_write_attribute_d(h_grp, "a_beg", c->a_begin);
  io_write_attribute_d(h_grp, "a_end", c->a_end);
  io_write_attribute_d(h_grp, "time_beg [internal units]", c->time_begin);
  io_write_attribute_d(h_grp, "time_end [internal units]", c->time_end);
803
804
805
  io_write_attribute_d(h_grp, "Universe age [internal units]", c->time);
  io_write_attribute_d(h_grp, "Lookback time [internal units]",
                       c->lookback_time);
806
807
  io_write_attribute_d(h_grp, "h", c->h);
  io_write_attribute_d(h_grp, "H0 [internal units]", c->H0);
808
  io_write_attribute_d(h_grp, "H [internal units]", c->H);
809
810
811
812
813
814
815
816
  io_write_attribute_d(h_grp, "Hubble time [internal units]", c->Hubble_time);
  io_write_attribute_d(h_grp, "Omega_m", c->Omega_m);
  io_write_attribute_d(h_grp, "Omega_r", c->Omega_r);
  io_write_attribute_d(h_grp, "Omega_b", c->Omega_b);
  io_write_attribute_d(h_grp, "Omega_k", c->Omega_k);
  io_write_attribute_d(h_grp, "Omega_lambda", c->Omega_lambda);
  io_write_attribute_d(h_grp, "w_0", c->w_0);
  io_write_attribute_d(h_grp, "w_a", c->w_a);
817
818
819
820
821
  io_write_attribute_d(h_grp, "w", c->w);
  io_write_attribute_d(h_grp, "Redshift", c->z);
  io_write_attribute_d(h_grp, "Scale-factor", c->a);
  io_write_attribute_d(h_grp, "Critical density [internal units]",
                       c->critical_density);
822
823
824
}
#endif

825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
/**
 * @brief Write a cosmology struct to the given FILE as a stream of bytes.
 *
 * @param cosmology the struct
 * @param stream the file stream
 */
void cosmology_struct_dump(const struct cosmology *cosmology, FILE *stream) {
  restart_write_blocks((void *)cosmology, sizeof(struct cosmology), 1, stream,
                       "cosmology", "cosmology function");
}

/**
 * @brief Restore a cosmology struct from the given FILE as a stream of
 * bytes.
 *
840
 * @param enabled whether cosmology is enabled.
841
842
843
 * @param cosmology the struct
 * @param stream the file stream
 */
844
845
void cosmology_struct_restore(int enabled, struct cosmology *cosmology,
                              FILE *stream) {
846
847
  restart_read_blocks((void *)cosmology, sizeof(struct cosmology), 1, stream,
                      NULL, "cosmology function");
848

849
850
  /* Re-initialise the tables if using a cosmology. */
  if (enabled) cosmology_init_tables(cosmology);
851
}