From d6d8f6781ece0d70ab678ce239c836c29b383b88 Mon Sep 17 00:00:00 2001
From: Bert Vandenbroucke <bert.vandenbroucke@gmail.com>
Date: Wed, 28 Nov 2018 16:48:33 +0000
Subject: [PATCH] Redisabled total energy evolution in Gizmo. Added entropy
 evolution to GizmoMFM and temporarily overwrote energy evolution. Still need
 to add an appropriate switch.

---
 src/const.h                          |  2 +-
 src/hydro/GizmoMFM/hydro.h           | 14 ++++++++++++++
 src/hydro/GizmoMFM/hydro_gradients.h |  2 +-
 src/hydro/GizmoMFM/hydro_part.h      |  6 ++++++
 4 files changed, 22 insertions(+), 2 deletions(-)

diff --git a/src/const.h b/src/const.h
index 737a6870fc..f41034101b 100644
--- a/src/const.h
+++ b/src/const.h
@@ -45,7 +45,7 @@
 /* Try to keep cells regular by adding a correction velocity. */
 //#define GIZMO_STEER_MOTION
 /* Use the total energy instead of the thermal energy as conserved variable. */
-#define GIZMO_TOTAL_ENERGY
+//#define GIZMO_TOTAL_ENERGY
 
 /* Options to control handling of unphysical values (GIZMO_SPH only). */
 /* In GIZMO, mass and energy (and hence density and pressure) can in principle
diff --git a/src/hydro/GizmoMFM/hydro.h b/src/hydro/GizmoMFM/hydro.h
index d57fd61ed9..dad33486ed 100644
--- a/src/hydro/GizmoMFM/hydro.h
+++ b/src/hydro/GizmoMFM/hydro.h
@@ -156,6 +156,10 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
      is to have a way of remembering that we need more neighbours for this
      particle */
   p->geometry.wcorr = 1.0f;
+
+  /* set the entropy to a non-zero initial value (to make sure we
+     don't zero out the density before we reach hydro_convert_quantities) */
+  p->conserved.entropy = 1.0f;
 }
 
 /**
@@ -334,6 +338,8 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
   /* energy contains the total thermal energy, we want the specific energy.
      this is why we divide by the volume, and not by the density */
   p->P = hydro_gamma_minus_one * energy * volume_inv;
+  p->A = p->conserved.entropy / m;
+  p->P = p->A * pow_gamma(p->rho);
 #endif
 
   /* sanity checks */
@@ -518,6 +524,11 @@ __attribute__((always_inline)) INLINE static void hydro_convert_quantities(
     const struct hydro_props* hydro_props) {
 
   p->conserved.energy /= cosmo->a_factor_internal_energy;
+
+  /* initialize the entropy */
+  p->conserved.entropy =
+    p->conserved.energy * pow_minus_gamma_minus_one(p->rho) *
+      hydro_gamma_minus_one;
 }
 
 /**
@@ -570,6 +581,7 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
     const float u = (p->conserved.energy + p->flux.energy * dt_drift) * m_inv;
 #endif
     p->P = hydro_gamma_minus_one * u * p->rho;
+    p->P = p->A * pow_gamma(p->rho);
 #endif
   }
 
@@ -635,12 +647,14 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
   p->conserved.momentum[0] += p->flux.momentum[0] * dt_therm;
   p->conserved.momentum[1] += p->flux.momentum[1] * dt_therm;
   p->conserved.momentum[2] += p->flux.momentum[2] * dt_therm;
+
 #if defined(EOS_ISOTHERMAL_GAS)
   /* We use the EoS equation in a sneaky way here just to get the constant u */
   p->conserved.energy =
       p->conserved.mass * gas_internal_energy_from_entropy(0.0f, 0.0f);
 #else
   p->conserved.energy += p->flux.energy * dt_therm;
+  p->conserved.energy = hydro_one_over_gamma_minus_one * p->conserved.entropy * pow_gamma_minus_one(p->rho);
 #endif
 
 #ifndef HYDRO_GAMMA_5_3
diff --git a/src/hydro/GizmoMFM/hydro_gradients.h b/src/hydro/GizmoMFM/hydro_gradients.h
index 6f751d9702..805cebf68d 100644
--- a/src/hydro/GizmoMFM/hydro_gradients.h
+++ b/src/hydro/GizmoMFM/hydro_gradients.h
@@ -95,7 +95,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
     struct part* restrict pi, struct part* restrict pj, float hi, float hj,
     const float* dx, float r, const float* xij_i, float* Wi, float* Wj) {
 
-  /* perform gradient reconstruction in space and time */
+  /* perform gradient reconstruction in space */
   /* Compute interface position (relative to pj, since we don't need the actual
    * position) eqn. (8) */
   const float xij_j[3] = {xij_i[0] + dx[0], xij_i[1] + dx[1], xij_i[2] + dx[2]};
diff --git a/src/hydro/GizmoMFM/hydro_part.h b/src/hydro/GizmoMFM/hydro_part.h
index 8097b1b256..0adfed9952 100644
--- a/src/hydro/GizmoMFM/hydro_part.h
+++ b/src/hydro/GizmoMFM/hydro_part.h
@@ -77,6 +77,9 @@ struct part {
   /* Pressure. */
   float P;
 
+  /* Entropic function. */
+  float A;
+
   union {
     /* Quantities used during the volume (=density) loop. */
     struct {
@@ -169,6 +172,9 @@ struct part {
     /* Fluid thermal energy (not per unit mass!). */
     float energy;
 
+    /* Fluid entropy. */
+    float entropy;
+
   } conserved;
 
   /* Geometrical quantities used for hydro. */
-- 
GitLab