From c62fe635952097cae7ca15a0e89c4aed031b4574 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Wed, 14 Sep 2016 14:27:35 +0100
Subject: [PATCH] Added a 'sound speed' variable to the Minimal SPH particle to
save redundant calculations.
---
src/const.h | 4 ++--
src/hydro/Minimal/hydro.h | 34 +++++++++++++++++++++++++++++++++-
src/hydro/Minimal/hydro_iact.h | 8 ++++----
src/hydro/Minimal/hydro_part.h | 3 +++
4 files changed, 42 insertions(+), 7 deletions(-)
diff --git a/src/const.h b/src/const.h
index 7a8bcbf6a8..a02e0b4437 100644
--- a/src/const.h
+++ b/src/const.h
@@ -63,8 +63,8 @@
//#define WENDLAND_C6_KERNEL
/* SPH variant to use */
-//#define MINIMAL_SPH
-#define GADGET2_SPH
+#define MINIMAL_SPH
+//#define GADGET2_SPH
//#define HOPKINS_PE_SPH
//#define DEFAULT_SPH
//#define GIZMO_SPH
diff --git a/src/hydro/Minimal/hydro.h b/src/hydro/Minimal/hydro.h
index 16939d49ac..d4dfc68fd3 100644
--- a/src/hydro/Minimal/hydro.h
+++ b/src/hydro/Minimal/hydro.h
@@ -141,6 +141,12 @@ __attribute__((always_inline)) INLINE static void hydro_set_internal_energy(
/* Compute the pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+
+ /* Compute the sound speed from the pressure*/
+ const float rho_inv = 1.f / p->rho;
+ const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
+
+ p->force.soundspeed = soundspeed;
p->force.pressure = pressure;
}
@@ -160,6 +166,12 @@ __attribute__((always_inline)) INLINE static void hydro_set_entropy(
/* Compute the pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+
+ /* Compute the sound speed from the pressure*/
+ const float rho_inv = 1.f / p->rho;
+ const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
+
+ p->force.soundspeed = soundspeed;
p->force.pressure = pressure;
}
@@ -277,10 +289,15 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
struct part *restrict p, struct xpart *restrict xp, int ti_current,
double timeBase) {
+ /* Compute the pressure */
const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
const float pressure = hydro_get_pressure(p, half_dt);
+
const float rho_inv = 1.f / p->rho;
+ /* Compute the sound speed */
+ const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
+
/* Compute the "grad h" term */
const float grad_h_term =
1.f / (1.f + hydro_dimension_inv * p->h * p->density.rho_dh * rho_inv);
@@ -288,6 +305,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Update variables. */
p->force.f = grad_h_term;
p->force.pressure = pressure;
+ p->force.soundspeed = soundspeed;
}
/**
@@ -347,7 +365,15 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
/* Drift the pressure */
const float dt_entr = (t1 - (p->ti_begin + p->ti_end) / 2) * timeBase;
- p->force.pressure = hydro_get_pressure(p, dt_entr);
+ const float pressure = hydro_get_pressure(p, dt_entr);
+
+ const float rho_inv = 1.f / p->rho;
+
+ /* Compute the new sound speed */
+ const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
+
+ p->force.pressure = pressure;
+ p->force.soundspeed = soundspeed;
}
/**
@@ -392,6 +418,12 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
/* Compute the pressure */
const float pressure = gas_pressure_from_internal_energy(p->rho, p->u);
+
+ /* Compute the sound speed from the pressure*/
+ const float rho_inv = 1.f / p->rho;
+ const float soundspeed = sqrtf(hydro_gamma * pressure * rho_inv);
+
+ p->force.soundspeed = soundspeed;
p->force.pressure = pressure;
}
diff --git a/src/hydro/Minimal/hydro_iact.h b/src/hydro/Minimal/hydro_iact.h
index 3e0770df6a..169947b99e 100644
--- a/src/hydro/Minimal/hydro_iact.h
+++ b/src/hydro/Minimal/hydro_iact.h
@@ -165,8 +165,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_force(
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Compute sound speeds and signal velocity */
- const float ci = sqrtf(hydro_gamma * pressurei / rhoi);
- const float cj = sqrtf(hydro_gamma * pressurej / rhoj);
+ const float ci = pi->force.soundspeed;
+ const float cj = pj->force.soundspeed;
const float v_sig = ci + cj - 3.f * mu_ij;
/* Construct the full viscosity term */
@@ -276,8 +276,8 @@ __attribute__((always_inline)) INLINE static void runner_iact_nonsym_force(
const float mu_ij = fac_mu * r_inv * omega_ij; /* This is 0 or negative */
/* Compute sound speeds and signal velocity */
- const float ci = sqrtf(hydro_gamma * pressurei / rhoi);
- const float cj = sqrtf(hydro_gamma * pressurej / rhoj);
+ const float ci = pi->force.soundspeed;
+ const float cj = pj->force.soundspeed;
const float v_sig = ci + cj - 3.f * mu_ij;
/* Construct the full viscosity term */
diff --git a/src/hydro/Minimal/hydro_part.h b/src/hydro/Minimal/hydro_part.h
index cc8086ae64..8542177278 100644
--- a/src/hydro/Minimal/hydro_part.h
+++ b/src/hydro/Minimal/hydro_part.h
@@ -131,6 +131,9 @@ struct part {
/*! Particle pressure. */
float pressure;
+ /*! Particle soundspeed. */
+ float soundspeed;
+
/*! Particle signal velocity */
float v_sig;
--
GitLab