Merge branch 'DOPAIR2_fix' into 'master'

Dopair2 fix

Large set of improvements to the hydro interactions functions:
 - Correct version of DOPAIR2() (Fixes #361). The symmetric condition was not correctly taken into account when computing the distance on the axis.
 - Removed the old vectorization scheme from the interaction functions. This is superseded by @jwillis' work.
 - Added the check for active/inactive in the *_NAIVE() functions. These can now correctly be used throughout a run.
 - Made all the interactions take place in the frame of the cell `cj` for PAIRs and `c` for the SELFs.
 - Added a new particle type. The particles do nothing apart from recording who they interact with in the density and force loops. This is solely designed for debugging. 

See merge request !419

configure.ac

@@ -627,7 +627,7 @@ fi
 # Hydro scheme.
 AC_ARG_WITH([hydro],
    [AS_HELP_STRING([--with-hydro=<scheme>],
-      [Hydro dynamics to use @<:@gadget2, minimal, hopkins, default, gizmo, shadowfax default: gadget2@:>@]
+      [Hydro dynamics to use @<:@gadget2, minimal, hopkins, default, gizmo, shadowfax debug default: gadget2@:>@]
    )],
    [with_hydro="$withval"],
    [with_hydro="gadget2"]
@@ -639,6 +639,9 @@ case "$with_hydro" in
    minimal)
       AC_DEFINE([MINIMAL_SPH], [1], [Minimal SPH])
    ;; 
+   debug)
+      AC_DEFINE([DEBUG_INTERACTIONS_SPH], [1], [Debug SPH])
+   ;; 
    hopkins)
       AC_DEFINE([HOPKINS_PE_SPH], [1], [Pressure-Entropy SPH])
    ;; 

examples/SedovBlast_3D/sedov.yml

@@ -32,4 +32,5 @@ SPH:
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  ./sedov.hdf5          # The file to read
+  h_scaling:           3.33
 

src/debug.c

@@ -21,16 +21,15 @@
  ******************************************************************************/
 
 /* Config parameters. */
-#include "../config.h"
+
+/* This object's header. */
+#include "debug.h"
 
 /* Some standard headers. */
 #include <float.h>
 #include <stdio.h>
 #include <unistd.h>
 
-/* This object's header. */
-#include "debug.h"
-
 /* Local includes. */
 #include "active.h"
 #include "cell.h"
@@ -41,7 +40,9 @@
 #include "space.h"
 
 /* Import the right hydro definition */
-#if defined(MINIMAL_SPH)
+#if defined(DEBUG_INTERACTIONS_SPH)
+#include "./hydro/DebugInteractions/hydro_debug.h"
+#elif defined(MINIMAL_SPH)
 #include "./hydro/Minimal/hydro_debug.h"
 #elif defined(GADGET2_SPH)
 #include "./hydro/Gadget2/hydro_debug.h"

src/debug.h

@@ -19,6 +19,9 @@
 #ifndef SWIFT_DEBUG_H
 #define SWIFT_DEBUG_H
 
+/* Config parameters. */
+#include "../config.h"
+
 /* Includes. */
 #include "cell.h"
 #include "part.h"

src/hydro.h

@@ -27,7 +27,11 @@
 #include "kernel_hydro.h"
 
 /* Import the right functions */
-#if defined(MINIMAL_SPH)
+#if defined(DEBUG_INTERACTIONS_SPH)
+#include "./hydro/DebugInteractions/hydro.h"
+#include "./hydro/DebugInteractions/hydro_iact.h"
+#define SPH_IMPLEMENTATION "Debug SELF/PAIR"
+#elif defined(MINIMAL_SPH)
 #include "./hydro/Minimal/hydro.h"
 #include "./hydro/Minimal/hydro_iact.h"
 #define SPH_IMPLEMENTATION "Minimal version of SPH (e.g. Price 2010)"

src/hydro/DebugInteractions/hydro.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * 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_DEBUG_INTERACTIONS_HYDRO_H
+#define SWIFT_DEBUG_INTERACTIONS_HYDRO_H
+
+/**
+ * @file DebugInteractions/hydro.h
+ * @brief Empty SPH implementation used solely to test the SELF/PAIR routines.
+ */
+
+#include "hydro_properties.h"
+#include "hydro_space.h"
+#include "part.h"
+
+#include <float.h>
+
+/**
+ * @brief Returns the internal energy of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
+    const struct part *restrict p) {
+
+  return 0.f;
+}
+
+/**
+ * @brief Returns the pressure of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_pressure(
+    const struct part *restrict p) {
+
+  return 0.f;
+}
+
+/**
+ * @brief Returns the entropy of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_entropy(
+    const struct part *restrict p) {
+
+  return 0.f;
+}
+
+/**
+ * @brief Returns the sound speed of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
+    const struct part *restrict p) {
+
+  return 0.f;
+}
+
+/**
+ * @brief Returns the density of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_density(
+    const struct part *restrict p) {
+
+  return 0.f;
+}
+
+/**
+ * @brief Returns the mass of a particle
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_mass(
+    const struct part *restrict p) {
+  return 0.f;
+}
+
+/**
+ * @brief Returns the time derivative of internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy_dt(
+    const struct part *restrict p) {
+  return 0.f;
+}
+
+/**
+ * @brief Returns the time derivative of internal energy of a particle
+ *
+ * We assume a constant density.
+ *
+ * @param p The particle of interest.
+ * @param du_dt The new time derivative of the internal energy.
+ */
+__attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
+    struct part *restrict p, float du_dt) {}
+
+/**
+ * @brief Computes the hydro time-step of a given particle
+ *
+ * @param p Pointer to the particle data
+ * @param xp Pointer to the extended particle data
+ *
+ */
+__attribute__((always_inline)) INLINE static float hydro_compute_timestep(
+    const struct part *restrict p, const struct xpart *restrict xp,
+    const struct hydro_props *restrict hydro_properties) {
+
+  return FLT_MAX;
+}
+
+/**
+ * @brief Does some extra hydro operations once the actual physical time step
+ * for the particle is known.
+ *
+ * @param p The particle to act upon.
+ * @param dt Physical time step of the particle during the next step.
+ */
+__attribute__((always_inline)) INLINE static void hydro_timestep_extra(
+    struct part *p, float dt) {}
+
+/**
+ * @brief Prepares a particle for the density calculation.
+ *
+ * Zeroes all the relevant arrays in preparation for the sums taking place in
+ * the variaous density tasks
+ *
+ * @param p The particle to act upon
+ * @param hs #hydro_space containing hydro specific space information.
+ */
+__attribute__((always_inline)) INLINE static void hydro_init_part(
+    struct part *restrict p, const struct hydro_space *hs) {
+
+  for (int i = 0; i < 256; ++i) p->ids_ngbs_density[i] = -1;
+  p->num_ngb_density = 0;
+
+  p->density.wcount = 0.f;
+  p->density.wcount_dh = 0.f;
+}
+
+/**
+ * @brief Finishes the density calculation.
+ *
+ * Multiplies the density and number of neighbours by the appropiate constants
+ * and add the self-contribution term.
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_density(
+    struct part *restrict p) {
+
+  /* Some smoothing length multiples. */
+  const float h = p->h;
+  const float h_inv = 1.0f / h;                       /* 1/h */
+  const float h_inv_dim = pow_dimension(h_inv);       /* 1/h^d */
+  const float h_inv_dim_plus_one = h_inv_dim * h_inv; /* 1/h^(d+1) */
+
+  /* Final operation on the density (add self-contribution). */
+  p->density.wcount += kernel_root;
+  p->density.wcount_dh -= hydro_dimension * kernel_root;
+
+  /* Finish the calculation by inserting the missing h-factors */
+  p->density.wcount *= h_inv_dim;
+  p->density.wcount_dh *= h_inv_dim_plus_one;
+}
+
+/**
+ * @brief Sets all particle fields to sensible values when the #part has 0 ngbs.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_part_has_no_neighbours(
+    struct part *restrict p, struct xpart *restrict xp) {
+
+  /* Some smoothing length multiples. */
+  const float h = p->h;
+  const float h_inv = 1.0f / h;                 /* 1/h */
+  const float h_inv_dim = pow_dimension(h_inv); /* 1/h^d */
+
+  /* Re-set problematic values */
+  p->density.wcount = kernel_root * kernel_norm * h_inv_dim;
+  p->density.wcount_dh = 0.f;
+}
+
+/**
+ * @brief Prepare a particle for the force calculation.
+ *
+ * Computes viscosity term, conduction term and smoothing length gradient terms.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ * @param ti_current The current time (on the timeline)
+ * @param timeBase The minimal time-step size
+ */
+__attribute__((always_inline)) INLINE static void hydro_prepare_force(
+    struct part *restrict p, struct xpart *restrict xp) {}
+
+/**
+ * @brief Reset acceleration fields of a particle
+ *
+ * Resets all hydro acceleration and time derivative fields in preparation
+ * for the sums taking place in the variaous force tasks
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
+    struct part *restrict p) {
+
+  for (int i = 0; i < 256; ++i) p->ids_ngbs_force[i] = -1;
+  p->num_ngb_force = 0;
+
+  p->force.h_dt = 0.f;
+}
+
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param p The particle.
+ * @param xp The extended data of this particle.
+ */
+__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values(
+    struct part *restrict p, const struct xpart *restrict xp) {}
+
+/**
+ * @brief Predict additional particle fields forward in time when drifting
+ *
+ * @param p The particle
+ * @param xp The extended data of the particle
+ * @param dt The drift time-step.
+ * @param t0 The time at the start of the drift
+ * @param t1 The time at the end of the drift
+ * @param timeBase The minimal time-step size
+ */
+__attribute__((always_inline)) INLINE static void hydro_predict_extra(
+    struct part *restrict p, const struct xpart *restrict xp, float dt) {}
+
+/**
+ * @brief Finishes the force calculation.
+ *
+ * Multiplies the forces and accelerationsby the appropiate constants
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_end_force(
+    struct part *restrict p) {}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * @param p The particle to act upon
+ * @param xp The particle extended data to act upon
+ * @param dt The time-step for this kick
+ */
+__attribute__((always_inline)) INLINE static void hydro_kick_extra(
+    struct part *restrict p, struct xpart *restrict xp, float dt) {}
+
+/**
+ *  @brief Converts hydro quantity of a particle at the start of a run
+ *
+ * Requires the density to be known
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
+    struct part *restrict p, struct xpart *restrict xp) {}
+
+/**
+ * @brief Initialises the particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * @param p The particle to act upon
+ * @param xp The extended particle data to act upon
+ */
+__attribute__((always_inline)) INLINE static void hydro_first_init_part(
+    struct part *restrict p, struct xpart *restrict xp) {
+
+  p->time_bin = 0;
+  xp->v_full[0] = p->v[0];
+  xp->v_full[1] = p->v[1];
+  xp->v_full[2] = p->v[2];
+
+  hydro_reset_acceleration(p);
+  hydro_init_part(p, NULL);
+}
+
+#endif /* SWIFT_DEBUG_INTERACTIONS_HYDRO_H */

src/hydro/DebugInteractions/hydro_debug.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * 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_DEBUG_INTERACTIONS_HYDRO_DEBUG_H
+#define SWIFT_DEBUG_INTERACTIONS_HYDRO_DEBUG_H
+
+/**
+ * @file DebugInteractions/hydro_debug.h
+ * @brief Empty SPH implementation used solely to test the SELF/PAIR routines.
+ */
+
+__attribute__((always_inline)) INLINE static void hydro_debug_particle(
+    const struct part* p, const struct xpart* xp) {
+  printf(
+      "x=[%.3e,%.3e,%.3e], "
+      "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
+      "h=%.3e, wcount=%.3f, wcount_dh=%.3e, dh/dt=%.3e time_bin=%d\n",
+      p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
+      xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
+      p->h, p->density.wcount, p->density.wcount_dh, p->force.h_dt,
+      p->time_bin);
+}
+
+#endif /* SWIFT_DEBUG_INTERACTIONS_HYDRO_DEBUG_H */

src/hydro/DebugInteractions/hydro_iact.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * 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_DEBUG_INTERACTIONS_HYDRO_IACT_H
+#define SWIFT_DEBUG_INTERACTIONS_HYDRO_IACT_H
+
+/**
+ * @file DebugInteractions/hydro_iact.h
+ * @brief Empty SPH implementation used solely to test the SELF/PAIR routines.
+ */
+
+/**
+ * @brief Density loop
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_density(
+    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+  float wi, wi_dx;
+  float wj, wj_dx;
+
+  /* Get r and r inverse. */
+  const float r = sqrtf(r2);
+
+  /* Compute the kernel function for pi */
+  const float hi_inv = 1.f / hi;
+  const float ui = r * hi_inv;
+  kernel_deval(ui, &wi, &wi_dx);
+
+  /* Compute contribution to the number of neighbours */
+  pi->density.wcount += wi;
+  pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+  /* Compute the kernel function for pj */
+  const float hj_inv = 1.f / hj;
+  const float uj = r * hj_inv;
+  kernel_deval(uj, &wj, &wj_dx);
+
+  /* Compute contribution to the number of neighbours */
+  pj->density.wcount += wj;
+  pj->density.wcount_dh -= (hydro_dimension * wj + uj * wj_dx);
+
+  /* Update ngb counters */
+  pi->ids_ngbs_density[pi->num_ngb_density] = pj->id;
+  ++pi->num_ngb_density;
+  pj->ids_ngbs_density[pj->num_ngb_density] = pi->id;
+  ++pj->num_ngb_density;
+}
+
+/**
+ * @brief Density loop (non-symmetric version)
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_density(
+    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+  float wi, wi_dx;
+
+  /* Get r and r inverse. */
+  const float r = sqrtf(r2);
+
+  /* Compute the kernel function */
+  const float hi_inv = 1.0f / hi;
+  const float ui = r * hi_inv;
+  kernel_deval(ui, &wi, &wi_dx);
+
+  /* Compute contribution to the number of neighbours */
+  pi->density.wcount += wi;
+  pi->density.wcount_dh -= (hydro_dimension * wi + ui * wi_dx);
+
+  /* Update ngb counters */
+  pi->ids_ngbs_density[pi->num_ngb_density] = pj->id;
+  ++pi->num_ngb_density;
+}
+
+/**
+ * @brief Force loop
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_force(
+    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+  /* Update ngb counters */
+  pi->ids_ngbs_force[pi->num_ngb_force] = pj->id;
+  ++pi->num_ngb_force;
+  pj->ids_ngbs_force[pj->num_ngb_force] = pi->id;
+  ++pj->num_ngb_force;
+}
+
+/**
+ * @brief Force loop (non-symmetric version)
+ */
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
+    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
+
+  /* Update ngb counters */
+  pi->ids_ngbs_force[pi->num_ngb_force] = pj->id;
+  ++pi->num_ngb_force;
+}
+
+#endif /* SWIFT_DEBUG_INTERACTIONS_HYDRO_IACT_H */

src/hydro/DebugInteractions/hydro_io.h

+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2017 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * 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_DEBUG_INTERACTIONS_HYDRO_IO_H
+#define SWIFT_DEBUG_INTERACTIONS_HYDRO_IO_H
+
+/**
+ * @file DebugInteractions/hydro_io.h
+ * @brief Empty SPH implementation used solely to test the SELF/PAIR routines.
+ */
+
+#include "adiabatic_index.h"
+#include "hydro.h"
+#include "io_properties.h"
+#include "kernel_hydro.h"
+
+/**
+ * @brief Specifies which particle fields to read from a dataset
+ *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to read.
+ * @param num_fields The number of i/o fields to read.
+ */
+void hydro_read_particles(struct part* parts, struct io_props* list,
+                          int* num_fields) {
+
+  *num_fields = 4;
+
+  /* List what we want to read */
+  list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
+                                UNIT_CONV_LENGTH, parts, x);
+  list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
+                                UNIT_CONV_SPEED, parts, v);
+  list[2] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
+                                UNIT_CONV_LENGTH, parts, h);
+  list[3] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
+                                UNIT_CONV_NO_UNITS, parts, id);
+}
+
+float convert_u(struct engine* e, struct part* p) {
+
+  return hydro_get_internal_energy(p);
+}
+
+float convert_P(struct engine* e, struct part* p) {
+
+  return hydro_get_pressure(p);
+}
+
+/**
+ * @brief Specifies which particle fields to write to a dataset
+ *