SLIDE 1
Central Schemes: a Powerful Black-Box-Solver for Nonlinear - - PowerPoint PPT Presentation
Central Schemes: a Powerful Black-Box-Solver for Nonlinear - - PowerPoint PPT Presentation
Central Schemes: a Powerful Black-Box-Solver for Nonlinear Hyperbolic PDEs Alexander Kurganov Southern University of Science and Technology, China and Tulane University, USA www.math.tulane.edu/ kurganov Supported by NSF and NSFC joint
SLIDE 2
SLIDE 3
Southern University of Science and Technology
Shenzhen (Special Economic Zone) Population: 20,000 (1980), 15,000,000 (2018) GDP per capita: Shenzhen 27,000$; China 16,000$ Nanshan District 50,000$ (higher than in Germany, Canada, Japan) SUSTech is the first “Western” major research university in China founded in 2011 Department of Mathematics was founded in 2015 http://math.sustc.edu.cn/?lang=en Size: 25 faculty (2018), 45-50 (within 3-5 years) Many tenure/tenure-track, long and short term visiting, post-doc and PhD positions are available
SLIDE 4
Finite-Volume Methods
1-D System:
qt + f(q)x = 0
Integrate it over the space-time control volumes [a, b] × [tn, tn+1] to
- btain the integral form:
b
- a
q(x, tn+1) dx =
b
- a
q(x, tn) dx −
tn+1
- tn
- f(q(b, t)) − f(q(a, t))
- dt
tn tn+1 a b x t
4
SLIDE 5
q n
j ≈
1 ∆x
- Cj
q(x, tn) dx : cell averages over Cj := (xj−1
2, xj+1 2)
This solution is approximated by a piecewise linear (conservative, second-order accurate, non-oscillatory) reconstruction:
- q n(x) = q n
j + (qx)n j (x − xj)
for x ∈ Cj
x
j
x
j−1
x
j+1
x
j+2
5
SLIDE 6
x
j
x
j−1
x
j+1
x
j+2
For example, (qx)n
j = minmod
- θ
q n
j −q n j−1
∆x ,
q n
j+1 −q n j−1
2∆x , θ
q n
j+1 −q n j
∆x
- θ ∈ [1, 2]
where the minmod function is defined as: minmod(z1, z2, ...) :=
minj{zj}, if zj > 0 ∀j maxj{zj}, if zj < 0 ∀j 0,
- therwise
6
SLIDE 7
Godunov-type upwind schemes are designed by integrating
qt + f(q)x = 0
- ver the space-time control volumes [xj−1
2, xj+1 2] × [tn, tn+1]
tn tn+1 x
j
x
j+1
x
j−1/2
x
j+1/2
x
j−1
7
SLIDE 8
tn tn+1 x
j
x
j+1
x
j−1/2
x
j+1/2
x
j−1
q n+1
j
= q n
j −
1 ∆x
tn+1
- tn
- f(q(xj+1
2, t)) − f(q(xj−1 2, t))
- dt
In order to evaluate the flux integrals on the RHS, one has to (approximately) solve the generalized Riemann problem. This may be hard or even impossible...
8
SLIDE 9
- tn+1
x
j+1
x
j−1/2
x
j+1/2
x
j−1
tn x
j
x
j+1
x
j−1/2
x
j+1/2
x
j−1
x
j
x t x q(x,t )
n
SLIDE 10
Nessyahu-Tadmor Scheme
The Nessyahu-Tadmor [Nessyahu, Tadmor; 1990] scheme is a Godunov- type central scheme. It is designed by integrating
qt + f(q)x = 0
- ver
the different set
- f
staggered space-time control volumes [xj, xj+1] × [tn, tn+1] containing the Riemann fans
x
j+1
tn tn+1 x
j
x
j+1/2
✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✄ ✄☎ ☎ ✆ ✆ ✝ ✝10
SLIDE 11
- tn+1
x
j+1
x
j−1/2
x
j+1/2
x
j−1
tn x
j
x
j+1
x
j−1/2
x
j+1/2
x
j−1
x
j
x t x q(x,t )
n
SLIDE 12
x
j+1
tn tn+1 x
j
x
j+1/2
✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✁✁✁✁✁✁✁✁✁✁✁✁✁✁ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂✁✂ ✄ ✄☎ ☎ ✆ ✆ ✝ ✝q n+1
j+1
2
= 1 ∆x
xj+1
- xj
- q n(x) dx −
1 ∆x
tn+1
- tn
- f(q(xj+1, t)) − f(q(xj, t))
- dt
Due to the finite speed of propagation, this can be reduced to:
q n+1
j+1
2
=
q n
j +q n j+1
2 + ∆x 8
- (qx)n
j − (qx)n j+1
- − ∆t
∆x
- f(q
n+1
2
j+1 ) − f(q n+1
2
j
)
SLIDE 13
Values of q at t = tn+1
2 are approximated using the Taylor expansion:
q
n+1
2
j
≈
q n(xj) + ∆t
2 qt(xj, tn)
q n(x) = q n
j + (qx)n j (x − xj) =
⇒
- q n(xj) = q n
j
- qt(xj, tn) = −f(q n
j )x
The space derivatives fx are computed using the (minmod) limiter:
f(q n
j )x = minmod
- θ
f(q n
j ) − f(q n j−1)
∆x ,
f(q n
j+1) − f(q n j−1)
2∆x , θ
f(q n
j+1) − f(q n j )
∆x
- 13
SLIDE 14
Higher-Order and Multidimensional Staggered Central Schemes
[Arminjon, Viallon, Madrane; 1997] [Jiang, Tadmor; 1998] [Liu, Tadmor; 1998] [Bianco, Puppo, Russo; 1999] [Levy, Puppo, Russo; 1999, 2000, 2002] [Lie, Noelle; 2000]
14
SLIDE 15
Central-Upwind Schemes
Goal: to reduce numerical dissipation of central schemes Example — Numerical Dissipation of the Staggered LxF Scheme qn+1
j+1
2
= qn
j+1 + qn j
2 − ∆t ∆x
- f(qn
j+1) − f(qn j )
- qn+1
j+1
2
− qn
j+1
2
+ ∆t ∆x
- f(qn
j+1) − f(qn j )
- =
qn
j+1 − 2qn j+1
2
+ qn
j
2 qn+1
j+1
2
− qn
j+1
2
∆t + f(qn
j+1) − f(qn j )
∆x = (∆x)2 8∆t · qn
j+1 − 2qn j+1
2
+ qn
j
(∆x/2)2
- As ∆t decreases, the numerical dissipation increases
- As ∆t ∼ (∆x)2, the LxF scheme is inconsistent
- As ∆t → 0, the numerical dissipation blows up
15
SLIDE 16
Central-Upwind Schemes
- Godunov-type finite-volume methods
- Central:
Riemann-problem-solver-free methods designed without tracking complicated nonlinear waves
- Upwind:
Use some information on wave propagation to reduce numerical dissipation and thus enhance the resolution of nonsmooth waves
- Can be applied as a “black-box” solver to hyperbolic (systems of)
PDEs
- Robust, efficient and highly accurate
16
SLIDE 17
Selected Scientific Applications
- Gas dynamics and granular hydrodynamics
- Compressible multi-fluid and multi-phase flows
- Dusty gas flows
- Shallow water equations and related models
- Propagation of pollutant models
- Traffic and pedestrian flow models
- Incompressible Navier-Stokes equations
- Buckley-Leverett equation and polymer systems
- Glacier growth models
- Combustion and stiff detonation waves
- Chemotaxis and haptotaxis models
- Elastic waves in heterogeneous media
- Optimal control problems governed by nonlinear hyperbolic PDEs
SLIDE 18
Some Industrial Application
- Central-upwind schemes have been implemented in the industrial open
source CFD software package OpenFOAM [Gasparini, Greenshields, Weller; 2008] OpenCFD Ltd – http://www.opencfd.co.uk 3-D, polyhedral meshes with collocated variables, complex geometry, complex physics Widely used in industry, consultancy and academia worldwide including the European Space Agency (ESA)
- Central-upwind schemes have been implemented at the National
Center of Computational Hydroscience and Engineering at the University
- f Mississippi, where they became the core part of the GIS-based
decision support systems for flood mitigation and management. These systems are currently being used by several US federal agencies
SLIDE 19
x
j
x
j−1
x
j+1
x
j+2
- q n(x) = q n
j + (qx)n j (x − xj)
for x ∈ Cj
q−
j+1
2
:= lim
x→xj+1
2
−
- q(x, tn) = qn
j + ∆x
2 (qx)n
j
q+
j+1
2
:= lim
x→xj+1
2
+
- q(x, tn) = qn
j+1 − ∆x
2 (qx)n
j+1
SLIDE 20
x
j
x
j−1
x
j+1
x
j+2
The discontinuities appearing at the reconstruction step at the interface points {xj+1
2} propagate at finite speeds estimated by:
a+
j+1
2
:= max
- λN
∂f
∂q (q−
j+1
2
)
- , λN
∂f
∂q (q+
j+1
2
)
- , 0
- a−
j+1
2
:= min
- λ1
∂f
∂q (q−
j+1
2
)
- , λ1
∂f
∂q (q+
j+1
2
)
- , 0
- λ1 < λ2 < . . . < λN: N eigenvalues of the Jacobian ∂f
∂q
SLIDE 21
Idea: Select control volumes according to the size
- f
each Riemann fan
xj−1/2 xj+1/2 xj tn tn+1 xj+1 xj−1
j−1/2 j+1/2
q q
int int
- xj−1
2 + a−
j−1
2
∆t, xj−1
2 + a+
j−1
2
∆t
- × [tn, tn+1]
- xj+1
2 + a−
j+1
2
∆t, xj+1
2 + a+
j+1
2
∆t
- × [tn, tn+1]
21
SLIDE 22
xj−1/2 xj+1/2 xj tn tn+1 xj+1 xj−1
qj
int
- xj−1
2 + a+
j−1
2
∆t, xj+1
2 + a−
j−1
2
∆t
- × [tn, tn+1]
22
SLIDE 23
Final Step: Projection onto the Original Grid A piecewise linear interpolant,
q int(x), reconstructed from the evolved
intermediate cell averages {q int
j
} and {q int
j+1
2
}, is projected back onto the
- riginal grid by averaging it over the intervals [xj−1
2, xj+1 2].
int j
q
j−1/2 int
q
j+1/2 int
q
n+1 n
t t xj xj+1/2 xj−1/2 xj+1 xj−1
23
SLIDE 24
int j
q
j−1/2 int
q
j+1/2 int
q xj xj+1/2 xj−1/2 xj+1 xj−1
New projected cell averages:
q n+1
j
= a+
j−1
2
∆t ∆x
q int
j−1
2
+
1 +
- a−
j−1
2
− a+
j+1
2
- ∆t
∆x
q int
j
− a−
j+1
2
∆t ∆x
q int
j+1
2
+(∆t)2 2∆x
(qx)int
j+1
2
a+
j+1
2
a−
j+1
2
− (qx)int
j−1
2
a+
j−1
2
a−
j−1
2
SLIDE 25
1-D Semi-Discrete Central-Upwind Scheme
d dtqj(tn) = lim
∆t→0
q n+1
j
−q n
j
∆t = a+
j−1
2
∆x lim
∆t→0q int j−1
2
− a−
j+1
2
∆x lim
∆t→0q int j+1
2
+ a−
j−1
2
− a+
j+1
2
∆x lim
∆t→0q int j
+ lim
∆t→0
q int
j
−q n
j
∆t
+ 1 2∆x lim
∆t→0
- ∆t
- (qx)int
j+1
2
a+
j+1
2
a−
j+1
2
− (qx)int
j−1
2
a+
j−1
2
a−
j−1
2
- We then substitute q int
j±1
2
, q int
j
and (qx)int
j±1
2
into here to obtain the 1-D semi-discrete central-upwind scheme (for details see [Kurganov, Lin; 2007])
25
SLIDE 26
d dtqj(t) = −
Hj+1
2(t) − Hj−1 2(t)
∆x The central-upwind numerical flux is:
Hj+1
2 =
a+
j+1
2
f(q−
j+1
2
) − a−
j+1
2
f(q+
j+1
2
) a+
j+1
2
− a−
j+1
2
+ a+
j+1
2
a−
j+1
2
q+
j+1
2
− q−
j+1
2
a+
j+1
2
− a−
j+1
2
− dj+1
2
The built-in “anti-diffusion” term is:
dj+1
2 = 1
2 lim
∆t→0
- ∆t(qx)int
j+1
2
- = minmod
q+
j+1
2
− q∗
j+1
2
a+
j+1
2
− a−
j+1
2
,
q∗
j+1
2
− q−
j+1
2
a+
j+1
2
− a−
j+1
2
The intermediate values q∗
j+1
2
are:
q∗
j+1
2
= lim
∆t→0q int j+1
2
= a+
j+1
2
q+
j+1
2
− a−
j+1
2
q−
j+1
2
−
- f(q+
j+1
2
) − f(q−
j+1
2
)
- a+
j+1
2
− a−
j+1
2 26
SLIDE 27
Remarks
- 1. dj+1
2 ≡ 0 corresponds to the central-upwind scheme from [Kurganov,
Noelle, Petrova; 2001]
- 2. dj+1
2 ≡ 0 and a+
j+1
2
≡ −a−
j+1
2
correspond to the original scheme from [Kurganov, Tadmor; 2000]
- 3. The (formal) order of the scheme is determined only by the order of
the piecewise polynomial reconstruction
q, used to compute the values q±
j+1
2
, and the order of the ODE solver
- 4. Recommended ODE solver for nonlinear hyperbolic problems is the 3-
stage third-order strong stability preserving (SSP) Runge-Kutta method [Shu, Osher; 1988, 1989] [Gottlieb, Shu, Tadmor; 2001] [Gottlieb, Ketcheson, Shu; 2011]
27
SLIDE 28
Central schemes have been successfully extended to more complicated problems arising in a wide variety of applications:
- Convection-diffusion equations
qt + f(q)x = νqxx
Chertock, Kurganov, Levy, Petrova, Tadmor, ...
- Balance laws
qt + f(q)x = S(q)
Castro, Chertock, Kurganov, Levy, Morales de Luna, Noelle, Pareschi, Petrova, Russo, ...
- Nonconservative hyperbolic systems
qt + f(q)x = B(q)qx
Castro, Cheng, Chertock, Karni, Kirr, Kurganov, Petrova, ...
28
SLIDE 29
1-D Gas Dynamics
∂ ∂t
ρ ρu E
+ ∂
∂x
ρu ρu2 + p u(E + p)
= 0
p = (γ − 1)
- E − ρu2
2
- : equation of state
ρ: density u: velocity p: pressure E: total energy γ = 1.4
29
SLIDE 30
Example — Moving Contact Wave Initial data: (ρ, u, p)(x, 0) =
- (1.4, 0.1, 1),
x < 0.5 (1.0, 0.1, 1), x > 0.5
30
SLIDE 31
Final time: t = 2
1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 0.4 0.5 0.6 0.7 0.8 0.9 1 OLD, 1-ORDER NEW, 1-ORDER 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 0.4 0.5 0.6 0.7 0.8 0.9 1 OLD, 2-ORDER NEW, 2-ORDER
31
SLIDE 32
Example — Stationary Contact Wave and Traveling Shock and Rarefaction Initial data: (ρ, u, p)(x, 0) =
- (1, −19.59745, 1000),
x < 0.8 (1, −19.59745, 0.01), x > 0.8 Final time: t = 0.012
32
SLIDE 33
1 2 3 4 5 6 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 REFERENCE OLD, 1-ORDER NEW, 1-ORDER 1 2 3 4 5 6 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 REFERENCE OLD, 2-ORDER NEW, 2-ORDER 1 2 3 4 5 6 0.7 0.75 0.8 0.85 REFERENCE OLD, 1-ORDER NEW, 1-ORDER 1 2 3 4 5 6 0.7 0.75 0.8 0.85 REFERENCE OLD, 2-ORDER NEW, 2-ORDER
SLIDE 34
2-D Finite-Volume Methods
qt + f(q)x + g(q)y = 0
Consider Cartesian cells Cj,k := [xj−1
2, xj+1 2] × [yk−1 2, yk+1 2] centered at
xj := j∆x, yk := k∆y. Consider a uniform grid in time with the time levels tn := n∆t.
q n
j,k ≈
1 ∆x∆y
- Cj,k
q(x, y, tn) dy dx : cell averages
These cell averages are then evolved to the next time level t = tn+1 in three consecutive steps: reconstruction, evolution and projection.
34
SLIDE 35
Reconstruction
Second-order piecewise linear interpolant:
- q(x, y, tn) =
- j,k
- q n
j,k + (qx)n j,k(x − xj) + (qy)n j,k(y − yk)
χj,k(x, y)
χj,k(x, y): the characteristic function of the interval Cj,k
(qx)n
j,k ≈ qx(xj, yk, tn),
(qy)n
j,k ≈ qy(xj, yk, tn)
(qx)n
j,k = minmod
- θ
q n
j+1,k −q n j,k
∆x ,
q n
j+1,k −q n j−1,k
2∆x , θ
q n
j,k −q n j−1,k
∆x
- (qy)n
j,k = minmod
- θ
q n
j,k+1 −q n j,k
∆y ,
q n
j,k+1 −q n j,k−1
2∆y , θ
q n
j,k −q n j,k−1
∆y
- θ ∈ [1, 2]
35
SLIDE 36
Reconstructed Point Values
qN
j,k := q n j,k + ∆y
2 (qy)n
j,k
qW
j,k := q n j,k − ∆x
2 (qx)n
j,k
qE
j,k := q n j,k + ∆x
2 (qx)n
j,k
qS
j,k := q n j,k − ∆y
2 (qy)n
j,k
36
SLIDE 37
Evolution
Consider a general quadrilateral
q n+1
D
= 1 |D|
- D
- q(x, y, tn) dx dy − 1
|D|
tn+1
- tn
- ∂D
[ηxf(q) + ηyg(q)] ds dt
η = (ηx, ηy)T: outer unit normal vector to ∂D
37
SLIDE 38
2-D upwinding is quite problematic even for 1-order schemes on Cartesian grid! Easy (robust, accurate, efficient) alternatives: staggered central and central-upwind schemes Multidimensional Staggered Central Schemes: [Arminjon, Viallon, Madrane; 1997] [Jiang, Tadmor; 1998] [Levy, Puppo, Russo; 1999, 2000, 2002] [Lie, Noelle; 2000]
38
SLIDE 39
2-D Central-Upwind Schemes
[Kurganov, Prugger, Wu; 2017] Consider a Cartesian cell Cj,k. The one-sided local speeds at the midpoint of the cell edges can be estimated by a+
j+1
2,k := max
- λN
∂f
∂q
- qW
j+1,k
- , λN
∂f
∂q
- qE
j,k
- , 0
- a−
j+1
2,k := min
- λ1
∂f
∂q
- qW
j+1,k
- , λ1
∂f
∂q
- qE
j,k
- , 0
- b+
j,k+1
2
:= max
- λN
∂g
∂q
- qS
j,k+1
- , λN
∂g
∂q
- qN
j,k
- , 0
- b−
j,k+1
2
:= min
- λ1
∂g
∂q
- qS
j,k+1
- , λ1
∂g
∂q
- qN
j,k
- , 0
- λ1 < λ2 < . . . < λN: eigenvalues of the Jacobians ∂f/∂q and ∂g/∂q
39
SLIDE 40
Using the one-sided speed bounds, we split the computational domain
Cj,k into the nonsymmetric subdomains:
40
SLIDE 41
Thanks to the finite speed of propagation, the evolved solution will remain smooth in the central subdomains Dj,k for all t ∈ [tn, tn+1).
41
SLIDE 42
Central Nessyahu-Tadmor-Type Evolution
Recall
q int
D
= 1 |D|
- D
- q(x, y, tn) dx dy − 1
|D|
tn+1
- tn
- ∂D
[ηxf(q) + ηyg(q)] ds dt
H12 :=
tn+1
- tn
A2
- A1
- ηx
12f(q) + ηy 12g(q)
- ds dt
≈ ∆t |A1A2| 2
- ηx
12
- f
- q(A1, tn+1
2)
- + f
- q(A2, tn+1
2)
- +ηy
12
- g
- q(A1, tn+1
2)
- + g
- q(A2, tn+1
2)
SLIDE 43
The solution at the points Ai, i = 1, 2, 3, 4 is smooth for t ∈ [tn, tn+1). The midpoint values are obtained using the Taylor expansion in time:
q(Ai, tn+1
2) = q(Ai, tn) − ∆t
2 [f(q(Ai, tn))x + g(q(Ai, tn))y] + O((∆t)2) ≈
q(Ai, tn) − ∆t
2
- (f(q)x)n
Ai + (g(q)y)n Ai
- (f(q)x)n
Ai = minmod
f(q n
j+1,k) − f(q n j,k)
∆x ,
f(q n
j,k) − f(q n j−1,k)
∆x
- (g(q)y)n
Ai = minmod
g(q n
j,k+1) − g(q n j,k)
∆y ,
g(q n
j,k) − g(q n j,k−1)
∆y
- 43
SLIDE 44
q int
D
= 1 |D|
- D
- q(x, y, tn) dx dy − 1
|D|
tn+1
- tn
- ∂D
[ηxf(q) + ηyg(q)] ds dt
q int
D
= q n
D − 1
|D| [H12 + H23 + H34 + H41]
44
SLIDE 45
Projection
We project the intermediate solution, realized in terms of: the cell averages {q int
D } and
the point values at A1, A2, A3 and A4:
qn+1
j±1
4,k±1 4
=
q(zj±1
4,k±1 4, tn) − ∆t
- (f(q)x)n
j,k + (g(q)y)n j,k
- nto the original uniform mesh.
To this end, we first use the evolved data to construct a conservative piecewise linear interpolant
- q int(x, y), and then integrate it over the
- riginal cells Cj,k to obtain the cell averages of q at the new time level
t = tn+1:
q n+1
j,k
= 1 ∆x∆y
- Cj,k
- q int(x, y) dx dy
45
SLIDE 46
Since Dj,k ⊂ Cj,k, we take the constant pieces in the central subdomains:
- q int(x, y) = q int
Dj,k
for (x, y) ∈ Dj,k
46
SLIDE 47
In the domains Dα,β with (α, β) = (j + 1
2, k + 1 2), (j, k + 1 2) or (j + 1 2, k),
the interpolant
q int(x, y) consists of four linear pieces that continuously
match along the segments connecting (xα, yβ) with the vertices of Dα,β.
47
SLIDE 48
The value of
q int at (xα, yβ), which we define by qn+1
α,β , is determined
from the conservation requirement: 1 |Dα,β|
- Dα,β
- q int(x, y) dx dy = q int
Dα,β
which guarantees the second order of accuracy and results in
qn+1
α,β
= 3q int
Dα,β −
|DN
α,β| + |DE α,β|
|Dα,β|
qn+1
α+1
4,β+1 4
− |DE
α,β| + |DS α,β|
|Dα,β|
qn+1
α+1
4,β−1 4
− |DS
α,β| + |DW α,β|
|Dα,β|
qn+1
α−1
4,β−1 4
− |DW
α,β| + |DN α,β|
|Dα,β|
qn+1
α−1
4,β+1 4 48
SLIDE 49
The constructed interpolant
q int(x, y) may be oscillatory.
In order to avoid appearance on new local extrema at (xα, yβ), we check whether (q(i))n+1
α,β
> max
- (q (i)) int
Dα,β, (q (i)) int Dα±1
2,β, (q (i)) int
Dα,β±1
2
- r
(q(i))n+1
α,β
< min
- (q (i)) int
Dα,β, (q (i)) int Dα±1
2,β, (q (i)) int
Dα,β±1
2
- for some component i of q, and then replace the corresponding piecewise
linear approximations with ( q(i)) int(x, y) = (q (i)) int
Dα,β
for (x, y) ∈ Dα,β
- Remark. This reconstruction correction procedure locally reduces the
- rder of the interpolant
q int to the first one, but this is the same clipping
effect as in the case of the minmod limiter.
49
SLIDE 50
2-D Semi-Discrete Central-Upwind Scheme
Similar to the 1-D case, a 2-D semi-discrete central-upwind scheme can be rigorously derived from the fully discrete scheme.
Rectangular Grid
[Kurganov, Petrova; 2001] [Kurganov, Noelle, Petrova; 2001] [Kurganov, Tadmor; 2002] [Kurganov, Lin; 2007]
Triangular Grid
[Kurganov, Petrova; 2005]
Quadrilateral and Polygonal Grids
Beljadid, Kurganov, Mohammadian, Seidou, Shirkhani
50
SLIDE 51
Scheme from [Kurganov, Tadmor; 2002]
d dtqj,k(t) = −
Hx
j+1
2,k(t) − Hx
j−1
2,k(t)
∆x −
Hy
j,k+1
2
(t) − Hy
j,k−1
2
(t) ∆y
Hx
j+1
2,k(t) =
a+
j+1
2,kf(qE
j,k) − a− j+1
2,kf(qW
j+1,k)
a+
j+1
2,k − a−
j+1
2,k
+ a+
j+1
2,ka−
j+1
2,k
a+
j+1
2,k − a−
j+1
2,k
- qW
j+1,k − qE j,k
- Hy
j,k+1
2
(t) = b+
j,k+1
2
g(qN
j,k) − b− j,k+1
2
g(qS
j,k+1)
b+
j,k+1
2
− b−
j,k+1
2
+ b+
j,k+1
2
b−
j,k+1
2
b+
j,k+1
2
− b−
j,k+1
2
- qS
j,k+1 − qN j,k
- 51
SLIDE 52
Gas Dynamics
∂ ∂t
ρ ρu ρv E
+ ∂
∂x
ρu ρu2 + p ρuv u(E + p)
+ ∂
∂y
ρv ρuv ρv2 + p v(E + p)
= 0
p = (γ − 1)
- E − ρ
2(u2 + v2)
- : equation of state
ρ: density u, v: x- and y-velocities p: pressure E: total energy, respectively γ = 1.4
52
SLIDE 53
Example — Explosion (ρ, u, v, p)
- (x,y,0)
=
(1.000, 0, 0, 1.0), x2 + y2 < 0.16 (0.125, 0, 0, 0.1),
- therwise
By t = 3.2 the circular contact curve typically develops instabilities. Therefore, this example is a good test for the amount of numerical dissipation present in the studied schemes.
53
SLIDE 54
0.3 0.6 0.9 1.2 0.3 0.6 0.9 1.2
0.1 0.12 0.14 0.16 0.18 0.2 0.22
0.3 0.6 0.9 1.2 0.3 0.6 0.9 1.2
0.1 0.12 0.14 0.16 0.18 0.2 0.22
Semi-discrete (left) vs. fully discrete (right)
54
SLIDE 55
Example — Implosion (ρ(x, y, 0), u(x, y, 0), v(x, y, 0), p(x, y, 0)) =
(0.125, 0, 0, 0.14), |x| + |y| < 0.15
(1.000, 0, 0, 1.00), otherwise. We set the solid wall boundary conditions. A jet of liquid is expected to emerge. However, the numerical dissipation present in many second-order schemes may smear the jet.
55
SLIDE 56
0.05 0.1 0.15 0.2 0.25 0.05 0.1 0.15 0.2 0.25
0.4 0.5 0.6 0.7 0.8 0.9 1
0.05 0.1 0.15 0.2 0.25 0.05 0.1 0.15 0.2 0.25
0.4 0.5 0.6 0.7 0.8 0.9 1
Semi-discrete (left) vs. fully discrete (right)
56
SLIDE 57
Back to 1-D Gas Dynamics
[Kurganov; in preparation] ∂ ∂t
ρ ρu E
+ ∂
∂x
ρu ρu2 + p u(E + p)
= 0
p = (γ − 1)
- E − ρu2
2
- : equation of state
ρ: density u: velocity p: pressure E: total energy γ = 1.4
57
SLIDE 58
int j
q
j−1/2 int
q
j+1/2 int
q
n+1 n
t t xj xj+1/2 xj−1/2 xj+1 xj−1
58
SLIDE 59
Idea: Modify the projection step int j
q
j−1/2 int
q
j+1/2 int
q xj xj+1/2 xj−1/2 xj+1 xj−1
- The solution is averaged over the Riemann fans
- All conservative variables remain continuous in the cell (xℓ
j+1
2
, xr
j+1
2
) This brings excessive numerical dissipation!
59
SLIDE 60
In principle, a better approximation of qint in the cell (xℓ
j+1
2
, xr
j+1
2
) can be obtained by incorporating the wave propagation information into the interpolant. However, this approach will require the solution of the (generalized) Riemann problem. Alternatively, recall that our goal is to derive a semi-discrete scheme, that is, to pass to the ∆t → 0 limit, in which case the interval (xℓ
j+1
2
, xr
j+1
2
) shrinks into a point where qint may have at most two
- ne-sided values.
We therefore replace the intermediate piece q int
j+1
2
with two constant pieces, q int,L
j+1
2
and q int,R
j+1
2
.
60
SLIDE 61
That is, instead of
int j
q
j−1/2 int
q
j+1/2 int
q xj xj+1/2 xj−1/2 xj+1 xj−1
61
SLIDE 62
we perform the projection as follows:
int j
q
j−1/2 int
q
j+1/2 int
q
j+1/2 int,L
q
j+1/2 int,R
q xj xj+1/2 xj−1/2 xj+1 xj−1
62
SLIDE 63
For the compressible Euler equations:
q int,L
j+1
2
=
- ρ int,L
j+1
2
,m int,L
j+1
2
,E int,L
j+1
2
T
and q int,R
j+1
2
=
- ρ int,R
j+1
2
,m int,R
j+1
2
,E int,R
j+1
2
T
represent six pieces of information, which can be used to adjust the projection step. For instance, one can enforce continuity of the velocity and pressure (which are continuous across contact discontinuities!) by setting uint,L
j+1
2
= uint,R
j+1
2
, pint,L
j+1
2
= pint,R
j+1
2
This together with three conservation requirements,
q int,R
j+1
2
a+
j+1
2
−q int,L
j+1
2
a−
j+1
2
= q int
j+1
2
(a+
j+1
2
− a−
j+1
2
) result in five equations to be satisfied. The remaining degree of freedom can be used for obtaining a sharper approximation of qint.
63
SLIDE 64
For example, one can make the value of ρ int,R
j+1
2
−ρ int,L
j+1
2
as close as possible to ρ int
j+1 −ρ int j
without sacrificing the monotonicity of ρ:
xj xj+1/2 xj−1/2 xj+1 xj−1
j−1/2 int int j
ρ
j+1/2 int
ρ ρ
64
SLIDE 65
The new projection procedure leads to the same semi-discrete central- upwind scheme d dtqj(t) = −
Hj+1
2(t) − Hj−1 2(t)
∆x
Hj+1
2 =
a+
j+1
2
f(q−
j+1
2
) − a−
j+1
2
f(q+
j+1
2
) a+
j+1
2
− a−
j+1
2
+ a+
j+1
2
a−
j+1
2
a+
j+1
2
− a−
j+1
2
- q+
j+1
2
− q−
j+1
2
- − dj+1
2 65
SLIDE 66
but with the modified “anti-diffusion” term:
dj+1
2 = −minmod
- a+
j+1
2
(ρ+
j+1
2
− ρ∗
j+1
2
) , −a−
j+1
2
(ρ∗
j+1
2
− ρ−
j+1
2
)
-
1 u∗
j+1
2
(u∗
j+1 2
)2 2
q∗
j+1
2
= a+
j+1
2
q+
j+1
2
− a−
j+1
2
q−
j+1
2
−
- f(q+
j+1
2
) − f(q−
j+1
2
)
- a+
j+1
2
− a−
j+1
2 66
SLIDE 67
Example — Moving Contact Wave Initial data: (ρ, u, p)(x, 0) =
- (1.4, 0.1, 1),
x < 0.3 (1.0, 0.1, 1), x > 0.3
1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 NEW OLD
67
SLIDE 68
Example — Stationary Contact Wave, Traveling Shock and Rarefaction ∂ ∂t
ρ ρu E
+ ∂
∂x
ρu ρu2 + p u(E + p)
= 0
p = (γ − 1)
- E − ρu2
2
- Initial data:
(ρ, u, p)(x, 0) =
- (1, −19.59745, 1000),
x < 0.8 (1, −19.59745, 0.01), x > 0.8 Final time: t = 0.03
68
SLIDE 69
1 2 3 4 5 6
- 1
- 0.8
- 0.6
- 0.4
- 0.2
0.2 0.4 0.6 0.8 1 NEW OLD
69
SLIDE 70
1 2 3 4 5 6 0.75 0.8 0.85 0.9 REF NEW OLD
70
SLIDE 71
Example — “Shock-Bubble” Interaction The initial data correspond to a left-moving shock, initially located at x = 0.75, and a “bubble” with radius 0.25, initially located at the origin: (ρ, u, p)(x, y, 0) =
(13.1538, 0, 1), |x| < 0.25 (1.3333, −0.3535, 1.5), x > 0.75 (1, 0, 1),
- therwise
Solid wall boundary condition on the right. Final time t = 3.
71
SLIDE 72
Density, ∆x = 1/100
5 10 15 20
- 0.7
- 0.6
- 0.5
- 0.4
- 0.3
- 0.2
- 0.1
REF NEW OLD
72
SLIDE 73
Velocity, ∆x = 1/100
- 0.3
- 0.25
- 0.2
- 0.15
- 0.1
- 0.05
- 1
- 0.5
0.5 1 REF NEW OLD
73
SLIDE 74
Pressure, ∆x = 1/100
1.4 1.6 1.8 2 2.2
- 1
- 0.5
0.5 1 REF NEW OLD
74
SLIDE 75
Density, ∆x = 1/200
5 10 15 20
- 0.7
- 0.6
- 0.5
- 0.4
- 0.3
- 0.2
- 0.1
REF NEW OLD
75
SLIDE 76
Velocity, ∆x = 1/200
- 0.3
- 0.25
- 0.2
- 0.15
- 0.1
- 0.05
- 1
- 0.5
0.5 1 REF NEW OLD
76
SLIDE 77
Pressure, ∆x = 1/200
1.4 1.6 1.8 2 2.2
- 1
- 0.5
0.5 1 REF NEW OLD
77
SLIDE 78
Velocity, ∆x = 1/200
- 0.3
- 0.25
- 0.2
- 0.15
- 0.1
- 0.05
- 1
- 0.5
0.5 1 REF NEW
78
SLIDE 79
Pressure, ∆x = 1/200
1.4 1.6 1.8 2 2.2
- 1
- 0.5
0.5 1 REF NEW
79