Commit 5b76e0a9 authored by Bert Vandenbroucke's avatar Bert Vandenbroucke
Browse files

Moved around the energy gravity update so that it uses the old momentum....

Moved around the energy gravity update so that it uses the old momentum. Removed a superfluous time step in mflux.
parent 30482afb
......@@ -646,15 +646,9 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
/* 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];
#if !defined(EOS_ISOTHERMAL_GAS)
/* This part still needs to be tested! */
/* If the energy needs to be updated, we need to do it before the momentum
is updated, as the old value of the momentum enters the equations. */
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]);
......@@ -663,6 +657,14 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
a_grav[1] * p->gravity.mflux[1] +
a_grav[2] * p->gravity.mflux[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] += 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];
}
/* reset fluxes */
......
......@@ -346,8 +346,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
}
dvdotdx = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
(Wi[3] - Wj[3]) * dx[2];
dvdotdx = min(dvdotdx,
(vi[0] - vj[0])*dx[0] + (vi[1]-vj[1])*dx[1] + (vi[2]-vj[2])*dx[2]);
if (dvdotdx < 0.) {
vmax -= dvdotdx / r;
/* the magical factor 3 also appears in Gadget2 */
vmax -= 3.*dvdotdx / r;
}
dt_min = r / vmax;
pi->timestepvars.dt_min = min(pi->timestepvars.dt_min, dt_min);
......@@ -517,7 +520,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
totflux[4] *= flux_limit_factor;
/* Store mass flux */
float mflux = mindt * Anorm * totflux[0];
float mflux = Anorm * totflux[0];
pi->gravity.mflux[0] += mflux * dx[0];
pi->gravity.mflux[1] += mflux * dx[1];
pi->gravity.mflux[2] += mflux * dx[2];
......@@ -555,7 +558,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
if (mode == 1 || pj->force.active == 0) {
/* Store mass flux */
mflux = mindt * Anorm * totflux[0];
mflux = Anorm * totflux[0];
pj->gravity.mflux[0] -= mflux * dx[0];
pj->gravity.mflux[1] -= mflux * dx[1];
pj->gravity.mflux[2] -= mflux * dx[2];
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment