No spurious factor of two in the radiateed energy statistics calculation

examples/CoolingBox/energy_plot.py

@@ -94,7 +94,7 @@ initial_energy_cgs = initial_energy/total_mass[0] * unit_length**2 / (unit_time)
 n_H_cgs = X_H * rho_cgs / m_p
 du_dt_cgs = -cooling_lambda * n_H_cgs**2 / rho_cgs
 cooling_time_cgs = (initial_energy_cgs/(-du_dt_cgs))[0]
-analytic_time_cgs = np.linspace(0, cooling_time_cgs*1.5,1000)
+analytic_time_cgs = np.linspace(0, cooling_time_cgs * 1.8, 1000)
 u_analytic_cgs = du_dt_cgs*analytic_time_cgs + initial_energy_cgs
 u_analytic_cgs[u_analytic_cgs < u_floor_cgs] = u_floor_cgs
 
@@ -103,11 +103,11 @@ print "Cooling time:", cooling_time_cgs, "[s]"
 figure()
 #plot(time_cgs, internal_energy_cgs, 'b-', lw=1.5, label="Gas internal energy")
 #plot(time_cgs, kinetic_energy_cgs, 'y-', lw=1.5, label="Gas kinetic energy")
-plot(time_cgs, total_energy_cgs, 'r-', lw=1.5, label="Gas total energy")
-plot(time_cgs, radiated_energy_cgs, 'g-', lw=1.5, label="Radiated energy")
-plot(time_cgs, total_energy_cgs + radiated_energy_cgs, 'k-', lw=0.6, label="Gas+radiated")
+plot(time_cgs, total_energy_cgs, 'r-', lw=1.6, label="Gas total energy")
+plot(time_cgs, radiated_energy_cgs, 'g-', lw=1.6, label="Radiated energy")
+plot(time_cgs, total_energy_cgs + radiated_energy_cgs, 'b-', lw=0.6, label="Gas total + radiated")
 
-plot(analytic_time_cgs, u_analytic_cgs, '--', color='k', alpha=0.8, lw=1.2, label="Analytic solution")
+plot(analytic_time_cgs, u_analytic_cgs, '--', color='k', alpha=0.8, lw=1.0, label="Analytic solution")
 
 legend(loc="upper right", fontsize=8)
 xlabel("${\\rm{Time~[s]}}$", labelpad=0)

src/cooling/const_lambda/cooling.h

@@ -89,8 +89,10 @@ __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 + cooling_du_dt * dt < u_floor) {
-    cooling_du_dt = (u_old - u_floor) / dt;
+  if (u_old + hydro_du_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;
   }
 
   /* Update the internal energy time derivative */

src/hydro/Gadget2/hydro.h

@@ -130,9 +130,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy_dt(
 __attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
     struct part *restrict p, float du_dt) {
 
-  const float ds_dt = gas_entropy_from_internal_energy(p->rho, du_dt);
-
-  p->entropy_dt = ds_dt;
+  p->entropy_dt = gas_entropy_from_internal_energy(p->rho, du_dt);
 }
 
 /**
@@ -379,12 +377,10 @@ __attribute__((always_inline)) INLINE static void hydro_end_force(
 __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*/
+  /* Do not decrease the entropy 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;
+  xp->entropy_full =
+      max(xp->entropy_full + entropy_change, 0.5f * xp->entropy_full);
 
   /* Compute the pressure */
   const float pressure = gas_pressure_from_entropy(p->rho, xp->entropy_full);

src/statistics.c

@@ -163,7 +163,7 @@ void stats_collect_part_mapper(void *map_data, int nr_parts, void *extra_data) {
     /* Collect energies. */
     stats.E_kin += 0.5f * m * (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]);
     stats.E_int += m * hydro_get_internal_energy(p);
-    stats.E_rad += cooling_get_radiated_energy(xp) / 2.;
+    stats.E_rad += cooling_get_radiated_energy(xp);
     stats.E_pot_self += 0.f;
     if (gp != NULL)
       stats.E_pot_ext += m * external_gravity_get_potential_energy(