diff --git a/src/hydro/GizmoMFV/hydro.h b/src/hydro/GizmoMFV/hydro.h index 6916fe33272692316354385b723ce9969606b6a2..281c66337cf761cd98bc82238460eca859bc6e88 100644 --- a/src/hydro/GizmoMFV/hydro.h +++ b/src/hydro/GizmoMFV/hydro.h @@ -591,12 +591,28 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( p->conserved.flux.momentum[2] * dt_drift / p->conserved.mass; #if !defined(EOS_ISOTHERMAL_GAS) - const float u = p->conserved.energy + p->conserved.flux.energy * dt_therm; - p->primitives.P = - hydro_gamma_minus_one * u * p->primitives.rho / p->conserved.mass; +#ifdef GIZMO_TOTAL_ENERGY + const float Etot = + p->conserved.energy + p->conserved.flux.energy * dt_therm; + const float v2 = (p->primitives.v[0] * p->primitives.v[0] + + p->primitives.v[1] * p->primitives.v[1] + + p->primitives.v[2] * p->primitives.v[2]); + const float u = (Etot / p->conserved.mass - 0.5 * v2); +#else + const float u = + (p->conserved.energy + p->conserved.flux.energy * dt_therm) / + p->conserved.mass; +#endif + p->primitives.P = hydro_gamma_minus_one * u * p->primitives.rho; #endif } + /* we use a sneaky way to get the gravitational contribtuion to the + velocity update */ + p->primitives.v[0] += p->v[0] - xp->v_full[0]; + p->primitives.v[1] += p->v[1] - xp->v_full[1]; + p->primitives.v[2] += p->v[2] - xp->v_full[2]; + #ifdef SWIFT_DEBUG_CHECKS if (p->h <= 0.) { error("Zero or negative smoothing length (%g)!", p->h); @@ -646,6 +662,33 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( float a_grav[3]; + /* Add gravity. We only do this if we have gravity activated. */ + if (p->gpart) { + /* Retrieve the current value of the gravitational acceleration from the + gpart. We are only allowed to do this because this is the kick. We still + need to check whether gpart exists though.*/ + a_grav[0] = p->gpart->a_grav[0]; + a_grav[1] = p->gpart->a_grav[1]; + a_grav[2] = p->gpart->a_grav[2]; + +#ifdef GIZMO_TOTAL_ENERGY + p->conserved.energy += dt * (p->conserved.momentum[0] * a_grav[0] + + p->conserved.momentum[1] * a_grav[1] + + p->conserved.momentum[2] * a_grav[2]); +#endif + + /* Kick the momentum for half a time step */ + /* Note that this also affects the particle movement, as the velocity for + the particles is set after this. */ + p->conserved.momentum[0] += p->conserved.mass * a_grav[0] * dt; + p->conserved.momentum[1] += p->conserved.mass * a_grav[1] * dt; + p->conserved.momentum[2] += p->conserved.mass * a_grav[2] * dt; + + p->conserved.energy -= 0.5f * dt * (p->gravity.mflux[0] * a_grav[0] + + p->gravity.mflux[1] * a_grav[1] + + p->gravity.mflux[2] * a_grav[2]); + } + /* Update conserved variables. */ p->conserved.mass += p->conserved.flux.mass * dt; p->conserved.momentum[0] += p->conserved.flux.momentum[0] * dt; @@ -688,28 +731,9 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( } #endif - /* Add gravity. We only do this if we have gravity activated. */ if (p->gpart) { - /* Retrieve the current value of the gravitational acceleration from the - gpart. We are only allowed to do this because this is the kick. We still - need to check whether gpart exists though.*/ - a_grav[0] = p->gpart->a_grav[0]; - a_grav[1] = p->gpart->a_grav[1]; - a_grav[2] = p->gpart->a_grav[2]; - /* Make sure the gpart knows the mass has changed. */ p->gpart->mass = p->conserved.mass; - - /* Kick the momentum for half a time step */ - /* Note that this also affects the particle movement, as the velocity for - the particles is set after this. */ - p->conserved.momentum[0] += dt * p->conserved.mass * a_grav[0]; - p->conserved.momentum[1] += dt * p->conserved.mass * a_grav[1]; - p->conserved.momentum[2] += dt * p->conserved.mass * a_grav[2]; - - p->conserved.energy += dt * (p->gravity.mflux[0] * a_grav[0] + - p->gravity.mflux[1] * a_grav[1] + - p->gravity.mflux[2] * a_grav[2]); } hydro_velocities_set(p, xp); diff --git a/src/hydro/GizmoMFV/hydro_velocities.h b/src/hydro/GizmoMFV/hydro_velocities.h index ea30d5a6c9c74d34ffd73fa6ab941640f37b02e4..a61a482683291cfc0662116ada99b94a4ccfe68d 100644 --- a/src/hydro/GizmoMFV/hydro_velocities.h +++ b/src/hydro/GizmoMFV/hydro_velocities.h @@ -103,7 +103,6 @@ __attribute__((always_inline)) INLINE static void hydro_velocities_set( p->v[2] = p->conserved.momentum[2] * inverse_mass; #ifdef GIZMO_STEER_MOTION - /* Add a correction to the velocity to keep particle positions close enough to the centroid of their mesh-free "cell". */