Tried implementing primitive variable drift. Does not work (yet).

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;

src/hydro/Gizmo/hydro.h

@@ -19,6 +19,7 @@
 
 #include <float.h>
 #include "adiabatic_index.h"
+#include "approx_math.h"
 #include "hydro_gradients.h"
 
 /**
@@ -219,6 +220,9 @@ __attribute__((always_inline)) INLINE static void hydro_reset_acceleration(
   p->a_hydro[0] = 0.0f;
   p->a_hydro[1] = 0.0f;
   p->a_hydro[2] = 0.0f;
+
+  /* Reset the time derivatives. */
+  p->force.h_dt = 0.0f;
 }
 
 /**
@@ -273,10 +277,17 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
 
   return;
   float dt = (t1 - t0) * timeBase;
-
-  p->primitives.rho =
-      (p->conserved.mass + p->conserved.flux.mass * dt / p->force.dt) /
-      p->geometry.volume;
+  float h_inv = 1.0f / p->h;
+  float w = -hydro_dimension * p->force.h_dt * h_inv * dt;
+
+  // p->primitives.rho =
+  //    (p->conserved.mass + p->conserved.flux.mass * dt / p->force.dt) /
+  //    p->geometry.volume;
+  if (fabsf(w) < 0.2f) {
+    p->primitives.rho *= approx_expf(w);
+  } else {
+    p->primitives.rho *= expf(w);
+  }
   p->primitives.v[0] += p->a_hydro[0] * dt;
   p->primitives.v[1] += p->a_hydro[1] * dt;
   p->primitives.v[2] += p->a_hydro[2] * dt;
@@ -299,6 +310,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
 __attribute__((always_inline)) INLINE static void hydro_end_force(
     struct part* p) {
 
+  p->force.h_dt *= p->h * hydro_dimension_inv;
+
   /* Set the hydro acceleration, based on the new momentum and mass */
   /* NOTE: the momentum and mass are only correct for active particles, since
            only active particles have received flux contributions from all their

src/hydro/Gizmo/hydro_iact.h

@@ -265,6 +265,20 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
   xj = r * hj_inv;
   kernel_deval(xj, &wj, &wj_dx);
 
+  /* Compute h_dt */
+  // float dvdr = (Wi[1] - Wj[1]) * dx[0] + (Wi[2] - Wj[2]) * dx[1] +
+  // (Wi[3]-Wj[3])*dx[2];
+  // float ri = 1.0f/r;
+  // float hidp1 = pow_dimension_plus_one(hi_inv);
+  // float hjdp1 = pow_dimension_plus_one(hj_inv);
+  // float wi_dr = hidp1 * wi_dx;
+  // float wj_dr = hjdp1 * wj_dx;
+  // dvdr *= ri;
+  // pi->force.h_dt -= pj->conserved.mass * dvdr / pj->primitives.rho * wi_dr;
+  // if(mode == 1){
+  //  pj->force.h_dt -= pi->conserved.mass * dvdr / pi->primitives.rho * wj_dr;
+  //}
+
   /* Compute area */
   /* eqn. (7) */
   Anorm = 0.0f;

src/hydro/Gizmo/hydro_part.h

@@ -167,7 +167,7 @@ struct part {
   /* Quantities used during the force loop. */
   struct {
 
-    /* Needed to compile the code. */
+    /* Needed to drift the primitive variables. */
     float h_dt;
 
     /* Physical time step of the particle. */