Merge branch 'moving_mesh' into derijcke_cooling

src/hydro/Shadowswift/hydro.h

@@ -339,16 +339,10 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
 
 #ifdef SHADOWSWIFT_EXTRAPOLATE_TIME
   /* Extrapolate primitive quantities in time */
-  float W[5], dW[5];
+  float W[6], dW[6];
   hydro_part_get_primitive_variables(p, W);
-  hydro_gradients_extrapolate_in_time(p, W, 0.5f * dt_therm, /*return*/ dW);
-
-  /* Update primitive quantities with extrapolations */
-  p->dW_time[0] += dW[0];
-  p->dW_time[1] += dW[1];
-  p->dW_time[2] += dW[2];
-  p->dW_time[3] += dW[3];
-  p->dW_time[4] += dW[4];
+  hydro_gradients_extrapolate_in_time(p, W, 0.5f * dt_therm, p->dW_time);
+
   // MATTHIEU: Apply the entropy floor here.
 #endif
 
@@ -412,7 +406,7 @@ hydro_convert_conserved_to_primitive(struct part *p, struct xpart *xp,
                 hydro_get_comoving_psize(p);
   float thermal_energy = Q[4] - Ekin;
   if (thermal_energy > 1e-2 * (Ekin + Egrav)) {
-    /* Recover pressure, thermal energy and entropy from total energy */
+    /* Recover thermal energy and entropy from total energy */
     p->thermal_energy = thermal_energy;
     W[5] = gas_entropy_from_internal_energy(W[0], thermal_energy * m_inv);
     p->conserved.entropy = Q[0] * W[5];
@@ -425,9 +419,10 @@ hydro_convert_conserved_to_primitive(struct part *p, struct xpart *xp,
   } else {
     /* Use evolved thermal energy to set entropy and total energy */
     p->conserved.energy = Ekin + p->thermal_energy;
-    W[5] = gas_entropy_from_internal_energy(W[0], thermal_energy * m_inv);
+    W[5] = gas_entropy_from_internal_energy(W[0], p->thermal_energy * m_inv);
     p->conserved.entropy = Q[0] * W[5];
   }
+  /* Calculate pressure from thermal energy */
   W[4] = gas_pressure_from_internal_energy(W[0], p->thermal_energy * m_inv);
 #endif
   /* reset the primitive variables if we are using Lloyd's algorithm */
@@ -654,6 +649,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
     p->dW_time[2] = 0.0f;
     p->dW_time[3] = 0.0f;
     p->dW_time[4] = 0.0f;
+    p->dW_time[5] = 0.0f;
 #endif
 
     /* Apply the minimal energy limit */

src/hydro/Shadowswift/hydro_gradients.h

@@ -73,7 +73,8 @@ hydro_gradients_extrapolate_in_time(const struct part* p, const float* W,
   }
   dW[4] = -dt * (hydro_gamma * W[4] * div_v + W[1] * dP[0] + W[2] * dP[1] +
                  W[3] * dP[2]);
-  dW[5] = -dt * (W[1] * dA[0] + W[2] * dA[1] + W[3] * dA[1]);
+  /* See eq. 51 in springel 2010 */
+  dW[5] = -dt * (W[1] * dA[0] + W[2] * dA[1] + W[3] * dA[2]);
 }
 
 /**

src/hydro/Shadowswift/hydro_gradients_shadowswift.h

@@ -125,12 +125,14 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
   dWi[2] += pi->dW_time[2];
   dWi[3] += pi->dW_time[3];
   dWi[4] += pi->dW_time[4];
+  dWi[5] += pi->dW_time[5];
 
   dWj[0] += pj->dW_time[0];
   dWj[1] += pj->dW_time[1];
   dWj[2] += pj->dW_time[2];
   dWj[3] += pj->dW_time[3];
   dWj[4] += pj->dW_time[4];
+  dWj[5] += pj->dW_time[5];
 #endif
 
   /* Apply the slope limiter at this interface */

src/hydro/Shadowswift/hydro_gradients_wls.h

@@ -154,12 +154,14 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
   dWi[2] += pi->dW_time[2];
   dWi[3] += pi->dW_time[3];
   dWi[4] += pi->dW_time[4];
+  dWi[5] += pi->dW_time[5];
 
   dWj[0] += pj->dW_time[0];
   dWj[1] += pj->dW_time[1];
   dWj[2] += pj->dW_time[2];
   dWj[3] += pj->dW_time[3];
   dWj[4] += pj->dW_time[4];
+  dWj[5] += pj->dW_time[5];
 #endif
 
   /* Apply the slope limiter at this interface */

src/hydro/Shadowswift/hydro_part.h

@@ -191,7 +191,7 @@ struct part {
 
 #ifdef SHADOWSWIFT_EXTRAPOLATE_TIME
   /* Time extrapolations of primitive variables (cumulative over timestep) */
-  float dW_time[5];
+  float dW_time[6];
 #endif
 
   /* The conserved hydrodynamical variables. */

src/hydro/Shadowswift/hydro_setters.h

@@ -80,7 +80,10 @@ __attribute__((always_inline)) INLINE static void
 hydro_set_comoving_internal_energy_dt(struct part* restrict p,
                                       const float du_dt) {
   const float old_du_dt = hydro_get_comoving_internal_energy_dt(p);
-  p->flux.energy += p->conserved.mass * (du_dt - old_du_dt) * p->flux.dt;
+  const float du = (du_dt - old_du_dt) * p->flux.dt;
+  p->flux.energy += p->conserved.mass * du;
+  p->flux.entropy +=
+      p->conserved.mass * gas_entropy_from_internal_energy(p->rho, du);
 }
 
 /**
@@ -119,19 +122,20 @@ hydro_set_comoving_internal_energy(struct part* p, const float u) {
     return;
   }
 
-  /* thermal_energy is NOT the specific energy (u), but the total thermal
-     energy (u*m) */
-  p->thermal_energy = u * p->conserved.mass;
-
   const float Ekin = 0.5f *
                      (p->conserved.momentum[0] * p->conserved.momentum[0] +
                       p->conserved.momentum[1] * p->conserved.momentum[1] +
                       p->conserved.momentum[2] * p->conserved.momentum[2]) /
                      mass;
 
-  p->conserved.energy = p->thermal_energy + Ekin;
   p->P = gas_pressure_from_internal_energy(p->rho, u);
   p->A = gas_entropy_from_internal_energy(p->rho, u);
+
+  /* thermal_energy is NOT the specific energy (u), but the total thermal
+     energy (u*m) */
+  p->thermal_energy = p->conserved.mass * u;
+  p->conserved.energy = p->thermal_energy + Ekin;
+  p->conserved.entropy = p->conserved.mass * p->A;
 }
 
 /**
@@ -200,31 +204,6 @@ hydro_diffusive_feedback_reset(struct part* restrict p) {
   /* Purposefully left empty */
 }
 
-/**
- * @brief Modifies the thermal state of a particle to the imposed entropy
- *
- * This overrides the current state of the particle but does *not* change its
- * time-derivatives
- *
- * @param p The particle
- * @param A The new entropy
- */
-__attribute__((always_inline)) INLINE static void hydro_set_entropy(
-    struct part* restrict p, float A) {
-
-  p->thermal_energy = A * pow_gamma_minus_one(p->rho) *
-                      hydro_one_over_gamma_minus_one * p->conserved.mass;
-  /* add the kinetic energy */
-  p->conserved.energy =
-      p->thermal_energy + 0.5f * p->conserved.mass *
-                              (p->conserved.momentum[0] * p->v[0] +
-                               p->conserved.momentum[1] * p->v[1] +
-                               p->conserved.momentum[2] * p->v[2]);
-
-  p->P = gas_pressure_from_entropy(p->rho, A);
-  p->A = A;
-}
-
 /**
  * @brief Overwrite the initial internal energy of a particle.
  *

src/hydro/Shadowswift/hydro_slope_limiters.h

@@ -86,6 +86,8 @@ __attribute__((always_inline)) INLINE static void hydro_slope_limiter_prepare(
   p->limiter.v[2][1] = p->v[2];
   p->limiter.P[0] = p->P;
   p->limiter.P[1] = p->P;
+  p->limiter.A[0] = p->A;
+  p->limiter.A[1] = p->A;
 
   p->limiter.extrapolations.rho[0] = 0.f;
   p->limiter.extrapolations.rho[1] = 0.f;
@@ -97,6 +99,8 @@ __attribute__((always_inline)) INLINE static void hydro_slope_limiter_prepare(
   p->limiter.extrapolations.v[2][1] = 0.f;
   p->limiter.extrapolations.P[0] = 0.f;
   p->limiter.extrapolations.P[1] = 0.f;
+  p->limiter.extrapolations.A[0] = 0.f;
+  p->limiter.extrapolations.A[1] = 0.f;
 }
 
 #endif  // SWIFTSIM_SHADOWSWIFT_HYDRO_SLOPE_LIMITERS_H