Add a new task to do the resyncing of feedback-hit particles.

examples/EAGLE_low_z/EAGLE_12/eagle_12.yml

@@ -141,7 +141,7 @@ EAGLEFeedback:
   IMF_max_mass_Msun:              100.0             # Maximal stellar mass considered for the Chabrier IMF in solar masses.
   SNII_min_mass_Msun:               6.0             # Minimal mass considered for SNII feedback (not SNII enrichment!) in solar masses.
   SNII_max_mass_Msun:             100.0             # Maximal mass considered for SNII feedback (not SNII enrichment!) in solar masses.
-  SNII_wind_delay_Gyr:              0.03            # Time in Gyr between a star's birth and the SNII thermal feedback event.
+  SNII_wind_delay_Gyr:              0.001           # Time in Gyr between a star's birth and the SNII thermal feedback event.
   SNII_delta_T_K:                   3.16228e7       # Change in temperature to apply to the gas particle in a SNII thermal feedback event in Kelvin.
   SNII_energy_erg:                  1.0e51          # Energy of one SNII explosion in ergs.
   SNII_energy_fraction_min:         0.3             # Minimal fraction of energy applied in a SNII feedback event.

examples/IsolatedGalaxy/IsolatedGalaxy_feedback/isolated_galaxy.yml

@@ -125,7 +125,7 @@ EAGLEFeedback:
   IMF_max_mass_Msun:              100.0             # Maximal stellar mass considered for the Chabrier IMF in solar masses.
   SNII_min_mass_Msun:               6.0             # Minimal mass considered for SNII feedback (not SNII enrichment!) in solar masses.
   SNII_max_mass_Msun:             100.0             # Maximal mass considered for SNII feedback (not SNII enrichment!) in solar masses.
-  SNII_wind_delay_Gyr:              0.001            # Time in Gyr between a star's birth and the SNII thermal feedback event.
+  SNII_wind_delay_Gyr:              0.0001           # Time in Gyr between a star's birth and the SNII thermal feedback event.
   SNII_delta_T_K:                   3.16228e7       # Change in temperature to apply to the gas particle in a SNII thermal feedback event in Kelvin.
   SNII_energy_erg:                  1.0e51          # Energy of one SNII explosion in ergs.
   SNII_energy_fraction_min:         1.0             # Minimal fraction of energy applied in a SNII feedback event.

examples/IsolatedGalaxy/IsolatedGalaxy_feedback/plotSolution.py

@@ -119,18 +119,18 @@ SF_thresh = KS_thresh_norm * (EAGLE_Z / KS_thresh_Z0) ** (KS_thresh_slope)
 gas_pos = f["/PartType0/Coordinates"][:, :]
 gas_mass = f["/PartType0/Masses"][:]
 gas_rho = f["/PartType0/Densities"][:]
-gas_T = f["/PartType0/Temperature"][:]
-gas_SFR = f["/PartType0/SFR"][:]
+gas_T = f["/PartType0/Temperatures"][:]
+gas_SFR = f["/PartType0/StarFormationRates"][:]
 gas_XH = f["/PartType0/ElementMassFractions"][:, 0]
-gas_Z = f["/PartType0/Metallicities"][:]
+gas_Z = f["/PartType0/MetalMassFractions"][:]
 gas_hsml = f["/PartType0/SmoothingLengths"][:]
 gas_sSFR = gas_SFR / gas_mass
 
 # Read the Star properties
 stars_pos = f["/PartType4/Coordinates"][:, :]
-stars_BirthDensity = f["/PartType4/BirthDensity"][:]
-stars_BirthTime = f["/PartType4/BirthTime"][:]
-stars_XH = f["/PartType4/ElementAbundance"][:, 0]
+stars_BirthDensity = f["/PartType4/BirthDensities"][:]
+stars_BirthTime = f["/PartType4/BirthTimes"][:]
+stars_XH = f["/PartType4/ElementMassFractions"][:, 0]
 
 # Centre the box
 gas_pos[:, 0] -= centre[0]

src/cell.h

@@ -293,6 +293,8 @@ enum cell_flags {
   cell_flag_do_bh_sub_drift = (1UL << 12),
   cell_flag_do_stars_resort = (1UL << 13),
   cell_flag_has_tasks = (1UL << 14),
+  cell_flag_do_hydro_sync = (1UL << 15),
+  cell_flag_do_hydro_sub_sync = (1UL << 16)
 };
 
 /**
@@ -764,6 +766,9 @@ struct cell {
   /*! The task to limit the time-step of inactive particles */
   struct task *timestep_limiter;
 
+  /*! The task to synchronize the time-step of inactive particles hit by feedback */
+  struct task *timestep_sync;
+  
 #ifdef WITH_LOGGER
   /*! The logger task */
   struct task *logger;

src/engine.c

@@ -2167,11 +2167,13 @@ void engine_step(struct engine *e) {
 
     /* Print some information to the screen */
     printf(
-        "  %6d %14e %12.7f %12.7f %14e %4d %4d %12lld %12lld %12lld %12lld "
+        "  %6d %14e %lld %12.7f %12.7f %14e %4d %4d %12lld %12lld %12lld "
+        "%12lld "
         "%21.3f %6d\n",
-        e->step, e->time, e->cosmology->a, e->cosmology->z, e->time_step,
-        e->min_active_bin, e->max_active_bin, e->updates, e->g_updates,
-        e->s_updates, e->b_updates, e->wallclock_time, e->step_props);
+        e->step, e->time, e->ti_current, e->cosmology->a, e->cosmology->z,
+        e->time_step, e->min_active_bin, e->max_active_bin, e->updates,
+        e->g_updates, e->s_updates, e->b_updates, e->wallclock_time,
+        e->step_props);
 #ifdef SWIFT_DEBUG_CHECKS
     fflush(stdout);
 #endif

src/engine_maketasks.c

@@ -849,8 +849,13 @@ void engine_make_hierarchical_tasks_common(struct engine *e, struct cell *c) {
         c->timestep_limiter = scheduler_addtask(
             s, task_type_timestep_limiter, task_subtype_none, 0, 0, c, NULL);
 
+        c->timestep_sync = scheduler_addtask(
+            s, task_type_timestep_sync, task_subtype_none, 0, 0, c, NULL);
+	
         scheduler_addunlock(s, c->timestep, c->timestep_limiter);
         scheduler_addunlock(s, c->timestep_limiter, c->kick1);
+        scheduler_addunlock(s, c->timestep, c->timestep_sync);
+        scheduler_addunlock(s, c->timestep_sync, c->kick1);
       }
     }
   } else { /* We are above the super-cell so need to go deeper */
@@ -1182,6 +1187,8 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c,
         scheduler_addunlock(s, c->hydro.limiter_out, c->super->timestep);
         scheduler_addunlock(s, c->hydro.limiter_out,
                             c->super->timestep_limiter);
+        scheduler_addunlock(s, c->hydro.limiter_out,
+                            c->super->timestep_sync);
 
         if (with_star_formation && c->hydro.count > 0) {
           scheduler_addunlock(s, c->top->hydro.star_formation,

src/engine_marktasks.c

@@ -965,8 +965,6 @@ int engine_marktasks(struct engine *e) {
   const ticks tic = getticks();
   int rebuild_space = 0;
 
-  message("MARKTASKS!");
-
   /* Run through the tasks and mark as skip or not. */
   size_t extra_data[3] = {(size_t)e, (size_t)rebuild_space, (size_t)&e->sched};
   threadpool_map(&e->threadpool, engine_marktasks_mapper, s->tasks, s->nr_tasks,

src/feedback/EAGLE/feedback.c

@@ -155,6 +155,8 @@ INLINE static void compute_SNII_feedback(
       return;
     }
 
+    message("Doing SNII feedback for star %lld", sp->id);
+
     /* Properties of the model (all in internal units) */
     const double delta_T =
         eagle_feedback_temperature_change(sp, feedback_props);

src/feedback/EAGLE/feedback_iact.h

@@ -34,6 +34,7 @@
  * @param pj Second particle (not updated).
  * @param xpj Extra particle data (not updated).
  * @param cosmo The cosmological model.
+ * @param e The #engine (for particle synchronization).
  * @param ti_current Current integer time value
  */
 __attribute__((always_inline)) INLINE static void
@@ -82,6 +83,7 @@ runner_iact_nonsym_feedback_density(const float r2, const float *dx,
  * @param pj Second (gas) particle.
  * @param xpj Extra particle data
  * @param cosmo The cosmological model.
+ * @param e The #engine (for particle synchronization).
  * @param ti_current Current integer time used value for seeding random number
  * generator
  */
@@ -298,7 +300,9 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx,
           pj->id, si->id, prob, rand, delta_u, delta_u / u_init);
 
       /* Synchronize the particle on the timeline */
-      timestep_sync_part(pj, xpj, e, cosmo);
+      timestep_sync_part(pj);
+
+      /* Recompute the signal velocity */
     }
   }
 }

src/feedback/none/feedback_iact.h

@@ -32,6 +32,7 @@
  * @param pj Second particle (not updated).
  * @param xp Extra particle data (not updated).
  * @param cosmo The cosmological model.
+ * @param e The #engine (for particle synchronization).
  * @param ti_current Current integer time value
  */
 __attribute__((always_inline)) INLINE static void
@@ -57,6 +58,7 @@ runner_iact_nonsym_feedback_density(const float r2, const float *dx,
  * @param pj Second (gas) particle.
  * @param xp Extra particle data
  * @param cosmo The cosmological model.
+ * @param e The #engine (for particle synchronization).
  * @param ti_current Current integer time used value for seeding random number
  * generator
  */

src/hydro/Minimal/hydro.h

@@ -870,6 +870,7 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
   p->time_bin = 0;
   p->min_ngb_time_bin = num_time_bins + 1;
   p->wakeup = time_bin_not_awake;
+  p->to_be_synchronized = 0;
   xp->v_full[0] = p->v[0];
   xp->v_full[1] = p->v[1];
   xp->v_full[2] = p->v[2];

src/hydro/Minimal/hydro_part.h

@@ -182,10 +182,13 @@ struct part {
 
   /* Need waking-up ? */
   timebin_t wakeup;
-
+  
   /*! Minimal time-bin across all neighbours */
   timebin_t min_ngb_time_bin;
 
+  /* Do we want this particle to be synched back on the time-line? */
+  char to_be_synchronized;
+  
 #ifdef SWIFT_DEBUG_CHECKS
 
   /* Time of the last drift */

src/runner.h

@@ -107,6 +107,7 @@ void runner_do_end_grav_force(struct runner *r, struct cell *c, int timer);
 void runner_do_init(struct runner *r, struct cell *c, int timer);
 void runner_do_cooling(struct runner *r, struct cell *c, int timer);
 void runner_do_limiter(struct runner *r, struct cell *c, int force, int timer);
+void runner_do_sync(struct runner *r, struct cell *c, int force, int timer);
 void runner_do_grav_mesh(struct runner *r, struct cell *c, int timer);
 void runner_do_grav_external(struct runner *r, struct cell *c, int timer);
 void runner_do_grav_fft(struct runner *r, int timer);

src/runner_main.c

@@ -367,6 +367,9 @@ void *runner_main(void *data) {
         case task_type_timestep_limiter:
           runner_do_limiter(r, ci, 0, 1);
           break;
+        case task_type_timestep_sync:
+          runner_do_sync(r, ci, 0, 1);
+          break;
 #ifdef WITH_MPI
         case task_type_send:
           if (t->subtype == task_subtype_tend_part) {

src/runner_time_integration.c

@@ -34,6 +34,7 @@
 #include "timers.h"
 #include "timestep.h"
 #include "timestep_limiter.h"
+#include "timestep_sync.h"
 #include "tracers.h"
 
 /**
@@ -1077,3 +1078,77 @@ void runner_do_limiter(struct runner *r, struct cell *c, int force, int timer) {
 
   if (timer) TIMER_TOC(timer_do_limiter);
 }
+
+void runner_do_sync(struct runner *r, struct cell *c, int force, int timer) {
+  
+  const struct engine *e = r->e;
+  const struct cosmology *cosmo = e->cosmology;
+  const int count = c->hydro.count;
+  struct part *restrict parts = c->hydro.parts;
+  struct xpart *restrict xparts = c->hydro.xparts;
+
+  TIMER_TIC;
+
+#ifdef SWIFT_DEBUG_CHECKS
+  /* Check that we only sync local cells. */
+  if (c->nodeID != engine_rank) error("Syncing of a foreign cell is nope.");
+#endif
+
+  /* integertime_t ti_hydro_end_min = max_nr_timesteps, ti_hydro_end_max = 0, */
+  /*               ti_hydro_beg_max = 0; */
+  /* integertime_t ti_gravity_end_min = max_nr_timesteps, ti_gravity_end_max = 0, */
+  /*               ti_gravity_beg_max = 0; */
+
+  /* Limit irrespective of cell flags? */
+  force = (force || cell_get_flag(c, cell_flag_do_hydro_sync));
+
+  /* Early abort? */
+  if (c->hydro.count == 0) {
+
+    /* Clear the sync flags. */
+    cell_clear_flag(c, cell_flag_do_hydro_sync | cell_flag_do_hydro_sub_sync);
+    return;
+  }
+
+  /* Loop over the progeny ? */
+  if (c->split && (force || cell_get_flag(c, cell_flag_do_hydro_sub_sync))) {
+    for (int k = 0; k < 8; k++) {
+      if (c->progeny[k] != NULL) {
+        struct cell *restrict cp = c->progeny[k];
+
+	/* Recurse */
+	runner_do_sync(r, cp, force, 0);
+      }
+    }
+    
+  } else if (!c->split && force) {
+
+    /* Loop over the gas particles in this cell. */
+    for (int k = 0; k < count; k++) {
+
+      /* Get a handle on the part. */
+      struct part *restrict p = &parts[k];
+      struct xpart *restrict xp = &xparts[k];
+      
+      /* Avoid inhibited particles */
+      if (part_is_inhibited(p, e)) continue;
+      
+      /* If the particle is active no need to sync it */
+      if (part_is_active(p, e) && p->to_be_synchronized)
+	p->to_be_synchronized = 0;
+	      
+      if (p->to_be_synchronized) {
+
+	timestep_process_sync_part(p, xp, e, cosmo);
+
+      }     
+    }    
+  }
+
+  /* Clear the sync flags. */
+  cell_clear_flag(c,
+                  cell_flag_do_hydro_sync | cell_flag_do_hydro_sub_sync);
+
+  if (timer) TIMER_TOC(timer_do_sync);
+
+}

src/space.c

@@ -254,6 +254,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
     c->hydro.limiter_in = NULL;
     c->hydro.limiter_out = NULL;
     c->timestep_limiter = NULL;
+    c->timestep_sync = NULL;
     c->hydro.end_force = NULL;
     c->hydro.drift = NULL;
     c->stars.drift = NULL;
@@ -5297,9 +5298,9 @@ void space_check_limiter_mapper(void *map_data, int nr_parts,
     if (parts[k].wakeup == time_bin_awake)
       error("Particle still woken up! id=%lld", parts[k].id);
 
-    if (parts[k].synchronized != 0)
+    if (parts[k].to_be_synchronized != 0)
       error("Synchronized particle not treated! id=%lld synchronized=%d",
-            parts[k].id, parts[k].synchronized);
+            parts[k].id, parts[k].to_be_synchronized);
 
     if (parts[k].gpart != NULL)
       if (parts[k].time_bin != parts[k].gpart->time_bin)

src/task.c

@@ -71,6 +71,7 @@ const char *taskID_names[task_type_count] = {"none",
                                              "kick2",
                                              "timestep",
                                              "timestep_limiter",
+					     "timestep_sync",
                                              "limiter_in",
                                              "limiter_out",
                                              "send",
@@ -433,6 +434,7 @@ void task_unlock(struct task *t) {
     case task_type_extra_ghost:
     case task_type_end_hydro_force:
     case task_type_timestep_limiter:
+    case task_type_timestep_sync:
       cell_unlocktree(ci);
       break;
 
@@ -585,6 +587,7 @@ int task_lock(struct task *t) {
     case task_type_extra_ghost:
     case task_type_end_hydro_force:
     case task_type_timestep_limiter:
+    case task_type_timestep_sync:
       if (ci->hydro.hold) return 0;
       if (cell_locktree(ci) != 0) return 0;
       break;

src/task.h

@@ -65,6 +65,7 @@ enum task_types {
   task_type_kick2,
   task_type_timestep,
   task_type_timestep_limiter,
+  task_type_timestep_sync,
   task_type_limiter_in,  /* Implicit */
   task_type_limiter_out, /* Implicit */
   task_type_send,

src/timeline.h

@@ -33,7 +33,7 @@ typedef long long integertime_t;
 typedef int8_t timebin_t;
 
 /*! The number of time bins */
-#define num_time_bins 56
+#define num_time_bins 36
 
 /*! The maximal number of timesteps in a simulation */
 #define max_nr_timesteps (1LL << (num_time_bins + 1))

src/timestep.h

@@ -222,7 +222,7 @@ __attribute__((always_inline)) INLINE static integertime_t get_spart_timestep(
   new_dt *= e->cosmology->time_step_factor;
 
   /* Limit timestep within the allowed range */
-  new_dt = min(new_dt, e->dt_max);
+  new_dt = min(new_dt, e->dt_max) / 8.;
   if (new_dt < e->dt_min) {
     error("spart (id=%lld) wants a time-step (%e) below dt_min (%e)", sp->id,
           new_dt, e->dt_min);