diff --git a/src/engine.c b/src/engine.c index ac4ac50d49f244646c762e8efa6534cf9cb407f5..7d5027237222ad3bfbfb832d7040a180158e8f4d 100644 --- a/src/engine.c +++ b/src/engine.c @@ -2720,7 +2720,7 @@ void engine_step(struct engine *e) { mask |= 1 << task_type_sub_pair; mask |= 1 << task_type_ghost; mask |= 1 << task_type_extra_ghost; /* Adding unnecessary things to the mask - does not harm */ + does not harm */ submask |= 1 << task_subtype_density; submask |= 1 << task_subtype_gradient; diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h index a4d9e58087879645d5c98912ec2039015bcc145a..b90632b986b16816eb4005963ee0bdf300340b6f 100644 --- a/src/hydro/Gizmo/hydro.h +++ b/src/hydro/Gizmo/hydro.h @@ -202,13 +202,6 @@ __attribute__((always_inline)) INLINE static void hydro_end_gradient( struct part* p) { hydro_gradients_finalize(p); - - /* Initialize fluxes */ - p->conserved.flux.mass = 0.0f; - p->conserved.flux.momentum[0] = 0.0f; - p->conserved.flux.momentum[1] = 0.0f; - p->conserved.flux.momentum[2] = 0.0f; - p->conserved.flux.energy = 0.0f; } /** @@ -339,7 +332,23 @@ __attribute__((always_inline)) INLINE static void hydro_end_force( * @param half_dt Half the physical time step. */ __attribute__((always_inline)) INLINE static void hydro_kick_extra( - struct part* p, struct xpart* xp, float dt, float half_dt) {} + struct part* p, struct xpart* xp, float dt, float half_dt) { + + p->conserved.mass += p->conserved.flux.mass; + p->conserved.momentum[0] += p->conserved.flux.momentum[0]; + p->conserved.momentum[1] += p->conserved.flux.momentum[1]; + p->conserved.momentum[2] += p->conserved.flux.momentum[2]; + p->conserved.energy += p->conserved.flux.energy; + + /* reset fluxes */ + /* we can only do this here, since we need to keep the fluxes for inactive + particles */ + p->conserved.flux.mass = 0.0f; + p->conserved.flux.momentum[0] = 0.0f; + p->conserved.flux.momentum[1] = 0.0f; + p->conserved.flux.momentum[2] = 0.0f; + p->conserved.flux.energy = 0.0f; +} /** * @brief Returns the internal energy of a particle diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h index d4e7bcb6bdefebdf279ddbdba526b238a7c6cfa3..f19082547d3c22fff689d44c2ae7805194611dd3 100644 --- a/src/hydro/Gizmo/hydro_iact.h +++ b/src/hydro/Gizmo/hydro_iact.h @@ -250,8 +250,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common( /* The flux will be exchanged using the smallest time step of the two * particles */ mindt = fminf(dti, dtj); - // dti = mindt; - // dtj = mindt; + dti = mindt; + dtj = mindt; /* Compute kernel of pi. */ hi_inv = 1.0 / hi; @@ -365,12 +365,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common( /* Update conserved variables */ /* eqn. (16) */ - pi->conserved.mass -= dti * Anorm * totflux[0]; - pi->conserved.momentum[0] -= dti * Anorm * totflux[1]; - pi->conserved.momentum[1] -= dti * Anorm * totflux[2]; - pi->conserved.momentum[2] -= dti * Anorm * totflux[3]; - pi->conserved.energy -= dti * Anorm * totflux[4]; - pi->conserved.flux.mass -= dti * Anorm * totflux[0]; pi->conserved.flux.momentum[0] -= dti * Anorm * totflux[1]; pi->conserved.flux.momentum[1] -= dti * Anorm * totflux[2]; @@ -380,11 +374,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common( float ekin = 0.5f * (pi->primitives.v[0] * pi->primitives.v[0] + pi->primitives.v[1] * pi->primitives.v[1] + pi->primitives.v[2] * pi->primitives.v[2]); - pi->conserved.energy += dti * Anorm * totflux[1] * pi->primitives.v[0]; - pi->conserved.energy += dti * Anorm * totflux[2] * pi->primitives.v[1]; - pi->conserved.energy += dti * Anorm * totflux[3] * pi->primitives.v[2]; - pi->conserved.energy -= dti * Anorm * totflux[0] * ekin; - pi->conserved.flux.energy += dti * Anorm * totflux[1] * pi->primitives.v[0]; pi->conserved.flux.energy += dti * Anorm * totflux[2] * pi->primitives.v[1]; pi->conserved.flux.energy += dti * Anorm * totflux[3] * pi->primitives.v[2]; @@ -402,14 +391,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common( UPDATE particle j. ==> we update particle j if (MODE IS 1) OR (j IS INACTIVE) */ - // if (mode == 1 || pj->ti_end > pi->ti_end) { - if (mode == 1) { - pj->conserved.mass += dtj * Anorm * totflux[0]; - pj->conserved.momentum[0] += dtj * Anorm * totflux[1]; - pj->conserved.momentum[1] += dtj * Anorm * totflux[2]; - pj->conserved.momentum[2] += dtj * Anorm * totflux[3]; - pj->conserved.energy += dtj * Anorm * totflux[4]; - + if (mode == 1 || pj->ti_end > pi->ti_end) { pj->conserved.flux.mass += dtj * Anorm * totflux[0]; pj->conserved.flux.momentum[0] += dtj * Anorm * totflux[1]; pj->conserved.flux.momentum[1] += dtj * Anorm * totflux[2]; @@ -419,11 +401,6 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common( ekin = 0.5f * (pj->primitives.v[0] * pj->primitives.v[0] + pj->primitives.v[1] * pj->primitives.v[1] + pj->primitives.v[2] * pj->primitives.v[2]); - pj->conserved.energy -= dtj * Anorm * totflux[1] * pj->primitives.v[0]; - pj->conserved.energy -= dtj * Anorm * totflux[2] * pj->primitives.v[1]; - pj->conserved.energy -= dtj * Anorm * totflux[3] * pj->primitives.v[2]; - pj->conserved.energy += dtj * Anorm * totflux[0] * ekin; - pj->conserved.flux.energy -= dtj * Anorm * totflux[1] * pj->primitives.v[0]; pj->conserved.flux.energy -= dtj * Anorm * totflux[2] * pj->primitives.v[1]; pj->conserved.flux.energy -= dtj * Anorm * totflux[3] * pj->primitives.v[2];