Split stars

examples/DwarfGalaxy/README

+This example is a galaxy extracted from the example "ZoomIn". It allows
+to test SWIFT on a smaller problem. See the README in "ZoomIn" for more
+information.
+
+
+MD5 check-sum of the ICS: 
+ae2af84d88f30011b6a8af3f37d140cf  dwarf_galaxy.hdf5
\ No newline at end of file

examples/DwarfGalaxy/dwarf_galaxy.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e21 # kpc in centimeters
+  UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
+  UnitCurrent_in_cgs:  1             # Amperes
+  UnitTemp_in_cgs:     1             # Kelvin
+
+# Structure finding options
+StructureFinding:
+  config_file_name:     stf_input.cfg # Name of the STF config file.
+  basename:             ./stf         # Common part of the name of output files.
+  output_time_format:   0             # Specifies the frequency format of structure finding. 0 for simulation steps (delta_step) and 1 for simulation time intervals (delta_time).
+  scale_factor_first:   0.92          # Scale-factor of the first snaphot (cosmological run)
+  time_first:           0.01        # Time of the first structure finding output (in internal units).
+  delta_step:           1000          # Time difference between consecutive structure finding outputs (in internal units) in simulation steps.
+  delta_time:           1.10          # Time difference between consecutive structure finding outputs (in internal units) in simulation time intervals.
+
+# Cosmological parameters
+Cosmology:
+  h:              0.673        # Reduced Hubble constant
+  a_begin:        0.9873046739     # Initial scale-factor of the simulation
+  a_end:          1.0           # Final scale factor of the simulation
+  Omega_m:        0.315         # Matter density parameter
+  Omega_lambda:   0.685         # Dark-energy density parameter
+  Omega_b:        0.0486        # Baryon density parameter
+  
+Scheduler:
+  max_top_level_cells:    8
+  cell_split_size:           400       # (Optional) Maximal number of particles per cell (this is the default value).
+  
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   1.  # The end time of the simulation (in internal units).
+  dt_min:     1e-10 # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e-3  # The maximal time-step size of the simulation (in internal units).
+  
+# Parameters governing the snapshots
+Snapshots:
+  basename:            dwarf_galaxy # Common part of the name of output files
+  time_first:          0.  # Time of the first output (non-cosmological run) (in internal units)
+  delta_time:          0.02  # Time difference between consecutive outputs (in internal units)
+  compression:         1
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  scale_factor_first:  0.987345 # Scale-factor of the first stat dump (cosmological run)
+  time_first:          0.01 # Time of the first stat dump (non-cosmological run) (in internal units)
+  delta_time:          1.05 # Time between statistics output
+
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:                    0.025    # Constant dimensionless multiplier for time integration.
+  theta:                  0.7      # Opening angle (Multipole acceptance criterion)
+  comoving_softening:     0.05 # Comoving softening length (in internal units).
+  max_physical_softening: 0.01    # Physical softening length (in internal units).
+  mesh_side_length:       16
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature:   100      # (internal units)
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./dwarf_galaxy.hdf5     # The file to read
+  cleanup_h_factors: 1               # Remove the h-factors inherited from Gadget
+  cleanup_velocity_factors: 1        # Remove the sqrt(a) factor in the velocities inherited from Gadget
+

examples/DwarfGalaxy/getIC.sh

+#!/bin/bash
+wget https://obswww.unige.ch/~lhausamm/swift/IC/DwarfGalaxy/dwarf_galaxy.hdf5

examples/DwarfGalaxy/run.sh

+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e dwarf_galaxy.hdf5 ]
+then
+    echo "Fetching initial conditions for the dwarf galaxy example..."
+    ./getIC.sh
+fi
+
+../swift -b -G -s -S -t 8 dwarf_galaxy.yml 2>&1 | tee output.log
+

examples/ZoomIn/README

+Initial conditions for a zoom in cosmological simulation of dwarf
+galaxies. These have been generated by MUSIC and ran up to z=0 with
+GEAR (see Revaz and Jablonka 2018 for more details on the simulation).
+
+The cosmology is taken from Planck 2015.
+
+The initial conditions have been cleaned to contain only the required
+fields. The ICs have been created for Gadget and the positions and box
+size are hence expressed in h-full units (e.g. box size of 32 / h Mpc).
+Similarly, the peculiar velocitites contain an extra sqrt(a) factor. 
+
+We will use SWIFT to cancel the h- and a-factors from the ICs. Gas
+particles will be generated at startup.
+
+MD5 check-sum of the ICS: 
+9aafe154438478ed435e88664c1c5dba zoom_in.hdf5

examples/ZoomIn/getIC.sh

+#!/bin/bash
+wget https://obswww.unige.ch/~lhausamm/swift/IC/ZoomIn/zoom_in.hdf5

examples/ZoomIn/run.sh

+#!/bin/bash
+
+ # Generate the initial conditions if they are not present.
+if [ ! -e zoom_in.hdf5 ]
+then
+    echo "Fetching initial conditions for the zoom in example..."
+    ./getIC.sh
+fi
+
+../swift -b -c -G -s -S -t 8 zoom_in.yml 2>&1 | tee output.log
+

examples/ZoomIn/zoom_in.yml

+# Define the system of units to use internally. 
+InternalUnitSystem:
+  UnitMass_in_cgs:     1.98848e43    # 10^10 M_sun in grams
+  UnitLength_in_cgs:   3.08567758e21 # kpc in centimeters
+  UnitVelocity_in_cgs: 1e5           # km/s in centimeters per second
+  UnitCurrent_in_cgs:  1             # Amperes
+  UnitTemp_in_cgs:     1             # Kelvin
+
+# Cosmological parameters
+Cosmology:
+  h:              0.673        # Reduced Hubble constant
+  a_begin:        0.9873046739     # Initial scale-factor of the simulation
+  a_end:          1.0           # Final scale factor of the simulation
+  Omega_m:        0.315         # Matter density parameter
+  Omega_lambda:   0.685         # Dark-energy density parameter
+  Omega_b:        0.0486        # Baryon density parameter
+  
+Scheduler:
+  max_top_level_cells:    8
+  
+# Parameters governing the time integration
+TimeIntegration:
+  time_begin: 0.    # The starting time of the simulation (in internal units).
+  time_end:   1e-2  # The end time of the simulation (in internal units).
+  dt_min:     1e-10 # The minimal time-step size of the simulation (in internal units).
+  dt_max:     1e-3  # The maximal time-step size of the simulation (in internal units).
+  
+# Parameters governing the snapshots
+Snapshots:
+  basename:            zoom_in # Common part of the name of output files
+  scale_factor_first:  0.987345  # Scale-factor of the first snaphot (cosmological run)
+  time_first:          0.01  # Time of the first output (non-cosmological run) (in internal units)
+  delta_time:          1.01  # Time difference between consecutive outputs (in internal units)
+  compression:         1
+
+# Parameters governing the conserved quantities statistics
+Statistics:
+  scale_factor_first:  0.987345 # Scale-factor of the first stat dump (cosmological run)
+  time_first:          0.01 # Time of the first stat dump (non-cosmological run) (in internal units)
+  delta_time:          1.05 # Time between statistics output
+
+# Parameters for the self-gravity scheme
+Gravity:
+  eta:                    0.025    # Constant dimensionless multiplier for time integration.
+  theta:                  0.7      # Opening angle (Multipole acceptance criterion)
+  comoving_softening:     0.05 # Comoving softening length (in internal units).
+  max_physical_softening: 0.01    # Physical softening length (in internal units).
+  mesh_side_length:       16
+
+# Parameters for the hydrodynamics scheme
+SPH:
+  resolution_eta:        1.2348   # Target smoothing length in units of the mean inter-particle separation (1.2348 == 48Ngbs with the cubic spline kernel).
+  CFL_condition:         0.1      # Courant-Friedrich-Levy condition for time integration.
+  minimal_temperature:   100      # (internal units)
+
+# Parameters related to the initial conditions
+InitialConditions:
+  file_name:  ./zoom_in.hdf5     # The file to read
+  cleanup_h_factors: 1               # Remove the h-factors inherited from Gadget
+  cleanup_velocity_factors: 1        # Remove the sqrt(a) factor in the velocities inherited from Gadget
+

examples/plot_task_level.py

+#!/usr/bin/python
+"""
+Usage:
+  ./plot_task_level.py task_level.txt
+
+Description:
+  Plot the number of tasks for each depth level and each type of task.
+"""
+
+
+import pandas as pd
+import matplotlib.pyplot as plt
+import sys
+
+# get filename
+filename = sys.argv[-1]
+
+# Column names
+names = ["type", "subtype", "depth", "count"]
+
+# read file
+data = pd.read_csv(filename, sep=' ', comment="#", names=names)
+
+# generate color map
+cmap = plt.get_cmap("hsv")
+N = data["depth"].max() + 5
+
+# plot data
+for i in range(data["depth"].max()):
+    ind = data["depth"] == i
+    label = "depth = %i" % i
+    c = cmap(i / N)
+    plt.plot(data["type"][ind] + "_" + data["subtype"][ind],
+             data["count"][ind], ".", label=label,
+             color=c)
+
+# modify figure parameters and show it
+plt.gca().set_yscale("log")
+plt.xticks(rotation=45)
+plt.ylabel("Number of Tasks")
+plt.gcf().subplots_adjust(bottom=0.15)
+plt.legend()
+plt.show()

examples/plot_tasks.py

@@ -157,11 +157,11 @@ for task in SUBTYPES:
 
 #  For fiddling with colours...
 if args.verbose:
-    print "#Selected colours:"
+    print("#Selected colours:")
     for task in sorted(TASKCOLOURS.keys()):
-        print "# " + task + ": " + TASKCOLOURS[task]
+        print("# " + task + ": " + TASKCOLOURS[task])
     for task in sorted(SUBCOLOURS.keys()):
-        print "# " + task + ": " + SUBCOLOURS[task]
+        print("# " + task + ": " + SUBCOLOURS[task])
 
 #  Read input.
 data = pl.loadtxt( infile )
@@ -169,11 +169,11 @@ data = pl.loadtxt( infile )
 #  Do we have an MPI file?
 full_step = data[0,:]
 if full_step.size == 13:
-    print "# MPI mode"
+    print("# MPI mode")
     mpimode = True
     if ranks == None:
         ranks = range(int(max(data[:,0])) + 1)
-    print "# Number of ranks:", len(ranks)
+    print("# Number of ranks:", len(ranks))
     rankcol = 0
     threadscol = 1
     taskcol = 2
@@ -181,7 +181,7 @@ if full_step.size == 13:
     ticcol = 5
     toccol = 6
 else:
-    print "# non MPI mode"
+    print("# non MPI mode")
     ranks = [0]
     mpimode = False
     rankcol = -1
@@ -194,10 +194,10 @@ else:
 #  Get CPU_CLOCK to convert ticks into milliseconds.
 CPU_CLOCK = float(full_step[-1]) / 1000.0
 if args.verbose:
-    print "# CPU frequency:", CPU_CLOCK * 1000.0
+    print("# CPU frequency:", CPU_CLOCK * 1000.0)
 
 nthread = int(max(data[:,threadscol])) + 1
-print "# Number of threads:", nthread
+print("# Number of threads:", nthread)
 
 # Avoid start and end times of zero.
 sdata = data[data[:,ticcol] != 0]
@@ -224,24 +224,24 @@ if delta_t == 0:
         dt = toc_step - tic_step
         if dt > delta_t:
             delta_t = dt
-    print "# Data range: ", delta_t / CPU_CLOCK, "ms"
+    print("# Data range: ", delta_t / CPU_CLOCK, "ms")
 
 # Once more doing the real gather and plots this time.
 for rank in ranks:
-    print "# Processing rank: ", rank
+    print("# Processing rank: ", rank)
     if mpimode:
         data = sdata[sdata[:,rankcol] == rank]
         full_step = data[0,:]
     tic_step = int(full_step[ticcol])
     toc_step = int(full_step[toccol])
-    print "# Min tic = ", tic_step
+    print("# Min tic = ", tic_step)
     data = data[1:,:]
     typesseen = []
     nethread = 0
 
     #  Dummy image for ranks that have no tasks.
     if data.size == 0:
-        print "# Rank ", rank, " has no tasks"
+        print("# Rank ", rank, " has no tasks")
         fig = pl.figure()
         ax = fig.add_subplot(1,1,1)
         ax.set_xlim(-delta_t * 0.01 / CPU_CLOCK, delta_t * 1.01 / CPU_CLOCK)
@@ -367,6 +367,6 @@ for rank in ranks:
     else:
         outpng = outbase + ".png"
     pl.savefig(outpng)
-    print "Graphics done, output written to", outpng
+    print("Graphics done, output written to", outpng)
 
 sys.exit(0)

src/cell.h

@@ -144,6 +144,9 @@ struct pcell {
     /*! Number of #spart in this cell. */
     int count;
 
+    /*! Maximal smoothing length. */
+    double h_max;
+
   } stars;
 
   /*! Relative indices of the cell's progeny. */
@@ -188,6 +191,10 @@ struct pcell_step {
 
   /*! Stars variables */
   struct {
+
+    /*! Maximal distance any #part has travelled since last rebuild */
+    float dx_max_part;
+
   } stars;
 };
 
@@ -426,6 +433,18 @@ struct cell {
     /*! Nr of #spart in this cell. */
     int count;
 
+    /*! Max smoothing length in this cell. */
+    double h_max;
+
+    /*! Values of h_max before the drifts, used for sub-cell tasks. */
+    float h_max_old;
+
+    /*! Maximum part movement in this cell since last construction. */
+    float dx_max_part;
+
+    /*! Values of dx_max before the drifts, used for sub-cell tasks. */
+    float dx_max_part_old;
+
     /*! Dependency implicit task for the star ghost  (in->ghost->out)*/
     struct task *ghost_in;
 
@@ -662,8 +681,12 @@ cell_can_recurse_in_self_hydro_task(const struct cell *c) {
 __attribute__((always_inline)) INLINE static int
 cell_can_recurse_in_pair_stars_task(const struct cell *c) {
 
-  // LOIC: To implement
-  return 0;
+  /* Is the cell split ? */
+  /* If so, is the cut-off radius plus the max distance the parts have moved */
+  /* smaller than the sub-cell sizes ? */
+  /* Note: We use the _old values as these might have been updated by a drift */
+  return c->split && ((kernel_gamma * c->stars.h_max_old +
+                       c->stars.dx_max_part_old) < 0.5f * c->dmin);
 }
 
 /**
@@ -675,8 +698,8 @@ cell_can_recurse_in_pair_stars_task(const struct cell *c) {
 __attribute__((always_inline)) INLINE static int
 cell_can_recurse_in_self_stars_task(const struct cell *c) {
 
-  // LOIC: To implement
-  return 0;
+  /* Is the cell split and not smaller than the smoothing length? */
+  return c->split && (kernel_gamma * c->stars.h_max_old < 0.5f * c->dmin);
 }
 
 /**
@@ -724,8 +747,13 @@ __attribute__((always_inline)) INLINE static int cell_can_split_self_hydro_task(
 __attribute__((always_inline)) INLINE static int cell_can_split_pair_stars_task(
     const struct cell *c) {
 
-  // LOIC: To implement
-  return 0;
+  /* Is the cell split ? */
+  /* If so, is the cut-off radius with some leeway smaller than */
+  /* the sub-cell sizes ? */
+  /* Note that since tasks are create after a rebuild no need to take */
+  /* into account any part motion (i.e. dx_max == 0 here) */
+  return c->split &&
+         (space_stretch * kernel_gamma * c->stars.h_max < 0.5f * c->dmin);
 }
 
 /**
@@ -737,8 +765,13 @@ __attribute__((always_inline)) INLINE static int cell_can_split_pair_stars_task(
 __attribute__((always_inline)) INLINE static int cell_can_split_self_stars_task(
     const struct cell *c) {
 
-  // LOIC: To implement
-  return 0;
+  /* Is the cell split ? */
+  /* If so, is the cut-off radius with some leeway smaller than */
+  /* the sub-cell sizes ? */
+  /* Note: No need for more checks here as all the sub-pairs and sub-self */
+  /* tasks will be created. So no need to check for h_max */
+  return c->split &&
+         (space_stretch * kernel_gamma * c->stars.h_max < 0.5f * c->dmin);
 }
 
 /**

src/drift.h

@@ -137,6 +137,12 @@ __attribute__((always_inline)) INLINE static void drift_spart(
   sp->x[0] += sp->v[0] * dt_drift;
   sp->x[1] += sp->v[1] * dt_drift;
   sp->x[2] += sp->v[2] * dt_drift;
+
+  /* Compute offsets since last cell construction */
+  for (int k = 0; k < 3; k++) {
+    const float dx = sp->v[k] * dt_drift;
+    sp->x_diff[k] -= dx;
+  }
 }
 
 #endif /* SWIFT_DRIFT_H */

src/engine.c

@@ -3329,14 +3329,11 @@ void engine_make_extra_hydroloop_tasks_mapper(void *map_data, int num_elements,
 }
 
 /**
- * @brief Duplicates the first hydro loop and construct all the
- * dependencies for the stars
+ * @brief Creates all the task dependencies for the stars
  *
- * This is done by looping over all the previously constructed tasks
- * and adding another task involving the same cells but this time
- * corresponding to the second hydro loop over neighbours.
- * With all the relevant tasks for a given cell available, we construct
- * all the dependencies for that cell.
+ * @param map_data The tasks
+ * @param num_elements number of tasks
+ * @param extra_data The #engine
  */
 void engine_link_stars_tasks_mapper(void *map_data, int num_elements,
                                     void *extra_data) {
@@ -3358,8 +3355,9 @@ void engine_link_stars_tasks_mapper(void *map_data, int num_elements,
 
       /* Now, build all the dependencies for the stars */
       engine_make_stars_loops_dependencies(sched, t, t->ci);
-      scheduler_addunlock(sched, t->ci->stars.ghost_out,
-                          t->ci->super->end_force);
+      if (t->ci == t->ci->super)
+        scheduler_addunlock(sched, t->ci->super->stars.ghost_out,
+                            t->ci->super->end_force);
     }
 
     /* Otherwise, pair interaction? */
@@ -5062,6 +5060,7 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
 #endif
 
   if (e->nodeID == 0) scheduler_write_dependencies(&e->sched, e->verbose);
+  if (e->nodeID == 0) scheduler_write_task_level(&e->sched);
 
   /* Run the 0th time-step */
   TIMER_TIC2;
@@ -5143,6 +5142,20 @@ void engine_init_particles(struct engine *e, int flag_entropy_ICs,
     }
   }
 
+  if (s->cells_top != NULL && s->nr_sparts > 0) {
+    for (int i = 0; i < s->nr_cells; i++) {
+      struct cell *c = &s->cells_top[i];
+      if (c->nodeID == engine_rank && c->stars.count > 0) {
+        float spart_h_max = c->stars.parts[0].h;
+        for (int k = 1; k < c->stars.count; k++) {
+          if (c->stars.parts[k].h > spart_h_max)
+            spart_h_max = c->stars.parts[k].h;
+        }
+        c->stars.h_max = max(spart_h_max, c->stars.h_max);
+      }
+    }
+  }
+
   clocks_gettime(&time2);
 
 #ifdef SWIFT_DEBUG_CHECKS

src/map.c

@@ -191,6 +191,13 @@ void map_h_max(struct part *p, struct cell *c, void *data) {
   if (p->h > (*p2)->h) *p2 = p;
 }
 
+void map_stars_h_max(struct spart *p, struct cell *c, void *data) {
+
+  struct spart **p2 = (struct spart **)data;
+
+  if (p->h > (*p2)->h) *p2 = p;
+}
+
 /**
  * @brief Mapping function for neighbour count.
  */

src/map.h

@@ -34,6 +34,7 @@ void map_wcount_min(struct part *p, struct cell *c, void *data);
 void map_wcount_max(struct part *p, struct cell *c, void *data);
 void map_h_min(struct part *p, struct cell *c, void *data);
 void map_h_max(struct part *p, struct cell *c, void *data);
+void map_stars_h_max(struct spart *p, struct cell *c, void *data);
 void map_icount(struct part *p, struct cell *c, void *data);
 void map_dump(struct part *p, struct cell *c, void *data);
 

src/runner.c

@@ -283,7 +283,7 @@ void runner_do_stars_ghost(struct runner *r, struct cell *c, int timer) {
         for (struct cell *finger = c; finger != NULL; finger = finger->parent) {
 
           /* Run through this cell's density interactions. */
-          for (struct link *l = finger->hydro.density; l != NULL; l = l->next) {
+          for (struct link *l = finger->stars.density; l != NULL; l = l->next) {
 
 #ifdef SWIFT_DEBUG_CHECKS
             if (l->t->ti_run < r->e->ti_current)
@@ -2322,6 +2322,8 @@ void runner_do_recv_spart(struct runner *r, struct cell *c, int timer) {
   const size_t nr_sparts = c->stars.count;
   const integertime_t ti_current = r->e->ti_current;
 
+  error("Need to add h_max computation");
+
   TIMER_TIC;
 
   integertime_t ti_gravity_end_min = max_nr_timesteps;

src/runner_doiact_stars.h

@@ -963,7 +963,7 @@ void runner_doself_branch_stars_density(struct runner *r, struct cell *c) {
   if (!cell_is_active_stars(c, e)) return;
 
   /* Did we mess up the recursion? */
-  if (c->hydro.h_max_old * kernel_gamma > c->dmin)
+  if (c->stars.h_max_old * kernel_gamma > c->dmin)
     error("Cell smaller than smoothing length");
 
   runner_doself_stars_density(r, c, 1);

src/scheduler.c

@@ -876,8 +876,6 @@ static void scheduler_splittask_hydro(struct task *t, struct scheduler *s) {
  */
 static void scheduler_splittask_stars(struct task *t, struct scheduler *s) {
 
-  // LOIC: This is un-tested. Need to verify that it works.
-
   /* Iterate on this task until we're done with it. */
   int redo = 1;
   while (redo) {
@@ -1829,6 +1827,7 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
     const float gcount_i = (t->ci != NULL) ? t->ci->grav.count : 0.f;
     const float gcount_j = (t->cj != NULL) ? t->cj->grav.count : 0.f;
     const float scount_i = (t->ci != NULL) ? t->ci->stars.count : 0.f;
+    const float scount_j = (t->cj != NULL) ? t->cj->stars.count : 0.f;
 
     switch (t->type) {
       case task_type_sort:
@@ -1837,22 +1836,30 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
         break;
 
       case task_type_self:
-        // LOIC: Need to do something for stars here
         if (t->subtype == task_subtype_grav)
           cost = 1.f * (wscale * gcount_i) * gcount_i;
         else if (t->subtype == task_subtype_external_grav)
           cost = 1.f * wscale * gcount_i;
+        else if (t->subtype == task_subtype_stars_density)
+          cost = 1.f * wscale * scount_i * count_i;
         else
           cost = 1.f * (wscale * count_i) * count_i;
         break;
 
       case task_type_pair:
-        // LOIC: Need to do something for stars here
         if (t->subtype == task_subtype_grav) {
           if (t->ci->nodeID != nodeID || t->cj->nodeID != nodeID)
             cost = 3.f * (wscale * gcount_i) * gcount_j;
           else
             cost = 2.f * (wscale * gcount_i) * gcount_j;
+        } else if (t->subtype == task_subtype_stars_density) {
+          if (t->ci->nodeID != nodeID)
+            cost = 3.f * wscale * count_i * scount_j * sid_scale[t->flags];
+          else if (t->cj->nodeID != nodeID)
+            cost = 3.f * wscale * scount_i * count_j * sid_scale[t->flags];
+          else
+            cost = 2.f * wscale * (scount_i * count_j + scount_j * count_i) *
+                   sid_scale[t->flags];
         } else {
           if (t->ci->nodeID != nodeID || t->cj->nodeID != nodeID)
             cost = 3.f * (wscale * count_i) * count_j * sid_scale[t->flags];
@@ -1862,23 +1869,34 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
         break;
 
       case task_type_sub_pair:
-        // LOIC: Need to do something for stars here
-        if (t->ci->nodeID != nodeID || t->cj->nodeID != nodeID) {
 #ifdef SWIFT_DEBUG_CHECKS
-          if (t->flags < 0) error("Negative flag value!");
+        if (t->flags < 0) error("Negative flag value!");
 #endif
-          cost = 3.f * (wscale * count_i) * count_j * sid_scale[t->flags];
+        if (t->subtype == task_subtype_stars_density) {
+          if (t->ci->nodeID != nodeID) {
+            cost = 3.f * (wscale * count_i) * scount_j * sid_scale[t->flags];
+          } else if (t->cj->nodeID != nodeID) {
+            cost = 3.f * (wscale * scount_i) * count_j * sid_scale[t->flags];
+          } else {
+            cost = 2.f * wscale * (scount_i * count_j + scount_j * count_i) *
+                   sid_scale[t->flags];
+          }