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