Pass the xparts to the i/o convertion function. Extrapolate the part...

Pass the xparts to the i/o convertion function. Extrapolate the part velocities to the current time when writting them out.

src/common_io.c

@@ -457,6 +457,7 @@ void io_convert_part_f_mapper(void* restrict temp, int N,
 
   const struct io_props props = *((const struct io_props*)extra_data);
   const struct part* restrict parts = props.parts;
+  const struct xpart* restrict xparts = props.xparts;
   const struct engine* e = props.e;
   const size_t dim = props.dimension;
 
@@ -465,7 +466,8 @@ void io_convert_part_f_mapper(void* restrict temp, int N,
   const ptrdiff_t delta = (temp_f - props.start_temp_f) / dim;
 
   for (int i = 0; i < N; i++)
-    props.convert_part_f(e, parts + delta + i, &temp_f[i * dim]);
+    props.convert_part_f(e, parts + delta + i, xparts + delta + i,
+                         &temp_f[i * dim]);
 }
 
 /**
@@ -477,6 +479,7 @@ void io_convert_part_d_mapper(void* restrict temp, int N,
 
   const struct io_props props = *((const struct io_props*)extra_data);
   const struct part* restrict parts = props.parts;
+  const struct xpart* restrict xparts = props.xparts;
   const struct engine* e = props.e;
   const size_t dim = props.dimension;
 
@@ -485,7 +488,8 @@ void io_convert_part_d_mapper(void* restrict temp, int N,
   const ptrdiff_t delta = (temp_d - props.start_temp_d) / dim;
 
   for (int i = 0; i < N; i++)
-    props.convert_part_d(e, parts + delta + i, &temp_d[i * dim]);
+    props.convert_part_d(e, parts + delta + i, xparts + delta + i,
+                         &temp_d[i * dim]);
 }
 
 /**

src/hydro/Gadget2/hydro_io.h

@@ -55,18 +55,20 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
                                 UNIT_CONV_DENSITY, parts, rho);
 }
 
-void convert_u(const struct engine* e, const struct part* p, float* ret) {
+void convert_part_u(const struct engine* e, const struct part* p,
+                    const struct xpart* xp, float* ret) {
 
   ret[0] = hydro_get_comoving_internal_energy(p);
 }
 
-void convert_P(const struct engine* e, const struct part* p, float* ret) {
+void convert_part_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]);
@@ -79,6 +81,34 @@ 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;
+  }
+
+  hydro_get_drifted_velocities(p, xp, dt_kick_hydro, dt_kick_grav, ret);
+}
+
 /**
  * @brief Specifies which particle fields to write to a dataset
  *
@@ -86,16 +116,17 @@ 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(const 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,
@@ -108,9 +139,9 @@ void hydro_write_particles(const struct part* parts, struct io_props* list,
       io_make_output_field("Density", FLOAT, 1, UNIT_CONV_DENSITY, parts, rho);
   list[7] = io_make_output_field_convert_part("InternalEnergy", FLOAT, 1,
                                               UNIT_CONV_ENERGY_PER_UNIT_MASS,
-                                              parts, convert_u);
+                                              parts, xparts, convert_part_u);
   list[8] = io_make_output_field_convert_part(
-      "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, convert_P);
+      "Pressure", FLOAT, 1, UNIT_CONV_PRESSURE, parts, xparts, convert_part_P);
 
 #ifdef DEBUG_INTERACTIONS_SPH
 

src/io_properties.h

@@ -34,9 +34,11 @@ enum DATA_IMPORTANCE { COMPULSORY = 1, OPTIONAL = 0, UNUSED = -1 };
 
 /* Helper typedefs */
 typedef void (*conversion_func_part_float)(const struct engine*,
-                                           const struct part*, float*);
+                                           const struct part*,
+                                           const struct xpart*, float*);
 typedef void (*conversion_func_part_double)(const struct engine*,
-                                            const struct part*, double*);
+                                            const struct part*,
+                                            const struct xpart*, double*);
 typedef void (*conversion_func_gpart_float)(const struct engine*,
                                             const struct gpart*, float*);
 typedef void (*conversion_func_gpart_double)(const struct engine*,
@@ -78,6 +80,7 @@ struct io_props {
 
   /* The particle arrays */
   const struct part* parts;
+  const struct xpart* xparts;
   const struct gpart* gparts;
 
   /* Are we converting? */
@@ -125,6 +128,7 @@ INLINE static struct io_props io_make_input_field_(
   r.field = field;
   r.partSize = partSize;
   r.parts = NULL;
+  r.xparts = NULL;
   r.gparts = NULL;
   r.conversion = 0;
   r.convert_part_f = NULL;
@@ -179,10 +183,10 @@ INLINE static struct io_props io_make_output_field_(
 /**
  * @brief Constructs an #io_props (with conversion) from its parameters
  */
-#define io_make_output_field_convert_part(name, type, dim, units, part, \
-                                          convert)                      \
-  io_make_output_field_convert_part_##type(name, type, dim, units,      \
-                                           sizeof(part[0]), part, convert)
+#define io_make_output_field_convert_part(name, type, dim, units, part, xpart, \
+                                          convert)                             \
+  io_make_output_field_convert_part_##type(                                    \
+      name, type, dim, units, sizeof(part[0]), part, xpart, convert)
 
 /**
  * @brief Construct an #io_props from its parameters
@@ -193,6 +197,7 @@ INLINE static struct io_props io_make_output_field_(
  * @param units The units of the dataset
  * @param partSize The size in byte of the particle
  * @param parts The particle array
+ * @param xparts The xparticle array
  * @param functionPtr The function used to convert a particle to a float
  *
  * Do not call this function directly. Use the macro defined above.
@@ -200,7 +205,8 @@ INLINE static struct io_props io_make_output_field_(
 INLINE static struct io_props io_make_output_field_convert_part_FLOAT(
     const char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension,
     enum unit_conversion_factor units, size_t partSize,
-    const struct part* parts, conversion_func_part_float functionPtr) {
+    const struct part* parts, const struct xpart* xparts,
+    conversion_func_part_float functionPtr) {
 
   struct io_props r;
   strcpy(r.name, name);
@@ -211,6 +217,7 @@ INLINE static struct io_props io_make_output_field_convert_part_FLOAT(
   r.field = NULL;
   r.partSize = partSize;
   r.parts = parts;
+  r.xparts = xparts;
   r.gparts = NULL;
   r.conversion = 1;
   r.convert_part_f = functionPtr;
@@ -230,6 +237,7 @@ INLINE static struct io_props io_make_output_field_convert_part_FLOAT(
  * @param units The units of the dataset
  * @param partSize The size in byte of the particle
  * @param parts The particle array
+ * @param xparts The xparticle array
  * @param functionPtr The function used to convert a particle to a float
  *
  * Do not call this function directly. Use the macro defined above.
@@ -237,7 +245,8 @@ INLINE static struct io_props io_make_output_field_convert_part_FLOAT(
 INLINE static struct io_props io_make_output_field_convert_part_DOUBLE(
     const char name[FIELD_BUFFER_SIZE], enum IO_DATA_TYPE type, int dimension,
     enum unit_conversion_factor units, size_t partSize,
-    const struct part* parts, conversion_func_part_double functionPtr) {
+    const struct part* parts, const struct xpart* xparts,
+    conversion_func_part_double functionPtr) {
 
   struct io_props r;
   strcpy(r.name, name);
@@ -248,6 +257,7 @@ INLINE static struct io_props io_make_output_field_convert_part_DOUBLE(
   r.field = NULL;
   r.partSize = partSize;
   r.parts = parts;
+  r.xparts = xparts;
   r.gparts = NULL;
   r.conversion = 1;
   r.convert_part_f = NULL;
@@ -293,6 +303,7 @@ INLINE static struct io_props io_make_output_field_convert_gpart_FLOAT(
   r.field = NULL;
   r.partSize = gpartSize;
   r.parts = NULL;
+  r.xparts = NULL;
   r.gparts = gparts;
   r.conversion = 1;
   r.convert_part_f = NULL;
@@ -330,6 +341,7 @@ INLINE static struct io_props io_make_output_field_convert_gpart_DOUBLE(
   r.field = NULL;
   r.partSize = gpartSize;
   r.parts = NULL;
+  r.xparts = NULL;
   r.gparts = gparts;
   r.conversion = 1;
   r.convert_part_f = NULL;

src/parallel_io.c

@@ -821,9 +821,10 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
                   const struct unit_system* internal_units,
                   const struct unit_system* snapshot_units) {
 
-  struct part* parts = e->s->parts;
-  struct gpart* gparts = e->s->gparts;
-  struct spart* sparts = e->s->sparts;
+  const struct part* parts = e->s->parts;
+  const struct xpart* xparts = e->s->xparts;
+  const struct gpart* gparts = e->s->gparts;
+  const struct spart* sparts = e->s->sparts;
   FILE* xmfFile = 0;
   int periodic = e->s->periodic;
   int numFiles = 1;
@@ -980,7 +981,7 @@ void prepare_file(struct engine* e, const char* baseName, long long N_total[6],
     switch (ptype) {
 
       case swift_type_gas:
-        hydro_write_particles(parts, list, &num_fields);
+        hydro_write_particles(parts, xparts, list, &num_fields);
         num_fields += chemistry_write_particles(parts, list + num_fields);
         break;
 
@@ -1045,10 +1046,11 @@ void write_output_parallel(struct engine* e, const char* baseName,
   const size_t Ngas = e->s->nr_parts;
   const size_t Nstars = e->s->nr_sparts;
   const size_t Ntot = e->s->nr_gparts;
-  struct part* parts = e->s->parts;
-  struct gpart* gparts = e->s->gparts;
+  const struct part* parts = e->s->parts;
+  const struct xpart* xparts = e->s->xparts;
+  const struct gpart* gparts = e->s->gparts;
   struct gpart* dmparts = NULL;
-  struct spart* sparts = e->s->sparts;
+  const struct spart* sparts = e->s->sparts;
 
   /* Number of unassociated gparts */
   const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0;
@@ -1208,7 +1210,7 @@ void write_output_parallel(struct engine* e, const char* baseName,
 
       case swift_type_gas:
         Nparticles = Ngas;
-        hydro_write_particles(parts, list, &num_fields);
+        hydro_write_particles(parts, xparts, list, &num_fields);
         num_fields += chemistry_write_particles(parts, list + num_fields);
         break;
 

src/serial_io.c

@@ -712,10 +712,11 @@ void write_output_serial(struct engine* e, const char* baseName,
   const size_t Ntot = e->s->nr_gparts;
   int periodic = e->s->periodic;
   int numFiles = 1;
-  struct part* parts = e->s->parts;
-  struct gpart* gparts = e->s->gparts;
+  const struct part* parts = e->s->parts;
+  const struct xpart* xparts = e->s->xparts;
+  const struct gpart* gparts = e->s->gparts;
   struct gpart* dmparts = NULL;
-  struct spart* sparts = e->s->sparts;
+  const struct spart* sparts = e->s->sparts;
   FILE* xmfFile = 0;
 
   /* Number of unassociated gparts */
@@ -962,7 +963,7 @@ void write_output_serial(struct engine* e, const char* baseName,
 
           case swift_type_gas:
             Nparticles = Ngas;
-            hydro_write_particles(parts, list, &num_fields);
+            hydro_write_particles(parts, xparts, list, &num_fields);
             num_fields += chemistry_write_particles(parts, list + num_fields);
             break;
 

src/single_io.c

@@ -578,10 +578,11 @@ void write_output_single(struct engine* e, const char* baseName,
   const size_t Ntot = e->s->nr_gparts;
   int periodic = e->s->periodic;
   int numFiles = 1;
-  struct part* parts = e->s->parts;
-  struct gpart* gparts = e->s->gparts;
+  const struct part* parts = e->s->parts;
+  const struct xpart* xparts = e->s->xparts;
+  const struct gpart* gparts = e->s->gparts;
   struct gpart* dmparts = NULL;
-  struct spart* sparts = e->s->sparts;
+  const struct spart* sparts = e->s->sparts;
 
   /* Number of unassociated gparts */
   const size_t Ndm = Ntot > 0 ? Ntot - (Ngas + Nstars) : 0;
@@ -779,7 +780,7 @@ void write_output_single(struct engine* e, const char* baseName,
 
       case swift_type_gas:
         N = Ngas;
-        hydro_write_particles(parts, list, &num_fields);
+        hydro_write_particles(parts, xparts, list, &num_fields);
         num_fields += chemistry_write_particles(parts, list + num_fields);
         break;
 

src/stars/Default/star_io.h

@@ -52,7 +52,7 @@ void star_read_particles(struct spart* sparts, struct io_props* list,
  * @param list The list of i/o properties to write.
  * @param num_fields The number of i/o fields to write.
  */
-void star_write_particles(struct spart* sparts, struct io_props* list,
+void star_write_particles(const struct spart* sparts, struct io_props* list,
                           int* num_fields) {
 
   /* Say how much we want to read */