Howard Elman University of Maryland Vicki Howles John Shadid - PowerPoint PPT Presentation

Howard Elman University of Maryland

• Vicki Howles • John Shadid • David Kay • David Silvester • Daniel Loghin • Ray Tuminaro • Alison Ramage • Andy Wathen 2

Incompressible Navier-Stokes Equations α − ν ∇ + ⋅ + = 2 ( grad ) grad u u u u p f t − = div 0 u α =0 ! steady state problem α =1 ! evolutionary problem Discretization and linearization Matrix equation x=b       T u f F B       =       −    0    p B C Goal: Robust general solution algorithms Easy to implement Derived from subsidiary building blocks Adaptible to a variety of scenarios (steady / evolutionary / Stokes / Boussinesq ) 3

Outline 1. General approach preconditioning for saddle point problems 2. Relation to traditional approaches projection methods SIMPLE 3. Details for Navier - Stokes equations 4. Analytic / experimental results 5. Potential for more general problems 4

Preliminary: Steady Stokes Equations − ∇ + = 2 grad Rusten & Winther, 1992 u p f − = div 0 Silvester & Wathen , 1993 u       Algebraic equation: A = discrete vector T u f A B     =   ,      0  Laplacian −       B C p Symmetric indefinite ! MINRES algorithm is applicable   0 Preconditioning operator: A      0  Q S         0 T Generalized eigenvalue problem: A B u A u       = λ           −        0  p Q p B C S Au+ B T p = λ Au , λ ≠ 1 ) u= 1/( λ - 1) A -1 B T p Case C=0: Bu = λ Q S p BA -1 B T p= λ ( λ - 1) Q S p

BA -1 B T p= µ Q S p, µ = λ ( λ - 1) Verf ü rth , 1984: For Q S = pressure mass matrix, µ 2 [a s ,b s ] independent of discretization parameter h µ λ =1 § (1+4 µ ) 1/2 , Under mapping λ ½ [ ] [ ] 0 d -a -b c / 2 k   − 1 ( ) /( ) bc ad   ≤ Convergence bound for MINRES : || || 2 || || r r   0 k + 1 ( ) /( )  bc ad  − 1   0 Computational requirements, for times a vector A      0  Q S Poisson solve: can be approximated, e.g. with multigrid Mass matrix solve: cheap 6

Generalize to Navier-Stokes Equations Linearization via Picard iteration (slightly inaccurate notation) : − ν ∇ + ⋅ + = α + + + + ( 1 ) 2 ( 1 ) ( ) ( 1 ) ( 1 ) ( ) m m ( m grad ) m grad m m u u u u p f ∆ t − = + ( 1 ) div m 0 u       T u f F B = Discretization           0  0  p B   0 Analogue of Stokes strategy: preconditioner F      0  Q S BF -1 B T p= µ Q S p , µ = λ ( λ - 1) Same analysis ! ) seek approximation Q S to Schur complement N.B. Same question arises for other strategies for linearization 7

Suppose Q S ¼ BF -1 B T so that eigenvalues of BF -1 B T p = µ Q S p are tightly clustered. Under mapping µ λ =1 § (1+4 µ ) 1/2 , eigenvalues λ are clustered in two regions, one on each side of imaginary axis Can improve this: Observation: symmetry is important for Stokes solver MINRES: optimal with fixed cost per step ) need block diagonal preconditioner For Navier - Stokes: don’t have symmetry need Krylov subspace method for nonsymmetric matrices (e.g. GMRES) 8

  Alternative: block triangular preconditioner T F B     −  0  Q S         ! Generalized eigenvalue problem T T u F B u F B         = µ         −     0   0   B p Q p S Fu+ B T p = µ Fu , µ ≠ 1 ) u= -F -1 B T p Bu = - µ Q S p BF -1 B T p= µ Q S p µ ½ {1} [ Theorem (Fischer, Ramage, Silvester, Wathen ): For preconditioned GMRES iteration, let p 0 be arbitrary and u 0 = F -1 (f -B T u 0 ) ( ) r 0 =(0,w 0 )) (u k ,p k ) T be generated with block triangular preconditioner , (u k ,p k ) D be generated with block diagonal preconditioner . Then (u 2k ,p 2k ) D = (u 2k+1 ,p 2k+1 ) D = ( u k ,p k ) T . 9

Computational requirements, to implement block triangular preconditioner : − 1       T w F B v     =   Compute       −      0  r q Q S Solve Q s r = - q, then solve Fw = v -B T r The only difference from block diagonal solve: matrix - vector product B T r (negligible) For second step: convection - diffusion solve: can be approximated, e.g. with multigrid For first step: something new needed: Q S 10

One more interpretation:       0 T T F B I F B       =       − − + − −  1  1    0 ( T )  B C BF I BF B C − 1       T T 0 I F B F B     =   )       − − + − − 1    1    0 ( T ) BF I B C BF B C − 1 ¼   T Q B   F   −  0  Q S Shows what is needed for stabilized discretization : Q S ¼ BF -1 B T + C 11

Relation to Projection Methods “Classical” O( ∆ t ) projection method ( Chorin 1967, Temam 1969): − (*) ( ) m u u Step 1: − ν ∇ + ⋅ = 2 * ( ) ( ) ( m grad ) m u u u f ( ) ∆ t − ν ∇ = − ⋅ + 1 1 2 * ( ) ( ) ( ) ( m grad ) m m I u f u u u ∆ ∆ In matrix form: ( t ) t +ν = − + 1 1 * ( ) ( ) m m M A u f Nu Mu ∆ ∆ t t ∇     +   ( 1 ) m 1 1 * I u Step 2: u     =   ∆ ∆ t t       − ∇ + ( 1 )   m   0 0   p     +   ( 1 ) In matrix form: T m 1 * Mu 1 M B u     =   ∆ ∆ t t       + ( 1 ) m   0  0    B p Performed via pressure - Poisson solve 12

Substitute u * from Step 1 into Step 2 : ( )( )   − − +     + + ν + ν 1 ( 1 ) m 1 1 1 u f Nu ( ) Mu ( ) T m m 1 M A M A M B       = ∆ ∆ ∆ ∆ t t t   t     + ( 1 )    m  0   p 0 B ( ) + ν   −   ( ) 1 0 1 M A 1 T I M B     ∆ ∆ t t Perot, 1993     − − 1 1 T   B B M B   0 I ∆ t Contrast: update derived purely from linearization & discretization : − + M ν +       + ( 1 ) T m 1 1 A B u f Nu ( ) Mu ( ) m m       = ∆ ∆ t t       + ( 1 )    0   m  0 B p ( ) + ν   −   ( ) + ν 1 0 1 M A 1 T I M A B     ∆ ∆ t t     − − + ν 1 1 T   B B M A B   0 I ( )( ) ∆ t − − + ν 1 = − ∆ ν = ∆ − 1 1 1 Error : T T ( ) T ( ) B M A M B t M AB O t ∆ ∆ t t 13

For higher order accuracy in time and related approaches: Dukowicz & Dvinsky 1992 Perot 1993 Quarteroni, Saleri & Veneziani 2000 Henriksen & Holman 2002 14

Relation to SIMPLE Patankar & Spaulding, 1972       − 0 1 T T F B F I F B       =       − − 1    T    0 0 B B BF B I     ˆ ¼ − 0 Q 1 T I F B     F     − ˆ −  1  T B B F B  0  I Q F : approximate convection -diffusion solve ^ F: diagonal part of F Many variants 15

Perspective of New Approach • Take on Schur complement directly • Separate time discretization from algebraic algorithm • Enable flexible treatment of time discretization, linearization Allow choice of linearization Allow large time steps / CFL numbers if circumstances warrant 16

Kay, Loghin , Approximation for the Schur Complement (I) Wathen 2002 Suppose the gradient and convection - diffusion operators approximately commute (w=u (m) ): ∇ − ν ∇ + ⋅ ∇ ≈ − ν ∇ + ⋅ ∇ ∇ 2 2 ( ) ( ) w w p u Require pressure convection - diffusion operator − − − − = 1 1 1 1 T T Discrete analogue: M B M F M F M B u p p u u ⇒ − = − − 1 1 1 T T BF B BM B F M u p p A p ≡ − 1 Q A F M In practice: don’t have equality, instead S p p p − 1 Requirements: Poisson solve A − 1 p for Q − 1 S M Mass matrix solve p − 1 + Convection -diffusion solve F 17