From 816954ee65c25d0cbcd92320018aafa07ea12af4 Mon Sep 17 00:00:00 2001
From: Matthieu Schaller <matthieu.schaller@durham.ac.uk>
Date: Sat, 6 Aug 2016 15:39:17 +0100
Subject: [PATCH] Define thermodynamical relations for the current EoS and use
them for the Gadget scheme
---
src/Makefile.am | 1 +
src/adiabatic_index.h | 4 ++
src/const.h | 4 ++
src/equation_of_state.h | 96 ++++++++++++++++++++++++++++++++
src/hydro/Gadget2/hydro.h | 113 +++++++++++++++++++-------------------
src/timestep.h | 2 +-
6 files changed, 163 insertions(+), 57 deletions(-)
create mode 100644 src/equation_of_state.h
diff --git a/src/Makefile.am b/src/Makefile.am
index 4aa7278ba5..b0d58ebe49 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -56,6 +56,7 @@ AM_SOURCES = space.c runner.c queue.c task.c cell.c engine.c \
nobase_noinst_HEADERS = approx_math.h atomic.h cycle.h error.h inline.h kernel_hydro.h kernel_gravity.h \
kernel_long_gravity.h vector.h runner_doiact.h runner_doiact_grav.h runner_doiact_fft.h \
units.h intrinsics.h minmax.h kick.h timestep.h drift.h adiabatic_index.h io_properties.h \
+ equation_of_state.h
gravity.h gravity_io.h \
gravity/Default/gravity.h gravity/Default/gravity_iact.h gravity/Default/gravity_io.h \
gravity/Default/gravity_debug.h gravity/Default/gravity_part.h \
diff --git a/src/adiabatic_index.h b/src/adiabatic_index.h
index ba7a1a64d5..50d24e8815 100644
--- a/src/adiabatic_index.h
+++ b/src/adiabatic_index.h
@@ -49,6 +49,10 @@
#define hydro_gamma_minus_one 1.f
#define hydro_one_over_gamma_minus_one 1.f
+#else
+
+#error "An adiabatic index needs to be chosen in const.h !"
+
#endif
/**
diff --git a/src/const.h b/src/const.h
index 736ba0155e..98143f0e53 100644
--- a/src/const.h
+++ b/src/const.h
@@ -41,6 +41,10 @@
//#define HYDRO_GAMMA_4_3
//#define HYDRO_GAMMA_2_1
+/* Equation of state choice */
+#define IDEAL_GAS
+//#define ISOTHERMAL_GAS
+
/* Kernel function to use */
#define CUBIC_SPLINE_KERNEL
//#define QUARTIC_SPLINE_KERNEL
diff --git a/src/equation_of_state.h b/src/equation_of_state.h
new file mode 100644
index 0000000000..410b515d11
--- /dev/null
+++ b/src/equation_of_state.h
@@ -0,0 +1,96 @@
+/*******************************************************************************
+ * This file is part of SWIFT.
+ * Copyright (c) 2016 Matthieu Schaller (matthieu.schaller@durham.ac.uk).
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU Lesser General Public License as published
+ * by the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ ******************************************************************************/
+#ifndef SWIFT_EQUATION_OF_STATE_H
+#define SWIFT_EQUATION_OF_STATE_H
+
+/* Config parameters. */
+#include "../config.h"
+
+/* Some standard headers. */
+#include <math.h>
+
+/* Local headers. */
+#include "adiabatic_index.h"
+#include "const.h"
+#include "debug.h"
+#include "inline.h"
+
+#if defined(IDEAL_GAS)
+
+/**
+ * @brief Returns the internal energy given density and entropy
+ *
+ * @param density The density
+ * @param entropy The entropy
+ */
+__attribute__((always_inline)) INLINE static float
+gas_internal_energy_from_entropy(float density, float entropy) {
+
+ return entropy * pow_gamma_minus_one(density) *
+ hydro_one_over_gamma_minus_one;
+}
+
+/**
+ * @brief Returns the pressure given density and entropy
+ *
+ * @param density The density
+ * @param entropy The entropy
+ */
+__attribute__((always_inline)) INLINE static float gas_pressure_from_entropy(
+ float density, float entropy) {
+
+ return entropy * pow_gamma(density);
+}
+
+/**
+ * @brief Returns the entropy given density and internal energy
+ *
+ * @param density The density
+ * @param u The internal energy
+ */
+__attribute__((always_inline)) INLINE static float
+gas_entropy_from_internal_energy(float density, float u) {
+
+ return hydro_gamma_minus_one * u * pow_minus_gamma_minus_one(density);
+}
+
+/**
+ * @brief Returns the pressure given density and internal energy
+ *
+ * @param density The density
+ * @param u The internal energy
+ */
+__attribute__((always_inline)) INLINE static float
+gas_pressure_from_internal_energy(float density, float u) {
+
+ return hydro_gamma_minus_one * u * density;
+}
+
+
+#elif defined(ISOTHERMAL_GAS)
+
+#error "Missing definitions !"
+
+#else
+
+#error "An Equation of state needs to be chosen in const.h !"
+
+#endif
+
+#endif /* SWIFT_EQUATION_OF_STATE_H */
diff --git a/src/hydro/Gadget2/hydro.h b/src/hydro/Gadget2/hydro.h
index 5efce17d78..813d655b77 100644
--- a/src/hydro/Gadget2/hydro.h
+++ b/src/hydro/Gadget2/hydro.h
@@ -18,6 +18,59 @@
******************************************************************************/
#include "adiabatic_index.h"
+#include "equation_of_state.h"
+
+/**
+ * @brief Returns the internal energy of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
+ const struct part *restrict p, float dt) {
+
+ const float entropy = p->entropy + p->entropy_dt * dt;
+
+ return gas_internal_energy_from_entropy(p->rho, entropy);
+}
+
+/**
+ * @brief Returns the pressure of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_pressure(
+ const struct part *restrict p, float dt) {
+
+ const float entropy = p->entropy + p->entropy_dt * dt;
+
+ return gas_pressure_from_entropy(p->rho, entropy);
+}
+
+/**
+ * @brief Returns the entropy of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_entropy(
+ const struct part *restrict p, float dt) {
+
+ return p->entropy + p->entropy_dt * dt;
+}
+
+/**
+ * @brief Returns the sound speed of a particle
+ *
+ * @param p The particle of interest
+ * @param dt Time since the last kick
+ */
+__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
+ const struct part *restrict p, float dt) {
+
+ return p->force.soundspeed;
+}
/**
* @brief Computes the hydro time-step of a given particle
@@ -140,8 +193,7 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
/* Compute the pressure */
const float half_dt = (ti_current - (p->ti_begin + p->ti_end) / 2) * timeBase;
- const float pressure =
- (p->entropy + p->entropy_dt * half_dt) * pow_gamma(p->rho);
+ const float pressure = hydro_get_pressure(p, half_dt);
const float irho = 1.f / p->rho;
@@ -200,8 +252,7 @@ __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;
- const float pressure =
- (p->entropy + p->entropy_dt * dt_entr) * pow_gamma(p->rho);
+ const float pressure = hydro_get_pressure(p, dt_entr);
const float irho = 1.f / p->rho;
@@ -265,56 +316,6 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
__attribute__((always_inline)) INLINE static void hydro_convert_quantities(
struct part *restrict p) {
- p->entropy =
- hydro_gamma_minus_one * p->entropy * pow_minus_gamma_minus_one(p->rho);
-}
-
-/**
- * @brief Returns the internal energy of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
- const struct part *restrict p, float dt) {
-
- const float entropy = p->entropy + p->entropy_dt * dt;
-
- return entropy * pow_gamma_minus_one(p->rho) * hydro_one_over_gamma_minus_one;
-}
-
-/**
- * @brief Returns the pressure of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_pressure(
- const struct part *restrict p, float dt) {
-
- return p->force.P_over_rho2 * p->rho * p->rho / p->rho_dh;
-}
-
-/**
- * @brief Returns the entropy of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_entropy(
- const struct part *restrict p, float dt) {
-
- return p->entropy + p->entropy_dt * dt;
-}
-
-/**
- * @brief Returns the sound speed of a particle
- *
- * @param p The particle of interest
- * @param dt Time since the last kick
- */
-__attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
- const struct part *restrict p, float dt) {
-
- return p->force.soundspeed;
+ /* We read u in the entropy field. We now get S from u */
+ p->entropy = gas_entropy_from_internal_energy(p->rho, p->entropy);
}
diff --git a/src/timestep.h b/src/timestep.h
index 99747484b2..569120cf9c 100644
--- a/src/timestep.h
+++ b/src/timestep.h
@@ -108,7 +108,7 @@ __attribute__((always_inline)) INLINE static int get_part_timestep(
e->external_potential, e->physical_constants, p->gpart);
/* const float new_dt_self = */
/* gravity_compute_timestep_self(e->physical_constants, p->gpart); */
- const float new_dt_self = FLT_MAX; // MATTHIEU
+ const float new_dt_self = FLT_MAX; // MATTHIEU
new_dt_grav = fminf(new_dt_external, new_dt_self);
}
--
GitLab