diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 7cfb29fa196bf40ed83bc1340cae807f9ccb4ec9..d322c493e01d7f89db5d7eec3db9974bdc5750c8 100644
--- a/src/hydro/Gizmo/hydro.h
+++ b/src/hydro/Gizmo/hydro.h
@@ -109,8 +109,12 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part(
p->conserved.momentum[1] = mass * p->primitives.v[1];
p->conserved.momentum[2] = mass * p->primitives.v[2];
- /* and the thermal energy */
+/* and the thermal energy */
+#if defined(EOS_ISOTHERMAL_GAS)
+ p->conserved.energy = mass * const_isothermal_internal_energy;
+#else
p->conserved.energy *= mass;
+#endif
#if defined(GIZMO_FIX_PARTICLES)
p->v[0] = 0.;
@@ -350,21 +354,24 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
else
p->h *= expf(w1);
- const float w2 = -hydro_dimension * w1;
- if (fabsf(w2) < 0.2f) {
- p->primitives.rho *= approx_expf(w2);
- } else {
- p->primitives.rho *= expf(w2);
- }
-
- p->primitives.v[0] += (p->a_hydro[0] + p->gravity.old_a[0]) * dt;
- p->primitives.v[1] += (p->a_hydro[1] + p->gravity.old_a[1]) * dt;
- p->primitives.v[2] += (p->a_hydro[2] + p->gravity.old_a[2]) * dt;
- if (p->conserved.mass > 0.) {
- const float u = p->conserved.energy + p->du_dt * dt;
- p->primitives.P =
- hydro_gamma_minus_one * u * p->primitives.rho / p->conserved.mass;
- }
+// const float w2 = -hydro_dimension * w1;
+// if (fabsf(w2) < 0.2f) {
+// p->primitives.rho *= approx_expf(w2);
+// } else {
+// p->primitives.rho *= expf(w2);
+// }
+
+// p->primitives.v[0] += (p->a_hydro[0] + p->gravity.old_a[0]) * dt;
+// p->primitives.v[1] += (p->a_hydro[1] + p->gravity.old_a[1]) * dt;
+// p->primitives.v[2] += (p->a_hydro[2] + p->gravity.old_a[2]) * dt;
+
+//#if !defined(EOS_ISOTHERMAL_GAS)
+// if (p->conserved.mass > 0.) {
+// const float u = p->conserved.energy + p->du_dt * dt;
+// p->primitives.P =
+// hydro_gamma_minus_one * u * p->primitives.rho / p->conserved.mass;
+// }
+//#endif
#ifdef SWIFT_DEBUG_CHECKS
if (p->h <= 0.) {
@@ -460,7 +467,11 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->conserved.momentum[0] += p->conserved.flux.momentum[0];
p->conserved.momentum[1] += p->conserved.flux.momentum[1];
p->conserved.momentum[2] += p->conserved.flux.momentum[2];
+#if defined(EOS_ISOTHERMAL_GAS)
+ p->conserved.energy = p->conserved.mass * const_isothermal_internal_energy;
+#else
p->conserved.energy += p->conserved.flux.energy;
+#endif
/* Add gravity. We only do this if we have gravity activated. */
if (p->gpart) {
@@ -475,6 +486,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
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)
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]);
@@ -482,6 +494,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->conserved.energy += dt * (a_grav[0] * p->gravity.mflux[0] +
a_grav[1] * p->gravity.mflux[1] +
a_grav[2] * p->gravity.mflux[2]);
+#endif
}
/* reset fluxes */
@@ -492,6 +505,11 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->conserved.flux.momentum[1] = 0.0f;
p->conserved.flux.momentum[2] = 0.0f;
p->conserved.flux.energy = 0.0f;
+
+ /* Set particle movement */
+ xp->v_full[0] = p->conserved.momentum[0] / p->conserved.mass;
+ xp->v_full[1] = p->conserved.momentum[1] / p->conserved.mass;
+ xp->v_full[2] = p->conserved.momentum[2] / p->conserved.mass;
}
/**
diff --git a/src/riemann/riemann_exact_isothermal.h b/src/riemann/riemann_exact_isothermal.h
index 8e8a6a4dcce456b55f85148bed7342ad4e651c1b..34471593e2293e653ee130835aaf3548a9020562 100644
--- a/src/riemann/riemann_exact_isothermal.h
+++ b/src/riemann/riemann_exact_isothermal.h
@@ -112,7 +112,7 @@ __attribute__((always_inline)) INLINE static float riemann_guess_rho(float* WL,
/* Currently three possibilities and not really an algorithm to decide which
one to choose: */
/* just the average */
- return 0.5f * (WL[0] + WR[0]);
+ // return 0.5f * (WL[0] + WR[0]);
/* two rarefaction approximation */
return sqrtf(WL[0] * WR[0] * expf((vL - vR) / const_isothermal_soundspeed));
@@ -276,11 +276,24 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
vL = WL[1] * n_unit[0] + WL[2] * n_unit[1] + WL[3] * n_unit[2];
vR = WR[1] * n_unit[0] + WR[2] * n_unit[1] + WR[3] * n_unit[2];
- /* VACUUM... */
+/* VACUUM... */
+#ifdef SWIFT_DEBUG_CHECKS
+ if (WL[0] == 0. || WL[4] == 0. || WR[0] == 0. || WR[4] == 0. ||
+ (4.0f * const_isothermal_soundspeed / hydro_gamma_minus_one <= vR - vL)) {
+ error("Vacuum not handled (yet)!");
+ }
+#endif
rho = 0.;
/* obtain a first guess for p */
rhoguess = riemann_guess_rho(WL, WR, vL, vR);
+
+#ifdef SWIFT_DEBUG_CHECKS
+ if (rhoguess <= 0.) {
+ error("Zero or negative initial density guess.");
+ }
+#endif
+
frho = riemann_f(rho, WL, WR, vL, vR);
frhoguess = riemann_f(rhoguess, WL, WR, vL, vR);
/* ok, rhostar is close to 0, better use Brent's method... */
@@ -289,14 +302,18 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
/* Newton-Raphson until convergence or until suitable interval is found
to use Brent's method */
unsigned int counter = 0;
- while (fabs(rho - rhoguess) > 1.e-6f * 0.5f * (rho + rhoguess) &&
+ while (fabs(rho - rhoguess) > 5.e-7f * (rho + rhoguess) &&
frhoguess < 0.0f) {
rho = rhoguess;
rhoguess = rhoguess - frhoguess / riemann_fprime(rhoguess, WL, WR);
frhoguess = riemann_f(rhoguess, WL, WR, vL, vR);
counter++;
if (counter > 1000) {
- error("Stuck in Newton-Raphson!\n");
+ error(
+ "Stuck in Newton-Raphson (rho: %g, rhoguess: %g, frhoguess: %g, "
+ "fprime: %g, rho-rhoguess: %g, WL: %g %g %g, WR: %g %g %g)!\n",
+ rho, rhoguess, frhoguess, riemann_fprime(rhoguess, WL, WR),
+ (rho - rhoguess), WL[0], vL, WL[4], WR[0], vR, WR[4]);
}
}
}