Merge branch 'madvice' into 'master'

Madvice

See merge request !19

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
  *

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

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

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;

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");