Parallel time-stepping for fluid-structure interactions

We present a parallel time-stepping method for fluid-structure interactions. The interaction between the incompressible Navier-Stokes equations and a hyperelastic solid is formulated in a fully monolithic framework. Discretization in space is based on equal order finite element for all variables and a variant of the Crank-Nicolson scheme is used as second order time integrator. To accelerate the solution of the systems, we analyze a parallel-in time method. For different numerical test cases in 2d and in 3d we present the efficiency of the resulting solution approach. We also discuss some special challenges and limitations that are connected to the special structure of fluid-structure interaction problem.

suited numerical schemes can be used for each of the subproblems. There are however application classes where partitioned approaches either fail or lack efficiency. The added mass effect [12] exactly describes this special stiffness connected to fluid-structure interactions. It is typical for problems with similar densities in the fluid and the solid -as it happens in the interaction of blood and tissue or in the interaction of water and the solid structure of a vessel. Here, monolithic approaches are considered to be favourable.
Monolithic approaches all give rise to strongly coupled, usually very large and nonlinear algebraic systems of equations. Although there has been substantial progress in designing efficient numerical schemes for tackling the nonlinear problems [17,24,27] (usually by Newton's method) and the resulting linear systems [3,21,31,38,48], the computational effort is still immense. Numerically accurate results and efficient approaches for 3d problems are still very rare [16,39].
In this contribution we exploit the perspectives (and limitations) of parallel time-stepping schemes for fluidstructure interaction systems. We will have to face various difficulties that come from the special type of equations and coupling, such as the hyperbolic property of the solid problem and the saddle-point structure of the fluid system. Parallel time-stepping methods are well established in the literature [7,11,13,18,23,30,44], but up to now there is little experience with fluid-structure interactions. However Parareal has shown to be unstable for hyperbolic problems, where [43] found the phase error between the fine and coarse propagators to be responsible, which already points to possible complications for our use case. In Section 4.2 we will analyse a benchmark configuration where exactly such a phase error will prevent convergence of the Parareal scheme. We also find that the moving domain formulation, which is characteristic for fluid-structure interactions, gives rise to further issues when an incompressible fluid is considered. Our discrete formulation will be based on Arbitrary Eulerian Lagrangian coordinates by mapping the moving fluid domain onto a fixed reference domain. Hereby geometric nonlinearities are introduced and when solutions from different discretizations, i.e. fine and coarse propagators, are linearly combined, incompressibility is disturbed.
Parareal has been investigated for the solution of the Navier-Stokes equations [13,18,46]. In [46] it has been shown that high Reynolds numbers lead to poor convergence of the algorithm. For low Reynolds numbers, as will be tested here, we expect that the viscous dissipation is dominant. The dissipation properties of coupled fluid-structure interactions are not fully understand and they depend on the material parameters but also on the geometry. Stability can be increased by adding damping within the structure or by damping the coupling between fluid and structure. The analysis is mostly restricted to linear model problems where the Stokes equation is coupled to linear elasticity, see [4][5][6].
In [13], the potential speedup of parareal beyond the saturation of spatial parallelization for the Navier-Stokes equations is investigated. Since we don't use other means of parallelization a comparison to the present work seems difficult.
In the following section we will shortly describe a monolithic Arbitrary Lagrangian Eulerian formulation for fluid-structure interaction problems. Further, we detail on the discretization of this system in space (with continuous finite elements) and time (with classical time-stepping methods). The following 3rd section will focus on describing a parallel time-stepping scheme and the special requirements for a realization in terms of fluid-structure interactions. Finally, in Section 4 we present numerical results and discuss possible applications but we also address some limitations.

Fluid-structure interactions
The presentation within this section mainly follows [39]. We consider fluid-structure interaction problems coupling an incompressible fluid with a hyperelastic solid. ByŜ we denote the Lagrangian reference framework of the solid, by S(t) its current configuration in Eulerian coordinates. By F(t) we denote the fluid domain at time t matching the solid at the common interface I(t) = ∂S(t) ∩ ∂F(t). By Ω(t) := F(t) ∪ I(t) ∪ S(t) we denote the Eulerian fluid-structure interaction domain. The domains F(t), S(t) and Ω(t) are all either two-dimensional or three-dimensional. The boundary of the fluid domain ∂F(t) = I(t) ∪ Γ D f (t) ∪ Γ out f (t) is split into the interface, a Dirichlet part Γ D f (t) (usually inflow or rigid walls) and an outflow part Γ out f (t), where we ask for the do-nothing condition [25]. For simplicity we assume Dirichlet conditions at the solid boundary apart from the interface Finally, by I = [0, T ] we denote the time interval. With the density ρ f the incompressible Navier-Stokes equations are given by with the right hand side field f f , boundary data v D f (t) and the interface velocity d t u s coming from the coupling to the solid equation.
By v 0 f we denote the initial velocity. The solid problem in term is given in Lagrangian reference formulation onŜ aŝ where we denote byF s = I +∇û s the deformation gradient, byρ s the reference density, byf s the right hand side vector, boundary data byû D s and by σ f (v f , p f ) n f the normal stresses from the coupling to the fluid equations. The hats "· " are added to distinguish Lagrangian variables form their counterparts in the reference framework. The attribution of the kinematic interface condition v f = d t u s to the fluid problem and the dynamic condition to the solid problem is artificial, as the coupled system of (2.1) and (2.2) must be considered as one entity. Finally, as material models we consider a Newtonian fluid and a St. Venant Kirchhoff solid where we denote byÊ s = 1 2 (F T sFs − I) the Green-Lagrange strain tensor and by µ s , λ s the Lamé parameters, by ν f the kinematic viscosity. By (2.1), (2.2) and (2.3) we denote the fluid-structure interaction problem.

Arbitrary Lagrangian Eulerian coordinates
To overcome the discrepancy between the Eulerian fluid framework and the Lagrangian solid framework, we map the flow problem to a fixed reference domainF that fits the Lagrangian solid domainŜ. Here we assume thatF = F(0) is just the known fluid domain at initial time t = 0. ByT f (t) :F → F(t) we denote the reference map, byF f :=∇T f its gradient and byĴ f : we denote the ALE representation of the Eulerian variable v f (same for the pressure). The Arbitrary Lagrangian Eulerian formulation (ALE) goes back to the 70s [14,26,29], and usually consists of adding convective terms with respect to the motion of the domain. We use a strict mapping to the fixed reference system and formulate the coupled variational formulation on this arbitrary frameworkF. Details are given in [39,40]. To close the ALE-formulation we construct the mapT f := id +û f by means of an extension of the solid deformation to the fluid domain, denoted byû f . For small deformations of the fluid domain a simple harmonic extension is sufficient −∆û f = 0 inF,û f =û s onÎ,û f = 0 on ∂F \Î, while problems with large changes in the fluid domain require more care in extending the deformation. Details are given in Section 5.3.5 of [39] and the references therein. We finally give the complete system. From hereon, all problems are given in the reference system such that we skip the hat for better readability Note that n f and n s are the outward facing normals in the reference framework such that the fluid stresses are given in terms of the Piola transform. The fluid stress tensor in ALE formulation reads In (2.4) we have split the hyperbolic solid problem into a system of first order (in time) equations by introducing the solid velocity v s that -on the interface -matches the fluid velocity.

Variational formulation and finite element discretization
To prepare for a discretization with finite elements we briefly sketch the variational formulation of the fluidstructure interaction problem (2.4) that also embeds the coupling conditions in a variational sense. Kinematic and geometric coupling conditions are taken care of by choosing global function spaces for the velocity and deformation, i.e.
where d ∈ {2, 3} is the spatial dimension, Ω = F ∪ I ∪ S and Ω D = Γ D f ∪ Γ D s the combined Dirichlet boundary and u D , v D are extensions of the Dirichlet data into the domain. Solid and fluid velocity (and deformation) are defined as the restrictions of v and u to the respective domain. The dynamic condition is realized by testing the momentum equations for both subproblems by one common and continuous (in the H 1 -sense) test functions φ ∈ H 1 0 (Ω; Ω D ) d and summing up both equations.
By ·, · Γ we denote the L 2 -inner product on a boundary segment Γ. The boundary term is introduced to correct the normal stresses to comply with the do-nothing outflow condition [25]. Dirichlet boundary values are embedded in the trial spaces, the initial conditions will be realized in the time stepping schemes. Next, let Ω h be a triangulation (or more general a mesh) of the domain Ω satisfying the usual conditions of structure-, shape-and form-regularity. In addition we assume that the mesh Ω h matches the interface I in the following sense: for each element K ∈ Ω h it holds K ∩ I = ∅. We refer to Section 5.3.1 of [39] for details on the construction of such meshes. Details on the specific implementation in the finite element library Gascoigne 3D [9] are given in Section 4.2 of [39]. The assumption of matching meshes is important to guarantee good approximation properties, as fluid-structure interactions are interface problems, where the solution has limited regularity across the interface and non-fitted meshes give rise to a breakdown in accuracy [19,42].
The finite element discretization of (2.6) is straightforward. We choose continuous polynomial spaces assembled on the mesh Ω h as where Q (r) is the space of bi-or tri-polynomial of degree r on quadrilateral or hexahedral meshes and where the iso-parametric reference mapT T ∈ [Q (r) ] d comes from this same space. The iso-parametric setup is able to give optimal order approximations on domains with curved boundaries, see [32] or Section 4.2.3 of [39]. Throughout this paper we choose piecewise quadratic discrete trial-spaces for the discretization of velocity, deformation and pressure in (2.6) and and also piecewise quadratic test-spaces for all equations with the necessary modifications to implement Dirichlet values. As the equal order finite element pair V h × Q h does not fulfil the inf-sup condition we add additional stabilization terms of local projection [8] type is the local stabilization parameter. Required small modifications in terms of the ALE formulation in fluid-structure interactions are discussed in [35] or Section 5.3.3 of [39]. For simplicity we assume that no mechanisms for stabilizing dominant transport are required.

Time discretization
For temporal discretization in the time interval I = [0, T ], we start by a splitting into discrete time-steps 0 = t 0 < t 1 < · · · < T N = T, k := t n − t n−1 .
To simplify notation we assume that this distribution is uniform with k being the same in all time steps t n−1 → t n . Modifications are however straightforward. At time t n we denote by v n = v h (t n ), u n = u h (t n ), p n = p h (t n ) the discrete approximations and by F n := I + ∇u n , J n := det F n , the deformation gradient and its determinant. By the index n − 1 2 we denote the mean values on I n = [t n−1 , t n ], e.g. v n := (v n−1 + v n )/2 or J n := (J n−1 + J n )/2.
As a compromise between accuracy (second order), very good stability properties (globally A-stable), simplicity and efficiency (simple one-step scheme) we consider an implicitly shifted version of the Crank-Nicolson method [34,36] and Section 5.1.2 of [39].
Applied to the finite element approximation to (2.6), each time step t n−1 → t n of the fully discrete problem reads: where for simplicity of notation we introduced the notation For θ = 1 2 this is the typical Crank-Nicolson scheme, θ = 1 would give the backward Euler method of first order. We usually consider θ = 1 2 + θ 0 k with a small fixed parameter θ 0 , which results in a second order accurate scheme, see [34,37] and [41] for a numerical study applied to fluid-structure interactions. A proper combination of three substeps with different choices of θ would give the fractional step theta method which is of second order and strongly A-stable, see [49] or [15] for an application to fluid-structure interactions. The discretization of the nonlinear terms including time derivatives (e.g. coming from the domain convection) is still discussed in literature. However, different choices give similar stability and accuracy results, see [41] and Section 5.1.2 of [39].

Parallel time-stepping
The motivation for increased parallelism in numerical algorithms arises from the architecture of modern computers with an ever-increasing number of processing units. Most numerical solvers for differential equations are formulated as sequential algorithms, which denies exploitation of parallelism.
The Parareal algorithm, which is probably the best known parallel time stepping method, promises to bypass this problem by basically providing a wrapper around common sequential algorithms. Although it is known in its fundamentals since 2001, it hasn't been widely applied to fluid structure interactions.
The Parareal-algorithm can be derived in multiple ways, see [20] for an overview. The first approach introduced in [33], is based on a predictor corrector scheme. First the time domain I is divided into L subintervals of equal length. Predictions with coarse timesteps h C are used to provide initial values for parallel computations with fine timesteps h C h F . These results are then used together with old coarse predictions for correcting the solution, providing new initial values. Timesteps are chosen such that T L = m h F = M h C . The algorithm is defined by using two different propagators for the solution over one subinterval. While C is the propagator using the coarse time step size, F is the fine counterpart.
With these definitions and the remarks from the last chapter in mind, the Parareal-algorithm for fluidstructure-interactions is defined by a simple recursive formula where l ∈ {0, . . . , L − 1} is the index of the subinterval and i ∈ N is the iteration count. While the predictor part is sequential (particularly the first iteration for initialization), the fine propagations on each subinterval can be parallelized. Due to the sequential part the parallel speedup will never be optimal, but by reducing the computational effort for the predictor the sequential part can be minimized. However to achieve fast convergence of the algorithm, the predictors accuracy shouldn't differ too much from the fine propagators accuracy. This criterion is hard to quantify since there is no way to predict the optimal corrector for a given fine propagator. In Figure 1 a shared memory implementation with OpenMP is sketched, following the distributed task scheduling from [2]. This results in a (theoretically) nearly optimal runtime since the sequential part of the algorithm is minimized. The theoretically achievable speedup with r where h F and h C are the timestep lengths of the fine and coarse propagator, respectively. The algorithm requires very little communication, so there is little overhead and the estimate holds well in practice. The main bottleneck are the sequential computations which can be minimized by optimal task scheduling and reducing the ratio r. In (3.2) the first term in the denominator corresponds to the initialization phase with the corrector and K N r corresponds to the sequential part of each further iteration.
From (3.2), we conclude that r and K should be minimized. In the impossible scenario that r = 0 and K = 1 we would get optimal speedup. Upon reducing r, the algorithm will require more iterations to converge, reducing K requires a more accurate predictor, resulting in a higher r. Since r and K interact this way, optimizing both at the same time is the main difficulty in applying the parareal-method effectively. As we will see in our examples, we cannot rely on the fact that r is actually the ratio between the computational costs of the coarse and fine propagator, but it still provides a useful guidance.
To our knowledge, there is no previous application of the Parareal idea to fluid-structure interactions. But from the experience with the Navier-Stokes equations we must expect to meet several problems: The nonlinearity of the will have an impact on the ratio between the effort on fine and coarse problem f , since larger time-steps will make the effect of the nonlinearity more severe. Considering higher Reynolds numbers a deterioration of the Parareal convergence rate is reported [46]. This is attributed to dominant imaginary Eigenvalues. The spectrum of the coupled fluid-structure interaction problem including the hyperbolic structure problem is even more diverse and less understood. Nevertheless, the numerical results presented below will show speedups up to about 4 (for a 2d example using 20 cores), which is comparable to a study on the Navier-Stokes equations [47] (showing a speedup up to 4.1 on 15 cores).
Recently a weighted version of the Parareal scheme, the θ-Parareal scheme has been developed [1], introducing weights for the contributions of the coarse propagators. To avoid confusion with the unrelated θ-timestepping scheme, we name the weights Λ:  The idea here was to replace the predictor term by yielding the formula above. One idea for obtaining Λ l i is to minimize the discrepancy of fine and coarse solution: In which case Λ l i can be computed as F, C C, C assuming the same arguments as above. We will use F, C C, C F, F , which is not a least squares solution, but can be interpreted as angle penalization. The factors are computed for each component separately and averaged to obtain a single scaling factor. This way we make sure that Λ l i ≤ 1 by construction. Later on we will shortly investigate the different approaches. In Section 4 the speedup and the convergence of Parareal will be investigated.

Numerical examples
In this section we will apply the Parareal algorithm to two FSI problems. Further we discuss various issues that we were facing when applying the Parareal algorithm to different configurations. Based on these examples we will measure the error introduced by the Parareal algorithm as well as the speedup and efficiency. Measurements are taken on a two socket machine with two Intel Xeon E5-2640 v4 each having 10 cores. The number of subintervals in the Parareal algorithm are chosen as the total number of cores, 20. Other discretization parameters are depicted in Table 2 Figures 2 and 3 show the geometry of the test cases in the 2d and 3d configuration. Both test cases assemble the flow around a wall-mounted elastic obstacle that will undergo a deformation. The problem is driven by a parabolic (bi-parabolic in 3d) inflow profile on Γ in which is oscillating in time  Table 1. The configurations yield the maximum Reynolds number Re 2d = 30 in 2d and Re 3d = 20 in 3d, where L 2d = 0.5 and L 3d = 0.2, the height of the obstacles is chosen as characteristic length scale. On the wall boundary Γ wall we prescribe homogeneous Dirichlet conditions v = 0 and on the outflow boundary Γ out the do-nothing condition, which, in ALE formulation, looks like     , where we subtract the transposed part of the gradient such that (4.1) remains. In 3d, we split the domain and introduce a symmetry boundary, where we prescribe a free-slip, no-penetration condition v · n = 0 and σ n · t = 0, where t are the tangential vectors.   Table 3. Parareal error (2d example). We show the maximum Parareal error over all 20 subintervals. In addition we indicate the discretization error using the fine step size h F = 0.001 on a mesh with 5 360 unknowns and mark the first step, where the maximum Parareal error is below the maximum discretization error. The problem comprises a large deformation such that the nonlinearity of the ALE mapping and the structure model plays an important part in the complexity of this test case. In Figure 5 we show the convergence behavior of the Parareal method for four choices of the large time step h C . For each iteration k of the Parareal method we plot the relative error between the Parareal approximation and a sequential solution (using the same small step size h F = 0.001). The time interval is split into 20 subintervals and we clearly recognize that after the k-th iteration, the solution on the k-th interval does not get any better. All plots also include a dashed line (this is exactly the same line in every plot) which indicates the discretization error, e.g. the error between the sequential solution with stepsize k and a refined reference solution. If the Parareal error is below this line, the Parareal approximation can be considered sufficiently accurate, e.g. considering h C = 0.01, two iterations k = 2 already yield a Parareal approximation error that is smaller than the discretization error. Table 3 shows the maximum (taken over all 20 subintervals) Parareal approximation for the first five Parareal iterations. We also indicate the discretization error for a sequential iteration using the time step size h F , which could serve as stopping criteria for the Parareal iteration. Parareal approximation errors far below the discretization error cannot reduce the overall error such that the simulation should optimally be stopped, once the Parareal error is below this threshold.

2d example
To test the robustness of the Parareal algorithm with respect to the discretization we indicate further results in Table 4. These show that the convergence of the Parareal method is not affected by the choice of the discretization, neither the fine time step h F nor the spatial mesh resolution has a significant impact on the convergence rate.
Furthermore, in Figure 6 we give an overview of the speedup that is obtained in comparison to the sequential simulation. In this example the optimal choice would be h C = 0.05, stopping after 2 iterations such that the discretization error is dominant and resulting in a speedup of 4.44. That is below the theoretical value of 5.78 estimated by from (3.2). The sub optimal performance is mostly due to the strong nonlinearity of the fluid-structure interaction problem. The nonlinear problems are approximated with a Newton's method that benefits from small time-steps as they limit the role of the ALE map. Using coarser time steps for the prediction (that should theoretically improve the speedup) we require more Newton iterations such that the efficiency is reduced. The theoretical predictions can only be reached when the computational costs of each time step are exactly the same, which we cannot guarantee for a nonlinear problem. Table 4. Parareal error (2d example) for fixed coarse time step h C = 0.05 and variations of the fine time step h F on two different spatial discretizations. We show the maximum Parareal error over all 20 subintervals. The tests show that the Parareal convergence is basically independent from the spatial discretization.

fsi-3 benchmark problem. A case where the algorithm fails
In the first experiments, we considered the fsi-3 benchmark problem as introduced by Hron and Turek [28]. In this case Parareal did not converge. Instead we observed pressure oscillations as shown in Figure 7. One possible explanation for these instabilities that are also known from the simulation of incompressible flow problems on moving meshes [10] is the violation of the weak divergence freeness after the Parareal update (u, v) final m+k = (u, v) C m+k−1 + (u, v) F m+k − −(u, v) C m+k , as the nonlinear ALE formulation (JF −1 : ∇v T , ξ) F with separate updates in velocity and deformation does not conserve weak incompressibility and gives rise to pressure oscillations Different approaches to cure this problem, e.g. applying a Stokes projection to the initial solutions on the subintervals in every time step, see [10], did not sufficiently reduce these oscillations. Another and more likely reason for the cause of these instabilities is the dependency of the damping property of the shifted Crank-Nicolson time-stepping scheme on both the spatial and temporal discretization parameters. Considering the Navier-Stokes flow around a rigid body, the frequency of the vortex shedding depends on the time step size. One observes higher frequencies for smaller time step sizes. In Figure 7 of [22] we have documented a similar dependency of the frequency on the spatial discretization. The Parareal update consists of combining solutions from two different time meshes which, due to the different damping, carry an offset. This leads to a systematic error at each macro step that cannot be reduced in further iterations. It would be a remedy to use a time stepping scheme whose damping properties do not depend on the time step size, such as the Crank-Nicolson scheme for θ = 1 2 . This however does not have sufficient stability for fluid-structure interactions, see [41,Fig. 4]. We repeated the fsi-3 test case with the BDF2 time stepping scheme, which however showed the same oscillations. In [43], Ruprecht suggests modifying the update in the parareal algorithm, such that information about the dissipation properties is taken into account.
The fsi-3 benchmark problem is driven by the inner dynamics of the system and the oscillatory motion of the solid is induced by the von Kármán vortex street arising from the flow around the obstacle. In this respect, the example discussed in Section 4.1 is simpler, since it is driven by an oscillating right hand side. For such problems we could not construct a convergent parareal method that simultaneously provides parallel acceleration. We think that the problems stem from a combination of different aspects of the interaction between the incompressible Navier-Stokes equations with an elastic structure and its discretization, and are not easily remedied.

3d example
Next we consider the 3d configuration shown in Figure 3. This test case corresponds to the 2d case discussed in Section 4.1 and since the dynamics is driven by the oscillatory inflow profile, it should be better suited for the Parareal algorithm than the previously discussed fsi-3 problem. Figure 8 shows the solution of the 3d problem t 17 = 6.8, t 18 = 7.2, t 19 = 7.6, t 20 = 8, analogous to the 2d case. The test case leads to a substantial displacement of the solid. Further, the Reynolds number is sufficiently high such that a complex vortex structure appears. The computational effort for the solution of such 3d fluidstructure interaction problems is very high, which is due to the size of the coupled problem (7 coupled unknowns describing pressure, velocity field and deformation field), but also due to the stiffness introduced by the coupling of two problems of different type (parabolic fluid problem vs. hyperbolic solid problem). The test case in this paper has been kept small with less than 100 000 unknowns, compare Table 2, such that a direct solver could be used for all linear problems. If finer meshes are considered, more elaborate approaches are required. Parallel solvers of Newton-Krylow type with multigrid preconditioning have shown to be most effective [3,16,38].
In Figure 9 we show the relative errors of the Parareal approximation in comparison to a sequential simulation for each iteration. Again, a dashed line represents the discretization error of the problem by comparison to a refined reference solution. We observe a similar behavior like in the 2d case: The Parareal error is quickly below Figure 8. Solutions on the last 4 subintervals. The deformation is magnified by a factor of 10. Table 5. Parareal error (3d example). We show the maximum Parareal error over all 20 subintervals. In addition we indicate the discretization error using the fine step size h F and mark the first step, where the maximum Parareal error is below the maximum discretization error. the level of the discretization error and initial intervals do not benefit from further steps. Table 5 collects these findings and here we indicate the maximum Parareal approximation error over all 20 subintervals. In addition we indicate the level of the discretization error. Figure 10 showing the speedups indicates a substantially reduced performance of the Parareal scheme when applied to this 3d test case. The main cause for this deterioration is the linearization of the problem: Newton's method takes longer to converge for big timesteps and the corrections in the Parareal algorithm only seem to aggravate this problem. Although the instabilities noted in Section 4.2 could not be observed here, they start to rise: Larger timesteps influence the dynamics and if they get too large, the problems are not feasible in a stable manner anymore. By choosing the timesteps carefully a speedup of 2.18 can still be obtained.

Conclusion
We applied the Parareal algorithm to fluid-structure interaction problems in 2d and 3d. Two cases were analyzed, in which speedup and convergence in very few iterations were achieved. In these two configurations the dynamics of the coupled problem was driven by oscillatory boundary data. Despite the hyperbolic character of the solid problem that prevents the use of too coarse timesteps, parameter sets are found where the Parareal Figure 9. Velocity error from the 3d case. We see the error between sequential (computed with the fine stepsize h F ) and Parareal solution with different coarse step sizes. They are calculated at the subinterval boundaries as the relative error in the Frobenius norm vp−vs F vs F . Where v s is the sequential solution and v p the Parareal solution. The bright grey line is the difference between the sequential solution and a reference solution. method is both convergent and efficient. In the future it might be beneficial to combine the Parareal approach with multirate time stepping, where different step sizes are chosen for fluid and solid [45].
However, we also saw notable limitations of the algorithm when applied to configurations, where oscillations arise from the internal dynamics of the flow around an obstacle. The combination of different step sizes for coarse and fine discretization can result in a different dynamics and an offset between predictor and corrector. The fsi-3 benchmark problem by Hron and Turek [28] showed to be too challenging and it was not possible to pick discretization parameters that gave speedup and a stable solution at the same time. Another issue is the severe nonlinearity of fluid-structure interactions. The design on nonlinear (and linear) solvers is already a challenge and coarse timesteps lead to an increased effort that limits the possible speedup.
The Parareal algorithm needs careful investigation of the problem before it can be successfully applied. Choosing the right time step sizes is crucial and even then it could still fail. If the problem is suitable on the other hand, it offers a simple way to speedup the solution.