diff --git a/src/feedback/EAGLE/feedback_iact.h b/src/feedback/EAGLE/feedback_iact.h
index 2576e0f05947d0c33230b3140885b0a4fbcaa58e..1cf8e42995dd87bc6c9dbf109fc877f73ea024cc 100644
--- a/src/feedback/EAGLE/feedback_iact.h
+++ b/src/feedback/EAGLE/feedback_iact.h
@@ -21,6 +21,7 @@
/* Local includes */
#include "random.h"
+#include "timestep_sync.h"
/**
* @brief Density interaction between two particles (non-symmetric).
@@ -36,13 +37,11 @@
* @param ti_current Current integer time value
*/
__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_feedback_density(const float r2, const float *dx,
- const float hi, const float hj,
- struct spart *restrict si,
- const struct part *restrict pj,
- const struct xpart *restrict xpj,
- const struct cosmology *restrict cosmo,
- const integertime_t ti_current) {
+runner_iact_nonsym_feedback_density(
+ const float r2, const float *dx, const float hi, const float hj,
+ struct spart *restrict si, const struct part *restrict pj,
+ const struct xpart *restrict xpj, const struct cosmology *restrict cosmo,
+ const struct engine *e, const integertime_t ti_current) {
/* Get the gas mass. */
const float mj = hydro_get_mass(pj);
@@ -85,13 +84,11 @@ runner_iact_nonsym_feedback_density(const float r2, const float *dx,
* generator
*/
__attribute__((always_inline)) INLINE static void
-runner_iact_nonsym_feedback_apply(const float r2, const float *dx,
- const float hi, const float hj,
- const struct spart *restrict si,
- struct part *restrict pj,
- struct xpart *restrict xpj,
- const struct cosmology *restrict cosmo,
- const integertime_t ti_current) {
+runner_iact_nonsym_feedback_apply(
+ const float r2, const float *dx, const float hi, const float hj,
+ const struct spart *restrict si, struct part *restrict pj,
+ struct xpart *restrict xpj, const struct cosmology *restrict cosmo,
+ const struct engine *e, const integertime_t ti_current) {
/* Get r and 1/r. */
const float r_inv = 1.0f / sqrtf(r2);
@@ -292,10 +289,13 @@ runner_iact_nonsym_feedback_apply(const float r2, const float *dx,
/* Impose maximal viscosity */
hydro_diffusive_feedback_reset(pj);
- /* message( */
- /* "We did some heating! id %llu star id %llu probability %.5e " */
- /* "random_num %.5e du %.5e du/ini %.5e", */
- /* pj->id, si->id, prob, rand, delta_u, delta_u / u_init); */
+ message(
+ "We did some heating! id %llu star id %llu probability %.5e "
+ "random_num %.5e du %.5e du/ini %.5e",
+ pj->id, si->id, prob, rand, delta_u, delta_u / u_init);
+
+ /* Synchronize the particle on the timeline */
+ timestep_sync_part(pj, xpj, e, cosmo);
}
}
}
diff --git a/src/runner_doiact_functions_stars.h b/src/runner_doiact_functions_stars.h
index d452fba01b38de09d71ca4a121dcd5d92388bc08..bb5e7ecf30c420d27e07646efadee5fbacaa7370 100644
--- a/src/runner_doiact_functions_stars.h
+++ b/src/runner_doiact_functions_stars.h
@@ -110,9 +110,9 @@ void DOSELF1_STARS(struct runner *r, struct cell *c, int timer) {
IACT_STARS(r2, dx, hi, hj, si, pj, a, H);
#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_feedback_density(r2, dx, hi, hj, si, pj, NULL, cosmo,
- ti_current);
+ e, ti_current);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
- runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, si, pj, xpj, cosmo,
+ runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, si, pj, xpj, cosmo, e,
ti_current);
#endif
}
@@ -219,9 +219,9 @@ void DO_NONSYM_PAIR1_STARS_NAIVE(struct runner *r, struct cell *restrict ci,
#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_feedback_density(r2, dx, hi, hj, si, pj, NULL, cosmo,
- ti_current);
+ e, ti_current);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
- runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, si, pj, xpj, cosmo,
+ runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, si, pj, xpj, cosmo, e,
ti_current);
#endif
}
@@ -392,10 +392,10 @@ void DO_SYM_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_feedback_density(r2, dx, hi, hj, spi, pj, NULL,
- cosmo, ti_current);
+ cosmo, e, ti_current);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, xpj, cosmo,
- ti_current);
+ e, ti_current);
#endif
}
} /* loop over the parts in cj. */
@@ -521,10 +521,10 @@ void DO_SYM_PAIR1_STARS(struct runner *r, struct cell *ci, struct cell *cj,
#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_feedback_density(r2, dx, hj, hi, spj, pi, xpi,
- cosmo, ti_current);
+ cosmo, e, ti_current);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
runner_iact_nonsym_feedback_apply(r2, dx, hj, hi, spj, pi, xpi, cosmo,
- ti_current);
+ e, ti_current);
#endif
}
} /* loop over the parts in ci. */
@@ -649,10 +649,10 @@ void DOPAIR1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_feedback_density(r2, dx, hi, hj, spi, pj, NULL,
- cosmo, ti_current);
+ cosmo, e, ti_current);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, xpj, cosmo,
- ti_current);
+ e, ti_current);
#endif
}
} /* loop over the parts in cj. */
@@ -711,10 +711,10 @@ void DOPAIR1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_feedback_density(r2, dx, hi, hj, spi, pj, NULL,
- cosmo, ti_current);
+ cosmo, e, ti_current);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, xpj, cosmo,
- ti_current);
+ e, ti_current);
#endif
}
} /* loop over the parts in cj. */
@@ -812,10 +812,10 @@ void DOPAIR1_SUBSET_STARS_NAIVE(struct runner *r, struct cell *restrict ci,
#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_feedback_density(r2, dx, hi, hj, spi, pj, NULL,
- cosmo, ti_current);
+ cosmo, e, ti_current);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
runner_iact_nonsym_feedback_apply(r2, dx, hi, hj, spi, pj, xpj, cosmo,
- ti_current);
+ e, ti_current);
#endif
}
} /* loop over the parts in cj. */
@@ -903,10 +903,10 @@ void DOSELF1_SUBSET_STARS(struct runner *r, struct cell *restrict ci,
IACT_STARS(r2, dx, hi, pj->h, spi, pj, a, H);
#if (FUNCTION_TASK_LOOP == TASK_LOOP_DENSITY)
runner_iact_nonsym_feedback_density(r2, dx, hi, pj->h, spi, pj, NULL,
- cosmo, ti_current);
+ cosmo, e, ti_current);
#elif (FUNCTION_TASK_LOOP == TASK_LOOP_FEEDBACK)
runner_iact_nonsym_feedback_apply(r2, dx, hi, pj->h, spi, pj, xpj,
- cosmo, ti_current);
+ cosmo, e, ti_current);
#endif
}
} /* loop over the parts in cj. */
diff --git a/src/timestep_sync.h b/src/timestep_sync.h
new file mode 100644
index 0000000000000000000000000000000000000000..aee0cde666a8c2f4f26fdebfafec38a16a200018
--- /dev/null
+++ b/src/timestep_sync.h
@@ -0,0 +1,128 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2018 Matthieu Schaller (matthieu.schaller@durham.ac.uk)
+ *
+ * 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_TIMESTEP_SYNC_H
+#define SWIFT_TIMESTEP_SYNC_H
+
+/* Config parameters. */
+#include "../config.h"
+
+#include "kick.h"
+
+INLINE static void timestep_sync_part(struct part *p, struct xpart *xp, const struct engine *e,
+ const struct cosmology *cosmo) {
+
+ const int with_cosmology = (e->policy & engine_policy_cosmology);
+ const integertime_t ti_current = e->ti_current;
+ const timebin_t max_active_bin = e->max_active_bin;
+ const double time_base = e->time_base;
+
+ /* This particle is already active. Nothing to do here... */
+ if (p->time_bin <= max_active_bin) return;
+
+ message("Synchronizing particle %lld", p->id);
+
+ /* We want to make the particle finish it's time-step now. */
+
+ /* Start by recovering the start and end point of the particle's time-step. */
+ const integertime_t old_ti_beg = get_integer_time_begin(ti_current, p->time_bin);
+ const integertime_t old_ti_end = get_integer_time_end(ti_current, p->time_bin);
+
+ /* Old time-step length on the time-line */
+ const integertime_t old_dti = old_ti_end - old_ti_beg;
+
+ /* The actual time-step size this particle will use */
+ const integertime_t new_ti_beg = old_ti_beg;
+ const integertime_t new_ti_end = ti_current;
+ const integertime_t new_dti = new_ti_end - new_ti_beg;
+
+#ifdef SWIFT_DEBUG_CHECKS
+ /* Some basic safety checks */
+ if (old_ti_beg >= ti_current)
+ error(
+ "Incorrect value for old time-step beginning ti_current=%lld, "
+ "old_ti_beg=%lld",
+ ti_current, old_ti_beg);
+
+ if (old_ti_end <= ti_current)
+ error(
+ "Incorrect value for old time-step end ti_current=%lld, "
+ "old_ti_end=%lld",
+ ti_current, old_ti_end);
+
+ if (new_ti_end > old_ti_end) error("New end of time-step after the old one");
+
+ if (new_dti > old_dti) error("New time-step larger than old one");
+#endif
+
+
+ double dt_kick_grav = 0., dt_kick_hydro = 0., dt_kick_therm = 0.,
+ dt_kick_corr = 0.;
+
+ /* Now we need to reverse the kick1... (the dt are negative here) */
+ if (with_cosmology) {
+ dt_kick_hydro = -cosmology_get_hydro_kick_factor(cosmo, old_ti_beg,
+ old_ti_beg + old_dti / 2);
+ dt_kick_grav = -cosmology_get_grav_kick_factor(cosmo, old_ti_beg,
+ old_ti_beg + old_dti / 2);
+ dt_kick_therm = -cosmology_get_therm_kick_factor(cosmo, old_ti_beg,
+ old_ti_beg + old_dti / 2);
+ dt_kick_corr = -cosmology_get_corr_kick_factor(cosmo, old_ti_beg,
+ old_ti_beg + old_dti / 2);
+ } else {
+ dt_kick_hydro = -(old_dti / 2) * time_base;
+ dt_kick_grav = -(old_dti / 2) * time_base;
+ dt_kick_therm = -(old_dti / 2) * time_base;
+ dt_kick_corr = -(old_dti / 2) * time_base;
+ }
+
+ kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm, dt_kick_corr,
+ e->cosmology, e->hydro_properties, e->entropy_floor,
+ old_ti_beg + old_dti / 2, old_ti_beg);
+
+ /* ...and apply the new one (dt is positiive) */
+ if (with_cosmology) {
+ dt_kick_hydro = cosmology_get_hydro_kick_factor(cosmo, new_ti_beg,
+ new_ti_beg + new_dti );
+ dt_kick_grav = cosmology_get_grav_kick_factor(cosmo, new_ti_beg,
+ new_ti_beg + new_dti );
+ dt_kick_therm = cosmology_get_therm_kick_factor(cosmo, new_ti_beg,
+ new_ti_beg + new_dti );
+ dt_kick_corr = cosmology_get_corr_kick_factor(cosmo, new_ti_beg,
+ new_ti_beg + new_dti );
+ } else {
+ dt_kick_hydro = (new_dti ) * time_base;
+ dt_kick_grav = (new_dti ) * time_base;
+ dt_kick_therm = (new_dti ) * time_base;
+ dt_kick_corr = (new_dti ) * time_base;
+ }
+
+ kick_part(p, xp, dt_kick_hydro, dt_kick_grav, dt_kick_therm, dt_kick_corr,
+ e->cosmology, e->hydro_properties, e->entropy_floor, new_ti_beg,
+ new_ti_beg + new_dti);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (p->ti_kick != ti_current) error("Particle has not been synchronized correctly.");
+#endif
+
+
+ /* The particle is now ready to compute its new time-step size and for the next kick */
+}
+
+
+#endif /* SWIFT_TIMESTEP_SYNC_H */