Updated the cooling to the new integration scheme

examples/CoolingBox/coolingBox.yml

@@ -10,8 +10,8 @@ InternalUnitSystem:
 TimeIntegration:
   time_begin: 0.    # The starting time of the simulation (in internal units).
   time_end:   4.    # 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_max:     1e-4  # The maximal time-step size of the simulation (in internal units).
+  dt_min:     1e-5  # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e-2  # The maximal time-step size of the simulation (in internal units).
 
 # Parameters governing the snapshots
 Snapshots:
@@ -35,7 +35,7 @@ InitialConditions:
 
 # Dimensionless pre-factor for the time-step condition
 LambdaCooling:
-  lambda_cgs:                 1.0e-22    # Cooling rate (in cgs units)
+  lambda_cgs:                 0. #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)

examples/CoolingBox/run.sh

@@ -5,7 +5,7 @@ echo "Generating initial conditions for the cooling box example..."
 
 python makeIC.py 10
 
-../swift -s -C -t 16 coolingBox.yml 
+../swift -s -C -t 1 coolingBox.yml
 
 #-C 2>&1 | tee output.log
 

src/cooling/const_du/cooling.h

@@ -63,7 +63,7 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
     struct part* restrict p, struct xpart* restrict xp, float dt) {
 
   /* Get current internal energy (dt=0) */
-  const float u_old = hydro_get_internal_energy(p, 0.f);
+  const float u_old = hydro_get_internal_energy(p);
 
   /* Get cooling function properties */
   const float du_dt = -cooling->cooling_rate;
@@ -102,7 +102,7 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
     const struct UnitSystem* restrict us, const struct part* restrict p) {
 
   const float cooling_rate = cooling->cooling_rate;
-  const float internal_energy = hydro_get_internal_energy(p, 0);
+  const float internal_energy = hydro_get_internal_energy(p);
   return cooling->cooling_tstep_mult * internal_energy / fabsf(cooling_rate);
 }
 

src/cooling/const_lambda/cooling.h

@@ -76,31 +76,28 @@ __attribute__((always_inline)) INLINE static void cooling_cool_part(
     const struct cooling_function_data* restrict cooling,
     struct part* restrict p, struct xpart* restrict xp, float dt) {
 
-  /* Get current internal energy (dt=0) */
-  const float u_old = hydro_get_internal_energy(p, 0.f);
-
   /* Internal energy floor */
   const float u_floor = cooling->min_energy;
 
-  /* Calculate du_dt */
-  const float du_dt = cooling_rate(phys_const, us, cooling, p);
+  /* Current energy */
+  const float u_old = hydro_get_internal_energy(p);
 
-  /* 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;
-  } else {
-    u_new = u_floor;
-  }
+  /* Current du_dt */
+  const float hydro_du_dt = hydro_get_internal_energy_dt(p);
+
+  /* Calculate cooling du_dt */
+  float cooling_du_dt = cooling_rate(phys_const, us, cooling, p);
 
-  /* Don't allow particle to cool too much in one timestep */
-  if (u_new < 0.5f * u_old) u_new = 0.5f * u_old;
+  /* Integrate cooling equation to enforce energy floor */
+  if (u_old + cooling_du_dt * dt < u_floor) {
+    cooling_du_dt = (u_old - u_floor) / dt;
+  }
 
-  /* Update the internal energy */
-  hydro_set_internal_energy(p, u_new);
+  /* Update the internal energy time derivative */
+  hydro_set_internal_energy_dt(p, hydro_du_dt + cooling_du_dt);
 
   /* Store the radiated energy */
-  xp->cooling_data.radiated_energy += hydro_get_mass(p) * (u_old - u_new);
+  xp->cooling_data.radiated_energy += hydro_get_mass(p) * cooling_du_dt * dt;
 }
 
 /**
@@ -117,11 +114,10 @@ __attribute__((always_inline)) INLINE static float cooling_timestep(
     const struct UnitSystem* restrict us, const struct part* restrict p) {
 
   /* Get current internal energy (dt=0) */
-  const float u = hydro_get_internal_energy(p, 0.f);
+  const float u = hydro_get_internal_energy(p);
   const float du_dt = cooling_rate(phys_const, us, cooling, p);
 
-  /* If we are close to (or below) the energy floor, we ignore cooling timestep
-   */
+  /* If we are close to (or below) the energy floor, we ignore the condition */
   if (u < 1.01f * cooling->min_energy)
     return FLT_MAX;
   else

src/engine.c

@@ -1833,20 +1833,11 @@ void engine_make_extra_hydroloop_tasks(struct engine *e) {
 #endif
     }
 
-    /* /\* Cooling tasks should depend on kick and unlock sourceterms *\/ */
-    /* else if (t->type == task_type_cooling) { */
-    /*   scheduler_addunlock(sched, t->ci->kick, t); */
-    /* } */
-    /* /\* source terms depend on cooling if performed, else on kick. It is the
-     * last */
-    /*    task *\/ */
-    /* else if (t->type == task_type_sourceterms) { */
-    /*   if (e->policy == engine_policy_cooling) */
-    /*     scheduler_addunlock(sched, t->ci->cooling, t); */
-    /*   else */
-    /*     scheduler_addunlock(sched, t->ci->kick, t); */
-    /* } */
-    // MATTHIEU
+    /* Cooling tasks take place after the second kick */
+    else if (t->type == task_type_cooling) {
+      scheduler_addunlock(sched, t->ci->kick2, t);
+      scheduler_addunlock(sched, t, t->ci->timestep);
+    }
   }
 }
 

src/hydro/Gadget2/hydro.h

@@ -43,40 +43,33 @@
  * @brief Returns the internal energy of a particle
  *
  * @param p The particle of interest
- * @param dt Time since the last kick
  */
 __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
-    const struct part *restrict p, float dt) {
-
-  const float entropy = p->entropy + dt * p->entropy_dt;
+    const struct part *restrict p) {
 
-  return gas_internal_energy_from_entropy(p->rho, entropy);
+  return gas_internal_energy_from_entropy(p->rho, p->entropy);
 }
 
 /**
  * @brief Returns the pressure of a particle
  *
  * @param p The particle of interest
- * @param dt Time since the last kick
  */
 __attribute__((always_inline)) INLINE static float hydro_get_pressure(
-    const struct part *restrict p, float dt) {
-
-  const float entropy = p->entropy + dt * p->entropy_dt;
+    const struct part *restrict p) {
 
-  return gas_pressure_from_entropy(p->rho, entropy);
+  return gas_pressure_from_entropy(p->rho, p->entropy);
 }
 
 /**
  * @brief Returns the entropy of a particle
  *
  * @param p The particle of interest
- * @param dt Time since the last kick
  */
 __attribute__((always_inline)) INLINE static float hydro_get_entropy(
-    const struct part *restrict p, float dt) {
+    const struct part *restrict p) {
 
-  return p->entropy + dt * p->entropy_dt;
+  return p->entropy;
 }
 
 /**
@@ -86,7 +79,7 @@ __attribute__((always_inline)) INLINE static float hydro_get_entropy(
  * @param dt Time since the last kick
  */
 __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
-    const struct part *restrict p, float dt) {
+    const struct part *restrict p) {
 
   return p->force.soundspeed;
 }
@@ -114,70 +107,32 @@ __attribute__((always_inline)) INLINE static float hydro_get_mass(
 }
 
 /**
- * @brief Modifies the thermal state of a particle to the imposed internal
- * energy
+ * @brief Returns the time derivative of internal energy of a particle
  *
- * This overwrites the current state of the particle but does *not* change its
- * time-derivatives. Entropy, pressure, sound-speed and signal velocity will be
- * updated.
+ * We assume a constant density.
  *
- * @param p The particle
- * @param u The new internal energy
+ * @param p The particle of interest
  */
-__attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
-    struct part *restrict p, float u) {
-
-  p->entropy = gas_entropy_from_internal_energy(p->rho, u);
-
-  /* Compute the new pressure */
-  const float pressure = gas_pressure_from_internal_energy(p->rho, u);
-
-  /* Compute the new sound speed */
-  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
-
-  /* Update the signal velocity */
-  const float v_sig_old = p->force.v_sig;
-  const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
-  const float v_sig = max(v_sig_old, v_sig_new);
-
-  const float rho_inv = 1.f / p->rho;
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy_dt(
+    const struct part *restrict p) {
 
-  p->force.soundspeed = soundspeed;
-  p->force.P_over_rho2 = pressure * rho_inv * rho_inv;
-  p->force.v_sig = v_sig;
+  return gas_internal_energy_from_entropy(p->rho, p->entropy_dt);
 }
 
 /**
- * @brief Modifies the thermal state of a particle to the imposed entropy
+ * @brief Returns the time derivative of internal energy of a particle
  *
- * This overwrites the current state of the particle but does *not* change its
- * time-derivatives. Entropy, pressure, sound-speed and signal velocity will be
- * updated.
+ * We assume a constant density.
  *
- * @param p The particle
- * @param S The new entropy
+ * @param p The particle of interest.
+ * @param du_dt The new time derivative of the internal energy.
  */
-__attribute__((always_inline)) INLINE static void hydro_set_entropy(
-    struct part *restrict p, float S) {
-
-  p->entropy = S;
-
-  /* Compute the pressure */
-  const float pressure = gas_pressure_from_entropy(p->rho, p->entropy);
-
-  /* Compute the new sound speed */
-  const float soundspeed = gas_soundspeed_from_pressure(p->rho, pressure);
-
-  /* Update the signal velocity */
-  const float v_sig_old = p->force.v_sig;
-  const float v_sig_new = p->force.v_sig - p->force.soundspeed + soundspeed;
-  const float v_sig = max(v_sig_old, v_sig_new);
+__attribute__((always_inline)) INLINE static void hydro_set_internal_energy_dt(
+    struct part *restrict p, float du_dt) {
 
-  const float rho_inv = 1.f / p->rho;
+  const float ds_dt = gas_entropy_from_internal_energy(p->rho, du_dt);
 
-  p->force.soundspeed = soundspeed;
-  p->force.P_over_rho2 = pressure * rho_inv * rho_inv;
-  p->force.v_sig = v_sig;
+  p->entropy_dt = ds_dt;
 }
 
 /**

src/hydro/Gadget2/hydro_debug.h

@@ -31,7 +31,7 @@ __attribute__((always_inline)) INLINE static void hydro_debug_particle(
       p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
       xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
       p->h, p->density.wcount, p->density.wcount_dh, p->mass, p->density.rho_dh,
-      p->rho, hydro_get_pressure(p, 0.), p->force.P_over_rho2, p->entropy,
+      p->rho, hydro_get_pressure(p), p->force.P_over_rho2, p->entropy,
       p->entropy_dt, p->force.soundspeed, p->density.div_v, p->density.rot_v[0],
       p->density.rot_v[1], p->density.rot_v[2], p->force.balsara,
       p->force.v_sig, p->force.h_dt, p->time_bin);

src/hydro/Gadget2/hydro_io.h

@@ -57,12 +57,12 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
 
 float convert_u(struct engine* e, struct part* p) {
 
-  return hydro_get_internal_energy(p, 0);
+  return hydro_get_internal_energy(p);
 }
 
 float convert_P(struct engine* e, struct part* p) {
 
-  return hydro_get_pressure(p, 0);
+  return hydro_get_pressure(p);
 }
 
 /**

src/hydro/Gizmo/hydro.h

@@ -378,7 +378,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
     struct part* p, struct xpart* xp, float dt) {
 
   float oldm, oldp[3], anew[3];
-  const float half_dt = 0.5f * dt; //MATTHIEU
+  const float half_dt = 0.5f * dt;  // MATTHIEU
 
   /* Retrieve the current value of the gravitational acceleration from the
      gpart. We are only allowed to do this because this is the kick. We still

src/runner.c

@@ -206,14 +206,17 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
   struct part *restrict parts = c->parts;
   struct xpart *restrict xparts = c->xparts;
   const int count = c->count;
-  // const int ti_current = r->e->ti_current;
-  const struct cooling_function_data *cooling_func = r->e->cooling_func;
-  const struct phys_const *constants = r->e->physical_constants;
-  const struct UnitSystem *us = r->e->internalUnits;
-  const double timeBase = r->e->timeBase;
+  const struct engine *e = r->e;
+  const struct cooling_function_data *cooling_func = e->cooling_func;
+  const struct phys_const *constants = e->physical_constants;
+  const struct UnitSystem *us = e->internalUnits;
+  const double timeBase = e->timeBase;
 
   TIMER_TIC;
 
+  /* Anything to do here? */
+  if (!cell_is_active(c, e)) return;
+
   /* Recurse? */
   if (c->split) {
     for (int k = 0; k < 8; k++)
@@ -227,13 +230,10 @@ void runner_do_cooling(struct runner *r, struct cell *c, int timer) {
       struct part *restrict p = &parts[i];
       struct xpart *restrict xp = &xparts[i];
 
-      /* Kick has already updated ti_end, so need to check ti_begin */
-      // if (p->ti_begin == ti_current) {
-      {
-
-        // const double dt = (p->ti_end - p->ti_begin) * timeBase;
-        const double dt = 1. * timeBase;  // MATTHIEU
+      if (part_is_active(p, e)) {
 
+        /* Let's cool ! */
+        const double dt = get_timestep(p->time_bin, timeBase);
         cooling_cool_part(constants, us, cooling_func, p, xp, dt);
       }
     }

src/statistics.c

@@ -162,14 +162,14 @@ 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, dt);
+    stats.E_int += m * hydro_get_internal_energy(p);
     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(
                                  time, potential, phys_const, gp);
     /* Collect entropy */
-    stats.entropy += m * hydro_get_entropy(p, dt);
+    stats.entropy += m * hydro_get_entropy(p);
   }
 
   /* Now write back to memory */