Idea metody Dirichleta-Neumanna

Metody dekompozycji obszaru (ang. Domain Decomposition Method), patrz np. [30], [32], [33], są sposobem równoległego rozwiązywania zadań powsta-łych w wyniku dyskretyzacji zagadnień różniczkowych cząstkowych. Ich idea opiera się na rozbiciu zagadnienia różniczkowego cząstkowego określonego na pewnym obszarze Ω na możliwie słabo ze sobą powiązane zagadnienia mniej-szych wymiarów, a w związku z tym łatwiejmniej-szych do rozwiązania. Każde takie zagadnienie odpowiada niewielkiej – na ogół jednemu – liczbie podobszarów Ωi, na jakie podzielony jest wyjściowy obszar, przy czym podobszary mogą się przecinać niepusto. Jeżeli części wspólne są podzbiorami brzegów pod-obszarów (wierzchołki, krawędzie, ściany), to mówimy o dekompozycji „bez zakładek” (patrz [33, rozdz. 1 i 4]), a metody oparte na takiej dekompo-zycji nazywamy strukturalnymi (ang. substructuring methods). W przeciw-nym razie mamy do czynienia z dekompozycją „z zakładkami” (patrz [33, rozdz. 1 i 3]), a metody dekompozycji nazywamy niestrukturalnymi. W nie-których książkach (porównaj np. [32, str. 101]) metody strukturalne dzielone są na iteracyjne i bezpośrednie. Pod nazwą bezpośrednich metod struktural-nych lub po prostu metod strukturalstruktural-nych współcześnie rozumie się algorytmy bezpośredniego rozwiązywania układów równań bądź oparte na faktoryzacji macierzy Schura. W odróżnieniu od nich iteracyjne metody strukturalne (ang.

iterative substructuring methods) inaczej zwane metodami uzupełnień Schura (ang. Schur complement methods) związane są z podejściem iteracyjnym.

Jedna z pierwszych prac poświęcona tym metodom, to [28]. Klasa tych metod obejmuje w szczególności metody przekształcające wyjściowe zagadnienie do równoważnego liniowego układu równań algebraicznych mniejszego wymiaru lub prostszego do rozwiązania oraz rozwiązanie go jedną z metod opartych na przestrzeniach Kryłowa (patrz [31]). Warto dodać, iż dekompozycja może być wprowadzona już na poziomie zagadnień różniczkowych. Często jednak jest wprowadzana dopiero na etapie rozwiązywania układów równań algebraicz-nych powstałych z aproksymacji rozpatrywaalgebraicz-nych dyskretyzacji. W praktyce spotyka się różne kombinacje powyższych podejść.

Zgodnie z koncepcją dekompozycji wprowadźmy na obszarze Ω triangula-cję T^H(Ω), konforemną i quasi-jednostajną, na N rozłącznych, wielokątnych w R² a wielościennych w R³ podobszarów Ω_i, i = 1, 2, . . . , N. Nazwijmy ją grubą siatką w odróżnieniu od drobnej siatki T^h(Ω). Zakładamy, że punkty grubej siatki stanowią podzbiór punktów siatki drobnej (tak zwane siatki zgodne). Oznaczmy maksymalną ze średnic podobszarów przez H.

Oczywi-ście na ogół H ≫ h. Warto dodać, iż w praktyce wygodniej jest wprowadzić najpierw grubą siatkę, a następnie dokonać dalszego podziału poszczególnych podobszarów, spełniającego warunki triangulacji przedstawione w poprzed-nim podrozdziale, w celu uzyskania drobnej siatki.

Wykorzystanie metod dekompozycji obszaru w większości przypadków sprowadza się do zastosowania procesu iteracyjnego (dla zagadnień syme-trycznych i dodatnio określonych na ogół metody gradientów sprzężonych) opartego o niezależne obliczenia na podobszarach. W przypadku zastosowa-nia komputerów równoległych do obliczeń opartych na metodach dekompo-zycji obszaru, każdemu procesorowi przypisujemy jedno lub więcej zagadnień lokalnych i wykonujemy na nich obliczenia niezależnie od siebie. Im więcej dostępnych procesorów, tym lepsze możliwości zrównoleglenia i tym lepsza wydajność całej metody. W idealnym przypadku:

• metoda dekompozycji obszaru powinna być optymalna (ang. optimal ), to znaczy szybkość zbieżności procesu iteracyjnego zastosowanego do znalezienia dobrego przybliżenia szukanego rozwiązania ma być nieza-leżna od wymiaru zadania, czyli od parametru triangulacji h (patrz [33, str. 9]),

• szybkość zbieżności procesu iteracyjnego zastosowanego do znalezienia dobrego przybliżenia szukanego rozwiązania powinna być niezależna od liczby podobszarów, na jakie podzielony został wyjściowy obszar Ω, a tym samym niezależna od parametru H; ang. scalability (porównaj def. 1.3 w [33]).

Niestety wymagania takie są na ogół zbyt rygorystyczne. Zamiast tego wpro-wadza się określenie metody prawie optymalnej (ang. quasi-optimal, zob. [33, str. 18]), to znaczy takiej, której wskaźnik uwarunkowania rośnie jak potęga logarytmu wraz ze wzrostem H/h. Potęga przy logarytmie powinna być oczy-wiście jak najmniejsza – w praktyce równa jeden lub dwa.

Spośród wielu iteracyjnych metod strukturalnych do najbardziej znanych należy z pewnością zaliczyć metody Neumanna-Neumanna opisane między innymi w [33, str. 10], metody Dirichleta-Dirichleta znane również pod na-zwą metod FETI w tej samej monografii na stronach 12-15 oraz metody Neumanna-Dirichleta, które zostały opisane w [32] na stronach 112-117 (przy-padek dwóch podobszarów) oraz stronach 124-125 (przy(przy-padek wielu podob-szarów). W pracy tej zajmować się będziemy metodami Dirichleta-Neumanna iteracyjnego rozwiązywania dyskretyzacji równań różniczkowych cząstko-wych. Były one rozważane już od połowy lat osiemdziesiątych ubiegłego stulecia. Do jednych z najwcześniejszych z pewnością należy zaliczyć prace [13] oraz [12]. Metody te opierają się na takim podziale wyjściowego obszaru

Ω, który w dwóch wymiarach wizualnie przypomina wzór szachownicy. Przy tym podziale każda iteracja w procesie iteracyjnym szukania rozwiązania li-niowego układu równań algebraicznych (1.10) sprowadza się do niezależnych obliczeń na polach (obszarach) koloru białego (typu Dirichleta) i czarnego (typu Neumanna). Dla dowolnego obszaru wprowadzenie takiego podziału może nie być możliwe. Problem ten zostanie dokładniej omówiony w dal-szej części pracy. Ponieważ tematem pracy są metody Dirichleta-Neumanna, celowe wydaje się przedstawienie najpierw idei stojącej za tymi metodami.

Ideę przedstawimy na przykładzie dwóch podobszarów w dwóch wymiarach (patrz rys. 1.5) dla zagadnienia modelowego, to jest:

 −∆u = f w Ω

u = 0 na ∂Ω . (1.11)

Przez Γ oznaczyliśmy tutaj szkielet (ang. interface):

Ω Ω

Γ

D N

n

n

Rysunek 1.5: Dekompozycja na dwa rozłączne podobszary.

Γ = (ΩD∪ ΩN) \ ∂Ω, (1.12)

czyli zbiór wszystkich punktów leżących na brzegach podobszarów, ale nie leżących na brzegu obszaru Ω. nD i nN natomiast oznaczają pochodne nor-malne zewnętrznie do podobszaru odpowiednio ΩD i ΩN. Podany w tym podrozdziale sposób przedstawienia metody jest wzorowany na opisie algo-rytmu dla dwóch podobszarów z [33, roz. 1.3.3]. Uwagi na temat uogólnienia

na większą liczbę podobszarów znajdują się w tej samej książce na stro-nach 18-20.

W ujęciu zagadnień różniczkowych n-ty krok procesu iteracyjnego metody Dirichleta-Neumanna możemy zapisać następująco (patrz [33, str. 8]):

Iteracja n-ta na podobszarze typu Dirichleta:

(D)

Iteracja n-ta na podobszarze typu Neumanna: (1.13)

(N ) τ jest parametrem z przedziału [0; 1].

W podrozdziale 1.2 przedstawione zostały etapy dyskretyzacji metodą MES prowadzące do zagadnienia dyskretnego (1.7) oraz macierzowego (1.10).

Po pogrupowaniu niewiadomych u na należące do szkieletu uΓ i do wnętrz podobszarów uD oraz uN możemy zapisać liniowy układ równań algebraicz-nych z tak zwaną macierzą sztywności A w postaci:



Podzielmy dodatkowo AΓΓ oraz fΓ na części związane z odpowiednimi pod-obszarami:

A_ΓΓ= A^D_ΓΓ+ A^N_ΓΓ, f_Γ = f_Γ^D+ f_Γ^N. (1.15) Stosując wprowadzone oznaczenia możemy zapisać przybliżenie λ^n+1/2_Γ po-chodnej normalnej ^∂u_∂n^n+1/2D

D jako (porównaj [33, str. 3]):

λ^n+1/2_Γ = AΓDu^n+1/2_D + A^D_ΓΓuⁿ_Γ− f_Γ^D. (1.16) W ten sposób n-ta iteracja metody Dirichleta-Neumanna może zostać

przed-stawiona w postaci macierzowej:

Iteracja n-ta na podobszarze typu Dirichleta:

(D) ADDu^n+1/2_D + ADΓuⁿ_Γ = fD

Iteracja n-ta na podobszarze typu Neumanna: (1.17) (N )

gdzie, jak w (1.14) i (1.15), fii f_Γⁱ oznaczają obcięcie wektora prawej strony f do odpowiednio wnętrza oraz części szkieletu związanego z podobszarem Ω_i dla i ∈ {D, N}.

Postać metody przedstawiona powyżej jest niestety związana ze wszyst-kimi niewiadomymi. Standardową praktyką przy projektowaniu algorytmów opartych na metodach dekompozycji obszaru jest eliminacja niewiadomych z wnętrz podobszarów. Operacja ta doprowadzi nas do otrzymania zwię-złej formuły na macierz preconditionera. Z równania (D) wyliczamy zatem zmienne odpowiadające niewiadomym z wnętrza podobszaru typu Dirichleta:

u^n+1/2_D = A⁻_DD¹ f_D− ADΓuⁿ_Γ

. (1.18)

Ze względu na postać zagadnienia wyjściowego (1.11) oraz sposób dyskretyza-cji mamy gwarancję, że macierz ADDjest nieosobliwa. Wstawiając wyliczone zmienne do wzoru na przybliżenie pochodnej normalnej (1.16) otrzymujemy

λ^n+1/2_Γ = AΓDA⁻_DD¹ f_D− A_DΓuⁿ_Γ

+ A^D_ΓΓuⁿ_Γ− f_Γ^D, (1.19) co można zapisać w prostszej postaci:

λ^n+1/2_Γ = − g_Γ^D− S_Duⁿ_Γ

, (1.20)

gdzie g^D_Γ i SD definiujemy (porównaj uzupełnienie Schura w [33, str. 5]) jako:

S_D = A^D_ΓΓ− AΓDA_DD⁻¹A_DΓ,

g_Γ^D = f_Γ^D− A_ΓDA_DD⁻¹f_D. (1.21) Stosując blokową eliminację Gaussa możemy wyeliminować niewiadome również z wnętrz podobszarów typu Neumanna otrzymując:

 A_NN A_NΓ

gdzie:

S_N = A^N_ΓΓ− A_ΓNA⁻_NN¹ A_NΓ, (1.23) g_Γ^N = f_Γ^N− AΓNA⁻_NN¹ f_N. (1.24) W ten sposób sprowadziliśmy zagadnienie (1.17) do rozwiązania w każdej iteracji następującego układu równań:

S_Nu^n+1/2_Γ = g_Γ^N− λ^n+1/2_Γ (1.25) oraz wyliczenia nowego przybliżenia uⁿ⁺¹_Γ z równania (Γ) w (1.17). Otrzy-mujemy w ten sposób:

S_N uⁿ⁺¹_Γ − (1 − τ )uⁿ_Γ

= τ 

g_Γ^N− λ^n+1/2_Γ 

. (1.26)

Podstawiając wyliczone wyrażenie na λ^n+1/2_Γ dostajemy:

S_N uⁿ⁺¹_Γ − (1 − τ )uⁿ_Γ

= τ g^N_Γ + g^D_Γ − S_Duⁿ_Γ

, (1.27) a po odjęciu od obu stron wyrazu τ S_Nuⁿ_Γ:

S_N uⁿ⁺¹_Γ − uⁿ_Γ

= τ gΓ− Suⁿ_Γ

, (1.28)

gdzie

g_Γ = g_Γ^D+ g_Γ^N, S = SD+ SN. (1.29) Pokazaliśmy tym samym, że algorytm Dirichleta-Neumanna jest po prostu iteracją Richardsona (patrz dodatek A) z parametrem τ dla układu równań algebraicznych

Su_Γ = gΓ (1.30)

z preconditionerem SN. W n-tej iteracji procesu iteracyjnego Richardsona musimy zatem:

• wykonać mnożenie wektora uⁿ_Γ przez macierz S,

• rozwiązać układ równań algebraicznych z macierzą S_N.

Mnożenie przez macierz S można sprowadzić do rozwiązania niezależnych zagadnień na poszczególnych podobszarach:

Suⁿ_Γ = S_Nuⁿ_Γ+ S_Duⁿ_Γ. (1.31)

Z definicji macierzy Sidla i ∈ {D, N} każde takie lokalne zagadnienie spro-wadza się do wykonania mnożeń przez macierze A_iΓ, A_Γi, Aⁱ_ΓΓ oraz rozwią-zania zagadnienia Dirichleta na podobszarze Ωi:

 A_iiv_iⁿ= AiΓuⁿ_Γ,

S_iuⁿ_Γ = Aⁱ_ΓΓuⁿ_Γ− A_Γiv_iⁿ. (1.32) Wszystkie te lokalne macierze są rozrzedzone, co ogranicza koszt pamięciowy ich przechowywania. Do rozwiązania zagadnienia Dirichleta zaś wystarczy za-pamiętać czynniki rozkładu trójkątno-trójkątnego⁵ macierzy A_ii. Natomiast rozwiązanie układu równań z macierzą SN i prawą stroną bΓ:

S_Nv_Γ = bΓ

jest równoważne rozwiązaniu zagadnienia Neumanna na podobszarze Ω_N:

 A_NN A_NΓ A_ΓN A^N_ΓΓ

  v_N v_Γ



=

 0 b_Γ



. (1.33)

Każdą iterację zatem można sprowadzić do równoległych obliczeń wykonywa-nych lokalnie na poszczególwykonywa-nych podobszarach. Warto podkreślić, że metoda jest dla dwóch podobszarów optymalna w sensie definicji podanej w tym podrozdziale, to jest szybkość zbieżności procesu iteracyjnego do dokładnego rozwiązania jest niezależna od wymiaru zadania (porównaj [33, str. 19-20]).

Wskaźnik uwarunkowania układu ze skonstruowanym powyżej preconditio-nerem:

S_N⁻¹Su_Γ = S_N⁻¹g_Γ (1.34) spełnia zależność (patrz [33, str. 10]):

cond S_N⁻¹S

≤ C, (1.35)

gdzie stała C nie zależy od wymiaru zadania (1.14). W istocie dla przypadku symetrycznego podziału na dwa podobszary, do rozwiązania wystarczy tylko jedna iteracja i w związku z tym stała C = 1. W dalszej części pracy zostanie pokazane, iż w przypadku uogólnienia algorytmu na większą liczbę podob-szarów dostajemy prawie optymalność. Udowodnimy bowiem, iż w ogólnym przypadku wskaźnik uwarunkowania otrzymanego układu rośnie jak druga potęga logarytmu wraz ze wzrostem H/h.

5Ze względu na symetryczność macierzy Aii zamiast rozkładu LU można dokonać rozkładu Choleskiego LL^T lub Banachiewicza-Choleskiego LDL^T (patrz [14, str. 50])

1.4 Abstrakcyjna teoria addytywnych metod