1 // Copyright (c) 2017-2022, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 /// @file 9 /// Helper functions for solving the Riemann problem. 10 // The left and right states are specified from the perspective of an outward-facing normal vector pointing left to right: 11 // 12 // (domain) 13 // / (outward facing normal) 14 // |------------| / 15 // | | / 16 // | Left |----> Right 17 // | (Interior) | (Exterior) 18 // |------------| 19 // 20 // The right state is exterior to the domain and the left state is the interior to the domain. 21 // Much of the work references Eleuterio F. Toro's "Riemann Solvers and Numerical Methods for Fluid Dynamics", 2009 22 23 #ifndef riemann_solver_h 24 #define riemann_solver_h 25 26 #include "newtonian_state.h" 27 #include "newtonian_types.h" 28 29 enum RiemannFluxType_ { RIEMANN_HLL, RIEMANN_HLLC }; 30 typedef enum RiemannFluxType_ RiemannFluxType; 31 32 typedef struct { 33 CeedScalar left, right; 34 } RoeWeights; 35 36 CEED_QFUNCTION_HELPER RoeWeights RoeSetup(CeedScalar rho_left, CeedScalar rho_right) { 37 CeedScalar sqrt_left = sqrt(rho_left), sqrt_right = sqrt(rho_right); 38 RoeWeights w = {sqrt_left / (sqrt_left + sqrt_right), sqrt_right / (sqrt_left + sqrt_right)}; 39 return w; 40 } 41 42 CEED_QFUNCTION_HELPER RoeWeights RoeSetup_fwd(CeedScalar rho_left, CeedScalar rho_right, CeedScalar drho_left, CeedScalar drho_right) { 43 CeedScalar sqrt_left = sqrt(rho_left), sqrt_right = sqrt(rho_right); 44 CeedScalar square_sum_root = Square(sqrt_left + sqrt_right); 45 CeedScalar r_right = (sqrt_left / (2 * sqrt_right * square_sum_root)) * drho_right - (sqrt_right / (2 * sqrt_left * square_sum_root)) * drho_left; 46 CeedScalar r_left = (sqrt_right / (2 * sqrt_left * square_sum_root)) * drho_left - (sqrt_left / (2 * sqrt_right * square_sum_root)) * drho_right; 47 RoeWeights dw = {r_left, r_right}; 48 return dw; 49 } 50 51 CEED_QFUNCTION_HELPER CeedScalar RoeAverage(RoeWeights r, CeedScalar q_left, CeedScalar q_right) { return r.left * q_left + r.right * q_right; } 52 53 CEED_QFUNCTION_HELPER CeedScalar RoeAverage_fwd(RoeWeights r, RoeWeights dr, CeedScalar q_left, CeedScalar q_right, CeedScalar dq_left, 54 CeedScalar dq_right) { 55 return q_right * dr.right + q_left * dr.left + r.right * dq_right + r.left * dq_left; 56 } 57 58 CEED_QFUNCTION_HELPER StateConservative Flux_HLL(State left, State right, StateConservative flux_left, StateConservative flux_right, 59 CeedScalar s_left, CeedScalar s_right) { 60 CeedScalar U_left[5], U_right[5], F_right[5], F_left[5], F_hll[5]; 61 UnpackState_U(left.U, U_left); 62 UnpackState_U(right.U, U_right); 63 UnpackState_U(flux_left, F_left); 64 UnpackState_U(flux_right, F_right); 65 for (int i = 0; i < 5; i++) { 66 F_hll[i] = (s_right * F_left[i] - s_left * F_right[i] + s_left * s_right * (U_right[i] - U_left[i])) / (s_right - s_left); 67 } 68 StateConservative F = { 69 F_hll[0], 70 {F_hll[1], F_hll[2], F_hll[3]}, 71 F_hll[4], 72 }; 73 return F; 74 } 75 76 CEED_QFUNCTION_HELPER StateConservative Flux_HLL_fwd(State left, State right, State dleft, State dright, StateConservative flux_left, 77 StateConservative flux_right, StateConservative dflux_left, StateConservative dflux_right, 78 CeedScalar S_l, CeedScalar S_r, CeedScalar dS_l, CeedScalar dS_r) { 79 CeedScalar U_l[5], U_r[5], F_r[5], F_l[5]; 80 UnpackState_U(left.U, U_l); 81 UnpackState_U(right.U, U_r); 82 UnpackState_U(flux_left, F_l); 83 UnpackState_U(flux_right, F_r); 84 85 CeedScalar dU_l[5], dU_r[5], dF_r[5], dF_l[5], dF_hll[5] = {0.}; 86 UnpackState_U(dleft.U, dU_l); 87 UnpackState_U(dright.U, dU_r); 88 UnpackState_U(dflux_left, dF_l); 89 UnpackState_U(dflux_right, dF_r); 90 for (int i = 0; i < 5; i++) { 91 const CeedScalar U_diff = U_r[i] - U_l[i]; 92 const CeedScalar S_diff = S_r - S_l; 93 const CeedScalar F_hll_denom = S_r * F_l[i] - S_l * F_r[i] + S_l * S_r * U_diff; 94 95 dF_hll[i] += ((F_l[i] + S_r * U_diff) * S_diff - F_hll_denom) / Square(S_diff) * dS_r; 96 dF_hll[i] += ((-F_r[i] + S_r * U_diff) * S_diff + F_hll_denom) / Square(S_diff) * dS_l; 97 dF_hll[i] += (S_r * dF_l[i] - S_l * dF_r[i] + S_r * S_l * dU_r[i] - S_r * S_l * dU_l[i]) / S_diff; 98 } 99 StateConservative dF = { 100 dF_hll[0], 101 {dF_hll[1], dF_hll[2], dF_hll[3]}, 102 dF_hll[4], 103 }; 104 return dF; 105 } 106 107 CEED_QFUNCTION_HELPER void ComputeHLLSpeeds_Roe(NewtonianIdealGasContext gas, State left, CeedScalar u_left, State right, CeedScalar u_right, 108 CeedScalar *s_left, CeedScalar *s_right) { 109 const CeedScalar gamma = HeatCapacityRatio(gas); 110 111 RoeWeights r = RoeSetup(left.U.density, right.U.density); 112 // Speed estimate 113 // Roe average eigenvalues for left and right non-linear waves. 114 // Stability requires that these speed estimates are *at least* as fast as the physical wave speeds. 115 CeedScalar u_roe = RoeAverage(r, u_left, u_right); 116 117 // TODO: revisit this for gravity 118 CeedScalar H_left = TotalSpecificEnthalpy(gas, left); 119 CeedScalar H_right = TotalSpecificEnthalpy(gas, right); 120 CeedScalar H_roe = RoeAverage(r, H_left, H_right); 121 CeedScalar a_roe = sqrt((gamma - 1) * (H_roe - 0.5 * Square(u_roe))); 122 123 // Einfeldt (1988) justifies (and Toro's book repeats) that Roe speeds can be used here. 124 *s_left = u_roe - a_roe; 125 *s_right = u_roe + a_roe; 126 } 127 128 CEED_QFUNCTION_HELPER void ComputeHLLSpeeds_Roe_fwd(NewtonianIdealGasContext gas, State left, State dleft, CeedScalar u_left, CeedScalar du_left, 129 State right, State dright, CeedScalar u_right, CeedScalar du_right, CeedScalar *s_left, 130 CeedScalar *ds_left, CeedScalar *s_right, CeedScalar *ds_right) { 131 const CeedScalar gamma = HeatCapacityRatio(gas); 132 133 RoeWeights r = RoeSetup(left.U.density, right.U.density); 134 RoeWeights dr = RoeSetup_fwd(left.U.density, right.U.density, dleft.U.density, dright.U.density); 135 // Speed estimate 136 // Roe average eigenvalues for left and right non-linear waves. 137 // Stability requires that these speed estimates are *at least* as fast as the physical wave speeds. 138 CeedScalar u_roe = RoeAverage(r, u_left, u_right); 139 CeedScalar du_roe = RoeAverage_fwd(r, dr, u_left, u_right, du_left, du_right); 140 141 CeedScalar H_left = TotalSpecificEnthalpy(gas, left); 142 CeedScalar H_right = TotalSpecificEnthalpy(gas, right); 143 CeedScalar dH_left = TotalSpecificEnthalpy_fwd(gas, left, dleft); 144 CeedScalar dH_right = TotalSpecificEnthalpy_fwd(gas, right, dright); 145 146 CeedScalar H_roe = RoeAverage(r, H_left, H_right); 147 CeedScalar dH_roe = RoeAverage_fwd(r, dr, H_left, H_right, dH_left, dH_right); 148 CeedScalar a_roe = sqrt((gamma - 1) * (H_roe - 0.5 * Square(u_roe))); 149 CeedScalar da_roe = 0.5 * (gamma - 1) / sqrt(H_roe) * dH_roe - 0.5 * sqrt(gamma - 1) * u_roe / sqrt(H_roe - 0.5 * Square(u_roe)) * du_roe; 150 151 *s_left = u_roe - a_roe; 152 *ds_left = du_roe - da_roe; 153 *s_right = u_roe + a_roe; 154 *ds_right = du_roe + da_roe; 155 } 156 157 // ***************************************************************************** 158 // @brief Harten Lax VanLeer (HLL) approximate Riemann solver. 159 // Taking in two states (left, right) and returns RiemannFlux_HLL. 160 // The left and right states are specified from the perspective of an outward-facing normal vector pointing left to right. 161 // 162 // @param[in] gas NewtonianIdealGasContext for the fluid 163 // @param[in] left Fluid state of the domain interior (the current solution) 164 // @param[in] right Fluid state of the domain exterior (free stream conditions) 165 // @param[in] normal Normalized, outward facing boundary normal vector 166 // 167 // @return StateConservative with HLL Riemann Flux 168 // ***************************************************************************** 169 CEED_QFUNCTION_HELPER StateConservative RiemannFlux_HLL(NewtonianIdealGasContext gas, State left, State right, const CeedScalar normal[3]) { 170 StateConservative flux_left = FluxInviscidDotNormal(gas, left, normal); 171 StateConservative flux_right = FluxInviscidDotNormal(gas, right, normal); 172 173 CeedScalar u_left = Dot3(left.Y.velocity, normal); 174 CeedScalar u_right = Dot3(right.Y.velocity, normal); 175 176 CeedScalar s_left, s_right; 177 ComputeHLLSpeeds_Roe(gas, left, u_left, right, u_right, &s_left, &s_right); 178 179 // Compute HLL flux 180 if (0 <= s_left) { 181 return flux_left; 182 } else if (s_right <= 0) { 183 return flux_right; 184 } else { 185 return Flux_HLL(left, right, flux_left, flux_right, s_left, s_right); 186 } 187 } 188 189 // ***************************************************************************** 190 // @brief Forward-mode Derivative of Harten Lax VanLeer (HLL) approximate Riemann solver. 191 // 192 // @param gas NewtonianIdealGasContext for the fluid 193 // @param left Fluid state of the domain interior (the current solution) 194 // @param right Fluid state of the domain exterior (free stream conditions) 195 // @param dleft Derivative of fluid state of the domain interior (the current solution) 196 // @param dright Derivative of fluid state of the domain exterior (free stream conditions) 197 // @param normal Normalized, outward facing boundary normal vector 198 // 199 // @return StateConservative with derivative of HLL Riemann Flux 200 // ***************************************************************************** 201 CEED_QFUNCTION_HELPER StateConservative RiemannFlux_HLL_fwd(NewtonianIdealGasContext gas, State left, State dleft, State right, State dright, 202 const CeedScalar normal[3]) { 203 StateConservative flux_left = FluxInviscidDotNormal(gas, left, normal); 204 StateConservative flux_right = FluxInviscidDotNormal(gas, right, normal); 205 StateConservative dflux_left = FluxInviscidDotNormal_fwd(gas, left, dleft, normal); 206 StateConservative dflux_right = FluxInviscidDotNormal_fwd(gas, right, dright, normal); 207 208 CeedScalar u_left = Dot3(left.Y.velocity, normal); 209 CeedScalar u_right = Dot3(right.Y.velocity, normal); 210 CeedScalar du_left = Dot3(dleft.Y.velocity, normal); 211 CeedScalar du_right = Dot3(dright.Y.velocity, normal); 212 213 CeedScalar s_left, ds_left, s_right, ds_right; 214 ComputeHLLSpeeds_Roe_fwd(gas, left, dleft, u_left, du_left, right, dright, u_right, du_right, &s_left, &ds_left, &s_right, &ds_right); 215 216 if (0 <= s_left) { 217 return dflux_left; 218 } else if (s_right <= 0) { 219 return dflux_right; 220 } else { 221 return Flux_HLL_fwd(left, right, dleft, dright, flux_left, flux_right, dflux_left, dflux_right, s_left, s_right, ds_left, ds_right); 222 } 223 } 224 225 CEED_QFUNCTION_HELPER StateConservative RiemannFlux_HLLC_Star(NewtonianIdealGasContext gas, State side, StateConservative F_side, 226 const CeedScalar normal[3], CeedScalar u_side, CeedScalar s_side, CeedScalar s_star) { 227 CeedScalar fact = side.U.density * (s_side - u_side) / (s_side - s_star); 228 CeedScalar denom = side.U.density * (s_side - u_side); 229 // U_* = fact * star 230 StateConservative star = { 231 1.0, 232 { 233 side.Y.velocity[0] + (s_star - u_side) * normal[0], 234 side.Y.velocity[1] + (s_star - u_side) * normal[1], 235 side.Y.velocity[2] + (s_star - u_side) * normal[2], 236 }, 237 side.U.E_total / side.U.density // 238 + (s_star - u_side) * (s_star + side.Y.pressure / denom) 239 }; 240 return StateConservativeAXPBYPCZ(1, F_side, s_side * fact, star, -s_side, side.U); 241 } 242 243 CEED_QFUNCTION_HELPER StateConservative RiemannFlux_HLLC_Star_fwd(NewtonianIdealGasContext gas, State side, State dside, StateConservative F_side, 244 StateConservative dF_side, const CeedScalar normal[3], CeedScalar u_side, 245 CeedScalar du_side, CeedScalar s_side, CeedScalar ds_side, CeedScalar s_star, 246 CeedScalar ds_star) { 247 CeedScalar fact = side.U.density * (s_side - u_side) / (s_side - s_star); 248 CeedScalar dfact = (side.U.density * (ds_side - du_side) + dside.U.density * (s_side - u_side)) / (s_side - s_star) // 249 - fact / (s_side - s_star) * (ds_side - ds_star); 250 CeedScalar denom = side.U.density * (s_side - u_side); 251 CeedScalar ddenom = side.U.density * (ds_side - du_side) + dside.U.density * (s_side - u_side); 252 253 StateConservative star = { 254 1.0, 255 { 256 side.Y.velocity[0] + (s_star - u_side) * normal[0], 257 side.Y.velocity[1] + (s_star - u_side) * normal[1], 258 side.Y.velocity[2] + (s_star - u_side) * normal[2], 259 }, 260 side.U.E_total / side.U.density // 261 + (s_star - u_side) * (s_star + side.Y.pressure / denom) 262 }; 263 StateConservative dstar = { 264 0., 265 { 266 dside.Y.velocity[0] + (ds_star - du_side) * normal[0], 267 dside.Y.velocity[1] + (ds_star - du_side) * normal[1], 268 dside.Y.velocity[2] + (ds_star - du_side) * normal[2], 269 }, 270 dside.U.E_total / side.U.density - side.U.E_total / Square(side.U.density) * dside.U.density // 271 + (ds_star - du_side) * (s_star + side.Y.pressure / denom) // 272 + (s_star - u_side) * (ds_star + dside.Y.pressure / denom - side.Y.pressure / Square(denom) * ddenom) // 273 }; 274 275 const CeedScalar a[] = {1, ds_side * fact + s_side * dfact, s_side * fact, -ds_side, -s_side}; 276 const StateConservative U[] = {dF_side, star, dstar, side.U, dside.U}; 277 return StateConservativeMult(5, a, U); 278 } 279 280 // HLLC Riemann solver (from Toro's book) 281 CEED_QFUNCTION_HELPER StateConservative RiemannFlux_HLLC(NewtonianIdealGasContext gas, State left, State right, const CeedScalar normal[3]) { 282 StateConservative flux_left = FluxInviscidDotNormal(gas, left, normal); 283 StateConservative flux_right = FluxInviscidDotNormal(gas, right, normal); 284 285 CeedScalar u_left = Dot3(left.Y.velocity, normal); 286 CeedScalar u_right = Dot3(right.Y.velocity, normal); 287 CeedScalar s_left, s_right; 288 ComputeHLLSpeeds_Roe(gas, left, u_left, right, u_right, &s_left, &s_right); 289 290 // Contact wave speed; Toro (10.37) 291 CeedScalar rhou_left = left.U.density * u_left, rhou_right = right.U.density * u_right; 292 CeedScalar numer = right.Y.pressure - left.Y.pressure + rhou_left * (s_left - u_left) - rhou_right * (s_right - u_right); 293 CeedScalar denom = left.U.density * (s_left - u_left) - right.U.density * (s_right - u_right); 294 CeedScalar s_star = numer / denom; 295 296 // Compute HLLC flux 297 if (0 <= s_left) { 298 return flux_left; 299 } else if (0 <= s_star) { 300 return RiemannFlux_HLLC_Star(gas, left, flux_left, normal, u_left, s_left, s_star); 301 } else if (0 <= s_right) { 302 return RiemannFlux_HLLC_Star(gas, right, flux_right, normal, u_right, s_right, s_star); 303 } else { 304 return flux_right; 305 } 306 } 307 308 CEED_QFUNCTION_HELPER StateConservative RiemannFlux_HLLC_fwd(NewtonianIdealGasContext gas, State left, State dleft, State right, State dright, 309 const CeedScalar normal[3]) { 310 StateConservative flux_left = FluxInviscidDotNormal(gas, left, normal); 311 StateConservative flux_right = FluxInviscidDotNormal(gas, right, normal); 312 StateConservative dflux_left = FluxInviscidDotNormal_fwd(gas, left, dleft, normal); 313 StateConservative dflux_right = FluxInviscidDotNormal_fwd(gas, right, dright, normal); 314 315 CeedScalar u_left = Dot3(left.Y.velocity, normal); 316 CeedScalar u_right = Dot3(right.Y.velocity, normal); 317 CeedScalar du_left = Dot3(dleft.Y.velocity, normal); 318 CeedScalar du_right = Dot3(dright.Y.velocity, normal); 319 320 CeedScalar s_left, ds_left, s_right, ds_right; 321 ComputeHLLSpeeds_Roe_fwd(gas, left, dleft, u_left, du_left, right, dright, u_right, du_right, &s_left, &ds_left, &s_right, &ds_right); 322 323 // Contact wave speed; Toro (10.37) 324 CeedScalar rhou_left = left.U.density * u_left, drhou_left = left.U.density * du_left + dleft.U.density * u_left; 325 CeedScalar rhou_right = right.U.density * u_right, drhou_right = right.U.density * du_right + dright.U.density * u_right; 326 CeedScalar numer = right.Y.pressure - left.Y.pressure // 327 + rhou_left * (s_left - u_left) // 328 - rhou_right * (s_right - u_right); 329 CeedScalar dnumer = dright.Y.pressure - dleft.Y.pressure // 330 + rhou_left * (ds_left - du_left) + drhou_left * (s_left - u_left) // 331 - rhou_right * (ds_right - du_right) - drhou_right * (s_right - u_right); 332 CeedScalar denom = left.U.density * (s_left - u_left) - right.U.density * (s_right - u_right); 333 CeedScalar ddenom = left.U.density * (ds_left - du_left) + dleft.U.density * (s_left - u_left) // 334 - right.U.density * (ds_right - du_right) - dright.U.density * (s_right - u_right); 335 CeedScalar s_star = numer / denom; 336 CeedScalar ds_star = dnumer / denom - numer / Square(denom) * ddenom; 337 338 // Compute HLLC flux 339 if (0 <= s_left) { 340 return dflux_left; 341 } else if (0 <= s_star) { 342 return RiemannFlux_HLLC_Star_fwd(gas, left, dleft, flux_left, dflux_left, normal, u_left, du_left, s_left, ds_left, s_star, ds_star); 343 } else if (0 <= s_right) { 344 return RiemannFlux_HLLC_Star_fwd(gas, right, dright, flux_right, dflux_right, normal, u_right, du_right, s_right, ds_right, s_star, ds_star); 345 } else { 346 return dflux_right; 347 } 348 } 349 350 #endif // riemann_solver_h 351