Project topics and materials
Get Complete Project Material File(s) Now! »
CHAPTER 2. NEWTON-KRYLOV METHOD
Introduction
In industrial applications, in order to avoid the inversion of very large matri-ces, time periodic states are often computed as the asymptotic limit solution of an initial value boundary value problem with arbitrarily chosen initial data. This kind of problems can be faced, for instance, in the cardiac contractions modeling. Another example concerns the steady rolling of a viscoelastic tyre [Le Tallec and Rahier 1994, Govindjee et al. 2014] with a periodic sculpture. In this case, the steady state satis es a « rolling » periodicity condition (cyclic steady state [Govin-djee et al. 2014]), including shifts both in time and space: the state u(t; X) at any point X is the same that at the corresponding point observed at the next sculpture one time period T ago u(t; X) = u(t T; R!T X):
Above, R!T denotes the rotation of angle !T , and ! is the rotation speed. So, the angle !T is the size of one sculpture.
In this chapter, we consider the Newton-Krylov method [Chan and Jackson 1984, Telichevesky et al. 1995, Knoll and Keyes 2004], considered as a shooting method for calculating the space-time periodic solution (steady cyclic state) of evolution problems [Govindjee et al. 2014, Brandstetter and Govindjee 2017]. The main idea consists in considering the periodicity condition as an equation for the unknown initial state, which would lead to the periodic solution. Indeed, once we have found the true initial state, the periodic solution is obtained solving the associated initial value problem during only one time period. The equation on the unknown initial state includes an evolution function, requiring to solve the non-linear initial value problem during one time period. The periodicity condition, comparing the resulting state after one time period to the initial one, provides a non-linear equation to solve. So, it is natural to solve this non-linear equation by a classical Newton-Raphson technique. The most challenging part is to construct the global Jacobian. We propose a way to construct it using the tangent matrices of the evolution problem, calculated and stored during each time step of the forward time integration of the problem. Then, the global Jacobian matrix is never assembled in practice, but de ned implicitly by its action during multiplication operation. With that in hand, one of the matrix-free Krylov methods [Saad 2003, Balay et al. 2016] (we use GMRES [Saad and Schultz 1986]) can be applied to solve the linear system, which appears at each global Newton step. This leads to the Jacobian-free Newton-Krylov method [Knoll and Keyes 2004]. Application of this technique, as a shooting method searching for the initial state of the space-time periodic solution (cyclic steady state), to a steady rolling threaded wheel is already proposed and discussed in [Govindjee et al. 2014]. In our work, we propose another point of view on this method, namely, to consider it as a correction after each time period of a standard rolling solution. This modi ed algorithm has been implemented in Michelin industrial code, applied to a full 3D tyre model and compared to the standard asymptotic convergence.
The organization of the chapter is as follows. The classical Newton-Krylov algorithm is described in Section 2.3. It looks for the « correct » initial state, which provides the periodic solution. On each global Newton iteration, given an initial state (starting with an arbitrary guess), it solves the non-linear initial value prob-lem during the rst period, computing the periodicity error, which is used as right hand side of the Newton linear system to correct the initial state for the next step. In this formulation, Newton algorithm is reminiscent to shooting methods.
Section 2.6 is dedicated to an alternative point of view of the same Newton algorithm, where a change of variable leads to the convenient way to implement the method. The working variable will now be the initial state of the current time period instead of the initial state of the rst period. In this way, the method can be easily derived from the classic time evolution by adding the correction step in the end of each period. Therefore, it can be compared with the methods, which look for the periodic state as an asymptotic limit of the solution, in particular those accelerating convergence to the limit state, like Delayed Feedback Control method (see the next chapter). Indeed, it may be considered as an observer-controller process, where correction, performed in the end of each period, is based on the observation during this period.
Sections 2.8 and 2.9 are devoted to implementation of the proposed technique, applied to the model problems, described in Chapter 1. First, the method has been tested on the linear 2D heat problem from Section 1.2, where a planar disk is heated with a source periodically moving along a circular path. As expected, the method only requires one Newton iteration to converge, but does need several Krylov iterations to solve the linear system. Finally, the method has been applied to our main model problem, that is the rolling 3D tyre in presence of stick-slip frictional contact with the soil (Section 1.1.1). This problem is three dimensional and non-linear. The Newton-Krylov method has been implemented along with the non-modi ed time evolution, in order to compare the rates of convergence to the limit space-time periodic solution. We study rst a simpli ed model using MatLab R2016b. Then, we present results of implementation of the Newton-Krylov shooting method into the Michelin industrial code.
Abstract problem formulation
Time continuous problem
Let H be a Hilbert space and let us consider a nonlinear evolution problem (2.1)
on u : R ! H with the space-time periodicity condition (2.2):
F (t; u(t); u(t)) = 0 8t 2 R; (2.1)
u(t + T ) = Su(t); (2.2)
where u = @tu denotes the time derivative of the state variable u, T > 0 is the time period, and the space shift S is a unitary linear operator on H (isometry).
A typical example of such problem is the rolling tyre model presented in Sec-tion 1.1.1. In this case, the function F is given by (1.11), H = f v 2 H1( )3 j vj 0 =
0g, and the space shift operator S is de ned through rotation of the coordinate:
8u 2 H S : u(X) 7!u(R!T X)
with R!T rotation of angle !T .
It is obvious in fact to see that all examples of Chapter 1 enter in the present framework (2.1)-(2.2), the simplest situation being the cardiac model where the space shift is the identity.
Back to the abstract formulation (2.1)-(2.2), let us denote by ’(t; t ; u ) solu-tion at time t (integral line or trajectory) of the di erential equation (2.1), passing through the point (t ; u ). Then, we can de ne the monodromy operator:
t : u 7!’(t + T ; t; u)
which represents the solution after one period of the initial value problem starting from the point (t; u). Under this notation, our problem written with the space-time periodicity condition (2.2) takes the form t(u(t)) = Su(t):
Thus, the space-time periodic solution of (2.1)-(2.2) is given by u(t) = ’(t; 0; u0);
where the initial data u0 is de ned as the solution to the space-time periodicity equation 0(u0) Su0 = 0: (2.3)
Indeed, if we know the « true » initial data (i.e. corresponding to the periodic solution), then the solution of (2.1)-(2.2) is obtained for the rst period by solving
F (t; u; u) = 0; t 2 [0; T ]
u(0) = u0 (2.4) and for any further period it is given by (2.2), provided that we have the space-time shift invariance property 8t 2 R; 8u; w 2 H : F (t + T; Su; Sw) = 0; when F (t; u; w) = 0: (2.5)
This space-time shift invariance has been proved to hold in Lemmas 1-3 for all systems studied in Chapter 1. We will take this assumption (2.5) as a general assumption to be satis ed by the abstract problem (2.1) under study.
So, to nd the solution of the problem (2.1)-(2.2) we have to solve the equa-tion (2.3) on the unknown u0, which is simply the value of u(t) at the point t = 0. As for the existence of a solution of (2.3), it has to be checked on a case by case basis. The theoretical framework of Chapter 3 introduces a general class of linear problems for which existence and stability of periodic solutions can be proved. For nonlinear problems, existence and uniqueness results can be proved when S 1 0 is contracting by a simple application of the xed point theorem.
Time discrete problem
Time discretization of evolution problem (2.1) consists in introducing a se-quence of calculation times t0 < ti < tm = t0 + T with time steps ti = ti+1 ti, and in replacing the implicit time evolution (2.1) by the induction [Quarteroni et al. 2010] ui+1 = ui + ti t(ti; ui; ti; F );
where the increment t is de ned by the time scheme. For an implicit Euler scheme, the increment t is solution of the non-linear equation F (ti + ti; ui + ti t; t) = 0: (2.6)
After such a discretization, the non-linear evolution problem (2.1) with space-time periodicity condition (2.2) becomes ui+1 = ui + ti t(ti; ui; ti; F );i = 0 : : : m 1 (2.7) um = Su0:
If we introduce the discrete monodromy operator operator 0 t : u0 7!um obtained by solving (2.7) from i = 0 to i = m 1,
the time discrete equivalent of (2.3) becomes 0 t(u0) Su0 = 0: (2.8)
In what follows, we will consider the present time discrete problem to be dis-cretized in space too, with solution living in the nite space Rd. In this case, the matrix, representing the space shift S, depends on the mesh and does not change until the mesh is modi ed. For a general mesh, the space shift is de ned through interpolation. It is strictly sparse and contains no more than Nnodes non-zeros by line, where Nnodes denotes the maximal number of nodes by element. If the mesh is periodic with respect to the sculpture size, the space shift is simply a permutation matrix.
Newton shooting technique
We are looking for the initial data u0 solution of the non-linear equation (2.3). This is a nonlinear problem for the unknown u0. We apply the Newton-Raphson method to solve it, which leads to the iterative correction process: un0+1 = un0 + un0+1
with the correction un0 being solution of the linear problem
(D 0(un0) S) (un0+1 un0) = 0(un0) Sun0
where D 0 denotes the gradient of 0 with respect to the initial data u0. On the right part we have the periodicity error (denoted by « nper) calculated during the n-th iteration. Thus the correction un0+1 is de ned by solving
(S D 0(u0n)) u0n+1 = « pern; (2.9)
which leads to the iteration process: u0n+1 = u0n + (S D 0(u0n)) 1 0(u0n) Su0n : (2.10)
The key question here is how to compute D 0, the gradient of the monodromy operator (which we call the tangent monodromy matrix). An approach is proposed and discussed in the following section.
Note that, during each Newton iteration n, computing the periodicity error « nper requires to perform all the evolution from t = 0 to t = T , i.e. to solve the initial value problem (2.4) with initial value un0, in order to obtain 0(un0). So, the residual computation for the « global » Newton loop, solving (2.3) or its time approximation (2.8), includes solving the nonlinear evolution problem (2.1), or respectively (2.7), with « local » Newton loops at each time step during one period.
This technique can be also considered as a kind of shooting method. Indeed, starting with an arbitrary left boundary value, we go through one period evolution to obtain the right boundary value and use it for correction of the left boundary value and so on.
Table of contents :
Introduction
Motivation
Issues and state of the art
Outline of the work
Chapter 1 Reference problems
1.1 Rolling tyre
1.1.1 Mechanical model
1.1.2 Space-time periodic solution (cyclic steady state)
1.1.3 Rotational invariance of the problem
1.1.4 Simplied academic model
1.2 Linear 2D heat problem
1.3 Cardiac model
1.3.1 General bio-mechanical model
1.3.2 Linearized mechanical model
Chapter 2 Newton-Krylov method
2.1 Introduction
2.2 Abstract problem formulation
2.2.1 Time continuous problem
2.2.2 Time discrete problem
2.3 Newton shooting technique
2.4 Approximation of the Jacobian
2.5 Matrix-free Krylov technique: GMRES method
2.6 Alternative realization
2.7 Detailed Newton-Krylov algorithm
2.8 Application. Linear 2D heat problem
2.9 Application. Rolling tyre
2.9.1 Simplied academic model
2.9.2 Industrial model
2.10 Conclusion
Chapter 3 Delayed feedback control
3.1 Introduction
3.2 Linear problem
3.2.1 Abstract problem formulation
3.2.2 Controlled problem. Explicit solution
3.2.3 Optimal gain operator
3.2.4 Eciency analysis
3.2.5 Quasi-optimal practical approach
3.2.6 Time-discrete feedback
3.2.7 Predictor-corrector form
3.2.8 Gain coecients
3.3 Non-linear problem
3.3.1 Predictor-corrector form
3.3.2 Adaptive tuning of the control parameter
3.4 Application. Linear 2D heat problem
3.5 Application. Rolling tyre
3.5.1 Simplied academic model
3.5.2 Industrial model
3.6 Conclusion
Chapter 4 DFC for dynamic systems
4.1 Delayed feedback control for linear dynamics
4.2 Practical approach: control applied to principal modes
4.3 Application. Linear 3D beating heart model
General conclusion and perspectives
Appendices
A Introduction to Lambert W function
B Maximum of the acceleration term
Bibliography
