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
RoeSetup(CeedScalar rho_left,CeedScalar rho_right)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
RoeSetup_fwd(CeedScalar rho_left,CeedScalar rho_right,CeedScalar drho_left,CeedScalar drho_right)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
RoeAverage(RoeWeights r,CeedScalar q_left,CeedScalar q_right)43 CEED_QFUNCTION_HELPER CeedScalar RoeAverage(RoeWeights r, CeedScalar q_left, CeedScalar q_right) { return r.left * q_left + r.right * q_right; }
44
RoeAverage_fwd(RoeWeights r,RoeWeights dr,CeedScalar q_left,CeedScalar q_right,CeedScalar dq_left,CeedScalar dq_right)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
Flux_HLL(State left,State right,StateConservative flux_left,StateConservative flux_right,CeedScalar s_left,CeedScalar s_right)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
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)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
ComputeHLLSpeeds_Roe(NewtonianIGProperties gas,State left,CeedScalar u_left,State right,CeedScalar u_right,CeedScalar * s_left,CeedScalar * s_right)97 CEED_QFUNCTION_HELPER void ComputeHLLSpeeds_Roe(NewtonianIGProperties 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
ComputeHLLSpeeds_Roe_fwd(NewtonianIGProperties 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)117 CEED_QFUNCTION_HELPER void ComputeHLLSpeeds_Roe_fwd(NewtonianIGProperties 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 NewtonianIGProperties 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 // *****************************************************************************
RiemannFlux_HLL(NewtonianIGProperties gas,State left,State right,const CeedScalar normal[3])159 CEED_QFUNCTION_HELPER StateConservative RiemannFlux_HLL(NewtonianIGProperties 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 NewtonianIGProperties 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 // *****************************************************************************
RiemannFlux_HLL_fwd(NewtonianIGProperties gas,State left,State dleft,State right,State dright,const CeedScalar normal[3])191 CEED_QFUNCTION_HELPER StateConservative RiemannFlux_HLL_fwd(NewtonianIGProperties 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
RiemannFlux_HLLC_Star(NewtonianIGProperties gas,State side,StateConservative F_side,const CeedScalar normal[3],CeedScalar u_side,CeedScalar s_side,CeedScalar s_star)215 CEED_QFUNCTION_HELPER StateConservative RiemannFlux_HLLC_Star(NewtonianIGProperties 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
RiemannFlux_HLLC_Star_fwd(NewtonianIGProperties 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)233 CEED_QFUNCTION_HELPER StateConservative RiemannFlux_HLLC_Star_fwd(NewtonianIGProperties 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)
RiemannFlux_HLLC(NewtonianIGProperties gas,State left,State right,const CeedScalar normal[3])271 CEED_QFUNCTION_HELPER StateConservative RiemannFlux_HLLC(NewtonianIGProperties 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
RiemannFlux_HLLC_fwd(NewtonianIGProperties gas,State left,State dleft,State right,State dright,const CeedScalar normal[3])298 CEED_QFUNCTION_HELPER StateConservative RiemannFlux_HLLC_fwd(NewtonianIGProperties 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