Commit 554ea68f authored by Loic Hausammann's avatar Loic Hausammann

Add drift/kick/timestep function for sink

parent 19bfdbab
......@@ -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 */
......@@ -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 */
......@@ -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 */
......@@ -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 */
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment