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