Updated i/o routines to be compatible with cosmology

src/hydro/PressureEnergy/hydro_io.h

@@ -68,18 +68,20 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
                                 UNIT_CONV_DENSITY, parts, rho);
 }
 
-void convert_S(const struct engine* e, const struct part* p, float* ret) {
+void convert_u(const struct engine* e, const struct part* p,
+	       const struct xpart* xp, float* ret) {
 
   ret[0] = hydro_get_comoving_entropy(p);
 }
 
-void convert_P(const struct engine* e, const struct part* p, float* ret) {
+void convert_P(const struct engine* e, const struct part* p, 
+	       const struct xpart* xp, float* ret) {
 
   ret[0] = hydro_get_comoving_pressure(p);
 }
 
 void convert_part_pos(const struct engine* e, const struct part* p,
-                      double* ret) {
+                      const struct xpart* xp, double* ret) {
 
   if (e->s->periodic) {
     ret[0] = box_wrap(p->x[0], 0.0, e->s->dim[0]);
@@ -92,6 +94,40 @@ void convert_part_pos(const struct engine* e, const struct part* p,
   }
 }
 
+void convert_part_vel(const struct engine* e, const struct part* p,
+                      const struct xpart* xp, float* ret) {
+
+  const int with_cosmology = (e->policy & engine_policy_cosmology);
+  const struct cosmology* cosmo = e->cosmology;
+  const integertime_t ti_current = e->ti_current;
+  const double time_base = e->time_base;
+
+  const integertime_t ti_beg = get_integer_time_begin(ti_current, p->time_bin);
+  const integertime_t ti_end = get_integer_time_end(ti_current, p->time_bin);
+
+  /* Get time-step since the last kick */
+  float dt_kick_grav, dt_kick_hydro;
+  if (with_cosmology) {
+    dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+    dt_kick_grav -=
+        cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+    dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, ti_beg, ti_current);
+    dt_kick_hydro -=
+        cosmology_get_hydro_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+  } else {
+    dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+    dt_kick_hydro = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+  }
+
+  /* Extrapolate the velocites to the current time */
+  hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
+
+  /* Conversion from internal units to peculiar velocities */
+  ret[0] *= cosmo->a2_inv;
+  ret[1] *= cosmo->a2_inv;
+  ret[2] *= cosmo->a2_inv;
+}
+
 /**
  * @brief Specifies which particle fields to write to a dataset
  *
@@ -99,28 +135,29 @@ void convert_part_pos(const struct engine* e, const struct part* p,
  * @param list The list of i/o properties to write.
  * @param num_fields The number of i/o fields to write.
  */
-void hydro_write_particles(struct part* parts, struct io_props* list,
-                           int* num_fields) {
+void hydro_write_particles(const struct part* parts,
+		           const struct xpart* xparts,
+		           struct io_props* list, int* num_fields) {
 
   *num_fields = 9;
 
   /* List what we want to write */
-  list[0] = io_make_output_field_convert_part(
-      "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, parts, convert_part_pos);
-  list[1] =
-      io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, v);
+  list[0] = io_make_output_field_convert_part("Coordinates", DOUBLE, 3,
+		                              UNIT_CONV_LENGTH, parts, xparts,
+					      convert_part_pos);
+  list[1] = io_make_output_field_convert_part(
+      "Velocities", FLOAT, 3, UNIT_CONV_SPEED, parts, xparts, convert_part_vel);
   list[2] =
       io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, parts, mass);
   list[3] = io_make_output_field("SmoothingLength", FLOAT, 1, UNIT_CONV_LENGTH,
                                  parts, h);
-  list[4] = io_make_output_field("InternalEnergy", FLOAT, 1,
-                                 UNIT_CONV_ENERGY_PER_UNIT_MASS, parts, u);
+  list[4] = io_make_output_field_convert_part(
+      "InternalEnergy", FLOAT, 1, UNIT_CONV_ENERGY_PER_UNIT_MASS, parts,
+      xparts, convert_u);
   list[5] = io_make_output_field("ParticleIDs", ULONGLONG, 1,
                                  UNIT_CONV_NO_UNITS, parts, id);
   list[6] =
       io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
-  list[7] = io_make_output_field_convert_part(
-      "Entropy", FLOAT, 1, UNIT_CONV_ENTROPY_PER_UNIT_MASS, parts, convert_S);
   list[8] = io_make_output_field(
       "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, pressure_bar);
 }