中科院計(jì)算流體力學(xué)最新講義CFD2011-第14講-MPI并行程序設(shè)計(jì)初步2

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

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

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

4、據(jù)可使用Call MPI_Recv(sum1,)Sum=sum0+sum1 7非阻塞發(fā)送啟動(dòng)發(fā)送立即返回計(jì)算通信完成釋放發(fā)送緩沖區(qū)發(fā)送消息非阻塞接收啟動(dòng)接收立即返回計(jì)算通信完成引用接收數(shù)據(jù)接收消息計(jì)算與通信重疊非阻塞消息發(fā)送與接收8非阻塞消息發(fā)送MPI_ISend(buf,count,datatype,dest,tag,comm,request,ierr)In buf,count,datatype,dest,tag,commOut request,ierr Request (返回的非阻塞通信對(duì)象, 整數(shù))非阻塞消息接收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) 等待消息收發(fā)完成 MPI_Test(request, flag,stutus,ierr) MPI_Waitall(const,request_array,status,ierr) 等待多個(gè)消息完成 In requestOut status， flag (logical型)9非阻塞通信調(diào)用后立即返回，緩沖區(qū)不能立即使用Sum= 計(jì)算某變量MPI_Isend(sum .) 發(fā)送

6、該變量 sum= 不能給變量重新賦值 （發(fā)送可能尚未完成）MPI_Irecv(sum1, )sum=sum0+sum1 數(shù)據(jù)不能立即使用 （接收可能未完成）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利用通信與計(jì)算重疊技術(shù)提高效率例： 計(jì)算差分串行程序 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并行程序 以兩個(gè)進(jìn)程并行為例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特點(diǎn)： 先收發(fā)邊界信息 再進(jìn)行計(jì)算缺點(diǎn)： 通信過(guò)程中CPU 空閑“內(nèi)邊界”13通信與計(jì)算重疊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特點(diǎn)： 傳遞邊界信息 同時(shí)進(jìn)行計(jì)算內(nèi)點(diǎn)讀取系統(tǒng)時(shí)間 doubleprecision time time=MPI_Wtime( ) 15二、 如何收發(fā)非連續(xù)數(shù)據(jù)例如： 發(fā)送數(shù)組的一行A(100,50)發(fā)送 A(1,1)

12、,A(1,2) ,A(1,3)A(1,1), A(1,2), A(1,3) 方法1. 多次發(fā)送 通信開銷大、效率低A(1,1), A(2,1), A(1,2), A(2,2) . A(1,3).16方法2. 將發(fā)送的數(shù)據(jù)拷貝到連續(xù)的數(shù)組中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不足： 額外的內(nèi)存占用 額外的拷貝操作通信不復(fù)雜的情況，內(nèi)存拷貝工作量不大，該方法也可以采用。效果還可以17方法3： 構(gòu)建新的數(shù)據(jù)結(jié)構(gòu) Count: 塊的數(shù)量； blocklength: 每塊的元素個(gè)數(shù)Stride: 跨度 （各塊起始元素之間的距離）Oldtype: 舊數(shù)據(jù)類型， Newtype: 新數(shù)據(jù)類型 （整數(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

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固定間隔（跨度）的非連續(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個(gè)元素，跨度為3（個(gè)元素）Fortran 數(shù)組的一行Real A(3,4). A(1,:

15、)在內(nèi)存中的排列次序18例： 發(fā)送三維數(shù)組中的一個(gè)面 (Fortran) 數(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,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） 連續(xù)分布，無(wú)需構(gòu)造新類型 19MPI_TYPE_INDEXED(count, array_of_blocklengths, array_of_displacements, oldtype,newtype，ierr)構(gòu)造數(shù)據(jù)類型更靈活的函數(shù) 直接指定每塊的元素個(gè)數(shù)及偏移量塊的數(shù)量（整數(shù)）每塊元素的個(gè)數(shù)（整形數(shù)組）每塊的偏移量（整形數(shù)組）例： 數(shù)組 real A(N,N), 欲將其上三角元素作為消息發(fā)送，試構(gòu)造其數(shù)據(jù)類型 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)內(nèi)存中的存儲(chǔ)次序(Fortran)N列N行注意： Fortran 行優(yōu)先次序存儲(chǔ)； C為列優(yōu)先次序存儲(chǔ)觀察規(guī)律： N塊； 第k塊有k個(gè)元素；第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的通信域和組預(yù)定義通訊域 MPI_Comm_World : 包含所有進(jìn)程的組通訊域的分割 MPI_Comm_S

19、plit(comm,color, key,New_Comm ) 02143576891011Color 相同的進(jìn)程在同一組根據(jù)key的大小排序 （key相同時(shí)按原ID排序）例如： 12個(gè)進(jìn)程， 分成 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, 分組標(biāo)準(zhǔn)Key, 排序依據(jù)如相同，按原ID排 提交新定義的組（否則新組無(wú)效，不要忘記）計(jì)算行號(hào)、列號(hào)21例： 計(jì)算差分 三維分割 A(M1,N1,P1)(M1=M/NM, N1=N/NN, P1=P/NP)基本思路：1） “擴(kuò)大”的數(shù)組 A(0: M1+1, 0: N1+1,0:P1+1)2）分割成三個(gè)組 Comm_X, Comm_Y, Comm_Z

21、得到組內(nèi)編號(hào) 建立三個(gè)方向通訊的數(shù)據(jù)結(jié)構(gòu)4） 通信 ， 計(jì)算內(nèi)點(diǎn)差分5） 計(jì)算邊界差分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)定義三個(gè)方向的通信域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) 定義新的數(shù)據(jù)結(jié)構(gòu)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 .內(nèi)點(diǎn)邊

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

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

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

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

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

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) 等待消息收發(fā)完成 MPI_Test(request, flag,stutus,ierr) MPI_Waitall(const,request_array,status,ierr) 等待多個(gè)消息完成 In requestOut status， flag (logical型)32發(fā)送非連續(xù)數(shù)據(jù)構(gòu)建新的數(shù)

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

33、).33通訊域的分割 MPI_Comm_Split(comm,color, key,New_Comm ) 02143576891011Color 相同的進(jìn)程在同一組根據(jù)key的大小排序例如： 12個(gè)進(jìn)程， 分成 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實(shí)例 1. 用（擬）譜方法求解二維不可壓N-S方程2p物理模型周期性邊界條件按照給定能譜布置初始流動(dòng) 研究流動(dòng)的演化規(guī)律35Fourier 變換 (1D)Fourier 變換 的特點(diǎn)： 求導(dǎo)數(shù) - 乘積困難： 非線性項(xiàng)卷積計(jì)算量巨大在物理空間計(jì)算Fourier 變換的快速算法FFT36二維 Fourier 變換兩次一維 Fourier 變換37求解步驟： 1） 讀入初值 2） 調(diào)用FFT 得到 3） 計(jì)算 調(diào)用FFT 得到 4） 計(jì)算 調(diào)用FFT 得到 5） 計(jì)算 6） 積分 求出下一時(shí)間步的值 7） 調(diào)用 FFT 得到 8）

35、循環(huán) 3）-7） 直到給定的時(shí)間 實(shí)際計(jì)算中，要采用抑制混淆誤差的措施38程序的并行化： 二維 FFT二維FFT: 調(diào)用兩次一維FFT一維 FFT 算法復(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(:

36、) enddo end39數(shù)據(jù)重分布的實(shí)現(xiàn)A1(M/P,N)A2 (M,N/P)1234abcd對(duì)等式編程思想 “我”需要完成的工作 1) 將數(shù)據(jù) A1(M/P,N) 切割成P塊 ，存入數(shù)組B1(M/P, N/P,P) 2) 將數(shù)據(jù) B1(:,:,k) 發(fā)到 進(jìn)程 k (k=0,1.P-1) 3) 從進(jìn)程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 問(wèn)題： 全部發(fā)送， 發(fā)送成功后再啟動(dòng)接收。 容易死鎖 按行分布 - 按列分布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 問(wèn)題： 按順序發(fā)送、接收，不易死鎖42數(shù)據(jù)全交換：MPI_AlltoAll(sendbuf,sendcount,sendtype,recvbuf,recvcount,recvtype,comm,ierr) sendbuf 發(fā)送緩沖區(qū)（首地址） recvbuf 接收緩沖區(qū)（首地址） sendcount 發(fā)送數(shù)目 recvcount 接收數(shù)目 sendtype 發(fā)送類型 recvtype 接收類型 Comm 通信域 ierr 整數(shù)，返回錯(cuò)誤值（0為成功） To 0To 1To 2To 3Sendbuf

40、的數(shù)據(jù)格式sendcountFrom 0From 1From 2From 3Recvbuf 的數(shù)據(jù)格式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 問(wèn)題： 無(wú)法做到計(jì)算與通信重疊 44 二維 并行FFT 的實(shí)現(xiàn) （輸入數(shù)據(jù)、輸出數(shù)據(jù)均為按列分布）1） 調(diào)用一維FFT實(shí)現(xiàn) i- 方向的變換 u - u12） 重新分布數(shù)據(jù) （按列- 按行） u1 - u2 調(diào)用一維FFT 實(shí)現(xiàn)j- 方向的變換 u2- Fu2 重新分布數(shù)據(jù) （按行 - 按列） Fu2- Fu45實(shí)例 （2） 利用流水線 實(shí)現(xiàn)緊致差分的并行化緊致型差分格式： 相同網(wǎng)格點(diǎn)上引入更

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

43、 i=1,N1 (N1=N/P) 計(jì)算 F(k,i); 如果myid P-1 向 myid+1 發(fā)送 F(k,N1) 缺點(diǎn)： 通信次數(shù)過(guò)多48通信次數(shù)過(guò)于頻繁解決方法： 分塊流水線步驟： 1) 計(jì)算 d(:,:) 2) for kp=1,MP 如果 myid=0, 計(jì)算 F(kp,0), 否則 從myid-1接收 F(kp,0); for j=1,N1 (N1=N/P) 計(jì)算 F(kp,j); 如果myid P-1 向 myid+1 發(fā)送 F(kp,N1) F(kp,i) 表示第kp 塊 49對(duì)稱緊致格式追趕法令則代入（1） 得對(duì)比（2） 得邊界處導(dǎo)數(shù)可由邊界條件或邊界格式給出： 則步驟： 1

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

45、時(shí)，可提高效率1, 2 N N+152 長(zhǎng)消息切割成多個(gè)短消息發(fā)送、接收 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 長(zhǎng)消息： 非緩沖； 短消息： 緩沖 緩沖區(qū)MPI_Send緩沖區(qū)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 優(yōu)化通信方式 例： 數(shù)據(jù)散發(fā)0號(hào) 進(jìn)程： 數(shù)據(jù) A(100)， 散發(fā)給 0-99方式1） 0 進(jìn)程執(zhí)行100次 MPI_Send 其他進(jìn)程執(zhí)行 MPI_Recv MPI_Scatter() 采用該算法方式 2） 0 進(jìn)程 把 A(100) 切割成10份 ， 發(fā)送給10個(gè)進(jìn)程 10個(gè)進(jìn)程接收A1(10) 后再散發(fā) 55OpenMP并行編程入門一、 特點(diǎn) 1. 針對(duì)共享內(nèi)存計(jì)算機(jī)結(jié)構(gòu) 全部CPU/線程均可訪問(wèn)內(nèi)存 2. 程序改動(dòng)量小、實(shí)現(xiàn)方便 （以編譯指示符為主） 3. 適用于小規(guī)模并行或與MPI配合 進(jìn)行大規(guī)模并行 內(nèi)

47、存CPU（核心）CPU（核心）CPU（核心）1臺(tái)PC機(jī) / 1個(gè)計(jì)算節(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)2Cluster結(jié)構(gòu)， 分布內(nèi)存構(gòu)架56 print*, code 1!$OMP PARALLEL print*, code 2!$OMP END PARALLEL print*, code 3 end例1 (test1.f90) ：編譯 （在深騰7000)運(yùn)行結(jié)果（屏幕截圖）ifort test1.f90 -openmp添加 openmp 選項(xiàng)運(yùn)行： 1. 設(shè)置線程數(shù)（并行執(zhí)行的數(shù)目） export OMP_NUM_THREADS=4 （例如，4個(gè)） 2.

48、 執(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*, code 2“!$OMP PARALLEL print*, “code 3”!$OMP END PARALLEL!$OMP END PARALLEL print*, code 4 end57DO 循環(huán)分解 （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簡(jiǎn)寫運(yùn)行結(jié)果（屏幕截圖）運(yùn)行結(jié)果：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屏幕截圖采用單線程執(zhí)行： 耗時(shí)2.15秒采用2線程執(zhí)行：耗時(shí) 1.43秒采用

51、4線程執(zhí)行： 耗時(shí)1.28秒59三、 OpenMP的數(shù)據(jù)結(jié)構(gòu)： 共享與私有!$OMP PARALLEL DO do k=1,6 print*, k enddo!$OMP END PARALLEL DOend線程0k線程1k循環(huán)變量k 在兩個(gè)線程中的值是不同的；K是一個(gè)進(jìn)程私有變量（PRIVATE)共享變量： 全體進(jìn)程均可訪問(wèn)的公共變量私有變量：各個(gè)進(jìn)程私有的變量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關(guān)鍵： 分析哪些是共享變量，哪些是私有變量。 顯然： 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計(jì)算流體力學(xué)課程2011 習(xí)題匯總1推導(dǎo)無(wú)量綱的Navier-Stokes方程1.2 對(duì)于一維Euler方程組 推導(dǎo)Jocabian矩陣 以及 中 的表達(dá)式。 要求： 給出具體推導(dǎo)過(guò)程，切忌從書上抄錄公式 （越詳細(xì)越好）642.1 公式推導(dǎo)（1） 一激波從左向右傳播。激波左側(cè)物理量為 ; 激

60、波右側(cè)壓力為 ， 試計(jì)算激波右側(cè)的速度 。（2） 有一扇膨脹波從左向右傳播。 膨脹波左側(cè)物理量為 ； 膨脹波右側(cè)壓力為 ， 試計(jì)算膨脹波右側(cè)的速度 。膨脹波要求： 務(wù)必寫出詳細(xì)的步驟推導(dǎo)（越詳細(xì)越好）。切忌照抄書上的公式。 652.2 如下Sod 激波管問(wèn)題:求出理論解， 并分別畫出t=0.14時(shí)刻 的分布曲線。663.1 對(duì)如下單波方程 構(gòu)建的差分格式如下：試?yán)肍ourier方法，分析其穩(wěn)定性674.1 構(gòu)造高分辨率差分格式，并進(jìn)行理論分析及數(shù)值實(shí)驗(yàn) 針對(duì)單波方程： 對(duì)于空間導(dǎo)數(shù)，構(gòu)造出一種不超過(guò)6點(diǎn)格式；并進(jìn)行Fourier誤差分析，畫出kr,ki的曲線。 要求：精度不限； 網(wǎng)格基架點(diǎn)數(shù)