Multigrid Methods for Solving the Convection-Diffusion Equations

1. Multigrid Methods for Solving the Convection-Diffusion Equations Chin-Tien Wu National Center for Theoretical Sciences Mathematics Division Tsing-Hua University

2. Motivation 1. Incompressible Navier-Stokes equation models many flows in practices. Using linearized backward Euler scheme, one has to solve an Oseen problem at each time levelUsing linearized backward Euler scheme, one has to solve an Oseen problem at each time level

4. Convection-Diffusion Equation

5. What are the foundations for building an efficient and accurate solver? Stable discretization methods for FEM. A good mesh. Reliable error estimation. Fast iterative linear solver.

6. I. Discretization Galerkin FEM Streamline diffusion upwinding FEM by Hughes and Brooks 1979 Discontinuous Galerkin FEM by Reed and Hill 1973 Edge-averaged FEM by Xu and Zikatanov 1999 (the linear system form this discretization is an M-matrix) In this thesis, we only consider GK and SDFEMIn this thesis, we only consider GK and SDFEM

7. Epson=10-3Epson=10-3

8. II. Methods for producing a good mesh? Generally, isotropic mesh is desired (DT). For problems with sharp gradient in some regions, anisotropic mesh, such as Shishkin mesh, is prefered. However, it is hard to know where to construct Shishkin-type of mesh without a priori knowledge of the exact solution especially when the domain is complicated. (Advancing front is able to produce long-thin element along boundary through proper chosen control function.) Generally, isotropic mesh is desired (DT). For problems with sharp gradient in some regions, anisotropic mesh, such as Shishkin mesh, is prefered. However, it is hard to know where to construct Shishkin-type of mesh without a priori knowledge of the exact solution especially when the domain is complicated. (Advancing front is able to produce long-thin element along boundary through proper chosen control function.)

9. III. A posteriori error indicators are used to pinpoint where the errors are large. A posteriori error estimation based on residual is proposed by Babuška and Rheinboldt 1978 A posteriori error estimation based on solving a local problem is proposed by Bank and Weiser 1985 A posteriori error estimation based on recovered gradient is proposed by Zienkiewicz and Zhu 1987 A posteriori error estimations for the convection-diffusion problems are proposed by Verfürth 1998, Kay and Silvester 2001.

10. Adaptive Solution Strategies

11. Improve Solution Accuracy by Mesh Refinement

12. IV. Solving Linear System Ax=b by Iterative Methods Methods: Stationary Methods: Jacobi, Gauss Seidel (GS), SOR. Krylov Subspace Methods: GMRES, by Saad and Schultz 1986, and MINRES, by Paige and Saunders 1975. Multigrid Methods: Geometric multigrid (MG), by Fedorenko 1961, and algebraic multigrid (AMG), by Ruge and Stüben 1985. 1) The convergence of the stationary iterative methods requires that the matrix from discretization is an M-matrix. 2) (i) The matrix from SDFEM is only positive definite. (ii) Convergence of stationary methods for convection-diffusion problem is shown for h << epsilon. (iii) MG with operator dependent prolongation is proved by Reusken. (iv) In practice, MG with standard prolongation still converges but no proof has been shown even for Problem2. 3) AMG convergence also requires the matrix from discretization is an M-matrix. However, AMG seems to be a robust solver for a wide range of applications. 4) To achieve fast convergence of GMRES, the eigenvalues of the matrix should be cluster to 1 as mentioned in earlier work of Elman. => GMRES need to be preconditioned. 5) (i) What is the best solver for solving the convection-diffusion equation discretized by SDFEM (on uniform mesh or on adaptive mesh)? (ii) What is a proper stopping criterion for the iterative solver during the mesh refinement process?1) The convergence of the stationary iterative methods requires that the matrix from discretization is an M-matrix. 2) (i) The matrix from SDFEM is only positive definite. (ii) Convergence of stationary methods for convection-diffusion problem is shown for h << epsilon. (iii) MG with operator dependent prolongation is proved by Reusken. (iv) In practice, MG with standard prolongation still converges but no proof has been shown even for Problem2. 3) AMG convergence also requires the matrix from discretization is an M-matrix. However, AMG seems to be a robust solver for a wide range of applications. 4) To achieve fast convergence of GMRES, the eigenvalues of the matrix should be cluster to 1 as mentioned in earlier work of Elman. => GMRES need to be preconditioned. 5) (i) What is the best solver for solving the convection-diffusion equation discretized by SDFEM (on uniform mesh or on adaptive mesh)? (ii) What is a proper stopping criterion for the iterative solver during the mesh refinement process?

13. Why multigrid works? 1. Relaxation methods converge slowly but smooth the error quickly. Consider

14. Multigrid (I)

15. Multigrid (II)

16. Multigrid (GMG) and Algebraic Multigrid (AMG)

17. Multigrid Components

18. GMG Convergence In unstructure grids, it remains challenge to generate semi-coarsening grids and operator-dependent prolongations. In unstructure grids, it remains challenge to generate semi-coarsening grids and operator-dependent prolongations.

19. GMG Convergence for Problem 1 1) ||E|| -> 0 when Epsilon -> 0: suggest HGS and MG are good preconditioner for h>>sqrt(?). 2) Theorem 4.3.2: ||ASk||< ? *(1+3* ? /h^2 min{h/sqrt(?),1})||Sk-1|| shows that HGS is a good smoother. The reason is: |||A|||<h-1 (from upper bound of Bsd ) => ||A||< ?-1/2 h-1 => ? 1/2h < ||A||-1. Since ? < ? 1/2h , for ?>> sqrt(h), Theorem 4.3.2 implies ||AS||< const ||A||-1. => smoothing property holds. 3) Using Theorem 4.3.2, a priori estimation (theorem 2.2.6) and regularity estimation (lemma 2.0.4), mg convergence can be shown.1) ||E|| -> 0 when Epsilon -> 0: suggest HGS and MG are good preconditioner for h>>sqrt(?). 2) Theorem 4.3.2: ||ASk||< ? *(1+3* ? /h^2 min{h/sqrt(?),1})||Sk-1|| shows that HGS is a good smoother. The reason is: |||A|||<h-1 (from upper bound of Bsd ) => ||A||< ?-1/2 h-1 => ? 1/2h < ||A||-1. Since ? < ? 1/2h , for ?>> sqrt(h), Theorem 4.3.2 implies ||AS||< const ||A||-1. => smoothing property holds. 3) Using Theorem 4.3.2, a priori estimation (theorem 2.2.6) and regularity estimation (lemma 2.0.4), mg convergence can be shown.

20. AMG convergence

22. AMG works for Problem 1

24. AMG Coarsening Algorithm (I)

25. AMG Coarsening Algorithm (II)

26. AMG Coarsening

27. AMG Coarsening on Uniform Meshes

29. GMG and AMG as a Solver (I) Adaptive meshes are generated from an initial 8x8 grid followed by 4 refinement with threshold 0.1, 0.01, 0.001 in cases of epsilon = 0.1, 0.01 and 0.001 Poor efficiency of GMG is due to the present of exponential boundary layer. Adaptive meshes are generated from an initial 8x8 grid followed by 4 refinement with threshold 0.1, 0.01, 0.001 in cases of epsilon = 0.1, 0.01 and 0.001 Poor efficiency of GMG is due to the present of exponential boundary layer.

30. GMG and AMG as a Solver (II) Adaptive meshes are generated from an initial 8x8 grid followed by 4 refinement with threshold 0.1 for epsilon = 0.1, 0.01 and 0.001 AMG convergence degrade as epsilon become smaller. GMG also perform poorly on uniform meshes. However, unlike in problem 1, GMG convergence on adaptive meshes is much faster than on the uniform mesh. The reasons are not only more coarse grid points are used for interpolation but also only characteristics present in this example. On the other hand, AMG converges much slower than GMGAdaptive meshes are generated from an initial 8x8 grid followed by 4 refinement with threshold 0.1 for epsilon = 0.1, 0.01 and 0.001 AMG convergence degrade as epsilon become smaller. GMG also perform poorly on uniform meshes. However, unlike in problem 1, GMG convergence on adaptive meshes is much faster than on the uniform mesh. The reasons are not only more coarse grid points are used for interpolation but also only characteristics present in this example. On the other hand, AMG converges much slower than GMG

31. GMG and AMG as a Preconditioner of GMRES

32. Stopping Criterion for Iterative Solver u : the weak solution of the convection-diffusion equation, uh : the finite element solution uh,n : the approximate iterative solution of uh. Given a prescribed tolerance ?. If || u-uh|| < 0.5?, clearly, we only have to ask || uh,n - uh || < 0.5? for large enough n. From a posteriori estimation || u-uh || < c?, where c is some constant, it is natural to acquire the iterative solution satisfies a stopping tolerance such that C1. || uh,n - uh || < c?, For adaptive mesh refinement, it is also desirable that C2. ? ? ?n (?n is the error indicator computed from uh,n ) Question : What stopping criteria should be imposed for iterative solutions? What solver requires least iterative steps to satisfy the stopping criteria?

33. Basic Ideas on Deriving the Stopping Criteria

34. Stopping criterion for satisfying C1.

35. Stopping criterion for satisfying C2.

36. Problem 1: Characteristic and downstream layers (ii)

37. Problem 2: Flow with Closed Characteristics (ii)

38. Problem 1: Characteristic and downstream layers (i)

39. Problem 2: Flow with Closed Characteristics (i)

40. Recent Numerical Results of Multigrid Methods on Solving Convection-Diffusion Equations GMRES and BiCGSTAB accelerated W-cycle MG with alternative zebra line GS smoother and upwind prolongation achieves h-independent convergence on problems with closed characteristics in which upwind discretization is used (by Oosterlee and Washio 1998). BICGSTAB accelerated V- and F- cycle AMG with symmetric GS smoother shows very slightly h-dependent convergence for problems with closed characteristics in which upwind discretization is used (Trottenberg, Oosterlee and Schüller: Multigrid p.519 2001). GMRES accelerated V-cycle MG with line GS smoother and bilinear, upwind or matrix-dependent prolongation achieves h-independent convergence for the model problems in which SDFEM discretization is used (by Ramage 1999).

41. Thank You

42. Conclusion SDFEM discretization is more stable than Galerkin discretization. Our error-adapted mesh refinement is able to produce a good mesh for resolving the boundary layers. On adaptive meshes, MG is a robust solver only for problems with only characteristics and AMG is robust for problems with only exponential layers. Both MG and AMG are good preconditioner for GMRES. Fewer iterative steps are required for the MG solver to satisfy our stopping criteria ( in Theorem 5.1.6 and Theorem 5.2.6) than to satisfy the heuristic tolerance (residual less than 10-6). No such saving can be seen if GMRES is used. The total saving of computation works is significant (can be more than half of the total works with heuristic tolerance). Robust in the sense of faster convergence and less dependent with mesh size and epsilonRobust in the sense of faster convergence and less dependent with mesh size and epsilon

43. Furture works

44. In this talk, we will discuss:

45. Galerking and SDFEM

46. It is well known that the solution from SDFEM discretization does not suffer from severe oscillation as the solution from Galerkin discretization due to improved coercivity from stablization. In the following, we present two of our test problems to show this phenomenon. It is well known that the solution from SDFEM discretization does not suffer from severe oscillation as the solution from Galerkin discretization due to improved coercivity from stablization. In the following, we present two of our test problems to show this phenomenon.

47. A Posteriori Error Estimation Based on Residual Proposed by Verfürth Alpha is from scaling between energy norm and H1 norm in the proof. For the test problem 2, alpha=h/sqrt(epsilon) Local lower bound blow up in the order of h/epsilpon in problem 2. |||u-u_h|||<h^{1/2}|u|_1<{h^{1/2}}/{\epsilon^{1/2}}||f||_0Alpha is from scaling between energy norm and H1 norm in the proof. For the test problem 2, alpha=h/sqrt(epsilon) Local lower bound blow up in the order of h/epsilpon in problem 2. |||u-u_h|||<h^{1/2}|u|_1<{h^{1/2}}/{\epsilon^{1/2}}||f||_0

48. A Posteriori Error Estimation Base on Solving a Local Problem Proposed by Kay and Silvester QT is the quadratic edge bubble function, BT is the cubic interior bubble function. Again, local lower bound blow up in the order of h/epsilon. ||b.?u|| is sqrt(epsilon/h) better than ||?u|| in the global a priori upper bound. => the error indicator blow up in the order of sqrt(h/epsilon).QT is the quadratic edge bubble function, BT is the cubic interior bubble function. Again, local lower bound blow up in the order of h/epsilon. ||b.?u|| is sqrt(epsilon/h) better than ||?u|| in the global a priori upper bound. => the error indicator blow up in the order of sqrt(h/epsilon).

49. Comparison of VR and KS Error Indicators (i) In the figure, VR is scaled by a factor of 32=sqrt(1024)In the figure, VR is scaled by a factor of 32=sqrt(1024)

50. Comparison of VR and KS Error Indicators (ii) E_T blow up in a rate of O(Pe) and E_{Omega} blow up in a rate of O(sqrt(Pe)) when h>> epsilon || b . grad(u) || is sqrt(epsilon/h) better than || grad(u) ||E_T blow up in a rate of O(Pe) and E_{Omega} blow up in a rate of O(sqrt(Pe)) when h>> epsilon || b . grad(u) || is sqrt(epsilon/h) better than || grad(u) ||

51. Moving Mesh and Error-Adapted Mesh Refinement

52. Mesh Moving by Equidistribution

53. Moving Mesh (I)

54. Moving Mesh (II)

55. Error-Adapted Mesh Refinement (i) On moving mesh, 1) meshes are no longer nested. 2) only very small moving step is allowed on anisotropic mesh. 3) Relaxation parameter r in moving mesh algorithm is important. (ii) On error-adapted mesh refinement, 1) meshes are nested. 2) highly stretched elements can be generated during earlier refinement stage. 3) Carefully chosen sensitivity parameter is needed. (iii) In 1992. Rippa points out that long-thin triangles can be good for linear interpolation ( triangle should be long in the direction where 2nd derivative is small and thin in the direction where 2nd derivative is large). !!! The error indicator represents grad(e_h). The external force in our algorithm represent the changes of gradient(e_h).(i) On moving mesh, 1) meshes are no longer nested. 2) only very small moving step is allowed on anisotropic mesh. 3) Relaxation parameter r in moving mesh algorithm is important. (ii) On error-adapted mesh refinement, 1) meshes are nested. 2) highly stretched elements can be generated during earlier refinement stage. 3) Carefully chosen sensitivity parameter is needed. (iii) In 1992. Rippa points out that long-thin triangles can be good for linear interpolation ( triangle should be long in the direction where 2nd derivative is small and thin in the direction where 2nd derivative is large). !!! The error indicator represents grad(e_h). The external force in our algorithm represent the changes of gradient(e_h).

56. Choice of error sensitivity parameter ?

57. Error-Adapted Mesh Refinement v.s. Regular Mesh Refinement (I)Problem 2 with e=10-4 Total 16 refinement steps are performed with threshold value 0.5, alpha=1/3 Initial mesh 4x4Total 16 refinement steps are performed with threshold value 0.5, alpha=1/3 Initial mesh 4x4

58. Variant IAHR/CEGB problem

59. Error-Adapted Mesh Refinement v.s. Regular Mesh Refinement (II)variant IAHR/CEGB problem with e=10-4 Initial mesh 4x4, error sensitivity parameter: 1/3Initial mesh 4x4, error sensitivity parameter: 1/3

60. Driven Cavity Flow for Re=100

61. GMRES

62. FEM Discretization