Extended the modifications to the other schemes.

examples/CoolingHaloWithSpin/cooling_halo.yml

@@ -9,8 +9,8 @@ InternalUnitSystem:
 # Parameters governing the time integration
 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-8 # The minimal time-step size of the simulation (in internal units).
+  time_end:   10.   # The end time of the simulation (in internal units).
+  dt_min:     1e-5  # 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
@@ -40,7 +40,7 @@ IsothermalPotential:
   position_z:      0.	
   vrot:            200.   # Rotation speed of isothermal potential in internal units
   timestep_mult:   0.03   # Controls time step
-  epsilon:         0.1    # Softening for the isothermal potential
+  epsilon:         1.0    # Softening for the isothermal potential
 
 # Cooling parameters
 LambdaCooling:

src/hydro/Minimal/hydro.h

@@ -360,8 +360,10 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
     struct part *restrict p, struct xpart *restrict xp, float dt) {
 
   /* Do not decrease the energy by more than a factor of 2*/
-  const float u_change = p->u_dt * dt;
-  xp->u_full = max(xp->u_full + u_change, 0.5f * xp->u_full);
+  if (p->u_dt < -0.5f * xp->u_full / dt) {
+    p->u_dt = -0.5f * xp->u_full / dt;
+  }
+  xp->u_full += p->u_dt * dt;
 
   /* Compute the pressure */
   const float pressure = gas_pressure_from_internal_energy(p->rho, xp->u_full);

src/hydro/PressureEntropy/hydro.h

@@ -394,11 +394,10 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
     struct part *restrict p, struct xpart *restrict xp, float dt) {
 
   /* Do not decrease the entropy (temperature) by more than a factor of 2*/
-  const float entropy_change = p->entropy_dt * dt;
-  if (entropy_change > -0.5f * xp->entropy_full)
-    xp->entropy_full += entropy_change;
-  else
-    xp->entropy_full *= 0.5f;
+  if (p->entropy_dt < -0.5f * xp->entropy_full / dt) {
+    p->entropy_dt = -0.5f * xp->entropy_full / dt;
+  }
+  xp->entropy_full += p->entropy_dt * dt;
 
   /* Compute the pressure */
   const float pressure =