Merged Stefan's corrections to the const-lambda cooling time-step criterion.

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-4  # The minimal time-step size of the simulation (in internal units).
+  dt_min:     1e-7  # 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
@@ -32,9 +32,6 @@ SPH:
 # Parameters related to the initial conditions
 InitialConditions:
   file_name:  CoolingHalo.hdf5       # The file to read
-  shift_x:    0.                  # A shift to apply to all particles read from the ICs (in internal units).
-  shift_y:    0.
-  shift_z:    0.
  
 # External potential parameters
 SoftenedIsothermalPotential:
@@ -43,12 +40,12 @@ SoftenedIsothermalPotential:
   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:
-  lambda_cgs:                  1.0e-22    # Cooling rate (in cgs units)
+  lambda_cgs:                  1.0e-22   # Cooling rate (in cgs units)
   minimum_temperature:         1.0e4  # Minimal temperature (Kelvin)
   mean_molecular_weight:       0.59   # Mean molecular weight
   hydrogen_mass_abundance:     0.75   # Hydrogen mass abundance (dimensionless)
-  cooling_tstep_mult:          1.0    # Dimensionless pre-factor for the time-step condition
+  cooling_tstep_mult:          0.1    # Dimensionless pre-factor for the time-step condition

src/cooling/const_lambda/cooling.h

@@ -24,6 +24,7 @@
 #define SWIFT_COOLING_CONST_LAMBDA_H
 
 /* Some standard headers. */
+#include <float.h>
 #include <math.h>
 
 /* Local includes. */
@@ -84,7 +85,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
   /* Calculate du_dt */
   const float du_dt = cooling_rate(phys_const, us, cooling, p);
 
-  /* Intergrate cooling equation, but enforce energy floor */
+  /* Integrate cooling equation, but enforce energy floor */
   float u_new;
   if (u_old + du_dt * dt > u_floor) {
     u_new = u_old + du_dt * dt;
@@ -92,6 +93,9 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
     u_new = u_floor;
   }
 
+  /* Don't allow particle to cool too much in one timestep */
+  if (u_new < 0.5f * u_old) u_new = 0.5f * u_old;
+
   /* Update the internal energy */
   hydro_set_internal_energy(p, u_new);
 
@@ -112,13 +116,16 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
     const struct phys_const* restrict phys_const,
     const struct UnitSystem* restrict us, const struct part* restrict p) {
 
-  /* Get du_dt */
-  const float du_dt = cooling_rate(phys_const, us, cooling, p);
-
   /* Get current internal energy (dt=0) */
   const float u = hydro_get_internal_energy(p, 0.f);
+  const float du_dt = cooling_rate(phys_const, us, cooling, p);
 
-  return u / fabsf(du_dt);
+  /* If we are close to (or below) the energy floor, we ignore cooling timestep
+   */
+  if (u < 1.01f * cooling->min_energy)
+    return FLT_MAX;
+  else
+    return cooling->cooling_tstep_mult * u / fabsf(du_dt);
 }
 
 /**