diff --git a/src/hydro/PressureEnergy/hydro_io.h b/src/hydro/PressureEnergy/hydro_io.h index 0d06f6fde0ba4f82fd3b37ff70fb0b5d5aa80ade..2bd919b8ddab52adb28ca3490b549d1e467b991a 100644 --- a/src/hydro/PressureEnergy/hydro_io.h +++ b/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); }