Merge branch 'output_bh_gas_velocities' into 'master'

Output the velocities of the gas surrounding the BHs

See merge request !1063

src/black_holes/EAGLE/black_holes.h

@@ -234,6 +234,10 @@ black_holes_bpart_has_no_neighbours(struct bpart* bp,
   /* Re-set problematic values */
   bp->density.wcount = kernel_root * h_inv_dim;
   bp->density.wcount_dh = 0.f;
+
+  bp->velocity_gas[0] = FLT_MAX;
+  bp->velocity_gas[1] = FLT_MAX;
+  bp->velocity_gas[2] = FLT_MAX;
 }
 
 /**
@@ -387,7 +391,7 @@ __attribute__((always_inline)) INLINE static void black_holes_prepare_feedback(
     const struct phys_const* constants, const struct cosmology* cosmo,
     const double time, const int with_cosmology, const double dt) {
 
-  if (dt == 0.) return;
+  if (dt == 0. || bp->rho_gas == 0.) return;
 
   /* Gather some physical constants (all in internal units) */
   const double G = constants->const_newton_G;

src/black_holes/EAGLE/black_holes_io.h

@@ -100,6 +100,37 @@ INLINE static void convert_bpart_vel(const struct engine* e,
   ret[2] *= cosmo->a_inv;
 }
 
+INLINE static void convert_bpart_gas_vel(const struct engine* e,
+                                         const struct bpart* bp, float* ret) {
+
+  const struct cosmology* cosmo = e->cosmology;
+
+  /* Convert velocities to peculiar velocities */
+  const double gas_v_peculiar[3] = {bp->velocity_gas[0] * cosmo->a_inv,
+                                    bp->velocity_gas[1] * cosmo->a_inv,
+                                    bp->velocity_gas[2] * cosmo->a_inv};
+
+  const double bh_v_peculiar[3] = {bp->v[0] * cosmo->a_inv,
+                                   bp->v[1] * cosmo->a_inv,
+                                   bp->v[2] * cosmo->a_inv};
+
+  ret[0] = gas_v_peculiar[0] - bh_v_peculiar[0];
+  ret[1] = gas_v_peculiar[1] - bh_v_peculiar[1];
+  ret[2] = gas_v_peculiar[2] - bh_v_peculiar[2];
+}
+
+INLINE static void convert_bpart_gas_circular_vel(const struct engine* e,
+                                                  const struct bpart* bp,
+                                                  float* ret) {
+
+  const struct cosmology* cosmo = e->cosmology;
+
+  /* Conversion from internal units to peculiar velocities */
+  ret[0] = bp->circular_velocity_gas[0] * cosmo->a_inv;
+  ret[1] = bp->circular_velocity_gas[1] * cosmo->a_inv;
+  ret[2] = bp->circular_velocity_gas[2] * cosmo->a_inv;
+}
+
 /**
  * @brief Specifies which b-particle fields to write to a dataset
  *
@@ -114,7 +145,7 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
                                                int with_cosmology) {
 
   /* Say how much we want to write */
-  *num_fields = 18;
+  *num_fields = 20;
 
   /* List what we want to write */
   list[0] = io_make_output_field_convert_bpart(
@@ -232,6 +263,20 @@ INLINE static void black_holes_write_particles(const struct bpart* bparts,
       "Physical angular momenta that the black holes have accumulated by "
       "swallowing gas particles.");
 
+  list[18] = io_make_output_field_convert_bpart(
+      "GasRelativeVelocities", FLOAT, 3, UNIT_CONV_SPEED, 0.f, bparts,
+      convert_bpart_gas_vel,
+      "Peculiar relative velocities of the gas particles around the black "
+      "holes. This is a * dx/dt where x is the co-moving position of the "
+      "particles.");
+
+  list[19] = io_make_output_field_convert_bpart(
+      "GasCircularVelocities", FLOAT, 3, UNIT_CONV_SPEED, 0.f, bparts,
+      convert_bpart_gas_circular_vel,
+      "Peculiar circular velocities of the gas particles around the black "
+      "holes. This is the curl of a * dx/dt where x is the co-moving position "
+      "of the particles.");
+
 #ifdef DEBUG_INTERACTIONS_BLACK_HOLES
 
   list += *num_fields;