Augmented Lagrangian Preconditioner for Linear Stability Analysis of incompressible fluid flows on large configurations J.Moulin 1 , J-L. Pfister 1 , O.Marquet 1 , P. Jolivet 2 1 Office National d’Etudes et de Recherches Aérospatiales 2 ENSEEIHT – Institut de Recherche en Informatique de Toulouse Funded by ERC Starting Grant FreeFem++ Days, Paris, 14-15 december 2017
Introduction What is Linear Stability Analysis ? One wants to know if some steady solution of equation (1) is temporally stable or unstable : �� �� = � � (1) Method : Linear Stability Analysis � � � = � Step 1 : compute a steady solution � � � �� around the Step 2 : test its stability for small monochromatic perturbations � steady solution � � ��� � = �(� � )� � Growth rate Frequency Mass matrix � = ℜ[�] � = ℑ[�] Jacobian matrix : (spatial discretization) �� �� (� � ) 2
Introduction A classical example From [Goharzadeh & Molki, 2015] Typical question : What is the critical Reynolds number above which the von Karman vortex street appears ? 3
Introduction A classical example Step 1 : compute a steady solution Steady Navier-Stokes solution ( �� = 50 ) [Sipp et al, 2010] 4
Introduction A classical example Step 1 : compute a steady solution Steady Navier-Stokes solution ( �� = 50 ) [Sipp et al, 2010] � � � �� around the Step 2 : test its stability for small monochromatic perturbations � steady solution � � ��� � = �(� � )� � � = ℜ[�] � = ℑ[�] 5
Introduction A classical example Step 1 : compute a steady solution Steady Navier-Stokes solution ( �� = 50 ) [Sipp et al, 2010] � � � �� around the Step 2 : test its stability for small monochromatic perturbations � steady solution � � ��� � = �(� � )� � � = ℜ[�] � = ℑ[�] Unstable eigenmode at �� = 50 [Sipp et al, 2010] 6
Introduction Why Linear Stability Analysis ? Some nice features : o Easy to determine a threshold value (sign of ℜ[�] ) o Less expensive than nonlinear time-integration But some computational burdens : o Find a (not necessarily stable) steady solution : Newton method � Multiple inversions of � o Find internal eigenvalues of generalized EV problems : Krylov-Schur + shift-and-invert � Multiple inversions of J − � � , where � is the shift 7
Introduction How to invert matrix of the type − ! " efficiently ? � For reasonably small configurations : direct sparse solvers (MUMPS, SUPERLU, etc) � For large configurations : iterative method (GMRES, BiCGSTAB, …) + good preconditioner 8
Introduction Linearized incompressible Navier-Stokes operator ( i.e. J − � � ) : � # + (# � ⋅ &)# + (# ⋅ &)# � + &' − 1 �� &²# = * −& ⋅ # = + Complex shift ( � = 0 in Newton) Once discretized with FE : classical saddle-point problem / 0 - . , = - 0 ' 1 How to precondition this ? SIMPLE [Patankar 1980] o Stokes Preconditioner [Tuckerman, 1989] (based on adaptation of existing time-stepping code) o Pressure Convection Diffusion [Silvester et al. 2001] o Least-Squares Commutator [Elman et al. 2006] o Augmentated Lagrangian [Benzi and Olshanskii 2006], [Heister and Rapin 2013] o 9
Introduction Linearized incompressible Navier-Stokes operator ( i.e. J − � � ) : � # + (# � ⋅ &)# + (# ⋅ &)# � + &' − 1 �� &²# = * −& ⋅ # = + Complex shift ( � = 0 in Newton) Once discretized with FE : classical saddle-point problem / 0 - . , = - 0 ' 1 How to precondition this ? SIMPLE [Patankar 1980] o Stokes Preconditioner [Tuckerman, 1989] (based on adaptation of existing time-stepping code) o Pressure Convection Diffusion [Silvester et al. 2001] o Least-Squares Commutator [Elman et al. 2006] o Augmentated Lagrangian [Benzi and Olshanskii 2006], [Heister and Rapin 2013] o 10
Overview 1 – Augmentation-based preconditioners 2 – Performances 3 – FreeFem++ parallel implementation 4 – Parallel 3D numerical examples 5 – Some further refinement … 11
1- Augmentation-based preconditioners Augmented problems Augmented Lagrangian (algebraic augmentation) / 0 / 0 - . - . , 2 = , + 3- . 4 56 - , 2 , 2 = = ' 1 ' 1 - 0 - 0 2 = 0 + 3- . 4 56 1 0 12
1- Augmentation-based preconditioners Augmented problems Augmented Lagrangian (algebraic augmentation) / 0 / 0 - . - . , 2 = , + 3- . 4 56 - , 2 , 2 = = ' 1 ' 1 - 0 - 0 2 = 0 + 3- . 4 56 1 0 Grad-Div augmentation (variational augmentation) 8 + �� 56 9#: 9# 7 � # ⋅ # 8 + (# ⋅ &)# ⋅ # 8 − p9 ⋅ # 8 < + 7 3(9 ⋅ #)(9 ⋅ # 8) = � < − 7 9 ⋅ # � 8 = � < 13
1- Augmentation-based preconditioners Augmented problems Augmented Lagrangian (algebraic augmentation) / 0 / 0 - . - . , 2 = , + 3- . 4 56 - , 2 , 2 = = ' 1 ' 1 - 0 - 0 2 = 0 + 3- . 4 56 1 0 Grad-Div augmentation (variational augmentation) 8 + �� 56 9#: 9# 7 � # ⋅ # 8 + (# ⋅ &)# ⋅ # 8 − p9 ⋅ # 8 < + 7 3(9 ⋅ #)(9 ⋅ # 8) = � < − 7 9 ⋅ # � 8 = � < Augmented Lagrangian leaves the discrete solution unchanged Grad-Div leaves the continuous solution unchanged 14
1- Augmentation-based preconditioners Augmented problems Augmented Lagrangian (algebraic augmentation) / 0 / 0 - . - . , 2 = , + 3- . 4 56 - , 2 , 2 = = ' 1 ' 1 - 0 - 0 2 = 0 + 3- . 4 56 1 0 Grad-Div augmentation (variational augmentation) 8 + �� 56 9#: 9# 7 � # ⋅ # 8 + (# ⋅ &)# ⋅ # 8 − p9 ⋅ # 8 < + 7 3(9 ⋅ #)(9 ⋅ # 8) = � < − 7 9 ⋅ # � 8 = � < Augmented Lagrangian leaves the discrete solution unchanged Grad-Div leaves the continuous solution unchanged 15
1- Augmentation-based preconditioners Classical vs. modified version In both cases, the same block structure arises : E 0 - . 56 - . , 2 0 , 2 E , 2 = 56 B M D = −BA 2 56 -, 2 E 0 D - 0 0 E N B C Classical preconditioner - . = >?@AA = BC = , 2 0 D with D 56 ≃ Re 56 + 3 � I 56 − � -� J - . 56 56 ≃ it’s complicated … , 2 Main features : Mesh optimality � Reynolds optimality � 56 ouch !) The higher 3 , the less iterations ( , 2 � 16
1- Augmentation-based preconditioners Classical vs. modified version In both cases, the same block structure arises : E 0 - . 56 - . , 2 0 , 2 E , 2 = 56 B M D = −BA 2 56 -, 2 E 0 D - 0 0 E N B C Classical preconditioner Modified preconditioner , 66,2 , 6U,2 - . - . = >?@AA = BC = , 2 = OPQRS = 0 , UU,2 0 D 0 D with with D 56 ≃ Re 56 + 3 � I 56 − � -� J - . 56 D 56 ≃ Re 56 + 3 � I 56 − � -� J - . 56 56 ≃ off-the-shelf algebraic multigrid 56 ≃ it’s complicated … , RR,2 , 2 Main features : Main features : Mesh optimality � Mesh optimality � Reynolds optimality � Reynolds dependent � 56 ouch !) The higher 3 , the less iterations ( , 2 � Exists an optimal and case dependent 3 � 17
1- Augmentation-based preconditioners Classical vs. modified version In both cases, the same block structure arises : E 0 - . 56 - . , 2 0 , 2 E , 2 = 56 B M D = −BA 2 56 -, 2 E 0 D - 0 0 E N B C Classical preconditioner Modified preconditioner , 66,2 , 6U,2 - . - . = >?@AA = BC = , 2 = OPQRS = 0 , UU,2 0 D 0 D with with D 56 ≃ Re 56 + 3 � I 56 − � -� J - . 56 D 56 ≃ Re 56 + 3 � I 56 − � -� J - . 56 56 ≃ off-the-shelf algebraic multigrid 56 ≃ it’s complicated … , RR,2 , 2 Main features : Main features : Mesh optimality � Mesh optimality � Reynolds optimality � Reynolds dependent � 56 ouch !) The higher 3 , the less iterations ( , 2 � Exists an optimal and case dependent 3 � 18
2- Performances Choice of 3 The choice of a good 3 is determinant for the preconditioning efficiency ! Bright side : since the preconditioner is �� ↑ independent of the mesh Optimal 3 can be found on a � coarse mesh Dark side : Optimal 3 is problem and �� – dependent Figure : Influence of �� ∈ [10,120] on optimal 3 for modified Grad-Div preconditioner 19
2- Performances CPU time in Newton method Full MUMPS Modified Grad-Div Velocity Pressure Mesh Facto Reso tot/ndof Facto Reso tot/ndof DOFs DOFs ( Y! ) ( Y! ) ( Z! ) ( Y! ) ( Y! ) ( Z! ) 32x32 9900 1300 140 0 20 30 50 14 64x64 39000 5000 810 10 27 320 250 20 96x96 88000 11000 2250 40 33 840 580 21 256x256 623400 78200 34480 290 62 8090 4780 25 Averaged timings for 1 Newton iteration ( 2D lid-driven cavity, �� = 100 , 3 = 0,1 ) 20 *All sub-systems are solved with MUMPS
Recommend
More recommend
