Exploring Alternative Approaches to CFD

1. Exploring Alternative Approaches to CFD Dimitri Mavriplis Dept.of Mechanical Engineering University of Wyoming Laramie, WY

2. Workshop Statement • Barriers to Progressand Technology Drivers • Misguided Government and Industry Policies • NASA: No base R&D budget, mission/vehicle emphasis • Boeing: Lack of interest in technology/internal expertise • Dell Computer Corp. • Multiple Scales • Anisotropy, Solver Stiffness, Mesh Refinement, Turbulence • Software Complexity • MDO, multiphase flow front tracking, education

3. Technology Drivers • High-Order Accuracy (Discontinuous Galerkin) • Attempt to mitigate multiple-scale issue • Large gains possible through more efficient approximations (temporal and spatial) • Requires efficient solvers to realize gains • More complex software • Low-Order Accuracy: Simpler Models (Lattice Boltzmann) • Simpler software • More complex physics: Mesoscales • Also requires efficient solvers

4. CFD Accuracy for Drag Prediction (AIAA workshop results) • Current simulations suffer from : • Discretization errors • Physical Modeling errors • Apparently, one does not dominate the other (still)

5. Apparently, one does not dominate the other (still) • Grid convergence (or lack thereof) from: AIAA-2003-3400: CFD Sensitivity Analysis of a Drag Prediction Wing/Body Transport Configuration, E.M. Lee-Rausch, P. G. Buning, J. H Morrison, M. A. Park, C. L. Rumsey and D. J Mavriplis

6. Higher-Order Methods • Simple asymptotic arguments indicate benefit of higher-order discretizations • Most beneficial for: • High accuracy requirements • Smooth functions

7. Motivation • Higher-order methods successes • Acoustics • Large Eddy Simulation (structured grids) • Other areas • High-order methods not demonstrated in: • Aerodynamics, Hydrodynamics • Unstructured mesh LES • Industrial CFD • Cost effectiveness not demonstrated: • Cost of discretization • Efficient solution of complex discrete equations

8. Motivation • Discretizations well developed • Spectral Methods, Spectral Elements • Streamwise Upwind Petrov Galerkin (SUPG) • Discontinuous Galerkin • Most implementations employ explicit or semi-implicit time stepping • e.g. Multi-Stage Runge Kutta ( ) • Need efficient solvers for: • Steady-State Problems • Time-Implicit Problems ( )

9. General Idea • Use element Jacobi as Smoother for p-multigrid solver (AIAA-2003-3989: Atkins, Helenbrook and Mavriplis) Solution of linear (2D) wave equation using 2 grid p-MG system for p=3 (Mavriplis and Atkins, 2002)

10. Multigrid Solver for Euler Equations • Develop efficient solvers (O(N)) for steady-state and time-implicit high-order spatial discretizations • Discontinuous Galerkin • Well suited for hyperbolic problems • Compact-element-based stencil • Use of Riemann solver at inter-element boundaries • Reduces to 1st order finite-volume at p=0 • Natural extension of FV unstructured mesh techniques • Closely related to spectral element methods

11. Discontinuous Galerkin (DG) Mass Matrix Spatial (convective or Stiffness) Matrix Element Based-Matrix Element-Boundary (Edge) Matrix

12. Steady-State Solver • Kijui=0 (Ignore Mass matrix) • Block form of Kij: • Eij = Block Diagonals (coupling of all modes within an element) • Fij = 3 Block Off-Diagonals (coupling between neighboring elements) Solve iteratively as: Eij (ui n+1 – uin ) = Kij uin

13. Steady-State Solver: Element Jacobi Solve iteratively as: Eij (ui n+1– uin) = Kij uin Duin+1 = E-1ij Kij uin Obtain E-1ij by Gaussian Elimination (LU Decomposition) 10X10 for p=3 on triangles

14. DG for Euler Equations • Mach = 0.5 over 10% sin bump • Cubic basis functions (p=3), 4406 elements

15. Entropy as Measure of Error • S  0.0 for exact solution • S is smaller for higher order accuracy

16. Error Convergence for Bump Case • P=1: Final Slope: 1.26 • P=2: Final Slope: 2.8 • P=3: Final Slope: 3.2 • Based on L2 norm of entropy • H-convergence • 382 elements • 1500 elements • 2800 elements • 4406 elements

17. Element Jacobi Convergence • P-Independent Convergence • H-dependence

18. Improving Convergence H-Dependence • Requires implicitness between grid elements • Multigrid methods based on use of coarser meshes for accelerating solution on fine mesh

19. Spectral Multigrid • Form coarse “grids” by reducing order of approximation on same grid • Simple implementation using hierarchical basis functions • When reach 1st order, agglomerate (h-coarsen) grid levels • Perform element Jacobi on each MG level

20. Hierarchical Basis Functions • Low order basis functions are subset of higher order basis functions • Low order expansion (linear in 2D): • U= a1F1 + a2F2 + a3F3 • Higher order (quadratic in 2D) • U=a1F1 + a2F2 + a3F3 + a4F4 + a5F5 + a6F6 • To project high order solution onto low order space: • Set a4=0, a5=0, a6=0

21. Hierarchical Basis on Triangles • Linear (p=1): F1=z1, F2=z2, F3=z3 • Quadratic (p=2): • F4=-0.5z1z2, F5=-0.5z2z3, F6=z3z1 • Cubic (p=3): F7=1/12z1z2(z1-z2), F8=1/12z2z3(z2-z3), F9=1/12z3z1(z3-z1), F10=z1z2z3

22. Spectral Multigrid • Fine/Coarse Grids contain same elements • Transfer operators almost trivial for hierarchical basis functions • Restriction: Fine to Coarse • Transfer low order (resolvable) modes to coarse level exactly • Omit higher order modes • Prolongation: Coarse to Fine • Transfer low order modes exactly • Zero out higher order modes

23. Euler Equations • System of Non-Linear Equations • Element Jacobi and Spectral MG Implemented for Euler Equations • FAS Multigrid (non-linear MG) • Possible tradeoffs using Linear MG • Requires storing entire Jacobian • Reduced cpu time (freeze Jacobian) • Bump and Airfoil Test Cases • Includes Curved Surface Elements • Necessary for maintaining accuracy

24. Element Jacobi Convergence • P-Independent Convergence • H-dependence

25. Multigrid Convergence • Nearly h-independent

26. Future Multigrid Work • Reduced complexity of element Jacobi • Develop implicit time integration scheme • High-order time discretization • Multigrid driven implicit solver • H-P Adaptivity in Multigrid Context • Towards a General Multi-resolution Framework

27. Using Alternate Models: The Lattice Boltzmann Approach • Lattice Boltzmann methods increasingly popular • Origins from Lattice Gas Automata (LGA) methods • LBM abandons particle formulation in favor of particle distribution function approach • Mesoscale Method (neither macro, nor molecular) • Shown to converge to Navier-Stokes equations in macroscopic limit • Facilitates physical modeling (and software complexity) • Front capturing in place of fitting • Physics incorporated at mesoscale level (collision term)

28. Front Capturing with LBM • c/o M. Krafczyk, TU Braunschweig, Germany

29. Governing Equations of Fluid Dynamics • Statistical Mechanics Description (Meso-Scale) • Boltzmann Equation (BGK approximation) • f=f(x, ,t): Distribution function in velocity space • = speed of gas molecules • Linear convection terms • Non-Linear Collision (source) term (local)

30. Lattice Boltzmann • 9 velocity model: • Velocities allow exact jump to neighboring grid point • Proven mathematically to reduce to Navier-Stokes equations in asymptotic limit • Obtain macroscopic values as:

31. 2 Step LBM Implementation • Implemented as Collision step followed by advection step (reminiscent of LGA method) • Collision step is entirely local : is a non-linear function of the • Advection: Shift in memory with no flops!

32. Motivation • Formally: 2nd order accurate in space • Formally: 1st order accurate in time • Remarkable delivered accuracy in many cases • Numerically efficient • Local collision term • Linear convection term: shift operator with no flops! • Corresponds to explicit time-stepping approach • Notoriously inefficient solution strategy

33. LBE for 2D Driven Cavity • Streamlines for 129 x 129 grid (Re=100)

34. Solution of Steady-State LBE • 16,000 time steps to machine zero on 33x33 grid

35. Solution of Steady-State LBE • 16,000 time steps to machine zero on 33x33 grid • 4 to 5 times slower on 129x129 grid

36. Solution of Steady-State LBE • 16,000 time steps to machine zero on 33x33 grid • 4 to 5 times slower on 129x129 grid • Even slower on 513x513 grid

37. Motivation • Develop efficient (optimal) solution strategies for LBE • LBE simplicity + Asymptotic efficiency • Steady-state problems • Decouple space and time • Preserve exact LBE spatial discretization • Obtain identical final result with faster solution method

38. Steady-State LBE Discretization • Drop time-level index n: • In residual form: • Linearize and solve with Newton’s method: • Iterative inversion of Jacobian matrix using Jacobi, GG, Multigrid…

39. Linear Multigrid Driven Newton Scheme • Efficiency of linear system • 150 multigrid cycles • Quadratic convergence of non-linear system

40. Linear vs. Non-Linear Approach • Linearized LBE is expensive • Block 9x9 matrices for collision term • Looses efficient 2 step non-linear evaluation

41. Operation Count • Linear system iteration (Jacobi): • 9x9 matrix-vector product: 2x9x9= 162 ops • Other terms : 63 ops • Total: 225 ops • Non-Linear LBE Time Step: • Macro variables: 18 ops • Collision: 75 ops • Advection: 0 ops • Total : 93 ops

42. Non-Linear Multigrid Approach • Use LBE time-step to drive non-linear multigrid algorithm (non-linear equivalent of Jacobi) • Jacobi has poor smoothing properties • Implement non-linear under-relaxed Jacobi in 3 steps:

43. Non-Linear Multigrid (FAS)

44. Non-Linear Multigrid Approach • On coarse grids, need to include defect-correction : • Implement as modification to under-relaxation stage • Preserve first two steps • Implemented by calling existing LBE code on all grid levels

45. Multigrid Cost • 1 MG cycle = 24 LBE steps • 4 pre-smoothing and 4 post-smoothing steps (factor 8) • Work of coarse grids (W-cycle in 2D: factor 2) • Remaining factor: 1.5: Under-relaxation, intergrid transfers • Memory Usage (factor 3) • 2 field arrays on fine grid (f + residuals) • Remainder is coarse grid variables (factor 1.5)

46. Solution of Steady-State LBE • 16,000 time steps to machine zero on 33x33 grid • 4 to 5 times slower on 129x129 grid • Even slower on 513x513 grid

47. Non-Linear Multigrid Solution of LBE • 2 to 3 orders of magnitude faster in terms of cycles • Multigrid cycle is 24 times cost of LBE time-step

48. Non-Linear Multigrid Solution of LBE • Work unit: CPU time for 1 LBE time-step on relevant grid • MG 1 to 2 orders of magnitude more efficient • Benefit increases with finer grids • Asymptotic property of MG: O(N)

49. Conclusions • Steady-state LBE can be solved efficiently O(N) using multigrid • Implemented by calling existing LBE code on each grid level • Plus defect correction and under-relaxation step • Better smoothers are available at the cost of decreased modularity • Linearizing LBE is a poor strategy • Destroys favorable operation count

50. Future Challenges • For steady-state cases, apply known PDE strategies • Adaptive meshing • Adjoint solution techniques for MDO sensitivities • Simple (but costly) linearization • Extend steady-state solver to time-implicit solver • Exact convection (shift to neighbor) is lost • Effect on temporal accuracy (?)