diff --git a/examples/CoolingHaloWithSpin/cooling_halo.yml b/examples/CoolingHaloWithSpin/cooling_halo.yml index 2e5f09b9cae199090f29ae2341c232218f8099d0..1501c1e1dc00987dee7bdf25b008f393affe8cff 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 757c5e6600d85e19910d2fd75e955e5342fb63ee..ae13b49c9561096b11d021f710c687d153f2d269 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 0bfbf7ffef18f22cac9440ccc8adde0ba1899616..202099839227dffb0244dee09af0c87be14d6ab2 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);