1、計算流體力學講義 第六講 MPI并行程序設計 （2） ；力學所主樓219； 82543801 知識點： 阻塞通信與非阻塞通信 非連續數據的發送與接收 OpenMP并行程序設計初步 講義、課件上傳至 （流體中文網） - “流體論壇” -“ CFD基礎理論”也可到如下網址下載：1服務器/前端機計算節點a.exea.exea.exeMPI 程序的運行原理： 服務器（前端機）編譯 可執行代碼復制 N 份，每個節點運行一份 調用MPI庫函數 得到每個節點號 my_id 根據my_id 不同，程序執行情況不同 調用MPI 庫函數進行通訊MPI 編程的基本思想： 主從式，對等式重點：對等式程序設計知識回顧2

2、計算節點a.exea.exea.exea.exe對等式設計“對等式”程序設計思想如果我是其中一個進程；我應當做完成我需要完成的任務站在其中一個進程的角度思考3基本的MPI函數（6個）MPI初始化 MPI_Init(ierr) ； MPI結束 MPI_Finalize(ierr)得到當前進程標識 MPI_Comm_rank(MPI_COMM_WORLD,myid,ierr) 得到通信域包含的進程數 MPI_Comm_size(MPI_COMM_WORLD,numprocs,ierr) 消息發送MPI_Send(buf,count,datatype,dest,tag,comm, ierr)消息接收

3、MPI_Recv(buf,count,datatype,source,tag,comm,status,ierr) 4MPI的消息發送機制 兩步進行MPI_Send( A, ) 發送MPI_Recv( B, ) 接收 發送 變量A接收 到變量B配合使用5阻塞發送開始結束消息成功發出緩沖區可釋放阻塞接收開始結束消息成功接收緩沖區數據可使用一、 阻塞式通信與非阻塞式通信阻塞式發送與接收MPI_Send( A, )MPI_Recv( B , )6 MPI_Send( ) 返回后緩沖區可釋放 sum= call MPI_Send(sum,) sum= 變量可重復利用 MPI_Recv() 返回后緩沖區數

4、據可使用Call MPI_Recv(sum1,)Sum=sum0+sum1 7非阻塞發送啟動發送立即返回計算通信完成釋放發送緩沖區發送消息非阻塞接收啟動接收立即返回計算通信完成引用接收數據接收消息計算與通信重疊非阻塞消息發送與接收8非阻塞消息發送MPI_ISend(buf,count,datatype,dest,tag,comm,request,ierr)In buf,count,datatype,dest,tag,commOut request,ierr Request (返回的非阻塞通信對象, 整數)非阻塞消息接收MPI_IRecv(buf,count,datatype,source,ta

5、g,comm,request,ierr)In buf,count,datatype,source,tag,commOut request,ierr非阻塞通信的完成 MPI_Wait(request,status,ierr) 等待消息收發完成 MPI_Test(request, flag,stutus,ierr) MPI_Waitall(const,request_array,status,ierr) 等待多個消息完成 In requestOut status， flag (logical型)9非阻塞通信調用后立即返回，緩沖區不能立即使用Sum= 計算某變量MPI_Isend(sum .) 發送

6、該變量 sum= 不能給變量重新賦值 （發送可能尚未完成）MPI_Irecv(sum1, )sum=sum0+sum1 數據不能立即使用 （接收可能未完成）MPI_Isend(sum, , request, )Call MPI_Wait(request,status,ierr)Sum= MPI_Irecv(sum1, , request, )Call MPI_Wait(request,status,ierr)Sum=sum0+sum1 10利用通信與計算重疊技術提高效率例： 計算差分串行程序 real A(N,N),B(N,N),h.Do i=1,NB(I,1)=(A(I,2)-A(I,1)/

7、h B(I,N)=(A(I,N)-A(I,N-1)/henddoDo j=2,N-1Do i=1,NB(i,j)=(A(i,j+1)-A(i,j-1)/(2.*h)EnddoEnddo0J=1,2,3 . N-1, Ni=1i=2 i=N11并行程序 以兩個進程并行為例real A(N,N/2),B(N,N/2),A1(N),hIf(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

8、,ierr)Else call MPI_recv(A1,N,MPI_real,0,99,MPI_Comm_World,status,ierr) call MPI_send(A(1,1),N,MPI_real,0,99,MPI_Comm_world,ierr)endif01J=1,2 N/2A(1,N/2)A(2,N/2)A(3,N/2)A(N,N/2)12If(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) EnddoElse Do i=1,N B(i,1)=(A(i,2)-A1(

9、i)/(2.*h) B(i,N)=(A(i,N)-A(i,N-1)/h EnddoendifDo j=2,N-1Do i=1,NB(i,j)=(A(i,j+1)-A(i,j-1)/(2.*h)EnddoEnddo01J=1,2 N/2特點： 先收發邊界信息 再進行計算缺點： 通信過程中CPU 空閑“內邊界”13通信與計算重疊real A(N,N/2),B(N,N/2),A1(N),hinteger myid,ierr, req1, req2,status()If(myid .eq. 0) then call MPI_ISend(A(1,N/2),N,MPI_real,1,99,MPI_Comm

10、_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)endif01J=1,2 N/214Do j=2,N-1Do i=1,NB(i,j)=(A(i,j+1)-A(i,j-1)/(2.*h)EnddoEnddoCall MPI_wait(re

11、q2,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) EnddoElse Do i=1,N B(I,1)=(A(I,2)-A1(i)/(2.*h) B(I,N)=(A1(i)-A(I,N-1)/h Enddoendif01J=1,2 N/2特點： 傳遞邊界信息 同時進行計算內點讀取系統時間 doubleprecision time time=MPI_Wtime( ) 15二、 如何收發非連續數據例如： 發送數組的一行A(100,50)發送 A(1,1)

12、,A(1,2) ,A(1,3)A(1,1), A(1,2), A(1,3) 方法1. 多次發送 通信開銷大、效率低A(1,1), A(2,1), A(1,2), A(2,2) . A(1,3).16方法2. 將發送的數據拷貝到連續的數組中dimension A(100,50), B(50)If(myid .eq. 0) then 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(

13、1,i)=B(i) Enddoendif不足： 額外的內存占用 額外的拷貝操作通信不復雜的情況，內存拷貝工作量不大，該方法也可以采用。效果還可以17方法3： 構建新的數據結構 Count: 塊的數量； blocklength: 每塊的元素個數Stride: 跨度 （各塊起始元素之間的距離）Oldtype: 舊數據類型， Newtype: 新數據類型 （整數）例：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

14、(3,1), A(1,2), A(2,2) , A(3,2), A(1,3), A(2,3), A(3,3), A(1,4), A(2,4), A(3,4)Stride=3固定間隔（跨度）的非連續數據 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 數組的一行Real A(3,4). A(1,:

15、)在內存中的排列次序18例： 發送三維數組中的一個面 (Fortran) 數組： 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,2,1),A(3,2,1) , A(1,1,2),A(2,1,2),A(3,1,2)

16、 , MPI_Type_Vector(P,M,M*N,MPI_Real,My_Type,ierr)通信3） 連續分布，無需構造新類型 19MPI_TYPE_INDEXED(count, array_of_blocklengths, array_of_displacements, oldtype,newtype，ierr)構造數據類型更靈活的函數 直接指定每塊的元素個數及偏移量塊的數量（整數）每塊元素的個數（整形數組）每塊的偏移量（整形數組）例： 數組 real A(N,N), 欲將其上三角元素作為消息發送，試構造其數據類型 A(1,1)A(1,2)A(1,3)A(1,4)A(2,2)A(2,3

17、)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)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)內存中的存儲次序(Fortran)N列N行注意： Fortran 行優先次序存儲； C為列優先次序存儲觀察規律： N塊； 第k塊有k個元素；第k塊的偏移為(k-1)*N （從0算起）Integer: count, blocklengths(N), displacements(N)Int

18、eger: Newtype,ierr count=N do k=1,N blocklengthes(k)=k displacements(k)=(k-1)*N enddocall MPI_TYPE_INDEXED(count, blocklengths, & displacements,MPI_REAL,newtype，ierr)Call MPI_TYPE_Commit(Newtype, ierr) call MPI_Send (A(1,1),1,Newtype, )N20三、 MPI的通信域和組預定義通訊域 MPI_Comm_World : 包含所有進程的組通訊域的分割 MPI_Comm_S

19、plit(comm,color, key,New_Comm ) 02143576891011Color 相同的進程在同一組根據key的大小排序 （key相同時按原ID排序）例如： 12個進程， 分成 3行4列Integer myid, Comm_Raw,Comm_column,myid_raw,myid_line,ierr,raw,columnRaw=mod(myid,3); column=int(myid/3)MPI_Comm_Split(MPI_Comm_World, raw, 0,Comm_Raw)MPI_Comm_Split(MPI_Comm_World,column,0,Comm_c

20、olumn)Call MPI_Comm_rank(Comm_Raw,myid_raw,ierr)Call MPI_Comm_rank(Comm_line, myid_line,ierr)MPI_Comm_WorldRAWColumnColor, 分組標準Key, 排序依據如相同，按原ID排 提交新定義的組（否則新組無效，不要忘記）計算行號、列號21例： 計算差分 三維分割 A(M1,N1,P1)(M1=M/NM, N1=N/NN, P1=P/NP)基本思路：1） “擴大”的數組 A(0: M1+1, 0: N1+1,0:P1+1)2）分割成三個組 Comm_X, Comm_Y, Comm_Z

21、得到組內編號 建立三個方向通訊的數據結構4） 通信 ， 計算內點差分5） 計算邊界差分02143576891011MPI_Comm_World22Parameter(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,myid,ierr)Call MPI_Comm_Split(MPI_Comm_World, mod(myid,NM),0,Comm_X,

22、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(Comm_Z,id_z,ierr)定義三個方向的通信域23Call MPI_Type_Vector(N1+2)*(P1+2),1,M1+2,MPI_

23、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_X_Pre+NMId_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,TY

24、PE_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(A(M1+1,0,0),1,TYPE_X,id_X_Pre,99,Comm_X,request(4),ierr) 定義新的數據結構24 Do k=2,P1-1 Do j=2,N1-1 Do i=2,M1-

25、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)EnddoEnddoEnddo call MPI_Wait_All(12,request,status,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 .內點邊

26、界點25四、分布數組的文件存儲 分布數組 real A(M/m1,N/n1) 存儲方式1. 每個進程存儲到獨立的文件 real A(M/m1,N/n1) character(len=50) write(,”(file-I4.4.dat)”) myid open(55,unformatted) write(55) A close(55) - 優點：程序簡單缺點： 數據文件多，不易處理； 改變處理器數目時需特殊處理 012326 分布數組 real A(M/m1,N/n1) 存儲方式2： 收集到0節點存儲 存儲到 一個文件 缺點： 改變處理器規模時，需要處理存儲方式3： 收集到0節點，重新裝配成大

27、數組 收集 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 Write(33) A0 else call MPI_Send(A,) endif 01230123027存儲方式4. 按列搜集后存儲 Real Aj(M)If( myid .eq. 0) then

28、open(33,file=“A.dat”,form= “binary”) do j=1,N 收集矩陣A0 的第 j 列存儲到 Aj(:) write(33) Aj enddoElse endif第 1列 第 2列 第 3列優點： 存儲的數據形式與內存中A0的存放格式一致。 存儲的文件串行程序可直接讀取 real A(M,N) open(55,file=“A.dat”,form=“binary”) read(55) A close(55)28存儲方式5 并行IO (MPI 2.0) 打開文件： MPI_(Comm,) mode 打開類型： MPI_Mode_RDONLY, MPI_Mode_RD

29、WR, fileno 文件號， info 整數 （信息） 關閉文件 ： MPI_() 指定偏移位置讀寫 MPI_() MPI_() offset 偏移， buff 緩沖區，const 數目 29Part 3 實例教學 CFD程序的MPI實現實例 （1） 用擬譜方法求解不可壓N-S方程 實例（2） 用流水線方法計算緊致差分 常用的優化方法30回顧 基本的MPI函數（6個）MPI初始化 MPI_Init(ierr) ； MPI結束 MPI_Finalize(ierr)得到當前進程標識 MPI_Comm_rank(MPI_COMM_WORLD,myid,ierr) 得到通信域包含的進程數 MPI_C

30、omm_size(MPI_COMM_WORLD,numprocs,ierr) 消息發送MPI_Send(buf,count,datatype,dest,tag,comm, ierr)消息接收MPI_Recv(buf,count,datatype,source,tag,comm,status,ierr) 31非阻塞消息發送MPI_ISend(buf,count,datatype,dest,tag,comm,request,ierr)In buf,count,datatype,dest,tag,commOut request,ierr Request (返回的非阻塞通信對象, 整數)非阻塞消息接收

31、MPI_IRecv(buf,count,datatype,source,tag,comm,request,ierr)In buf,count,datatype,source,tag,commOut request,ierr非阻塞通信的完成 MPI_Wait(request,status,ierr) 等待消息收發完成 MPI_Test(request, flag,stutus,ierr) MPI_Waitall(const,request_array,status,ierr) 等待多個消息完成 In requestOut status， flag (logical型)32發送非連續數據構建新的數

32、據結構MPI_TYPE_VECTOR(count,blocklength,stride,oldtype,newtype,ierr)Count: 塊的數量； blocklength: 每塊的元素個數Stride: 跨度 （各塊起始元素之間的距離）Oldtype: 舊數據類型， Newtype: 新數據類型 （整數）例：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

33、).33通訊域的分割 MPI_Comm_Split(comm,color, key,New_Comm ) 02143576891011Color 相同的進程在同一組根據key的大小排序例如： 12個進程， 分成 3行4列Line=mod(myid,3); raw=myid/3MPI_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

34、, myid_line,ierr)MPI_Comm_World34實例 1. 用（擬）譜方法求解二維不可壓N-S方程2p物理模型周期性邊界條件按照給定能譜布置初始流動 研究流動的演化規律35Fourier 變換 (1D)Fourier 變換 的特點： 求導數 - 乘積困難： 非線性項卷積計算量巨大在物理空間計算Fourier 變換的快速算法FFT36二維 Fourier 變換兩次一維 Fourier 變換37求解步驟： 1） 讀入初值 2） 調用FFT 得到 3） 計算 調用FFT 得到 4） 計算 調用FFT 得到 5） 計算 6） 積分 求出下一時間步的值 7） 調用 FFT 得到 8）

35、循環 3）-7） 直到給定的時間 實際計算中，要采用抑制混淆誤差的措施38程序的并行化： 二維 FFT二維FFT: 調用兩次一維FFT一維 FFT 算法復雜，并行化難度大二維 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(:

36、) enddo end39數據重分布的實現A1(M/P,N)A2 (M,N/P)1234abcd對等式編程思想 “我”需要完成的工作 1) 將數據 A1(M/P,N) 切割成P塊 ，存入數組B1(M/P, N/P,P) 2) 將數據 B1(:,:,k) 發到 進程 k (k=0,1.P-1) 3) 從進程k 接收 B2(:,:,k) 4) 組合 B2(:,:,k) 成 A240程序：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

37、/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, .) Enddodo 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) Enddoend 問題： 全部發送， 發送成功后再啟動接收。 容易死鎖 按行分布 - 按列分布41Subroutine Redistibute_ItoJ(

38、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/(P*P),MPI_Real, id_send, .) call MPI_Recv(B2,M*N/(P*P),MPI_Real, id_recv, .)

39、 A2(k-1)*M/P+1: k*M/P) , : )=B2(:,:,P) Enddoend 問題： 按順序發送、接收，不易死鎖42數據全交換：MPI_AlltoAll(sendbuf,sendcount,sendtype,recvbuf,recvcount,recvtype,comm,ierr) sendbuf 發送緩沖區（首地址） recvbuf 接收緩沖區（首地址） sendcount 發送數目 recvcount 接收數目 sendtype 發送類型 recvtype 接收類型 Comm 通信域 ierr 整數，返回錯誤值（0為成功） To 0To 1To 2To 3Sendbuf

40、的數據格式sendcountFrom 0From 1From 2From 3Recvbuf 的數據格式recvcount43程序：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) ) enddo call MPI_AlltoAll (B1,M*N/(P*P),MPI_Real,

41、 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) Enddoend 問題： 無法做到計算與通信重疊 44 二維 并行FFT 的實現 （輸入數據、輸出數據均為按列分布）1） 調用一維FFT實現 i- 方向的變換 u - u12） 重新分布數據 （按列- 按行） u1 - u2 調用一維FFT 實現j- 方向的變換 u2- Fu2 重新分布數據 （按行 - 按列） Fu2- Fu45實例 （2） 利用流水線 實現緊致差分的并行化緊致型差分格式： 相同網格點上引入更

42、多信息。 性能更優化。 是 的差分逼近普通差分格式： 顯式給出 Fi 的表達式緊致型差分格式： 隱式給出 Fi 的表達式 6 階中心6 階對稱緊致 (Lele)5 階迎風緊致 (Fu) j-2 j-1 j j+1 j+246 普通差分格式： 直接計算導數，并行容易緊致格式的計算： 遞推遞推公式：計算出 （由邊界條件或邊界格式給出）2） 由 遞推計算 出全部導數 后面的數據必須等待前一步計算完成，無法并行47二維問題： 流水線法求解 流水線示意圖步驟： 1) 計算 d(:,:) 2) for k=1,M 如果 myid=0, 計算 F(k,0), 否則 從myid-1接收 F(k,0); for

43、 i=1,N1 (N1=N/P) 計算 F(k,i); 如果myid P-1 向 myid+1 發送 F(k,N1) 缺點： 通信次數過多48通信次數過于頻繁解決方法： 分塊流水線步驟： 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(kp,N1) F(kp,i) 表示第kp 塊 49對稱緊致格式追趕法令則代入（1） 得對比（2） 得邊界處導數可由邊界條件或邊界格式給出： 則步驟： 1

44、） 2） 由 （3）式遞推，得到 3) 4) 由 （2）式遞推，得到特點： 兩次遞推。 并行方法與前文類似50常用的并行優化方法 1） 通信與計算重疊 采用非阻塞通信 Isend, Irecv 2） 用重復計算代替通信 3） 拆分長消息、合并短消息 4） 優化通信方式 51 用重復計算代替通信 例如： 計算差分 u 分布存儲, f(u) 為 u 的函數 01方法 1） 計算出 v=f(u) 通信得到 uN+1, vN+1 計算差分方法 2） 計算出 v=f(u) 通信得到 uN+1 （邊界外） 計算出 vN+1 =f(uN+1) 計算差分 方法 2） 計算量大，通信量小 當函數 f(u)不復雜

45、時，可提高效率1, 2 N N+152 長消息切割成多個短消息發送、接收 call MPI_Send(A(1),100000, MPI_Real, 1, ) 改為： do m=1,10 call MPI_Send(A(m-1)*10000+1),10000,MPI_real,1 ) enddo 長消息： 非緩沖； 短消息： 緩沖 緩沖區MPI_Send緩沖區MPI_SendMPI_RecvMPI_Recv53合并短消息 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 cal

46、l MPI_Send(B(1),100, MPI_Real, 1, ) 54 優化通信方式 例： 數據散發0號 進程： 數據 A(100)， 散發給 0-99方式1） 0 進程執行100次 MPI_Send 其他進程執行 MPI_Recv MPI_Scatter() 采用該算法方式 2） 0 進程 把 A(100) 切割成10份 ， 發送給10個進程 10個進程接收A1(10) 后再散發 55OpenMP并行編程入門一、 特點 1. 針對共享內存計算機結構 全部CPU/線程均可訪問內存 2. 程序改動量小、實現方便 （以編譯指示符為主） 3. 適用于小規模并行或與MPI配合 進行大規模并行 內

47、存CPU（核心）CPU（核心）CPU（核心）1臺PC機 / 1個計算節點 （共享內存構架）CPU內存CPU內存CPU內存外部網絡節點1節點2Cluster結構， 分布內存構架56 print*, code 1!$OMP PARALLEL print*, code 2!$OMP END PARALLEL print*, code 3 end例1 (test1.f90) ：編譯 （在深騰7000)運行結果（屏幕截圖）ifort test1.f90 -openmp添加 openmp 選項運行： 1. 設置線程數（并行執行的數目） export OMP_NUM_THREADS=4 （例如，4個） 2.

48、 執行： ./a.out顯示結果： code 1 code 2 code 2 code 2 code 2 code 3并行域中的代碼執行了4次Test2.f90 : print*, code 1!$OMP PARALLEL print*, code 2“!$OMP PARALLEL print*, “code 3”!$OMP END PARALLEL!$OMP END PARALLEL print*, code 4 end57DO 循環分解 （openMP最常用的并行方法）!$OMP PARALLEL!$OMP DO do k=1,12 print*, k enddo!$OMP END DO!

49、$OMP END PARALLELend示例：線程0k=1,2,3線程1k=4,5,6線程2k=7,8,9線程2k=10,11,12!$OMP PARALLEL!$OMP DO!$OMP PARALLEL DO簡寫運行結果（屏幕截圖）運行結果：123789456101112線程0線程2線程1線程358 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),

50、z(N)!$ time1=OMP_get_wtime()!$OMP PARALLEL 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屏幕截圖采用單線程執行： 耗時2.15秒采用2線程執行：耗時 1.43秒采用

51、4線程執行： 耗時1.28秒59三、 OpenMP的數據結構： 共享與私有!$OMP PARALLEL DO do k=1,6 print*, k enddo!$OMP END PARALLEL DOend線程0k線程1k循環變量k 在兩個線程中的值是不同的；K是一個進程私有變量（PRIVATE)共享變量： 全體進程均可訪問的公共變量私有變量：各個進程私有的變量x=8.0 ; y=x+2.0; . !$OMP PARALLEL DO SHARED(x,y) PRIVATE (k, z) do k=1,6 z=k*x+y print*, x, y, z enddo!$OMP END PARALL

52、EL DOend線程0k , z線程1k, zx, y 私有變量公共變量60例： 將下面代碼并行化Integer, parameter: N=1024Real,dimension(N): x,y,zReal r. （給x, y賦值） Do k=1,N r=sqrt(x(k)*x(k)+y(k)*y(k) z(k)=1.0/(1.0+r) Enddo關鍵： 分析哪些是共享變量，哪些是私有變量。 顯然： r,k 是私有變量，其他均為共享變量!$OMP PARALLEL DO SHARED(DEFAULT) PREATE (r,k) Do k=1,N r=sqrt(x(k)*x(k)+y(k)*y(

53、k) z(k)=1.0/(1.0+r) Enddo!$OMP END PARALLEL DO61四、 OpenCFD-EC 3D 的OpenMP并行化舉例!$OMP PARALLEL DO do k=-1,nz+1 do j=-1,ny+1 do i=-1,nx+1 d(i,j,k)= B%U(1,i,j,k) uu(i,j,k)=B%U(2,i,j,k)/d(i,j,k) v(i,j,k)= B%U(3,i,j,k)/d(i,j,k) w(i,j,k)= B%U(4,i,j,k)/d(i,j,k) T(i,j,k)=(B%U(5,i,j,k)-0.5d0*d(i,j,k)*(uu(i,j,k

54、)*uu(i,j,k)+v(i,j,k)*v(i,j,k) & +w(i,j,k)*w(i,j,k)/(Cv*d(i,j,k) vt(i,j,k)=B%U(6,i,j,k) . . enddo enddo enddo!$OMP END PARALLEL DO62!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,k,s1x,s1y,s1z,xi,yi,zi,xj,yj,zj,xk,yk,zk,Jac,ix,iy,iz,jx,jy,jz,kx,ky,kz, &!$ ui,vi,wi,Ti,uj,vj,wj,Tj,uk,vk,wk,Tk,ux,uy,uz,v

55、x,vy,vz,wx,wy,wz,Tx,Ty,Tz,s11,s12,s13,s22,s23,s33,u1,v1,w1,E1,E2,E3,&!$ mu0,k0,v0,vti,vtj,vtk,vtx,vty,vtz,vn1,vn2,vn0,vfi) do k=1,nz-1 do j=1,ny-1 do i=1,nx s1x=B%ni1(i,j,k); s1y=B%ni2(i,j,k) ; s1z= B%ni3(i,j,k) ! xi=B%xc(i,j,k)-B%xc(i-1,j,k) yi=B%yc(i,j,k)-B%yc(i-1,j,k) zi=B%zc(i,j,k)-B%zc(i-1,j,k)

56、 xj=(B%x(i,j+1,k)-B%x(i,j,k) +B%x(i,j+1,k+1)-B%x(i,j,k+1)* yj=(B%y(i,j+1,k)-B%y(i,j,k) +B%y(i,j+1,k+1)-B%y(i,j,k+1)*0.5d0 zj=(B%z(i,j+1,k)-B%z(i,j,k) +B%z(i,j+1,k+1)-B%z(i,j,k+1)*0.5d0 xk=(B%x(i,j,k+1)-B%x(i,j,k) +B%x(i,j+1,k+1)-B%x(i,j+1,k)*0.5d0 yk=(B%y(i,j,k+1)-B%y(i,j,k) +B%y(i,j+1,k+1)-B%y(i,j+

57、1,k)*0.5d0 zk=(B%z(i,j,k+1)-B%z(i,j,k) +B%z(i,j+1,k+1)-B%z(i,j+1,k)*0.5d0 Jac=1.d0/(xi*yj*zk+yi*zj*xk+zi*xj*yk-xi*zj*yk-yi*xj*zk-zi*yj*xk) ix=Jac*(yj*zk-zj*yk) iy=Jac*(zj*xk-xj*zk) iz=Jac*(xj*yk-yj*xk) jx=Jac*(yk*zi-zk*yi) jy=Jac*(zk*xi-xk*zi) jz=Jac*(xk*yi-yk*xi) kx=Jac*(yi*zj-zi*yj) ky=Jac*(zi*xj-

58、xi*zj) kz=Jac*(xi*yj-yi*xj) ui=uu(i,j,k)-uu(i-1,j,k) vi=v(i,j,k)-v(i-1,j,k) wi=w(i,j,k)-w(i-1,j,k) Ti=T(i,j,k)-T(i-1,j,k) vti=vt(i,j,k)-vt(i-1,j,k) uj=0.25d0*(uu(i,j+1,k)-uu(i,j-1,k)+uu(i-1,j+1,k)-uu(i-1,j-1,k) vj=0.25d0*(v(i,j+1,k)-v(i,j-1,k)+v(i-1,j+1,k)-v(i-1,j-1,k) wj=0.25d0*(w(i,j+1,k)-w(i,j-1,

59、k)+w(i-1,j+1,k)-w(i-1,j-1,k) Tj=0.25d0*(T(i,j+1,k)-T(i,j-1,k)+T(i-1,j+1,k)-T(i-1,j-1,k) vtj=0.25d0*(vt(i,j+1,k)-vt(i,j-1,k)+vt(i-1,j+1,k)-vt(i-1,j-1,k) 63計算流體力學課程2011 習題匯總1推導無量綱的Navier-Stokes方程1.2 對于一維Euler方程組 推導Jocabian矩陣 以及 中 的表達式。 要求： 給出具體推導過程，切忌從書上抄錄公式 （越詳細越好）642.1 公式推導（1） 一激波從左向右傳播。激波左側物理量為 ; 激

60、波右側壓力為 ， 試計算激波右側的速度 。（2） 有一扇膨脹波從左向右傳播。 膨脹波左側物理量為 ； 膨脹波右側壓力為 ， 試計算膨脹波右側的速度 。膨脹波要求： 務必寫出詳細的步驟推導（越詳細越好）。切忌照抄書上的公式。 652.2 如下Sod 激波管問題:求出理論解， 并分別畫出t=0.14時刻 的分布曲線。663.1 對如下單波方程 構建的差分格式如下：試利用Fourier方法，分析其穩定性674.1 構造高分辨率差分格式，并進行理論分析及數值實驗 針對單波方程： 對于空間導數，構造出一種不超過6點格式；并進行Fourier誤差分析，畫出kr,ki的曲線。 要求：精度不限； 網格基架點數