14.Juni2007 SebastianPlitzko ADI-Verfahren

ADI-Verfahren

Sebastian Plitzko

Seminar Numerik und wissenschaftliches Rechnen Technische Universität Darmstadt

14. Juni 2007

Inhaltsverzeichnis

1 Einführung

2 ADI-Verfahren

3 Sequentieller ADI-Algorithmus

4 Paralleler ADI-Algorithmus

5 Zusammenfassung

Inhaltsverzeichnis

1 Einführung

2 ADI-Verfahren

3 Sequentieller ADI-Algorithmus

4 Paralleler ADI-Algorithmus

5 Zusammenfassung

Inhaltsverzeichnis

1 Einführung

2 ADI-Verfahren

3 Sequentieller ADI-Algorithmus

4 Paralleler ADI-Algorithmus

5 Zusammenfassung

Inhaltsverzeichnis

1 Einführung

2 ADI-Verfahren

3 Sequentieller ADI-Algorithmus

4 Paralleler ADI-Algorithmus

5 Zusammenfassung

Inhaltsverzeichnis

1 Einführung

2 ADI-Verfahren

3 Sequentieller ADI-Algorithmus

4 Paralleler ADI-Algorithmus

5 Zusammenfassung

Übersicht

1 Einführung

2 ADI-Verfahren

3 Sequentieller ADI-Algorithmus

4 Paralleler ADI-Algorithmus

5 Zusammenfassung

Grundlegendes

implizite Behandlung mehrdimensionaler Probleme nicht einfach ADI≡alternate direction implicit

speziell für das Lösen parabolische DGL entwickelt Algorithmus ist dimensionsabhängig

Matrix K wird in eine Summe von Matrizen zerlegt, die in tridiagonalform gebracht werden können

Übersicht

1 Einführung

2 ADI-Verfahren

3 Sequentieller ADI-Algorithmus

4 Paralleler ADI-Algorithmus

5 Zusammenfassung

ADI-Verfahren

ADI-Verfahren in 2D (u

= u

+ u

)

1. Halbschritt u⁽ⁿ⁺

1 2)

k,l =u_k,l⁽ⁿ⁾+ ^δt₂

"

u⁽ⁿ⁺

12) k+1,l +2u⁽ⁿ⁺

12) k,l +u⁽ⁿ⁺

12) k−1,l

h²_x + ^u

(n)

k+1,l+2u_k,l⁽ⁿ⁾+u⁽ⁿ⁾_k−1,l h_y²

#

2. Halbschritt u⁽ⁿ⁺¹⁾_k,l =u⁽ⁿ⁺

1 2) k,l + ^δt₂

"

u⁽ⁿ⁺

1 2) k+1,l +2u⁽ⁿ⁺

1 2) k,l +u⁽ⁿ⁺

1 2) k−1,l

h²_x +^u

(n+1)

k+1,l+2u_k,l⁽ⁿ⁺¹⁾+u⁽ⁿ⁺¹⁾_k−1,l h_y²

#

zwei Teilschritte ^δt₂

zunächst wird nur x-Richtung implizit integriert danach wird die y-Richtung weiter integriert ein voller Schritt ist vollzogen

Verfahren

Verfahrensvorschrift

u_k,l⁽ⁿ⁺¹⁾ =u_k,l⁽ⁿ⁾+δt^u

(n+1 2) k+1,l +2u⁽ⁿ⁺

1 2) k,l +u⁽ⁿ⁺

1 2) k−1,l

h²_x

+^δt₂

u⁽ⁿ⁺¹⁾_k+1,l+2u_k,l⁽ⁿ⁺¹⁾+u⁽ⁿ⁺¹⁾_k−1,l

h²_y +^u

(n)

k+1,l+2u⁽ⁿ⁾_k,l+u_k−1,l⁽ⁿ⁾ h²_y

Integrationen:

1 x-Richtung Tangententrapezregel: Fehlerordnungδt²

2 y-Richtung Sehnentrapezregel: Fehlerordnungδt²

=⇒ Gesamtfehlerordnung δt³

Stabilitätsanalyse für Wärmeleitungsproblem

Courant-Friedrichs-Levy Stabilitätsanalyse mit dem Ansatz u_k,l⁽ⁿ⁾ =λⁿe^{i(qxkhx+qylhy)}

Ansatz folgt aus Eigensystem des Laplaceoperators für den ersten Halbschritt

λ

1 2

1 =1+δt 2

λ

1 2

1

e^iqxhx +e^−iqxhx −2

h²_x +e^iqyhy +e^−iqyhy −2 h²_y

=⇒ λ

1 2

1 =1+2δt

"

λ

1 2

1

sin^{2 qx}₂^hx

h²_x +sin^{2 qy}₂^hy h²_y

#

Stabilitätsanalyse für Wärmeleitungsproblem

λ

1 2

1 =

1−2δt^sin2

qy hy 2

h²_y

1+2δt^sin2_h^{qx hx}2² x

λ

1 2

2 = 1−2δt^sin2_h^{qx hx}2² x

1+2δt^sin

2qy hy 2

h²_y

numerische Diffusionsprobleme δt >0 für~q 6=0−→ |λ₁|¹²|λ₂|¹² <1

für positive Schritte unbedingt stabil

unbedingte Stabilität ist spezielle Eigenschaft des 2D ADI-Verfahrens

Übersicht

1 Einführung

2 ADI-Verfahren

3 Sequentieller ADI-Algorithmus

4 Paralleler ADI-Algorithmus

5 Zusammenfassung

Sequentieller Algorithmus

Algorithmus

Chooseu⁰ r :=f −K ·u⁰ σ:=σ₀ := (r,r) k :=0

whileσ > tol²·σ₀ do

(H +ρ_k+1I)·u^k∗ = (ρ_k+1I −V)·u^k +f (V +ρk+1I)·u^k+1 = (ρk+1I −H)·u^k∗+f

r :=f −K ·u^k+1 σ = (r,r)

k :=k+1 end

Konvergenzbeschleunigung

bekannt sind die Eigenwerte bzw. Eigenvektoren der Matrizen Kµjm =λjmµjm

µ_jm =sin(jπy)sin(^mπx₂ ) Hµjm =λmµjm

Vµ_jm =λ_jµ_jm λjm =λj +λm

Konvergenzbeschleunigung

Zur Beschleunigung der Konvergenz wird eine Zeitentwicklungsmatrix betrachtet

Aus dem Algerithmus:

Tρ= (V +ρI)⁻¹(H −ρI)(H +ρI)⁻¹(V −ρI) somitTρµ_jm = ^(λ_(λ^{j−ρ)(λm−ρ)}

j+ρ)(λm+ρ)µ_jm

wir erhalten ein Minimierungsproblem mit Variable ρ

Konvergenzbeschleunigung

0< α < λ_j , λ_m ≤β c = ^α_β , δ = (√

2−1)² , n = [logδc] +1 ρ_j =βc^n−1j−1

Für ein(Nx +1)x(Ny +1)-Gitter erhalten wir:

α= _h¹2(2−2 cos(πh))≈π β = _h¹2(2+2 cos(πh))≈ _h⁴2

ρ_j ≈ _h⁴2δ^j−1

Übersicht

1 Einführung

2 ADI-Verfahren

3 Sequentieller ADI-Algorithmus

4 Paralleler ADI-Algorithmus

5 Zusammenfassung

Paralleler Algorithmus

Aufteilung der Knoten zu den 4 Prozessoren in Streifen

K =

KE KEI

KIE KI

E für Randknoten und I für innere Knoten

K-Matrix für unseren Fall von 4 Prozessoren





























KE01 0 0 0 0 KE01I1 0 0 0

0 KE12 0 0 0 KE12I1 KE12I2 0 0 0 0 K_E23 0 0 0 K_E23I2 K_E23I3 0

0 0 0 KE34 0 0 0 KE34I3 KE34I4

0 0 0 0 KE40 0 0 0 KE40I4

K_I1E01 K_I1E12 0 0 0 K_I1 0 0 0 0 KI2E12 KI2E23 0 0 0 KI2 0 0 0 0 KI3E23 KI3E34 0 0 0 KI3 0

0 0 0 K_I4E34 K_I4E40 0 0 0 K_I4





























Parallleler Algorithmus

mit Richtungsunterscheidung

neue Aufteilung inK=V+H V=

K_E,y K_IE KEI KI,y

und H=

K_E,x 0 0 KI,x

K_E,y Diagonalmatrix

KI,y,KE,x,KI,x parallel invertierbar erstes Ansatz für den Algorithmus

Erster Algorithmus(einfache Übersetzung)

Deklarationsteil

Chooseu⁰ r:=

fE − KE ·u_E⁰ − PP

i=1A^T_i KEI,i ·u_l⁰_,i f_I − K_IE ·u_E⁰ − K_IE ·u_E⁰

σ:=σ₀ := (r,R⁻¹r) k :=0

Erster Algorithmus(einfache Übersetzung)

Hauptteil

whileσ > tol²·σ₀ do T_x :=H+ρ_k+1I u^k∗ :=

T⁻¹_E,x[fE + ρ_k+1u_E^k − K_E,y ·u_E^k − PP

i=1A^T_i K_EI,i ·u_l,i^k T⁻¹_I,x[fI + ρk+1u_I^k − KIE ·u_E^k − KI,y ·u_I^k

T_y :=V+ρ_k+1I g:=

fE − ρk+1u_E^k∗ − KE,x·u_E^k∗

fI − ρ_k+1u_I^k∗ − K_I,x·u_I^k∗

Solve T_yu^k+1 =g r:=

fE − KE ·u_E^k+1 − PP

i=1A^T_i KEI,i ·u^k_l,i⁺¹ fI − K_IE ·u_E^k+1 − K_IE ·u_E^k+1

σ:=σ0 := (r,R⁻¹r) k :=k+1

end

Verbesserter Algorithmus

Deklarationsteil

Versuch eines verbesserten Algorithmus Reduzierung der Kommunikationsschritte

Chooseu⁰ v:=

fE − PP

i=1A^T_i KEI,i ·u_l,i⁰ f_I − K_IE ·u_E⁰

r:=v−K

K_E 0 0 K_I

·u⁰ σ:=σ₀ := (r,R⁻¹r)

Verbesserter Algorithmus

whileσ > tol²·σ₀ do T_x :=H+ρ_k+1I u^k∗ :=T⁻¹_x ·

v+ρ_k+1u^k −

KE,y 0 0 K_I,y

·u^k

Ty :=V+ρk+1I g:=f +ρ_k+1u^k∗−

K_E,x 0 0 K_I,x

·u^k∗

Solve T_yu^k+1 =g v:=

fE − PP

i=1A^T_i KEI,i ·u_l,i^k+1 f_I − K_IE ·u_E^k+1

r:=v−K

K_E 0 0 KI

·u^k+1 σ:=σ₀ := (r,R⁻¹r)

k :=k+1 end

Paralleler Algorithmus

Einführung eines Hilfsvektors v

Somit nur ein Kommunikationsschritt pro Iteration Inversion T⁻¹_x mit einem sequentiellen Invertierer

Einsetzen eines parallelen Algorithmus für das zu lösende GLS (parallel Gauß-Seidel)

Übersicht

1 Einführung

2 ADI-Verfahren

3 Sequentieller ADI-Algorithmus

4 Paralleler ADI-Algorithmus

5 Zusammenfassung

Zusammenfassung

sehr gutes Verfahren für parabolische Probleme in zwei Dimensionen unbedingt stabil

nur ein Kommunikationsschritt pro Iteration anwendbar auf Schrödinger-Probleme aus der QM (imaginärzeitliche Diffusionsprobleme)