diff --git a/src/csds_cosmology.c b/src/csds_cosmology.c
index 6c99f48c04974d8e344f3a3981d31b9056aec648..6f1a4e167cc371bf3a2d502b3c6451968a27f525 100644
--- a/src/csds_cosmology.c
+++ b/src/csds_cosmology.c
@@ -26,131 +26,6 @@
/* Include SWIFT */
#include "parser.h"
-/* Some standard headers */
-#include <math.h>
-
-/**
- * @brief Add the cosmological factors for the position.
- * See Hausammann et al. 2021 for details.
- *
- * @param cosmo The cosmological model.
- * @param scale The scale factor.
- * @param vel The velocity to modify.
- * @param acc The acceleration to modify.
- */
-void csds_cosmology_add_factor_position(const struct csds_cosmology *cosmo,
- const double scale, float *vel,
- float *acc) {
-
- const double a_dot = csds_cosmology_get_scale_factor_derivative(cosmo, scale);
- const double scale2 = scale * scale;
- if (acc) {
- const double a_ddot =
- csds_cosmology_get_scale_factor_second_derivative(cosmo, scale);
- *acc = (*acc - 2 * a_dot * *vel) * scale;
- *acc -= scale2 * a_ddot * *vel / a_dot;
- *acc /= scale2 * scale2 * a_dot * a_dot;
- }
- *vel /= scale2 * a_dot;
-}
-
-/**
- * @brief Add the cosmological factors for the velocities.
- * \f$ \frac{dv'}{da} = \frac{g'}{\dot{a}a} \f$.
- *
- * @param cosmo The cosmological model.
- * @param scale The scale factor.
- * @param acc The acceleration to modify.
- */
-void csds_cosmology_add_factor_velocity(const struct csds_cosmology *cosmo,
- const double scale, float *acc) {
-
- const double a_dot = csds_cosmology_get_scale_factor_derivative(cosmo, scale);
- *acc /= scale * a_dot;
-}
-
-/**
- * @brief Compute the dark energy equation of state.
- *
- * @param cosmo The cosmological model
- * @param scale The scale factor
- *
- * @return The equation of state.
- */
-double csds_cosmology_get_w_tilde(const struct csds_cosmology *cosmo,
- const double scale) {
- const double w_tilde =
- (scale - 1.) * cosmo->w_a - (1. + cosmo->w_0 + cosmo->w_a) * log(scale);
- return exp(3. * w_tilde);
-}
-
-/**
- * @brief Compute E(a) for the derivatives of the scale factor.
- *
- * @param cosmo The cosmological model
- * @param scale The scale factor
- *
- * @return The first derivative of the scale factor.
- */
-double csds_cosmology_get_E(const struct csds_cosmology *cosmo,
- const double scale) {
- const double a_inv = 1. / scale;
- const double a_inv2 = a_inv * a_inv;
- const double Omega_m = cosmo->Omega_b + cosmo->Omega_cdm;
- const double dark_eos = csds_cosmology_get_w_tilde(cosmo, scale);
- const double E = sqrt(cosmo->Omega_r * a_inv2 * a_inv2 + /* Radiation */
- Omega_m * a_inv * a_inv2 + /* Matter */
- cosmo->Omega_k * a_inv2 + /* Curvature */
- cosmo->Omega_lambda * dark_eos);
- return E;
-}
-
-/**
- * @brief Compute the first time derivative of the scale factor
- *
- * @param cosmo The cosmological model
- * @param scale The scale factor
- *
- * @return The first derivative of the scale factor.
- */
-double csds_cosmology_get_scale_factor_derivative(
- const struct csds_cosmology *cosmo, const double scale) {
- const double E = csds_cosmology_get_E(cosmo, scale);
- return cosmo->H0 * scale * E;
-}
-
-/**
- * @brief Compute the first time derivative of the scale factor
- *
- * @param cosmo The cosmological model
- * @param scale The scale factor
- *
- * @return The first derivative of the scale factor.
- */
-double csds_cosmology_get_scale_factor_second_derivative(
- const struct csds_cosmology *cosmo, const double scale) {
-
- /* Compute some powers of a */
- const double a_inv = 1. / scale;
- const double a_inv2 = a_inv * a_inv;
- const double a_inv4 = a_inv2 * a_inv2;
-
- /* Compute E */
- const double E = csds_cosmology_get_E(cosmo, scale);
-
- /* Compute the derivative of E */
- double dark_eos_der = (cosmo->w_a) - (1 + cosmo->w_0 + cosmo->w_a) * a_inv;
- dark_eos_der *= csds_cosmology_get_w_tilde(cosmo, scale);
- const double Omega_m = cosmo->Omega_b + cosmo->Omega_cdm;
- const double E_der =
- (-4 * cosmo->Omega_r * a_inv4 * a_inv - 3 * Omega_m * a_inv4 -
- 2 * cosmo->Omega_k * a_inv2 * a_inv + dark_eos_der);
-
- /* Compute the full expression */
- const double a_dot = csds_cosmology_get_scale_factor_derivative(cosmo, scale);
- return cosmo->H0 * (a_dot * E + 0.5 * scale * E_der / E);
-}
-
/**
* @brief Compute the first time derivative of the scale factor
*
diff --git a/src/csds_cosmology.h b/src/csds_cosmology.h
index 03afcd648b7644564f2db29351a53dc93bbed96e..7fe7be47a56dfe3f16d5c5b1b21883bc7466c953 100644
--- a/src/csds_cosmology.h
+++ b/src/csds_cosmology.h
@@ -25,6 +25,12 @@
#ifndef CSDS_COSMOLOGY_H
#define CSDS_COSMOLOGY_H
+/* SWIFT includes */
+#include "inline.h"
+
+/* Some standard headers */
+#include <math.h>
+
struct swift_params;
/**
@@ -57,15 +63,134 @@ struct csds_cosmology {
double H0;
};
-double csds_cosmology_get_scale_factor_derivative(
- const struct csds_cosmology *cosmo, const double scale);
-double csds_cosmology_get_scale_factor_second_derivative(
- const struct csds_cosmology *cosmo, const double scale);
+/**
+ * @brief Compute the dark energy equation of state.
+ *
+ * @param cosmo The cosmological model
+ * @param scale The scale factor
+ *
+ * @return The equation of state.
+ */
+__attribute__((always_inline)) INLINE static double csds_cosmology_get_w_tilde(
+ const struct csds_cosmology *cosmo, const double scale) {
+ const double w_tilde =
+ (scale - 1.) * cosmo->w_a - (1. + cosmo->w_0 + cosmo->w_a) * log(scale);
+ return exp(3. * w_tilde);
+}
+
+/**
+ * @brief Compute E(a) for the derivatives of the scale factor.
+ *
+ * @param cosmo The cosmological model
+ * @param scale The scale factor
+ *
+ * @return The first derivative of the scale factor.
+ */
+__attribute__((always_inline)) INLINE static double csds_cosmology_get_E(
+ const struct csds_cosmology *cosmo, const double scale) {
+
+ const double a_inv = 1. / scale;
+ const double a_inv2 = a_inv * a_inv;
+ const double Omega_m = cosmo->Omega_b + cosmo->Omega_cdm;
+ const double dark_eos = csds_cosmology_get_w_tilde(cosmo, scale);
+ const double E = sqrt(cosmo->Omega_r * a_inv2 * a_inv2 + /* Radiation */
+ Omega_m * a_inv * a_inv2 + /* Matter */
+ cosmo->Omega_k * a_inv2 + /* Curvature */
+ cosmo->Omega_lambda * dark_eos);
+ return E;
+}
+
+/**
+ * @brief Compute the first time derivative of the scale factor
+ *
+ * @param cosmo The cosmological model
+ * @param scale The scale factor
+ *
+ * @return The first derivative of the scale factor.
+ */
+__attribute__((always_inline)) INLINE static double
+csds_cosmology_get_scale_factor_derivative(const struct csds_cosmology *cosmo,
+ const double scale) {
+ const double E = csds_cosmology_get_E(cosmo, scale);
+ return cosmo->H0 * scale * E;
+}
+
+/**
+ * @brief Compute the first time derivative of the scale factor
+ *
+ * @param cosmo The cosmological model
+ * @param scale The scale factor
+ *
+ * @return The first derivative of the scale factor.
+ */
+__attribute__((always_inline)) INLINE static double
+csds_cosmology_get_scale_factor_second_derivative(
+ const struct csds_cosmology *cosmo, const double scale) {
+
+ /* Compute some powers of a */
+ const double a_inv = 1. / scale;
+ const double a_inv2 = a_inv * a_inv;
+ const double a_inv4 = a_inv2 * a_inv2;
+
+ /* Compute E */
+ const double E = csds_cosmology_get_E(cosmo, scale);
+
+ /* Compute the derivative of E */
+ double dark_eos_der = (cosmo->w_a) - (1 + cosmo->w_0 + cosmo->w_a) * a_inv;
+ dark_eos_der *= csds_cosmology_get_w_tilde(cosmo, scale);
+ const double Omega_m = cosmo->Omega_b + cosmo->Omega_cdm;
+ const double E_der =
+ (-4 * cosmo->Omega_r * a_inv4 * a_inv - 3 * Omega_m * a_inv4 -
+ 2 * cosmo->Omega_k * a_inv2 * a_inv + dark_eos_der);
+
+ /* Compute the full expression */
+ const double a_dot = csds_cosmology_get_scale_factor_derivative(cosmo, scale);
+ return cosmo->H0 * (a_dot * E + 0.5 * scale * E_der / E);
+}
+
+/**
+ * @brief Add the cosmological factors for the position.
+ * See Hausammann et al. 2021 for details.
+ *
+ * @param cosmo The cosmological model.
+ * @param scale The scale factor.
+ * @param vel The velocity to modify.
+ * @param acc The acceleration to modify.
+ */
+__attribute__((always_inline)) INLINE static void
+csds_cosmology_add_factor_position(const struct csds_cosmology *cosmo,
+ const double scale, float *vel, float *acc) {
+
+ const double a_dot = csds_cosmology_get_scale_factor_derivative(cosmo, scale);
+ const double scale2 = scale * scale;
+ if (acc) {
+ const double a_dot_inv = 1. / a_dot;
+ const double scale2_inv = 1. / scale2;
+ const double a_ddot =
+ csds_cosmology_get_scale_factor_second_derivative(cosmo, scale);
+ *acc = (*acc - 2 * a_dot * *vel) * scale;
+ *acc -= scale2 * a_ddot * *vel * a_dot_inv;
+ *acc *= scale2_inv * scale2_inv * a_dot_inv * a_dot_inv;
+ }
+ *vel /= scale2 * a_dot;
+}
+
+/**
+ * @brief Add the cosmological factors for the velocities.
+ * \f$ \frac{dv'}{da} = \frac{g'}{\dot{a}a} \f$.
+ *
+ * @param cosmo The cosmological model.
+ * @param scale The scale factor.
+ * @param acc The acceleration to modify.
+ */
+__attribute__((always_inline)) INLINE static void
+csds_cosmology_add_factor_velocity(const struct csds_cosmology *cosmo,
+ const double scale, float *acc) {
+
+ const double a_dot = csds_cosmology_get_scale_factor_derivative(cosmo, scale);
+ *acc /= scale * a_dot;
+}
+
void csds_cosmology_init(struct csds_cosmology *cosmo,
struct swift_params *params);
-void csds_cosmology_add_factor_position(const struct csds_cosmology *cosmo,
- const double scale, float *vel,
- float *acc);
-void csds_cosmology_add_factor_velocity(const struct csds_cosmology *cosmo,
- const double scale, float *acc);
#endif // CSDS_CSDS_PARAMETERS_H
diff --git a/src/csds_loader_io.h b/src/csds_loader_io.h
index 55ec7ea984d8538232573d2e0b6ffcdbac27ed5a..666c9901581b3fef99207551dcf49de7aaefc51e 100644
--- a/src/csds_loader_io.h
+++ b/src/csds_loader_io.h
@@ -26,8 +26,10 @@
#include "csds_header.h"
#include "csds_tools.h"
+/* Standard headers */
#include <stdio.h>
#include <stdlib.h>
+#include <sys/mman.h>
#define CSDS_VERBOSE_TIMERS -1
@@ -106,4 +108,37 @@ __attribute__((always_inline)) INLINE static char *csds_loader_io_write_data(
return data + size;
};
+/* Defines a few advice for the mmap */
+#define CSDS_ADVICE_SEQUENTIAL(log) \
+ { \
+ /* Warn the OS that we will read in a sequential way */ \
+ int test_ret = madvise(log.map, log.mmap_size, MADV_SEQUENTIAL); \
+ if (test_ret != 0) { \
+ error("Failed to advise the mmap"); \
+ } \
+ }
+
+#define CSDS_ADVICE_NORMAL(log) \
+ { \
+ /* Warn the OS that we will read in a normal way */ \
+ int test_ret = madvise(log.map, log.mmap_size, MADV_NORMAL); \
+ if (test_ret != 0) { \
+ error("Failed to advise the mmap"); \
+ } \
+ }
+
+#define CSDS_ADVICE_WILL_NEED(mmap, addr, length) \
+ { \
+ /* Adjust the size to be at least the page size. */ \
+ const size_t page_size = sysconf(_SC_PAGE_SIZE); \
+ const size_t mem_align = (addr + (page_size - 1)) & ~(page_size - 1); \
+ \
+ /* Warn the OS that we will read soon this address */ \
+ int test_ret = \
+ madvise(mmap + mem_align, length + page_size, MADV_WILLNEED); \
+ if (test_ret != 0) { \
+ error("Failed to advise the mmap"); \
+ } \
+ }
+
#endif // CSDS_CSDS_LOADER_IO_H
diff --git a/src/csds_reader.c b/src/csds_reader.c
index cf5efab39ab7aca6c822c14e4519fbc8b562fd17..41e3fd0b0895a98efa210396d82621f87dc8ed5d 100644
--- a/src/csds_reader.c
+++ b/src/csds_reader.c
@@ -104,6 +104,13 @@ void csds_reader_init_index(struct csds_reader *reader, int number_index,
}
}
+ /* Ensures that we do not have an unexpected situation */
+ if (number_index == 1 || (count == 1 && !restart)) {
+ error(
+ "The CSDS cannot run with only one index file."
+ " Did you forget to restart their generation?");
+ }
+
/* Check if need to generate the index files or not. */
if ((count != number_index && restart) || count == 0) {
csds_reader_generate_index_files(reader, number_index, count);
@@ -168,10 +175,15 @@ void csds_reader_set_time(struct csds_reader *reader, double time) {
/* Set the time */
reader->time.time = time;
- /* Find the correct index */
unsigned int left = 0;
unsigned int right = reader->index.n_files - 1;
+ /* Check if the time is meaningful */
+ if (reader->index.times[left] > time || reader->index.times[right] < time)
+ error("The requested time (%f) is not within the index limits [%f, %f]",
+ time, reader->index.times[left], reader->index.times[right]);
+
+ /* Find the correct index */
while (left != right) {
/* Do a ceil - division */
unsigned int m = (left + right + 1) / 2;
diff --git a/src/csds_reader_generate_index.c b/src/csds_reader_generate_index.c
index 854ad91cdaf80279445d5469c4e8cd9536d4a45c..7030be3ae25da9b1e3e694758a246de8321118f9 100644
--- a/src/csds_reader_generate_index.c
+++ b/src/csds_reader_generate_index.c
@@ -296,12 +296,14 @@ size_t csds_reader_get_initial_state(const struct csds_reader *reader,
struct csds_hashmap **current_state,
struct time_record *time_record) {
const ticks tic = getticks();
-
/* Get a few variables. */
const struct csds_logfile *log = &reader->log;
const struct header *h = &log->header;
const int size_record_header = CSDS_MASK_SIZE + CSDS_OFFSET_SIZE;
+ /* Warn the OS that we will read in a sequential way */
+ CSDS_ADVICE_SEQUENTIAL(log->log);
+
/* Get the offset after the dump of all the particles
and the time information of the first time record. */
if (log->times.size < 2) {
@@ -370,6 +372,9 @@ size_t csds_reader_get_initial_state(const struct csds_reader *reader,
message("took %.3f %s", clocks_from_ticks(getticks() - tic),
clocks_getunit());
+ /* Go back to normal */
+ CSDS_ADVICE_NORMAL(log->log);
+
return offset_max;
}
@@ -440,8 +445,11 @@ void csds_reader_update_particles_to_next_index_mapper(void *map_data,
csds_loader_io_read_mask(h, log->log.map + current_offset, &full_mask,
&h_offset);
- /* Remove the special mask */
- full_mask = full_mask & !SPECIAL_FLAGS_MASK;
+ /* Ensures that a special flag is present in the mask */
+ full_mask |= SPECIAL_FLAGS_MASK;
+
+ /* Now remove it */
+ full_mask = full_mask ^ SPECIAL_FLAGS_MASK;
/* Find the last offset before the current time */
size_t last_full_offset = current_offset;
@@ -536,6 +544,9 @@ size_t csds_reader_update_state_to_next_index(
/* Record the initial time */
const ticks tic = getticks();
+ /* Warn the OS that we will read in a sequential way */
+ CSDS_ADVICE_SEQUENTIAL(log->log);
+
while (offset < time_record.offset) {
/* Print status */
@@ -622,6 +633,9 @@ size_t csds_reader_update_state_to_next_index(
}
}
+ /* Go back to normal */
+ CSDS_ADVICE_NORMAL(log->log);
+
/* Cleanup output */
if (reader->verbose > 0) {
printf("\n");