Commit e91f9b38 authored by Matthieu Schaller's avatar Matthieu Schaller
Browse files

Added cosmology terms in the drift for parts.

parent 673c75f2
......@@ -2343,7 +2343,6 @@ int cell_has_tasks(struct cell *c) {
void cell_drift_part(struct cell *c, const struct engine *e, int force) {
const float hydro_h_max = e->hydro_properties->h_max;
const double timeBase = e->timeBase;
const integertime_t ti_old_part = c->ti_old_part;
const integertime_t ti_current = e->ti_current;
struct part *const parts = c->parts;
......@@ -2384,7 +2383,19 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
} else if (!c->split && force && ti_current > ti_old_part) {
/* Drift from the last time the cell was drifted to the current time */
const double dt = (ti_current - ti_old_part) * timeBase;
double dt_drift, dt_kick_grav, dt_kick_hydro;
if (e->policy & engine_policy_cosmology) {
dt_drift =
cosmology_get_drift_factor(e->cosmology, ti_old_part, ti_current);
dt_kick_grav =
cosmology_get_grav_kick_factor(e->cosmology, ti_old_part, ti_current);
dt_kick_hydro = cosmology_get_hydro_kick_factor(e->cosmology, ti_old_part,
ti_current);
} else {
dt_drift = (ti_current - ti_old_part) * e->time_base;
dt_kick_grav = (ti_current - ti_old_part) * e->time_base;
dt_kick_hydro = (ti_current - ti_old_part) * e->time_base;
}
/* Loop over all the gas particles in the cell */
const size_t nr_parts = c->count;
......@@ -2395,7 +2406,8 @@ void cell_drift_part(struct cell *c, const struct engine *e, int force) {
struct xpart *const xp = &xparts[k];
/* Drift... */
drift_part(p, xp, dt, timeBase, ti_old_part, ti_current);
drift_part(p, xp, dt_drift, dt_kick_hydro, dt_kick_grav, ti_old_part,
ti_current);
/* Limit h to within the allowed range */
p->h = min(p->h, hydro_h_max);
......
......@@ -68,14 +68,16 @@ __attribute__((always_inline)) INLINE static void drift_gpart(
*
* @param p The #part to drift.
* @param xp The #xpart of the particle.
* @param dt The drift time-step
* @param timeBase The minimal allowed time-step size.
* @param ti_old Integer start of time-step
* @param ti_current Integer end of time-step
* @param dt_drift The drift time-step
* @param dt_kick_grav The kick time-step for gravity accelerations.
* @param dt_kick_hydro The kick time-step for hydro accelerations.
* @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_part(
struct part *restrict p, struct xpart *restrict xp, double dt,
double timeBase, integertime_t ti_old, integertime_t ti_current) {
struct part *restrict p, struct xpart *restrict xp, double dt_drift,
double dt_kick_hydro, double dt_kick_grav, integertime_t ti_old,
integertime_t ti_current) {
#ifdef SWIFT_DEBUG_CHECKS
if (p->ti_drift != ti_old)
......@@ -88,21 +90,25 @@ __attribute__((always_inline)) INLINE static void drift_part(
#endif
/* Drift... */
p->x[0] += xp->v_full[0] * dt;
p->x[1] += xp->v_full[1] * dt;
p->x[2] += xp->v_full[2] * dt;
p->x[0] += xp->v_full[0] * dt_drift;
p->x[1] += xp->v_full[1] * dt_drift;
p->x[2] += xp->v_full[2] * dt_drift;
/* Predict velocities (for hydro terms) */
p->v[0] += p->a_hydro[0] * dt;
p->v[1] += p->a_hydro[1] * dt;
p->v[2] += p->a_hydro[2] * dt;
p->v[0] += p->a_hydro[0] * dt_kick_hydro;
p->v[1] += p->a_hydro[1] * dt_kick_hydro;
p->v[2] += p->a_hydro[2] * dt_kick_hydro;
p->v[0] += xp->a_grav[0] * dt_kick_grav;
p->v[1] += xp->a_grav[1] * dt_kick_grav;
p->v[2] += xp->a_grav[2] * dt_kick_grav;
/* Predict the values of the extra fields */
hydro_predict_extra(p, xp, dt);
hydro_predict_extra(p, xp, dt_drift);
/* Compute offsets since last cell construction */
for (int k = 0; k < 3; k++) {
const float dx = xp->v_full[k] * dt;
const float dx = xp->v_full[k] * dt_drift;
xp->x_diff[k] -= dx;
xp->x_diff_sort[k] -= dx;
}
......
......@@ -513,6 +513,9 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
xp->v_full[0] = p->v[0];
xp->v_full[1] = p->v[1];
xp->v_full[2] = p->v[2];
xp->a_grav[0] = 0.f;
xp->a_grav[1] = 0.f;
xp->a_grav[2] = 0.f;
xp->entropy_full = p->entropy;
hydro_reset_acceleration(p);
......
......@@ -46,6 +46,9 @@ struct xpart {
/* Velocity at the last full step. */
float v_full[3];
/* Gravitational acceleration at the last full step. */
float a_grav[3];
/* Entropy at the last full step. */
float entropy_full;
......
Supports Markdown
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