Projekt 6: Równanie Poissona - rozwiązanie metodą algebraiczną.

(1)

Projekt 6: Równanie Poissona - rozwiązanie metodą algebraiczną.

Tomasz Chwiej 29 sierpnia 2018

1 Wstęp

1.1 Dyskretyzacja

i= 0 1

n_x 2

0 1 2 n_y

j=

V₁

V₂

V₃

V₄

Rysunek 1: Geometria układu i schemat siatki obliczeniowej (docelowej). Potencjały V₁, V₂, V₃, V₄ określają warunki brzegowe Dirichleta, ρ1, ρ2to gęstości ładunku, ε1, ε2- stałe dielektryczne w obszarze lewym i prawym. Potencjał w 4 narożnych czerwonych węzłach nie wpływa na rozwiązanie i może mieć w nich dowolną wartość.

Na zajęciach rozwiążemy równanie Poissona dla układu pokazanego na Rys.1 postępując następu- jąco: i) zdyskretyzujemy równanie na regularnej siatce przy użyciu ilorazów różnicowych, ii) zapiszemy równanie jako układ równań liniowych z macierzą rzadką, iii) układ równań rozwiążemy stosując me- tody dla macierzy rzadkich.

Ogólne rów. Poissona dla obszaru obejmującego różne wartości stałej dielektrycznej

∇ε∇V = −ρ (1)

zapisujemy jako

ε∇²V +∇ε · ∇V = −ρ (2)

Dyskretyzujemy rów. Poissona stosując ilorazy różnicowe: ^d_dx²^f2 = ^fⁱ⁺¹^−2f_2∆xⁱ^+fⁱ⁻¹ (dokładność O(∆x²)) oraz _dx^df = ^fⁱ⁺¹_∆x^−fⁱ (gorsza dokładność O(∆x) ale skok ε dokonuje się na jednym oczku siatki). Przyj- mujemy oznaczenia

∆x = ∆y = ∆ (3)

xi= ∆· i, i = 0, 1, 2, . . . , nx (4) y_i = ∆· j, j = 0, 1, 2, . . . , ny (5)

V (x, y)→ V (xi, y_j)→ Vi,j (6)

ρ(x, y)→ ρ(xi, yj)→ ρi,j (7)

ε(x, y)→ ε(xi, yj)→ εi,j (8)

(2)

wówczas równanie zdyskretyzowane ma postać εi,j

(Vi+1,j − 2Vi,j+ Vi−1,j

∆² +Vi,j+1− 2Vi,j+ Vi,j−1

∆²

)

+

(εi+1,j− εi,j

∆

Vi+1,j − Vi,j

∆ + εi,j+1− εi,j

∆

Vi,j+1− Vi,j

∆ )

= −ρi,j (9)

1.2 Algebraizacja równania

Równanie (9) zapiszemy w postaci macierzowej. Najpierw dokonamy reindekscaji węzłów, węzły po- numeryjemy tak aby każdemu odpowiadał jeden indeks (l) zamiast dwóch (i,j)

l = i + j· (nx+ 1), (i = 0, 1, . . . , nx; j = 0, 1, 2, . . . , ny) (10) l = 0, 1, 2, . . . , N− 1, N = (nx+ 1)· (ny + 1) (11) przejście w drugą stronę l→ (i, j) (przydatne np. przy sprządzaniu map potencjału)

j = f loor ( l

n_x+ 1 )

(wyznaczamy część całkowitą) (12)

i = l− j · (nx+ 1) (13)

Równanie (9) zapiszemy w postaci macierzowej (algebraicznej)

A⃗V = ⃗b (14)

w której mnożenie l-tego wiersza A przez ⃗V przebiega jak poniżej

a_l,l_−n_x₋₁· Vl−n^x−1+ a_l,l₋₁· Vl−1+ a_l,l· Vl+ a_l,l+1· Vl+1+ a_l,l+n_x₊₁· Vl+nx+1 = b_l (15) a elementy macierzowe są zdefiniowane następująco

a_l,l_−n_x₋₁ = ε_l

∆² (16)

al,l−1 = εl

∆² (17)

al,l = −2ε_l+ ε_l+1+ ε_l+n_x₊₁

∆² (18)

a_l,l+1 = ε_l+1

∆² (19)

a_l,l+n_x₊₁ = ε_l+n_x₊₁

∆² (20)

Macierz układu jest macierzą rzadką, 5-przekątniową (zob. Rys.2).

Zmianę stałej dielektrycznej w przestrzeni określa poniższa relacja ε_l=

{

ε1, i¬ nx/2

ε₂, i > n_x/2 (21)

1.3 Warunki brzegowe Dirichleta

Zastosujemy warunki brzegowe Dirichelta tj. ustalone wartości potencjału w węzłach leżących na brzegu obszaru (rys. 1). Dla węzła leżącego na brzegu (indeks k na rys.2) WB wprowadzamy następująco:

a_k,k_−n_x₋₁= a_k,k₋₁ = a_k,k+1= a_k,k+n_x₊₁= 0 (zerowanie wiersza poza diagonalą) (22)

(3)

V_k = V_B

k

Rysunek 2: Macierz układu jest macierzą rzadką 5-przekątniową.

2 Zadania do wykonania

Macierz układu zapiszemy w formacie CSR (compressed sparse row) indeksując elementy zaczynając od 0. Układ rozwiążemy stosując metodę GMRES używając procedury numerycznej napisanej w języku C. Aby sprawnie rozwiązać problem proszę postępować według poniższych wskazówek

1. Ustalamy początkowe wartości parametrów: ∆ = 0.1, nx = ny = 4, ε1 = ε2 = 1, V1 = V3 = 10, V₂ = V₄ =−10, ρ1(x, y) = ρ_x,y = 0, N = (n_x+ 1)(n_y+ 1)

2. Tworzymy 3 wektory definiujące macierz układu: aaa,iaiaia,jajaja oraz wektor wyrazów wolnych bbb i wektor rozwiązań VVV gdzie:

• aaa to wektor typu double o 5 · N elementach (niezerowe wartości el. macierzowych)

• jajaja to wektor typu int o 5· N elementach (przechowuje informacje o numerach kolumn)

• iaiaia - wektor typu int, ilość elementów N +1 (wskaźniki do elementów rozpoczynających dany wiersz). Uwaga: po utworzeniu wektora iaiaia wszystkie jego elementy wypełniamy wartością

−1.

• bbb i VVV to wektory typu double o N elementach

3. Wyznaczamy niezerowe elementy macierzy A oraz wektora wyrazów wolnych bbb uwzględniając WB Dirichleta zgodnie zalgorytmem przedstawionym w sekcji (3). W celu sprawdzenia popraw- ności wypełnienia macierzy A i wektora bbb należy zapisać je do pliku, podając niezerowe elementy l, i_l, j_l, a[l] (dla macierzy) oraz l, i_l, j_l, b[l] (dla wektora). 20 pkt.

4. Rozwiązujemy układ równań z macierzą rzadką stosując procedurę (dołączyć pliki: mgmres.c i mgmres.h)

vo i d p m g m r e s _ i l u _ c r ( int N , int nz_num , int ia [] , int ja [] , d o u b l e a [] , d o u b l e V [] , d o u b l e b [] , int itr_max , int mr , d o u b l e tol_abs ,

d o u b l e t o l _ r e l )

Opis znaczenia argumentów znaleźć można w plikumgmres.c. Wraz ze zmianą n_x i n_y zmieniają się: N - rozmiar układu (wzór 11) oraz nz num - ilość elementów niezerowych w macierzy A (przedostatnia instrukcja w algorytmie w sekcji 3), dla pozostałych parametrów można przyjąć następujące wartości: itr max = 500, mr = 500, tol abs = 10⁻⁸, tol rel = 10⁻⁸.

5. Sporządzić mapy potencjału dla ε₁ = ε₂ = 1, V₁ = V₃ = −V2 = −V4 = 10, ρ⁽¹⁾ = ρ⁽²⁾ = 0 i trzech siatek: a) nx = ny = 50, b) nx = ny = 100, c) nx= ny = 200. 50 pkt.

6. Ustalamy wartości parametrów: n_x = n_y = 100, V₁ = V₃ = V₂ = V₄ = 0, x_max = ∆ · nx,

(4)

y_max = ∆· ny, σ = x_max/10 oraz ρ⁽¹⁾(x, y) = (+1)· exp

(

−(x− 0.25 · xmax)²

σ² −(y− 0.5 · ymax)² σ²

)

(25)

ρ⁽²⁾(x, y) = (−1) · exp (

−(x− 0.75 · xmax)²

σ² −(y− 0.5 · ymax)² σ²

)

(26) Znaleźć rozkłady potencjału dla 3 przypadków: a) ε₁ = 1 i ε₂ = 1, b) ε₁ = 1 i ε₂ = 2 oraz c) ε1 = 1 i ε2 = 10. Za przedstawienie 3 map potencjału 30 pkt.

3 Algorytm wypełniania macierzy rzadkiej w formacie CSR + WB Dirichleta

Po utworzeniu tablic: aaa, iaiaia, jajaja oraz wektora wyrazów wolnych bbb, ich elementy wypełniamy zgodnie z poniższym algorytmem

i n i c j a l i z a c j a : k = -1 // n u m e r u j e n i e z e r o w e e l e m e n t y A FOR l =0 TO l < N S T E P 1 DO

int b r z e g =0 // w s k a ź n i k p o ł o ż e n i a : 0 - ś r o d e k o b s z a r u ; 1 - b r z e g d o u b l e vb =0. // p o t e n c j a l na b r z e g u

IF ( i = = 0 ) T H E N // l e w y b r z e g b r z e g =1

vb = v1 END IF

IF ( j == ny ) T H E N // g ó r n y b r z e g b r z e g =1

vb = v2 END IF

IF ( i == nx ) T H E N // p r a w y b r z e g b r z e g =1

vb = v3 END IF

IF ( j = = 0 ) T H E N // d o l n y b r z e g b r z e g =1

vb = v4 END IF

// w y p e ł n i a m y od r a z u w e k t o r w y r a z ó w w o l n y c h b [ l ]= -(ρ⁽¹⁾_l +ρ⁽²⁾_l ); // jeśli w środku jest g ę s to ś ć IF ( b r z e g = = 1 ) T H E N

b [ l ]= vb ; // w y m u s z a m y w a r t o ś ć p o t e n c j a ł u na b r z e g u END IF

(5)

ia [ l ]= -1; // w s k a ź n i k dla p i e r w s z e g o el . w w i e r s z u // l e w a s k r a j n a p r z e k a t n a

IF ( l - nx -1 >=0 AND b r z e g ==0 ) T H E N k ++

if ( ia [ l ] <0) ia [ l ]= k a [ k ]=al,l−nx−1

ja [ k ]=l− nx− 1 E N D I F

// p o d d i a g o n a l a

IF ( l -1 >=0 AND b r z e g ==0 ) T H E N k ++

if ( ia [ l ] <0) ia [ l ]= k a [ k ]=a_l,l₋₁

ja [ k ]=l− 1 END IF

// d i a g o n a l a k ++

if ( ia [ l ] <0) ia [ l ]= k IF ( b r z e g = = 0 ) T H E N

a [ k ]=al,l

E L S E a [ k ]=1 END IF ja [ k ]=l // n a d d i a g o n a l a

IF ( l < N AND b r z e g ==0 ) T H E N k ++

a [ k ]=al,l+1

ja [ k ]=l + 1 END IF

// p r a w a s k r a j n a p r z e k ą t n a

IF ( l <N− nx− 1 AND brzeg ==0 ) THEN k ++

a [ k ]=a_l,l+n_x₊₁ ja [ k ]=l + nx+ 1 END IF

END DO

n z _ n u m = k +1 // i l o s c n i e z e r o w y c h e l e m e n t o w (1 e l e m e n t ma i n d e k s 0) ia [ N ]= n z _ n u m

4 Przykładowe rozwiązania

(6)

nx=ny=50

’v50.dat’ u 1:2:3

0 10 20 30 40 50

-10 -5 0 5 10

example

Rysunek 3: nx= ny = 50, ρ⁽¹⁾ = ρ⁽²⁾= 0, ε1 = ε2= 1

nx=ny=100, e1=1, e2=2

’g_e12.dat’ u 1:2:3

0 20 40 60 80 100

-0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8

example

Rysunek 4: n_x = n_y = 100, ρ⁽¹⁾ = ρ⁽²⁾̸= 0, ε1 = 1, ε₂ = 2