diff --git a/src/cosmology.c b/src/cosmology.c
index bd80c9809c4f21cb183c01fded6d65271030f483..40e86b3f7948c1c84f23f0ccfb3571f9c75aae29 100644
--- a/src/cosmology.c
+++ b/src/cosmology.c
@@ -309,7 +309,7 @@ void cosmology_init_tables(struct cosmology *c) {
c->time_interp_table =
(double *)malloc(cosmology_table_length * sizeof(double));
c->scale_factor_interp_table =
- (double *)malloc(cosmology_table_length * sizeof(double));
+ (double *)malloc(cosmology_table_length * sizeof(double));
/* Prepare a table of scale factors for the integral bounds */
const double delta_a =
@@ -376,7 +376,8 @@ void cosmology_init_tables(struct cosmology *c) {
/* Inverse t(a) */
const double time_init = c->time_interp_table_offset;
- const double delta_t = (c->universe_age_at_present_day - time_init) / cosmology_table_length;
+ const double delta_t =
+ (c->universe_age_at_present_day - time_init) / cosmology_table_length;
int i_prev = 0;
for (int i = 0; i < cosmology_table_length; i++) {
@@ -384,7 +385,8 @@ void cosmology_init_tables(struct cosmology *c) {
double time_interp = delta_t * i;
/* Find next time in time_interp_table */
- while (i_prev < cosmology_table_length && c->time_interp_table[i_prev] <= time_interp) {
+ while (i_prev < cosmology_table_length &&
+ c->time_interp_table[i_prev] <= time_interp) {
i_prev++;
}
@@ -394,7 +396,9 @@ void cosmology_init_tables(struct cosmology *c) {
scale += i_prev;
/* Compute interpolated scale factor */
- double log_a = c->log_a_begin + scale * (c->log_a_end - c->log_a_begin) / cosmology_table_length;
+ double log_a =
+ c->log_a_begin +
+ scale * (c->log_a_end - c->log_a_begin) / cosmology_table_length;
c->scale_factor_interp_table[i] = exp(log_a) - c->a_begin;
}
@@ -665,16 +669,6 @@ double cosmology_get_delta_time(const struct cosmology *c,
return t2 - t1;
}
-
-/**
- * @brief Compute the scale factor from the redshift.
- *
- * @brief redshift The redshift to compute.
- */
-double cosmology_get_a_from_z(double redshift) {
- return 1. / (1. + redshift);
-}
-
/**
* @brief Compute scale factor from time since big bang (in internal units).
*
@@ -689,7 +683,7 @@ 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,
- c->universe_age_at_present_day);
+ c->universe_age_at_present_day);
return a + c->a_begin;
}
diff --git a/src/cosmology.h b/src/cosmology.h
index 0b7ddfaa75cc89520ea30fa0d6706a17a617c40d..7136b65667195953971060b76ddfd447a5fdf500 100644
--- a/src/cosmology.h
+++ b/src/cosmology.h
@@ -183,8 +183,6 @@ double cosmology_get_therm_kick_factor(const struct cosmology *cosmo,
double cosmology_get_delta_time(const struct cosmology *c,
integertime_t ti_start, integertime_t ti_end);
-double cosmology_get_a_from_z(double redshift);
-
double cosmology_get_scale_factor(const struct cosmology *cosmo, double t);
double cosmology_get_time_since_big_bang(const struct cosmology *c, double a);
@@ -196,7 +194,6 @@ void cosmology_init_no_cosmo(struct cosmology *c);
void cosmology_print(const struct cosmology *c);
void cosmology_clean(struct cosmology *c);
-
#ifdef HAVE_HDF5
void cosmology_write_model(hid_t h_grp, const struct cosmology *c);
#endif
diff --git a/src/engine.c b/src/engine.c
index 7c1648b25f5cbeb2a5327c51ab2aa8cb493b9302..ea3426f487e1d997e567651f10c27fb224d1e7c4 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -6569,6 +6569,12 @@ void engine_compute_next_statistics_time(struct engine *e) {
return;
}
+ /* Do outputlist file case */
+ if (e->outputlist_stats) {
+ engine_read_next_statistics_time(e);
+ return;
+ }
+
/* Find upper-bound on last output */
double time_end;
if (e->policy & engine_policy_cosmology)
@@ -6677,6 +6683,66 @@ void engine_compute_next_stf_time(struct engine *e) {
}
}
+/**
+ * @brief Read the next time (on the time line) for a dump
+ *
+ * @param e The #engine.
+ */
+void engine_read_next_snapshot_time(struct engine *e) {
+
+ int is_cosmo = e->policy & engine_policy_cosmology;
+ const struct outputlist *t = e->outputlist_snapshots;
+
+ /* Find upper-bound on last output */
+ double time_end;
+ if (is_cosmo)
+ time_end = e->cosmology->a_end;
+ else
+ time_end = e->time_end;
+
+ /* Find next snasphot above current time */
+ double time;
+ time = t->times[0];
+ size_t ind = 0;
+ while (time < time_end) {
+
+ /* Output time on the integer timeline */
+ if (is_cosmo)
+ e->ti_next_snapshot = log(time / e->cosmology->a_begin) / e->time_base;
+ else
+ e->ti_next_snapshot = (time - e->time_begin) / e->time_base;
+
+ /* Found it? */
+ if (e->ti_next_snapshot > e->ti_current) break;
+
+ ind += 1;
+ if (ind == t->size) break;
+
+ time = t->times[ind];
+ }
+
+ /* Deal with last snapshot */
+ if (e->ti_next_snapshot >= max_nr_timesteps || ind == t->size ||
+ time >= time_end) {
+ e->ti_next_snapshot = -1;
+ if (e->verbose) message("No further output time.");
+ } else {
+
+ /* Be nice, talk... */
+ if (is_cosmo) {
+ const double next_snapshot_time =
+ exp(e->ti_next_snapshot * e->time_base) * e->cosmology->a_begin;
+ if (e->verbose)
+ message("Next snapshot time set to a=%e.", next_snapshot_time);
+ } else {
+ const double next_snapshot_time =
+ e->ti_next_snapshot * e->time_base + e->time_begin;
+ if (e->verbose)
+ message("Next snapshot time set to t=%e.", next_snapshot_time);
+ }
+ }
+}
+
/**
* @brief Read the next time (on the time line) for a statistics dump
*
@@ -6708,16 +6774,19 @@ void engine_read_next_statistics_time(struct engine *e) {
if (e->ti_next_stats > e->ti_current) break;
ind += 1;
- if (ind == t->size)
- break;
-
+ if (ind == t->size) break;
+
time = t->times[ind];
}
/* Deal with last statistics */
- if (e->ti_next_stats >= max_nr_timesteps ||
- ind == t->size || time >= time_end) {
+ if (e->ti_next_stats >= max_nr_timesteps || ind == t->size ||
+ time >= time_end) {
e->ti_next_stats = -1;
+ if (e->verbose) message("No further output time.");
+ } else {
+
+ /* Be nice, talk... */
if (is_cosmo) {
const double next_statistics_time =
exp(e->ti_next_stats * e->time_base) * e->cosmology->a_begin;
@@ -6735,7 +6804,7 @@ void engine_read_next_statistics_time(struct engine *e) {
}
/**
- * @brief Initialize all the outputlist required by the engine
+ * @brief Initialize all the outputlist required by the engine
*
* @param e The #engine.
* @param params The #swift_params.
@@ -6746,20 +6815,20 @@ void engine_init_outputlists(struct engine *e, struct swift_params *params) {
/* get cosmo */
struct cosmology *cosmo = NULL;
- if (e->policy & engine_policy_cosmology)
- cosmo = e->cosmology;
-
+ if (e->policy & engine_policy_cosmology) cosmo = e->cosmology;
+
/* Deal with snapshots */
- int outputlist_on = parser_get_opt_param_int(params, "Snapshots:output_list_on", 0);
+ int outputlist_on =
+ parser_get_opt_param_int(params, "Snapshots:output_list_on", 0);
/* Read outputlist for snapshots */
if (outputlist_on) {
- e->outputlist_snapshots = (struct outputlist*) malloc(sizeof(struct outputlist));
+ e->outputlist_snapshots =
+ (struct outputlist *)malloc(sizeof(struct outputlist));
list = e->outputlist_snapshots;
- parser_get_param_string(params, "Snapshots:output_list",
- filename);
-
+ parser_get_param_string(params, "Snapshots:output_list", filename);
+
message("Reading snapshots output file.");
outputlist_read_file(list, filename, cosmo);
@@ -6770,25 +6839,24 @@ void engine_init_outputlists(struct engine *e, struct swift_params *params) {
if (cosmo) {
e->delta_time_snapshot = list->times[1] / list->times[0];
e->a_first_snapshot = list->times[0];
- }
- else {
+ } else {
e->delta_time_snapshot = list->times[1] - list->times[0];
e->time_first_snapshot = list->times[0];
}
}
- /* Deal with stats */
-
- outputlist_on = parser_get_opt_param_int(params, "Statistics:output_list_on", 0);
+ /* Deal with stats */
+ outputlist_on =
+ parser_get_opt_param_int(params, "Statistics:output_list_on", 0);
/* Read outputlist for stats */
if (outputlist_on) {
- e->outputlist_stats = (struct outputlist*) malloc(sizeof(struct outputlist));
+ e->outputlist_stats =
+ (struct outputlist *)malloc(sizeof(struct outputlist));
list = e->outputlist_stats;
- parser_get_param_string(params, "Statistics:output_list",
- filename);
-
+ parser_get_param_string(params, "Statistics:output_list", filename);
+
message("Reading statistics output file.");
outputlist_read_file(list, filename, cosmo);
@@ -6799,8 +6867,7 @@ void engine_init_outputlists(struct engine *e, struct swift_params *params) {
if (cosmo) {
e->delta_time_statistics = list->times[1] / list->times[0];
e->a_first_statistics = list->times[0];
- }
- else {
+ } else {
e->delta_time_statistics = list->times[1] - list->times[0];
e->time_first_statistics = list->times[0];
}
@@ -6992,8 +7059,7 @@ void engine_struct_dump(struct engine *e, FILE *stream) {
parser_struct_dump(e->parameter_file, stream);
if (e->outputlist_snapshots)
outputlist_struct_dump(e->outputlist_snapshots, stream);
- if (e->outputlist_stats)
- outputlist_struct_dump(e->outputlist_stats, stream);
+ if (e->outputlist_stats) outputlist_struct_dump(e->outputlist_stats, stream);
}
/**
diff --git a/src/outputlist.c b/src/outputlist.c
index 95a031ab8a4b9977eaa0f855456311a36f2095e4..37ebbcb5ea6200602b048493014187213e117b20 100644
--- a/src/outputlist.c
+++ b/src/outputlist.c
@@ -1,7 +1,6 @@
/*******************************************************************************
* This file is part of SWIFT.
- * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
- * Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ * Copyright (c) 2018 Loic Hausamman (loic.hausammann@epfl.ch)
*
* 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
@@ -39,11 +38,11 @@
* @param filename The file to read.
* @param cosmo The #cosmology model.
*/
-void outputlist_read_file(struct outputlist *outputlist, const char* filename,
- struct cosmology *cosmo) {
+void outputlist_read_file(struct outputlist *outputlist, const char *filename,
+ struct cosmology *cosmo) {
/* Open file */
- FILE* file = fopen(filename, "r");
+ FILE *file = fopen(filename, "r");
if (file == NULL) error("Error opening file '%s'", filename);
/* Declare reading variables */
@@ -60,13 +59,15 @@ void outputlist_read_file(struct outputlist *outputlist, const char* filename,
/* Return to start of file and initialize time array */
fseek(file, 0, SEEK_SET);
- outputlist->times = (double *) malloc(sizeof(double) * nber_line);
+ outputlist->times = (double *)malloc(sizeof(double) * nber_line);
if (!outputlist->times)
- error("Unable to malloc outputlist time array. "
- "Try reducing the number of lines in %s", filename);
+ error(
+ "Unable to malloc outputlist time array. "
+ "Try reducing the number of lines in %s",
+ filename);
/* Read header */
- if((read = getline(&line, &len, file)) == -1)
+ if ((read = getline(&line, &len, file)) == -1)
error("Unable to read header in file '%s'", filename);
/* Remove end of line character */
@@ -81,14 +82,14 @@ void outputlist_read_file(struct outputlist *outputlist, const char* filename,
else if (!strcmp(line, "# Scale Factor"))
type = OUTPUTLIST_SCALE_FACTOR;
else
- error("Unable to interpret the header (%s) in file '%s'",
- line, filename);
+ error("Unable to interpret the header (%s) in file '%s'", line, filename);
if (!cosmo &&
(type == OUTPUTLIST_SCALE_FACTOR || type == OUTPUTLIST_REDSHIFT))
- error("Unable to compute a redshift or a scale factor without cosmology. "
- "Please change the header in '%s'", filename);
-
+ error(
+ "Unable to compute a redshift or a scale factor without cosmology. "
+ "Please change the header in '%s'",
+ filename);
/* Read file */
size_t ind = 0;
@@ -97,14 +98,13 @@ void outputlist_read_file(struct outputlist *outputlist, const char* filename,
/* Write data to outputlist */
if (sscanf(line, "%lf", time) != 1) {
error(
- "Tried parsing double but found '%s' with illegal double "
- "characters in file '%s'.",
- line, filename);
+ "Tried parsing double but found '%s' with illegal double "
+ "characters in file '%s'.",
+ line, filename);
}
/* Transform input into correct time (e.g. ages or scale factor) */
- if (type == OUTPUTLIST_REDSHIFT)
- *time = cosmology_get_a_from_z(*time);
+ if (type == OUTPUTLIST_REDSHIFT) *time = 1. / (1. + *time);
if (cosmo && type == OUTPUTLIST_AGE)
*time = cosmology_get_scale_factor(cosmo, *time);
@@ -123,7 +123,7 @@ void outputlist_print(const struct outputlist *outputlist) {
printf("/*\t Time Array\t */\n");
printf("Number of Line: %lu\n", outputlist->size);
- for(size_t ind = 0; ind < outputlist->size; ind++) {
+ for (size_t ind = 0; ind < outputlist->size; ind++) {
printf("\t%lf\n", outputlist->times[ind]);
}
}
@@ -139,19 +139,21 @@ void outputlist_clean(struct outputlist *outputlist) {
* @brief Dump an #outputlist in a restart file
*/
void outputlist_struct_dump(struct outputlist *list, FILE *stream) {
- restart_write_blocks(list, sizeof(struct outputlist), 1, stream, "outputlist", "outputlist struct");
+ restart_write_blocks(list, sizeof(struct outputlist), 1, stream, "outputlist",
+ "outputlist struct");
- restart_write_blocks(list->times, list->size, sizeof(double), stream, "outputlist", "times");
+ restart_write_blocks(list->times, list->size, sizeof(double), stream,
+ "outputlist", "times");
}
-
/**
* @brief Restore an #outputlist from a restart file
*/
-void outputlist_struct_restore(struct outputlist * list, FILE *stream) {
- restart_read_blocks(list, sizeof(struct outputlist), 1, stream, NULL, "outputlist struct");
-
- list->times = (double *) malloc(sizeof(double) * list->size);
- restart_read_blocks(list->times, list->size, sizeof(double), stream, NULL, "times");
+void outputlist_struct_restore(struct outputlist *list, FILE *stream) {
+ restart_read_blocks(list, sizeof(struct outputlist), 1, stream, NULL,
+ "outputlist struct");
+ list->times = (double *)malloc(sizeof(double) * list->size);
+ restart_read_blocks(list->times, list->size, sizeof(double), stream, NULL,
+ "times");
}
diff --git a/src/outputlist.h b/src/outputlist.h
index 829e6e4a7e2499197c4066936d909dc4e20a8f0e..9067f584141472fee9bb340580716be60bae8365 100644
--- a/src/outputlist.h
+++ b/src/outputlist.h
@@ -1,7 +1,6 @@
/*******************************************************************************
* This file is part of SWIFT.
- * Copyright (c) 2012 Pedro Gonnet (pedro.gonnet@durham.ac.uk),
- * Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ * Copyright (c) 2018 Loic Hausamman (loic.hausammann@epfl.ch)
*
* 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
@@ -37,11 +36,11 @@ struct outputlist {
size_t size;
};
-
-void outputlist_read_file(struct outputlist *outputlist, const char* filename, struct cosmology *cosmo);
+void outputlist_read_file(struct outputlist *outputlist, const char *filename,
+ struct cosmology *cosmo);
void outputlist_print(const struct outputlist *outputlist);
void outputlist_clean(struct outputlist *outputlist);
void outputlist_struct_dump(struct outputlist *list, FILE *stream);
-void outputlist_struct_restore(struct outputlist * list, FILE *stream);
+void outputlist_struct_restore(struct outputlist *list, FILE *stream);
-#endif // SWIFT_OUTPUTLIST_H
+#endif /* SWIFT_OUTPUTLIST_H */
diff --git a/tests/testCosmology.c b/tests/testCosmology.c
index f24b920df7027071df9613b701bb8bef03354208..698351ad952e7d0b5f7d8e354c45a1a2dd53f968 100644
--- a/tests/testCosmology.c
+++ b/tests/testCosmology.c
@@ -1,6 +1,6 @@
/*******************************************************************************
* This file is part of SWIFT.
- * Copyright (C) 2015 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ * Copyright (c) 2018 Loic Hausamman (loic.hausammann@epfl.ch)
*
* 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
@@ -48,7 +48,7 @@ int main(int argc, char *argv[]) {
struct unit_system us;
units_init_cgs(&us);
- /* initialization of phys_const */
+ /* initialization of phys_const */
struct phys_const phys_const;
phys_const_init(&us, ¶ms, &phys_const);
@@ -58,7 +58,7 @@ int main(int argc, char *argv[]) {
message("Start checking computation...");
- for(int i=0; i < N_CHECK; i++) {
+ for (int i = 0; i < N_CHECK; i++) {
double a = 0.1 + 0.9 * i / (N_CHECK - 1.);
/* Compute a(t(a)) and check if same results */
double tmp = cosmology_get_time_since_big_bang(&cosmo, a);