diff --git a/logger/Makefile.am b/logger/Makefile.am
index c724ebb85360fa6543b66d27adf3bf1bdaffcf57..9623ffc9e203d7d0645aed89400c7183b5b6501a 100644
--- a/logger/Makefile.am
+++ b/logger/Makefile.am
@@ -80,12 +80,13 @@ endif
SUBDIRS = tests
# List required headers
-include_HEADERS = logger_header.h logger_loader_io.h logger_particle.h logger_time.h logger_tools.h logger_reader.h \
- logger_logfile.h logger_index.h quick_sort.h logger_python_tools.h logger_interpolation.h
+include_HEADERS = logger_header.h logger_loader_io.h logger_particle.h logger_time.h logger_tools.h
+include_HEADERS += logger_reader.h logger_logfile.h logger_index.h quick_sort.h logger_python_tools.h
+include_HEADERS += logger_interpolation.h logger_parameters.h
# Common source files
AM_SOURCES = logger_header.c logger_loader_io.c logger_particle.c logger_time.c logger_tools.c logger_reader.c
-AM_SOURCES += logger_logfile.c logger_index.c quick_sort.c
+AM_SOURCES += logger_logfile.c logger_index.c quick_sort.c logger_parameters.c
AM_SOURCES += $(GRAVITY_SRC) $(STARS_SRC) $(HYDRO_SRC) $(CHEMISTRY_SRC) $(STAR_FORMATION_SRC)
if HAVEPYTHON
diff --git a/logger/chemistry/GEAR/logger_chemistry.h b/logger/chemistry/GEAR/logger_chemistry.h
index f0fbdd6119fc9689e7dfd301e31f97378f621ea9..635e0d65d5cd4d104f4de2cc052b590cf9583639 100644
--- a/logger/chemistry/GEAR/logger_chemistry.h
+++ b/logger/chemistry/GEAR/logger_chemistry.h
@@ -25,6 +25,7 @@
#include "chemistry_logger.h"
#include "logger_interpolation.h"
#include "logger_loader_io.h"
+#include "logger_parameters.h"
#include "logger_python_tools.h"
/* Index of the mask in the header mask array */
@@ -70,12 +71,14 @@ chemistry_logger_reader_link_derivatives_spart(struct header *head) {}
* @param t Requested time.
* @param field The field to reconstruct (follows the order of
* #chemistry_logger_fields).
+ * @param params The simulation's #logger_parameters.
*/
__attribute__((always_inline)) INLINE static void
chemistry_logger_interpolate_field_part(
const double t_before, const struct logger_field *restrict before,
const double t_after, const struct logger_field *restrict after,
- void *restrict output, const double t, const int field) {
+ void *restrict output, const double t, const int field,
+ const struct logger_parameters *params) {
#ifdef SWIFT_DEBUG_CHECKS
/* Check the times */
if (t_before > t || t_after < t) {
@@ -112,12 +115,14 @@ chemistry_logger_interpolate_field_part(
* @param t Requested time.
* @param field The field to reconstruct (follows the order of
* #chemistry_logger_fields).
+ * @param params The simulation's #logger_parameters.
*/
__attribute__((always_inline)) INLINE static void
chemistry_logger_interpolate_field_spart(
const double t_before, const struct logger_field *restrict before,
const double t_after, const struct logger_field *restrict after,
- void *restrict output, const double t, const int field) {
+ void *restrict output, const double t, const int field,
+ const struct logger_parameters *params) {
#ifdef SWIFT_DEBUG_CHECKS
/* Check the times */
@@ -132,7 +137,8 @@ chemistry_logger_interpolate_field_spart(
switch (field) {
case chemistry_logger_field_spart_metal_mass_fractions:
interpolate_linear_double_ND(t_before, before, t_after, after, output, t,
- GEAR_CHEMISTRY_ELEMENT_COUNT);
+ GEAR_CHEMISTRY_ELEMENT_COUNT,
+ /* periodic= */ 0, params);
break;
default:
diff --git a/logger/chemistry/none/logger_chemistry.h b/logger/chemistry/none/logger_chemistry.h
index e776cb352d514e3bbf0a50fdaaed682941b0f18d..8a2a60d161ba77d0ba10ae07f0012aa2507aa294 100644
--- a/logger/chemistry/none/logger_chemistry.h
+++ b/logger/chemistry/none/logger_chemistry.h
@@ -25,6 +25,7 @@
#include "chemistry_logger.h"
#include "logger_interpolation.h"
#include "logger_loader_io.h"
+#include "logger_parameters.h"
#include "logger_python_tools.h"
/* Index of the mask in the header mask array */
@@ -70,12 +71,14 @@ chemistry_logger_reader_link_derivatives_spart(struct header *head) {}
* @param t Requested time.
* @param field The field to reconstruct (follows the order of
* #chemistry_logger_fields).
+ * @param params The simulation's #logger_parameters.
*/
__attribute__((always_inline)) INLINE static void
chemistry_logger_interpolate_field_part(
const double t_before, const struct logger_field *restrict before,
const double t_after, const struct logger_field *restrict after,
- void *restrict output, const double t, const int field) {}
+ void *restrict output, const double t, const int field,
+ const struct logger_parameters *params) {}
/**
* @brief Interpolate a field of the #spart at the given time.
@@ -92,12 +95,14 @@ chemistry_logger_interpolate_field_part(
* @param t Requested time.
* @param field The field to reconstruct (follows the order of
* #chemistry_logger_fields).
+ * @param params The simulation's #logger_parameters.
*/
__attribute__((always_inline)) INLINE static void
chemistry_logger_interpolate_field_spart(
const double t_before, const struct logger_field *restrict before,
const double t_after, const struct logger_field *restrict after,
- void *restrict output, const double t, const int field) {}
+ void *restrict output, const double t, const int field,
+ const struct logger_parameters *params) {}
#ifdef HAVE_PYTHON
/**
diff --git a/logger/gravity/MultiSoftening/logger_gravity.h b/logger/gravity/MultiSoftening/logger_gravity.h
index c5cfcdc07e476f8c04191fe139b446ed281fd107..07154fc24a3a55e33df8b3269a8a13493f120209 100644
--- a/logger/gravity/MultiSoftening/logger_gravity.h
+++ b/logger/gravity/MultiSoftening/logger_gravity.h
@@ -25,6 +25,7 @@
#include "gravity_logger.h"
#include "logger_interpolation.h"
#include "logger_loader_io.h"
+#include "logger_parameters.h"
#include "logger_python_tools.h"
/* Index of the mask in the header mask array */
@@ -72,6 +73,7 @@ gravity_logger_reader_link_derivatives(struct header *head) {
* @param t Requested time.
* @param field The field to reconstruct (follows the order of
* #gravity_logger_fields).
+ * @param params The simulation's #logger_parameters.
*/
__attribute__((always_inline)) INLINE static void
gravity_logger_interpolate_field(const double t_before,
@@ -79,7 +81,8 @@ gravity_logger_interpolate_field(const double t_before,
const double t_after,
const struct logger_field *restrict after,
void *restrict output, const double t,
- const int field) {
+ const int field,
+ const struct logger_parameters *params) {
#ifdef SWIFT_DEBUG_CHECKS
/* Check the times */
@@ -94,18 +97,21 @@ gravity_logger_interpolate_field(const double t_before,
switch (field) {
case gravity_logger_field_coordinates:
interpolate_quintic_double_float_ND(t_before, before, t_after, after,
- output, t, /* dimension= */ 3);
+ output, t, /* dimension= */ 3,
+ params->periodic, params);
break;
case gravity_logger_field_velocities:
interpolate_cubic_float_ND(t_before, before, t_after, after, output, t,
- /* dimension= */ 3);
+ /* dimension= */ 3, /* periodic= */ 0, params);
break;
case gravity_logger_field_accelerations:
interpolate_linear_float_ND(t_before, before, t_after, after, output, t,
- /* dimension= */ 3);
+ /* dimension= */ 3, /* periodic= */ 0,
+ params);
break;
case gravity_logger_field_masses:
- interpolate_linear_float(t_before, before, t_after, after, output, t);
+ interpolate_linear_float(t_before, before, t_after, after, output, t,
+ /* periodic= */ 0, params);
break;
case gravity_logger_field_particle_ids:
interpolate_ids(t_before, before, t_after, after, output, t);
diff --git a/logger/hydro/Gadget2/logger_hydro.h b/logger/hydro/Gadget2/logger_hydro.h
index affa61e8ec7bb6491f99ffdb134a0fa3bd0bf84e..c8f5efe2c4897c78ff81979aab1d7c6ac89ec2e0 100644
--- a/logger/hydro/Gadget2/logger_hydro.h
+++ b/logger/hydro/Gadget2/logger_hydro.h
@@ -25,6 +25,7 @@
#include "hydro_logger.h"
#include "logger_interpolation.h"
#include "logger_loader_io.h"
+#include "logger_parameters.h"
#include "logger_python_tools.h"
/* Index of the mask in the header mask array */
@@ -72,6 +73,7 @@ hydro_logger_reader_link_derivatives(struct header *head) {
* @param t Requested time.
* @param field The field to reconstruct (follows the order of
* #hydro_logger_fields).
+ * @param params The simulation's #logger_parameters.
*/
__attribute__((always_inline)) INLINE static void
hydro_logger_interpolate_field(const double t_before,
@@ -79,7 +81,8 @@ hydro_logger_interpolate_field(const double t_before,
const double t_after,
const struct logger_field *restrict after,
void *restrict output, const double t,
- const int field) {
+ const int field,
+ const struct logger_parameters *params) {
#ifdef SWIFT_DEBUG_CHECKS
/* Check the times */
@@ -95,23 +98,25 @@ hydro_logger_interpolate_field(const double t_before,
/* Do the position */
case hydro_logger_field_coordinates:
interpolate_quintic_double_float_ND(t_before, before, t_after, after,
- output, t, /* dimesion= */ 3);
+ output, t, /* dimesion= */ 3,
+ params->periodic, params);
break;
/* Do the velocity */
case hydro_logger_field_velocities:
interpolate_cubic_float_ND(t_before, before, t_after, after, output, t,
- /* dimesion= */ 3);
+ /* dimesion= */ 3, /* periodic= */ 0, params);
break;
case hydro_logger_field_accelerations:
interpolate_linear_float_ND(t_before, before, t_after, after, output, t,
- /* dimesion= */ 3);
+ /* dimesion= */ 3, /* periodic= */ 0, params);
break;
/* Do the linear interpolation of float. */
case hydro_logger_field_smoothing_lengths:
case hydro_logger_field_entropies:
case hydro_logger_field_densities:
case hydro_logger_field_masses:
- interpolate_linear_float(t_before, before, t_after, after, output, t);
+ interpolate_linear_float(t_before, before, t_after, after, output, t,
+ /* periodic= */ 0, params);
break;
/* Check the ids */
case hydro_logger_field_particle_ids:
diff --git a/logger/hydro/SPHENIX/logger_hydro.h b/logger/hydro/SPHENIX/logger_hydro.h
index 7f4288f24e269ced02d4d066fde43a1e01155538..bb60efd67e58ab4c54594543d00eb6e14a55496e 100644
--- a/logger/hydro/SPHENIX/logger_hydro.h
+++ b/logger/hydro/SPHENIX/logger_hydro.h
@@ -26,6 +26,7 @@
#include "hydro_logger.h"
#include "logger_interpolation.h"
#include "logger_loader_io.h"
+#include "logger_parameters.h"
#include "logger_python_tools.h"
/* Index of the mask in the header mask array */
@@ -73,6 +74,7 @@ hydro_logger_reader_link_derivatives(struct header *head) {
* @param t Requested time.
* @param field The field to reconstruct (follows the order of
* #hydro_logger_fields).
+ * @param params The simulation's #logger_parameters.
*/
__attribute__((always_inline)) INLINE static void
hydro_logger_interpolate_field(const double t_before,
@@ -80,7 +82,8 @@ hydro_logger_interpolate_field(const double t_before,
const double t_after,
const struct logger_field *restrict after,
void *restrict output, const double t,
- const int field) {
+ const int field,
+ const struct logger_parameters *params) {
#ifdef SWIFT_DEBUG_CHECKS
/* Check the times */
@@ -96,23 +99,26 @@ hydro_logger_interpolate_field(const double t_before,
/* Do the position */
case hydro_logger_field_coordinates:
interpolate_quintic_double_float_ND(t_before, before, t_after, after,
- output, t, /* dimension= */ 3);
+ output, t, /* dimension= */ 3,
+ params->periodic, params);
break;
/* Do the velocity */
case hydro_logger_field_velocities:
interpolate_cubic_float_ND(t_before, before, t_after, after, output, t,
- /* dimension= */ 3);
+ /* dimension= */ 3, /* periodic= */ 0, params);
break;
case hydro_logger_field_accelerations:
interpolate_linear_float_ND(t_before, before, t_after, after, output, t,
- /* dimension= */ 3);
+ /* dimension= */ 3, /* periodic= */ 0,
+ params);
break;
/* Do the linear interpolation of float. */
case hydro_logger_field_masses:
case hydro_logger_field_smoothing_lengths:
case hydro_logger_field_internal_energies:
case hydro_logger_field_densities:
- interpolate_linear_float(t_before, before, t_after, after, output, t);
+ interpolate_linear_float(t_before, before, t_after, after, output, t,
+ /* periodic= */ 0, params);
break;
/* Check the ids */
case hydro_logger_field_particle_ids:
@@ -136,10 +142,12 @@ hydro_logger_interpolate_field(const double t_before,
/* Div v */
x[n_linear] = interpolate_cubic_hermite_spline(
t_before, bef[n_linear], bef[n_linear + 1], t_after, aft[n_linear],
- aft[n_linear + 1], t);
+ aft[n_linear + 1], t, /* periodic= */ 0, params);
/* d Div v / dt */
x[n_linear + 1] = wa * aft[n_linear + 1] + wb * bef[n_linear + 1];
+ /* Use the linear interpolation */
+ x[1] = wa * div_aft[1] + wb * div_bef[1];
break;
}
diff --git a/logger/logger_index.c b/logger/logger_index.c
index b4bfe7adb5f62b70cc46cf73e07a09e97bdb13d3..606f7cd8c059b313ac71efab51a79c8bc0aaaa45 100644
--- a/logger/logger_index.c
+++ b/logger/logger_index.c
@@ -208,7 +208,7 @@ struct index_data *logger_index_get_removed_history(struct logger_index *index,
*/
int logger_index_contains_time_array(struct logger_index *index) {
/* Only the first index file should have a time array */
- if (index->time != 0.) {
+ if (index->integer_time != 0) {
error_python("Only the first index file can have a time array.");
}
diff --git a/logger/logger_interpolation.h b/logger/logger_interpolation.h
index 7bc125e1c38fe742d19e7c1a34d9d295d5b793c1..775059aa3086d9ab288425e925331a7c82cb0d7d 100644
--- a/logger/logger_interpolation.h
+++ b/logger/logger_interpolation.h
@@ -20,6 +20,7 @@
#ifndef LOGGER_LOGGER_INTERPOLATION_H
#define LOGGER_LOGGER_INTERPOLATION_H
+#include "logger_parameters.h"
#include "logger_tools.h"
/**
@@ -151,6 +152,8 @@ interpolate_cubic_hermite_spline(const double t0, const float v0,
* @param t_before Time of field_before (< t).
* @param t_after Time of field_after (> t).
* @param t Requested time.
+ * @param periodic Should the periodic boundary be applied?
+ * @param params The simulation's #logger_parameters.
*/
__attribute__((always_inline)) INLINE static void
interpolate_quintic_double_float_ND(const double t_before,
@@ -158,35 +161,72 @@ interpolate_quintic_double_float_ND(const double t_before,
const double t_after,
const struct logger_field *restrict after,
void *restrict output, const double t,
- const int dimension) {
- /* Compute the interpolation scaling. */
- const double wa = (t - t_before) / (t_after - t_before);
- const double wb = 1. - wa;
+ const int dimension, int periodic,
+ const struct logger_parameters *params) {
- for (int i = 0; i < dimension; i++) {
- double *x = (double *)output;
- const double *x_bef = (double *)before->field;
- const float *v_bef = (float *)before->first_deriv;
- const float *a_bef = (float *)before->second_deriv;
+ /* Get the arrays */
+ const double *x_bef = (double *)before->field;
+ const float *v_bef = (float *)before->first_deriv;
+ const float *a_bef = (float *)before->second_deriv;
+
+ const double *x_aft = (double *)after->field;
+ const float *v_aft = (float *)after->first_deriv;
+ const float *a_aft = (float *)after->second_deriv;
- const double *x_aft = (double *)after->field;
- const float *v_aft = (float *)after->first_deriv;
- const float *a_aft = (float *)after->second_deriv;
+ double *x = (double *)output;
+
+ /* Loop over the dimensions */
+ for (int i = 0; i < dimension; i++) {
+ /* Get the current variables */
+ double x0 = x_bef[i];
+ double x1 = x_aft[i];
+ float v0 = v_bef ? v_bef[i] : -1;
+ float v1 = v_aft ? v_aft[i] : -1;
+ float a0 = a_bef ? a_bef[i] : -1;
+ float a1 = a_aft ? a_aft[i] : -1;
+
+ /* Apply the periodic boundaries */
+ if (periodic) {
+ /* Move towards the origin for more accuracy */
+ const double half = 0.5 * params->box_size[i];
+ if (x0 > x1 + half)
+ x0 -= params->box_size[i];
+ else if (x1 > x0 + half)
+ x1 -= params->box_size[i];
+ }
/* Use quintic hermite spline. */
if (v_bef && v_aft && a_bef && a_aft) {
- x[i] = interpolate_quintic_hermite_spline(t_before, x_bef[i], v_bef[i],
- a_bef[i], t_after, x_aft[i],
- v_aft[i], a_aft[i], t);
+ x[i] = interpolate_quintic_hermite_spline(t_before, x0, v0, a0, t_after,
+ x1, v1, a1, t);
}
/* Use cubic hermite spline. */
else if (v_bef && v_aft) {
- x[i] = interpolate_cubic_hermite_spline(t_before, x_bef[i], v_bef[i],
- t_after, x_aft[i], v_aft[i], t);
+ x[i] = interpolate_cubic_hermite_spline(t_before, x0, v0, t_after, x1, v1,
+ t);
}
/* Use linear interpolation. */
else {
- x[i] = wa * x_aft[i] + wb * x_bef[i];
+ const double wa = (t - t_before) / (t_after - t_before);
+ const double wb = 1. - wa;
+ x[i] = wa * x1 + wb * x0;
+ }
+
+ /* Apply the periodic boundaries to the output */
+ if (periodic) {
+ /* We suppose that only 1 wrapping is enough
+ otherwise something really bad happened */
+ box_wrap(x[i], 0, params->box_size[i]);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check if the periodic box is correctly applied */
+ if (x[i] >= params->box_size[i] || x[i] < 0) {
+ error_python(
+ "A particle is outside the periodic box (%g):"
+ "before=%g, after=%g interpolation=%g",
+ params->box_size[i], x0, x1, x[i]);
+ }
+#endif
}
}
}
@@ -203,11 +243,21 @@ interpolate_quintic_double_float_ND(const double t_before,
* @param t_before Time of field_before (< t).
* @param t_after Time of field_after (> t).
* @param t Requested time.
+ * @param periodic Should the periodic boundary be applied?
+ * @param params The simulation's #logger_parameters.
*/
__attribute__((always_inline)) INLINE static void interpolate_cubic_float_ND(
const double t_before, const struct logger_field *restrict before,
const double t_after, const struct logger_field *restrict after,
- void *restrict output, const double t, const int dimension) {
+ void *restrict output, const double t, const int dimension, int periodic,
+ const struct logger_parameters *params) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* The position is expected to be a double => error with periodic */
+ if (periodic) {
+ error_python("Not implemented yet");
+ }
+#endif
/* Compute the interpolation scaling. */
const float wa = (t - t_before) / (t_after - t_before);
@@ -244,11 +294,21 @@ __attribute__((always_inline)) INLINE static void interpolate_cubic_float_ND(
* @param t_before Time of field_before (< t).
* @param t_after Time of field_after (> t).
* @param t Requested time.
+ * @param periodic Should the periodic boundary be applied?
+ * @param params The simulation's #logger_parameters.
*/
__attribute__((always_inline)) INLINE static void interpolate_linear_float_ND(
const double t_before, const struct logger_field *restrict before,
const double t_after, const struct logger_field *restrict after,
- void *restrict output, const double t, const int dimension) {
+ void *restrict output, const double t, const int dimension, int periodic,
+ const struct logger_parameters *params) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* The position is expected to be a double => error with periodic */
+ if (periodic) {
+ error_python("Not implemented yet");
+ }
+#endif
/* Compute the interpolation scaling. */
const float wa = (t - t_before) / (t_after - t_before);
@@ -274,11 +334,22 @@ __attribute__((always_inline)) INLINE static void interpolate_linear_float_ND(
* @param t_before Time of field_before (< t).
* @param t_after Time of field_after (> t).
* @param t Requested time.
+ * @param periodic Should the periodic boundary be applied?
+ * @param params The simulation's #logger_parameters.
*/
__attribute__((always_inline)) INLINE static void interpolate_linear_double_ND(
const double t_before, const struct logger_field *restrict before,
const double t_after, const struct logger_field *restrict after,
- void *restrict output, const double t, const int dimension) {
+ void *restrict output, const double t, const int dimension, int periodic,
+ const struct logger_parameters *params) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* The velocity and acceleration are supposed to be provided => quintic =>
+ * error here */
+ if (periodic) {
+ error_python("Not implemented yet");
+ }
+#endif
/* Compute the interpolation scaling. */
const double wa = (t - t_before) / (t_after - t_before);
@@ -302,11 +373,21 @@ __attribute__((always_inline)) INLINE static void interpolate_linear_double_ND(
* @param t_before Time of field_before (< t).
* @param t_after Time of field_after (> t).
* @param t Requested time.
+ * @param periodic Should the periodic boundary be applied?
+ * @param params The simulation's #logger_parameters.
*/
__attribute__((always_inline)) INLINE static void interpolate_linear_float(
const double t_before, const struct logger_field *restrict before,
const double t_after, const struct logger_field *restrict after,
- void *restrict output, const double t) {
+ void *restrict output, const double t, int periodic,
+ const struct logger_parameters *params) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* The position is expected to be a double => error with periodic */
+ if (periodic) {
+ error_python("Not implemented yet");
+ }
+#endif
/* Compute the interpolation scaling. */
const float wa = (t - t_before) / (t_after - t_before);
@@ -337,4 +418,5 @@ __attribute__((always_inline)) INLINE static void interpolate_ids(
}
*(long long *)output = *(long long *)after->field;
}
+
#endif // LOGGER_LOGGER_INTERPOLATION_H
diff --git a/logger/logger_parameters.c b/logger/logger_parameters.c
new file mode 100644
index 0000000000000000000000000000000000000000..cb78d7cea44a750a2dc2cfd43faee51b6bd22905
--- /dev/null
+++ b/logger/logger_parameters.c
@@ -0,0 +1,70 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 Loic Hausammann (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 <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* Include corresponding header */
+#include "logger_parameters.h"
+
+/* Include the logger */
+#include "logger_reader.h"
+#include "logger_tools.h"
+
+/* Include swift */
+#include "parser.h"
+
+/**
+ * @brief Initialize the parameters from the yaml file.
+ *
+ * @param reader The #logger_reader.
+ * @param basename The basename of the logger files.
+ * @param verbose The verbose level.
+ */
+void logger_parameters_init(struct logger_parameters *params,
+ const struct logger_reader *reader) {
+
+ /* Generate filename */
+ char filename[STRING_SIZE];
+ strcpy(filename, reader->basename);
+ size_t len = strlen(filename);
+ /* Remove the rank number */
+ filename[len - 5] = '\0';
+ strcat(filename, ".yml");
+
+ /* Initialize the parser */
+ struct swift_params swift_params;
+ parser_read_file(filename, &swift_params);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Print the parameters */
+ if (reader->verbose > 0) parser_print_params(&swift_params);
+#endif
+
+ /* Read the number of dimension */
+ params->dimension = parser_get_param_int(&swift_params, "Header:Dimension");
+
+ /* Read the box size */
+ parser_get_param_double_array(&swift_params, "Header:BoxSize", 3,
+ params->box_size);
+
+ /* Read the periodicity */
+ params->periodic = parser_get_param_int(&swift_params, "Header:Periodic");
+
+ /* Read if we are running with cosmology */
+ params->with_cosmology =
+ parser_get_param_int(&swift_params, "Policy:cosmological integration");
+}
diff --git a/logger/logger_parameters.h b/logger/logger_parameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..8bf8bd2f55a7f85e570b48fb68f9e8ef1095e702
--- /dev/null
+++ b/logger/logger_parameters.h
@@ -0,0 +1,51 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2019 Loic Hausammann (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 <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/**
+ * @file logger_parameters.h
+ * @brief This file contains the functions that deals with the yaml file.
+ */
+
+#ifndef LOGGER_LOGGER_PARAMETERS_H
+#define LOGGER_LOGGER_PARAMETERS_H
+
+struct logger_reader;
+
+/**
+ * @brief Structure containing the simulation's parameters.
+ */
+struct logger_parameters {
+
+ /* Number of dimension */
+ int dimension;
+
+ /* Simulation volume */
+ double box_size[3];
+
+ /* Is the volume periodic? */
+ int periodic;
+
+ /* Are we doing a cosmological simulation? */
+ int with_cosmology;
+};
+
+void logger_parameters_init(struct logger_parameters *params,
+ const struct logger_reader *reader);
+
+#endif // LOGGER_LOGGER_PARAMETERS_H
diff --git a/logger/logger_particle.c b/logger/logger_particle.c
index 6a5084763ea90320f68a2a3ae8332bef00e0185d..e3b02d61e17313eb9e304e5a8d7dbe815045bde5 100644
--- a/logger/logger_particle.c
+++ b/logger/logger_particle.c
@@ -147,12 +147,14 @@ logger_particle_read_special_flag(const struct logger_reader *reader,
* @param time_after Time of field_after (> t).
* @param time Requested time.
* @param field The field to reconstruct.
+ * @param params The simulation's #logger_parameters.
*/
void logger_particle_interpolate_field(
const double time_before, const struct logger_field *restrict before,
const double time_after, const struct logger_field *restrict after,
void *restrict output, const double time,
- const struct field_information *field, enum part_type type) {
+ const struct field_information *field, enum part_type type,
+ const struct logger_parameters *params) {
/* Select the correct interpolation */
switch (type) {
@@ -162,13 +164,14 @@ void logger_particle_interpolate_field(
switch (field->module) {
case field_module_default:
hydro_logger_interpolate_field(time_before, before, time_after, after,
- output, time, field->local_index);
+ output, time, field->local_index,
+ params);
break;
case field_module_chemistry:
- chemistry_logger_interpolate_field_part(time_before, before,
- time_after, after, output,
- time, field->local_index);
+ chemistry_logger_interpolate_field_part(
+ time_before, before, time_after, after, output, time,
+ field->local_index, params);
break;
default:
@@ -183,7 +186,8 @@ void logger_particle_interpolate_field(
if (field->module != field_module_default)
error("Module not implemented");
gravity_logger_interpolate_field(time_before, before, time_after, after,
- output, time, field->local_index);
+ output, time, field->local_index,
+ params);
break;
/* Stars */
@@ -191,19 +195,20 @@ void logger_particle_interpolate_field(
switch (field->module) {
case field_module_default:
stars_logger_interpolate_field(time_before, before, time_after, after,
- output, time, field->local_index);
+ output, time, field->local_index,
+ params);
break;
case field_module_chemistry:
- chemistry_logger_interpolate_field_spart(time_before, before,
- time_after, after, output,
- time, field->local_index);
+ chemistry_logger_interpolate_field_spart(
+ time_before, before, time_after, after, output, time,
+ field->local_index, params);
break;
case field_module_star_formation:
- star_formation_logger_interpolate_field(time_before, before,
- time_after, after, output,
- time, field->local_index);
+ star_formation_logger_interpolate_field(
+ time_before, before, time_after, after, output, time,
+ field->local_index, params);
break;
default:
diff --git a/logger/logger_particle.h b/logger/logger_particle.h
index ce73e9e05a40550d3d624a94bd07d11acc9682d4..c20aabcf82f4d46a9af0d27082dd783d4723bb5d 100644
--- a/logger/logger_particle.h
+++ b/logger/logger_particle.h
@@ -27,6 +27,7 @@
#include "logger_gravity.h"
#include "logger_header.h"
#include "logger_hydro.h"
+#include "logger_parameters.h"
#include "logger_star_formation.h"
#include "logger_stars.h"
#include "logger_time.h"
@@ -53,7 +54,8 @@ void logger_particle_interpolate_field(
const double t_before, const struct logger_field *restrict before,
const double t_after, const struct logger_field *restrict after,
void *restrict output, const double t,
- const struct field_information *field, enum part_type type);
+ const struct field_information *field, enum part_type type,
+ const struct logger_parameters *params);
enum logger_special_flags logger_particle_read_special_flag(
const struct logger_reader *reader, size_t offset, size_t *mask,
diff --git a/logger/logger_python_wrapper.c b/logger/logger_python_wrapper.c
index 6a69b8b95dfcf0bc3afd6b09f941f41941d9b72b..b696022712ef5428f8719c52c8f9b8887d80cb26 100644
--- a/logger/logger_python_wrapper.c
+++ b/logger/logger_python_wrapper.c
@@ -181,6 +181,30 @@ static PyObject *getTimeLimits(PyObject *self, PyObject *Py_UNUSED(ignored)) {
return (PyObject *)out;
}
+/**
+ * @brief Read the box size.
+ *
+ * <b>returns</b> tuple containing the box size.
+ */
+static PyObject *get_box_size(PyObject *self, PyObject *Py_UNUSED(ignored)) {
+ if (!((PyObjectReader *)self)->ready) {
+ error_python(
+ "The logger is not ready yet."
+ "Did you forget to open it with \"with\"?");
+ }
+
+ /* initialize the reader. */
+ struct logger_reader *reader = &((PyObjectReader *)self)->reader;
+
+ /* Create the output */
+ PyObject *out = PyTuple_New(reader->params.dimension);
+ for (int i = 0; i < reader->params.dimension; i++) {
+ PyTuple_SetItem(out, i, PyFloat_FromDouble(reader->params.box_size[i]));
+ }
+
+ return (PyObject *)out;
+}
+
#define find_field_in_module_internal(MODULE, PART) \
for (int local = 0; local < MODULE##_logger_field##PART##_count; local++) { \
const int global = MODULE##_logger_local_to_global##PART[local]; \
@@ -871,6 +895,12 @@ static PyMethodDef libloggerReaderMethods[] = {
"-------\n\n"
"times: tuple\n"
" time min, time max\n"},
+ {"get_box_size", get_box_size, METH_NOARGS,
+ "Gives the box size of the simulation.\n\n"
+ "Returns\n"
+ "-------\n\n"
+ "box_size: tuple\n"
+ " The box size.\n"},
{"get_data", (PyCFunction)pyGetData, METH_VARARGS | METH_KEYWORDS,
"Read some fields from the logfile at a given time.\n\n"
"Parameters\n"
diff --git a/logger/logger_reader.c b/logger/logger_reader.c
index 8875860878730d57199675ec3b61b5912fe5e083..b7713f7e158799b752a68b4dc39f71f8f231d655 100644
--- a/logger/logger_reader.c
+++ b/logger/logger_reader.c
@@ -20,6 +20,8 @@
/* Include corresponding header */
#include "logger_reader.h"
+#include "logger_parameters.h"
+
/* Include standard library */
#include <sys/sysinfo.h>
#include <unistd.h>
@@ -50,6 +52,7 @@ void logger_reader_init(struct logger_reader *reader, const char *basename,
/* Initialize the reader variables. */
reader->verbose = verbose;
+ logger_parameters_init(&reader->params, reader);
/* Generate the logfile filename */
char logfile_name[STRING_SIZE];
@@ -539,7 +542,8 @@ int logger_reader_read_field(struct logger_reader *reader, double time,
/* Interpolate the data. */
logger_particle_interpolate_field(time_before, &before, time_after, &after,
- output, time, field_wanted, type);
+ output, time, field_wanted, type,
+ &reader->params);
}
return 0;
}
diff --git a/logger/logger_reader.h b/logger/logger_reader.h
index d770d8abf41949081559c354302ed391a8969808..8189589d7d47480471a8383163e5c498577300f9 100644
--- a/logger/logger_reader.h
+++ b/logger/logger_reader.h
@@ -50,6 +50,7 @@
#include "logger_index.h"
#include "logger_loader_io.h"
#include "logger_logfile.h"
+#include "logger_parameters.h"
#include "logger_particle.h"
/**
@@ -101,6 +102,9 @@ struct logger_reader {
size_t index;
} time;
+ /* Information from the yaml file */
+ struct logger_parameters params;
+
/* Level of verbosity. */
int verbose;
};
diff --git a/logger/star_formation/GEAR/logger_star_formation.h b/logger/star_formation/GEAR/logger_star_formation.h
index bcdd80e4dd8d16b39cb771cd31e87a663c1d324d..729d914cec8cb6ef70fa229d26e34b2d359c5785 100644
--- a/logger/star_formation/GEAR/logger_star_formation.h
+++ b/logger/star_formation/GEAR/logger_star_formation.h
@@ -24,6 +24,7 @@
/* local includes */
#include "logger_interpolation.h"
#include "logger_loader_io.h"
+#include "logger_parameters.h"
#include "logger_python_tools.h"
#include "star_formation_particle_logger.h"
@@ -58,12 +59,14 @@ star_formation_logger_reader_link_derivatives(struct header *head) {}
* @param t Requested time.
* @param field The field to reconstruct (follows the order of
* #star_formation_logger_fields).
+ * @param params The simulation's #logger_parameters.
*/
__attribute__((always_inline)) INLINE static void
star_formation_logger_interpolate_field(
const double t_before, const struct logger_field *restrict before,
const double t_after, const struct logger_field *restrict after,
- void *restrict output, const double t, const int field) {
+ void *restrict output, const double t, const int field,
+ const struct logger_parameters *params) {
#ifdef SWIFT_DEBUG_CHECKS
/* Check the times */
@@ -78,7 +81,7 @@ star_formation_logger_interpolate_field(
switch (field) {
case star_formation_logger_field_all:
interpolate_linear_float_ND(t_before, before, t_after, after, output, t,
- /* Dimension */ 2);
+ /* Dimension */ 2, /* periodic= */ 0, params);
/* Deal with the progenitor id */
long long *proj_id =
diff --git a/logger/star_formation/none/logger_star_formation.h b/logger/star_formation/none/logger_star_formation.h
index ad5721b94b6cf758b76abcb6df36186590f28920..2c1ed94887f4b9f0abc98cf607edfdd867b4d6c1 100644
--- a/logger/star_formation/none/logger_star_formation.h
+++ b/logger/star_formation/none/logger_star_formation.h
@@ -24,6 +24,7 @@
/* local includes */
#include "logger_interpolation.h"
#include "logger_loader_io.h"
+#include "logger_parameters.h"
#include "logger_python_tools.h"
#include "star_formation_particle_logger.h"
@@ -58,12 +59,14 @@ star_formation_logger_reader_link_derivatives(struct header *head) {}
* @param t Requested time.
* @param field The field to reconstruct (follows the order of
* #star_formation_logger_fields).
+ * @param params The simulation's #logger_parameters.
*/
__attribute__((always_inline)) INLINE static void
star_formation_logger_interpolate_field(
const double t_before, const struct logger_field *restrict before,
const double t_after, const struct logger_field *restrict after,
- void *restrict output, const double t, const int field) {}
+ void *restrict output, const double t, const int field,
+ const struct logger_parameters *params) {}
#ifdef HAVE_PYTHON
__attribute__((always_inline)) INLINE static void
diff --git a/logger/stars/Basic/logger_stars.h b/logger/stars/Basic/logger_stars.h
index 65d6b37dc58b069d1aa44e4be59d7ba34a7489af..a979d74f6406392e4c7e3acc652c2f996e6a3bd3 100644
--- a/logger/stars/Basic/logger_stars.h
+++ b/logger/stars/Basic/logger_stars.h
@@ -24,6 +24,7 @@
/* local includes */
#include "logger_interpolation.h"
#include "logger_loader_io.h"
+#include "logger_parameters.h"
#include "logger_python_tools.h"
#include "stars_logger.h"
@@ -72,6 +73,7 @@ stars_logger_reader_link_derivatives(struct header *head) {
* @param t Requested time.
* @param field The field to reconstruct (follows the order of
* #stars_logger_fields).
+ * @param params The simulation's #logger_parameters.
*/
__attribute__((always_inline)) INLINE static void
stars_logger_interpolate_field(const double t_before,
@@ -79,7 +81,8 @@ stars_logger_interpolate_field(const double t_before,
const double t_after,
const struct logger_field *restrict after,
void *restrict output, const double t,
- const int field) {
+ const int field,
+ const struct logger_parameters *params) {
#ifdef SWIFT_DEBUG_CHECKS
/* Check the times */
@@ -94,19 +97,22 @@ stars_logger_interpolate_field(const double t_before,
switch (field) {
case stars_logger_field_coordinates:
interpolate_quintic_double_float_ND(t_before, before, t_after, after,
- output, t, /* dimension= */ 3);
+ output, t, /* dimension= */ 3,
+ params->periodic, params);
break;
case stars_logger_field_velocities:
interpolate_cubic_float_ND(t_before, before, t_after, after, output, t,
- /* dimension= */ 3);
+ /* dimension= */ 3, /* periodic= */ 0, params);
break;
case stars_logger_field_accelerations:
interpolate_linear_float_ND(t_before, before, t_after, after, output, t,
- /* dimension= */ 3);
+ /* dimension= */ 3, /* periodic= */ 0,
+ params);
break;
case stars_logger_field_smoothing_lengths:
case stars_logger_field_masses:
- interpolate_linear_float(t_before, before, t_after, after, output, t);
+ interpolate_linear_float(t_before, before, t_after, after, output, t,
+ /* periodic= */ 0, params);
break;
case stars_logger_field_particle_ids:
interpolate_ids(t_before, before, t_after, after, output, t);
diff --git a/logger/stars/GEAR/logger_stars.h b/logger/stars/GEAR/logger_stars.h
index 2e7f8ad19253fbe435b7a383f31d700c7d471525..70ad1938ac867d3493c78a15228880effde67846 100644
--- a/logger/stars/GEAR/logger_stars.h
+++ b/logger/stars/GEAR/logger_stars.h
@@ -24,6 +24,7 @@
/* local includes */
#include "logger_interpolation.h"
#include "logger_loader_io.h"
+#include "logger_parameters.h"
#include "logger_python_tools.h"
#include "stars_logger.h"
@@ -72,6 +73,7 @@ stars_logger_reader_link_derivatives(struct header *head) {
* @param t Requested time.
* @param field The field to reconstruct (follows the order of
* #stars_logger_fields).
+ * @param params The simulation's #logger_parameters.
*/
__attribute__((always_inline)) INLINE static void
stars_logger_interpolate_field(const double t_before,
@@ -79,7 +81,8 @@ stars_logger_interpolate_field(const double t_before,
const double t_after,
const struct logger_field *restrict after,
void *restrict output, const double t,
- const int field) {
+ const int field,
+ const struct logger_parameters *params) {
#ifdef SWIFT_DEBUG_CHECKS
/* Check the times */
@@ -94,20 +97,23 @@ stars_logger_interpolate_field(const double t_before,
switch (field) {
case stars_logger_field_coordinates:
interpolate_quintic_double_float_ND(t_before, before, t_after, after,
- output, t, /* dimension= */ 3);
+ output, t, /* dimension= */ 3,
+ params->periodic, params);
break;
case stars_logger_field_velocities:
interpolate_cubic_float_ND(t_before, before, t_after, after, output, t,
- /* dimension= */ 3);
+ /* dimension= */ 3, /* periodic= */ 0, params);
break;
case stars_logger_field_accelerations:
interpolate_linear_float_ND(t_before, before, t_after, after, output, t,
- /* dimension= */ 3);
+ /* dimension= */ 3, /* periodic= */ 0,
+ params);
break;
case stars_logger_field_smoothing_lengths:
case stars_logger_field_masses:
case stars_logger_field_birth_scale_factors:
- interpolate_linear_float(t_before, before, t_after, after, output, t);
+ interpolate_linear_float(t_before, before, t_after, after, output, t,
+ /* periodic= */ 0, params);
break;
case stars_logger_field_particle_ids:
interpolate_ids(t_before, before, t_after, after, output, t);
diff --git a/src/gravity/MultiSoftening/gravity_logger.h b/src/gravity/MultiSoftening/gravity_logger.h
index 0923e6127c52ec141f26c04649a33187cc47ecd5..a29418179adf784e754ecbf06e202c602277fb5e 100644
--- a/src/gravity/MultiSoftening/gravity_logger.h
+++ b/src/gravity/MultiSoftening/gravity_logger.h
@@ -148,8 +148,13 @@ INLINE static char *gravity_logger_write_particle(
/* Write the acceleration. */
if (logger_should_write_field(mask_data[gravity_logger_field_accelerations],
mask)) {
- memcpy(buff, p->a_grav, 3 * sizeof(float));
- buff += 3 * sizeof(float);
+ float acc[3] = {
+ p->a_grav[0] + p->a_grav_mesh[0],
+ p->a_grav[1] + p->a_grav_mesh[1],
+ p->a_grav[2] + p->a_grav_mesh[2],
+ };
+ memcpy(buff, acc, sizeof(acc));
+ buff += sizeof(acc);
}
/* Write the mass. */
diff --git a/src/logger.c b/src/logger.c
index b4dd2901c1f0c6dc640003a82537cb964e05e44d..ff0fc9a467412f433c9504240fdaeade505b75e6 100644
--- a/src/logger.c
+++ b/src/logger.c
@@ -300,6 +300,12 @@ void logger_log_parts(struct logger_writer *log, const struct part *p,
size_t offset_new;
char *buff = (char *)dump_get(&log->dump, size_total, &offset_new);
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Save the buffer position in order to test if the requested buffer was
+ * really used */
+ const char *buff_before = buff;
+#endif
+
/* Write the particles */
for (int i = 0; i < count; i++) {
/* Get the masks */
@@ -343,6 +349,15 @@ void logger_log_parts(struct logger_writer *log, const struct part *p,
xp[i].logger_data.last_offset);
}
}
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Ensure that the buffer was fully used */
+ const int diff = buff - buff_before;
+ if (diff != (int)size_total) {
+ error("It seems that the requested buffer was not totally used: %i != %zi",
+ diff, size_total);
+ }
+#endif
}
/**
@@ -462,6 +477,11 @@ void logger_log_sparts(struct logger_writer *log, struct spart *sp, int count,
/* Allocate a chunk of memory in the dump of the right size. */
size_t offset_new;
char *buff = (char *)dump_get(&log->dump, size_total, &offset_new);
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Save the buffer position in order to test if the requested buffer was
+ * really used */
+ const char *buff_before = buff;
+#endif
for (int i = 0; i < count; i++) {
/* Get the masks */
@@ -509,6 +529,14 @@ void logger_log_sparts(struct logger_writer *log, struct spart *sp, int count,
sp[i].logger_data.last_offset);
}
}
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Ensure that the buffer was fully used */
+ const int diff = buff - buff_before;
+ if (diff != (int)size_total) {
+ error("It seems that the requested buffer was not totally used: %i != %zi",
+ diff, size_total);
+ }
+#endif
}
/**
@@ -595,34 +623,32 @@ void logger_log_gparts(struct logger_writer *log, struct gpart *p, int count,
const uint32_t special_flags = logger_pack_flags_and_data(flag, flag_data);
/* Compute the size of the buffer. */
+ /* As we might have some non DM particles, we cannot log_all_fields blindly */
size_t size_total = 0;
- if (log_all_fields) {
- size_t size = log->max_size_record_gpart + logger_header_bytes;
+ for (int i = 0; i < count; i++) {
+ /* Log only the dark matter */
+ if (p[i].type != swift_type_dark_matter &&
+ p[i].type != swift_type_dark_matter_background)
+ continue;
+
+ unsigned int mask = 0;
+ size_t size = 0;
+ gravity_logger_compute_size_and_mask(log->mask_data_pointers.gravity, &p[i],
+ log_all_fields, &size, &mask);
if (flag != 0) {
size += size_special_flag;
}
- size_total = count * size;
- } else {
- for (int i = 0; i < count; i++) {
- /* Log only the dark matter */
- if (p[i].type != swift_type_dark_matter &&
- p[i].type != swift_type_dark_matter_background)
- continue;
-
- unsigned int mask = 0;
- size_t size = 0;
- gravity_logger_compute_size_and_mask(log->mask_data_pointers.gravity,
- &p[i], log_all_fields, &size, &mask);
- if (flag != 0) {
- size += size_special_flag;
- }
- size_total += size + logger_header_bytes;
- }
+ size_total += size + logger_header_bytes;
}
/* Allocate a chunk of memory in the dump of the right size. */
size_t offset_new;
char *buff = (char *)dump_get(&log->dump, size_total, &offset_new);
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Save the buffer position in order to test if the requested buffer was
+ * really used */
+ const char *buff_before = buff;
+#endif
for (int i = 0; i < count; i++) {
/* Log only the dark matter */
@@ -668,6 +694,14 @@ void logger_log_gparts(struct logger_writer *log, struct gpart *p, int count,
p[i].id_or_neg_offset, p[i].logger_data.last_offset);
}
}
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Ensure that the buffer was fully used */
+ const int diff = buff - buff_before;
+ if (diff != (int)size_total) {
+ error("It seems that the requested buffer was not totally used: %i != %zi",
+ diff, size_total);
+ }
+#endif
}
/**
@@ -946,6 +980,8 @@ void logger_init_masks(struct logger_writer *log, const struct engine *e) {
}
/* Compute the maximal size of the records. */
+
+ /* Hydro */
log->max_size_record_part = 0;
for (int i = 0; i < hydro_logger_field_count; i++) {
log->max_size_record_part += log->mask_data_pointers.hydro[i].size;
@@ -954,11 +990,13 @@ void logger_init_masks(struct logger_writer *log, const struct engine *e) {
log->max_size_record_part += log->mask_data_pointers.chemistry_part[i].size;
}
+ /* Gravity */
log->max_size_record_gpart = 0;
for (int i = 0; i < gravity_logger_field_count; i++) {
log->max_size_record_gpart += log->mask_data_pointers.gravity[i].size;
}
+ /* Stars */
log->max_size_record_spart = 0;
for (int i = 0; i < stars_logger_field_count; i++) {
log->max_size_record_spart += log->mask_data_pointers.stars[i].size;
@@ -976,10 +1014,12 @@ void logger_init_masks(struct logger_writer *log, const struct engine *e) {
log->logger_count_mask = num_fields;
#ifdef SWIFT_DEBUG_CHECKS
- message("The logger contains the following masks:");
- for (int i = 0; i < log->logger_count_mask; i++) {
- message("%20s:\t mask=%03u\t size=%i", log->logger_mask_data[i].name,
- log->logger_mask_data[i].mask, log->logger_mask_data[i].size);
+ if (e->nodeID == 0) {
+ message("The logger contains the following masks:");
+ for (int i = 0; i < log->logger_count_mask; i++) {
+ message("%20s:\t mask=%03u\t size=%i", log->logger_mask_data[i].name,
+ log->logger_mask_data[i].mask, log->logger_mask_data[i].size);
+ }
}
#endif
}
diff --git a/src/logger_io.c b/src/logger_io.c
index 6e444aeae904a6511244fcdc4f55760a044471c8..5f68be24dd605b5c3fd3a4c1a1256e5a91e18414 100644
--- a/src/logger_io.c
+++ b/src/logger_io.c
@@ -417,12 +417,11 @@ void logger_write_index_file(struct logger_writer* log, struct engine* e) {
* @params e The #engine.
*/
void logger_write_description(struct logger_writer* log, struct engine* e) {
+
/* Only the master writes the description */
- if (engine_rank != 0) {
+ if (e->nodeID != 0) {
return;
}
- /* const struct unit_system *internal_units = e->internal_units; */
- /* const struct unit_system *snapshot_units = e->snapshot_units; */
/* File name */
char fileName[FILENAME_BUFFER_SIZE];
@@ -430,13 +429,53 @@ void logger_write_description(struct logger_writer* log, struct engine* e) {
/* Open file */
FILE* f = NULL;
- f = fopen(fileName, "wb");
+ f = fopen(fileName, "w");
if (f == NULL) {
error("Failed to open file %s", fileName);
}
- /* TODO Write stuff */
+ /* Write the header section */
+ fprintf(f, "Header:\n");
+ fprintf(f, " BoxSize: [%g, %g, %g]\n", e->s->dim[0], e->s->dim[1],
+ e->s->dim[2]);
+ fprintf(f, " Dimension: %i\n", (int)hydro_dimension);
+ fprintf(f, " RunName: %s\n", e->run_name);
+ fprintf(f, " Periodic: %i\n", e->s->periodic);
+ fprintf(f, "\n");
+
+ /* Write unit system */
+ const struct unit_system* us = e->internal_units;
+ fprintf(f, "InternalUnitSystem:\n");
+ fprintf(f, " UnitMass_in_cgs: %g\n", units_get_base_unit(us, UNIT_MASS));
+ fprintf(f, " UnitLength_in_cgs: %g\n", units_get_base_unit(us, UNIT_LENGTH));
+ fprintf(f, " UnitTime_in_cgs: %g\n", units_get_base_unit(us, UNIT_TIME));
+ fprintf(f, " UnitCurrent_in_cgs: %g\n",
+ units_get_base_unit(us, UNIT_CURRENT));
+ fprintf(f, " UnitTemp_in_cgs: %g\n",
+ units_get_base_unit(us, UNIT_TEMPERATURE));
+ fprintf(f, "\n");
+
+ /* Write the code section */
+ fprintf(f, "Code:\n");
+ fprintf(f, " Code: SWIFT\n");
+ fprintf(f, " CodeVersion: %s\n", package_version());
+ fprintf(f, " CompilerName: %s\n", compiler_name());
+ fprintf(f, " CompilerVersion: %s\n", compiler_version());
+ fprintf(f, " GitBranch: %s\n", git_branch());
+ fprintf(f, " GitRevision: %s\n", git_revision());
+ fprintf(f, " GitDate: %s\n", git_date());
+ fprintf(f, " ConfigurationOptions: %s\n", configuration_options());
+ fprintf(f, " RandomSeed: %i\n", SWIFT_RANDOM_SEED_XOR);
+ fprintf(f, "\n");
+
+ /* Write the policy section */
+ fprintf(f, "Policy:\n");
+ for (int i = 1; i < engine_maxpolicy; i++) {
+ const int value = e->policy & (1 << i) ? 1 : 0;
+ fprintf(f, " %s: %i\n", engine_policy_names[i + 1], value);
+ }
+ fprintf(f, "\n");
/* Close file */
fclose(f);
diff --git a/tools/read_index_file.py b/tools/read_index_file.py
index 51fdec59208e10e62ba5293f23e72674866cffe1..944d6aa4306bfe22772440823afe66af0225ffe3 100644
--- a/tools/read_index_file.py
+++ b/tools/read_index_file.py
@@ -2,64 +2,231 @@
import sys
import numpy as np
-
-filename = sys.argv[-1]
-n_type = 6
-
-
-# dtype for the particle's data
-dt = np.dtype([("ids", np.ulonglong),
- ("offset", np.uint64)])
-
-# Read the file
-with open(filename, "rb") as f:
- # read the time
- time = np.fromfile(f, dtype=float, count=1)
- time_int = np.fromfile(f, dtype=np.longlong, count=1)
- print("Time: {}, integer time: {}".format(
- time[0], time_int[0]))
-
- # read the number of particles
- nparts = np.fromfile(f, dtype=np.uint64, count=n_type)
-
- print("Number of particles:", nparts)
-
- # read if the file is sorted
- sort = np.fromfile(f, dtype=np.bool, count=1)
- print("File is sorted?", sort[0])
-
- # read the memory alignment garbage
- n = ((f.tell() + 7) & ~7) - f.tell()
- f.read(n)
-
- # read the particles
- print("Particles data (ids / offset):")
- for n in nparts:
- if n == 0:
- continue
-
- data = np.fromfile(f, dtype=dt, count=n)
-
- print("\t", data)
-
- # print the history of new particles
- n_new = np.fromfile(f, dtype=np.uint64, count=n_type)
- print("New particles: ", n_new)
-
- for n in n_new:
- if n == 0:
- continue
-
- data = np.fromfile(f, dtype=dt, count=n)
- print("\t", data)
-
- # print the history of particles removed
- n_rem = np.fromfile(f, dtype=np.uint64, count=n_type)
- print("Particles removed: ", n_rem)
-
- for n in n_rem:
- if n == 0:
- continue
-
- data = np.fromfile(f, dtype=dt, count=n)
- print("\t", data)
+import mmap
+import struct
+import os
+
+
+class DarkMatterParticle():
+ def __init__(self, pos, vel, acc, mass, ids):
+ self.pos = pos
+ self.vel = vel
+ self.acc = acc
+ self.mass = mass
+ self.ids = ids
+
+class LogfileReader():
+ n_type = 6
+ swift_type_dark_matter = 1
+
+ # dtype for the particle's data
+ index_dtype = np.dtype([("ids", np.ulonglong),
+ ("offset", np.uint64)])
+
+ def __init__(self, basename):
+ self.basename = basename
+ logfile = basename + ".dump"
+ self.f = open(logfile, "rb")
+ self.m = mmap.mmap(self.f.fileno(), 0, prot=mmap.PROT_READ)
+ self.read_logfile_header()
+
+ def __del__(self):
+ self.m.close()
+ self.f.close()
+
+ def read_int(self):
+ val = self.m.read(4)
+ return int.from_bytes(val, byteorder="little")
+
+ def read_uint(self):
+ val = self.m.read(4)
+ return struct.unpack("I", val)[0]
+
+ def read_offset(self):
+ offset = self.m.read(6)
+ return int.from_bytes(offset, byteorder="little")
+
+ def read_string(self):
+ string = self.m.read(self.string_length)
+ string = struct.unpack("%is" % self.string_length, string)[0]
+ i = string.find(b"\x00")
+ return string[:i]
+
+ def read_logfile_header(self):
+ self.m.seek(0)
+ self.vers = self.m.read(20)
+ self.major = self.read_int()
+ self.minor = self.read_int()
+ self.offset_direction = self.read_int()
+ self.first_record = self.read_offset()
+ self.string_length = self.read_uint()
+ self.n_masks = self.read_uint()
+
+ mask_dtype = np.dtype([("name", "U%i" % self.string_length),
+ ("value", "int32"), ("size", "uint32")])
+
+ self.masks = np.empty(self.n_masks, dtype=mask_dtype)
+ for i in range(self.n_masks):
+ self.masks[i]["name"] = self.read_string()
+ self.masks[i]["value"] = 1 << i
+ self.masks[i]["size"] = self.read_uint()
+
+ print("Fields found (name, mask, size): ", self.masks)
+
+
+ def read_record_header(self):
+ mask = self.m.read(2)
+ mask = struct.unpack("h", mask)[0]
+
+ if (mask != 316):
+ print(mask)
+
+ offset = self.read_offset()
+ return mask, offset
+
+ def read_index(self, index_number):
+ filename = basename + "_%04i.index" % index_number
+
+ # Read the file
+ with open(filename, "rb") as f:
+ # read the time
+ time = np.fromfile(f, dtype=float, count=1)
+ time_int = np.fromfile(f, dtype=np.longlong, count=1)
+ print("Time: {}, integer time: {}".format(
+ time[0], time_int[0]))
+
+ # read the number of particles
+ nparts = np.fromfile(f, dtype=np.uint64, count=LogfileReader.n_type)
+
+ print("Number of particles:", nparts)
+
+ # read if the file is sorted
+ sort = np.fromfile(f, dtype=np.bool, count=1)
+ print("File is sorted?", sort[0])
+
+ # read the memory alignment garbage
+ n = ((f.tell() + 7) & ~7) - f.tell()
+ f.read(n)
+
+ # read the particles
+ print("Particles data (ids / offset):")
+ data = []
+ for n in nparts:
+ if n == 0:
+ data.append(None)
+ continue
+
+ tmp = np.fromfile(f, dtype=LogfileReader.index_dtype, count=n)
+ data.append(tmp)
+
+ print("\t", data)
+
+ # print the history of new particles
+ n_new = np.fromfile(f, dtype=np.uint64, count=LogfileReader.n_type)
+ print("New particles: ", n_new)
+
+ new = []
+ for n in n_new:
+ if n == 0:
+ new.append(None)
+ continue
+
+ new.append(np.fromfile(f, dtype=LogfileReader.index_dtype, count=n))
+ print("\t", new)
+
+ # print the history of particles removed
+ n_rem = np.fromfile(f, dtype=np.uint64, count=LogfileReader.n_type)
+ print("Particles removed: ", n_rem)
+
+ removed = []
+ for n in n_rem:
+ if n == 0:
+ removed.append(None)
+ continue
+
+ removed.append(np.fromfile(f, dtype=LogfileReader.index_dtype, count=n))
+ print("\t", removed)
+ return data, new, removed
+
+ def seek(self, new_offset):
+ self.m.seek(new_offset)
+
+ def read_particle(self, mask):
+ flag = None
+ ref = self.masks["name"] == "SpecialFlags"
+ if (mask & self.masks[ref]["value"]):
+ flag = self.m.read(4)
+ print("Flag: ", struct.unpack("i", flag)[0])
+
+ pos = None
+ ref = self.masks["name"] == "Coordinates"
+ if (mask & self.masks[ref]["value"]):
+ pos = self.m.read(24)
+ pos = struct.unpack("ddd", pos)[0]
+
+ vel = None
+ ref = self.masks["name"] == "Velocities"
+ if (mask & self.masks[ref]["value"]):
+ vel = self.m.read(12)
+ vel = struct.unpack("fff", vel)[0]
+
+ acc = None
+ ref = self.masks["name"] == "Accelerations"
+ if (mask & self.masks[ref]["value"]):
+ acc = self.m.read(12)
+ acc = struct.unpack("fff", acc)[0]
+
+ # print(pos, vel, acc)
+ mass = None
+ ref = self.masks["name"] == "Masses"
+ if (mask & self.masks[ref]["value"]):
+ mass = self.m.read(4)
+ mass = struct.unpack("f", mass)[0]
+
+ ids = None
+ ref = self.masks["name"] == "ParticleIDs"
+ if (mask & self.masks[ref]["value"]):
+ ids = self.m.read(8)
+ ids = struct.unpack("q", ids)[0]
+
+ part = DarkMatterParticle(pos, vel, acc, mass, ids)
+ return part
+
+
+if __name__ == "__main__":
+ logfile = sys.argv[-1]
+ basename = logfile[:-5]
+
+ # Open the logfile
+ reader = LogfileReader(basename)
+
+ # Read the first index file
+ data, _, _ = reader.read_index(0)
+
+ # Keep only the DM particles
+ data = data[LogfileReader.swift_type_dark_matter]
+
+ for i in range(len(data)):
+ if i > 10:
+ print("Early exit due to low performances")
+ break
+
+ off = data["offset"][i]
+ off = int(off)
+ index_ids = data["ids"][i]
+ reader.seek(off)
+
+ for i in range(10):
+ mask, offset = reader.read_record_header()
+ print(mask, offset, off)
+
+ part = reader.read_particle(mask)
+ if (index_ids != part.ids):
+ raise Exception("Failed", index_ids, part.ids)
+
+ if offset == 0:
+ print("End of particle")
+ break
+
+ off += offset
+ reader.seek(off)