From 01a6209f4829f13cf28319604493116f936ff7dc Mon Sep 17 00:00:00 2001
From: Loic Hausammann <loic.hausammann@protonmail.ch>
Date: Wed, 14 Jul 2021 17:24:02 +0000
Subject: [PATCH] Add test for interpolation

---
 Jenkinsfile                  |   8 +-
 configure.ac                 |   5 +
 coverage.sh                  |   5 +-
 src/csds_cosmology.h         |  12 +-
 src/csds_error.h             |   2 +-
 src/csds_fields.h            |  18 +--
 src/csds_hashmap.h           |   6 +-
 src/csds_header.h            |   6 +-
 src/csds_inline.h            |  16 ++-
 src/csds_interpolation.h     |  14 +-
 src/csds_loader_io.h         |   6 +-
 src/csds_logfile_writer.h    |   6 +-
 src/csds_openmp.h            |   2 +-
 src/csds_particle.h          |   6 +-
 src/csds_python_tools.h      |   8 +-
 src/csds_python_wrapper.c    |   6 +-
 tests/Makefile.am            |   5 +-
 tests/testInterpolation.c    | 245 +++++++++++++++++++++++++++++++++++
 tests/test_interpolation.yml |  66 ++++++++++
 19 files changed, 379 insertions(+), 63 deletions(-)
 create mode 100644 tests/testInterpolation.c
 create mode 100644 tests/test_interpolation.yml

diff --git a/Jenkinsfile b/Jenkinsfile
index 82547eb..b2a4c84 100644
--- a/Jenkinsfile
+++ b/Jenkinsfile
@@ -67,6 +67,7 @@ pipeline {
 
              module purge
              module load gnu_comp/7.3.0
+             module load python/3.8.7-C8  # For coverage
              module load llvm/10.0.1
              export CLANG_FORMAT_CMD="clang-format"
 
@@ -98,12 +99,7 @@ pipeline {
              export CLANG_FORMAT_CMD="clang-format"
 
              git clean -fdx
-             ./autogen.sh
-
-             ./configure --enable-debugging-checks --enable-coverage
-             make -j 2
-
-             ./coverage.sh -r -i
+             ./coverage.sh -c -r -i
              make clean
              '''
 
diff --git a/configure.ac b/configure.ac
index ba1ec7e..2e73d2b 100644
--- a/configure.ac
+++ b/configure.ac
@@ -180,9 +180,14 @@ AC_ARG_ENABLE([coverage],
 if test "$enable_coverage" = "yes"; then
    if test "$ax_cv_c_compiler_vendor" = "gnu"; then
       CFLAGS="$CFLAGS -fprofile-arcs -ftest-coverage"
+      # Remove inline in order to profile headers
+      CFLAGS="$CFLAGS -fno-inline -fno-inline-small-functions -fno-default-inline"
+      # unused function is an issue without inlining
+      CFLAGS="$CFLAGS -Wno-unused-function"
    else
       AC_MSG_ERROR(Cannot use coverage without gcc)
    fi
+   AC_DEFINE([CSDS_NO_INLINE], 1, [Disable forced inlining])
 fi
 
 # Autoconf stuff.
diff --git a/coverage.sh b/coverage.sh
index 4326fb4..8f66d7b 100755
--- a/coverage.sh
+++ b/coverage.sh
@@ -12,13 +12,12 @@ function run() {
 
 function compile() {
     ./autogen.sh
-    ./configure --enable-coverage
+    ./configure --enable-coverage --enable-debugging-checks
     make -j 10
 }
 
 function generate() {
-    gcov tests/*c --object-directory src/.libs/
-    lcov -c --directory src/ --output-file test
+    lcov -c --directory src/ --directory tests/ --output-file test
     genhtml test --output-directory coverage
 }
 
diff --git a/src/csds_cosmology.h b/src/csds_cosmology.h
index 62027d9..b0372d4 100644
--- a/src/csds_cosmology.h
+++ b/src/csds_cosmology.h
@@ -71,7 +71,7 @@ struct csds_cosmology {
  *
  * @return The equation of state.
  */
-__attribute__((always_inline)) INLINE static double csds_cosmology_get_w_tilde(
+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);
@@ -86,7 +86,7 @@ __attribute__((always_inline)) INLINE static double csds_cosmology_get_w_tilde(
  *
  * @return The first derivative of the scale factor.
  */
-__attribute__((always_inline)) INLINE static double csds_cosmology_get_E(
+INLINE static double csds_cosmology_get_E(
     const struct csds_cosmology *cosmo, const double scale) {
 
   const double a_inv = 1. / scale;
@@ -108,7 +108,7 @@ __attribute__((always_inline)) INLINE static double csds_cosmology_get_E(
  *
  * @return The first derivative of the scale factor.
  */
-__attribute__((always_inline)) INLINE static double
+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);
@@ -123,7 +123,7 @@ csds_cosmology_get_scale_factor_derivative(const struct csds_cosmology *cosmo,
  *
  * @return The first derivative of the scale factor.
  */
-__attribute__((always_inline)) INLINE static double
+INLINE static double
 csds_cosmology_get_scale_factor_second_derivative(
     const struct csds_cosmology *cosmo, const double scale) {
 
@@ -157,7 +157,7 @@ csds_cosmology_get_scale_factor_second_derivative(
  * @param vel The velocity to modify.
  * @param acc The acceleration to modify.
  */
-__attribute__((always_inline)) INLINE static void
+INLINE static void
 csds_cosmology_add_factor_position(const struct csds_cosmology *cosmo,
                                    const double scale, float *vel, float *acc) {
 
@@ -183,7 +183,7 @@ csds_cosmology_add_factor_position(const struct csds_cosmology *cosmo,
  * @param scale The scale factor.
  * @param acc The acceleration to modify.
  */
-__attribute__((always_inline)) INLINE static void
+INLINE static void
 csds_cosmology_add_factor_velocity(const struct csds_cosmology *cosmo,
                                    const double scale, float *acc) {
 
diff --git a/src/csds_error.h b/src/csds_error.h
index 6c85460..7b8bb47 100644
--- a/src/csds_error.h
+++ b/src/csds_error.h
@@ -74,7 +74,7 @@
 /**
  * @brief Print the python trace back
  */
-__attribute__((always_inline)) INLINE static void csds_loader_print_traceback(
+INLINE static void csds_loader_print_traceback(
     void) {
 
   PyGILState_STATE gstate = PyGILState_Ensure();
diff --git a/src/csds_fields.h b/src/csds_fields.h
index 923b793..1a04ac3 100644
--- a/src/csds_fields.h
+++ b/src/csds_fields.h
@@ -95,7 +95,7 @@ struct field_information {
 /**
  * @brief Get the #field_enum from the field name.
  */
-__attribute__((always_inline)) INLINE static enum field_enum field_enum_get(
+INLINE static enum field_enum field_enum_get(
     const char *name) {
   for (enum field_enum i = field_enum_special_flags; i < field_enum_count;
        i++) {
@@ -108,7 +108,7 @@ __attribute__((always_inline)) INLINE static enum field_enum field_enum_get(
 /**
  * @brief Get the name of the field.
  */
-__attribute__((always_inline)) INLINE static const char *field_enum_get_name(
+INLINE static const char *field_enum_get_name(
     enum field_enum field) {
   return fields_names[field];
 }
@@ -116,7 +116,7 @@ __attribute__((always_inline)) INLINE static const char *field_enum_get_name(
 /**
  * @brief Get the #field_enum for the first derivative.
  */
-__attribute__((always_inline)) INLINE static enum field_enum
+INLINE static enum field_enum
 field_enum_get_first_derivative(enum field_enum field) {
   switch (field) {
     case field_enum_coordinates:
@@ -132,7 +132,7 @@ field_enum_get_first_derivative(enum field_enum field) {
 /**
  * @brief Get the #field_enum for the second derivative.
  */
-__attribute__((always_inline)) INLINE static enum field_enum
+INLINE static enum field_enum
 field_enum_get_second_derivative(enum field_enum field) {
   switch (field) {
     case field_enum_coordinates:
@@ -146,7 +146,7 @@ field_enum_get_second_derivative(enum field_enum field) {
 /**
  * @brief Get the size of the field.
  */
-__attribute__((always_inline)) INLINE static unsigned int field_enum_get_size(
+INLINE static unsigned int field_enum_get_size(
     enum field_enum field) {
   switch (field) {
 
@@ -192,7 +192,7 @@ __attribute__((always_inline)) INLINE static unsigned int field_enum_get_size(
 /**
  * @brief Set a #csds_python_field for the field.
  */
-__attribute__((always_inline)) INLINE static struct csds_python_field
+INLINE static struct csds_python_field
 field_enum_get_python_field(enum field_enum field) {
   switch (field) {
     case field_enum_coordinates:
@@ -311,7 +311,7 @@ struct mask_data {
  * @param all_data All the #mask_data from the header.
  * @param size_all_data Size of all_data.
  */
-__attribute__((always_inline)) INLINE static void
+INLINE static void
 field_information_set_derivative(struct single_field *der,
                                  enum field_enum field,
                                  struct mask_data *all_data,
@@ -355,7 +355,7 @@ field_information_set_derivative(struct single_field *der,
  * @param all_data All the #mask_data from the header.
  * @param size_all_data Size of all_data.
  */
-__attribute__((always_inline)) INLINE static void field_information_init(
+INLINE static void field_information_init(
     struct field_information *info, struct mask_data *data,
     struct mask_data *all_data, int size_all_data) {
 
@@ -388,7 +388,7 @@ __attribute__((always_inline)) INLINE static void field_information_init(
  * @param info An array of #field_information.
  * @param size The size of the array.
  */
-__attribute__((always_inline)) INLINE static void
+INLINE static void
 field_information_set_positions(struct field_information *info,
                                 const unsigned int size) {
   for (unsigned int i = 0; i < size; i++) {
diff --git a/src/csds_hashmap.h b/src/csds_hashmap.h
index b44be08..3288c76 100644
--- a/src/csds_hashmap.h
+++ b/src/csds_hashmap.h
@@ -88,12 +88,12 @@ size_t csds_hashmap_get_number_buckets(struct csds_hashmap *map);
 struct index_data *csds_hashmap_get_from_position(struct csds_hashmap *map,
                                                   size_t i);
 
-__attribute__((always_inline)) INLINE static struct bucket *bucket_at(
+INLINE static struct bucket *bucket_at(
     struct csds_hashmap *map, size_t index) {
   return (struct bucket *)(((char *)map->buckets) + (map->bucketsz * index));
 }
 
-__attribute__((always_inline)) INLINE static struct index_data *bucket_item(
+INLINE static struct index_data *bucket_item(
     struct bucket *entry) {
   char *out = ((char *)entry) + sizeof(struct bucket);
   return (struct index_data *)out;
@@ -214,7 +214,7 @@ static uint64_t SIP64(const id_type *key) {
 }
 
 // hashmap_sip returns a hash value for `data` using SipHash-2-4.
-__attribute__((always_inline)) INLINE static uint64_t get_hash(id_type x) {
+INLINE static uint64_t get_hash(id_type x) {
   return SIP64(&x) & ((1LU << 48) - 1);
 }
 
diff --git a/src/csds_header.h b/src/csds_header.h
index d301c43..4f2aff6 100644
--- a/src/csds_header.h
+++ b/src/csds_header.h
@@ -81,7 +81,7 @@ const struct field_information *header_get_field_from_name(
  * @brief Check if the offset are forward.
  * @param h The #header.
  */
-__attribute__((always_inline)) INLINE static int header_is_forward(
+INLINE static int header_is_forward(
     const struct header *h) {
   return h->offset_direction == csds_offset_forward;
 }
@@ -90,7 +90,7 @@ __attribute__((always_inline)) INLINE static int header_is_forward(
  * @brief Check if the offset are backward.
  * @param h The #header.
  */
-__attribute__((always_inline)) INLINE static int header_is_backward(
+INLINE static int header_is_backward(
     const struct header *h) {
   return h->offset_direction == csds_offset_backward;
 }
@@ -99,7 +99,7 @@ __attribute__((always_inline)) INLINE static int header_is_backward(
  * @brief Check if the offset are corrupted.
  * @param h The #header.
  */
-__attribute__((always_inline)) INLINE static int header_is_corrupted(
+INLINE static int header_is_corrupted(
     const struct header *h) {
   return h->offset_direction == csds_offset_corrupted;
 }
diff --git a/src/csds_inline.h b/src/csds_inline.h
index e2a41ff..30b1662 100644
--- a/src/csds_inline.h
+++ b/src/csds_inline.h
@@ -26,16 +26,20 @@
  * @brief Defines inline
  */
 #ifndef INLINE
+#ifdef CSDS_NO_INLINE
+#define INLINE
+#else
 #ifdef __cplusplus
-#define INLINE inline
+#define INLINE __attribute__((always_inline)) inline
 #else
 #if __GNUC__ && !__GNUC_STDC_INLINE__
-#define INLINE extern inline
+#define INLINE __attribute__((always_inline)) extern inline
 #else
-#define INLINE inline
-#endif
-#endif
-#endif
+#define INLINE __attribute__((always_inline)) inline
+#endif // GNU
+#endif // C++
+#endif // CSDS_NO_INLINE
+#endif // INLINE
 
 /**
  * @brief Defines unused
diff --git a/src/csds_interpolation.h b/src/csds_interpolation.h
index d522731..a870f96 100644
--- a/src/csds_interpolation.h
+++ b/src/csds_interpolation.h
@@ -61,7 +61,7 @@
  *
  * @return The function evaluated at t.
  */
-__attribute__((always_inline)) INLINE static double
+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,
@@ -122,7 +122,7 @@ interpolate_quintic_hermite_spline(const double t0, const double x0,
  *
  * @return The function evaluated at t.
  */
-__attribute__((always_inline)) INLINE static float
+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,
@@ -169,7 +169,7 @@ interpolate_cubic_hermite_spline(const double t0, const float v0,
  * @param periodic Should the periodic boundary be applied?
  * @param params The simulation's #csds_parameters.
  */
-__attribute__((always_inline)) INLINE static void
+INLINE static void
 interpolate_quintic_double_float_ND(
     const double t_before, const struct csds_reader_field *restrict before,
     const double t_after, const struct csds_reader_field *restrict after,
@@ -266,7 +266,7 @@ interpolate_quintic_double_float_ND(
  * @param periodic Should the periodic boundary be applied?
  * @param params The simulation's #csds_parameters.
  */
-__attribute__((always_inline)) INLINE static void interpolate_cubic_float_ND(
+INLINE static void interpolate_cubic_float_ND(
     const double t_before, const struct csds_reader_field *restrict before,
     const double t_after, const struct csds_reader_field *restrict after,
     void *restrict output, const double t, const int dimension,
@@ -325,7 +325,7 @@ __attribute__((always_inline)) INLINE static void interpolate_cubic_float_ND(
  * @param periodic Should the periodic boundary be applied?
  * @param params The simulation's #csds_parameters.
  */
-__attribute__((always_inline)) INLINE static void interpolate_linear_float_ND(
+INLINE static void interpolate_linear_float_ND(
     const double t_before, const struct csds_reader_field *restrict before,
     const double t_after, const struct csds_reader_field *restrict after,
     void *restrict output, const double t, const int dimension,
@@ -366,7 +366,7 @@ __attribute__((always_inline)) INLINE static void interpolate_linear_float_ND(
  * @param periodic Should the periodic boundary be applied?
  * @param params The simulation's #csds_parameters.
  */
-__attribute__((always_inline)) INLINE static void interpolate_linear_double_ND(
+INLINE static void interpolate_linear_double_ND(
     const double t_before, const struct csds_reader_field *restrict before,
     const double t_after, const struct csds_reader_field *restrict after,
     void *restrict output, const double t, const int dimension,
@@ -406,7 +406,7 @@ __attribute__((always_inline)) INLINE static void interpolate_linear_double_ND(
  * @param periodic Should the periodic boundary be applied?
  * @param params The simulation's #csds_parameters.
  */
-__attribute__((always_inline)) INLINE static void interpolate_linear_float(
+INLINE static void interpolate_linear_float(
     const double t_before, const struct csds_reader_field *restrict before,
     const double t_after, const struct csds_reader_field *restrict after,
     void *restrict output, const double t, ATTR_UNUSED int periodic,
diff --git a/src/csds_loader_io.h b/src/csds_loader_io.h
index e546bb5..c352a2d 100644
--- a/src/csds_loader_io.h
+++ b/src/csds_loader_io.h
@@ -58,7 +58,7 @@ void csds_loader_io_munmap_file(struct mapped_file *map);
  *
  * @return memory after the record header.
  */
-__attribute__((always_inline)) INLINE static char *csds_loader_io_read_mask(
+INLINE static char *csds_loader_io_read_mask(
     char *data, mask_type *mask, size_t *diff_offset) {
   /* read mask */
   if (mask) {
@@ -86,7 +86,7 @@ __attribute__((always_inline)) INLINE static char *csds_loader_io_read_mask(
 
  * @return memory after the data written.
  */
-__attribute__((always_inline)) INLINE static char *csds_loader_io_read_data(
+INLINE static char *csds_loader_io_read_data(
     char *data, const size_t size, void *p) {
   memcpy(p, data, size);
   return data + size;
@@ -101,7 +101,7 @@ __attribute__((always_inline)) INLINE static char *csds_loader_io_read_data(
  *
  * @return memory after the data written.
  */
-__attribute__((always_inline)) INLINE static char *csds_loader_io_write_data(
+INLINE static char *csds_loader_io_write_data(
     char *data, const size_t size, const void *p) {
   memcpy(data, p, size);
 
diff --git a/src/csds_logfile_writer.h b/src/csds_logfile_writer.h
index 4c5d6fe..1df1501 100644
--- a/src/csds_logfile_writer.h
+++ b/src/csds_logfile_writer.h
@@ -82,7 +82,7 @@ char *csds_logfile_writer_write_begining_header(
  *
  * @return updated buff
  */
-__attribute__((always_inline)) INLINE static char *csds_write_record_header(
+INLINE static char *csds_write_record_header(
     char *buff, const unsigned int *mask, const size_t *offset,
     const size_t offset_new) {
   /* write mask. */
@@ -106,7 +106,7 @@ __attribute__((always_inline)) INLINE static char *csds_write_record_header(
  * csds_logfile_writer file.
  * @return A pointer to the memory-mapped chunk of data.
  */
-__attribute__((always_inline)) INLINE static void *csds_logfile_writer_get(
+INLINE static void *csds_logfile_writer_get(
     struct csds_logfile_writer *d, size_t count, size_t *offset) {
   size_t local_offset = 0;
 
@@ -158,7 +158,7 @@ INLINE static uint32_t csds_pack_flags_and_data(enum csds_special_flags flag,
  * @param size number of bytes to write
  * @param p pointer to the data
  */
-__attribute__((always_inline)) INLINE static void csds_write_data(
+INLINE static void csds_write_data(
     struct csds_logfile_writer *d, size_t *offset, size_t size, const void *p) {
   /* get buffer. */
   char *buff = csds_logfile_writer_get(d, size, offset);
diff --git a/src/csds_openmp.h b/src/csds_openmp.h
index 119c4c5..972da5b 100644
--- a/src/csds_openmp.h
+++ b/src/csds_openmp.h
@@ -32,7 +32,7 @@
 /**
  * @brief Wrapper for omp_get_thread_num
  */
-__attribute__((always_inline)) INLINE static int csds_get_thread_num(void) {
+INLINE static int csds_get_thread_num(void) {
 #ifdef CSDS_WITH_OPENMP
   return omp_get_thread_num();
 #else
diff --git a/src/csds_particle.h b/src/csds_particle.h
index 7d69b91..15c2f0a 100644
--- a/src/csds_particle.h
+++ b/src/csds_particle.h
@@ -56,7 +56,7 @@ enum csds_reader_type {
  *
  * @return position after the record.
  */
-__attribute__((always_inline)) INLINE static size_t csds_particle_read_field(
+INLINE static size_t csds_particle_read_field(
     size_t offset, void *output, const struct field_information *field_read,
     const int derivative, mask_type *mask, size_t *h_offset,
     const struct field_information *fields, int number_fields, char *log_map) {
@@ -156,7 +156,7 @@ INLINE static enum csds_special_flags csds_unpack_flags_and_data(
  *
  * @return The special flag.
  */
-__attribute__((always_inline)) INLINE static enum csds_special_flags
+INLINE static enum csds_special_flags
 csds_particle_read_special_flag(size_t offset, mask_type *mask,
                                 size_t *h_offset, int *data, int *part_type,
                                 char *map) {
@@ -193,7 +193,7 @@ csds_particle_read_special_flag(size_t offset, mask_type *mask,
  * @param field The field to reconstruct.
  * @param params The simulation's #csds_parameters.
  */
-__attribute__((always_inline)) INLINE static void
+INLINE static void
 csds_particle_interpolate_field(const double time_before,
                                 const struct csds_reader_field *restrict before,
                                 const double time_after,
diff --git a/src/csds_python_tools.h b/src/csds_python_tools.h
index 868ca90..2e18fa7 100644
--- a/src/csds_python_tools.h
+++ b/src/csds_python_tools.h
@@ -72,7 +72,7 @@ struct csds_python_field {
  *
  * @return The initialized structure.
  */
-__attribute__((always_inline)) INLINE static struct csds_python_field
+INLINE static struct csds_python_field
 csds_loader_python_field(int dimension, int numpy_type) {
   struct csds_python_field ret;
 
@@ -104,7 +104,7 @@ csds_loader_python_field(int dimension, int numpy_type) {
  * ...)
  * @param name The name of the subfield.
  */
-__attribute__((always_inline)) INLINE static void
+INLINE static void
 csds_loader_python_field_add_subfield(struct csds_python_field *field,
                                       const char *name, size_t offset,
                                       int dimension, int numpy_type) {
@@ -134,7 +134,7 @@ csds_loader_python_field_add_subfield(struct csds_python_field *field,
  *
  * @param field The #csds_python_field to clean.
  */
-__attribute__((always_inline)) INLINE static void csds_loader_python_field_free(
+INLINE static void csds_loader_python_field_free(
     struct csds_python_field *field) {
   if (field->typenum == CUSTOM_NPY_TYPE) {
     free(field->subfields);
@@ -149,7 +149,7 @@ __attribute__((always_inline)) INLINE static void csds_loader_python_field_free(
  * @param typenum The C typenum (e.g. NPY_FLOAT32, NPY_DOUBLE, ...).
  * @param dimension The number of dimensions.
  */
-__attribute__((always_inline)) INLINE static PyObject *
+INLINE static PyObject *
 csds_python_tools_get_format_string(int typenum, int dimension) {
 
   const char *type_format = NULL;
diff --git a/src/csds_python_wrapper.c b/src/csds_python_wrapper.c
index 81d8123..b7ac857 100644
--- a/src/csds_python_wrapper.c
+++ b/src/csds_python_wrapper.c
@@ -82,7 +82,7 @@ typedef struct {
 /**
  * @brief Deallocate the memory of a particle reader.
  */
-__attribute__((always_inline)) INLINE static void particle_reader_dealloc(
+INLINE static void particle_reader_dealloc(
     PyParticleReader *self) {
   Py_TYPE(self)->tp_free((PyObject *)self);
 }
@@ -101,7 +101,7 @@ static void Reader_dealloc(PyObjectReader *self) {
 /**
  * @brief Allocate the memory of a particle reader.
  */
-__attribute__((always_inline)) INLINE static PyObject *particle_reader_new(
+INLINE static PyObject *particle_reader_new(
     PyTypeObject *type, ATTR_UNUSED PyObject *args,
     ATTR_UNUSED PyObject *kwds) {
   PyParticleReader *self = (PyParticleReader *)type->tp_alloc(type, 0);
@@ -229,7 +229,7 @@ static PyObject *get_box_size(PyObject *self, ATTR_UNUSED PyObject *args) {
  *
  * @return The python list of numpy array.
  */
-__attribute__((always_inline)) INLINE static PyObject *
+INLINE static PyObject *
 csds_loader_create_output(void **output, const enum field_enum *fields,
                           const int n_fields, uint64_t n_tot) {
 
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 215adbd..cbc09da 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -28,11 +28,11 @@ AM_CFLAGS += $(PYTHON_EXTRA_COMPILER_FLAG)
 
 # List of programs and scripts to run in the test suite
 TESTS = testLogfileHeader testLogfileReader testTimeArray testQuickSort testVirtualReality
-TESTS += testHashmap
+TESTS += testHashmap testInterpolation
 
 # List of test programs to compile
 check_PROGRAMS = testLogfileHeader testLogfileReader testTimeArray testQuickSort testVirtualReality
-check_PROGRAMS += testHashmap
+check_PROGRAMS += testHashmap testInterpolation
 
 # Rebuild tests when CSDS is updated.
 $(check_PROGRAMS): ../src/.libs/libcsds.a ../src/.libs/libcsds_writer.a
@@ -44,6 +44,7 @@ testTimeArray_SOURCES = testTimeArray.c
 testQuickSort_SOURCES = testQuickSort.c
 testVirtualReality_SOURCES = testVirtualReality.c
 testHashmap_SOURCES = testHashmap.c
+testInterpolation_SOURCES = testInterpolation.c
 
 # Files necessary for distribution
 EXTRA_DIST = 
diff --git a/tests/testInterpolation.c b/tests/testInterpolation.c
new file mode 100644
index 0000000..ba94e60
--- /dev/null
+++ b/tests/testInterpolation.c
@@ -0,0 +1,245 @@
+/*******************************************************************************
+ * This file is part of CSDS.
+ * Copyright (C) 2019 Loic Hausammann (loic.hausammann@epfl.ch)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+
+/* Include local headers */
+#include "csds_header.h"
+#include "csds_logfile.h"
+#include "csds_parser.h"
+#include "csds_reader.h"
+
+/* Include standard headers */
+#include <assert.h>
+
+
+#define err_abs(a, b) fabs(a - b)
+
+void pos_func(double t, double *p) {
+  const double t2 = t * t;
+  const double f = exp(t2 - 5 * t - 2);
+  p[0] = f;
+  p[1] = f;
+  p[2] = f;
+}
+
+void vel_func(double t, float *v) {
+  const double t2 = t * t;
+  const float f = (2 * t - 5) * exp(t2 - 5 * t - 2);
+  v[0] = f;
+  v[1] = f;
+  v[2] = f;
+}
+
+void acc_func(double t, float *a) {
+  const double t2 = t * t;
+  const float f = (4 * t2 - 20 * t + 27) *
+    exp(t2 - 5 * t - 2);
+  a[0] = f;
+  a[1] = f;
+  a[2] = f;
+}
+
+long long id_func(void) {
+  return 1698241543;
+}
+
+int main(void) {
+  /* Fake initialization of the reader */
+  struct csds_reader reader;
+  strcpy(reader.basename, "test_interpolation_0000");
+  reader.verbose = 2;
+
+  /* Initialize the parameters */
+  struct csds_parameters params;
+  csds_parameters_init(&params, &reader);
+
+  /* Construct the fields */
+  const int n_fields = 4;
+  struct mask_data data[n_fields];
+
+  /* Set the coordinate */
+  data[0] = (struct mask_data) {
+      .size = field_enum_get_size(field_enum_coordinates),
+      .mask = 1 << 2,
+      .name = "Coordinates"
+  };
+
+  /* Set the velocity */
+  data[1] = (struct mask_data) {
+      .size = field_enum_get_size(field_enum_velocities),
+      .mask = 1 << 3,
+      .name = "Velocities"
+  };
+
+  /* Set the acceleration */
+  data[2] = (struct mask_data) {
+      .size = field_enum_get_size(field_enum_accelerations),
+      .mask = 1 << 4,
+      .name = "Accelerations"
+  };
+
+  /* Set the IDs */
+  data[3] = (struct mask_data) {
+      .size = field_enum_get_size(field_enum_particles_ids),
+      .mask = 1 << 5,
+      .name = "ParticleIDs",
+  };
+
+  /* Construct the field information */
+  struct field_information info[n_fields];
+  for(int i = 0; i < n_fields; i++) {
+    field_information_init(&info[i], &data[i], data, n_fields);
+  }
+  field_information_set_positions(info, n_fields);
+
+  /* Now test */
+  const int n_evaluations = 10;
+  const double t_beg = 0.2;
+  const double t_end = 0.8;
+  const double t_eval = 0.5 * (t_beg + t_end);
+  double last_error[n_fields];
+  for(int i = 0; i < n_fields; i++) {
+    last_error[i] = 0;
+  }
+
+  /* Check the convergence */
+  for(int i = 1; i <= n_evaluations; i++) {
+    /* Get the time */
+    const double dt = 0.5 * (t_end - t_beg) / i;
+    const double t1 = t_eval - dt;
+    const double t2 = t_eval + dt;
+    message("Interpolating at t=%g from t=[%g, %g]",
+            t_eval, t1, t2);
+
+    /* Set the state at t1 and t2 */
+    double p1[3];
+    pos_func(t1, p1);
+    double p2[3];
+    pos_func(t2, p2);
+
+    float v1[3];
+    vel_func(t1, v1);
+    float v2[3];
+    vel_func(t2, v2);
+
+    float a1[3];
+    acc_func(t1, a1);
+    float a2[3];
+    acc_func(t2, a2);
+
+    long long id1 = id_func();
+    long long id2 = id_func();
+
+    /* Create the output variables */
+    double out_p[3];
+    float out_v[3];
+    float out_a[3];
+    long long out_id = 0;
+
+    /* Do the interpolations */
+    /* Coordinates */
+    struct csds_reader_field before = {
+        .field = p1, .first_deriv = v1,
+        .second_deriv = a1,
+   };
+   struct csds_reader_field after = {
+       .field = p2, .first_deriv = v2,
+       .second_deriv = a2,
+   };
+
+   csds_particle_interpolate_field(
+       t1, &before, t2, &after, &out_p,
+       t_eval, &info[0], &params);
+
+   /* Velocities */
+   before = (struct csds_reader_field) {
+       .field = v1, .first_deriv = a1,
+       .second_deriv = NULL,
+   };
+   after = (struct csds_reader_field) {
+       .field = v2, .first_deriv = a2,
+       .second_deriv = NULL,
+   };
+
+   csds_particle_interpolate_field(
+       t1, &before, t2, &after, &out_v,
+       t_eval, &info[1], &params);
+
+   /* Accelerations */
+   before = (struct csds_reader_field) {
+       .field = a1, .first_deriv = NULL,
+       .second_deriv = NULL,
+   };
+   after = (struct csds_reader_field) {
+       .field = a2, .first_deriv = NULL,
+       .second_deriv = NULL,
+   };
+
+   csds_particle_interpolate_field(
+       t1, &before, t2, &after, &out_a,
+       t_eval, &info[2], &params);
+
+   /* IDs */
+   before = (struct csds_reader_field) {
+       .field = &id1, .first_deriv = NULL,
+       .second_deriv = NULL,
+   };
+   after = (struct csds_reader_field) {
+       .field = &id2, .first_deriv = NULL,
+       .second_deriv = NULL,
+   };
+
+   csds_particle_interpolate_field(
+       t1, &before, t2, &after, &out_id,
+       t_eval, &info[3], &params);
+
+   /* Now check the results */
+   double current_error[n_fields];
+   double p_ref[3];
+   pos_func(t_eval, p_ref);
+   current_error[0] = err_abs(out_p[0], p_ref[0]);
+
+   float v_ref[3];
+   vel_func(t_eval, v_ref);
+   current_error[1] = err_abs(out_v[0], v_ref[0]);
+
+   float a_ref[3];
+   acc_func(t_eval, a_ref);
+   current_error[2] = err_abs(out_a[0], a_ref[0]);
+
+   assert(out_id == id_func());
+
+   /* For the first element, we cannot check the convergence */
+   if (i != 1) {
+     const double fact = (i - 1.) / (double) i;
+     /* Quintic convergence for coordinates */
+     assert(current_error[0] <= fmax(last_error[0] * pow(fact, 5), 1e-10));
+     /* Cubic convergence for velocities */
+     assert(current_error[1] <= fmaxf(last_error[1] * pow(fact, 3), 1e-6));
+     /* Linear convergence for accelerations */
+     assert(current_error[2] <= fmaxf(last_error[2] * fact, 1e-6));
+   }
+
+   /* Update the last error */
+   for(int field = 0; field < n_fields; field++) {
+     last_error[field] = current_error[field];
+   }
+ }
+
+ return 0;
+}
diff --git a/tests/test_interpolation.yml b/tests/test_interpolation.yml
new file mode 100644
index 0000000..e75ad3b
--- /dev/null
+++ b/tests/test_interpolation.yml
@@ -0,0 +1,66 @@
+Header:
+  BoxSize: [100000, 100000, 100000]
+  Dimension: 3
+  RunName: Untitled SWIFT simulation
+  Periodic: 1
+  NumberParts: 0
+  NumberSParts: 0
+  NumberGParts: 65536
+
+Cosmology:
+  Omega_cdm: 0.2305
+  Omega_lambda: 0.724
+  Omega_b: 0.0455
+  Omega_r: 0
+  Omega_k: 0
+  Omega_nu_0: 0
+  w_0: -1
+  w_a: 0
+  Hubble0: 70.3
+
+InternalUnitSystem:
+  UnitMass_in_cgs: 1.98841e+43
+  UnitLength_in_cgs: 3.08568e+24
+  UnitTime_in_cgs: 3.08568e+19
+  UnitCurrent_in_cgs: 1
+  UnitTemp_in_cgs: 1
+
+Code:
+  Code: SWIFT
+  CodeVersion: 0.9.0
+  CompilerName: GCC
+  CompilerVersion: 9.3.0
+  GitBranch: master
+  GitRevision: v0.9.0-565-g0122ba6f-dirty
+  GitDate: 2021-07-14 09:33:30 +0000
+  ConfigurationOptions: '--enable-csds'
+  RandomSeed: 0
+
+Policy:
+  steal: 1
+  keep: 1
+  block: 0
+  cpu tight: 0
+  mpi: 1
+  numa affinity: 1
+  hydro: 0
+  self gravity: 1
+  external gravity: 0
+  cosmological integration: 0
+  drift everything: 0
+  reconstruct multi-poles: 0
+  temperature: 0
+  cooling: 0
+  stars: 0
+  structure finding: 0
+  star formation: 0
+  feedback: 0
+  black holes: 0
+  fof search: 0
+  time-step limiter: 0
+  time-step sync: 0
+  csds: 1
+  line of sight: 0
+  sink: 0
+  rt: 0
+
-- 
GitLab