Document Sample

Iterative linear solvers for PDE-constrained optimization involving ﬂuid ﬂow Andy Wathen Oxford University joint work with Tyrone Rees, Martin Stoll and Sue Thorne Dave Grifﬁths Meeting, 2010 – p.1/32 Let A ∈ Rn×n be symmetric and positive deﬁnite 1 T minx 2 x Ax − xT b equivalent to solving Ax = b Dave Grifﬁths Meeting, 2010 – p.2/32 Constrained Optimization and Saddle-point problems: 1 T minx 2 x Ax − xT b subject to Bx = c Dave Grifﬁths Meeting, 2010 – p.3/32 Constrained Optimization and Saddle-point problems: 1 T minx 2 x Ax − xT b subject to Bx = c Lagrangian: L(x, λ) = 1 xTAx − xT b + (Bx − c)T λ 2 λ : Lagrange multipliers. minx ⇒ Ax + B ⋆ λ =b maxλ ⇒ Bx =c A BT x b ⇒ = B 0 λ c (Benzi, Golub & Liesen (2005)) Dave Grifﬁths Meeting, 2010 – p.3/32 Classic problem of this type: the Stokes problem ‘Minimise energy subject to conserving mass’ ⇒ −∇2 y + ∇p = u subject to ∇·y = 0 y: velocity, p: pressure, u: body forces Dave Grifﬁths Meeting, 2010 – p.4/32 Classic problem of this type: the Stokes problem ‘Minimise energy subject to conserving mass’ ⇒ −∇2 y + ∇p = u subject to ∇·y = 0 y: velocity, p: pressure, u: body forces Mixed ﬁnite elements/MAC ﬁnite difference approximation: K BT y u K 0 = , K= in R2 B 0 p g 0 K K : discrete Laplacian Dave Grifﬁths Meeting, 2010 – p.4/32 K BT B 0 is symmetric indeﬁnite → MINRES Krylov subspace iterative solver (Paige & Saunders (1975)) Preconditioning? Dave Grifﬁths Meeting, 2010 – p.5/32 K BT B 0 is symmetric indeﬁnite → MINRES Krylov subspace iterative solver (Paige & Saunders (1975)) Preconditioning? Block diagonal/triangular preconditioners Dave Grifﬁths Meeting, 2010 – p.5/32 Block diagonal/triangular preconditioners: based on observation (Murphy, Golub & W (2000), Korzak(1999)) A BT B 0 preconditioned by √ A 0 • has 3 distinct eigenvalues (1, 1 ± 2 5 2 ) 0 S A BT • has 2 distinct eigenvalues 0 S where S = BA−1 B T (Schur Complement) ⇒ MINRES /GMRES terminates in 3 / 2 iterations ⇒ want approximations A, S ⇒ 3 / 2 clusters ⇒ fast convergence Dave Grifﬁths Meeting, 2010 – p.6/32 For Stokes problem: A = K is discrete Laplacians: use A =multigrid cycles • geometric mutigrid: relaxed Jacobi smoothing, standard grid transfers • algebraic multigrid: HSL routine HSL MI20 (Boyle, Mihajlovic & Scott (2007)) Dave Grifﬁths Meeting, 2010 – p.7/32 S for Stokes problem: Babuska-Brezzi stability ⇒ Schur Complement spectrally equivalent to the ﬁnite element identity matrix, the mass matrix, M ie. the Gram matrix of the ﬁnite element basis functions {φj , j = 1, . . . , n} in L2 (Ω), Ω ⊂ Rd : BA−1 B T ≈ ∇ · (∇2 )−1 ∇ ≈ Id → M √ Speciﬁcally ∃γ > 0, Γ ≤ d such that xT BA−1 B T x γ≤ ≤Γ xT M x where M = {mi,j }, mi,j = φi φj Ω (Bank, Welfert & Yserentant (’90), Silvester & W (’94), Grifﬁths (’95)) Dave Grifﬁths Meeting, 2010 – p.8/32 Mass matrix is effectively preconditioned by its diagonal D =diag(M ) (W (1987)) eg. for Q1 (bilinear) elements in 2-dimensions (rectangles) eigenvalues of D −1 M all in [1/4, 9/4] eg. for Q1 (trilinear) ﬁnite elements in 3-dimensions (bricks) eigenvalues of D −1 M all in [1/8, 27/8] independently of mesh size h (ie. independently of discrete problem dimension) Could use S = D but better a few iterations of a diagonally preconditioned iteration for M : diagonally scaled Conjugate Gradients: leads to a nonlinear preconditioner diagonally scaled Chebyshev (semi-)iteration: is linear and we have precise eigenvalue inclusion intervals (W & Rees (2008)) Dave Grifﬁths Meeting, 2010 – p.9/32 4 10 CG Chebyshev 2 10 −1 A 0 || r || 10 k −2 10 −4 10 0 5 10 15 20 Number of iterations Dave Grifﬁths Meeting, 2010 – p.10/32 For Stokes problem: A 0 AAM G 0 P = = 0 S 0 T20 (D −1 M ) Here: Q2-Q1 mixed ﬁnite element (Taylor-Hood), A = AAM G : 1 AMG V-cycle h Iterations CPU time (s) 2−2 (187) 19 0.015 2−3 (659) 24 0.073 2−4 (2,467) 26 0.082 2−5 (9,539) 28 0.21 2−6 (37,507) 29 3.80 2−7 (148,739) 29 15.5 Dave Grifﬁths Meeting, 2010 – p.11/32 PDE-constrained Optimization General problem: Given Ω ⊂ R2 or R3 , y ∈ L2 (Ω) as some desired state then for some (regularisation) parameter β 1 2 β 2 min y−y L2 (Ω) + u L2 (Ω) y,u 2 2 subject to Ly = u in Ω, y=y on ∂Ω where L represents a partial differential operator Dave Grifﬁths Meeting, 2010 – p.12/32 can also include simple bounds on the control: u≤u≤u and even bounds on the state (via Moreau-Yoshida regulrization) y≤y≤y also boundary control. . . Dave Grifﬁths Meeting, 2010 – p.13/32 Simple sample problem: desirable y, controlable body force u 1 2 β 2 min y−y L2 (Ω) + u L2 (Ω) y,u 2 2 subject to −∇2 y = u in Ω, y=y on ∂Ω Dave Grifﬁths Meeting, 2010 – p.14/32 1 2 β 2 min y−y + u y,u 2 2 subject to − ∇2 y = u in Ω, y=y on ∂Ω Discretisation: ﬁnite elements yh = yj φj , y = (y1 , y2 , . . . , yn )T uh = uj φj , u = (u1 , u2 , . . . , un )T 1 T T β min y My + y b + uT M u y,u 2 2 subject to Ky = M u + d M = {mi,j }, mi,j = Ω φi φj — mass matrix K = {ki,j }, ki,j = Ω ∇φi · ∇φj — stiffness matrix as before Dave Grifﬁths Meeting, 2010 – p.15/32 so, Lagrangian: 1 T T β y My + y b + uT M u + λT (Ky − M u − d) 2 2 stationarity ⇒ Saddle point system M 0 K T y −b 0 βM −M u = 0 K −M 0 λ d M 0 Note B= K −M and A= 0 βM in usual saddle point form A BT B 0 Dave Grifﬁths Meeting, 2010 – p.16/32 M 0 Recall B= K −M and A= 0 βM so S ? S = BA−1 B T (Schur Complement) M −1 0 KT = K −M 1 0 β M −1 −M 1 = M + KM −1 K T β Unless approx β < 10−6 dominant part is S = KM −1 K T Dave Grifﬁths Meeting, 2010 – p.17/32 Hence preconditioner for M 0 K T M 0 0 A = 0 βM −M is P = 0 βM 0 K −M 0 0 0 KM −1 K T Eigenvalues ν of P −1 A ν = 1, 1 2α1 h4 1 2α2 1 + 5+ ≤ ν ≤ 1+ 5+ 2 β 2 β 1 2α2 1 2α1 h4 or 1− 5+ ≤ ν ≤ 1 − 5+ , 2 β 2 β where α1 , α2 are positive constants independent of h. Dave Grifﬁths Meeting, 2010 – p.18/32 M 0 0 P = 0 βM 0 , β = 10−2 0 0 KM −1 K T 1.5 1 0.5 0 −0.5 0 5 10 15 20 25 Dave Grifﬁths Meeting, 2010 – p.19/32 But M 0 0 P = 0 βM 0 0 0 KM −1 K T still expensive to use in practice so employ approximations M ≃M and K ≃ K giving M 0 0 A 0 P = 0 βM 0 = 0 S 0 0 KM −1 K T Dave Grifﬁths Meeting, 2010 – p.20/32 Important subtlety: K≃K does not imply that KM −1 K T ≃ KM −1 K T or indeed that K K T ≃ KK T without further conditions which are satisﬁed in this case (Braess & Peisker (1986)) Dave Grifﬁths Meeting, 2010 – p.21/32 M 0 0 A 0 P = 0 βM 0 = 0 S 0 0 KM −1 K T so M ?: T20 (D −1 M ) as before For L = −∇2 , K is a discrete Laplacian: use multigrid cycles as before In our examples: K is the action of 2 V-cycles M is the action of 20 Chebyshev semi-iterative steps Dave Grifﬁths Meeting, 2010 – p.22/32 Example problem: Ω = [0, 1]d , d = 2, 3 1 2 β 2 min y−y L2 (Ω) + u L2 (Ω) y,u 2 2 subject to −∇2 y = u in Ω, y=y on ∂Ω Q1 (bilinear) ﬁnite elements, β = 10−2 1 y: 0.8 0.6 0.4 0.2 0 1 1 0.5 0.8 0.6 0.4 0.2 0 0 Dave Grifﬁths Meeting, 2010 – p.23/32 CPU times (MINRES iterations) in 2D, tol 10−6 h 3n backslash MINRES MINRES (PAM G ) (PM G ) 2−2 27 0.0003 0.02 (7) 0.13 (7) 2−3 147 0.002 0.03 (9) 0.16 (9) 2−4 675 0.01 0.05 (9) 0.21 (9) 2−5 2883 0.08 0.14 (9) 0.41 (9) 2−6 11907 0.46 0.61 (9) 1.29 (9) 2−7 48387 3.10 2.61 (9) 5.09 (9) 2−8 195075 15.5 15.0 (11) 23.6 (9) 2−9 783363 — 75.6 (11) 136 (9) Dave Grifﬁths Meeting, 2010 – p.24/32 Control: 0 −0.05 −0.1 −0.15 −0.2 1 0.8 1 0.6 0.8 0.4 0.6 0.4 0.2 0.2 0 0 Dave Grifﬁths Meeting, 2010 – p.25/32 CPU times (MINRES iterations) in 3D, tol 10−6 h 3n backslash MINRES MINRES (PAM G ) (PM G ) 2−2 81 0.001 0.02 (7) 0.14 (8) 2−3 1029 0.013 0.13 (9) 0.26 (8) 2−4 10125 25.1 1.89 (8) 1.69 (8) 2−5 89373 — 22.1 (8) 15.9 (8) 2−6 750141 — — 230 (10) Dave Grifﬁths Meeting, 2010 – p.26/32 Stokes Control 1 2 1 2 β 2 min ˆ y−y L2 (Ω) + ˆ p−p L2 (Ω) + u L2 (Ω) y,p,u 2 2 2 subject to − ∇2 y + ∇p = u ∇·y = 0 y: velocity, p: pressure. Mixed ﬁnite elements for (forward) Stokes problem: K BT y u K 0 = , K= in R2 B 0 p g 0 K Dave Grifﬁths Meeting, 2010 – p.27/32 Cost functional 1 T T 1 T T β y My y − y b + p Mp p − p d + u T Mu u 2 2 2 combined with constraint via the Lagrangian ⇒ My 0 0 K BT y b 0 Mp 0 B 0 p d 0 0 βMu −Mu 0 u = 0 . K B T −M u 0 0 λ h B 0 0 0 0 µ k Dave Grifﬁths Meeting, 2010 – p.28/32 Schur Complement: K BT My−1 0 K BT 1 −1 + Mu B 0 0 Mp B 0 β and again ignore 2nd term for moderate β So overall block diagonal preconditioner requires: My , Mp , Mu → Chebyshev and 2 Stokes approximations Dave Grifﬁths Meeting, 2010 – p.29/32 Stokes preconditioners: K BT K 0 = B 0 B Mp and K BT K BT = B 0 0 Mp on left and right respectively where K is multigrid cycles for each discrete scalar Laplacian as before (Silvester & W (1993), Klawonn (1998)) Gives symmetric Schur complement approximation Dave Grifﬁths Meeting, 2010 – p.30/32 Here: Q2-Q1 mixed ﬁnite elements for cavity ﬂow 4 AMG V-cycles to approx each K 20 Chebyshev semi-iterations for each M h Iterations CPU time (s) 2−2 (344) 26 0.48 2−3 (1,512) 31 1.05 2−4 (6,344) 33 3.69 2−5 (25,992) 33 18.0 2−6 (105,224) 34 84.2 2−7 (423,432) 34 342 recall: 29 MINRES iterations for forward Stokes solve (but only 1 AMG V-cycle) ⇒ approx 10 times more work to solve the control problem than a single PDE solve. Dave Grifﬁths Meeting, 2010 – p.31/32 References Rees, T., Dollar, H.S. & Wathen, A.J., 2009, ‘Optimal solvers for PDE-constrained Optimization’, to appear in SIAM J. Sci. Comput., Oxford University NA report 08/10. Murphy, M.F., Golub, G.H. & Wathen, A.J., 2000, ‘A note on preconditioning for indeﬁnite linear systems’, SIAM J. Sci. Comput. 21(6), pp. 1969-1972. Silvester, D.J. & Wathen, A.J., 1994, ‘Fast iterative solution of stabilised Stokes systems Part II: Using general block preconditioners’, SIAM J. Numer Anal. 31, 1352-1367. Wathen, A.J. & Rees, T., 2009, ‘Chebyshev semi-iteration in preconditioning’, Victor Pereyra special issue of ETNA. Benzi, M, Golub, G.H. & Liesen, J., 2005, ‘Numerical Solution of Saddle Point Problems’, Acta Numerica 2005. Dave Grifﬁths Meeting, 2010 – p.32/32

DOCUMENT INFO

Shared By:

Categories:

Tags:
iterative linear solvers, iterative solvers, iterative methods, linear solver, linear systems, the matrix, scientific computing, Finite Elements, linear solvers, Andrew J. Wathen

Stats:

views: | 13 |

posted: | 4/4/2011 |

language: | English |

pages: | 35 |

OTHER DOCS BY nyut545e2

Docstoc is the premier online destination to start and grow small businesses. It hosts the best quality and widest selection of professional documents (over 20 million) and resources including expert videos, articles and productivity tools to make every small business better.

Search or Browse for any specific document or resource you need for your business. Or explore our curated resources for Starting a Business, Growing a Business or for Professional Development.

Feel free to Contact Us with any questions you might have.