diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h index beb6f98b8c0d781aa709fb6ee3ca564a52704db2..801e864709bbf20e20b476d93e3ce69f47a4a1b0 100644 --- a/src/hydro/Minimal/hydro.h +++ b/src/hydro/Minimal/hydro.h @@ -210,26 +210,6 @@ __attribute__((always_inline)) INLINE static float hydro_compute_timestep( return dt_cfl; } -/** - * @brief Initialises the particles for the first time - * - * This function is called only once just after the ICs have been - * read in to do some conversions or assignments between the particle - * and extended particle fields. - * - * @param p The particle to act upon - * @param xp The extended particle data to act upon - */ -__attribute__((always_inline)) INLINE static void hydro_first_init_part( - struct part *restrict p, struct xpart *restrict xp) { - - p->ti_begin = 0; - p->ti_end = 0; - xp->v_full[0] = p->v[0]; - xp->v_full[1] = p->v[1]; - xp->v_full[2] = p->v[2]; -} - /** * @brief Prepares a particle for the density calculation. * @@ -292,16 +272,12 @@ __attribute__((always_inline)) INLINE static void hydro_end_density( * * @param p The particle to act upon * @param xp The extended particle data to act upon - * @param ti_current The current time (on the timeline) - * @param timeBase The minimal time-step size */ __attribute__((always_inline)) INLINE static void hydro_prepare_force( - struct part *restrict p, struct xpart *restrict xp, int ti_current, - double timeBase) { + struct part *restrict p, struct xpart *restrict xp) { /* Compute the pressure */ - const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase; - const float pressure = hydro_get_pressure(p, half_dt); + const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); /* Compute the sound speed */ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); @@ -339,6 +315,25 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( p->force.v_sig = 0.0f; } +/** + * @brief Sets the values to be predicted in the drifts to their values at a + * kick time + * + * @param p The particle. + * @param xp The extended data of this particle. + */ +__attribute__((always_inline)) INLINE static void hydro_reset_predicted_values( + struct part *restrict p, const struct xpart *restrict xp) { + + /* Re-set the predicted velocities */ + p->v[0] = xp->v_full[0]; + p->v[1] = xp->v_full[1]; + p->v[2] = xp->v_full[2]; + + /* Re-set the entropy */ + p->u = xp->u_full; +} + /** * @brief Predict additional particle fields forward in time when drifting * @@ -348,13 +343,9 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration( * @param p The particle. * @param xp The extended data of the particle. * @param dt The drift time-step. - * @param t0 The time at the start of the drift (on the timeline). - * @param t1 The time at the end of the drift (on the timeline). - * @param timeBase The minimal time-step size. */ __attribute__((always_inline)) INLINE static void hydro_predict_extra( - struct part *restrict p, const struct xpart *restrict xp, float dt, int t0, - int t1, double timeBase) { + struct part *restrict p, const struct xpart *restrict xp, float dt) { const float h_inv = 1.f / p->h; @@ -372,9 +363,11 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( else p->rho *= expf(w2); - /* Drift the pressure */ - const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase; - const float pressure = hydro_get_pressure(p, dt_entr); + /* Predict the internal energy */ + p->u += p->u_dt * dt; + + /* Compute the new pressure */ + const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); /* Compute the new sound speed */ const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure); @@ -407,24 +400,19 @@ __attribute__((always_inline)) INLINE static void hydro_end_force( * @param p The particle to act upon * @param xp The particle extended data to act upon * @param dt The time-step for this kick - * @param half_dt The half time-step for this kick */ __attribute__((always_inline)) INLINE static void hydro_kick_extra( - struct part *restrict p, struct xpart *restrict xp, float dt, - float half_dt) { + struct part *restrict p, struct xpart *restrict xp, float dt) { /* Do not decrease the energy by more than a factor of 2*/ const float u_change = p->u_dt * dt; - if (u_change > -0.5f * p->u) - p->u += u_change; + if (u_change > -0.5f * xp->u_full) + xp->u_full += u_change; else - p->u *= 0.5f; - - /* Do not 'overcool' when timestep increases */ - if (p->u + p->u_dt * half_dt < 0.5f * p->u) p->u_dt = -0.5f * p->u / half_dt; + xp->u_full *= 0.5f; /* Compute the pressure */ - const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); + const float pressure = gas_pressure_from_internal_energy(p->rho, xp->u_full); /* Compute the sound speed */ const float soundspeed = gas_soundspeed_from_internal_energy(p->rho, p->u); @@ -442,9 +430,10 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( * This can be used to convert internal energy into entropy for instance. * * @param p The particle to act upon + * @param xp The extended particle to act upon */ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( - struct part *restrict p) { + struct part *restrict p, struct xpart *restrict xp) { /* Compute the pressure */ const float pressure = gas_pressure_from_internal_energy(p->rho, p->u); @@ -456,4 +445,27 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities( p->force.soundspeed = soundspeed; } +/** + * @brief Initialises the particles for the first time + * + * This function is called only once just after the ICs have been + * read in to do some conversions or assignments between the particle + * and extended particle fields. + * + * @param p The particle to act upon + * @param xp The extended particle data to act upon + */ +__attribute__((always_inline)) INLINE static void hydro_first_init_part( + struct part *restrict p, struct xpart *restrict xp) { + + p->time_bin = 0; + xp->v_full[0] = p->v[0]; + xp->v_full[1] = p->v[1]; + xp->v_full[2] = p->v[2]; + xp->u_full = p->u; + + hydro_reset_acceleration(p); + hydro_init_part(p); +} + #endif /* SWIFT_MINIMAL_HYDRO_H */ diff --git a/src/hydro/Minimal/hydro_debug.h b/src/hydro/Minimal/hydro_debug.h index 16ae62413a0d76b7bf871e615fe5684219752fee..876ce148824489d4c43358c2c519aa3b90dcf002 100644 --- a/src/hydro/Minimal/hydro_debug.h +++ b/src/hydro/Minimal/hydro_debug.h @@ -40,12 +40,11 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle( "v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e], " "u=%.3e, du/dt=%.3e v_sig=%.3e, P=%.3e\n" "h=%.3e, dh/dt=%.3e wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, " - "t_begin=%d, t_end=%d\n", + "time_bin=%d\n", p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0], xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2], p->u, p->u_dt, p->force.v_sig, p->force.pressure, p->h, p->force.h_dt, - (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho, p->ti_begin, - p->ti_end); + (int)p->density.wcount, p->mass, p->density.rho_dh, p->rho, p->time_bin); } #endif /* SWIFT_MINIMAL_HYDRO_DEBUG_H */ diff --git a/src/hydro/Minimal/hydro_part.h b/src/hydro/Minimal/hydro_part.h index 8542177278998d5e0b830dc164988611549ef24d..dabae1a546d66f61db4f9796c21b71817ca20aac 100644 --- a/src/hydro/Minimal/hydro_part.h +++ b/src/hydro/Minimal/hydro_part.h @@ -49,6 +49,9 @@ struct xpart { /*! Velocity at the last full step. */ float v_full[3]; + /*! Internal energy at the last full step. */ + float u_full; + /*! Additional data used to record cooling information */ struct cooling_xpart_data cooling_data; @@ -63,6 +66,12 @@ struct xpart { */ struct part { + /*! Particle unique ID. */ + long long id; + + /*! Pointer to corresponding gravity part. */ + struct gpart* gpart; + /*! Particle position. */ double x[3]; @@ -78,12 +87,6 @@ struct part { /*! Particle smoothing length. */ float h; - /*! Time at the beginning of time-step. */ - int ti_begin; - - /*! Time at the end of time-step. */ - int ti_end; - /*! Particle internal energy. */ float u; @@ -143,11 +146,18 @@ struct part { } force; }; - /*! Particle unique ID. */ - long long id; + /*! Time-step length */ + timebin_t time_bin; - /*! Pointer to corresponding gravity part. */ - struct gpart* gpart; +#ifdef SWIFT_DEBUG_CHECKS + + /* Time of the last drift */ + integertime_t ti_drift; + + /* Time of the last kick */ + integertime_t ti_kick; + +#endif } SWIFT_STRUCT_ALIGN;