/******************************************************************************* * This file is part of SWIFT. * Copyright (c) 2016 Bert Vandenbroucke (bert.vandenbroucke@gmail.com) * Copyright (c) 2021 Mladen Ivkovic (mladen.ivkovic@hotmail.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 * 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 . * ******************************************************************************/ #ifndef SWIFT_RT_SLOPE_LIMITERS_FACE_GEAR_H #define SWIFT_RT_SLOPE_LIMITERS_FACE_GEAR_H #include "sign.h" /** * @file src/rt/GEAR/rt_slope_limiters_face.h * @brief File containing routines concerning the face slope * limiter for the GEAR RT scheme, i.e. limiting the actual * interparticle flux. */ /** * @brief Slope limit a single quantity at the interface. CURRENTLY NOT IN USE. * * @param phi_i Value of the quantity at the particle position. * @param phi_j Value of the quantity at the neighbouring particle position. * @param phi_mid0 Extrapolated value of the quantity at the interface position. * @param xij_norm Distance between the particle position and the interface * position. * @param r_inv Inverse distance between the particle and its neighbour. * @return The slope limited difference between the quantity at the particle * position and the quantity at the interface position. */ __attribute__((always_inline)) INLINE static float rt_slope_limit_face_quantity( float phi_i, float phi_j, float phi_mid0, float xij_norm, float r_inv) { const float dphi = phi_i - phi_j; const float gamma1 = 0.5f; const float gamma2 = 0.25f; const float delta1 = gamma1 * fabsf(dphi); const float delta2 = gamma2 * fabsf(dphi); const float phimin = min(phi_i, phi_j); const float phimax = max(phi_i, phi_j); const float phibar = phi_i + xij_norm * r_inv * dphi; float phiplus, phiminus, phi_mid; if (same_signf(phimax + delta1, phimax)) { phiplus = phimax + delta1; } else { phiplus = phimax / (1.0f + delta1 / (fabsf(phimax) + FLT_MIN)); } if (same_signf(phimin - delta1, phimin)) { phiminus = phimin - delta1; } else { phiminus = phimin / (1.0f + delta1 / (fabsf(phimin) + FLT_MIN)); } if (phi_i < phi_j) { const float temp = min(phibar + delta2, phi_mid0); phi_mid = max(phiminus, temp); } else { const float temp = max(phibar - delta2, phi_mid0); phi_mid = min(phiplus, temp); } return phi_mid - phi_i; } /** * the minmod slope limiter. * * @param dQi left slope * @param dQj right slope */ __attribute__((always_inline)) INLINE static void rt_limiter_minmod( float* dQi, float* dQj) { if (*dQi * *dQj > 0.f) { if (fabsf(*dQi) < fabsf(*dQj)) { *dQj = *dQi; } else { *dQi = *dQj; } } else { *dQi = 0.f; *dQj = 0.f; } } /** * the monotonized cenetral limiter. * * @param dQi left slope * @param dQj right slope * @return factor to slope limit the slope dQi */ __attribute__((always_inline)) INLINE static float rt_limiter_mc( const float dQi, const float dQj) { const float r = dQj == 0.f ? dQi * 1e6 : dQi / dQj; const float minterm = min3(0.5 * (1. + r), 2.f, 2.f * r); return max(0.f, minterm); } /** * the van Leer limiter. * * @param dQi left slope * @param dQj right slope * @return factor to slope limit the slope dQi */ __attribute__((always_inline)) INLINE static float rt_limiter_vanLeer( const float dQi, const float dQj) { const float r = dQj == 0.f ? dQi * 1e6 : dQi / dQj; const float absr = fabsf(r); return (r + absr) / (1.f + absr); } /** * the superbee limiter. * * @param dQi left slope * @param dQj right slope * @return factor to slope limit the slope dQi */ __attribute__((always_inline)) INLINE static float rt_limiter_superbee( const float dQi, const float dQj) { const float r = dQj == 0.f ? dQi * 1e6 : dQi / dQj; const float minterm1 = min(1.f, 2.f * r); const float minterm2 = min(2.f, r); return max3(0.f, minterm1, minterm2); } /** * @brief Slope limit the slopes at the interface between two particles * * @param Qi RT quantities of particle i (energy density + fluxes in 3 dim) * @param Qj RT quantities of particle j (energy density + fluxes in 3 dim) * @param dQi Difference between the RT quantities of particle i at the * position of particle i and at the interface position. * @param dQj Difference between the RT quantities of particle j at the * position of particle j and at the interface position. * @param dx Comoving distance vector between the particles (dx = pi->x - * pj->x). * @param r Comoving distance between particle i and particle j. * @param xij_i Position of the "interface" w.r.t. position of particle i * @param xij_j Position of the "interface" w.r.t. position of particle j */ __attribute__((always_inline)) INLINE static void rt_slope_limit_face( const float Qi[4], const float Qj[4], float dQi[4], float dQj[4], const float* dx, const float r, const float xij_i[3], const float xij_j[3]) { /* In 1D advection tests, any limiter works with the * GLF solver. If you also activate (uncomment) the * cell slope limiting, you should avoid the superbee * and MC limiters. The Gizmo-style clever limiter * that is used for hydro leads to a more diffusive * behaviour. * For the HLL solver, avoid superbee and MC limiters, * as they lead to unstable results. */ /* const float xij_i_norm = */ /* sqrtf(xij_i[0] * xij_i[0] + xij_i[1] * xij_i[1] + xij_i[2] * xij_i[2]); */ /* */ /* const float xij_j_norm = */ /* sqrtf(xij_j[0] * xij_j[0] + xij_j[1] * xij_j[1] + xij_j[2] * xij_j[2]); */ /* */ /* const float r_inv = 1.f / r; */ for (int i = 0; i < 4; i++) { rt_limiter_minmod(&dQi[i], &dQj[i]); /* const float alphai = rt_limiter_mc(dQi[i], dQj[i]); */ /* const float alphaj = rt_limiter_mc(dQj[i], dQi[i]); */ /* dQi[i] *= alphai; */ /* dQj[i] *= alphaj; */ /* const float alphai = rt_limiter_superbee(dQi[i], dQj[i]); */ /* const float alphaj = rt_limiter_superbee(dQj[i], dQi[i]); */ /* dQi[i] *= alphai; */ /* dQj[i] *= alphaj; */ /* const float alphai = rt_limiter_vanLeer(dQi[i], dQj[i]); */ /* const float alphaj = rt_limiter_vanLeer(dQj[i], dQi[i]); */ /* dQi[i] *= alphai; */ /* dQj[i] *= alphaj; */ /* dQi[i] = rt_slope_limit_face_quantity(Qi[i], Qj[i], Qi[i] + dQi[i], */ /* xij_i_norm, r_inv); */ /* */ /* dQj[i] = rt_slope_limit_face_quantity(Qj[i], Qi[i], Qj[i] + dQj[i], */ /* xij_j_norm, r_inv); */ } } #endif /* SWIFT_RT_SLOPE_LIMITERS_FACE_GEAR_H */