diff --git a/src/drift.h b/src/drift.h
index a4bdf9be74aade4fe0f1349544cf472363c81c99..7e874fe0ceabe5b091cc7c5bb53adbef2c9a3efd 100644
--- a/src/drift.h
+++ b/src/drift.h
@@ -28,6 +28,7 @@
#include "dimension.h"
#include "hydro.h"
#include "part.h"
+#include "stars.h"
/**
* @brief Perform the 'drift' operation on a #gpart.
@@ -138,6 +139,9 @@ __attribute__((always_inline)) INLINE static void drift_spart(
sp->x[1] += sp->v[1] * dt_drift;
sp->x[2] += sp->v[2] * dt_drift;
+ /* Predict the values of the extra fields */
+ stars_predict_extra(sp, dt_drift);
+
/* Compute offsets since last cell construction */
for (int k = 0; k < 3; k++) {
const float dx = sp->v[k] * dt_drift;
diff --git a/src/engine_maketasks.c b/src/engine_maketasks.c
index cff4ba42f88fcc3a55506f4ebd5822e6db6e883c..ff0baeff6991437a48206e7fd3f358f9b05c5264 100644
--- a/src/engine_maketasks.c
+++ b/src/engine_maketasks.c
@@ -269,6 +269,9 @@ void engine_addtasks_send_stars(struct engine *e, struct cell *ci,
t_feed = scheduler_addtask(s, task_type_send, task_subtype_spart,
ci->mpi.tag, 0, ci, cj);
+ /* The send_stars task should unlock the super_cell's kick task. */
+ scheduler_addunlock(s, t_feed, ci->super->end_force);
+
/* Ghost before you send */
scheduler_addunlock(s, ci->super->stars.ghost_out, t_feed);
}
@@ -492,10 +495,19 @@ void engine_addtasks_recv_stars(struct engine *e, struct cell *c,
scheduler_addunlock(s, l->t, t_feed);
}
+ struct task *recv_rho = NULL;
+ if (c->mpi.hydro.recv_rho != NULL)
+ recv_rho = c->mpi.hydro.recv_rho;
+
for (struct link *l = c->stars.feedback; l != NULL; l = l->next) {
scheduler_addunlock(s, t_feed, l->t);
+
+ /* Need gas density before feedback */
+ if (recv_rho != NULL)
+ scheduler_addunlock(s, c->mpi.hydro.recv_rho, l->t);
}
+
/* Recurse? */
if (c->split)
for (int k = 0; k < 8; k++)
@@ -964,6 +976,16 @@ void engine_make_hierarchical_tasks_stars(struct engine *e, struct cell *c) {
scheduler_addtask(s, task_type_stars_ghost_out, task_subtype_none, 0,
/* implicit = */ 1, c, NULL);
engine_add_stars_ghosts(e, c, c->stars.ghost_in, c->stars.ghost_out);
+
+ /* Need to compute the gas density before moving to the feedback */
+ struct task *hydro_ghost = NULL;
+ if (c->hydro.super)
+ hydro_ghost = c->hydro.super->hydro.ghost_out;
+
+ if (hydro_ghost) {
+ scheduler_addunlock(s, hydro_ghost,
+ c->super->stars.ghost_out);
+ }
}
} else { /* We are above the super-cell so need to go deeper */
diff --git a/src/runner.c b/src/runner.c
index bf3c457fff7222c631540ce4dbebec9217b76f47..a7cdd142920fb9cf8d3d4498acce025dbfc99b53 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -2664,7 +2664,7 @@ void runner_do_end_force(struct runner *r, struct cell *c, int timer) {
if (spart_is_active(sp, e)) {
/* Finish the force loop */
- stars_end_force(sp);
+ stars_end_feedback(sp);
}
}
}
diff --git a/src/scheduler.c b/src/scheduler.c
index 9763d6f9a238dab407b95afce83a5964cd995fda..f443ec8724d01fa41542bb89c414a8204da4f38e 100644
--- a/src/scheduler.c
+++ b/src/scheduler.c
@@ -1147,7 +1147,6 @@ static void scheduler_splittask_stars(struct task *t, struct scheduler *s) {
/* Otherwise, split it. */
} else {
-
/* Take a step back (we're going to recycle the current task)... */
redo = 1;
@@ -1710,7 +1709,6 @@ void scheduler_splittasks_mapper(void *map_data, int num_elements,
* @param s The #scheduler.
*/
void scheduler_splittasks(struct scheduler *s) {
-
/* Call the mapper on each current task. */
threadpool_map(s->threadpool, scheduler_splittasks_mapper, s->tasks,
s->nr_tasks, sizeof(struct task), 0, s);
diff --git a/src/stars/Default/stars.h b/src/stars/Default/stars.h
index 7a60ed6339ea8428376d63a7897ef708bb7f6dec..deb0f671496de84d9b450667a0e862a73fce3f93 100644
--- a/src/stars/Default/stars.h
+++ b/src/stars/Default/stars.h
@@ -65,6 +65,26 @@ __attribute__((always_inline)) INLINE static void stars_init_spart(
sp->density.wcount_dh = 0.f;
}
+/**
+ * @brief Predict additional particle fields forward in time when drifting
+ *
+ * @param p The particle
+ * @param dt_drift The drift time-step for positions.
+ */
+__attribute__((always_inline)) INLINE static void stars_predict_extra(
+ struct spart *restrict sp, float dt_drift) {
+
+ const float h_inv = 1.f / sp->h;
+
+ /* Predict smoothing length */
+ const float w1 = sp->feedback.h_dt * h_inv * dt_drift;
+ if (fabsf(w1) < 0.2f)
+ sp->h *= approx_expf(w1); /* 4th order expansion of exp(w) */
+ else
+ sp->h *= expf(w1);
+
+}
+
/**
* @brief Sets the values to be predicted in the drifts to their values at a
* kick time
@@ -75,12 +95,14 @@ __attribute__((always_inline)) INLINE static void stars_reset_predicted_values(
struct spart* restrict sp) {}
/**
- * @brief Finishes the calculation of the feedback
+ * @brief Finishes the calculation of (non-gravity) forces acting on stars
*
* @param sp The particle to act upon
*/
__attribute__((always_inline)) INLINE static void stars_end_feedback(
- struct spart* sp) {}
+ struct spart* sp) {
+ sp->feedback.h_dt *= sp->h * hydro_dimension_inv;
+}
/**
* @brief Kick the additional variables
@@ -155,6 +177,10 @@ __attribute__((always_inline)) INLINE static void stars_evolve_spart(
*/
__attribute__((always_inline)) INLINE static void stars_reset_acceleration(
struct spart* restrict p) {
+
+ /* Reset time derivative */
+ p->feedback.h_dt = 0.f;
+
#ifdef DEBUG_INTERACTIONS_STARS
for (int i = 0; i < MAX_NUM_OF_NEIGHBOURS_STARS; ++i)
p->ids_ngbs_force[i] = -1;
diff --git a/src/stars/Default/stars_iact.h b/src/stars/Default/stars_iact.h
index 09c6f16c7972506a3a762e0a92bf3c85c4c66345..a99b4ade72be208f991f800594dacb2a276b7100 100644
--- a/src/stars/Default/stars_iact.h
+++ b/src/stars/Default/stars_iact.h
@@ -58,6 +58,27 @@ runner_iact_nonsym_stars_feedback(float r2, const float *dx, float hi, float hj,
struct spart *restrict si,
struct part *restrict pj, float a, float H) {
+ const float mj = pj->mass;
+ const float rhoj = pj->rho;
+ const float r = sqrtf(r2);
+ const float ri = 1.f / r;
+
+ /* Get the kernel for hi. */
+ float hi_inv = 1.0f / hi;
+ float hid_inv = pow_dimension_plus_one(hi_inv); /* 1/h^(d+1) */
+ float xi = r * hi_inv;
+ float wi, wi_dx;
+ kernel_deval(xi, &wi, &wi_dx);
+ float wi_dr = hid_inv * wi_dx;
+
+ /* Compute dv dot r */
+ float dvdr = (si->v[0] - pj->v[0]) * dx[0] +
+ (si->v[1] - pj->v[1]) * dx[1] +
+ (si->v[2] - pj->v[2]) * dx[2];
+
+ /* Get the time derivative for h. */
+ si->feedback.h_dt -= mj * dvdr * ri / rhoj * wi_dr;
+
#ifdef DEBUG_INTERACTIONS_STARS
/* Update ngb counters */
if (si->num_ngb_force < MAX_NUM_OF_NEIGHBOURS_STARS)
diff --git a/src/stars/Default/stars_part.h b/src/stars/Default/stars_part.h
index 25342ae9a93e94616536bad1ff9b5f6dff81cc17..49a95a2d406e7d0b741daaac1bd3d8235d3a4d87 100644
--- a/src/stars/Default/stars_part.h
+++ b/src/stars/Default/stars_part.h
@@ -65,6 +65,12 @@ struct spart {
} density;
+ struct {
+ /* Change in smoothing length over time. */
+ float h_dt;
+
+ } feedback;
+
#ifdef SWIFT_DEBUG_CHECKS
/* Time of the last drift */
diff --git a/src/timestep.h b/src/timestep.h
index e9943a41a0536b65944f0256c827d43386aadd88..94301046682cd068cd9ab3ce8434be0d08238b93 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -201,8 +201,15 @@ __attribute__((always_inline)) INLINE static integertime_t get_spart_timestep(
new_dt_self = gravity_compute_timestep_self(
sp->gpart, a_hydro, e->gravity_properties, e->cosmology);
+ /* Limit change in smoothing length */
+ const float dt_h_change =
+ (sp->feedback.h_dt != 0.0f)
+ ? fabsf(e->stars_properties->log_max_h_change * sp->h / sp->feedback.h_dt)
+ : FLT_MAX;
+
/* Take the minimum of all */
float new_dt = min3(new_dt_stars, new_dt_self, new_dt_ext);
+ new_dt = min(new_dt, dt_h_change);
/* Apply the maximal displacement constraint (FLT_MAX if non-cosmological)*/
new_dt = min(new_dt, e->dt_max_RMS_displacement);
@@ -212,9 +219,10 @@ __attribute__((always_inline)) INLINE static integertime_t get_spart_timestep(
/* Limit timestep within the allowed range */
new_dt = min(new_dt, e->dt_max);
- if (new_dt < e->dt_min)
+ if (new_dt < e->dt_min) {
error("spart (id=%lld) wants a time-step (%e) below dt_min (%e)", sp->id,
new_dt, e->dt_min);
+ }
/* Convert to integer time */
const integertime_t new_dti = make_integer_timestep(