From aaf537a92b080fbd17f610f4ca2db62d3a282b1c Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Sat, 13 Feb 2016 17:07:29 +0000
Subject: [PATCH] Soundspeed gets pre-computed before the force loop
---
src/hydro/Gadget2/hydro.h | 10 +++++++---
src/hydro/Gadget2/hydro_debug.h | 5 +++--
src/hydro/Gadget2/hydro_iact.h | 8 ++++----
src/hydro/Gadget2/hydro_part.h | 3 +++
4 files changed, 17 insertions(+), 9 deletions(-)
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 0c359cce21..88984d49ea 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -126,6 +126,9 @@ INLINE static void hydro_prepare_force(struct part* p, struct xpart* xp, float t
const float dt = time - 0.5f * (p->t_begin + p->t_end);
p->force.pressure =
(p->entropy + p->entropy_dt * dt) * powf(p->rho, const_hydro_gamma);
+
+ /* Compute the sound speed */
+ p->force.soundspeed = sqrtf(const_hydro_gamma * p->force.pressure / p->rho);
}
/**
@@ -164,12 +167,13 @@ __attribute__((always_inline))
__attribute__((always_inline)) INLINE static void hydro_predict_extra(
struct part* p, struct xpart* xp, float t0, float t1) {
- // p->rho *= expf(-p->div_v * dt);
- // p->h *= expf(0.33333333f * p->div_v * dt)
-
+ /* Drift the pressure */
const float dt_entr = t1 - 0.5f * (p->t_begin + p->t_end);
p->force.pressure =
(p->entropy + p->entropy_dt * dt_entr) * powf(p->rho, const_hydro_gamma);
+
+ /* Compute the new sound speed */
+ p->force.soundspeed = sqrtf(const_hydro_gamma * p->force.pressure / p->rho);
}
/**
diff --git a/src/hydro/Gadget2/hydro_debug.h b/src/hydro/Gadget2/hydro_debug.h
index 9948d8f0f4..0f3f3ac945 100644
--- a/src/hydro/Gadget2/hydro_debug.h
+++ b/src/hydro/Gadget2/hydro_debug.h
@@ -24,13 +24,14 @@ __attribute__((always_inline))
"v=[%.3e,%.3e,%.3e],v_full=[%.3e,%.3e,%.3e] \n a=[%.3e,%.3e,%.3e],\n "
"h=%.3e, "
"wcount=%d, m=%.3e, dh_drho=%.3e, rho=%.3e, P=%.3e, S=%.3e, "
- "dS/dt=%.3e,\n"
+ "dS/dt=%.3e, c=%.3e\n"
"divV=%.3e, curlV=%.3e, rotV=[%.3e,%.3e,%.3e] \n "
"v_sig=%e dh/dt=%.3e t_begin=%.3e, t_end=%.3e\n",
p->x[0], p->x[1], p->x[2], p->v[0], p->v[1], p->v[2], xp->v_full[0],
xp->v_full[1], xp->v_full[2], p->a_hydro[0], p->a_hydro[1], p->a_hydro[2],
p->h, (int)p->density.wcount, p->mass, p->rho_dh, p->rho, p->force.pressure,
- p->entropy, p->entropy_dt, p->div_v, p->force.curl_v, p->density.rot_v[0],
+ p->entropy, p->entropy_dt, p->force.soundspeed,
+ p->div_v, p->force.curl_v, p->density.rot_v[0],
p->density.rot_v[1], p->density.rot_v[2], p->force.v_sig, p->force.h_dt,
p->t_begin, p->t_end);
}
diff --git a/src/hydro/Gadget2/hydro_iact.h b/src/hydro/Gadget2/hydro_iact.h
index 80111a759e..0919b14427 100644
--- a/src/hydro/Gadget2/hydro_iact.h
+++ b/src/hydro/Gadget2/hydro_iact.h
@@ -223,8 +223,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
/* Compute sound speeds */
- const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
- const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
+ const float ci = pi->force.soundspeed;
+ const float cj = pj->force.soundspeed;
float v_sig = ci + cj;
/* Compute dv dot r. */
@@ -334,8 +334,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float P_over_rho_j = pressurej / (rhoj * rhoj) * pj->rho_dh;
/* Compute sound speeds */
- const float ci = sqrtf(const_hydro_gamma * pressurei / rhoi);
- const float cj = sqrtf(const_hydro_gamma * pressurej / rhoj);
+ const float ci = pi->force.soundspeed;
+ const float cj = pj->force.soundspeed;
float v_sig = ci + cj;
/* Compute dv dot r. */
diff --git a/src/hydro/Gadget2/hydro_part.h b/src/hydro/Gadget2/hydro_part.h
index 3fbf5dfb82..cda0123b5e 100644
--- a/src/hydro/Gadget2/hydro_part.h
+++ b/src/hydro/Gadget2/hydro_part.h
@@ -95,6 +95,9 @@ struct part {
/* Particle pressure */
float pressure;
+ /* Particle sound speed */
+ float soundspeed;
+
} force;
};
--
GitLab