IV (2)

COMPUTATIONAL FLUID DYNAMICS (CFD) IV (2) The Analysis of Numerical Schemes (2) 11. Iterative methods for algebraic systems Reima Iwatsu, e-mail : iwatsu@cck.dendai.ac.jp Winter Semester 2007, Graduate School, Tokyo Denki University 19, 2007 11 16F, 1601 1

IV (2) 10. 10.1 10.2 10.3 11. 11.1 11.2 11.3 11.4 11.5 2

11. Iterative methods for algebraic systems CFD 3D RANS (7 100 = 700 ),,, :, Varga (1962), Dahlquist & Bjork (1974), Saad (2003) Dongarra 3

11.1 11.1.1 2 5 u = f 0 x, y L (1) u = g on the bounradies (2) (u i+1,j 2u ij + u i 1,j ) + (u i,j+1 2u ij + u i 1,j ) = f ij x 2 (3) M j 1 0 y 6 6 6 6 0 1 i M 6 x 4

U = (u 11, u 21,, u M,1, u 1,2, u 2,2,, u 1,j, u 2,j,, u M,j, ) (4) 4 1 1 1 4 1 1 1 4 1 1 (i) S SU = F x 2 + G Q (5) (ii) S (, ) N = M 2, I = i + (j 1)M N J=1 u 11. u i1. = x2 f 11. f i1. g 10 + g 01. g i0. s IJ u J = Q I I = 1,, N (7) (6) 5

11.1.2 / ( ) u n ij u n+1 un+1 ij (n ) ij = 1 4 (un i+1,j + un i 1,j + un i,j+1 + un i,j 1 ) + 1 4 qn ij (8) u n+1 i = 1 N q n s i s ij u n j (9) ij j=1,j i 6 6 6 6 : Unknown at iteration n + 1 : Known from iteration n 6

D:, E:, F : S = D + E + F (10) DU n+1 = Q n (E + F )U n (11) (D + E + F )U = Q (12) - R n R n R(U n ) = SU n + Q n = (D + E + F )U n + Q n (13) D U n = R n (14) U n U n+1 U n (15) 4(u n+1 ij u n ij ) = (un i+1,j + un i 1,j + un i,j+1 + un i,j 1 4un ij ) + qn ij (16) 7

u n+1 ij u n+1 I = 1 s II = 1 4 (un+1 i+1,j + un i 1,j + un i,j+1 + un+1 i,j 1 ) + 1 4 qn ij I 1 qi n N s IJ u n+1 J s IJ u n J J=1 J=I+1 (17a) (17b) (D + E)U n+1 = Q n F U n (18) (D + E) U n = R n (19) 6 6 : Unknown at iteration n + 1 : Known from iteration n + 1 : Known from iteration n 8

11.1.3 n e n 0, e n = U n Ū Ū : (20) S S = P + A (21) P U n+1 = Q n AU n (22) P U n = R n (23) P : ( (22) ) SU n+1 = Q n or S U n = R n (24) (25) (22) SŪ = (P + A)Ū = Q (25) 9

P e n+1 = Ae n (26) e n+1 = (P 1 A)e n = (P 1 A) n e 1 = (1 P 1 S) n e 1 (27) n e n+1 0 P 1 A ( ) G G = 1 P 1 S (28) ρ(g) 1 (29) λ J (G) 1 for all J (30) S P ˆP 1 O(N) ˆP S 10

R e R n = Se n (31) 0 0 ( ) R R U n Ū SU 2 x = 10 2 U 10 4 U n+1 = (1 P 1 S)U n P 1 Q = GU n P 1 Q (32) e n+1 = Ge n = G n+1 e 0 (33) n 11

G = D 1 (E + F ) = 1 D 1 S G = (D + E) 1 F = 1 (D + E) 1 S, : S S : S S, (2D S) (S 2D) : S S, (2D S) 12

n (33) Varga (1962) e n ( e n ) 1/n e 0 (34) e n G n 1/n (35) log s n 1 (10 ) s lim n Gn 1/n = ρ(g) (36) ρ(g) < 1 13 ( ) 1 n s = ρ(g) (37) 10 n 1/ log ρ(g) (38)

11.1.4 G (P 1 S ) U n = J=1 (λ J ) n U 0 V (J) J=1 Q J Ω J [1 (λ J ) n ]V (J) (39) V (J) : S, U 0 :, Ω J : S, S, G, (λ J : G ) lim n U n = J=1 U = S 1 Q Q J Ω J V (J) (40) Ū = J=1 U J V (J) (41) SŪ = J=1 Ω J U J V (J) = J=1 Q J V (J) (42) U J = Q J Ω J (43) 14

n (e 0J : ) e n = J=1(λ J ) n e 0J V (J) (44) V (J) ( ) R n = J=1 (λ J ) n e 0J SV (J) = J=1(λ J ) n Ω J e 0J V (J) (45) L 2 R n L2 = J=1[(λ J ) n Ω J e 0J ] 2 (46) n (λ J ) (λ J ) 1 R b a n 15

G S λ J = λ(ω J ) D G J = 1 D 1 S (47) λ J = 1 1 d J Ω J (S) (48) 16

11.1.5 S, G S Ω = 4(sin 2 φ x /2 + sin 2 φ y /2) φ x = lπ/m φ y = mπ/m J = l+(m 1)M l, m = 1,, M (49) G J λ(g J ) = 1 (sin 2 φ x /2 + sin 2 φ y /2) = 1 2 (cos φ x + cos φ y ) (50) (φ x, φ y 0) (λ 1) (Ω- ) x = y ρ(g J ) = cos πm (51) ρ(g J ) 1 π2 2M 2 = 1 π2 2N = 1 O( x2 ) (52) 17

(4 e iφ x e iφ y )λ GS = e iφ x + e iφ y (53) 2 (Young, 1971) λ(g GS ) = 1 4 (cos φ x + cos φ y ) 2 = λ(g J ) 2 (54) ρ(g GS ) = ρ(g J ) 2 = cos 2 π/m (55) M ρ(g GS ) 1 π2 M 2 = 1 π2 N 2 1 N/(0.43π 2 ) Varga (1962) (56) 18

10.2 U n = U n+1 U n (Frankel & Young, 1950) U n+1 ω U n+1 = ωū n+1 + (1 ω)u n (57) U n = U n+1 U n with U n = ω U n (58) ω 19

10.2.1 u n+1 ij = ω 4 (un i+1,j + un i 1,j + un i,j+1 + un i,j 1 + qn ij ) + (1 ω)un ij (59) ( ) DU n+1 = ωq n ω(e+f )U n +(1 ω)du n = ω(su n +Q n )+DU n (60) D U n = ωr n (61) G J (ω) = (1 ω)i + ωg J = 1 ωd 1 S (62) ρ(g J (ω)) 1 ω + ωρ(g J ) < 1 (63) 20

0 < ω < 2 1 + ρ(g J ) (64) λ(g J (ω)) = (1 ω) + ωλ(g J ) (65) ω ω = 1/(1 λ(g J )) λ(g J ) = 0 λ 1, ω < 1 λ 1, ω ω max 21

ω ω opt 1 + ω opt (1 λ min ) = 1 ω opt (1 λ max ) (66) λ min λ max ω opt = (50) ω opt = 1 2 2 (λ min + λ max ) (67) λ max = λ min = cos π/m (68) λ(g J (ω)) Unstable 1 1 1 λ min ω opt 1 λ max ω 22

10.2.2 : SOR (Successive Overrelaxation) ( ) ū n+1 ij = ω 4 (un i+1,j + un+1 i 1,j + un i,j+1 + un+1 i,j 1 + qn ij ) + (1 ω)un ij u n+1 ij = ωū n+1 ij + (1 ω)u n ij (69) D U n = R n E U n (70) - ( (D + ωe)) (D + ωe) U n = ωr n (71) G SOR (ω) = 1 ω(d + ωe) 1 S = (D + ωe) 1 [(1 ω)d ωf ] (72) 23

det G SOR (ω) = det(i + ωd 1 E) 1 det[(1 ω)i ωd 1 F ] = I det[(1 ω)i ωd 1 F ] = (1 ω) N (73) ρ(g SOR (ω)) N det G SOR (ω) = (1 ω) N (74) 1 (1 ω) ρ(g SOR (ω)) < 1 (75) SOR (Young, 1971) 0 < ω < 0 < ω < 2 (76) 2 1 + ρ(g GS ) (77) 24

( ) D S = 1 F F T (78) D 2 D 1, D 2 :, λ(g SOR (ω)) = λ(ω) i) SOR λ(ω) = 1 ω + ωλ 1/2 (ω)λ(g J ) (79) ii) ω = 1 λ(ω = 1) λ(g GS ) = λ 2 (G J ) (80) iii) ω opt = 2 1 + 1 ρ 2 (G J ) (81) ρ(g SOR (ω opt )) = ω opt 1 (82) 25

ω opt = ρ(g SOR (ω opt )) = ρ(g J ) = cos π/m (83) 2 1 + sin π/m 2(1 π M + π2 M 2) (84) 1 sin π/m 1 + sin π/m 1 2π M + O( 1 M 2) (85) 1 2.3M/π N SOR k kn N ρ(g) = lim n e n+1 e n n (79) ρ(g J ) (81) ω (86) 26

10.2.3 SSOR (Symmetric Successive Overrelaxation) SOR U n+1 SOR (D + ωe) U n+1/2 = ωr n (87) (D + ωf ) U n = ωr(ū n+1/2 ) (88) U n+1 = Ū n+1/2 + Ū n+1/2 = U n + U n = U n + U n+1/2 U n+1/2 + U n (90) SSOR SOR ω opt 2 1 + 2(1 ρ 2 (G J )) (91) SOR 2 ;2 27

10.2.4 SLOR (Successive Line Overrelaxation) U n SLOR (VLOR) ū n+1 ij = ω 4 (un i+1,j + un+1 i 1,j + ūn i,j+1 + ūn i,j 1 + 1 4 qn ij ) u n+1 ij = ωū n+1 ij + (1 ω)u n ij (92) 4 Uij n U i,j 1 n U i,j+1 n ω U i 1,j n = ωrn ij (93) SLOR ω opt 2 6 : Unknown at iteration n + 1 : Known from iteration n : Known from iteration n + 1 28

- : I = i + (j 1)M, : I : Red point : Black point 29

10.3 ( -form) n P U n τ = ω(su n + Q n ) = ωr n (94) τ P du dt = ω(su + Q) (95) ( ) P 1 G = 1 ωτp 1 S (96) ρ(g) 1 ( ωτp 1 )S ( ) P = (ωτs), G = 0 1 P/ωτ S P 30

10.3.1 P U n τ = ω(su n + Q n ) = ωr n P/τ = 1 ( S ) U n+1 = U n + ω(su n + Q n ) = (1 + ωs)u n + ωq n G R U n + ωq n (97) ( ) G R = 1 + ωs G R S Ω J λ J = 1 + ωω J (98) 0 < ω < 2 = 2 Ω J max ρ(s) (99) 31

λ 1 Ω ω 1 max opt ω opt = Unstable Ω min ( ) S ρ(g R ) = 1 2 Ω J max + Ω J min (100) 2 Ω J max = κ(s) 1 Ω J max + Ω J min κ(s) + 1 ω (101) κ(s) = Ω J max Ω J min (102) Ω J max = 8 Ω J min = 8 sin 2 π 2M 2π2 N 32 (103)

κ(s) = 1 sin 2 π 4N 2M π 2 (104) M 50 50 κ(s) 1000 ( ) ρ(g R ) = 1 2π2 N = 1 2 κ(s) (105) ω 1/ Ω J max 1/ Ω J min non-stationary cf. stationary 33

10.3.2 ADI (Alternating Direction Implicit) P ADI (1 τs x )(1 τs y )(1 τs z ) = τω(su n + Q n ) (106) (1 τs x ) U = τω(su n + Q n ) (1 τs y ) U = U (1 τs z ) U = U (107) ω, τ ADI N log N ( ) ω opt 2 τ opt = τ 1 ΩJ min Ω J max (108) λ(g ADI ) = (1 4τ 1 sin 2 φ x /2)(1 4τ 1 sin 2 φ y /2) (1 + 4τ 1 sin 2 φ x /2)(1 + 4τ 1 sin 2 φ y /2) (109) 34

10.3.3 P P U n τ = ω(su n + Q n ) = ωr n U n+1 = U n τωp 1 (SU n + Q n ) B = τp 1 S G = 1 τp 1 S, G ω opt = λ(g) = 1 ωλ(τp 1 S) (110) 2 λ min (τp 1 S) + λ max (τp 1 S) ρ(g) = κ(τp 1 s) 1 κ(τp 1 s) + 1 (111) = 10 5 10 6 : τp 1 S P 35

(Choleski) Meijerink & Van der Vorst U n+1 = U n τωp 1 (SU n + Q n ) (112) Strongly Implicit Procedure (SIP) Stone, Schneider & Zedan λ(g) = 1 ωλ(τp 1 S) (113) CG (Conjugate Gradient) Reid, Concus, Kershaw ω opt = 2 λ min (τp 1 S) + λ max (τp 1 S) GMRES (Generalized Minimum Residual) (114) Saad & Schultz ρ(g) = κ(τp 1 s) 1 κ(τp 1 s) + 1 (115) 36

10.4 S(U) = Q (116) S(U n+1 ) = S(U n + U) = S(U n ) + (Jacobian) ( ) J(U) Ū ( ) S U ( ) S U = Q (115) U (116) J(U) U n = R n (117) Ū = U n J(U) 1 R n (118) 37

e n = U n Ū J(U)e n = R n (119) (117) J(U) U n = R n P/τ P τ U n = R n (120) e n G e n+1 = U n+1 Ū = en + U = e n τp 1 R n = (1 τp 1 J)e n Ge n P 1 (121) G = 1 τp 1 J (122) 38

10.5 Brandt (1972,1977,1982), Thomas et al. (2003), Briggs (2000) = = G 1 = S Ω J min φ x = φ y = π/m S (λ J ) n (G ) λ J = 0.5 n 2 n λ J = 0.999 1 10 23000 Error after smoothing x 39

10.5.1 (smoothing) µ = max π/2 φ π : φ x = φ y = π 1 λ(g J ) = 1, +1/2 λ(g) (123) µ(g J (ω)) = max[ 1 2ω, 1 ω/2 ] (124) ω = 4/5 µ = 0.6 0.6 0.5 0.447 0.25 0.25 0.048 40

10.5.2 CGC (Coarse Grid Correction) h:, H:, ex. H = 2h (, residual form) S h U h = Q h (125) P U h = R h (126) 1. (restriction) 2. R H = I H h R h (127) S H U H = Q H (128) S H U H = R H (129) 3. U H (prolongation) U H = I h H U H (130) 41

10.5.3 Two-grid Iteration Method h, H 1. Uh n S 1 n 1 2. U n+1 h = Uh n + U h 3. U n+1 h S 2 n 2 S 1 6 6 n 1 S 2 n 2 6 coarse grid correstion h ( N) N ( N) 42

10.5.4 FAS (Full Approximation Scheme) Ih H, Îh H S h (U h ) = Q h (131) S h (U h + U h ) S h (U h ) = R h (132) R H = Ih H R h (133) U H = ÎH h U H (134) S H (U H + U H ) S H (U H ) = R H (135) 43

( ) P G = 1 P 1 S G S O(N) 44

n ˆ 1 ˆ 2 n ˆ ( 6 ), I, Toro ˆ Chorin (1967,1968), Harlow & Welch (1965) Patankar & Spalding (1972) 45