Tìm nghiệm của phương trình Poisson ba chiều bằng hai phương pháp TFQMR và Gmres(m)

TÌM NGHIỆM CỦA PHƯƠNG TRÌNH POISSON BA CHIỀU<br /> BẰNG HAI PHƯƠNG PHÁP TFQMR VÀ GMRES(m)<br /> ĐÀO HỮU HÀ<br /> Trường Phổ thông Dân tộc Nội trú Tỉnh Kon Tum<br /> ĐINH NHƯ THẢO<br /> Trường Đại học Sư phạm – Đại học Huế<br /> Tóm tắt: Bài báo trình bày việc xây dựng chương trình giải phương trình<br /> Poisson ba chiều dựa trên hai thuật toán TFQMR và GMRES(m) để sử dụng<br /> trong chương trình mô phỏng linh kiện na-nô bán dẫn bằng phương pháp<br /> Monte – Carlo tập hợp tự hợp. Để kiểm tra hiệu năng, các chương trình mô<br /> phỏng tương ứng được áp dụng để mô phỏng các đi-ốt p-i-n bán dẫn GaAs.<br /> Các kết quả chỉ ra rằng chương trình giải phương trình Poisson dựa trên<br /> thuật toán TFQMR có tốc độ hội tụ nhanh hơn nhiều so với chương trình sử<br /> dụng thuật toán GMRES(m). Cả hai thuật toán chạy chậm hơn so với thuật<br /> toán BICGSTAB(3) nhưng bù lại có tính ổn định cao hơn nhiều.<br /> <br /> 1. GIỚI THIỆU<br /> Nghiên cứu và phát triển các linh kiện na-nô bán dẫn đang thu hút sự quan tâm mạnh<br /> mẽ của giới khoa học do tính ứng dụng cao của nó [1], [2]. Nghiên cứu thực nghiệm các<br /> linh kiện na-nô nói chung là rất tốn kém, đòi hỏi phải sử dụng công nghệ cao và mất<br /> nhiều thời gian. Các phương pháp nghiên cứu lý thuyết có thể giúp khắc phục được các<br /> hạn chế nêu trên [3], đặc biệt phương pháp mô phỏng Monte – Carlo tập hợp tự hợp với<br /> các ưu điểm nổi trội là tính chính xác và tính ổn định.<br /> Trong quá trình mô phỏng, phương pháp Monte – Carlo tập hợp tự hợp cần cập nhật<br /> phân bố của điện thế trong linh kiện thông qua việc giải phương trình Poisson, thông<br /> thường bằng phương pháp sai phân hữu hạn [3]. Khi đó việc giải phương trình Poisson<br /> chuyển thành việc giải một hệ phương trình tuyến tính thưa cực lớn với hàng triệu<br /> phương trình và hàng triệu ẩn. Thông thường để giải hệ phương trình trên người ta phải<br /> sử dụng các phương pháp số chạy trên một siêu máy tính với bộ nhớ cực lớn mà Việt<br /> Nam hiện nay chưa có. May mắn là các phương pháp không gian con Krylov có thể hỗ<br /> trợ cách tính toán không cần lưu trữ các số liệu tính toán trung gian [4], [5], [6]. Một số<br /> tác giả đã sử dụng các phương pháp BICGSTAB, BICGSTAB tiền điều kiện,<br /> BICGSTAB2, BICGSTAB(3) và GPBICG để giải phương trình Poisson và đã thu được<br /> các kết quả chính xác với thời gian tính toán được rút ngắn nhiều lần [7], [8], [9], [10].<br /> Đó là động lực để chúng tôi tiến hành tìm nghiệm của phương trình Poisson bằng hai<br /> phương pháp TFQMR và GMRES(m) [5] với mục đích tìm ra những phương pháp tối<br /> ưu, hoạt động ổn định hơn và cho kết quả nhanh hơn.<br /> Chúng tôi đã xây dựng hai chương trình mô phỏng mới và thực hiện tính toán trên máy<br /> tính Dell Inspiron 14R-N4010 (Intel(R) Core(TM) i3 CPU M370 @ 2.4GHz DDR 4GB).<br /> Tạp chí Khoa học và Giáo dục, Trường Đại học Sư phạm Huế<br /> ISSN 1859-1612, Số 04(20)/2011: tr. 5-12<br /> <br /> 6<br /> <br /> ĐÀO HỮU HÀ - ĐINH NHƯ THẢO<br /> <br /> Kết quả chỉ ra rằng chương trình mô phỏng sử dụng thuật toán TFQMR có nhiều ưu điểm<br /> còn chương trình mô phỏng sử dụng thuật toán GMRES(m) thực tế không hiệu quả.<br /> 2. GIẢI PHƯƠNG TRÌNH POISSON BA CHIỀU BẰNG HAI THUẬT TOÁN<br /> TFQMR VÀ GMRES(m)<br /> Giả sử vật liệu là đồng nhất thì phương trình Poisson trong trường hợp ba chiều có dạng<br /> <br /> ∂ 2ϕ ∂ 2ϕ ∂ 2ϕ<br /> ρ<br /> + 2 + 2 =− ,<br /> 2<br /> ∂x<br /> ∂y<br /> ∂z<br /> εS<br /> <br /> (1)<br /> <br /> ở đây ϕ là điện thế, ρ là mật độ điện tích, ε S là hằng số điện môi tĩnh trong linh kiện;<br /> x , y , z là ba biến không gian. Để có thể dễ dàng thực hiện sai phân hữu hạn ta chia<br /> mô hình linh kiện thành các ô lưới và giả sử khoảng cách giữa các nút lưới theo các<br /> chiều không gian là bằng nhau, Δx = Δy = Δz . Tiến hành lấy sai phân hữu hạn phương<br /> trình (1) ta thu được hệ phương trình sau.<br /> <br /> ϕi −1, j ,k + ϕi , j −1,k + ϕi , j ,k −1 − 6ϕi , j ,k + ϕi +1, j ,k + ϕi , j +1,k + ϕi , j ,k +1 = −<br /> <br /> ρi , j , k 2<br /> Δx , (2)<br /> εS<br /> <br /> ở đây i = 1, N x , j = 1, N y , k = 1, N z với N x , N y , N z lần lượt là số nút lưới theo các<br /> chiều không gian Ox , Oy , Oz . Đây chính là hệ phương trình tuyến tính thưa cực lớn<br /> mà ta cần giải. Hai giải thuật TFQMR và GMRES(m) [5] để tìm nghiệm của phương<br /> trình Poisson được khai triển trong Bảng 1.<br /> Bảng 1. Thuật toán TFQMR và GMRES(m) để tìm nghiệm của phương trình Poisson.<br /> 1<br /> <br /> Thuật toán GMRES(m)<br /> Chọn ϕ0 ban đầu<br /> <br /> 1<br /> <br /> Thuật toán TFQMR<br /> Chọn ϕ0 ban đầu<br /> <br /> 2<br /> <br /> r0 = b – Aϕ0, β = ||r0||2, v1 = r0/β,<br /> <br /> 2<br /> <br /> w0 = u0 = r0 = b – Aϕ0, v0 = Au0, d0 = 0<br /> <br /> 3<br /> <br /> Do j = 1, 2,…, m<br /> <br /> 3<br /> <br /> τ0 = ||r0||2, θ0 = η0 = 0<br /> <br /> 4<br /> <br /> wj = Avj<br /> <br /> 4<br /> <br /> Chọn rg ≠ 0 bất kỳ, ρ0 = (rg,r0);<br /> <br /> 5<br /> <br /> Doi: Do I = 1,…, j<br /> <br /> 5<br /> <br /> Do m = 0, 1, 2,…<br /> <br /> 6<br /> <br /> hi,j = (vi,wj), wj = wj – hi,jvi<br /> <br /> 6<br /> <br /> If m chẵn thì tính<br /> <br /> 7<br /> <br /> Enddo Doi<br /> <br /> 7<br /> <br /> αm+1 = αm = ρm/(vm,rg)<br /> <br /> 8<br /> <br /> hj+1,j = ||wj||2<br /> <br /> 8<br /> <br /> um+1 = um - αmvm<br /> <br /> 9<br /> <br /> vj+1 = wj/hj+1,j<br /> <br /> 9<br /> <br /> Endif<br /> <br /> 10<br /> <br /> Enddo<br /> <br /> 10<br /> <br /> 11<br /> <br /> ym (là tối thiểu của ||βe1 – Hmy||2)<br /> <br /> 11<br /> <br /> wm+1 = wm - αmAum<br /> dm+1 = um + ( /αm)ηmdm;<br /> <br /> 12<br /> <br /> ϕ m = ϕ 0 + v my m<br /> <br /> 12<br /> <br /> 13<br /> <br /> If ϕm thoả mãn thì stop<br /> <br /> 13<br /> <br /> θm+1 = ||wm+1||/τm; cm+1 = 1/(1+θm+1)1/2<br /> τm+1 = τmθm+1cm+1; ηm+1 =<br /> α m;<br /> <br /> TÌM NGHIỆM CỦA PHƯƠNG TRÌNH POISSON BA CHIỀU BẰNG HAI PHƯƠNG PHÁP...<br /> <br /> 14<br /> <br /> Else<br /> <br /> 14<br /> <br /> ϕm+1 = ϕm ηm+1dm+1;<br /> <br /> 15<br /> <br /> ϕ0 = ϕm<br /> <br /> 15<br /> <br /> If m lẻ thì tính:<br /> <br /> 16<br /> <br /> Goto 2<br /> <br /> 16<br /> <br /> ρm+1 = (wm+1,rg); βm-1 = ρm+1/ρm-1<br /> <br /> 17<br /> <br /> um+1 = wm+1 + βm-1um<br /> <br /> 18<br /> <br /> vm+1 = Aum+1 + βm-1(Aum + βm-1vm-1)<br /> <br /> 19<br /> <br /> Endif<br /> <br /> 20<br /> <br /> Enddo<br /> <br /> 7<br /> <br /> 3. KẾT QUẢ MÔ PHỎNG VÀ THẢO LUẬN<br /> Để kiểm tra hiệu năng của chương trình giải phương trình Poisson ba chiều chúng tôi<br /> tích hợp chương trình này vào chương trình mô phỏng bằng phương pháp Monte –<br /> Carlo tập hợp tự hợp để mô phỏng động lực học của hạt tải trong các đi-ốt p-i-n bán dẫn<br /> GaAs. Mô hình cấu trúc của đi-ốt p-i-n bán dẫn GaAs gồm một lớp bán dẫn thuần (i)<br /> kẹp giữa hai lớp bán dẫn pha tạp loại p và loại n như được chỉ ra trong Hình 1, trong đó<br /> mỗi lớp có độ dày tương ứng là di = 340 nm , d p = d n = 50 nm . Mật độ pha tạp acceptor<br /> và donor tương ứng là N A = 0.5 ×1017 cm −3 và N D = 2.5 × 1017 cm −3 . Các hạt tải kích<br /> thích quang được tạo ra trong linh kiện bằng cách chiếu một xung laser với chiều dài<br /> xung là 12 fs và năng lượng photon là 1.49 eV , mật độ hạt tải quang là<br /> N ex = 5 × 1016 cm −3 sau thời gian 1 ps . Kích thước theo ba chiều không gian của đi-ốt là<br /> Lx × Ly × Lz = 440 nm ×100 nm ×100 nm , giả sử đi-ốt được nuôi cấy theo phương Ox .<br /> Linh kiện được chia thành các ô lưới không gian với Δx = Δy = Δz = 50 ×10−10 m . Như<br /> vậy ta sẽ có N x = 89 nút<br /> lưới theo phương Ox ,<br /> N y = 21 nút lưới theo<br /> phương Oy và N z = 21 nút<br /> lưới theo phương Oz . Điện<br /> trường ngoài được đặt vào<br /> linh kiện dọc theo phương<br /> Ox và đi-ốt được phân cực<br /> nghịch, xem Hình 1.<br /> Hình 1. Mô hình đi-ốt p-i-n GaAs<br /> Hình 2 mô tả sự thay đổi<br /> vận tốc trôi dạt của điện tử theo các phương Ox , Oy và Oz và vận tốc trôi dạt toàn<br /> phần ứng với điện trường ngoài Eext = 100 kV cm , được tính toán với hai giải thuật giải<br /> phương trình Poisson khác nhau: Hình 2a) là kết quả tính toán với giải thuật<br /> GMRES(m) còn Hình 2b) là kết quả tính toán với giải thuật TFQMR. Cả hai giải thuật<br /> đều cho những kết quả tương tự nhau.<br /> <br /> 8<br /> <br /> ĐÀO HỮU HÀ - ĐINH NHƯ THẢO<br /> <br /> Hình 3 mô tả sự phụ thuộc của vận tốc trôi dạt của điện tử theo thời gian ứng với các giá<br /> trị điện trường ngoài Eext = 100 kV cm và 150kV cm , cũng được tính toán với hai giải<br /> thuật giải phương trình Poisson khác nhau và hai giải thuật cho những kết quả gần như<br /> trùng khớp.<br /> <br /> Hình 2. Vận tốc trôi dạt của điện tử theo các phương khác nhau và vận tốc trôi dạt toàn phần<br /> như là hàm của thời gian ứng với Eext = 100 kV cm được tính toán với hai giải thuật giải<br /> phương trình Poisson khác nhau: a) GMRES(m), b) TFQMR<br /> <br /> Hình 3. Vận tốc trôi dạt toàn phần của điện tử<br /> như là hàm của thời gian ứng với các điện<br /> trường ngoài khác nhau và được tính toán với<br /> hai giải thuật khác nhau<br /> <br /> Hình 4. Vận tốc trôi dạt của điện tử theo<br /> phương Ox như là hàm của thời gian thu<br /> được bằng ba chương trình mô phỏng khác<br /> nhau ứng với Eext = 100 kV cm<br /> <br /> Hình 4 cho kết quả so sánh vận tốc của điện tử thu được bằng ba chương trình mô<br /> phỏng ba chiều sử dụng thuật toán GMRES(m), thuật toán TFQMR và thuật toán<br /> BICGSTAB(3) [8], ba đồ thị gần như trùng nhau hoàn toàn, đặc biệt chương trình dùng<br /> thuật toán TFQMR cho kết quả trùng hoàn toàn với kết quả thu được khi sử dụng thuật<br /> <br /> TÌM NGHIỆM CỦA PHƯƠNG TRÌNH POISSON BA CHIỀU BẰNG HAI PHƯƠNG PHÁP...<br /> <br /> 9<br /> <br /> toán BICGSTAB(3). Điều đó cho thấy phương pháp TFQMR hoạt động hiệu quả hơn<br /> phương pháp GMRES(m) do chương trình ba chiều sử dụng thuật toán BICGSTAB(3)<br /> đã được chứng minh là cho kết quả phù hợp với thực nghiệm [8].<br /> Hình 5 mô tả sự phân bố điện thế không gian trong đi-ốt p-i-n bán dẫn GaAs theo hai<br /> phương Ox , Oy tại mặt cắt z = 50 nm ứng với điện trường ngoài Eext = 100 kV cm được<br /> tính toán bằng ba chương trình mô phỏng ba chiều sử dụng ba thuật toán GMRES(m),<br /> TFQMR và BICGSTAB(3). Ba đồ thị trong hình 5 đều gần như trùng nhau hoàn toàn.<br /> Kết quả này cũng phù hợp với các kết quả đã được công bố trước đây [1], [8].<br /> <br /> Hình 5. Phân bố điện thế không gian trong đi-ốt<br /> p-i-n bán dẫn GaAs theo hai phương Ox, Oy tại<br /> mặt cắt z = 50 nm ứng với Eext = 100 kV cm ,<br /> được tính toán bằng ba chương trình mô phỏng<br /> khác nhau<br /> <br /> Hình 6. Sự phụ thuộc của chuẩn Euclid của<br /> vectơ thặng dư vào số vòng lặp của chương<br /> trình con Poisson ứng với Eext = 100 kV cm<br /> <br /> Để so sánh tốc độ hội tụ và tính ổn định của hai thuật toán TFQMR và GMRES(m) với<br /> các đặc trưng tương ứng của thuật toán BICGSTAB(3), chúng tôi đã tiến hành khảo sát<br /> sự phụ thuộc của chuẩn Euclid của vectơ thặng dư vào số vòng lặp của chương trình con<br /> Poisson, Hình 6. Chuẩn Euclid của vectơ thặng dư được tính theo công thức [4].<br /> <br /> r 2 = rT r ,<br /> <br /> (7)<br /> <br /> với r = b − Aϕ là vectơ thặng dư. Từ đồ thị ta thấy rằng, cả hai thuật toán TFQMR và<br /> GMRES(m) đều cần số vòng lặp lớn hơn rất nhiều để tìm ra nghiệm có cùng chuẩn<br /> Euclid của vector thặng dư so với nghiệm tìm được bằng thuật toán BICGSTAB(3).<br /> Điều này có nghĩa là cả hai thuật toán này có tốc độ hội tụ chậm hơn thuật toán<br /> BICGSTAB(3). Tuy nhiên, đồ thị tương ứng với hai thuật toán TFQMR và GMRES(m)<br /> trơn hơn đồ thị tương ứng với thuật toán BICGSTAB(3), hàm ý rằng hai thuật toán này<br /> cho kết quả ổn định hơn thuật toán BICGSTAB(3).<br /> <br />