Cleaned up material for merge

src/adiabatic_index.h

@@ -432,40 +432,4 @@ __attribute__((always_inline)) INLINE static float pow_one_over_gamma(float x) {
 #endif
 }
 
-/**
- * @brief Return the argument to the power one minus two over the adiabatic
- * index
- *
- * Computes \f$x^{1 - \frac{2}{\gamma}}\f$.
- *
- * @param x Argument
- * @return Argument to the power one minus two over the adiabatic index
- */
-__attribute__((always_inline)) INLINE static float pow_one_minus_two_over_gamma(
-    float x) {
-
-#if defined(HYDRO_GAMMA_5_3)
-
-  return powf(x, -0.2f);
-
-#elif defined(HYDRO_GAMMA_7_5)
-
-  return powf(x, -0.428571429f);
-
-#elif defined(HYDRO_GAMMA_4_3)
-
-  return 1.f / sqrtf(x);
-
-#elif defined(HYDRO_GAMMA_2_1)
-
-  return 1.f;
-
-#else
-
-  error("The adiabatic index is not defined !");
-  return 0.f;
-
-#endif
-}
-
 #endif /* SWIFT_ADIABATIC_INDEX_H */

src/const.h

@@ -40,9 +40,9 @@
 #define const_isothermal_internal_energy 20.2615290634f
 
 /* Dimensionality of the problem */
-//#define HYDRO_DIMENSION_3D
+#define HYDRO_DIMENSION_3D
 //#define HYDRO_DIMENSION_2D
-#define HYDRO_DIMENSION_1D
+//#define HYDRO_DIMENSION_1D
 
 /* Hydrodynamical adiabatic index. */
 #define HYDRO_GAMMA_5_3
@@ -64,8 +64,8 @@
 
 /* SPH variant to use */
 //#define MINIMAL_SPH
-//#define GADGET2_SPH
-#define HOPKINS_PE_SPH
+#define GADGET2_SPH
+//#define HOPKINS_PE_SPH
 //#define DEFAULT_SPH
 //#define GIZMO_SPH
 

src/hydro/PressureEntropy/hydro_io.h

@@ -56,9 +56,6 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
                                 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,
-   * entropy); */
   list[4] =
       io_make_input_field("InternalEnergy", FLOAT, 1, COMPULSORY,
                           UNIT_CONV_ENTROPY_PER_UNIT_MASS, parts, entropy);
@@ -127,8 +124,11 @@ void writeSPHflavour(hid_t h_grpsph) {
   /* Viscosity and thermal conduction */
   /* Nothing in this minimal model... */
   writeAttribute_s(h_grpsph, "Thermal Conductivity Model", "No treatment");
-  writeAttribute_s(h_grpsph, "Viscosity Model",
-                   "Minimal treatment as in Monaghan (1992)");
+  writeAttribute_s(
+      h_grpsph, "Viscosity Model",
+      "as in Springel (2005), i.e. Monaghan (1992) with Balsara (1995) switch");
+  writeAttribute_f(h_grpsph, "Viscosity alpha", const_viscosity_alpha);
+  writeAttribute_f(h_grpsph, "Viscosity beta", 3.f);
 
   /* Time integration properties */
   writeAttribute_f(h_grpsph, "Maximal Delta u change over dt",

tests/test125cells.c

@@ -95,6 +95,8 @@ void set_energy_state(struct part *part, enum pressure_field press, float size,
 
 #if defined(GADGET2_SPH)
   part->entropy = pressure / pow_gamma(density);
+#elif defined(HOPKINS_PE_SPH)
+  part->entropy = pressure / pow_gamma(density);
 #elif defined(DEFAULT_SPH)
   part->u = pressure / (hydro_gamma_minus_one * density);
 #elif defined(MINIMAL_SPH)

tests/test27cells.c

@@ -104,11 +104,18 @@ struct cell *make_cell(size_t n, double *offset, double size, double h,
         }
         part->h = size * h / (float)n;
         part->id = ++(*partId);
+
 #ifdef GIZMO_SPH
         part->conserved.mass = density * volume / count;
 #else
         part->mass = density * volume / count;
 #endif
+
+#if defined(HOPKINS_PE_SPH)
+        part->entropy = 1.f;
+        part->entropy_one_over_gamma = 1.f;
+#endif
+
         part->ti_begin = 0;
         part->ti_end = 1;
         ++part;
@@ -192,17 +199,14 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
             hydro_get_density(&main_cell->parts[pid]),
 #if defined(GIZMO_SPH)
             0.f,
+#elif defined(HOPKINS_PE_SPH)
+            main_cell->parts[pid].density.rho_dh,
 #else
             main_cell->parts[pid].rho_dh,
 #endif
             main_cell->parts[pid].density.wcount,
             main_cell->parts[pid].density.wcount_dh,
-#if defined(GADGET2_SPH)
-            main_cell->parts[pid].density.div_v,
-            main_cell->parts[pid].density.rot_v[0],
-            main_cell->parts[pid].density.rot_v[1],
-            main_cell->parts[pid].density.rot_v[2]
-#elif defined(DEFAULT_SPH)
+#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH)
             main_cell->parts[pid].density.div_v,
             main_cell->parts[pid].density.rot_v[0],
             main_cell->parts[pid].density.rot_v[1],
@@ -234,14 +238,13 @@ void dump_particle_fields(char *fileName, struct cell *main_cell,
               cj->parts[pjd].v[2], hydro_get_density(&cj->parts[pjd]),
 #if defined(GIZMO_SPH)
               0.f,
+#elif defined(HOPKINS_PE_SPH)
+              main_cell->parts[pjd].density.rho_dh,
 #else
               main_cell->parts[pjd].rho_dh,
 #endif
               cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh,
-#if defined(GADGET2_SPH)
-              cj->parts[pjd].density.div_v, cj->parts[pjd].density.rot_v[0],
-              cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2]
-#elif defined(DEFAULT_SPH)
+#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH)
               cj->parts[pjd].density.div_v, cj->parts[pjd].density.rot_v[0],
               cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2]
 #else

tests/testPair.c

@@ -138,11 +138,13 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
             ci->parts[pid].v[2], hydro_get_density(&ci->parts[pid]),
 #if defined(GIZMO_SPH)
             0.f,
+#elif defined(HOPKINS_PE_SPH)
+            ci->parts[pid].density.rho_dh,
 #else
-            cj->parts[pid].rho_dh,
+            ci->parts[pid].rho_dh,
 #endif
             ci->parts[pid].density.wcount, ci->parts[pid].density.wcount_dh,
-#if defined(GADGET2_SPH) || defined(DEFAULT_SPH)
+#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH)
             ci->parts[pid].density.div_v, ci->parts[pid].density.rot_v[0],
             ci->parts[pid].density.rot_v[1], ci->parts[pid].density.rot_v[2]
 #else
@@ -162,11 +164,13 @@ void dump_particle_fields(char *fileName, struct cell *ci, struct cell *cj) {
             cj->parts[pjd].v[2], hydro_get_density(&cj->parts[pjd]),
 #if defined(GIZMO_SPH)
             0.f,
+#elif defined(HOPKINS_PE_SPH)
+            cj->parts[pjd].density.rho_dh,
 #else
             cj->parts[pjd].rho_dh,
 #endif
             cj->parts[pjd].density.wcount, cj->parts[pjd].density.wcount_dh,
-#if defined(GADGET2_SPH) || defined(DEFAULT_SPH)
+#if defined(GADGET2_SPH) || defined(DEFAULT_SPH) || defined(HOPKINS_PE_SPH)
             cj->parts[pjd].density.div_v, cj->parts[pjd].density.rot_v[0],
             cj->parts[pjd].density.rot_v[1], cj->parts[pjd].density.rot_v[2]
 #else