Multigrid preconditioners for linear systems arising in PDE - PowerPoint PPT Presentation

Model problems Unconstrained problems with linear PDE constraints Nonlinear constraints, control constraints Summary Multigrid preconditioners for linear systems arising in PDE constrained optimization Andrei Draganescu Department of Mathematics and Statistics University of Maryland, Baltimore County March 5, 2013 Andrei Draganescu, UMBC NIST, March 5, 2013

Model problems Unconstrained problems with linear PDE constraints Nonlinear constraints, control constraints Summary Acknowledgments Former collaborators: Todd Dupont (U Chicago), Volkan Akçelik (Exxon), George Biros (U Texas), Omar Ghattas (U Texas), Judith Hill (ORNL), Cosmin Petra (ANL), Bart van Bloemen Waanders (Sandia). Current collaborators: (UMBC) Mona Hajghassem , Jyoti Saraswat , Ana Maria Soane Grants: NSF awards DMS-1016177 and DMS-0821311. DOE contract no: DE-SC0005455. Andrei Draganescu, UMBC NIST, March 5, 2013

Model problems Unconstrained problems with linear PDE constraints Nonlinear constraints, control constraints Summary Outline Model problems 1 Unconstrained problems with linear PDE constraints 2 Nonlinear constraints, control constraints 3 A semilinear elliptic constrained problem Control-constrained problems Optimal control problems constrained by the Stokes equations Andrei Draganescu, UMBC NIST, March 5, 2013

Model problems Unconstrained problems with linear PDE constraints Nonlinear constraints, control constraints Summary Abstract problem formulation  J ( y , u ) = 1 2 || y − y d || 2 L 2 (Ω) + R ( u , y ) , minimize    (1) u ∈ U ad ⊂ U , y ∈ Y ad ⊂ Y , subj . to    e ( y , u ) = 0 . U ad and Y ad – sets of admissible controls resp. states (convex, closed, non-empty). Ex.: U ad = { u ∈ U : u ≤ u ≤ u } , Y ad = { y ∈ Y : y ≤ y ≤ y } . Equality constraint is a well-posed PDE: for all u ∈ U there is a unique y ∈ Y (depending continuously on u ), so that def e ( y , u ) = 0 , y = K ( u ) . Andrei Draganescu, UMBC NIST, March 5, 2013

Model problems Unconstrained problems with linear PDE constraints Nonlinear constraints, control constraints Summary Reduced problem formulation If U ad = U and Y ad = Y , problem can be reformulated as unconstrained: J ( u ) = 1 | 2 + β | 2 , u ∈ U ad . 2 | | K ( u ) − y d | 2 | | Lu | min (2) u ∈ U If β ≪ 1, essentially we want solve K ( u ) = y d . However, problems of interest are ill-posed, need regularization: L = I ⇒ find u of smallest norm ; L = ∇ ⇒ find u of smallest variation. Andrei Draganescu, UMBC NIST, March 5, 2013

Model problems Unconstrained problems with linear PDE constraints Nonlinear constraints, control constraints Summary Motivating applications 1. Reverse advection-diffusion problems (source inversion): T > 0 fixed “end-time”, y d end-time state, u initial state z ( · , t ) transported quantity subjected to:  ∂ t z − ∇ · ( a ∇ z + bz ) + cz = 0 on Ω  z ( x , t ) = 0 for x ∈ ∂ Ω , t ∈ [ 0 , T ]  z ( x , 0 ) = u ( x ) for x ∈ Ω K = S ( T ) : initial - to - final def K u = S ( T ) u = z ( · , T ) Andrei Draganescu, UMBC NIST, March 5, 2013

Model problems Unconstrained problems with linear PDE constraints Nonlinear constraints, control constraints Summary Further motivating applications 2. Elliptic optimal control problem: PDE-constrained optimal control problem  | 2 + β | 2 , 1 minimize 2 | | y − y d | 2 | | u |   subj to: − ∆ y = u , u | ∂ Ω = 0 ,   u ≤ u ≤ u . If unconstrained, then K = ( − ∆) − 1 . Andrei Draganescu, UMBC NIST, March 5, 2013

Model problems Unconstrained problems with linear PDE constraints Nonlinear constraints, control constraints Summary The case of linear constraints Assume K linear, U ad = U : u J ( u ) = 1 | 2 + β | 2 2 | | Ku − y d | 2 | | u | min Newton’s method gives the solution explicitly in one step: u min = u guess − G − 1 ∇ J ( u guess ) , 0 0 where G = G ( β ) = I + β − 1 K ∗ · K , ∇ J ( u ) = u + β − 1 K ∗ ( Ku − y d ) . Formulation is equivalent to the regularized normal equations ( β I + K ∗ · K ) u = K ∗ y d . Andrei Draganescu, UMBC NIST, March 5, 2013

Model problems Unconstrained problems with linear PDE constraints Nonlinear constraints, control constraints Summary Strategy: discretize-then-optimize Natural FE discretization for the operator K : u J ( u ) = 1 | 2 + β | 2 . min 2 | | K h u − y d | 2 | | u | Solution of discrete problem: min = u guess − 1 ∇ J h ( u guess u h − G h 0 , h ) , 0 , h where ∗ · K h , G h = G h ( β ) = I + β − 1 K h ∇ J h ( u ) = u + β − 1 K h ∗ ( K h u − π h y d ) , π h is the orthogonal projection onto the finite element space V h Main problem: need to invert the operator G h efficiently. Andrei Draganescu, UMBC NIST, March 5, 2013

Model problems Unconstrained problems with linear PDE constraints Nonlinear constraints, control constraints Summary Main issues The matrix representing the linear operator G h is dense, potentially large, and not available. Matrix-vector product cost is comparable to two forward computations (expensive, but feasible): G h · u = u + β − 1 K ∗ h · K h u . Gradient computation also costs as much as two forward computations (only done once): ∇ J h ( u ) = u + β − 1 K ∗ h ( K h u − π h y d ) . Need iterative methods. Andrei Draganescu, UMBC NIST, March 5, 2013

Model problems Unconstrained problems with linear PDE constraints Nonlinear constraints, control constraints Summary Solution using conjugate gradient Eigenvalues of G h cluster around 1 ⇒ CG is a good choice for solving inverting G h : the number of iterations is independent of the resolution; grows only logarithmically with β → 0. A measure of success: speedup over CG. Andrei Draganescu, UMBC NIST, March 5, 2013

Model problems Unconstrained problems with linear PDE constraints Nonlinear constraints, control constraints Summary Multigrid strategies major differences between G h and an elliptic operator A h : G h A h smoothing roughening nonlocal local cond ( G h ) bounded cond ( A h ) → ∞ as h → 0 Related multigrid work: Hackbusch (1981), King (1992), Rieder (1997), Hanke and Vogel (1999), Kaltenbacher (2003), Donatelli (2005), Biros and Dogan (2008), Draganescu and Dupont (2008), Borzi and Kunisch (2005). Lewis and Nash (2000) overview: Borzi and Schultz (SIAM review, 2009) more recent: Wathen, Stoll, Rees, Dollar, Draganescu and Soane, etc Andrei Draganescu, UMBC NIST, March 5, 2013

Model problems Unconstrained problems with linear PDE constraints Nonlinear constraints, control constraints Summary Two-grid approximation (heuristics) “ rough “ functions “ smooth “ functions ���� ���� ⊕ V h = V 2 h W denote π = π 2 h , ρ = I − π 2 h G h = π G h π + ρ G h π + π G h ρ + ρ G h ρ � �� � � �� � � �� � � �� � M 1 M 2 M 3 M 4 � � I + β − 1 K h ∗ K h since G h ρ = ρ ≈ ρ M 2 ≈ 0 M 3 ≈ 0 M 1 ≈ G 2 h π M 4 ≈ ρ conclusion: def G h ≈ M h = G 2 h π 2 h ⊕ ( I − π 2 h ) . Andrei Draganescu, UMBC NIST, March 5, 2013

Model problems Unconstrained problems with linear PDE constraints Nonlinear constraints, control constraints Summary Two-grid approximation (heuristics) “ rough “ functions “ smooth “ functions ���� ���� ⊕ V h = V 2 h W denote π = π 2 h , ρ = I − π 2 h G h = π G h π + ρ G h π + π G h ρ + ρ G h ρ � �� � � �� � � �� � � �� � M 1 M 2 M 3 M 4 � � I + β − 1 K h ∗ K h since G h ρ = ρ ≈ ρ M 2 ≈ 0 M 3 ≈ 0 M 1 ≈ G 2 h π M 4 ≈ ρ conclusion: def G h ≈ M h = G 2 h π 2 h ⊕ ( I − π 2 h ) . Andrei Draganescu, UMBC NIST, March 5, 2013

Model problems Unconstrained problems with linear PDE constraints Nonlinear constraints, control constraints Summary Two-grid approximation (heuristics) “ rough “ functions “ smooth “ functions ���� ���� ⊕ V h = V 2 h W denote π = π 2 h , ρ = I − π 2 h G h = π G h π + ρ G h π + π G h ρ + ρ G h ρ � �� � � �� � � �� � � �� � M 1 M 2 M 3 M 4 � � I + β − 1 K h ∗ K h since G h ρ = ρ ≈ ρ M 2 ≈ 0 M 3 ≈ 0 M 1 ≈ G 2 h π M 4 ≈ ρ conclusion: def G h ≈ M h = G 2 h π 2 h ⊕ ( I − π 2 h ) . Andrei Draganescu, UMBC NIST, March 5, 2013

Model problems Unconstrained problems with linear PDE constraints Nonlinear constraints, control constraints Summary Two-grid approximation (heuristics) “ rough “ functions “ smooth “ functions ���� ���� ⊕ V h = V 2 h W denote π = π 2 h , ρ = I − π 2 h G h = π G h π + ρ G h π + π G h ρ + ρ G h ρ � �� � � �� � � �� � � �� � M 1 M 2 M 3 M 4 � � I + β − 1 K h ∗ K h since G h ρ = ρ ≈ ρ M 2 ≈ 0 M 3 ≈ 0 M 1 ≈ G 2 h π M 4 ≈ ρ conclusion: def G h ≈ M h = G 2 h π 2 h ⊕ ( I − π 2 h ) . Andrei Draganescu, UMBC NIST, March 5, 2013