diff --git a/configure.ac b/configure.ac
index 8202e6c405e97dde2e8a788b2ecfebde0ddc378f..d119cfa3a99c1c85456056eab5d73b334ec86d15 100644
--- a/configure.ac
+++ b/configure.ac
@@ -2382,6 +2382,12 @@ AM_CONDITIONAL([HAVEGRACKLECOOLING], [test "${with_cooling:0:7}" = "grackle"])
# check if using gear feedback
AM_CONDITIONAL([HAVEGEARFEEDBACK], [test "$with_feedback" = "GEAR"])
+# check if using SPHENIX
+AM_CONDITIONAL([HAVE_SPHENIX], [test "$with_hydro" = "sphenix"])
+
+# check if using GADGET2 SPH
+AM_CONDITIONAL([HAVE_GADGET2], [test "$with_hydro" = "gadget2"])
+
# Handle .in files.
AC_CONFIG_FILES([Makefile src/Makefile examples/Makefile examples/Cooling/CoolingRates/Makefile doc/Makefile doc/Doxyfile tests/Makefile])
AC_CONFIG_FILES([argparse/Makefile tools/Makefile logger/Makefile logger/tests/Makefile])
diff --git a/examples/main.c b/examples/main.c
index 9f59e9d955f955ff0d8fe96b21c24d382e269cbd..1b9c21bdb1e9f9effccd29ece728d7bfb60df3bf 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -1638,8 +1638,13 @@ int main(int argc, char *argv[]) {
engine_dump_index(&e);
/* Write a sentinel timestamp */
- logger_log_timestamp(e.logger, e.ti_current, e.time,
- &e.logger->timestamp_offset);
+ if (e.policy & engine_policy_cosmology) {
+ logger_log_timestamp(e.logger, e.ti_current, e.cosmology->a,
+ &e.logger->timestamp_offset);
+ } else {
+ logger_log_timestamp(e.logger, e.ti_current, e.time,
+ &e.logger->timestamp_offset);
+ }
}
#endif
diff --git a/logger/Makefile.am b/logger/Makefile.am
index af9c040e9c8fb5dbaf01cb8acf8ec671a2bac8f6..441aaab2a9be434b80b7903e290ff3535f1c0cf6 100644
--- a/logger/Makefile.am
+++ b/logger/Makefile.am
@@ -45,14 +45,20 @@ lib_LTLIBRARIES = liblogger.la
GRAVITY_SRC = gravity/MultiSoftening/logger_gravity.c
STARS_SRC = stars/Default/logger_stars.c
+
+if HAVE_SPHENIX
+HYDRO_SRC = hydro/SPHENIX/logger_hydro.c
+endif
+if HAVE_GADGET2
HYDRO_SRC = hydro/Gadget2/logger_hydro.c
+endif
# subdirectories
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_logfile.h logger_index.h quick_sort.h logger_python_tools.h logger_interpolation.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 \
@@ -70,7 +76,8 @@ endif
# Include files for distribution, not installation.
nobase_noinst_HEADERS = logger_hydro.h hydro/Gadget2/logger_hydro.h \
logger_stars.h stars/Default/logger_stars.h \
- logger_gravity.h gravity/MultiSoftening/logger_gravity.h
+ logger_gravity.h gravity/MultiSoftening/logger_gravity.h \
+ hydro/SPHENIX/logger_hydro.h
# Sources and flags for regular library
liblogger_la_SOURCES = $(AM_SOURCES)
diff --git a/logger/examples/reader_example.py b/logger/examples/reader_example.py
index 6fe65399cc70de43b6638820c19f41486ee2d9df..e3518cf61914ba6916075433c90b1f0877f0e4f8 100644
--- a/logger/examples/reader_example.py
+++ b/logger/examples/reader_example.py
@@ -6,12 +6,31 @@ Example: ./reader_example.py ../../examples/SedovBlast_3D/index 0.1
import sys
import numpy as np
import matplotlib.pyplot as plt
+from glob import glob
from mpl_toolkits.mplot3d import Axes3D
+import argparse
sys.path.append("../.libs/")
import liblogger as logger
+def parse_arguments():
+ parser = argparse.ArgumentParser(
+ description='Read a logfile and plots some basic properties')
+
+ default_files = "../../examples/HydroTests/SedovBlast_3D/index_*dump"
+ default_files = glob(default_files)
+
+ parser.add_argument("-t", '--time', dest='time',
+ type=float, default=0.01,
+ help='Simulation time')
+ parser.add_argument('files', metavar='filenames', type=str, nargs="*",
+ help='The filenames of the logfiles')
+ args = parser.parse_args()
+ if len(args.files) == 0:
+ args.files = default_files
+ return args
+
def plot3D(pos):
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
@@ -31,36 +50,34 @@ def plot2D(pos, entropies):
plt.ylabel("Entropy")
-basename = "../../examples/HydroTests/SedovBlast_3D/index_0000"
-time = 0.03
-if len(sys.argv) >= 2:
- basename = sys.argv[1]
-else:
- print("No basename supplied (first argument), using default.")
-if len(sys.argv) >= 3:
- time = float(sys.argv[2])
-else:
- print("No time supplied (second argument), using default.")
-if len(sys.argv) > 3:
- print("Ignoring excess arguments '%s'." % sys.argv[3:])
-print("basename: %s" % basename)
-print("time: %g" % time)
+args = parse_arguments()
+print("basename: %s" % args.files)
+print("time: %g" % args.time)
# read the logger
-gas_type = 0
-with logger.Reader(basename, verbose=0) as reader:
- t = reader.get_time_limits()
-
- fields = reader.get_list_fields(gas_type)
- if ("Coordinates" not in fields or
- "Entropies" not in fields):
- raise Exception("Field not found in the logfile")
-
- pos, ent = reader.get_particle_data(
- ["Coordinates", "Entropies"], time, gas_type)
-
-print("Min/Max of the position:", pos.min(), pos.max())
-print("Min/Max of the entropy:", ent.min(), ent.max())
-plot3D(pos)
-plot2D(pos, ent)
+positions = np.empty((0, 3))
+entropies = np.empty(0)
+for f in args.files:
+ if f.endswith(".dump"):
+ filename = f[:-5]
+ else:
+ raise Exception("It seems that you are not providing a logfile (.dump)")
+ with logger.Reader(filename, verbose=0) as reader:
+ t = reader.get_time_limits()
+ out = reader.get_particle_data(
+ ["Coordinates", "Entropies"], args.time)
+
+ fields = reader.get_list_fields(gas_type)
+ if ("Coordinates" not in fields or
+ "Entropies" not in fields):
+ raise Exception("Field not found in the logfile")
+
+ # add the data to the list
+ positions = np.append(positions, out[0], axis=0)
+ entropies = np.append(entropies, out[1], axis=0)
+
+print("Min/Max of the position:", positions.min(), positions.max())
+print("Min/Max of the entropy:", entropies.min(), entropies.max())
+plot3D(positions)
+plot2D(positions, entropies)
plt.show()
diff --git a/logger/gravity/MultiSoftening/logger_gravity.h b/logger/gravity/MultiSoftening/logger_gravity.h
index e89acb9afadda60a8c4d42620c0b67ad09340dd0..c5cfcdc07e476f8c04191fe139b446ed281fd107 100644
--- a/logger/gravity/MultiSoftening/logger_gravity.h
+++ b/logger/gravity/MultiSoftening/logger_gravity.h
@@ -23,6 +23,7 @@
/* local includes */
#include "gravity_logger.h"
+#include "logger_interpolation.h"
#include "logger_loader_io.h"
#include "logger_python_tools.h"
@@ -90,77 +91,24 @@ gravity_logger_interpolate_field(const double t_before,
}
#endif
- /* Compute the interpolation scaling. */
- const double wa = (t - t_before) / (t_after - t_before);
- const double wb = 1. - wa;
-
switch (field) {
case gravity_logger_field_coordinates:
- for (int i = 0; i < 3; 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;
-
- const double *x_aft = (double *)after->field;
- const float *v_aft = (float *)after->first_deriv;
- const float *a_aft = (float *)after->second_deriv;
-
- /* Use quintic hermite spline. */
- if (v_bef && v_aft && a_bef && a_aft) {
- x[i] = logger_tools_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);
- }
- /* Use cubic hermite spline. */
- else if (v_bef && v_aft) {
- x[i] = logger_tools_cubic_hermite_spline(
- t_before, x_bef[i], v_bef[i], t_after, x_aft[i], v_aft[i], t);
- }
- /* Use linear interpolation. */
- else {
- x[i] = wa * x_aft[i] + wb * x_bef[i];
- }
- }
+ interpolate_quintic_double_float_ND(t_before, before, t_after, after,
+ output, t, /* dimension= */ 3);
break;
case gravity_logger_field_velocities:
- for (int i = 0; i < 3; i++) {
- float *v = (float *)output;
- const float *v_bef = (float *)before->field;
- const float *a_bef = (float *)before->first_deriv;
-
- const float *v_aft = (float *)after->field;
- const float *a_aft = (float *)after->first_deriv;
-
- /* Use a cubic hermite spline. */
- if (a_bef && a_aft) {
- v[i] = logger_tools_cubic_hermite_spline(
- t_before, v_bef[i], a_bef[i], t_after, v_aft[i], a_aft[i], t);
- }
- /* Use linear interpolation. */
- else {
- v[i] = wa * v_aft[i] + wb * v_bef[i];
- }
- }
+ interpolate_cubic_float_ND(t_before, before, t_after, after, output, t,
+ /* dimension= */ 3);
break;
case gravity_logger_field_accelerations:
- /* interpolate vectors. */
- for (int i = 0; i < 3; i++) {
- float *a = (float *)output;
- const float *a_bef = (float *)before->field;
- const float *a_aft = (float *)after->field;
- a[i] = wa * a_aft[i] + wb * a_bef[i];
- }
+ interpolate_linear_float_ND(t_before, before, t_after, after, output, t,
+ /* dimension= */ 3);
break;
case gravity_logger_field_masses:
- ((float *)output)[0] =
- wa * ((float *)after->field)[0] + wb * ((float *)before->field)[0];
+ interpolate_linear_float(t_before, before, t_after, after, output, t);
break;
case gravity_logger_field_particle_ids:
- if (*(long long *)after->field != *(long long *)before->field) {
- error("Interpolating different particles");
- }
- *(long long *)output = *(long long *)after->field;
+ interpolate_ids(t_before, before, t_after, after, output, t);
break;
default:
error("Not implemented");
diff --git a/logger/hydro/Gadget2/logger_hydro.c b/logger/hydro/Gadget2/logger_hydro.c
index 57c8cb35c9d46a14bc296a764fe7520abf173e46..9f252c4152179da33dc4c477db54e8c81bc00e41 100644
--- a/logger/hydro/Gadget2/logger_hydro.c
+++ b/logger/hydro/Gadget2/logger_hydro.c
@@ -19,6 +19,8 @@
#include "logger_hydro.h"
+#include "config.h"
+
int hydro_logger_local_to_global[hydro_logger_field_count];
/* Define the size of all the fields. */
#define member_size(type, member) sizeof(((type *)0)->member)
diff --git a/logger/hydro/Gadget2/logger_hydro.h b/logger/hydro/Gadget2/logger_hydro.h
index 96431e35735c88186c6aad1692f675f03814989c..affa61e8ec7bb6491f99ffdb134a0fa3bd0bf84e 100644
--- a/logger/hydro/Gadget2/logger_hydro.h
+++ b/logger/hydro/Gadget2/logger_hydro.h
@@ -23,6 +23,7 @@
/* local includes */
#include "hydro_logger.h"
+#include "logger_interpolation.h"
#include "logger_loader_io.h"
#include "logger_python_tools.h"
@@ -90,84 +91,31 @@ hydro_logger_interpolate_field(const double t_before,
}
#endif
- /* Compute the interpolation scaling. */
- const double wa = (t - t_before) / (t_after - t_before);
- const double wb = 1. - wa;
-
switch (field) {
/* Do the position */
case hydro_logger_field_coordinates:
- for (int i = 0; i < 3; 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;
-
- const double *x_aft = (double *)after->field;
- const float *v_aft = (float *)after->first_deriv;
- const float *a_aft = (float *)after->second_deriv;
-
- /* Use quintic hermite spline. */
- if (v_bef && v_aft && a_bef && a_aft) {
- x[i] = logger_tools_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);
- }
- /* Use cubic hermite spline. */
- else if (v_bef && v_aft) {
- x[i] = logger_tools_cubic_hermite_spline(
- t_before, x_bef[i], v_bef[i], t_after, x_aft[i], v_aft[i], t);
- }
- /* Use linear interpolation. */
- else {
- x[i] = wa * x_aft[i] + wb * x_bef[i];
- }
- }
+ interpolate_quintic_double_float_ND(t_before, before, t_after, after,
+ output, t, /* dimesion= */ 3);
break;
/* Do the velocity */
case hydro_logger_field_velocities:
- for (int i = 0; i < 3; i++) {
- float *v = (float *)output;
- const float *v_bef = (float *)before->field;
- const float *a_bef = (float *)before->first_deriv;
-
- const float *v_aft = (float *)after->field;
- const float *a_aft = (float *)after->first_deriv;
-
- /* Use a cubic hermite spline. */
- if (a_bef && a_aft) {
- v[i] = logger_tools_cubic_hermite_spline(
- t_before, v_bef[i], a_bef[i], t_after, v_aft[i], a_aft[i], t);
- }
- /* Use linear interpolation. */
- else {
- v[i] = wa * v_aft[i] + wb * v_bef[i];
- }
- }
+ interpolate_cubic_float_ND(t_before, before, t_after, after, output, t,
+ /* dimesion= */ 3);
break;
case hydro_logger_field_accelerations:
- /* interpolate vectors. */
- for (int i = 0; i < 3; i++) {
- float *a = (float *)output;
- const float *a_bef = (float *)before->field;
- const float *a_aft = (float *)after->field;
- a[i] = wa * a_aft[i] + wb * a_bef[i];
- }
+ interpolate_linear_float_ND(t_before, before, t_after, after, output, t,
+ /* dimesion= */ 3);
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:
- ((float *)output)[0] =
- wa * ((float *)after->field)[0] + wb * ((float *)before->field)[0];
+ interpolate_linear_float(t_before, before, t_after, after, output, t);
break;
/* Check the ids */
case hydro_logger_field_particle_ids:
- if (*(long long *)after->field != *(long long *)before->field) {
- error("Interpolating different particles");
- }
- *(long long *)output = *(long long *)after->field;
+ interpolate_ids(t_before, before, t_after, after, output, t);
break;
default:
error("Not implemented");
diff --git a/logger/hydro/SPHENIX/logger_hydro.c b/logger/hydro/SPHENIX/logger_hydro.c
new file mode 100644
index 0000000000000000000000000000000000000000..d129cb51c06e12f2a59214c2c47ddb8b9589b6fc
--- /dev/null
+++ b/logger/hydro/SPHENIX/logger_hydro.c
@@ -0,0 +1,42 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 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 "logger_hydro.h"
+
+#include "config.h"
+
+int hydro_logger_local_to_global[hydro_logger_field_count];
+
+/* Define the size of all the fields. */
+#define member_size(type, member) sizeof(((type *)0)->member)
+
+const int hydro_logger_field_size[hydro_logger_field_count] = {
+ member_size(struct part, x), // coordinates
+ member_size(struct part, v), // velocities
+ member_size(struct part, a_hydro), // accelerations
+ member_size(struct part, mass), // massses
+ member_size(struct part, h), // Smoothing Length
+ member_size(struct part, u), // InternalEnergies
+ member_size(struct part, id), // IDs
+ member_size(struct part, rho), // density
+ sizeof(float), // Entropy
+ sizeof(float), // Pressure
+ 3 * sizeof(float), // Viscosity / diffusion
+ 2 * sizeof(float), // Velocity divergence + deriv
+};
diff --git a/logger/hydro/SPHENIX/logger_hydro.h b/logger/hydro/SPHENIX/logger_hydro.h
new file mode 100644
index 0000000000000000000000000000000000000000..0565db40a5a6a681460bb2641cb2d540864965cb
--- /dev/null
+++ b/logger/hydro/SPHENIX/logger_hydro.h
@@ -0,0 +1,177 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_SPHENIX_LOGGER_HYDRO_H
+#define SWIFT_SPHENIX_LOGGER_HYDRO_H
+
+#include "../config.h"
+
+/* local includes */
+#include "hydro.h"
+#include "hydro_logger.h"
+#include "logger_interpolation.h"
+#include "logger_loader_io.h"
+#include "logger_python_tools.h"
+
+/* Index of the mask in the header mask array */
+extern int hydro_logger_local_to_global[hydro_logger_field_count];
+
+/* Size for each mask */
+extern const int hydro_logger_field_size[hydro_logger_field_count];
+
+/**
+ * @brief Create the link between the fields and their derivatives.
+ *
+ * @param head The #header.
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_logger_reader_link_derivatives(struct header *head) {
+
+ /* Set the first and second derivatives */
+ const int pos_id =
+ hydro_logger_local_to_global[hydro_logger_field_coordinates];
+ const int vel_id =
+ hydro_logger_local_to_global[hydro_logger_field_velocities];
+ const int acc_id =
+ hydro_logger_local_to_global[hydro_logger_field_accelerations];
+
+ /* Coordinates */
+ header_set_first_derivative(head, pos_id, vel_id);
+ header_set_second_derivative(head, pos_id, acc_id);
+
+ /* Velocities */
+ header_set_first_derivative(head, vel_id, acc_id);
+}
+
+/**
+ * @brief Interpolate a field of the particle at the given time.
+ * Here we use a linear interpolation for most of the fields.
+ * For the position (velocity), we use a quintic (cubic) hermite interpolation
+ * based on the positions, velocities and accelerations at the time of the two
+ * particles.
+ *
+ * @param before Pointer to the #logger_field at a time < t.
+ * @param after Pointer to the #logger_field at a time > t.
+ * @param otuput Pointer to the output value.
+ * @param t_before Time of field_before (< t).
+ * @param t_after Time of field_after (> t).
+ * @param t Requested time.
+ * @param field The field to reconstruct (follows the order of
+ * #hydro_logger_fields).
+ */
+__attribute__((always_inline)) INLINE static void
+hydro_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) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check the times */
+ if (t_before > t || t_after < t) {
+ error(
+ "The times for the interpolation are not "
+ "monotonically increasing %g < %g < %g.",
+ t_before, t, t_after);
+ }
+#endif
+
+ switch (field) {
+ /* Do the position */
+ case hydro_logger_field_coordinates:
+ interpolate_quintic_double_float_ND(t_before, before, t_after, after,
+ output, t, /* dimension= */ 3);
+ break;
+ /* Do the velocity */
+ case hydro_logger_field_velocities:
+ interpolate_cubic_float_ND(t_before, before, t_after, after, output, t,
+ /* dimension= */ 3);
+ break;
+ case hydro_logger_field_accelerations:
+ case hydro_logger_field_viscosity_diffusion:
+ interpolate_linear_float_ND(t_before, before, t_after, after, output, t,
+ /* dimension= */ 3);
+ 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:
+ case hydro_logger_field_entropies:
+ case hydro_logger_field_pressures:
+ interpolate_linear_float(t_before, before, t_after, after, output, t);
+ break;
+ /* Check the ids */
+ case hydro_logger_field_particle_ids:
+ interpolate_ids(t_before, before, t_after, after, output, t);
+ break;
+
+ case hydro_logger_field_velocity_divergences: {
+ /* Get some variables */
+ const float wa = (t - t_before) / (t_after - t_before);
+ const float wb = 1. - wa;
+ float *x = (float *)output;
+ const float *div_bef = (float *)before->field;
+ const float *div_aft = (float *)after->field;
+
+ /* Use cubic hermite spline. */
+ x[0] = interpolate_cubic_hermite_spline(
+ t_before, div_bef[0], div_bef[1], t_after, div_aft[0], div_aft[1], t);
+ /* Use the linear interpolation */
+ x[1] = wa * div_aft[1] + wb * div_bef[1];
+ break;
+ }
+
+ default:
+ error("Not implemented");
+ }
+}
+
+#ifdef HAVE_PYTHON
+__attribute__((always_inline)) INLINE static void hydro_logger_generate_python(
+ struct logger_python_field *fields) {
+
+ fields[hydro_logger_field_coordinates] =
+ logger_loader_python_field(/* Dimension */ 3, NPY_DOUBLE);
+ fields[hydro_logger_field_velocities] =
+ logger_loader_python_field(/* Dimension */ 3, NPY_FLOAT32);
+ fields[hydro_logger_field_accelerations] =
+ logger_loader_python_field(/* Dimension */ 3, NPY_FLOAT32);
+ fields[hydro_logger_field_masses] =
+ logger_loader_python_field(/* Dimension */ 1, NPY_FLOAT32);
+ fields[hydro_logger_field_smoothing_lengths] =
+ logger_loader_python_field(/* Dimension */ 1, NPY_FLOAT32);
+ fields[hydro_logger_field_internal_energies] =
+ logger_loader_python_field(/* Dimension */ 1, NPY_FLOAT32);
+ fields[hydro_logger_field_particle_ids] =
+ logger_loader_python_field(/* Dimension */ 1, NPY_LONGLONG);
+ fields[hydro_logger_field_densities] =
+ logger_loader_python_field(/* Dimension */ 1, NPY_FLOAT32);
+ fields[hydro_logger_field_entropies] =
+ logger_loader_python_field(/* Dimension */ 1, NPY_FLOAT32);
+ fields[hydro_logger_field_pressures] =
+ logger_loader_python_field(/* Dimension */ 1, NPY_FLOAT32);
+ fields[hydro_logger_field_viscosity_diffusion] =
+ logger_loader_python_field(/* Dimension */ 3, NPY_FLOAT32);
+ fields[hydro_logger_field_velocity_divergences] =
+ logger_loader_python_field(/* Dimension */ 2, NPY_FLOAT32);
+}
+
+#endif // HAVE_PYTHON
+#endif /* SWIFT_SPHENIX_LOGGER_HYDRO_H */
diff --git a/logger/logger_header.c b/logger/logger_header.c
index ec0d1453df757e4f128d637c9f08b8d0e7438f9b..52bc8132f7f0e218299139146ec01223ee024c5d 100644
--- a/logger/logger_header.c
+++ b/logger/logger_header.c
@@ -171,6 +171,11 @@ void header_read(struct header *h, struct logger_logfile *log) {
map = logger_loader_io_read_data(map, sizeof(unsigned int),
&h->masks[i].size);
+ /* Print the information. */
+ if (reader->verbose > 1) {
+ message("Field %s found in the logfile", h->masks[i].name);
+ }
+
/* Keep the timestamp mask in memory */
if (strcmp(h->masks[i].name, "Timestamp") == 0) {
h->timestamp_mask = h->masks[i].mask;
@@ -182,6 +187,16 @@ void header_read(struct header *h, struct logger_logfile *log) {
error_python("Unable to find the timestamp mask.");
}
+ /* check the special flag. */
+ if (strcmp(h->masks[logger_index_special_flags].name, "SpecialFlags") != 0) {
+ error("The special flag is expected to be at position %i",
+ logger_index_special_flags);
+ }
+
+ if (h->masks[logger_index_special_flags].size != sizeof(uint32_t)) {
+ error("The special flag is expected to be 32 bits");
+ }
+
/* Check the logfile header's size. */
if (map != (char *)log->log.map + h->offset_first_record) {
size_t offset = (char *)map - (char *)log->log.map;
diff --git a/logger/logger_hydro.h b/logger/logger_hydro.h
index f9729c55c4bc1e7b6b85e0977a517a23d80305d3..c012ad5962482c18cffa48d0efbd89a3e8c7710e 100644
--- a/logger/logger_hydro.h
+++ b/logger/logger_hydro.h
@@ -50,7 +50,6 @@
#error TODO
#include "./hydro/Planetary/logger_hydro.h"
#elif defined(SPHENIX_SPH)
-#error TODO
#include "./hydro/SPHENIX/logger_hydro.h"
#elif defined(ANARCHY_PU_SPH)
#error TODO
diff --git a/logger/logger_index.c b/logger/logger_index.c
index a62ff168bea2e21af4bd0998c37c186d020bd49c..50a8545b64a943bbd2446c7a3e16a9bf71abbf7f 100644
--- a/logger/logger_index.c
+++ b/logger/logger_index.c
@@ -72,15 +72,40 @@ void logger_index_read_header(struct logger_index *index,
(char *)index->index.map + logger_index_integer_time_offset,
logger_index_integer_time_size);
- /* Read the number of particles */
+ /* Read the number of particles. */
memcpy(index->nparts, (char *)index->index.map + logger_index_npart_offset,
logger_index_npart_size);
- /* Read if the file is sorted */
+ /* Read whether the file is sorted. */
memcpy(&index->is_sorted,
(char *)index->index.map + logger_index_is_sorted_offset,
logger_index_is_sorted_size);
+ /* Compute the position of the history of new particles. */
+ size_t count = 0;
+ for (int i = 0; i < swift_type_count; i++) {
+ count += index->nparts[i];
+ }
+ count *= sizeof(struct index_data);
+
+ /* Read the number of new particles. */
+ memcpy(&index->nparts_created,
+ (char *)index->index.map + logger_index_data_offset + count,
+ logger_index_npart_size);
+
+ /* Compute the position of the history of particles removed. */
+ size_t count_new = 0;
+ for (int i = 0; i < swift_type_count; i++) {
+ count_new += index->nparts_created[i];
+ }
+ count_new *= sizeof(struct index_data);
+
+ /* Read the number of particles removed. */
+ memcpy(&index->nparts_removed,
+ (char *)index->index.map + logger_index_data_offset + count +
+ logger_index_npart_size + count_new,
+ logger_index_npart_size);
+
/* Close the file */
logger_index_free(index);
}
@@ -104,6 +129,9 @@ void logger_index_write_sorted(struct logger_index *index) {
/**
* @brief Get the #index_data of a particle type.
*
+ * The #index_data contains the ids and offset of the particle.
+ * This function should be used when looking for the last known particles.
+ *
* @param index The #logger_index.
* @param type The particle type.
*/
@@ -121,6 +149,60 @@ struct index_data *logger_index_get_data(struct logger_index *index, int type) {
return (struct index_data *)ret;
}
+/**
+ * @brief Get the #index_data for the particles created.
+ *
+ * The #index_data contains the ids and offset of the particle.
+ * This function should be used when looking for the particles created.
+ *
+ * @param index The #logger_index.
+ * @param type The particle type.
+ */
+struct index_data *logger_index_get_created_history(struct logger_index *index,
+ int type) {
+
+ /* Get the position after the previous array. */
+ char *ret = (char *)logger_index_get_data(index, swift_type_count);
+
+ /* Get the offset of particles of this type by skipping entries of lower
+ * types. */
+ size_t count = 0;
+ for (int i = 0; i < type; i++) {
+ count += index->nparts_created[i];
+ }
+ const size_t type_offset = count * sizeof(struct index_data);
+
+ ret += type_offset + logger_index_npart_size;
+ return (struct index_data *)ret;
+}
+
+/**
+ * @brief Get the #index_data for the particles removed.
+ *
+ * The #index_data contains the ids and offset of the particle.
+ * This function should be used when looking for the particles removed.
+ *
+ * @param index The #logger_index.
+ * @param type The particle type.
+ */
+struct index_data *logger_index_get_removed_history(struct logger_index *index,
+ int type) {
+
+ /* Get the position after the previous array. */
+ char *ret = (char *)logger_index_get_created_history(index, swift_type_count);
+
+ /* Get the offset of particles of this type by skipping entries of lower
+ * types. */
+ size_t count = 0;
+ for (int i = 0; i < type; i++) {
+ count += index->nparts_removed[i];
+ }
+ const size_t type_offset = count * sizeof(struct index_data);
+
+ ret += type_offset + logger_index_npart_size;
+ return (struct index_data *)ret;
+}
+
/**
* @brief Map the file and if required sort it.
*
diff --git a/logger/logger_index.h b/logger/logger_index.h
index 8419e0aaa18e00846e2c3ad9f1218780bf44f431..1a50ed812acd3b317c52f4ba7c040e55f131a03d 100644
--- a/logger/logger_index.h
+++ b/logger/logger_index.h
@@ -68,6 +68,12 @@ struct logger_index {
/* The mapped file */
struct mapped_file index;
+
+ /* Number of particles created. */
+ uint64_t nparts_created[swift_type_count];
+
+ /* Number of particles removed. */
+ uint64_t nparts_removed[swift_type_count];
};
void logger_index_write_sorted(struct logger_index *index);
@@ -81,5 +87,9 @@ size_t logger_index_get_particle_offset(struct logger_index *index,
void logger_index_free(struct logger_index *index);
void logger_index_sort_file(struct logger_index *index);
struct index_data *logger_index_get_data(struct logger_index *index, int type);
+struct index_data *logger_index_get_created_history(struct logger_index *index,
+ int type);
+struct index_data *logger_index_get_removed_history(struct logger_index *index,
+ int type);
#endif // LOGGER_LOGGER_INDEX_H
diff --git a/logger/logger_interpolation.h b/logger/logger_interpolation.h
new file mode 100644
index 0000000000000000000000000000000000000000..c573e8395217a44049c215174514fcd4b57e2d4d
--- /dev/null
+++ b/logger/logger_interpolation.h
@@ -0,0 +1,310 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 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/>.
+ *
+ ******************************************************************************/
+
+#ifndef LOGGER_LOGGER_INTERPOLATION_H
+#define LOGGER_LOGGER_INTERPOLATION_H
+
+#include "logger_tools.h"
+
+/**
+ * @brief Compute the quintic hermite spline interpolation.
+ * See https://en.wikipedia.org/wiki/Hermite_interpolation for more information.
+ * @f$ p(x) = f(x_0) + f'(x_0) (x - x_0) + \frac{1}{2}f''(x_0) (x - x_0)^2 +
+ \frac{f(x_1) - f(x_0) - f'(x_0) (x_1 - x_0) - \frac{1}{2} f''(x_0) (x_1 -
+ x_0)^2}{(x_1 - x_0)^3} (x - x_0)^3
+ + \frac{3 f(x_0) - 3 f(x_1) + \left( 2 f'(x_0) + f'(x_1) \right) (x_1 - x_0) +
+ \frac{1}{2} f''(x_0) (x_1 - x_0)^2}{(x_1 - x_0)^4} (x - x_0)^3 (x - x_1)
+ + \frac{6 f(x_1) - 6 f(x_0) - 3 \left( f'(x_0) + f'(x_1) \right) (x_1 - x_0) +
+ \frac{1}{2}\left( f''(x_1) - f''(x_0) \right) (x_1 - x_0)^2}{(x_1 - x_0)^5} (x
+ - x_0)^3 (x - x_1)^2 @f$
+ *
+ * @param t0 The time at the left of the interval.
+ * @param x0 The function at the left of the interval.
+ * @param v0 The first derivative at the left of the interval.
+ * @param a0 The second derivative at the left of the interval.
+ * @param t1 The time at the right of the interval.
+ * @param x1 The function at the right of the interval.
+ * @param v1 The first derivative at the right of the interval.
+ * @param a1 The second derivative at the right of the interval.
+ * @param t The time of the interpolation.
+ *
+ * @return The function evaluated at t.
+ */
+__attribute__((always_inline)) INLINE static double
+interpolate_quintic_hermite_spline(const double t0, const double x0,
+ const float v0, const float a0,
+ const double t1, const double x1,
+ const float v1, const float a1,
+ const double t) {
+
+ /* Generates recurring variables */
+ /* Time differences */
+ const double dt = t1 - t0;
+ const double dt2 = dt * dt;
+ const double dt3 = dt2 * dt;
+ const double dt4 = dt3 * dt;
+ const double dt5 = dt4 * dt;
+
+ const double t_t0 = t - t0;
+ const double t_t0_2 = t_t0 * t_t0;
+ const double t_t0_3 = t_t0_2 * t_t0;
+ const double t_t1 = t - t1;
+ const double t_t1_2 = t_t1 * t_t1;
+
+ /* Derivatives */
+ const double v0_dt = v0 * dt;
+ const double a0_dt2 = 0.5 * a0 * dt2;
+ const double v1_dt = v1 * dt;
+ const double a1_dt2 = 0.5 * a1 * dt2;
+
+ /* Do the first 3 terms of the hermite spline */
+ double x = x0 + v0 * t_t0 + 0.5 * a0 * t_t0_2;
+
+ /* Cubic term */
+ x += (x1 - x0 - v0_dt - a0_dt2) * t_t0_3 / dt3;
+
+ /* Quartic term */
+ x += (3. * x0 - 3. * x1 + v1_dt + 2. * v0_dt + a0_dt2) * t_t0_3 * t_t1 / dt4;
+
+ /* Quintic term */
+ x += (6. * x1 - 6. * x0 - 3. * v0_dt - 3. * v1_dt + a1_dt2 - a0_dt2) *
+ t_t0_3 * t_t1_2 / dt5;
+
+ return x;
+}
+
+/**
+ * @brief Compute the cubic hermite spline interpolation.
+ * See https://en.wikipedia.org/wiki/Hermite_interpolation for more information.
+ * @f$ p(x) = f(x_0) + f'(x_0) (x - x_0) + \frac{f(x_1) - f(x_0) - f'(x_0) (x_1
+ - x_0)}{(x_1 - x_0)^2} (x - x_0)^2
+ + \frac{2 f(x_0) - 2 f(x_1) + f'(x_0) (x_1 - x_0) + f'(x_1) (x_1 - x_0)}{(x_1
+ - x_0)^3} (x - x_0)^2 (x - x_1) @f$
+ *
+ * @param t0 The time at the left of the interval.
+ * @param v0 The first derivative at the left of the interval.
+ * @param a0 The second derivative at the left of the interval.
+ * @param t1 The time at the right of the interval.
+ * @param v1 The first derivative at the right of the interval.
+ * @param a1 The second derivative at the right of the interval.
+ * @param t The time of the interpolation.
+ *
+ * @return The function evaluated at t.
+ */
+__attribute__((always_inline)) INLINE static float
+interpolate_cubic_hermite_spline(const double t0, const float v0,
+ const float a0, const double t1,
+ const float v1, const float a1,
+ const double t) {
+
+ /* Generates recurring variables */
+ /* Time differences */
+ const float dt = t1 - t0;
+ const float dt2 = dt * dt;
+ const float dt3 = dt2 * dt;
+
+ const float t_t0 = t - t0;
+ const float t_t0_2 = t_t0 * t_t0;
+ const float t_t1 = t - t1;
+
+ /* Derivatives */
+ const float a0_dt = a0 * dt;
+ const float a1_dt = a1 * dt;
+
+ /* Do the first 2 terms of the hermite spline */
+ float x = v0 + a0 * t_t0;
+
+ /* Square term */
+ x += (v1 - v0 - a0_dt) * t_t0_2 / dt2;
+
+ /* Cubic term */
+ x += (2. * v0 - 2. * v1 + a1_dt + a0_dt) * t_t0_2 * t_t1 / dt3;
+
+ return x;
+}
+
+/**
+ * @brief Interpolates a field in N dimension using the first and second
+ * derivatives if possible.
+ *
+ * The field is in double and both derivatives in float.
+ *
+ * @param before Pointer to the #logger_field at a time < t.
+ * @param after Pointer to the #logger_field at a time > t.
+ * @param otuput Pointer to the output value.
+ * @param t_before Time of field_before (< t).
+ * @param t_after Time of field_after (> t).
+ * @param t Requested time.
+ */
+__attribute__((always_inline)) INLINE static void
+interpolate_quintic_double_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) {
+ /* Compute the interpolation scaling. */
+ const double wa = (t - t_before) / (t_after - t_before);
+ const double wb = 1. - wa;
+
+ 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;
+
+ const double *x_aft = (double *)after->field;
+ const float *v_aft = (float *)after->first_deriv;
+ const float *a_aft = (float *)after->second_deriv;
+
+ /* 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);
+ }
+ /* 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);
+ }
+ /* Use linear interpolation. */
+ else {
+ x[i] = wa * x_aft[i] + wb * x_bef[i];
+ }
+ }
+}
+
+/**
+ * @brief Interpolates a field in N dimension using the first derivative if
+ * possible.
+ *
+ * The field and the first derivatives are floats.
+ *
+ * @param before Pointer to the #logger_field at a time < t.
+ * @param after Pointer to the #logger_field at a time > t.
+ * @param otuput Pointer to the output value.
+ * @param t_before Time of field_before (< t).
+ * @param t_after Time of field_after (> t).
+ * @param t Requested time.
+ */
+__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) {
+
+ /* Compute the interpolation scaling. */
+ const float wa = (t - t_before) / (t_after - t_before);
+ const float wb = 1. - wa;
+
+ for (int i = 0; i < dimension; i++) {
+ float *v = (float *)output;
+ const float *v_bef = (float *)before->field;
+ const float *a_bef = (float *)before->first_deriv;
+
+ const float *v_aft = (float *)after->field;
+ const float *a_aft = (float *)after->first_deriv;
+
+ /* Use a cubic hermite spline. */
+ if (a_bef && a_aft) {
+ v[i] = interpolate_cubic_hermite_spline(t_before, v_bef[i], a_bef[i],
+ t_after, v_aft[i], a_aft[i], t);
+ }
+ /* Use linear interpolation. */
+ else {
+ v[i] = wa * v_aft[i] + wb * v_bef[i];
+ }
+ }
+}
+
+/**
+ * @brief Interpolates a field in N dimension.
+ *
+ * The field is in float.
+ *
+ * @param before Pointer to the #logger_field at a time < t.
+ * @param after Pointer to the #logger_field at a time > t.
+ * @param otuput Pointer to the output value.
+ * @param t_before Time of field_before (< t).
+ * @param t_after Time of field_after (> t).
+ * @param t Requested time.
+ */
+__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) {
+
+ /* Compute the interpolation scaling. */
+ const float wa = (t - t_before) / (t_after - t_before);
+ const float wb = 1. - wa;
+
+ /* interpolate vectors. */
+ for (int i = 0; i < dimension; i++) {
+ float *a = (float *)output;
+ const float *a_bef = (float *)before->field;
+ const float *a_aft = (float *)after->field;
+ a[i] = wa * a_aft[i] + wb * a_bef[i];
+ }
+}
+
+/**
+ * @brief Interpolates a field stored as a float.
+ *
+ * @param before Pointer to the #logger_field at a time < t.
+ * @param after Pointer to the #logger_field at a time > t.
+ * @param otuput Pointer to the output value.
+ * @param t_before Time of field_before (< t).
+ * @param t_after Time of field_after (> t).
+ * @param t Requested time.
+ */
+__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) {
+
+ /* Compute the interpolation scaling. */
+ const float wa = (t - t_before) / (t_after - t_before);
+ const float wb = 1. - wa;
+ ((float *)output)[0] =
+ wa * ((float *)after->field)[0] + wb * ((float *)before->field)[0];
+}
+
+/**
+ * @brief Interpolation function for the ids.
+ *
+ * As the IDs should not change during the simulation, we ensure
+ * that the IDs are the same between two records.
+ *
+ * @param before Pointer to the #logger_field at a time < t.
+ * @param after Pointer to the #logger_field at a time > t.
+ * @param otuput Pointer to the output value.
+ * @param t_before Time of field_before (< t).
+ * @param t_after Time of field_after (> t).
+ * @param t Requested time.
+ */
+__attribute__((always_inline)) INLINE static void interpolate_ids(
+ 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) {
+ if (*(long long *)after->field != *(long long *)before->field) {
+ error("Interpolating different particles");
+ }
+ *(long long *)output = *(long long *)after->field;
+}
+#endif // LOGGER_LOGGER_INTERPOLATION_H
diff --git a/logger/logger_particle.c b/logger/logger_particle.c
index 331fd74190d603ffc379df8d7a6abf4712fbfa50..4e01c38fd942a2ec356207940730aa7bbf730a5d 100644
--- a/logger/logger_particle.c
+++ b/logger/logger_particle.c
@@ -66,7 +66,12 @@ __attribute__((always_inline)) INLINE size_t logger_particle_read_field(
/* Check if it is not a time record. */
if (*mask == h->timestamp_mask) {
- error_python("Cannot particle from timestep record.");
+ error_python("Cannot read a particle from timestep record.");
+ }
+
+ /* Skip the special field. */
+ if (*mask & h->masks[logger_index_special_flags].mask) {
+ map += h->masks[logger_index_special_flags].size;
}
/* Read the record and copy it to a particle. */
@@ -86,5 +91,50 @@ __attribute__((always_inline)) INLINE size_t logger_particle_read_field(
map += h->masks[global].size;
}
}
+
return map - reader->log.log.map;
}
+
+/**
+ * @brief Read the special flag of a particle (of any type) in the log file.
+ *
+ * @param reader The #logger_reader.
+ * @param offset offset of the record to read.
+ * @param mask (out) The mask of the record.
+ * @param h_offset (out) Difference of offset with the next record.
+ * @param data (out) The data of the flag.
+ *
+ * @return The special flag.
+ */
+__attribute__((always_inline)) INLINE enum logger_special_flags
+logger_particle_read_special_flag(const struct logger_reader *reader,
+ size_t offset, size_t *mask, size_t *h_offset,
+ int *data) {
+
+ /* Get a few pointers. */
+ const struct header *h = &reader->log.header;
+ void *map = reader->log.log.map;
+
+ *mask = 0;
+ *h_offset = 0;
+
+ /* Read the record's mask. */
+ map = logger_loader_io_read_mask(h, (char *)map + offset, mask, h_offset);
+
+ /* Check that the mask is within the limits. */
+ if (*mask > (unsigned int)(1 << h->masks_count)) {
+ error_python("Found an unexpected mask %zi", *mask);
+ }
+
+ /* Check if it is not a time record. */
+ if (*mask == h->timestamp_mask) {
+ error_python("Cannot read a particle from timestep record.");
+ }
+
+ /* Read the special flag */
+ uint32_t packed_data = 0;
+ map = logger_loader_io_read_data(
+ map, h->masks[logger_index_special_flags].size, &packed_data);
+
+ return logger_unpack_flags_and_data(packed_data, data);
+}
diff --git a/logger/logger_particle.h b/logger/logger_particle.h
index d76dfe917190d3aa3c8c98b7aa13768b79654358..22dc6460b05275ae50419d93371ae827641cafc2 100644
--- a/logger/logger_particle.h
+++ b/logger/logger_particle.h
@@ -45,6 +45,11 @@ size_t logger_particle_read_field(const struct logger_reader *reader,
const int *local_to_global,
const int local_count, int field,
size_t *mask, size_t *h_offset);
+
+enum logger_special_flags logger_particle_read_special_flag(
+ const struct logger_reader *reader, size_t offset, size_t *mask,
+ size_t *h_offset, int *data);
+
/**
* @brief Generate the data for the special flags.
*
diff --git a/logger/logger_python_wrapper.c b/logger/logger_python_wrapper.c
index 9c7336dacd09e868cab2f879a398f38e2ec6a4bd..29d6852c4c908a9192471bc7e9b75f7fd38ab5ea 100644
--- a/logger/logger_python_wrapper.c
+++ b/logger/logger_python_wrapper.c
@@ -428,14 +428,45 @@ static PyObject *pyGetParticleData(__attribute__((unused)) PyObject *self,
/* input variables. */
PyObject *fields = NULL;
double time = 0;
+ PyObject *types = Py_None;
+
/* parse the arguments. */
- if (!PyArg_ParseTuple(args, "Od", &fields, &time)) return NULL;
+ if (!PyArg_ParseTuple(args, "Od|O", &fields, &time, &types)) return NULL;
/* Check the inputs. */
if (!PyList_Check(fields)) {
error_python("Expecting a list of fields");
}
+ /* Get the type of particles to read. */
+ int read_types[swift_type_count] = {0};
+ /* By default, we read everything */
+ if (types == Py_None) {
+ for (int i = 0; i < swift_type_count; i++) {
+ read_types[i] = 1;
+ }
+ }
+ /* Deal with the case of a single int. */
+ else if (PyLong_Check(types)) {
+ const size_t type = PyLong_AsSize_t(types);
+ if (type >= swift_type_count) {
+ error_python("Unexpected particle type %zi", type);
+ }
+ read_types[type] = 1;
+ }
+ /* Deal with the case of a list */
+ else if (PyList_Check(types)) {
+ const size_t size = PyList_Size(types);
+ for (size_t i = 0; i < size; i++) {
+ PyObject *cur = PyList_GetItem(types, i);
+ const size_t type = PyLong_AsSize_t(cur);
+ if (type >= swift_type_count) {
+ error_python("Unexpected particle type %zi", type);
+ }
+ read_types[type] = 1;
+ }
+ }
+
/* initialize the reader. */
struct logger_reader *reader = &self_reader->reader;
const struct header *h = &reader->log.header;
@@ -474,10 +505,10 @@ static PyObject *pyGetParticleData(__attribute__((unused)) PyObject *self,
logger_reader_set_time(reader, time);
/* Get the number of particles. */
- int n_type = 0;
- const uint64_t *n_part = logger_reader_get_number_particles(reader, &n_type);
+ uint64_t n_part[swift_type_count];
+ logger_reader_get_number_particles(reader, n_part, read_types);
uint64_t n_tot = 0;
- for (int i = 0; i < n_type; i++) {
+ for (int i = 0; i < swift_type_count; i++) {
n_tot += n_part[i];
}
diff --git a/logger/logger_reader.c b/logger/logger_reader.c
index feb757fd264cba983134754d873839c57dddd38b..493ba6d6d81ae781baf5ee61b1cd07b88d20b56f 100644
--- a/logger/logger_reader.c
+++ b/logger/logger_reader.c
@@ -72,7 +72,8 @@ void logger_reader_init(struct logger_reader *reader, const char *basename,
*/
void logger_reader_init_index(struct logger_reader *reader) {
/* Initialize the logger_index */
- logger_index_init(&reader->index.index, reader);
+ logger_index_init(&reader->index.index_prev, reader);
+ logger_index_init(&reader->index.index_next, reader);
/* Count the number of files */
int count = 0;
@@ -106,11 +107,11 @@ void logger_reader_init_index(struct logger_reader *reader) {
sprintf(filename, "%s_%04i.index", reader->basename, i);
/* Read the header */
- logger_index_read_header(&reader->index.index, filename);
+ logger_index_read_header(&reader->index.index_prev, filename);
/* Save the required information */
- reader->index.times[i] = reader->index.index.time;
- reader->index.int_times[i] = reader->index.index.integer_time;
+ reader->index.times[i] = reader->index.index_prev.time;
+ reader->index.int_times[i] = reader->index.index_prev.integer_time;
}
}
@@ -124,7 +125,7 @@ void logger_reader_free(struct logger_reader *reader) {
logger_logfile_free(&reader->log);
if (reader->time.time != -1.) {
- logger_index_free(&reader->index.index);
+ logger_index_free(&reader->index.index_prev);
}
free(reader->index.int_times);
@@ -156,17 +157,27 @@ void logger_reader_set_time(struct logger_reader *reader, double time) {
}
/* Generate the filename */
- char filename[STRING_SIZE + 50];
- sprintf(filename, "%s_%04u.index", reader->basename, left);
+ char filename_prev[STRING_SIZE + 50];
+ sprintf(filename_prev, "%s_%04u.index", reader->basename, left);
+ char filename_next[STRING_SIZE + 50];
+ sprintf(filename_next, "%s_%04u.index", reader->basename, left + 1);
/* Check if the file is already mapped */
- if (reader->index.index.index.map != NULL) {
- logger_index_free(&reader->index.index);
+ if (reader->index.index_prev.index.map != NULL) {
+ logger_index_free(&reader->index.index_prev);
+ }
+ if (reader->index.index_next.index.map != NULL) {
+ logger_index_free(&reader->index.index_next);
}
/* Read the file */
- logger_index_read_header(&reader->index.index, filename);
- logger_index_map_file(&reader->index.index, filename, /* sorted */ 1);
+ logger_index_read_header(&reader->index.index_prev, filename_prev);
+ logger_index_map_file(&reader->index.index_prev, filename_prev,
+ /* sorted */ 1);
+
+ logger_index_read_header(&reader->index.index_next, filename_next);
+ logger_index_map_file(&reader->index.index_next, filename_next,
+ /* sorted */ 1);
/* Get the offset of the time chunk */
size_t ind = time_array_get_index_from_time(&reader->log.times, time);
@@ -189,17 +200,133 @@ void logger_reader_set_time(struct logger_reader *reader, double time) {
}
/**
- * @brief Provides the number of particle (per type) from the index file.
+ * @brief Count the number of new particles at the time set since last index
+ * file.
+ *
+ * @param reader The #logger_reader.
+ * @param part_type The index corresponding to the particle type.
+ *
+ * @param The number of new particles.
+ */
+size_t logger_reader_count_number_new_particles(struct logger_reader *reader,
+ enum part_type part_type) {
+
+ const size_t threshold = reader->time.time_offset;
+
+ /* Get the history of created particles. */
+ struct index_data *data =
+ logger_index_get_created_history(&reader->index.index_next, part_type);
+
+ /* Do we have any new particle? */
+ if (reader->index.index_next.nparts_created[part_type] == 0 ||
+ threshold < data[0].offset) {
+ return 0;
+ }
+
+ /* Perform a binary search */
+ size_t right = reader->index.index_next.nparts_created[part_type] - 1;
+ size_t left = 0;
+ while (left <= right) {
+ size_t m = (left + right) / 2;
+ if (data[m].offset < threshold) {
+ left = m + 1;
+ } else if (data[m].offset > threshold) {
+ right = m - 1;
+ } else {
+ return m;
+ }
+ }
+
+ /* Compute the return value. */
+ size_t ret = (left + right) / 2;
+
+ /* Check if the binary search is correct. */
+ if (data[ret].offset > threshold) {
+ error_python("The binary search has failed.");
+ }
+
+ /* We need the count, not the index. */
+ return ret + 1;
+}
+
+/**
+ * @brief Count the number of removed particles at the time set since last index
+ * file.
*
* @param reader The #logger_reader.
* @param n_type (output) The number of particle type possible.
+ * @param part_type The index corresponding to the particle type.
*
- * @return For each type possible, the number of particle.
+ * @param The number of removed particles.
*/
-const uint64_t *logger_reader_get_number_particles(struct logger_reader *reader,
- int *n_type) {
- *n_type = swift_type_count;
- return reader->index.index.nparts;
+size_t logger_reader_count_number_removed_particles(
+ struct logger_reader *reader, enum part_type part_type) {
+
+ const size_t threshold = reader->time.time_offset;
+
+ /* Get the history of created particles. */
+ struct index_data *data =
+ logger_index_get_removed_history(&reader->index.index_next, part_type);
+
+ /* Do we have any new particle? */
+ if (reader->index.index_next.nparts_removed[part_type] == 0 ||
+ threshold < data[0].offset) {
+ return 0;
+ }
+
+ /* Perform a binary search */
+ size_t right = reader->index.index_next.nparts_removed[part_type] - 1;
+ size_t left = 0;
+ while (left <= right) {
+ size_t m = (left + right) / 2;
+ if (data[m].offset < threshold) {
+ left = m + 1;
+ } else if (data[m].offset > threshold) {
+ right = m - 1;
+ } else {
+ return m;
+ }
+ }
+
+ /* Compute the return value. */
+ size_t ret = (left + right) / 2;
+
+ /* Check if the binary search is correct. */
+ if (data[ret].offset > threshold) {
+ error_python("The binary search has failed.");
+ }
+
+ /* We need the count, not the index. */
+ return ret + 1;
+}
+
+/**
+ * @brief Provides the number of particle (per type) from the index file.
+ *
+ * @param reader The #logger_reader.
+ * @param n_parts (out) Number of particles at the time set in the reader.
+ * @param read_types Should we read this type of particles?
+ */
+void logger_reader_get_number_particles(struct logger_reader *reader,
+ uint64_t *n_parts,
+ const int *read_types) {
+ for (int i = 0; i < swift_type_count; i++) {
+ /* Should we skip this type of particles? */
+ if (read_types[i] == 0) {
+ n_parts[i] = 0;
+ continue;
+ }
+
+ /* Count the number of particles present in the last index file. */
+ const uint64_t count_prev = reader->index.index_prev.nparts[i];
+ /* Count the number of particles that have been created since last index. */
+ const uint64_t count_new =
+ logger_reader_count_number_new_particles(reader, i);
+ /* Count the number of particles that have been removed since last index. */
+ const uint64_t count_removed =
+ logger_reader_count_number_removed_particles(reader, i);
+ n_parts[i] = (count_prev + count_new) - count_removed;
+ }
}
/**
@@ -227,18 +354,24 @@ const uint64_t *logger_reader_get_number_particles(struct logger_reader *reader,
(-1 if none)
* @param output Pointers to the output array
* @param type The #part_type.
+ *
+ * @return Is the particle removed from the logfile?
*/
-void logger_reader_read_field(struct logger_reader *reader, double time,
- size_t offset_time,
- enum logger_reader_type interp_type,
- const size_t offset_last_full_record,
- const int field, const int first,
- const int second, void *output,
- enum part_type type) {
+int logger_reader_read_field(struct logger_reader *reader, double time,
+ size_t offset_time,
+ enum logger_reader_type interp_type,
+ const size_t offset_last_full_record,
+ const int field, const int first, const int second,
+ void *output, enum part_type type) {
const struct header *h = &reader->log.header;
size_t offset = offset_last_full_record;
+ /* Check if the offsets are correct. */
+ if (offset > offset_time) {
+ error_python("Last offset is larger than the requested time.");
+ }
+
/* Get the particle type variables. */
int *local_to_global = NULL;
int field_count = 0;
@@ -248,6 +381,7 @@ void logger_reader_read_field(struct logger_reader *reader, double time,
field_count = hydro_logger_field_count;
break;
case swift_type_dark_matter:
+ case swift_type_dark_matter_background:
local_to_global = gravity_logger_local_to_global;
field_count = gravity_logger_field_count;
break;
@@ -285,6 +419,22 @@ void logger_reader_read_field(struct logger_reader *reader, double time,
logger_loader_io_read_mask(h, reader->log.log.map + offset, &mask,
&h_offset);
+ /* Is the particle removed from the logfile? */
+ if (mask & h->masks[logger_index_special_flags].mask) {
+ int data = 0;
+ enum logger_special_flags flag = logger_particle_read_special_flag(
+ reader, offset, &mask, &h_offset, &data);
+ if (flag == logger_flag_change_type || flag == logger_flag_mpi_exit ||
+ flag == logger_flag_delete) {
+ return 1;
+ }
+ }
+
+ /* Check if there is a next record (avoid infinite loop). */
+ if (h_offset == 0) {
+ error_python("No next record found.");
+ }
+
/* Check if the field is present */
if (mask & mask_field->mask) {
offset_before = offset;
@@ -405,6 +555,7 @@ void logger_reader_read_field(struct logger_reader *reader, double time,
output, time, field);
break;
case swift_type_dark_matter:
+ case swift_type_dark_matter_background:
gravity_logger_interpolate_field(time_before, &before, time_after,
&after, output, time, field);
break;
@@ -416,6 +567,7 @@ void logger_reader_read_field(struct logger_reader *reader, double time,
error_python("Particle type not implemented");
}
}
+ return 0;
}
/**
@@ -443,22 +595,19 @@ void logger_reader_global_to_local(
/* Get the correct variables. */
int n_max = 0;
int *local_to_global = NULL;
- const char **local_names = NULL;
switch (type) {
case swift_type_gas:
n_max = hydro_logger_field_count;
local_to_global = hydro_logger_local_to_global;
- local_names = hydro_logger_field_names;
break;
case swift_type_dark_matter:
+ case swift_type_dark_matter_background:
n_max = gravity_logger_field_count;
local_to_global = gravity_logger_local_to_global;
- local_names = gravity_logger_field_names;
break;
case swift_type_stars:
n_max = stars_logger_field_count;
local_to_global = stars_logger_local_to_global;
- local_names = stars_logger_field_names;
break;
default:
error_python("Particle type not implemented yet.");
@@ -496,7 +645,9 @@ void logger_reader_global_to_local(
/* Check that we found the fields */
for (int local = 0; local < n_fields_wanted; local++) {
if (local_fields_wanted[local] < 0) {
- error_python("Field %s not found in particle type %s", local_names[local],
+ const int global_field = global_fields_wanted[local];
+ const char *name = h->masks[global_field].name;
+ error_python("Field %s not found in particle type %s", name,
part_type_names[type]);
}
}
@@ -546,61 +697,185 @@ void logger_reader_read_all_particles(struct logger_reader *reader, double time,
/* Do the hydro. */
if (n_part[swift_type_gas] != 0) {
struct index_data *data =
- logger_index_get_data(&reader->index.index, swift_type_gas);
+ logger_index_get_data(&reader->index.index_prev, swift_type_gas);
+ struct index_data *data_created = logger_index_get_created_history(
+ &reader->index.index_next, swift_type_gas);
/* Sort the fields in order to read the correct bits. */
logger_reader_global_to_local(
reader, global_fields_wanted, local_fields_wanted, local_first_deriv,
local_second_deriv, n_fields_wanted, swift_type_gas);
+ size_t current_in_index = 0;
+ int reading_history = 0;
+ const size_t size_index = reader->index.index_prev.nparts[swift_type_gas];
+ const size_t size_history =
+ logger_reader_count_number_new_particles(reader, swift_type_gas);
+
/* Read the particles */
for (size_t i = 0; i < n_part[swift_type_gas]; i++) {
- /* Get the offset */
- size_t offset = data[i].offset;
-
- /* Sort the output into output_single. */
- for (int field = 0; field < n_fields_wanted; field++) {
- const int global = global_fields_wanted[field];
- const int local = local_fields_wanted[field];
- const int first = local_first_deriv[field];
- const int second = local_second_deriv[field];
- void *output_single = output[field] + i * h->masks[global].size;
-
- /* Read the field. */
- logger_reader_read_field(reader, time, reader->time.time_offset,
- interp_type, offset, local, first, second,
- output_single, swift_type_gas);
+ int particle_removed = 1;
+ /* Do it until finding a particle not removed. */
+ while (particle_removed) {
+ /* Should we start to read the history? */
+ if (!reading_history && current_in_index == size_index) {
+ current_in_index = 0;
+ reading_history = 1;
+ }
+
+ /* Check if we still have some particles available. */
+ if (reading_history && current_in_index == size_history) {
+ error_python("The logger was not able to find enough particles.");
+ }
+
+ /* Get the offset */
+ size_t offset = reading_history ? data_created[current_in_index].offset
+ : data[current_in_index].offset;
+
+ /* Sort the output into output_single. */
+ for (int field = 0; field < n_fields_wanted; field++) {
+ const int global = global_fields_wanted[field];
+ const int local = local_fields_wanted[field];
+ const int first = local_first_deriv[field];
+ const int second = local_second_deriv[field];
+ void *output_single = output[field] + i * h->masks[global].size;
+
+ /* Read the field. */
+ particle_removed = logger_reader_read_field(
+ reader, time, reader->time.time_offset, interp_type, offset,
+ local, first, second, output_single, swift_type_gas);
+
+ /* Should we continue to read the fields of this particle? */
+ if (particle_removed) {
+ break;
+ }
+ }
+ current_in_index++;
}
}
}
/* Do the dark matter. */
if (n_part[swift_type_dark_matter] != 0) {
- struct index_data *data =
- logger_index_get_data(&reader->index.index, swift_type_dark_matter);
+ struct index_data *data = logger_index_get_data(&reader->index.index_prev,
+ swift_type_dark_matter);
+ struct index_data *data_created = logger_index_get_created_history(
+ &reader->index.index_next, swift_type_dark_matter);
/* Sort the fields in order to read the correct bits. */
logger_reader_global_to_local(
reader, global_fields_wanted, local_fields_wanted, local_first_deriv,
local_second_deriv, n_fields_wanted, swift_type_dark_matter);
+ size_t current_in_index = 0;
+ int reading_history = 0;
+ const size_t size_index =
+ reader->index.index_prev.nparts[swift_type_dark_matter];
+ const size_t size_history = logger_reader_count_number_new_particles(
+ reader, swift_type_dark_matter);
+
/* Read the particles */
for (size_t i = 0; i < n_part[swift_type_dark_matter]; i++) {
- /* Get the offset */
- size_t offset = data[i].offset;
-
- /* Sort the output into output_single. */
- for (int field = 0; field < n_fields_wanted; field++) {
- const int global = global_fields_wanted[field];
- const int local = local_fields_wanted[field];
- const int first = local_first_deriv[field];
- const int second = local_second_deriv[field];
- void *output_single = output[field] + i * h->masks[global].size;
-
- /* Read the field. */
- logger_reader_read_field(reader, time, reader->time.time_offset,
- interp_type, offset, local, first, second,
- output_single, swift_type_dark_matter);
+ int particle_removed = 1;
+ /* Do it until finding a particle not removed. */
+ while (particle_removed) {
+ /* Should we start to read the history? */
+ if (!reading_history && current_in_index == size_index) {
+ current_in_index = 0;
+ reading_history = 1;
+ }
+
+ /* Check if we still have some particles available. */
+ if (reading_history && current_in_index == size_history) {
+ error_python("The logger was not able to find enough particles.");
+ }
+
+ /* Get the offset */
+ size_t offset = reading_history ? data_created[current_in_index].offset
+ : data[current_in_index].offset;
+
+ /* Sort the output into output_single. */
+ for (int field = 0; field < n_fields_wanted; field++) {
+ const int global = global_fields_wanted[field];
+ const int local = local_fields_wanted[field];
+ const int first = local_first_deriv[field];
+ const int second = local_second_deriv[field];
+ void *output_single = output[field] + i * h->masks[global].size;
+
+ /* Read the field. */
+ particle_removed = logger_reader_read_field(
+ reader, time, reader->time.time_offset, interp_type, offset,
+ local, first, second, output_single, swift_type_dark_matter);
+
+ /* Should we continue to read the fields of this particle? */
+ if (particle_removed) {
+ break;
+ }
+ }
+ current_in_index++;
+ }
+ }
+ }
+
+ /* Do the dark matter background. */
+ if (n_part[swift_type_dark_matter_background] != 0) {
+ struct index_data *data = logger_index_get_data(
+ &reader->index.index_prev, swift_type_dark_matter_background);
+ struct index_data *data_created = logger_index_get_created_history(
+ &reader->index.index_next, swift_type_dark_matter_background);
+
+ /* Sort the fields in order to read the correct bits. */
+ logger_reader_global_to_local(
+ reader, global_fields_wanted, local_fields_wanted, local_first_deriv,
+ local_second_deriv, n_fields_wanted, swift_type_dark_matter_background);
+
+ size_t current_in_index = 0;
+ int reading_history = 0;
+ const size_t size_index =
+ reader->index.index_prev.nparts[swift_type_dark_matter_background];
+ const size_t size_history = logger_reader_count_number_new_particles(
+ reader, swift_type_dark_matter_background);
+
+ /* Read the particles */
+ for (size_t i = 0; i < n_part[swift_type_dark_matter_background]; i++) {
+ int particle_removed = 1;
+ /* Do it until finding a particle not removed. */
+ while (particle_removed) {
+ /* Should we start to read the history? */
+ if (!reading_history && current_in_index == size_index) {
+ current_in_index = 0;
+ reading_history = 1;
+ }
+
+ /* Check if we still have some particles available. */
+ if (reading_history && current_in_index == size_history) {
+ error_python("The logger was not able to find enough particles.");
+ }
+
+ /* Get the offset */
+ size_t offset = reading_history ? data_created[current_in_index].offset
+ : data[current_in_index].offset;
+
+ /* Sort the output into output_single. */
+ for (int field = 0; field < n_fields_wanted; field++) {
+ const int global = global_fields_wanted[field];
+ const int local = local_fields_wanted[field];
+ const int first = local_first_deriv[field];
+ const int second = local_second_deriv[field];
+ void *output_single = output[field] + i * h->masks[global].size;
+
+ /* Read the field. */
+ particle_removed = logger_reader_read_field(
+ reader, time, reader->time.time_offset, interp_type, offset,
+ local, first, second, output_single,
+ swift_type_dark_matter_background);
+
+ /* Should we continue to read the fields of this particle? */
+ if (particle_removed) {
+ break;
+ }
+ }
+ current_in_index++;
}
}
}
@@ -608,30 +883,60 @@ void logger_reader_read_all_particles(struct logger_reader *reader, double time,
/* Do the stars. */
if (n_part[swift_type_stars] != 0) {
struct index_data *data =
- logger_index_get_data(&reader->index.index, swift_type_stars);
+ logger_index_get_data(&reader->index.index_prev, swift_type_stars);
+ struct index_data *data_created = logger_index_get_created_history(
+ &reader->index.index_next, swift_type_stars);
/* Sort the fields in order to read the correct bits. */
logger_reader_global_to_local(
reader, global_fields_wanted, local_fields_wanted, local_first_deriv,
local_second_deriv, n_fields_wanted, swift_type_stars);
+ size_t current_in_index = 0;
+ int reading_history = 0;
+ const size_t size_index = reader->index.index_prev.nparts[swift_type_stars];
+ const size_t size_history =
+ logger_reader_count_number_new_particles(reader, swift_type_stars);
+
/* Read the particles */
for (size_t i = 0; i < n_part[swift_type_stars]; i++) {
- /* Get the offset */
- size_t offset = data[i].offset;
-
- /* Sort the output into output_single. */
- for (int field = 0; field < n_fields_wanted; field++) {
- const int global = global_fields_wanted[field];
- const int local = local_fields_wanted[field];
- const int first = local_first_deriv[field];
- const int second = local_second_deriv[field];
- void *output_single = output[field] + i * h->masks[global].size;
-
- /* Read the field. */
- logger_reader_read_field(reader, time, reader->time.time_offset,
- interp_type, offset, local, first, second,
- output_single, swift_type_stars);
+ int particle_removed = 1;
+ /* Do it until finding a particle not removed. */
+ while (particle_removed) {
+ /* Should we start to read the history? */
+ if (!reading_history && current_in_index == size_index) {
+ current_in_index = 0;
+ reading_history = 1;
+ }
+
+ /* Check if we still have some particles available. */
+ if (reading_history && current_in_index == size_history) {
+ error_python("The logger was not able to find enough particles.");
+ }
+
+ /* Get the offset */
+ size_t offset = reading_history ? data_created[current_in_index].offset
+ : data[current_in_index].offset;
+
+ /* Sort the output into output_single. */
+ for (int field = 0; field < n_fields_wanted; field++) {
+ const int global = global_fields_wanted[field];
+ const int local = local_fields_wanted[field];
+ const int first = local_first_deriv[field];
+ const int second = local_second_deriv[field];
+ void *output_single = output[field] + i * h->masks[global].size;
+
+ /* Read the field. */
+ particle_removed = logger_reader_read_field(
+ reader, time, reader->time.time_offset, interp_type, offset,
+ local, first, second, output_single, swift_type_stars);
+
+ /* Should we continue to read the fields of this particle? */
+ if (particle_removed) {
+ break;
+ }
+ }
+ current_in_index++;
}
}
}
diff --git a/logger/logger_reader.h b/logger/logger_reader.h
index dc74850188a2ba3b3def6b7adccda3a2734778fd..44bb5e6b62d02916611641f9759a757f4f2f66af 100644
--- a/logger/logger_reader.h
+++ b/logger/logger_reader.h
@@ -67,8 +67,11 @@ struct logger_reader {
char basename[STRING_SIZE];
struct {
- /* Information contained in the index file */
- struct logger_index index;
+ /* Information contained in the previous index file. */
+ struct logger_index index_prev;
+
+ /* Information contained in the next index file. */
+ struct logger_index index_next;
/* Number of index files */
int n_files;
@@ -119,8 +122,9 @@ double logger_reader_get_time_begin(struct logger_reader *reader);
double logger_reader_get_time_end(struct logger_reader *reader);
size_t logger_reader_get_next_offset_from_time(struct logger_reader *reader,
double time);
-const uint64_t *logger_reader_get_number_particles(struct logger_reader *reader,
- int *n_type);
+void logger_reader_get_number_particles(struct logger_reader *reader,
+ uint64_t *n_parts,
+ const int *read_types);
void logger_reader_read_all_particles(struct logger_reader *reader, double time,
enum logger_reader_type interp_type,
diff --git a/logger/logger_tools.c b/logger/logger_tools.c
index 34fb0b5185bd3220ada2f657a0ef4f6cf23b8cd6..bef86f7207d70bf069b6de2db17ba559209a15ee 100644
--- a/logger/logger_tools.c
+++ b/logger/logger_tools.c
@@ -139,7 +139,6 @@ size_t tools_reverse_offset(const struct header *h, void *file_map,
/* first records do not have a previous partner. */
if (prev_offset == cur_offset) return after_current_record;
-
if (prev_offset > cur_offset)
error_python("Unexpected offset: header %lu, current %lu.", prev_offset,
cur_offset);
@@ -187,9 +186,21 @@ size_t tools_check_record_consistency(const struct logger_reader *reader,
size_t mask;
size_t pointed_offset;
+ const size_t mask_special_flag =
+ h->masks[header_get_field_index(h, "SpecialFlags")].mask;
+
/* read mask + offset. */
map = logger_loader_io_read_mask(h, map, &mask, &pointed_offset);
+ /* set offset after current record. */
+ map = (char *)map + header_get_record_size_from_mask(h, mask);
+ const size_t offset_ret = (size_t)((char *)map - (char *)file_init);
+
+ /* If something happened, skip the check. */
+ if (mask & mask_special_flag) {
+ return offset_ret;
+ }
+
/* get absolute offset. */
if (header_is_forward(h))
pointed_offset += offset;
@@ -202,11 +213,7 @@ size_t tools_check_record_consistency(const struct logger_reader *reader,
error_python("Offset are corrupted.");
}
- /* set offset after current record. */
- map = (char *)map + header_get_record_size_from_mask(h, mask);
-
- if (pointed_offset == offset || pointed_offset == 0)
- return (size_t)((char *)map - (char *)file_init);
+ if (pointed_offset == offset || pointed_offset == 0) return offset_ret;
/* read mask of the pointed record. */
size_t pointed_mask = 0;
@@ -219,103 +226,5 @@ size_t tools_check_record_consistency(const struct logger_reader *reader,
error_python("Error in the offset (mask %lu at %lu != %lu at %lu).", mask,
offset, pointed_mask, pointed_offset);
- return (size_t)((char *)map - (char *)file_init);
-}
-
-/**
- * @brief Compute the quintic hermite spline interpolation.
- *
- * @param t0 The time at the left of the interval.
- * @param x0 The function at the left of the interval.
- * @param v0 The first derivative at the left of the interval.
- * @param a0 The second derivative at the left of the interval.
- * @param t1 The time at the right of the interval.
- * @param x1 The function at the right of the interval.
- * @param v1 The first derivative at the right of the interval.
- * @param a1 The second derivative at the right of the interval.
- * @param t The time of the interpolation.
- *
- * @return The function evaluated at t.
- */
-double logger_tools_quintic_hermite_spline(double t0, double x0, float v0,
- float a0, double t1, double x1,
- float v1, float a1, double t) {
-
- /* Generates recurring variables */
- /* Time differences */
- const double dt = t1 - t0;
- const double dt2 = dt * dt;
- const double dt3 = dt2 * dt;
- const double dt4 = dt3 * dt;
- const double dt5 = dt4 * dt;
-
- const double t_t0 = t - t0;
- const double t_t0_2 = t_t0 * t_t0;
- const double t_t0_3 = t_t0_2 * t_t0;
- const double t_t1 = t - t1;
- const double t_t1_2 = t_t1 * t_t1;
-
- /* Derivatives */
- const double v0_dt = v0 * dt;
- const double a0_dt2 = 0.5 * a0 * dt2;
- const double v1_dt = v1 * dt;
- const double a1_dt2 = 0.5 * a1 * dt2;
-
- /* Do the first 3 terms of the hermite spline */
- double x = x0 + v0 * t_t0 + 0.5 * a0 * t_t0_2;
-
- /* Cubic term */
- x += (x1 - x0 - v0_dt - a0_dt2) * t_t0_3 / dt3;
-
- /* Quartic term */
- x += (3. * x0 - 3. * x1 + v1_dt + 2. * v0_dt + a0_dt2) * t_t0_3 * t_t1 / dt4;
-
- /* Quintic term */
- x += (6. * x1 - 6. * x0 - 3. * v0_dt - 3. * v1_dt + a1_dt2 - a0_dt2) *
- t_t0_3 * t_t1_2 / dt5;
-
- return x;
-}
-
-/**
- * @brief Compute the cubic hermite spline interpolation.
- *
- * @param t0 The time at the left of the interval.
- * @param v0 The first derivative at the left of the interval.
- * @param a0 The second derivative at the left of the interval.
- * @param t1 The time at the right of the interval.
- * @param v1 The first derivative at the right of the interval.
- * @param a1 The second derivative at the right of the interval.
- * @param t The time of the interpolation.
- *
- * @return The function evaluated at t.
- */
-float logger_tools_cubic_hermite_spline(double t0, float v0, float a0,
- double t1, float v1, float a1,
- double t) {
-
- /* Generates recurring variables */
- /* Time differences */
- const float dt = t1 - t0;
- const float dt2 = dt * dt;
- const float dt3 = dt2 * dt;
-
- const float t_t0 = t - t0;
- const float t_t0_2 = t_t0 * t_t0;
- const float t_t1 = t - t1;
-
- /* Derivatives */
- const float a0_dt = a0 * dt;
- const float a1_dt = a1 * dt;
-
- /* Do the first 2 terms of the hermite spline */
- float x = v0 + a0 * t_t0;
-
- /* Square term */
- x += (v1 - v0 - a0_dt) * t_t0_2 / dt2;
-
- /* Cubic term */
- x += (2. * v0 - 2. * v1 + a1_dt + a0_dt) * t_t0_2 * t_t1 / dt3;
-
- return x;
+ return offset_ret;
}
diff --git a/logger/logger_tools.h b/logger/logger_tools.h
index 1b5a031019d2cd0dbdf05d46fc56d52615588ec7..1e1037c016d11b7929bd49dd58b5452216510840 100644
--- a/logger/logger_tools.h
+++ b/logger/logger_tools.h
@@ -71,13 +71,6 @@ size_t tools_reverse_offset(const struct header *h, void *map, size_t offset);
size_t tools_check_record_consistency(const struct logger_reader *reader,
size_t offset);
-double logger_tools_quintic_hermite_spline(double t0, double x0, float v0,
- float a0, double t1, double x1,
- float v1, float a1, double t);
-float logger_tools_cubic_hermite_spline(double t0, float v0, float a0,
- double t1, float v1, float a1,
- double t);
-
#ifndef HAVE_PYTHON
#define error_python(s, ...) error(s, ##__VA_ARGS__);
#else
diff --git a/logger/stars/Default/logger_stars.h b/logger/stars/Default/logger_stars.h
index 63ad7d7fffc32a06ce40834cd847d8efb9d91eda..8f35039183cc975e51df091372ea514bfa7a0a69 100644
--- a/logger/stars/Default/logger_stars.h
+++ b/logger/stars/Default/logger_stars.h
@@ -22,6 +22,7 @@
#include "../config.h"
/* local includes */
+#include "logger_interpolation.h"
#include "logger_loader_io.h"
#include "logger_python_tools.h"
#include "stars_logger.h"
@@ -90,78 +91,25 @@ stars_logger_interpolate_field(const double t_before,
}
#endif
- /* Compute the interpolation scaling. */
- const double wa = (t - t_before) / (t_after - t_before);
- const double wb = 1. - wa;
-
switch (field) {
case stars_logger_field_coordinates:
- for (int i = 0; i < 3; 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;
-
- const double *x_aft = (double *)after->field;
- const float *v_aft = (float *)after->first_deriv;
- const float *a_aft = (float *)after->second_deriv;
-
- /* Use quintic hermite spline. */
- if (v_bef && v_aft && a_bef && a_aft) {
- x[i] = logger_tools_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);
- }
- /* Use cubic hermite spline. */
- else if (v_bef && v_aft) {
- x[i] = logger_tools_cubic_hermite_spline(
- t_before, x_bef[i], v_bef[i], t_after, x_aft[i], v_aft[i], t);
- }
- /* Use linear interpolation. */
- else {
- x[i] = wa * x_aft[i] + wb * x_bef[i];
- }
- }
+ interpolate_quintic_double_float_ND(t_before, before, t_after, after,
+ output, t, /* dimension= */ 3);
break;
case stars_logger_field_velocities:
- for (int i = 0; i < 3; i++) {
- float *v = (float *)output;
- const float *v_bef = (float *)before->field;
- const float *a_bef = (float *)before->first_deriv;
-
- const float *v_aft = (float *)after->field;
- const float *a_aft = (float *)after->first_deriv;
-
- /* Use a cubic hermite spline. */
- if (a_bef && a_aft) {
- v[i] = logger_tools_cubic_hermite_spline(
- t_before, v_bef[i], a_bef[i], t_after, v_aft[i], a_aft[i], t);
- }
- /* Use linear interpolation. */
- else {
- v[i] = wa * v_aft[i] + wb * v_bef[i];
- }
- }
+ interpolate_cubic_float_ND(t_before, before, t_after, after, output, t,
+ /* dimension= */ 3);
break;
case stars_logger_field_accelerations:
- /* interpolate vectors. */
- for (int i = 0; i < 3; i++) {
- float *a = (float *)output;
- const float *a_bef = (float *)before->field;
- const float *a_aft = (float *)after->field;
- a[i] = wa * a_aft[i] + wb * a_bef[i];
- }
+ interpolate_linear_float_ND(t_before, before, t_after, after, output, t,
+ /* dimension= */ 3);
break;
case stars_logger_field_smoothing_lengths:
case stars_logger_field_masses:
- ((float *)output)[0] =
- wa * ((float *)after->field)[0] + wb * ((float *)before->field)[0];
+ interpolate_linear_float(t_before, before, t_after, after, output, t);
break;
case stars_logger_field_particle_ids:
- if (*(long long *)after->field != *(long long *)before->field) {
- error("Interpolating different particles");
- }
- *(long long *)output = *(long long *)after->field;
+ interpolate_ids(t_before, before, t_after, after, output, t);
break;
default:
error("Not implemented");
diff --git a/src/Makefile.am b/src/Makefile.am
index 300b844a8240dea2ba8ae746e612d3f0d61a53d6..b9d9ff5cf46a568b538321687ec6844a198aab04 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -174,7 +174,7 @@ nobase_noinst_HEADERS += hydro/AnarchyPU/hydro.h hydro/AnarchyPU/hydro_iact.h hy
nobase_noinst_HEADERS += hydro/AnarchyPU/hydro_debug.h hydro/AnarchyPU/hydro_part.h
nobase_noinst_HEADERS += hydro/AnarchyPU/hydro_parameters.h
nobase_noinst_HEADERS += hydro/SPHENIX/hydro.h hydro/SPHENIX/hydro_iact.h hydro/SPHENIX/hydro_io.h
-nobase_noinst_HEADERS += hydro/SPHENIX/hydro_debug.h hydro/SPHENIX/hydro_part.h
+nobase_noinst_HEADERS += hydro/SPHENIX/hydro_debug.h hydro/SPHENIX/hydro_part.h hydro/SPHENIX/hydro_logger.h
nobase_noinst_HEADERS += hydro/SPHENIX/hydro_parameters.h
nobase_noinst_HEADERS += hydro/Gizmo/hydro_parameters.h
nobase_noinst_HEADERS += hydro/Gizmo/hydro_io.h hydro/Gizmo/hydro_debug.h
diff --git a/src/engine.c b/src/engine.c
index 064a21b6679cb320954225de753670deddd91556..f2e5c3d4c6629a662488012dcbfa742f38c7b104 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1742,8 +1742,13 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
#ifdef WITH_LOGGER
if (e->policy & engine_policy_logger) {
/* Mark the first time step in the particle logger file. */
- logger_log_timestamp(e->logger, e->ti_current, e->time,
- &e->logger->timestamp_offset);
+ if (e->policy & engine_policy_cosmology) {
+ logger_log_timestamp(e->logger, e->ti_current, e->cosmology->a,
+ &e->logger->timestamp_offset);
+ } else {
+ logger_log_timestamp(e->logger, e->ti_current, e->time,
+ &e->logger->timestamp_offset);
+ }
/* Make sure that we have enough space in the particle logger file
* to store the particles in current time step. */
logger_ensure_size(e->logger, s->nr_parts, s->nr_gparts, s->nr_sparts);
@@ -2117,8 +2122,13 @@ void engine_step(struct engine *e) {
#ifdef WITH_LOGGER
if (e->policy & engine_policy_logger) {
/* Mark the current time step in the particle logger file. */
- logger_log_timestamp(e->logger, e->ti_current, e->time,
- &e->logger->timestamp_offset);
+ if (e->policy & engine_policy_cosmology) {
+ logger_log_timestamp(e->logger, e->ti_current, e->cosmology->a,
+ &e->logger->timestamp_offset);
+ } else {
+ logger_log_timestamp(e->logger, e->ti_current, e->time,
+ &e->logger->timestamp_offset);
+ }
/* Make sure that we have enough space in the particle logger file
* to store the particles in current time step. */
logger_ensure_size(e->logger, e->s->nr_parts, e->s->nr_gparts,
diff --git a/src/hydro/Gadget2/hydro_logger.h b/src/hydro/Gadget2/hydro_logger.h
index e7a12a5b0b3defec623aca3101ef879e6fba5ba6..907259cb699e86f7646966349ccd1d0bbcd5cd2b 100644
--- a/src/hydro/Gadget2/hydro_logger.h
+++ b/src/hydro/Gadget2/hydro_logger.h
@@ -19,6 +19,7 @@
#ifndef SWIFT_GADGET2_HYDRO_LOGGER_H
#define SWIFT_GADGET2_HYDRO_LOGGER_H
+#include "hydro.h"
#include "logger_io.h"
#ifdef WITH_LOGGER
@@ -180,9 +181,14 @@ INLINE static char *hydro_logger_write_particle(
/* Compute the acceleration due to hydro and gravity */
float *acc = (float *)buff;
- acc[0] = p->a_hydro[0] + p->gpart->a_grav[0];
- acc[1] = p->a_hydro[1] + p->gpart->a_grav[1];
- acc[2] = p->a_hydro[2] + p->gpart->a_grav[2];
+ acc[0] = p->a_hydro[0];
+ acc[1] = p->a_hydro[1];
+ acc[2] = p->a_hydro[2];
+ if (p->gpart) {
+ acc[0] += p->gpart->a_grav[0];
+ acc[1] += p->gpart->a_grav[1];
+ acc[2] += p->gpart->a_grav[2];
+ }
memcpy(buff, acc, 3 * sizeof(float));
buff += 3 * sizeof(float);
diff --git a/src/hydro/SPHENIX/hydro_logger.h b/src/hydro/SPHENIX/hydro_logger.h
new file mode 100644
index 0000000000000000000000000000000000000000..44f9274e6a6a352aea789d3eaab2f24fc86825b5
--- /dev/null
+++ b/src/hydro/SPHENIX/hydro_logger.h
@@ -0,0 +1,309 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2020 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/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_SPHENIX_HYDRO_LOGGER_H
+#define SWIFT_SPHENIX_HYDRO_LOGGER_H
+
+#include "hydro.h"
+#include "logger_io.h"
+
+#ifdef WITH_LOGGER
+
+/*
+ * List of all possible mask.
+ * Outside the module, only hydro_logger_field_count is used.
+ */
+enum hydro_logger_fields {
+ hydro_logger_field_coordinates = 0,
+ hydro_logger_field_velocities,
+ hydro_logger_field_accelerations,
+ hydro_logger_field_masses,
+ hydro_logger_field_smoothing_lengths,
+ hydro_logger_field_internal_energies,
+ hydro_logger_field_particle_ids,
+ hydro_logger_field_densities,
+ hydro_logger_field_entropies,
+ hydro_logger_field_pressures,
+ hydro_logger_field_viscosity_diffusion,
+ hydro_logger_field_velocity_divergences,
+ hydro_logger_field_count,
+};
+
+/* Name of each possible mask. */
+static const char *hydro_logger_field_names[hydro_logger_field_count] = {
+ "Coordinates", "Velocities", "Accelerations",
+ "Masses", "SmoothingLengths", "InternalEnergies",
+ "ParticleIDs", "Densities", "Entropies",
+ "Pressures", "ViscosityDiffusion", "VelocityDivergences",
+};
+
+/**
+ * @brief Initialize the logger.
+ *
+ * WARNING: The order should be the same in all the functions and
+ * #hydro_logger_fields!
+ *
+ * @param mask_data Data for each type of mask.
+ *
+ * @return Number of masks used.
+ */
+INLINE static int hydro_logger_writer_populate_mask_data(
+ struct mask_data *mask_data) {
+ mask_data[hydro_logger_field_coordinates] = logger_create_mask_entry(
+ hydro_logger_field_names[hydro_logger_field_coordinates],
+ 3 * sizeof(double));
+
+ mask_data[hydro_logger_field_velocities] = logger_create_mask_entry(
+ hydro_logger_field_names[hydro_logger_field_velocities],
+ 3 * sizeof(float));
+
+ mask_data[hydro_logger_field_accelerations] = logger_create_mask_entry(
+ hydro_logger_field_names[hydro_logger_field_accelerations],
+ 3 * sizeof(float));
+
+ mask_data[hydro_logger_field_masses] = logger_create_mask_entry(
+ hydro_logger_field_names[hydro_logger_field_masses], sizeof(float));
+
+ mask_data[hydro_logger_field_smoothing_lengths] = logger_create_mask_entry(
+ hydro_logger_field_names[hydro_logger_field_smoothing_lengths],
+ sizeof(float));
+
+ mask_data[hydro_logger_field_internal_energies] = logger_create_mask_entry(
+ hydro_logger_field_names[hydro_logger_field_internal_energies],
+ sizeof(float));
+
+ mask_data[hydro_logger_field_particle_ids] = logger_create_mask_entry(
+ hydro_logger_field_names[hydro_logger_field_particle_ids],
+ sizeof(long long));
+
+ mask_data[hydro_logger_field_densities] = logger_create_mask_entry(
+ hydro_logger_field_names[hydro_logger_field_densities], sizeof(float));
+
+ mask_data[hydro_logger_field_entropies] = logger_create_mask_entry(
+ hydro_logger_field_names[hydro_logger_field_entropies], sizeof(float));
+
+ mask_data[hydro_logger_field_pressures] = logger_create_mask_entry(
+ hydro_logger_field_names[hydro_logger_field_pressures], sizeof(float));
+
+ mask_data[hydro_logger_field_viscosity_diffusion] = logger_create_mask_entry(
+ hydro_logger_field_names[hydro_logger_field_viscosity_diffusion],
+ 3 * sizeof(float));
+
+ mask_data[hydro_logger_field_velocity_divergences] = logger_create_mask_entry(
+ hydro_logger_field_names[hydro_logger_field_velocity_divergences],
+ 2 * sizeof(float));
+
+ return hydro_logger_field_count;
+}
+
+/**
+ * @brief Generates the mask and compute the size of the record.
+ *
+ * WARNING: The order should be the same in all the functions and
+ * #hydro_logger_fields!
+ *
+ * @param masks The list of masks (same order than in #hydro_logger_init).
+ * @param part The #part that will be written.
+ * @param xpart The #xpart that will be written.
+ * @param write_all Are we forcing to write all the fields?
+ *
+ * @param buffer_size (out) The requested size for the buffer.
+ * @param mask (out) The mask that will be written.
+ */
+INLINE static void hydro_logger_compute_size_and_mask(
+ const struct mask_data *masks, const struct part *part,
+ const struct xpart *xpart, const int write_all, size_t *buffer_size,
+ unsigned int *mask) {
+
+ /* Here you can decide your own writing logic */
+
+ /* Add the coordinates. */
+ *mask |= logger_add_field_to_mask(masks[hydro_logger_field_coordinates],
+ buffer_size);
+
+ /* Add the velocities. */
+ *mask |= logger_add_field_to_mask(masks[hydro_logger_field_velocities],
+ buffer_size);
+
+ /* Add the accelerations. */
+ *mask |= logger_add_field_to_mask(masks[hydro_logger_field_accelerations],
+ buffer_size);
+
+ /* Add the masses. */
+ *mask |=
+ logger_add_field_to_mask(masks[hydro_logger_field_masses], buffer_size);
+
+ /* Add the smoothing lengths. */
+ *mask |= logger_add_field_to_mask(masks[hydro_logger_field_smoothing_lengths],
+ buffer_size);
+
+ /* Add the energies. */
+ *mask |= logger_add_field_to_mask(masks[hydro_logger_field_internal_energies],
+ buffer_size);
+
+ /* Add the ID. */
+ *mask |= logger_add_field_to_mask(masks[hydro_logger_field_particle_ids],
+ buffer_size);
+
+ /* Add the density. */
+ *mask |= logger_add_field_to_mask(masks[hydro_logger_field_densities],
+ buffer_size);
+
+ /* Add the entropies. */
+ *mask |= logger_add_field_to_mask(masks[hydro_logger_field_entropies],
+ buffer_size);
+
+ /* Add the pressures. */
+ *mask |= logger_add_field_to_mask(masks[hydro_logger_field_pressures],
+ buffer_size);
+
+ /* Add the viscosity / diffusion. */
+ *mask |= logger_add_field_to_mask(
+ masks[hydro_logger_field_viscosity_diffusion], buffer_size);
+
+ /* Add the velocity divergences + derivative. */
+ *mask |= logger_add_field_to_mask(
+ masks[hydro_logger_field_velocity_divergences], buffer_size);
+}
+
+/**
+ * @brief Write a particle to the logger.
+ *
+ * WARNING: The order should be the same in all the functions and
+ * #hydro_logger_fields!
+ *
+ * @param masks The list of masks (same order than in #hydro_logger_init).
+ * @param p The #part to write.
+ * @param xp The #xpart to write.
+ * @param mask The mask to use for this record.
+ * @param buff The buffer where to write the particle.
+ *
+ * @return The buffer after the data.
+ */
+INLINE static char *hydro_logger_write_particle(
+ const struct mask_data *mask_data, const struct part *p,
+ const struct xpart *xp, unsigned int *mask, char *buff) {
+
+ /* Write the coordinate. */
+ if (logger_should_write_field(mask_data[hydro_logger_field_coordinates],
+ mask)) {
+ memcpy(buff, p->x, 3 * sizeof(double));
+ buff += 3 * sizeof(double);
+ }
+
+ /* Write the velocity. */
+ if (logger_should_write_field(mask_data[hydro_logger_field_velocities],
+ mask)) {
+ memcpy(buff, p->v, 3 * sizeof(float));
+ buff += 3 * sizeof(float);
+ }
+
+ /* Write the acceleration. */
+ if (logger_should_write_field(mask_data[hydro_logger_field_accelerations],
+ mask)) {
+
+ /* Compute the acceleration due to hydro and gravity */
+ float *acc = (float *)buff;
+ acc[0] = p->a_hydro[0];
+ acc[1] = p->a_hydro[1];
+ acc[2] = p->a_hydro[2];
+ if (p->gpart) {
+ acc[0] += p->gpart->a_grav[0];
+ acc[1] += p->gpart->a_grav[1];
+ acc[2] += p->gpart->a_grav[2];
+ }
+
+ memcpy(buff, acc, 3 * sizeof(float));
+ buff += 3 * sizeof(float);
+ }
+
+ /* Write the mass. */
+ if (logger_should_write_field(mask_data[hydro_logger_field_masses], mask)) {
+ memcpy(buff, &p->mass, sizeof(float));
+ buff += sizeof(float);
+ }
+
+ /* Write the smoothing length. */
+ if (logger_should_write_field(mask_data[hydro_logger_field_smoothing_lengths],
+ mask)) {
+ memcpy(buff, &p->h, sizeof(float));
+ buff += sizeof(float);
+ }
+
+ /* Write the energy. */
+ if (logger_should_write_field(mask_data[hydro_logger_field_internal_energies],
+ mask)) {
+ memcpy(buff, &p->u, sizeof(float));
+ buff += sizeof(float);
+ }
+
+ /* Write the Id. */
+ if (logger_should_write_field(mask_data[hydro_logger_field_particle_ids],
+ mask)) {
+ memcpy(buff, &p->id, sizeof(long long));
+ buff += sizeof(long long);
+ }
+
+ /* Write the density. */
+ if (logger_should_write_field(mask_data[hydro_logger_field_densities],
+ mask)) {
+ memcpy(buff, &p->rho, sizeof(float));
+ buff += sizeof(float);
+ }
+
+ /* Write the entropy. */
+ if (logger_should_write_field(mask_data[hydro_logger_field_entropies],
+ mask)) {
+ const float entropy = hydro_get_comoving_entropy(p, xp);
+ memcpy(buff, &entropy, sizeof(float));
+ buff += sizeof(float);
+ }
+
+ /* Write the pressures. */
+ if (logger_should_write_field(mask_data[hydro_logger_field_pressures],
+ mask)) {
+ const float pressure = hydro_get_comoving_pressure(p);
+ memcpy(buff, &pressure, sizeof(float));
+ buff += sizeof(float);
+ }
+
+ /* Write the viscosity / diffusion. */
+ if (logger_should_write_field(
+ mask_data[hydro_logger_field_viscosity_diffusion], mask)) {
+ const float coef[3] = {p->viscosity.alpha * p->force.balsara,
+ p->diffusion.alpha, p->diffusion.laplace_u};
+ memcpy(buff, coef, sizeof(coef));
+ buff += sizeof(coef);
+ }
+
+ /* Write the velocity divergence. */
+ if (logger_should_write_field(
+ mask_data[hydro_logger_field_velocity_divergences], mask)) {
+ const float div[2] = {
+ p->viscosity.div_v,
+ p->viscosity.div_v_dt,
+ };
+ memcpy(buff, div, sizeof(div));
+ buff += sizeof(div);
+ }
+
+ return buff;
+}
+
+#endif // WITH_LOGGER
+#endif // SWIFT_SPHENIX_HYDRO_LOGGER_H
diff --git a/src/hydro/SPHENIX/hydro_part.h b/src/hydro/SPHENIX/hydro_part.h
index 5904dfadeb1403181a6b6040795b3fc49d9bc433..836d5166c679b08085ba74fc69b39cba165e7d5c 100644
--- a/src/hydro/SPHENIX/hydro_part.h
+++ b/src/hydro/SPHENIX/hydro_part.h
@@ -30,6 +30,7 @@
#include "chemistry_struct.h"
#include "cooling_struct.h"
#include "feedback_struct.h"
+#include "logger.h"
#include "particle_splitting_struct.h"
#include "pressure_floor_struct.h"
#include "star_formation_struct.h"
@@ -75,6 +76,11 @@ struct xpart {
/* Additional data used by the feedback */
struct feedback_part_data feedback_data;
+#ifdef WITH_LOGGER
+ /* Additional data for the particle logger */
+ struct logger_part_data logger_data;
+#endif
+
} SWIFT_STRUCT_ALIGN;
/**
diff --git a/src/hydro_io.h b/src/hydro_io.h
index 60934acb41a7da337efc68b2f2d77de1f35533c5..75bb63cb9b7dba157b2069347046842e733d91d2 100644
--- a/src/hydro_io.h
+++ b/src/hydro_io.h
@@ -27,7 +27,6 @@
#include "./hydro/Minimal/hydro_io.h"
#elif defined(GADGET2_SPH)
#include "./hydro/Gadget2/hydro_io.h"
-#include "./hydro/Gadget2/hydro_logger.h"
#elif defined(HOPKINS_PE_SPH)
#include "./hydro/PressureEntropy/hydro_io.h"
#elif defined(HOPKINS_PU_SPH)
diff --git a/src/hydro_logger.h b/src/hydro_logger.h
index 503a3b4acab6dae349694009943088e4056945b0..a4fef3617ef6c76fc7bdfcee37ad492392a351ca 100644
--- a/src/hydro_logger.h
+++ b/src/hydro_logger.h
@@ -48,7 +48,7 @@
#elif defined(PLANETARY_SPH)
#error TODO
#elif defined(SPHENIX_SPH)
-#error TODO
+#include "./hydro/SPHENIX/hydro_logger.h"
#elif defined(ANARCHY_PU_SPH)
#error TODO
#else
diff --git a/src/logger.c b/src/logger.c
index be462bcfe58bf78c3f571bb1c03648604789cbd6..f9c0b407ac3e5d521f092bb213b17ddbfc74b3d1 100644
--- a/src/logger.c
+++ b/src/logger.c
@@ -38,13 +38,14 @@
#include "logger.h"
/* Local headers. */
+#include "active.h"
#include "atomic.h"
#include "dump.h"
#include "engine.h"
#include "error.h"
#include "gravity_logger.h"
-#include "hydro_io.h"
-#include "stars_io.h"
+#include "hydro_logger.h"
+#include "stars_logger.h"
#include "units.h"
/*
@@ -127,6 +128,7 @@ void logger_write_data(struct dump *d, size_t *offset, size_t size,
/**
* @brief log all particles in the engine.
*
+ * TODO use threadpool + logger function for multiple particles.
* @param log The #logger_writer
* @param e The #engine
*/
@@ -140,18 +142,37 @@ void logger_log_all_particles(struct logger_writer *log,
const struct space *s = e->s;
/* log the parts. */
- logger_log_parts(log, s->parts, s->xparts, s->nr_parts, e,
- /* log_all_fields= */ 1, /* flag= */ 0, /* flag_data= */ 0);
+ for (size_t i = 0; i < s->nr_parts; i++) {
+ struct part *p = &s->parts[i];
+ struct xpart *xp = &s->xparts[i];
+ if (!part_is_inhibited(p, e) && p->time_bin != time_bin_not_created) {
+ logger_log_part(log, p, xp, e,
+ /* log_all_fields */ 1, /* Special flags */ 0,
+ /* data */ 0);
+ }
+ }
/* log the gparts */
- logger_log_gparts(log, s->gparts, s->nr_gparts, e,
- /* log_all_fields= */ 1, /* flag= */ 0,
- /* flag_data= */ 0);
+ for (size_t i = 0; i < s->nr_gparts; i++) {
+ struct gpart *gp = &s->gparts[i];
+ if (!gpart_is_inhibited(gp, e) && gp->time_bin != time_bin_not_created &&
+ (gp->type == swift_type_dark_matter ||
+ gp->type == swift_type_dark_matter_background)) {
+ logger_log_gpart(log, gp, e,
+ /* log_all_fields */ 1, /* Special flags */ 0,
+ /* data */ 0);
+ }
+ }
/* log the parts */
- logger_log_sparts(log, s->sparts, s->nr_sparts, e,
- /* log_all_fields= */ 1, /* flag= */ 0,
- /* flag_data= */ 0);
+ for (size_t i = 0; i < s->nr_sparts; i++) {
+ struct spart *sp = &s->sparts[i];
+ if (!spart_is_inhibited(sp, e) && sp->time_bin != time_bin_not_created) {
+ logger_log_spart(log, sp, e,
+ /* log_all_fields */ 1, /* Special flags */ 0,
+ /* data */ 0);
+ }
+ }
if (e->total_nr_bparts > 0) error("Not implemented");
}
@@ -184,10 +205,6 @@ void logger_copy_part_fields(const struct logger_writer *log,
/* Write the header. */
buff = logger_write_record_header(buff, &mask, offset, offset_new);
- /* Write the hydro fields */
- buff = hydro_logger_write_particle(log->mask_data_pointers.hydro, p, xp,
- &mask, buff);
-
/* Special flags */
if (mask & log->logger_mask_data[logger_index_special_flags].mask) {
memcpy(buff, &special_flags,
@@ -196,6 +213,10 @@ void logger_copy_part_fields(const struct logger_writer *log,
mask &= ~log->logger_mask_data[logger_index_special_flags].mask;
}
+ /* Write the hydro fields */
+ buff = hydro_logger_write_particle(log->mask_data_pointers.hydro, p, xp,
+ &mask, buff);
+
#ifdef SWIFT_DEBUG_CHECKS
if (mask) {
error("Requested logging of values not present in parts. %u", mask);
@@ -243,18 +264,27 @@ void logger_log_parts(struct logger_writer *log, const struct part *p,
const int flag_data) {
/* Build the special flag */
+ const int size_special_flag =
+ log->logger_mask_data[logger_index_special_flags].size;
const uint32_t special_flags = logger_pack_flags_and_data(flag, flag_data);
/* Compute the size of the buffer. */
size_t size_total = 0;
if (log_all_fields) {
- size_total = count * (log->max_size_record_part + logger_header_bytes);
+ size_t size = log->max_size_record_part + logger_header_bytes;
+ if (flag != 0) {
+ size += size_special_flag;
+ }
+ size_total = count * size;
} else {
for (int i = 0; i < count; i++) {
unsigned int mask = 0;
size_t size = 0;
hydro_logger_compute_size_and_mask(log->mask_data_pointers.hydro, &p[i],
&xp[i], log_all_fields, &size, &mask);
+ if (flag != 0) {
+ size += size_special_flag;
+ }
size_total += size + logger_header_bytes;
}
}
@@ -272,8 +302,14 @@ void logger_log_parts(struct logger_writer *log, const struct part *p,
&xp[i], log_all_fields, &size, &mask);
size += logger_header_bytes;
- if (special_flags != 0) {
+ /* Add the special flag. */
+ if (flag != 0) {
mask |= log->logger_mask_data[logger_index_special_flags].mask;
+ size += size_special_flag;
+ /* reset the offset of the previous log */
+ if (flag == logger_flag_create || flag == logger_flag_mpi_enter) {
+ xp[i].logger_data.last_offset = 0;
+ }
}
/* Copy everything into the buffer */
@@ -281,22 +317,21 @@ void logger_log_parts(struct logger_writer *log, const struct part *p,
&xp[i].logger_data.last_offset, offset_new, buff,
special_flags);
- /* Write the particle into the history if needed. */
- if (flag & logger_flag_create || flag & logger_flag_mpi_enter) {
- logger_history_log(&log->history_new[swift_type_gas], p->id,
- xp->logger_data.last_offset);
-
- } else if (flag & logger_flag_change_type || flag & logger_flag_delete ||
- flag & logger_flag_mpi_exit) {
- logger_history_log(&log->history_removed[swift_type_gas], p->id,
- xp->logger_data.last_offset);
- }
-
/* Update the pointers */
xp[i].logger_data.last_offset = offset_new;
xp[i].logger_data.steps_since_last_output = 0;
buff += size;
offset_new += size;
+
+ /* Write the particle into the history if needed. */
+ if (flag == logger_flag_create || flag == logger_flag_mpi_enter) {
+ logger_history_log(&log->history_new[swift_type_gas], p[i].id,
+ xp[i].logger_data.last_offset);
+ } else if (flag == logger_flag_change_type || flag == logger_flag_delete ||
+ flag == logger_flag_mpi_exit) {
+ logger_history_log(&log->history_removed[swift_type_gas], p[i].id,
+ xp[i].logger_data.last_offset);
+ }
}
}
@@ -327,10 +362,6 @@ void logger_copy_spart_fields(const struct logger_writer *log,
/* Write the header. */
buff = logger_write_record_header(buff, &mask, offset, offset_new);
- /* Write the stellar fields */
- buff = stars_logger_write_particle(log->mask_data_pointers.stars, sp, &mask,
- buff);
-
/* Special flags */
if (mask & log->logger_mask_data[logger_index_special_flags].mask) {
memcpy(buff, &special_flags,
@@ -339,6 +370,9 @@ void logger_copy_spart_fields(const struct logger_writer *log,
mask &= ~log->logger_mask_data[logger_index_special_flags].mask;
}
+ /* Write the stellar fields */
+ buff = stars_logger_write_particle(log->mask_data_pointers.stars, sp, &mask,
+ buff);
#ifdef SWIFT_DEBUG_CHECKS
if (mask) {
error("Requested logging of values not present in sparts. %u", mask);
@@ -380,18 +414,27 @@ void logger_log_sparts(struct logger_writer *log, struct spart *sp, int count,
const enum logger_special_flags flag,
const int flag_data) {
/* Build the special flag */
+ const int size_special_flag =
+ log->logger_mask_data[logger_index_special_flags].size;
const uint32_t special_flags = logger_pack_flags_and_data(flag, flag_data);
/* Compute the size of the buffer. */
size_t size_total = 0;
if (log_all_fields) {
- size_total = count * (log->max_size_record_spart + logger_header_bytes);
+ size_t size = log->max_size_record_spart + logger_header_bytes;
+ if (flag != 0) {
+ size += size_special_flag;
+ }
+ size_total = count * size;
} else {
for (int i = 0; i < count; i++) {
unsigned int mask = 0;
size_t size = 0;
stars_logger_compute_size_and_mask(log->mask_data_pointers.stars, &sp[i],
log_all_fields, &size, &mask);
+ if (flag != 0) {
+ size += size_special_flag;
+ }
size_total += size + logger_header_bytes;
}
}
@@ -408,8 +451,15 @@ void logger_log_sparts(struct logger_writer *log, struct spart *sp, int count,
log_all_fields, &size, &mask);
size += logger_header_bytes;
- if (special_flags != 0) {
+ /* Add the special flag. */
+ if (flag != 0) {
mask |= log->logger_mask_data[logger_index_special_flags].mask;
+ size += size_special_flag;
+
+ /* reset the offset of the previous log */
+ if (flag == logger_flag_create || flag == logger_flag_mpi_enter) {
+ sp[i].logger_data.last_offset = 0;
+ }
}
/* Copy everything into the buffer */
@@ -417,22 +467,21 @@ void logger_log_sparts(struct logger_writer *log, struct spart *sp, int count,
&sp[i].logger_data.last_offset, offset_new, buff,
special_flags);
- /* Write the particle into the history if needed. */
- if (flag & logger_flag_create || flag & logger_flag_mpi_enter) {
- logger_history_log(&log->history_new[swift_type_stars], sp->id,
- sp->logger_data.last_offset);
-
- } else if (flag & logger_flag_change_type || flag & logger_flag_delete ||
- flag & logger_flag_mpi_exit) {
- logger_history_log(&log->history_removed[swift_type_stars], sp->id,
- sp->logger_data.last_offset);
- }
-
/* Update the pointers */
sp[i].logger_data.last_offset = offset_new;
sp[i].logger_data.steps_since_last_output = 0;
buff += size;
offset_new += size;
+
+ /* Write the particle into the history if needed. */
+ if (flag == logger_flag_create || flag == logger_flag_mpi_enter) {
+ logger_history_log(&log->history_new[swift_type_stars], sp[i].id,
+ sp[i].logger_data.last_offset);
+ } else if (flag == logger_flag_change_type || flag == logger_flag_delete ||
+ flag == logger_flag_mpi_exit) {
+ logger_history_log(&log->history_removed[swift_type_stars], sp[i].id,
+ sp[i].logger_data.last_offset);
+ }
}
}
@@ -463,10 +512,6 @@ void logger_copy_gpart_fields(const struct logger_writer *log,
/* Write the header. */
buff = logger_write_record_header(buff, &mask, offset, offset_new);
- /* Write the hydro fields */
- buff = gravity_logger_write_particle(log->mask_data_pointers.gravity, gp,
- &mask, buff);
-
/* Special flags */
if (mask & log->logger_mask_data[logger_index_special_flags].mask) {
memcpy(buff, &special_flags,
@@ -475,6 +520,10 @@ void logger_copy_gpart_fields(const struct logger_writer *log,
mask &= ~log->logger_mask_data[logger_index_special_flags].mask;
}
+ /* Write the gravity fields */
+ buff = gravity_logger_write_particle(log->mask_data_pointers.gravity, gp,
+ &mask, buff);
+
#ifdef SWIFT_DEBUG_CHECKS
if (mask) {
error("Requested logging of values not present in gparts. %u", mask);
@@ -515,21 +564,32 @@ void logger_log_gparts(struct logger_writer *log, struct gpart *p, int count,
const enum logger_special_flags flag,
const int flag_data) {
/* Build the special flag */
+ const int size_special_flag =
+ log->logger_mask_data[logger_index_special_flags].size;
const uint32_t special_flags = logger_pack_flags_and_data(flag, flag_data);
/* Compute the size of the buffer. */
size_t size_total = 0;
if (log_all_fields) {
- size_total = count * (log->max_size_record_gpart + logger_header_bytes);
+ size_t size = log->max_size_record_gpart + logger_header_bytes;
+ 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) continue;
+ 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;
}
}
@@ -540,7 +600,9 @@ void logger_log_gparts(struct logger_writer *log, struct gpart *p, int count,
for (int i = 0; i < count; i++) {
/* Log only the dark matter */
- if (p[i].type != swift_type_dark_matter) continue;
+ if (p[i].type != swift_type_dark_matter &&
+ p[i].type != swift_type_dark_matter_background)
+ continue;
/* Get the masks */
size_t size = 0;
@@ -549,30 +611,36 @@ void logger_log_gparts(struct logger_writer *log, struct gpart *p, int count,
log_all_fields, &size, &mask);
size += logger_header_bytes;
- if (special_flags != 0) {
+ /* Add the special flag. */
+ if (flag != 0) {
mask |= log->logger_mask_data[logger_index_special_flags].mask;
+ size += size_special_flag;
+
+ /* reset the offset of the previous log */
+ if (flag == logger_flag_create || flag == logger_flag_mpi_enter) {
+ p[i].logger_data.last_offset = 0;
+ }
}
/* Copy everything into the buffer */
logger_copy_gpart_fields(log, &p[i], e, mask, &p[i].logger_data.last_offset,
offset_new, buff, special_flags);
- /* Write the particle into the history if needed. */
- if (flag & logger_flag_create || flag & logger_flag_mpi_enter) {
- logger_history_log(&log->history_new[swift_type_dark_matter],
- p->id_or_neg_offset, p->logger_data.last_offset);
-
- } else if (flag & logger_flag_change_type || flag & logger_flag_delete ||
- flag & logger_flag_mpi_exit) {
- logger_history_log(&log->history_removed[swift_type_dark_matter],
- p->id_or_neg_offset, p->logger_data.last_offset);
- }
-
/* Update the pointers */
p[i].logger_data.last_offset = offset_new;
p[i].logger_data.steps_since_last_output = 0;
buff += size;
offset_new += size;
+
+ /* Write the particle into the history if needed. */
+ if (flag == logger_flag_create || flag == logger_flag_mpi_enter) {
+ logger_history_log(&log->history_new[swift_type_dark_matter],
+ p[i].id_or_neg_offset, p[i].logger_data.last_offset);
+ } else if (flag == logger_flag_change_type || flag == logger_flag_delete ||
+ flag == logger_flag_mpi_exit) {
+ logger_history_log(&log->history_removed[swift_type_dark_matter],
+ p[i].id_or_neg_offset, p[i].logger_data.last_offset);
+ }
}
}
@@ -873,7 +941,7 @@ void logger_init(struct logger_writer *log, const struct engine *e,
const float max_memory_size =
parser_get_opt_param_float(params, "Logger:maximal_memory_size", 1.);
log->maximal_size_history =
- max_memory_size / sizeof(struct logger_index_data);
+ 1e9 * max_memory_size / sizeof(struct logger_index_data);
if (e->nodeID == 0) {
message("Maximal memory size for the logger history: %g GB",
diff --git a/src/logger.h b/src/logger.h
index c2810219f7f8dc4ba695ab1df6edf433a442be65..646fc901842176882761292d88fa000a01553fad 100644
--- a/src/logger.h
+++ b/src/logger.h
@@ -44,6 +44,15 @@ struct engine;
/* Size of the strings. */
#define logger_string_length 200
+/*
+ * The two following defines need to correspond to the list's order
+ * in logger_init_masks.
+ */
+/* Index of the special flags in the list of masks */
+#define logger_index_special_flags 0
+/* Index of the timestamp in the list of masks */
+#define logger_index_timestamp 1
+
/**
* Logger entries contain messages representing the particle data at a given
* point in time during the simulation.
diff --git a/src/logger_io.c b/src/logger_io.c
index 178df2a9b7384c8e64f789368390ffdf4a005a02..6e444aeae904a6511244fcdc4f55760a044471c8 100644
--- a/src/logger_io.c
+++ b/src/logger_io.c
@@ -167,6 +167,8 @@ void logger_write_history(struct logger_history* history, struct engine* e,
*/
void logger_write_index_file(struct logger_writer* log, struct engine* e) {
+ const int with_DM_background = e->s->with_DM_background;
+
struct part* parts = e->s->parts;
struct xpart* xparts = e->s->xparts;
struct gpart* gparts = e->s->gparts;
@@ -178,6 +180,15 @@ void logger_write_index_file(struct logger_writer* log, struct engine* e) {
const size_t Nstars = e->s->nr_sparts;
/* const size_t Nblackholes = e->s->nr_bparts; */
+ if (e->s->nr_sinks != 0) {
+ error("Sink particles are not implemented");
+ }
+
+ size_t Ndm_background = 0;
+ if (with_DM_background) {
+ Ndm_background = io_count_dm_background_gparts(gparts, Ntot);
+ }
+
/* Number of particles that we will write */
const size_t Ntot_written =
e->s->nr_gparts - e->s->nr_inhibited_gparts - e->s->nr_extra_gparts;
@@ -190,11 +201,12 @@ void logger_write_index_file(struct logger_writer* log, struct engine* e) {
const size_t Nbaryons_written =
Ngas_written + Nstars_written + Nblackholes_written;
const size_t Ndm_written =
- Ntot_written > 0 ? Ntot_written - Nbaryons_written : 0;
+ Ntot_written > 0 ? Ntot_written - Nbaryons_written - Ndm_background : 0;
/* Format things in a Gadget-friendly array */
uint64_t N_total[swift_type_count] = {
- (uint64_t)Ngas_written, (uint64_t)Ndm_written, 0, 0,
+ (uint64_t)Ngas_written, (uint64_t)Ndm_written,
+ (uint64_t)Ndm_background, 0,
(uint64_t)Nstars_written, (uint64_t)Nblackholes_written};
/* File name */
@@ -212,7 +224,11 @@ void logger_write_index_file(struct logger_writer* log, struct engine* e) {
}
/* Write double time */
- fwrite(&e->time, sizeof(double), 1, f);
+ if (e->policy & engine_policy_cosmology) {
+ fwrite(&e->cosmology->a, sizeof(double), 1, f);
+ } else {
+ fwrite(&e->time, sizeof(double), 1, f);
+ }
/* Write integer time */
fwrite(&e->ti_current, sizeof(integertime_t), 1, f);
@@ -314,6 +330,27 @@ void logger_write_index_file(struct logger_writer* log, struct engine* e) {
}
break;
+ case swift_type_dark_matter_background: {
+
+ /* Ok, we need to fish out the particles we want */
+ N = Ndm_background;
+
+ /* Allocate temporary array */
+ if (swift_memalign("gparts_written", (void**)&gparts_written,
+ gpart_align,
+ Ndm_background * sizeof(struct gpart)) != 0)
+ error("Error while allocating temporary memory for gparts");
+
+ /* Collect the non-inhibited DM particles from gpart */
+ const int with_stf = 0;
+ io_collect_gparts_background_to_write(
+ gparts, e->s->gpart_group_data, gparts_written,
+ gpart_group_data_written, Ntot, Ndm_background, with_stf);
+
+ /* Select the fields to write */
+ num_fields += darkmatter_write_index(gparts_written, list);
+ } break;
+
case swift_type_stars:
if (Nstars == Nstars_written) {
diff --git a/src/runner_others.c b/src/runner_others.c
index 479da517a7e435d91161ba4397d7acda9dbd76fe..e0420324d5b9c2ccb9505918bc6667cfbe66c821 100644
--- a/src/runner_others.c
+++ b/src/runner_others.c
@@ -278,14 +278,6 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
if (star_formation_should_convert_to_star(p, xp, sf_props, e,
dt_star)) {
-#ifdef WITH_LOGGER
- /* Write the particle */
- /* Logs all the fields request by the user */
- // TODO select only the requested fields
- logger_log_part(e->logger, p, xp, e, /* log_all */ 1,
- logger_flag_change_type, swift_type_stars);
-#endif
-
/* Convert the gas particle to a star particle */
struct spart *sp = NULL;
const int spawn_spart =
@@ -301,6 +293,13 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
} else {
/* Convert the gas particle to a star particle */
sp = cell_convert_part_to_spart(e, c, p, xp);
+#ifdef WITH_LOGGER
+ /* Write the particle */
+ /* Logs all the fields request by the user */
+ // TODO select only the requested fields
+ logger_log_part(e->logger, p, xp, e, /* log_all */ 1,
+ logger_flag_change_type, swift_type_stars);
+#endif
}
} else {
@@ -343,8 +342,13 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
}
#ifdef WITH_LOGGER
- /* Copy the properties back to the stellar particle */
- sp->logger_data = xp->logger_data;
+ if (spawn_spart) {
+ /* Set to zero the logger data. */
+ logger_part_data_init(&sp->logger_data);
+ } else {
+ /* Copy the properties back to the stellar particle */
+ sp->logger_data = xp->logger_data;
+ }
/* Write the s-particle */
logger_log_spart(e->logger, sp, e, /* log_all */ 1,
@@ -769,7 +773,9 @@ void runner_do_logger(struct runner *r, struct cell *c, int timer) {
struct gpart *restrict gp = &gparts[k];
/* Write only the dark matter particles */
- if (gp->type != swift_type_dark_matter) continue;
+ if (gp->type != swift_type_dark_matter &&
+ gp->type != swift_type_dark_matter_background)
+ continue;
/* If particle needs to be log */
if (gpart_is_active(gp, e)) {
diff --git a/src/stars_io.h b/src/stars_io.h
index 6cf92eac91c2eca09070864a8410ac720cc2ee5b..18224b1ef9a189c719d3674a037eb4b26cd14d4e 100644
--- a/src/stars_io.h
+++ b/src/stars_io.h
@@ -27,7 +27,6 @@
#include "./stars/const/stars_io.h"
#elif defined(STARS_NONE)
#include "./stars/Default/stars_io.h"
-#include "./stars/Default/stars_logger.h"
#elif defined(STARS_EAGLE)
#include "./stars/EAGLE/stars_io.h"
#elif defined(STARS_GEAR)
diff --git a/tools/read_index_file.py b/tools/read_index_file.py
index 96ab1df2e7e8426b5c84e3822b8fb6a05aabfd46..51fdec59208e10e62ba5293f23e72674866cffe1 100644
--- a/tools/read_index_file.py
+++ b/tools/read_index_file.py
@@ -51,7 +51,6 @@ with open(filename, "rb") as f:
continue
data = np.fromfile(f, dtype=dt, count=n)
-
print("\t", data)
# print the history of particles removed
@@ -63,5 +62,4 @@ with open(filename, "rb") as f:
continue
data = np.fromfile(f, dtype=dt, count=n)
-
print("\t", data)