Stability and Algorithmic Efficiency of a High Order Finite Volume Method for Compressible Flows

1. Stability and Algorithmic Efficiency of a High Order Finite Volume Method for Compressible Flows Jean-Marie Le Gouez, Onera CFD Department ECCOMAS CFD 2016

2. Overview of the presentation • Context and objectives of the project • Main features of the H.O. compressible solver numerical scheme : reconstruction of conservative variables fields, flux models • Analysis of the accuracy and stability of the scheme on arbitrary grids • Spectrum of linear convection and diffusion discrete operators • Error associated to low order flux integration • Programming on hybrid Hardware architectures : CPU, GPU • Validation in terms of accuracy and algorithmic efficiency : H.O. CFD workshops • Unsteady laminar cases • Direct Numerical Simulation of the TAYLOR-GREEN vortex • Extension to Finite Volumes with curved faces

3. Context and objectives of the project • A prototype solver of the Onera CFD department • High order Reconstructed Finite Volume method for unstructured grids with mixed type elements, non compact scheme • Unlimited projection of the reconstructed polynomial on the interfaces : leads to types of fluxes schemes Characteristic upwind (state upwind, Roe-like), • Centered scheme, stabilized by a difference between left and right high-order extrapolated C.V. • Loworder or high order flux integration • Inviscid or laminar, (RANS and URANS), ALE • Participation to the H.O. (compressible) CFD solvers workshops : Ringleb flow, flapping Naca12 in the wake of a cylinder (Rey=1000), Heaving and pitching Naca12 (Rey=1000, 5000), convection of an isentropic vortex, TAYLOR-GREEN vortex : Accuracy, algorithmic and parallel efficiency • Also computing on GPU clusters, programmed in CUDA

4. Basics of the NXO scheme on arbitrary unstructured grids • Euler, Navier-Stokes for perfect gas law of state • 1 dof per cell and per equation (Volume average) • Non compact scheme : up to 20 cells to compute one flux (2D), up to 70 in 3D (4th to 5th order) • Parallel programming on node in shared memory using • Open-MP, loop-based on CPUs, or Cuda on GPUs • MPI & coarse partitioning at the internode level • (big partitions up to 10th of millions cells  Many-core, GPU) • Polynomial Reconstruction algorithm for the conservative variables fields • Preprocessor phase : Weighted Least-Square polynomial fit adapts • to the “quality” of the stencils • Gives the interpolation coefficients of conservative variable fields • from volume averages to surface averages • Different stencils and polynomial degrees for the convective and the • diffusive operators • Characteristic Upwind-biased convective scheme or centered scheme, • Explicit RK time integration or Dual Time-stepping , Target interface ,

5. L C Z X,Y l r R H.O. Full 3D Volume to face interpolation : Reconstruction and projection Reconstruction error functional ns = Stencil size : nb of monomials + 50% 1/ Reconstruction in a stencil : L, R or centred : Volume moment of order ijk

6. L C Z X,Y l r R FV NXO method : Reconstruction and projection 2/ Projection on the interface

7. : Volume average : NXO scheme : surface average 2 options for the inviscid fluxes Characteristic Upwind or centred, low order integration Higher order integration only used for high order geometry FV NXO method : Inviscid fluxes options , Upwind scheme : one average flux evaluation from the left and right extrapolated average conservative variables, characteristic splitting ‘state upwind’ Centred scheme : interpolation of the cell-average flux density tensor in all cells of the stencil to the interface + a stabilization term , Main inacuracy sources : asymptotically 2nd order Upwind scheme Centred scheme

8. Eigenvalues of the linear convective operator, depending on weights distribution and meshing method Uniform velocity orientation pi/8 , Periodic square, different meshing methods (gmsh) Unstable : too low weight decay with distance from the stencil center , Frontal Automatic Delaunay

9. Eigenvalues and Eigenvectors if the discrete linear convective operator, depending on the reconstruction degree , Typical stencils Mapping of one numerical eigenvector with higher Wave number (real part) onto the mesh Optimal weights distribution for the WLSQ, for each reconstruction degree k1  k4 Higher concentration of eigenvalues near the complex axis for K4 Enhanced damping of the higher wave number eigenmodes Delaunay mesh, 4th degree reconstruction Real part of the Eigenvector nearest Eigenvalue

10. Eigenvalues and Eigenvectors if the discrete diffusion operator (Laplacian) weight decay rate = 0.25 weight decay rate = 0.65 Insufficient weight decay rate as a function of the topological distance from the stencil center leads to an unstable scheme (positive real part for some eigenvalues) also for the diffusive operator

11. Eigenvalues and Eigenvectors if the discrete diffusion operator (Laplacian) , Delaunay mesh, 4th degree reconstruction Example of eigenvector of the periodic Laplacian, near to cos6x sin5y For all 4000 eigenvalues, nearest –(n2+m2), offset from it For reconstructions k1(purple)  k4 (red) : highest wave numbers Cell size H = 2pi/50

12. Non linear effects in flux integration Total error over the grid due to the flux evaluation from single projected values (interface integral of the reconstructed polynomials for conservative variables) x-axis : position along an interface y-axis : flux distributions for each equation from the reconstructed polynomials k1k5 Field of y-momentum Test case : Initial field for the isentropic vortex transport on the interfaces Upwind fluxes Comparison between :  Analytical expressions for the conservative variables and the normal flux densities  High order polynomial reconstructions in the cells and along the interfaces from exact cell averages

13. Non linear effects in flux integration Total error over the grid due to the flux evaluation from one single projected value (interface integral of the reconstructed polynomial) Mesh convergence index for reconstruction k5 : Optimal order for the fluxes of continuity equation (linear in the conservative variables), loss for other equations Gain with respect to the second order process of replacing the space flux integrand by its mean value (dotted lines for the analytical evaluations : exact averages, solid for the NXO method) Total error over the cell interfaces of the fine grid Upwind scheme Convergence with the reconstruction degree

14. Addressng the algorithmic efficiency barrier : accuracy versus number of dofsUnsteady laminar and Euler vortical flows, turbulent RANS ALE Spatially High-Order Finite Volume method for RANS / LES, coarse partitoning, wide halos , ,

15. Participation to the HO Workshop Workshop 1 Laminar / DNS Taylor-Green Vortex computed on a regular grid of tetraedra • Workshop 2 • Unsteady Laminar : dual local time stepping Heaving Naca12 in the unsteady wake of a cylinder, computation of the frequency lock-in of the vortex shedding Euler Isentropic vortex transport (right figure), results improved recently by using the flux reconstruction method • Workshop 3 • Unsteay Laminar Heaving and pitching naca12 at Reynolds numbers 1000 and 5000 demonstration of the grid convergence for both cases, moderate computation costs Euler, High order geometry Ringleb Flow, with a novel formulation of the method on Finite Volumes of H.O. geometry. Convergence diffiulties on the finest grids, but good results in the coarser ones (even with low count of dof, accurate CR-BC on curved walls) • Workshop 3 • Unsteay Laminar, BL3 Heaving and pitching Naca12, 3 different motions of the wing Grid / time step / polynomial order convergence (to be confirmed), moderate computation costs Grid and order convergence for the overset grid caseComparison in the conditions of the HO workshop : 50 tc , ,

16. High Order CFD Workshop Case 1.6 Vortex transport by uniform flow , , K convergence 2nd High Order CFD Workshop Köln May 2013

17. Case 1.6 Vortex transport by uniform flow Decomposition of the reconstructed solution on the modes 7 Mode 7 , Mode 6 Mode 4 Mode 5 Mode 3 Mode 2 , Mode 1 Mode 0 2nd High Order CFD Workshop Köln May 2013

18. NXO : FV on high order grids for the Euler equations RINGLEB Flow 1/ Characteristic upwind fluxes : “upwinded state” computed with the Jacobian eigenvectors 2/ The WLSQ polynolial reconstruction was extended to high order geometry Evaluation of the high order geometric moments of cells and interface (integrals of each monomial) by curvilinear integrals 3/ Higher order integration of the normal fluxes along the interface (Riemann solver on multiple points of the interface) Multiple evaluation of the projection of the same polynomial for each conservative variable on each sub-face element (averages on curvilinear sub-interfaces) • Parallel programming on node in shared memory using Open-MP (64 threads on SGI UV) • MPI & coarser partitioning only for inter-node communication Polynomial degrees : k4 to k5 inside the domain, k3 for slip wall b.c. stencils, k4 for non reflective far field bc 3rd High Order CFD Workshop Kissimmee january 2015

19. High Order CFD Workshop Case 1.1 Transonic Ringleb flow • Reconstructed field of Mach number • Exact and converged solutions plotted on each graph • Left : coarse grid 16*48 • Right medium grid 32*96 , , 3rd High Order CFD Workshop Kissimmee january 2015

20. High Order CFD Workshop Case 1.1 Transonic Ringleb flowHigh order flux integration Computations on 3 grids with p4 geometry: 16*48, 32*96, 64*192 • Effect of the high order flux integration : error divided by 11 to 30 for different fields Results on the finest grid, Nflx number of flux evaluation points along a curvilinear interface , , 3rd High Order CFD Workshop Kissimmee january 2015

21. High Order CFD Workshop Case 1.1 Transonic Ringleb flow , , • Spectral content of the reconstructed Y_momentum on grid2 (32*96) Left to right : mode 0, mode1, mode 2, mode 3 and full reconstructed solution 3rd High Order CFD Workshop Kissimmee january 2015

22. High Order CFD Workshop Case 3.4 2D Laminar Flapping wing Description of meshes used for the case. Generated by gmsh : 10500, 20000, 39000 triangles, external boundary at 32 chords ALE fully conservative formulation Non-linear implicit by dual time-stepping was used for this test-case, with RK 4 stages, CFL 3 for the dual local time stepping iteration. Time differencing scheme : 3rd order finite differences , , The computation is done in the absolute reference frame with fluid at rest at infinity, the free stream velocity is transferred in the grid motion and in the moving wall boundary condition 2nd High Order CFD Workshop Köln May 2013

23. High Order CFD Workshop Case 3.4 2D Laminar Flapping wingoccurrence of frequency lock-in , , Lift on cylinder and spectral analysis 2nd High Order CFD Workshop Köln May 2013

24. NXO : Laminar Heaving and Flapping wing Case setup • Computation done in the absolute reference frame with flow at rest at infinity • Translating and deforming grid • Grids by gmsh, isotropic triangles, shrunk in the normal direction to the profile by a factor 3 near the wall • External boundary at 6 chords 3rd High Order CFD Workshop Kissimmee january 2015

25. NXO : case setup • Results for Reynolds = 1000 • Results for Reynolds = 5000 3rd High Order CFD Workshop Kissimmee january 2015

26. High Order CFD Workshop Case 2.3 Heaving and Pitching profile , , Reynolds 5000 3rd High Order CFD Workshop Kissimmee january 2015

27. HO Workshop 4 : Case BL3 Energy extracting Grid and order convergence for the overset grid caseComparison in the conditions of the HO workshop : 50 tc , ,

28. High Order CFD Workshop Case 3.5 Taylor-Green Vortex Comparison of time derivatives of enstrophy and scaled Dissipation of kinetic energy Occurrence of an acoustic phenomenon , , High Order CFD Workshop Nashville Jan. 2012

29. GPU implementation of the NextFlow solver • Performance on each K20Xm GPU : • in k3 1,8e-8 s per RHS, 0.36s for 20 000 000 cells • in k4 2,5e-8 s per RHS, 0.50s for 20 000 000 cells • Taylor Green vortex 256**3 - wall-clock = 12 hours on 16 IVY-Bridge processors (total 128 cores) : 1600 hours CPU Intel core • 25 minutes on 16 Tesla K20M GPU • By comparison, at the 1st HO CFD workshop , this case requestedbetween 1100 and 33000 Intel coreCpuhours, depending on the numericalmethod • Taylor Green vortex 512**3 - wall-clock : 4 hours on 16 Tesla K20M GPUs Taylor-Green Vortex Rey = 1600 Computations on wedges