Coverage for dragonfly_sim/utils/Euler_utils.py: 76%
95 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-17 03:37 +0000
1import ufl
3# TODO: Return the velocity components all in one tuple, so that only the length of that tuple changes
5def primitives(U):
6 dim = U.ufl_shape[0] - 2
7 rho = U[0]
8 u = ufl.as_vector([U[i] / rho for i in range(1,1+dim)])
9 E = U[1 + dim]
10 return rho, u, E
12def pressure(U, gamma):
13 rho, u, E = primitives(U)
14 return (gamma - 1.0) * (E - 0.5 * rho * ufl.dot(u, u))
16def energy_density(p_in, rho_in, u_in, gamma):
17 kinetic = 0.5 * rho_in * ufl.dot(u_in, u_in)
18 E_new = p_in/(gamma - 1.0) + kinetic
19 return E_new
21def flux(U, gamma):
22 rho, u, E = primitives(U)
23 p = pressure(U, gamma)
25 d = u.ufl_shape[0]
26 I = ufl.Identity(d)
28 # Build block structure
29 # Shape: (d+2, d)
30 rows = []
32 # mass
33 rows.append(rho * u)
35 # momentum (d rows)
36 for i in range(d):
37 rows.append(rho * u[i] * u + p * I[i, :])
39 # energy
40 rows.append((E + p) * u)
42 return ufl.as_tensor(rows)
43 # F = ufl.as_vector([rho*u[0], rho*u[0]*u[0] + p, rho*u[0]*u[1], (E+p)*u[0]])
44 # G = ufl.as_vector([rho*u[1], rho*u[0]*u[1], rho*u[1]*u[1] + p, (E+p)*u[1]])
45 # return ufl.as_vector([F, G]).T
47def normal_flux(U, n_normal, gamma):
48 flux_U = flux(U, gamma)
49 return ufl.dot(flux_U, n_normal)
52def slip_wall_state(U, n):
53 rho, u, E = primitives(U)
55 u_normal = ufl.dot(u, n)
57 u_ref = u - 2 * u_normal * n
59 return ufl.as_vector([
60 U[0],
61 *[rho*u_ref[i] for i in range(u_ref.ufl_shape[0])],
62 U[1 + u_ref.ufl_shape[0]]
63 ])
65def llf_penalty_term(U_p, U_m, n, gamma):
66 rho_p, u_p, _ = primitives(U_p)
67 p_p = pressure(U_p, gamma)
68 c_p = ufl.sqrt(gamma * p_p / rho_p)
69 lam_p = [ufl.dot(u_p, n) - c_p, ufl.dot(u_p, n), ufl.dot(u_p, n) + c_p]
71 rho_m, u_m, _ = primitives(U_m)
72 p_m = pressure(U_m, gamma)
73 c_m = ufl.sqrt(gamma * p_m / rho_m)
74 lam_m = [ufl.dot(u_m, n) - c_m, ufl.dot(u_m, n), ufl.dot(u_m, n) + c_m]
76 def max_reduce(items):
77 res = items[0]
78 for item in items[1:]:
79 res = ufl.max_value(res, item)
80 return res
82 eigen_vals_max_p = max_reduce([abs(l) for l in lam_p])
83 eigen_vals_max_m = max_reduce([abs(l) for l in lam_m])
84 return ufl.max_value(eigen_vals_max_p, eigen_vals_max_m) / 2.0
87def llf_flux(U_p, U_m, n, gamma):
88 """
89 Local Lax-Friedrichs (LLF) numerical flux, as implemented in dolfin_dg.
90 Mathematically, this evaluates to the same expression as the Rusanov flux
91 for the compressible Euler equations.
92 """
93 alpha = llf_penalty_term(U_p, U_m, n, gamma)
95 return 0.5*(normal_flux(U_p, n, gamma) + normal_flux(U_m, n, gamma)) \
96 + alpha*(U_p - U_m)
98def signed_wave_speeds(U, n, gamma):
99 """
100 The signed extremal acoustic wave speeds (u.n - c, u.n + c) of the
101 normal flux Jacobian.
103 The Lax-Friedrichs family only ever needs the largest wave speed in
104 magnitude, |u.n| + c, whereas an approximate Riemann solver needs to
105 know which *direction* the fan is travelling in -- that is precisely
106 what lets it become fully one-sided in supersonic flow instead of
107 always dissipating both ways. So this returns the signed quantity that
108 the magnitude throws away.
110 Floors rho/p at 1e-8 before c = sqrt(gamma*p/rho). Worth being
111 explicit about why this floor is not merely cosmetic here: sqrt of a
112 negative argument returns NaN, not just an implausibly large number,
113 so a single transiently negative pointwise pressure -- which a DG
114 discretization can readily produce in early, far-from-converged
115 Newton iterations -- would otherwise poison the entire assembled
116 residual rather than just the facet it occurred on.
118 Works for facet-restricted traces (U('+')/U('-')) exactly like
119 primitives()/pressure().
120 """
121 rho, u, _ = primitives(U)
122 p = pressure(U, gamma)
123 rho_safe = ufl.max_value(rho, 1e-8)
124 p_safe = ufl.max_value(p, 1e-8)
125 c = ufl.sqrt(gamma * p_safe / rho_safe)
126 un = ufl.dot(u, n)
127 return un - c, un + c
129def hll_flux(U_p, U_m, n, gamma):
130 """
131 HLL (Harten-Lax-van Leer) approximate Riemann solver, a drop-in
132 alternative to llf_flux for flows containing shocks.
134 The Riemann fan is bounded by two waves with Davis speed estimates
136 SL = min(un_p - c_p, un_m - c_m)
137 SR = max(un_p + c_p, un_m + c_m)
139 and the flux is the corresponding single averaged intermediate state
141 SL >= 0 : F(U_p).n (fully upwind)
142 SR <= 0 : F(U_m).n (fully upwind)
143 else : (SR F(U_p).n - SL F(U_m).n + SL SR (U_m - U_p))
144 / (SR - SL)
146 Written branchlessly below via SL_m = min(SL, 0), SR_p = max(SR, 0),
147 which reproduces all three cases *exactly*: SL >= 0 sends SL_m to 0
148 and collapses the quotient to F(U_p).n, SR <= 0 sends SR_p to 0 and
149 collapses it to F(U_m).n. This is not merely a tidier spelling. It
150 needs only min_value/max_value rather than nested ufl.conditional
151 nodes, so the expression is genuinely continuous across the branch
152 boundaries instead of only piecewise defined -- which is what keeps
153 the assembled Jacobian consistent with the residual for a scheme
154 whose branch boundaries (un = +/- c) sit exactly on the sonic lines
155 of a transonic solution. That matters twice over here, since this
156 model's Jacobian also feeds the adjoint used for shape optimization.
158 Verified numerically against an explicit three-branch reference in
159 test_algorithms/test_hll_flux.py -- branchless vs branched agree to
160 ~1e-16 with states swept across sub-, trans- and supersonic so that
161 every branch is exercised, in both 2D and 3D.
163 Why HLL rather than Roe, and a caveat about dissipation
164 -------------------------------------------------------
165 HLL is the robust member of the characteristic-scheme family: it is
166 positivity-preserving given sound wave-speed estimates, and it
167 admits neither the carbuncle phenomenon nor entropy-violating
168 expansion shocks -- both of which a Roe flux needs an explicit
169 entropy fix to suppress. It also needs no eigenvector decomposition,
170 no Roe average and no entropy fix at all, which makes it markedly
171 cheaper to assemble and, more importantly here, to differentiate.
173 It is NOT more dissipative than the llf_flux it is offered as an
174 alternative to -- the opposite. Rusanov/LLF damps every
175 characteristic field at the speed of the fastest one (|u.n| + c),
176 which upper-bounds any consistent Riemann solver, so LLF is the
177 dissipation ceiling of this family and nothing swapped in here can
178 exceed it. Measured mean facet dissipation |F_hat - {F.n}| over
179 random state pairs, relative to LLF:
181 M~0.3 M~0.85 M~1.3
182 LLF 1.00x 1.00x 1.00x
183 HLL (Davis) 0.86x 0.75x 0.72x
184 Roe + entropy fix 0.81x 0.67x 0.67x
185 Van Leer FVS 0.52x 0.54x 0.61x
187 The exact ratios depend on how the trace pairs are sampled -- the
188 test's own sampling reports HLL at 0.89x/0.81x/0.80x, for instance --
189 so it is the *ordering* that is the invariant here, and that is what
190 test_hll_flux.py asserts rather than the individual numbers.
192 The consequence to watch for: with HLL the sensor-driven volume
193 viscosity in this model is doing relatively more of the shock
194 capturing than it was under LLF, and those constants were tuned
195 against LLF. If oscillations appear near a shock after switching,
196 the viscosity constants are the thing to revisit, not this flux.
198 Wave-speed estimates: Davis. The Einfeldt (HLLE) variant, which
199 additionally folds in the Roe-average speeds, measured within 3% of
200 Davis on the same states (0.85x/0.72x/0.70x) -- not worth pulling
201 the Roe-average machinery in for. It remains the first thing to try
202 if a near-vacuum positivity failure ever surfaces, since Einfeldt's
203 estimate is the one carrying the positivity proof.
205 Signature and orientation convention match llf_flux (n
206 is the "+"-side outward normal, so U_p is the upstream state for
207 u.n > 0), so this can be handed straight to
208 CompressibleEulerModel(..., flux_function=hll_flux). Dimension
209 agnostic: written against n alone, the expression is identical in 2D
210 and 3D.
211 """
212 sl_p, sr_p = signed_wave_speeds(U_p, n, gamma)
213 sl_m, sr_m = signed_wave_speeds(U_m, n, gamma)
215 SL = ufl.min_value(sl_p, sl_m)
216 SR = ufl.max_value(sr_p, sr_m)
218 SL_minus = ufl.min_value(SL, 0.0)
219 SR_plus = ufl.max_value(SR, 0.0)
221 # Non-negative by construction, and zero only in the degenerate
222 # SL = SR = 0 vacuum case -- the floor is essentially never active,
223 # it is here so that case yields a finite number instead of a 0/0.
224 denom = ufl.max_value(SR_plus - SL_minus, 1e-8)
226 return (
227 SR_plus * normal_flux(U_p, n, gamma)
228 - SL_minus * normal_flux(U_m, n, gamma)
229 + SL_minus * SR_plus * (U_m - U_p)
230 ) / denom
232def hllc_flux(U_p, U_m, n, gamma):
233 """
234 HLLC approximate Riemann solver: HLL with the contact wave restored.
236 hll_flux collapses the entire Riemann fan into a single averaged
237 intermediate state, which throws the contact wave away entirely.
238 HLLC reinstates it as a third wave of speed S*, so the two
239 intermediate states either side of the contact are resolved
240 separately. Wave speeds SL/SR come from the same
241 signed_wave_speeds/Davis estimate hll_flux uses -- deliberately, so
242 the two schemes differ only in the contact treatment -- and the
243 contact speed is
245 S* = (p_m - p_p + rho_p un_p (SL - un_p)
246 - rho_m un_m (SR - un_m))
247 / (rho_p (SL - un_p) - rho_m (SR - un_m))
249 with the star state on side K (Toro eq. 10.39; E is the volumetric
250 total energy, as primitives() returns it)
252 U*_K = rho_K (S_K - un_K)/(S_K - S*) *
253 [ 1,
254 u_K + (S* - un_K) n,
255 E_K/rho_K + (S* - un_K)(S* + p_K/(rho_K (S_K - un_K))) ]
257 What this buys over hll_flux
258 ----------------------------
259 Exact resolution of isolated contact discontinuities -- the fields
260 HLL still smears. On a stationary contact (density jump, uniform
261 pressure, zero normal velocity) this returns the exact physical
262 flux to machine precision (~4e-15), where hll_flux and llf_flux are
263 both O(1) wrong on the same states. In practice that means sharper
264 contacts, shear layers and boundary layers.
266 Branch structure: one conditional, and it is continuous
267 ------------------------------------------------------
268 Textbook HLLC is four-branched. Absorbing the two outer branches
269 into min(SL, 0) / max(SR, 0) -- the same trick hll_flux uses to
270 become fully branchless -- leaves exactly one select, on the sign
271 of S*:
273 F*_L = F(U_p).n + min(SL, 0) (U*_L - U_p)
274 F*_R = F(U_m).n + max(SR, 0) (U*_R - U_m)
275 F = conditional(S* >= 0, F*_L, F*_R)
277 SL >= 0 forces S* >= 0 (since SL <= S* <= SR) and collapses F*_L to
278 F(U_p).n; SR <= 0 forces S* <= 0 and collapses F*_R to F(U_m).n.
280 That the remaining conditional is *continuous* is what makes it an
281 acceptable concession here. At S* = 0 the two star fluxes are
282 equal: the mass flux rho* S*, the tangential momentum and the
283 energy all vanish, and the normal momentum is p*, which is equal on
284 both sides by construction. So unlike a generic branch, this one
285 introduces no jump in the assembled Jacobian -- which matters
286 because that Jacobian also feeds the adjoint used for shape
287 optimization. Verified by bisecting states onto S* = 0
288 (|F*_L - F*_R| ~ 5e-15) in test_algorithms/test_hllc_flux.py, along
289 with exact agreement against an explicit four-branch reference.
291 Cost relative to hll_flux
292 -------------------------
293 Roughly 1.5x the flux arithmetic. Counting arithmetic nodes in the
294 lowered expression graph (unique-node traversal, so common
295 subexpressions are collapsed as code generation would) and weighting
296 by approximate operation latency, the 2D interior facet term costs
297 1.53x for the residual and 1.49x for the Jacobian; the ratio moves
298 by under 5% across weightings from "all ops equal" to
299 "div=10, sqrt=20", so it is not an artifact of the weights. 3D is
300 the same picture (1.49x / 1.43x).
302 The extra cost is divisions, not the conditional. Division count in
303 the 2D residual goes 7 (LLF) -> 7 (HLL) -> 13 (HLLC), and in the
304 Jacobian 16 -> 16 -> 28: the star-state construction divides once
305 for S*, then per side for rho(S_K - un)/(S_K - S*),
306 p/(rho(S_K - un)) and E/rho. Square roots stay at 2 throughout (the
307 two sound speeds). For scale, hll_flux costs only 1.03x/1.05x over
308 llf_flux -- HLLC is the first flux here with a real price tag.
310 Read that 1.5x as a flux-arithmetic upper bound, NOT as an
311 end-to-end slowdown. Diluted into the full inviscid residual
312 (volume term plus facet term) it is already 1.44x/1.41x in 2D, and
313 in the assembled kernel it dilutes further against geometry, basis
314 tabulation and the quadrature loop -- all identical between the two
315 -- while this model's stabilization terms are unchanged and dilute
316 it again. Per Newton step the linear solve typically dominates
317 assembly outright. This is a static arithmetic proxy, produced
318 without running anything; it models neither memory traffic,
319 vectorization, FFCx's own optimizations, nor wall time.
321 Caveat worth knowing
322 --------------------
323 HLL's contact smearing is exactly what makes it carbuncle-free.
324 Restoring the contact reintroduces HLLC's standard susceptibility
325 to the carbuncle/odd-even instability at strong grid-aligned bow
326 shocks. Unlikely on a transonic airfoil c-grid, but it is the
327 mechanism to suspect if a strong normal shock ever develops a
328 ragged, grid-locked profile -- and hll_flux stays available as the
329 immediate fallback.
331 Signature and orientation convention match hll_flux/llf_flux (n is
332 the "+"-side outward normal, so U_p is the upstream state for
333 u.n > 0), so this can be handed straight to
334 CompressibleEulerModel(..., flux_function=hllc_flux). Dimension
335 agnostic: written against n alone, identical in 2D and 3D.
336 """
337 sl_p, sr_p = signed_wave_speeds(U_p, n, gamma)
338 sl_m, sr_m = signed_wave_speeds(U_m, n, gamma)
340 SL = ufl.min_value(sl_p, sl_m)
341 SR = ufl.max_value(sr_p, sr_m)
343 rho_p, u_p, E_p = primitives(U_p)
344 rho_m, u_m, E_m = primitives(U_m)
345 p_p = pressure(U_p, gamma)
346 p_m = pressure(U_m, gamma)
347 un_p = ufl.dot(u_p, n)
348 un_m = ufl.dot(u_m, n)
350 # SL - un_p <= -c_p < 0 and SR - un_m >= c_m > 0, so this denominator
351 # is strictly negative analytically. Clamp it *away from zero on the
352 # negative side* rather than flooring toward zero, so a degenerate
353 # state cannot flip the sign of S* -- which would select the wrong
354 # star flux rather than merely perturb its magnitude.
355 S_star = (
356 p_m - p_p + rho_p * un_p * (SL - un_p) - rho_m * un_m * (SR - un_m)
357 ) / ufl.min_value(rho_p * (SL - un_p) - rho_m * (SR - un_m), -1e-8)
359 def star_state(rho, u, E, p, un, S_K, denom):
360 factor = rho * (S_K - un) / denom
361 momentum = u + (S_star - un) * n
362 e_star = E / rho + (S_star - un) * (S_star + p / (rho * (S_K - un)))
363 return ufl.as_vector([
364 factor,
365 *[factor * momentum[i] for i in range(u.ufl_shape[0])],
366 factor * e_star,
367 ])
369 # SL - S* <= 0 <= SR - S* by construction; clamp each on its own side
370 # for the same sign-preserving reason as above.
371 U_star_p = star_state(rho_p, u_p, E_p, p_p, un_p, SL,
372 ufl.min_value(SL - S_star, -1e-8))
373 U_star_m = star_state(rho_m, u_m, E_m, p_m, un_m, SR,
374 ufl.max_value(SR - S_star, 1e-8))
376 # min/max absorb the two fully-upwind outer branches: SL >= 0 zeroes
377 # the first correction, SR <= 0 zeroes the second.
378 F_star_p = normal_flux(U_p, n, gamma) + ufl.min_value(SL, 0.0) * (U_star_p - U_p)
379 F_star_m = normal_flux(U_m, n, gamma) + ufl.max_value(SR, 0.0) * (U_star_m - U_m)
381 return ufl.conditional(ufl.ge(S_star, 0.0), F_star_p, F_star_m)
383def boundary_flux(U, U_ext, n, gamma):
384 alpha = llf_penalty_term(U, U_ext, n, gamma)
385 return normal_flux(U, n, gamma) + alpha * (U - U_ext)