Time Evolution Time-evolution problems are widely solved in - PowerPoint PPT Presentation

・・・ ・・・ ・・・ Aug. 30 - Sep. 2, 2011 ParCo2011, Ghent, Belgium Time Evolution Time-evolution problems are widely solved in scientific Parareal Acceleration of simulations described by discretized differential equations. Matrix Multiplication Parallel technique is usually applied through domain decomposition in the space direction , where quantity Toshiya Takami and Akira Nishida on the surface of each domain must be shared with its Kyushu University, Japan neighbors. On the other hand, efficient parallelism by the time- domain decomposition seems difficult because of its severe dependency on the previous state . 1 3 Contents Time-domain Parallelism Introduction: Time-domain Decomposition Time-evolution is usually defined by strictly dependent relations, which is difficult to be parallelized What is Parareal? x k +1 = F k ( x k ) Parareal-in-Time Algorithm as a Perturbation ‘’Parareal-in-Time’’ is one of the time-domain methods that Application: Series by Matrix-Vector Multiplications can be used in spite of such strict dependency. J-L. Lions, Y. Maday, and G. Turinici, C. R. Acad. Sci., Ser. I: Math. 232, 661–668 (2001). Convergence Property Domain Decomposition in Time Direction Speed-up Ratio and Efficiency x 0 x 1 x 2 x 3 x 4 x 5 x k − 3 x k − 2 x k − 1 x k Discussion: Applicability to Other Linear Calculations F 1 F 2 F k − 2 F k − 1 F 0 F 3 F 4 F 5 F k − 4 F k − 3 Conclusion t = 0 1 3 4 5 k − 3 k − 2 k − 1 k 2 2 4

Sequence Defined by What is “Parareal”? Parareal-in-Time “Its main goal concerns real time problems, hence the Instead of the sequence defined by , { x k } x k +1 = F k ( x k ) { x ( r ) consider an approximated one by an iterative relation, proposed terminology of ‘parareal’ algorithm. k } x ( r +1) = G k ( x ( r +1) ) + F k ( x ( r ) k ) − G k ( x ( r ) k ) J-L. Lions, Y. Maday, and G. Turinici, C. R. Acad. Sci., Ser. I: Math. 232, 661–668 (2001). k +1 k where is a coarse evolution that approximates the “Parareal is not the first algorithm to propose the solution of G k ( · ) original one , and r is the order of the approximation. F k ( · ) evolution problems in a time-parallel fashion” M. J. Gander and S. Vandewalle, SIAM J. Sci. Comput. 29, 556-578 (2007). The exact operators can be calculated in parallel. F k ( · ) Multiple Shooting Methods with Newton’s Iterations In the limit , it has been shown that the approximate r → ∞ sequence converges to the original one: { x ( r ) k } → { x k } Space-Time Multigrid Method which is called the Parareal-in-Time method. J-L. Lions, Y. Maday, and G. Turinici, C. R. Acad. Sci., Ser. I: Math. 232, 661–668 (2001). 5 7 Parareal-in-Time for Parareal-in-Time Scientific Applications as a Perturbation (1) Since the first proposal by J-L. Lion, et al ., various time- While convergence of the Parareal-in-Time method is evolution problems has been analyzed. already shown as the multiple-shooting method, a simple explanation is given here using a linear problem defined by PDE with a fluid-structure coupling: bounded operators: C. Farhat and M. Chandesris, Int. J. Numer. Meth. Engng. 58 , 1397-1434 (2003). consider a sequence of vectors { x k } Molecular Dynamics with Quantum Iteraction: defined by x k +1 = F k x k L. Baffico, S. Bernard, Y. Maday, G. Turinici, and G. Zérah, Phys. Rev. E 66 , 057701 (2002). = [ G k + ( F k − G k )] x k Quantum Control Problem: where we assume and are bounded. G k F k − G k Y. Maday, G. Turinici, Int. J. Quant. Chem. 93 , 223-228 (2003). Symplectic Integrator: Convergence of the sequence above can be explained as a G. Bal and Q. Wu, LNCSE 60, 401–408 (Springer, 2008). perturbation expansion and the spectral radius of operators. 6 8

Parareal-in-Time Series defined by Matrix- as a Perturbation (2) Vector Multiplication For linear problems, the operator is separated into a A linear time-evolution defined by matrix multiplications, F j coarse operator and its correction , x k = ( G + ε F ) x k − 1 = ( G + ε F ) k x 0 G j F j − G j k − 1 can be represented as an expanded sum in the order of . ε � x k = [ G k − 1 + ( F k − 1 − G k − 1 )] x k − 1 = [ G j + ( F j − G j )] x 0 j =0 { x ( r ) If we introduce as r-th order approximation of k } { x k } The third order approximation: � G k + ε x ( r +1) = [ G k + ( F k − G k )] x ( r +1) x (3) � = (Terms with a F ) k +1 k k = G k x ( r +1) + ( F k − G k ) x ( r +1) + ε 2 � k k (Terms with two F ’s) ≈ G k x ( r +1) + ( F k − G k ) x ( r ) k k � + ε 3 � (Terms with three F ’s) x 0 Thus, the Parareal-in-Time calculations for linear problems can be analyzed as a higher order perturbation . 9 11 How to Calculate Residual Errors the Parareal Sequence from Spectral Radius The parallel procedure of this calculation is: Errors can be estimated from the rest of expanded terms: k are parallel in k , but are sequantial. � ε j � F k − G k G k (Terms with F j ) x 0 j = r +1 k -direction is first, r -direction is sequential Error from the exact sequence is bounded by ... x 0 x 1 x 2 x 3 x 4 x k − 3 x k − 2 x k − 1 x k � � � x ( r ) Time − x 0 k � j � � ... � � � ερ ( F ) k r = 0 k � � G G G G ρ ( G k ) | x 0 | ≤ j ρ ( G ) F − G F − G F − G F − G F − G x (1) j = r +1 ... r = 1 k − 2 G G G G G By the Stirling’s formula, the magnitude of r -th term is F − G F − G F − G F − G F − G x (2) ... � r = 2 k − 1 � � � ερ ( F ) � r � ek � r � ερ ( F ) � r k k G G G G G = r ρ ( G ) 2 π ( k − r ) r r ρ ( G ) F − G F − G F − G F − G F − G x (3) ... r = 3 G G G k G G Order of Perturbation 10 12

Example 1: Implementation by MPI Real-Symmetric Matrix Consider a sequence of the real vector defined by a real- We can implement the Parareal-in-Time iteration on parallel { x k } computers by the use of MPI_Send/Recv. symmetrix matrix , where G is diagonal and F is a G + ε F 1 1 random matrix with elements satisfying � | F ii | 2 � = � | F ij | 2 � = 2 N , The original number of calculations, , is compared K ( T g + T f ) 4 N with the parareal one, . KT g + T c ( P − 1) + rK ( T g + T f ) /P N =1024, � =0.01 and in this case, (a) Iterate r times we can set 10 -5 ρ ( F ) ≈ ρ ( G ) ≈ 1 G F G F G F G -7 r=5 Normalized Error 10 Process 1: P resources -9 10 G F G F G F G dashed: error of ( r +1)-th term Process 2: -11 10 10 -13 G F G F G F G Process 3: eigenvalues r=10 -15 10 of F G F G F G F G -17 10 Process 4: r=15 KT g + T c ( P − 1) rK ( T f + T g ) /P 10 -19 Computational Time -1 1 0 50 100 150 200 E. Aubanel, Parallel Computing 37 , 172-182 (2011). Length of Sequence 13 15 Example 2: Estimation of Unitary Matrix Speed-up Ratio Unitary time-evolution is often used in quantum mechanics: Speed-up Ratio will be represented in the function, � ∆ t � x k − 1 ≈ ( I − i ε H ) k x 0 K ( T g + T f ) P x k = exp i � H S ( r, K, P ) = = KT g + T c ( P − 1) + rK ( T f + T g ) r + TP + P ( P − 1) t 1 where H is a Hermitian matrix satisfying P K � | H ii | 2 � = � | H ij | 2 � = 4 N where and . T ≡ T g / ( T g + T f ) t ≡ T c / ( T g + T f ) (b) 10 -5 N =1024, � =0.01 and in this case, -6 10 r=10 Property of the Speed-up: we can set Normalized Error -7 10 ρ ( H ) ≈ 1 Limit of Speedup -8 1/T 10 r=20 Speedup Ratio Max. Efficiency: small P Max. Efficiency (1/r) dashed: error of ( r +1)-th term K>>P>>1 K~P -9 10 10 -10 r=30 eigenvalues -11 Max. Speed-up: 10 r=40 K ≫ P ≫ 1 of H 10 -12 r=50 -13 10 0 0 200 400 600 800 1000 0 K Total Number of Ranks (P) Length of Sequence -1 1 14 16