4. Computation of steps for optimal control problems 87
4.5. Approximation of operators
4.5.2. Approximation of the stiffness matrix
Regarding the discretizations of the differential operators, occurring in the con-straint and the adjoint equation, hierarchical multigrid solvers and preconditioners are among the most promising methods for adaptive algorithms. The two main reasons are the h-independent convergence rate and linear complexity with respect to the number of unknowns. Both properties are not only provable, but also can be verified for reasonable implementations of multigrid algorithms [260].
These mainly exploit the observation that the classification of high- and low-frequency error components in finite element (FE)-computations is dependent on the resolu-tion of the spatial domain. High-frequency errors on a fine grid are not captured on coarser grids. In contrast low-frequency errors on fine grids may appear high-frequency on coarser ones, whereas smooth error components essentially can be represented on significantly coarser grids.
In order to exploit this insight, multigrid solvers use cheap smoothers, such as damped Jacobi- or Gauss-Seidel-iterations [79], to eliminate high-frequency error components. Repeated application on grid with different spatial resolution then admits to significantly reduce the oscillatory components of the algebraic error.
Eventually on the coarsest grid the remaining error can be eliminated by a direct solver.
To express this formally we consider the operator equation Au = b in X∗ and a sequence of hierarchical grids, resp. nested FE-spaces
S0 ⊂ · · · ⊂Sj ⊂X
with corresponding projected differential operators Ak = A|Sk. We denote the pro-jection operators by Ik−1k : Sk−1 → Sk and the corresponding restriction operators byIkk−1 : Sk →Sk−1. The essential ingredients of the multigrid method are captured in Alg. 4.8.
For adequate choices of the algorithmic parameters, ν1 applications of the smoother eliminate the high-frequency error components in Sk. Then, the remaining error is projected to a coarser grid and eliminated there. Assuming that the remaining error is “smooth” with respect to Sk we expect relatively small loss due to restriction and interpolation operators. After the computation of the coarse grid correction one
Chapter 4 Computation of steps for optimal control problems
Require: given smoothing parametersν1, ν2, grid levelk, operatorAkand right hand side bk, initial value uk.
1: mgCycle:
2: if k = 0 do
3: use direct solver to solveA0δu0 =b0
4: u0 ←u0+δu0
5: else do
6: apply ν1 steps of a smoother to the system Akuk=bk
7: compute the restriction of the residual to the coarse spacerk−1 ←Ikk−1rk= Ikk−1(bk−Akuk)
8: compute Ak−1 =Ik−1k AkIkk−1, set δuk−1 = 0
9: apply mgCycle(ν1, ν2, k−1, Ak−1, rk−1, δuk−1)
10: correct the fine grid solutionuk ←uk+Ik−1k δuk−1
11: apply ν2 steps of a smoother to the system Akuk=bk Algorithm 4.8: Two-grid correction scheme.
further applies a smoother ν2 times to eliminate possibly remaining high-frequency error contributions. On the coarsest grid a direct solver can be employed to compute a highly accurate coarse grid correction. Alg. 4.9 is essentially the V-cycle multigrid method. For the realization of a linear solver, that can satisfy prescribed accu-racy requirements, we repeatedly apply the above algorithm on the defect equation.
See Briggs et al. [45], Trottenberg et al. [260] for more details and the W-cycle as well as the full multigrid (FMG) scheme. Same as the Chebyshev semi-iteration the multigrid method provides a linear preconditioner if employed with a fixed number of iterations.
Require: given smoothing parametersν1, ν2, operator A and right hand side b.
1: while convergence test failed do
2: setr←b−Au and δu= 0
3: apply mgCycle(ν1, ν2,0, A, r, δu)
4: set u←u+δu
Algorithm 4.9: V-cycle multigrid method.
On optimized academic examples multigrid solvers are extremely fast. This is il-lustrated in Tab. 4.2 where the iteration numbers are given for a simple example of linear heat transfer on the unit cube, with right hand side f = 1. On real-world problems, where already small and coarse geometries may contain larger numbers of degrees of freedoms in a linear finite element space, their performance is less out-standing, but still impressive. This is illustrated in Tab. 4.3 for equations of linear elasticity with the same right hand side. On the left of this table the iterations
112
4.5 Approximation of operators
numbers and computational time on the unit cube are given. On the right we give the corresponding numbers on the geometry used in Sec. 6.2.2.1. Thus, we will have to apply several V-cycles to get a reasonable preconditioner for our KKT-systems.
If used to approximate the differential operators in the constraint preconditioner (4.2.4) we even need to actually solve the constraint and adjoint equation. For these reason the repeated application of multigrid solvers and preconditioners for problems in constraint preconditioners for PPCG methods still is quite expensive.
Laplace (2D, unit cube) Laplace (3D, unit cube)
dof #iter. time dof #iter. time
545 7 2.2 ms 369 8 2.6 ms
8 321 7 32 ms 2 465 9 21 ms
33 025 7 0.13 s 17 985 10 0.18 s
131 585 7 0.59 s 137 345 10 1.7 s
2 099 201 7 11 s 1 073 409 11 16 s
Table 4.2.: Computation times for simple test problems (rel. acc.: 10−9, smoothing steps: 10).
Linearized elasticity (unit cube) Linearized elasticity (real world geometry)
dof #iter. time dof #iter. time
53 955 43 3.6 s 213 687 25 18 s
412 035 49 35 s 1 523 604 48 4.4 min
3 220 227 53 5 min 11 461 518 72 49 min
Table 4.3.: Computation times for simple test problems from linearized elasticity (both 3D, rel. acc.: 10−9, smoothing steps: 10, resp. 20 for the second problem)
We note that in the context of pressure-type boundary conditions the lower left and upper right block in the KKT-system not only contain the differential operator Wϕϕ, but also non-symmetric contributions from the Piola-transformed pressure-type boundary conditions gcof(∇ϕ)n. Moreover, due to the polyconvexity of the stored energy function, the corresponding differential operator may not be elliptic.
Consequently, in our setting, the application of multigrid solvers and preconditioners is not backed by a solid theoretical basis. However, in our computations they seemed to work well.
Chapter 4 Computation of steps for optimal control problems
4.6. Summary
Let us summarize the step computations within the different settings. For the computations of δn, px and δs we can always assume positive definiteness ofM on kerc0(x), and thus unique solvability of the corresponding system. For moderately sized problems, or a low dimensional control space, the solution can be found by direct elimination methods. Otherwise a PPCG method, i.e. a conjugate gradient method combined with a constraint preconditioner, can be used.
The situation is different for tangential steps δt. As Lxx is in general indefinite we have to use one of the modifications from Sec. 4.3, the HCG method, to compute descent directions for the cost functional. The restriction to kerc0(x) is again in-corporated with the help of a constraint preconditioner. For problems of moderate size or low dimensional control space we can reuse the direct factorization which was computed for the determination of the normal step as preconditioner. If this approach is not admissible we will use the same block triangular preconditioner
Qsc =
as in the computation of the normal step. To efficiently evaluate this preconditioner we have to efficiently solve the state, variational and adjoint equation. For this we replace the inversion of bothA and AT by a multigrid solver. Since the constraint preconditioner has to project onto kerc0(x0), and in the absence of further analysis, it is necessary to solve the arising systems A0(y0)δny =rp+Bδnu to high accuracy.
Relaxing this condition on Qsc is subject to ongoing work. In contrast Mu−1 can be replaced by a fixed number of Chebyshev semi-iterations [110, 124, 274], which needs not to be overly accurate. The required spectral bounds for the Chebyshev semi-iterations for the preconditioned matrixQjacMu, whereQjacrepresents one step of the Jacobi iteration, are taken from [273].
For the purpose of error estimation and adaptive mesh refinement in an affine co-variant setting, a hierarchical estimator for the error of the primal variables in the Lagrange-Newton step was proposed. Both the block triangular preconditioner and the error indicator require the solution of equations involving the mass matrix resp. the differential operators of the state and adjoint equation. For the first again the Chebyshev semi-iteration is employed, whereas the differential operators are treated with multigrid preconditioners (25 V-cycles and 20 pre- and post-smoothing steps).
114