diff --git a/src/engine.c b/src/engine.c
index 4a91c6be94125999a23db50aa5873830974aaa0e..dcdd956250045d6274ebb1ee2cbd5c43be216e83 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 60f7cf55f35213ba3d570c355412f6587db7a219..1e3aef0652e25a062f517deaed6b96c6831b8d89 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -318,15 +318,14 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
            neighbours. Since this method is only called for active particles,
            this is indeed the case. */
   if (p->force.dt) {
-    p->a_hydro[0] =
-        (p->conserved.momentum[0] / p->conserved.mass - p->primitives.v[0]) /
-        p->force.dt;
-    p->a_hydro[1] =
-        (p->conserved.momentum[1] / p->conserved.mass - p->primitives.v[1]) /
-        p->force.dt;
-    p->a_hydro[2] =
-        (p->conserved.momentum[2] / p->conserved.mass - p->primitives.v[2]) /
-        p->force.dt;
+    float mnew = p->conserved.mass + p->conserved.flux.mass;
+    float pnew[3];
+    pnew[0] = p->conserved.momentum[0] + p->conserved.flux.momentum[0];
+    pnew[1] = p->conserved.momentum[1] + p->conserved.flux.momentum[1];
+    pnew[2] = p->conserved.momentum[2] + p->conserved.flux.momentum[2];
+    p->a_hydro[0] = (pnew[0] / mnew - p->primitives.v[0]) / p->force.dt;
+    p->a_hydro[1] = (pnew[1] / mnew - p->primitives.v[1]) / p->force.dt;
+    p->a_hydro[2] = (pnew[2] / mnew - p->primitives.v[2]) / p->force.dt;
   } else {
     p->a_hydro[0] = 0.0f;
     p->a_hydro[1] = 0.0f;