Coverage for dragonfly_sim/utils/NS_utils.py: 99%
78 statements
« prev ^ index » next coverage.py v7.16.2, created at 2026-10-17 03:37 +0000
1"""
2Navier-Stokes / RANS building blocks on top of Euler_utils (used by
3core/RANS_model.py and by CompressibleEulerModel whenever the state carries
4extra transported scalars).
6The state is the conservative mean flow U_m = (rho, rho*u, rho*E), n_mean =
72 + dim components, optionally followed by conservative transported scalars
8(rho*phi_k, e.g. rho*nu_tilde for Spalart-Allmaras). Euler_utils only ever
9sees the mean-flow slice; the scalars are convected passively here, with the
10mean flow's wave speeds.
12Every flux is written against an explicit n_mean, and reduces to the
13corresponding Euler_utils expression when there are no extra components.
15Nondimensionalization is the solver's: rho_inf = p_inf = 1, reference length
161, so T is only ever needed as the ratio (p/rho)/(p_inf/rho_inf). The heat
17flux is written in terms of p/rho, which removes the gas constant:
19 q = -kappa grad T = -(gamma/(gamma-1)) (mu/Pr + mu_t/Pr_t) grad(p/rho)
20"""
21import ufl
23from dragonfly_sim.utils.Euler_utils import (primitives, flux,
24 signed_wave_speeds, llf_penalty_term)
27def mean_flow(U, n_mean):
28 """The mean-flow slice of the state expression U."""
29 if U.ufl_shape[0] == n_mean:
30 return U
31 return ufl.as_vector([U[i] for i in range(n_mean)])
34def extra_components(U, n_mean):
35 """The transported scalars rho*phi_k of U, as a list."""
36 return [U[i] for i in range(n_mean, U.ufl_shape[0])]
39def concat_state(U_mean, extras):
40 """A full state from a mean-flow vector and a list of scalars."""
41 if len(extras) == 0:
42 return U_mean
43 return ufl.as_vector([U_mean[i] for i in range(U_mean.ufl_shape[0])] + list(extras))
46# ----------------------------------------------------------------------
47# Inviscid fluxes
48# ----------------------------------------------------------------------
49def extended_flux(U, gamma, n_mean):
50 """The inviscid flux tensor, (n_state, dim): Euler rows, then rho*phi*u."""
51 if U.ufl_shape[0] == n_mean:
52 return flux(U, gamma)
53 Um = mean_flow(U, n_mean)
54 _, u, _ = primitives(Um)
55 Fm = flux(Um, gamma)
56 rows = [Fm[i, :] for i in range(n_mean)]
57 rows += [s * u for s in extra_components(U, n_mean)]
58 return ufl.as_tensor(rows)
61def extended_normal_flux(U, n, gamma, n_mean):
62 return ufl.dot(extended_flux(U, gamma, n_mean), n)
65def extended_hll_flux(U_p, U_m, n, gamma, n_mean):
66 """
67 hll_flux on the full state. The wave speeds come from the mean flow
68 alone, and the same single-intermediate-state formula is applied to
69 every component, so the mean-flow rows are exactly hll_flux and each
70 scalar is upwinded with the fan of the flow carrying it.
71 """
72 sl_p, sr_p = signed_wave_speeds(mean_flow(U_p, n_mean), n, gamma)
73 sl_m, sr_m = signed_wave_speeds(mean_flow(U_m, n_mean), n, gamma)
75 SL_minus = ufl.min_value(ufl.min_value(sl_p, sl_m), 0.0)
76 SR_plus = ufl.max_value(ufl.max_value(sr_p, sr_m), 0.0)
77 denom = ufl.max_value(SR_plus - SL_minus, 1e-8)
79 return (
80 SR_plus * extended_normal_flux(U_p, n, gamma, n_mean)
81 - SL_minus * extended_normal_flux(U_m, n, gamma, n_mean)
82 + SL_minus * SR_plus * (U_m - U_p)
83 ) / denom
86def extended_boundary_flux(U, U_ext, n, gamma, n_mean):
87 """Euler_utils.boundary_flux on the full state."""
88 alpha = llf_penalty_term(mean_flow(U, n_mean), mean_flow(U_ext, n_mean), n, gamma)
89 return extended_normal_flux(U, n, gamma, n_mean) + alpha * (U - U_ext)
92# ----------------------------------------------------------------------
93# Gradients of the primitive quantities
94# ----------------------------------------------------------------------
95def primitive_gradients(U, G, gamma, n_mean):
96 """
97 grad(u), grad(p/rho) and grad(phi_k) from the state U and a tensor G
98 standing in for grad(U), shape (n_state, dim).
100 G is an argument rather than ufl.grad(U) because the interior-penalty
101 terms evaluate the viscous flux with the jump tensor [U] (x) n in
102 place of the gradient. Everything returned is linear in G.
103 """
104 d = n_mean - 2
105 rho = U[0]
106 u = ufl.as_vector([U[1 + i] / rho for i in range(d)])
107 rhoE = U[1 + d]
108 p = (gamma - 1.0) * (rhoE - 0.5 * rho * ufl.dot(u, u))
110 grad_rho = ufl.as_vector([G[0, j] for j in range(d)])
111 grad_m = [ufl.as_vector([G[1 + i, j] for j in range(d)]) for i in range(d)]
112 grad_rhoE = ufl.as_vector([G[1 + d, j] for j in range(d)])
114 grad_u = ufl.as_tensor([(grad_m[i] - u[i] * grad_rho) / rho for i in range(d)])
115 grad_p = (gamma - 1.0) * (grad_rhoE - sum(u[i] * grad_m[i] for i in range(d))
116 + 0.5 * ufl.dot(u, u) * grad_rho)
117 grad_p_over_rho = (grad_p - (p / rho) * grad_rho) / rho
119 grad_phi = []
120 for k in range(n_mean, U.ufl_shape[0]):
121 phi = U[k] / rho
122 grad_s = ufl.as_vector([G[k, j] for j in range(d)])
123 grad_phi.append((grad_s - phi * grad_rho) / rho)
125 return grad_u, grad_p_over_rho, grad_phi, grad_rho
128def vorticity_magnitude(grad_u, eps=1e-16):
129 """
130 sqrt(2 W:W), W the antisymmetric part of grad(u). Regularized by eps
131 inside the root so its derivative stays finite where the flow is
132 irrotational (the whole free stream).
133 """
134 W = 0.5 * (grad_u - grad_u.T)
135 return ufl.sqrt(2.0 * ufl.inner(W, W) + eps)
138# ----------------------------------------------------------------------
139# Viscous flux
140# ----------------------------------------------------------------------
141def viscous_stress(grad_u, mu_eff):
142 d = grad_u.ufl_shape[0]
143 return mu_eff * (grad_u + grad_u.T - (2.0 / 3.0) * ufl.tr(grad_u) * ufl.Identity(d))
146def viscous_flux(U, G, gamma, n_mean, mu, mu_t, Pr, Pr_t,
147 extra_fluxes=None, heat_flux=True):
148 """
149 The viscous flux tensor F_v(U, G), shape (n_state, dim), linear in G:
151 mass 0
152 momentum tau = (mu + mu_t)(grad u + grad u^T - 2/3 div u I)
153 energy tau u - q
154 scalars extra_fluxes(grad_phi, grad_rho), one row per scalar
156 `mu` and `mu_t` are expressions of the state (not of G). heat_flux=False
157 drops the conduction term, which is how the adiabatic wall imposes
158 q.n = 0.
159 """
160 grad_u, grad_T, grad_phi, grad_rho = primitive_gradients(U, G, gamma, n_mean)
161 d = n_mean - 2
162 rho = U[0]
163 u = ufl.as_vector([U[1 + i] / rho for i in range(d)])
165 tau = viscous_stress(grad_u, mu + mu_t)
166 energy = ufl.dot(tau, u)
167 if heat_flux:
168 energy = energy + gamma / (gamma - 1.0) * (mu / Pr + mu_t / Pr_t) * grad_T
170 rows = [ufl.as_vector([0.0] * d)]
171 rows += [ufl.as_vector([tau[i, j] for j in range(d)]) for i in range(d)]
172 rows += [energy]
173 if U.ufl_shape[0] > n_mean:
174 rows += list(extra_fluxes(grad_phi, grad_rho))
175 return ufl.as_tensor(rows)
178def sutherland_viscosity(U, gamma, n_mean, mu_inf, T_ratio_inf, S_ratio):
179 """
180 Sutherland's law in the solver's units:
182 mu/mu_inf = theta^(3/2) (1 + S) / (theta + S),
183 theta = T/T_inf = (p/rho) / (p_inf/rho_inf), S = S_dim / T_inf_dim
185 `T_ratio_inf` is p_inf/rho_inf. theta is floored so a transiently
186 negative pressure cannot hand theta^(3/2) a negative base.
187 """
188 Um = mean_flow(U, n_mean)
189 rho, u, E = primitives(Um)
190 p = (gamma - 1.0) * (E - 0.5 * rho * ufl.dot(u, u))
191 theta = ufl.max_value((p / rho) / T_ratio_inf, 1e-8)
192 return mu_inf * theta ** 1.5 * (1.0 + S_ratio) / (theta + S_ratio)