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 = f(x), x ∈ Rd Most algorithms are sequential in nature: xn+1 = Φ∆T (xn). 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 and G∆T = (ΦdT )∆T/dT 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: first, coarse propagation that yields

• xk=0

n

• n:

∀n, xk=0

n+1 = G∆T (xk=0 n

) Then, iterate over k ≥ 0: compute jumps (in parallel): Jk

n = F∆T (xk n) − G∆T (xk n)

sequential update to obtain

• xk+1

n

• n:

∀n, xk+1

n+1 = G∆T (xk+1 n

) + Jk

n

The fine solver is called only in the parallel part of the algorithm.

The University of Chicago – p. 3

The parareal paradigm - 2

xk+1

n+1 = G∆T (xk+1 n

) + F∆T (xk

n) − G∆T (xk n)

If F ≡ G, then xk+1

n+1 = F∆T (xk+1 n

): similar to a fine sequential integration. It turns out that, at iteration k, xk

k = (F∆T )(k) (x0),

as if the fine scheme was used sequentially all the way from x0 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)

1e-35 1e-30 1e-25 1e-20 1e-15 1e-10 1e-05 1 10 20 30 40 50 60 t fine integrator (sequential) parareal iteration 0 (coarse integrator) parareal iteration 3 parareal iteration 6 parareal iteration 10

Trajectory xk(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

• ver 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 (q, p), ˙ p = −∂H ∂q (q, p), q ∈ Rd, p ∈ Rd. Introducing x = (q, p) ∈ R2d, recast the dynamics as ˙ x = f(x) = J∇H(x), J =   Id −Id   . Separated case: H(q, p) = 1 2pT p + V (q). Then ˙ q = p, ˙ p = −∇V (q). Verlet algorithm (explicit, second order, symmetric, symplectic): pn+1/2 = pn − δt 2 ∇V (qn), qn+1 = qn + δt pn+1/2, pn+1 = pn+1/2 − δt 2 ∇V (qn+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) dt = J∇H(ϕt(x)), ϕt=0(x) = x. we will focus on the case H(q, p) = 1 2pT p + V (q), for which the flow is reversible: ∀t, ρ ◦ ϕt = (ϕt)−1 ◦ ρ 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) dt = J∇H(ϕt(x)), ϕt=0(x) = x. we will focus on the case H(q, p) = 1 2pT p + V (q), for which the flow is reversible: ∀t, ρ ◦ ϕt = (ϕt)−1 ◦ ρ 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 = 0, ˙ θ = ω(a), a ∈ Rd, θ ∈ Td. Thus a(q, p) ∈ Rd is left invariant by the dynamics.

The University of Chicago – p. 7

Backward error analysis for ODEs

exact flow exact flow numerical scheme ˙ x = f(x) ˙ x = fδt(x) ϕδt(x) Φδ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 = f(x) = J∇H(x), x = (q, p) ∈ R2d Standard parareal algorithm: xk+1

n+1 = G∆T (xk+1 n

) + F∆T (xk

n) − G∆T (xk n)

with F∆T = (Φδt)∆T/δt and G∆T = (ΦdT )∆T/dT 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-10 1e-08 1e-06 1e-04 0.01 1 100 10000 1e+06 1 10 100 1000 t Sequential coarse scheme Parareal it. #1 Parareal it. #4 Parareal it. #10 Sequential fine scheme

Error on the energy 1 2q2 + 1 2p2

The University of Chicago – p. 11

Outline

Adapt the parareal algorithm to the Hamiltonian framework by using symmetrization projection on the constant energy manifold M =

• (q, p) ∈ R2d; H(q, p) = H0
• The University of Chicago – p. 12

Symmetric parareal algorithm

Symmetrizing by hand a given integrator is easy: consider a numerical flow Φ∆T , that reads xn+1 = Φ∆T (xn) consider the adjoint of Φ∆T , defined by Φ⋆

∆T = Φ−1 −∆T

then the scheme xn+1 = Φ∆T/2 ◦ Φ⋆

∆T/2(xn) is symmetric:

xn+1/2 = Φ⋆

∆T/2(xn),

xn+1 = Φ∆T/2(xn+1/2) e.g. Φ−∆T/2(xn+1/2) = xn, xn+1 = Φ∆T/2(xn+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 xk+1

n+1 = G∆T (xk+1 n

) + F∆T (xk

n) − G∆T (xk n)

(⋆) It is NOT a one-step integrator on xn!

The University of Chicago – p. 14

Symmetric parareal

The standard form of the parareal algorithm writes xk+1

n+1 = G∆T (xk+1 n

) + F∆T (xk

n) − G∆T (xk n)

(⋆) It is NOT a one-step integrator on xn! Define Un = (x0

n, x1 n, . . . , xK n )

Knowing Un, the parareal algorithm defines Un+1: we write Un+1 = Φ∆T (Un). The symmetrized form of Φ, that is Φ∆T/2 ◦ Φ⋆

∆T/2, reads

Un = Φ−∆T/2(Un+1/2), Un+1 = Φ∆T/2(Un+1/2).

The University of Chicago – p. 14

Symmetric parareal

The standard form of the parareal algorithm writes xk+1

n+1 = G∆T (xk+1 n

) + F∆T (xk

n) − G∆T (xk n)

(⋆) It is NOT a one-step integrator on xn! Define Un = (x0

n, x1 n, . . . , xK n )

Knowing Un, the parareal algorithm defines Un+1: we write Un+1 = Φ∆T (Un). The symmetrized form of Φ, that is Φ∆T/2 ◦ Φ⋆

∆T/2, reads

Un = Φ−∆T/2(Un+1/2), Un+1 = Φ∆T/2(Un+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: xk=0

n+1/2 = G−1 −∆T/2(xk=0 n

), xk=0

n+1 = G∆T/2(xk=0 n+1/2).

we next iterate over k:      xk+1

n+1/2

= G−1

−∆T/2

• xk+1

n

− F−∆T/2(xk

n+1/2) + G−∆T/2(xk n+1/2)

• ,

xk+1

n+1

= G∆T/2(xk+1

n+1/2) + F∆T/2(xk n+1/2) − G∆T/2(xk n+1/2)

The flow is symmetric in the following sense: if, for any k and any n, (x0

n, . . . , xk n, xk+1 n

) − → (x0

n+1, . . . , xk n+1, xk+1 n+1)

by the previous integrator, then (x0

n+1, . . . , xk n+1, xk+1 n+1) −

→ (x0

n, . . . , xk n, xk+1 n

) by the exact same algorithm reversing the time.

The University of Chicago – p. 15

Symmetric parareal: some basic remarks - 2

     xk+1

n+1/2

= G−1

−∆T/2

• xk+1

n

− F−∆T/2(xk

n+1/2) + G−∆T/2(xk n+1/2)

• ,

xk+1

n+1

= G∆T/2(xk+1

n+1/2) + F∆T/2(xk n+1/2) − G∆T/2(xk n+1/2)

If F ≡ G, then the symmetrized form reads xk+1

n+1 = F∆T/2 ◦ F−1 −∆T/2(xk+1 n

) standard symmetrized version of the fine (= coarse) propagator. Formally taking the limit k − → +∞ yields x∞

n+1 = F∆T/2 ◦ F−1 −∆T/2(x∞ n ).

Limit of the algorithm in terms of parareal iterations: standard symmetrized form of F∆T . Note also that we never call F−1. All the expensive computations (involving F) are performed in parallel.

The University of Chicago – p. 16

Harmonic oscillator case: ˙ q = p, ˙ p = −q, symmetric algorithm

1e-08 1e-06 1e-04 0.01 1 100 1 10 100 1000 10000 t Sequential coarse scheme Parareal it. #1 Parareal it. #4 Parareal it. #10 Sequential fine scheme

Error on the energy 1 2q2 + 1 2p2

The University of Chicago – p. 17

Parareal integrator interpretation

A parareal integrator may be seen, at parareal iteration k, as an integrator of a system consisting of k + 1 identical replicas of the original system. Consider the symmetric parareal integrator, used on an integrable reversible Hamiltonian system: The first replica (k = 0) is integrated by the symmetric algorithm G∆T/2 ◦ G−1

−∆T/2: energy preservation.

Replicas are noninteracting: the system of replicas is an integrable reversible system, with energy = k

j=0 H(qj, pj)

At iteration k + 1, we have integrated an integrable reversible system with a symmetric algorithm: its energy is preserved. Hence, for all k, k

j=0 H(qj, pj) is preserved: the energy of each replica

H(qj, pj) is preserved. However, the non resonant condition is NOT satisfied: replicas . . . !

The University of Chicago – p. 18

Resonances

Consider a set of k identical replicas of the original system: before time discretization, they are resonant (identical!), but uncoupled: energy preservation, . . . OK the parareal time discretization introduces coupling. these resonances are present in the standard version as well as in the symmetric version of the algorithm. even if the symmetric algorithm is used, these resonances may impede preservation of the invariants (energy, . . . ) in the long-time limit. One possibility to avoid resonances: use a timestep dT (for G∆T ) that depends on k.    xk=0

n+1

= G(k=0)

∆T

(xk=0

n

) xk=1

n+1

= G(k=1)

∆T

(xk=1

n

) + F∆T (xk=0

n

) − G(k=1)

∆T

(xk=0

n

) + symmetrization.

The University of Chicago – p. 19

Fix resonances by using dTk

1e-10 1e-09 1e-08 1e-07 1e-06 1e-05 1e-04 0.001 0.01 0.1 10 100 1000 Error on energy, symetric parareal (shift), Verlet $t$

Kepler, symmetric parareal based on Verlet, dTk Error on the energy H(q, p) = pT p 2 − 1 |q|

The University of Chicago – p. 20

Going from dTk to a smaller dT with no k dependence

1e-10 1e-09 1e-08 1e-07 1e-06 1e-05 1e-04 0.001 0.01 0.1 10 100 1000 Error on energy, symetric parareal (shift), Verlet $t$ 1e-10 1e-09 1e-08 1e-07 1e-06 1e-05 1e-04 0.001 0.01 0.1 10 100 1000 Error on energy, symetric parareal (without shift), Verlet $t$ Initial coarse scheme Parareal iteration #1 Parareal iteration #4 Parareal iteration #10 At convergence

Kepler, symmetric parareal based on Verlet Error on the energy Left: dTk Right: dT smaller Both pictures with equal cost

The University of Chicago – p. 21

Another idea

On top of symmetrizing, let us project on the constant energy manifold M =

• x = (q, p) ∈ R2d; H(q, p) = H0
• ,

in a way that preserves symmetry.

The University of Chicago – p. 22

Another idea

On top of symmetrizing, let us project on the constant energy manifold M =

• x = (q, p) ∈ R2d; H(q, p) = H0
• ,

in a way that preserves symmetry. Take any symmetric integrator Ψ∆T , and consider xn → xn+1 defined by       

• xn

= xn + µ∇H(xn),

• xn+1

= Ψ∆T ( xn), xn+1 =

• xn+1 + µ∇H(xn+1),

with µ chosen such that xn+1 ∈ M. As the same Lagrange multiplier µ is used in the first and third lines, and as Ψ∆T is symmetric, the integrator xn → xn+1 is symmetric.

The University of Chicago – p. 22

Symmetric projection

      

• xn

= xn + µ∇H(xn),

• xn+1

= Ψ∆T ( xn), xn+1 =

• xn+1 + µ∇H(xn+1),

u

3

u

1

u

M

2

u

Need to solve a nonlinear problem to compute xn+1 and µ. In practice: use a Newton-like algorithm, and stop after a few iterations, or when the residu is small enough. In the following examples, stop after 2 iterations.

The University of Chicago – p. 23

Application to the parareal setting - 1

We need to start from a symmetric algorithm. We again introduce Un := (x0

n, x1 n, · · · , xK n ),

and consider the symmetric parareal algorithm previously obtained: Un+1 = Ψsym

∆T (Un),

where the map Ψsym

∆T is symmetric in the classical sense.

Introduce Hk(Un) := H(xk

n),

1 ≤ k ≤ K and consider a symmetric algorithm, based on the symmetric map Ψsym

∆T and

• n a symmetric projection on the manifold where all the energies Hk are

preserved.

The University of Chicago – p. 24

Application to the parareal setting - 2

• Un

= Un +

K

• k=1

µk∇Hk(Un)

• Un+1

= Ψsym

∆T (

Un) Un+1 =

• Un+1 +

K

• k=1

µk∇Hk(Un+1), with the Lagrange multipliers µk such that Hk(Un+1) = H0 for any 1 ≤ k ≤ K. Resulting algorithm: Initialization: at k = 0, use the symmetrized version of G: x0

n+1/2 = G−1 −∆T/2(x0 n),

x0

n+1 = G∆T/2(x0 n+1/2).

Set x0

n+1/2 = x0 n+1/2.

The University of Chicago – p. 25

Symmetric parareal integrator with projection

Assume iteration k is completed. We compute the solution at iteration k + 1 by: Set xk+1 = u0; For 0 ≤ n ≤ N − 1, compute xk+1

n+1 from xk+1 n

as             

• xk+1

n

= xk+1

n

+ µk+1∇H(xk+1

n

)

• xk+1

n+1/2

= G−1

−∆T/2

• xk+1

n

− F−∆T/2( xk

n+1/2) + G−∆T/2(

xk

n+1/2)

• xk+1

n+1

= G∆T/2( xk+1

n+1/2) + F∆T/2(

xk

n+1/2) − G∆T/2(

xk

n+1/2)

xk+1

n+1

=

• xk+1

n+1 + µk+1∇H(xk+1 n+1)

where µk+1 is such that H(xk+1

n+1) = H0.

As before, µk+1 and xk+1

n+1 are determined iteratively.

Jumps are precomputed in parallel, before the nonlinear iterations at step k + 1. All the expensive computations (fine propagator) are performed in parallel.

The University of Chicago – p. 26

Kepler problem (Nmax

iter = 2): energy preservation

1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 1e-04 0.01 0.01 0.1 1 10 100 1000 10000 t Sequential coarse scheme Parareal it. #1 Parareal it. #4 Parareal it. #10 Sequential fine scheme

Relative error on the energy H(q, p) (δt = 10−4, dT = 0.01, ∆T = 0.2). Energy is well preserved at any parareal iteration k ≥ 1, although we limit the number of iterations to solve the nonlinear projection step.

The University of Chicago – p. 27

Kepler problem (Nmax

iter = 2): angular momentum preservation

1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 1e-04 0.01 1 10 100 1000 10000 t Sequential coarse scheme Parareal it. #1 Parareal it. #4 Parareal it. #7 Parareal it. #10 Sequential fine scheme

Relative error on the angular momentum q ∧ p For all parareal iterations k ≥ 7, the angular momentum is preserved with a relative accuracy of 10−4 over the complete time range.

The University of Chicago – p. 28

Kepler problem (Nmax

iter = 2): trajectory accuracy

1e-08 1e-07 1e-06 1e-05 1e-04 0.001 0.01 0.1 1 10 1 10 100 1000 10000 t Sequential coarse scheme Parareal it. #1 Parareal it. #3 Parareal it. #5 Parareal it. #6 Sequential fine scheme

Error on the trajectory After only 6 iterations, errors on the trajectory are comparable to the errors on the trajectory of the fine scheme used sequentially on the complete time range.

The University of Chicago – p. 29

Kepler problem: summary of the results

energy is well preserved at any parareal iteration k ≥ 1 good preservation of the angular momentum (although we project on the manifold where only the energy is preserved)

• nly 6 parareal iterations are needed to reach the best possible accuracy
• n the trajectory.

These results are better than those obtained with the standard parareal algorithm or its symmetrized version: energy was badly preserved the standard parareal algorithm coupled with a standard, not symmetrized, projection step: better preservation of energy (remember N max

iter < ∞), of the angular momentum, and trajectory convergence is

reached for fewer parareal iterations. Both symmetrization and projection are useful.

The University of Chicago – p. 30

The outer solar system

H(q, p) = 1 2pT M −1p + V (q), q ∈ R18, p ∈ R18, with V (q) = −

• 1≤i<j≤6

G mi mj qi − qj, M = diag(mi). First particle ≡ Sun, with mass m1 ≥ 1000mi for any planet (i ≥ 2). We show here that the symmetric parareal algorithm with symmetric projection on the constant energy manifold again behaves well. that the coarse solver can be driven by a simplified dynamics without damaging the good properties of the algorithm.

The University of Chicago – p. 31

A simplified coarse integrator

Consider the simplified dynamics H(q, p) = 1 2pT M −1p + Vsimp(q), q ∈ R18, p ∈ R18, with Vsimp(q) = −

• 2≤j≤6

G m1 mj q1 − qj. In Vsimp, we only take into account the gravitational interaction between the Sun (at position q1) and the other planets (at position qj, 2 ≤ j ≤ 6), and we ignore the interaction between pairs of planets. m1 ≫ mi → The exact dynamics is a perturbation of the simplified dynamics. The simplified system is less expensive to simulate than the exact one, since Vsimp is a sum of 5 terms, whereas V is a sum of 15 terms.

The University of Chicago – p. 32

Numerical results (Nmax

iter = 2): energy preservation

1e-14 1e-12 1e-10 1e-08 1e-06 1e-04 100 1000 10000 100000 t Sequential coarse scheme Parareal it. #1 Parareal it. #4 Parareal it. #8 Parareal it. #10 Sequential fine scheme

Relative error on the energy (δt = 10−2, dT = 50, ∆T = 200). For k ≥ 8, energy is preserved up to the prescribed tolerance.

The University of Chicago – p. 33

Numerical results: trajectory accuracy

1e-10 1e-08 1e-06 1e-04 0.01 1 100 100 1000 10000 100000 t Sequential coarse scheme Parareal it. #1 Parareal it. #5 Parareal it. #8 Parareal it. #10 Parareal it. #14 Sequential fine scheme

Error on the trajectory For k ≤ 5, large trajectory error. At k = 15, the trajectory error is comparable to the error made by the fine propagator used sequentially.

The University of Chicago – p. 34

Numerical results: angular momentum preservation

1e-16 1e-14 1e-12 1e-10 1e-08 1e-06 1e-04 0.01 1000 10000 100000 t Sequential coarse scheme Parareal it. #1 Parareal it. #5 Parareal it. #10 Parareal it. #15 Sequential fine scheme

Relative error on the angular momentum Good preservation of the angular momentum (error < 1 %), although the trajectory may be wrong.

The University of Chicago – p. 35

Numerical results: qualitative trajectory

• 15
• 10
• 5

5 10 15

• 30
• 20
• 10

10 20 30 40 Sun Jupiter Saturn Uranus Neptune Pluto

• 15
• 10
• 5

5 10 15

• 30
• 20
• 10

10 20 30 40 Sun Jupiter Saturn Uranus Neptune Pluto

• 15
• 10
• 5

5 10 15

• 30
• 20
• 10

10 20 30 40 Sun Jupiter Saturn Uranus Neptune Pluto

• 15
• 10
• 5

5 10 15

• 30
• 20
• 10

10 20 30 40 Sun Jupiter Saturn Uranus Neptune Pluto

k = 1 (top left), k = 5 (top right), k = 10 (bottom left), k = 15 (bottom right). For k = 5, the trajectory is quantitatively wrong but qualitatively correct! For k ≥ 10, the trajectory is quantitatively correct.

The University of Chicago – p. 36

Computational speed-up

With our parameters, the computation cost is mostly due to the fine-scale

• integrations. Hence, the speed-up is

G = T K∆T = number of processors number of parareal iterations For the outer-solar system test-case, G ≈ 66 provided we have 1000 processors.

The University of Chicago – p. 37

Symplectic parareal algorithm

Unclear (no standard way to make a given algorithm symplectic). A first possible approach by G. Bal and Q. Wu (2008): build a generating function S by interpolating parareal results: in practice, for some chosen φi(q, P), S(q, P) =

• i

ai φi(q, P), ai to determine. derive from this generating function a symplectic algorithm: p = P + ∂S ∂q (q, P), Q = q + ∂S ∂P (q, P). How to go further?

The University of Chicago – p. 38

Conclusions

the standard parareal algorithm is neither symplectic nor symmetric, even if the building blocks are: this creates issues in the long-time behaviour. easy to symmetrize, not enough due to resonances . . . symmetrize + symmetric projection → good results! Both ingredients are needed. Some open questions: Symplectic parareal algorithm? Understand better the behaviour of the integrators through their modified equations (work in progress).

• X. Dai, CLB, F

. Legoll, Y. Maday, Symmetric parareal algorithms for Hamiltonian

systems, arXiv preprint 1011.6222, to appear in M2AN.

Related work: M. Gander, E. Hairer, Analysis for parareal algorithms applied to

Hamiltonian differential equations, preprint

The University of Chicago – p. 39