From c4640e79dbbfc4bc8441651c3148fd40fb60c508 Mon Sep 17 00:00:00 2001
From: Bert Vandenbroucke <bert.vandenbroucke@ugent.be>
Date: Wed, 1 Mar 2017 15:53:29 +0000
Subject: [PATCH] Improved initial density guess in isothermal Riemann solver,
which leads to (faster) convergence. Added some extra isothermal eos flags in
Gizmo hydro.h. Temporarily disabled primitive variable drift and changed the
way the actual particle movement is calculated.
---
src/hydro/Gizmo/hydro.h | 50 +++++++++++++++++---------
src/riemann/riemann_exact_isothermal.h | 25 ++++++++++---
2 files changed, 55 insertions(+), 20 deletions(-)
diff --git a/src/hydro/Gizmo/hydro.h b/src/hydro/Gizmo/hydro.h
index 7cfb29fa19..d322c493e0 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 8e8a6a4dcc..34471593e2 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]);
}
}
}
--
GitLab