Added some extra isothermal eos flags in Gizmo/hydro.h. Made sure gpart mass...

Added some extra isothermal eos flags in Gizmo/hydro.h. Made sure gpart mass is updated when conserved mass changes. Fixed a number of possible issues with the isothermal Riemann solver.

src/const.h

@@ -58,7 +58,7 @@
 /* Options to control the movement of particles for GIZMO_SPH. */
 /* This option disables particle movement */
 //#define GIZMO_FIX_PARTICLES
-#define GIZMO_TOTAL_ENERGY
+//#define GIZMO_TOTAL_ENERGY
 
 /* Source terms */
 #define SOURCETERMS_NONE

src/hydro/Gizmo/hydro.h

@@ -23,6 +23,7 @@
 #include "approx_math.h"
 #include "equation_of_state.h"
 #include "hydro_gradients.h"
+#include "riemann.h"
 #include "minmax.h"
 
 /**
@@ -238,6 +239,10 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
   p->primitives.v[1] = momentum[1] / m;
   p->primitives.v[2] = momentum[2] / m;
 
+#ifdef EOS_ISOTHERMAL_GAS
+  p->primitives.P = const_isothermal_soundspeed*const_isothermal_soundspeed*p->primitives.rho;
+#else
+
   float energy = p->conserved.energy;
 
 #ifdef GIZMO_TOTAL_ENERGY
@@ -247,6 +252,7 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
 #endif
 
   p->primitives.P = hydro_gamma_minus_one * energy / volume;
+#endif
 
   /* sanity checks */
   if (p->primitives.rho < 0.0f || p->primitives.P < 0.0f) {
@@ -278,13 +284,9 @@ __attribute__((always_inline)) INLINE static void hydro_prepare_force(
   p->timestepvars.vmax = 0.0f;
 
   /* Set the actual velocity of the particle */
-  /* This should be v and not v_full. The reason is that v is updated with the
-     acceleration times the full time step, while v_full is updated using a
-     mixture of the half time steps that don't necessarily correspond to the
-     full time step. */
-  p->force.v_full[0] = p->v[0];
-  p->force.v_full[1] = p->v[1];
-  p->force.v_full[2] = p->v[2];
+  p->force.v_full[0] = xp->v_full[0];
+  p->force.v_full[1] = xp->v_full[1];
+  p->force.v_full[2] = xp->v_full[2];
 }
 
 /**
@@ -506,12 +508,20 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
     a_grav[1] = p->gpart->a_grav[1];
     a_grav[2] = p->gpart->a_grav[2];
 
+    /* Store the gravitational acceleration for later use. */
+    p->gravity.old_a[0] = a_grav[0];
+    p->gravity.old_a[1] = a_grav[1];
+    p->gravity.old_a[2] = a_grav[2];
+
+    /* Make sure the gpart knows the mass has changed. */
+    p->gpart->mass = p->conserved.mass;
+
     /* Kick the momentum for half a time step */
     p->conserved.momentum[0] += dt * p->conserved.mass * a_grav[0];
     p->conserved.momentum[1] += dt * p->conserved.mass * a_grav[1];
     p->conserved.momentum[2] += dt * p->conserved.mass * a_grav[2];
 
-#if !defined(EOS_ISOTHERMAL_GAS)
+#if !defined(EOS_ISOTHERMAL_GAS) && defined(GIZMO_TOTAL_ENERGY)
     p->conserved.energy += dt * (p->conserved.momentum[0] * a_grav[0] +
                                  p->conserved.momentum[1] * a_grav[1] +
                                  p->conserved.momentum[2] * a_grav[2]);
@@ -578,7 +588,11 @@ __attribute__((always_inline)) INLINE static float hydro_get_internal_energy(
     const struct part* restrict p) {
 
   if (p->primitives.rho > 0.) {
+#ifdef EOS_ISOTHERMAL_GAS
     return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
+#else
+    return const_isothermal_internal_energy;
+#endif
   } else {
     return 0.;
   }
@@ -608,7 +622,11 @@ __attribute__((always_inline)) INLINE static float hydro_get_soundspeed(
     const struct part* restrict p) {
 
   if (p->primitives.rho > 0.) {
+#ifdef EOS_ISOTHERMAL_GAS
+    return sqrtf(const_isothermal_internal_energy * hydro_gamma * hydro_gamma_minus_one);
+#else
     return sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
+#endif
   } else {
     return 0.;
   }

src/hydro/Gizmo/hydro_gradients.h

@@ -22,6 +22,7 @@
 #define SWIFT_HYDRO_GRADIENTS_H
 
 #include "hydro_slope_limiters.h"
+#include "riemann.h"
 
 #if defined(GRADIENTS_SPH)
 
@@ -142,6 +143,27 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
 
   /* time */
   if (Wi[0] > 0.0f) {
+#ifdef EOS_ISOTHERMAL_GAS
+    dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
+                             Wi[2] * pi->primitives.gradients.rho[1] +
+                             Wi[3] * pi->primitives.gradients.rho[2] +
+                             Wi[0] * (pi->primitives.gradients.v[0][0] +
+                                      pi->primitives.gradients.v[1][1] +
+                                      pi->primitives.gradients.v[2][2]));
+    dWi[1] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[0][0] +
+                             Wi[2] * pi->primitives.gradients.v[0][1] +
+                             Wi[3] * pi->primitives.gradients.v[0][2] +
+                             const_isothermal_soundspeed*const_isothermal_soundspeed*pi->primitives.gradients.rho[0] / Wi[0]);
+    dWi[2] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[1][0] +
+                             Wi[2] * pi->primitives.gradients.v[1][1] +
+                             Wi[3] * pi->primitives.gradients.v[1][2] +
+                             const_isothermal_soundspeed*const_isothermal_soundspeed*pi->primitives.gradients.rho[1] / Wi[0]);
+    dWi[3] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.v[2][0] +
+                             Wi[2] * pi->primitives.gradients.v[2][1] +
+                             Wi[3] * pi->primitives.gradients.v[2][2] +
+                             const_isothermal_soundspeed*const_isothermal_soundspeed*pi->primitives.gradients.rho[2] / Wi[0]);
+    /* we don't care about P in this case */
+#else
     dWi[0] -= 0.5 * mindt * (Wi[1] * pi->primitives.gradients.rho[0] +
                              Wi[2] * pi->primitives.gradients.rho[1] +
                              Wi[3] * pi->primitives.gradients.rho[2] +
@@ -167,9 +189,30 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
                hydro_gamma * Wi[4] * (pi->primitives.gradients.v[0][0] +
                                       pi->primitives.gradients.v[1][1] +
                                       pi->primitives.gradients.v[2][2]));
+#endif
   }
 
   if (Wj[0] > 0.0f) {
+#ifdef EOS_ISOTHERMAL_GAS
+    dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
+                             Wj[2] * pj->primitives.gradients.rho[1] +
+                             Wj[3] * pj->primitives.gradients.rho[2] +
+                             Wj[0] * (pj->primitives.gradients.v[0][0] +
+                                      pj->primitives.gradients.v[1][1] +
+                                      pj->primitives.gradients.v[2][2]));
+    dWj[1] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[0][0] +
+                             Wj[2] * pj->primitives.gradients.v[0][1] +
+                             Wj[3] * pj->primitives.gradients.v[0][2] +
+                             const_isothermal_soundspeed*const_isothermal_soundspeed*pj->primitives.gradients.rho[0] / Wj[0]);
+    dWj[2] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[1][0] +
+                             Wj[2] * pj->primitives.gradients.v[1][1] +
+                             Wj[3] * pj->primitives.gradients.v[1][2] +
+                             const_isothermal_soundspeed*const_isothermal_soundspeed*pj->primitives.gradients.rho[1] / Wj[0]);
+    dWj[3] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.v[2][0] +
+                             Wj[2] * pj->primitives.gradients.v[2][1] +
+                             Wj[3] * pj->primitives.gradients.v[2][2] +
+                             const_isothermal_soundspeed*const_isothermal_soundspeed*pj->primitives.gradients.rho[2] / Wj[0]);
+#else
     dWj[0] -= 0.5 * mindt * (Wj[1] * pj->primitives.gradients.rho[0] +
                              Wj[2] * pj->primitives.gradients.rho[1] +
                              Wj[3] * pj->primitives.gradients.rho[2] +
@@ -195,6 +238,7 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
                hydro_gamma * Wj[4] * (pj->primitives.gradients.v[0][0] +
                                       pj->primitives.gradients.v[1][1] +
                                       pj->primitives.gradients.v[2][2]));
+#endif
   }
 
   if (-dWi[0] > Wi[0]) {

src/hydro/Gizmo/hydro_iact.h

@@ -232,8 +232,12 @@ __attribute__((always_inline)) INLINE static void runner_iact_fluxes_common(
 
   /* calculate the maximal signal velocity */
   if (Wi[0] > 0.0f && Wj[0] > 0.0f) {
+#ifdef EOS_ISOTHERMAL_GAS
+    vmax = 2.*const_isothermal_soundspeed;
+#else
     vmax =
         sqrtf(hydro_gamma * Wi[4] / Wi[0]) + sqrtf(hydro_gamma * Wj[4] / Wj[0]);
+#endif
   } else {
     vmax = 0.0f;
   }

src/hydro/Gizmo/hydro_io.h

@@ -70,7 +70,15 @@ void hydro_read_particles(struct part* parts, struct io_props* list,
  * @return Internal energy of the particle
  */
 float convert_u(struct engine* e, struct part* p) {
-  return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
+  if(p->primitives.rho > 0.){
+#ifdef EOS_ISOTHERMAL_GAS
+    return const_isothermal_internal_energy;
+#else
+    return p->primitives.P / hydro_gamma_minus_one / p->primitives.rho;
+#endif
+  } else {
+    return 0.;
+  }
 }
 
 /**
@@ -81,7 +89,11 @@ float convert_u(struct engine* e, struct part* p) {
  * @return Entropic function of the particle
  */
 float convert_A(struct engine* e, struct part* p) {
+  if(p->primitives.rho > 0.){
   return p->primitives.P / pow_gamma(p->primitives.rho);
+  } else  {
+    return 0.;
+  }
 }
 
 /**
@@ -92,6 +104,9 @@ float convert_A(struct engine* e, struct part* p) {
  * @return Total energy of the particle
  */
 float convert_Etot(struct engine* e, struct part* p) {
+#ifdef GIZMO_TOTAL_ENERGY
+  return p->conserved.energy;
+#else
   float momentum2;
 
   momentum2 = p->conserved.momentum[0] * p->conserved.momentum[0] +
@@ -99,6 +114,7 @@ float convert_Etot(struct engine* e, struct part* p) {
               p->conserved.momentum[2] * p->conserved.momentum[2];
 
   return p->conserved.energy + 0.5f * momentum2 / p->conserved.mass;
+#endif
 }
 
 /**
@@ -111,7 +127,7 @@ float convert_Etot(struct engine* e, struct part* p) {
 void hydro_write_particles(struct part* parts, struct io_props* list,
                            int* num_fields) {
 
-  *num_fields = 13;
+  *num_fields = 14;
 
   /* List what we want to write */
   list[0] = io_make_output_field("Coordinates", DOUBLE, 3, UNIT_CONV_LENGTH,
@@ -142,6 +158,9 @@ void hydro_write_particles(struct part* parts, struct io_props* list,
   list[12] =
       io_make_output_field_convert_part("TotEnergy", FLOAT, 1, UNIT_CONV_ENERGY,
                                         parts, conserved.energy, convert_Etot);
+  list[13] = io_make_output_field("GravAcceleration", FLOAT, 3, UNIT_CONV_ACCELERATION,
+                                 parts, gravity.old_a);
+
 }
 
 /**

src/riemann/riemann_exact_isothermal.h

@@ -292,7 +292,8 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
   }
 #endif
 
-  frho = riemann_f(rho, WL, WR, vL, vR);
+  /* we know that the value of riemann_f for rho=0 is negative (it is -inf). */
+  frho = -1.;
   frhoguess = riemann_f(rhoguess, WL, WR, vL, vR);
   /* ok, rhostar is close to 0, better use Brent's method... */
   /* we use Newton-Raphson until we find a suitable interval */
@@ -319,7 +320,7 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
   if (1.e6 * fabs(rho - rhoguess) > 0.5f * (rho + rhoguess) &&
       frhoguess > 0.0f) {
     rho = 0.0f;
-    frho = riemann_f(rho, WL, WR, vL, vR);
+    frho = -1.;
     /* use Brent's method to find the zeropoint */
     rho = riemann_solve_brent(rho, rhoguess, frho, frhoguess, 1.e-6, WL, WR, vL,
                               vR);