Fixed enforcement of energy floor. entropy_dt is limited in hydro_kick_extra

examples/CoolingHaloWithSpin/cooling_halo.yml

@@ -10,7 +10,7 @@ InternalUnitSystem:
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
   time_end:   10.    # The end time of the simulation (in internal units).
-  dt_min:     1e-10 # The minimal time-step size of the simulation (in internal units).
+  dt_min:     1e-8 # The minimal time-step size of the simulation (in internal units).
   dt_max:     1e-1  # The maximal time-step size of the simulation (in internal units).
 
 # Parameters governing the conserved quantities statistics

src/cooling/const_lambda/cooling.h

@@ -89,10 +89,8 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
   float cooling_du_dt = cooling_rate(phys_const, us, cooling, p);
 
   /* Integrate cooling equation to enforce energy floor */
-  if (u_old + hydro_du_dt * dt < u_floor) {
-    cooling_du_dt = 0.f;
-  } else if (u_old + (hydro_du_dt + cooling_du_dt) * dt < u_floor) {
-    cooling_du_dt = (u_old + dt * hydro_du_dt - u_floor) / dt;
+  if (u_old + (hydro_du_dt + cooling_du_dt) * dt < u_floor) {
+    cooling_du_dt = -(u_old + dt * hydro_du_dt - u_floor) / dt;
   }
 
   /* Update the internal energy time derivative */

src/hydro/Gadget2/hydro.h

@@ -332,6 +332,8 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
   else
     p->rho *= expf(w2);
 
+  if(p->entropy + p->entropy_dt*dt < 0.0)
+    printf("Negative entropy: Old entropy = %g, New entropy = %g ", p->entropy,p->entropy + p->entropy_dt*dt);
   /* Predict the entropy */
   p->entropy += p->entropy_dt * dt;
 
@@ -377,9 +379,8 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
     struct part *restrict p, struct xpart *restrict xp, float dt) {
 
   /* Do not decrease the entropy by more than a factor of 2 */
-  const float entropy_change = p->entropy_dt * dt;
-  xp->entropy_full =
-      max(xp->entropy_full + entropy_change, 0.5f * xp->entropy_full);
+  p->entropy_dt = max(-0.5f * xp->entropy_full / dt , p->entropy_dt); 
+  xp->entropy_full += p->entropy_dt * dt;
 
   /* Compute the pressure */
   const float pressure = gas_pressure_from_entropy(p->rho, xp->entropy_full);