From 844e06eab773241d8df069861b1fc11f687ecd46 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Tue, 9 Feb 2016 21:21:25 +0000
Subject: [PATCH] Smoothing length time derivative computed correctly in the
Gadget-2 case
---
src/hydro/Default/hydro.h | 29 ++++++++++++++++++++++++++---
src/hydro/Gadget2/hydro.h | 2 +-
src/hydro/Gadget2/hydro_iact.h | 9 ++++++++-
src/runner.c | 5 +++++
4 files changed, 40 insertions(+), 5 deletions(-)
diff --git a/src/hydro/Default/hydro.h b/src/hydro/Default/hydro.h
index fd18233b12..3215390e83 100644
--- a/src/hydro/Default/hydro.h
+++ b/src/hydro/Default/hydro.h
@@ -70,9 +70,10 @@ __attribute__((always_inline))
* and add the self-contribution term.
*
* @param p The particle to act upon
+ * @param time The current time
*/
__attribute__((always_inline))
- INLINE static void hydro_end_density(struct part* p) {
+ INLINE static void hydro_end_density(struct part* p, float time) {
/* Some smoothing length multiples. */
const float h = p->h;
@@ -166,12 +167,15 @@ __attribute__((always_inline))
*
* @param p The particle
* @param xp The extended data of the particle
- * @param dt The time-step over which to drift
+ * @param t0 The time at the start of the drift
+ * @param t1 The time at the end of the drift
*/
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
- struct part* p, struct xpart* xp, float dt) {
+ struct part* p, struct xpart* xp, float t0, float t1) {
float u, w;
+ const float dt = t1 - t0;
+
/* Predict internal energy */
w = p->force.u_dt / p->u * dt;
if (fabsf(w) < 0.01f) /* 1st order expansion of exp(w) */
@@ -196,3 +200,22 @@ __attribute__((always_inline))
p->force.h_dt *= p->h * 0.333333333f;
}
+
+/**
+ * @brief Kick the additional variables
+ *
+ * @param p The particle to act upon
+ * @param dt The time-step for this kick
+ */
+__attribute__((always_inline))
+ INLINE static void hydro_kick_extra(struct part* p, float dt) {}
+
+/**
+ * @brief Converts hydro quantity of a particle
+ *
+ * Requires the density to be known
+ *
+ * @param p The particle to act upon
+ */
+__attribute__((always_inline))
+ INLINE static void hydro_convert_quantities(struct part* p) {}
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 53366bd84d..30d714e253 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -203,7 +203,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
*/
__attribute__((always_inline))
INLINE static void hydro_end_force(struct part* p) {
-
+
p->entropy_dt *=
(const_hydro_gamma - 1.f) * powf(p->rho, -(const_hydro_gamma - 1.f));
}
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 31528b63f5..260f78ba0e 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -197,7 +197,7 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float r_inv = 1.0f / r;
/* Get some values in local variables. */
- // const float mi = pi->mass;
+ const float mi = pi->mass;
const float mj = pj->mass;
const float rhoi = pi->rho;
const float rhoj = pj->rho;
@@ -280,6 +280,10 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
pj->a[1] += acc * dx[1];
pj->a[2] += acc * dx[2];
+ /* Get the time derivative for h. */
+ pi->force.h_dt -= mj * dvdr / rhoj * wi_dr;
+ pj->force.h_dt -= mi * dvdr / rhoi * wj_dr;
+
/* Update the signal velocity. */
pi->v_sig = fmaxf(pi->v_sig, v_sig);
pj->v_sig = fmaxf(pj->v_sig, v_sig);
@@ -366,6 +370,9 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
pi->a[1] -= acc * dx[1];
pi->a[2] -= acc * dx[2];
+ /* Get the time derivative for h. */
+ pi->force.h_dt -= mj * dvdr / rhoj * wi_dr;
+
/* Update the signal velocity. */
pi->v_sig = fmaxf(pi->v_sig, v_sig);
diff --git a/src/runner.c b/src/runner.c
index 66c1e02deb..eeee039344 100644
--- a/src/runner.c
+++ b/src/runner.c
@@ -851,8 +851,13 @@ void runner_dokick(struct runner *r, struct cell *c, int timer) {
if (is_fixdt || p->t_end <= t_current) {
/* First, finish the force loop */
+ p->force.h_dt *= p->h * 0.333333333f;
+
+ /* And do the same of the extra variable */
hydro_end_force(p);
+ /* Now we are ready to compute the next time-step size */
+
if (is_fixdt) {
/* Now we have a time step, proceed with the kick */
--
GitLab