Parallel-in-time integrators for Hamiltonian systems Claude Le Bris - PowerPoint PPT Presentation

Parallel-in-time integrators for Hamiltonian systems Claude Le Bris ENPC and INRIA Visiting Professor, The University of Chicago joint work with X. Dai (Paris 6 and Chinese Academy of Sciences), F. Legoll (ENPC and INRIA) and Y. Maday (Paris 6 and Brown) SSC Seminar, October 2012 The University of Chicago – p. 1

Parallel in time algorithm for ODEs x ∈ R d x = f ( x ) , ˙ Most algorithms are sequential in nature: x n +1 = Φ ∆ T ( x n ) . The parareal algorithm (Lions, Maday and Turinici, 2001) is the first efficient version of a time integrator that is parallel in time. It is based upon two integrators to propagate the system over a time ∆ T : a fine integrator F ∆ T , to advance of ∆ T by many tiny timesteps; a coarse integrator G ∆ T , to advance of ∆ T by a few large timesteps (hence much cheaper). F ∆ T = (Φ δt ) ∆ T/δt G ∆ T = (Φ dT ) ∆ T/dT and with δt ≪ dT Φ δt ≡ a single step of any one-step integrator. The University of Chicago – p. 2

The parareal paradigm - 1 The parareal algorithm combines the two integrators as follows: � � x k =0 first, coarse propagation that yields n : n x k =0 n +1 = G ∆ T ( x k =0 ∀ n, ) n Then, iterate over k ≥ 0 : compute jumps (in parallel): J k n = F ∆ T ( x k n ) − G ∆ T ( x k n ) � � x k +1 sequential update to obtain n : n x k +1 n +1 = G ∆ T ( x k +1 ) + J k ∀ n, n n The fine solver is called only in the parallel part of the algorithm. The University of Chicago – p. 3

The parareal paradigm - 2 x k +1 n +1 = G ∆ T ( x k +1 ) + F ∆ T ( x k n ) − G ∆ T ( x k n ) n If F ≡ G , then x k +1 n +1 = F ∆ T ( x k +1 ) : similar to a fine sequential integration. n It turns out that, at iteration k , k = ( F ∆ T ) ( k ) ( x 0 ) , x k as if the fine scheme was used sequentially all the way from x 0 up to time k ∆ T . In practice, convergence occurs much faster: only a few iterations are needed to obtain a good accuracy on a large range. The approach has been successfully tested, in the past few years, on many equations (ODEs or PDEs, e.g. heat equation). The University of Chicago – p. 4

A simple test-case: ˙ x = − x (coarse and fine propagators: forward Euler) 1 fine integrator (sequential) parareal iteration 0 (coarse integrator) parareal iteration 3 parareal iteration 6 1e-05 parareal iteration 10 1e-10 1e-15 1e-20 1e-25 1e-30 1e-35 0 10 20 30 40 50 60 t Trajectory x k ( t ) ( ∆ T = 0 . 5 , dT = ∆ T , δt = 0 . 005 ) At iteration k = 10 , trajectory is exact up to time 10∆ T = 5 , and is accurate over a much longer range. The University of Chicago – p. 5

Hamiltonian dynamics Consider a Hamiltonian function H ( q, p ) , and the Hamiltonian dynamics q = ∂H p = − ∂H q ∈ R d , p ∈ R d . ˙ ∂p ( q, p ) , ˙ ∂q ( q, p ) , Introducing x = ( q, p ) ∈ R 2 d , recast the dynamics as   0 Id   . x = f ( x ) = J ∇ H ( x ) , ˙ J = − Id 0 Separated case: H ( q, p ) = 1 2 p T p + V ( q ) . Then q = p, ˙ p = −∇ V ( q ) . ˙ Verlet algorithm (explicit, second order, symmetric, symplectic): p n − δt = 2 ∇ V ( q n ) , p n +1 / 2 = q n + δt p n +1 / 2 , q n +1 p n +1 / 2 − δt = 2 ∇ V ( q n +1 ) . p n +1 The University of Chicago – p. 6

Reversibility x = f ( x ) = J ∇ H ( x ) ˙ Let ϕ t ( x ) be the flow of the equation, namely the solution at time t with I.C. x : dϕ t ( x ) = J ∇ H ( ϕ t ( x )) , ϕ t =0 ( x ) = x. dt we will focus on the case H ( q, p ) = 1 2 p T p + V ( q ) , for which the flow is reversible: ρ ◦ ϕ t = ( ϕ t ) − 1 ◦ ρ ∀ t, where ρ ( x ) = ρ ( q, p ) = ( q, − p ) . The University of Chicago – p. 7

Reversibility x = f ( x ) = J ∇ H ( x ) ˙ Let ϕ t ( x ) be the flow of the equation, namely the solution at time t with I.C. x : dϕ t ( x ) = J ∇ H ( ϕ t ( x )) , ϕ t =0 ( x ) = x. dt we will focus on the case H ( q, p ) = 1 2 p T p + V ( q ) , for which the flow is reversible: ρ ◦ ϕ t = ( ϕ t ) − 1 ◦ ρ ∀ t, where ρ ( x ) = ρ ( q, p ) = ( q, − p ) . a subset of such systems is the case of integrable reversible systems: up to a reversible change of variables ( q, p ) �→ ( a, θ ) , the dynamics reads ˙ a ∈ R d , θ ∈ T d . a = 0 , ˙ θ = ω ( a ) , Thus a ( q, p ) ∈ R d is left invariant by the dynamics. The University of Chicago – p. 7

Backward error analysis for ODEs ϕ δt ( x ) exact flow x = f ( x ) ˙ numerical scheme Φ δt ( x ) exact flow x = f δt ( x ) ˙ Exact flow of the modified equation at time δt ≡ Flow Φ δt ( x ) given by the numerical scheme after one-time step. The University of Chicago – p. 8

Classical analysis for Hamiltonian systems Symmetric integrators on reversible systems: ∀ δt, Φ δt ◦ Φ − δt = Id Symmetric integrator ⇒ Reversible modified equation Reversible modified equation + integrable reversible system + non resonant condition ⇒ Conservation of invariants (hence of energy). Another class of efficient schemes: Symplectic integrators The University of Chicago – p. 9

Parareal scheme for Hamiltonian dynamics x = ( q, p ) ∈ R 2 d x = f ( x ) = J ∇ H ( x ) , ˙ Standard parareal algorithm: x k +1 n +1 = G ∆ T ( x k +1 ) + F ∆ T ( x k n ) − G ∆ T ( x k n ) n F ∆ T = (Φ δt ) ∆ T/δt G ∆ T = (Φ dT ) ∆ T/dT with and with δt ≪ dT, where Φ δt ≡ one step of (say) Verlet algorithm with δt time-step. Symplectic or symmetric algorithms are known to be efficient for Hamiltonian dynamics. PROBLEM: Even if F ∆ T and G ∆ T are symplectic (or symmetric), the resulting parareal algorithm is NOT symplectic (or symmetric) . . . Long-term (and geometric) properties of the approach are to be clarified! The University of Chicago – p. 10

Harmonic oscillator case: q = p , ˙ p = − q ˙ 1e+06 Sequential coarse scheme Parareal it. #1 Parareal it. #4 10000 Parareal it. #10 Sequential fine scheme 100 1 0.01 1e-04 1e-06 1e-08 1e-10 1 10 100 1000 t Error on the energy 1 2 q 2 + 1 2 p 2 The University of Chicago – p. 11

Outline Adapt the parareal algorithm to the Hamiltonian framework by using symmetrization projection on the constant energy manifold � � ( q, p ) ∈ R 2 d ; H ( q, p ) = H 0 M = The University of Chicago – p. 12

Symmetric parareal algorithm Symmetrizing by hand a given integrator is easy: consider a numerical flow Φ ∆ T , that reads x n +1 = Φ ∆ T ( x n ) ∆ T = Φ − 1 consider the adjoint of Φ ∆ T , defined by Φ ⋆ − ∆ T then the scheme x n +1 = Φ ∆ T/ 2 ◦ Φ ⋆ ∆ T/ 2 ( x n ) is symmetric: x n +1 / 2 = Φ ⋆ ∆ T/ 2 ( x n ) , x n +1 = Φ ∆ T/ 2 ( x n +1 / 2 ) e.g. Φ − ∆ T/ 2 ( x n +1 / 2 ) = x n , x n +1 = Φ ∆ T/ 2 ( x n +1 / 2 ) Use this observation to make a symmetric version of parareal: Φ ∆ T ≡ the standard parareal algorithm; we can write explicitly what is Φ ⋆ ∆ T . Remark: it is much more difficult to make symplectic a given integrator . . . The University of Chicago – p. 13

Symmetric parareal The standard form of the parareal algorithm writes x k +1 n +1 = G ∆ T ( x k +1 ) + F ∆ T ( x k n ) − G ∆ T ( x k n ) ( ⋆ ) n It is NOT a one-step integrator on x n ! The University of Chicago – p. 14

Symmetric parareal The standard form of the parareal algorithm writes x k +1 n +1 = G ∆ T ( x k +1 ) + F ∆ T ( x k n ) − G ∆ T ( x k n ) ( ⋆ ) n It is NOT a one-step integrator on x n ! Define U n = ( x 0 n , x 1 n , . . . , x K n ) Knowing U n , the parareal algorithm defines U n +1 : we write U n +1 = Φ ∆ T ( U n ) . The symmetrized form of Φ , that is Φ ∆ T/ 2 ◦ Φ ⋆ ∆ T/ 2 , reads U n = Φ − ∆ T/ 2 ( U n +1 / 2 ) , U n +1 = Φ ∆ T/ 2 ( U n +1 / 2 ) . The University of Chicago – p. 14

Symmetric parareal The standard form of the parareal algorithm writes x k +1 n +1 = G ∆ T ( x k +1 ) + F ∆ T ( x k n ) − G ∆ T ( x k n ) ( ⋆ ) n It is NOT a one-step integrator on x n ! Define U n = ( x 0 n , x 1 n , . . . , x K n ) Knowing U n , the parareal algorithm defines U n +1 : we write U n +1 = Φ ∆ T ( U n ) . The symmetrized form of Φ , that is Φ ∆ T/ 2 ◦ Φ ⋆ ∆ T/ 2 , reads U n = Φ − ∆ T/ 2 ( U n +1 / 2 ) , U n +1 = Φ ∆ T/ 2 ( U n +1 / 2 ) . In the specific context of the parareal algorithm ( ⋆ ), this two-step iteration yields the symmetric parareal integrator. The University of Chicago – p. 14

Symmetric parareal: some basic remarks - 1 At the initial iteration ( k = 0 ), we use the symmetric version of G : n +1 / 2 = G − 1 x k =0 − ∆ T/ 2 ( x k =0 x k =0 n +1 = G ∆ T/ 2 ( x k =0 ) , n +1 / 2 ) . n we next iterate over k :  � �  x k +1 G − 1 x k +1 − F − ∆ T/ 2 ( x k n +1 / 2 ) + G − ∆ T/ 2 ( x k  = n +1 / 2 ) , n n +1 / 2 − ∆ T/ 2   x k +1 G ∆ T/ 2 ( x k +1 n +1 / 2 ) + F ∆ T/ 2 ( x k n +1 / 2 ) − G ∆ T/ 2 ( x k = n +1 / 2 ) n +1 The flow is symmetric in the following sense: if, for any k and any n , n +1 , x k +1 ( x 0 n , . . . , x k n , x k +1 → ( x 0 n +1 , . . . , x k ) − n +1 ) n by the previous integrator, then ( x 0 n +1 , . . . , x k n +1 , x k +1 → ( x 0 n , . . . , x k n , x k +1 n +1 ) − ) n by the exact same algorithm reversing the time. The University of Chicago – p. 15