From 7894d3b5cf73ed6c74e7333b9d0f130923f7eb61 Mon Sep 17 00:00:00 2001
From: Stefan Arridge <stefan.arridge@durham.ac.uk>
Date: Thu, 26 Jan 2017 18:27:45 +0000
Subject: [PATCH] Fixed enforcement of energy floor. entropy_dt is limited in
 hydro_kick_extra

---
 examples/CoolingHaloWithSpin/cooling_halo.yml | 2 +-
 src/cooling/const_lambda/cooling.h            | 6 ++----
 src/hydro/Gadget2/hydro.h                     | 7 ++++---
 3 files changed, 7 insertions(+), 8 deletions(-)

diff --git a/examples/CoolingHaloWithSpin/cooling_halo.yml b/examples/CoolingHaloWithSpin/cooling_halo.yml
index 2e5f09b9ca..1501c1e1dc 100644
--- a/examples/CoolingHaloWithSpin/cooling_halo.yml
+++ b/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
diff --git a/src/cooling/const_lambda/cooling.h b/src/cooling/const_lambda/cooling.h
index 757c5e6600..ae13b49c95 100644
--- a/src/cooling/const_lambda/cooling.h
+++ b/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 */
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 0bfbf7ffef..2020998392 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/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);
-- 
GitLab