diff --git a/src/const.h b/src/const.h
index 5c95bd27c13487208e2b12837ec49e6480633363..010cad4167c28c93a69d9e23f7dc8612462c20bf 100644
--- a/src/const.h
+++ b/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
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 165d4b54fb4fa55853479dc6a081c0830473c6fc..42a10cae7f6c8a42db629d51ce3975a799fe1b84 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/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.;
}
diff --git a/src/hydro/Gizmo/hydro_gradients.h b/src/hydro/Gizmo/hydro_gradients.h
index e52c39a1f43defc72a2f42d2bb742e262f5b701d..ac99f072c88846a17d8a746a9659ca829f1a91eb 100644
--- a/src/hydro/Gizmo/hydro_gradients.h
+++ b/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]) {
diff --git a/src/hydro/Gizmo/hydro_iact.h b/src/hydro/Gizmo/hydro_iact.h
index 342b59c695df304b53deb52bb33c27ad375fd6f5..c0ce7eca0fa4955b9f921adb9cbca5f4fa6e1d4c 100644
--- a/src/hydro/Gizmo/hydro_iact.h
+++ b/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;
}
diff --git a/src/hydro/Gizmo/hydro_io.h b/src/hydro/Gizmo/hydro_io.h
index 818217baee18c2128c8f004e3f2b14ecd454acfc..3e33d6c69876b06d1c830dcfceb1ac3a44284913 100644
--- a/src/hydro/Gizmo/hydro_io.h
+++ b/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);
+
}
/**
diff --git a/src/riemann/riemann_exact_isothermal.h b/src/riemann/riemann_exact_isothermal.h
index f8c6837bcf7938e757b45d4d828464b692426210..26fb649f02ddd62c102e31c9191c20f9003c0133 100644
--- a/src/riemann/riemann_exact_isothermal.h
+++ b/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);