Add sink module

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
 

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

src/sink.h

+/*******************************************************************************
+ * 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

src/sink/Default/sink.h

+/*******************************************************************************
+ * 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 */

src/sink/Default/sink_io.h

+/*******************************************************************************
+ * 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 */

src/sink/Default/sink_part.h

+/*******************************************************************************
+ * 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 */

src/sink_io.h

+/*******************************************************************************
+ * 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 */