Merge branch 'add-sphenix-output-#647' into 'master'

Added new outputs to SPHENIX scheme

Closes #647

See merge request !1032

src/hydro/SPHENIX/hydro.h

@@ -774,6 +774,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
 
   /* Set our old div_v to the one for the next loop */
   p->viscosity.div_v_previous_step = p->viscosity.div_v;
+  p->viscosity.div_v_dt = div_v_dt;
 
   /* Now for the diffusive alpha */
 
@@ -1046,6 +1047,7 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
   p->viscosity.alpha = hydro_props->viscosity.alpha;
   /* Initialise this here to keep all the AV variables together */
   p->viscosity.div_v_previous_step = 0.f;
+  p->viscosity.div_v_dt = 0.f;
 
   /* Set the initial values for the thermal diffusion */
   p->diffusion.alpha = hydro_props->diffusion.alpha;

src/hydro/SPHENIX/hydro_io.h

@@ -140,7 +140,7 @@ INLINE static void convert_part_potential(const struct engine* e,
 INLINE static void convert_viscosity(const struct engine* e,
                                      const struct part* p,
                                      const struct xpart* xp, float* ret) {
-  ret[0] = p->viscosity.alpha;
+  ret[0] = p->viscosity.alpha * p->force.balsara;
 }
 
 INLINE static void convert_diffusion(const struct engine* e,
@@ -161,7 +161,7 @@ INLINE static void hydro_write_particles(const struct part* parts,
                                          struct io_props* list,
                                          int* num_fields) {
 
-  *num_fields = 11;
+  *num_fields = 14;
   /* List what we want to write */
   list[0] = io_make_output_field_convert_part(
       "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, 1.f, parts, xparts,
@@ -203,11 +203,34 @@ INLINE static void hydro_write_particles(const struct part* parts,
 
   list[9] = io_make_output_field_convert_part(
       "ViscosityParameters", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, parts, xparts,
-      convert_viscosity, "Visosity coefficient (alpha_visc) of the particles");
+      convert_viscosity,
+      "Visosity coefficient (alpha_visc) of the particles, multiplied by the "
+      "balsara switch");
 
   list[10] = io_make_output_field_convert_part(
       "DiffusionParameters", FLOAT, 1, UNIT_CONV_NO_UNITS, 0.f, parts, xparts,
       convert_diffusion, "Diffusion coefficient (alpha_diff) of the particles");
+
+  list[11] = io_make_output_field(
+      "LaplacianInternalEnergies", FLOAT, 1, UNIT_CONV_FREQUENCY_SQUARED,
+      1.f - 3.f * hydro_gamma, parts, diffusion.laplace_u,
+      "Laplacian (del squared) of the Internal Energy per "
+      "unit mass of the particles");
+
+  list[12] = io_make_output_field(
+      "VelocityDivergences", FLOAT, 1, UNIT_CONV_FREQUENCY, 0.f, parts,
+      viscosity.div_v,
+      "Local velocity divergence field around the particles. Provided without "
+      "cosmology, as this includes the Hubble flow. To return to a peculiar "
+      "velocity divergence, div . v_pec = a^2 (div . v - n_D H)");
+
+  list[13] = io_make_output_field(
+      "VelocityDivergenceTimeDifferentials", FLOAT, 1,
+      UNIT_CONV_FREQUENCY_SQUARED, 0.f, parts, viscosity.div_v_dt,
+      "Time differential (over the previous step) of the "
+      "velocity divergence field around the particles. Again, provided without "
+      "cosmology as this includes a Hubble flow term. To get back to a peculiar
+      velocity divergence time differential, x_pec = a^4 (x - a^{-2} n_D dH / dt)");
 }
 
 /**

src/hydro/SPHENIX/hydro_part.h

@@ -118,6 +118,9 @@ struct part {
     /*! Particle velocity divergence */
     float div_v;
 
+    /*! Time differential of velocity divergence */
+    float div_v_dt;
+
     /*! Particle velocity divergence from previous step */
     float div_v_previous_step;
 

src/units.c

@@ -258,6 +258,10 @@ void units_get_base_unit_exponents_array(float baseUnitsExp[5],
       baseUnitsExp[UNIT_TIME] = -1.f;
       break;
 
+    case UNIT_CONV_FREQUENCY_SQUARED:
+      baseUnitsExp[UNIT_TIME] = -2.f;
+      break;
+
     case UNIT_CONV_DENSITY:
       baseUnitsExp[UNIT_MASS] = 1.f;
       baseUnitsExp[UNIT_LENGTH] = -3.f;

src/units.h

@@ -86,6 +86,7 @@ enum unit_conversion_factor {
   UNIT_CONV_POWER,
   UNIT_CONV_PRESSURE,
   UNIT_CONV_FREQUENCY,
+  UNIT_CONV_FREQUENCY_SQUARED,
   UNIT_CONV_ELECTRIC_CHARGE,
   UNIT_CONV_ELECTRIC_VOLTAGE,
   UNIT_CONV_ELECTRIC_CAPACITANCE,