diff --git a/src/gravity/Default/gravity_io.h b/src/gravity/Default/gravity_io.h index 26f7a6d15e11cf2e2d9ccd9978a1c67bab17d852..06d8099a88b7c7753c4eabbf7345ad0a59ebb792 100644 --- a/src/gravity/Default/gravity_io.h +++ b/src/gravity/Default/gravity_io.h @@ -35,6 +35,32 @@ void convert_gpart_pos(const struct engine* e, const struct gpart* gp, } } +void convert_gpart_vel(const struct engine* e, const struct gpart* gp, + 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, gp->time_bin); + const integertime_t ti_end = get_integer_time_end(ti_current, gp->time_bin); + + /* Get time-step since the last kick */ + float dt_kick_grav; + 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); + } else { + dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base; + } + + ret[0] = gp->v_full[0] + gp->a_grav[0] * dt_kick_grav; + ret[1] = gp->v_full[1] + gp->a_grav[1] * dt_kick_grav; + ret[2] = gp->v_full[2] + gp->a_grav[2] * dt_kick_grav; +} + /** * @brief Specifies which g-particle fields to read from a dataset * @@ -69,20 +95,20 @@ void darkmatter_read_particles(struct gpart* gparts, struct io_props* list, void darkmatter_write_particles(const struct gpart* gparts, struct io_props* list, int* num_fields) { - /* Say how much we want to read */ + /* Say how much we want to write */ *num_fields = 5; - /* List what we want to read */ + /* List what we want to write */ list[0] = io_make_output_field_convert_gpart( "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, gparts, convert_gpart_pos); - list[1] = io_make_output_field("Velocities", FLOAT, 3, UNIT_CONV_SPEED, - gparts, v_full); + list[1] = io_make_output_field_convert_gpart( + "Velocities", FLOAT, 3, UNIT_CONV_SPEED, gparts, convert_gpart_vel); list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, gparts, mass); list[3] = io_make_output_field("ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, gparts, id_or_neg_offset); - list[4] = io_make_output_field("Acceleration", FLOAT, 3, - UNIT_CONV_ACCELERATION, gparts, a_grav); + list[4] = io_make_output_field("Potential", FLOAT, 1, UNIT_CONV_POTENTIAL, + gparts, potential); } #endif /* SWIFT_DEFAULT_GRAVITY_IO_H */