PARALLEL IN TIME ALGORITHMS FOR NONLINEAR ITERATIVE METHODS

In this paper we investigate the feasibility of applying the Parareal algorithm [5, 6] for quasi-static nonlinear structural analysis problems. We describe how this proposal has been realized and present some preliminary numerical results of applying this algorithm to a beam undergoing nonlinear deflection with a contact boundary condition. Further numerical experiments are needed to provide an evidence for the efficiency of the method. Résumé. Dans cet article on cherche à determiner si l’algorithme de Parareel [5, 6] est applicable pour des structures quasi-statiques non linéaires. On décrit comment cette approche a été réalisée puis on présente quelques résultats numériques préliminaires de notre algorithme pour un problème poutre non linéaire dont une condition au limite modélise un contact. D’autres expériences numériques sont nécessaires afin de prouver l’efficacité de la méthode.


Introduction
The simulation of complex nonlinear structures via the finite element method is generally based on quasi-static incremental loading procedure leading to successive nonlinear problems solved by Newton-like methods [7,9] which themselves entail solving multiple ill-conditioned large linear systems.For complex structures, the total number of steps may be very large and, since the procedure is recursive, the computational cost is very large.Reducing the time for such simulations is of great interest for numerical engineering in industry.A recently developed parallel in time method for time dependent problems is the Parareal algorithm [5,6] that allows the parallelization of the computation in the temporal domain making use of HPC facilities.
The objective of this work is to present a proposal for the application of the Parareal algorithm for quasistatic nonlinear processes in order to perform these quasi-static steps in parallel by analogy with the Parareal method for time-dependent systems of ordinary or partial differential equations.In this work we describe how this proposal has been realized on a common example of a beam undergoing nonlinear elastic deformation with additional boundary nonlinearity.Numerical results demonstrate the possibility of using the Parareal algorithm on elementary quasi-static problems, although more conclusive numerical tests are needed to provide a substantial evidence for the efficiency of the approach.
Plan of the paper.We present the model problem, its spacial discretization and linearization in Section 1. Basic idea and precise desription of the Parareal algorithm for the case of quasi-static nonlinear problems are laid out in Section 2. Numerical experiments for several test cases are reported in Section 3. We discuss some reoccurring themes and questions that faced us during this work in Section 4. Finally, we draw some concluding remarks in Section 5.

Nonlinear structures and incremental loading
We treat the boundary value problem for the equations of nonlinear elasticity [1,8] as a minimization problem.For a domain Ω = [0, 1] × [0, 1] ⊂ R 2 we want to find the displacement field (u 1 , u 2 ) : Ω → R 2 that minimizes the functional of total potential energy J(u 1 , u 2 ): min where V is a suitable function space that satisfies some Dirichlet boundary conditions: where ū1 and ū2 are some given constants, Γ v denotes the part of the boundary on which the displacements are given.The total potential energy is given by the following formula: where ) is a bilinear symmetric positive form with respect to matrices X and Y , f is a given C 2 function corresponding to hyperelastic constitutive law, P α is some external force which we apply on the boundary surface Γ P ⊂ ∂Ω (for the sake of simplicity we suppose that the force act in the x 2 -direction) and is the Green-Saint-Venant strain tensor.
The finite element method is the standard modeling approach to simulate and analyze the behavior of solids [2].Let us approximate the space V by the finite element space V h .We thus introduce a regular mesh T h on Ω, a triangle K of the mesh T h and a standard P 1 finite element space on it V h ⊂ V : Thus we obtain the finite element problem of finding the field (u h , v h ) that satisfy (2) and min The solution of the large strain nonlinear problem will require an iterative process and Newton's method forms the basis of practical schemes [9].We start with an initial guess (u 0 , v 0 ) which is reasonably close to the true correct solution.Let us suppose that we already know the solution (u n−1 , v n−1 ) from the step n − 1 then we attempt to correct this guess to bring it closer to the proper solution by setting where the correction term (du, dv) is given by where DJ(u, v) and D 2 J(u, v) are the Jacobian and the Hessian of J(u, v) respectively.Then we check if the magnitude of the correction (du, dv) is small enough, and if not, we go to the next iteration n + 1 and compute the new correction term ( du, dv).  1. Convergence of Newton's method depending on the number of loading steps.Cantilever rubber beam loaded at the upper bound with a uniformly distributed vertical force F = 200 P a with a contact surface, mesh 60 × 20 points, N I -average number of Newton iterations, N IC -number of Newton iterations at the loading step when the contact appears for the first time.

Incremental loading
There is no guarantee that Newton's method will converge.It will converge quadratically to the exact solution only if the initial guess is sufficiently close to the correct solution.A common way to avoid this problem is to apply the load in a series of increments instead of all at once [7,9].For example, we would like to solve our problem for some certain boundary force F , but numerical test shows us that Newton's method does not converge.We can divide the force into 2 loads: at the beginning we solve our problem with force F/2 and after we solve the problem once more but with force F using the solution at the end of the preceding increment as the initial guess.The quality of this guess could be improved by reducing the increment once more.In general, small load increments are essential for a better accuracy of the solution.In Table 1 we show the case of a cantilever rubber beam loaded with a uniformly distributed force F = 200 P a and a contact surface.In the case of applying the force all at once, we can't achieve the convergence of Newton's method.The situation is persistent when using 2 and 3 loading steps, whereas we achieve convergence when we start using 4 loading steps and so on.Moreover, this example demonstrates that the number of incremental loads can affect the number of iterations required by Newton's method to achieve the convergence.Indeed, when the force is divided into 7 loading steps we only need 3 Newton iterations to achieve the same tolerance like in the case of 4 loading steps with 4 Newton iterations.This example shows that the incremental loading procedure is required in the case when the initial guess for Newton's method is too far from the exact solution.Note that with incremental loads we obtain essentially a quasi-static nonlinear problem instead of a static one; from now on we have a so-called "pseudo-time" direction which is in fact an additional dimension where it is possible to apply the idea of Parareal method.But before discussing this topic we will introduce a boundary nonlinearity to our problem.

Boundary nonlinearity
The incremental loading procedure allows us to insert a different form of nonlinearity in our model: a contact problem [4].Let us consider the cantilever beam: under a vertical load P α and in one moment the beam touches a solid surface X 2 = X sf with border Γ S where X sf and (X 1 , X 2 ) are the Eulerian coordinates.We model the appearance of the contact surface with the following algorithm: at every incremental loading step we look for the nodes of our mesh at the border of the beam and check if the vertical displacements of those nodes are equal to or greater than the vertical coordinate of the surface X sf (See Figure 1 a), in what follows we call them violated nodes [3].On the next incremental loading step for every violated node we set the sliding contact boundary condition (a vertex of the finite element mesh can slide over the surface, but cannot move away from it), in other words we take the vertical displacement equal to the displacement corresponding to the contact surface with Dirichlet boundary conditions (See Figure 1 b): The number of incremental loadings is very important for detecting the contact surface.From Table 1 we see that for the cantilever beam configuration loaded with an uniformly distributed force we achieve the maximum number of the Newton iterations in the load when the contact appears, because we are converging in presence of such strong nonlinear effect like a contact.Thus the accuracy of the solution depends on the nonlinearity involved in the model and the number of taken load increments.

Parareal method for ODEs
In order to understand the mechanism of the parallelization in time we will briefly explain the basic idea of the parareal method with a simple first order ordinary differential equation [5]: where a is a constant.Let us introduce the temporal mesh on the time interval [0, T ]: with time steps ∆T = T /N .The Parareal algorithm utilizes a coarse solver G to quickly step through time by computing relatively cheap approximate solutions for all time intervals of interest, and then simultaneously refines all of these approximate solutions using an accurate fine solver F .Let G ∆T (u n−1 , T n−1 ) for all n = 1, . . ., N be a rough approximation of u(T n ) with an initial condition u( ) is a more accurate approximation of the solution u(T n ) with an initial condition u(T n−1 ) = u n−1 .Then starting with a coarse approximation U 0 n at the time points T 1 , T 2 , . . .T N , Parareal performs for k = 0, 1, . . . the correction iteration: The application of the fine solver to each time interval [T n , T n+1 ] is independent of the other time intervals, and thus parallelizable.From (8) we see that the refined solutions are then fed back to the coarse solver, and the iterative cycle continues until all time intervals are converged.The dependencies between the time intervals are carried through the coarse solver which involves stepping sequentially through time.This, of course, represents a sequential process, but it leads to a more rapid solution, assuming that the coarse solver is much faster than the fine solver, and hence, computationally negligible.In the simple numerical example for equation (7) we set T = 10 and perform 5 coarse steps each parareal iteration, thus the coarse solution is obtained by the implicit Euler scheme with step ∆T = 2.After that, we introduce a finer mesh with a smaller time step ∆t = 0.02 on each interval [T n , T n−1 ], thus we have 100 fine time steps on each interval [T n , T n−1 ].We use the same implicit Euler scheme as a fine solver.
We present the error in log scale between the parareal solution and the analytical solution of equation ( 7) for each of the 5 parareal iterations in Figure 2. In order to show the convergence of the parareal method at each parareal iteration we also present the error between the solution of Euler's method with 500 time steps and the analytical solution of equation (7).
The solution of the first coarse time step is processed by the fine solver using the exact initial state u(0).For every further parareal iteration, the converged accurate solution is computed with initial conditions from the previous iteration, thus the second iteration will give us convergence at least for the first two steps, the third iteration for the first three steps, and so on.Thus the algorithm is guaranteed to converge in at most K = N iterations (see Figure 2), though in practice it is possible to achieve much faster convergence with judicious choice of the coarse propagator.The computational complexity of the parareal algorithm is given by the following formula: where C C is the cost to perform the coarse propagator over one time step, C F is the complexity of one time step for the fine propagator, K is the number of iterations required to achieve the required convergence of the parareal algorithm and N is the number of coarse steps (and the number of processors as well).
Looking at the Parareal algorithm from a different perspective, it is reasonable to presume that the idea of parallelizing the time domain for time-dependent problems can be used for the case of quasi-static problems in order to reduce the computation time.

Parareal method for nonlinear structures
In what follows we shall propose a "time" parallel iterative method for solving the minimization problem (5).At first, we introduce N M load increment steps and assume that we compute the coarse solution at points 0, M, 2M, . . ., N M and the fine solution at points 0, 1, 2, . . ., N M .Thus we introduce a pseudo-time parameter t ∈ [0, t1 , t2 , . . ., tNM ] which is used to describe our algorithm via an analogy with a time-dependent problem.For the sake of simplicity we assume that all increments are the same and denote points for the coarse solution as [0, T1 , T2 , . . ., TN ], in other words tM = T1 , t2M = T2 , . . ., tMN = TMN .At every point [0, t1 , t2 , . . ., tNM ] we introduce the exact solution of problem (5) (u( tn ), v( tn )) for n ∈ [0, 1, . . ., N M ].We want to find approximate solution where index k is the number of parareal iterations.Like in the case of a partial differential equation, let G(U k nM , V k nM ) for all n = 0, . . ., N − 1 be a rough approximation of (u( Tn+1 ), v( Tn+1 )) which we compute starting from initial guess (U k nM , V k nM ) by using the coarse propagator that is described below.Let F (U k nM , V k nM ) be a more accurate approximation of the solution (u( Tn+1 ), v( Tn+1 )) which we compute starting from displacements (U k nM , V k nM ) by using the fine propagator that is described below.We start from the initial guess (U 0 0 , V 0 0 ) = (u 0 , v 0 ) for Newton's method at the first increment step and obtain an initial coarse solution via That is we obtained Parareal solution of iteration k = 0. Hereafter on every Parareal iteration k ≥ 1 we determine our numerical solution by the formula: The coarse propagator Once the solution is known at iteration k − 1, the coarse problem at iteration k + 1 is Newton's method for N incremental loads.For all coarse increment steps n = 1, . . ., N starting with the initial guess for Newton's method as displacements from the previous parareal iteration coarse solution at the Newton iteration j as: where the correction term (du j n , dv j n ) by analogy with (6) follows from Let's assume that at the Newton iteration number j * we achieved the required accuracy for approximate solution, then the coarse propagator is

The fine propagator
We use a procedure similar to (11) and ( 12) to obtain the fine solution F (U k nM , V k nM ) on parareal iteration k.Inside every "coarse" pseudo-time interval [ Tn , Tn+1 ] we introduce a local fine solver on fine mesh F (U k nM +m , V k nM +m ) to emphasize the fact that this solver propagates the solution on fine increment steps.We compute the solution with the same Newton's method for M incremental loads using the displacements from the previous Parareal iteration as an initial guess for the first incremental load 10) Correction End for End for Table 2. Parareal algorithm for nonlinear quasi-static structure.
For all the next fine increment steps m = 1, . . ., M − 1 we use displacements from the previous increment step as an initial guess for Newton's method, in other words, we start by setting the initial guess as of the Newton iteration j − 1 then the solution of iteration j is: where the correction term (du j m , dv j m ) like in the case of the coarse solver ( 12) is defined by (14) Let's imagine that at the Newton iteration number j * we achieve the required tolerance then we determine the fine propagator as F U k nM +m , V k nM +m = F j * U k nM +m , V k nM +m for all fine increment steps m = 1, . . ., M − 1.Thus we determine the fine propagator at coarse time points as

A numerical study of the Parareal algorithm for a nonlinear beam
We now turn to numerical examples for our algorithm.What we report here are preliminary results that have to be extended to more complex cases.We run our experiments for four different cases with different boundary conditions: a beam anchored at one or two ends, with and without a contact; see Figure 3.However, here we report numerical results only for the case of a one-end fixed beam with a contact, see Figure 3 d 3. Nonlinear beam with one fixed end and a contact.Convergence for fixed mesh 60 × 20 , F = 120 P a, 6 coarse steps, 6 fine steps.The absolute error is calculated between the parareal solution and the direct method in L ∞ norm .
it is an interesting configuration in the sense of the degree of nonlinearity involved.All the results reported hereafter correspond to a vertically uniform distributed external load applied to upper bound of a rubber beam.In the numerical experiments, FreeFem++ is used for the finite element formulation and the linear algebra implementation, the fine solver is parallelized using MPI.
The first experiment considered here is the Parareal algorithm parallelized on 4 cores, 6 coarse steps and 6 fine steps are used, see Table 3.In order to measure the quality of the solution we are looking for L ∞ and L 2 norms of displacements.The first line of Table 3 corresponds to displacements obtained by applying 36 loading steps for Newton's method.We observe that we have a convergence from the first Parareal iteration, which is reasonable since our numerical example is performed with a sufficiently fine mesh.Indeed, the coarse solver converges and it means that 6 loading steps are enough to reach an accurate solution under the final load and since our process is quasi-static we do not care about the solution at every "time" point.We would like to underline that our aim here is to demonstrate the possibility of using Parareal algorithm on elementary quasi-static problems in order to inspire ourselves for further more realistic and strongly nonlinear models.
The problem formulation for the second experiment is the same as that of the previous experiment except that here we have used a different space mesh size for the coarse and fine propagation.Due to the sequential nature of the coarse propagator we would like it to be as fast and cheap as possible.Since the coarse propagator will be corrected with the fine propagator we can use rough space mesh coarse steps.Results with different coarse and fine mesh sizes are summarized in Table 4. First line of Table 4 corresponds to displacements obtained with loading steps for direct Newton's method.Note that as in the previous example we have a convergence from the first Parareal iteration.

Discussion
The results presented in Section 3 support the proposal for the feasibility of applying the parareal method for nonlinear structural analysis problems.That said, there is a crucial point that needs to be clarified; the pseudo time stepping is only there in order to ensure the convergence of Newton's method.It is thus more like a stability problem: once we have a converged result then it is the good one.Thus a legitimate question is why

1 Figure 1 .
Figure 1.Contact surface modulation: we detecte the violated nodes at displacements corresponding to load increment n, then we impose Dirichlet boundary conditions for the violated nodes before the start of load increment n + 1 .

Figure 3 .
Figure 3. Beam configuration: a) beam fixed at two ends without a contact, b) beam fixed at two ends with a contact, c) beam fixed at one end without a contact, d) beam fixed at one end with a contact.