Preconditioners for computations of subduction zone modelling - PowerPoint PPT Presentation

Preconditioners for computations of subduction zone modelling Sander Rhebergen, Richard Katz, Andrew Wathen FEniCS’13, University of Cambridge

Overview March, 2013 · Preconditioners for Stokes · PETSc · Wedge flow benchmark · Conclusions

Stokes equations March, 2013 Let Ω ∈ R d be the domain of interest and denote the boundary of Ω by ∂ Ω. The Stokes equations are then given by: −∇ 2 u + ∇ p = f , ∇ · u = 0 in Ω . with Dirichlet and Neumann boundary conditions: � u = w on Γ D , ∇ u − p I ) · n = s on Γ N , where ∂ Ω = Γ D ∪ Γ N and Γ D ∩ Γ N = ∅ . u ∈ R d is the velocity vector p ∈ R the pressure f ∈ R d a given vector body force

FEM for Stokes March, 2013 E 0 and M h ⊂ L 2 (Ω). The weak formulation: Let X h 0 ⊂ H 1 E and p h ∈ M h such that Find a u h ∈ X h ∀ v h ∈ X h 0 and ∀ q h ∈ M h , a ( u h , v h )+ b ( v h , p h )+ b ( u h , q h ) = l ( v h ) where � a ( u , v ) := ∇ u : ∇ v , Ω � b ( v , q ) := − q ∇ · v , Ω � � l ( v ) := v · f + s · v Ω Γ N

Discrete linear system of Stokes March, 2013 E and p h ∈ M h such that Find a u h ∈ X h ∀ v h ∈ X h a ( u h , v h ) + b ( v h , p h ) = l ( v h ) 0 ∀ q h ∈ M h b ( u h , q h ) =0 with ( u h , p h ) an approximation to ( u , p ). The discrete linear system has the form:       B T  A  u  f  =  , A U = b ⇐ ⇒  0 B p g in which A corresponds to the viscous term, B T to the gradient operator and B to minus the divergence operator.

Solving the discrete linear system March, 2013 Solving A U = b · Direct solvers · Iterative methods

Conjugate Gradient Theory March, 2013 Solve AU = b using a Conjugate Gradient method: One can show that � √ κ − 1 � k | e ( k ) | | | A κ = λ max ( A ) √ κ + 1 ≤ λ ∈ σ ( A ) | p k ( λ ) | ≤ 2 p k ∈ Π k ,p k (0)=1 max min , | e (0) | | | A λ min ( A ) For fast convergence you want: · a few distinct eigenvalues or clustered eigenvalues → it will be easier to find a good polynomial p k ( z ) · reduce κ → if the condition number is small, there will be rapid convergence of the CG method

Preconditioners March, 2013 Instead of solving A U = b solve P − 1 A U = P − 1 b where the preconditioning matrix P is such that P − 1 A has · a few distinct eigenvalues · clustered eigenvalues · κ ( P − 1 A ) is small

Preconditioners for Stokes March, 2013 Stokes:       B T  A  u  f  = A U = b ⇐ ⇒   0 B p g LDU decomposition of A :         A − 1 B T B T 0 0  A I  A  I  =  , A = LDU =    BA − 1 B 0 I 0 S 0 I in which S = − BA − 1 B T is the Schur-complement matrix.

Preconditioners for Stokes March, 2013 Examples of “impractical” preconditioners: · P = LDU (direct inverse → convergence in 1 iteration) · P = D (three distinct eigenvalues, MINRES converges in 3 iterations) · P = DU (two distinct eigenvalues, BiCGStab converges in 2 iterations) Required: the inverse of the Schur-complement S = − BA − 1 B T which is a full matrix. Can we approximate S such that we can obtain good “practical” preconditioners?

Preconditioners for Stokes March, 2013 Practical preconditioners:     0 0  A MG �  A � P = ≈ D = ,   0 diag( Q ) 0 S     B T B T  A MG �  A � P = ≈ DU = ,   0 diag( Q ) 0 S � where Q = Ω pq is the pressure mass-matrix. · diag( Q ) is spectrally equivalent to S · multigrid is optimal for A Combined: iterative methods for solving P − 1 A U = P − 1 b will be optimal, that is independent of the problem size!

Overview March, 2013 · Preconditioners for Stokes · PETSc · Wedge flow benchmark · Conclusions

PETSc March, 2013 · Large collection of iterative solvers and preconditioners. · Freely available at http://mcs.anl.gov/petsc

PETSc, P ≈ D March, 2013 -ksp type minres -pc type fieldsplit -pc fieldsplit detect saddle point -pc fieldsplit type schur -pc fieldsplit schur fact type diag -pc fieldsplit schur precondition user -fieldsplit 0 ksp type preonly -fieldsplit 0 pc type hypre -fieldsplit 1 ksp type preonly -fieldsplit 1 pc type jacobi

PETSc, P ≈ DU March, 2013 -ksp type bcgs -pc type fieldsplit -pc fieldsplit detect saddle point -pc fieldsplit type schur -pc fieldsplit schur fact type upper -pc fieldsplit schur precondition user -fieldsplit 0 ksp type preonly -fieldsplit 0 pc type hypre -fieldsplit 1 ksp type preonly -fieldsplit 1 pc type jacobi

FEniCS+PETSc results March, 2013 Stokes equations, lid-driven cavity flow, 2D. Grid P ≈ D P ≈ DU 16 2 44 17 32 2 42 16 64 2 42 16 128 2 40 18 256 2 40 18

Overview March, 2013 · Preconditioners for Stokes · PETSc · Wedge flow benchmark · Conclusions

Wedge flow benchmark March, 2013 660 km Rigid overriding plate (to 50 km) 600 km Fixed plate velocity 5 cm/yr (a) van Keken et al. 2008 (b) http://www.geosci.usyd.edu.au Cornerflow model for subduction zone dynamics (van Keken et al. 2008).

Velocity + pressure solution March, 2013 u = (0 , 0) ∇ · ǫ − ∇ p = 0 ∇ · u = 0 √ u = (1 , − 1) / 2 The deviatoric stress tensor ǫ = η ( ∇ u + ( ∇ u ) T ) with viscosity: � η 0 � E diff � − 1 � + η 0 η diff , effective = η diff ( T ) = A diff exp , . η diff η max RTT 0 Dimensionless viscosity: η ∈ [38 , 10 5 ]. Dimensional: multiply by 10 21 .

Boundary conditions Stokes March, 2013 u = (0 , 0) √ u = (1 , − 1) / 2 ( ǫ − p I ) · n = 0

Temperature solution March, 2013 u · T = Pe − 1 ∇ 2 T Here Pe ≈ 1310 is the Peclet number.

Boundary conditions Temperature March, 2013 T = 273 / 1573 T = f ( y ) T = g ( y ) T = 1 on inflow ∇ T · n = 0 on outflow ∇ T · n = 0

FEM weak formulation on Ω 1 March, 2013 E and p h ∈ M h such that Find a u h ∈ X h ∀ v h ∈ X h 0 , and ∀ q h ∈ M h a ( u h , v h ; T h )+ b ( v h , p h )+ b ( u h , q h ) = l ( v h ) where � a ( u , v ; T ) := 2 η ( T ) ǫ ( u ) : ǫ ( v ) , Ω � b ( v , q ) := − q ∇ · v , Ω � l ( v ) := v · f Ω

FEM weak formulation on Ω March, 2013 Find a T h ∈ S h E such that ∀ W h ∈ S h a t ( T h , W h ; u h ) = l t ( W h ) 0 where � � � 1 a t ( T, W ; u ) := ( u · ∇ T ) W + P e ∇ W · ∇ T , Ω � l t ( W ) := Wf t Ω

Solving iteratively March, 2013 We solve the two systems iteratively, i.e., given T (0) h , we iterate over k ≥ 1: a ( u ( k ) h , v h ; T ( k − 1) ) + b ( v h , p ( k ) h ) + b ( u ( k ) h , q h ) = l ( v h ) a t ( T ( k ) h , W h ; u ( k ) h ) = l t ( W h ) until some convergence criterium has been met. We consider the | T ( k ) − T ( k − 1) | 2 < δ , δ = 10 − 2 . system solved if | | h h

Solving FEM discretization of Stokes March, 2013 To solve Stokes, we use the same preconditioner as before:     0 0  A MG �  A � P = ≈ D =   0 diag( Q ) 0 S � Ω η ( T ) − 1 pq dx is the scaled pressure mass-matrix. where Q =

Solving FEM discretization of Stokes March, 2013 -ksp_type minres -pc_type fieldsplit -pc_fieldsplit_detect_saddle_point -pc_fieldsplit_type schur -pc_fieldsplit_schur_fact_type diag -pc_fieldsplit_schur_precondition user -fieldsplit_0_ksp_type preonly -fieldsplit_0_pc_type hypre -fieldsplit_1_ksp_type preonly -fieldsplit_1_pc_type jacobi

Solving FEM discretization of Temperature eq. March, 2013 -adv_pc_type asm -adv_sub_pc_type lu -adv_pc_asm_blocks 128 -adv_pc_asm_overlap 2

Plots of solution March, 2013 (c) Temperature field. (d) Velocity field.

Comparison with University of Michigan March, 2013 | | T slab | | | | T wedge | | Code (elements) T 60 OX (50,398) 571.52 602.24 1001.96 OX (200,836) 576.75 604.87 1002.27 OX (804,060) 578.08 605.62 1002.03 UM 580.66 607.11 1003.20 All 570.30-581.30 606.94-614.09 1002.15-1007.31 Temperature (in ◦ C ) comparison between the OX code, the UM code and the interval of results of all codes in van Keken et al. (2008). · T 60 is the temperature (in ◦ C ) at ( x, y ) = (60 km, − 60 km ). · | | T slab | | L2 norm of the slabwedge interface temperature · | | T wedge | | L2 norm of the temperature in the triangular part of the tip of the wedge

Optimal solver? March, 2013 Code (elements) Outer Its. Stokes Adv.Dif. OX (50,398) 16 127.94 (11s) 55.75 (2s) OX (200,836) 17 135.47 (46s) 83.35 (9s) OX (804,060) 18 138.00 (190s) 107.61 (50s) Optimal solver for Stokes part, but the advection diffusion equation for the temperature is dependent on the grid size.

Overview March, 2013 · Preconditioners for Stokes · PETSc · Wedge flow benchmark · Conclusions

Conclusions March, 2013 · Optimal preconditioners for Stokes · Implementation using FEniCS for the discretization and PETSc for solving · Wedge flow benchmark results compare well with literature results. Optimal preconditioner for Stokes part, good solver for advection-diffusion part.