diff --git a/doc/RTD/source/GettingStarted/parameter_file.rst b/doc/RTD/source/GettingStarted/parameter_file.rst
index cecf2fa085b23946b8bcbebd324372b01f4940c2..1bf4cbd3940a104c126aa256759edc24ca996338 100644
--- a/doc/RTD/source/GettingStarted/parameter_file.rst
+++ b/doc/RTD/source/GettingStarted/parameter_file.rst
@@ -12,6 +12,35 @@ parameters. Each section in this file corresponds to a different option in
SWIFT and are not always required depending on the configuration options and
the run time parameters.
+Output List
+~~~~~~~~~~~
+
+In the sections ``Snapshots`` and ``Statistics``, you can specify the options ``output_list_on`` and ``output_list`` which receive an int and a filename.
+The ``output_list_on`` enable or not the output list and ``output_list`` is the filename containing the output times.
+With the file header, you can choose between writing redshifts, scale factors or times.
+
+Example of file containing with times (in internal units):
+::
+ # Time
+ 0.5
+ 1.5
+ 3.0
+ 12.5
+
+Example of file with scale factors:
+::
+ # Scale Factor
+ 0.1
+ 0.2
+ 0.3
+
+Example of file with redshift:
+::
+ # Redshift
+ 20
+ 15
+ 10
+ 5
Output Selection
~~~~~~~~~~~~~~~~
diff --git a/examples/parameter_example.yml b/examples/parameter_example.yml
index 783b6030be74802d455034420d39ef9c741896cb..ddb71c594122a3e8d6ddbd7c5b73e0474b404a75 100644
--- a/examples/parameter_example.yml
+++ b/examples/parameter_example.yml
@@ -81,6 +81,8 @@ Snapshots:
UnitVelocity_in_cgs: 1 # (Optional) Unit system for the outputs (Centimeters per second)
UnitCurrent_in_cgs: 1 # (Optional) Unit system for the outputs (Amperes)
UnitTemp_in_cgs: 1 # (Optional) Unit system for the outputs (Kelvin)
+ output_list_on: 0 # (Optional) Enable the output list
+ output_list: snaplist.txt # (Optional) File containing the output times (see documentation in "Parameter File" section)
# Parameters governing the conserved quantities statistics
Statistics:
@@ -89,6 +91,8 @@ Statistics:
time_first: 0. # (Optional) Time of the first stats output if non-cosmological time-integration (in internal units)
energy_file_name: energy # (Optional) File name for energy output
timestep_file_name: timesteps # (Optional) File name for timing information output. Note: No underscores "_" allowed in file name
+ output_list_on: 0 # (Optional) Enable the output list
+ output_list: statlist.txt # (Optional) File containing the output times (see documentation in "Parameter File" section)
# Parameters related to the initial conditions
InitialConditions:
@@ -220,3 +224,14 @@ EAGLEChemistry:
CalciumOverSilicon: 0.0941736 # Constant ratio of Calcium over Silicon abundance
SulphurOverSilicon: 0.6054160 # Constant ratio of Sulphur over Silicon abundance
+# Structure finding options (requires velociraptor)
+StructureFinding:
+ config_file_name: stf_input.cfg # Name of the STF config file.
+ basename: ./stf # Common part of the name of output files.
+ output_time_format: 0 # Specifies the frequency format of structure finding. 0 for simulation steps (delta_step) and 1 for simulation time intervals (delta_time).
+ scale_factor_first: 0.92 # Scale-factor of the first snaphot (cosmological run)
+ time_first: 0.01 # Time of the first structure finding output (in internal units).
+ delta_step: 1000 # Time difference between consecutive structure finding outputs (in internal units) in simulation steps.
+ delta_time: 1.10 # Time difference between consecutive structure finding outputs (in internal units) in simulation time intervals.
+ output_list_on: 0 # (Optional) Enable the output list
+ output_list: stflist.txt # (Optional) File containing the output times (see documentation in "Parameter File" section)
diff --git a/src/Makefile.am b/src/Makefile.am
index 7e54ecab60a5b97a585bd560a04836e4e6f78336..fbc1d9fdc025dca2a498208ddce98e716fa2fd03 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -48,7 +48,7 @@ include_HEADERS = space.h runner.h queue.h task.h lock.h cell.h part.h const.h \
dump.h logger.h active.h timeline.h xmf.h gravity_properties.h gravity_derivatives.h \
gravity_softened_derivatives.h vector_power.h collectgroup.h hydro_space.h sort_part.h \
chemistry.h chemistry_io.h chemistry_struct.h cosmology.h restart.h space_getsid.h utilities.h \
- mesh_gravity.h cbrt.h velociraptor_interface.h swift_velociraptor_part.h
+ mesh_gravity.h cbrt.h velociraptor_interface.h swift_velociraptor_part.h outputlist.h
# Common source files
AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
@@ -60,7 +60,8 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
statistics.c runner_doiact_vec.c profiler.c dump.c logger.c \
part_type.c xmf.c gravity_properties.c gravity.c \
collectgroup.c hydro_space.c equation_of_state.c \
- chemistry.c cosmology.c restart.c mesh_gravity.c velociraptor_interface.c
+ chemistry.c cosmology.c restart.c mesh_gravity.c velociraptor_interface.c \
+ outputlist.c
# Include files for distribution, not installation.
nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
diff --git a/src/cosmology.c b/src/cosmology.c
index e0ea1e9cd30610504821926720a73cb839e37312..40e86b3f7948c1c84f23f0ccfb3571f9c75aae29 100644
--- a/src/cosmology.c
+++ b/src/cosmology.c
@@ -308,6 +308,8 @@ void cosmology_init_tables(struct cosmology *c) {
(double *)malloc(cosmology_table_length * sizeof(double));
c->time_interp_table =
(double *)malloc(cosmology_table_length * sizeof(double));
+ c->scale_factor_interp_table =
+ (double *)malloc(cosmology_table_length * sizeof(double));
/* Prepare a table of scale factors for the integral bounds */
const double delta_a =
@@ -372,6 +374,34 @@ void cosmology_init_tables(struct cosmology *c) {
GSL_INTEG_GAUSS61, space, &result, &abserr);
c->universe_age_at_present_day = result;
+ /* 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;
+
+ int i_prev = 0;
+ for (int i = 0; i < cosmology_table_length; i++) {
+ /* Current time */
+ 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) {
+ i_prev++;
+ }
+
+ /* Find linear interpolation scaling */
+ double scale = time_interp - c->time_interp_table[i_prev - 1];
+ scale /= c->time_interp_table[i_prev] - c->time_interp_table[i_prev - 1];
+ 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;
+ c->scale_factor_interp_table[i] = exp(log_a) - c->a_begin;
+ }
+
/* Free the workspace and temp array */
gsl_integration_workspace_free(space);
free(a_table);
@@ -639,6 +669,24 @@ double cosmology_get_delta_time(const struct cosmology *c,
return t2 - t1;
}
+/**
+ * @brief Compute scale factor from time since big bang (in internal units).
+ *
+ * WARNING: This method has a low accuracy at high redshift.
+ * The relative error is around 1e-3 (testCosmology.c is measuring it).
+ *
+ * @param c The current #cosmology.
+ * @param t time since the big bang
+ * @return The scale factor.
+ */
+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);
+ return a + c->a_begin;
+}
+
/**
* @brief Prints the #cosmology model to stdout.
*/
@@ -660,6 +708,7 @@ void cosmology_clean(struct cosmology *c) {
free(c->grav_kick_fac_interp_table);
free(c->hydro_kick_fac_interp_table);
free(c->time_interp_table);
+ free(c->scale_factor_interp_table);
}
#ifdef HAVE_HDF5
diff --git a/src/cosmology.h b/src/cosmology.h
index 97d62ec98156caf4cf70275798dd8e9da806b4b2..7136b65667195953971060b76ddfd447a5fdf500 100644
--- a/src/cosmology.h
+++ b/src/cosmology.h
@@ -156,6 +156,9 @@ struct cosmology {
/*! Time interpolation table */
double *time_interp_table;
+ /*! Scale factor interpolation table */
+ double *scale_factor_interp_table;
+
/*! Time between Big Bang and first entry in the table */
double time_interp_table_offset;
@@ -180,6 +183,9 @@ 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_scale_factor(const struct cosmology *cosmo, double t);
+
+double cosmology_get_time_since_big_bang(const struct cosmology *c, double a);
void cosmology_init(struct swift_params *params, const struct unit_system *us,
const struct phys_const *phys_const, struct cosmology *c);
diff --git a/src/engine.c b/src/engine.c
index 369dba72b4c9edac550a5726caba4c8311e3ccb0..28dd94ca26f369dfde1cebb03cc020d13dbec3e4 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -69,6 +69,7 @@
#include "map.h"
#include "memswap.h"
#include "minmax.h"
+#include "outputlist.h"
#include "parallel_io.h"
#include "part.h"
#include "partition.h"
@@ -5829,6 +5830,8 @@ void engine_init(struct engine *e, struct space *s, struct swift_params *params,
e->time_base_inv = e->cosmology->time_base_inv;
e->ti_current = 0;
}
+
+ engine_init_output_lists(e, params);
}
/**
@@ -5902,6 +5905,9 @@ void engine_config(int restart, struct engine *e, struct swift_params *params,
error(
"Invalid flag (%d) set for output time format of structure finding.",
e->stf_output_freq_format);
+
+ /* overwrite input if outputlist */
+ if (e->output_list_stf) e->stf_output_freq_format = TIME;
}
/* Get the number of queues */
@@ -6435,6 +6441,12 @@ void engine_print_policy(struct engine *e) {
* @param e The #engine.
*/
void engine_compute_next_snapshot_time(struct engine *e) {
+ /* Do outputlist file case */
+ if (e->output_list_snapshots) {
+ output_list_read_next_time(e->output_list_snapshots, e, "snapshots",
+ &e->ti_next_snapshot);
+ return;
+ }
/* Find upper-bound on last output */
double time_end;
@@ -6493,6 +6505,12 @@ void engine_compute_next_snapshot_time(struct engine *e) {
* @param e The #engine.
*/
void engine_compute_next_statistics_time(struct engine *e) {
+ /* Do output_list file case */
+ if (e->output_list_stats) {
+ output_list_read_next_time(e->output_list_stats, e, "stats",
+ &e->ti_next_stats);
+ return;
+ }
/* Find upper-bound on last output */
double time_end;
@@ -6553,6 +6571,11 @@ void engine_compute_next_statistics_time(struct engine *e) {
* @param e The #engine.
*/
void engine_compute_next_stf_time(struct engine *e) {
+ /* Do output_list file case */
+ if (e->output_list_stf) {
+ output_list_read_next_time(e->output_list_stf, e, "stf", &e->ti_nextSTF);
+ return;
+ }
/* Find upper-bound on last output */
double time_end;
@@ -6602,6 +6625,53 @@ void engine_compute_next_stf_time(struct engine *e) {
}
}
+/**
+ * @brief Initialize all the output_list required by the engine
+ *
+ * @param e The #engine.
+ * @param params The #swift_params.
+ */
+void engine_init_output_lists(struct engine *e, struct swift_params *params) {
+ /* Deal with snapshots */
+ double snaps_time_first;
+ e->output_list_snapshots = NULL;
+ output_list_init(&e->output_list_snapshots, e, "Snapshots",
+ &e->delta_time_snapshot, &snaps_time_first);
+
+ if (e->output_list_snapshots) {
+ if (e->policy & engine_policy_cosmology)
+ e->a_first_snapshot = snaps_time_first;
+ else
+ e->time_first_snapshot = snaps_time_first;
+ }
+
+ /* Deal with stats */
+ double stats_time_first;
+ e->output_list_stats = NULL;
+ output_list_init(&e->output_list_stats, e, "Statistics",
+ &e->delta_time_statistics, &stats_time_first);
+
+ if (e->output_list_stats) {
+ if (e->policy & engine_policy_cosmology)
+ e->a_first_statistics = stats_time_first;
+ else
+ e->time_first_statistics = stats_time_first;
+ }
+
+ /* Deal with stf */
+ double stf_time_first;
+ e->output_list_stf = NULL;
+ output_list_init(&e->output_list_stf, e, "StructureFinding", &e->deltaTimeSTF,
+ &stf_time_first);
+
+ if (e->output_list_stf) {
+ if (e->policy & engine_policy_cosmology)
+ e->a_first_stf = stf_time_first;
+ else
+ e->timeFirstSTFOutput = stf_time_first;
+ }
+}
+
/**
* @brief Computes the maximal time-step allowed by the max RMS displacement
* condition.
@@ -6735,6 +6805,18 @@ void engine_clean(struct engine *e) {
}
free(e->runners);
free(e->snapshot_units);
+ if (e->output_list_snapshots) {
+ output_list_clean(e->output_list_snapshots);
+ free(e->output_list_snapshots);
+ }
+ if (e->output_list_stats) {
+ output_list_clean(e->output_list_stats);
+ free(e->output_list_stats);
+ }
+ if (e->output_list_stf) {
+ output_list_clean(e->output_list_stf);
+ free(e->output_list_stf);
+ }
free(e->links);
free(e->cell_loc);
scheduler_clean(&e->sched);
@@ -6777,6 +6859,11 @@ void engine_struct_dump(struct engine *e, FILE *stream) {
chemistry_struct_dump(e->chemistry, stream);
sourceterms_struct_dump(e->sourceterms, stream);
parser_struct_dump(e->parameter_file, stream);
+ if (e->output_list_snapshots)
+ output_list_struct_dump(e->output_list_snapshots, stream);
+ if (e->output_list_stats)
+ output_list_struct_dump(e->output_list_stats, stream);
+ if (e->output_list_stf) output_list_struct_dump(e->output_list_stf, stream);
}
/**
@@ -6872,6 +6959,27 @@ void engine_struct_restore(struct engine *e, FILE *stream) {
parser_struct_restore(parameter_file, stream);
e->parameter_file = parameter_file;
+ if (e->output_list_snapshots) {
+ struct output_list *output_list_snapshots =
+ (struct output_list *)malloc(sizeof(struct output_list));
+ output_list_struct_restore(output_list_snapshots, stream);
+ e->output_list_snapshots = output_list_snapshots;
+ }
+
+ if (e->output_list_stats) {
+ struct output_list *output_list_stats =
+ (struct output_list *)malloc(sizeof(struct output_list));
+ output_list_struct_restore(output_list_stats, stream);
+ e->output_list_stats = output_list_stats;
+ }
+
+ if (e->output_list_stf) {
+ struct output_list *output_list_stf =
+ (struct output_list *)malloc(sizeof(struct output_list));
+ output_list_struct_restore(output_list_stf, stream);
+ e->output_list_stf = output_list_stf;
+ }
+
/* Want to force a rebuild before using this engine. Wait to repartition.*/
e->forcerebuild = 1;
e->forcerepart = 0;
diff --git a/src/engine.h b/src/engine.h
index a0cd9eab3587e6460f16a3976042b9246dde67f6..aeb57c65ac36ff5ddbf4b74185adeb94f3d460da 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -215,6 +215,9 @@ struct engine {
double time_first_snapshot;
double delta_time_snapshot;
+ /* Output_List for the snapshots */
+ struct output_list *output_list_snapshots;
+
/* Integer time of the next snapshot */
integertime_t ti_next_snapshot;
@@ -231,6 +234,9 @@ struct engine {
double deltaTimeSTF;
int deltaStepSTF;
+ /* Output_List for the structure finding */
+ struct output_list *output_list_stf;
+
/* Integer time of the next stf output */
integertime_t ti_nextSTF;
@@ -241,6 +247,9 @@ struct engine {
double time_first_statistics;
double delta_time_statistics;
+ /* Output_List for the stats */
+ struct output_list *output_list_stats;
+
/* Integer time of the next statistics dump */
integertime_t ti_next_stats;
@@ -373,6 +382,7 @@ void engine_drift_top_multipoles(struct engine *e);
void engine_reconstruct_multipoles(struct engine *e);
void engine_print_stats(struct engine *e);
void engine_dump_snapshot(struct engine *e);
+void engine_init_output_lists(struct engine *e, struct swift_params *params);
void engine_init(struct engine *e, struct space *s, struct swift_params *params,
long long Ngas, long long Ngparts, long long Nstars,
int policy, int verbose, struct repartition *reparttype,
diff --git a/src/outputlist.c b/src/outputlist.c
new file mode 100644
index 0000000000000000000000000000000000000000..d141c9f0b632704dc4e1b245645b44b35419bdc3
--- /dev/null
+++ b/src/outputlist.c
@@ -0,0 +1,276 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * 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
+ * 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 .
+ *
+ ******************************************************************************/
+
+/* Config parameters. */
+#include "../config.h"
+
+/* This object's header. */
+#include "outputlist.h"
+
+/* Local includes. */
+#include "cosmology.h"
+#include "engine.h"
+#include "error.h"
+#include "restart.h"
+
+/* Some standard headers. */
+#include
+
+/**
+ * @brief Read a file containing a list of time
+ *
+ * @param outputlist The #output_list to fill.
+ * @param filename The file to read.
+ * @param cosmo The #cosmology model.
+ */
+void output_list_read_file(struct output_list *outputlist, const char *filename,
+ struct cosmology *cosmo) {
+
+ /* Open file */
+ FILE *file = fopen(filename, "r");
+ if (file == NULL) error("Error opening file '%s'", filename);
+
+ /* Declare reading variables */
+ ssize_t read;
+ size_t len = 0;
+ char *line = NULL;
+
+ /* Count number of lines */
+ size_t nber_line = -1; /* Do not count header */
+ while (getline(&line, &len, file) != -1) {
+ nber_line++;
+ }
+ outputlist->size = nber_line;
+
+ /* Return to start of file and initialize time array */
+ fseek(file, 0, SEEK_SET);
+ outputlist->times = (double *)malloc(sizeof(double) * nber_line);
+ if (!outputlist->times)
+ error(
+ "Unable to malloc output_list. "
+ "Try reducing the number of lines in %s",
+ filename);
+
+ /* Read header */
+ if ((read = getline(&line, &len, file)) == -1)
+ error("Unable to read header in file '%s'", filename);
+
+ /* Remove end of line character */
+ line[strcspn(line, "\n")] = 0;
+
+ /* Find type of data in file */
+ int type = -1;
+ if (!strcmp(line, "# Redshift"))
+ type = OUTPUT_LIST_REDSHIFT;
+ else if (!strcmp(line, "# Time"))
+ type = OUTPUT_LIST_AGE;
+ else if (!strcmp(line, "# Scale Factor"))
+ type = OUTPUT_LIST_SCALE_FACTOR;
+ else
+ error("Unable to interpret the header (%s) in file '%s'", line, filename);
+
+ if (!cosmo &&
+ (type == OUTPUT_LIST_SCALE_FACTOR || type == OUTPUT_LIST_REDSHIFT))
+ 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;
+ while ((read = getline(&line, &len, file)) != -1) {
+ double *time = &outputlist->times[ind];
+ /* 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);
+ }
+
+ /* Transform input into correct time (e.g. ages or scale factor) */
+ if (type == OUTPUT_LIST_REDSHIFT) *time = 1. / (1. + *time);
+
+ if (cosmo && type == OUTPUT_LIST_AGE)
+ *time = cosmology_get_scale_factor(cosmo, *time);
+
+ /* Update size */
+ ind += 1;
+ }
+
+ /* set current indice to 0 */
+ outputlist->cur_ind = 0;
+
+ fclose(file);
+}
+
+/**
+ * @brief Read the next time for an output
+ *
+ * @param t The #output_list
+ * @param e The #engine.
+ * @param name The name of the output (e.g. 'stats', 'snapshots', 'stf')
+ * @param ti_next updated to the next output time
+ */
+void output_list_read_next_time(struct output_list *t, const struct engine *e,
+ const char *name, integertime_t *ti_next) {
+ int is_cosmo = e->policy & engine_policy_cosmology;
+
+ /* 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 = t->times[t->cur_ind];
+ size_t ind = t->cur_ind;
+ while (time < time_end) {
+
+ /* Output time on the integer timeline */
+ if (is_cosmo)
+ *ti_next = log(time / e->cosmology->a_begin) / e->time_base;
+ else
+ *ti_next = (time - e->time_begin) / e->time_base;
+
+ /* Found it? */
+ if (*ti_next > e->ti_current) break;
+
+ ind += 1;
+ if (ind == t->size) break;
+
+ time = t->times[ind];
+ t->cur_ind = ind;
+ }
+
+ /* Deal with last statistics */
+ if (*ti_next >= max_nr_timesteps || ind == t->size || time >= time_end) {
+ *ti_next = -1;
+ if (e->verbose) message("No further output time for %s.", name);
+ } else {
+
+ /* Be nice, talk... */
+ if (is_cosmo) {
+ const double next_time =
+ exp(*ti_next * e->time_base) * e->cosmology->a_begin;
+ if (e->verbose)
+ message("Next output time for %s set to a=%e.", name, next_time);
+ } else {
+ const double next_time = *ti_next * e->time_base + e->time_begin;
+ if (e->verbose)
+ message("Next output time for %s set to t=%e.", name, next_time);
+ }
+ }
+}
+
+/**
+ * @brief initialize an output list
+ *
+ * @param list The output list to initialize
+ * @param e The #engine
+ * @param name The name of the section in params
+ * @param delta_time updated to the initial delta time
+ * @param time_first updated to the time of first output (scale factor or cosmic
+ * time)
+ */
+void output_list_init(struct output_list **list, const struct engine *e,
+ char *name, double *delta_time, double *time_first) {
+ struct swift_params *params = e->parameter_file;
+
+ /* get cosmo */
+ struct cosmology *cosmo = NULL;
+ if (e->policy & engine_policy_cosmology) cosmo = e->cosmology;
+
+ /* Read output on/off */
+ char param_name[PARSER_MAX_LINE_SIZE];
+ sprintf(param_name, "%s:output_list_on", name);
+ int outputlist_on = parser_get_opt_param_int(params, param_name, 0);
+
+ /* Check if read outputlist */
+ if (!outputlist_on) return;
+
+ /* Read outputlist for snapshots */
+ *list = (struct output_list *)malloc(sizeof(struct output_list));
+
+ /* Read filename */
+ char filename[PARSER_MAX_LINE_SIZE];
+ sprintf(param_name, "%s:output_list", name);
+ parser_get_param_string(params, param_name, filename);
+
+ message("Reading %s output file.", name);
+ output_list_read_file(*list, filename, cosmo);
+
+ if ((*list)->size < 2)
+ error("You need to provide more snapshots in '%s'", filename);
+
+ /* Set data for later checks */
+ if (cosmo) {
+ *delta_time = (*list)->times[1] / (*list)->times[0];
+ *time_first = (*list)->times[0];
+ } else {
+ *delta_time = (*list)->times[1] - (*list)->times[0];
+ *time_first = (*list)->times[0];
+ }
+}
+
+/**
+ * @brief Print an #output_list
+ */
+void output_list_print(const struct output_list *outputlist) {
+
+ printf("/*\t Time Array\t */\n");
+ printf("Number of Line: %lu\n", outputlist->size);
+ for (size_t ind = 0; ind < outputlist->size; ind++) {
+ if (ind == outputlist->cur_ind)
+ printf("\t%lf <-- Current\n", outputlist->times[ind]);
+ else
+ printf("\t%lf\n", outputlist->times[ind]);
+ }
+}
+
+/**
+ * @brief Clean an #output_list
+ */
+void output_list_clean(struct output_list *outputlist) {
+ free(outputlist->times);
+}
+
+/**
+ * @brief Dump an #output_list in a restart file
+ */
+void output_list_struct_dump(struct output_list *list, FILE *stream) {
+ restart_write_blocks(list, sizeof(struct output_list), 1, stream,
+ "output_list", "output_list struct");
+
+ restart_write_blocks(list->times, list->size, sizeof(double), stream,
+ "output_list", "times");
+}
+
+/**
+ * @brief Restore an #output_list from a restart file
+ */
+void output_list_struct_restore(struct output_list *list, FILE *stream) {
+ restart_read_blocks(list, sizeof(struct output_list), 1, stream, NULL,
+ "output_list 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
new file mode 100644
index 0000000000000000000000000000000000000000..ddc08c84d3739097e990c175cd83e809231845a9
--- /dev/null
+++ b/src/outputlist.h
@@ -0,0 +1,65 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * 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
+ * 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 .
+ *
+ ******************************************************************************/
+#ifndef SWIFT_OUTPUT_LIST_H
+#define SWIFT_OUTPUT_LIST_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Local includes */
+#include "cosmology.h"
+
+struct engine;
+
+/**
+ * @brief the different output_list type
+ */
+enum output_list_type {
+ OUTPUT_LIST_AGE,
+ OUTPUT_LIST_REDSHIFT,
+ OUTPUT_LIST_SCALE_FACTOR,
+};
+
+/**
+ * @brief the array containing the output times
+ */
+struct output_list {
+
+ /* Time array */
+ double *times;
+
+ /* Size of the time array */
+ size_t size;
+
+ /* Current index */
+ size_t cur_ind;
+};
+
+void output_list_read_file(struct output_list *outputlist, const char *filename,
+ struct cosmology *cosmo);
+void output_list_read_next_time(struct output_list *t, const struct engine *e,
+ const char *name, integertime_t *ti_next);
+void output_list_init(struct output_list **list, const struct engine *e,
+ char *name, double *delta_time, double *time_first);
+void output_list_print(const struct output_list *outputlist);
+void output_list_clean(struct output_list *outputlist);
+void output_list_struct_dump(struct output_list *list, FILE *stream);
+void output_list_struct_restore(struct output_list *list, FILE *stream);
+
+#endif /* SWIFT_OUTPUT_LIST_H */
diff --git a/src/swift.h b/src/swift.h
index 3222376888bada072114f58ac5294f4567e8659a..e10938addb99956c202b3e4dd2b0592b580fa948 100644
--- a/src/swift.h
+++ b/src/swift.h
@@ -49,6 +49,7 @@
#include "map.h"
#include "mesh_gravity.h"
#include "multipole.h"
+#include "outputlist.h"
#include "parallel_io.h"
#include "parser.h"
#include "part.h"
diff --git a/tests/Makefile.am b/tests/Makefile.am
index ad7dcf92d6e4133e80bdd077447841626e59d722..2744164e2c794ef4471a51ffb66db3049d9fa048 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -28,7 +28,7 @@ TESTS = testGreetings testMaths testReading.sh testSingle testKernel testSymmetr
testVoronoi1D testVoronoi2D testVoronoi3D testGravityDerivatives \
testPeriodicBC.sh testPeriodicBCPerturbed.sh testPotentialSelf \
testPotentialPair testEOS testUtilities testSelectOutput.sh \
- testCbrt
+ testCbrt testCosmology
# List of test programs to compile
check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
@@ -39,7 +39,7 @@ check_PROGRAMS = testGreetings testReading testSingle testTimeIntegration \
testRiemannHLLC testMatrixInversion testDump testLogger \
testVoronoi1D testVoronoi2D testVoronoi3D testPeriodicBC \
testGravityDerivatives testPotentialSelf testPotentialPair testEOS testUtilities \
- testSelectOutput testCbrt
+ testSelectOutput testCbrt testCosmology
# Rebuild tests when SWIFT is updated.
$(check_PROGRAMS): ../src/.libs/libswiftsim.a
@@ -53,6 +53,8 @@ testReading_SOURCES = testReading.c
testSelectOutput_SOURCES = testSelectOutput.c
+testCosmology_SOURCES = testCosmology.c
+
testSymmetry_SOURCES = testSymmetry.c
# Added because of issues using memcmp on clang 4.x
diff --git a/tests/testCosmology.c b/tests/testCosmology.c
new file mode 100644
index 0000000000000000000000000000000000000000..698351ad952e7d0b5f7d8e354c45a1a2dd53f968
--- /dev/null
+++ b/tests/testCosmology.c
@@ -0,0 +1,76 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * 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
+ * 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 .
+ *
+ ******************************************************************************/
+
+/* Some standard headers. */
+#include
+
+/* Includes. */
+#include "swift.h"
+
+#define N_CHECK 20
+#define TOLERANCE 1e-3
+
+void test_params_init(struct swift_params *params) {
+ parser_init("", params);
+ parser_set_param(params, "Cosmology:Omega_m:0.3075");
+ parser_set_param(params, "Cosmology:Omega_lambda:0.6910");
+ parser_set_param(params, "Cosmology:Omega_b:0.0486");
+ parser_set_param(params, "Cosmology:h:0.6774");
+ parser_set_param(params, "Cosmology:a_begin:0.1");
+ parser_set_param(params, "Cosmology:a_end:1.0");
+}
+
+int main(int argc, char *argv[]) {
+
+ message("Initialization...");
+
+ /* pseudo initialization of params */
+ struct swift_params params;
+ test_params_init(¶ms);
+
+ /* initialization of unit system */
+ struct unit_system us;
+ units_init_cgs(&us);
+
+ /* initialization of phys_const */
+ struct phys_const phys_const;
+ phys_const_init(&us, ¶ms, &phys_const);
+
+ /* initialization of cosmo */
+ struct cosmology cosmo;
+ cosmology_init(¶ms, &us, &phys_const, &cosmo);
+
+ message("Start checking computation...");
+
+ 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);
+ tmp = cosmology_get_scale_factor(&cosmo, tmp);
+
+ /* check accuracy */
+ tmp = (tmp - a) / a;
+ message("Accuracy of %g at a=%g", tmp, a);
+ assert(fabs(tmp) < TOLERANCE);
+ }
+
+ message("Everything seems fine with cosmology.");
+
+ return 0;
+}