Commit 35fa1dc3 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Smoothing length time derivative computed correctly in the Gadget-2 case

parent 29bf9e61
......@@ -203,7 +203,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
*/
__attribute__((always_inline))
INLINE static void hydro_end_force(struct part* p) {
p->entropy_dt *=
(const_hydro_gamma - 1.f) * powf(p->rho, -(const_hydro_gamma - 1.f));
}
......
......@@ -197,7 +197,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float r_inv = 1.0f / r;
/* Get some values in local variables. */
// const float mi = pi->mass;
const float mi = pi->mass;
const float mj = pj->mass;
const float rhoi = pi->rho;
const float rhoj = pj->rho;
......@@ -280,6 +280,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj->a[1] += acc * dx[1];
pj->a[2] += acc * dx[2];
/* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdr / rhoj * wi_dr;
pj->force.h_dt -= mi * dvdr / rhoi * wj_dr;
/* Update the signal velocity. */
pi->v_sig = fmaxf(pi->v_sig, v_sig);
pj->v_sig = fmaxf(pj->v_sig, v_sig);
......@@ -366,6 +370,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->a[1] -= acc * dx[1];
pi->a[2] -= acc * dx[2];
/* Get the time derivative for h. */
pi->force.h_dt -= mj * dvdr / rhoj * wi_dr;
/* Update the signal velocity. */
pi->v_sig = fmaxf(pi->v_sig, v_sig);
......
......@@ -851,8 +851,13 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
if (is_fixdt || p->t_end <= t_current) {
/* First, finish the force loop */
p->force.h_dt *= p->h * 0.333333333f;
/* And do the same of the extra variable */
hydro_end_force(p);
/* Now we are ready to compute the next time-step size */
if (is_fixdt) {
/* Now we have a time step, proceed with the kick */
......
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