/*******************************************************************************
* 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
/* This object's header. */
#include "output_list.h"
/* Local includes. */
#include "cosmology.h"
#include "engine.h"
#include "error.h"
#include "restart.h"
#include "tools.h"
/* Some standard headers. */
#include
/**
* @brief Read a file containing a list of time
*
* @param output_list The #output_list to fill.
* @param filename The file to read.
* @param cosmo The #cosmology model.
*/
void output_list_read_file(struct output_list *output_list,
const char *filename, struct cosmology *cosmo) {
/* Open file */
FILE *file = fopen(filename, "r");
if (file == NULL) error("Error opening file '%s'", filename);
/* Count number of lines */
size_t len = 0;
char *line = NULL;
size_t nber_line = 0;
while (getline(&line, &len, file) != -1) nber_line++;
output_list->size = nber_line - 1; /* Do not count header */
/* Return to start of file and initialize time array */
fseek(file, 0, SEEK_SET);
output_list->times = (double *)malloc(sizeof(double) * output_list->size);
output_list->snapshot_labels = (int *)malloc(sizeof(int) * output_list->size);
output_list->select_output_indices =
(int *)malloc(sizeof(int) * output_list->size);
if ((!output_list->times) || (!output_list->select_output_indices)) {
error(
"Unable to malloc output_list. "
"Try reducing the number of lines in %s",
filename);
}
/* Read header */
if (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;
output_list->select_output_on = 0;
output_list->select_output_number_of_names = 0;
output_list->alternative_labels_on = 0;
trim_trailing(line);
if (strcasecmp(line, "# Redshift") == 0) {
type = OUTPUT_LIST_REDSHIFT;
} else if (strcasecmp(line, "# Time") == 0) {
type = OUTPUT_LIST_AGE;
} else if (strcasecmp(line, "# Scale Factor") == 0) {
type = OUTPUT_LIST_SCALE_FACTOR;
} else if (strcasecmp(line, "# Comoving Distance") == 0) {
type = OUTPUT_LIST_COMOVING_DISTANCE;
} else if (strcasecmp(line, "# Redshift, Select Output") == 0) {
type = OUTPUT_LIST_REDSHIFT;
output_list->select_output_on = 1;
} else if (strcasecmp(line, "# Time, Select Output") == 0) {
type = OUTPUT_LIST_AGE;
output_list->select_output_on = 1;
} else if (strcasecmp(line, "# Scale Factor, Select Output") == 0) {
type = OUTPUT_LIST_SCALE_FACTOR;
output_list->select_output_on = 1;
} else if (strcasecmp(line, "# Comoving Distance, Select Output") == 0) {
type = OUTPUT_LIST_COMOVING_DISTANCE;
output_list->select_output_on = 1;
} else if (strcasecmp(line, "# Redshift, Select Output, Label") == 0) {
type = OUTPUT_LIST_REDSHIFT;
output_list->select_output_on = 1;
output_list->alternative_labels_on = 1;
} else if (strcasecmp(line, "# Time, Select Output, Label") == 0) {
type = OUTPUT_LIST_AGE;
output_list->select_output_on = 1;
output_list->alternative_labels_on = 1;
} else if (strcasecmp(line, "# Scale Factor, Select Output, Label") == 0) {
type = OUTPUT_LIST_SCALE_FACTOR;
output_list->select_output_on = 1;
output_list->alternative_labels_on = 1;
} else if (strcasecmp(line, "# Comoving Distance, Select Output, Label") ==
0) {
type = OUTPUT_LIST_COMOVING_DISTANCE;
output_list->select_output_on = 1;
output_list->alternative_labels_on = 1;
} else {
error("Unable to interpret the header (%s) in file '%s'", line, filename);
}
if (!cosmo &&
(type == OUTPUT_LIST_SCALE_FACTOR || type == OUTPUT_LIST_REDSHIFT ||
type == OUTPUT_LIST_COMOVING_DISTANCE))
error(
"Unable to compute a redshift or a scale factor without cosmology. "
"Please change the header in '%s'",
filename);
if (!output_list->select_output_on && output_list->alternative_labels_on)
error(
"Found an output list with alternative labels but not individual "
"output selections");
/* Read file */
size_t ind = 0;
int read_successfully = 0;
int found_select_output = 0;
char select_output_buffer[FIELD_BUFFER_SIZE] =
select_output_header_default_name;
while (getline(&line, &len, file) != -1) {
double *time = &output_list->times[ind];
int *label = &output_list->snapshot_labels[ind];
/* Write data to output_list */
if (output_list->select_output_on && output_list->alternative_labels_on) {
read_successfully = sscanf(line, "%lf, %[^,], %d", time,
select_output_buffer, label) == 3;
} else if (output_list->select_output_on) {
read_successfully =
sscanf(line, "%lf, %s", time, select_output_buffer) == 2;
} else {
read_successfully = sscanf(line, "%lf", time) == 1;
}
if (!read_successfully) {
error(
"Tried parsing output_list but found '%s' with illegal "
"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);
if (cosmo && type == OUTPUT_LIST_COMOVING_DISTANCE)
*time = cosmology_scale_factor_at_comoving_distance(cosmo, *time);
/* Search to find index for select output - select_output_index is the index
* in the select_output_names array that corresponds to this select output
* name. */
found_select_output = 0;
for (int i = 0; i < output_list->select_output_number_of_names; i++) {
if (strcmp(select_output_buffer, output_list->select_output_names[i]) ==
0) {
/* We already have this select output list string in the buffer! */
output_list->select_output_indices[ind] = i;
found_select_output = 1;
}
}
/* If we did not assign it above, we haven't encountered this name before
* and we need to create this name in the array */
if (!found_select_output) {
strcpy(
output_list
->select_output_names[output_list->select_output_number_of_names],
select_output_buffer);
output_list->select_output_indices[ind] =
output_list->select_output_number_of_names;
output_list->select_output_number_of_names += 1;
}
/* Update size */
ind += 1;
}
/* Cleanup */
free(line);
if (ind != output_list->size)
error("Did not read the correct number of output times.");
/* Check that the list is in monotonic order */
for (size_t i = 1; i < output_list->size; ++i) {
if ((type == OUTPUT_LIST_REDSHIFT) &&
(output_list->times[i] <= output_list->times[i - 1]))
error("Output list not having monotonically decreasing redshifts.");
if ((type == OUTPUT_LIST_AGE) &&
(output_list->times[i] <= output_list->times[i - 1]))
error("Output list not having monotonically increasing ages.");
if ((type == OUTPUT_LIST_SCALE_FACTOR) &&
(output_list->times[i] <= output_list->times[i - 1]))
error("Output list not having monotonically increasing scale-factors.");
if ((type == OUTPUT_LIST_COMOVING_DISTANCE) &&
(output_list->times[i] <= output_list->times[i - 1]))
error(
"Output list not having monotonically decreasing comoving distance.");
}
/* set current indice to 0 */
output_list->cur_ind = 0;
/* We check if this is true later */
output_list->final_step_dump = 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 snapshot 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;
}
/* Do we need to do a dump at the end of the last timestep?
* Note that what this really does is given that a t=t_max, z=0,
* or a=1 is found in output_list.txt set the flag `final_step_dump`
* to 1 - this is not special behaviour that is controlled by a
* parameter file flag. */
if (time == time_end ||
(time > time_end && time - time_end < OUTPUT_LIST_EPS_TIME_END)) {
t->final_step_dump = 1;
if (e->verbose) {
if (is_cosmo) {
message("Next output time for %s set to a=%e.", name, time_end);
} else {
message("Next output time for %s set to t=%e.", name, time_end);
}
}
}
/* Deal with last statistics */
if (*ti_next >= max_nr_timesteps || ind == t->size || time >= time_end) {
*ti_next = -1;
if (e->verbose && t->final_step_dump != 1) {
message("No further output time for %s.", name);
/* Do not print anything about the output style (below) */
return;
}
} 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);
}
}
/* Finally, talk if we are a snapshot and we are using SelectOutput */
if (e->verbose) {
if (t->select_output_number_of_names > 1) {
char select_output_style[FIELD_BUFFER_SIZE];
output_list_get_current_select_output(t, select_output_style);
message("Next output style for %s set to %s.", name, select_output_style);
}
}
}
/**
* @brief Copys the string describing the current output name into the
* buffer described by select_output_name.
*
* @param t The #output_list
* @param select_output_name updated to the current Select Output choice.
**/
void output_list_get_current_select_output(struct output_list *t,
char *select_output_name) {
strcpy(select_output_name,
t->select_output_names[t->select_output_indices[t->cur_ind]]);
}
/**
* @brief initialize an output list
*
* @param list The output list to initialize.
* @param e The #engine.
* @param name The name of the section in the param file.
* @param delta_time (return) The delta between the first two outputs
*/
void output_list_init(struct output_list **list, const struct engine *e,
const char *name, double *const delta_time) {
struct swift_params *params = e->parameter_file;
if (*list != NULL) error("Output list already allocated!");
/* 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);
const int output_list_on = parser_get_opt_param_int(params, param_name, 0);
/* Check if read output_list */
if (!output_list_on) return;
/* Read output_list for snapshots */
*list = (struct output_list *)malloc(sizeof(struct output_list));
bzero(*list, sizeof(struct output_list));
(*list)->output_list_on = output_list_on;
/* Read filename */
char filename[PARSER_MAX_LINE_SIZE];
sprintf(param_name, "%s:output_list", name);
parser_get_param_string(params, param_name, filename);
if (e->verbose) message("Reading %s output file.", name);
output_list_read_file(*list, filename, cosmo);
if ((*list)->size < 2)
error("You need to provide more entries in '%s'", filename);
if (strcmp(name, "Snapshots") == 0) {
if ((*list)->select_output_on &&
e->output_options->select_output->paramCount == 0)
error(
"Cannot use an output list with individual output selection entries "
"if 'Snapshots:select_output_on' is set to 0.");
}
/* Set data for later checks */
if (cosmo) {
*delta_time = (*list)->times[1] / (*list)->times[0];
} else {
*delta_time = (*list)->times[1] - (*list)->times[0];
}
}
/**
* @brief Verify that a given output list only contains output selections
* that have actually been read in.
*
* @param list the #output_list.
* @param output_options The #output_options we read in earlier.
*/
void output_list_check_selection(const struct output_list *list,
const struct output_options *output_options) {
for (size_t i = 0; i < list->size; ++i) {
const int list_id = list->select_output_indices[i];
const int selection_id = parser_get_section_id(
output_options->select_output, list->select_output_names[list_id]);
if (selection_id < 0)
error("Invalid output selection found in output list: '%s'",
list->select_output_names[list_id]);
}
}
/**
* @brief Print an #output_list
*/
void output_list_print(const struct output_list *output_list) {
printf("/*\t Total number of Select Output options: %d \t */\n",
output_list->select_output_number_of_names);
printf("/*\t Select Output options found in output list: \t */\n");
for (int i = 0; i < output_list->select_output_number_of_names; i++) {
printf("\t Index: %d, Name: %s\n", i, output_list->select_output_names[i]);
}
printf("\n/*\t Time Array, Select Output \t */\n");
printf("Number of Lines: %zu\n", output_list->size);
for (size_t ind = 0; ind < output_list->size; ind++) {
if (ind == output_list->cur_ind)
printf(
"\t %lf, %s <-- Current\n", output_list->times[ind],
output_list
->select_output_names[output_list->select_output_indices[ind]]);
else
printf(
"\t %lf, %s\n", output_list->times[ind],
output_list
->select_output_names[output_list->select_output_indices[ind]]);
}
}
/**
* @brief Clean an #output_list
*/
void output_list_clean(struct output_list **output_list) {
if (*output_list) {
free((*output_list)->times);
free((*output_list)->snapshot_labels);
free((*output_list)->select_output_indices);
free(*output_list);
*output_list = NULL;
}
}
/**
* @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");
restart_write_blocks(list->select_output_indices, list->size, sizeof(int),
stream, "output_list", "select_output_indices");
restart_write_blocks(list->snapshot_labels, list->size, sizeof(int), stream,
"output_list", "snapshot_labels");
}
/**
* @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");
list->select_output_indices = (int *)malloc(sizeof(int) * list->size);
restart_read_blocks(list->select_output_indices, list->size, sizeof(int),
stream, NULL, "select_output_indices");
list->snapshot_labels = (int *)malloc(sizeof(int) * list->size);
restart_read_blocks(list->snapshot_labels, list->size, sizeof(int), stream,
NULL, "snapshot_labels");
}