Parameter‐robust preconditioning for unsteady Stokes control problems

We propose a saddle‐point preconditioner for an optimization problem constrained by the time‐dependent Stokes equations, discretized using the backward Euler method in time. The key ingredients are an inner iteration for the (1, 1)‐block, accelerated by a known preconditioner for the heat control problem, and an approximation of the Schur complement involving a commutator argument applied to a block matrix. Numerical results demonstrate the efficacy and robustness of this approach.


Problem statement and resulting linear system
The rapid iterative solution of time-dependent partial (PDE) [9] or fractional (FDE) [10] differential equation constrained optimization problems is a key challenge in numerical linear algebra. Here, we consider the unsteady Stokes control problem: The state variables (⃗ v, p) denote velocity and pressure in space-time domain Q := Ω×(0, T ), with ⃗ v equipped with boundary data on ∂Ω and an initial condition, ⃗ u and ⃗ v d denote the control variable and desired state, and β is a regularization parameter.
We discretize the PDE constraints above and the corresponding adjoint equations (with adjoint variables ⃗ λ, µ to ⃗ v, p), using the gradient equation β⃗ u − ⃗ λ = ⃗ 0 to eliminate ⃗ u. Applying backward Euler in time leads to the 'saddle-point type' system: The vectors b 1 , b 2 , b 3 arise from the boundary data, v d corresponds to the desired state, and with K discretizing the (negative) vector Laplacian, M the identity operator in a vectorial sense, and B the (negative) divergence. We take a time-step size τ = T /n t , with n t time-steps. Each block matrix in (1) contains n t + 1 blocks, with the MATLAB notation blkdiag defining a block diagonal matrix. For the spatial discretization, we elect to use Taylor-Hood (Q 2 -Q 1 ) finite elements. We now present a preconditioner for (1), based on work in [6] for Navier-Stokes control problems.

Preconditioned iterative method for linear system
To solve the linear system (1) we employ a block-triangular preconditioner P which may be summarized as follows: which may be applied within a Krylov subspace method such as GMRES [11]. The matrix A is related to a heat control problem, solved over each spatial dimension, and A is known to be a potent and parameter-robust preconditioner for A [9]. Here, ϵ > 0 is a 'perturbation parameter' which guarantees invertibility of M v , and M = blkdiag (0, τ M, ..., τ M, √ ϵ τ M) en- To ensure a sufficiently accurate approximation of A within the Stokes control preconditioner P, we apply A −1 within an inner iterative method: either GMRES, or an Uzawa iteration (see [2]).
The matrix S approximates the (negative) Schur complement S = BA −1 B T of the matrix in (1). We apply a commutator argument [3,Ch. 9], which allows one to build approximations of the form BL −1 B T ≈ K p L −1 p M p , where L is a matrix derived on the velocity space, with K p , M p , L p analogous operators on the pressure space to K, M, L. Based on this we take See [6] for a complete description of this argument, which involves applying the commutator argument to each block of A.

Numerical results for two solution strategies
We test our preconditioner on a problem posed on the space-time domain (−1, 1) 2 × (0, 2), with spatial coordinates (x 1 , x 2 ), the boundary condition ⃗ f (⃗ x, t) = [min{t, 1}, 0] T on x 2 = 1 with ⃗ f (⃗ x, t) = 0 elsewhere, the initial condition ⃗ g(⃗ x) = ⃗ 0, and Tables 1-2 present results for the following two approaches, also reporting the degrees of freedom of each system solved: ▷ Table 1 shows results with flexible GMRES [12] as the outer method, and an inner GMRES solver for the matrix A. ▷ Table 2 shows results with flexible GMRES as the outer method, and an inner Uzawa method for the matrix A. We run all tests using MATLAB R2018b, on a 1.70GHz Intel quad-core i5 processor and 8 GB RAM on an Ubuntu 18.04.1 LTS operating system. We restart the outer GMRES method every 10 iterations, and solve to a relative residual-norm tolerance of 10 −6 . The flexible GMRES routine is based on that in [8], using MATLAB's inbuilt routine for the standard GMRES method. We take 5 inner GMRES or Uzawa iterations to approximate A, approximate M and M p with 20 steps of Chebyshev semi-iteration [4], each diagonal block of L + 1 √ β M and its transpose with 4 V-cycles of the algebraic multigrid routine AGMG [7], K p with 2 V-cycles of the HSL_MI20 solver [1] (with 2 pre-/post-smoothing steps), and take ϵ = 10 −3 . In Tables 1-2, spatial mesh level ℓ indicates a constant step of 2 −ℓ (2 1−ℓ ) in each spatial dimension for Q 2 (Q 1 ) basis functions.