Well-balanced DG scheme for Euler equations with gravity Praveen - - PowerPoint PPT Presentation
Well-balanced DG scheme for Euler equations with gravity Praveen Chandrashekar praveen@tifrbng.res.in Center for Applicable Mathematics Tata Institute of Fundamental Research Bangalore 560065 Higher Order Numerical Methods for Evolutionary
Praveen Chandrashekar praveen@tifrbng.res.in
Center for Applicable Mathematics Tata Institute of Fundamental Research Bangalore 560065
Higher Order Numerical Methods for Evolutionary PDEs: Applied Mathematics Meets Astrophysical Applications (15w5134) BIRS, Banff Centre, 10–15 May, 2015 Supported by Airbus Foundation Chair at TIFR-CAM, Bangalore
1 / 55
∂q ∂t + ∂f ∂x = − ρ ρu ∂Φ ∂x where q = ρ ρu E , f = ρu p + ρu2 (E + p)u , Φ = gravitational potential Calorically ideal gas p = (γ − 1)
2ρu2
γ = cp cv > 1
2 / 55
ue = 0
∂tρ + ∂x(ρu) = 0, ∂tE + ∂x[(E + p)u] = −ρu∂xΦ
dpe dx = −ρe dΦ dx (1)
3 / 55
p = ρRT, R = gas constant
pe(x) exp Φ(x) RTe
(2) Density ρe(x) = pe(x) RTe
4 / 55
peρ−ν
e
= α = const, ν > 1 constant (3) From (1) and (3), we obtain ναρν−1
e
(x) ν − 1 + Φ(x) = β = const (4)
5 / 55
◮ Gauss-Lobatto-Legendre nodes ◮ Arbitrary quadrilateral cells in 2-D
◮ isothermal solutions under ideal gas ◮ polytropic solutions (with Markus Zenk)
◮ Discontinuous density: contact preserving flux 6 / 55
Conservation law with source term stationary solution qe qt + f(q)x = s(q) f(qe)x = s(qe) Weak formulation: Find q(t) ∈ V such that d dt (q, ϕ) + a(q, ϕ) = (s(q), ϕ), ∀ϕ ∈ V FEM with quadrature: Find qh(t) ∈ Vh such that d dt (qh, ϕh)h + ah(qh, ϕh) = (sh(qh), ϕh)h, ∀ϕh ∈ Vh
7 / 55
Stationary solution is not a polynomial, qe / ∈ Vh. Let qe,h = Πh(qe), Πh : V → Vh, interpolation or projection FEM is well-balanced if ah(qe,h, ϕh) = (sh(qe,h), ϕh)h, ∀ϕh ∈ Vh because if qh(0) = qe,h = ⇒ qh(t) = qe,h ∀t
8 / 55
Ci = (xi− 1
2 , xi+ 1 2 ),
∆xi = xi+ 1
2 − xi− 1 2
Ii−1 Ii Ii−1
9 / 55
x = ξ∆xi + xi− 1
2
(5)
roots of (1 − ξ2)P ′
N(ξ) in [−1, +1]
ℓj(ξk) = δjk, 0 ≤ j, k ≤ N
ϕj(x) = ℓj(ξ), 0 ≤ j ≤ N
d dxϕj(x) = ℓ′
j(ξ) dξ
dx = 1 ∆xi ℓ′
j(ξ)
10 / 55
xj = ξj∆xi + xi− 1
2 ,
0 ≤ j ≤ N
11 / 55
Consider the single conservation law with source term ∂q ∂t + ∂f ∂x = s(q, x) Solution inside cell Ci is polynomial of degree N qh(x, t) =
N
qj(t)ϕj(x), qj(t) = qh(xj, t) Approximate the flux f(qh) ≈ fh(x, t) =
N
f(qh(xj, t))ϕj(x) =
N
fj(t)ϕj(x)
12 / 55
Gauss-Lobatto-Legendre quadrature
φ(x)ψ(x)dx ≈ (φ, ψ)h = ∆xi
N
ωqφ(xq)ψ(xq)
Semi-discrete DG: For 0 ≤ j ≤ N
d dt (qh, ϕj)h + (∂xfh, ϕj)h +[ ˆ fi+ 1
2 − fh(x−
i+ 1
2 )]ϕj(x−
i+ 1
2 )
−[ ˆ fi− 1
2 − fh(x+
i− 1
2 )]ϕj(x+
i− 1
2 ) = (sh, ϕj)h
(6) where ˆ fi+ 1
2 = ˆ
f(q−
i+ 1
2 , q+
i+ 1
2 ) is a numerical flux function. This is also
called a DG Spectral Element Method.
13 / 55
Consistency of numerical flux
ˆ f(q, q) = f(q)
Def: Contact property
The numerical flux ˆ f is said to satisfy contact property if for any two states wL = [ρL, 0, p] and wR = [ρR, 0, p] we have ˆ f(q(wL), q(wR)) = [0, p, 0]⊤
contact discontinuity.
⇒ exactly supports stationary contact
14 / 55
Let ¯ Ti = temperature corresponding to the cell average value in cell Ci Rewrite the source term in the momentum equation as (Xing & Shu) s = −ρ∂Φ ∂x = ρR ¯ Ti exp Φ R ¯ Ti ∂ ∂xexp
R ¯ Ti
sh(x) = ρh(x)R ¯ Ti exp Φ(x) R ¯ Ti ∂ ∂x
N
exp
R ¯ Ti
(7) Source term in the energy equation 1 ρh (ρu)hsh
15 / 55
Define H(x) inside each cell Ci as H(x) = ν ν − 1 ln ν − 1 ναi (βi − Φ(x))
x ∈ Ci αi, βi: constants to be chosen. Rewrite source term s(x) = −ρ∂Φ ∂x = ν − 1 ν ρ(βi − Φ(x)) exp(−H(x)) ∂ ∂xexp(H(x)), x ∈ Ci The source term is approximated as sh(x) = ν − 1 ν ρh(x)(βi − Φ(x)) exp(−H(x)) ∂ ∂x
N
exp(H(xj))ϕj(x) βi = max
0≤j≤N
ν − 1 pj ρj + Φ(xj)
αi = pj∗ρ−ν
j∗
16 / 55
Well-balanced property
Let the initial condition to the DG scheme (6), (7) be obtained by interpolating the isothermal/polytropic hydrostatic solution corresponding to a continuous gravitational potential Φ. Then the scheme (6), (7) preserves the initial condition under any time integration scheme. Proof: For continuous hydrostatic solution, qh(0) is continuous. By flux consistency ˆ fi+ 1
2 − fh(x−
i+ 1
2 ) = 0,
ˆ fi− 1
2 − fh(x+
i− 1
2 ) = 0
Above is true even if density is discontinuous, provided flux satisfies contact property. = ⇒ density and energy equations are well-balanced
17 / 55
Momentum eqn: flux fh has the form fh(x, t) =
N
pj(t)ϕj(x), pj = pressure at the GLL point xj Isothermal initial condition, ¯ Ti = Te = const. The source term evaluated
18 / 55
at any GLL node xk is given by sh(xk) = ρh(xk)RTe
exp Φ(xk) RTe N
exp
RTe ∂ ∂xϕj(xk) = pk exp Φ(xk) RTe N
exp
RTe ∂ ∂xϕj(xk) =
N
pk exp Φ(xk) RTe
RTe ∂ ∂xϕj(xk) =
N
pj exp Φ(xj) RTe
RTe ∂ ∂xϕj(xk) =
N
pj ∂ ∂xϕj(xk) = ∂ ∂xfh(xk)
19 / 55
Since ∂xfh(xk) = sh(xk) at all the GLL nodes xk we can conclude that (∂xfh, ϕj)h = (sh, ϕj)h , 0 ≤ j ≤ N = ⇒ scheme is well-balanced for the momentum equation The proof for polytropic case is similar.
20 / 55
∂q ∂t + ∂f ∂x + ∂g ∂y = s q = vector of conserved variables, (f, g) = flux vector and s = source term, given by q = ρ ρu ρv E , f = ρu p + ρu2 ρuv (E + p)u , g = ρv ρuv p + ρv2 (E + p)v s = −ρ ∂Φ
∂x
−ρ ∂Φ
∂y
−
∂x + v ∂Φ ∂y
21 / 55
Momentum equation ∇pe = −ρe∇Φ Assuming the ideal gas equation of state p = ρRT and a constant temperature T = Te = const, we get pe(x, y) exp Φ(x, y) RTe
(8) We will exploit the above property of the hydrostatic state to construct the well-balanced scheme.
22 / 55
A quadrilateral cell K and reference cell ˆ K = [0, 1] × [0, 1]
1 2 3 4 K x y ξ η 1 2 3 4 ˆ K TK
1-D GLL points ξr ∈ [0, 1], 0 ≤ r ≤ N
23 / 55
Tensor product of GLL points N = 1 N = 2 N = 3 Basis functions in cell K: i ↔ (r, s) ϕK
i (x, y) = ℓr(ξ)ℓs(η),
(x, y) = TK(ξ, η) Computation of ∂xfh requires derivatives of the basis functions ∂ ∂xϕK
i (x, y) = ℓ′ r(ξ)ℓs(η) ∂ξ
∂x + ℓr(ξ)ℓ′
s(η)∂η
∂x
24 / 55
Consider a single conservation law with source term ∂q ∂t + ∂f ∂x + ∂g ∂x = s Solution inside cell K qh(x, y, t) =
M
qK
i (t)ϕK i (x, y),
M = (N + 1)2 Approximate fluxes inside the cell by interpolation at the same GLL nodes fh(x, y) =
M
f(qh(xi, yi))ϕK
i (x, y) = M
fiϕK
i (x, y)
gh(x, y) =
M
g(qh(xi, yi))ϕK
i (x, y) = M
giϕK
i (x, y)
25 / 55
Quadrature on element K using GLL points (φ, ψ)K =
N
N
φ(ξr, ξs)ψ(ξr, ξs)ωrωs|JK(ξr, ξs)| Quadrature on the faces of the cell K (φ, ψ)∂K =
(φ, ψ)e where (·, ·)e are one dimensional quadrature rules using the subset of GLL nodes located on the boundary of the cell.
26 / 55
Semi-discrete DG scheme: For 1 ≤ i ≤ M
d dt
i
i
+
Fh − F −
h , ϕK i
i
Numerical flux function ˆ Fh = ˆ F(q−
h , q+ h , n),
n = (nx, ny), unit outward normal to ∂K Trace of flux on ∂K from cell K F −
h = f− h nx + g− h ny
27 / 55
Let ¯ TK = temperature corresponding to the cell average value in cell K Rewrite source term in the x momentum equation s = −ρ∂Φ ∂x = ρR ¯ TK exp
R ¯ TK ∂ ∂x exp
Φ R ¯ TK
sh(x, y) = ρh(x, y)R ¯ TK exp Φ(x, y) R ¯ TK ∂ ∂x
M
exp
R ¯ TK
j (x, y)
(10)
28 / 55
Well-balanced property
Let the initial condition to the DG scheme (9), (10) be obtained by interpolating the isothermal/polytropic hydrostatic solution corresponding to a continuous gravitational potential Φ. Then the scheme preserves the initial condition under any time integration scheme.
29 / 55
◮ If residual in cell K is zero, then dont apply limiter in that cell. 30 / 55
dt = R(t, q)
Second order accurate SSP Runge-Kutta scheme q(1) = qn + ∆tR(tn, qn) q(2) = 1 2qn + 1 2[q(1) + ∆tR(tn + ∆t, q(1))] qn+1 = q(2) Third order accurate SSP Runge-Kutta scheme q(1) = qn + ∆tR(tn, qn) q(2) = 3 4qn + 1 4[q(1) + ∆tR(tn + ∆t, q(1))] q(3) = 1 3qn + 2 3[q(2) + ∆tR(tn + ∆t/2, q(2))] qn+1 = q(3)
31 / 55
32 / 55
Potential Φ = x, Φ = sin(2πx) Initial state u = 0, ρ = p = exp(−Φ(x)) Mesh ρu ρ E 25x25 1.03822e-13 2.72604e-14 9.53913e-14 50x50 1.04783e-13 2.67559e-14 9.36725e-14 100x100 1.05019e-13 2.66323e-14 9.34503e-14 200x200 1.05088e-13 2.66601e-14 9.33861e-14
Table: Well-balanced test on Cartesian mesh using Q1 and potential Φ = y
33 / 55
Mesh ρu ρ E 25x25 1.04518e-13 2.7548e-14 9.64205e-14 50x50 1.04983e-13 2.69317e-14 9.43158e-14 100x100 1.05069e-13 2.69998e-14 9.39126e-14 200x200 1.05089e-13 2.68828e-14 9.462e-14
Table: Well-balanced test on Cartesian mesh using Q2 and potential Φ = x
Mesh ρu ρ E 25x25 9.23424e-13 2.31405e-13 8.16645e-13 50x50 9.36459e-13 2.28315e-13 8.04602e-13 100x100 9.39613e-13 2.28001e-13 8.03005e-13 200x200 9.40422e-13 2.2792e-13 8.02653e-13
Table: Well-balanced test on Cartesian mesh using Q1 and Φ = sin(2πx)
34 / 55
Mesh ρu ρ E 25x25 9.34536e-13 2.35173e-13 8.30316e-13 50x50 9.39556e-13 2.29908e-13 8.10055e-13 100x100 9.40442e-13 2.28538e-13 8.04638e-13 200x200 9.40668e-13 2.28051e-13 8.0357e-13
Table: Well-balanced test on Cartesian mesh using Q2 and Φ = sin(2πx)
35 / 55
Potential Φ(x) = x Initial condition u = 0, ρ = exp(−x), p = exp(−x) + η exp(−100(x − 1/2)2)
36 / 55
0.0 0.2 0.4 0.6 0.8 1.0
x
−0.002 0.000 0.002 0.004 0.006 0.008 0.010
Pressure perturbation Initial 100 cells 1000 cells
0.0 0.2 0.4 0.6 0.8 1.0
x
−0.002 0.000 0.002 0.004 0.006 0.008 0.010
Pressure perturbation Initial 50 cells 500 cells
(a) (b)
Figure: Evolution of perturbations for η = 10−2: (a) Q1 (b) Q2
37 / 55
0.0 0.2 0.4 0.6 0.8 1.0
x
−0.00002 0.00000 0.00002 0.00004 0.00006 0.00008 0.00010
Pressure perturbation Initial 100 cells 1000 cells
0.0 0.2 0.4 0.6 0.8 1.0
x
−0.00002 0.00000 0.00002 0.00004 0.00006 0.00008 0.00010
Pressure perturbation Initial 50 cells 500 cells
(a) (b)
Figure: Evolution of perturbations for η = 10−4: (a) Q1 (b) Q2
38 / 55
Potential Φ(x) = x and initial condition (ρ, u, p) =
x < 1
2
(0.125, 0, 0.1) x > 1
2
0.0 0.2 0.4 0.6 0.8 1.0
x
0.0 0.2 0.4 0.6 0.8 1.0 1.2
Density Q1, 100 cells Q2, 100 cells FVM, 2000 cells
0.0 0.2 0.4 0.6 0.8 1.0
x
0.0 0.2 0.4 0.6 0.8 1.0 1.2
Density Q1, 200 cells Q2, 200 cells FVM, 2000 cells
39 / 55
Potential Φ = x + y Hydrostatic solution is given by ρ = ρ0 exp
p0 (x + y)
p = p0 exp
p0 (x + y)
u = v = 0 ρu ρv ρ E Q1, 25 × 25 9.85926e-14 9.85855e-14 5.32357e-14 1.55361e-13 Q1, 50 × 50 9.94493e-14 9.94451e-14 5.37084e-14 1.56669e-13 Q1, 100 × 100 9.96481e-14 9.96474e-14 5.38404e-14 1.57062e-13 Q2, 25 × 25 9.9256e-14 9.92682e-14 5.39863e-14 1.57435e-13 Q2, 50 × 50 9.961e-14 9.96538e-14 5.41091e-14 1.57521e-13 Q2, 100 × 100 9.95889e-14 9.97907e-14 5.43145e-14 1.57728e-13
40 / 55
p = p0 exp
p0 (x + y)
p0 [(x − 0.3)2 + (y − 0.3)2]
x y
0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8
x y
0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8
(a) (b)
Figure: Pressure perturbation on 50 × 50 mesh at time t = 0.15 (a) Q1 (b) Q2. Showing 20 contours lines between −0.0002 to +0.0002
42 / 55
x y
0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8
x y
0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8
(a) (b)
Figure: Pressure perturbation on 200 × 200 mesh at time t = 0.15 (a) Q1 (b) Q2. Showing 20 contours lines between −0.0002 to +0.0002
43 / 55
Potential Φ(x, y) = y Temperature is discontinuous T =
y < 0 Tu, y > 0 Hydrostatic pressure and density p =
y ≤ 0 p0e−y/RTu, y > 0 , ρ =
y < 0 p/RTu, y > 0 Density is discontinuous at y = 0.
44 / 55
Case 1, Tl = 1, Tu = 2 Q1, 25x100 5.13108e-13 1.10971e-13 2.56836e-13 8.91784e-13 Q1, 50x200 5.15744e-13 1.12234e-13 2.84309e-13 8.97389e-13 Q2, 25x100 5.15726e-13 1.1341e-13 3.22713e-13 9.01183e-13 Q2, 50x200 5.16397e-13 1.13707e-13 3.93623e-13 9.02503e-13
Table: Well-balanced test for Rayleigh-Taylor problem
Case 2, Tl = 2, Tu = 1 Q1, 25x100 3.487e-13 1.0989e-13 1.98949e-13 7.42265e-13 Q1, 50x200 3.50153e-13 1.1121e-13 2.34991e-13 7.4747e-13 Q2, 25x100 3.50149e-13 1.12159e-13 2.80645e-13 7.5104e-13 Q2, 50x200 3.50548e-13 1.12533e-13 3.63806e-13 7.52151e-13
Table: Well-balanced test for Rayleigh-Taylor problem
45 / 55
Exact solution of Euler equations with gravity given by (Xing & Shu) ρ = 1 + 0.2 sin(π(x + y − t(u0 + v0))), u = u0, v = v0 p = p0 + t(u0 + v0) − x − y + 0.2 cos(π(x + y − t(u0 + v0)))/π Compute solution error in L2 norm at time t = 0.1
1/h ρu ρv ρ E Error Rate Error Rate Error Rate Error Rate 50 0.00134154 – 0.00134154 – 0.0012837 – 0.00161287 – 100 0.000335446 1.99 0.000335446 1.99 0.00032044 2.00 0.000411141 1.97 200 8.35627e-05 2.00 8.35627e-05 2.00 7.97842e-05 2.00 0.00010335 1.99 400 2.08348e-05 2.00 2.08348e-05 2.00 1.98754e-05 2.00 2.58109e-05 2.00
Table: Convergence of error for degree N = 1
1/h ρu ρv ρ E Error Rate Error Rate Error Rate Error Rate 25 7.7019e-05 – 7.7019e-05 – 7.80868e-05 – 9.32865e-05 – 50 9.68863e-06 2.99 9.68863e-06 2.99 9.76471e-06 2.99 1.16849e-05 2.99 100 1.21506e-06 2.99 1.21506e-06 2.99 1.22031e-06 3.00 1.46256e-06 2.99 200 1.52134e-07 2.99 1.52134e-07 2.99 1.52503e-07 3.00 1.8247e-07 3.00
Table: Convergence of error for degree N = 2
46 / 55
30 × 30 cells 956 cells
47 / 55
Radial potential Φ(x, y) = r =
Hydrostatic solution ρ = p = exp(−r) ρu ρv ρ E Q1, 30x30 6.08113e-13 6.06081e-13 2.11611e-14 4.13127e-14 Q1, 50x50 5.45634e-13 5.44973e-13 1.29319e-14 2.66376e-14 Q1, 100x100 5.4918e-13 5.49012e-13 9.63433e-15 2.11463e-14 Q2, 30x30 6.29776e-13 6.27278e-13 1.9777e-14 5.08689e-14 Q2, 50x50 5.5376e-13 5.52903e-13 1.5635e-14 3.85134e-14 Q2, 100x100 5.51645e-13 5.5142e-13 2.173e-14 5.25578e-14
Table: Well balanced test for radial Rayleigh-Taylor problem on Cartesian mesh
48 / 55
ρu ρv ρ E Q1, 956 3.03033e-16 3.23738e-16 6.66245e-16 2.11449e-16 Q1, 2037 5.20653e-16 5.10865e-16 9.82565e-16 4.06402e-16 Q1, 10710 2.00037e-15 1.57078e-15 2.21284e-15 4.66515e-15 Q2, 956 2.69429e-14 3.31591e-14 4.50499e-14 1.61455e-13 Q2, 2037 7.68549e-14 1.1834e-13 1.01632e-13 3.84886e-13 Q2, 10710 2.92633e-13 2.44323e-13 2.9449e-13 4.10069e-13
Table: Well balanced test for radial Rayleigh-Taylor problem on unstructured mesh
49 / 55
Initial pressure and density p =
r ≤ r0 e− r
α +r0 (1−α) α
r > r0 , ρ =
r ≤ ri
1 αe− r
α +r0 (1−α) α
r > ri where ri = r0(1 + η cos(kθ)), α = exp(−r0)/(exp(−r0) + ∆ρ) The density jumps by an amount ∆ρ > 0 at the interface defined by r = ri whereas the pressure is continuous. Following [2], we take ∆ρ = 0.1, η = 0.02, k = 20
50 / 55
Cartesian mesh of 240 × 240 and Q1 basis
51 / 55
v = [ρ exp(Φ/R ¯ Ti), u, p exp(Φ/R ¯ Ti)]⊤
vh = Πhv(qh)
d dt (qh, ϕj)h + (∂xfh(qh), ϕj)h +[{ ˆ f(q(v−
h ), q(v+ h )) − f(q(v− h ))}ϕ− j ]i+ 1
2
−[{ ˆ f(q(v−
h ), q(v+ h )) − f(q(v+ h ))}ϕ+ j ]i− 1
2
= (sh(qh), ϕj)h
52 / 55
53 / 55
54 / 55
[1] R. K¨ appeli and S. Mishra, Well-balanced schemes for the Euler equations with gravitation, J. Comput. Phys., 259 (2014),
[2] Randall. J. LeVeque and Derek. S. Bale, Wave propagation methods for conservation laws with source terms, in Hyperbolic Problems: Theory, Numerics, Applications, Rolf Jeltsch and Michael Fey, eds., vol. 130 of International Series of Numerical Mathematics, Birkh¨ auser Basel, 1999, pp. 609–618. [3] Yulong Xing and Chi-Wang Shu, High order well-balanced WENO scheme for the gas dynamics equations under gravitational fields, J. Sci. Comput., 54 (2013), pp. 645–662.
55 / 55