diff --git a/configure.ac b/configure.ac
index df1cec11823702f98675385c01de80a7c5c1980f..ade84b9410992c77af1a9f377b22fdd013b2a98b 100644
--- a/configure.ac
+++ b/configure.ac
@@ -1476,6 +1476,7 @@ with_subgrid_stars=none
with_subgrid_star_formation=none
with_subgrid_feedback=none
with_subgrid_task_order=default
+with_subgrid_sink=none
case "$with_subgrid" in
yes)
@@ -1492,6 +1493,7 @@ case "$with_subgrid" in
with_subgrid_feedback=GEAR
with_subgrid_black_holes=none
with_subgrid_task_order=GEAR
+ with_subgrid_sink=none
enable_fof=no
;;
QLA)
@@ -1504,6 +1506,7 @@ case "$with_subgrid" in
with_subgrid_feedback=none
with_subgrid_black_holes=none
with_subgrid_task_order=default
+ with_subgrid_sink=none
enable_fof=no
;;
EAGLE)
@@ -1516,6 +1519,7 @@ case "$with_subgrid" in
with_subgrid_feedback=EAGLE
with_subgrid_black_holes=EAGLE
with_subgrid_task_order=EAGLE
+ with_subgrid_sink=none
enable_fof=yes
;;
*)
@@ -2021,6 +2025,32 @@ case "$with_black_holes" in
;;
esac
+# sink model.
+AC_ARG_WITH([sink],
+[AS_HELP_STRING([--with-sink=<model>],
+[Sink particle model to use @<:@none, default: none@:>@]
+)],
+[with_sink="$withval"],
+[with_sink="none"]
+)
+
+if test "$with_subgrid" != "none"; then
+if test "$with_sink" != "none"; then
+AC_MSG_ERROR([Cannot provide with-subgrid and with-sink together])
+else
+with_sink="$with_subgrid_sink"
+fi
+fi
+
+case "$with_sink" in
+none)
+AC_DEFINE([SINK_NONE], [1], [No sink particle model])
+;;
+*)
+AC_MSG_ERROR([Unknown sink particle model model: $with_sink])
+;;
+esac
+
# Task order
AC_ARG_WITH([task-order],
[AS_HELP_STRING([--with-task-order=<model>],
@@ -2314,6 +2344,7 @@ AC_MSG_RESULT([
Stellar model : $with_stars
Star formation model : $with_star_formation
Star feedback model : $with_feedback
+ Sink particle model : $with_sink
Black holes model : $with_black_holes
Task dependencies : $with_task_order
diff --git a/src/Makefile.am b/src/Makefile.am
index b5390e12166b5fa500ca78fd5226c61c80926676..b3b095f089d14ffedf4f298cc11d61118b6a72bd 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -279,7 +279,9 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
pressure_floor/GEAR/pressure_floor.h pressure_floor/none/pressure_floor.h \
pressure_floor/GEAR/pressure_floor_iact.h pressure_floor/none/pressure_floor_iact.h \
pressure_floor/GEAR/pressure_floor_struct.h pressure_floor/none/pressure_floor_struct.h \
- task_order/GEAR/task_order.h task_order/EAGLE/task_order.h task_order/default/task_order.h
+ task_order/GEAR/task_order.h task_order/EAGLE/task_order.h task_order/default/task_order.h \
+ sink/Default/sink.h sink/Default/sink_io.h sink/Default/sink_part.h \
+ sink.h sink_io.h
# Sources and special flags for the gravity library
diff --git a/src/sink.h b/src/sink.h
new file mode 100644
index 0000000000000000000000000000000000000000..991bfd27c0b9f829359b90354864acab1d26ea23
--- /dev/null
+++ b/src/sink.h
@@ -0,0 +1,32 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2020 Loic Hausammann (loic.hausammann@epfl.ch)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_SINK_H
+#define SWIFT_SINK_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Select the correct sink model */
+#if defined(SINK_NONE)
+#include "./sink/Default/sink.h"
+#else
+#error "Invalid choice of sink model"
+#endif
+
+#endif
diff --git a/src/sink/Default/sink.h b/src/sink/Default/sink.h
new file mode 100644
index 0000000000000000000000000000000000000000..b2c9e2a240443e08f4dfa72903a50d33c6023b45
--- /dev/null
+++ b/src/sink/Default/sink.h
@@ -0,0 +1,91 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2020 Loic Hausammann (loic.hausammann@epfl.ch)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_DEFAULT_SINK_H
+#define SWIFT_DEFAULT_SINK_H
+
+#include <float.h>
+
+/* Local includes */
+#include "dimension.h"
+#include "kernel_hydro.h"
+#include "minmax.h"
+#include "sink_part.h"
+
+/**
+ * @brief Computes the time-step of a given sink particle.
+ *
+ * @param sp Pointer to the sink-particle data.
+ */
+__attribute__((always_inline)) INLINE static float sink_compute_timestep(
+ const struct sink* const sp) {
+
+ return FLT_MAX;
+}
+
+/**
+ * @brief Initialises the sink-particles for the first time
+ *
+ * This function is called only once just after the ICs have been
+ * read in to do some conversions.
+ *
+ * @param sp The particle to act upon
+ * @param props The properties of the sink model.
+ */
+__attribute__((always_inline)) INLINE static void sink_first_init_sink(
+ struct sink* sp) {
+
+ sp->time_bin = 0;
+}
+
+/**
+ * @brief Prepares a sink-particle for its interactions
+ *
+ * @param sp The particle to act upon
+ */
+__attribute__((always_inline)) INLINE static void sink_init_sink(
+ struct sink* sp) {}
+
+/**
+ * @brief Predict additional particle fields forward in time when drifting
+ *
+ * @param sp The particle
+ * @param dt_drift The drift time-step for positions.
+ */
+__attribute__((always_inline)) INLINE static void sink_predict_extra(
+ struct sink* restrict sp, float dt_drift) {}
+
+/**
+ * @brief Sets the values to be predicted in the drifts to their values at a
+ * kick time
+ *
+ * @param sp The particle.
+ */
+__attribute__((always_inline)) INLINE static void sink_reset_predicted_values(
+ struct sink* restrict sp) {}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * @param sp The particle to act upon
+ * @param dt The time-step for this kick
+ */
+__attribute__((always_inline)) INLINE static void sink_kick_extra(
+ struct sink* sp, float dt) {}
+
+#endif /* SWIFT_DEFAULT_SINK_H */
diff --git a/src/sink/Default/sink_io.h b/src/sink/Default/sink_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..a1d4a8f8f20efdaa537a0301eb8d410094ee75be
--- /dev/null
+++ b/src/sink/Default/sink_io.h
@@ -0,0 +1,130 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2020 Loic Hausammann (loic.hausammann@epfl.ch)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_DEFAULT_SINK_IO_H
+#define SWIFT_DEFAULT_SINK_IO_H
+
+#include "io_properties.h"
+#include "sink_part.h"
+
+/**
+ * @brief Specifies which sink-particle fields to read from a dataset
+ *
+ * @param sinks The sink-particle array.
+ * @param list The list of i/o properties to read.
+ * @param num_fields The number of i/o fields to read.
+ */
+INLINE static void sink_read_particles(struct sink* sinks,
+ struct io_props* list, int* num_fields) {
+
+ /* Say how much we want to read */
+ *num_fields = 4;
+
+ /* List what we want to read */
+ list[0] = io_make_input_field("Coordinates", DOUBLE, 3, COMPULSORY,
+ UNIT_CONV_LENGTH, sinks, x);
+ list[1] = io_make_input_field("Velocities", FLOAT, 3, COMPULSORY,
+ UNIT_CONV_SPEED, sinks, v);
+ list[2] = io_make_input_field("Masses", FLOAT, 1, COMPULSORY, UNIT_CONV_MASS,
+ sinks, mass);
+ list[3] = io_make_input_field("ParticleIDs", LONGLONG, 1, COMPULSORY,
+ UNIT_CONV_NO_UNITS, sinks, id);
+}
+
+INLINE static void convert_sink_pos(const struct engine* e,
+ const struct sink* sp, double* ret) {
+
+ const struct space* s = e->s;
+ if (s->periodic) {
+ ret[0] = box_wrap(sp->x[0] - s->pos_dithering[0], 0.0, s->dim[0]);
+ ret[1] = box_wrap(sp->x[1] - s->pos_dithering[1], 0.0, s->dim[1]);
+ ret[2] = box_wrap(sp->x[2] - s->pos_dithering[2], 0.0, s->dim[2]);
+ } else {
+ ret[0] = sp->x[0];
+ ret[1] = sp->x[1];
+ ret[2] = sp->x[2];
+ }
+}
+
+INLINE static void convert_sink_vel(const struct engine* e,
+ const struct sink* sp, float* ret) {
+
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const struct cosmology* cosmo = e->cosmology;
+ const integertime_t ti_current = e->ti_current;
+ const double time_base = e->time_base;
+
+ const integertime_t ti_beg = get_integer_time_begin(ti_current, sp->time_bin);
+ const integertime_t ti_end = get_integer_time_end(ti_current, sp->time_bin);
+
+ /* Get time-step since the last kick */
+ float dt_kick_grav;
+ if (with_cosmology) {
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, ti_beg, ti_current);
+ dt_kick_grav -=
+ cosmology_get_grav_kick_factor(cosmo, ti_beg, (ti_beg + ti_end) / 2);
+ } else {
+ dt_kick_grav = (ti_current - ((ti_beg + ti_end) / 2)) * time_base;
+ }
+
+ /* Extrapolate the velocites to the current time */
+ const struct gpart* gp = sp->gpart;
+ ret[0] = gp->v_full[0] + gp->a_grav[0] * dt_kick_grav;
+ ret[1] = gp->v_full[1] + gp->a_grav[1] * dt_kick_grav;
+ ret[2] = gp->v_full[2] + gp->a_grav[2] * dt_kick_grav;
+
+ /* Conversion from internal units to peculiar velocities */
+ ret[0] *= cosmo->a_inv;
+ ret[1] *= cosmo->a_inv;
+ ret[2] *= cosmo->a_inv;
+}
+
+/**
+ * @brief Specifies which sink-particle fields to write to a dataset
+ *
+ * @param sinks The sink-particle array.
+ * @param list The list of i/o properties to write.
+ * @param num_fields The number of i/o fields to write.
+ * @param with_cosmology Are we running a cosmological simulation?
+ */
+INLINE static void sink_write_particles(const struct sink* sinks,
+ struct io_props* list, int* num_fields,
+ int with_cosmology) {
+
+ /* Say how much we want to write */
+ *num_fields = 4;
+
+ /* List what we want to write */
+ list[0] = io_make_output_field_convert_sink(
+ "Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH, 1.f, sinks, convert_sink_pos,
+ "Co-moving position of the particles");
+
+ list[1] = io_make_output_field_convert_sink(
+ "Velocities", FLOAT, 3, UNIT_CONV_SPEED, 0.f, sinks, convert_sink_vel,
+ "Peculiar velocities of the particles. This is a * dx/dt where x is the "
+ "co-moving position of the particles.");
+
+ list[2] = io_make_output_field("Masses", FLOAT, 1, UNIT_CONV_MASS, 0.f, sinks,
+ mass, "Masses of the particles");
+
+ list[3] =
+ io_make_output_field("ParticleIDs", ULONGLONG, 1, UNIT_CONV_NO_UNITS, 0.f,
+ sinks, id, "Unique ID of the particles");
+}
+
+#endif /* SWIFT_DEFAULT_SINK_IO_H */
diff --git a/src/sink/Default/sink_part.h b/src/sink/Default/sink_part.h
new file mode 100644
index 0000000000000000000000000000000000000000..a969d2f44ed9e06121d214a6cd8bdf1ed3f3e425
--- /dev/null
+++ b/src/sink/Default/sink_part.h
@@ -0,0 +1,64 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2020 Loic Hausammann (loic.hausammann@epfl.ch)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_DEFAULT_SINK_PART_H
+#define SWIFT_DEFAULT_SINK_PART_H
+
+#include "timeline.h"
+
+/**
+ * @brief Particle fields for the sink particles.
+ *
+ * All quantities related to gravity are stored in the associate #gpart.
+ */
+struct sink {
+
+ /*! Particle ID. */
+ long long id;
+
+ /*! Pointer to corresponding gravity part. */
+ struct gpart* gpart;
+
+ /*! Particle position. */
+ double x[3];
+
+ /* Offset between current position and position at last tree rebuild. */
+ float x_diff[3];
+
+ /*! Particle velocity. */
+ float v[3];
+
+ /*! Sink particle mass */
+ float mass;
+
+ /*! Particle time bin */
+ timebin_t time_bin;
+
+#ifdef SWIFT_DEBUG_CHECKS
+
+ /* Time of the last drift */
+ integertime_t ti_drift;
+
+ /* Time of the last kick */
+ integertime_t ti_kick;
+
+#endif
+
+} SWIFT_STRUCT_ALIGN;
+
+#endif /* SWIFT_DEFAULT_SINK_PART_H */
diff --git a/src/sink_io.h b/src/sink_io.h
new file mode 100644
index 0000000000000000000000000000000000000000..0f3be8c7278dd3f2ab0f344e46ba0db5af1a442b
--- /dev/null
+++ b/src/sink_io.h
@@ -0,0 +1,31 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Coypright (c) 2020 Loic Hausammann (loic.hausammann@epfl.ch)
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_SINK_IO_H
+#define SWIFT_SINK_IO_H
+
+#include "../config.h"
+
+/* Load the correct sink type */
+#if defined(SINK_NONE)
+#include "./sink/Default/sink_io.h"
+#else
+#error "Invalid choice of sink model"
+#endif
+
+#endif /* SWIFT_SINK_IO_H */