Robust discretization and solvers for elliptic optimal control problems with energy regularization
RRobust discretization and solvers for elliptic optimalcontrol problems with energy regularization
Ulrich Langer ∗ , Olaf Steinbach † , Huidong Yang ‡ Abstract
We consider the finite element discretization and the iterative solution ofsingularly perturbed elliptic reaction-diffusion equations in three-dimensionalcomputational domains. These equations arise from the optimality condi-tions for elliptic distributed optimal control problems with energy regulariza-tion that were recently studied by M. Neum¨uller and O. Steinbach (2020).We provide quasi-optimal a priori finite element error estimates which de-pend both on the mesh size h and on the regularization parameter (cid:37) . Thechoice (cid:37) = h ensures optimal convergence which only depends on the reg-ularity of the target function. For the iterative solution, we employ an al-gebraic multigrid preconditioner and a balancing domain decomposition byconstraints (BDDC) preconditioner. We numerically study robustness andefficiency of the proposed algebraic preconditioners with respect to the meshsize h , the regularization parameter (cid:37) , and the number of subdomains (cores) p . Furthermore, we investigate the parallel performance of the BDDC precon-ditioned conjugate gradient solver. Keywords:
Elliptic optimal control problems, energy regularization,finite element discretization, a priori error estimates, fast solvers
In the recent work [22], M. Neum¨uller and O. Steinbach have investigated regu-larization error estimates for the solution of distributed optimal control problemsin energy spaces on the basis of the following model problem: Minimize the costfunctional J ( u (cid:37) , z (cid:37) ) = 12 (cid:90) Ω [ u (cid:37) ( x ) − u ( x )] dx + (cid:37) (cid:107) z (cid:37) (cid:107) H − (Ω) (1)with respect to the state u (cid:37) ∈ H (Ω) and the control z (cid:37) ∈ H − (Ω) subject to theconstraints − ∆ u (cid:37) = z (cid:37) in Ω , u (cid:37) = 0 on Γ , (2)where u ∈ L (Ω) represents the given desired state (target), and (cid:37) ∈ R + denotesthe regularization parameter. Here, the computational domain Ω ⊂ R is assumed ∗ Institute for Computational Mathematics, Johannes Kepler University Linz, AltenbergerStraße 69, 4040 Linz, Austria, Email: [email protected] † Institut f¨ur Angewandte Mathematik, Technische Universit¨at Graz, Steyrergaße 30, 8010 Graz,Austria, Email: [email protected] ‡ Johann Radon Institute for Computational and Applied Mathematics, Austrian Academy ofSciences, Altenberger Straße 69, 4040 Linz, Austria, Email: [email protected] a r X i v : . [ m a t h . NA ] F e b o be a bounded Lipschitz domain with boundary Γ = ∂ Ω. Note that (cid:107) z (cid:37) (cid:107) H − (Ω) = (cid:107)∇ u (cid:37) (cid:107) L (Ω) . The associated adjoint state equation reads: Find the adjoint state p (cid:37) ∈ H (Ω)such that − ∆ p (cid:37) = u (cid:37) − u in Ω , p (cid:37) = 0 on Γ . (3)Further, the associated gradient equation is given by the relation p (cid:37) + (cid:37) u (cid:37) = 0 in Ω . (4)Using the gradient equation (4) for eliminating the adjoint state p (cid:37) from (3), wefinally obtain the singularly perturbed reaction-diffusion equation − (cid:37) ∆ u (cid:37) + u (cid:37) = u in Ω , u (cid:37) = 0 on Γ , (5)for defining the state u (cid:37) . Once the state u (cid:37) is computed from (5), the adjoint state p (cid:37) is given by (4), and the control z (cid:37) = − ∆ u (cid:37) = (cid:37) − ( u − u (cid:37) ) by (2) and (5).The variational formulation of the Dirichlet boundary value problem (5) is tofind u (cid:37) ∈ H (Ω) such that (cid:37) (cid:90) Ω ∇ u (cid:37) ( x ) · ∇ v ( x ) dx + (cid:90) Ω u (cid:37) ( x ) v ( x ) dx = (cid:90) Ω u ( x ) v ( x ) dx (6)is satisfied for all v ∈ H (Ω). When assuming some regularity of the given target u ∈ H s (Ω) := [ L (Ω) , H (Ω)] s for s ∈ [0 ,
1] or u ∈ H (Ω) ∩ H s (Ω) for s ∈ (1 , (cid:37) has beenshown in [22]: (cid:107) u (cid:37) − u (cid:107) L (Ω) ≤ c (cid:37) s/ (cid:107) u (cid:107) H s (Ω) . (7)Moreover, in the case u ∈ H (Ω) ∩ H s (Ω) for s ∈ (1 , (cid:107)∇ ( u (cid:37) − u ) (cid:107) L (Ω) ≤ c (cid:37) ( s − / (cid:107) u (cid:107) H s (Ω) . (8)For (cid:37) →
0, the boundary value problem (5) belongs to a class of singularly perturbedproblems, see, e.g., [7]. For robust numerical methods to treat such problems, werefer to the monograph [8].In this paper, we investigate the errors (cid:107) u (cid:37)h − u (cid:107) L (Ω) and (cid:107)∇ ( u (cid:37)h − ¯ u ) (cid:107) L (Ω) of the finite element approximation u (cid:37)h to the target u in terms of both the reg-ularization parameter (cid:37) and the discretization parameter h . The error estimate isbased on the estimates (7) and (8), and well-known finite element discretizationerror estimates. It turns out that the choice (cid:37) = h gives the optimal convergencerate O ( h ) for continuous, piecewise linear finite elements. Moreover, we consideriterative methods for the solution of the linear system of algebraic equations thatarise from the Galerkin finite element discretization of the variational problem (6),which are robust with respect to the mesh size h and the regularization parame-ter (cid:37) . Such solvers have been addressed in many works. For example, geometricmultigrid methods, which are robust with respect to the mesh size h and the pa-rameter (cid:37) in the case of standard Galerkin discretization methods, were analyzed in[25], based on the convergence rate min { , h /(cid:37) } in the L norm as provided in [28].The parameter independent contraction number of the multigrid method is achievedby a combination of such a deteriorate approximation property with an improvedsmoothing property. Robust and optimal algebraic multilevel (AMLI) methods forreaction-diffusion type problems were studied in [14]. Therein, a uniformly con-vergent AMLI method with optimal complexity was shown using the constant inthe strengthened Cauchy-Bunyakowski-Schwarz inequality computed for both mass2nd stiffness matrices. Robust block-structured preconditioning on both piecewiseuniform meshes and graded meshes with boundary layers has been considered in[19, 23] in comparison with standard robust multigrid preconditioners. The newblock diagonal preconditioner was based on the partitioning of degrees of freedominto those on the corners, edges, and interior points, respectively. The perturbationparameter independent condition number of the preconditioned linear system usingdiagonal and incomplete Cholesky preconditioning methods for singularly perturbedproblems on layer-adapted meshes has been recently shown in [24]. Considering theadaptive finite element discretization on simplicial meshes with local refinements,we employ an algebraic multigrid (AMG) preconditioner [4, 5, 27] and the balancingdomain decomposition by constraints (BDDC) preconditioner [6, 20, 21] for solvingthe discrete system. We make performance studies of the preconditioned conjugategradient (PCG) method that uses these AMG and BDDC preconditioners. In par-ticular, the BDDC preconditioned conjugate gradient solver shows excellent strongscalability.The reminder of the paper is organized as follows: Section 2 deals with the finiteelement discretization of the equation (5). In Section 3, we describe both the AMGand BDDC preconditioners that are used in solving the system of finite elementequations. Numerical results are presented and discussed in Section 4. Finally,some conclusions are drawn in Section 5. To perform the finite element discretization of the variational form (6), we introduceconforming finite element spaces V h ⊂ H (Ω). Particularly, we use the standardfinite element space V h = S h (Ω h ) ∩ H (Ω) spanned by continuous and piecewiselinear basis functions. These functions are defined with respect to some admissibledecomposition T h (Ω) of the domain Ω into shape regular simplicial finite elements τ (cid:96) , and are zero on Γ. Here, Ω h = (cid:83) (cid:96) τ (cid:96) , and h denotes a suitable mesh-sizeparameter, see, e.g., [3, 29]. Then the finite element approximation of (6) is to find u (cid:37)h ∈ V h such that (cid:37) (cid:90) Ω ∇ u (cid:37)h ( x ) · ∇ v h ( x ) dx + (cid:90) Ω u (cid:37)h ( x ) v h ( x ) dx = (cid:90) Ω u ( x ) v h ( x ) dx (9)for all test functions v h ∈ V h .In [22], the convergence of u (cid:37) towards the target u , and the error estimates(7) and (8) were proved. In the numerical experiments, presented in [22], u (cid:37) wascomputed on a very fine grid to minimize the influence of the discretization. Nowwe are going to investigate the effects of the finite element discretization, and thefinal aim is to provide estimates of the errors (cid:107) u (cid:37)h − u (cid:107) L (Ω) and (cid:107)∇ ( u (cid:37)h − u ) (cid:107) L (Ω) in terms of both (cid:37) and h for sufficiently smooth target functions u . Theorem 1.
Let us assume that Ω ⊂ R is convex, and that the target functionsatisfies u ∈ H (Ω) ∩ H (Ω) . Then there are positive constants C , C , C and C such that the error estimates (cid:107) u (cid:37)h − u (cid:107) L (Ω) ≤ (cid:16) C h + C (cid:37)h + C (cid:37) (cid:17) (cid:107) u (cid:107) H (Ω) (10) and (cid:107)∇ ( u (cid:37)h − u ) (cid:107) L (Ω) ≤ (cid:16) C h (cid:37) − + C h + C (cid:37) (cid:17) (cid:107) u (cid:107) H (Ω) (11) hold. For the choice (cid:37) = h , we therefore have (cid:107) u (cid:37)h − u (cid:107) L (Ω) ≤ c h (cid:107) u (cid:107) H (Ω) and (cid:107)∇ ( u (cid:37)h − u ) (cid:107) L (Ω) ≤ ˜ c h (cid:107) u (cid:107) H (Ω) . roof. When subtracting the Galerkin variational formulation (9) from (6) for testfunctions v h ∈ V h this gives the Galerkin orthogonality (cid:37) (cid:90) Ω ∇ ( u (cid:37) − u (cid:37)h ) · ∇ v h dx + (cid:90) Ω ( u (cid:37) − u (cid:37)h ) v h dx = 0 ∀ v h ∈ V h . Therefore, (cid:37) (cid:107)∇ ( u (cid:37) − u (cid:37)h ) (cid:107) L (Ω) + (cid:107) u (cid:37) − u (cid:37)h (cid:107) L (Ω) = (cid:37) (cid:90) Ω ∇ ( u (cid:37) − u (cid:37)h ) · ∇ ( u (cid:37) − u (cid:37)h ) dx + (cid:90) Ω ( u (cid:37) − u (cid:37)h ) ( u (cid:37) − u (cid:37)h ) dx = (cid:37) (cid:90) Ω ∇ ( u (cid:37) − u (cid:37)h ) · ∇ ( u (cid:37) − v h ) dx + (cid:90) Ω ( u (cid:37) − u (cid:37)h ) ( u (cid:37) − v h ) dx ≤ (cid:37) (cid:107)∇ ( u (cid:37) − u (cid:37)h ) (cid:107) L (Ω) (cid:107)∇ ( u (cid:37) − v h ) (cid:107) L (Ω) + (cid:107) u (cid:37) − u (cid:37)h (cid:107) L (Ω) (cid:107) u (cid:37) − v h (cid:107) L (Ω) ≤ (cid:113) (cid:37) (cid:107)∇ ( u (cid:37) − u (cid:37)h ) (cid:107) L (Ω) + (cid:107) u (cid:37) − u (cid:37),h (cid:107) L (Ω) · (cid:113) (cid:37) (cid:107)∇ ( u (cid:37) − v h ) (cid:107) L (Ω) + (cid:107) u (cid:37) − v h (cid:107) L (Ω) follows, i.e., we have Cea’s estimate (cid:37) (cid:107)∇ ( u (cid:37) − u (cid:37)h ) (cid:107) L (Ω) + (cid:107) u (cid:37) − u (cid:37)h (cid:107) L (Ω) ≤ (cid:37) (cid:107)∇ ( u (cid:37) − v h ) (cid:107) L (Ω) + (cid:107) u (cid:37) − v h (cid:107) L (Ω) for all v h ∈ V h . Inserting the Lagrangian interpolation v h = I h u (cid:37) of u (cid:37) , and usingstandard interpolation error estimates, we arrive at the error estimate (cid:37) (cid:107)∇ ( u (cid:37) − u (cid:37)h ) (cid:107) L (Ω) + (cid:107) u (cid:37) − u (cid:37)h (cid:107) L (Ω) ≤ c (cid:37) h | u (cid:37) | H (Ω) + c h | u (cid:37) | H (Ω) = (cid:2) c (cid:37) + c h (cid:3) h | u (cid:37) | H (Ω) , (12)where the positive constants c and c are nothing but the constants in the H and L interpolation error estimates; see, e.g., [3, 29]. Note that, under the assumptionsmade, we have u (cid:37) ∈ H (Ω) ∩ H (Ω). Due to − (cid:37) ∆ u (cid:37) = u − u (cid:37) , we now conclude | u (cid:37) | H (Ω) ≤ c (cid:37) (cid:107) u − u (cid:37) (cid:107) L (Ω) (13)with the positive H coercivity constant c . Since we assume u ∈ H (Ω) ∩ H (Ω),we can use (7) for s = 2, i.e., (cid:107) u − u (cid:37) (cid:107) L (Ω) ≤ c (cid:37) (cid:107) u (cid:107) H (Ω) . (14)For less regular u , we have a reduced order in (cid:37) ; see Theorem 3.2 in [22] as well as(7). Combining (13) and (14), we get | u (cid:37) | H (Ω) ≤ c c (cid:107) u (cid:107) H (Ω) . (15)Inserting (15) into (12) this gives (cid:37) (cid:107)∇ ( u (cid:37) − u (cid:37)h ) (cid:107) L (Ω) + (cid:107) u (cid:37) − u (cid:37)h (cid:107) L (Ω) ≤ c c (cid:2) c (cid:37) + c h (cid:3) h (cid:107) u (cid:107) H (Ω) . Using the triangle inequality, and again (14), we arrive at (cid:107) u (cid:37)h − u (cid:107) L (Ω) ≤ (cid:104) (cid:107) u (cid:37)h − u (cid:37) (cid:107) L (Ω) + (cid:107) u (cid:37) − u (cid:107) L (Ω) (cid:105) ≤ (cid:104) (cid:107) u (cid:37)h − u (cid:37) (cid:107) L (Ω) + (cid:107) u (cid:37) − u (cid:107) L (Ω) (cid:105) ≤ c c (cid:2) c (cid:37) + c h (cid:3) h (cid:107) u (cid:107) H (Ω) + 2 c (cid:37) (cid:107) u (cid:107) H (Ω) , C = 2 c c c , C = 2 c c c , and C = 2 c .Moreover, we also have (cid:107)∇ ( u (cid:37) − u (cid:37)h ) (cid:107) L (Ω) ≤ c c (cid:2) c (cid:37) + c h (cid:3) h (cid:37) − (cid:107) u (cid:107) H (Ω) . Now, using (8) for s = 2, we finally obtain (cid:107)∇ ( u (cid:37)h − u ) (cid:107) L (Ω) ≤ (cid:104) (cid:107)∇ ( u (cid:37)h − u (cid:37) ) (cid:107) L (Ω) + (cid:107)∇ ( u (cid:37) − u ) (cid:107) L (Ω) (cid:105) ≤ (cid:104) (cid:107)∇ ( u (cid:37)h − u (cid:37) ) (cid:107) L (Ω) + (cid:107)∇ ( u (cid:37) − u ) (cid:107) L (Ω) (cid:105) ≤ c c (cid:2) c (cid:37) + c h (cid:3) h (cid:37) − (cid:107) u (cid:107) H (Ω) + 2 c (cid:37) (cid:107) u (cid:107) H (Ω) , which is (11) with C = 2 c .The finite element space V h = span { ϕ , . . . , ϕ N h } is spanned by the standard nodalCourant basis functions { ϕ , . . . , ϕ N h } . Once the basis is chosen, the finite elementscheme (9) is equivalent to the system of finite element equations Au = f , (16)where the system matrix A = (cid:37)K + M consists of the scaled stiffness matrix K and the mass matrix M , respectively. The matrix entries and the entries of theright-hand side f are defined by K ij = (cid:90) Ω ∇ ϕ j · ∇ ϕ i dx, M ij = (cid:90) Ω ϕ j ϕ i dx, f i = (cid:90) Ω u ϕ i dx for i, j = 1 , . . . , N h . The solution u = ( u j ) ∈ R N h of (16) contains the unknowncoefficients of the finite element solution u (cid:37)h = N h (cid:88) j =1 u j ϕ j ∈ V h of (9). Robust and fast iterative solvers for (16) will be considered in the nextsection. In this section, we present two precondioners that lead to robust and efficient solversfor the system (16) of finite element equations when applying a preconditionedconjugate gradient algorithm. The BDDC preconditioner is especially suited forparallel computers.
In constrast to geometric multigrid methods, e.g, [10], usually relying on underlyinghierarchical meshes, algebraic multigrid methods [4, 5, 27, 34] are more flexible withrespect to complex geometries, adaptive mesh refinements, and so on. We refer to[9] for a comparison between these two methods. In this work, the particular AMGpreconditioner developed in [13] in combination with the conjugate gradient (CG)method is adopted to solve the discrete system. This AMG method was originallydeveloped for second-order elliptic problems. It has been extended to act as arobust preconditioner to systems arising from the discretization of coupled vectorfield problems, e.g., to fluid and elasticity problems in fluid-structure interaction5olvers [17, 35]. More recently, in [31], by choosing a proper smoother [26], it wasalso used as a preconditioner for the non-symmetric and positive definite systemthat arises from the Petrov-Galerkin space-time finite element discretization of theheat equation [30].In this AMG method, the coarsening strategy is based on a simple red-blackcolouring algorithm [13, Algorithm 2]: Pick up an unmarked degree of freedom(Dof), and mark it black (coarse); all neighboring Dofs are marked in red (fine);repeat this procedure until all Dofs are visited. Afterwards, the fine Dofs (red)are interpolated by an average over their neighboring coarse Dofs (black), whichdefines a linear interpolation operator P . The coarse grid operator is constructedby Galerkin’s method, i.e., A c = P (cid:62) AP . Assume that we apply m pre- and post-smoothing steps, the iteration operator for the two grid method is given by M amg := S m ( I − P A − c P (cid:62) A ) S m , where S is the iteration matrix for a smoothing step. For such a symmetric andpositive definite system, one sweep of the symmetric Gauss–Seidel method is used,that is, one forward on the down cycle together with one backward on the up cycle.Then, we may formulate the two-grid preconditioner P AMG as follows (with m = 1): P − AMG := S P A − c P (cid:62) S . (17)Note that, in the coarse grid correction step, one may apply AMG recursively; see[4, 5, 27], i.e., we replace A − c by AMG iterations starting with a zero initial guess.In connection with optimal control problems as considered in this paper, we areespecially interested in small (cid:37) ∈ (0 , (cid:37) →
0, the mass matrix M becomesdominant in A . For the mass matrix, the condition number is O (1), which evenmakes the system easier to solve as observed from our numerical tests. For the (two-level) BDDC preconditioner, we first make a reordering of the system(16), i.e., ˜ A ˜ u = ˜ f , (18)where˜ A := (cid:20) ˜ A II ˜ A IC ˜ A CI ˜ A CC (cid:21) , ˜ u := (cid:20) ˜ u I ˜ u C (cid:21) , ˜ f := (cid:20) ˜ f I ˜ f C (cid:21) , and ˜ A II = diag (cid:104) ˜ A II , ..., ˜ A pII (cid:105) . Here, p represents the number of polyhedral subdomains Ω i , as well as the number ofcores. These subdomains Ω i are obtained from the graph partitioning tool [12] anda non-overlapping domain decomposition of Ω [32]. Following common notations,the degrees of freedom in (18) are respectively decomposed into internal ( I ) andinterface ( C ) Dofs. Then, we arrive at the following Schur complement systemwhich is defined on the interface Γ C : S C ˜ u C = g, (19)where S C := ˜ A CC − ˜ A CI ˜ A − II ˜ A IC , g := ˜ f C − ˜ A CI ˜ A − II ˜ f I . Now, following similar notations as in [6, 11, 20, 21], the BDDC preconditioner P BDDC for (19) is formed as P − BDDC = R (cid:62) C ( T sub + T ) R C , R C represents the direct sum of restriction operators R iC thatmap the global interface vector to its components on a local interface Γ i := ∂ Ω i ∩ Γ C with a proper scaling. Furthermore, the coarse level correction operator T is definedby T = Φ(Φ (cid:62) S C Φ) − Φ (cid:62) , (20)where Φ = (cid:2) (Φ ) (cid:62) , . . . , (Φ N ) (cid:62) (cid:3) (cid:62) is the matrix of the coarse level basis functions.Each basis function matrix Φ i on the subdomain interface Γ i is computed by solvingthe augmented system: (cid:20) S iC (cid:0) C i (cid:1) (cid:62) C i (cid:21) (cid:20) Φ i Λ i (cid:21) = (cid:20) R i Π (cid:21) , (21)where S iC is nothing but the local Schur complement, C i denotes the given primalconstraints of the subdomain Ω i , and each column of Λ i contains the vector ofLagrange multipliers. The number of columns of each Φ i is equal to the numberof global coarse level degrees of freedom, usually living on the subdomain corners,and/or interface edges, and/or faces. Moreover, the restriction operator R i Π mapsthe global interface vector in the continuous primal variable space on the coarselevel to its component on Γ i . We mention that these corners/edges/faces on thecoarse space may be characterized as objects from a pure algebraic manner; see,e.g., [1]. We have recently used this abstract method in [18] for constructing BDDCpreconditioners in solving the linear system that arises from the space-time finiteelement discretization of parabolic equations [15].The subdomain correction operator T sub is then defined as T sub = N (cid:88) i =1 (cid:2) ( R iC ) (cid:62) (cid:3) (cid:20) S iC (cid:0) C i (cid:1) (cid:62) C i (cid:21) − (cid:20) R iC (cid:21) , with vanishing primal variables on all the coarse levels. Here, the restriction oper-ator R iC maps global interface vectors to their components on Γ i . Note that thistwo-level method can be extended to a multi-level method when the coarse problembecomes too large; see, e.g., [2, 36]. This means that we may apply the BDDCmethod to the inversion of Φ (cid:62) S C Φ in the coarse problem (20) recursively.
In all numerical examples presented in this section, we consider the unit cubeΩ = (0 , as computational domain. In the first example, we choose a smoothtarget, whereas a discontinuous target is used in the second example. The firstexample is covered by Theorem 1, the second not. Finally, we consider an examplewith a smooth target that violates the homogeneous Dirichlet boundary conditions.The system (16) of finite element equations is solved by a preconditioned conjugategradient (PCG) method with both AMG and BDDC preconditioners. The PCGiteration is stopped when the relative residual error of the preconditioned systemreaches ε = 10 − . In the AMG method, we have applied one V -cycle multigridpreconditioner with one pre- and post-smoothing step, respectively. We run theAMG preconditioned CG on the shared memory supercomputer MACH-2 . TheBDDC preconditioned CG is performed on the high-performance distributed mem-ory computing cluster RADON-1 . .1 Smooth target In the first example, we consider the smooth function¯ u ( x , x , x ) = − (cid:37) ∆ u (cid:37) + u (cid:37) = sin( πx ) sin( πx ) sin( πx )as target, which was computed from the manufactured solution u (cid:37) ( x , x , x ) = (3 (cid:37)π + 1) − sin( πx ) sin( πx ) sin( πx )of the reduced optimality equation (5). Then the adjoint state p (cid:37) ( x , x , x ) = − u (cid:37) ( x , x , x ) (cid:37) (3 (cid:37)π + 1) = − (cid:37) (3 (cid:37)π + 1) sin( πx ) sin( πx ) sin( πx )results from the gradient equation, and the optimal control z (cid:37) ( x , x , x ) = − ∆ u (cid:37) ( x , x , x ) = 3 π (cid:37)π + 1 sin( πx ) sin( πx ) sin( πx )from the state equation. To check the convergence in the L (Ω) and H (Ω) norms,we first compute the discretization errors (cid:107) u (cid:37) − u (cid:37)h (cid:107) L (Ω) and (cid:107)∇ ( u (cid:37) − u (cid:37)h ) (cid:107) L (Ω) for finer and finer mesh sizes h and 6 different values of the regularization parameter (cid:37) ; see Table 1 and Table 2, respectively. Since the given solution u (cid:37) is smooth, weobserve optimal convergence rates with respect to the mesh size h . Furthermore,the convergence does not deteriorate with respect to the regularization parameter (cid:37) .In Table 3, we numerically study the convergence of the approximation u (cid:37)h to thetarget ¯ u in the L (Ω) norm for decreasing (cid:37) by choosing h = (cid:37) / as predicted byTheorem 1. From Table 3, we clearly see a second-order convergence. We furtherobserve a second-order convergence in the H − (Ω) norm; see also [22, Table 5]. h (cid:37) eoc 10 − eoc 10 − eoc2 − . − . − . − − . − .
82 8 . − .
08 8 . − . − . − .
98 2 . − .
03 2 . − . − . − .
00 5 . − .
99 5 . − . − . − .
00 1 . − .
98 1 . − . − . − .
00 3 . − .
98 3 . − . − . − .
00 8 . − .
99 8 . − . h (cid:37) − eoc 10 − eoc 10 − eoc2 − . − . − . − − . − .
19 8 . − .
19 8 . − . − . − .
08 2 . − .
08 2 . − . − . − .
04 4 . − .
04 4 . − . − . − .
02 1 . − .
03 1 . − . − . − .
01 3 . − .
01 3 . − . − . − .
00 7 . − .
01 7 . − . (cid:107) u (cid:37) − u (cid:37)h (cid:107) L (Ω) with respect to h (columns), and fixed,but different (cid:37) (rows).Table 4 provides the number of AMG preconditioned CG iterations, and the totalcoarsening and solving time measured in second (s) with respect to h and (cid:37) . Fromthese numerical results, it is clear to see the relative robustness of the AMG pre-conditioner with respect to the mesh refinement and the decreasing regularization8 (cid:37) eoc 10 − eoc 10 − eoc2 − . − . − . − − . − .
92 4 . − .
00 5 . − . − . − .
99 2 . − .
02 2 . − . − . − .
00 1 . − .
01 1 . − . − . − .
00 5 . − .
00 6 . − . − . − .
00 2 . − .
00 3 . − . − . − .
00 1 . − .
00 1 . − . h (cid:37) − eoc 10 − eoc 10 − eoc2 − . − . − . − − . − .
09 5 . − .
09 5 . − . − . − .
04 2 . − .
04 2 . − . − . − .
01 1 . − .
01 1 . − . − . − .
01 6 . − .
00 6 . − . − . − .
00 3 . − .
00 3 . − . − . − .
00 1 . − .
00 1 . − . (cid:107)∇ ( u (cid:37) − u (cid:37)h ) (cid:107) L (Ω) with respect to h (columns), andfixed, but different (cid:37) (rows). (cid:37) (cid:107) ¯ u − u (cid:37)h (cid:107) H − (Ω) eoc (cid:107) ¯ u − u (cid:37)h (cid:107) L (Ω) eoc10 . − . − − . − .
22 6 . − . − . − .
03 6 . − . − . − .
80 1 . − . − . − .
98 1 . − . (cid:107) ¯ u − u (cid:37)h (cid:107) H − (Ω) and (cid:107) ¯ u − u (cid:37)h (cid:107) L (Ω) for ¯ u ∈ H (Ω) ∩ H (Ω).parameter (cid:37) . The computational time scales well with respect to the number ofunknowns. In Table 5, we study the robustness with respect to (cid:37) and parallel per-formance (strong scaling) of the BDDC preconditioned CG solver. So, we fix thetotal number of unknowns to Dof s = 2 , ,
689 that is related to h = 1 / p , and the varying regularization parameter (cid:37) . Further, weobserve perfect strong scalability (checking each row). Over-scaling from 32 to 64cores may result from heterogeneity of the cluster, larger memory consumption forsmall number of cores, or the shared node resources with other users. In the second example, we consider a discontinuous target function¯ u ( x ) = (cid:40) x ∈ (0 . , . x ∈ Ω \ (0 . , . , (22)that is similar to the first example in [22]. In this case, we solve the equationon adaptively refined meshes that are driven by a residual based error indicator[33]. This leads to a lot of local refinements near the interface where the targetis discontinuous. The energy regularization leads to a control that concentrates9 h
18 116 132 164 1128 1256 − − − − − − Table 4: Example 4.1: Number of AMG preconditioned CG iterations, total coars-ening and solving time measured in second (s), with respect to h and (cid:37) . (cid:37) p
32 64 128 25610
34 (46.43 s) 34 (16.04 s) 35 (5.92 s) 36 (2.98 s)10 −
32 (43.84 s) 32 (15.20 s) 33 (5.58 s) 34 (2.70 s)10 −
27 (37.20 s) 27 (13.12 s) 28 (4.85 s) 30 (2.39 s)10 −
20 (28.01 s) 20 (9.77 s) 22 (3.83 s) 22 (1.87 s)10 −
20 (28.15 s) 20 (9.75 s) 21 (3.65 s) 21 (1.76 s)10 −
20 (28.02 s) 20 (9.69 s) 21 (3.63 s) 21 (1.74 s)10 −
20 (28.01 s) 20 (9.71 s) 21 (3.69 s) 21 (1.67 s)Table 5: Example 4.1: Number of iterations, and time measured in second (s) forthe BDDC preconditioned CG solver, where h = 1 /
128 and
Dof s = 2 , , u (cid:37)h , and the control z (cid:37)h in Fig. 1, which are comparable to the results in [22,Fig. 3]. Moreover, we provide the convergence of the approximation u (cid:37)h to thetarget ¯ u with respect to the regularization parameter (cid:37) in Table 6, i.e., an order √ . H − (Ω) and √ . L (Ω). In this case, since¯ u ∈ H / − ε (Ω), we obtain reduced convergence of the finite element approxima-tion to the target; see also [22, Table 1]. Due to adaptive refinements for differentregularization parameters (cid:37) = 10 k , k = 0 , − , . . . , −
6, the final meshes may con-tain different number of degrees of freedom. In Fig. 2, we visualize the numberof AMG preconditioned CG iterations (left plot) and the computational time inseconds (right plot) for different values of (cid:37) during the adaptive refinement proce-dure. We observe that the AMG preconditioned CG iterations stay in a small rangebetween 5 and 25; the computational time scales well with respect to the numberof unknowns. To compare the performance of BDDC preconditioners with respectto an adaptive refinement, we select similar numbers of Dofs for all regulariza-tion parameters, i.e., (cid:37) = 10 ), (cid:37) = 10 − ), (cid:37) = 10 − ), (cid:37) = 10 − ), (cid:37) = 10 − ), (cid:37) = 10 − ), and (cid:37) = 10 − ). InTable 7, we present the number of BDDC preconditioned CG iterations, and thecorresponding solving time, measured in second (s), with respect to the number ofsubdomains p and the regularization parameter (cid:37) . For such adaptive meshes, weobserve relatively good performance in both the number of BDDC preconditionedCG iterations and the computational time, as well as good strong scalability withrespect to the number of subdomains p . In the third example, we use the smooth target¯ u ( x , x , x ) = 1 + sin( πx ) sin( πx ) sin( πx )10igure 1: Example 4.2: Adaptive meshes at the 21th refinement step with 1 , , u (cid:37)h in Ω and on the cutting plane x = 0 . , . , .
5] and [1 , . , .
5] (down), (cid:37) = 10 − . (cid:37) (cid:107) ¯ u − u (cid:37)h (cid:107) H − (Ω) eoc (cid:107) ¯ u − u (cid:37)h (cid:107) L (Ω) eoc10 . − . − − . − .
20 9 . − . − . − .
85 3 . − . − . − .
33 1 . − . − . − .
45 3 . − . − . − .
48 1 . − . − . −
10 1 .
50 3 . − . (cid:107) ¯ u − u (cid:37)h (cid:107) H − (Ω) and (cid:107) ¯ u − u (cid:37)h (cid:107) L (Ω) of the approxi-mations u (cid:37)h to the target ¯ u ∈ H / − ε (Ω).considered in [22]. This target violates the homogeneous Dirichlet boundary con-ditions. As in the previous example, we solve the equation on adaptively refinedmeshes, which results in many local refinements near the boundary, where the con-trol concentrates. For an illustration, we visualize the adaptive mesh, the state u (cid:37)h ,and the control z (cid:37)h in Fig. 3. Further, we observe the convergence order √ . H − (Ω), and √ . L (Ω) with respect to (cid:37) ; see Table 8.Since ¯ u ∈ C ∞ (Ω), but ¯ u / ∈ H (Ω), we only expect a reduced order of convergencefor this example; see also [22, Table 7]. In Fig. 4, for each (cid:37) , we visualize the num-ber of AMG preconditioned CG iterations (left plot) and the computational timein seconds (right plot) during the adaptive refinement levels. The AMG precondi-tioned CG iterations are in a small range between 1 and 20. The computationaltime scales well with respect to the number of unknowns. To see the performance ofthe BDDC preconditioner with respect to an adaptive refinement, we select similarnumbers of Dofs for all regularization parameters, i.e., (cid:37) = 10 ), (cid:37) = 10 − ), (cid:37) = 10 − ), (cid:37) = 10 − ), (cid:37) = 10 − ), (cid:37) = 10 − ), and11igure 2: Example 4.2: Number of AMG preconditioned CG iterations (left)and computational time in seconds (right) with respect to different (cid:37) ∈{ , − , ... − } . (cid:37) p
32 64 128 256
37 (60.77 s) 46 (24.34 s) 42 (8.20 s) 52 (5.02 s) 2 , , −
42 (66.85 s) 45 (24.85 s) 46 (9.28 s) 43 (4.64 s) 2 , , −
34 (60.06 s) 38 (23.86 s) 42 (9.75 s) 43 (4.79 s) 2 , , −
35 (40.92 s) 39 (15.35 s) 35 (5.12 s) 40 (3.23 s) 1 , , −
33 (16.79 s) 31 (6.30 s) 39 (3.17 s) 36 (1.53 s) 1 , , −
33 (17.63 s) 33 (7.65 s) 33 (3.28 s) 38 (1.88 s) 1 , , −
36 (18.30 s) 37 (8.42 s) 35 (3.31 s) 44 (1.93 s) 1 , , p (= number ofcores), regularization parameter (cid:37) , and different adaptive meshes ( (cid:37) = 10 − ). Table 9 presents the number of BDDC precondi-tioned CG iterations, and solving times measured in second (s) with respect to thenumber of subdomains p and the regularization parameter (cid:37) . We again observe arelatively good performance in both the number of BDDC preconditioned CG iter-ations and the computational time, as well as good strong scalability with respectto the number of subdomains p . (cid:37) (cid:107) ¯ u − u (cid:37)h (cid:107) H − (Ω) eoc (cid:107) ¯ u − u (cid:37)h (cid:107) L (Ω) eoc10 . − . − − . − .
21 1 . − . − . − .
93 3 . − . − . − .
46 9 . − . − . − .
51 3 . − . − . − .
50 9 . − . (cid:107) ¯ u − u (cid:37)h (cid:107) H − (Ω) and (cid:107) ¯ u − u (cid:37)h (cid:107) L (Ω) of the approxi-mations u (cid:37)h to the target ¯ u ∈ C ∞ (Ω, ¯ u / ∈ H (Ω).12igure 3: Example 4.3: Adaptive meshes at the 12th refinement step with 3 , , u (cid:37)h in Ω and z (cid:37)h in Ω (up); u (cid:37)h and z (cid:37)h along the line between [0 , . , . , . , .
5] (down), (cid:37) = 10 − . Using estimates for the error between the exact solution u (cid:37) of the optimal controlproblem (1) and the target u in terms of the regularization parameter (cid:37) , derivedin [22], we have investigated the error between the finite element solution u (cid:37)h andthe target u in terms of (cid:37) and h . Furthermore, we have studied AMG and BDDCpreconditioned CG solvers for (adaptive) finite element equations arising from theelliptic optimal control problem with energy regularization. Both preconditionershave shown good performance with respect to (adaptive) mesh refinements and a(decreasing) regularization parameter. Moreover, the BDDC preconditioner hasshown its strong scalability with respect to the number of subdomains on a dis-tributed memory computer.While in this paper we have considered a distributed control problem subject tothe Poisson equation, a related analysis can be done in the case of a parabolic heatequation. These results will be published elsewhere, but see [16] for a space-timediscretization approach. Acknowledgments
The third author was partly supported by the Austrian Science Fund (FWF) viathe grant NFN S117-03.
References [1] S. Badia, A. Martin, and J. Principe. Implementation and scalability analysisof balancing domain decomposition methods.
Arch. Computat. Methods Eng. ,20:239–262, 2013. 13igure 4: Example 4.3: Number of AMG preconditioned CG iterations (left) andcomputational time in seconds (right) with respect to different (cid:37) ∈ { , − , ... − } . (cid:37) p
32 64 128 256
33 (45.61 s) 32 (15.37 s) 32 (5.62 s) 33 (2.78 s) 2 , , −
30 (41.51 s) 31 (15.16 s) 32 (5.71 s) 32 (2.84 s) 2 , , −
30 (42.13 s) 31 (15.46 s) 30 (5.14 s) 32 (2.59 s) 2 , , −
29 (32.46 s) 34 (14.30 s) 32 (7.20 s) 39 (3.27 s) 2 , , −
29 (15.87 s) 31 (6.98 s) 31 (2.94 s) 29 (1.40 s) 1 , , −
30 (21.03 s) 30 (9.22 s) 33 (4.34 s) 37 (2.73 s) 3 , , −
31 (22.57 s) 33 (10.06 s) 31 (4.05 s) 38 (2.57 s) 3 , , p (= number ofcores), regularization parameter (cid:37) , and different adaptive meshes ( SIAM J. Sci. Comput. , 38(1):C22–C52, 2016.[3] D. Braess.
Finite Elements: Theory, Fast Solvers, and Applications in SolidMechanics . Cambridge University Press, Cambridge, 2007.[4] A. Brandt, S. McCormick, and J. Ruge. Algebraic multigrid (AMG) for sparsematrix equations. In D. Evans, editor,
Sparsity and its Applications , pages257–284. Cambridge University Press, Cambridge, 1985.[5] W. Briggs, V. Henson, and S. McCormick.
A Multigrid Tutorial . SIAM,Philadelphia, second edition, 2000.[6] C. Dohrmann. A preconditioner for substructuring based on constrained energyminimization.
SIAM J. Sci. Comput. , 25(1):246–258, 2003.[7] H. Goering, A. Felgenhauer, G. Lube, H. Roos, and L. Tobiska.
SingularlyPerturbed Differential Equations . Akademie Verlag, Berlin, 1983.[8] H. Goering, M. Stynes, and L. Tobiska.
Robust Numerical Methods for Singu-laryly Perturbed Differential Equations . Springer Verlag, Berlin, 2008.[9] G. Haase and U. Langer. Multigrid methods: from geometrical to algebraicversions. In A. Bourlioux, M. Gander, and G. Sabidussi, editors,
ModernMethods in Scientific Computing and Applications , pages 103–153. SpringerNetherlands, Dordrecht, 2002. 1410] W. Hackbusch.
Multi-Grid Methods and Applications . Springer, Heidelberg,1985.[11] L. Jing and O. B. Widlund. FETI-DP, BDDC, and block Cholesky methods.
Int. J. Numer. Meth. Engng. , 66(2):250–271, 2006.[12] G. Karypis and V. Kumar. A fast and high quality multilevel scheme forpartitioning irregular graphs.
SIAM J. Sci. Comput. , 20(1):359–392, 1998.[13] F. Kickinger. Algebraic multi-grid for discrete elliptic second-order problems.In W. Hackbusch and G. Wittum, editors,
Multigrid Methods V , pages 157–172,Berlin, Heidelberg, 1998. Springer Verlag.[14] J. Kraus and M. Wolfmayr. On the robustness and optimality of algebraicmultilevel methods for reaction-diffusion type problems.
Comput. Visual Sci. ,16:15–32, 2013.[15] U. Langer, M. Neum¨uller, and A. Schafelner. Space-time finite element meth-ods for parabolic evolution problems with variable coefficients. In T. Apel,U. Langer, A. Meyer, and O. Steinbach, editors,
Advanced Finite ElementMethods with Applications: Selected Papers from the 30th Chemnitz FiniteElement Symposium 2017 , pages 247–275. Springer International Publishing,Cham, 2019.[16] U. Langer, O. Steinbach, F. Tr¨oltzsch, and H. Yang. Space-time finite elementdiscretization of parabolic optimal control problems with energy regularization.
SIAM J. Numer. Anal. , 2020. accepted.[17] U. Langer and H. Yang. Robust and efficient monolithic fluid-structure-interaction solvers.
Int. J. Numer. Meth. Engng. , 108(4):303–325, 2016.[18] U. Langer and H. Yang. BDDC preconditioners for a space-time finite elementdiscretization of parabolic problems. In R. Haynes, S. MacLachlan, X. Cai,L. Halpern, H. Kim, A. Klawonn, and O. Widlund, editors,
Domain Decompo-sition Methods in Science and Engineering XXV , pages 367–374, Cham, 2020.Springer.[19] S. MacLachlan and N. Madden. Robust solution of singularly perturbed prob-lems using multigrid methods.
SIAM J. Sci. Comput. , 35(5):A2225–A2254,2013.[20] J. Mandel and C. Dohrmann. Convergence of a balancing domain decompo-sition by constraints and energy minimization.
Numer. Linear Algebra Appl. ,10(7):639–659, 2003.[21] J. Mandel, C. Dohrmann, and R. Tezaur. An algebraic theory for primal anddual substructuring methods by constraints.
Appl. Numer. Math. , 54(2):167–193, 2005.[22] M. Neum¨uller and O. Steinbach. Regularization error estimates for distributedcontrol problems in energy spaces.
Math. Meth. Appl. Sci. , 2020. publishedonline.[23] T. A. Nhan, S. MacLachlan, and N. Madden. Boundary layer precondition-ers for finite-element discretizations of singularly perturbed reaction-diffusionproblems.
Numer. Algor. , 79:281–310, 2018.1524] T. A. Nhan and N. Madden. An analysis of diagonal and incomplete Choleskypreconditioners for singularly perturbed problems on layer-adapted meshes.
J.Appl. Math. Computing , 65:245–272, 2021.[25] M. Olshanskii and A. Reusken. On the convergence of a multigrid method forlinear reaction-diffusion problems.
Computing , 65:193–202, 2000.[26] C. Popa. Algebraic multigrid smoothing property of kaczmarz’s relaxation forgeneral rectangular linear systems.
Electron. Trans. Numer. Anal. , 29:150–162,2007-2008.[27] J. Ruge and K. St¨uben. Algebraic multigrid. In S. McCormick, editor,
MultigridMethods , chapter 4, pages 73–130. SIAM, Philadelphia, 1987.[28] A. H. Schatz and L. B. Wahlbin. On the finite element method for singu-larly perturbed reaction-diffusion problems in two and one dimensions.
Math.Comp. , 40(161):47–89, 1983.[29] O. Steinbach.
Numerical Approximation Methods for Elliptic Boundary ValueProblems: Finite and Boundary Elements . Springer, 2008.[30] O. Steinbach. Space-time finite element methods for parabolic problems.
Com-put. Methods Appl. Math. , 15:551–566, 2015.[31] O. Steinbach and H. Yang. Comparison of algebraic multigrid methods for anadaptive space-time finite-element discretization of the heat equation in 3d and4d.
Numer. Linear Algebra Appl. , 25(3):e2143, 2018.[32] A. Toselli and O. B. Widlund.
Domain Decomposition Methods — Algorithmsand Theory . Springer, Berlin, Heidelberg, 2005.[33] R. Verf¨urth.
A posteriori error esimtation techniques for finite element meth-ods . Oxford University Press, Oxford, 2013.[34] J. Xu and L. Zikatanov. Algebraic multigrid methods.
Acta Numer. ,26:591–721, 2017.[35] H. Yang and W. Zulehner. Numerical simulation of fluid-structure interactionproblems on hybrid meshes with algebraic multigrid methods.
J. Comput.Appl. Math. , 235(18):5367–5379, 2011.[36] S. Zampini and X. Tu. Multilevel balancing domain decomposition by con-straints deluxe algorithms with adaptive coarse spaces for flow in porous media.