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ạp chí Khoa học và Giáo dục, Trường Đại học Sư phạm Huế ISSN 1859-1612, Số 04(20)/2011: tr. 5-12 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) ĐÀO HỮU HÀ Trường Phổ thông Dân tộc Nội trú Tỉnh Kon Tum ĐINH NHƯ THẢO Trường Đại học Sư phạm – Đại học Huế 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 Poisson ba chiều dựa trên hai thuật toán TFQMR và GMRES(m) để sử dụng trong chương trình mô phỏng linh kiện na-nô bán dẫn bằng phương pháp Monte – Carlo tập hợp tự hợp. Để kiểm tra hiệu năng, các chương trình mô phỏng tương ứng được áp dụng để mô phỏng các đi-ốt p-i-n bán dẫn GaAs. Các kết quả chỉ ra rằng chương trình giải phương trình Poisson dựa trên thuật toán TFQMR có tốc độ hội tụ nhanh hơn nhiều so với chương trình sử 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 toán BICGSTAB(3) nhưng bù lại có tính ổn định cao hơn nhiều. 1. GIỚI THIỆU 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 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 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 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 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 các ưu điểm nổi trội là tính chính xác và tính ổn định. 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 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 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 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 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 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 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ỗ 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ố tác giả đã sử dụng các phương pháp BICGSTAB, BICGSTAB tiền điều kiện, BICGSTAB2, BICGSTAB(3) và GPBICG để giải phương trình Poisson và đã thu được 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]. Đó 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 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 ưu, hoạt động ổn định hơn và cho kết quả nhanh hơn. 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 tính Dell Inspiron 14R-N4010 (Intel(R) Core(TM) i3 CPU M370 @ 2.4GHz DDR 4GB). ĐÀO HỮU HÀ - ĐINH NHƯ THẢO 6 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 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ả. 2. GIẢI PHƯƠNG TRÌNH POISSON BA CHIỀU BẰNG HAI THUẬT TOÁN TFQMR VÀ GMRES(m) 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 2 2 2 2 2 2 , Sx y z ϕ ϕ ϕ ρ ε ∂ ∂ ∂ + + = − ∂ ∂ ∂ ở đâ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; 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 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 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 trình (1) ta thu được hệ phương trình sau. , , 2 1, , , 1, , , 1 , , 1, , , 1, , , 16 , i j k i j k i j k i j k i j k i j k i j k i j k S x ρ ϕ ϕ ϕ ϕ ϕ ϕ ϕ ε− − − + + + + + − + + + = − Δ ở đây 1, xi N= , 1, yj N= , 1, zk N= với xN , yN , zN lần lượt là số nút lưới theo các 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 mà ta cần giải. Hai giải thuật TFQMR và GMRES(m) [5] để tìm nghiệm của phương trình Poisson được khai triển trong Bảng 1. Bảng 1. Thuật toán TFQMR và GMRES(m) để tìm nghiệm của phương trình Poisson. Thuật toán GMRES(m) Thuật toán TFQMR 1 Chọn ϕ0 ban đầu 1 Chọn ϕ0 ban đầu 2 r0 = b – Aϕ0, β = ||r0||2, v1 = r0/β, 2 w0 = u0 = r0 = b – Aϕ0, v0 = Au0, d0 = 0 3 Do j = 1, 2,, m 3 τ0 = ||r0||2, θ0 = η0 = 0 4 wj = Avj 4 Chọn rg ≠ 0 bất kỳ, ρ0 = (rg,r0); 5 Doi: Do I = 1,, j 5 Do m = 0, 1, 2, 6 hi,j = (vi,wj), wj = wj – hi,jvi 6 If m chẵn thì tính 7 Enddo Doi 7 αm+1 = αm = ρm/(vm,rg) 8 hj+1,j = ||wj||2 8 um+1 = um - αmvm 9 vj+1 = wj/hj+1,j 9 Endif 10 Enddo 10 wm+1 = wm - αmAum 11 ym (là tối thiểu của ||βe1 – Hmy||2) 11 dm+1 = um + ( /αm)ηmdm; 12 ϕm = ϕ0 + vmym 12 θm+1 = ||wm+1||/τm; cm+1 = 1/(1+θm+1)1/2 13 If ϕm thoả mãn thì stop 13 τm+1 = τmθm+1cm+1; ηm+1 = αm; (1) (2) TÌM NGHIỆM CỦA PHƯƠNG TRÌNH POISSON BA CHIỀU BẰNG HAI PHƯƠNG PHÁP... 7 14 Else 14 ϕm+1 = ϕm ηm+1dm+1; 15 ϕ0 = ϕm 15 If m lẻ thì tính: 16 Goto 2 16 ρm+1 = (wm+1,rg); βm-1 = ρm+1/ρm-1 17 um+1 = wm+1 + βm-1um 18 vm+1 = Aum+1 + βm-1(Aum + βm-1vm-1) 19 Endif 20 Enddo 3. KẾT QUẢ MÔ PHỎNG VÀ THẢO LUẬN Để 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 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 – 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 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) 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 đó mỗi lớp có độ dày tương ứng là 340id nm= , 50p nd d nm= = . Mật độ pha tạp acceptor và donor tương ứng là 17 30.5 10AN cm −= × và 17 32.5 10DN cm −= × . Các hạt tải kích 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 xung là 12 sf và năng lượng photon là 1.49 eV , mật độ hạt tải quang là 16 35 10exN cm −= × sau thời gian 1 ps . Kích thước theo ba chiều không gian của đi-ốt là 440 100 100x y zL L L nm nm nm× × = × × , giả sử đi-ốt được nuôi cấy theo phương Ox . Linh kiện được chia thành các ô lưới không gian với 1050 10x y z m−Δ = Δ = Δ = × . Như vậy ta sẽ có 89xN = nút lưới theo phương Ox , 21yN = nút lưới theo phương Oy và 21zN = nút lưới theo phương Oz . Điện trường ngoài được đặt vào linh kiện dọc theo phương Ox và đi-ốt được phân cực nghịch, xem Hình 1. Hình 2 mô tả sự thay đổi 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 phần ứng với điện trường ngoài 100extE kV cm= , được tính toán với hai giải thuật giải 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 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 đều cho những kết quả tương tự nhau. Hình 1. Mô hình đi-ốt p-i-n GaAs ĐÀO HỮU HÀ - ĐINH NHƯ THẢO 8 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á trị điện trường ngoài 100extE kV cm= và 150kV cm , cũng được tính toán với hai giải 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ư trùng khớp. 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 như là hàm của thời gian ứng với 100extE kV cm= được tính toán với hai giải thuật giải phương trình Poisson khác nhau: a) GMRES(m), b) TFQMR 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ô phỏng ba chiều sử dụng thuật toán GMRES(m), thuật toán TFQMR và thuật toán BICGSTAB(3) [8], ba đồ thị gần như trùng nhau hoàn toàn, đặc biệt chương trình dùng 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 Hình 3. Vận tốc trôi dạt toàn phần của điện tử như là hàm của thời gian ứng với các điện trường ngoài khác nhau và được tính toán với hai giải thuật khác nhau Hình 4. Vận tốc trôi dạt của điện tử theo phương Ox như là hàm của thời gian thu được bằng ba chương trình mô phỏng khác nhau ứng với 100extE kV cm= TÌM NGHIỆM CỦA PHƯƠNG TRÌNH POISSON BA CHIỀU BẰNG HAI PHƯƠNG PHÁP... 9 toán BICGSTAB(3). Điều đó cho thấy phương pháp TFQMR hoạt động hiệu quả hơn phương pháp GMRES(m) do chương trình ba chiều sử dụng thuật toán BICGSTAB(3) đã được chứng minh là cho kết quả phù hợp với thực nghiệm [8]. 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 phương Ox , Oy tại mặt cắt 50z nm= ứng với điện trường ngoài 100extE kV cm= được 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), TFQMR và BICGSTAB(3). Ba đồ thị trong hình 5 đều gần như trùng nhau hoàn toàn. 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]. Hình 5. Phân bố điện thế không gian trong đi-ốt p-i-n bán dẫn GaAs theo hai phương Ox, Oy tại mặt cắt 50z nm= ứng với 100extE kV cm= , được tính toán bằng ba chương trình mô phỏng khác nhau Hình 6. 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 Poisson ứng với 100extE kV cm= Để 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 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 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 Poisson, Hình 6. Chuẩn Euclid của vectơ thặng dư được tính theo công thức [4]. 2 Tr r r= , (7) 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à 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 Euclid của vector thặng dư so với nghiệm tìm được bằng thuật toán BICGSTAB(3). Đ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 BICGSTAB(3). Tuy nhiên, đồ thị tương ứng với hai thuật toán TFQMR và GMRES(m) 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 cho kết quả ổn định hơn thuật toán BICGSTAB(3). ĐÀO HỮU HÀ - ĐINH NHƯ THẢO 10 Trong Bảng 2 chúng tôi so sánh chi phí tính toán trong một vòng lặp của chương trình con Poisson gồm có số phép nhân ma trận-vectơ (MVS), số lần cập nhật vectơ (AXPYS), số lần tính tích nội của hai vectơ (DOTS) [8] của ba thuật toán GMRES(m), TFQMR và BiCGSTAB(3). Bảng 2. Bảng so sánh chi phí tính toán trong một vòng lặp của chương trình con Poisson Thuật toán MVS AXPYS DOTS GMRES(m) 3 10 6 TFQMR 2 10 4 BiCGSTAB(3) 6 24 16 Ta thấy rằng số phép nhân ma trận-vectơ của thuật toán GMRES(m) chỉ bằng một nửa so với thuật toán BiCGSTAB(3) và gấp 1.5 lần so với thuật toán TFQMR, còn số lần tính tích nội thì ít hơn khoảng 2.7 lần so với thuật toán BiCGSTAB(3) và chỉ hơn 1.5 lần so với thuật toán TFQMR. Riêng số lần cập nhật vectơ thì bằng với thuật toán TFQMR và ít hơn 2.4 lần so với thuật toán BiCGSTAB(3). Kết quả cho thấy thuật toán TFQMR yêu cầu chi phí tính toán ít hơn các phương pháp khác. Bảng 3. Bảng so sánh số vòng lặp trung bình của chương trình con Poisson, thời gian trung bình trong một vòng lặp của chương trình mô phỏng và tổng thời gian mô phỏng Thuật toán loop T/loops T GMRES(m) 48706 264.4 giây 53.1 giờ TFQMR 317 1.43 giây 17.2 phút BiCGSTAB(3) 62 2.75 giây 33.1 phút Bảng 3 so sánh trung bình số vòng lặp của chương trình con Poisson (loop), thời gian trung bình trong một vòng lặp của chương trình mô phỏng (T/loops), tổng thời gian mô phỏng (T). Số liệu cho thấy số vòng lặp của chương trình con Poisson và thời gian mô phỏng đối với thuật toán GMRES(m) là rất lớn so với các giải thuật khác và là một hạn chế của thuật toán này. Ngược lại, mặc dù số vòng lặp của chương trình con Poisson sử dụng thuật toán TFQMR gấp khoảng 5 lần so với chương trình sử dụng thuật toán BiCGSTAB(3) nhưng thời gian mô phỏng trung bình trên một vòng lặp chỉ bằng một nửa so với thuật toán BiCGSTAB(3), dẫn đến việc tổng thời gian mô phỏng được rút ngắn. Đây chính là ưu thế của thuật toán TFQMR trong việc giải phương trình Poisson ba chiều. 4. KẾT LUẬN Chúng tôi đã xây dựng thành công hai chương trình giải phương trình Poisson ba chiều dựa trên hai thuật toán TFQMR và GMRES(m) dùng để tích hợp trong chương trình mô phỏng linh kiện na-nô bán dẫn bằng phương pháp Monte – Carlo tập hợp tự hợp. Để khảo sát các đặc trưng của phương pháp chúng tôi đã áp dụng để mô phỏng động lực học ba chiều của hạt tải trong các đi-ốt p-i-n bán dẫn GaAs và so sánh với các kết quả mô phỏng đã được công bố trước đây. Các kết quả chỉ ra rằng chương trình giải phương trình Poisson dựa trên hai thuật toán TFQMR và GMRES(m) có tốc độ hội tụ chậm hơn nhưng cho kết quả ổn định hơn chương trình giải phương trình Poisson dựa trên thuật toán BICGSTAB(3). TÌM NGHIỆM CỦA PHƯƠNG TRÌNH POISSON BA CHIỀU BẰNG HAI PHƯƠNG PHÁP... 11 LỜI CẢM ƠN Các tác giả xin gửi lời cảm ơn chân thành tới PGS. TS. Nguyễn Hồng Quang, Ban Hợp tác Quốc tế, Viện Khoa học và Công nghệ Việt Nam, chủ nhiệm đề tài NAFOSTED (103.02.109.09) đã hỗ trợ mọi mặt cho các tác giả trong quá trình hoàn thành bài báo này. TÀI LIỆU THAM KHẢO [1] D. N. Thao, S. Katayama, and K. Tomizawa (2004). Numerical simulation of THz radiation by coherent LO phonons in GaAs p-i-n diodes under high electric fields. Journal of the Physical Society of Japan 73, 3177–3181. [2] G. Klatt et al. (2011). Photo-Dember terahertz emitter excited with an Er: fiber laser, Appl. Phys. Lett. 98, 021114–021114-3. [3] K. Tomizawa (1993)., Numerical simulation of submicron semiconductor devices. Artech House, Boston London. [4] H. A. Vorst (2003). Iterative Krylov methods for large linear systems. Cambridge University. [5] Y. Saad (2000). Iterative Methods for Sparse Linear Systems. SIAM. [6] Shao-Liang Zhang (1997). GPBi-CG: Generalized product-type methods based on Bi-CG for solving nonsymmetric linear systems. SIAM J. Sci. Comput., Vol 18, No. 2, 537 – 551. [7] G. Speyer, D. Vasileska and S. M. Goodnick (2001). Efficient Poisson equation solvers for large scale 3D simulations. Technical Proceedings of the 2001 International Conference on Modeling and Simulation of Microsystems, Nanotech 2001, Vol. 1, 23 – 26. [8] D. N. Thao and L. H. Hai (2010). 3D simulation of semiconductor devices using BICGSTAB (3) for the solution of Poisson’s equation. Journal of Science and Education 15, 19-26. [9] D. N. Thao and N. T. Ngoc (2010). 3D simulation of semiconductor devices using preconditioned BICGSTAB algorithm with Jacobi preconditioner for the solution of Poisson’s equation. Journal of Science and Education, 16, 34-41. [10] D. N. Thao, D. T. D. My, N. C. P. Thi and N. T. Thuy (2011). Three- dimensional simulation of nano semiconductor devices using GPBICG algorithm for the solution of the Poisson's equation. Journal of Science 65, 215-223. ĐÀO HỮU HÀ - ĐINH NHƯ THẢO 12 Title: OBTAINING THE SOLUTION OF THE 3D POISSON’S EQUATION BY MEANS OF THE TFQMR AND GMRES(m) METHODS Abstract: The paper presents building approach for 3D Poisson solvers based on GMRES(m) and TFQMR algorithms for incorporating into self-consistent ensemble Monte Carlo simulation programs of nano semiconductor devices. In order to test the efficiency, the programs are used to simulate the well-known GaAs p-i-n diodes. The obtained results show that the TFQMR based Poisson solver owns much higher convergent rate compared to the GMRES(m) based solver. Both algorithms run more slowly than BICGSTAB(3) algorithm but in return, have much higher stability. TS. ĐINH NHƯ THẢO Khoa Vật lý, Trường Đại học Sư phạm – Đại học Huế ĐT: 0996.867.668. Email: dnthao@gmail.com ThS. ĐÀO HỮU HÀ Trường Phổ thông Dân tộc Nội trú Tỉnh Kon Tum ĐT: 0984.797.370. Email: dhha80@gmail.com