MPI communications for the SF split

examples/SmallCosmoVolume/SmallCosmoVolume_cooling/small_cosmo_volume.yml

@@ -51,6 +51,8 @@ Statistics:
 Scheduler:
   max_top_level_cells: 8
   cell_split_size:     50
+  cell_extra_gparts:         100         # (Optional) Number of spare gparts per top-level allocated at rebuild time for on-the-fly creation.
+  cell_extra_sparts:         100       # (Optional) Number of spare sparts per top-level allocated at rebuild time for on-the-fly creation.
   
 # Parameters related to the initial conditions
 InitialConditions:
@@ -96,22 +98,40 @@ EAGLEEntropyFloor:
   Cool_over_density_threshold:     10.       # Overdensity above which the EAGLE Cool limiter entropy floor can kick in.
   Cool_temperature_norm_K:         8000      # Temperature of the EAGLE Cool limiter entropy floor at the density threshold expressed in Kelvin.
   Cool_gamma_effective:            1.        # Slope the of the EAGLE Cool limiter entropy floor
-  
+
+GEARPressureFloor:
+  jeans_factor: 10.       # Number of particles required to suppose a resolved clump and avoid the pressure floor.
+
+
 # Cooling with Grackle 3.0
 GrackleCooling:
-  cloudy_table: CloudyData_UVB=HM2012.h5 # Name of the Cloudy Table (available on the grackle bitbucket repository)
-  with_UV_background: 1                   # Enable or not the UV background
-  redshift: -1                           # Redshift to use (-1 means time based redshift)
-  with_metal_cooling: 1                   # Enable or not the metal cooling
-  provide_volumetric_heating_rates: 0      # (optional) User provide volumetric heating rates
-  provide_specific_heating_rates: 0        # (optional) User provide specific heating rates
-  self_shielding_method: 0                # (optional) Grackle (<= 3) or Gear self shielding method
-  max_steps: 10000                       # (optional) Max number of step when computing the initial composition
-  convergence_limit: 1e-2                # (optional) Convergence threshold (relative) for initial composition
-  thermal_time_myr: 5
+  cloudy_table: CloudyData_UVB=HM2012.h5       # Name of the Cloudy Table (available on the grackle bitbucket repository)
+  with_UV_background: 1                        # Enable or not the UV background
+  redshift: 0                                  # Redshift to use (-1 means time based redshift)
+  with_metal_cooling: 1                        # Enable or not the metal cooling
+  provide_volumetric_heating_rates: 0          # (optional) User provide volumetric heating rates
+  provide_specific_heating_rates: 0            # (optional) User provide specific heating rates
+  max_steps: 10000                             # (optional) Max number of step when computing the initial composition
+  convergence_limit: 1e-2                      # (optional) Convergence threshold (relative) for initial composition
+  thermal_time_myr: 5                          # (optional) Time (in Myr) for adiabatic cooling after a feedback event.
+  self_shielding_method: -1                    # (optional) Grackle (1->3 for Grackle's ones, 0 for none and -1 for GEAR)
+  self_shielding_threshold_atom_per_cm3: 0.007 # Required only with GEAR's self shielding. Density threshold of the self shielding
+
 
+# GEAR chemistry model (Revaz and Jablonka 2018)
 GEARChemistry:
-  initial_metallicity: 0.01295
+  initial_metallicity: 1         # Initial metallicity of the gas (mass fraction)
+  scale_initial_metallicity: 1   # Should we scale the initial metallicity with the solar one?
 
-GEARPressureFloor:
-  jeans_factor: 10.       # Number of particles required to suppose a resolved clump and avoid the pressure floor.
+  # GEAR star formation model (Revaz and Jablonka 2018)
+GEARStarFormation:
+  star_formation_efficiency: 0.01   # star formation efficiency (c_*)
+  maximal_temperature:  3e4         # Upper limit to the temperature of a star forming particle
+  n_stars_per_particle: 4
+  min_mass_frac: 0.5
+
+# GEAR feedback model
+GEARFeedback:
+  supernovae_energy_erg: 0.1e51                            # Energy released by a single supernovae.
+  yields_table: chemistry-AGB+OMgSFeZnSrYBaEu-16072013.h5  # Table containing the yields.
+  discrete_yields: 0                                       # Should we use discrete yields or the IMF integrated one?

src/cell.c

@@ -233,6 +233,7 @@ int cell_link_gparts(struct cell *c, struct gpart *gparts) {
 #endif
 
   c->grav.parts = gparts;
+  c->grav.parts_rebuild = gparts;
 
   /* Fill the progeny recursively, depth-first. */
   if (c->split) {
@@ -1154,7 +1155,12 @@ int cell_pack_sf_counts(struct cell *restrict c,
   pcells[0].stars.count = c->stars.count;
   pcells[0].stars.dx_max_part = c->stars.dx_max_part;
 
+  /* Pack this cell's data. */
+  pcells[0].grav.delta_from_rebuild = c->grav.parts - c->grav.parts_rebuild;
+  pcells[0].grav.count = c->grav.count;
+
 #ifdef SWIFT_DEBUG_CHECKS
+  /* Stars */
   if (c->stars.parts_rebuild == NULL)
     error("Star particles array at rebuild is NULL! c->depth=%d", c->depth);
 
@@ -1163,6 +1169,16 @@ int cell_pack_sf_counts(struct cell *restrict c,
 
   if (pcells[0].stars.delta_from_rebuild > 0 && c->depth == 0)
     error("Shifting the top-level pointer is not allowed!");
+
+  /* Grav */
+  if (c->grav.parts_rebuild == NULL)
+    error("Grav. particles array at rebuild is NULL! c->depth=%d", c->depth);
+
+  if (pcells[0].grav.delta_from_rebuild < 0)
+    error("Grav part pointer moved in the wrong direction!");
+
+  if (pcells[0].grav.delta_from_rebuild > 0 && c->depth == 0)
+    error("Shifting the top-level pointer is not allowed!");
 #endif
 
   /* Fill in the progeny, depth-first recursion. */
@@ -1198,6 +1214,8 @@ int cell_unpack_sf_counts(struct cell *restrict c,
 #ifdef SWIFT_DEBUG_CHECKS
   if (c->stars.parts_rebuild == NULL)
     error("Star particles array at rebuild is NULL!");
+  if (c->grav.parts_rebuild == NULL)
+    error("Grav particles array at rebuild is NULL!");
 #endif
 
   /* Unpack this cell's data. */
@@ -1205,6 +1223,9 @@ int cell_unpack_sf_counts(struct cell *restrict c,
   c->stars.parts = c->stars.parts_rebuild + pcells[0].stars.delta_from_rebuild;
   c->stars.dx_max_part = pcells[0].stars.dx_max_part;
 
+  c->grav.count = pcells[0].grav.count;
+  c->grav.parts = c->grav.parts_rebuild + pcells[0].grav.delta_from_rebuild;
+
   /* Fill in the progeny, depth-first recursion. */
   int count = 1;
   for (int k = 0; k < 8; k++)
@@ -1967,6 +1988,7 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
     c->progeny[k]->grav.count = bucket_count[k];
     c->progeny[k]->grav.count_total = c->progeny[k]->grav.count;
     c->progeny[k]->grav.parts = &c->grav.parts[bucket_offset[k]];
+    c->progeny[k]->grav.parts_rebuild = c->progeny[k]->grav.parts;
   }
 }
 
@@ -5593,11 +5615,12 @@ struct spart *cell_add_spart(struct engine *e, struct cell *const c) {
   lock_lock(&top->stars.star_formation_lock);
 
   /* Are there any extra particles left? */
-  if (top->stars.count == top->stars.count_total - 1) {
+  if (top->stars.count >= top->stars.count_total - 1) {
     /* Release the local lock before exiting. */
     if (lock_unlock(&top->stars.star_formation_lock) != 0)
       error("Failed to unlock the top-level cell.");
-    message("We ran out of star particles!");
+    if (top->stars.count == top->stars.count_total - 1)
+      message("We ran out of star particles!");
     atomic_inc(&e->forcerebuild);
     return NULL;
   }
@@ -5725,11 +5748,12 @@ struct gpart *cell_add_gpart(struct engine *e, struct cell *c) {
   lock_lock(&top->grav.star_formation_lock);
 
   /* Are there any extra particles left? */
-  if (top->grav.count == top->grav.count_total - 1) {
+  if (top->grav.count >= top->grav.count_total - 1) {
     /* Release the local lock before exiting. */
     if (lock_unlock(&top->grav.star_formation_lock) != 0)
       error("Failed to unlock the top-level cell.");
-    message("We ran out of gravity particles!");
+    if (top->grav.count == top->grav.count_total - 1)
+      message("We ran out of gravity particles!");
     atomic_inc(&e->forcerebuild);
     return NULL;
   }

src/cell.h

@@ -273,6 +273,17 @@ struct pcell_sf {
     float dx_max_part;
 
   } stars;
+
+  /*! Grav. variables */
+  struct {
+
+    /* Distance by which the gpart pointer has moved since the last rebuild */
+    ptrdiff_t delta_from_rebuild;
+
+    /* Number of particles in the cell */
+    int count;
+
+  } grav;
 };
 
 /**
@@ -461,6 +472,9 @@ struct cell {
     /*! Pointer to the #gpart data. */
     struct gpart *parts;
 
+    /*! Pointer to the #spart data at rebuild time. */
+    struct gpart *parts_rebuild;
+
     /*! This cell's multipole. */
     struct gravity_tensors *multipole;
 

src/runner_others.c

@@ -267,6 +267,14 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
 
         /* Update current cell using child cells */
         star_formation_logger_add(&c->stars.sfh, &cp->stars.sfh);
+
+        /* Update the dx_max */
+        if (star_formation_need_update_dx_max) {
+          c->hydro.dx_max_part =
+              max(cp->hydro.dx_max_part, c->hydro.dx_max_part);
+          c->hydro.dx_max_sort =
+              max(cp->hydro.dx_max_sort, c->hydro.dx_max_sort);
+        }
       }
   } else {
 
@@ -359,6 +367,25 @@ void runner_do_star_formation(struct runner *r, struct cell *c, int timer) {
               /* Update the Star formation history */
               star_formation_logger_log_new_spart(sp, &c->stars.sfh);
 
+              /* Update the displacement information */
+              if (star_formation_need_update_dx_max) {
+                const float dx2_part = xp->x_diff[0] * xp->x_diff[0] +
+                                       xp->x_diff[1] * xp->x_diff[1] +
+                                       xp->x_diff[2] * xp->x_diff[2];
+                const float dx2_sort = xp->x_diff_sort[0] * xp->x_diff_sort[0] +
+                                       xp->x_diff_sort[1] * xp->x_diff_sort[1] +
+                                       xp->x_diff_sort[2] * xp->x_diff_sort[2];
+
+                const float dx_part = sqrtf(dx2_part);
+                const float dx_sort = sqrtf(dx2_sort);
+
+                /* No need to climb up the tree as the star formation starts
+                   higher in the hierarchy than the hydro and goes to the bottom
+                 */
+                c->hydro.dx_max_part = max(c->hydro.dx_max_part, dx_part);
+                c->hydro.dx_max_sort = max(c->hydro.dx_max_sort, dx_sort);
+              }
+
 #ifdef WITH_LOGGER
               /* Copy the properties back to the stellar particle */
               sp->logger_data = xp->logger_data;

src/space.c

@@ -265,6 +265,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
     c->hydro.parts = NULL;
     c->hydro.xparts = NULL;
     c->grav.parts = NULL;
+    c->grav.parts_rebuild = NULL;
     c->stars.parts = NULL;
     c->stars.parts_rebuild = NULL;
     c->black_holes.parts = NULL;
@@ -1954,6 +1955,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
 
       /* Store the state at rebuild time */
       c->stars.parts_rebuild = c->stars.parts;
+      c->grav.parts_rebuild = c->grav.parts;
 
       c->hydro.count_total = c->hydro.count + space_extra_parts;
       c->grav.count_total = c->grav.count + space_extra_gparts;

src/star_formation/EAGLE/star_formation.h

@@ -41,6 +41,8 @@
  * @brief Star formation model used in the EAGLE model
  */
 
+#define star_formation_need_update_dx_max 0
+
 /**
  * @brief Properties of the EAGLE star formation model.
  */

src/star_formation/GEAR/star_formation.h

@@ -34,6 +34,8 @@
 #include "star_formation_struct.h"
 #include "units.h"
 
+#define star_formation_need_update_dx_max 1
+
 /**
  * @brief Calculate if the gas has the potential of becoming
  * a star.
@@ -207,6 +209,77 @@ INLINE static void star_formation_update_part_not_SFR(
     struct part* p, struct xpart* xp, const struct engine* e,
     const struct star_formation* starform, const int with_cosmology) {}
 
+/**
+ * @brief Separate the #spart and #part by randomly moving both of them.
+ *
+ * @param e The #engine.
+ * @param p The #part generating a star.
+ * @param xp The #xpart generating a star.
+ * @param sp The new #spart.
+ */
+INLINE static void star_formation_separate_particles(const struct engine* e,
+                                                     struct part* p,
+                                                     struct xpart* xp,
+                                                     struct spart* sp) {
+#ifdef SWIFT_DEBUG_CHECKS
+  if (p->x[0] != sp->x[0] || p->x[1] != sp->x[1] || p->x[2] != sp->x[2]) {
+    error(
+        "Moving particles that are not at the same location."
+        " (%g, %g, %g) - (%g, %g, %g)",
+        p->x[0], p->x[1], p->x[2], sp->x[0], sp->x[1], sp->x[2]);
+  }
+#endif
+
+  /* Move a bit the particle in order to avoid
+     division by 0.
+  */
+  const float max_displacement = 0.2;
+  const double delta_x =
+      2.f * random_unit_interval(p->id, e->ti_current,
+                                 (enum random_number_type)0) -
+      1.f;
+  const double delta_y =
+      2.f * random_unit_interval(p->id, e->ti_current,
+                                 (enum random_number_type)1) -
+      1.f;
+  const double delta_z =
+      2.f * random_unit_interval(p->id, e->ti_current,
+                                 (enum random_number_type)2) -
+      1.f;
+
+  sp->x[0] += delta_x * max_displacement * p->h;
+  sp->x[1] += delta_y * max_displacement * p->h;
+  sp->x[2] += delta_z * max_displacement * p->h;
+
+  /* Copy the position to the gpart */
+  sp->gpart->x[0] = sp->x[0];
+  sp->gpart->x[1] = sp->x[1];
+  sp->gpart->x[2] = sp->x[2];
+
+  /* Do the gas particle. */
+  const double mass_ratio = sp->mass / hydro_get_mass(p);
+  const double dx[3] = {mass_ratio * delta_x * max_displacement * p->h,
+                        mass_ratio * delta_y * max_displacement * p->h,
+                        mass_ratio * delta_z * max_displacement * p->h};
+
+  p->x[0] -= dx[0];
+  p->x[1] -= dx[1];
+  p->x[2] -= dx[2];
+
+  /* Compute offsets since last cell construction */
+  xp->x_diff[0] += dx[0];
+  xp->x_diff[1] += dx[1];
+  xp->x_diff[1] += dx[2];
+  xp->x_diff_sort[0] += dx[0];
+  xp->x_diff_sort[1] += dx[1];
+  xp->x_diff_sort[2] += dx[2];
+
+  /* Copy the position to the gpart */
+  p->gpart->x[0] = p->x[0];
+  p->gpart->x[1] = p->x[1];
+  p->gpart->x[2] = p->x[2];
+}
+
 /**
  * @brief Copies the properties of the gas particle over to the
  * star particle.
@@ -225,10 +298,9 @@ INLINE static void star_formation_update_part_not_SFR(
  * @param convert_part Did we convert a part (or spawned one)?
  */
 INLINE static void star_formation_copy_properties(
-    struct part* p, const struct xpart* xp, struct spart* sp,
-    const struct engine* e, const struct star_formation* starform,
-    const struct cosmology* cosmo, const int with_cosmology,
-    const struct phys_const* phys_const,
+    struct part* p, struct xpart* xp, struct spart* sp, const struct engine* e,
+    const struct star_formation* starform, const struct cosmology* cosmo,
+    const int with_cosmology, const struct phys_const* phys_const,
     const struct hydro_props* restrict hydro_props,
     const struct unit_system* restrict us,
     const struct cooling_function_data* restrict cooling,
@@ -247,6 +319,8 @@ INLINE static void star_formation_copy_properties(
     /* Update the part */
     hydro_set_mass(p, new_mass_gas);
     p->gpart->mass = new_mass_gas;
+
+    star_formation_separate_particles(e, p, xp, sp);
   } else {
     sp->mass = mass_gas;
   }

src/star_formation/QLA/star_formation.h

@@ -34,6 +34,8 @@
  * @brief Star formation model used in the EAGLE model
  */
 
+#define star_formation_need_update_dx_max 0
+
 /**
  * @brief Properties of the EAGLE star formation model.
  */

src/star_formation/none/star_formation.h

@@ -32,6 +32,11 @@
 /* Starformation struct */
 struct star_formation {};
 
+/* Does the star formation model move the hydro particles?
+   This will update the c->hydro.dx_max_part and
+   c->hydro.dx_max_sort after forming a star. */
+#define star_formation_need_update_dx_max 0
+
 /**
  * @brief Calculate if the gas has the potential of becoming
  * a star.