Commit b403c0fa authored by Bert Vandenbroucke's avatar Bert Vandenbroucke Committed by Matthieu Schaller
Browse files

Removed unused isothermal Riemann solver. Added generic input and output debug...

Removed unused isothermal Riemann solver. Added generic input and output debug checks for all remaining Riemann solvers. Made sure HLLC Riemann solver works for 0 input pressure. Made sure all quantities are set to vacuum if density/mass or pressure/energy is zero.
parent 35609c96
......@@ -110,6 +110,7 @@ nobase_noinst_HEADERS = align.h approx_math.h atomic.h barrier.h cycle.h error.h
hydro/Shadowswift/voronoi_cell.h \
riemann/riemann_hllc.h riemann/riemann_trrs.h \
riemann/riemann_exact.h riemann/riemann_vacuum.h \
riemann/riemann_checks.h \
stars.h stars_io.h \
stars/Default/star.h stars/Default/star_iact.h stars/Default/star_io.h \
stars/Default/star_debug.h stars/Default/star_part.h \
......
......@@ -362,8 +362,9 @@ __attribute__((always_inline)) INLINE static void hydro_end_density(
#endif
/* sanity checks */
gizmo_check_physical_quantity("density", p->primitives.rho);
gizmo_check_physical_quantity("pressure", p->primitives.P);
gizmo_check_physical_quantities("density", "pressure", p->primitives.rho,
p->primitives.v[0], p->primitives.v[1],
p->primitives.v[2], p->primitives.P);
#ifdef GIZMO_LLOYD_ITERATION
/* overwrite primitive variables to make sure they still have safe values */
......@@ -563,7 +564,6 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
p->primitives.v[2] +=
p->conserved.flux.momentum[2] * dt_drift / p->conserved.mass;
// MATTHIEU: Bert is this correct?
#if !defined(EOS_ISOTHERMAL_GAS)
const float u = p->conserved.energy + p->conserved.flux.energy * dt_therm;
p->primitives.P =
......@@ -576,6 +576,10 @@ __attribute__((always_inline)) INLINE static void hydro_predict_extra(
error("Zero or negative smoothing length (%g)!", p->h);
}
#endif
gizmo_check_physical_quantities("density", "pressure", p->primitives.rho,
p->primitives.v[0], p->primitives.v[1],
p->primitives.v[2], p->primitives.P);
}
/**
......@@ -632,13 +636,14 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
/* Apply the minimal energy limit */
const float min_energy =
hydro_props->minimal_internal_energy * cosmo->a_factor_internal_energy;
if (p->conserved.energy < min_energy) {
p->conserved.energy = min_energy;
if (p->conserved.energy < min_energy * p->conserved.mass) {
p->conserved.energy = min_energy * p->conserved.mass;
p->conserved.flux.energy = 0.f;
}
gizmo_check_physical_quantity("mass", p->conserved.mass);
gizmo_check_physical_quantity("energy", p->conserved.energy);
gizmo_check_physical_quantities(
"mass", "energy", p->conserved.mass, p->conserved.momentum[0],
p->conserved.momentum[1], p->conserved.momentum[2], p->conserved.energy);
#ifdef SWIFT_DEBUG_CHECKS
/* Note that this check will only have effect if no GIZMO_UNPHYSICAL option
......@@ -675,6 +680,10 @@ __attribute__((always_inline)) INLINE static void hydro_kick_extra(
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];
p->conserved.energy += dt * (p->gravity.mflux[0] * a_grav[0] +
p->gravity.mflux[1] * a_grav[1] +
p->gravity.mflux[2] * a_grav[2]);
}
hydro_velocities_set(p, xp);
......
......@@ -150,11 +150,10 @@ __attribute__((always_inline)) INLINE static void hydro_gradients_predict(
Wj[3] += dWj[3];
Wj[4] += dWj[4];
/* Check that we don't have problematic densities or pressures */
gizmo_check_physical_quantity("density", Wi[0]);
gizmo_check_physical_quantity("pressure", Wi[4]);
gizmo_check_physical_quantity("density", Wj[0]);
gizmo_check_physical_quantity("pressure", Wj[4]);
gizmo_check_physical_quantities("density", "pressure", Wi[0], Wi[1], Wi[2],
Wi[3], Wi[4]);
gizmo_check_physical_quantities("density", "pressure", Wj[0], Wj[1], Wj[2],
Wj[3], Wj[4]);
}
#endif /* SWIFT_HYDRO_GIZMO_GRADIENTS_H */
......@@ -46,9 +46,23 @@
quantity = 0.f; \
}
#define gizmo_check_physical_quantities( \
mass_name, energy_name, mass, momentum_x, momentum_y, momentum_z, energy) \
gizmo_check_physical_quantity(mass_name, mass); \
gizmo_check_physical_quantity(energy_name, energy); \
/* now check for vacuum and make sure we have a real vacuum */ \
if (mass == 0.f || energy == 0.f) { \
mass = 0.f; \
momentum_x = 0.f; \
momentum_y = 0.f; \
momentum_z = 0.f; \
energy = 0.f; \
}
#else // defined(GIZMO_UNPHYSICAL_ERROR) || defined(GIZMO_UNPHYSICAL_RESCUE)
#define gizmo_check_physical_quantity(name, quantity)
#define gizmo_check_physical_quantities( \
mass_name, energy_name, mass, momentum_x, momentum_y, momentum_z, energy)
#endif // defined(GIZMO_UNPHYSICAL_ERROR) || defined(GIZMO_UNPHYSICAL_RESCUE)
......
......@@ -94,10 +94,13 @@ __attribute__((always_inline)) INLINE static void hydro_velocities_set(
#else // GIZMO_FIX_PARTICLES
if (p->conserved.mass > 0.f && p->primitives.rho > 0.f) {
const float inverse_mass = 1.f / p->conserved.mass;
/* Normal case: set particle velocity to fluid velocity. */
p->v[0] = p->conserved.momentum[0] / p->conserved.mass;
p->v[1] = p->conserved.momentum[1] / p->conserved.mass;
p->v[2] = p->conserved.momentum[2] / p->conserved.mass;
p->v[0] = p->conserved.momentum[0] * inverse_mass;
p->v[1] = p->conserved.momentum[1] * inverse_mass;
p->v[2] = p->conserved.momentum[2] * inverse_mass;
#ifdef GIZMO_STEER_MOTION
......@@ -118,8 +121,7 @@ __attribute__((always_inline)) INLINE static void hydro_velocities_set(
const float soundspeed =
sqrtf(hydro_gamma * p->primitives.P / p->primitives.rho);
/* We only apply the correction if the offset between centroid and position
is
too large. */
is too large. */
if (d > 0.9f * etaR) {
float fac = xi * soundspeed / d;
if (d < 1.1f * etaR) {
......
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2018 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_RIEMANN_CHECKS_H
#define SWIFT_RIEMANN_CHECKS_H
#include "error.h"
#ifdef SWIFT_DEBUG_CHECKS
/**
* @brief Check if the given state vector is physically valid.
*
* A valid state vector has 5 finite elements (not NaN), and a positive density
* and pressure (first and fifth element).
*
* @param W State vector.
* @return 0 if the state vector is valid, 1 otherwise.
*/
__attribute__((always_inline)) INLINE static int riemann_check_state(
const float *W) {
int errorFlag = 0;
/* check the density: should be finite and positive */
if (W[0] != W[0]) {
message("NaN density!");
errorFlag = 1;
}
if (W[0] < 0.f) {
message("Negative density!");
errorFlag = 1;
}
/* check the velocities: should be finite */
if (W[1] != W[1]) {
message("NaN x velocity!");
errorFlag = 1;
}
if (W[2] != W[2]) {
message("NaN y velocity!");
errorFlag = 1;
}
if (W[3] != W[3]) {
message("NaN z velocity!");
errorFlag = 1;
}
/* check the pressure: should be positive and finite */
if (W[4] != W[4]) {
message("NaN pressure!");
errorFlag = 1;
}
if (W[4] < 0.f) {
message("Negative pressure!");
errorFlag = 1;
}
return errorFlag;
}
/**
* @brief Check that the given vector is physically valid.
*
* A valid vector has 3 finite elements (not NaN).
*
* @param x Vector to check.
* @return 0 if the vector is valid, 1 otherwise.
*/
__attribute__((always_inline)) INLINE static int riemann_check_vector(
const float *x) {
int errorFlag = 0;
/* check that all components are finite */
if (x[0] != x[0]) {
message("NaN x component!");
errorFlag = 1;
}
if (x[1] != x[1]) {
message("NaN y component!");
errorFlag = 1;
}
if (x[2] != x[2]) {
message("NaN z component!");
errorFlag = 1;
}
return errorFlag;
}
/**
* @brief Check that the given input for a Riemann solver is physically valid.
*
* Physically valid input consists of a valid left and right state, and valid
* surface normal and surface velocity vectors.
* If no valid input is provided, an error is thrown.
*
* @param WL Left state vector.
* @param WR Right state vector.
* @param n Surface normal vector.
* @param vij Surface velocity vector.
*/
__attribute__((always_inline)) INLINE static void riemann_check_input(
const float *WL, const float *WR, const float *n, const float *vij) {
int errorFlag = 0;
if (riemann_check_state(WL)) {
message("Invalid left state!");
errorFlag = 1;
}
if (riemann_check_state(WR)) {
message("Invalid right state!");
errorFlag = 1;
}
if (riemann_check_vector(n)) {
message("Invalid surface normal vector!");
errorFlag = 1;
}
if (riemann_check_vector(vij)) {
message("Invalid face velocity vector!");
errorFlag = 1;
}
if (errorFlag) {
message("WL: %g %g %g %g %g", WL[0], WL[1], WL[2], WL[3], WL[4]);
message("WR: %g %g %g %g %g", WR[0], WR[1], WR[2], WR[3], WR[4]);
message("n: %g %g %g", n[0], n[1], n[2]);
message("vij: %g %g %g", vij[0], vij[1], vij[2]);
error("Invalid Riemann solver input!");
}
}
/**
* @brief Check that the given output from a Riemann solver is physically valid.
*
* Physically valid output consists of 5 finite (not NaN) flux components.
* If no valid output is provided, an error is thrown.
*
* @param WL Left state vector.
* @param WR Right state vector.
* @param n Surface normal vector.
* @param vij Surface velocity vector.
* @param totflux Riemann solver flux result.
*/
__attribute__((always_inline)) INLINE static void riemann_check_output(
const float *WL, const float *WR, const float *n, const float *vij,
const float *totflux) {
int errorFlag = 0;
/* check that all the fluxes are finite */
if (totflux[0] != totflux[0]) {
message("NaN mass flux!");
errorFlag = 1;
}
if (totflux[1] != totflux[1]) {
message("NaN x momentum flux!");
errorFlag = 1;
}
if (totflux[2] != totflux[2]) {
message("NaN y momentum flux!");
errorFlag = 1;
}
if (totflux[3] != totflux[3]) {
message("NaN z momentum flux!");
errorFlag = 1;
}
if (totflux[4] != totflux[4]) {
message("NaN energy flux!");
errorFlag = 1;
}
if (errorFlag) {
message("WL: %g %g %g %g %g", WL[0], WL[1], WL[2], WL[3], WL[4]);
message("WR: %g %g %g %g %g", WR[0], WR[1], WR[2], WR[3], WR[4]);
message("n: %g %g %g", n[0], n[1], n[2]);
message("vij: %g %g %g", vij[0], vij[1], vij[2]);
message("totflux: %g %g %g %g %g", totflux[0], totflux[1], totflux[2],
totflux[3], totflux[4]);
error("Invalid Riemann solver output!");
}
}
#endif
#endif /* SWIFT_RIEMANN_CHECKS_H */
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
* Copyright (c) 2015 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
* 2018 Bert Vandenbroucke (bert.vandenbroucke@gmail.com)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
......@@ -38,6 +39,7 @@
#include "adiabatic_index.h"
#include "error.h"
#include "minmax.h"
#include "riemann_checks.h"
#include "riemann_vacuum.h"
/**
......@@ -318,30 +320,6 @@ __attribute__((always_inline)) INLINE static void riemann_solver_solve(
float SHR, STR;
float pdpL, SL;
float SHL, STL;
int errorFlag = 0;
/* sanity checks */
if (WL[0] != WL[0] || WL[4] != WL[4]) {
printf("NaN WL!\n");
errorFlag = 1;
}
if (WR[0] != WR[0] || WR[4] != WR[4]) {
printf("NaN WR!\n");
errorFlag = 1;
}
if (WL[0] < 0.0f || WL[4] < 0.0f) {
printf("Negative WL!\n");
errorFlag = 1;
}
if (WR[0] < 0.0f || WR[4] < 0.0f) {
printf("Negative WR!\n");
errorFlag = 1;
}
if (errorFlag) {
printf("WL: %g %g %g %g %g\n", WL[0], WL[1], WL[2], WL[3], WL[4]);
printf("WR: %g %g %g %g %g\n", WR[0], WR[1], WR[2], WR[3], WR[4]);
error("Riemman solver input error!\n");
}
/* calculate velocities in interface frame */
vL = WL[1] * n_unit[0] + WL[2] * n_unit[1] + WL[3] * n_unit[2];
......@@ -511,6 +489,10 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
const float* Wi, const float* Wj, const float* n_unit, const float* vij,
float* totflux) {
#ifdef SWIFT_DEBUG_CHECKS
riemann_check_input(Wi, Wj, n_unit, vij);
#endif
float Whalf[5];
float flux[5][3];
float vtot[3];
......@@ -555,6 +537,10 @@ __attribute__((always_inline)) INLINE static void riemann_solve_for_flux(
flux[3][0] * n_unit[0] + flux[3][1] * n_unit[1] + flux[3][2] * n_unit[2];
totflux[4] =
flux[4][0] * n_unit[0] + flux[4][1] * n_unit[1] + flux[4][2] * n_unit[2];
#ifdef SWIFT_DEBUG_CHECKS
riemann_check_output(Wi, Wj, n_unit, vij, totflux);
#endif
}
#endif /* SWIFT_RIEMANN_EXACT_H */
/*******************************************************************************
* This file is part of SWIFT.
* Coypright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@ugent.be)
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU Lesser General Public License as published
* by the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*
******************************************************************************/
#ifndef SWIFT_RIEMANN_EXACT_ISOTHERMAL_H
#define SWIFT_RIEMANN_EXACT_ISOTHERMAL_H
#include <float.h>
#include "adiabatic_index.h"
#include "minmax.h"
#include "riemann_vacuum.h"
#define const_isothermal_soundspeed \
sqrtf(hydro_gamma_minus_one* const_isothermal_internal_energy)
/**
* @brief Relative difference between the middle state velocity and the left or
* right state velocity used in the middle state density iteration.
*
* @param rho Current estimate of the middle state density.
* @param W Left or right state vector.
* @return Density dependent part of the middle state velocity.
*/
__attribute__((always_inline)) INLINE static float riemann_fb(float rho,
float* W) {
if (rho < W[0]) {
return const_isothermal_soundspeed * logf(rho / W[0]);
} else {
return const_isothermal_soundspeed *
(sqrtf(rho / W[0]) - sqrtf(W[0] / rho));
}
}
/**
* @brief Derivative w.r.t. rho of the function riemann_fb.
*
* @param rho Current estimate of the middle state density.
* @param W Left or right state vector.
* @return Derivative of riemann_fb.
*/
__attribute__((always_inline)) INLINE static float riemann_fprimeb(float rho,
float* W) {
if (rho < W[0]) {
return const_isothermal_soundspeed * W[0] / rho;
} else {
return 0.5f * const_isothermal_soundspeed *
(sqrtf(rho / W[0]) + sqrtf(W[0] / rho)) / rho;
}
}
/**
* @brief Difference between the left and right middle state velocity estimates.
*
* Since the middle state velocity takes on a constant value, we want to get
* this difference as close to zero as possible.
*
* @param rho Current estimate of the middle state density.
* @param WL Left state vector.
* @param WR Right state vector.
* @param vL Left state velocity along the interface normal.
* @param vR Right state velocity along the interface normal.
* @return Difference between the left and right middle state velocity
* estimates.
*/
__attribute__((always_inline)) INLINE static float riemann_f(
float rho, float* WL, float* WR, float vL, float vR) {
return riemann_fb(rho, WR) + riemann_fb(rho, WL) + vR - vL;
}
/**
* @brief Derivative of riemann_f w.r.t. rho.
*
* @param rho Current estimate of the middle state density.
* @param WL Left state vector.
* @param WR Right state vector.
* @return Derivative of riemann_f.
*/
__attribute__((always_inline)) INLINE static float riemann_fprime(float rho,
float* WL,
float* WR) {
return riemann_fprimeb(rho, WL) + riemann_fprimeb(rho, WR);
}
/**
* @brief Get a good first guess for the middle state density.
*
* @param WL The left state vector
* @param WR The right state vector
* @param vL The left velocity along the interface normal
* @param vR The right velocity along the interface normal
*/
__attribute__((always_inline)) INLINE static float riemann_guess_rho(float* WL,
float* WR,
float vL,
float vR) {
/* Currently three possibilities and not really an algorithm to decide which
one to choose: */
/* just the average */
// return 0.5f * (WL[0] + WR[0]);
/* two rarefaction approximation */
return sqrtf(WL[0] * WR[0] * expf((vL - vR) / const_isothermal_soundspeed));
/* linearized primitive variable approximation */
return 0.25f * (WL[0] + WR[0]) * (vL - vR) / const_isothermal_soundspeed +
0.5f * (WL[0] + WR[0]);
}
/**
* @brief Find the zeropoint of riemann_f(rho) using Brent's method.
*
* @param lower_limit Lower limit for the method (riemann_f(lower_limit) < 0)
* @param upper_limit Upper limit for the method (riemann_f(upper_limit) > 0)
* @param lowf Value of riemann_f(lower_limit).
* @param upf Value of riemann_f(upper_limit).
* @param error_tol Tolerance used to decide if the solution is converged
* @param WL Left state vector
* @param WR Right state vector
* @param vL The left velocity along the interface normal
* @param vR The right velocity along the interface normal
*/
__attribute__((always_inline)) INLINE static float riemann_solve_brent(
float lower_limit, float upper_limit, float lowf, float upf,
float error_tol, float* WL, float* WR, float vL, float vR) {
float a, b, c, d, s;
float fa, fb, fc, fs;
float tmp, tmp2;
int mflag;
int i;
a = lower_limit;
b = upper_limit;
c = 0.0f;
d = FLT_MAX;
fa = lowf;
fb = upf;
fc = 0.0f;
s = 0.0f;
fs = 0.0f;
/* if f(a) f(b) >= 0 then error-exit */
if (fa * fb >= 0.0f) {
error(
"Brent's method called with equal sign function values!\n"
"f(%g) = %g, f(%g) = %g\n",
a, fa, b, fb);
/* return NaN */
return 0.0f / 0.0f;
}
/* if |f(a)| < |f(b)| then swap (a,b) */
if (fabs(fa) < fabs(fb)) {
tmp = a;
a = b;
b = tmp;
tmp = fa;
fa = fb;
fb = tmp;
}
c = a;
fc = fa;
mflag = 1;
i = 0;
while (!(fb == 0.0f) && (fabs(a - b) > error_tol * 0.5f * (a + b))) {
if ((fa != fc) && (fb != fc)) /* Inverse quadratic interpolation */
s = a * fb * fc / (fa - fb) / (fa - fc) +
b * fa * fc / (fb - fa) / (fb - fc) +
c * fa * fb / (fc - fa) / (fc - fb);
else
/* Secant Rule */
s = b - fb * (b - a) / (fb - fa);
tmp2 = 0.25f * (3.0f * a + b);
if (!(((s > tmp2) && (s < b)) || ((s < tmp2) && (s > b))) ||