Added GIZMO_SPH option in const.h. Made sure the resulting code compiles. It...

Added GIZMO_SPH option in const.h. Made sure the resulting code compiles. It probably does not work correctly though.

src/const.h

@@ -51,8 +51,9 @@
 
 /* SPH variant to use */
 //#define MINIMAL_SPH
-#define GADGET2_SPH
+//#define GADGET2_SPH
 //#define DEFAULT_SPH
+#define GIZMO_SPH
 
 /* Self gravity stuff. */
 #define const_gravity_multipole_order 2

src/debug.c

@@ -43,6 +43,8 @@
 #include "./hydro/Gadget2/hydro_debug.h"
 #elif defined(DEFAULT_SPH)
 #include "./hydro/Default/hydro_debug.h"
+#elif defined(GIZMO_SPH)
+#include "./hydro/Gizmo/hydro_debug.h"
 #else
 #error "Invalid choice of SPH variant"
 #endif

src/hydro.h

@@ -38,6 +38,10 @@
 #include "./hydro/Default/hydro.h"
 #include "./hydro/Default/hydro_iact.h"
 #define SPH_IMPLEMENTATION "Default version of SPH"
+#elif defined(GIZMO_SPH)
+#include "./hydro/Gizmo/hydro.h"
+#include "./hydro/Gizmo/hydro_iact.h"
+#define SPH_IMPLEMENTATION "Mesh-free hydrodynamics as in GIZMO (Hopkins 2015)"
 #else
 #error "Invalid choice of SPH variant"
 #endif

src/hydro/Gizmo/hydro.h

@@ -17,6 +17,9 @@
  *
  ******************************************************************************/
 
+#include <float.h>
+#include "adiabatic_index.h"
+
 /**
  * @brief Computes the hydro time-step of a given particle
  *
@@ -25,9 +28,12 @@
  *
  */
 __attribute__((always_inline)) INLINE static float hydro_compute_timestep(
-    struct part* p, struct xpart* xp) {
+    const struct part* restrict p, const struct xpart* restrict xp,
+    const struct hydro_props* restrict hydro_properties) {
 
-  return const_cfl * p->h / fabs(p->timestepvars.vmax);
+  const float CFL_condition = hydro_properties->CFL_condition;
+
+  return CFL_condition * p->h / fabsf(p->timestepvars.vmax);
 }
 
 /**
@@ -59,7 +65,7 @@ __attribute__((always_inline)) INLINE static void hydro_init_part(
     p->primitives.v[1] = p->conserved.momentum[1] / p->conserved.mass;
     p->primitives.v[2] = p->conserved.momentum[2] / p->conserved.mass;
     p->primitives.P =
-        (const_hydro_gamma - 1.) *
+        hydro_gamma_minus_one *
         (p->conserved.energy -
          0.5 * (p->conserved.momentum[0] * p->conserved.momentum[0] +
                 p->conserved.momentum[1] * p->conserved.momentum[1] +
@@ -158,15 +164,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
 
   float detE, volume;
   float E[3][3];
-  GFLOAT m, momentum[3], energy;
+  float m, momentum[3], energy;
 
 #ifndef THERMAL_ENERGY
-  GFLOAT momentum2;
+  float momentum2;
 #endif
 
 #if defined(SPH_GRADIENTS) && defined(SLOPE_LIMITER)
-  GFLOAT gradrho[3], gradv[3][3], gradP[3];
-  GFLOAT gradtrue, gradmax, gradmin, alpha;
+  float gradrho[3], gradv[3][3], gradP[3];
+  float gradtrue, gradmax, gradmin, alpha;
 #endif
 
   /* Final operation on the geometry. */
@@ -374,9 +380,9 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_gradient(
     p->primitives.v[2] = momentum[2] / m;
 #ifndef THERMAL_ENERGY
     p->primitives.P =
-        (const_hydro_gamma - 1.) * (energy - 0.5 * momentum2 / m) / volume;
+        hydro_gamma_minus_one * (energy - 0.5 * momentum2 / m) / volume;
 #else
-    p->primitives.P = (const_hydro_gamma - 1.) * energy / volume;
+    p->primitives.P = hydro_gamma_minus_one * energy / volume;
 #endif
   }
 }
@@ -392,8 +398,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_gradient(
 #ifndef SPH_GRADIENTS
   float h, ih, ih2, ih3;
 #ifdef SLOPE_LIMITER
-  GFLOAT gradrho[3], gradv[3][3], gradP[3];
-  GFLOAT gradtrue, gradmax, gradmin, alpha;
+  float gradrho[3], gradv[3][3], gradP[3];
+  float gradtrue, gradmax, gradmin, alpha;
 #endif
 
   /* add kernel normalization to gradients */
@@ -575,10 +581,10 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
     struct part* p) {
 
   float volume;
-  GFLOAT m;
-  GFLOAT momentum[3];
+  float m;
+  float momentum[3];
 #ifndef THERMAL_ENERGY
-  GFLOAT momentum2;
+  float momentum2;
 #endif
   volume = p->geometry.volume;
 
@@ -587,7 +593,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
   p->primitives.v[1] = p->v[1];
   p->primitives.v[2] = p->v[2];
   /* P actually contains internal energy at this point */
-  p->primitives.P *= (const_hydro_gamma - 1.) * p->primitives.rho;
+  p->primitives.P *= hydro_gamma_minus_one * p->primitives.rho;
 
   p->conserved.mass = m = p->primitives.rho * volume;
   p->conserved.momentum[0] = momentum[0] = m * p->primitives.v[0];
@@ -597,9 +603,9 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
   momentum2 = momentum[0] * momentum[0] + momentum[1] * momentum[1] +
               momentum[2] * momentum[2];
   p->conserved.energy =
-      p->primitives.P / (const_hydro_gamma - 1.) * volume + 0.5 * momentum2 / m;
+      p->primitives.P / hydro_gamma_minus_one * volume + 0.5 * momentum2 / m;
 #else
-  p->conserved.energy = p->primitives.P / (const_hydro_gamma - 1.) * volume;
+  p->conserved.energy = p->primitives.P / hydro_gamma_minus_one * volume;
 #endif
 }
 
@@ -614,7 +620,51 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
     struct part* p) {}
 __attribute__((always_inline)) INLINE static void hydro_kick_extra(
     struct part* p, struct xpart* xp, float dt, float half_dt) {}
+
+/**
+ * @brief Returns the internal energy of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
 __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
-    struct part* p) {
-  return 0.f;
+    const struct part* restrict p, float dt) {
+
+  return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
+}
+
+/**
+ * @brief Returns the entropy of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_entropy(
+    const struct part* restrict p, float dt) {
+
+  return p->primitives.P / pow_gamma(p->primitives.rho);
+}
+
+/**
+ * @brief Returns the sound speed of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
+    const struct part* restrict p, float dt) {
+
+  return sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
+}
+
+/**
+ * @brief Returns the pressure of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_pressure(
+    const struct part* restrict p, float dt) {
+
+  return p->primitives.P;
 }

src/hydro/Gizmo/hydro_debug.h

@@ -18,7 +18,7 @@
  ******************************************************************************/
 
 __attribute__((always_inline)) INLINE static void hydro_debug_particle(
-    struct part* p, struct xpart* xp) {
+    const struct part* p, const struct xpart* xp) {
   printf(
       "x=[%.16e,%.16e,%.16e], "
       "v=[%.3e,%.3e,%.3e], a=[%.3e,%.3e,%.3e], volume=%.3e\n",

src/hydro/Gizmo/hydro_iact.h

@@ -19,6 +19,7 @@
  *
  ******************************************************************************/
 
+#include "adiabatic_index.h"
 #include "riemann.h"
 
 #define USE_GRADIENTS
@@ -26,7 +27,7 @@
 /* #define PRINT_ID 0 */
 
 /* this corresponds to task_subtype_hydro_loop1 */
-__attribute__((always_inline)) INLINE static void runner_iact_hydro_loop1(
+__attribute__((always_inline)) INLINE static void runner_iact_density(
     float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
 
   float r = sqrtf(r2);
@@ -194,9 +195,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_hydro_loop1(
 }
 
 /* this corresponds to task_subtype_hydro_loop1 */
-__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_hydro_loop1(float r2, float *dx, float hi, float hj,
-                               struct part *pi, struct part *pj) {
+__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 r;
   float xi;
@@ -296,7 +296,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_hydro_loop2(
   int k, l;
   float Bi[3][3];
   float Bj[3][3];
-  GFLOAT Wi[5], Wj[5];
+  float Wi[5], Wj[5];
 
   /* Initialize local variables */
   for (k = 0; k < 3; k++) {
@@ -497,7 +497,7 @@ runner_iact_nonsym_hydro_loop2(float r2, float *dx, float hi, float hj,
   float wi, wi_dx;
   int k, l;
   float Bi[3][3];
-  GFLOAT Wi[5], Wj[5];
+  float Wi[5], Wj[5];
 
   /* Initialize local variables */
   for (k = 0; k < 3; k++) {
@@ -619,13 +619,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   float xij_i[3], xfac, xijdotdx;
   float vmax, dvdotdx;
   float vi[3], vj[3], vij[3];
-  GFLOAT Wi[5], Wj[5];  //, Whalf[5];
+  float Wi[5], Wj[5];  //, Whalf[5];
 #ifdef USE_GRADIENTS
-  GFLOAT dWi[5], dWj[5];
+  float dWi[5], dWj[5];
   float xij_j[3];
 #endif
-  //    GFLOAT rhoe;
-  //    GFLOAT flux[5][3];
+  //    float rhoe;
+  //    float flux[5][3];
   float dti, dtj, mindt;
   float n_unit[3];
 
@@ -658,8 +658,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   //    }
 
   /* calculate the maximal signal velocity */
-  vmax = sqrtf(const_hydro_gamma * Wi[4] / Wi[0]) +
-         sqrtf(const_hydro_gamma * Wj[4] / Wj[0]);
+  vmax =
+      sqrtf(hydro_gamma * Wi[4] / Wi[0]) + sqrtf(hydro_gamma * Wj[4] / Wj[0]);
   dvdotdx = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
             (Wi[3] - Wj[3]) * dx[2];
   if (dvdotdx > 0.) {
@@ -790,14 +790,14 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
 #ifdef PER_FACE_LIMITER
 
   float xij_i_norm;
-  GFLOAT phi_i, phi_j;
-  GFLOAT delta1, delta2;
-  GFLOAT phiminus, phiplus;
-  GFLOAT phimin, phimax;
-  GFLOAT phibar;
+  float phi_i, phi_j;
+  float delta1, delta2;
+  float phiminus, phiplus;
+  float phimin, phimax;
+  float phibar;
   /* free parameters, values from Hopkins */
-  GFLOAT psi1 = 0.5, psi2 = 0.25;
-  GFLOAT phi_mid0, phi_mid;
+  float psi1 = 0.5, psi2 = 0.25;
+  float phi_mid0, phi_mid;
 
   for (k = 0; k < 10; k++) {
     if (k < 5) {
@@ -877,13 +877,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
                            Wi[2] * pi->primitives.gradients.v[2][1] +
                            Wi[3] * pi->primitives.gradients.v[2][2] +
                            pi->primitives.gradients.P[2] / Wi[0]);
-  dWi[4] -= 0.5 * mindt *
-            (Wi[1] * pi->primitives.gradients.P[0] +
-             Wi[2] * pi->primitives.gradients.P[1] +
-             Wi[3] * pi->primitives.gradients.P[2] +
-             const_hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
-                                          pi->primitives.gradients.v[1][1] +
-                                          pi->primitives.gradients.v[2][2]));
+  dWi[4] -=
+      0.5 * mindt * (Wi[1] * pi->primitives.gradients.P[0] +
+                     Wi[2] * pi->primitives.gradients.P[1] +
+                     Wi[3] * pi->primitives.gradients.P[2] +
+                     hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
+                                            pi->primitives.gradients.v[1][1] +
+                                            pi->primitives.gradients.v[2][2]));
 
   dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
                            Wj[2] * pj->primitives.gradients.rho[1] +
@@ -903,13 +903,13 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
                            Wj[2] * pj->primitives.gradients.v[2][1] +
                            Wj[3] * pj->primitives.gradients.v[2][2] +
                            pj->primitives.gradients.P[2] / Wj[0]);
-  dWj[4] -= 0.5 * mindt *
-            (Wj[1] * pj->primitives.gradients.P[0] +
-             Wj[2] * pj->primitives.gradients.P[1] +
-             Wj[3] * pj->primitives.gradients.P[2] +
-             const_hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
-                                          pj->primitives.gradients.v[1][1] +
-                                          pj->primitives.gradients.v[2][2]));
+  dWj[4] -=
+      0.5 * mindt * (Wj[1] * pj->primitives.gradients.P[0] +
+                     Wj[2] * pj->primitives.gradients.P[1] +
+                     Wj[3] * pj->primitives.gradients.P[2] +
+                     hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
+                                            pj->primitives.gradients.v[1][1] +
+                                            pj->primitives.gradients.v[2][2]));
 
   //    printf("WL: %g %g %g %g %g\n", Wi[0], Wi[1], Wi[2], Wi[3], Wi[4]);
   //    printf("WR: %g %g %g %g %g\n", Wj[0], Wj[1], Wj[2], Wj[3], Wj[4]);
@@ -969,7 +969,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
     error("Negative density or pressure!\n");
   }
 
-  GFLOAT totflux[5];
+  float totflux[5];
   riemann_solve_for_flux(Wi, Wj, n_unit, vij, totflux);
 
   /* Update conserved variables */
@@ -1014,16 +1014,32 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
 }
 
 /* this corresponds to task_subtype_fluxes */
-__attribute__((always_inline)) INLINE static void runner_iact_hydro_loop3(
+__attribute__((always_inline)) INLINE static void runner_iact_force(
     float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
 
   runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 1);
 }
 
 /* this corresponds to task_subtype_fluxes */
-__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_hydro_loop3(float r2, float *dx, float hi, float hj,
-                               struct part *pi, struct part *pj) {
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
+    float r2, float *dx, float hi, float hj, struct part *pi, struct part *pj) {
 
   runner_iact_fluxes_common(r2, dx, hi, hj, pi, pj, 0);
 }
+
+/* PLACEHOLDERS */
+__attribute__((always_inline)) INLINE static void runner_iact_vec_density(
+    float *R2, float *Dx, float *Hi, float *Hj, struct part **pi,
+    struct part **pj) {}
+
+__attribute__((always_inline)) INLINE static void
+runner_iact_nonsym_vec_density(float *R2, float *Dx, float *Hi, float *Hj,
+                               struct part **pi, struct part **pj) {}
+
+__attribute__((always_inline)) INLINE static void runner_iact_vec_force(
+    float *R2, float *Dx, float *Hi, float *Hj, struct part **pi,
+    struct part **pj) {}
+
+__attribute__((always_inline)) INLINE static void runner_iact_nonsym_vec_force(
+    float *R2, float *Dx, float *Hi, float *Hj, struct part **pi,
+    struct part **pj) {}

src/hydro/Gizmo/hydro_io.h

 /*******************************************************************************
  * This file is part of SWIFT.
- * Coypright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ * Coypright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
  *
  * 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
@@ -17,77 +17,75 @@
  *
  ******************************************************************************/
 
+#include "adiabatic_index.h"
+#include "io_properties.h"
+
 /**
- * @brief Reads the different particles to the HDF5 file
- *
- * @param h_grp The HDF5 group in which to read the arrays.
- * @param N The number of particles on that MPI rank.
- * @param N_total The total number of particles (only used in MPI mode)
- * @param offset The offset of the particles for this MPI rank (only used in MPI
- *mode)
- * @param parts The particle array
+ * @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.
  */
-__attribute__((always_inline)) INLINE static void hydro_read_particles(
-    hid_t h_grp, int N, long long N_total, long long offset,
-    struct part* parts) {
+void hydro_read_particles(struct part* parts, struct io_props* list,
+                          int* num_fields) {
+
+  *num_fields = 8;
 
-  /* Read arrays */
-  readArray(h_grp, "Coordinates", DOUBLE, N, 3, parts, N_total, offset, x,
-            COMPULSORY);
-  readArray(h_grp, "Velocities", FLOAT, N, 3, parts, N_total, offset, v,
-            COMPULSORY);
-  readArray(h_grp, "Masses", FLOAT, N, 1, parts, N_total, offset,
-            conserved.mass, COMPULSORY);
-  readArray(h_grp, "SmoothingLength", FLOAT, N, 1, parts, N_total, offset, h,
-            COMPULSORY);
-  readArray(h_grp, "InternalEnergy", FLOAT, N, 1, parts, N_total, offset,
-            primitives.P, COMPULSORY);
-  readArray(h_grp, "ParticleIDs", ULONGLONG, N, 1, parts, N_total, offset, id,
-            COMPULSORY);
-  readArray(h_grp, "Acceleration", FLOAT, N, 3, parts, N_total, offset, a_hydro,
-            OPTIONAL);
-  readArray(h_grp, "Density", FLOAT, N, 1, parts, N_total, offset,
-            primitives.rho, OPTIONAL);
+  /* 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("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
+                                parts, mass);
+  list[3] = io_make_input_field("SmoothingLength", FLOAT, 1, COMPULSORY,
+                                UNIT_CONV_LENGTH, parts, h);
+  list[4] =
+      io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
+                          UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, primitives.P);
+  list[5] = io_make_input_field("ParticleIDs", ULONGLONG, 1, COMPULSORY,
+                                UNIT_CONV_NO_UNITS, parts, id);
+  list[6] = io_make_input_field("Accelerations", FLOAT, 3, OPTIONAL,
+                                UNIT_CONV_ACCELERATION, parts, a_hydro);
+  list[7] = io_make_input_field("Density", FLOAT, 1, OPTIONAL,
+                                UNIT_CONV_DENSITY, parts, primitives.rho);
+}
+
+float convert_u(struct engine* e, struct part* p) {
+  return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
 }
 
 /**
- * @brief Writes the different particles to the HDF5 file
- *
- * @param h_grp The HDF5 group in which to write the arrays.
- * @param fileName The name of the file (unsued in MPI mode).
- * @param xmfFile The XMF file to write to (unused in MPI mode).
- * @param N The number of particles on that MPI rank.
- * @param N_total The total number of particles (only used in MPI mode)
- * @param mpi_rank The MPI rank of this node (only used in MPI mode)
- * @param offset The offset of the particles for this MPI rank (only used in MPI
- *mode)
- * @param parts The particle array
- * @param us The unit system to use
+ * @brief Specifies which particle fields to write to a dataset
  *
+ * @param parts The particle array.
+ * @param list The list of i/o properties to write.
+ * @param num_fields The number of i/o fields to write.
  */
-__attribute__((always_inline)) INLINE static void hydro_write_particles(
-    hid_t h_grp, char* fileName, FILE* xmfFile, int N, long long N_total,
-    int mpi_rank, long long offset, struct part* parts, struct UnitSystem* us) {
+void hydro_write_particles(struct part* parts, struct io_props* list,
+                           int* num_fields) {
 
-  /* Write arrays */
-  writeArray(h_grp, fileName, xmfFile, "Coordinates", DOUBLE, N, 3, parts,
-             N_total, mpi_rank, offset, x, us, UNIT_CONV_LENGTH);
-  writeArray(h_grp, fileName, xmfFile, "Velocities", FLOAT, N, 3, parts,
-             N_total, mpi_rank, offset, v, us, UNIT_CONV_SPEED);
-  writeArray(h_grp, fileName, xmfFile, "Masses", FLOAT, N, 1, parts, N_total,
-             mpi_rank, offset, conserved.mass, us, UNIT_CONV_MASS);
-  writeArray(h_grp, fileName, xmfFile, "SmoothingLength", FLOAT, N, 1, parts,
-             N_total, mpi_rank, offset, h, us, UNIT_CONV_LENGTH);
-  writeArray(h_grp, fileName, xmfFile, "InternalEnergy", FLOAT, N, 1, parts,
-             N_total, mpi_rank, offset, primitives.P, us,
-             UNIT_CONV_ENTROPY_PER_UNIT_MASS);
-  writeArray(h_grp, fileName, xmfFile, "ParticleIDs", ULONGLONG, N, 1, parts,
-             N_total, mpi_rank, offset, id, us, UNIT_CONV_NO_UNITS);
-  writeArray(h_grp, fileName, xmfFile, "Acceleration", FLOAT, N, 3, parts,
-             N_total, mpi_rank, offset, a_hydro, us, UNIT_CONV_ACCELERATION);
-  writeArray(h_grp, fileName, xmfFile, "Density", FLOAT, N, 1, parts, N_total,
-             mpi_rank, offset, primitives.rho, us, UNIT_CONV_DENSITY);
+  *num_fields = 8;
+
+  /* List what we want to write */
+  list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
+                                 parts, x);
+  list[1] = io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts,
+                                 primitives.v);
+  list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts,
+                                 conserved.mass);
+  list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
+                                 parts, h);
+  list[4] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
+                                              UNIT_CONV_ENERGY_PER_UNIT_MASS,
+                                              parts, primitives.P, convert_u);
+  list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
+                                 UNIT_CONV_NO_UNITS, parts, id);
+  list[6] = io_make_output_field("Acceleration", FLOAT, 3,
+                                 UNIT_CONV_ACCELERATION, parts, a_hydro);
+  list[7] = io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts,
+                                 primitives.rho);
 }
 
 /**
@@ -96,25 +94,27 @@ __attribute__((always_inline)) INLINE static void hydro_write_particles(
  */
 void writeSPHflavour(hid_t h_grpsph) {
 
-  /* Kernel function description */
-  writeAttribute_s(h_grpsph, "Kernel", kernel_name);
-  writeAttribute_f(h_grpsph, "Kernel eta", const_eta_kernel);
-  writeAttribute_f(h_grpsph, "Weighted N_ngb", kernel_nwneigh);
-  writeAttribute_f(h_grpsph, "Delta N_ngb", const_delta_nwneigh);
-  writeAttribute_f(h_grpsph, "Hydro gamma", const_hydro_gamma);
-
   /* Viscosity and thermal conduction */
   writeAttribute_s(h_grpsph, "Thermal Conductivity Model",
-                   "(No treatment) Legacy Gadget-2 as in Springel (2005)");
+                   "Price (2008) without switch");
+  writeAttribute_f(h_grpsph, "Thermal Conductivity alpha",
+                   const_conductivity_alpha);
   writeAttribute_s(h_grpsph, "Viscosity Model",
-                   "Legacy Gadget-2 as in Springel (2005)");
-  writeAttribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
-  writeAttribute_f(h_grpsph, "Viscosity beta", 3.f);
+                   "Morris & Monaghan (1997), Rosswog, Davies, Thielemann & "
+                   "Piran (2000) with additional Balsara (1995) switch");
+  writeAttribute_f(h_grpsph, "Viscosity alpha_min", const_viscosity_alpha_min);
+  writeAttribute_f(h_grpsph, "Viscosity alpha_max", const_viscosity_alpha_max);
+  writeAttribute_f(h_grpsph, "Viscosity beta", 2.f);
+  writeAttribute_f(h_grpsph, "Viscosity decay length", const_viscosity_length);
 
   /* Time integration properties */
-  writeAttribute_f(h_grpsph, "CFL parameter", const_cfl);
-  writeAttribute_f(h_grpsph, "Maximal ln(Delta h) change over dt",
-                   const_ln_max_h_change);
-  writeAttribute_f(h_grpsph, "Maximal Delta h change over dt",
-                   exp(const_ln_max_h_change));
+  writeAttribute_f(h_grpsph, "Maximal Delta u change over dt",
+                   const_max_u_change);
 }
+
+/**
+ * @brief Are we writing entropy in the internal energy field ?
+ *
+ * @return 1 if entropy is in 'internal energy', 0 otherwise.
+ */
+int writeEntropyFlag() { return 0; }

src/hydro/Gizmo/hydro_part.h

@@ -19,23 +19,15 @@
  *
  ******************************************************************************/
 
-/* Some standard headers. */
-#include <stdlib.h>
-
-#define GFLOAT float
-
 /* Extra particle data not needed during the computation. */
 struct xpart {
 
-  /* Old position, at last tree rebuild. */
-  double x_old[3];
+  /* Offset between current position and position at last tree rebuild. */
+  float x_diff[3];
 
   /* Velocity at the last full step. */
   float v_full[3];
 
-  /* Entropy at the half-step. */
-  float u_hdt;
-
   /* Old density. */
   float omega;
 
@@ -47,17 +39,12 @@ struct part {
   /* Particle position. */
   double x[3];
 
-  /* Particle velocity. */
+  /* Particle predicted velocity. */
   float v[3];
 
   /* Particle acceleration. */
   float a_hydro[3];
 
-  float mass;  // MATTHIEU
-  float h_dt;
-  float rho;
-  float rho_dh;
-
   /* Particle cutoff radius. */
   float h;
 
@@ -67,76 +54,86 @@ struct part {
   /* Particle time of end of time-step. */
   int ti_end;
 
-  /* The primitive hydrodynamical variables */
+  /* The primitive hydrodynamical variables. */
   struct {
 
-    /* fluid velocity */
-    GFLOAT v[3];
+    /* Fluid velocity. */
+    float v[3];
 
-    /* density */
-    GFLOAT rho;
+    /* Density. */
+    float rho;
 
-    /* pressure */
-    GFLOAT P;
+    /* Pressure. */
+    float P;
 
+    /* Gradients of the primitive variables. */
     struct {
 
-      GFLOAT rho[3];
+      /* Density gradients. */
+      float rho[3];
 
-      GFLOAT v[3][3];
+      /* Fluid velocity gradients. */
+      float v[3][3];
 
-      GFLOAT P[3];
+      /* Pressure gradients. */
+      float P[3];
 
     } gradients;
 
+    /* Quantities needed by the slope limiter. */
     struct {
 
-      /* extreme values among the neighbours */
-      GFLOAT rho[2];
+      /* Extreme values of the density among the neighbours. */
+      float rho[2];
 
-      GFLOAT v[3][2];
+      /* Extreme values of the fluid velocity among the neighbours. */
+      float v[3][2];
 
-      GFLOAT P[2];
+      /* Extreme values of the pressure among the neighbours. */
+      float P[2];
 
-      /* maximal distance to all neighbouring faces */
+      /* Maximal distance to all neighbouring faces. */
       float maxr;
 
     } limiter;
 
   } primitives;
 
-  /* The conserved hydrodynamical variables */
+  /* The conserved hydrodynamical variables. */
   struct {
 
-    /* fluid momentum */
-    GFLOAT momentum[3];
+    /* Fluid momentum. */
+    float momentum[3];
 
-    /* fluid mass */
-    GFLOAT mass;
+    /* Fluid mass (this field already exists outside of this struct as well). */
+    float mass;
 
-    /* fluid energy */
-    GFLOAT energy;
+    /* Fluid total energy (= internal energy + fluid kinetic energy). */
+    float energy;
 
   } conserved;
 
-  /* Geometrical quantities used for hydro */
+  /* Geometrical quantities used for hydro. */
   struct {
 
-    /* volume of the particle */
+    /* Volume of the particle. */
     float volume;
 
-    /* gradient matrix */
+    /* Geometrical shear matrix used to calculate second order accurate
+       gradients */
     float matrix_E[3][3];
 
   } geometry;
 
+  /* Variables used for timestep calculation (currently not used). */
   struct {
 
+    /* Maximum fluid velocity among all neighbours. */
     float vmax;
 
   } timestepvars;
 
-  /* Quantities used during the density loop */
+  /* Quantities used during the volume (=density) loop. */
   struct {
 
     /* Particle velocity divergence. */
@@ -153,10 +150,21 @@ struct part {
 
   } density;
 
+  /* Particle mass (this field is also part of the conserved quantities...). */
+  float mass;
+
   /* Particle ID. */
   unsigned long long id;
 
   /* Associated gravitas. */
   struct gpart *gpart;
 
+  /* Variables needed for the code to compile (should be removed/replaced). */
+  float rho;
+  struct {
+    float h_dt;
+    float v_sig;
+  } force;
+  float rho_dh;
+
 } __attribute__((aligned(part_align)));

src/hydro_io.h

@@ -28,6 +28,8 @@
 #include "./hydro/Gadget2/hydro_io.h"
 #elif defined(DEFAULT_SPH)
 #include "./hydro/Default/hydro_io.h"
+#elif defined(GIZMO_SPH)
+#include "./hydro/Gizmo/hydro_io.h"
 #else
 #error "Invalid choice of SPH variant"
 #endif

src/part.h

@@ -45,6 +45,8 @@
 #include "./hydro/Gadget2/hydro_part.h"
 #elif defined(DEFAULT_SPH)
 #include "./hydro/Default/hydro_part.h"
+#elif defined(GIZMO_SPH)
+#include "./hydro/Gizmo/hydro_part.h"
 #else
 #error "Invalid choice of SPH variant"
 #endif

src/riemann/riemann_exact.h

@@ -50,12 +50,11 @@
  * @param W The left or right state vector
  * @param a The left or right sound speed
  */
-__attribute__((always_inline)) INLINE static GFLOAT riemann_fb(GFLOAT p,
-                                                               GFLOAT* W,
-                                                               GFLOAT a) {
+__attribute__((always_inline)) INLINE static float riemann_fb(float p, float* W,
+                                                              float a) {
 
-  GFLOAT fval = 0.;
-  GFLOAT A, B;
+  float fval = 0.;
+  float A, B;
   if (p > W[4]) {
     A = const_riemann_tdgp1 / W[0];
     B = const_riemann_gm1dgp1 * W[4];
@@ -78,9 +77,8 @@ __attribute__((always_inline)) INLINE static GFLOAT riemann_fb(GFLOAT p,
  * @param aL The left sound speed
  * @param aR The right sound speed
  */
-__attribute__((always_inline)) INLINE static GFLOAT riemann_f(
-    GFLOAT p, GFLOAT* WL, GFLOAT* WR, GFLOAT vL, GFLOAT vR, GFLOAT aL,
-    GFLOAT aR) {
+__attribute__((always_inline)) INLINE static float riemann_f(
+    float p, float* WL, float* WR, float vL, float vR, float aL, float aR) {
 
   return riemann_fb(p, WL, aL) + riemann_fb(p, WR, aR) + (vR - vL);
 }
@@ -92,12 +90,12 @@ __attribute__((always_inline)) INLINE static GFLOAT riemann_f(
  * @param W The left or right state vector
  * @param a The left or right sound speed
  */
-__attribute__((always_inline)) INLINE static GFLOAT riemann_fprimeb(GFLOAT p,
-                                                                    GFLOAT* W,
-                                                                    GFLOAT a) {
+__attribute__((always_inline)) INLINE static float riemann_fprimeb(float p,
+                                                                   float* W,
+                                                                   float a) {
 
-  GFLOAT fval = 0.;
-  GFLOAT A, B;
+  float fval = 0.;
+  float A, B;
   if (p > W[4]) {
     A = const_riemann_tdgp1 / W[0];
     B = const_riemann_gm1dgp1 * W[4];
@@ -117,8 +115,8 @@ __attribute__((always_inline)) INLINE static GFLOAT riemann_fprimeb(GFLOAT p,
  * @param aL The left sound speed
  * @param aR The right sound speed
  */
-__attribute__((always_inline)) INLINE static GFLOAT riemann_fprime(
-    GFLOAT p, GFLOAT* WL, GFLOAT* WR, GFLOAT aL, GFLOAT aR) {
+__attribute__((always_inline)) INLINE static float riemann_fprime(
+    float p, float* WL, float* WR, float aL, float aR) {
 
   return riemann_fprimeb(p, WL, aL) + riemann_fprimeb(p, WR, aR);
 }
@@ -129,10 +127,10 @@ __attribute__((always_inline)) INLINE static GFLOAT riemann_fprime(
  * @param p The current guess for the pressure
  * @param W The left or right state vector
  */
-__attribute__((always_inline)) INLINE static GFLOAT riemann_gb(GFLOAT p,
-                                                               GFLOAT* W) {
+__attribute__((always_inline)) INLINE static float riemann_gb(float p,
+                                                              float* W) {
 
-  GFLOAT A, B;
+  float A, B;
   A = const_riemann_tdgp1 / W[0];
   B = const_riemann_gm1dgp1 * W[4];
   return sqrtf(A / (p + B));
@@ -151,11 +149,11 @@ __attribute__((always_inline)) INLINE static GFLOAT riemann_gb(GFLOAT p,
  * @param aL The left sound speed
  * @param aR The right sound speed
  */
-__attribute__((always_inline)) INLINE static GFLOAT riemann_guess_p(
-    GFLOAT* WL, GFLOAT* WR, GFLOAT vL, GFLOAT vR, GFLOAT aL, GFLOAT aR) {
+__attribute__((always_inline)) INLINE static float riemann_guess_p(
+    float* WL, float* WR, float vL, float vR, float aL, float aR) {
 
-  GFLOAT pguess, pmin, pmax, qmax;
-  GFLOAT ppv;
+  float pguess, pmin, pmax, qmax;
+  float ppv;
 
   pmin = fminf(WL[4], WR[4]);
   pmax = fmaxf(WL[4], WR[4]);
@@ -202,14 +200,14 @@ __attribute__((always_inline)) INLINE static GFLOAT riemann_guess_p(
  * @param aL The left sound speed
  * @param aR The right sound speed
  */
-__attribute__((always_inline)) INLINE static GFLOAT riemann_solve_brent(
-    GFLOAT lower_limit, GFLOAT upper_limit, GFLOAT lowf, GFLOAT upf,
-    GFLOAT error_tol, GFLOAT* WL, GFLOAT* WR, GFLOAT vL, GFLOAT vR, GFLOAT aL,
-    GFLOAT aR) {
-
-  GFLOAT a, b, c, d, s;
-  GFLOAT fa, fb, fc, fs;
-  GFLOAT tmp, tmp2;
+__attribute__((always_inline)) INLINE static float riemann_solve_brent(
+    float lower_limit, float upper_limit, float lowf, float upf,
+    float error_tol, float* WL, float* WR, float vL, float vR, float aL,
+    float aR) {
+
+  float a, b, c, d, s;
+  float fa, fb, fc, fs;
+  float tmp, tmp2;
   int mflag;
   int i;
 
@@ -309,11 +307,11 @@ __attribute__((always_inline)) INLINE static GFLOAT riemann_solve_brent(
  * @param n_unit Normal vector of the interface
  */
 __attribute__((always_inline)) INLINE static void riemann_solve_vacuum(
-    GFLOAT* WL, GFLOAT* WR, GFLOAT vL, GFLOAT vR, GFLOAT aL, GFLOAT aR,
-    GFLOAT* Whalf, float* n_unit) {
+    float* WL, float* WR, float vL, float vR, float aL, float aR, float* Whalf,
+    float* n_unit) {
 
-  GFLOAT SL, SR;
-  GFLOAT vhalf;
+  float SL, SR;
+  float vhalf;
 
   if (!WR[0] && !WL[0]) {
     /* if both states are vacuum, the solution is also vacuum */
@@ -456,20 +454,20 @@ __attribute__((always_inline)) INLINE static void riemann_solve_vacuum(
  * @param n_unit Normal vector of the interface
  */
 __attribute__((always_inline)) INLINE static void riemann_solver_solve(
-    GFLOAT* WL, GFLOAT* WR, GFLOAT* Whalf, float* n_unit) {
+    float* WL, float* WR, float* Whalf, float* n_unit) {
 
   /* velocity of the left and right state in a frame aligned with n_unit */
-  GFLOAT vL, vR, vhalf;
+  float vL, vR, vhalf;
   /* sound speeds */
-  GFLOAT aL, aR;
+  float aL, aR;
   /* variables used for finding pstar */
-  GFLOAT p, pguess, fp, fpguess;
+  float p, pguess, fp, fpguess;
   /* variables used for sampling the solution */
-  GFLOAT u;
-  GFLOAT pdpR, SR;
-  GFLOAT SHR, STR;
-  GFLOAT pdpL, SL;
-  GFLOAT SHL, STL;
+  float u;
+  float pdpR, SR;
+  float SHR, STR;
+  float pdpL, SL;
+  float SHL, STL;
   int errorFlag = 0;
 
   /* sanity checks */
@@ -661,10 +659,10 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
 }
 
 __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
-    GFLOAT* Wi, GFLOAT* Wj, float* n_unit, float* vij, GFLOAT* totflux) {
+    float* Wi, float* Wj, float* n_unit, float* vij, float* totflux) {
 
-  GFLOAT Whalf[5];
-  GFLOAT flux[5][3];
+  float Whalf[5];
+  float flux[5][3];
   float vtot[3];
   float rhoe;
 

src/riemann/riemann_hllc.h

@@ -20,13 +20,15 @@
 #ifndef SWIFT_RIEMANN_HLLC_H
 #define SWIFT_RIEMANN_HLLC_H
 
+#include "adiabatic_index.h"
+
 __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
-    GFLOAT *WL, GFLOAT *WR, float *n, float *vij, GFLOAT *totflux) {
+    float *WL, float *WR, float *n, float *vij, float *totflux) {
 
-  GFLOAT uL, uR, aL, aR;
-  GFLOAT rhobar, abar, pPVRS, pstar, qL, qR, SL, SR, Sstar;
-  GFLOAT v2, eL, eR;
-  GFLOAT UstarL[5], UstarR[5];
+  float uL, uR, aL, aR;
+  float rhobar, abar, pPVRS, pstar, qL, qR, SL, SR, Sstar;
+  float v2, eL, eR;
+  float UstarL[5], UstarR[5];
 
   /* Handle pure vacuum */
   if (!WL[0] && !WR[0]) {
@@ -41,14 +43,14 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
   /* STEP 0: obtain velocity in interface frame */
   uL = WL[1] * n[0] + WL[2] * n[1] + WL[3] * n[2];
   uR = WR[1] * n[0] + WR[2] * n[1] + WR[3] * n[2];
-  aL = sqrtf(const_hydro_gamma * WL[4] / WL[0]);
-  aR = sqrtf(const_hydro_gamma * WR[4] / WR[0]);
+  aL = sqrtf(hydro_gamma * WL[4] / WL[0]);
+  aR = sqrtf(hydro_gamma * WR[4] / WR[0]);
 
   /* Handle vacuum: vacuum does not require iteration and is always exact */
   if (!WL[0] || !WR[0]) {
     error("Vacuum not yet supported");
   }
-  if (2. * aL / (const_hydro_gamma - 1.) + 2. * aR / (const_hydro_gamma - 1.) <
+  if (2. * aL / hydro_gamma_minus_one + 2. * aR / hydro_gamma_minus_one <
       fabs(uL - uR)) {
     error("Vacuum not yet supported");
   }
@@ -64,14 +66,12 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
   qL = 1.;
   if (pstar > WL[4]) {
     qL = sqrtf(1. +
-               0.5 * (const_hydro_gamma + 1.) / const_hydro_gamma *
-                   (pstar / WL[4] - 1.));
+               0.5 * (hydro_gamma + 1.) / hydro_gamma * (pstar / WL[4] - 1.));
   }
   qR = 1.;
   if (pstar > WR[4]) {
     qR = sqrtf(1. +
-               0.5 * (const_hydro_gamma + 1.) / const_hydro_gamma *
-                   (pstar / WR[4] - 1.));
+               0.5 * (hydro_gamma + 1.) / hydro_gamma * (pstar / WR[4] - 1.));
   }
   SL = uL - aL * qL;
   SR = uR + aR * qR;
@@ -88,7 +88,7 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
     totflux[2] = WL[0] * WL[2] * uL + WL[4] * n[1];
     totflux[3] = WL[0] * WL[2] * uL + WL[4] * n[2];
     v2 = WL[1] * WL[1] + WL[2] * WL[2] + WL[3] * WL[3];
-    eL = WL[4] / (const_hydro_gamma - 1.) / WL[0] + 0.5 * v2;
+    eL = WL[4] / hydro_gamma_minus_one / WL[0] + 0.5 * v2;
     totflux[4] = WL[0] * eL * uL + WL[4] * uL;
     if (SL < 0.) {
       /* add flux FstarL */
@@ -118,7 +118,7 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
     totflux[2] = WR[0] * WR[2] * uR + WR[4] * n[1];
     totflux[3] = WR[0] * WR[3] * uR + WR[4] * n[2];
     v2 = WR[1] * WR[1] + WR[2] * WR[2] + WR[3] * WR[3];
-    eR = WR[4] / (const_hydro_gamma - 1.) / WR[0] + 0.5 * v2;
+    eR = WR[4] / hydro_gamma_minus_one / WR[0] + 0.5 * v2;
     totflux[4] = WR[0] * eR * uR + WR[4] * uR;
     if (SR > 0.) {
       /* add flux FstarR */

src/riemann/riemann_trrs.h

@@ -50,14 +50,14 @@
  * @param n_unit Normal vector of the interface
  */
 __attribute__((always_inline)) INLINE static void riemann_solver_solve(
-    GFLOAT* WL, GFLOAT* WR, GFLOAT* Whalf, float* n_unit) {
-  GFLOAT aL, aR;
-  GFLOAT PLR;
-  GFLOAT vL, vR;
-  GFLOAT ustar, pstar;
-  GFLOAT vhalf;
-  GFLOAT pdpR, SHR, STR;
-  GFLOAT pdpL, SHL, STL;
+    float* WL, float* WR, float* Whalf, float* n_unit) {
+  float aL, aR;
+  float PLR;
+  float vL, vR;
+  float ustar, pstar;
+  float vhalf;
+  float pdpR, SHR, STR;
+  float pdpL, SHL, STL;
 
   /* calculate the velocities along the interface normal */
   vL = WL[1] * n_unit[0] + WL[2] * n_unit[1] + WL[3] * n_unit[2];

src/timestep.h

@@ -108,7 +108,7 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
         e->external_potential, e->physical_constants, p->gpart);
     /* const float new_dt_self = */
     /*     gravity_compute_timestep_self(e->physical_constants, p->gpart); */
-    const float new_dt_self = FLT_MAX; // MATTHIEU
+    const float new_dt_self = FLT_MAX;  // MATTHIEU
 
     new_dt_grav = fminf(new_dt_external, new_dt_self);
   }

tests/test125cells.c

@@ -99,6 +99,8 @@ void set_energy_state(struct part *part, enum pressure_field press, float size,
   part->u = pressure / (hydro_gamma_minus_one * density);
 #elif defined(MINIMAL_SPH)
   part->u = pressure / (hydro_gamma_minus_one * density);
+#elif defined(GIZMO_SPH)
+  part->primitives.P = pressure;
 #else
   error("Need to define pressure here !");
 #endif