diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h index 461d2f99ad76eb566f33ad84f0eef629f40b6cc6..5d1f5aca8af5e906fa29d86b45b1105df6b6e411 100644 --- a/src/hydro/Gadget2/hydro_iact.h +++ b/src/hydro/Gadget2/hydro_iact.h @@ -158,8 +158,17 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density( const float ih3 = h_inv * h_inv * h_inv; const float ih4 = h_inv * h_inv * h_inv * h_inv; + const float fac = mj * wi_dx * ri; + + /* Compute dv dot r */ + dv[0] = pi->v[0] - pj->v[0]; + dv[1] = pi->v[1] - pj->v[1]; + dv[2] = pi->v[2] - pj->v[2]; + const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2]; + pi->div_v -= fac * dvdr; + if(pi->id == 1000 && pj->id == 1103) - message("Interacting with %lld. r=%e hi=%e u=%e W=%e dW/dx=%e dh_drho1=%e dh_drho2=%e %e\n", + message("Interacting with %lld. r=%e hi=%e u=%e W=%e dW/dx=%e dh_drho1=%e dh_drho2=%e %e\n fac=%e dvdr=%e", pj->id, r, hi, @@ -168,18 +177,11 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_density( wi_dx * ih4, -mj * (3.f * kernel_igamma * wi) * ih4, -mj * u * wi_dx * kernel_igamma * ih4, - kernel_igamma + fac, + dvdr ); - const float fac = mj * wi_dx * ri; - /* Compute dv dot r */ - dv[0] = pi->v[0] - pj->v[0]; - dv[1] = pi->v[1] - pj->v[1]; - dv[2] = pi->v[2] - pj->v[2]; - const float dvdr = dv[0] * dx[0] + dv[1] * dx[1] + dv[2] * dx[2]; - pi->div_v -= fac * dvdr; - /* Compute dv cross r */ curlvr[0] = dv[1] * dx[2] - dv[2] * dx[1]; curlvr[1] = dv[2] * dx[0] - dv[0] * dx[2];