A Multi-Vector Interface Quasi-Newton Method with Linear Complexity for Partitioned Fluid-Structure Interaction
AA Multi-Vector Interface Quasi-Newton Method with Linear Complexity forPartitioned Fluid-Structure Interaction
Thomas Spenke ∗ , Norbert Hosters, Marek Behr Chair for Computational Analysis of Technical Systems (CATS), Center for Computational Engineering Science (CCES),RWTH Aachen University, 52062 Aachen, Germany
Abstract
In recent years, interface quasi-Newton methods have gained growing attention in the fluid-structure interactioncommunity by significantly improving partitioned solution schemes: They not only help to control the inherent added-mass instability, but also prove to substantially speed up the coupling’s convergence. In this work, we present a novelvariant: The key idea is to build on the multi-vector Jacobian update scheme first presented by Bogaers et al. [1] andavoid any explicit representation of the (inverse) Jacobian approximation, since it slows down the solution for largesystems. Instead, all terms involving a quadratic complexity have been systematically eliminated. The result is a newmulti-vector interface quasi-Newton variant whose computational cost scales linearly with the problem size.
Keywords:
Fluid-Structure Interaction, Interface Quasi-Newton Methods, Partitioned Algorithm
1. Introduction
Partitioned solution schemes for fluid-structure interaction (FSI) are widely used in modern computational engi-neering science, as they o ff er great flexibility and modularity concerning the two solvers employed for the fluid and thestructure. Treated as black boxes, the solvers are coupled only via the exchange of interface data. The major downsideof partitioned schemes, however, is an inherent instability caused by the so-called added-mass e ff ect [2, 3]. Dependingon the application, its influence might be severe enough to impede a numerical solution – if no counter-measures aretaken.The simplest way to retain stability is a relaxation of the exchanged coupling data. The relaxation factor mighteither be constant or updated dynamically via Aitken’s method [4, 5]. Unfortunately, the high price to be paid in termsof the coupling’s convergence speed often renders this approach infeasible.As opposed to this, interface quasi-Newton (IQN) methods have proven capable of both stabilizing and accel-erating partitioned solution schemes [6, 7, 8]. Identifying the converged solution as a fixed point, their basic ideais to speed up the coupling iteration using Newton’s method. Since the required (inverse) Jacobian is typically notaccessible, interface quasi-Newton approaches settle for approximating it instead.Early work in the field has been done by Gerbeau et al. [9] as well as van Brummelen et al. [10]. A breakthrough,however, was the interface quasi-Newton inverse least-squares (IQN-ILS) method by Degroote et al. [11]: On the onehand, it directly solves for the inverse Jacobian required in the Newton linearization; on the other hand, the IQN-ILSvariant introduces the least-squares approximation based on input-output data pairs that is still common today.Research of the last decade has shown that a reutilization of data pairs from previous time steps is extremelyadvantageous. Unfortunately, an explicit incorporation in the IQN-ILS method su ff ers from numerical di ffi cultiessuch as rank deficiency; moreover, good choices for the number of incorporated time steps are in general problem-dependent [12, 13]. However, the works by Degroote and Vierendeels [14] as well as Haltermann et al. [15] haveshown that filtering out linearly dependent data pairs is an e ff ective way of alleviating these issues. ∗ Corresponding author
Email addresses: [email protected] (Thomas Spenke), [email protected] (Norbert Hosters), [email protected] (Marek Behr)
Preprint submitted to Computer Methods in Applied Mechanics and Engineering January 23, 2020 a r X i v : . [ c s . C E ] J a n s an alternative, Bogaers et al. [1] and Lindner et al. [13] formulated a very beneficial implicit reutilization ofpast information, yielding the interface quasi-Newton inverse multi-vector Jacobian (IQN-IMVJ) method. Its onlyreal drawback is the required explicit representation of the approximated inverse Jacobian: Since the related cost isincreasing quadratically with the number of degrees of freedom at the FSI boundary, this explicit form seriously slowsdown the numerical simulation for larger problem scales.In this work, we present an enhancement of the IQN-IMVJ concept tackling exactly this shortcoming: Retainingthe e ff ective implicit reutilization of past data, the new interface quasi-Newton implicit multi-vector least-squares(IQN-IMVLS) method completely avoids any explicit Jacobian or quadratic dependency on the interface resolution.While the advantages of the multi-vector approach are preserved, the resulting linear complexity ensures the negligi-bility of computational cost even for large problem scales.This paper is structured as follows: Section 2 presents the governing equations of the considered FSI problems,before Section 3 briefly covers the numerical methods applied to solve them. In Section 4, the IQN-ILS and theIQN-IMVJ approach are discussed in detail, before the IQN-IMVLS method is derived. The e ffi ciency of the newapproach is validated in Section 5 based on numerical test cases.
2. Governing Equations
In general, fluid-structure interaction considers a fluid domain Ω ft ⊂ R nsd and a structure Ω st ⊂ R nsd , that areconnected via a shared boundary Γ f st = ∂ Ω ft ∩ ∂ Ω st ⊂ R nsd − , the FSI interface . The subscript t refers to the time level,while nsd denotes the number of spatial dimensions. This section introduces the models and equations employed forthe two subproblems along with the boundary conditions interlinking their solutions at the shared interface. The velocity u f ( x , t ) and the pressure p f ( x , t ) of the fluid are governed by the unsteady Navier-Stokes equationsfor an incompressible fluid, reading ρ f (cid:32) ∂ u f ∂ t + u f · ∇ u f − f f (cid:33) − ∇ · T f = in Ω ft ∀ t ∈ [0 , T ] , (1a) ∇ · u f = Ω ft ∀ t ∈ [0 , T ] . (1b)Therein, ρ f is the constant fluid density, while f f denotes the resultant of all external body forces per unit mass offluid. For a Newtonian fluid with the dynamic viscosity µ f , the Cauchy stress tensor T f is modeled by Stokes law as T f ( u f , p f ) = − p f I + µ f (cid:16) ∇ u f + ( ∇ u f ) T (cid:17) . (2)The problem is closed by defining not only a divergence-free initial velocity field u f , but also a prescribed velocity g f on the Dirichlet boundary Γ fD , t and prescribed tractions h f on the Neumann boundary Γ fN , t with its outer normal n f : u f ( x , t = = u f ( x ) in Ω f , (3a) u f = g f on Γ fD , t ∀ t ∈ [0 , T ] , (3b) T f n f = h f on Γ fN , t ∀ t ∈ [0 , T ] . (3c) The response of the structure to external loads is expressed via the displacement field d s ( x , t ), which is governedby the equation of motion, stating a dynamic balance of inner and outer stresses: ρ s d d s dt = ∇ · T s + b s in Ω st ∀ t ∈ [0 , T ] . (4)2n this relation, ρ s denotes the material density and b s the resultant of all body forces per unit volume, whereas T s represents the Cauchy stress tensor.As constitutive relation, the St. Venant-Kirchho ff material model is used: It relates the 2nd Piola-Kirchho ff stresses S : = det( F ) F − T s F − T to the Green-Lagrange strains E : = (cid:16) F T F − I (cid:17) via a linear stress-strain law, reading S = λ s tr ( E ) + µ s E . (5)Therein, F denotes the deformation gradient, while λ s and µ s are the Lam´e constants. As the Green-Lagrange straindefinition forms a nonlinear kinematic relation, the structural model is geometrically nonlinear. Hence, it is capableof representing large displacements and rotations, but only small strains [16, 17].Collecting all this information, the equation of motion can be expressed in the (undeformed) reference configura-tion Ω s in a total Lagrangian fashion: ρ s d d s dt = ∇ · (cid:16) SF T (cid:17) + b s in Ω s ∀ t ∈ [0 , T ] . (6)Again, the problem is closed by defining an initial displacement d , which is typically zero, and a set of boundaryconditions on two complementary subsets of Γ s = ∂ Ω s : Prescribing the displacement g s on the Dirichlet part Γ sD , andthe tractions h s on the Neumann part Γ sN , , with the outer normal in the reference state n , the conditions read d s ( x , t = = d s ( x ) in Ω s , (7a) d s = g s on Γ sD , ∀ t ∈ [0 , T ] , (7b) F S n = h s on Γ sN , ∀ t ∈ [0 , T ] . (7c) The essence of fluid-structure interaction is that the two subproblems cannot be solved independently, since theirsolution fields are interlinked via the so-called coupling conditions arising at the shared interface Γ f st : • The kinematic coupling condition requires the continuity of displacements d f and d s , velocities u f and u s , andaccelerations a f and a s across the FSI boundary [18]: d f = d s on Γ f st ∀ t ∈ [0 , T ] , (8a) u f = u s on Γ f st ∀ t ∈ [0 , T ] , (8b) a f = a s on Γ f st ∀ t ∈ [0 , T ] . (8c) • In agreement with Newton’s third law, the dynamic coupling condition, T f n f = T s n s on Γ f st ∀ t ∈ [0 , T ] , (9)enforces the equality of stresses at the interface. Therein, n f and n s denote the associated normal vectors [19].In the continuous problem formulation, satisfying both these coupling conditions for every time t ∈ [0 , T ] ensures theconservation of mass, momentum, and energy over the FSI boundary [20].
3. Numerical Methods
The Navier-Stokes equations formulated in Section 2.1 are discretized in this work using P1P1 finite elements,i.e., linear interpolations for both velocity and pressure. Numerical instabilities, caused by P1P1 elements violatingthe LBB condition, are overcome using a
Galerkin / Least-Squares (GLS) stabilization [21, 22].Contrary to the usual practice, both space and time are discretized by finite elements. More precisely, we employthe deforming-spatial-domain / stabilized-space-time (DSD / SST) approach introduced by Tezduyar and Behr [23, 24].By formulating the variational form over the space-time domain, this approach inherently accounts for an evolvingspatial domain. While the finite-element interpolation functions are continuous in space, the linear discontinuousGalerkin ansatz in time allows to solve one space-time slab after another.The mesh is adjusted to moving boundaries, e.g., the FSI interface, by means of interface tracking [18, 25]. Theresulting mesh distortion is compensated by the elastic mesh update method (EMUM) [26].3 .2. Structural Solution
The structural subproblem is discretized in space by isogeometric analysis [27], while the time integration isperformed using the generalized- α scheme [28, 29]. Introduced by Hughes et al. [30] in 2005, isogeometric analysis is a finite-element variant aimed at achievinggeometrical accuracy and a closer linkage between CAD and numerical analysis . The essential idea is to express thesolution space via the same basis functions as used for the geometry description. Since CAD systems are commonlybased on non-uniform rational B-Splines (NURBS), we choose a NURBS ansatz for the unknown displacement field.Mathematically, a NURBS is a linear combination of its n control points x I ∈ R nsd and the rational basis functions N I ( ξ ), defined in the parameter space ξ ∈ Ξ ⊂ R d , with the spline dimension d [31]. A NURBS surface ¯x ( ξ ) ∈ R nsd hence has the form ¯x ( ξ , ξ ) = n (cid:88) I = N I ( ξ , ξ ) x I with ( ξ , ξ ) = ξ (10) Shell structures are thin-walled structures capable of providing lightweight, cost-e ffi cient and yet stable construc-tions for numerous engineering applications [32]. Mathematical shell models exploit the small thickness by reducingthe structure from a volumetric description to the midsurface plus an interpolation over the thickness h . Combiningisogeometric analysis with Reissner-Mindlin shell theory, a NURBS midsurface ¯x ( ξ , ξ ) and a linear interpolationare chosen; any position in the structural domain is hence expressed as x ( ξ , ξ , ξ ) = ¯x ( ξ , ξ ) + ξ b ( ξ , ξ ) . (11)The parametric coordinate ξ ∈ [ − h , + h ] changes along the thickness direction, defined by the director vector b ( ξ , ξ ).Consequently, the unknown displacement field is a combination of the midsurface displacement ¯d s ( ξ , ξ ) and a secondterm accounting for changes of the director vector ∆ b ( ξ , ξ ): d s ( ξ , ξ , ξ ) = ¯d s ( ξ , ξ ) + ξ ∆ b ( ξ , ξ ) . (12)For detailed information on nonlinear isogeometric Reissner-Mindlin shell elements, the works by Dornisch andKlinkel [33, 34, 35] are recommended. This work pursues a partitioned FSI coupling approach: Two distinct solvers are employed for the fluid and thestructure; they are connected via a coupling module, which handles the exchange of interface data in accordanceto the coupling conditions. While the strengths of this approach are its great flexibility and modularity regardingthe single-field solvers, these advantages come at the price of additional considerations to be made in terms of dataexchange: • In general, the meshes – or even the discretization techniques – do not match at the FSI boundary. Hence, aconservative projection is required to transfer relevant data between the two solvers. In this work, we employa spline-based variant of the finite interpolation elements (FIE) method; a detailed description can be found inHosters et al. [36]. For the sake of observability, however, any e ff ects of this spatial coupling are neglected inthe following, since they do not interfere with the presented methods. • Due to the (potentially strong) interdependency between the two subproblems, a consistent FSI solution ingeneral requires an iterative procedure, which is referred to as strong (temporal) coupling [18, 14].4 luid SolverT fk = F ( x kd ) Structural Solver˜x kd = S ( T fk ) k = k + Yes
Next TimeStep
Nok = PreviousTime Step ˜x kd ? ≈ x kd CFD Loads T fk Interface Deformation x k + d = ˜x kd Figure 1: Sketch of the Dirichlet-Neumann coupling scheme.
Nowadays, the de facto standard among strong coupling approaches for FSI problems is the Dirichlet-Neumanncoupling scheme illustrated in Figure 1: The first coupling iteration ( k =
0) of a new time step t n → t n + starts withthe fluid solver: Based on the interface deformation x kd of time level t n (or an extrapolated one), it computes a newflow field. In compliance with the dynamic coupling condition, the resulting fluid stresses T fk = F ( x kd ) acting on Γ f s are passed as a Neumann condition to the structural solver, which determines the corresponding interface deformation ˜x kd = S ( T fk ). If it matches (within a certain accuracy) the previous deformation , i.e., ˜x kd ≈ x kd , the coupled system hasconverged and we can go on to the next time step. Otherwise, the interface deformation x k + d = ˜x kd is passed back tothe flow solver as a Dirichlet condition, following kinematic continuity, and the next coupling iteration k + deformation at the interface x d [4, 13]; the fixed-point operator H ( x kd ) ≡ S ◦ F ( x kd ) corresponds to thetwo subsequent solver calls, i.e., to running one coupling iteration k → k + ˜x kd = H ( x kd ) , x k + d = ˜x kd . (13)Finding the converged solution of the next time level is hence equivalent to finding a fixed point x (cid:63) d = H ( x (cid:63) d ) of theinterface deformation. Defining the fixed-point residual R ( x kd ) and its inverse ˜R ( ˜x kd ) as R ( x kd ) : = H ( x kd ) − x kd = ˜x kd − x kd , (14a) ˜R ( ˜x kd ) : = ˜x kd − H − ( ˜x kd ) = ˜x kd − x kd , (14b)with ˜R ( ˜x kd ) ≡ R ( x kd ), the convergence criteria can analogously be expressed as a root of the residual, i.e., R ( x ∗ d ) = ˜R ( x ∗ d ) = x ∗ d − x ∗ d = .Note that although it is much less common, the fixed-point ansatz and the convergence criteria could analogouslybe formulated based on the fluid loads at the FSI interface. Unfortunately, partitioned algorithms for fluid-structure interaction involving incompressible fluids exhibit aninherent instability, which may be severe, depending on the simulated problem. This so-called (artificial) added-masse ff ect originates from the CFD forces inevitably being determined based on a defective structural deformation duringthe iterative procedure. As a consequence, they di ff er from the correct loads of the next time level. This deviationpotentially acts like an additional fluid mass on the structural degrees of freedom. More precisely, due to the kinematiccontinuity an overestimated structural deformation entails an exaggerated fluid acceleration – and hence excessiveinertia stress terms. In many cases the e ff ect is amplified throughout the coupling iterations, causing a divergence ofthe coupled problem. A major influencing factor are high density ratios between the fluid and the structure, ρ f /ρ s ,but dependencies on the viscous terms and the temporal discretization are observed, too. We recommend the worksby F¨orster [2], F¨orster et al. [3], and Causin et al. [37] for more details on the added-mass e ff ect. Moreover, vanBrummelen et al. [38] investigate similar di ffi culties for compressible flows.5 .3.3. Stabilizing and Accelerating the Coupling Scheme One way to deal with the added-mass instabilities is to adjust the computed interface deformation ˜x kd by an updatestep x k + d = U ( ˜x kd ) before passing it back to the flow solver. Depending on the chosen update technique, both thestability and e ffi ciency of the coupling scheme can be improved.The simplest way to increase stability is a constant under-relaxation of the interface deformation: x k + d = U Relax ( ˜x kd ) = ω ˜x kd + (1 − ω ) x kd , (15)with ω <
1. E ff ectively, it yields an interpolation between the current and the previous interface displacements ˜x kd and x kd . As the factor must be chosen small enough to avoid the coupling’s divergence for any time step, unfortunately thisapproach often comes at a very high price in terms of e ffi ciency [7]. Aitken’s dynamic relaxation [4, 5] tackles thisissue by dynamically adapting the relaxation factor in Equation (15) for each coupling iteration by ω k = − ω k − ( R ( x k − d )) T ( R ( x kd ) − R ( x k − d )) (cid:107) R ( x kd ) − R ( x k − d ) (cid:107) , (16)i.e., based on the two most recent fixed-point residuals R ( x kd ) and R ( x k − d ). Despite its rather simple implementation,in many cases Aitken’s relaxation provides significant speed-ups without perturbing the stability of the relaxation.Still, the performance of the coupling scheme can be pushed further by more sophisticated approaches like interfacequasi-Newton (IQN) methods. Remark 1:
In the following, we will express the interface deformation as the structural solution at the FSI interface.While its representation via the boundary displacement of the fluid mesh is likewise possible, the typically finerinterface resolution would result in higher numerical cost for the discussed update techniques.
4. Interface Quasi-Newton Methods
In a partitioned solution procedure for fluid-structure interaction, the single-field solvers are typically the mostexpensive part, rendering all further costs, e.g., for data exchange, negligible in comparison [1, 39].Since the converged time step solution has been identified as a root of the fixed-point residual R ( x d ) or its inverseform ˜R ( ˜x d ) defined in Equation (14), increasing the e ffi ciency essentially comes down to reducing the number ofcoupling iterations required to find such a root. Consequently, employing Newton’s method as an update step of theinterface deformation would be a promising approach to speed up convergence, yielding [11, 13] x k + d = U Newton ( ˜x kd ) = ˜x kd − J − R ( ˜x kd ) ˜R ( ˜x kd ) . (17)Unfortunately, an evaluation of the required inverse Jacobian J − R : = d ˜R / d ˜x d in general is not accessible, as it wouldinvolve the derivatives of both the fluid and the structural solver. The idea of interface quasi-Newton methods isto circumvent this issue by approximating the (inverse) Jacobian rather than evaluating it exactly. Using a Taylorexpansion of ˜R ( ˜x d ), we can approximate J − R via J − R ( ˜x kd ) ∆ R ( ∆ ˜x d ) ≈ ∆ ˜x d , (18)based on an input-output data pair. Such a pair is formed by some change in the interface deformation ∆ ˜x d ∈ R m andthe corresponding change of the fixed-point residual ∆ R ∈ R m . Therein, m denotes the number of structural degreesof freedom at the FSI interface . Essentially, this idea is an m -dimensional version of approximating a derivative withthe slope of a secant. To avoid additional solver calls, the required data pairs are formed from the intermediate results ˜x id and ˜R i of the k coupling iterations already performed for the current time step. More precisely, they are stored inthe input matrix V k ∈ R m × k and the output matrix W k ∈ R m × k [12]: V k = (cid:2) ∆ R , ∆ R , · · · , ∆ R kk − (cid:3) with ∆ R ji = ˜R ( ˜x jd ) − ˜R ( ˜x id ) = R ( x jd ) − R ( x id ) , (19a) W k = (cid:2) ∆ ˜x , ∆ ˜x , · · · , ∆ ˜x kk − (cid:3) with ∆ ˜x ji = ˜x jd − ˜x id . (19b)6ith the collected data, an approximation (cid:98) J − ≈ J − R ∈ R m × m of the inverse Jacobian can be formulated via (cid:98) J − V k = W k . (20)Since the number of data pairs stored in W k and V k typically is much smaller than the number of structural degreesof freedom at the FSI interface, i.e., k << m , the linear system of equations (20) is underdetermined. The existence ofa unique solution is ensured by demanding the minimization of the Frobenius norm [13]: (cid:107) (cid:98) J − (cid:107) F → min . (21)Together, equations (20) and (21) form a constrained optimization problem, which leads to an explicit form of theinverse Jacobian approximation [12, 13]: (cid:98) J − = W k (cid:16) V Tk V k (cid:17) − V Tk . (22)Inserting this inverse Jacobian approximation into the Newton update formula in Equation (17) yields a quasi-Newtonupdate scheme for the interface deformation, x k + d = U IQN ( ˜x kd ) = ˜x kd + (cid:98) J − (cid:16) − R ( x kd ) (cid:17) = ˜x kd + W k (cid:16) V Tk V k (cid:17) − V Tk (cid:16) − R ( x kd ) (cid:17)(cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) : = α = ˜x kd + W k α . (23)Defining the vector α = (cid:16) V Tk V k (cid:17) − V Tk (cid:16) − R ( x kd ) (cid:17) exploits that the inverse Jacobian is not needed explicitly here, butonly its product with the residual vector. An e ffi cient way of computing α ∈ R k is to solve the least-squares problem[11, 12] min α ∈ R k (cid:107) V k α + R ( x kd ) (cid:107) , (24)e.g., using a QR decomposition via Householder reflections. In the first coupling iteration, a relaxation step is used,as no data pairs are available yet.The approach presented so far is referred to as interface quasi-Newton inverse least-squares (IQN-ILS) method[11]; it can be interpreted as the basis of the other IQN variants discussed in this work. So far, the inverse Jacobian of the residual operator is approximated only based on information gathered in thecoupling iterations of the current time step. In principle, however, the e ffi ciency of IQN methods can be significantlyincreased by incorporating data from previous time steps as well. The most straightforward way to do so is to explicitlyinclude the data pairs of past time steps in the input and output matrices V k and W k [11].Unfortunately, apart from increasing costs for handling big data sets, a growing number of data pairs stored inthese matrices entails two major problems [1, 13, 40]: (1) The matrix V k might be very close to rank-deficient dueto (almost) linear dependent columns; related to that, the condition number of the least-squares problem quicklyincreases; (2) information gathered in di ff erent time steps might be contradictory. Together, these drawbacks carry therisk of preventing the system from being numerically solvable at all.An obvious remedy is to incorporate only the data of the q most recent time steps. While this approach is in factcapable of yielding a superior convergence speed, its major shortcomings are the risk of rank deficiency and in a goodchoice for the parameter q being problem-dependent [1, 12]. One way to mitigate these drawbacks is by employing afiltering technique [14, 15]. An alternative way of reusing data from previous time steps was introduced by Bogaers et al. [1] and developedfurther by Lindner et al. [13]: The interface quasi-Newton inverse multi-vector Jacobian (IQN-IMVJ) methodcombines the IQN approach with the idea of Broyden’s method: In a Newton iteration, one might expect the Jacobiansevaluated in subsequent iterations to be similar in some sense. Therefore, Broyden’s method limits the typically costly7acobian evaluation to the first iteration only; after that, the new Jacobian is instead approximated by a rank-one updateof the previous one. For the interface quasi-Newton framework, this concept is adopted by formulating the currentinverse Jacobian approximation n + (cid:98) J − as an update of the one determined in the previous time step n (cid:98) J − : n + (cid:98) J − = n (cid:98) J − + ∆ (cid:98) J − , (25)where ∆ (cid:98) J − denotes the update increment. Introducing this concept changes the constrained optimization problemdiscussed in Section 4.1 to (cid:107) ∆ (cid:98) J − (cid:107) F = (cid:107) n + (cid:98) J − − n (cid:98) J − (cid:107) F → min subject to ∆ (cid:98) J − V k = W k − n (cid:98) J − V k . (26)This system allows for a quite descriptive interpretation: Fitting Broyden’s idea, the di ff erence between the inverseJacobians of the current and the previous time step is minimized, rather than the norm of the approximation itself.Again, the result is an explicit approximation of the current inverse Jacobian [13]: n + (cid:98) J − = n (cid:98) J − + (cid:16) W k − n (cid:98) J − V k (cid:17) (cid:16) V Tk V k (cid:17) − V Tk (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) : = Z k = n (cid:98) J − + (cid:16) W k − n (cid:98) J − V k (cid:17) Z k . (27)In contrast to the IQN-ILS approach, the IQN-IMVJ method explicitly updates the Jacobian in every coupling iterationvia Equation (27). Therein, the matrix Z k = (cid:16) V Tk V k (cid:17) − V Tk ∈ R k × m is determined by solving the least-squares problemmin z j ∈ R k (cid:107) V k z j − ˆ e j (cid:107) (28)for each column z j with the j -th unit vector ˆ e j , again using a Householder QR decomposition. After n + (cid:98) J − has beenupdated, the quasi-Newton step is performed: x k + d = U IQN ( ˜x kd ) = ˜x kd − n + (cid:98) J − R k , where R k : = R ( x kd ) (29)Since the inverse Jacobian is initialized to zero, i.e., (cid:98) J − = , the first time step is equivalent to the IQN-ILS method.Note that the matrices W k and V k only contain data collected in the current time step; because, rather than ex-plicitly including data pairs collected in previous time steps, the IQN-IMVJ approach incorporates them implicitly inform of the previous Jacobian approximation n (cid:98) J − . This implicit reutilization of past information entails some signif-icant advantages [1, 13]: (1) Since the matrices W k and V k are not a ff ected, neither is the condition number of theleast-squares problem; (2) it renders any tuning of (strongly) problem-dependent parameters obsolete; (3) since pastinformation is matched in a minimum norm sense only, it is automatically less emphasized than that of the currenttime step. This avoids the risk of relying on outdated or contradicting data.These benefits come at the price of an explicit representation of the inverse Jacobian approximation, that is oftenvery expensive to store and handle, as the entailed complexity is growing quadratically with the problem scale. In this paper, we introduce an enhancement of the IQN-IMVJ approach tackling its main issue, i.e., the high costrelated to the explicit Jacobian approximation. By successively replacing the expensive parts, this section derives thenew interface quasi-Newton implicit multi-vector least-squares (IQN-IMVLS) method.As a first step, the explicit use of the current inverse Jacobian approximation n + (cid:98) J − is eliminated: By insertingthe Jacobian update (27) into the quasi-Newton step, the vector α = (cid:16) V Tk V k (cid:17) − V Tk ( − R k ) can be identified in analogyto the IQN-ILS approach, changing the update formula to x k + d = ˜x kd − n + (cid:98) J − R k = ˜x kd − n (cid:98) J − R k + (cid:16) W k − n (cid:98) J − V k (cid:17) (cid:16) V Tk V k (cid:17) − V Tk ( − R k ) (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) : = α = ˜x kd − n (cid:98) J − R k + (cid:16) W k − n (cid:98) J − V k (cid:17) α . (30)8y reducing the least-squares problem back to Equation (24), this drastically reduces the computational cost. Thepast inverse Jacobian n (cid:98) J − ∈ R m × m , however, is still explicitly required: Once a time step has converged, it is updatedto the next time level by Equation (27); the associated cost is quadratic in m . Beyond that, n (cid:98) J − is involved in thequasi-Newton formula in Equation (30). In particular, it is needed for the potentially very costly matrix product in B k : = ( W k − n (cid:98) J − V k ), which has a complexity of O ( m k ). As mentioned before, however, the matrices W k and V k aresuccessively built up only by appending new columns. Taking this into account, we can reformulate B k = W k − n (cid:98) J − V k = (cid:104) W k − , ∆ ˜x kk − (cid:105) − n (cid:98) J − (cid:104) V k − , R k − R k − (cid:105) = (cid:104) W k − − n (cid:98) J − V k − (cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) : = B k − , ∆ ˜x kk − − n (cid:98) J − R k (cid:124) (cid:32) (cid:123)(cid:122) (cid:32) (cid:125) : = b k + n (cid:98) J − R k − (cid:124) (cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32) (cid:125) : = b k − (cid:105) . = (cid:104) B k − , ∆ ˜x kk − − b k + b k − (cid:105) (31)As a consequence, restoring the terms B k − : = W k − − n (cid:98) J − V k − ∈ R m × k and b k − : = n (cid:98) J − R k − ∈ R m from the previousiteration in fact allows to reduce the number of matrix-vector products with the inverse Jacobian to one per couplingiteration, that is b k = n (cid:98) J − R k .This Jacobian-vector product remains the only operation in Equation (30) involving an O ( m ) complexity. For anincreasing number of structural degrees of freedom at the interface m , this term’s cost therefore may become dominantand slow down the overall procedure (see Remark 4).In order to tackle this issue, we introduce an alternative purely implicit formula for evaluating b k = n (cid:98) J − R k without any explicit previous Jacobian n (cid:98) J − . Note that the explicit Jacobian update will become obsolete as a directconsequence. Unraveling the recursion in Equation (27), the update of the inverse Jacobian can be reformulated to n + (cid:98) J − = (cid:98) J − n (cid:89) i = (cid:16) I − V ik Z ik (cid:17) + n (cid:88) i = W ik Z ik n (cid:89) j = i + (cid:16) I − V jk Z jk (cid:17) , (32)see Appendix A for the proof. Therein, the upper index of the matrices W ik , V ik , and Z ik refers to the time step t i → t i + they were determined in. Considering the initial Jacobian approximation, (cid:98) J − = , this expression simplifies to n + (cid:98) J − = n (cid:88) i = W ik Z ik n (cid:89) j = i + (cid:16) I − V jk Z jk (cid:17) . (33)While being mathematically equivalent, this formulation entails two main advantages: First, the matrix-vector product n (cid:98) J − R k can be determined in an implicit manner via n (cid:98) J − R k = n − (cid:88) i = W ik Z ik n − (cid:89) j = i + (cid:16) R k − V jk Z jk R k (cid:17) , (34)without an O ( m ) complexity, and hence potentially cheaper. However, the cost for evaluating the expression inEquation (34) obviously increases with the number of processed time steps. This issue is addressed using the secondadvantage: While data from past time steps is still incorporated in an implicit manner, with all the associated advan-tages discussed in Section 4.3, the new formulation explicitly identifies the contribution of each time step, in the formof the matrices V ik , W ik , and Z ik . As a consequence, this update technique allows to incorporate only the q most recenttime steps for the previous Jacobian approximation as a way of limiting the cost of the matrix-vector product, whichthen reads n (cid:98) J − R k ≈ n − (cid:88) i = n − q W ik Z ik n − (cid:89) j = i + (cid:16) R k − V jk Z jk R k (cid:17) . (35)In fact, by taking advantage of the repetition of terms in the product (see Appendix B), this expression can be evaluatedin O ( m ¯ k q ), where ¯ k denotes the average number of coupling iterations per time step. This step can be justified by thefact that a certain time step’s contribution to the previous Jacobian approximation is gradually becoming less and lessimportant the further the simulation progresses. 9t a first glance, one might argue that the parameter q reintroduces exactly the problems the multi-vector approachwas designed to avoid in the first place, i.e., the dependency of the IQN-ILS method on the number of reused timesteps. However, this issue arose from the explicit incorporation of past data, e.g., due to high condition numbers andlinear-dependent columns in V k . Since the IQN-IMVLS approach reuses data in an implicit manner, it does not su ff erfrom these drawbacks. Instead, the method’s quality in general benefits from increasing the number of past time stepstaken into account, as the limit of reusing all steps is analytically equivalent to the explicit Jacobian multiplication.For this variant, the matrix Z k = (cid:16) V Tk V k (cid:17) − V Tk ∈ R k × m still has to be determined and stored once after each timestep. Since the complexity of computing it via the m least-squares problems (28) is O ( m ¯ k ) for the Householder QRapproach, i.e., growing quadratically with m , we use a matrix inversion of V Tk V k via a LU decomposition using partialpivoting with row interchanges instead [41]. While being slightly less robust to bad conditioning, the big advantageis that it requires O ( m ¯ k ) operations and therefore scales linearly with the problem size.Combining all this, the modified multi-vector update completely avoids any O ( m )-terms that might slow downthe IQN-IMVLS approach for large systems, where m >> ¯ k q holds, see Table 1. A pseudo-code realization of thepurely implicit IQN-IMVLS method is outlined in Algorithm 1. Increment step: ∆ x = (cid:16) W k − n (cid:98) J − V k (cid:17) α Operation Expression Explicit Implicit
Least-squares problem α m ¯ k m ¯ k via Householder QRMatrix-vector product b k = n (cid:98) J − R k m m ¯ k q Compute new column B k = (cid:104) B k − , ∆ ˜x kk − − b k + b k − (cid:105) m m Matrix-vector product B k α m ¯ k m ¯ k Jacobian update: n + (cid:98) J − = (cid:16) W k − n (cid:98) J − V k (cid:17) Z k = B k Z k Operation Expression Explicit Implicit
Least-squares problem Z k m ¯ k – via Householder QRMatrix-matrix product B k Z k m ¯ k – Determine and store Z k Operation Expression Explicit Implicit
Matrix inversion via (cid:16) V Tk V k (cid:17) − – m ¯ k LU decompositionMatrix-matrix product (cid:16) V Tk V k (cid:17) − V Tk – m ¯ k Table 1: Complexities of the suboperations required for the IQN-IMVLS update method with an explicit previous Jacobian compared to the purelyimplicit version.
In direct analogy to the computational complexity, the memory requirements, too, are no longer scaling quadrati-cally but linearly with the problem size: Although the matrices V ik , W ik ∈ R m × ¯ k as well as Z ik ∈ R ¯ k × m have to be storedfor the q most recent time steps, the required amount of storage is much smaller than for the explicit Jacobian as longas m >> ¯ kq holds.Of course, the e ff ectiveness of the purely implicit IQN-IMVLS update strongly depends on this ratio of m and ¯ k q ,so that the explicit Jacobian approximation might still be the better option for small systems. However, the assumption m >> ¯ k q is not very restrictive for common application scales.10 ime Step Loop: for n = , · · · do Initialization: W n = [ ], V n = [ ], B = [ ]First iteration: ˜x d = H ( x d )Form residual: R = ˜x d − x d if n == then Under-relaxtion step: x d = ω ˜x d + (1 − ω ) x d else Previous inverse Jacobian times residual (implicitly): b = (cid:80) n − i = n − q W ik Z ik (cid:81) n − j = i + (cid:16) R − V jk Z jk R (cid:17) IQN-update with Jacobian product: x d = ˜x d − b endCoupling Loop: for k = , · · · until convergence do Perform coupling iteration: ˜x kd = H ( x kd )Form residual: R k = ˜x kd − x kd Append new column to input matrix: V nk = [ V nk − , ∆ R kk − ] with ∆ R ji = R j − R i Append new column to output matrix: W nk = [ W nk − , ∆ ˜x kk − ] with ∆ ˜x ji = ˜x jd − ˜x id Previous inverse Jacobian times residual (implicitly): b k = (cid:80) n − i = n − q W ik Z ik (cid:81) n − j = i + (cid:16) R k − V jk Z jk R k (cid:17) Build B k restoring terms from previous iteration: B k = [ B k − , ∆ ˜x kk − + b k − − b k ]Solve least-squares problem: min (cid:107) V nk α + R k (cid:107) Implicit update step: x k + d = ˜x kd − B k α end Determine Z nk via LU decomposition: Z nk = (cid:16) V Tk V k (cid:17) − V Tk Store for future time steps: Z nk , W nk , V nk → Store end
Algorithm 1: Pseudo-code of the interface quasi-Newton implicit multi-vector least-squares (IQN-IMVLS) method.
Remark 2:
Aside from suggesting improvements similar to using an explicit past Jacobian within the IQN-IMVLSmethod, Scheufele and Mehl [12] also derived a multi-vector variant with linear complexity. Therein, the past inverseJacobian is represented based on the matrices Z ik and B ik from past time steps, i.e., n (cid:98) J − = n − (cid:88) i = B ik Z ik . (36)While the idea to avoid the costly explicit Jacobian via a sum formulation is similar to the IQN-IMVLS method, thereis one essential di ff erence: Due to the recursive definition of the multi-vector Jacobian approximation, each matrix B ik = W ik − i (cid:98) J − V ik contains information from the first i time steps. By implication, this means that the contribution ofa time step j is contained in all subsequent matrices B jk , B j + k , B j + k , and so forth. As a consequence, it is impossibleto drop old time steps from the approximation while keeping more recent ones, as the IQN-IMVLS method does.Instead, Scheufele and Mehl divide the simulation into chunks of several time steps; after each of these chunks,the sum in equation (36) has to be reset. Since a plain erasure of old data at every restart dramatically impairsthe e ffi ciency of the multi-vector approach, focus is put on di ff erent restart options: In particular, they introduce the multi-vector Jacobian restart singular value decomposition (IQN-MVJ-RS-SVD) method, in which a truncated SVDaccounts for the dropped data. It introduces two main parameters: While the influence of the chunk size M is typicallyless marked, the method’s accuracy and e ffi ciency strongly depends on the threshold ε svd for cutting o ff insignificantsingular values. However, choices around ε svd ≈ .
01 were shown to be suitable for various test cases [12]. For anin-depth discussion of the method refer to Scheufele and Mehl [12, 42].In contrast to this restart-based approach, the IQN-IMVLS method always considers the q most recent time steps.In combination with its arguably simpler implementation, it therefore represents a straightforward way of achieving alinear complexity without any need for restart techniques.11 emark 3: In the beginning of a new time step, the input and output matrix of the IQN-IMVLS method rely onvery few data pairs only. To improve the significance of the least-squares problem, it can be beneficial to includethe information of the most recent time step explicitly in the matrices V k and W k – in addition to the implicit reuti-lization of past data. The e ff ect of this option will be included in the discussion of the numerical test cases in Section 5. Remark 4:
The computational cost of the interface quasi-Newton approaches discussed in this work increases withthe number of structural degrees of freedom at the FSI interface. To assess the significance of this cost in comparisonto the solver calls, the complexity of the structural solution is therefore of particular interest. Nowadays, the structuralsubproblem is typically solved using some finite-element variant, such as isogeometric analysis. These procedurescan be subdivided into two main tasks: (1) The assembly of the system matrices is done in an element-by-elementmanner and hence scales with the number of elements n el , i.e., O ( n el ). Moreover, each of the n do f degrees of freedomcauses b nonzero entries in the sparse system matrix [43]. As this number b << n do f primarily depends on theelement type this yields O ( n do f ). For common meshes, O ( n el ) ≈ O ( n do f ) is a reasonable estimate [44]. (2) Usually,an iterative solver is employed for solving the resulting sparse linear system of equations. In the ideal case, theassociated computational complexity is O ( n do f n iter ) [45, 46]. Since increasing the number of unknowns n do f typicallybrings along a moderately growing number of solver iterations n iter , the total complexity of the structural subproblemis expected to be superlinear, but significantly lower than quadratic.
5. Numerical Results
The first test case considers an elastic cylindrical tube that is filled with an incompressible fluid, as depicted inFigure 2a. Caused by a short excitation in the beginning, a pressure pulse propagates through the pipe structure. Whilethe configuration is inspired by similar test cases discussed, among others, by Degroote et al. [7] and Lindner et al.[13], the prescribed pressure pulse is increased by a factor of ten.
Incompressible FluidElastic Wall R = . m L = . ms = . m (a) Configuration of the elastic tube test case. i n l e t p r e ss u r e i n P a time in msp max =13322 Pa (b) Pressure pulse imposed on the inlet.Figure 2: Elastic tube test case The tube has a length of L = . m , an inner radius of R = . m , and a wall thickness of s = . m . Theelastic structure is characterized by the density ρ s = kgm , its Young’s modulus E = . · kgm s , and a Poissonratio of ν = .
3; the fluid’s density is ρ f = kgm and its dynamic viscosity µ f = . kgm s . Based on the densityratio of ρ f /ρ s = . . Pa is prescribedon one end, the other one is excited by a short pressure pulse with a peak of 13322 . Pa in the beginning of thesimulation, see Figure 2b. After that, its boundary pressure is held fixed at 0 . Pa as well.The fluid domain is discretized by 19193 tetrahedral elements with 4231 nodes per time level, i.e., a total of8462 nodes due to the space-time formulation. The mesh is locally refined in the vicinity of the walls and near the12rescribed pressure pulse. The elastic structure is clamped at both ends; it is represented by 16 × =
640 nonlinearisogeometric shell elements defined by a quadratic NURBS. This spline surface has 714 control points and 4284degrees of freedom. The simulation runs till T = ms in 80 time steps of size ∆ t = . ms . The convergence of thecoupling scheme is detected by a combination of an absolute bound ε abs = − and a relative criterion ε rel = − for the norm of the fixed-point residual R k . As stated above, the excitation causes a pressure pulse propagating through the elastic tube, which is depicted inFigure 3 for three sample time levels. The snapshots show that the moving pressure peak is accompanied by a largewidening of the structure. While retaining its basic profile on its way through the pipe, the pulse clearly exhibits somedi ff usive flattening. Qualitatively, the numerical results agree with both physical expectations and the discussionsof similar test cases in literature [1, 11, 13]. Note that this solution is independent (within the chosen convergencetolerance) from the update technique employed. (a) t = . s : The peak detaches from the inlet.(b) t = . s : The pressure wave propagates through the structure.(c) t = . s : Eventually, it reaches the outflow end of the tube.Figure 3: Illustration of the pressure peak running through the tube. Aside from the pressure field, the deformation of the fluid mesh is shown. The focus of this test case is on the e ffi ciency of di ff erent updating schemes. A comparison is provided in Table 2:Aside from the average number of coupling iterations per time step, it lists the relative speed-up in terms of couplingiterations and runtime with respect to a constant under-relaxation; lastly, the percentages of computation time spent forthe interface quasi-Newton update are indicated. The first observation is that, although Aitken’s dynamic relaxationconverges much faster than the constant under-relaxation, neither of them is really feasible for this test case, since theyrequire about 100 or 26 coupling iterations per time step, respectively. The problem is the high added-mass instability,requiring very small relaxation factors; in this case, ω ≈ .
05 was found to be the highest value precluding divergence.As opposed to this, all employed interface quasi-Newton variants entail a substantial increase in performance –both in terms of reducing the number of coupling iterations and the runtime. In particular, the incorporation of datafrom past time steps proves very beneficial. With this, using an IQN approach reduces the total runtime comparedto a constant relaxation by up to 92%; perhaps even more remarkable is the decrease by about 70% with respect toAitken’s method.The results for the IQN-ILS variant identify its dependency on the number of explicitly reused time steps q : Awayfrom the optimum of this test case, q ≈
5, including either too few or too many previous steps compromises thee ff ectiveness; in fact, the coupling scheme was even diverging for values of q higher than about 25. Note in thiscontext that no filtering technique is applied within this work.13ethod Average Coupling Relative Coupling Relative IQN-Time / Iterations Iterations Runtime RuntimeUnder-relaxation, ω = .
05 100.39 100.00 % 100.00 % -Aitken’s relaxation 26.38 26.28 % 26.31 % -IQN-ILS ( q =
0) 14.96 14.90 % 14.93 % 0.087 %IQN-ILS, q = q = q =
10 8.66 8.63 % 9.03 % 1.042 %IQN-ILS, q =
20 9.76 9.72 % 10.14 % 2.890 %IQN-IMVJ 8.11 8.08 % 9.53 % 11.736 %IQN-MVJ-RS-SVD, M = ε svd = .
01 8.09 8.06 % 8.41 % 0.370 %IQN-IMVLS, explicit Jacobian 8.14 8.10 % 8.50 % 1.263 %IQN-IMVLS, q = q =
10 8.28 8.24 % 8.68 % 0.078 %IQN-IMVLS, q =
20 7.95 7.92 % 8.28 % 0.089 %IQN-IMVLS, q =
50 8.05 8.02 % 8.34 % 0.113 %IQN-IMVLS, q =
80 8.05 8.02 % 8.35 % 0.120 %IQN-IMVLS, q =
80, 1 explicit step 7.76 7.73 % 8.02 % 0.253 %
Table 2: Comparison of di ff erent update schemes for the pressure tube test case. Against this backdrop, the implicit incorporation of past data via the previous inverse Jacobian deserves specialnotice: Despite not requiring any problem-dependent parameters, it enables the IQN-IMVJ variant to keep up withthe optimum of the IQN-ILS approach in terms of coupling iterations. However, the actual runtime indicates its majordrawback: While the computational cost of the IQN-ILS method barely exceeds 1 % of the total runtime for moderatechoices of q , handling the explicit Jacobian approximation raises the e ff ort for the IQN-IMVJ method to 11 .
74 %;hence, its cost is by no means negligible against the solver calls.That is exactly what the IQN-IMVLS variant is designed to provide a remedy for. Since the analytic form ofthe quasi-Newton update stays untouched, the required coupling iterations of the IQN-IMVJ and the IQN-IMVLSapproach with an explicit Jacobian show only an insignificant di ff erence arising from numerical rounding. The IQN-IMVLS method, however, employs more advantageous numerical techniques and a beneficial storage of already deter-mined quantities (as discussed in Section 4.4) to reduce the update’s computational cost. In particular, the number ofoperations involving an explicit Jacobian approximation is reduced. The results in Table 2 prove that these adjustmentsentail a significant speed-up of the multi-vector IQN approach – in this case by a factor of ten.Moreover, the IQN-IMVLS method also provides a purely implicit version avoiding any explicit representation ofthe Jacobian approximation. This variant reintroduces the number of past time steps q to be incorporated; in contrastto the explicit reutilization, however, this parameter is not very problem-dependent or critical, as high values of q donot entail any numerical problems. Instead, the quality of the update scheme in general benefits from an increasing q , as confirmed by Table 2: Although there seems to be a marginal optimum at around q =
20 for this test case, theconvergence speed afterwards aligns more and more with the one observed for including all past time steps ( q = q =
50. In combination with the computational e ff ort,which stays at around 0 . ff ectiveness, the restart-based IQN-MVJ-RS-SVD method significantly reduces the costof the multi-vector approach, too. Although being slightly more expensive than the new algorithm for this test case,the observed di ff erence is small enough to be blurred by implementation details and parameter choices to some extent;in relation to the overall simulation time, both methods come at negligible cost.When it comes to including the most recent time step explicitly in the IQN-IMVLS method as suggested in Re-mark 3, the last row in Table 2 confirms that this option does in fact slightly reduce the number of coupling iterationswithout being too expensive. The e ff ect, however, is not very marked.14s indicated before, the main advantage of the purely implicit update, i.e., its linear complexity, becomes moredistinct for larger problem scales. In general, the cost of the discussed interface quasi-Newton concept depends onthe number of structural degrees of freedom at the FSI interface. As the elastic tube is modeled by shell elements,however, in this test case the FSI interface is equivalent to the whole structural domain. For shell analysis, a distinctionbetween all structural degrees of freedom and those at the interface is redundant.Therefore, Figure 4 plots the absolute time spent for the multi-vector IQN variants as well as the IQN-ILS methodwith q = ff erent methods,both axes are scaled logarithmically. The first observation is the exploding cost of the IQN-IMVJ method, caused byits extensive usage of the explicit Jacobian approximation. While the IQN-IMVLS method with an explicit inverseJacobian again significantly reduces the numerical e ff ort, it does not keep the cost from growing quadratically withthe problem scale. In contrast, the purely implicit IQN-IMVLS version is not only much cheaper; the plot clearlypoints out its linear scaling with the interface resolution. Moreover, for this test setting it outperforms both the IQN-MVJ-RS-SVD and the IQN-ILS method, which also show a linear complexity. ti m e s p e n t f o r I QN upd a t e [ i n s ec ] structural degrees of freedom at the FSI interfaceIMVJIMVLS, expl. Jac.ILS, q=5MVJ-RS-SVDIMVLS, q=80 Figure 4: Double-logarithmic plot of the computational cost for di ff erent interface quasi-Newton variants in dependence on the interface resolution. While Figure 4 compares the e ffi ciency of di ff erent quasi-Newton schemes, it does not allow any judgment on towhich extent their computational cost is negligible compared to the solver calls. Therefore, Table 3 puts them intocontext: Taking into account that the fluid mesh does not have a direct influence on the IQN approach and that thestructural domain is equivalent to the FSI interface due to the usage of shell elements, its last column relates the timespent for the quasi-Newton methods to that required for the structural solution; it considers three di ff erent levels ofrefinement of the structural discretization, i.e., 9234, 22644, and 44228 degrees of freedom. Moreover, Table 3 liststhe average numbers of coupling iterations per time step, which stay close to 8 for all tests conducted. Nevertheless,the investigated methods show vast di ff erences regarding their computational cost.Again, the most striking observation is the very high cost of the IQN-IMVJ method: For the three refinementlevels, the runtime it accounts for ranges from 27 .
35% to 83 .
58% compared to the time required for the structuralsolution, and therefore seriously slows down the whole simulation. Moreover, this percentage is growing with theproblem scale, indicating that the quadratic complexity causes the cost of the IQN-IMVJ approach to be increasingfaster than that of solving the structural subproblem. As discussed before, the IQN-IMVLS method with an explicitJacobian approximation significantly reduces the numerical e ff ort, but not the method’s complexity; hence, its relative15ethod Average Coupling IQN / StructureIterations in %9324 Structural Degrees of Freedom (at the FSI Interface)IQN-ILS, q = M = ε svd = .
01 7.99 0.35 %IQN-IMVLS, explicit Jacobian 8.05 2.66 %IQN-IMVLS, q =
80 7.99 0.12 %IQN-IMVLS, q =
80, 1 explicit step 7.48 0.22 %22644 Structural Degrees of Freedom (at the FSI Interface)IQN-ILS, q = M = ε svd = .
01 8.09 0.32 %IQN-IMVLS, explicit Jacobian 7.99 5.08 %IQN-IMVLS, q =
80 7.86 0.09 %IQN-IMVLS, q =
80, 1 explicit step 7.79 0.19 %42228 Structural Degrees of Freedom (at the FSI Interface)IQN-ILS, q = M = ε svd = .
01 8.04 0.27 %IQN-IMVLS, explicit Jacobian 7.88 7.31 %IQN-IMVLS, q =
80 7.97 0.06 %IQN-IMVLS, q =
80, 1 explicit step 8.20 0.15 %
Table 3: Comparison of IQN variants for three refinement levels. Due to the shell model, all structural degrees of freedom are at the FSI interface. cost is still growing with the number of degrees of freedom, in this case from 2 .
66% to 7 . The second example is a cylindrical deformable tank partially filled with an incompressible fluid. While Figure 5illustrates the basic configuration, the test case parameters are listed in Table 4.The simulation is started from the steady state where the fluid is at rest and the deformation has adjusted to thehydrostatic pressure. During the simulation, the system is excited in horizontal direction by a periodic movementof the tank bottom with the prescribed displacement ∆ x ( t ) = A sin(2 π f t ). Due to inertia, this excitation results incomplex and highly intertwined motions of the deforming tank structure and the liquid sloshing inside.The evolving fluid domain is discretized by a space-time mesh of 16243 elements and 4024 nodes per time level;since it is adapted to the free-surface via interface tracking [18, 25], the tank wall and its bottom involve slip-conditionsfor both the fluid and the mesh. The cylindrical tank wall is modeled by 32 × =
960 nonlinear isogeometric shellelements of degree 2 defined by a NURBS with 1184 control points. The simulation runs for 1600 time steps of size ∆ t = . ε abs = − and ε rel = − .16 Fluid
H R
Elastic WallMoving Base
Figure 5: Sloshing tank test case.
Field Parameter Symbol Value
Fluid Density ρ f kgm Dynamic viscosity µ f . kgm s Structure Density ρ s kgm Young’s modulus E . + kgm s Poisson ratio ν . R m Tank height H m Filling level h m Wall thickness s . m Excitation Amplitude A . m Frequency f . Hz Table 4: Test case parameters.
To give an impression of the simulation results, Figure 6 depicts three di ff erent snapshots of the sloshing tank.Both the free surface and the tank structure exhibit large displacements. (a) t = . s (b) t = . s (c) t = . s Figure 6: Snapshots of the sloshing tank test case. The colored arrows visualize the velocity field.
The e ffi ciency of di ff erent update techniques is compared in Table 5, both in terms of coupling iterations and theircomputational cost. All in all, the results are in very good agreement with those of the elastic tube scenario: Whileneither a constant nor Aitken’s dynamic relaxation converge fast enough to be a reasonable option, all interface quasi-Newton variants drastically speed up the coupling scheme. Again, incorporating information from past time steps isextremely beneficial. The IQN-ILS method performs very well for this test case, showing a quite flat optimum ofchoosing the parameter q ; the cost of solving the least-squares problem, however, is clearly increasing with q .Although being slightly less e ffi cient for this test case, the implicit reutilization of the multi-vector methods almostkeeps up with the IQN-ILS approach in terms of the required coupling iterations. While the IQN-IMVJ variant paysfor this fast convergence with the expensive explicit Jacobian, the computational cost of both the IQN-IMVLS and theIQN-MVJ-RS-SVD method is comparable to that of the IQN-ILS method, thanks to the linear complexity. In directcomparison, the new IQN-IMVLS method outruns the restart-based variant for this test case; however, this statementis again best put into the perspective of the specific implementation and the chosen parameters.Increasing the value q brings the convergence speed closer to that of the explicit multi-vector approach. Althougha rather high value is required here to reach this speed, the implicit Jacobian product proves to be very cheap as thecomputational cost stays non-critical even for the maximum of reusing all time steps, i.e., q = / RuntimeUnder-relaxation, ω = . q =
0) 17.90 0.08 %IQN-ILS, q = q = q =
10 4.45 0.39 %IQN-ILS, q =
20 4.17 0.64 %IQN-ILS, q =
50 4.20 2.29 %IQN-ILS, q =
100 4.40 8.58 %IQN-IMVJ 5.10 11.33 %IQN-MVJ-RS-SVD, M = ε svd = .
01 5.18 1.28 %IQN-IMVLS, explicit Jacobian 5.12 1.67 %IQN-IMVLS, q =
50 7.14 0.08 %IQN-IMVLS, q =
100 6.23 0.12 %IQN-IMVLS, q =
200 5.56 0.20 %IQN-IMVLS, q =
500 5.05 0.37 %IQN-IMVLS, q = q = q = q = Table 5: Comparison of di ff erent update techniques for the sloshing tank test case. Altogether, the test case proves that the IQN-IMVLS method competes with the IQN-ILS approach not only interms of the convergence speed, but also regarding the computational cost.
6. Conclusion
This paper presents a novel interface quasi-Newton method improving the partitioned simulation of fluid-structureinteraction. Building on the inverse multi-vector Jacobian (IQN-IMVJ) variant, it overcomes the main drawback: Asthe cost of handling an explicit inverse Jacobian approximation increases quadratically with the problem size, theIQN-IMVJ method significantly slows down large-scale simulations. The new interface quasi-Newton implicit multi-vector least-squares (IQN-IMVLS) method, in contrast, systematically avoids any explicit Jacobian approximation,so that the multi-vector update is realized with linear complexity. As a consequence, its computational cost staysnegligible even and particularly for industrial-scale applications.The e ff ectiveness of the IQN-IMVLS method is confirmed by the numerical test cases in Section 5. Although thenumber of required coupling iterations stays virtually unchanged compared to the IQN-IMVJ variant, the results showa substantial reduction of computational cost. In accordance to the discussion above, this di ff erence becomes moreand more distinct for finer resolutions of the structural discretization. One major strength of the multi-vector conceptis that its implicit incorporation of past data does not rely on problem-dependent parameters – and yet keeps up withan optimal parameter choice of the interface quasi-Newton implicit least-squares (IQN-ILS) method, regarding therequired coupling iterations. While the IQN-IMVJ approach had to pay a high price in computational cost for thiscapability, the IQN-IMVLS method is comparable to the IQN-ILS variant concerning the runtime.Although the IQN-IMVLS variant reintroduces the number of reused time steps q , this parameter is not problem-atic at all as the quality of the quasi-Newton update in general benefits from including more steps.In conclusion, the new IQN-IMVLS method combines the advantages of the IQN-IMVJ and the IQN-ILS variantwithout their individual drawbacks. Coping without critical parameters, it succeeds in realizing the favorable implicitreutilization of past data with a computational complexity that grows linearly with the problem size.18 ppendix A. Proof of New Jacobian Update Via mathematical induction, this section proves the equality of the two Jacobian updates discussed in Section 4.4: n + (cid:98) J − = n (cid:98) J − + (cid:16) W nk − n (cid:98) J − V nk (cid:17) Z nk and n + (cid:98) J − = (cid:98) J − n (cid:89) i = (cid:16) I − V ik Z ik (cid:17) + n (cid:88) i = W ik Z ik n (cid:89) j = i + (cid:16) I − V jk Z jk (cid:17) .
1. Induction base case: For n = (cid:98) J − = (cid:98) J − (cid:89) i = (cid:16) I − V ik Z ik (cid:17) + (cid:88) i = W ik Z ik (cid:89) j = i + (cid:16) I − V jk Z jk (cid:17)(cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) (cid:81) qp = p > q = (cid:98) J − (cid:16) I − V k Z k (cid:17) + W k Z k = (cid:98) J − + (cid:16) W k − (cid:98) J − V k (cid:17) Z k .
2. Inductive step: If the equality is satisfied for any n , it can be shown to hold for n + n + (cid:98) J − = n (cid:98) J − + (cid:16) W nk − n (cid:98) J − V nk (cid:17) Z nk = n (cid:98) J − (cid:16) I − V nk Z nk (cid:17) + W nk Z nk = (cid:98) J − n − (cid:89) i = (cid:16) I − V ik Z ik (cid:17) + n − (cid:88) i = W ik Z ik n − (cid:89) j = i + (cid:16) I − V jk Z jk (cid:17) (cid:16) I − V nk Z nk (cid:17) + W nk Z nk = (cid:98) J − n (cid:89) i = (cid:16) I − V ik Z ik (cid:17) + n − (cid:88) i = W ik Z ik n (cid:89) j = i + (cid:16) I − V jk Z jk (cid:17) + W nk Z nk = (cid:98) J − n (cid:89) i = (cid:16) I − V ik Z ik (cid:17) + n − (cid:88) i = W ik Z ik n (cid:89) j = i + (cid:16) I − V jk Z jk (cid:17) + W nk Z nk n (cid:89) j = n + (cid:16) I − V jk Z jk (cid:17)(cid:124) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:123)(cid:122) (cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32)(cid:32) (cid:125) = = (cid:98) J − n (cid:89) i = (cid:16) I − V ik Z ik (cid:17) + n (cid:88) i = W ik Z ik n (cid:89) j = i + (cid:16) I − V jk Z jk (cid:17) . (cid:3) Appendix B. Implicit Jacobian Multiplication
Based on the new Jacobian update formulation introduced in Section 4.4 and proven in Appendix A, the productof the previous Jacobian approximation with a (residual) vector R k can be evaluated as n (cid:98) J − R k ≈ n − (cid:88) i = n − q W ik Z ik n − (cid:89) j = i + (cid:16) R k − V jk Z jk R k (cid:17) , (B.1)where q of the n past time steps are considered. Note that the initial choice (cid:98) J − = has been exploited. Thekey for an e ffi cient evaluation of this expression is to realize that the term (cid:16) I − V jk Z jk (cid:17) is needed for all i < j . As aconsequence, when looping over the previous time steps in reversed order, it is su ffi cient to update the product term (cid:81) nj = i + (cid:16) R k − V jk Z jk R k (cid:17) in every iteration i by one product – rather than evaluating it again and again. This way, theoverall computational complexity of the Jacobian-vector product evaluation reduces to O ( m ¯ kq ), where ¯ k is the averagenumber of coupling iterations per time step. For more details on the procedure and the cost of each step involved, apseudo-code realization is presented in Algorithm 2. Acknowledgement
We gratefully acknowledge the computing time granted by the J¨ulich-Aachen Research Alliance JARA-HPC.19mplicit Jacobian-Vector Product b k = n (cid:98) J − R k Initialization:
Auxiliary vector for Π -terms: a π = R k Result vector: b k = q time steps: for i = n − , · · · , n − q do Add contribution of time step i to result b k : b k = b k − W ik Z ik a π Update a π with Π -term of time step i : a π = a π − V ik Z ik a π end m → O ( m ) m → O ( m ) q iterations m (¯ k + → O ( m ¯ k )2 m (¯ k + → O ( m ¯ k )Total complexity: 4 m (¯ k + q → O ( m ¯ kq ) Algorithm 2: Pseudo-code for computing the product of the previous time step’s inverse Jacobian with a vector in an implicit manner. Thenumber of floating point operations and the resulting computational complexity are indicated in italic.
References [1] A. E. Bogaers, S. Kok, B. D. Reddy, T. Franz, Quasi-Newton Methods for Implicit Black-Box FSI Coupling, Computer Methods in AppliedMechanics and Engineering 279 (2014) 113–132.[2] C. F¨orster, Robust Methods for Fluid-Structure Interaction with Stabilized Finite Elements, Ph.D. thesis, University of Stuttgart, 2007.[3] C. F¨orster, W. A. Wall, E. Ramm, Artificial Added Mass Instabilities in Sequential Staggered Coupling of Nonlinear Structures and Incom-pressible Viscous Flows, Computer Methods in Applied Mechanics and Engineering 196 (2007) 1278–1293.[4] U. K¨uttler, W. A. Wall, Fixed-Point Fluid–Structure Interaction Solvers with Dynamic Relaxation, Computational Mechanics, Springer 43(2008) 61–72.[5] B. M. Irons, R. C. Tuck, A Version of the Aitken Accelerator for Computer Iteration, International Journal for Numerical Methods inEngineering 1 (1969) 275–277.[6] B. Gatzhammer, E ffi cient and Flexible Partitioned Simulation of Fluid-Structure Interactions, Ph.D. thesis, Technical University of Munich,2014.[7] J. Degroote, R. Haelterman, S. Annerel, P. Bruggeman, J. Vierendeels, Performance of Partitioned Procedures in Fluid–Structure Interaction,Computers & Structures 88 (2010) 446–457.[8] J. Degroote, A. Souto-Iglesias, W. Van Paepegem, S. Annerel, P. Bruggeman, J. Vierendeels, Partitioned Simulation of the Interaction betweenan Elastic Structure and Free Surface Flow, Computer Methods in Applied Mechanics and Engineering 199 (2010) 2085–2098.[9] J.-F. Gerbeau, M. Vidrascu, A Quasi-Newton Algorithm based on a Reduced Model for Fluid-Structure Interaction Problems in Blood Flows,ESAIM: Mathematical Modelling and Numerical Analysis 37 (2003) 631–647.[10] E. Van Brummelen, C. Michler, R. De Borst, Interface-GMRES (R) Acceleration of Subiteration for Fluid-Structure-Interaction Problems,Report DACS-05-001 (2005).[11] J. Degroote, K.-J. Bathe, J. Vierendeels, Performance of a New Partitioned Procedure Versus a Monolithic Procedure in Fluid-StructureInteraction, Computers & Structures 87 (2009) 793–801.[12] K. Scheufele, M. Mehl, Robust Multisecant Quasi-Newton Variants for Parallel Fluid-Structure Simulations—and Other Multiphysics Ap-plications, SIAM Journal on Scientific Computing 39 (2017) 404–433.[13] F. Lindner, M. Mehl, K. Scheufele, B. Uekermann, A Comparison of Various Quasi-Newton Schemes for Partitioned Fluid-Structure Inter-action, Proceedings of 6th International Conference on Computational Methods for Coupled Problems in Science and Engineering, Venice(2015) 1–12.[14] J. Degroote, J. Vierendeels, Multi-Solver Algorithms for the Partitioned Simulation of Fluid-Structure Interaction, Computer Methods inApplied Mechanics and Engineering 200 (2011) 2195–2210.[15] R. Haelterman, A. E. Bogaers, K. Scheufele, B. Uekermann, M. Mehl, Improving the Performance of the Partitioned QN-ILS Procedure forFluid–Structure Interaction Problems: Filtering, Computers & Structures 171 (2016) 9–17.[16] F. Yibin, R. Ogden, Nonlinear Elasticity: Theory and Applications, volume 283, 2001.[17] K.-J. Bathe, Finite Element Procedures, TBS, 1996.[18] N. Hosters, Spline-Based Methods for Fluid-Structure Interaction, Ph.D. thesis, RWTH Aachen University, 2018.[19] C. Braun, Ein modulares Verfahren f¨ur die numerische aeroelastische Analyse von Luftfahrzeugen, Ph.D. thesis, RWTH Aachen University,2007.[20] U. K¨uttler, C. F¨orster, W. A. Wall, A Solution for the Incompressibility Dilemma in Partitioned Fluid-Structure Interaction with Pure DirichletFluid Domains, Computational Mechanics, Springer 38 (2006) 417–429.[21] L. Pauli, Stabilized Finite Element Methods for Computational Design of Blood-Handling Devices, Ph.D. thesis, RWTH Aachen University,2016.[22] J. Donea, A. Huerta, Finite Element Methods for Flow Problems, WILEY, 2003.[23] T. E. Tezduyar, M. Behr, A New Strategy for Finite Element Computations Involving Moving Boundaries and Interfaces - The Deforming-Spatial-Domain / Space-Time Procedure: I. The Concept and the Preliminary Numerical Tests, Computer Methods in Applied Mechanics andEngineering 94 (1992) 353–371.[24] T. E. Tezduyar, M. Behr, A New Strategy for Finite Element Computations Involving Moving Boundaries and Interfaces - The Deforming- patial-Domain / Space-Time Procedure: II. Computation of Free-Surface Flows, Two-Liquid Flows, and Flows with Drifting Cylinders,Computer Methods in Applied Mechanics and Engineering 94 (1992) 353–371.[25] S. Elgeti, H. Sauerland, Deforming Fluid Domains within the Finite Element Method: Five Mesh-Based Tracking Methods in Comparison,Archives of Computational Methods in Engineering 23 (2016) 323–361.[26] M. Behr, F. Abraham, Free-Surface Flow Simulation in the Presence of Inclined Walls, Computer Methods in Applied Mechanics andEngineering 191 (2002) 5467–5483.[27] J. A. Cottrell, T. J. R. Hughes, Y. Bazilevs, Isogeometric Analysis - Toward Integration of CAD and FEA, WILEY, 2009.[28] J. Chung, G. M. Hulbert, A Time Integration Algorithm for Structural Dynamics with Improved Numerical Dissipation: The Generalized- α Method, Journal of Applied Mechanics 60 (1993) 371–375.[29] D. Kuhl, A. Crisfield, Energy-Conserving and Decaying Algorithms in Non-Linear Structural Dynamics, International Journal for NumericalMethods in Engineering 45 (1999) 569–599.[30] T. J. R. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric Analysis: CAD, Finite Elements, NURBS, Exact Geometry and Mesh Refinement,Computer Methods in Applied Mechanics and Engineering 194 (2005) 4135–4195.[31] L. Piegl, W. Tiller, The NURBS Book, Springer Science & Business Media, 2012.[32] E. Ramm, W. A. Wall, Shell Structures - A Sensitive Interrelation between Physics and Numerics, International Journal for NumericalMethods in Engineering 60 (2004) 381–427.[33] W. Dornisch, Interpolation of Rotations and Coupling of Patches in Isogeometric Reissner-Mindlin Shell Analysis, Ph.D. thesis, RWTHAachen University, 2015.[34] W. Dornisch, S. Klinkel, B. Simeon, Isogeometric Reissner-Mindlin Shell Analysis with Exactly Calculated Director Vectors, ComputerMethods in Applied Mechanics and Engineering 253 (2013) 491–504.[35] W. Dornisch, S. Klinkel, Treatment of Reissner-Mindlin Shells with Kinks without the Need for Drilling Rotation Stabilization in anIsogeometric Framework, Computer Methods in Applied Mechanics and Engineering 276 (2014) 35–66.[36] N. Hosters, J. Helmig, A. Stavrev, M. Behr, S. Elgeti, Fluid–Structure Interaction with NURBS-based Coupling, Computer Methods inApplied Mechanics and Engineering 332 (2018) 520–539.[37] P. Causin, J.-F. Gerbeau, F. Nobile, Added-Mass E ff ect in the Design of Partitioned Algorithms for Fluid-Structure Problems, ComputerMethods in Applied Mechanics and Engineering 194 (2005) 4506–4527.[38] E. H. Van Brummelen, Added Mass E ff ects of Compressible and Incompressible Flows in Fluid-Structure Interaction, Journal of AppliedMechanics 76 (2009) 021206.[39] B. W. Uekermann, Partitioned Fluid-Structure Interaction on Massively Parallel Systems, Ph.D. thesis, Technical University of Munich, 2016.[40] A. E. Bogaers, S. Kok, T. Franz, Strongly Coupled Partitioned FSI Using Proper Orthogonal Decomposition, Conference Paper, Eighth SouthAfrican Conference on Computational and Applied Mechanics (SACAM) (2012).[41] G. H. Golub, C. F. Van Loan, Matrix Computations, JHU Press, 2012.[42] K. Scheufele, Coupling Schemes and Inexact Newton for Multi-Physics and Coupled Optimization Problems, Ph.D. thesis, University ofStuttgart, 2018.[43] I. Farmaga, P. Shmigelskyi, P. Spiewak, L. Ciupinski, Evaluation of Computational Complexity of Finite Element Analysis, 11th InternationalConference The Experience of Designing and Application of CAD Systems in Microelectronics (CADSM) (2011) 213–214.[44] B. Graham, A. Adler, A Nodal Jacobian Inverse Solver for Reduced Complexity EIT Reconstructions, International Journal of Informationand Systems Sciences 2 (2006) 453–468.[45] B. Zhou, D. Jiao, A Linear Complexity Direct Finite Element Solver for Large-scale 3-D Electromagnetic Analysis, 2013 IEEE Antennasand Propagation Society International Symposium (APSURSI) (2013) 1684–1685.[46] P. Greisen, M. Runo, P. Guillet, S. Heinzle, A. Smolic, H. Kaeslin, M. Gross, Evaluation and FPGA Implementation of Sparse Linear Solversfor Video Processing Applications, IEEE Transactions on Circuits and Systems for Video Technology 23 (2013) 1402–1407.ects of Compressible and Incompressible Flows in Fluid-Structure Interaction, Journal of AppliedMechanics 76 (2009) 021206.[39] B. W. Uekermann, Partitioned Fluid-Structure Interaction on Massively Parallel Systems, Ph.D. thesis, Technical University of Munich, 2016.[40] A. E. Bogaers, S. Kok, T. Franz, Strongly Coupled Partitioned FSI Using Proper Orthogonal Decomposition, Conference Paper, Eighth SouthAfrican Conference on Computational and Applied Mechanics (SACAM) (2012).[41] G. H. Golub, C. F. Van Loan, Matrix Computations, JHU Press, 2012.[42] K. Scheufele, Coupling Schemes and Inexact Newton for Multi-Physics and Coupled Optimization Problems, Ph.D. thesis, University ofStuttgart, 2018.[43] I. Farmaga, P. Shmigelskyi, P. Spiewak, L. Ciupinski, Evaluation of Computational Complexity of Finite Element Analysis, 11th InternationalConference The Experience of Designing and Application of CAD Systems in Microelectronics (CADSM) (2011) 213–214.[44] B. Graham, A. Adler, A Nodal Jacobian Inverse Solver for Reduced Complexity EIT Reconstructions, International Journal of Informationand Systems Sciences 2 (2006) 453–468.[45] B. Zhou, D. Jiao, A Linear Complexity Direct Finite Element Solver for Large-scale 3-D Electromagnetic Analysis, 2013 IEEE Antennasand Propagation Society International Symposium (APSURSI) (2013) 1684–1685.[46] P. Greisen, M. Runo, P. Guillet, S. Heinzle, A. Smolic, H. Kaeslin, M. Gross, Evaluation and FPGA Implementation of Sparse Linear Solversfor Video Processing Applications, IEEE Transactions on Circuits and Systems for Video Technology 23 (2013) 1402–1407.