Merge branch 'sink_prepare_drift2' into 'master'

Add drift/kick/timestep function for sink

See merge request !1121

src/drift.h

@@ -31,6 +31,7 @@
 #include "hydro.h"
 #include "hydro_properties.h"
 #include "part.h"
+#include "sink.h"
 #include "stars.h"
 
 /**
@@ -199,4 +200,42 @@ __attribute__((always_inline)) INLINE static void drift_bpart(
   }
 }
 
+/**
+ * @brief Perform the 'drift' operation on a #sink
+ *
+ * @param sink The #sink to drift.
+ * @param dt_drift The drift time-step.
+ * @param ti_old Integer start of time-step (for debugging checks).
+ * @param ti_current Integer end of time-step (for debugging checks).
+ */
+__attribute__((always_inline)) INLINE static void drift_sink(
+    struct sink *restrict sink, double dt_drift, integertime_t ti_old,
+    integertime_t ti_current) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (sink->ti_drift != ti_old)
+    error(
+        "s-particle has not been drifted to the current time "
+        "sink->ti_drift=%lld, "
+        "c->ti_old=%lld, ti_current=%lld",
+        sink->ti_drift, ti_old, ti_current);
+
+  sink->ti_drift = ti_current;
+#endif
+
+  /* Drift... */
+  sink->x[0] += sink->v[0] * dt_drift;
+  sink->x[1] += sink->v[1] * dt_drift;
+  sink->x[2] += sink->v[2] * dt_drift;
+
+  /* Predict the values of the extra fields */
+  sink_predict_extra(sink, dt_drift);
+
+  /* Compute offsets since last cell construction */
+  for (int k = 0; k < 3; k++) {
+    const float dx = sink->v[k] * dt_drift;
+    sink->x_diff[k] -= dx;
+  }
+}
+
 #endif /* SWIFT_DRIFT_H */

src/kick.h

@@ -26,6 +26,7 @@
 #include "black_holes.h"
 #include "const.h"
 #include "debug.h"
+#include "sink.h"
 #include "stars.h"
 #include "timeline.h"
 
@@ -194,4 +195,42 @@ __attribute__((always_inline)) INLINE static void kick_bpart(
   black_holes_kick_extra(bp, dt_kick_grav);
 }
 
+/**
+ * @brief Perform the 'kick' operation on a #sink
+ *
+ * @param sink The #sink to kick.
+ * @param dt_kick_grav The kick time-step for gravity accelerations.
+ * @param ti_start The starting (integer) time of the kick (for debugging
+ * checks).
+ * @param ti_end The ending (integer) time of the kick (for debugging checks).
+ */
+__attribute__((always_inline)) INLINE static void kick_sink(
+    struct sink *restrict sink, double dt_kick_grav, integertime_t ti_start,
+    integertime_t ti_end) {
+
+#ifdef SWIFT_DEBUG_CHECKS
+  if (sink->ti_kick != ti_start)
+    error(
+        "sink-particle has not been kicked to the current time "
+        "sink->ti_kick=%lld, "
+        "ti_start=%lld, ti_end=%lld id=%lld",
+        sink->ti_kick, ti_start, ti_end, sink->id);
+
+  sink->ti_kick = ti_end;
+#endif
+
+  /* Kick particles in momentum space */
+  sink->v[0] += sink->gpart->a_grav[0] * dt_kick_grav;
+  sink->v[1] += sink->gpart->a_grav[1] * dt_kick_grav;
+  sink->v[2] += sink->gpart->a_grav[2] * dt_kick_grav;
+
+  /* Give the gpart friend the same velocity */
+  sink->gpart->v_full[0] = sink->v[0];
+  sink->gpart->v_full[1] = sink->v[1];
+  sink->gpart->v_full[2] = sink->v[2];
+
+  /* Kick extra variables */
+  sink_kick_extra(sink, dt_kick_grav);
+}
+
 #endif /* SWIFT_KICK_H */

src/minmax.h

@@ -103,4 +103,40 @@
     max(_temp1, _temp2);                      \
   })
 
+/**
+ * @brief Minimum of five numbers
+ *
+ * This macro evaluates its arguments exactly once.
+ */
+#define min5(x, y, z, w, u)                           \
+  ({                                                  \
+    const __typeof__(x) _x = (x);                     \
+    const __typeof__(y) _y = (y);                     \
+    const __typeof__(z) _z = (z);                     \
+    const __typeof__(w) _w = (w);                     \
+    const __typeof__(u) _u = (u);                     \
+    const __typeof__(x) _temp1 = min(_x, _y);         \
+    const __typeof__(x) _temp2 = min(_z, _w);         \
+    const __typeof__(x) _temp3 = min(_temp1, _temp2); \
+    min(_temp3, _u);                                  \
+  })
+
+/**
+ * @brief Maximum of five numbers
+ *
+ * This macro evaluates its arguments exactly once.
+ */
+#define max5(x, y, z, w, u)                           \
+  ({                                                  \
+    const __typeof__(x) _x = (x);                     \
+    const __typeof__(y) _y = (y);                     \
+    const __typeof__(z) _z = (z);                     \
+    const __typeof__(w) _w = (w);                     \
+    const __typeof__(u) _u = (u);                     \
+    const __typeof__(x) _temp1 = max(_x, _y);         \
+    const __typeof__(x) _temp2 = max(_z, _w);         \
+    const __typeof__(x) _temp3 = max(_temp1, _temp2); \
+    max(_temp3, _u);                                  \
+  })
+
 #endif /* SWIFT_MINMAX_H */

src/timestep.h

@@ -292,4 +292,51 @@ __attribute__((always_inline)) INLINE static integertime_t get_bpart_timestep(
   return new_dti;
 }
 
+/**
+ * @brief Compute the new (integer) time-step of a given #sink.
+ *
+ * @param sink The #sink.
+ * @param e The #engine (used to get some constants).
+ */
+__attribute__((always_inline)) INLINE static integertime_t get_sink_timestep(
+    const struct sink *restrict sink, const struct engine *restrict e) {
+
+  /* Sink time-step */
+  float new_dt_sink = sink_compute_timestep(sink);
+
+  /* Gravity time-step */
+  float new_dt_self = FLT_MAX, new_dt_ext = FLT_MAX;
+
+  if (e->policy & engine_policy_external_gravity)
+    new_dt_ext = external_gravity_timestep(e->time, e->external_potential,
+                                           e->physical_constants, sink->gpart);
+
+  const float a_hydro[3] = {0.f, 0.f, 0.f};
+  if (e->policy & engine_policy_self_gravity)
+    new_dt_self = gravity_compute_timestep_self(
+        sink->gpart, a_hydro, e->gravity_properties, e->cosmology);
+
+  /* Take the minimum of all */
+  float new_dt = min3(new_dt_sink, new_dt_self, new_dt_ext);
+
+  /* Apply the maximal dibslacement constraint (FLT_MAX  if non-cosmological)*/
+  new_dt = min(new_dt, e->dt_max_RMS_displacement);
+
+  /* Apply cosmology correction (This is 1 if non-cosmological) */
+  new_dt *= e->cosmology->time_step_factor;
+
+  /* Limit timestep within the allowed range */
+  new_dt = min(new_dt, e->dt_max);
+  if (new_dt < e->dt_min) {
+    error("sink (id=%lld) wants a time-step (%e) below dt_min (%e)", sink->id,
+          new_dt, e->dt_min);
+  }
+
+  /* Convert to integer time */
+  const integertime_t new_dti = make_integer_timestep(
+      new_dt, sink->time_bin, num_time_bins, e->ti_current, e->time_base_inv);
+
+  return new_dti;
+}
+
 #endif /* SWIFT_TIMESTEP_H */