A Newton multigrid framework for optimal control of fluid–structure interactions

In this paper we consider optimal control of nonlinear time-dependent fluid structure interactions. To determine a time-dependent control variable a BFGS algorithm is used, whereby gradient information is computed via a dual problem. To solve the resulting ill conditioned linear problems occurring in every time step of state and dual equation, we develop a highly efficient monolithic solver that is based on an approximated Newton scheme for the primal equation and a preconditioned Richardson iteration for the dual problem. The performance of the presented algorithms is tested for one 2d and one 3d example numerically.


Introduction
Fluid-structure interactions are part of various applications ranging from classical engineering problems like aeroelasticity or naval design to medical applications, e.g. the flow of blood in the heart or in blood vessels. More and more of these applications are regarded recently in combination with optimal control, shape-optimization, and parameter estimation. Especially in hemodynamical applications-in order to get a deeper understanding of the development of vascular diseases-patient specific properties have to be incorporated into the models. For example, in Bertoglio 1 3 et al. (2012,2013,2014), D'Elia et al. (2012), Lassila et al. (2013), Pant et al. (2014), Moireau et al. (2013) patient specific boundary conditions and vessel material parameters are determined to simulate arterial blood flow. Similar approaches using gradient information have been proposed in D' Elia et al. (2012), Bertagna et al. (2014), Perego et al. (2011) to estimate Young's modulus of an artery.
As computer tomography (CT) and magnetic resonance imaging (MRI) evolve rapidly, already very accurate measurements of the movement of the vessel wall are possible nowadays and even averaged flow profiles in blood vessels can be provided, see Asner et al. (2016), Bertoglio et al. (2014), Lamata et al. (2014). To incorporate the data in the vascular models, it is necessary to improve the available parameter estimation and optimal control algorithms for fluid-structure interaction applications, in particular since only few approaches in the literature take the sensitivity information of the full time-dependent nonlinear system into account. For example in Degroote et al. (2013), Martin et al. (2005), adjoint equations are derived for one-dimensional fluid-structure interaction configurations and in Richter and Wick (2013) for a stationary fluid-structure interaction problem. In contrast, the authors of Pant et al. (2014), Bertoglio et al. (2012Bertoglio et al. ( , 2014, Moireau et al. (2013) use a sequential reduced Kalman filter. Perego et al. (2011) compute sensitivity information to estimate the wall stiffness. To reduce the computational time, they solve in every time-point an optimal control problem. As the mesh motion is discretized via an explicit time-stepping scheme, no sensitivity information of the mesh motion equation has to be computed. Similar to the articles Bertoglio et al. (2012Bertoglio et al. ( , 2014, Moireau et al. (2013), the estimated parameters are updated in every time step and the forward simulation only runs once.
In this paper we are going to compute gradient information for the full timedependent nonlinear system for 3d applications. Thereby, the optimization algorithm takes the intrinsic property of fluid-structure interaction, transport over time, into account. In addition the here presented approach enables to regard tracking type functionals with observation at a singular time-point or on a specific time-interval. Furthermore a time-dependent parameter can be reconstructed. This would not be possible, if we would use a Kalman filter or would solve an optimization problem in every time step as in the literature cited above. The dual problem to compute sensitivities can be derived as in Failer (2017) or as in Failer and Wick (2018), where sensitivity information was used for a dual-weighted residual error estimator.
For various applications and a general overview on modeling and discretization techniques for fluid-structure interactions we refer to Bazilevs et al. (2013), Richter (2017). Mathematically, two challenges come together in fluid-structure interactions: First, fluid-structure interactions are free boundary value problems. The governing domains for the fluid-we will consider the incompressible Navier-Stokes equations-and the solid-we consider hyperelastic materials like the St. Venant Kirchhoff model-move and the motion is determined by the coupled dynamics, i.e. it is not known a priori. This geometric problem is treated by mapping onto a fixed domain, see for example Donea (1982), such that movement of the boundary is incorporated into the equation and we can derive the dual problem on the fixed reference domain. Second, the two problems that are coupled are of different type, the parabolic Navier-Stokes equation and the hyperbolic solid problem.
A Newton multigrid framework for optimal control of fluid-… On the common and moving interface, both systems are coupled by different conditions. This coupling gives rise to stability problems that can call for small time steps or many subiterations. Most prominently this problem shows itself in the so called added mass effect presented in Causin et al. (2005). The added mass effect is of particular relevance in hemodynamical applications, that are focus of this work, and calls according to Hron et al. (2010) for monolithic formulations and strongly coupled discretizations and solution techniques. This property is transmitted to the dual problem, such that we have to derive strongly coupled solution techniques for the dual problem. To compute dual information, we extend the Newton solver proposed for time-dependent fluid-structure interactions in Failer and Richter (2020) to the dual problem. Thereby, iterative solvers, preconditioned with geometric multigrid, can solve the resulting linear problems in every state and dual time step very robust and efficiently.
In Sect. 2 we present the optimal control problem, which is discretized in Sect. 3 in space and time. For the discretized system we derive optimality conditions. In Sect. 4 we discuss modifications for the Newton scheme presented first in Failer and Richter (2020) and extend the approach to the dual problem. Finally we test the proposed algorithm in Sect. 5 numerically to analyze the behavior of the Newton scheme. In addition we take a closer look on the convergence behavior of the iterative solvers.

Governing equations
Here, we present the optimal control problem of a tracking type functional subject to fluid structure interactions. We use a monolithic formulation for the fluid-structure interaction model coupling the incompressible Navier-Stokes equations and an hyperelastic solid, based on the St. Venant Kirchhoff material. For details we refer to Richter (2017). The here presented optimization approach can be directly extended to specific material laws used in hemodynamics. As control variable we chose exemplarily the mean pressure over time at the outflow boundary. In the following we restrict us to the control space Q = L 2 (I) , but the here presented optimization algorithm can as well be applied to determine material parameters (e.g. Q = ℝ n ) or space-distributed parameters (e.g. Q = L 2 (Ω) ) entering the fluid-or solid-problem.
On the d-dimensional domain, partitioned in reference configuration Ω = F ∪ I ∪ S , where F is the fluid domain, S the solid domain and I the fluid structure interface, we denote by the velocity field, split into fluid velocity f ∶= | F and solid velocity s ∶= | S , and by the deformation field, again with s ∶= | S and f ∶= | F . The boundary of the fluid domain Γ f ∶= F⧵I is split into inflow boundary Γ in f and wall boundary Γ wall f , where we usually assume Dirichlet conditions, Γ D f ∶= Γ in f ∪ Γ wall f , and a possible outflow boundary Γ out f , where we enforce the do-nothing outflow condition presented in Heywood et al. (1992), and 1 3 the control boundary Γ q . The solid boundary Γ s = S⧵I is split into Dirichlet part Γ D s and a Neumann part Γ N s . We formulate the coupled fluid-structure interaction problem in a strictly monolithic scheme by mapping the moving fluid domain onto the reference state via the ALE map T f (t) ∶ F → F(t) , constructed by a fluid domain deformation T f (t) = id + f (t) . In the solid domain, this map T s (t) = id + s (t) denotes the Lagrange-Euler mapping and as the deformation field will be defined globally on Ω we simply use the notation T(t) = id + (t) with the deformation gradient ∶= ∇T and its determinant J ∶= det( ).
For given desired states ̃ (t) ∈ L 2 (F) or ̃ (t) ∈ L 2 (S) , we find the global (in fluid and solid domain) velocity and deformation fields the pressure p ∈ L 2 (F) and the control parameter q ∈ Q satisfying the initial condition (0) = 0 and (0) = 0 , as solution to and subject to where the test functions are given in By 0 s we denote the solid's density, by D (t) ∈ H 1 (Ω) d and D (t) ∈ H 1 (Ω) d extensions of the Dirichlet data into the domain. The Cauchy stress tensor of the Navier-Stokes equations in ALE coordinates is given by with the kinematic viscosity f and the density f . In the solid we consider the St. Venant Kirchhoff material with the second Piola Kirchhoff tensor s based on the Green Lagrange strain tensor s and with the shear modulus s and the Lamé coefficient s . In (2) we construct the ALE extension f = | F by a simple harmonic extension. A detailed discussion and further literature on the construction of this extension is found in Yirgit et al. (2008), Richter (2017). For shorter notation, we denote by U ∶= ( , , p f ) ∈ X the solution variable and with X the corresponding ansatz space and by Φ ∶= ( , f , s , ) ∈ Y the test functions and the corresponding test space. For a control q ∈ L 2 (I) and the here given tracking-type functional constrained by linear-fluid structure interaction, we were able to proof in Failer et al. (2016) existence of a unique solution and H 1 (I) regularity of the optimal control. In addition an optimality system could be rigorously derived. Due to the missing regularity results for the here regarded nonlinear control to state mapping, no further theoretical conclusions are possible here.

Discretization
In the following we give a description of the discretization of the fluid-structure interaction system (2) in space and in time. While there exist many variants and different realizations, our choice of methods is based on the following principles • Since the FSI system is a constraint in the optimization process we base the discretization on Galerkin methods in space and time. This helps us to derive the discrete optimality system. As far as possible (up to quadrature error) we aim at permutability of discretization and optimization. • Aiming at three dimensional problems we consider methods of reasonable approximation error at feasible costs. In space we will use second order finite elements and in time a second order time stepping scheme. This approach is similar to Hron et al. (2010) or our previous work documented in Richter (2017). • Since the key component of the linear solver is a geometric multigrid method with Vanka type blocking in the smoother we choose equal-order finite elements for all unknowns, pressure, velocity and deformation adding stabilization terms for the inf-sup condition. This setup allows for efficient linear algebra and local blocking of the unknowns that is in favor of strong local couplings taking care of all nonlinearities, see also Braack and Richter (2006) for a detailed description of the realization in the context of reactive flows. • The temporal dynamics of fluid-structure interactions is governed by the parabolic/hyperbolic character of the coupling. In particular long term simulations give rise to stability problems. The Crank-Nicolson shows stability problems such that variants will be considered, see Richter and Wick (2015).

Temporal discretization
In Richter and Wick (2015) and Richter (2017, Section 4.1) many aspects of time discretization of monolithic fluid-structure interactions are discussed. It turns out that the standard Crank-Nicolson scheme is not sufficiently stable for long time simulations. Suitable variants are the fractional step theta method or shifted versions of the Crank-Nicolson scheme which we refer to as theta time stepping methods. Applied to the ode u � = f (t, u(t)) they take the form if 0 = t 0 < t 1 < ⋯ < t N = T are the discrete time steps with step size k n = t n − t n−1 . The choice = 1 2 + O(k) gives second order convergence and sufficient stability Richter and Wick (2015). Alternative approaches are the fractional step theta scheme that consists of three sub steps with specific choices for and the step size or the enrichment of the Crank-Nicolson scheme with occasional Euler steps, see Rannacher (1984).
In the context of optimization problems we aim at permutability of optimization and discretization such that Galerkin approaches are of a favor. In Richter (2014, 2015) we have demonstrated an interpretation of the general theta scheme and the fractional step theta scheme as Galerkin method with adapted function spaces: the solution is found in the space of continuous and piecewise (on I n = (t n−1 , t n ) ) linear functions, the test-space is a space rotated constant functions with jumps at the discrete time steps t n , namely The theta scheme is recovered exactly for linear problems and approximated by a suitable quadrature rule for nonlinear problems.
In case of fluid-structure interactions the domain motion term (J −1 t ⋅ ∇ , ) takes a special role since it couples temporal and spatial differential operators. In Richter and Wick (2015) various discretizations are analyzed and all found to give results in close agreement.
Here, we consider the Galerkin variant of the theta scheme and we approximate all temporal integrals by the quadrature rule [see Meidner and Richter (2014)] The resulting discrete scheme is-up to quadrature error-the standard theta time stepping scheme, which we use in our implementation for reasons of efficiency.
For the following we denote by U n ≈ U(t n ) the approximation at time t n . Further we introduce and the step t n−1 ↦ t n is given as � . (3) with J n = 1∕2(J n−1 + J n ) and ̄ n = 1∕2( n−1 + n ) . The divergence equation A div and the pressure coupling A p are fully implicit, which can be considered as a post processing step, see Meidner and Richter (2015). If the optimality system is first derived and then discretized using the Petrov-Galerkin discretization, we could observe that the control variable has to be in the theta dependent test space of the adjoint variable. As this space is very difficult to interpret the control variable q ∈ Q is approximated by piece-wise constant functions q n on every time-interval in the following. An alternative interpretation is to actually use the theta dependent test space for the adjoint variable but to approximate these integrals with the midpoint rule giving Given the choice = 1∕2 + O(k n ) this gives correct second order convergence. Numerical studies comparing both approaches did not result in a different behavior of the optimization algorithm.

Finite elements
Spatial discretization of the primal and adjoint problem is by means of quadratic finite elements in all variables on quadrilateral and hexahedral meshes. The interface I is resolved by the mesh such that no additional approximation error appears. To cope with the saddle point structure of the flow problem we use the local projection method for stabilization (Becker and Braack 2001;Frei 2016;Molnar 2015;Richter 2017). In the context of optimization problems this scheme has the advantage that stabilization and optimization commute, see Braack (2009). Further details on this and comparable approaches are found in the literature (Hron et al. 2010;Richter and Wick 2010;Richter 2017).
The use of equal order finite elements in all variables has the advantage that one set of scalar test functions { (1) h , … , (N) h } can be chosen for all variables. The discrete solution U h can then be written as h ) are small but dense local matrices of size (2d + 1) × (2d + 1) . All linear algebra routines act on these blocks, e.g. inversion of a matrix entry corresponds to the inversion of these blocks A −1 ij , which results in a better cache efficiency and reduced effort for indirect indexing of matrix and vector entries. The effect of this approach is described in Braack and Richter (2006).

Optimality system and adjoint equation
As gradient based algorithms for parameter estimation are not very common in the hemodynamics community, we shortly derive the Karush-Kuhn-Tucker system and show how gradient information can thereby be extracted. To derive the Karush-Kuhn-Tucker system, we define Lagrange multi- If the triplet U n = (p n , n , n ) ∈ X h is the solution of the discrete fluid-structure interaction system of (4) in every time step n = 0, … , N with the control parameter (q n ) N n=1 in the boundary condition, the useful identity ) the derivative of the state variable with respect to the control, we obtain via the Lagrange functional the representation of the derivative of the reduced functional. If we choose the Lagrange multiplier Z n ∈ Y h for n = N, … , 0 such that the dual problem is fulfilled, then we can evaluate the derivative of the reduced functional j((q n ) N n=1 ) in an arbitrary direction ( q) N n=1 by evaluating This enables us to apply any gradient based optimization algorithm. We will later use a limited memory version of the Broyden-Fletcher-Goldfarb-Shanno (BFGS) update formula [see for example Geiger and Kanzow (2013)] to find a local minima of the discretized optimization problem.
In every update step we first have to solve for the solution U n of the state equation (4) for n = 0, … , N and then compute the dual problem (7) for n = N, … , 0 . Thereby, the dual problem for n = N − 1, … , 1 consists of three equations with the test function Φ ∈ X h . Due to derivatives with respect to the velocity variable v n , we obtain: Due to derivatives of the Lagrangian with respect to the displacement u n , we obtain: Finally due to derivatives of the Lagrangian with respect to the pressure variable p n , we obtain: The first and last step of the discrete dual problem have a slightly different structure, but can be derived in a similar way. Since the monolithic formulation is a Petrov Galerkin formulation with different trial and test spaces, the adjoint coupling conditions differ from the primal ones. In the primal problem the solid displacement field enters as Dirichlet condition on the interface for the ALE extension problem. In the adjoint problem the shape derivatives of the adjoint ALE equation are coupled with the adjoint solid problem via a global test function which corresponds to a Neumann condition. As uf fulfills zero Dirichlet conditions on the interface, this corresponds to a back coupling of the shape derivatives into the adjoint solid problem via residuum terms. Similar to the primal problem the adjoint velocity v has to match on the interface and in addition an "adjoint dynamic" coupling condition is hidden in the test function . Therefore, block preconditioners suggested in the literature cannot be directly applied to the adjoint problem, but have to be adapted to the new structure.

Short notation for state and dual equation
Short notation of the state equation Key to the efficiency of the multigrid approach demonstrated in Failer and Richter (2020) is a condensation of the deformation unknown n from the solid problem. The last equation in (4) gives a relation for the new deformation at time t n (10) and we will use this representation to eliminate the unknown deformation and base the solid stresses purely on the last time step and the unknown velocity, i.e. by expressing the deformation gradient as Removing the solid deformation from the momentum equation will help to reduce the algebraic systems in Sect. 4. A similar technique within an Eulerian formulation and using a characteristics method is presented in Pironneau (2016Pironneau ( , 2019. For each time step t n−1 ↦ t n we introduce the following short notation for the system of algebraic equations that is based on the splitting of the solution into unknowns acting in the fluid domain ( f , f ) , on the interface ( i , i ) and those on the solid ( s , s ) . The pressure variable p acts in the fluid and on the interface.
D describes the divergence equation which acts in the fluid domain and on the interface, M the two momentum equations, acting in the fluid domain, on the interface and in the solid domain (which is indicated by a corresponding index), E describes the ALE extension in the fluid domain and U is the relation between solid velocity and solid deformation, acting on the interface degrees of freedom and in the solid. Note that M i and M s , the term describing the momentum equations, do not directly depend on the solid deformation s as we base the deformation gradient on the velocity, see (13).
Short notation dual equation We aim at applying a similar reduction scheme to the adjoint problem. Here, there is no direct counterpart to (12). Instead, we first introduce the new variable ̃ us n such that Thereby we can substitute all terms in (9), (10) and (11) which depend on us n by the new variable ̃ us n , such as (9). Now the adjoint terms M i and M s resulting from derivatives of the momentum equation with respect to the displacement variable do not depend on the adjoint velocity variable v anymore which will enable later to decouple the problem in three well conditioned subproblems. Furthermore the "adjoint dynamic" coupling conditions now corresponds to equivalents of adjoint boundary forces on the interface as in the state equation.
For each time step t n+1 ↦ t n we introduce again a short notation for the system of algebraic equations that is based on the splitting of the adjoint solution into unknowns acting in the fluid domain

Solution of the algebraic systems
In Failer and Richter (2020), we have derived an efficient approximated Newton scheme for the forward fluid-structure interaction problem. We briefly outline the main steps and then focus on transferring these ideas to the dual equations. The general idea is described by the following two steps 1. In the Jacobian, we omit the derivatives of the Navier-Stokes equations with respect to the fluid domain deformation, which results in an approximated Newton scheme. In Richter (2017, chapter 5) it is documented that this approximation will slightly increase the iteration counts of the Newton scheme. On the other hand, the overall computational time is nevertheless reduced, since assembly times for these neglected terms are especially high. Since the Newton residual is not changed, the resulting nonlinear solver is of an approximated Newton type.
We use the discretization of the relation t = between solid deformation and solid velocity, namely n+1 = n + k n+1 + (1 − )k n to reformulate the solid's deformation gradient based on the velocity instead of the deformation. This step has been explained in the previous section.
These two steps, the first one being an approximation, while the second is an equivalence transformation, allow to reduce the number of couplings in the Jacobian in such a way that each linear step falls apart into three successive linear systems. The first one describes the coupled momentum equation for fluid-and solid-velocity, the second realizes the solid's velocity-deformation relation and the third one stands for the ALE extension. We finally note that the approximations only involve the Jacobian. The residual of the systems is not altered such that we still solve the original problem and compute the exact discrete gradient. Then, in Sect. 4.2 we describe the extension of this solution mechanism to the adjoint system. Two major differences occur: first, the adjoint system is linear, such that we realize the solver in the framework of a preconditioned Richardson iteration. The preconditioner takes the place of the approximated Jacobian. Second, the adjoint interface coupling conditions differ from the primal conditions as outlined in the last paragraph of Sect. 3.3. This will call for a modification of the condensation procedure introduced as second reduction step in the primal solver.

Solution of the primal problem
In each time step of the forward problem we must solve a nonlinear problem. We employ an approximated Newton scheme where (l) is a line search parameter, U (0) an initial guess. By A � (U) we denote the Jacobian, by Ã � (U) an approximation. As outlined in Failer and Richter (2020) the Jacobian is modified in two essential steps: first, in the Navier-Stokes problem, we skip the derivatives with respect to the ALE discretization. These terms are computationally expensive and they further introduce the only couplings from the fluid problem to the deformation unknowns. In Richter (2017, chapter 5) it has been shown that while this approximation does slightly worsen Newton's convergence rate, the overall efficiency is nevertheless increased, as the number of additional Newton steps is very small in comparison to the savings in assembly time. Second, we employ the reduction step outlines in Sect. 3.4, which is a static condensation of the deformation unknowns from the solid's momentum equation. Taken together, both steps completely remove all deformation couplings from the com-

Dual
Due to the unsymmetrical structure of the fluid-structure interaction model the block collocation and coupling of the blocks in the transposed Jacobian A � (U) T in the dual problem differs to the Jacobian of the primal problem. This stays in strong relation to the adjoint coupling conditions, see Sect. 3.3. Hence, block preconditioners developed for the state problem can not be applied in a black box way to the linear systems arising in the dual problem, but have to be adjusted. Furthermore, the dual system is linear such that the approximated Newton scheme must be replaced by a different concept. We start by indicating the full system matrix of the dual problem

3
A Newton multigrid framework for optimal control of fluid-… given as the transposed of the primal Jacobian, A D = A � (U) T , see Failer and Richter (2020).
For solving the dual problem we want to mimic the primal approach: first, approximate the system matrix by neglecting couplings, second, use the static condensation as described in Sect. 3.4 in (15). As the problem is linear, a direct modification of the system matrix would alter the dual solution. Instead, we approximate the solution by a preconditioned Richardson iteration with an inexact matrix Ã � D ≈ A � (U) T as preconditioner (approximated by a geometric multigrid solver) where

̃ us s } and the update in every Richardson iteration is given by W
The residual is computed based on the full Jacobian A � (U) T (including the ALE derivatives) such that we still converge to the original adjoint problem. Since we never assemble the complete Jacobian A � (U) in the primal solver, the adjoint residual B d − A � (U) T Z (l−1) can be established in a matrix free setting.
Then, similar to the described approach in case of the primal system, we neglect the ALE terms (shaded entries). Finally, we reorder to reach the preconditioned iteration with the preconditioner Ã D that decomposes into three separate steps. First, the equation for the adjoint mesh deformation variable Usually a symmetric extension operator E f f can be chosen. This avoids re-assembly of this matrix and possible preparations for the linear solver. See Richter (2017, Section 5.3.5) for different efficient options for extension operators. Second, the update for the adjoint solid deformation, which only involves inversion of the mass matrix and finally the update for the adjoint velocity and adjoint pressure The numbering of the right hand side d 1 , … , d 7 is according to (23). As we do not modify the residuum, the derivatives with respect to the ALE transformation i ( v f , v i ) on the interface. These terms contain the adjoint geometric coupling condition. An explicit update by one vector-addition as for the corresponding primal equation is not possible. The "adjoint kinematic" and "adjoint dynamic" coupling conditions are fully incorporated in (28), similar as for the state equation, and thereby these coupling conditions are fully resolved in every Richardson iteration.

Solution of the linear problems
In each step of the Newton iteration for solving the state equation and in each step of the Richardson iteration in the case of the adjoint system, we must approximate three individual linear systems of equations. The mesh-update problems are usually of elliptic type, the vector Laplacian or a linear elasticity problem. Here, standard geometric multigrid solvers are highly efficient. Problem (27) and the primal counterpart correspond to zero order equations. Multigrid solvers or the CG method converge with optimal efficiency. It remains to approximate the coupled momentum equations, given by (28) in the dual case. Here, we are lacking any desirable structure. The matrices are not symmetric, they feature a saddle-point structure and involve different scaling of the fluid-and solid-problem. We approximate these equations by a GMRES iteration that is preconditioned with a geometric multigrid (27) solver. Within the multigrid iteration we employ a smoother of Vanka type, where we invert local patches exactly. These patches correspond to all degrees of freedom of one element (within the fluid) and to a union of 2 d elements (within the solid). For the sake of simplicity (in terms of implementational effort) we use this highly robust solver also for the other two problems, despite their simpler character. For details we refer to Failer and Richter (2020).

Algorithm
To get an overview how the final optimization routine works we summarized all the intermediate steps in the following algorithm:

Problem configuration 2d
We modify the well known FSI Benchmark from Hron and Turek (2006) by adding an additional boundary Γ q as in Fig. 1. The material parameters are chosen as for the FSI 1 Benchmark. In the solid the Lame parameters with = 2.0 × 10 6 kg/ms 2 and = 0.5 × 10 6 kg/ms 2 and a fluid viscosity f = 0.001 m 2 ∕s are chosen. The solid and fluid densities are given by s = 1000.0 kg/m 3 and f = 1000.0 kg/m 3 . The inflow velocity is increased slowly during the time interval I = [0, 2] as for the instationary FSI 2 and FSI 3 benchmark.
In the first example, we would like to determine the pressure profile q(t) ∈ L 2 (I) at the control boundary Γ q on the time interval I = [0, 6] , leading to the displacement profile over time in y-direction in the area Ω obs = {0.525 ≤ x ≤ 0.6, 0.19 ≤ y ≤ 0.21} at the tip of the flag (see Fig. 1). To do so, we minimize the functional constrained by the fluid-structure interaction problem. We discretize the system as presented in Sect. 3 in time using a shifted Crank-Nicolson time stepping scheme with time step size k = 0.01 and = 0.5 + 2k and quadrilateral meshes with varying mesh size, see Table 1. The control variable is chosen to be piece-wise constant on every time interval ( dim(Q) = 600 ). The Tikhonov regularization parameter is set to = 1.0 × 10 −17 .

Problem configuration 3d
In the second example we regard a pressure wave in straight cylinder as presented in Gerbeau and Vidrascu (2003). The cylinder has the length 5 cm and a radius of 0.5 cm . The cylinder is surrounded by an elastic structure with constant thickness h = 0.1 cm . The elastic structure is clamped at the inflow boundary and the outflow domain is fixed in x-direction and free to move in y-and z-direction. At the inlet we describe a pressure step function p in = 1.33 × 10 4 g cm −1 s = 10 mmHg for t ≤ 0.003 s , afterwards we set the pressure to zero. In the solid the Lame parameters with = 1.73 × 10 6 g/cms 2 and = 1.15 × 10 6 g/cms 2 and a fluid viscosity

3
A Newton multigrid framework for optimal control of fluid-… f = 0.03 cm 2 ∕s are chosen. The solid and fluid densities are given by s = 1.2 g∕cm 3 and f = 1.0 g/cm 3 . We plotted the solution at t = 0.006 s in Fig. 2.
If only a do-nothing condition is described with constant pressure at the outflow boundary, then pressure waves are going to be reflected at the outflow boundary. If the pressure along the outflow boundary is chosen appropriate the pressure wave will leave the cylindrical domain without any reflection such that the system will be at rest after some time. To determine the corresponding pressure profile q(t) ∈ L 2 (I) on the time interval I = [0, 0.04] , we minimize the kinetic energy in the fluid domain for t > 0.03 s . Hence we minimize the functional constrained by the fluid-structure interaction problem. In time we use, as already in the previous example, a shifted Crank-Nicolson time stepping scheme with time step size k = 0.0001 and = 0.5 + 2k and hexahedral meshes with varying mesh size, see Table 1. Only at the time points t = 0 and t = 0.003 , we use for four steps a time step-size of k = 0.00005 with = 1.0 . Thereby, no artificial effects occur due to the jump in the pressure at the inflow boundary. The Tikhonov regularization parameter is set to = 1.0 × 10 −8 .

Optimization algorithm
Given the computed gradient information using Formula (8), we apply a BFGS algorithm implemented in the optimization library RoDoBo (Becker et al. 2020b) to solve the optimization problem. We use a limited memory version as presented e.g. in Geiger and Kanzow (2013) such that there is no need to assemble and store the BFGS update matrix. To guarantee that the update matrix keeps symmetric and positive definite a Powell-Wolfe step size control should be used. As this step size criteria is in most cases very conservative, we only check if we have descend in the functional value. Control constraints could be added in the optimization algorithm via projection of the update in every optimization step. In this paper we only regard unconstrained examples. Fast convergence of the BFGS algorithm only can be expected close to the optimal solution. Hence, we take advantage of the mesh hierarchy, which is used in the geometric multigrid algorithm, and first solve the optimization problem on a coarse grid to have a good initial guess on finer meshes. As the computation time rises for 3d examples on finer meshes very fast, we can save a lot of computation time using this approach. In the following we reduce the norm of the gradient by a factor of 10 −1 in every optimization loop and then refine the mesh and restart the optimization with the control from the coarser mesh. To compute the gradient, we have to solve the state and dual problem for all time steps. The state solutions are stored on the hard disc and loaded during the computation of the dual problem if necessary. Thereby, only the current state and adjoint solutions have to be held in the memory. In every Newton step a GMRES iterative solver preconditioned with a geometric multigrid solver provided by the FEM software library Gascoigne (Becker et al. 2020a) is used. All linear systems are solved up to a relative accuracy of 10 −4 . In every time step the initial Newton residuum is reduced by a factor of 10 −6 . We use the same tolerances for the state and dual problem. The matrices are only reassembled, if the nonlinear convergence rate falls below = 0.05 as in Failer and Richter (2020). In every dual step, the matrices are assembled at least once at the beginning of every Richardson iteration.

Numerical results 2d example
In Fig. 3 we plot the value of the regularized functional j(q) and the norm of the gradient ‖j � (q)‖ in every optimization step. The optimization algorithm is started with q n = 0 for n = 1, … , N . The computed optimal control is given in Fig. 4. The functional value reduces in every optimization step and the gradient can be reduced by a factor of 10 −3 after less then 40 optimization cycles (see Fig. 3). In addition we Table 1 Spatial degrees of freedom for 2d and 3d configuration on every refinement level. In 3d the mesh on level 4 and 5 are locally refined along the interface. In time we use a uniform partitioning. In 3d, we add 4 initial backward Euler steps with reduced step size for smoothing at times t = 0 s and t = 0.03 s each plotted the optimal solution for the 2d example on the finest mesh at the point B in the center of the observation domain Ω obs and compare the solution to the reference solution in Fig. 4. Only around the kink of the desired state a mismatch between desired state and optimized solution can be seen. To evaluate our approach to solve the optimization algorithm first on coarser grids and then to refine systematically, we restarted the optimization algorithm directly on meshlevel 5. A similar behavior in functional and gradient to the previous approach can be observed in Fig. 3. To reduce the gradient to a tolerance of 10 −12 only 20 optimization loops are required. But since about 14000 s of computational time are needed to solve one cycle of the state and the adjoint system for all time steps on meshlevel 5, but only less then 4000 s on the coarser meshlevel 4 (even less time on still coarser grids), systematical refinement of the mesh is much more efficient. The computations were performed on a Intel(R) Xeon(R) Gold 6150 CPU @ 2.70GHz with 18 threads. We parallelized the assembling of matrix and the Vanka smoother as well as the matrix vector multiplication via OpenMP, see Failer and Richter (2020).

Numerical results 3d example
For the 3d example we used the same optimization algorithm as in the 2d case.
We can see in Fig. 5 that using the optimized pressure profile at the outflow boundary about 98.9% of the kinetic energy after t > 0.03s now leaves the cylinder. The jumps in the gradient after every refinement step indicate that a more accurate computation on the coarse grid would not result in better starting values on the finer meshes. The norm of the gradient could be reduced from 1.15 × 10 −4 to 3.68 × 10 −7 in 23 optimization steps, whereby only 9 optimization cycles were necessary on the computationally costly fine grid on meshlevel 4. To compare the controlled solution with the uncontrolled solution, we computed in addition the solution on a cylinder with length 10 cm . As the reflection on the outflow boundary occurs later the pressure and flow values in the center at x = 5 cm can be seen as reference values for optimal non reflective boundary conditions for t < 0.02 s . As we only control the pressure on the fluid domain, reflections on the solid boundary can still occur. Furthermore we can observe that the pressure is not constant along the virtual outflow boundary at x = 5 cm for the long cylinder. Thus, we can not expect the reference solution to fully match the optimized solution. Nevertheless, we can see in outflow boundary at x = 10 cm . This explains the different behavior of pressure and outflow after t = 0.025 s.

Test of the Newton scheme, of the Richardson iteration and of the iterative linear solver
How to evaluate the performance of the quasi Newton scheme or of the iterative linear solver is not so obvious. Due to the changing control in every optimization cycle and the nonlinear character of the problem, the condition numbers of the matrices occurring in the linear subproblems will vary in every time step and optimization cycle. Hence, we first compute only one optimization step with q = 0 on various meshlevels to analyze the h-dependence of our solution algorithm. Thereby, we compare the mean number of Newton/Richardson iterations and GMRES steps per time step on every meshlevel. Afterwards, we compute mean values in every optimization loop to analyze how the performance can vary during the optimization procedure. In Tables 2 and 3 we present the mean number of Newton/Richardson iterations and Matrix assemblies per time step for the 2d and 3d examples. In addition we present the mean number of GMRES steps needed per Newton/Richardson iteration to solve the linear subproblems (22), (21) and (20) and (26),  (27) and (28). We can observe that the number of Newton/Richardson iterations per time step ranges between 3 and 4 for state and dual problem. The coupled momentum equation remains the most complex problem with the highest number of steps required. Equations (21) and (27) (2020) a more detailed analysis of the smoother in the geometric multigrid algorithm can be found. As the dual equation is linear with respect to the adjoint variable, we would have expected to need only one Richardson iteration per time step. However, since we neglected terms occurring due to the ALE transformation, we loose the optimal convergence and need about 3 Richardson iterations per time step. On the other hand, the matrices occurring in cascade of subproblems in the dual system have the same condition number as the matrices in the state equation. Only by this approximation and splitting, iterative solvers can successfully be applied to solve the linear problems. The number of GMRES steps in the dual subproblems in Tables 2 and 3 are rather small and close to the values for the state problem. As (28) and (20) are still fully coupled problems of fluid and elastic solid, most of the computational time is spent in solving these two subproblems. The least number of GMRES steps is needed to invert the mass matrix in (21) and (27). In all subproblems the number of GMRES steps only increases slightly under mesh refinement.
In Fig. 7 we show the mean number of GMRES steps per Newton step in every timestep and the number of Newton steps per timestep. For the given examples the values only vary slightly over time. Therefore the mean values in Tables 2 and 3 represent the overall behavior very well.
To understand how the performance of the solution algorithm in the case of the 2d example varies during the optimization loop, we show in Fig. 8 the average number of Newton steps per time step and the average number of GMRES steps per Newton step for each optimization step. The computation was started with q = 0 on meshlevel 5. No further mesh was applied.
The dependency on the time step size of the presented quasi Newton scheme for the state equation was already analyzed in Failer and Richter (2020). Therein, we could observe an increasing Newton iteration count for larger time steps, in particular for very large time steps. A similar behavior can be expected for the dual problem. While the presented solution approach can be regarded as highly efficient and appropriate for optimal control of nonstationary fluid-structure interactions, stationary or quasi stationary applications will call for alternative approaches like the geometric multigrid solver presented by Aulisa et al. (2018).
In Failer and Richter (2020) more details regarding the computational time and savings due to parallelization can be found. As we have to evaluate state and adjoint variables, as well as additional terms due to linearization in every dual 1 3  (20) and (28) plotted over optimization steps . Right: mean Newton steps per time step plotted over optimization steps step, the computational time to assemble the matrix and to compute the Newton residuum is slightly larger for the dual equation then for the state equation.

Summary
We extended the Newton multigrid framework for monolithic fluid-structure interactions in ALE coordinates presented in Failer and Richter (2020) to the dual system of fluid-structure interaction problems. The solver is based on replacing the adjoint solid deformation by a new variable and on skipping the ALE derivatives within the adjoint Navier-Stokes equations. As we do not modify the residuum, state and dual solution in each time step still converge to the exact discrete solution. The adjoint coupling conditions incorporated in the monolithic formulation are still fulfilled. As we compute correct gradient information, we see fast convergence in our optimization algorithm. The coupled system is better conditioned (as compared to monolithic Jacobians) which allows to use very simple multigrid smoothers that are easy to parallelize. Only this makes gradient based algorithm feasible and efficient for 3d fluid-structure interaction applications, where memory consumption prevents the use of direct solvers. It would be straightforward to combine the presented algorithm with dualweighted residual error estimators for mesh and time step refinement. Instead of global refinement of the mesh after every optimization loop the error estimators indicate where to refine the mesh locally. The sensitivity information from the optimization algorithm can be directly used to evaluate the error estimators.