McGill University
McGill University 4th
High-order CFD Workshop
Farshad Navah PhD Candidate
Siva Nadarajah Associate Professor McGill University Computational Aerodynamics Group
June, 2016
Farshad Navah
4th High-order CFD Workshop
June, 2016
1 / 32
McGill University
Outline
1 High-order CFD Code
Code description 2 Code verification
NS solver verification via MMS Example of solution verification 3 BS1 - DNS of the Taylor-Green Vortex
Farshad Navah
4th High-order CFD Workshop
June, 2016
2 / 32
McGill University High-order CFD Code Code description
I
Governing equations: Compressible Navier-Stokes equations
I
Discretization scheme: high-order correction procedure via flux reconstruction (CPR)
I
Numerical Flux: Roe scheme for inviscid terms and BR2 for viscous terms
I
Divergence computation method: Chain rule (CR) for inviscid fluxes and the Lagrange polynomials (LP) for viscous fluxes
I
Solution method: Backward Euler → Full Newton
I
Nodes: GLL
I
Parallelization: Open-MPI
Farshad Navah
4th High-order CFD Workshop
June, 2016
3 / 32
McGill University High-order CFD Code Code description
Recent Developments (Farshad Navah) I
Governing equations: Compressible Reynolds-averaged Navier-Stokes equations closed by the negative Spalart-Allmaras (SA) turbulence model (ICCFD7-1902)
I
Solution method: 15-digits accurate analytical Jacobian of RANS-SA(pos/neg), verified via complex step
I
Code verification: Method of manufactured solutions (Euler, NS, RANS-SA)
Farshad Navah
4th High-order CFD Workshop
June, 2016
4 / 32
McGill University Code verification
Code verification in CFD
Farshad Navah
4th High-order CFD Workshop
June, 2016
5 / 32
McGill University Code verification
Simulation Step
Sources of Error in CFD
Reality ⇓
=⇒
Modelling error
Conceptual model ⇓
=⇒
Programming & Discretization errors
=⇒
Round-off & Iterative convergence errors
Numerical model ⇓ Numerical solution
Farshad Navah
4th High-order CFD Workshop
June, 2016
6 / 32
McGill University Code verification
V&V Model validation
Sources of Error in CFD ⇐⇒
Modelling error
⇑ Solution verification ⇐⇒
Discretization error
⇑ Code verification
⇐⇒
Programming errors
⇑ Solution process
Farshad Navah
⇐⇒ Round-off & Iterative convergence errors
4th High-order CFD Workshop
June, 2016
7 / 32
McGill University Code verification
Methods of Code Verification in CFD I
Method of analytical solutions Pros: -non-intrusive Cons: -limited range of models (no RANS solutions) -(often) over-simplified flows (ex: Couette flow)
I
Method of manufactured solutions (MMS) Pros: -Covers all possible models/flow regimes -Verifies targeted boundary conditions -Allows for debugging Cons: -Creation of a proper MS is delicate wrt to model validity/numerical stability, etc. -Deployment needs expertise
Farshad Navah
4th High-order CFD Workshop
June, 2016
8 / 32
McGill University Code verification NS solver verification via MMS
Examples of Code and Solution Verification
Focus: Discretization and Programming Errors ◦ Round-off error ◦ Iterative convergence error
Farshad Navah
)
−→ Residual norm is at least 3 orders of magnitude lower than error norm
4th High-order CFD Workshop
June, 2016
9 / 32
McGill University Code verification NS solver verification via MMS
Manufactured solution: ρM S = ρ0 + ρx sin(aρx πx/L)
+ ρy cos(aρy πy/L)
+ ρxy cos(aρxy πx/L) cos(aρxy πy/L)
UM S = u0 + ux sin(aux πx/L) + uy cos(auy πy/L) + uxy cos(auxy πx/L) cos(auxy πy/L) VM S =
v0 + vx cos(avx πx/L)
+ vy sin(avy πy/L)
+ vxy cos(avxy πx/L) cos(avxy πy/L)
PM S = p0 + px cos(apx πx/L)
+ py sin(apy πy/L)
+ pxy cos(apxy πx/L) cos(apxy πy/L)
EM S = PM S /((γ−1)ρM S ) +
1 2 (U +V 2 ) 2 MS MS
Domain: Ω = [0, 1]2
Grids: Series of doubling isotropic quads Farshad Navah
4th High-order CFD Workshop
June, 2016
10 / 32
McGill University Code verification NS solver verification via MMS
Viscous - (Full NS) Solution: subsonic
Viscosity: µ = 0.001
Boundary conditions: weak Dirichlet
Farshad Navah
4th High-order CFD Workshop
June, 2016
11 / 32
McGill University Code verification NS solver verification via MMS
ρE error distribution versus grid refinement for P 4
(a) 4 × 4
(b) 8 × 8
(d) 32 × 32 Farshad Navah
(c) 16 × 16
(e) 64 × 64
4th High-order CFD Workshop
June, 2016
12 / 32
McGill University Code verification NS solver verification via MMS
NS solver verification Order of accuracy - P1 3.0
3.0 kρ − ρex k1
kρu − ρuex k1
2.5
kρu − ρuex k2
kρ − ρex k∞
Order of accuracy
Order of accuracy
kρ − ρex k2
2.0
1.5
1.0
2.5
2.0
1.5
1.0 10−2
h= 3.0 kρv − ρvex k1
p
10−1
10−2
(1/N )
h= 3.0 kρE − ρEex k1
2.0
1.5
1.0
(1/N )
2.5
p
(1/N )
kρE − ρEex k∞
2.0
1.5
1.0 10−2
h=
Farshad Navah
10−1
p
kρE − ρEex k2
kρv − ρvex k∞
Order of accuracy
Order of accuracy
kρv − ρvex k2 2.5
kρu − ρuex k∞
p
10−1
(1/N )
10−2
h=
4th High-order CFD Workshop
10−1
June, 2016
13 / 32
McGill University Code verification NS solver verification via MMS
NS solver verification Order of accuracy - P2 4.0
4.0 kρ − ρex k1
kρu − ρuex k1
3.5
kρu − ρuex k2
kρ − ρex k∞
Order of accuracy
Order of accuracy
kρ − ρex k2
3.0
2.5
2.0
3.5
3.0
2.5
2.0 10−2
h= 4.0 kρv − ρvex k1
p
10−2
(1/N )
h= 4.0 kρE − ρEex k1
3.0
2.5
2.0
(1/N )
3.5
p
(1/N )
kρE − ρEex k∞
3.0
2.5
2.0 10−2
h=
Farshad Navah
p
kρE − ρEex k2
kρv − ρvex k∞
Order of accuracy
Order of accuracy
kρv − ρvex k2 3.5
kρu − ρuex k∞
p
10−2
(1/N )
h=
4
th
High-order CFD Workshop
June, 2016
14 / 32
McGill University Code verification NS solver verification via MMS
NS solver verification Order of accuracy - P3 5.0
5.0 kρ − ρex k1
kρu − ρuex k1
4.5
kρu − ρuex k2
kρ − ρex k∞
Order of accuracy
Order of accuracy
kρ − ρex k2
4.0
3.5
3.0
4.5
4.0
3.5
3.0 10−2
h= 5.0 kρv − ρvex k1
p
10−2
(1/N )
h= 5.0 kρE − ρEex k1
4.0
3.5
3.0
(1/N )
4.5
p
(1/N )
kρE − ρEex k∞
4.0
3.5
3.0 10−2
h=
Farshad Navah
p
kρE − ρEex k2
kρv − ρvex k∞
Order of accuracy
Order of accuracy
kρv − ρvex k2 4.5
kρu − ρuex k∞
p
10−2
(1/N )
h=
4
th
High-order CFD Workshop
June, 2016
15 / 32
McGill University Code verification NS solver verification via MMS
NS solver verification Order of accuracy - P4 6.0
6.0 kρ − ρex k1
kρu − ρuex k1
5.5
kρu − ρuex k2
kρ − ρex k∞
Order of accuracy
Order of accuracy
kρ − ρex k2
5.0
4.5
4.0
5.5
kρu − ρuex k∞
5.0
4.5
4.0 −2
h= 6.0 kρv − ρvex k1
−2
p10 (1/N )
h= 6.0 kρE − ρEex k1
5.5
kρE − ρEex k2
kρv − ρvex k∞
Order of accuracy
Order of accuracy
kρv − ρvex k2
5.0
4.5
4.0
5.5
kρE − ρEex k∞
5.0
4.5
4.0 −2
h=
Farshad Navah
p10 (1/N )
p10 (1/N )
−2
h=
4th High-order CFD Workshop
p10 (1/N )
June, 2016
16 / 32
McGill University Code verification NS solver verification via MMS
NS solver verification Order of accuracy - P5 7.0
7.0 kρ − ρex k1
kρu − ρuex k1
6.5
kρu − ρuex k2
kρ − ρex k∞
Order of accuracy
Order of accuracy
kρ − ρex k2
6.0
5.5
5.0
6.5
6.0
5.5
5.0 10−2
h= 7.0 kρv − ρvex k1
p
10−2
(1/N )
h= 7.0 kρE − ρEex k1
6.0
5.5
5.0
(1/N )
6.5
p
(1/N )
kρE − ρEex k∞
6.0
5.5
5.0 10−2
h=
Farshad Navah
p
kρE − ρEex k2
kρv − ρvex k∞
Order of accuracy
Order of accuracy
kρv − ρvex k2 6.5
kρu − ρuex k∞
p
10−2
(1/N )
h=
4
th
High-order CFD Workshop
June, 2016
17 / 32
McGill University Code verification NS solver verification via MMS
NS solver verification qx → 2.0 × qx
Order of accuracy - P3 1.4
2.5 kρ − ρex k1
kρu − ρuex k1
kρ − ρex k2
2.0
kρ − ρex k∞
1.0
Order of accuracy
Order of accuracy
1.2
0.8 0.6 0.4 0.2
10−2
h= 1.4 kρv − ρvex k1
1.0 0.5
p
−0.5
h= 2.5 kρE − ρEex k1
kρv − ρvex k2
2.0
kρv − ρvex k∞
0.8 0.6 0.4
0.0
Farshad Navah
kρE − ρEex k2
kρE − ρEex k∞
0.5
−0.5
p
(1/N )
10−2
h=
4
th
(1/N )
1.0
0.0
h=
p
1.5
0.2
10−2
10−2
(1/N )
Order of accuracy
Order of accuracy
1.0
kρu − ρuex k∞
0.0
0.0 −0.2
1.2
kρu − ρuex k2
1.5
High-order CFD Workshop
p
(1/N )
June, 2016
18 / 32
McGill University Code verification Example of solution verification
Example of solution verification
Farshad Navah
4th High-order CFD Workshop
June, 2016
19 / 32
McGill University Code verification Example of solution verification
Turbulent Boundary Layer from TMR 2D zero-pressure-gradient flat plate with Re = 5 × 106 , M a = 0.2, χ∞ = 0.3 and χw = 0:
Figure: Domain and boundary conditions description
Discretization: 5 levels of h refinement 3 levels of p refinement: P 1, P 2 and P 3 Farshad Navah
4th High-order CFD Workshop
June, 2016
20 / 32
McGill University Code verification Example of solution verification
−3 3.00 ×10
10−4
2.95
10−5
|Cd − Cdexact|
2.90
10−6
P1 P2 P3
10−8 10−9
Cd
2.85
10−7
2.80 2.75
O(h2 ) O(h3 )
10−10
2.70
O(h4 ) 10
−11
P1 P2 P3 CFL3D FUN3D
2.65 10−3
h=
p
10−2
(1/N )
h=
(a) Estimated Cd Error
Farshad Navah
10−3
4th High-order CFD Workshop
p
10−2
(1/N )
(b) Cd
June, 2016
21 / 32
McGill University BS1 - DNS of the Taylor-Green Vortex
BS1 - DNS of Taylor-Green Vortex
Farshad Navah
4th High-order CFD Workshop
June, 2016
22 / 32
McGill University BS1 - DNS of the Taylor-Green Vortex
Code optimization CPR on Tensor-products
−→ Very sparse D and L operators
BR2 on Tensor-products
−→ Interior Penalty.
TGV for P 3 − 64 is 5 times cheaper after optimization
Farshad Navah
4th High-order CFD Workshop
June, 2016
23 / 32
McGill University BS1 - DNS of the Taylor-Green Vortex
12 simulations Resolution: 643 , 1283 , 2563
(based on dofs)
Polynomial: P 3, P 4, P 5, P 9
Farshad Navah
4th High-order CFD Workshop
June, 2016
24 / 32
McGill University BS1 - DNS of the Taylor-Green Vortex
Kinetic Energy, Ek , vs t∗ 0.14
0.14 P3 – 64 P3 – 128 P3 – 256 Spectral DNS
0.13 0.12
0.12 0.11
Ek
Ek
0.11 0.10
0.10
0.09
0.09
0.08
0.08
0.07
0.07
0.06
0.06 0
2
4
t∗
6
8
10
0
0.14
2
4
t∗
6
8
10
0.14 P5 – 66 P5 – 126 P5 – 258 Spectral DNS
0.13 0.12
P9 – 60 P9 – 130 P9 – 260 Spectral DNS
0.13 0.12 0.11
Ek
0.11
Ek
P4 – 65 P4 – 130 P4 – 255 Spectral DNS
0.13
0.10
0.10
0.09
0.09
0.08
0.08
0.07
0.07
0.06
0.06 0
Farshad Navah
2
4
t∗
6
8
10
0
4th High-order CFD Workshop
2
4
t∗
6
8
10
June, 2016
25 / 32
McGill University BS1 - DNS of the Taylor-Green Vortex
Kinetic energy dissipation, −∂Ek /∂t, vs t∗ 0.014
0.014 P3 – 64 P3 – 128 P3 – 256 Spectral DNS
0.012
−∂Ek /∂t∗
0.010
−∂Ek /∂t∗
0.010
P4 – 65 P4 – 130 P4 – 255 Spectral DNS
0.012
0.008
0.008
0.006
0.006
0.004
0.004
0.002
0.002
0.000
0.000 0
2
4
t∗
6
8
10
0.014
0
2
4
t∗
6
8
10
6
8
10
0.014 P5 – 66 P5 – 126 P5 – 258 Spectral DNS
0.012
−∂Ek /∂t∗
0.010
−∂Ek /∂t∗
0.010
P9 – 60 P9 – 130 P9 – 260 Spectral DNS
0.012
0.008
0.008
0.006
0.006
0.004
0.004
0.002
0.002
0.000
0.000 0
Farshad Navah
2
4
t∗
6
8
10
0
2
4th High-order CFD Workshop
4
t∗
June, 2016
26 / 32
McGill University BS1 - DNS of the Taylor-Green Vortex
Enstrophy, , vs time, t∗ 12
12 P3 – 64 P3 – 128 P3 – 256 Spectral DNS
8
8
6
6
4
4
2
2
0
0 0
2
4
t∗
6
8
10
0
12
2
4
t∗
6
8
10
6
8
10
12 P5 – 66 P5 – 126 P5 – 258 Spectral DNS
10
P9 – 60 P9 – 130 P9 – 260 Spectral DNS
10
8
8
6
6
ε
ε
P4 – 65 P4 – 130 P4 – 255 Spectral DNS
10
ε
ε
10
4
4
2
2
0
0 0
Farshad Navah
2
4
t∗
6
8
10
0
2
4th High-order CFD Workshop
4
t∗
June, 2016
27 / 32
McGill University BS1 - DNS of the Taylor-Green Vortex
Vorticity isocontours at x/L0 = −π and t∗ = 8 663
P = P5 1263
2583
Farshad Navah
4th High-order CFD Workshop
June, 2016
28 / 32
McGill University BS1 - DNS of the Taylor-Green Vortex
Vorticity isocontours at x/L0 = −π and t∗ = 8
Farshad Navah
Res = 2563
P3
P4
P5
P9
4th High-order CFD Workshop
June, 2016
29 / 32
McGill University BS1 - DNS of the Taylor-Green Vortex
Energy spectrum
10−1 P3 – 64 P3 – 128
10−2
Ek
10−3 10−4 10−5 10−6 10−7 100
101
k
Farshad Navah
4th High-order CFD Workshop
June, 2016
30 / 32
McGill University BS1 - DNS of the Taylor-Green Vortex
Work units vs DOFs
Work Units
107
P3 P4 P5 P9
106
105
104 105
106
107 DOFs
Farshad Navah
4th High-order CFD Workshop
June, 2016
31 / 32
McGill University BS1 - DNS of the Taylor-Green Vortex
Thank you for your attention! Questions?
Farshad Navah
4th High-order CFD Workshop
June, 2016
32 / 32
