Skip to content
Snippets Groups Projects
Commit c4640e79 authored by Bert Vandenbroucke's avatar Bert Vandenbroucke
Browse files

Improved initial density guess in isothermal Riemann solver, which leads to...

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.
parent b7424e3b
No related branches found
No related tags found
1 merge request!317Cleaned up GIZMO code, added SineWavePotential tests.
...@@ -109,8 +109,12 @@ __attribute__((always_inline)) INLINE static void hydro_first_init_part( ...@@ -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[1] = mass * p->primitives.v[1];
p->conserved.momentum[2] = mass * p->primitives.v[2]; 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; p->conserved.energy *= mass;
#endif
#if defined(GIZMO_FIX_PARTICLES) #if defined(GIZMO_FIX_PARTICLES)
p->v[0] = 0.; p->v[0] = 0.;
...@@ -350,21 +354,24 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra( ...@@ -350,21 +354,24 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
else else
p->h *= expf(w1); p->h *= expf(w1);
const float w2 = -hydro_dimension * w1; // const float w2 = -hydro_dimension * w1;
if (fabsf(w2) < 0.2f) { // if (fabsf(w2) < 0.2f) {
p->primitives.rho *= approx_expf(w2); // p->primitives.rho *= approx_expf(w2);
} else { // } else {
p->primitives.rho *= expf(w2); // p->primitives.rho *= expf(w2);
} // }
p->primitives.v[0] += (p->a_hydro[0] + p->gravity.old_a[0]) * dt; // 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[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; // 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; //#if !defined(EOS_ISOTHERMAL_GAS)
p->primitives.P = // if (p->conserved.mass > 0.) {
hydro_gamma_minus_one * u * p->primitives.rho / p->conserved.mass; // 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 #ifdef SWIFT_DEBUG_CHECKS
if (p->h <= 0.) { if (p->h <= 0.) {
...@@ -460,7 +467,11 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( ...@@ -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[0] += p->conserved.flux.momentum[0];
p->conserved.momentum[1] += p->conserved.flux.momentum[1]; p->conserved.momentum[1] += p->conserved.flux.momentum[1];
p->conserved.momentum[2] += p->conserved.flux.momentum[2]; 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; p->conserved.energy += p->conserved.flux.energy;
#endif
/* Add gravity. We only do this if we have gravity activated. */ /* Add gravity. We only do this if we have gravity activated. */
if (p->gpart) { if (p->gpart) {
...@@ -475,6 +486,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( ...@@ -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[1] += dt * p->conserved.mass * a_grav[1];
p->conserved.momentum[2] += dt * p->conserved.mass * a_grav[2]; 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.energy += dt * (p->conserved.momentum[0] * a_grav[0] +
p->conserved.momentum[1] * a_grav[1] + p->conserved.momentum[1] * a_grav[1] +
p->conserved.momentum[2] * a_grav[2]); p->conserved.momentum[2] * a_grav[2]);
...@@ -482,6 +494,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( ...@@ -482,6 +494,7 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->conserved.energy += dt * (a_grav[0] * p->gravity.mflux[0] + p->conserved.energy += dt * (a_grav[0] * p->gravity.mflux[0] +
a_grav[1] * p->gravity.mflux[1] + a_grav[1] * p->gravity.mflux[1] +
a_grav[2] * p->gravity.mflux[2]); a_grav[2] * p->gravity.mflux[2]);
#endif
} }
/* reset fluxes */ /* reset fluxes */
...@@ -492,6 +505,11 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra( ...@@ -492,6 +505,11 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
p->conserved.flux.momentum[1] = 0.0f; p->conserved.flux.momentum[1] = 0.0f;
p->conserved.flux.momentum[2] = 0.0f; p->conserved.flux.momentum[2] = 0.0f;
p->conserved.flux.energy = 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;
} }
/** /**
......
...@@ -112,7 +112,7 @@ __attribute__((always_inline)) INLINE static float riemann_guess_rho(float* WL, ...@@ -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 /* Currently three possibilities and not really an algorithm to decide which
one to choose: */ one to choose: */
/* just the average */ /* just the average */
return 0.5f * (WL[0] + WR[0]); // return 0.5f * (WL[0] + WR[0]);
/* two rarefaction approximation */ /* two rarefaction approximation */
return sqrtf(WL[0] * WR[0] * expf((vL - vR) / const_isothermal_soundspeed)); 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( ...@@ -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]; 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]; 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.; rho = 0.;
/* obtain a first guess for p */ /* obtain a first guess for p */
rhoguess = riemann_guess_rho(WL, WR, vL, vR); 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); frho = riemann_f(rho, WL, WR, vL, vR);
frhoguess = riemann_f(rhoguess, WL, WR, vL, vR); frhoguess = riemann_f(rhoguess, WL, WR, vL, vR);
/* ok, rhostar is close to 0, better use Brent's method... */ /* ok, rhostar is close to 0, better use Brent's method... */
...@@ -289,14 +302,18 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve( ...@@ -289,14 +302,18 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
/* Newton-Raphson until convergence or until suitable interval is found /* Newton-Raphson until convergence or until suitable interval is found
to use Brent's method */ to use Brent's method */
unsigned int counter = 0; 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) { frhoguess < 0.0f) {
rho = rhoguess; rho = rhoguess;
rhoguess = rhoguess - frhoguess / riemann_fprime(rhoguess, WL, WR); rhoguess = rhoguess - frhoguess / riemann_fprime(rhoguess, WL, WR);
frhoguess = riemann_f(rhoguess, WL, WR, vL, vR); frhoguess = riemann_f(rhoguess, WL, WR, vL, vR);
counter++; counter++;
if (counter > 1000) { 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]);
} }
} }
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment