計算流體力學(xué)李新亮cfd第14講mpi并行程序設(shè)計初步PPT學(xué)習(xí)教案

1、會計學(xué)1 計算流體力學(xué)李新亮計算流體力學(xué)李新亮cfd第第14講講mpi并行并行 程序設(shè)計初步程序設(shè)計初步 2 計算節(jié)點(diǎn)計算節(jié)點(diǎn) a.exe a.exe a.exe a.exe 對等式對等式 設(shè)計設(shè)計 “對等式對等式”程序設(shè)計思想程序設(shè)計思想 如果我是其中一個進(jìn)程；如果我是其中一個進(jìn)程； 我應(yīng)當(dāng)做我應(yīng)當(dāng)做 完成我需要完成的任務(wù)完成我需要完成的任務(wù) 站在其中一個進(jìn)程的角度思考站在其中一個進(jìn)程的角度思考 第1頁/共62頁 3 第2頁/共62頁 發(fā)送發(fā)送 變量變量A 接收接收 到變量到變量B 配合使用配合使用 4 第3頁/共62頁 阻塞發(fā)送阻塞發(fā)送 開始開始 結(jié)束結(jié)束 消息成功發(fā)出消息成功發(fā)出 緩沖區(qū)

2、可釋放緩沖區(qū)可釋放 阻塞接收阻塞接收 開始開始 結(jié)束結(jié)束 消息成功接收消息成功接收 緩沖區(qū)數(shù)據(jù)可使緩沖區(qū)數(shù)據(jù)可使 用用 一、一、 阻塞式通信與非阻塞式通信阻塞式通信與非阻塞式通信 阻塞式發(fā)送與接收阻塞式發(fā)送與接收 MPI_Send( A, ) MPI_Recv( B , ) 5 第4頁/共62頁 MPI_SendMPI_Send( ) ( ) 返回后緩沖區(qū)可釋放返回后緩沖區(qū)可釋放 sum= sum= call call MPI_SendMPI_Send(sum,)(sum,) sum= sum= 變量可重復(fù)利用變量可重復(fù)利用 MPI_RecvMPI_Recv() () 返回后緩沖區(qū)數(shù)據(jù)可使用返

3、回后緩沖區(qū)數(shù)據(jù)可使用 Call Call MPI_RecvMPI_Recv(sum1,)(sum1,) Sum=sum0+sum1Sum=sum0+sum1 6 第5頁/共62頁 非阻塞發(fā)送非阻塞發(fā)送 啟動發(fā)送啟動發(fā)送 立即返回立即返回 計計 算算 通信完成通信完成 釋放發(fā)送緩沖區(qū)釋放發(fā)送緩沖區(qū) 發(fā)發(fā) 送送 消消 息息 非阻塞接收非阻塞接收 啟動接收啟動接收 立即返回立即返回 計計 算算 通信完成通信完成 引用接收數(shù)據(jù)引用接收數(shù)據(jù) 接接 收收 消消 息息 計算計算 與與 通信通信 重疊重疊 非阻塞消息發(fā)送與接收非阻塞消息發(fā)送與接收 7 第6頁/共62頁 8 第7頁/共62頁 非阻塞通信調(diào)用后立

4、即返回，緩沖區(qū)不能非阻塞通信調(diào)用后立即返回，緩沖區(qū)不能立即使用立即使用 Sum= 計算某變量計算某變量 MPI_Isend(sum .) 發(fā)送該變量發(fā)送該變量 sum= 不能給變量重新賦值不能給變量重新賦值 （發(fā)送可能尚未完成）（發(fā)送可能尚未完成） MPI_Irecv(sum1, ) sum=sum0+sum1 數(shù)據(jù)不能立即使用數(shù)據(jù)不能立即使用 （接收可能未完成）（接收可能未完成） MPI_Isend(sum, , request, ) Call MPI_Wait(request,status,ierr) Sum= MPI_Irecv(sum1, , request, ) Call MPI_W

5、ait(request,status,ierr) Sum=sum0+sum1 9 第8頁/共62頁 利用通信與計算重疊技術(shù)提高效率利用通信與計算重疊技術(shù)提高效率 例：例： 計算差分計算差分 串行程序串行程序 real A(N,N),B(N,N),h . Do i=1,N B(I,1)=(A(I,2)-A(I,1)/h B(I,N)=(A(I,N)-A(I,N-1)/h enddo Do j=2,N-1 Do i=1,N B(i,j)=(A(i,j+1)-A(i,j-1)/(2.*h) Enddo Enddo 0 )2/()()( 1,1, hfff jijijiy J=1,2,3 . N-1,

6、 N i=1 i=2 i=N 10 第9頁/共62頁 并行程序并行程序 以兩個進(jìn)程并行為例以兩個進(jìn)程并行為例 real A(N,N/2),B(N,N/2),A1(N),h If(myid .eq. 0) then call MPI_send(A(1,N/2),N,MPI_real,1,99,MPI_Comm_world,ierr) call MPI_recv(A1,N,MPI_real,1,99,MPI_Comm_World,status,ierr) Else call MPI_recv(A1,N,MPI_real,0,99,MPI_Comm_World,status,ierr) call M

7、PI_send(A(1,1),N,MPI_real,0,99,MPI_Comm_world,ierr) endif 01 J=1,2 N/2 A(1,N/2) A(2,N/2) A(3,N/2) A(N,N/ 2) 11 第10頁/共62頁 If(myid .eq. 0) then Do i=1,N B(i,1)=(A(i,2)-A(i,1)/h B(i,N)=(A1(i)-A(i,N-1)/(2.*h) Enddo Else Do i=1,N B(i,1)=(A(i,2)-A1(i)/(2.*h) B(i,N)=(A(i,N)-A(i,N-1)/h Enddo endif Do j=2,N-

8、1 Do i=1,N B(i,j)=(A(i,j+1)-A(i,j-1)/(2.*h) Enddo Enddo 01 J=1,2 N/2 特點(diǎn)：特點(diǎn)： 先收發(fā)邊界信息先收發(fā)邊界信息 再進(jìn)行計算再進(jìn)行計算 缺點(diǎn)：缺點(diǎn)： 通信過程中通信過程中CPU 空閑空閑 12 “內(nèi)邊界” 第11頁/共62頁 通信與計算重疊通信與計算重疊 real A(N,N/2),B(N,N/2),A1(N),h integer myid,ierr, req1, req2,status() If(myid .eq. 0) then call MPI_ISend(A(1,N/2),N,MPI_real,1,99,MPI_Com

9、m_world,req1, ierr) call MPI_Irecv(A1,N,MPI_real,1,99,MPI_Comm_World,req2,ierr) Else call MPI_Irecv(A1,N,MPI_real,0,99,MPI_Comm_World,req2,ierr) call MPI_Isend(A(1,1),N,MPI_real,0,99,MPI_Comm_world,req1,ierr) endif 01 J=1,2 N/2 13 第12頁/共62頁 Do j=2,N-1 Do i=1,N B(i,j)=(A(i,j+1)-A(i,j-1)/(2.*h) Enddo

10、Enddo Call MPI_wait(req2,statue,ierr) If(myid .eq. 0) then Do i=1,N B(I,1)=(A(I,2)-A(I,1)/h B(I,N)=(A1(i)-A(I,N-1)/(2.*h) Enddo Else Do i=1,N B(I,1)=(A(I,2)-A1(i)/(2.*h) B(I,N)=(A1(i)-A(I,N-1)/h Enddo endif 01 J=1,2 N/2 特點(diǎn)：特點(diǎn)： 傳遞邊界信息傳遞邊界信息 同時同時進(jìn)行計算進(jìn)行計算 內(nèi)點(diǎn)內(nèi)點(diǎn) 讀取系統(tǒng)時間讀取系統(tǒng)時間 doubleprecision time time=MPI

11、_Wtime( ) 14 第13頁/共62頁 二、二、 如何收發(fā)非連續(xù)數(shù)據(jù)如何收發(fā)非連續(xù)數(shù)據(jù) 例如：例如： 發(fā)送數(shù)組的一行發(fā)送數(shù)組的一行 A(100,50) 發(fā)送發(fā)送 A(1,1),A(1,2) ,A(1,3) A(1,1), A(1,2), A(1,3) 方法方法1. 多次發(fā)送多次發(fā)送 通信開銷大、效率低通信開銷大、效率低 A(1,1), A(2,1), A(1,2), A(2,2) . A(1,3). 15 第14頁/共62頁 方法方法2. 將發(fā)送的數(shù)據(jù)拷貝到連續(xù)的數(shù)組中將發(fā)送的數(shù)據(jù)拷貝到連續(xù)的數(shù)組中 dimension A(100,50), B(50) If(myid .eq. 0) t

12、hen Do i=1,50 B(i)=A(1,i) Enddo call MPI_Send(B,50,MPI_REAL,1,99,MPI_COMM_WORLD,ierr) Else call MPI_Recv(B,50,MPI_Real,0,99, ) Do i=1,50 A(1,i)=B(i) Enddo endif 不足：不足： 額外的內(nèi)存占用額外的內(nèi)存占用 額外的拷貝操作額外的拷貝操作 通信不復(fù)雜的情況，內(nèi)存拷貝工作量不大，該方法也可以采用。通信不復(fù)雜的情況，內(nèi)存拷貝工作量不大，該方法也可以采用。 效果還可以效果還可以 16 第15頁/共62頁 方法方法3： 構(gòu)建新的數(shù)據(jù)結(jié)構(gòu)構(gòu)建新的數(shù)據(jù)

13、結(jié)構(gòu) Count: 塊的數(shù)量；塊的數(shù)量； blocklength: 每塊的元素個數(shù)每塊的元素個數(shù) Stride: 跨度跨度 （各塊起始元素之間的距離）（各塊起始元素之間的距離） Oldtype: 舊數(shù)據(jù)類型，舊數(shù)據(jù)類型， Newtype: 新數(shù)據(jù)類型新數(shù)據(jù)類型 （整數(shù)）（整數(shù)） 例：例：integer MY_TYPE Call MPI_TYPE_VECTOR(4,1,3,MPI_REAL,MY_TYPE,ierr) Call MPI_TYPE_Commit(MY_TYPE,ierr) A(1,1), A(2,1) , A(3,1), A(1,2), A(2,2) , A(3,2), A(1,3

14、), A(2,3), A(3,3), A(1,4), A(2,4), A(3,4) Stride=3 固定間隔（跨度）的非連續(xù)數(shù)據(jù)固定間隔（跨度）的非連續(xù)數(shù)據(jù) MPI_TYPE_VECTOR(count ,blocklength, stride ,oldtype, newtype, ierr) A(1,1) A(1,2) A(1,3) A(1,4) A(2,1) A(2,2) A(2,3) A(2,4) A(3,1) A(3,2) A(3,3) A(3,4) 4塊，每塊塊，每塊1個元素，跨度為個元素，跨度為3（個元素）（個元素） Fortran 數(shù)組的一行數(shù)組的一行 Real A(3,4) .

15、 A(1,:) 在內(nèi)存中的排列次序 17 第16頁/共62頁 例：例： 發(fā)送三維數(shù)組中的一個面發(fā)送三維數(shù)組中的一個面 (Fortran) 數(shù)組：數(shù)組： real A(M,N,P) 通信通信 1) A(i,:,:) ; 2) A(:,j,:) ; 3) A(:,:,k) 通信通信1） A(1,1,1),A(2,1,1), A(3,1,1) ,A(M,1,1), A(1,2,1),A(2,2,1)., MPI_Type_Vector(N*P,1,M,MPI_Real, My_Type,ierr) 通信通信2) A(1,1,1),A(2,1,1), A(3,1,1) ., A(1,2,1),A(2,

16、2,1),A(3,2,1) , A(1,1,2),A(2,1,2),A(3,1,2) , MPI_Type_Vector(P,M,M*N,MPI_Real,My_Type,ierr) 通信通信3） 連續(xù)分布，無需構(gòu)造新類型連續(xù)分布，無需構(gòu)造新類型 18 第17頁/共62頁 MPI_TYPE_INDEXED(count, array_of_blocklengths, array_of_displacements, oldtype,newtype，ierr) 構(gòu)造數(shù)據(jù)類型更靈活的函數(shù)構(gòu)造數(shù)據(jù)類型更靈活的函數(shù) 直接指定每塊的元素個數(shù)及偏移量直接指定每塊的元素個數(shù)及偏移量 塊的數(shù)量（整數(shù)）塊的數(shù)量（整

17、數(shù)）每塊元素的個數(shù)每塊元素的個數(shù) （整形數(shù)組）（整形數(shù)組） 每塊的偏移量每塊的偏移量 （整形數(shù)組）（整形數(shù)組） 例：例： 數(shù)組數(shù)組 real A(N,N), 欲將其欲將其上三角元素上三角元素作為消息發(fā)送，試構(gòu)造其數(shù)據(jù)類型作為消息發(fā)送，試構(gòu)造其數(shù)據(jù)類型 A(1, 1) A(1, 2) A(1, 3) A(1, 4) A(2, 2) A(2, 3) A(2, 4) A(4, 4) A(3, 3) A(3, 4) A(2, 1) A(3, 1) A(3, 2) A(4, 1) A(4, 2) A(4, 3) A(1, 1) A(2, 1) A(1, 2) A(2, 2) A(3, 1) A(4, 1

18、) A(3, 2) A(4, 2) A(1, 3) A(2, 3) A(3, 3) A(4, 3) A(1, 4) A(2, 4) A(3, 4) A(4, 4) 內(nèi)存中的存儲次序 (Fortran) N列 N行 注意：注意： Fortran 行優(yōu)先次序存儲；行優(yōu)先次序存儲； C為列優(yōu)先次序存儲為列優(yōu)先次序存儲 觀察規(guī)律：觀察規(guī)律： N塊；塊； 第第k塊有塊有k個元素；第個元素；第k塊的偏移為塊的偏移為(k-1)*N （從（從0算起）算起） Integer: count, blocklengths(N), displacements(N) Integer: Newtype,ierr count

19、=N do k=1,N blocklengthes(k)=k displacements(k)=(k-1)*N enddo call MPI_TYPE_INDEXED(count, blocklengths, column=int(myid/3) MPI_Comm_Split(MPI_Comm_World, raw, 0,Comm_Raw) MPI_Comm_Split(MPI_Comm_World,column,0,Comm_column) Call MPI_Comm_rank(Comm_Raw,myid_raw,ierr) Call MPI_Comm_rank(Comm_line, myi

20、d_line,ierr) MPI_Comm_World RAW Column Color, 分組標(biāo)準(zhǔn) Key, 排序依據(jù) 如相同，按原ID排 提交新定義的組提交新定義的組 （否則新組無效，不要忘記）（否則新組無效，不要忘記） 計算行號、列號 20 第19頁/共62頁 例：例： 計算差分計算差分 三維分割三維分割 A(M1,N1,P1) (M1=M/NM, N1=N/NN, P1=P/NP) 基本思路：基本思路： 1） “擴(kuò)大擴(kuò)大”的數(shù)組的數(shù)組 A(0: M1+1, 0: N1+1,0:P1+1) 2）分割成三個組）分割成三個組 Comm_X, Comm_Y, Comm_Z 得到組內(nèi)編號得到組內(nèi)

21、編號 3)建立三個方向通訊的數(shù)據(jù)結(jié)構(gòu)建立三個方向通訊的數(shù)據(jù)結(jié)構(gòu) 4） 通信通信 ， 計算內(nèi)點(diǎn)差分計算內(nèi)點(diǎn)差分 5） 計算邊界差分計算邊界差分 z f y f x f , 0 2 14 3 5 7 6 8 9 10 11 MPI_Comm_World 21 第20頁/共62頁 Parameter(M1=M/NM,N1=N/NN,P1=P/NP) Real A(0:M1+1,0:N1+1,0:P1+1) Integer myid,Comm_X,Comm_Y,Comm_Z,id_X,id_Y,id_Z, request(12),. Call MPI_Comm_Rank(MPI_Comm_World,

22、myid,ierr) Call MPI_Comm_Split(MPI_Comm_World, mod(myid,NM),0,Comm_X,ierr) Call MPI_Comm_Split(MPI_Comm_World,mod(myid,NM*NN)/NM,0,Comm_Y,ierr) Call MPI_Comm_Split(MPI_Comm_World,myid/(NM*NN),0,Comm_Z,ierr) Call MPI_Comm_Rank(Comm_X,id_x,ierr) Call MPI_Comm_Rank(Comm_Y,id_y,ierr) Call MPI_Comm_Rank(

23、Comm_Z,id_z,ierr) 定義三個方向的通信域定義三個方向的通信域 22 第21頁/共62頁 Call MPI_Type_Vector(N1+2)*(P1+2),1,M1+2,MPI_real,Type_X,ierr) Call MPI_Type_Vector(P1+2,N1+2,(M1+2)*(N1+2),MPI_real,Type_Y,ierr) Call MPI_Type_Commit(Type_X,ierr) Call MPI_Type_Commit(Type_Y,ierr) . id_X_Pre=id_X-1, if(id_X_Pre .le. 0) id_X_pre=id

24、_X_Pre+NM Id_X_Next=id_X+1, if(id_X_Next .ge. NM) id_X_Next=id_X_Next-NM Call MPI_Isend(A(1,0,0) ,1,TYPE_X, id_X_Pre, 99,Comm_X,request(1),ierr) Call MPI_Isend(A(M1,0,0),1,TYPE_X,id_X_next,99,Comm_X,request(2),ierr) Call MPI_Irecv(A(0,0,0),1,TYPE_X,id_X_next,99,Comm_X,request(3),ierr) Call MPI_Irecv

25、(A(M1+1,0,0),1,TYPE_X,id_X_Pre,99,Comm_X,request(4),ierr) 定義新的數(shù)據(jù)結(jié)構(gòu)定義新的數(shù)據(jù)結(jié)構(gòu) 23 第22頁/共62頁 Do k=2,P1-1 Do j=2,N1-1 Do i=2,M1-1 Ax(I,j,k)=(A(i+1,j,k)-A(i-1,j,k)/(2.*hx) Ay(I,j,k)=(A(I,j+1,k)-A(I,j-1,k)/(2.*hy) Az(I,j,k)=(A(I,j,k+1)-A(I,j,k-1)/(2.*hz) Enddo Enddo Enddo call MPI_Wait_All(12,request,status

26、,ierr) do k=1,P1 do j=1,N1 Ax(1,j,k)=(A(2,j,k)-A(0,j,k)/(2.*hx) Ax(M1 ,j,k)=(A(M1+1,j,k)-A(M1-1,j,k)/(2.*hx) enddo Enddo . 內(nèi)點(diǎn)內(nèi)點(diǎn) 邊界點(diǎn)邊界點(diǎn) 24 第23頁/共62頁 四、分布數(shù)組的文件存儲四、分布數(shù)組的文件存儲 分布數(shù)組分布數(shù)組 real A(M/m1,N/n1)real A(M/m1,N/n1) 存儲方式存儲方式1. 每個進(jìn)程存儲到獨(dú)立的文件每個進(jìn)程存儲到獨(dú)立的文件 real A(M/m1,N/n1) character(len=50) filename writ

27、e(filename,”(file-I4.4.dat)”) myid open(55,file=filename,form=unformatted) write(55) A close(55) - file-0000.dat file-0001.dat file-0002.dat 優(yōu)點(diǎn)：程序簡單優(yōu)點(diǎn)：程序簡單 缺點(diǎn)：缺點(diǎn)： 數(shù)據(jù)文件多，不易處理；數(shù)據(jù)文件多，不易處理； 改變處理器數(shù)目時需特殊處理改變處理器數(shù)目時需特殊處理 0 1 2 3 25 第24頁/共62頁 分布數(shù)組分布數(shù)組 real A(M/m1,N/n1)real A(M/m1,N/n1) 存儲方式存儲方式2： 收集到收集到0節(jié)點(diǎn)存儲

28、節(jié)點(diǎn)存儲 存儲到存儲到 一個文件一個文件 缺點(diǎn)：缺點(diǎn)： 改變處理器規(guī)模時，需要處理改變處理器規(guī)模時，需要處理 存儲方式存儲方式3： 收集到收集到0節(jié)點(diǎn)，重新裝配成大數(shù)組節(jié)點(diǎn)，重新裝配成大數(shù)組 收集收集 A(M/m1,N/n1) 組成組成 A0(M,N) real A0(M,N), A(M/m1,N/n1), A1(M/m1,N/n1) if(myid.eq.0) then do k=0,m1*n1 call MPI_recv(A1, M/m1*N/n1,MPI_real,k,.) . A0( i_global, j_global ) = A1(i,j ) 把把A1 裝配到裝配到A0 enddo

29、 Write(33) A0 else call MPI_Send(A,) endif 0 1 2 3 0 1 2 3 0 26 第25頁/共62頁 存儲方式存儲方式4. 按列搜集后存儲按列搜集后存儲 Real Aj(M) If( myid .eq. 0) then open(33,file=“A.dat”,form= “binary”) do j=1,N 收集矩陣收集矩陣A0 的第的第 j 列存儲到列存儲到 Aj(:) write(33) Aj enddo Else endif 第第 1列列 第第 2列列 第第 3列列 優(yōu)點(diǎn)：優(yōu)點(diǎn)： 存儲的數(shù)據(jù)形式與內(nèi)存中存儲的數(shù)據(jù)形式與內(nèi)存中A0的存放格式一

30、致。的存放格式一致。 存儲的文件串行程序可直接讀取存儲的文件串行程序可直接讀取 real A(M,N) open(55,file=“A.dat”,form=“binary”) read(55) A close(55) 27 第26頁/共62頁 存儲方式存儲方式5 并行并行IO (MPI 2.0) 打開文件：打開文件： MPI_file_open(Comm,filename,mode,info,fileno,ierr) mode 打開類型：打開類型： MPI_Mode_RDONLY, MPI_Mode_RDWR, fileno 文件號，文件號， info 整數(shù)整數(shù) （信息）（信息） 關(guān)閉文件關(guān)閉

31、文件 ： MPI_file_close(fileno,ierr) 指定偏移位置讀寫指定偏移位置讀寫 MPI_file_read_at(fileno,offset,buff,const,datatype,status,ierr) MPI_file_write_at(fileno,offset,buff,const,datatype,status,ierr) offset 偏移，偏移， buff 緩沖區(qū)，緩沖區(qū)，const 數(shù)目數(shù)目 28 第27頁/共62頁 Part 3 實(shí)例教學(xué)實(shí)例教學(xué) CFD程序的程序的MPI實(shí)現(xiàn)實(shí)現(xiàn) 實(shí)例實(shí)例 （1） 用擬譜方法求解不可壓用擬譜方法求解不可壓N-S方程方程

32、實(shí)例（實(shí)例（2） 用流水線方法計算緊致差分用流水線方法計算緊致差分 常用的優(yōu)化方法常用的優(yōu)化方法 29 第28頁/共62頁 30 第29頁/共62頁 31 第30頁/共62頁 發(fā)送非連續(xù)數(shù)據(jù)發(fā)送非連續(xù)數(shù)據(jù)構(gòu)建新的數(shù)據(jù)結(jié)構(gòu)構(gòu)建新的數(shù)據(jù)結(jié)構(gòu) MPI_TYPE_VECTOR(count,blocklength,stride,oldtype,newtype,ierr) Count: 塊的數(shù)量；塊的數(shù)量； blocklength: 每塊的元素個數(shù)每塊的元素個數(shù) Stride: 跨度跨度 （各塊起始元素之間的距離）（各塊起始元素之間的距離） Oldtype: 舊數(shù)據(jù)類型，舊數(shù)據(jù)類型， Newtype: 新

33、數(shù)據(jù)類型新數(shù)據(jù)類型 （整數(shù)）（整數(shù)） 例：例：integer MY_TYPE Call MPI_TYPE_VECTOR(50,1,100,MPI_REAL,MY_TYPE,ierr) Call MPI_TYPE_Commit(MY_TYPE,ierr) A(1,1), A(2,1), A(1,2), A(2,2) . A(1,3). 32 第31頁/共62頁 通訊域的分割通訊域的分割 MPI_Comm_Split(comm,color, key,New_Comm ) 0 2 14 3 5 7 6 8 9 10 11 Color 相同的進(jìn)程在同一組相同的進(jìn)程在同一組 根據(jù)根據(jù)key的大小排序的大

34、小排序 例如：例如： 12個進(jìn)程，個進(jìn)程， 分成分成 3行行4列列 Line=mod(myid,3); raw=myid/3 MPI_Comm_Split(MPI_Comm_World, raw, 0,Comm_Raw) MPI_Comm_Split(MPI_Comm_World,line,Comm_Line) Call MPI_Comm_rank(Comm_Raw,myid_raw,ierr) Call MPI_Comm_rank(Comm_line, myid_line,ierr) MPI_Comm_World 33 第32頁/共62頁 實(shí)例實(shí)例 1. 用（擬）譜方法求解二維不可壓用（擬）

35、譜方法求解二維不可壓N-S方程方程 jj i jj i j i i i xx u x p x u u t u x u Re 1 0 2p 物理模型物理模型 周期性邊界條件周期性邊界條件 按照給定能譜布置初始流動按照給定能譜布置初始流動 研究流動的演化規(guī)律研究流動的演化規(guī)律 34 第33頁/共62頁 Fourier 變換變換 (1D) Njxeku N juuFu j N Nk kxi j /2,) ( 1 )( : ) ( 12/ 2/ p 1 0 1 )() ( : )( N j kxi j ejukuuFu Fourier 變換變換 的特點(diǎn)：的特點(diǎn)： 求導(dǎo)數(shù)求導(dǎo)數(shù) - 乘積乘積 )/( 1

36、 dxduFuik 0 i i x u 0 2211 ukuk 困難：困難： 非線性項非線性項 )( j i j x u uF )(*)( j i j x u FuF 卷積卷積 計算量巨大計算量巨大 在物理空間計算在物理空間計算 j i j x u u Fourier 變換的快速算法變換的快速算法FFT 35 第34頁/共62頁 二維 Fourier 變換變換 12/ 2/ 12/ 2/ 2 2 ), ( 1 ),( : ) ( N Nk myi kxi N Nm j ie emku N jiuuFu )()( ,2jixymk uFFuF 兩次一維兩次一維 Fourier 變換變換 jj i

37、 jj i j i i i xx u x p x u u t u x u Re 1 0 )2( Re 1 ) 1 (0 2 ii i ii ukpk t u uk )( j i j x u uF 0)2( iii ukconsiderk 2 /kkp i ii i ukkk t u Re 1 )/1 ( 22 36 第35頁/共62頁 ii i ukkk t u Re 1 )/1 ( 22 )( j i j x u uF 求解步驟：求解步驟： 1） 讀入初值讀入初值 2） 調(diào)用調(diào)用FFT 得到得到 3） 計算計算 調(diào)用調(diào)用FFT 得到得到 4） 計算計算 調(diào)用調(diào)用FFT 得到得到 5） 計算計

38、算 6） 積分積分 求出下一時間步的值求出下一時間步的值 7） 調(diào)用調(diào)用 FFT 得到得到 8） 循環(huán)循環(huán) 3）-7） 直到給定的時間直到給定的時間 ),(vuu i i u iju ki j i x u j i j x u u )( j i j x u uF iii ukkkQ Re 1 )/1 ( 22 i i Q t u )1( n i u )1(n i u 37 實(shí)際計算中，要采用抑制混淆誤差的措施 第36頁/共62頁 程序的并行化：程序的并行化： 二維二維 FFT )()( ,2jixymk uFFuF 二維二維FFT: 調(diào)用兩次一維調(diào)用兩次一維FFT 一維一維 FFT 算法復(fù)雜，并

39、行化難度大算法復(fù)雜，并行化難度大 二維二維 FFT 的并行：的并行： 重新分布重新分布 Subroutine FFT2d(nx,ny,u) integer nx,ny Complex u(nx,ny),Fu(nx,ny), u1(ny),u2(nx), do i=1,nx u1(:)= u(i,:) call FFT1d(ny,u1) Fu(i,:)=u1(:) enddo do j=1,ny u2(:)=Fu(:,j) call FFT1d(nx,u1) u(:,j)=u1(:) enddo end 38 第37頁/共62頁 數(shù)據(jù)重分布的實(shí)現(xiàn)數(shù)據(jù)重分布的實(shí)現(xiàn) A1(M/P,N) A2 (M,

40、N/P) 1234 a b c d 對等式編程思想對等式編程思想 “我我”需要完成的工作需要完成的工作 1) 將數(shù)據(jù)將數(shù)據(jù) A1(M/P,N) 切割成切割成P塊塊 ，存入數(shù)組，存入數(shù)組B1(M/P, N/P,P) 2) 將數(shù)據(jù)將數(shù)據(jù) B1(:,:,k) 發(fā)到發(fā)到 進(jìn)程進(jìn)程 k (k=0,1.P-1) 3) 從進(jìn)程從進(jìn)程k 接收接收 B2(:,:,k) 4) 組合組合 B2(:,:,k) 成成 A2 39 第38頁/共62頁 程序：程序： Subroutine Redistibute_ItoJ(A1,A2,M,N,P) Integer M,N,P,k,ierr,status(MPI_Status

41、_Size) real A1(M/P,N), A2(M,N/P), B1(M/P,N/P,P), B2(M/P,N/P,P) do k=1,P B1(:,:,P)=A1(:, (k-1)*N/P+1: k*N/P) ) call MPI_Send(B1,M*N/(P*P),MPI_Real, k-1, .) Enddo do k=1,P call MPI_Recv(B2,M*N/(P*P),MPI_Real, k-1, .) A2(k-1)*M/P+1: k*M/P) , : )=B2(:,:,P) Enddo end 問題：問題： 全部發(fā)送，全部發(fā)送， 發(fā)送成功發(fā)送成功后再啟動接收。后再啟動

42、接收。 容易死鎖容易死鎖 按行分布按行分布 - 按列分布按列分布 40 第39頁/共62頁 Subroutine Redistibute_ItoJ(A1,A2,M,N,P) Integer M,N,P,k,ierr,status(MPI_Status_Size) real A1(M/P,N), A2(M,N/P), B1(M/P,N/P,P), B2(M/P,N/P,P) do k=1,P B1(:,:,P)=A1(:, (k-1)*N/P+1: k*N/P) ) id_send=myid-k mod P id_recv= myid+k mod P call MPI_Send(B1,M*N/(

43、P*P),MPI_Real, id_send, .) call MPI_Recv(B2,M*N/(P*P),MPI_Real, id_recv, .) A2(k-1)*M/P+1: k*M/P) , : )=B2(:,:,P) Enddo end 問題：問題： 按順序發(fā)送、接收，不易死鎖按順序發(fā)送、接收，不易死鎖 41 第40頁/共62頁 數(shù)據(jù)全交換：數(shù)據(jù)全交換： MPI_AlltoAll(sendbuf,sendcount,sendtype,recvbuf,recvcount,recvtype,comm, ierr) sendbuf 發(fā)送緩沖區(qū)（首地址）發(fā)送緩沖區(qū)（首地址） recvbuf

44、接收緩沖區(qū)（首地址）接收緩沖區(qū)（首地址） sendcount 發(fā)送數(shù)目發(fā)送數(shù)目 recvcount 接收數(shù)目接收數(shù)目 sendtype 發(fā)送類型發(fā)送類型 recvtype 接收類型接收類型 Comm 通信域通信域 ierr 整數(shù)，返回錯誤值（整數(shù)，返回錯誤值（0為成功）為成功） To 0To 1To 2To 3 Sendbuf 的數(shù)據(jù)格式的數(shù)據(jù)格式 sendcount From 0From 1From 2From 3 Recvbuf 的數(shù)據(jù)格式的數(shù)據(jù)格式 recvcount 42 第41頁/共62頁 程序：程序： Subroutine Redistibute_ItoJ(A1,A2,M,N,P)

45、 Integer M,N,P,k,ierr,status(MPI_Status_Size) real A1(M/P,N), A2(M,N/P), B1(M/P,N/P,P), B2(M/P,N/P,P) do k=1,P B1(:,:,P)=A1(:, (k-1)*N/P+1: k*N/P) ) enddo call MPI_AlltoAll (B1,M*N/(P*P),MPI_Real, B2, M*N/(P*P),MPI_Real, MPI_Comm_World,ierr) do k=1,P A2(k-1)*M/P+1: k*M/P) , : )=B2(:,:,P) Enddo end 問

46、題：問題： 無法做到計算與通信重疊無法做到計算與通信重疊 43 第42頁/共62頁 二維二維 并行并行FFT 的實(shí)現(xiàn)的實(shí)現(xiàn) （輸入數(shù)據(jù)、輸出數(shù)據(jù)均為按列分布）（輸入數(shù)據(jù)、輸出數(shù)據(jù)均為按列分布） 1） 調(diào)用一維調(diào)用一維FFT實(shí)現(xiàn)實(shí)現(xiàn) i- 方向的變換方向的變換 u - u1 2） 重新分布數(shù)據(jù)重新分布數(shù)據(jù) （按列（按列- 按行）按行） u1 - u2 3) 調(diào)用一維調(diào)用一維FFT 實(shí)現(xiàn)實(shí)現(xiàn)j- 方向的變換方向的變換 u2- Fu2 4) 重新分布數(shù)據(jù)重新分布數(shù)據(jù) （按行（按行 - 按列）按列） Fu2- Fu 44 第43頁/共62頁 實(shí)例實(shí)例 （2） 利用流水線利用流水線 實(shí)現(xiàn)緊致差分的并行化

47、實(shí)現(xiàn)緊致差分的并行化 緊致型差分格式：緊致型差分格式： 相同網(wǎng)格點(diǎn)上引入更多信息。相同網(wǎng)格點(diǎn)上引入更多信息。 性能更優(yōu)化。性能更優(yōu)化。 0 x u a t u i F i x u 是是 的差分逼近的差分逼近 普通差分格式：普通差分格式： 顯式給出顯式給出 Fi 的表達(dá)式的表達(dá)式 )2/()( 11 huuF jjj )945459( 60 1 321123 jjjjjjj uuuuuu h F 緊致型差分格式：緊致型差分格式： 隱式給出隱式給出 Fi 的表達(dá)式的表達(dá)式 )1236443( 120 1 5 2 5 3 21121 jjjjjjj uuuuu x FF 6 階中心 6 階對稱緊致階

48、對稱緊致 (Lele) 5 階迎風(fēng)緊致階迎風(fēng)緊致 (Fu) j-2 j-1 j j+1 j+2 hffffFFF jjjjjjj 36/ )2828(3/13/1 211211 45 第44頁/共62頁 )945459( 60 1 321123 jjjjjjj uuuuuu h F 普通差分格式：普通差分格式： 直接計算導(dǎo)數(shù)，并行容易直接計算導(dǎo)數(shù)，并行容易 jjj dFF 1 5 2 5 3 緊致格式的計算：緊致格式的計算： 遞推遞推 )1236443( 120 1 2112 jjjjjj uuuuu x d 1 3/23/5 jjj FdF 遞推公式：遞推公式： 1)計算出計算出 （由邊界條

49、件或邊界格式給出）（由邊界條件或邊界格式給出） 2） 由由 遞推計算遞推計算 出全部導(dǎo)數(shù)出全部導(dǎo)數(shù) 1 F 1 3/23/5 jjj FdF 后面的數(shù)據(jù)必須等待前一步計算完成，無法并行后面的數(shù)據(jù)必須等待前一步計算完成，無法并行 46 第45頁/共62頁 1, 3/23/5 jijiji FdF 二維問題：二維問題： 流水線法求解流水線法求解 流水線示意圖流水線示意圖 步驟：步驟： 1) 計算計算 d(:,:) 2) for k=1,M 如果如果 myid=0, 計算計算 F(k,0), 否則否則 從從myid-1接收接收 F(k,0); for i=1,N1 (N1=N/P) 計算計算 F(k

50、,i); 如果如果myid P-1 向向 myid+1 發(fā)送發(fā)送 F(k,N1) 缺點(diǎn)：缺點(diǎn)： 通信次數(shù)過多通信次數(shù)過多 47 第46頁/共62頁 通信次數(shù)過于頻繁通信次數(shù)過于頻繁解決方法：解決方法： 分塊流水線分塊流水線 步驟：步驟： 1) 計算計算 d(:,:) 2) for kp=1,MP 如果如果 myid=0, 計算計算 F(kp,0), 否則否則 從從myid-1接收接收 F(kp,0); for j=1,N1 (N1=N/P) 計算計算 F(kp,j); 如果如果myid P-1 向向 myid+1 發(fā)送發(fā)送 F(kp,N1) F(kp,i) 表示第表示第kp 塊塊 48 第47

51、頁/共62頁 對稱緊致格式對稱緊致格式 jjjj dFFF 11 3/13/1hffffd jjjjj 36/ )2828( 2112 追趕法追趕法) 1 ( 11jjjjjj dFFF )2( 1jjjj BFAF 111 jjjj BFAF 令 則 代入（代入（1） 得得 111 ) 1( jjjjjjjj BdFFA 對比（對比（2） 得得 ) 3() 1/()(, ) 1/( 111 jjjjjjjjjj ABdBAA 邊界處導(dǎo)數(shù)可由邊界條件或邊界格式給出邊界處導(dǎo)數(shù)可由邊界條件或邊界格式給出： * 11NN FFFF 則 * 111 , 0FBA 步驟：步驟： 1） 2） 由由 （3）

52、式遞推，得到）式遞推，得到 3) 4) 由由 （2）式遞推，得到）式遞推，得到 jj BA , * NN FF j F 特點(diǎn)：特點(diǎn)： 兩次遞推。兩次遞推。 并行方法與前文類似并行方法與前文類似 49 第48頁/共62頁 3.常用的并行優(yōu)化方法常用的并行優(yōu)化方法 1） 通信與計算重疊通信與計算重疊 采用非阻塞通信采用非阻塞通信 Isend, Irecv 2） 用重復(fù)計算代替通信用重復(fù)計算代替通信 3） 拆分長消息、合并短消息拆分長消息、合并短消息 4） 優(yōu)化通信方式優(yōu)化通信方式 50 第49頁/共62頁 用重復(fù)計算代替通信用重復(fù)計算代替通信 例如：例如： 計算差分計算差分 u 分布存儲分布存儲,

53、 f(u) 為為 u 的函數(shù)的函數(shù) 01 )(,uf xx u 方法方法 1） 計算出計算出 v=f(u) 通信得到通信得到 uN+1, vN+1 計算差分計算差分 方法方法 2） 計算出計算出 v=f(u) 通信得到通信得到 uN+1 （邊界外）（邊界外） 計算出計算出 vN+1 =f(uN+1) 計算差分計算差分 方法方法 2） 計算量大，通信量小計算量大，通信量小 當(dāng)函數(shù)當(dāng)函數(shù) f(u)不復(fù)雜時，可提高效率不復(fù)雜時，可提高效率 1, 2 N N+1 51 第50頁/共62頁 長消息切割成多個短消息發(fā)送、接收長消息切割成多個短消息發(fā)送、接收 call MPI_Send(A(1),10000

54、0, MPI_Real, 1, ) 改為：改為： do m=1,10 call MPI_Send(A(m-1)*10000+1),10000,MPI_real,1 ) enddo 長消息：長消息： 非緩沖；非緩沖； 短消息：短消息： 緩沖緩沖 緩沖區(qū) MPI_Send 緩沖區(qū) MPI_Send MPI_Recv MPI_Recv 52 第51頁/共62頁 合并短消息合并短消息 do m=1,100 call MPI_Send(A(1,m),1,MPI_real,1 ) enddo 改為改為 do m=1,100 B(m)=A(1,m) enddo call MPI_Send(B(1),100,

55、 MPI_Real, 1, ) 53 第52頁/共62頁 優(yōu)化通信方式優(yōu)化通信方式 例：例： 數(shù)據(jù)散發(fā)數(shù)據(jù)散發(fā) 0號號 進(jìn)程：進(jìn)程： 數(shù)據(jù)數(shù)據(jù) A(100)， 散發(fā)給散發(fā)給 0-99 方式方式1） 0 進(jìn)程執(zhí)行進(jìn)程執(zhí)行100次次 MPI_Send 其他進(jìn)程執(zhí)行其他進(jìn)程執(zhí)行 MPI_Recv MPI_Scatter() 采用該算法采用該算法 方式方式 2） 0 進(jìn)程進(jìn)程 把把 A(100) 切割成切割成10份份 ， 發(fā)送給發(fā)送給10個進(jìn)程個進(jìn)程 10個進(jìn)程接收個進(jìn)程接收A1(10) 后再散發(fā)后再散發(fā) 54 第53頁/共62頁 OpenMP并行編程入門并行編程入門 一、一、 特點(diǎn)特點(diǎn) 1. 針對針

56、對共享內(nèi)存共享內(nèi)存計算機(jī)結(jié)構(gòu)計算機(jī)結(jié)構(gòu) 全部全部CPU/線程均可訪問內(nèi)存線程均可訪問內(nèi)存 2. 程序改動量小、實(shí)現(xiàn)方便程序改動量小、實(shí)現(xiàn)方便 （以編譯指示符為主）（以編譯指示符為主） 3. 適用于小規(guī)模并行或與適用于小規(guī)模并行或與MPI配合配合 進(jìn)行大規(guī)模并行進(jìn)行大規(guī)模并行 內(nèi)存 CPU（ 核心） CPU（ 核心） CPU（ 核心） 1臺PC機(jī) / 1個計算節(jié)點(diǎn) （共享內(nèi)存構(gòu)架） CPU 內(nèi)存 CPU 內(nèi)存 CPU 內(nèi)存 外部網(wǎng)絡(luò) 節(jié)點(diǎn)1節(jié)點(diǎn)2 Cluster結(jié)構(gòu)， 分布內(nèi)存構(gòu)架 第54頁/共62頁 print*, code 1 !$OMP PARALLEL print*, code 2 !$

57、OMP END PARALLEL print*, code 3 end 例1 (test1.f90) ： 編譯編譯 （在深騰（在深騰7000) 運(yùn)行結(jié)果（屏幕截圖） ifort test1.f90 -openmp 添加 openmp 選項 運(yùn)行： 1. 設(shè)置線程數(shù)（并行執(zhí)行的數(shù)目） export OMP_NUM_THREADS=4 （例如，4個） 2. 執(zhí)行： ./a.out 顯示結(jié)果： code 1 code 2 code 2 code 2 code 2 code 3 并行域中的代碼執(zhí)行了4次 Test2.f90 : print*, code 1 !$OMP PARALLEL print*,

58、 code 2“ !$OMP PARALLEL print*, “code 3” !$OMP END PARALLEL !$OMP END PARALLEL print*, code 4 end 第55頁/共62頁 3.DO 循環(huán)分解循環(huán)分解 （openMP最常用的并行方法）最常用的并行方法） !$OMP PARALLEL !$OMP DO do k=1,12 print*, k enddo !$OMP END DO !$OMP END PARALLEL end 示例： 線程線程0 k=1,2,3 線程線程1 k=4,5,6 線程線程2 k=7,8,9 線程線程2 k=10,11,12 !$O

59、MP PARALLEL !$OMP DO !$OMP PARALLEL DO 簡寫 運(yùn)行結(jié)果 （屏幕截圖） 運(yùn)行結(jié)果： 1 2 3 7 8 9 4 5 6 10 11 12 線程0 線程2 線程1 線程3 第56頁/共62頁 implicit none integer,parameter: N=100000000 integer: k real*8,dimension(:),allocatable: x,y,z real*8: time1,time2,OMP_get_wtime allocate(x(N),y(N),z(N) !$ time1=OMP_get_wtime() !$OMP PAR

60、ALLEL DO SHARED(x,y,z) PRIVATE(k) do k=1,N x(k)=(k-1.d0)/(N-1.d0) y(k)=(k+1.d0)/(N-1.d0) z(k)=x(k)+y(k) enddo !$OMP END PARALLEL DO !$ time2=OMP_get_wtime() deallocate(x,y,z) print*, Total Wall Time is , time2-time1 end 例：test 4屏幕截圖 采用單線程執(zhí)行：采用單線程執(zhí)行： 耗時耗時2.15秒秒 采用采用2線程執(zhí)行：耗時線程執(zhí)行：耗時 1.43秒秒 采用采用4線程執(zhí)行：線程