diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_sink/getIC.sh b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/getIC.sh
index a3d16db27aac06abda683a7bd75e72a275f8b9d4..f354eb1c253039e02b409c9562c7c0b0716b07f6 100755
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_sink/getIC.sh
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/getIC.sh
@@ -1,2 +1,3 @@
#!/bin/bash
-wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/3e11-star-only-DM-halo-galaxy.hdf5
+wget http://virgodb.cosma.dur.ac.uk/swift-webstorage/ICs/IsolatedGalaxies/$1
+
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_sink/isolated_galaxy.yml b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/isolated_galaxy.yml
index c83abbd55d207095655cca733fd77f596a5b6c15..7428181f3f17b5e90b216f597149fcbca3187cc7 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_sink/isolated_galaxy.yml
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/isolated_galaxy.yml
@@ -8,6 +8,9 @@ InternalUnitSystem:
# Parameters for the self-gravity scheme
Gravity:
+ MAC: adaptive # Choice of mulitpole acceptance criterion: 'adaptive' OR 'geometric'.
+ epsilon_fmm: 0.001 # Tolerance parameter for the adaptive multipole acceptance criterion.
+ theta_cr: 0.7 # Opening angle for the purely gemoetric criterion.
eta: 0.025 # Constant dimensionless multiplier for time integration.
theta: 0.7 # Opening angle (Multipole acceptance criterion).
max_physical_DM_softening: 0.7 # Physical softening length (in internal units).
@@ -39,4 +42,41 @@ InitialConditions:
file_name: isolated_galaxy.hdf5 # The file to read
periodic: 0 # Are we running with periodic ICs?
+# 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.
+ use_mass_weighted_num_ngb: 0 # (Optional) Are we using the mass-weighted definition of the number of neighbours?
+ h_tolerance: 1e-4 # (Optional) Relative accuracy of the Netwon-Raphson scheme for the smoothing lengths.
+ h_max: 10. # (Optional) Maximal allowed smoothing length in internal units. Defaults to FLT_MAX if unspecified.
+ h_min_ratio: 0. # (Optional) Minimal allowed smoothing length in units of the softening. Defaults to 0 if unspecified.
+ max_volume_change: 1.4 # (Optional) Maximal allowed change of kernel volume over one time-step.
+ max_ghost_iterations: 30 # (Optional) Maximal number of iterations allowed to converge towards the smoothing length.
+ particle_splitting: 1 # (Optional) Are we splitting particles that are too massive (default: 0)
+ particle_splitting_mass_threshold: 7e-4 # (Optional) Mass threshold for particle splitting (in internal units)
+ generate_random_ids: 0 # (Optional) When creating new particles via splitting, generate ids at random (1) or use new IDs beyond the current range (0) (default: 0)
+ initial_temperature: 0 # (Optional) Initial temperature (in internal units) to set the gas particles at start-up. Value is ignored if set to 0.
+ minimal_temperature: 0 # (Optional) Minimal temperature (in internal units) allowed for the gas particles. Value is ignored if set to 0.
+ H_mass_fraction: 0.755 # (Optional) Hydrogen mass fraction used for initial conversion from temp to internal energy. Default value is derived from the physical constants.
+ H_ionization_temperature: 1e4 # (Optional) Temperature of the transition from neutral to ionized Hydrogen for primoridal gas.
+ viscosity_alpha: 0.8 # (Optional) Override for the initial value of the artificial viscosity. In schemes that have a fixed AV, this remains as alpha throughout the run.
+ viscosity_alpha_max: 2.0 # (Optional) Maximal value for the artificial viscosity in schemes that allow alpha to vary.
+ viscosity_alpha_min: 0.1 # (Optional) Minimal value for the artificial viscosity in schemes that allow alpha to vary.
+ viscosity_length: 0.1 # (Optional) Decay length for the artificial viscosity in schemes that allow alpha to vary.
+ diffusion_alpha: 0.0 # (Optional) Override the initial value for the thermal diffusion coefficient in schemes with thermal diffusion.
+ diffusion_beta: 0.01 # (Optional) Override the decay/rise rate tuning parameter for the thermal diffusion.
+ diffusion_alpha_max: 1.0 # (Optional) Override the maximal thermal diffusion coefficient that is allowed for a given particle.
+ diffusion_alpha_min: 0.0 # (Optional) Override the minimal thermal diffusion coefficient that is allowed for a given particle.
+# Hernquist potential parameters
+HernquistPotential:
+ useabspos: 0 # 0 -> positions based on centre, 1 -> absolute positions
+ position: [0.,0.,0.] # Location of centre of isothermal potential with respect to centre of the box (if 0) otherwise absolute (if 1) (internal units)
+ idealizeddisk: 1 # Run with an idealized galaxy disk
+ M200: 137.0 # M200 of the galaxy disk
+ h: 0.704 # reduced Hubble constant (value does not specify the used units!)
+ concentration: 9.0 # concentration of the Halo
+ diskfraction: 0.040 # Disk mass fraction
+ bulgefraction: 0.014 # Bulge mass fraction
+ timestep_mult: 0.01 # Dimensionless pre-factor for the time-step condition, basically determines the fraction of the orbital time we use to do the time integration
+ epsilon: 0.2 # Softening size (internal units)
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_sink/makeIC.py b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/makeIC.py
index 3587c9d858ab1a72a39b1b243b6d58bae7fad383..61744122a760bb3afef1f2b61db816bdff766bfb 100644
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_sink/makeIC.py
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/makeIC.py
@@ -20,14 +20,14 @@
import h5py
import numpy as np
from shutil import copyfile
+import sys
# Add sink particles to the isolated galaxy
-fileName = "3e11-star-only-DM-halo-galaxy.hdf5"
+fileName = sys.argv[-1]
output = "isolated_galaxy.hdf5"
-L = 13756 # Size of the box (internal units)
-L_sink = 1000 # Size of the sink particle area (L_sink < L)
+L_sink = 50 # Size of the sink particle area
N = 1000 # Number of sink particles
max_vel = 100 # Maximal velocity of the sink particles (in km / s)
mass = 0.000142 # Mass of a sink particle (internal units)
@@ -40,6 +40,11 @@ copyfile(fileName, output)
# File
file = h5py.File(output, 'a')
+# Get the box size (suppose that the galaxy is at the center of the box)
+L = 2 * file["PartType0/Coordinates"][:].mean()
+if (L < L_sink):
+ raise Exception("Error the size of the sink box is too large")
+
pos = np.random.rand(N, 3) * L_sink + 0.5 * (L - L_sink)
vel = 2 * (np .random.rand(N, 3) - 0.5) * max_vel
m = mass * np.ones(N)
diff --git a/examples/IsolatedGalaxy/IsolatedGalaxy_sink/run.sh b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/run.sh
index 7d2cee8c29c0eadc886b3188d532567ea861bde5..8f3e976956e4b79a7d5a708eb65ed76284b19b7d 100755
--- a/examples/IsolatedGalaxy/IsolatedGalaxy_sink/run.sh
+++ b/examples/IsolatedGalaxy/IsolatedGalaxy_sink/run.sh
@@ -1,12 +1,20 @@
#!/bin/bash
-if [ ! -e 3e11-star-only-DM-halo-galaxy.hdf5 ]
+# filename=fid.hdf5
+# filename=f10.hdf5
+# filename=f90.hdf5
+filename=lowres8.hdf5
+# filename=lowres64.hdf5
+# filename=lowres512.hdf5
+# filename=highres6.hdf5
+
+if [ ! -e $filename ]
then
echo "Fetching initial conditons for the isolated galaxy with an external potential ..."
- ./getIC.sh
+ ./getIC.sh $filename
fi
echo "Generate initial conditions"
-python3 makeIC.py
+python3 makeIC.py $filename
-../../swift --sinks --external-gravity --self-gravity --threads=16 isolated_galaxy.yml 2>&1 | tee output.log
+../../swift --hydro --sinks --external-gravity --self-gravity --threads=1 isolated_galaxy.yml 2>&1 | tee output.log
diff --git a/examples/main.c b/examples/main.c
index 3a3bee4e51a021310d8a08b5ff92ca77e3d597c1..f8b645111310314575bdf4872ff81bbd11ca2267 100644
--- a/examples/main.c
+++ b/examples/main.c
@@ -1298,7 +1298,7 @@ int main(int argc, char *argv[]) {
if (with_fof) engine_policies |= engine_policy_fof;
if (with_logger) engine_policies |= engine_policy_logger;
if (with_line_of_sight) engine_policies |= engine_policy_line_of_sight;
- if (with_sink) engine_policies |= engine_policy_sink;
+ if (with_sink) engine_policies |= engine_policy_sinks;
if (with_rt) engine_policies |= engine_policy_rt;
/* Initialize the engine with the space and policies. */
@@ -1401,10 +1401,12 @@ int main(int argc, char *argv[]) {
/* Legend */
if (myrank == 0) {
- printf("# %6s %14s %12s %12s %14s %9s %12s %12s %12s %12s %16s [%s] %6s\n",
- "Step", "Time", "Scale-factor", "Redshift", "Time-step", "Time-bins",
- "Updates", "g-Updates", "s-Updates", "b-Updates", "Wall-clock time",
- clocks_getunit(), "Props");
+ printf(
+ "# %6s %14s %12s %12s %14s %9s %12s %12s %12s %12s %12s %16s [%s] "
+ "%6s\n",
+ "Step", "Time", "Scale-factor", "Redshift", "Time-step", "Time-bins",
+ "Updates", "g-Updates", "s-Updates", "sink-Updates", "b-Updates",
+ "Wall-clock time", clocks_getunit(), "Props");
fflush(stdout);
}
@@ -1548,19 +1550,21 @@ int main(int argc, char *argv[]) {
/* Print some information to the screen */
printf(
- " %6d %14e %12.7f %12.7f %14e %4d %4d %12lld %12lld %12lld %12lld"
+ " %6d %14e %12.7f %12.7f %14e %4d %4d %12lld %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.sink_updates, e.b_updates, e.wallclock_time, e.step_props);
fflush(stdout);
fprintf(e.file_timesteps,
" %6d %14e %12.7f %12.7f %14e %4d %4d %12lld %12lld %12lld %12lld"
- " %21.3f %6d\n",
+ " %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.s_updates, e.sink_updates, e.b_updates, e.wallclock_time,
+ e.step_props);
fflush(e.file_timesteps);
/* Print information to the SFH logger */
diff --git a/src/cell.c b/src/cell.c
index 30fd9031e159b2dab033623fc2adef34cf2d5bbf..2a947f6b6ed32080a19b08a87b71fdf5a12bb8b9 100644
--- a/src/cell.c
+++ b/src/cell.c
@@ -2065,6 +2065,7 @@ void cell_split(struct cell *c, ptrdiff_t parts_offset, ptrdiff_t sparts_offset,
/* Store the counts and offsets. */
for (int k = 0; k < 8; k++) {
c->progeny[k]->sinks.count = bucket_count[k];
+ c->progeny[k]->sinks.count_total = c->progeny[k]->sinks.count;
c->progeny[k]->sinks.parts = &c->sinks.parts[bucket_offset[k]];
}
@@ -2294,6 +2295,44 @@ void cell_check_gpart_drift_point(struct cell *c, void *data) {
#endif
}
+/**
+ * @brief Checks that the #sink in a cell are at the
+ * current point in time
+ *
+ * Calls error() if the cell is not at the current time.
+ *
+ * @param c Cell to act upon
+ * @param data The current time on the integer time-line
+ */
+void cell_check_sink_drift_point(struct cell *c, void *data) {
+#ifdef SWIFT_DEBUG_CHECKS
+
+ const integertime_t ti_drift = *(integertime_t *)data;
+
+ /* Only check local cells */
+ if (c->nodeID != engine_rank) return;
+
+ /* Only check cells with content */
+ if (c->sinks.count == 0) return;
+
+ if (c->sinks.ti_old_part != ti_drift)
+ error(
+ "Cell in an incorrect time-zone! c->sinks.ti_old_part=%lld "
+ "ti_drift=%lld",
+ c->sinks.ti_old_part, ti_drift);
+
+ for (int i = 0; i < c->sinks.count; ++i)
+ if (c->sinks.parts[i].ti_drift != ti_drift &&
+ c->sinks.parts[i].time_bin != time_bin_inhibited)
+ error(
+ "sink-part in an incorrect time-zone! sink->ti_drift=%lld "
+ "ti_drift=%lld",
+ c->sinks.parts[i].ti_drift, ti_drift);
+#else
+ error("Calling debugging code without debugging flag activated.");
+#endif
+}
+
/**
* @brief Checks that the #spart in a cell are at the
* current point in time
@@ -2633,7 +2672,8 @@ void cell_clear_drift_flags(struct cell *c, void *data) {
cell_flag_do_grav_drift | cell_flag_do_grav_sub_drift |
cell_flag_do_bh_drift | cell_flag_do_bh_sub_drift |
cell_flag_do_stars_drift |
- cell_flag_do_stars_sub_drift);
+ cell_flag_do_stars_sub_drift |
+ cell_flag_do_sink_drift | cell_flag_do_sink_sub_drift);
}
/**
@@ -3027,6 +3067,44 @@ void cell_activate_drift_bpart(struct cell *c, struct scheduler *s) {
}
}
+/**
+ * @brief Activate the #sink drifts on the given cell.
+ */
+void cell_activate_drift_sink(struct cell *c, struct scheduler *s) {
+
+ /* If this cell is already marked for drift, quit early. */
+ if (cell_get_flag(c, cell_flag_do_sink_drift)) return;
+
+ /* Mark this cell for drifting. */
+ cell_set_flag(c, cell_flag_do_sink_drift);
+
+ /* Set the do_sink_sub_drifts all the way up and activate the super
+ drift if this has not yet been done. */
+ if (c == c->hydro.super) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (c->sinks.drift == NULL)
+ error("Trying to activate un-existing c->sinks.drift");
+#endif
+ scheduler_activate(s, c->sinks.drift);
+ } else {
+ for (struct cell *parent = c->parent;
+ parent != NULL && !cell_get_flag(parent, cell_flag_do_sink_sub_drift);
+ parent = parent->parent) {
+ /* Mark this cell for drifting */
+ cell_set_flag(parent, cell_flag_do_sink_sub_drift);
+
+ if (parent == c->hydro.super) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (parent->sinks.drift == NULL)
+ error("Trying to activate un-existing parent->sinks.drift");
+#endif
+ scheduler_activate(s, parent->sinks.drift);
+ break;
+ }
+ }
+ }
+}
+
/**
* @brief Activate the drifts on the given cell.
*/
@@ -3514,6 +3592,102 @@ void cell_activate_subcell_black_holes_tasks(struct cell *ci, struct cell *cj,
} /* Otherwise, pair interation */
}
+/**
+ * @brief Traverse a sub-cell task and activate the sinks drift tasks that
+ * are required by a sinks task
+ *
+ * @param ci The first #cell we recurse in.
+ * @param cj The second #cell we recurse in.
+ * @param s The task #scheduler.
+ * @param with_timestep_sync Are we running with time-step synchronization on?
+ */
+void cell_activate_subcell_sinks_tasks(struct cell *ci, struct cell *cj,
+ struct scheduler *s,
+ const int with_timestep_sync) {
+ const struct engine *e = s->space->e;
+
+ /* Store the current dx_max and h_max values. */
+ ci->sinks.dx_max_part_old = ci->sinks.dx_max_part;
+ ci->hydro.dx_max_part_old = ci->hydro.dx_max_part;
+ ci->hydro.h_max_old = ci->hydro.h_max;
+
+ if (cj != NULL) {
+ cj->sinks.dx_max_part_old = cj->sinks.dx_max_part;
+ cj->hydro.dx_max_part_old = cj->hydro.dx_max_part;
+ cj->hydro.h_max_old = cj->hydro.h_max;
+ }
+
+ /* Self interaction? */
+ if (cj == NULL) {
+ /* Do anything? */
+ if (!cell_is_active_sinks(ci, e) || ci->hydro.count == 0 ||
+ ci->sinks.count == 0)
+ return;
+
+ /* Recurse? */
+ if (cell_can_recurse_in_self_sinks_task(ci)) {
+ /* Loop over all progenies and pairs of progenies */
+ for (int j = 0; j < 8; j++) {
+ if (ci->progeny[j] != NULL) {
+ cell_activate_subcell_sinks_tasks(ci->progeny[j], NULL, s,
+ with_timestep_sync);
+ for (int k = j + 1; k < 8; k++)
+ if (ci->progeny[k] != NULL)
+ cell_activate_subcell_sinks_tasks(ci->progeny[j], ci->progeny[k],
+ s, with_timestep_sync);
+ }
+ }
+ } else {
+ /* We have reached the bottom of the tree: activate drift */
+ cell_activate_drift_sink(ci, s);
+ cell_activate_drift_part(ci, s);
+ }
+ }
+
+ /* Otherwise, pair interation */
+ else {
+ error("Not implemented yet.");
+ /* /\* Should we even bother? *\/ */
+ /* if (!cell_is_active_sinks(ci, e) && */
+ /* !cell_is_active_sinks(cj, e)) */
+ /* return; */
+
+ /* /\* Get the orientation of the pair. *\/ */
+ /* double shift[3]; */
+ /* const int sid = space_getsid(s->space, &ci, &cj, shift); */
+
+ /* /\* recurse? *\/ */
+ /* if (cell_can_recurse_in_pair_sinks_task(ci, cj) && */
+ /* cell_can_recurse_in_pair_sinks_task(cj, ci)) { */
+ /* const struct cell_split_pair *csp = &cell_split_pairs[sid]; */
+ /* for (int k = 0; k < csp->count; k++) { */
+ /* const int pid = csp->pairs[k].pid; */
+ /* const int pjd = csp->pairs[k].pjd; */
+ /* if (ci->progeny[pid] != NULL && cj->progeny[pjd] != NULL) */
+ /* cell_activate_subcell_sinks_tasks( */
+ /* ci->progeny[pid], cj->progeny[pjd], s, with_timestep_sync); */
+ /* } */
+ /* } */
+
+ /* /\* Otherwise, activate the drifts. *\/ */
+ /* else if (cell_is_active_sinks(ci, e) || */
+ /* cell_is_active_sinks(cj, e)) { */
+
+ /* /\* Activate the drifts if the cells are local. *\/ */
+ /* if (ci->nodeID == engine_rank) cell_activate_drift_sinks(ci, s); */
+ /* if (cj->nodeID == engine_rank) cell_activate_drift_part(cj, s); */
+ /* if (cj->nodeID == engine_rank && with_timestep_sync) */
+ /* cell_activate_sync_part(cj, s); */
+
+ /* /\* Activate the drifts if the cells are local. *\/ */
+ /* if (ci->nodeID == engine_rank) cell_activate_drift_part(ci, s); */
+ /* if (cj->nodeID == engine_rank) cell_activate_drift_sinks(cj, s); */
+ /* if (ci->nodeID == engine_rank && with_timestep_sync) */
+ /* cell_activate_sync_part(ci, s); */
+ /* } */
+ } /* Otherwise, pair interation */
+}
+
/**
* @brief Traverse a sub-cell task and activate the gravity drift tasks that
* are required by a self gravity task.
@@ -4609,6 +4783,38 @@ int cell_unskip_black_holes_tasks(struct cell *c, struct scheduler *s) {
return rebuild;
}
+/**
+ * @brief Un-skips all the sinks tasks associated with a given cell and
+ * checks if the space needs to be rebuilt.
+ *
+ * @param c the #cell.
+ * @param s the #scheduler.
+ *
+ * @return 1 If the space needs rebuilding. 0 otherwise.
+ */
+int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s) {
+
+ struct engine *e = s->space->e;
+ const int nodeID = e->nodeID;
+ int rebuild = 0;
+
+ if (c->sinks.drift != NULL && cell_is_active_sinks(c, e)) {
+ cell_activate_drift_sink(c, s);
+ }
+
+ /* Unskip all the other task types. */
+ if (c->nodeID == nodeID && cell_is_active_sinks(c, e)) {
+ if (c->kick1 != NULL) scheduler_activate(s, c->kick1);
+ if (c->kick2 != NULL) scheduler_activate(s, c->kick2);
+ if (c->timestep != NULL) scheduler_activate(s, c->timestep);
+#ifdef WITH_LOGGER
+ if (c->logger != NULL) scheduler_activate(s, c->logger);
+#endif
+ }
+
+ return rebuild;
+}
+
/**
* @brief Set the super-cell pointers for all cells in a hierarchy.
*
@@ -5439,6 +5645,156 @@ void cell_drift_bpart(struct cell *c, const struct engine *e, int force) {
cell_clear_flag(c, cell_flag_do_bh_drift | cell_flag_do_bh_sub_drift);
}
+/**
+ * @brief Recursively drifts the #sinks in a cell hierarchy.
+ *
+ * @param c The #cell.
+ * @param e The #engine (to get ti_current).
+ * @param force Drift the particles irrespective of the #cell flags.
+ */
+void cell_drift_sink(struct cell *c, const struct engine *e, int force) {
+
+ const int periodic = e->s->periodic;
+ const double dim[3] = {e->s->dim[0], e->s->dim[1], e->s->dim[2]};
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const integertime_t ti_old_bpart = c->sinks.ti_old_part;
+ const integertime_t ti_current = e->ti_current;
+ struct sink *const sinks = c->sinks.parts;
+
+ float dx_max = 0.f, dx2_max = 0.f;
+
+ /* Drift irrespective of cell flags? */
+ force = (force || cell_get_flag(c, cell_flag_do_sink_drift));
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that we only drift local cells. */
+ if (c->nodeID != engine_rank) error("Drifting a foreign cell is nope.");
+
+ /* Check that we are actually going to move forward. */
+ if (ti_current < ti_old_bpart) error("Attempt to drift to the past");
+#endif
+
+ /* Early abort? */
+ if (c->sinks.count == 0) {
+
+ /* Clear the drift flags. */
+ cell_clear_flag(c, cell_flag_do_sink_drift | cell_flag_do_sink_sub_drift);
+
+ /* Update the time of the last drift */
+ c->sinks.ti_old_part = ti_current;
+
+ return;
+ }
+
+ /* Ok, we have some particles somewhere in the hierarchy to drift */
+
+ /* Are we not in a leaf ? */
+ if (c->split && (force || cell_get_flag(c, cell_flag_do_sink_sub_drift))) {
+
+ /* Loop over the progeny and collect their data. */
+ for (int k = 0; k < 8; k++) {
+ if (c->progeny[k] != NULL) {
+ struct cell *cp = c->progeny[k];
+
+ /* Recurse */
+ cell_drift_sink(cp, e, force);
+
+ /* Update */
+ dx_max = max(dx_max, cp->sinks.dx_max_part);
+ }
+ }
+
+ /* Store the values */
+ c->sinks.dx_max_part = dx_max;
+
+ /* Update the time of the last drift */
+ c->sinks.ti_old_part = ti_current;
+
+ } else if (!c->split && force && ti_current > ti_old_bpart) {
+
+ /* Drift from the last time the cell was drifted to the current time */
+ double dt_drift;
+ if (with_cosmology) {
+ dt_drift =
+ cosmology_get_drift_factor(e->cosmology, ti_old_bpart, ti_current);
+ } else {
+ dt_drift = (ti_current - ti_old_bpart) * e->time_base;
+ }
+
+ /* Loop over all the star particles in the cell */
+ const size_t nr_sinks = c->sinks.count;
+ for (size_t k = 0; k < nr_sinks; k++) {
+
+ /* Get a handle on the sink. */
+ struct sink *const sink = &sinks[k];
+
+ /* Ignore inhibited particles */
+ if (sink_is_inhibited(sink, e)) continue;
+
+ /* Drift... */
+ drift_sink(sink, dt_drift, ti_old_bpart, ti_current);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Make sure the particle does not drift by more than a box length. */
+ if (fabs(sink->v[0] * dt_drift) > e->s->dim[0] ||
+ fabs(sink->v[1] * dt_drift) > e->s->dim[1] ||
+ fabs(sink->v[2] * dt_drift) > e->s->dim[2]) {
+ error("Particle drifts by more than a box length!");
+ }
+#endif
+
+ /* In non-periodic BC runs, remove particles that crossed the border */
+ if (!periodic) {
+
+ /* Did the particle leave the box? */
+ if ((sink->x[0] > dim[0]) || (sink->x[0] < 0.) || // x
+ (sink->x[1] > dim[1]) || (sink->x[1] < 0.) || // y
+ (sink->x[2] > dim[2]) || (sink->x[2] < 0.)) { // z
+
+ lock_lock(&e->s->lock);
+
+ /* Re-check that the particle has not been removed
+ * by another thread before we do the deed. */
+ if (!sink_is_inhibited(sink, e)) {
+
+ /* Remove the particle entirely */
+ // cell_remove_sink(e, c, bp);
+ error("TODO: loic implement cell_remove_sink");
+ }
+
+ if (lock_unlock(&e->s->lock) != 0)
+ error("Failed to unlock the space!");
+
+ continue;
+ }
+ }
+
+ /* Compute (square of) motion since last cell construction */
+ const float dx2 = sink->x_diff[0] * sink->x_diff[0] +
+ sink->x_diff[1] * sink->x_diff[1] +
+ sink->x_diff[2] * sink->x_diff[2];
+ dx2_max = max(dx2_max, dx2);
+
+ /* Get ready for a density calculation */
+ if (sink_is_active(sink, e)) {
+ sink_init_sink(sink);
+ }
+ }
+
+ /* Now, get the maximal particle motion from its square */
+ dx_max = sqrtf(dx2_max);
+
+ /* Store the values */
+ c->sinks.dx_max_part = dx_max;
+
+ /* Update the time of the last drift */
+ c->sinks.ti_old_part = ti_current;
+ }
+
+ /* Clear the drift flags. */
+ cell_clear_flag(c, cell_flag_do_sink_drift | cell_flag_do_sink_sub_drift);
+}
+
/**
* @brief Recursively drifts all multipoles in a cell hierarchy.
*
@@ -5587,7 +5943,7 @@ void cell_check_timesteps(const struct cell *c, const integertime_t ti_current,
if (c->hydro.ti_end_min == 0 && c->grav.ti_end_min == 0 &&
c->stars.ti_end_min == 0 && c->black_holes.ti_end_min == 0 &&
- c->nr_tasks > 0)
+ c->sinks.ti_end_min == 0 && c->nr_tasks > 0)
error("Cell without assigned time-step");
if (c->split) {
diff --git a/src/cell.h b/src/cell.h
index c8c0d664a0798d310e8579a74a4000a0b54a18e9..55e2c89387ce42078f42c753c516da15542b8b96 100644
--- a/src/cell.h
+++ b/src/cell.h
@@ -304,10 +304,12 @@ enum cell_flags {
cell_flag_do_stars_sub_drift = (1UL << 10),
cell_flag_do_bh_drift = (1UL << 11),
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)
+ cell_flag_do_sink_drift = (1UL << 13),
+ cell_flag_do_sink_sub_drift = (1UL << 14),
+ cell_flag_do_stars_resort = (1UL << 15),
+ cell_flag_has_tasks = (1UL << 16),
+ cell_flag_do_hydro_sync = (1UL << 17),
+ cell_flag_do_hydro_sub_sync = (1UL << 18)
};
/**
@@ -317,7 +319,7 @@ enum cell_flags {
*/
struct cell {
- /*! The cell location on the grid. */
+ /*! The cell location on the grid (corner nearest to the origin). */
double loc[3];
/*! The cell dimensions. */
@@ -764,12 +766,21 @@ struct cell {
/*! Nr of #sink this cell can hold after addition of new one. */
int count_total;
+ /*! Number of #sink updated in this cell. */
+ int updated;
+
/*! Is the #sink data of this cell being used in a sub-cell? */
int hold;
/*! Spin lock for various uses (#sink case). */
swift_lock_type lock;
+ /*! 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;
+
/*! Last (integer) time the cell's sink were drifted forward in time. */
integertime_t ti_old_part;
@@ -783,6 +794,10 @@ struct cell {
* tasks.
*/
integertime_t ti_beg_max;
+
+ /*! The drift task for sinks */
+ struct task *drift;
+
} sinks;
#ifdef WITH_MPI
@@ -940,16 +955,19 @@ void cell_clean(struct cell *c);
void cell_check_part_drift_point(struct cell *c, void *data);
void cell_check_gpart_drift_point(struct cell *c, void *data);
void cell_check_spart_drift_point(struct cell *c, void *data);
+void cell_check_sink_drift_point(struct cell *c, void *data);
void cell_check_multipole_drift_point(struct cell *c, void *data);
void cell_reset_task_counters(struct cell *c);
int cell_unskip_hydro_tasks(struct cell *c, struct scheduler *s);
int cell_unskip_stars_tasks(struct cell *c, struct scheduler *s,
const int with_star_formation);
+int cell_unskip_sinks_tasks(struct cell *c, struct scheduler *s);
int cell_unskip_black_holes_tasks(struct cell *c, struct scheduler *s);
int cell_unskip_gravity_tasks(struct cell *c, struct scheduler *s);
void cell_drift_part(struct cell *c, const struct engine *e, int force);
void cell_drift_gpart(struct cell *c, const struct engine *e, int force);
void cell_drift_spart(struct cell *c, const struct engine *e, int force);
+void cell_drift_sink(struct cell *c, const struct engine *e, int force);
void cell_drift_bpart(struct cell *c, const struct engine *e, int force);
void cell_drift_multipole(struct cell *c, const struct engine *e);
void cell_drift_all_multipoles(struct cell *c, const struct engine *e);
@@ -977,6 +995,7 @@ void cell_activate_super_spart_drifts(struct cell *c, struct scheduler *s);
void cell_activate_drift_part(struct cell *c, struct scheduler *s);
void cell_activate_drift_gpart(struct cell *c, struct scheduler *s);
void cell_activate_drift_spart(struct cell *c, struct scheduler *s);
+void cell_activate_drift_sink(struct cell *c, struct scheduler *s);
void cell_activate_drift_bpart(struct cell *c, struct scheduler *s);
void cell_activate_sync_part(struct cell *c, struct scheduler *s);
void cell_activate_hydro_sorts(struct cell *c, int sid, struct scheduler *s);
@@ -1185,6 +1204,22 @@ cell_can_recurse_in_self_black_holes_task(const struct cell *c) {
(kernel_gamma * c->hydro.h_max_old < 0.5f * c->dmin);
}
+/**
+ * @brief Can a sub-self sinks task recurse to a lower level based
+ * on the status of the particles in the cell.
+ *
+ * @param c The #cell.
+ */
+__attribute__((always_inline)) INLINE static int
+cell_can_recurse_in_self_sinks_task(const struct cell *c) {
+
+ /* Is the cell split and not smaller than the smoothing length? */
+ // loic TODO: add cut off radius
+ return c->split &&
+ //(kernel_gamma * c->sinks.h_max_old < 0.5f * c->dmin) &&
+ (kernel_gamma * c->hydro.h_max_old < 0.5f * c->dmin);
+}
+
/**
* @brief Can a pair hydro task associated with a cell be split into smaller
* sub-tasks.
diff --git a/src/collectgroup.c b/src/collectgroup.c
index 6aa120a5708df486229773c8ff9dde4f24180ccf..1c75d8c152472ff52b7e706cf6fd50469e0e80b1 100644
--- a/src/collectgroup.c
+++ b/src/collectgroup.c
@@ -37,19 +37,22 @@
/* Local collections for MPI reduces. */
struct mpicollectgroup1 {
- long long updated, g_updated, s_updated, b_updated;
- long long inhibited, g_inhibited, s_inhibited, b_inhibited;
+ long long updated, g_updated, s_updated, sink_updated, b_updated;
+ long long inhibited, g_inhibited, s_inhibited, sink_inhibited, b_inhibited;
integertime_t ti_hydro_end_min;
integertime_t ti_gravity_end_min;
integertime_t ti_stars_end_min;
+ integertime_t ti_sinks_end_min;
integertime_t ti_black_holes_end_min;
integertime_t ti_hydro_end_max;
integertime_t ti_gravity_end_max;
integertime_t ti_stars_end_max;
+ integertime_t ti_sinks_end_max;
integertime_t ti_black_holes_end_max;
integertime_t ti_hydro_beg_max;
integertime_t ti_gravity_beg_max;
integertime_t ti_stars_beg_max;
+ integertime_t ti_sinks_beg_max;
integertime_t ti_black_holes_beg_max;
int forcerebuild;
long long total_nr_cells;
@@ -106,20 +109,30 @@ void collectgroup1_apply(const struct collectgroup1 *grp1, struct engine *e) {
e->ti_black_holes_end_min = grp1->ti_black_holes_end_min;
e->ti_black_holes_end_max = grp1->ti_black_holes_end_max;
e->ti_black_holes_beg_max = grp1->ti_black_holes_beg_max;
- e->ti_end_min = min4(e->ti_hydro_end_min, e->ti_gravity_end_min,
- e->ti_stars_end_min, e->ti_black_holes_end_min);
- e->ti_end_max = max4(e->ti_hydro_end_max, e->ti_gravity_end_max,
- e->ti_stars_end_max, e->ti_black_holes_end_max);
- e->ti_beg_max = max4(e->ti_hydro_beg_max, e->ti_gravity_beg_max,
- e->ti_stars_beg_max, e->ti_black_holes_beg_max);
+ e->ti_sinks_end_min = grp1->ti_sinks_end_min;
+ e->ti_sinks_end_max = grp1->ti_sinks_end_max;
+ e->ti_sinks_beg_max = grp1->ti_sinks_beg_max;
+
+ e->ti_end_min =
+ min5(e->ti_hydro_end_min, e->ti_gravity_end_min, e->ti_sinks_end_min,
+ e->ti_stars_end_min, e->ti_black_holes_end_min);
+ e->ti_end_max =
+ max5(e->ti_hydro_end_max, e->ti_gravity_end_max, e->ti_sinks_end_max,
+ e->ti_stars_end_max, e->ti_black_holes_end_max);
+ e->ti_beg_max =
+ max5(e->ti_hydro_beg_max, e->ti_gravity_beg_max, e->ti_sinks_beg_max,
+ e->ti_stars_beg_max, e->ti_black_holes_beg_max);
+
e->updates = grp1->updated;
e->g_updates = grp1->g_updated;
e->s_updates = grp1->s_updated;
e->b_updates = grp1->b_updated;
+ e->sink_updates = grp1->sink_updated;
e->nr_inhibited_parts = grp1->inhibited;
e->nr_inhibited_gparts = grp1->g_inhibited;
e->nr_inhibited_sparts = grp1->s_inhibited;
e->nr_inhibited_bparts = grp1->b_inhibited;
+ e->nr_inhibited_sinks = grp1->sink_inhibited;
e->forcerebuild = grp1->forcerebuild;
e->total_nr_cells = grp1->total_nr_cells;
e->total_nr_tasks = grp1->total_nr_tasks;
@@ -138,6 +151,8 @@ void collectgroup1_apply(const struct collectgroup1 *grp1, struct engine *e) {
* @param g_updated the number of updated gravity particles on this node this
* step.
* @param s_updated the number of updated star particles on this node this step.
+ * @param sink_updated the number of updated sink particles on this node this
+ * step.
* @param b_updated the number of updated black hole particles on this node this
* step.
* @param inhibited the number of inhibited hydro particles on this node this
@@ -146,6 +161,8 @@ void collectgroup1_apply(const struct collectgroup1 *grp1, struct engine *e) {
* this step.
* @param s_inhibited the number of inhibited star particles on this node this
* step.
+ * @param sink_inhibited the number of inhibited sink particles on this node
+ * this step.
* @param b_inhibited the number of inhibited black hole particles on this node
* this step.
* @param ti_hydro_end_min the minimum end time for next hydro time step after
@@ -166,6 +183,12 @@ void collectgroup1_apply(const struct collectgroup1 *grp1, struct engine *e) {
* after this step.
* @param ti_stars_beg_max the maximum begin time for next stars time step
* after this step.
+ * @param ti_sinks_end_min the minimum end time for next sinks time step
+ * after this step.
+ * @param ti_sinks_end_max the maximum end time for next sinks time step
+ * after this step.
+ * @param ti_sinks_beg_max the maximum begin time for next sinks time step
+ * after this step.
* @param ti_black_holes_end_min the minimum end time for next black holes time
* step after this step.
* @param ti_black_holes_end_max the maximum end time for next black holes time
@@ -181,25 +204,30 @@ void collectgroup1_apply(const struct collectgroup1 *grp1, struct engine *e) {
*/
void collectgroup1_init(
struct collectgroup1 *grp1, size_t updated, size_t g_updated,
- size_t s_updated, size_t b_updated, size_t inhibited, size_t g_inhibited,
- size_t s_inhibited, size_t b_inhibited, integertime_t ti_hydro_end_min,
+ size_t s_updated, size_t sink_updated, size_t b_updated, size_t inhibited,
+ size_t g_inhibited, size_t s_inhibited, size_t b_inhibited,
+ size_t sink_inhibited, integertime_t ti_hydro_end_min,
integertime_t ti_hydro_end_max, integertime_t ti_hydro_beg_max,
integertime_t ti_gravity_end_min, integertime_t ti_gravity_end_max,
integertime_t ti_gravity_beg_max, integertime_t ti_stars_end_min,
integertime_t ti_stars_end_max, integertime_t ti_stars_beg_max,
- integertime_t ti_black_holes_end_min, integertime_t ti_black_holes_end_max,
- integertime_t ti_black_holes_beg_max, int forcerebuild,
- long long total_nr_cells, long long total_nr_tasks, float tasks_per_cell,
- const struct star_formation_history sfh, float runtime) {
+ integertime_t ti_sinks_end_min, integertime_t ti_sinks_end_max,
+ integertime_t ti_sinks_beg_max, integertime_t ti_black_holes_end_min,
+ integertime_t ti_black_holes_end_max, integertime_t ti_black_holes_beg_max,
+ int forcerebuild, long long total_nr_cells, long long total_nr_tasks,
+ float tasks_per_cell, const struct star_formation_history sfh,
+ float runtime) {
grp1->updated = updated;
grp1->g_updated = g_updated;
grp1->s_updated = s_updated;
grp1->b_updated = b_updated;
+ grp1->sink_updated = sink_updated;
grp1->inhibited = inhibited;
grp1->g_inhibited = g_inhibited;
grp1->s_inhibited = s_inhibited;
grp1->b_inhibited = b_inhibited;
+ grp1->sink_inhibited = sink_inhibited;
grp1->ti_hydro_end_min = ti_hydro_end_min;
grp1->ti_hydro_end_max = ti_hydro_end_max;
grp1->ti_hydro_beg_max = ti_hydro_beg_max;
@@ -212,6 +240,9 @@ void collectgroup1_init(
grp1->ti_black_holes_end_min = ti_black_holes_end_min;
grp1->ti_black_holes_end_max = ti_black_holes_end_max;
grp1->ti_black_holes_beg_max = ti_black_holes_beg_max;
+ grp1->ti_sinks_end_min = ti_sinks_end_min;
+ grp1->ti_sinks_end_max = ti_sinks_end_max;
+ grp1->ti_sinks_beg_max = ti_sinks_beg_max;
grp1->forcerebuild = forcerebuild;
grp1->total_nr_cells = total_nr_cells;
grp1->total_nr_tasks = total_nr_tasks;
@@ -237,22 +268,27 @@ void collectgroup1_reduce(struct collectgroup1 *grp1) {
mpigrp11.updated = grp1->updated;
mpigrp11.g_updated = grp1->g_updated;
mpigrp11.s_updated = grp1->s_updated;
+ mpigrp11.sink_updated = grp1->sink_updated;
mpigrp11.b_updated = grp1->b_updated;
mpigrp11.inhibited = grp1->inhibited;
mpigrp11.g_inhibited = grp1->g_inhibited;
mpigrp11.s_inhibited = grp1->s_inhibited;
+ mpigrp11.sink_inhibited = grp1->sink_inhibited;
mpigrp11.b_inhibited = grp1->b_inhibited;
mpigrp11.ti_hydro_end_min = grp1->ti_hydro_end_min;
mpigrp11.ti_gravity_end_min = grp1->ti_gravity_end_min;
mpigrp11.ti_stars_end_min = grp1->ti_stars_end_min;
+ mpigrp11.ti_sinks_end_min = grp1->ti_sinks_end_min;
mpigrp11.ti_black_holes_end_min = grp1->ti_black_holes_end_min;
mpigrp11.ti_hydro_end_max = grp1->ti_hydro_end_max;
mpigrp11.ti_gravity_end_max = grp1->ti_gravity_end_max;
mpigrp11.ti_stars_end_max = grp1->ti_stars_end_max;
+ mpigrp11.ti_sinks_end_max = grp1->ti_sinks_end_max;
mpigrp11.ti_black_holes_end_max = grp1->ti_black_holes_end_max;
mpigrp11.ti_hydro_beg_max = grp1->ti_hydro_beg_max;
mpigrp11.ti_gravity_beg_max = grp1->ti_gravity_beg_max;
mpigrp11.ti_stars_beg_max = grp1->ti_stars_beg_max;
+ mpigrp11.ti_sinks_beg_max = grp1->ti_sinks_beg_max;
mpigrp11.ti_black_holes_beg_max = grp1->ti_black_holes_beg_max;
mpigrp11.forcerebuild = grp1->forcerebuild;
mpigrp11.total_nr_cells = grp1->total_nr_cells;
@@ -269,23 +305,28 @@ void collectgroup1_reduce(struct collectgroup1 *grp1) {
/* And update. */
grp1->updated = mpigrp12.updated;
grp1->g_updated = mpigrp12.g_updated;
+ grp1->sink_updated = mpigrp12.sink_updated;
grp1->s_updated = mpigrp12.s_updated;
grp1->b_updated = mpigrp12.b_updated;
grp1->inhibited = mpigrp12.inhibited;
grp1->g_inhibited = mpigrp12.g_inhibited;
grp1->s_inhibited = mpigrp12.s_inhibited;
+ grp1->sink_inhibited = mpigrp12.sink_inhibited;
grp1->b_inhibited = mpigrp12.b_inhibited;
grp1->ti_hydro_end_min = mpigrp12.ti_hydro_end_min;
grp1->ti_gravity_end_min = mpigrp12.ti_gravity_end_min;
grp1->ti_stars_end_min = mpigrp12.ti_stars_end_min;
+ grp1->ti_sinks_end_min = mpigrp12.ti_sinks_end_min;
grp1->ti_black_holes_end_min = mpigrp12.ti_black_holes_end_min;
grp1->ti_hydro_end_max = mpigrp12.ti_hydro_end_max;
grp1->ti_gravity_end_max = mpigrp12.ti_gravity_end_max;
grp1->ti_stars_end_max = mpigrp12.ti_stars_end_max;
+ grp1->ti_sinks_end_max = mpigrp12.ti_sinks_end_max;
grp1->ti_black_holes_end_max = mpigrp12.ti_black_holes_end_max;
grp1->ti_hydro_beg_max = mpigrp12.ti_hydro_beg_max;
grp1->ti_gravity_beg_max = mpigrp12.ti_gravity_beg_max;
grp1->ti_stars_beg_max = mpigrp12.ti_stars_beg_max;
+ grp1->ti_sinks_beg_max = mpigrp12.ti_sinks_beg_max;
grp1->ti_black_holes_beg_max = mpigrp12.ti_black_holes_beg_max;
grp1->forcerebuild = mpigrp12.forcerebuild;
grp1->total_nr_cells = mpigrp12.total_nr_cells;
@@ -312,12 +353,14 @@ static void doreduce1(struct mpicollectgroup1 *mpigrp11,
mpigrp11->updated += mpigrp12->updated;
mpigrp11->g_updated += mpigrp12->g_updated;
mpigrp11->s_updated += mpigrp12->s_updated;
+ mpigrp11->sink_updated += mpigrp12->sink_updated;
mpigrp11->b_updated += mpigrp12->b_updated;
/* Sum of inhibited */
mpigrp11->inhibited += mpigrp12->inhibited;
mpigrp11->g_inhibited += mpigrp12->g_inhibited;
mpigrp11->s_inhibited += mpigrp12->s_inhibited;
+ mpigrp11->sink_inhibited += mpigrp12->sink_inhibited;
mpigrp11->b_inhibited += mpigrp12->b_inhibited;
/* Minimum end time. */
@@ -327,6 +370,8 @@ static void doreduce1(struct mpicollectgroup1 *mpigrp11,
min(mpigrp11->ti_gravity_end_min, mpigrp12->ti_gravity_end_min);
mpigrp11->ti_stars_end_min =
min(mpigrp11->ti_stars_end_min, mpigrp12->ti_stars_end_min);
+ mpigrp11->ti_sinks_end_min =
+ min(mpigrp11->ti_sinks_end_min, mpigrp12->ti_sinks_end_min);
mpigrp11->ti_black_holes_end_min =
min(mpigrp11->ti_black_holes_end_min, mpigrp12->ti_black_holes_end_min);
@@ -337,6 +382,8 @@ static void doreduce1(struct mpicollectgroup1 *mpigrp11,
max(mpigrp11->ti_gravity_end_max, mpigrp12->ti_gravity_end_max);
mpigrp11->ti_stars_end_max =
max(mpigrp11->ti_stars_end_max, mpigrp12->ti_stars_end_max);
+ mpigrp11->ti_sinks_end_max =
+ max(mpigrp11->ti_sinks_end_max, mpigrp12->ti_sinks_end_max);
mpigrp11->ti_black_holes_end_max =
max(mpigrp11->ti_black_holes_end_max, mpigrp12->ti_black_holes_end_max);
@@ -347,6 +394,8 @@ static void doreduce1(struct mpicollectgroup1 *mpigrp11,
max(mpigrp11->ti_gravity_beg_max, mpigrp12->ti_gravity_beg_max);
mpigrp11->ti_stars_beg_max =
max(mpigrp11->ti_stars_beg_max, mpigrp12->ti_stars_beg_max);
+ mpigrp11->ti_sinks_beg_max =
+ max(mpigrp11->ti_sinks_beg_max, mpigrp12->ti_sinks_beg_max);
mpigrp11->ti_black_holes_beg_max =
max(mpigrp11->ti_black_holes_beg_max, mpigrp12->ti_black_holes_beg_max);
diff --git a/src/collectgroup.h b/src/collectgroup.h
index fcb49fe0596b23125525c944fa9bfa80418f5b80..99229a257f8544bac749dc330a6255698a12f3da 100644
--- a/src/collectgroup.h
+++ b/src/collectgroup.h
@@ -36,10 +36,10 @@ struct engine;
struct collectgroup1 {
/* Number of particles updated */
- long long updated, g_updated, s_updated, b_updated;
+ long long updated, g_updated, s_updated, b_updated, sink_updated;
/* Number of particles inhibited */
- long long inhibited, g_inhibited, s_inhibited, b_inhibited;
+ long long inhibited, g_inhibited, s_inhibited, b_inhibited, sink_inhibited;
/* SFH logger */
struct star_formation_history sfh;
@@ -50,6 +50,7 @@ struct collectgroup1 {
integertime_t ti_stars_end_min, ti_stars_end_max, ti_stars_beg_max;
integertime_t ti_black_holes_end_min, ti_black_holes_end_max,
ti_black_holes_beg_max;
+ integertime_t ti_sinks_end_min, ti_sinks_end_max, ti_sinks_beg_max;
/* Force the engine to rebuild? */
int forcerebuild;
@@ -69,16 +70,19 @@ void collectgroup_init(void);
void collectgroup1_apply(const struct collectgroup1 *grp1, struct engine *e);
void collectgroup1_init(
struct collectgroup1 *grp1, size_t updated, size_t g_updated,
- size_t s_updated, size_t b_updated, size_t inhibited, size_t g_inhibited,
- size_t s_inhibited, size_t b_inhibited, integertime_t ti_hydro_end_min,
+ size_t s_updated, size_t b_updated, size_t sink_updated, size_t inhibited,
+ size_t g_inhibited, size_t s_inhibited, size_t sink_inhibited,
+ size_t b_inhibited, integertime_t ti_hydro_end_min,
integertime_t ti_hydro_end_max, integertime_t ti_hydro_beg_max,
integertime_t ti_gravity_end_min, integertime_t ti_gravity_end_max,
integertime_t ti_gravity_beg_max, integertime_t ti_stars_end_min,
integertime_t ti_stars_end_max, integertime_t ti_stars_beg_max,
- integertime_t ti_black_holes_end_min, integertime_t ti_black_holes_end_max,
- integertime_t ti_black_holes_beg_max, int forcerebuild,
- long long total_nr_cells, long long total_nr_tasks, float tasks_per_cell,
- const struct star_formation_history sfh, float runtime);
+ integertime_t ti_sinks_end_min, integertime_t ti_sinks_end_max,
+ integertime_t ti_sinks_beg_max, integertime_t ti_black_holes_end_min,
+ integertime_t ti_black_holes_end_max, integertime_t ti_black_holes_beg_max,
+ int forcerebuild, long long total_nr_cells, long long total_nr_tasks,
+ float tasks_per_cell, const struct star_formation_history sfh,
+ float runtime);
void collectgroup1_reduce(struct collectgroup1 *grp1);
#ifdef WITH_MPI
void mpicollect_free_MPI_type(void);
diff --git a/src/common_io.c b/src/common_io.c
index 5db398c71e081ef34fd89240fbd66fe38266e923..8bc4cdc803709642d3e81c1e4b50a63ee3a2c554 100644
--- a/src/common_io.c
+++ b/src/common_io.c
@@ -2883,9 +2883,6 @@ int get_ptype_fields(const int ptype, struct io_props* list,
num_fields += velociraptor_write_gparts(NULL, list + num_fields);
break;
- case 3:
- break;
-
case swift_type_stars:
stars_write_particles(NULL, list, &num_fields, with_cosmology);
num_fields += chemistry_write_sparticles(NULL, list + num_fields);
@@ -2897,6 +2894,10 @@ int get_ptype_fields(const int ptype, struct io_props* list,
num_fields += velociraptor_write_sparts(NULL, list + num_fields);
break;
+ case swift_type_sink:
+ sink_write_particles(NULL, list, &num_fields, with_cosmology);
+ break;
+
case swift_type_black_hole:
black_holes_write_particles(NULL, list, &num_fields, with_cosmology);
num_fields += chemistry_write_bparticles(NULL, list + num_fields);
diff --git a/src/engine.c b/src/engine.c
index bfd5a542ec19eae11f24dbaa35569480cde81cce..662979c3cd49e214fb1b5a436d90626d67a90632 100644
--- a/src/engine.c
+++ b/src/engine.c
@@ -1591,6 +1591,10 @@ int engine_estimate_nr_tasks(const struct engine *e) {
n1 += 6;
#endif
}
+ if (e->policy & engine_policy_sinks) {
+ /* 1 drift, 2 kicks, 1 time-step */
+ n1 += 4;
+ }
if (e->policy & engine_policy_fof) {
n1 += 2;
}
@@ -1697,46 +1701,52 @@ void engine_rebuild(struct engine *e, const int repartitioned,
/* Report the number of particles and memory */
if (e->verbose)
message(
- "Space has memory for %zd/%zd/%zd/%zd part/gpart/spart/bpart "
- "(%zd/%zd/%zd/%zd MB)",
+ "Space has memory for %zd/%zd/%zd/%zd/%zd part/gpart/spart/sink/bpart "
+ "(%zd/%zd/%zd/%zd/%zd MB)",
e->s->size_parts, e->s->size_gparts, e->s->size_sparts,
- e->s->size_bparts,
+ e->s->size_sinks, e->s->size_bparts,
e->s->size_parts * sizeof(struct part) / (1024 * 1024),
e->s->size_gparts * sizeof(struct gpart) / (1024 * 1024),
e->s->size_sparts * sizeof(struct spart) / (1024 * 1024),
+ e->s->size_sinks * sizeof(struct sink) / (1024 * 1024),
e->s->size_bparts * sizeof(struct bpart) / (1024 * 1024));
if (e->verbose)
message(
- "Space holds %zd/%zd/%zd/%zd part/gpart/spart/bpart (fracs: "
- "%f/%f/%f/%f)",
- e->s->nr_parts, e->s->nr_gparts, e->s->nr_sparts, e->s->nr_bparts,
+ "Space holds %zd/%zd/%zd/%zd/%zd part/gpart/spart/sink/bpart (fracs: "
+ "%f/%f/%f/%f/%f)",
+ e->s->nr_parts, e->s->nr_gparts, e->s->nr_sparts, e->s->nr_sinks,
+ e->s->nr_bparts,
e->s->nr_parts ? e->s->nr_parts / ((double)e->s->size_parts) : 0.,
e->s->nr_gparts ? e->s->nr_gparts / ((double)e->s->size_gparts) : 0.,
e->s->nr_sparts ? e->s->nr_sparts / ((double)e->s->size_sparts) : 0.,
+ e->s->nr_sinks ? e->s->nr_sinks / ((double)e->s->size_sinks) : 0.,
e->s->nr_bparts ? e->s->nr_bparts / ((double)e->s->size_bparts) : 0.);
const ticks tic2 = getticks();
/* Update the global counters of particles */
- long long num_particles[4] = {
+ long long num_particles[5] = {
(long long)(e->s->nr_parts - e->s->nr_extra_parts),
(long long)(e->s->nr_gparts - e->s->nr_extra_gparts),
(long long)(e->s->nr_sparts - e->s->nr_extra_sparts),
+ (long long)(e->s->nr_sinks - e->s->nr_extra_sinks),
(long long)(e->s->nr_bparts - e->s->nr_extra_bparts)};
#ifdef WITH_MPI
- MPI_Allreduce(MPI_IN_PLACE, num_particles, 4, MPI_LONG_LONG, MPI_SUM,
+ MPI_Allreduce(MPI_IN_PLACE, num_particles, 5, MPI_LONG_LONG, MPI_SUM,
MPI_COMM_WORLD);
#endif
e->total_nr_parts = num_particles[0];
e->total_nr_gparts = num_particles[1];
e->total_nr_sparts = num_particles[2];
- e->total_nr_bparts = num_particles[3];
+ e->total_nr_sinks = num_particles[3];
+ e->total_nr_bparts = num_particles[4];
/* Flag that there are no inhibited particles */
e->nr_inhibited_parts = 0;
e->nr_inhibited_gparts = 0;
e->nr_inhibited_sparts = 0;
+ e->nr_inhibited_sinks = 0;
e->nr_inhibited_bparts = 0;
if (e->verbose)
@@ -1822,6 +1832,7 @@ void engine_rebuild(struct engine *e, const int repartitioned,
e->updates_since_rebuild = 0;
e->g_updates_since_rebuild = 0;
e->s_updates_since_rebuild = 0;
+ e->sink_updates_since_rebuild = 0;
e->b_updates_since_rebuild = 0;
/* Flag that a rebuild has taken place */
@@ -2031,8 +2042,8 @@ void engine_skip_force_and_kick(struct engine *e) {
/* Skip everything that updates the particles */
if (t->type == task_type_drift_part || t->type == task_type_drift_gpart ||
t->type == task_type_drift_spart || t->type == task_type_drift_bpart ||
- t->type == task_type_kick1 || t->type == task_type_kick2 ||
- t->type == task_type_timestep ||
+ t->type == task_type_drift_sink || t->type == task_type_kick1 ||
+ t->type == task_type_kick2 || t->type == task_type_timestep ||
t->type == task_type_timestep_limiter ||
t->type == task_type_timestep_sync ||
t->type == task_type_end_hydro_force || t->type == task_type_cooling ||
@@ -2061,6 +2072,7 @@ void engine_skip_force_and_kick(struct engine *e) {
t->subtype == task_subtype_tend_part ||
t->subtype == task_subtype_tend_gpart ||
t->subtype == task_subtype_tend_spart ||
+ t->subtype == task_subtype_tend_sink ||
t->subtype == task_subtype_tend_bpart ||
t->subtype == task_subtype_rho || t->subtype == task_subtype_sf_counts)
t->skip = 1;
@@ -2087,7 +2099,8 @@ void engine_skip_drift(struct engine *e) {
/* Skip everything that moves the particles */
if (t->type == task_type_drift_part || t->type == task_type_drift_gpart ||
- t->type == task_type_drift_spart || t->type == task_type_drift_bpart)
+ t->type == task_type_drift_spart || t->type == task_type_drift_bpart ||
+ t->type == task_type_drift_sink)
t->skip = 1;
}
@@ -2473,10 +2486,11 @@ 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 %21.3f %6d\n",
+ "%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->s_updates, e->sink_updates, e->b_updates, e->wallclock_time,
+ e->step_props);
#ifdef SWIFT_DEBUG_CHECKS
fflush(stdout);
#endif
@@ -2499,10 +2513,11 @@ void engine_step(struct engine *e) {
fprintf(
e->file_timesteps,
" %6d %14e %12.7f %12.7f %14e %4d %4d %12lld %12lld %12lld %12lld "
- "%21.3f %6d\n",
+ "%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->s_updates, e->sink_updates, e->b_updates, e->wallclock_time,
+ e->step_props);
#ifdef SWIFT_DEBUG_CHECKS
fflush(e->file_timesteps);
#endif
@@ -2731,6 +2746,7 @@ void engine_step(struct engine *e) {
e->updates_since_rebuild += e->collect_group1.updated;
e->g_updates_since_rebuild += e->collect_group1.g_updated;
e->s_updates_since_rebuild += e->collect_group1.s_updated;
+ e->sink_updates_since_rebuild += e->collect_group1.sink_updated;
e->b_updates_since_rebuild += e->collect_group1.b_updated;
#ifdef SWIFT_DEBUG_CHECKS
@@ -4323,12 +4339,12 @@ void engine_config(int restart, int fof, struct engine *e,
engine_step_prop_stf, engine_step_prop_fof,
engine_step_prop_logger_index);
- fprintf(
- e->file_timesteps,
- "# %6s %14s %12s %12s %14s %9s %12s %12s %12s %12s %16s [%s] %6s\n",
- "Step", "Time", "Scale-factor", "Redshift", "Time-step", "Time-bins",
- "Updates", "g-Updates", "s-Updates", "b-Updates", "Wall-clock time",
- clocks_getunit(), "Props");
+ fprintf(e->file_timesteps,
+ "# %6s %14s %12s %12s %14s %9s %12s %12s %12s %12s %12s %16s "
+ "[%s] %6s\n",
+ "Step", "Time", "Scale-factor", "Redshift", "Time-step",
+ "Time-bins", "Updates", "g-Updates", "s-Updates", "Sink-Updates",
+ "b-Updates", "Wall-clock time", clocks_getunit(), "Props");
fflush(e->file_timesteps);
}
@@ -5214,8 +5230,8 @@ void engine_recompute_displacement_constraint(struct engine *e) {
/* Start by reducing the minimal mass of each particle type */
float min_mass[swift_type_count] = {
- e->s->min_part_mass, e->s->min_gpart_mass, FLT_MAX, FLT_MAX,
- e->s->min_spart_mass, e->s->min_bpart_mass};
+ e->s->min_part_mass, e->s->min_gpart_mass, FLT_MAX,
+ e->s->min_sink_mass, e->s->min_spart_mass, e->s->min_bpart_mass};
#ifdef WITH_MPI
MPI_Allreduce(MPI_IN_PLACE, min_mass, swift_type_count, MPI_FLOAT, MPI_MIN,
@@ -5235,9 +5251,12 @@ void engine_recompute_displacement_constraint(struct engine *e) {
#endif
/* Do the same for the velocity norm sum */
- float vel_norm[swift_type_count] = {
- e->s->sum_part_vel_norm, e->s->sum_gpart_vel_norm, 0.f, 0.f,
- e->s->sum_spart_vel_norm, e->s->sum_spart_vel_norm};
+ float vel_norm[swift_type_count] = {e->s->sum_part_vel_norm,
+ e->s->sum_gpart_vel_norm,
+ 0.f,
+ e->s->sum_sink_vel_norm,
+ e->s->sum_spart_vel_norm,
+ e->s->sum_spart_vel_norm};
#ifdef WITH_MPI
MPI_Allreduce(MPI_IN_PLACE, vel_norm, swift_type_count, MPI_FLOAT, MPI_SUM,
MPI_COMM_WORLD);
@@ -5252,17 +5271,19 @@ void engine_recompute_displacement_constraint(struct engine *e) {
(float)e->total_nr_parts,
(float)total_nr_dm_gparts,
(float)e->total_nr_DM_background_gparts,
- 0.f,
+ (float)e->total_nr_sinks,
(float)e->total_nr_sparts,
(float)e->total_nr_bparts};
/* Count of particles for the two species */
const float N_dm = count_parts[1];
- const float N_b = count_parts[0] + count_parts[4] + count_parts[5];
+ const float N_b =
+ count_parts[0] + count_parts[3] + count_parts[4] + count_parts[5];
/* Peculiar motion norm for the two species */
const float vel_norm_dm = vel_norm[1];
- const float vel_norm_b = vel_norm[0] + vel_norm[4] + vel_norm[5];
+ const float vel_norm_b =
+ vel_norm[0] + vel_norm[3] + vel_norm[4] + vel_norm[5];
/* Mesh forces smoothing scale */
float r_s;
@@ -5293,7 +5314,8 @@ void engine_recompute_displacement_constraint(struct engine *e) {
if (N_b > 0.f) {
/* Minimal mass for the baryons */
- const float min_mass_b = min3(min_mass[0], min_mass[4], min_mass[5]);
+ const float min_mass_b =
+ min4(min_mass[0], min_mass[3], min_mass[4], min_mass[5]);
/* Inter-particle sepration for the baryons */
const float d_b = cbrtf(min_mass_b / (Ob * rho_crit0));
diff --git a/src/engine.h b/src/engine.h
index 9dae02eef75004b4afc1e347a8453b0a465e8307..9f9c04fa20d4029db726e66ec340a7f56ad08381 100644
--- a/src/engine.h
+++ b/src/engine.h
@@ -81,7 +81,7 @@ enum engine_policy {
engine_policy_timestep_sync = (1 << 22),
engine_policy_logger = (1 << 23),
engine_policy_line_of_sight = (1 << 24),
- engine_policy_sink = (1 << 25),
+ engine_policy_sinks = (1 << 25),
engine_policy_rt = (1 << 26),
};
#define engine_maxpolicy 27
@@ -224,6 +224,15 @@ struct engine {
/* Maximal black holes ti_beg for the next time-step */
integertime_t ti_black_holes_beg_max;
+ /* Minimal sinks ti_end for the next time-step */
+ integertime_t ti_sinks_end_min;
+
+ /* Maximal sinks ti_end for the next time-step */
+ integertime_t ti_sinks_end_max;
+
+ /* Maximal sinks ti_beg for the next time-step */
+ integertime_t ti_sinks_beg_max;
+
/* Minimal overall ti_end for the next time-step */
integertime_t ti_end_min;
@@ -234,12 +243,13 @@ struct engine {
integertime_t ti_beg_max;
/* Number of particles updated in the previous step */
- long long updates, g_updates, s_updates, b_updates;
+ long long updates, g_updates, s_updates, b_updates, sink_updates;
/* Number of updates since the last rebuild */
long long updates_since_rebuild;
long long g_updates_since_rebuild;
long long s_updates_since_rebuild;
+ long long sink_updates_since_rebuild;
long long b_updates_since_rebuild;
/* Star formation logger information */
@@ -252,6 +262,7 @@ struct engine {
long long total_nr_parts;
long long total_nr_gparts;
long long total_nr_sparts;
+ long long total_nr_sinks;
long long total_nr_bparts;
long long total_nr_DM_background_gparts;
@@ -265,6 +276,7 @@ struct engine {
long long nr_inhibited_parts;
long long nr_inhibited_gparts;
long long nr_inhibited_sparts;
+ long long nr_inhibited_sinks;
long long nr_inhibited_bparts;
#ifdef SWIFT_DEBUG_CHECKS
@@ -510,9 +522,6 @@ struct engine {
integertime_t ti_next_los;
int los_output_count;
- /* Total number of sink particles in the system. */
- long long total_nr_sinks;
-
#ifdef SWIFT_GRAVITY_FORCE_CHECKS
/* Run brute force checks only on steps when all gparts active? */
int force_checks_only_all_active;
diff --git a/src/engine_collect_end_of_step.c b/src/engine_collect_end_of_step.c
index 0973399c11146baef3b3be60b3635426021b23cb..4281379f98b70b8f3ce169822fa435aad2557708 100644
--- a/src/engine_collect_end_of_step.c
+++ b/src/engine_collect_end_of_step.c
@@ -35,11 +35,12 @@
*/
struct end_of_step_data {
- size_t updated, g_updated, s_updated, b_updated;
- size_t inhibited, g_inhibited, s_inhibited, b_inhibited;
+ size_t updated, g_updated, s_updated, sink_updated, b_updated;
+ size_t inhibited, g_inhibited, s_inhibited, sink_inhibited, b_inhibited;
integertime_t ti_hydro_end_min, ti_hydro_end_max, ti_hydro_beg_max;
integertime_t ti_gravity_end_min, ti_gravity_end_max, ti_gravity_beg_max;
integertime_t ti_stars_end_min, ti_stars_end_max, ti_stars_beg_max;
+ integertime_t ti_sinks_end_min, ti_sinks_end_max, ti_sinks_beg_max;
integertime_t ti_black_holes_end_min, ti_black_holes_end_max,
ti_black_holes_beg_max;
struct engine *e;
@@ -273,6 +274,61 @@ void engine_collect_end_of_step_recurse_black_holes(struct cell *c,
c->black_holes.updated = updated;
}
+/**
+ * @brief Recursive function gathering end-of-step data.
+ *
+ * We recurse until we encounter a timestep or time-step MPI recv task
+ * as the values will have been set at that level. We then bring these
+ * values upwards.
+ *
+ * @param c The #cell to recurse into.
+ * @param e The #engine.
+ */
+void engine_collect_end_of_step_recurse_sinks(struct cell *c,
+ const struct engine *e) {
+
+ /* Skip super-cells (Their values are already set) */
+ if (c->timestep != NULL) return;
+#ifdef WITH_MPI
+ if (cell_get_recv(c, task_subtype_tend_sink) != NULL) return;
+#endif /* WITH_MPI */
+
+#ifdef SWIFT_DEBUG_CHECKS
+ // if (!c->split) error("Reached a leaf without finding a time-step task!");
+#endif
+
+ /* Counters for the different quantities. */
+ size_t updated = 0;
+ integertime_t ti_sinks_end_min = max_nr_timesteps, ti_sinks_end_max = 0,
+ ti_sinks_beg_max = 0;
+
+ /* Collect the values from the progeny. */
+ for (int k = 0; k < 8; k++) {
+ struct cell *cp = c->progeny[k];
+ if (cp != NULL && cp->sinks.count > 0) {
+
+ /* Recurse */
+ engine_collect_end_of_step_recurse_sinks(cp, e);
+
+ /* And update */
+ ti_sinks_end_min = min(ti_sinks_end_min, cp->sinks.ti_end_min);
+ ti_sinks_end_max = max(ti_sinks_end_max, cp->sinks.ti_end_max);
+ ti_sinks_beg_max = max(ti_sinks_beg_max, cp->sinks.ti_beg_max);
+
+ updated += cp->sinks.updated;
+
+ /* Collected, so clear for next time. */
+ cp->sinks.updated = 0;
+ }
+ }
+
+ /* Store the collected values in the cell. */
+ c->sinks.ti_end_min = ti_sinks_end_min;
+ c->sinks.ti_end_max = ti_sinks_end_max;
+ c->sinks.ti_beg_max = ti_sinks_beg_max;
+ c->sinks.updated = updated;
+}
+
/**
* @brief Mapping function to collect the data from the end of the step
*
@@ -294,19 +350,23 @@ void engine_collect_end_of_step_mapper(void *map_data, int num_elements,
const int with_ext_grav = (e->policy & engine_policy_external_gravity);
const int with_grav = (with_self_grav || with_ext_grav);
const int with_stars = (e->policy & engine_policy_stars);
+ const int with_sinks = (e->policy & engine_policy_sinks);
const int with_black_holes = (e->policy & engine_policy_black_holes);
struct space *s = e->s;
int *local_cells = (int *)map_data;
struct star_formation_history *sfh_top = &data->sfh;
/* Local collectible */
- size_t updated = 0, g_updated = 0, s_updated = 0, b_updated = 0;
+ size_t updated = 0, g_updated = 0, s_updated = 0, sink_updated = 0,
+ b_updated = 0;
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;
integertime_t ti_stars_end_min = max_nr_timesteps, ti_stars_end_max = 0,
ti_stars_beg_max = 0;
+ integertime_t ti_sinks_end_min = max_nr_timesteps, ti_sinks_end_max = 0,
+ ti_sinks_beg_max = 0;
integertime_t ti_black_holes_end_min = max_nr_timesteps,
ti_black_holes_end_max = 0, ti_black_holes_beg_max = 0;
@@ -320,7 +380,7 @@ void engine_collect_end_of_step_mapper(void *map_data, int num_elements,
struct cell *c = &s->cells_top[local_cells[ind]];
if (c->hydro.count > 0 || c->grav.count > 0 || c->stars.count > 0 ||
- c->black_holes.count > 0) {
+ c->black_holes.count > 0 || c->sinks.count > 0) {
/* Make the top-cells recurse */
if (with_hydro) {
@@ -332,6 +392,9 @@ void engine_collect_end_of_step_mapper(void *map_data, int num_elements,
if (with_stars) {
engine_collect_end_of_step_recurse_stars(c, e);
}
+ if (with_sinks) {
+ engine_collect_end_of_step_recurse_sinks(c, e);
+ }
if (with_black_holes) {
engine_collect_end_of_step_recurse_black_holes(c, e);
}
@@ -352,6 +415,11 @@ void engine_collect_end_of_step_mapper(void *map_data, int num_elements,
ti_stars_end_max = max(ti_stars_end_max, c->stars.ti_end_max);
ti_stars_beg_max = max(ti_stars_beg_max, c->stars.ti_beg_max);
+ if (c->sinks.ti_end_min > e->ti_current)
+ ti_sinks_end_min = min(ti_sinks_end_min, c->sinks.ti_end_min);
+ ti_sinks_end_max = max(ti_sinks_end_max, c->sinks.ti_end_max);
+ ti_sinks_beg_max = max(ti_sinks_beg_max, c->sinks.ti_beg_max);
+
if (c->black_holes.ti_end_min > e->ti_current)
ti_black_holes_end_min =
min(ti_black_holes_end_min, c->black_holes.ti_end_min);
@@ -363,6 +431,7 @@ void engine_collect_end_of_step_mapper(void *map_data, int num_elements,
updated += c->hydro.updated;
g_updated += c->grav.updated;
s_updated += c->stars.updated;
+ sink_updated += c->sinks.updated;
b_updated += c->black_holes.updated;
/* Check if the cell was inactive and in that case reorder the SFH */
@@ -378,6 +447,7 @@ void engine_collect_end_of_step_mapper(void *map_data, int num_elements,
c->hydro.updated = 0;
c->grav.updated = 0;
c->stars.updated = 0;
+ c->sinks.updated = 0;
c->black_holes.updated = 0;
}
}
@@ -388,6 +458,7 @@ void engine_collect_end_of_step_mapper(void *map_data, int num_elements,
data->updated += updated;
data->g_updated += g_updated;
data->s_updated += s_updated;
+ data->sink_updated += sink_updated;
data->b_updated += b_updated;
/* Add the SFH information from this engine to the global data */
@@ -411,6 +482,11 @@ void engine_collect_end_of_step_mapper(void *map_data, int num_elements,
data->ti_stars_end_max = max(ti_stars_end_max, data->ti_stars_end_max);
data->ti_stars_beg_max = max(ti_stars_beg_max, data->ti_stars_beg_max);
+ if (ti_sinks_end_min > e->ti_current)
+ data->ti_sinks_end_min = min(ti_sinks_end_min, data->ti_sinks_end_min);
+ data->ti_sinks_end_max = max(ti_sinks_end_max, data->ti_sinks_end_max);
+ data->ti_sinks_beg_max = max(ti_sinks_beg_max, data->ti_sinks_beg_max);
+
if (ti_black_holes_end_min > e->ti_current)
data->ti_black_holes_end_min =
min(ti_black_holes_end_min, data->ti_black_holes_end_min);
@@ -446,12 +522,15 @@ void engine_collect_end_of_step(struct engine *e, int apply) {
struct space *s = e->s;
struct end_of_step_data data;
data.updated = 0, data.g_updated = 0, data.s_updated = 0, data.b_updated = 0;
+ data.sink_updated = 0;
data.ti_hydro_end_min = max_nr_timesteps, data.ti_hydro_end_max = 0,
data.ti_hydro_beg_max = 0;
data.ti_gravity_end_min = max_nr_timesteps, data.ti_gravity_end_max = 0,
data.ti_gravity_beg_max = 0;
data.ti_stars_end_min = max_nr_timesteps, data.ti_stars_end_max = 0,
data.ti_stars_beg_max = 0;
+ data.ti_sinks_end_min = max_nr_timesteps, data.ti_sinks_end_max = 0,
+ data.ti_sinks_beg_max = 0;
data.ti_black_holes_end_min = max_nr_timesteps,
data.ti_black_holes_end_max = 0, data.ti_black_holes_beg_max = 0;
data.e = e;
@@ -472,20 +551,22 @@ void engine_collect_end_of_step(struct engine *e, int apply) {
data.inhibited = s->nr_inhibited_parts;
data.g_inhibited = s->nr_inhibited_gparts;
data.s_inhibited = s->nr_inhibited_sparts;
+ data.sink_inhibited = s->nr_inhibited_sinks;
data.b_inhibited = s->nr_inhibited_bparts;
/* Store these in the temporary collection group. */
collectgroup1_init(
&e->collect_group1, data.updated, data.g_updated, data.s_updated,
- data.b_updated, data.inhibited, data.g_inhibited, data.s_inhibited,
- data.b_inhibited, data.ti_hydro_end_min, data.ti_hydro_end_max,
- data.ti_hydro_beg_max, data.ti_gravity_end_min, data.ti_gravity_end_max,
- data.ti_gravity_beg_max, data.ti_stars_end_min, data.ti_stars_end_max,
- data.ti_stars_beg_max, data.ti_black_holes_end_min,
- data.ti_black_holes_end_max, data.ti_black_holes_beg_max, e->forcerebuild,
- e->s->tot_cells, e->sched.nr_tasks,
- (float)e->sched.nr_tasks / (float)e->s->tot_cells, data.sfh,
- data.runtime);
+ data.sink_updated, data.b_updated, data.inhibited, data.g_inhibited,
+ data.s_inhibited, data.sink_inhibited, data.b_inhibited,
+ data.ti_hydro_end_min, data.ti_hydro_end_max, data.ti_hydro_beg_max,
+ data.ti_gravity_end_min, data.ti_gravity_end_max, data.ti_gravity_beg_max,
+ data.ti_stars_end_min, data.ti_stars_end_max, data.ti_stars_beg_max,
+ data.ti_sinks_end_min, data.ti_sinks_end_max, data.ti_sinks_beg_max,
+ data.ti_black_holes_end_min, data.ti_black_holes_end_max,
+ data.ti_black_holes_beg_max, e->forcerebuild, e->s->tot_cells,
+ e->sched.nr_tasks, (float)e->sched.nr_tasks / (float)e->s->tot_cells,
+ data.sfh, data.runtime);
/* Aggregate collective data from the different nodes for this step. */
#ifdef WITH_MPI
diff --git a/src/engine_drift.c b/src/engine_drift.c
index 0e447d6235b6e3f4b1f71dfea173daa8971d56f7..ccabb07575afc638e8162d96516c380619d95422 100644
--- a/src/engine_drift.c
+++ b/src/engine_drift.c
@@ -220,6 +220,54 @@ void engine_do_drift_all_bpart_mapper(void *map_data, int num_elements,
}
}
+/**
+ * @brief Mapper function to drift *all* the #sink to the current time.
+ *
+ * @param map_data An array of #cell%s.
+ * @param num_elements Chunk size.
+ * @param extra_data Pointer to an #engine.
+ */
+void engine_do_drift_all_sink_mapper(void *map_data, int num_elements,
+ void *extra_data) {
+
+ const struct engine *e = (const struct engine *)extra_data;
+ const int restarting = e->restarting;
+ struct space *s = e->s;
+ struct cell *cells_top;
+ int *local_cells_top;
+
+ if (restarting) {
+
+ /* When restarting, we loop over all top-level cells */
+ cells_top = (struct cell *)map_data;
+ local_cells_top = NULL;
+
+ } else {
+
+ /* In any other case, we use the list of local cells with tasks */
+ cells_top = s->cells_top;
+ local_cells_top = (int *)map_data;
+ }
+
+ for (int ind = 0; ind < num_elements; ind++) {
+
+ struct cell *c;
+
+ /* When restarting, the list of local cells does not
+ yet exist. We use the raw list of top-level cells instead */
+ if (restarting)
+ c = &cells_top[ind];
+ else
+ c = &cells_top[local_cells_top[ind]];
+
+ if (c->nodeID == e->nodeID) {
+
+ /* Drift all the particles */
+ cell_drift_sink(c, e, /* force the drift=*/1);
+ }
+ }
+}
+
/**
* @brief Mapper function to drift *all* the multipoles to the current time.
*
@@ -305,6 +353,11 @@ void engine_drift_all(struct engine *e, const int drift_mpoles) {
e->s->local_cells_top, e->s->nr_local_cells, sizeof(int),
threadpool_auto_chunk_size, e);
}
+ if (e->s->nr_sinks > 0) {
+ threadpool_map(&e->threadpool, engine_do_drift_all_sink_mapper,
+ e->s->local_cells_top, e->s->nr_local_cells, sizeof(int),
+ threadpool_auto_chunk_size, e);
+ }
if (e->s->nr_bparts > 0) {
threadpool_map(&e->threadpool, engine_do_drift_all_bpart_mapper,
e->s->local_cells_top, e->s->nr_local_cells, sizeof(int),
@@ -332,6 +385,11 @@ void engine_drift_all(struct engine *e, const int drift_mpoles) {
e->s->cells_top, e->s->nr_cells, sizeof(struct cell),
threadpool_auto_chunk_size, e);
}
+ if (e->s->nr_sinks > 0) {
+ threadpool_map(&e->threadpool, engine_do_drift_all_sink_mapper,
+ e->s->cells_top, e->s->nr_cells, sizeof(struct cell),
+ threadpool_auto_chunk_size, e);
+ }
if (e->s->nr_bparts > 0) {
threadpool_map(&e->threadpool, engine_do_drift_all_bpart_mapper,
e->s->cells_top, e->s->nr_cells, sizeof(struct cell),
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index b38cfca652f7685d943a772f334ed4e5636fff7a..f04d74261213b0290c06cceaab2a41d70192429a 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -1096,6 +1096,7 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c,
struct scheduler *s = &e->sched;
const int with_stars = (e->policy & engine_policy_stars);
+ const int with_sinks = (e->policy & engine_policy_sinks);
const int with_feedback = (e->policy & engine_policy_feedback);
const int with_cooling = (e->policy & engine_policy_cooling);
const int with_star_formation = (e->policy & engine_policy_star_formation);
@@ -1174,6 +1175,13 @@ void engine_make_hierarchical_tasks_hydro(struct engine *e, struct cell *c,
scheduler_addunlock(s, c->stars.drift, c->super->kick2);
}
+ /* Sinks */
+ if (with_sinks) {
+ c->sinks.drift = scheduler_addtask(s, task_type_drift_sink,
+ task_subtype_none, 0, 0, c, NULL);
+ scheduler_addunlock(s, c->sinks.drift, c->super->kick2);
+ }
+
/* Black holes */
if (with_black_holes) {
c->black_holes.drift = scheduler_addtask(
@@ -2878,6 +2886,7 @@ void engine_make_hydroloop_tasks_mapper(void *map_data, int num_elements,
const int periodic = e->s->periodic;
const int with_feedback = (e->policy & engine_policy_feedback);
const int with_stars = (e->policy & engine_policy_stars);
+ const int with_sinks = (e->policy & engine_policy_sinks);
const int with_black_holes = (e->policy & engine_policy_black_holes);
struct space *s = e->s;
@@ -2902,6 +2911,7 @@ void engine_make_hydroloop_tasks_mapper(void *map_data, int num_elements,
/* Skip cells without hydro or star particles */
if ((ci->hydro.count == 0) && (!with_stars || ci->stars.count == 0) &&
+ (!with_sinks || ci->sinks.count == 0) &&
(!with_black_holes || ci->black_holes.count == 0))
continue;
@@ -2933,6 +2943,7 @@ void engine_make_hydroloop_tasks_mapper(void *map_data, int num_elements,
if ((cid >= cjd) ||
((cj->hydro.count == 0) &&
(!with_feedback || cj->stars.count == 0) &&
+ (!with_sinks || cj->sinks.count == 0) &&
(!with_black_holes || cj->black_holes.count == 0)) ||
(ci->nodeID != nodeID && cj->nodeID != nodeID))
continue;
diff --git a/src/engine_marktasks.c b/src/engine_marktasks.c
index 3fb51038b6296d39c70ab740c69b88eca50b5659..115725602a93abecc89abba35a8ae57b5c7e6e54 100644
--- a/src/engine_marktasks.c
+++ b/src/engine_marktasks.c
@@ -878,11 +878,17 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
else if (t_type == task_type_kick1 || t_type == task_type_kick2) {
if (cell_is_active_hydro(t->ci, e) || cell_is_active_gravity(t->ci, e) ||
- cell_is_active_stars(t->ci, e) ||
+ cell_is_active_stars(t->ci, e) || cell_is_active_sinks(t->ci, e) ||
cell_is_active_black_holes(t->ci, e))
scheduler_activate(s, t);
}
+ /* Sink drift ? */
+ else if (t_type == task_type_drift_sink) {
+
+ if (cell_is_active_sinks(t->ci, e)) cell_activate_drift_sink(t->ci, s);
+ }
+
/* Hydro ghost tasks ? */
else if (t_type == task_type_ghost || t_type == task_type_extra_ghost ||
t_type == task_type_ghost_in || t_type == task_type_ghost_out) {
@@ -959,9 +965,10 @@ void engine_marktasks_mapper(void *map_data, int num_elements,
t->ci->hydro.updated = 0;
t->ci->grav.updated = 0;
t->ci->stars.updated = 0;
+ t->ci->sinks.updated = 0;
t->ci->black_holes.updated = 0;
if (cell_is_active_hydro(t->ci, e) || cell_is_active_gravity(t->ci, e) ||
- cell_is_active_stars(t->ci, e) ||
+ cell_is_active_stars(t->ci, e) || cell_is_active_sinks(t->ci, e) ||
cell_is_active_black_holes(t->ci, e))
scheduler_activate(s, t);
}
diff --git a/src/engine_redistribute.c b/src/engine_redistribute.c
index 35354ec848d8890b9157c5c3a6825f84094cdd89..a7d954655f6dce4b18d16814106dde3098785c30 100644
--- a/src/engine_redistribute.c
+++ b/src/engine_redistribute.c
@@ -521,7 +521,7 @@ void engine_redistribute(struct engine *e) {
#ifdef SWIFT_DEBUG_CHECKS
const int nr_sinks_new = 0;
#endif
- if (e->policy & engine_policy_sink) {
+ if (e->policy & engine_policy_sinks) {
error("Not implemented yet");
}
diff --git a/src/engine_unskip.c b/src/engine_unskip.c
index 9142ff447ac023d0428d466029110233cef37149..bff3bf842b746b8b618b27a42dcf0d3a6bd19f18 100644
--- a/src/engine_unskip.c
+++ b/src/engine_unskip.c
@@ -44,6 +44,7 @@ enum task_broad_types {
task_broad_types_hydro = 1,
task_broad_types_gravity,
task_broad_types_stars,
+ task_broad_types_sinks,
task_broad_types_black_holes,
task_broad_types_count,
};
@@ -174,6 +175,38 @@ static void engine_do_unskip_black_holes(struct cell *c, struct engine *e) {
if (forcerebuild) atomic_inc(&e->forcerebuild);
}
+/**
+ * @brief Unskip any sink tasks associated with active cells.
+ *
+ * @param c The cell.
+ * @param e The engine.
+ */
+static void engine_do_unskip_sinks(struct cell *c, struct engine *e) {
+
+ /* Early abort (are we below the level where tasks are)? */
+ if (!cell_get_flag(c, cell_flag_has_tasks)) return;
+
+ /* Ignore empty cells. */
+ if (c->sinks.count == 0) return;
+
+ /* Skip inactive cells. */
+ if (!cell_is_active_sinks(c, e)) return;
+
+ /* Recurse */
+ if (c->split) {
+ for (int k = 0; k < 8; k++) {
+ if (c->progeny[k] != NULL) {
+ struct cell *cp = c->progeny[k];
+ engine_do_unskip_sinks(cp, e);
+ }
+ }
+ }
+
+ /* Unskip any active tasks. */
+ const int forcerebuild = cell_unskip_sinks_tasks(c, &e->sched);
+ if (forcerebuild) atomic_inc(&e->forcerebuild);
+}
+
/**
* @brief Unskip any gravity tasks associated with active cells.
*
@@ -269,6 +302,13 @@ void engine_do_unskip_mapper(void *map_data, int num_elements,
#endif
engine_do_unskip_stars(c, e, with_star_formation);
break;
+ case task_broad_types_sinks:
+#ifdef SWIFT_DEBUG_CHECKS
+ if (!(e->policy & engine_policy_sinks))
+ error("Trying to unskip sink tasks in a non-sinks run!");
+#endif
+ engine_do_unskip_sinks(c, e);
+ break;
case task_broad_types_black_holes:
#ifdef SWIFT_DEBUG_CHECKS
if (!(e->policy & engine_policy_black_holes))
@@ -300,6 +340,7 @@ void engine_unskip(struct engine *e) {
const int with_self_grav = e->policy & engine_policy_self_gravity;
const int with_ext_grav = e->policy & engine_policy_external_gravity;
const int with_stars = e->policy & engine_policy_stars;
+ const int with_sinks = e->policy & engine_policy_sinks;
const int with_feedback = e->policy & engine_policy_feedback;
const int with_black_holes = e->policy & engine_policy_black_holes;
@@ -322,6 +363,7 @@ void engine_unskip(struct engine *e) {
cell_is_active_gravity(c, e)) ||
(with_feedback && cell_is_active_stars(c, e)) ||
(with_stars && c->nodeID == nodeID && cell_is_active_stars(c, e)) ||
+ (with_sinks && cell_is_active_sinks(c, e)) ||
(with_black_holes && cell_is_active_black_holes(c, e))) {
if (num_active_cells != k)
@@ -346,6 +388,10 @@ void engine_unskip(struct engine *e) {
data.task_types[multiplier] = task_broad_types_stars;
multiplier++;
}
+ if (with_sinks) {
+ data.task_types[multiplier] = task_broad_types_sinks;
+ multiplier++;
+ }
if (with_black_holes) {
data.task_types[multiplier] = task_broad_types_black_holes;
multiplier++;
diff --git a/src/gravity/MultiSoftening/gravity.h b/src/gravity/MultiSoftening/gravity.h
index 6ed0c38b7ee90807cc412efeb765167ded77e70c..e6dfccd81e129a8dfde316e76d3b54f119ed0a62 100644
--- a/src/gravity/MultiSoftening/gravity.h
+++ b/src/gravity/MultiSoftening/gravity.h
@@ -238,6 +238,9 @@ __attribute__((always_inline)) INLINE static void gravity_predict_extra(
case swift_type_dark_matter:
gp->epsilon = grav_props->epsilon_DM_cur;
break;
+ case swift_type_sink:
+ gp->epsilon = grav_props->epsilon_baryon_cur;
+ break;
case swift_type_stars:
gp->epsilon = grav_props->epsilon_baryon_cur;
break;
@@ -298,6 +301,9 @@ __attribute__((always_inline)) INLINE static void gravity_first_init_gpart(
case swift_type_stars:
gp->epsilon = grav_props->epsilon_baryon_cur;
break;
+ case swift_type_sink:
+ gp->epsilon = grav_props->epsilon_baryon_cur;
+ break;
case swift_type_gas:
gp->epsilon = grav_props->epsilon_baryon_cur;
break;
diff --git a/src/runner.h b/src/runner.h
index 6990dacd46938c4ff16de0e95e3d74155d0d5f72..0b22f1710301c43ca148a002864d3f35a660f237 100644
--- a/src/runner.h
+++ b/src/runner.h
@@ -98,6 +98,7 @@ void runner_do_all_stars_sort(struct runner *r, struct cell *c);
void runner_do_drift_part(struct runner *r, struct cell *c, int timer);
void runner_do_drift_gpart(struct runner *r, struct cell *c, int timer);
void runner_do_drift_spart(struct runner *r, struct cell *c, int timer);
+void runner_do_drift_sink(struct runner *r, struct cell *c, int timer);
void runner_do_drift_bpart(struct runner *r, struct cell *c, int timer);
void runner_do_kick1(struct runner *r, struct cell *c, int timer);
void runner_do_kick2(struct runner *r, struct cell *c, int timer);
diff --git a/src/runner_drift.c b/src/runner_drift.c
index 8c4376743cd50ffea4709cb471959864cedcc4b7..db1918de8c7012de6780e5e2afa2fd44f1c2dc16 100644
--- a/src/runner_drift.c
+++ b/src/runner_drift.c
@@ -94,3 +94,19 @@ void runner_do_drift_bpart(struct runner *r, struct cell *c, int timer) {
if (timer) TIMER_TOC(timer_drift_bpart);
}
+
+/**
+ * @brief Drift all sink in a cell.
+ *
+ * @param r The runner thread.
+ * @param c The cell.
+ * @param timer Are we timing this ?
+ */
+void runner_do_drift_sink(struct runner *r, struct cell *c, int timer) {
+
+ TIMER_TIC;
+
+ cell_drift_sink(c, r->e, 0);
+
+ if (timer) TIMER_TOC(timer_drift_sink);
+}
diff --git a/src/runner_main.c b/src/runner_main.c
index effcdf3cda94c5b74276af40f69c0e384aecadf4..6be1405e6a5503a0f197a45f489be25c0a1218d0 100644
--- a/src/runner_main.c
+++ b/src/runner_main.c
@@ -340,6 +340,9 @@ void *runner_main(void *data) {
case task_type_drift_spart:
runner_do_drift_spart(r, ci, 1);
break;
+ case task_type_drift_sink:
+ runner_do_drift_sink(r, ci, 1);
+ break;
case task_type_drift_bpart:
runner_do_drift_bpart(r, ci, 1);
break;
diff --git a/src/runner_others.c b/src/runner_others.c
index 4fd380be04bdc23b598bcceae8223107cc9c48ec..738353ab4f733e7120daf7786973cd8124e6cd4e 100644
--- a/src/runner_others.c
+++ b/src/runner_others.c
@@ -570,6 +570,8 @@ void runner_do_end_grav_force(struct runner *r, struct cell *c, int timer) {
id = e->s->parts[-gp->id_or_neg_offset].id;
else if (gp->type == swift_type_stars)
id = e->s->sparts[-gp->id_or_neg_offset].id;
+ else if (gp->type == swift_type_sink)
+ id = e->s->sinks[-gp->id_or_neg_offset].id;
else if (gp->type == swift_type_black_hole)
id = e->s->bparts[-gp->id_or_neg_offset].id;
else
@@ -604,6 +606,8 @@ void runner_do_end_grav_force(struct runner *r, struct cell *c, int timer) {
my_id = e->s->parts[-gp->id_or_neg_offset].id;
else if (gp->type == swift_type_stars)
my_id = e->s->sparts[-gp->id_or_neg_offset].id;
+ else if (gp->type == swift_type_sink)
+ my_id = e->s->sinks[-gp->id_or_neg_offset].id;
else if (gp->type == swift_type_black_hole)
error("Unexisting type");
else
@@ -660,6 +664,9 @@ void runner_do_logger(struct runner *r, struct cell *c, int timer) {
if (c->black_holes.count != 0) {
error("Black holes are not implemented in the logger.");
}
+ if (c->sinks.count != 0) {
+ error("Sink particles are not implemented in the logger.");
+ }
/* Anything to do here? */
if (!cell_is_active_hydro(c, e) && !cell_is_active_gravity(c, e) &&
diff --git a/src/runner_time_integration.c b/src/runner_time_integration.c
index 0e14524aedda6276bd3e280e2c50040e9210a7fb..16d7267fa55ccd8b7444f5c518766a4f237b4dab 100644
--- a/src/runner_time_integration.c
+++ b/src/runner_time_integration.c
@@ -91,10 +91,12 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
struct xpart *restrict xparts = c->hydro.xparts;
struct gpart *restrict gparts = c->grav.parts;
struct spart *restrict sparts = c->stars.parts;
+ struct sink *restrict sinks = c->sinks.parts;
struct bpart *restrict bparts = c->black_holes.parts;
const int count = c->hydro.count;
const int gcount = c->grav.count;
const int scount = c->stars.count;
+ const int sink_count = c->sinks.count;
const int bcount = c->black_holes.count;
const integertime_t ti_current = e->ti_current;
const double time_base = e->time_base;
@@ -103,7 +105,8 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
/* Anything to do here? */
if (!cell_is_starting_hydro(c, e) && !cell_is_starting_gravity(c, e) &&
- !cell_is_starting_stars(c, e) && !cell_is_starting_black_holes(c, e))
+ !cell_is_starting_stars(c, e) && !cell_is_starting_sinks(c, e) &&
+ !cell_is_starting_black_holes(c, e))
return;
/* Recurse? */
@@ -183,10 +186,7 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
/* If the g-particle has no counterpart and needs to be kicked */
if ((gp->type == swift_type_dark_matter ||
- gp->type == swift_type_dark_matter_background ||
- gp->type == swift_type_sink) &&
- // TODO loic remove this
-
+ gp->type == swift_type_dark_matter_background) &&
gpart_is_starting(gp, e)) {
const integertime_t ti_step = get_integer_timestep(gp->time_bin);
@@ -256,6 +256,44 @@ void runner_do_kick1(struct runner *r, struct cell *c, int timer) {
}
}
+ /* Loop over the sink particles in this cell. */
+ for (int k = 0; k < sink_count; k++) {
+
+ /* Get a handle on the s-part. */
+ struct sink *restrict sink = &sinks[k];
+
+ /* If particle needs to be kicked */
+ if (sink_is_starting(sink, e)) {
+
+ const integertime_t ti_step = get_integer_timestep(sink->time_bin);
+ const integertime_t ti_begin =
+ get_integer_time_begin(ti_current + 1, sink->time_bin);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ const integertime_t ti_end =
+ get_integer_time_end(ti_current + 1, sink->time_bin);
+
+ if (ti_begin != ti_current)
+ error(
+ "Sink-particle in wrong time-bin, ti_end=%lld, ti_begin=%lld, "
+ "ti_step=%lld time_bin=%d ti_current=%lld",
+ ti_end, ti_begin, ti_step, sink->time_bin, ti_current);
+#endif
+
+ /* Time interval for this half-kick */
+ double dt_kick_grav;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_begin,
+ ti_begin + ti_step / 2);
+ } else {
+ dt_kick_grav = (ti_step / 2) * time_base;
+ }
+
+ /* do the kick */
+ kick_sink(sink, dt_kick_grav, ti_begin, ti_begin + ti_step / 2);
+ }
+ }
+
/* Loop over the black hole particles in this cell. */
for (int k = 0; k < bcount; k++) {
@@ -317,11 +355,13 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
const int count = c->hydro.count;
const int gcount = c->grav.count;
const int scount = c->stars.count;
+ const int sink_count = c->sinks.count;
const int bcount = c->black_holes.count;
struct part *restrict parts = c->hydro.parts;
struct xpart *restrict xparts = c->hydro.xparts;
struct gpart *restrict gparts = c->grav.parts;
struct spart *restrict sparts = c->stars.parts;
+ struct sink *restrict sinks = c->sinks.parts;
struct bpart *restrict bparts = c->black_holes.parts;
const integertime_t ti_current = e->ti_current;
const double time_base = e->time_base;
@@ -330,7 +370,8 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
/* Anything to do here? */
if (!cell_is_active_hydro(c, e) && !cell_is_active_gravity(c, e) &&
- !cell_is_active_stars(c, e) && !cell_is_active_black_holes(c, e))
+ !cell_is_active_stars(c, e) && !cell_is_active_sinks(c, e) &&
+ !cell_is_active_black_holes(c, e))
return;
/* Recurse? */
@@ -410,10 +451,7 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
/* If the g-particle has no counterpart and needs to be kicked */
if ((gp->type == swift_type_dark_matter ||
- gp->type == swift_type_dark_matter_background ||
- gp->type == swift_type_sink) &&
- // TODO loic remove this
-
+ gp->type == swift_type_dark_matter_background) &&
gpart_is_active(gp, e)) {
const integertime_t ti_step = get_integer_timestep(gp->time_bin);
@@ -491,6 +529,48 @@ void runner_do_kick2(struct runner *r, struct cell *c, int timer) {
}
}
+ /* Loop over the sink particles in this cell. */
+ for (int k = 0; k < sink_count; k++) {
+
+ /* Get a handle on the part. */
+ struct sink *restrict sink = &sinks[k];
+
+ /* If particle needs to be kicked */
+ if (sink_is_active(sink, e)) {
+
+ const integertime_t ti_step = get_integer_timestep(sink->time_bin);
+ const integertime_t ti_begin =
+ get_integer_time_begin(ti_current, sink->time_bin);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (ti_begin + ti_step != ti_current)
+ error("Particle in wrong time-bin");
+#endif
+
+ /* Time interval for this half-kick */
+ double dt_kick_grav;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(
+ cosmo, ti_begin + ti_step / 2, ti_begin + ti_step);
+ } else {
+ dt_kick_grav = (ti_step / 2) * time_base;
+ }
+
+ /* Finish the time-step with a second half-kick */
+ kick_sink(sink, dt_kick_grav, ti_begin + ti_step / 2,
+ ti_begin + ti_step);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Check that kick and the drift are synchronized */
+ if (sink->ti_drift != sink->ti_kick)
+ error("Error integrating sink-part in time.");
+#endif
+
+ /* Prepare the values to be drifted */
+ sink_reset_predicted_values(sink);
+ }
+ }
+
/* Loop over the particles in this cell. */
for (int k = 0; k < bcount; k++) {
@@ -553,32 +633,39 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
const int count = c->hydro.count;
const int gcount = c->grav.count;
const int scount = c->stars.count;
+ const int sink_count = c->sinks.count;
const int bcount = c->black_holes.count;
struct part *restrict parts = c->hydro.parts;
struct xpart *restrict xparts = c->hydro.xparts;
struct gpart *restrict gparts = c->grav.parts;
struct spart *restrict sparts = c->stars.parts;
+ struct sink *restrict sinks = c->sinks.parts;
struct bpart *restrict bparts = c->black_holes.parts;
TIMER_TIC;
/* Anything to do here? */
if (!cell_is_active_hydro(c, e) && !cell_is_active_gravity(c, e) &&
- !cell_is_active_stars(c, e) && !cell_is_active_black_holes(c, e)) {
+ !cell_is_active_stars(c, e) && !cell_is_active_sinks(c, e) &&
+ !cell_is_active_black_holes(c, e)) {
c->hydro.updated = 0;
c->grav.updated = 0;
c->stars.updated = 0;
+ c->sinks.updated = 0;
c->black_holes.updated = 0;
return;
}
- int updated = 0, g_updated = 0, s_updated = 0, b_updated = 0;
+ int updated = 0, g_updated = 0, s_updated = 0, sink_updated = 0,
+ b_updated = 0;
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;
integertime_t ti_stars_end_min = max_nr_timesteps, ti_stars_end_max = 0,
ti_stars_beg_max = 0;
+ integertime_t ti_sinks_end_min = max_nr_timesteps, ti_sinks_end_max = 0,
+ ti_sinks_beg_max = 0;
integertime_t ti_black_holes_end_min = max_nr_timesteps,
ti_black_holes_end_max = 0, ti_black_holes_beg_max = 0;
@@ -678,9 +765,7 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
/* If the g-particle has no counterpart */
if (gp->type == swift_type_dark_matter ||
- gp->type == swift_type_dark_matter_background ||
- gp->type == swift_type_sink) {
- // Loic TODO remove sink
+ gp->type == swift_type_dark_matter_background) {
/* need to be updated ? */
if (gpart_is_active(gp, e)) {
@@ -712,13 +797,6 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
/* What is the next starting point for this cell ? */
ti_gravity_beg_max = max(ti_current, ti_gravity_beg_max);
- // Loic TODO remove this
- /* Synchronize the time step */
- if (gp->type == swift_type_sink) {
- struct sink *sink = &e->s->sinks[-gp->id_or_neg_offset];
- sink->time_bin = gp->time_bin;
- }
-
} else { /* gpart is inactive */
if (!gpart_is_inhibited(gp, e)) {
@@ -805,7 +883,67 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
}
}
- /* Loop over the star particles in this cell. */
+ /* Loop over the sink particles in this cell. */
+ for (int k = 0; k < sink_count; k++) {
+
+ /* Get a handle on the part. */
+ struct sink *restrict sink = &sinks[k];
+
+ /* need to be updated ? */
+ if (sink_is_active(sink, e)) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Current end of time-step */
+ const integertime_t ti_end =
+ get_integer_time_end(ti_current, sink->time_bin);
+
+ if (ti_end != ti_current)
+ error("Computing time-step of rogue particle.");
+#endif
+ /* Get new time-step */
+ const integertime_t ti_new_step = get_sink_timestep(sink, e);
+
+ /* Update particle */
+ sink->time_bin = get_time_bin(ti_new_step);
+ sink->gpart->time_bin = get_time_bin(ti_new_step);
+
+ /* Number of updated sink-particles */
+ sink_updated++;
+ g_updated++;
+
+ ti_sinks_end_min = min(ti_current + ti_new_step, ti_sinks_end_min);
+ ti_sinks_end_max = max(ti_current + ti_new_step, ti_sinks_end_max);
+ ti_gravity_end_min = min(ti_current + ti_new_step, ti_gravity_end_min);
+ ti_gravity_end_max = max(ti_current + ti_new_step, ti_gravity_end_max);
+
+ /* What is the next starting point for this cell ? */
+ ti_sinks_beg_max = max(ti_current, ti_sinks_beg_max);
+ ti_gravity_beg_max = max(ti_current, ti_gravity_beg_max);
+
+ /* sink particle is inactive but not inhibited */
+ } else {
+
+ if (!sink_is_inhibited(sink, e)) {
+
+ const integertime_t ti_end =
+ get_integer_time_end(ti_current, sink->time_bin);
+
+ const integertime_t ti_beg =
+ get_integer_time_begin(ti_current + 1, sink->time_bin);
+
+ ti_sinks_end_min = min(ti_end, ti_sinks_end_min);
+ ti_sinks_end_max = max(ti_end, ti_sinks_end_max);
+ ti_gravity_end_min = min(ti_end, ti_gravity_end_min);
+ ti_gravity_end_max = max(ti_end, ti_gravity_end_max);
+
+ /* What is the next starting point for this cell ? */
+ ti_sinks_beg_max = max(ti_beg, ti_sinks_beg_max);
+ ti_gravity_beg_max = max(ti_beg, ti_gravity_beg_max);
+ }
+ }
+ }
+
+ /* Loop over the black hole particles in this cell. */
for (int k = 0; k < bcount; k++) {
/* Get a handle on the part. */
@@ -880,6 +1018,7 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
/* And aggregate */
updated += cp->hydro.updated;
g_updated += cp->grav.updated;
+ sink_updated += cp->sinks.updated;
s_updated += cp->stars.updated;
b_updated += cp->black_holes.updated;
@@ -895,6 +1034,10 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
ti_stars_end_max = max(cp->stars.ti_end_max, ti_stars_end_max);
ti_stars_beg_max = max(cp->stars.ti_beg_max, ti_stars_beg_max);
+ ti_sinks_end_min = min(cp->sinks.ti_end_min, ti_sinks_end_min);
+ ti_sinks_end_max = max(cp->sinks.ti_end_max, ti_sinks_end_max);
+ ti_sinks_beg_max = max(cp->sinks.ti_beg_max, ti_sinks_beg_max);
+
ti_black_holes_end_min =
min(cp->black_holes.ti_end_min, ti_black_holes_end_min);
ti_black_holes_end_max =
@@ -909,6 +1052,7 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
c->hydro.updated = updated;
c->grav.updated = g_updated;
c->stars.updated = s_updated;
+ c->sinks.updated = sink_updated;
c->black_holes.updated = b_updated;
c->hydro.ti_end_min = ti_hydro_end_min;
@@ -920,6 +1064,9 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
c->stars.ti_end_min = ti_stars_end_min;
c->stars.ti_end_max = ti_stars_end_max;
c->stars.ti_beg_max = ti_stars_beg_max;
+ c->sinks.ti_end_min = ti_sinks_end_min;
+ c->sinks.ti_end_max = ti_sinks_end_max;
+ c->sinks.ti_beg_max = ti_sinks_beg_max;
c->black_holes.ti_end_min = ti_black_holes_end_min;
c->black_holes.ti_end_max = ti_black_holes_end_max;
c->black_holes.ti_beg_max = ti_black_holes_beg_max;
@@ -934,6 +1081,9 @@ void runner_do_timestep(struct runner *r, struct cell *c, int timer) {
if (c->stars.ti_end_min == e->ti_current &&
c->stars.ti_end_min < max_nr_timesteps)
error("End of next stars step is current time!");
+ if (c->sinks.ti_end_min == e->ti_current &&
+ c->sinks.ti_end_min < max_nr_timesteps)
+ error("End of next sinks step is current time!");
if (c->black_holes.ti_end_min == e->ti_current &&
c->black_holes.ti_end_min < max_nr_timesteps)
error("End of next black holes step is current time!");
diff --git a/src/scheduler.c b/src/scheduler.c
index bb380e5a51f478fb53762e805ef1bbbdb0fceec7..165c3d2fe99ce7dbb9734fd9f82d5805a244027e 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -545,6 +545,7 @@ static void scheduler_splittask_hydro(struct task *t, struct scheduler *s) {
access the engine... */
const int with_feedback = (s->space->e->policy & engine_policy_feedback);
const int with_stars = (s->space->e->policy & engine_policy_stars);
+ const int with_sinks = (s->space->e->policy & engine_policy_sinks);
const int with_black_holes =
(s->space->e->policy & engine_policy_black_holes);
@@ -558,6 +559,7 @@ static void scheduler_splittask_hydro(struct task *t, struct scheduler *s) {
const int is_self =
(t->type == task_type_self) && (t->ci != NULL) &&
((t->ci->hydro.count > 0) || (with_stars && t->ci->stars.count > 0) ||
+ (with_sinks && t->ci->sinks.count > 0) ||
(with_black_holes && t->ci->black_holes.count > 0));
/* Is this a non-empty pair-task? */
@@ -565,9 +567,11 @@ static void scheduler_splittask_hydro(struct task *t, struct scheduler *s) {
(t->cj != NULL) &&
((t->ci->hydro.count > 0) ||
(with_feedback && t->ci->stars.count > 0) ||
+ (with_sinks && t->ci->sinks.count > 0) ||
(with_black_holes && t->ci->black_holes.count > 0)) &&
((t->cj->hydro.count > 0) ||
(with_feedback && t->cj->stars.count > 0) ||
+ (with_sinks && t->cj->sinks.count > 0) ||
(with_black_holes && t->cj->black_holes.count > 0));
/* Empty task? */
@@ -1356,6 +1360,7 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
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;
+ const float sink_count_i = (t->ci != NULL) ? t->ci->sinks.count : 0.f;
const float bcount_i = (t->ci != NULL) ? t->ci->black_holes.count : 0.f;
const float bcount_j = (t->cj != NULL) ? t->cj->black_holes.count : 0.f;
@@ -1545,6 +1550,9 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
case task_type_drift_spart:
cost = wscale * scount_i;
break;
+ case task_type_drift_sink:
+ cost = wscale * sink_count_i;
+ break;
case task_type_drift_bpart:
cost = wscale * bcount_i;
break;
@@ -1576,13 +1584,16 @@ void scheduler_reweight(struct scheduler *s, int verbose) {
cost = wscale * (count_i + scount_i);
break;
case task_type_kick1:
- cost = wscale * (count_i + gcount_i + scount_i + bcount_i);
+ cost =
+ wscale * (count_i + gcount_i + scount_i + sink_count_i + bcount_i);
break;
case task_type_kick2:
- cost = wscale * (count_i + gcount_i + scount_i + bcount_i);
+ cost =
+ wscale * (count_i + gcount_i + scount_i + sink_count_i + bcount_i);
break;
case task_type_timestep:
- cost = wscale * (count_i + gcount_i + scount_i + bcount_i);
+ cost =
+ wscale * (count_i + gcount_i + scount_i + sink_count_i + bcount_i);
break;
case task_type_timestep_limiter:
cost = wscale * count_i;
diff --git a/src/space.c b/src/space.c
index da840cd83c437f55778cb3b6666d12c61bd6f888..d9a13d92922aa9b9bb2331f02a319d0db03dc419 100644
--- a/src/space.c
+++ b/src/space.c
@@ -206,6 +206,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->grav.mm = NULL;
c->hydro.dx_max_part = 0.0f;
c->hydro.dx_max_sort = 0.0f;
+ c->sinks.dx_max_part = 0.f;
c->stars.dx_max_part = 0.f;
c->stars.dx_max_sort = 0.f;
c->black_holes.dx_max_part = 0.f;
@@ -252,6 +253,7 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->timestep_sync = NULL;
c->hydro.end_force = NULL;
c->hydro.drift = NULL;
+ c->sinks.drift = NULL;
c->stars.drift = NULL;
c->stars.stars_in = NULL;
c->stars.stars_out = NULL;
@@ -285,6 +287,8 @@ void space_rebuild_recycle_mapper(void *map_data, int num_elements,
c->hydro.ti_end_max = -1;
c->grav.ti_end_min = -1;
c->grav.ti_end_max = -1;
+ c->sinks.ti_end_min = -1;
+ c->sinks.ti_end_max = -1;
c->stars.ti_end_min = -1;
c->stars.ti_end_max = -1;
c->black_holes.ti_end_min = -1;
@@ -621,6 +625,7 @@ void space_regrid(struct space *s, int verbose) {
c->hydro.count = 0;
c->grav.count = 0;
c->stars.count = 0;
+ c->sinks.count = 0;
c->top = c;
c->super = c;
c->hydro.super = c;
@@ -628,6 +633,7 @@ void space_regrid(struct space *s, int verbose) {
c->hydro.ti_old_part = ti_current;
c->grav.ti_old_part = ti_current;
c->stars.ti_old_part = ti_current;
+ c->sinks.ti_old_part = ti_current;
c->black_holes.ti_old_part = ti_current;
c->grav.ti_old_multipole = ti_current;
#ifdef WITH_MPI
@@ -1613,6 +1619,8 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
s->parts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
} else if (s->gparts[k].type == swift_type_stars) {
s->sparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
+ } else if (s->gparts[k].type == swift_type_sink) {
+ s->sparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
} else if (s->gparts[k].type == swift_type_black_hole) {
s->bparts[-s->gparts[k].id_or_neg_offset].gpart = &s->gparts[k];
}
@@ -1623,6 +1631,9 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
} else if (s->gparts[nr_gparts].type == swift_type_stars) {
s->sparts[-s->gparts[nr_gparts].id_or_neg_offset].gpart =
&s->gparts[nr_gparts];
+ } else if (s->gparts[nr_gparts].type == swift_type_sink) {
+ s->sparts[-s->gparts[nr_gparts].id_or_neg_offset].gpart =
+ &s->gparts[nr_gparts];
} else if (s->gparts[nr_gparts].type == swift_type_black_hole) {
s->bparts[-s->gparts[nr_gparts].id_or_neg_offset].gpart =
&s->gparts[nr_gparts];
@@ -2087,6 +2098,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
c->grav.ti_old_part = ti_current;
c->grav.ti_old_multipole = ti_current;
c->stars.ti_old_part = ti_current;
+ c->sinks.ti_old_part = ti_current;
c->black_holes.ti_old_part = ti_current;
#if defined(SWIFT_DEBUG_CHECKS) || defined(SWIFT_CELL_GRAPH)
@@ -2114,6 +2126,7 @@ void space_rebuild(struct space *s, int repartitioned, int verbose) {
c->hydro.count_total = c->hydro.count + space_extra_parts;
c->grav.count_total = c->grav.count + space_extra_gparts;
c->stars.count_total = c->stars.count + space_extra_sparts;
+ c->sinks.count_total = c->sinks.count + space_extra_sinks;
c->black_holes.count_total = c->black_holes.count + space_extra_bparts;
finger = &finger[c->hydro.count_total];
@@ -2991,7 +3004,7 @@ void space_sinks_get_cell_index_mapper(void *map_data, int nr_sinks,
if (cell_counts[k]) atomic_add(&data->cell_counts[k], cell_counts[k]);
free(cell_counts);
- /* Write the count of inhibited and extra bparts */
+ /* Write the count of inhibited and extra sinks */
if (count_inhibited_sink)
atomic_add(&data->count_inhibited_sink, count_inhibited_sink);
if (count_extra_sink) atomic_add(&data->count_extra_sink, count_extra_sink);
@@ -3167,6 +3180,10 @@ void space_sinks_get_cell_index(struct space *s, int *sink_ind,
const ticks tic = getticks();
+ /* Re-set the counters */
+ s->min_sink_mass = FLT_MAX;
+ s->sum_sink_vel_norm = 0.f;
+
/* Pack the extra information */
struct index_data data;
data.s = s;
@@ -3757,6 +3774,8 @@ void space_split_recursive(struct space *s, struct cell *c,
ti_gravity_beg_max = 0;
integertime_t ti_stars_end_min = max_nr_timesteps, ti_stars_end_max = 0,
ti_stars_beg_max = 0;
+ integertime_t ti_sinks_end_min = max_nr_timesteps, ti_sinks_end_max = 0,
+ ti_sinks_beg_max = 0;
integertime_t ti_black_holes_end_min = max_nr_timesteps,
ti_black_holes_end_max = 0, ti_black_holes_beg_max = 0;
struct part *parts = c->hydro.parts;
@@ -3895,6 +3914,7 @@ void space_split_recursive(struct space *s, struct cell *c,
cp->grav.ti_old_part = c->grav.ti_old_part;
cp->grav.ti_old_multipole = c->grav.ti_old_multipole;
cp->stars.ti_old_part = c->stars.ti_old_part;
+ cp->sinks.ti_old_part = c->sinks.ti_old_part;
cp->black_holes.ti_old_part = c->black_holes.ti_old_part;
cp->loc[0] = c->loc[0];
cp->loc[1] = c->loc[1];
@@ -3914,6 +3934,7 @@ void space_split_recursive(struct space *s, struct cell *c,
cp->stars.h_max = 0.f;
cp->stars.dx_max_part = 0.f;
cp->stars.dx_max_sort = 0.f;
+ cp->sinks.dx_max_part = 0.f;
cp->black_holes.h_max = 0.f;
cp->black_holes.dx_max_part = 0.f;
cp->nodeID = c->nodeID;
@@ -3981,6 +4002,9 @@ void space_split_recursive(struct space *s, struct cell *c,
ti_stars_end_min = min(ti_stars_end_min, cp->stars.ti_end_min);
ti_stars_end_max = max(ti_stars_end_max, cp->stars.ti_end_max);
ti_stars_beg_max = max(ti_stars_beg_max, cp->stars.ti_beg_max);
+ ti_sinks_end_min = min(ti_sinks_end_min, cp->sinks.ti_end_min);
+ ti_sinks_end_max = max(ti_sinks_end_max, cp->sinks.ti_end_max);
+ ti_sinks_beg_max = max(ti_sinks_beg_max, cp->sinks.ti_beg_max);
ti_black_holes_end_min =
min(ti_black_holes_end_min, cp->black_holes.ti_end_min);
ti_black_holes_end_max =
@@ -4204,6 +4228,30 @@ void space_split_recursive(struct space *s, struct cell *c,
sparts[k].x_diff[2] = 0.f;
}
+ /* sinks: Get dt_min/dt_max */
+ for (int k = 0; k < sink_count; k++) {
+#ifdef SWIFT_DEBUG_CHECKS
+ if (sinks[k].time_bin == time_bin_not_created)
+ error("Extra sink-particle present in space_split()");
+ if (sinks[k].time_bin == time_bin_inhibited)
+ error("Inhibited sink-particle present in space_split()");
+#endif
+
+ /* When does this particle's time-step start and end? */
+ const timebin_t time_bin = sinks[k].time_bin;
+ const integertime_t ti_end = get_integer_time_end(ti_current, time_bin);
+ const integertime_t ti_beg = get_integer_time_begin(ti_current, time_bin);
+
+ ti_sinks_end_min = min(ti_sinks_end_min, ti_end);
+ ti_sinks_end_max = max(ti_sinks_end_max, ti_end);
+ ti_sinks_beg_max = max(ti_sinks_beg_max, ti_beg);
+
+ /* Reset x_diff */
+ sinks[k].x_diff[0] = 0.f;
+ sinks[k].x_diff[1] = 0.f;
+ sinks[k].x_diff[2] = 0.f;
+ }
+
/* bparts: Get dt_min/dt_max */
for (int k = 0; k < bcount; k++) {
#ifdef SWIFT_DEBUG_CHECKS
@@ -4274,6 +4322,9 @@ void space_split_recursive(struct space *s, struct cell *c,
c->stars.ti_end_max = ti_stars_end_max;
c->stars.ti_beg_max = ti_stars_beg_max;
c->stars.h_max = stars_h_max;
+ c->sinks.ti_end_min = ti_sinks_end_min;
+ c->sinks.ti_end_max = ti_sinks_end_max;
+ c->sinks.ti_beg_max = ti_sinks_beg_max;
c->black_holes.ti_end_min = ti_black_holes_end_min;
c->black_holes.ti_end_max = ti_black_holes_end_max;
c->black_holes.ti_beg_max = ti_black_holes_beg_max;
@@ -5460,10 +5511,12 @@ void space_init(struct space *s, struct swift_params *params,
s->sinks = sinks;
s->min_part_mass = FLT_MAX;
s->min_gpart_mass = FLT_MAX;
+ s->min_sink_mass = FLT_MAX;
s->min_spart_mass = FLT_MAX;
s->min_bpart_mass = FLT_MAX;
s->sum_part_vel_norm = 0.f;
s->sum_gpart_vel_norm = 0.f;
+ s->sum_sink_vel_norm = 0.f;
s->sum_spart_vel_norm = 0.f;
s->sum_bpart_vel_norm = 0.f;
s->nr_queues = 1; /* Temporary value until engine construction */
diff --git a/src/space.h b/src/space.h
index f597db6421fe39db1fad17e0201597efdba060f8..3b942a03d653f2b0e7d6255698b2e8c1439f5bf2 100644
--- a/src/space.h
+++ b/src/space.h
@@ -262,6 +262,9 @@ struct space {
/*! Minimal mass of all the #spart */
float min_spart_mass;
+ /*! Minimal mass of all the #sink */
+ float min_sink_mass;
+
/*! Minimal mass of all the #bpart */
float min_bpart_mass;
@@ -274,6 +277,9 @@ struct space {
/*! Sum of the norm of the velocity of all the #spart */
float sum_spart_vel_norm;
+ /*! Sum of the norm of the velocity of all the #sink */
+ float sum_sink_vel_norm;
+
/*! Sum of the norm of the velocity of all the #bpart */
float sum_bpart_vel_norm;
diff --git a/src/task.c b/src/task.c
index a766d5b36392df4b3a0a472be0e1dcee77cf5a28..e978d6f23e0526524dd83e3a22b04b9302e6a603 100644
--- a/src/task.c
+++ b/src/task.c
@@ -63,6 +63,7 @@ const char *taskID_names[task_type_count] = {"none",
"extra_ghost",
"drift_part",
"drift_spart",
+ "drift_sink",
"drift_bpart",
"drift_gpart",
"drift_gpart_out",
@@ -104,15 +105,37 @@ const char *taskID_names[task_type_count] = {"none",
"fof_pair"};
/* Sub-task type names. */
-const char *subtaskID_names[task_subtype_count] = {
- "none", "density", "gradient", "force",
- "limiter", "grav", "external_grav", "tend_part",
- "tend_gpart", "tend_spart", "tend_bpart", "xv",
- "rho", "part_swallow", "bpart_merger", "gpart",
- "multipole", "spart", "stars_density", "stars_feedback",
- "sf_count", "bpart_rho", "bpart_swallow", "bpart_feedback",
- "bh_density", "bh_swallow", "do_gas_swallow", "do_bh_swallow",
- "bh_feedback"};
+const char *subtaskID_names[task_subtype_count] = {"none",
+ "density",
+ "gradient",
+ "force",
+ "limiter",
+ "grav",
+ "external_grav",
+ "tend_part",
+ "tend_gpart",
+ "tend_spart",
+ "tend_sink",
+ "tend_bpart",
+ "xv",
+ "rho",
+ "part_swallow",
+ "bpart_merger",
+ "gpart",
+ "multipole",
+ "spart",
+ "stars_density",
+ "stars_feedback",
+ "sf_count",
+ "bpart_rho",
+ "bpart_swallow",
+ "bpart_feedback",
+ "bh_density",
+ "bh_swallow",
+ "do_gas_swallow",
+ "do_bh_swallow",
+ "bh_feedback",
+ "sink"};
const char *task_category_names[task_category_count] = {
"drift", "sort", "hydro", "gravity", "feedback",
@@ -152,6 +175,7 @@ MPI_Comm subtaskMPI_comms[task_subtype_count];
TASK_CELL_OVERLAP(part, hydro.parts, hydro.count);
TASK_CELL_OVERLAP(gpart, grav.parts, grav.count);
TASK_CELL_OVERLAP(spart, stars.parts, stars.count);
+TASK_CELL_OVERLAP(sink, sinks.parts, sinks.count);
TASK_CELL_OVERLAP(bpart, black_holes.parts, black_holes.count);
/**
@@ -187,6 +211,10 @@ __attribute__((always_inline)) INLINE static enum task_actions task_acts_on(
return task_action_spart;
break;
+ case task_type_drift_sink:
+ return task_action_sink;
+ break;
+
case task_type_drift_bpart:
case task_type_bh_density_ghost:
case task_type_bh_swallow_ghost3:
@@ -314,6 +342,7 @@ float task_overlap(const struct task *restrict ta,
(ta_act == task_action_gpart || ta_act == task_action_all);
const int ta_spart =
(ta_act == task_action_spart || ta_act == task_action_all);
+ const int ta_sink = (ta_act == task_action_sink || ta_act == task_action_all);
const int ta_bpart =
(ta_act == task_action_bpart || ta_act == task_action_all);
const int tb_part = (tb_act == task_action_part || tb_act == task_action_all);
@@ -321,6 +350,7 @@ float task_overlap(const struct task *restrict ta,
(tb_act == task_action_gpart || tb_act == task_action_all);
const int tb_spart =
(tb_act == task_action_spart || tb_act == task_action_all);
+ const int tb_sink = (tb_act == task_action_sink || tb_act == task_action_all);
const int tb_bpart =
(tb_act == task_action_bpart || tb_act == task_action_all);
@@ -387,6 +417,27 @@ float task_overlap(const struct task *restrict ta,
return ((float)size_intersect) / (size_union - size_intersect);
}
+ /* In the case where both tasks act on sink */
+ else if (ta_sink && tb_sink) {
+
+ /* Compute the union of the cell data. */
+ size_t size_union = 0;
+ if (ta->ci != NULL) size_union += ta->ci->sinks.count;
+ if (ta->cj != NULL) size_union += ta->cj->sinks.count;
+ if (tb->ci != NULL) size_union += tb->ci->sinks.count;
+ if (tb->cj != NULL) size_union += tb->cj->sinks.count;
+
+ if (size_union == 0) return 0.f;
+
+ /* Compute the intersection of the cell data. */
+ const size_t size_intersect = task_cell_overlap_spart(ta->ci, tb->ci) +
+ task_cell_overlap_sink(ta->ci, tb->cj) +
+ task_cell_overlap_sink(ta->cj, tb->ci) +
+ task_cell_overlap_sink(ta->cj, tb->cj);
+
+ return ((float)size_intersect) / (size_union - size_intersect);
+ }
+
/* In the case where both tasks act on bparts */
else if (ta_bpart && tb_bpart) {
@@ -1354,6 +1405,7 @@ enum task_categories task_get_category(const struct task *t) {
case task_type_drift_part:
case task_type_drift_spart:
+ case task_type_drift_sink:
case task_type_drift_bpart:
case task_type_drift_gpart:
return task_category_drift;
diff --git a/src/task.h b/src/task.h
index c8b7e587c8afe24875fb4568de6d8963920f680c..2fdc4ba2b1774256eca225c0e1997a66dc332109 100644
--- a/src/task.h
+++ b/src/task.h
@@ -57,6 +57,7 @@ enum task_types {
task_type_extra_ghost,
task_type_drift_part,
task_type_drift_spart,
+ task_type_drift_sink,
task_type_drift_bpart,
task_type_drift_gpart,
task_type_drift_gpart_out, /* Implicit */
@@ -113,6 +114,7 @@ enum task_subtypes {
task_subtype_tend_part,
task_subtype_tend_gpart,
task_subtype_tend_spart,
+ task_subtype_tend_sink,
task_subtype_tend_bpart,
task_subtype_xv,
task_subtype_rho,
@@ -132,6 +134,7 @@ enum task_subtypes {
task_subtype_do_gas_swallow,
task_subtype_do_bh_swallow,
task_subtype_bh_feedback,
+ task_subtype_sink,
task_subtype_count
} __attribute__((packed));
@@ -143,6 +146,7 @@ enum task_actions {
task_action_part,
task_action_gpart,
task_action_spart,
+ task_action_sink,
task_action_bpart,
task_action_all,
task_action_multipole,
diff --git a/src/timers.c b/src/timers.c
index cc4ebc90d491832178f1e48e96960f489c73f4c1..23f2b9e8b69d201104d0f6354990db83488e9d4b 100644
--- a/src/timers.c
+++ b/src/timers.c
@@ -123,6 +123,7 @@ const char* timers_names[timer_count] = {
"do_stars_resort",
"fof_self",
"fof_pair",
+ "drift_sink",
};
/* File to store the timers */
diff --git a/src/timers.h b/src/timers.h
index 96742567f549c156944333fafd4342dba14a4bee..68cea16ff9629e95760325d3b333745b732f27d5 100644
--- a/src/timers.h
+++ b/src/timers.h
@@ -124,6 +124,7 @@ enum {
timer_do_stars_resort,
timer_fof_self,
timer_fof_pair,
+ timer_drift_sink,
timer_count,
};