中科院计算流体力学最新讲义CFD1114讲MPI并行程序设计初步2

《中科院计算流体力学最新讲义CFD1114讲MPI并行程序设计初步2》由会员分享，可在线阅读，更多相关《中科院计算流体力学最新讲义CFD1114讲MPI并行程序设计初步2（79页珍藏版）》请在金锄头文库上搜索。

1、中科院计算流体力学最新讲义CFD1114讲MPI并行程序设计初步2Stillwatersrundeep.流静水深流静水深,人静心深人静心深Wherethereislife,thereishope。有生命必有希望。有生命必有希望服务器服务器/前端机前端机计算节点计算节点a.exea.exea.exeMPI 程序的运行原理：程序的运行原理： 服务器（前端机）编译服务器（前端机）编译 可执行代码可执行代码复制复制 N 份，份，每个节点运行一份每个节点运行一份 调用调用MPI库函数库函数 得到每个节点号得到每个节点号 my_id 根据根据my_id 不同，程序执行情况不同不同，程序执行情况不同 调用调

2、用MPI 库函数进行通讯库函数进行通讯MPI 编程的基本思想：编程的基本思想： 主从式，主从式，对等式对等式2Copyright by Li Xinliang重点：对等式程序设计重点：对等式程序设计知识回顾知识回顾Copyright by Li Xinliang3计算节点计算节点a.exea.exea.exea.exe对等式对等式设计设计“对等式对等式”程序设计思想程序设计思想如果我是其中一个进程；如果我是其中一个进程；我应当做我应当做完成我需要完成的任务完成我需要完成的任务站在其中一个进程的角度思考站在其中一个进程的角度思考基本的基本的MPI函数（函数（6个）个）MPI初始化初始化 MPI_

3、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)消息接收消息接收MPI_Recv(buf,count,datatype,source,tag,comm,status,ierr) 4Copyright by Li

4、 XinliangMPI的消息发送机制的消息发送机制 两步进行两步进行MPI_Send( A, ) 发送发送MPI_Recv( B, ) 接收接收 发送发送 变量变量A接收接收 到变量到变量B配合使用配合使用5Copyright by Li Xinliang阻塞发送阻塞发送开始开始结束结束消息成功发出消息成功发出缓冲区可释放缓冲区可释放阻塞接收阻塞接收开始开始结束结束消息成功接收消息成功接收缓冲区数据可使用缓冲区数据可使用一、一、 阻塞式通信与非阻塞式通信阻塞式通信与非阻塞式通信阻塞式发送与接收阻塞式发送与接收MPI_Send( A, )MPI_Recv( B , )6Copyright by

5、 Li Xinliang MPI_Send( ) MPI_Send( ) 返回后缓冲区可释放返回后缓冲区可释放 sum= sum= call MPI_Send(sum,) call MPI_Send(sum,) sum= sum= 变量可重复利用变量可重复利用 MPI_Recv() MPI_Recv() 返回后缓冲区数据可使用返回后缓冲区数据可使用Call MPI_Recv(sum1,)Call MPI_Recv(sum1,)Sum=sum0+sum1Sum=sum0+sum1 7Copyright by Li Xinliang非阻塞发送非阻塞发送启动发送启动发送立即返回立即返回计计算算通信完

6、成通信完成释放发送缓冲区释放发送缓冲区发发送送消消息息非阻塞接收非阻塞接收启动接收启动接收立即返回立即返回计计算算通信完成通信完成引用接收数据引用接收数据接接收收消消息息计算计算与与通信通信重叠重叠非阻塞消息发送与接收非阻塞消息发送与接收8Copyright by Li Xinliang非阻塞消息发送非阻塞消息发送MPI_ISend(buf,count,datatype,dest,tag,comm,request,ierr)In buf,count,datatype,dest,tag,commOut request,ierr Request (返回的非阻塞通信对象返回的非阻塞通信对象, 整数整

7、数)非阻塞消息接收非阻塞消息接收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 requestO

8、ut status， flag (logical型型)9Copyright by Li Xinliang非阻塞通信调用后立即返回，缓冲区不能非阻塞通信调用后立即返回，缓冲区不能立即使用立即使用Sum= 计算某变量计算某变量MPI_Isend(sum .) 发送该变量发送该变量 sum= 不能给变量重新赋值不能给变量重新赋值 （发送可能尚未完成）（发送可能尚未完成）MPI_Irecv(sum1, )sum=sum0+sum1 数据不能立即使用数据不能立即使用 （接收可能未完成）（接收可能未完成）MPI_Isend(sum, , request, )Call MPI_Wait(request,st

9、atus,ierr)Sum= MPI_Irecv(sum1, , request, )Call MPI_Wait(request,status,ierr)Sum=sum0+sum1 10Copyright by Li Xinliang利用通信与计算重叠技术提高效率利用通信与计算重叠技术提高效率例：例： 计算差分计算差分串行程序串行程序 real A(N,N),B(N,N),h.Do i=1,NB(I,1)=(A(I,2)-A(I,1)/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.*

10、h)EnddoEnddo0J=1,2,3 . N-1, Ni=1i=2 i=N11Copyright by Li Xinliang并行程序并行程序 以两个进程并行为例以两个进程并行为例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,ierr)Else call MPI_recv(A1,N,MPI_real,0,99

11、,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)12Copyright by Li XinliangIf(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)=(A

12、(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 空闲空闲13Copyright by Li Xinliang“内边界”通信与计算重叠通信与计算重叠real A(N,N/2),B(N,N/2),A1(N),hinteger myid,ierr, req1, req2,status()If(myid .eq. 0) then call MPI_IS

13、end(A(1,N/2),N,MPI_real,1,99,MPI_Comm_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/214Copyright by Li XinliangDo j=2,N-1Do i=

14、1,NB(i,j)=(A(i,j+1)-A(i,j-1)/(2.*h)EnddoEnddoCall 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) 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特点：特点： 传递边界信息传递边界信息 同时同时进行计算进行计算内点内点读取系统时间读取系统时

15、间 doubleprecision time time=MPI_Wtime( ) 15Copyright by Li Xinliang二、二、 如何收发非连续数据如何收发非连续数据例如：例如： 发送数组的一行发送数组的一行A(100,50)发送发送 A(1,1),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).16Copyright by Li Xinliang方法方法2. 将发送的数据拷贝到连续的数组中将发送的数据拷贝到

16、连续的数组中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(1,i)=B(i) Enddoendif不足：不足： 额外的内存占用额外的内存占用 额外的拷贝操作额外的拷贝操作通信不复杂的情况，内存拷贝工作量不大，该方法也可以采用。通信不复杂的情况，内存拷贝工作量不大，该方法也可以采用。效

17、果还效果还可以可以17Copyright by Li Xinliang方法方法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

18、,1), A(2,1) , A(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个元素，跨度为个元素，

19、跨度为3（个元素）（个元素）Fortran 数组的一行数组的一行Real A(3,4). A(1,:)在内存中的排列次序18Copyright by Li Xinliang例：例： 发送三维数组中的一个面发送三维数组中的一个面 (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,ier

20、r) 通信通信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) , MPI_Type_Vector(P,M,M*N,MPI_Real,My_Type,ierr)通信通信3） 连续分布，无需构造新类型连续分布，无需构造新类型 19Copyright by Li XinliangMPI_TYPE_INDEXED(count, array_of_blocklengths, array_of_displacements, oldtype,newtype，ierr)构造数据类

21、型更灵活的函数构造数据类型更灵活的函数 直接指定每块的元素个数及偏移量直接指定每块的元素个数及偏移量块的数量块的数量（整数）（整数）每块元素的个数每块元素的个数（整形数组）（整形数组）每块的偏移量每块的偏移量（整形数组）（整形数组）例：例： 数组数组 real A(N,N), 欲将其欲将其上三角元素上三角元素作为消息发送，试构造其数据类型作为消息发送，试构造其数据类型 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)

22、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)Integer: Newtype,ierr count=N do k=1,N

23、 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, )20Copyright by Li XinliangN三、三、 MPI的通信域和组的通信域和组1.预定义通讯域预定义通讯域 MPI_Comm_World : 包含所有进程的组包含所有进程的组2.通讯

24、域的分割通讯域的分割 MPI_Comm_Split(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,

25、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, myid_line,ierr)MPI_Comm_WorldRAWColumnColor, 分组标准Key, 排序依据如相同，按原ID排 提交新定义的组提交新定义的组（否则新组无效，不要忘记）（否则新组无效，不要忘记）计算行号、列号21Copyright by Li Xinliang例：例： 计算差分计算差分 三维分割三维分割 A(

26、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 得到组内编号得到组内编号 3)建立三个方向通讯的数据结构建立三个方向通讯的数据结构4） 通信通信 ， 计算内点差分计算内点差分5） 计算边界差分计算边界差分02143576891011MPI_Comm_World22Copyright by Li XinliangParameter(M1=M/NM,N1=N/NN,P1=P/NP)Real A(

27、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,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/(

28、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)定义三个方向的通信域定义三个方向的通信域23Copyright by Li XinliangCall 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,ie

29、rr)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,TYPE_X, id_X_Pre, 99,Comm_X,request(1),ierr) Call MPI_Isend(A(M1,0,0),1,TYPE_X,id_

30、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) 定义新的数据结构定义新的数据结构24Copyright by Li Xinliang 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,

31、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 .内点内点边界点边界点25Copyright by Li Xinliang四、分布数组的文件存储四、分布

32、数组的文件存储 分布数组分布数组 real A(M/m1,N/n1)real A(M/m1,N/n1) 存储方式存储方式1. 每个进程存储到独立的文件每个进程存储到独立的文件 real A(M/m1,N/n1) character(len=50) filename write(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 优点：程序简单优点：程序简单缺点：缺点： 数

33、据文件多，不易处理；数据文件多，不易处理； 改变处理器数目时需特殊处理改变处理器数目时需特殊处理 012326Copyright by Li Xinliang 分布数组分布数组 real A(M/m1,N/n1)real A(M/m1,N/n1) 存储方式存储方式2： 收集到收集到0节点存储节点存储 存储到存储到 一个文件一个文件 缺点：缺点： 改变处理器规模时，需要处理改变处理器规模时，需要处理存储方式存储方式3： 收集到收集到0节点，重新装配成大数组节点，重新装配成大数组 收集收集 A(M/m1,N/n1) 组成组成 A0(M,N) real A0(M,N), A(M/m1,N/n1),

34、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 01230123027Copyright by Li Xinliang存储方式存储方式4. 按列搜集后存储按列搜集后存储 Real Aj(M)If( myid .eq. 0) then open(33,file=“A.dat

35、”,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)28Copyright by Li Xinliang存储方式存储方式5 并行并行IO (MPI

36、 2.0) 打开文件：打开文件： MPI_file_open(Comm,filename,mode,info,fileno,ierr) mode 打开类型：打开类型： MPI_Mode_RDONLY, MPI_Mode_RDWR, fileno 文件号，文件号， info 整数整数 （信息）（信息） 关闭文件关闭文件 ： MPI_file_close(fileno,ierr) 指定偏移位置读写指定偏移位置读写 MPI_file_read_at(fileno,offset,buff,const,datatype,status,ierr) MPI_file_write_at(fileno,offs

37、et,buff,const,datatype,status,ierr) offset 偏移，偏移， buff 缓冲区，缓冲区，const 数目数目 29Copyright by Li XinliangPart 3 实例教学实例教学 CFD程序的程序的MPI实现实现实例实例 （1） 用拟谱方法求解不可压用拟谱方法求解不可压N-S方程方程 实例（实例（2） 用流水线方法计算紧致差分用流水线方法计算紧致差分 常用的优化方法常用的优化方法30Copyright by Li Xinliang回顾回顾 基本的基本的MPI函数（函数（6个）个）MPI初始化初始化 MPI_Init(ierr) ； MPI结束

38、结束 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)消息接收消息接收MPI_Recv(buf,count,datatype,source,tag,comm,status,ierr) 31Copyright by Li Xinliang非阻塞消息发送非

39、阻塞消息发送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,tag,comm,request,ierr)In buf,count,datatype,source,tag,commOut request,ierr非阻塞通信的完成非阻塞通信的完成 MPI_Wait

40、(request,status,ierr) 等待消息收发完成等待消息收发完成 MPI_Test(request, flag,stutus,ierr) MPI_Waitall(const,request_array,status,ierr) 等待多个消息完成等待多个消息完成 In requestOut status， flag (logical型型)32Copyright by Li Xinliang发送非连续数据发送非连续数据构建新的数据结构构建新的数据结构MPI_TYPE_VECTOR(count,blocklength,stride,oldtype,newtype,ierr)Count:

41、块的数量；块的数量； 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).33Copyright by Li Xinliang

42、通讯域的分割通讯域的分割 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

43、)Call MPI_Comm_rank(Comm_line, myid_line,ierr)MPI_Comm_World34Copyright by Li Xinliang实例实例 1. 用（拟）谱方法求解二维不可压用（拟）谱方法求解二维不可压N-S方程方程2p物理模型物理模型周期性边界条件周期性边界条件按照给定能谱布置初始流动按照给定能谱布置初始流动 研究流动的演化规律研究流动的演化规律35Copyright by Li XinliangFourier 变换变换 (1D)Fourier 变换变换 的特点：的特点： 求导数求导数 - 乘积乘积困难：困难： 非线性项非线性项卷积卷积计算量巨大计算

44、量巨大在物理空间计算在物理空间计算Fourier 变换的快速算法变换的快速算法FFT36Copyright by Li Xinliang二维 Fourier 变换变换两次一维两次一维 Fourier 变换变换37Copyright by Li Xinliang求解步骤：求解步骤： 1） 读入初值读入初值 2） 调用调用FFT 得到得到 3） 计算计算 调用调用FFT 得到得到 4） 计算计算 调用调用FFT 得到得到 5） 计算计算 6） 积分积分 求出下一时间步的值求出下一时间步的值 7） 调用调用 FFT 得到得到 8） 循环循环 3）-7） 直到给定的时间直到给定的时间 38Copyri

45、ght by Li Xinliang实际计算中，要采用抑制混淆误差的措施程序的并行化：程序的并行化： 二维二维 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(:

46、)=Fu(:,j) call FFT1d(nx,u1) u(:,j)=u1(:) enddo end39Copyright by Li Xinliang数据重分布的实现数据重分布的实现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) 成成 A

47、240Copyright by Li Xinliang程序：程序：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) ) call MPI_Send(B1,M*N/(P*P),MPI_Real, k-1, .) Enddodo k=1,P call MPI_Recv(B2,M*N

48、/(P*P),MPI_Real, k-1, .) A2(k-1)*M/P+1: k*M/P) , : )=B2(:,:,P) Enddoend 问题：问题： 全部发送，全部发送， 发送成功发送成功后再启动接收。后再启动接收。 容易死锁容易死锁 按行分布按行分布 - 按列分布按列分布41Copyright by Li XinliangSubroutine 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

49、,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, .) A2(k-1)*M/P+1: k*M/P) , : )=B2(:,:,P) Enddoend 问题：问题： 按顺序发送、接收，不易死锁按顺序发送、接收，不易死锁42Copyright by Li Xinlia

50、ng数据全交换：数据全交换：MPI_AlltoAll(sendbuf,sendcount,sendtype,recvbuf,recvcount,recvtype,comm,ierr) sendbuf 发送缓冲区（首地址）发送缓冲区（首地址） recvbuf 接收缓冲区（首地址）接收缓冲区（首地址） sendcount 发送数目发送数目 recvcount 接收数目接收数目 sendtype 发送类型发送类型 recvtype 接收类型接收类型 Comm 通信域通信域 ierr 整数，返回错误值（整数，返回错误值（0为成功）为成功） To 0To 1To 2To 3Sendbuf 的数据格式的数

51、据格式sendcountFrom 0From 1From 2From 3Recvbuf 的数据格式的数据格式recvcount43Copyright by Li Xinliang程序：程序：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_Al

52、ltoAll (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) Enddoend 问题：问题： 无法做到计算与通信重叠无法做到计算与通信重叠 44Copyright by Li Xinliang 二维二维 并行并行FFT 的实现的实现 （输入数据、输出数据均为按列分布）（输入数据、输出数据均为按列分布）1） 调用一维调用一维FFT实现实现 i- 方向的变换方向的变换 u - u12） 重新分布数据重新分布数据 （按

53、列（按列- 按行）按行） u1 - u23) 调用一维调用一维FFT 实现实现j- 方向的变换方向的变换 u2- Fu24) 重新分布数据重新分布数据 （按行（按行 - 按列）按列） Fu2- Fu45Copyright by Li Xinliang实例实例 （2） 利用流水线利用流水线 实现紧致差分的并行化实现紧致差分的并行化紧致型差分格式：紧致型差分格式： 相同网格点上引入更多信息。相同网格点上引入更多信息。 性能更优化。性能更优化。 是是 的差分逼近的差分逼近普通差分格式：普通差分格式： 显式给出显式给出 Fi 的表达式的表达式紧致型差分格式：紧致型差分格式： 隐式给出隐式给出 Fi 的

54、表达式的表达式 6 阶中心6 阶对称紧致阶对称紧致 (Lele)5 阶迎风紧致阶迎风紧致 (Fu) j-2 j-1 j j+1 j+246Copyright by Li Xinliang 普通差分格式：普通差分格式： 直接计算导数，并行容易直接计算导数，并行容易紧致格式的计算：紧致格式的计算： 递推递推递推公式：递推公式：1)计算出计算出 （由边界条件或边界格式给出）（由边界条件或边界格式给出）2） 由由 递推计算递推计算 出全部导数出全部导数 后面的数据必须等待前一步计算完成，无法并行后面的数据必须等待前一步计算完成，无法并行47Copyright by Li Xinliang二维问题：二维

55、问题： 流水线法求解流水线法求解 流水线示意图流水线示意图步骤：步骤： 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,i); 如果如果myid P-1 向向 myid+1 发送发送 F(k,N1) 缺点：缺点： 通信次数过多通信次数过多48Copyright by Li Xinliang通信次数过于频繁通信次数过于频繁解决方法：解决方法： 分块流水线分块流水线步骤：步骤： 1) 计算计算 d(:,:) 2) for kp=1

56、,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 块块 49Copyright by Li Xinliang对称紧致格式对称紧致格式追赶法追赶法令则代入（代入（1） 得得对比（对比（2） 得得边界处导数可由边界条件或边界格式给出边界处导数可由边界条件或边界格式给出： 则步骤：步骤： 1） 2） 由由 （3）式递推，得到）式递推，得到 3) 4) 由由 （2

57、）式递推，得到）式递推，得到特点：特点： 两次递推。两次递推。 并行方法与前文类似并行方法与前文类似50Copyright by Li Xinliang3.常用的并行优化方法常用的并行优化方法 1） 通信与计算重叠通信与计算重叠 采用非阻塞通信采用非阻塞通信 Isend, Irecv 2） 用重复计算代替通信用重复计算代替通信 3） 拆分长消息、合并短消息拆分长消息、合并短消息 4） 优化通信方式优化通信方式 51Copyright by Li Xinliang 用重复计算代替通信用重复计算代替通信 例如：例如： 计算差分计算差分 u 分布存储分布存储, f(u) 为为 u 的函数的函数 01

58、方法方法 1） 计算出计算出 v=f(u) 通信得到通信得到 uN+1, vN+1 计算差分计算差分方法方法 2） 计算出计算出 v=f(u) 通信得到通信得到 uN+1 （边界外）（边界外） 计算出计算出 vN+1 =f(uN+1) 计算差分计算差分 方法方法 2） 计算量大，通信量小计算量大，通信量小 当函数当函数 f(u)不复杂时，可提高效率不复杂时，可提高效率1, 2 N N+152Copyright by Li Xinliang 长消息切割成多个短消息发送、接收长消息切割成多个短消息发送、接收 call MPI_Send(A(1),100000, MPI_Real, 1, ) 改为：

59、改为： do m=1,10 call MPI_Send(A(m-1)*10000+1),10000,MPI_real,1 ) enddo 长消息：长消息： 非缓冲；非缓冲； 短消息：短消息： 缓冲缓冲 缓冲区MPI_Send缓冲区MPI_SendMPI_RecvMPI_Recv53Copyright by Li Xinliang合并短消息合并短消息 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, MPI_Real, 1,

60、) 54Copyright by Li Xinliang 优化通信方式优化通信方式 例：例： 数据散发数据散发0号号 进程：进程： 数据数据 A(100)， 散发给散发给 0-99方式方式1） 0 进程执行进程执行100次次 MPI_Send 其他进程执行其他进程执行 MPI_Recv MPI_Scatter() 采用该算法采用该算法方式方式 2） 0 进程进程 把把 A(100) 切割成切割成10份份 ， 发送给发送给10个进程个进程 10个进程接收个进程接收A1(10) 后再散发后再散发 55Copyright by Li XinliangOpenMP并行编程入门并行编程入门一、一、 特点

61、特点 1. 针对针对共享内存共享内存计算机结构计算机结构 全部全部CPU/线程均可访问内存线程均可访问内存 2. 程序改动量小、实现方便程序改动量小、实现方便 （以编译指示符为主）（以编译指示符为主） 3. 适用于小规模并行或与适用于小规模并行或与MPI配合配合 进行大规模并行进行大规模并行 内存CPU（核心）CPU（核心）CPU（核心）1台PC机 / 1个计算节点 （共享内存构架）CPU内存CPU内存CPU内存外部网络节点1节点2Cluster结构， 分布内存构架 print*, code 1!$OMP PARALLEL print*, code 2!$OMP END PARALLEL pr

62、int*, code 3 end例1 (test1.f90) ：编译编译 （在深腾（在深腾7000)运行结果（屏幕截图）ifort test1.f90 -openmp添加 openmp 选项运行： 1. 设置线程数（并行执行的数目） export OMP_NUM_THREADS=4 （例如，4个） 2. 执行： ./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*

63、, “code 3”!$OMP END PARALLEL!$OMP END PARALLEL print*, code 4 end3.DO 循环分解循环分解 （openMP最常用的并行方法）最常用的并行方法）!$OMP PARALLEL!$OMP DO do k=1,12 print*, k enddo!$OMP END DO!$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简写运行结果（屏幕截图）运行结果：12378

64、9456101112线程0线程2线程1线程3 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 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

65、)=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秒秒采用采用4线程执行：线程执行： 耗时耗时1.28秒秒三、三、 OpenMP的数据结构：的数据结构： 共享与私有共享与私有!$OMP PARALLEL DO do k=1,6 print*, k enddo!$OMP EN

66、D 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 PARALLEL DOend线程线程0k , z线程线程1k, zx, y 私有变量公共变量例： 将下面代码并行化

67、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(k) z(k)=1.0

68、/(1.0+r) Enddo!$OMP END PARALLEL DO四、四、 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)*uu(

69、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 DO!$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,vx,vy,vz

70、,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) xj=(B%

71、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+1,k)*0.

72、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-xi*zj)

73、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,k)+w(i-

74、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) 64Copyright by Li Xinliang计算流体力学计算流体力学课程课程2011 习题汇总习题汇总1.1推导无量纲的推导无量纲的Navier-Stokes方程方程1.2 对于一维对于一维Euler方程组方程组 推导推导Jocabian矩阵矩阵 以及以及 中中 的表达式。的表达式。 要求：要求： 给

75、出具体推导过程，切忌从书上抄录公式给出具体推导过程，切忌从书上抄录公式 （越详细越好）（越详细越好）Copyright by Li Xinliang652.1 公式推导公式推导（1） 一激波从左向右传播。激波左侧物理量为一激波从左向右传播。激波左侧物理量为 ; 激波右侧压力为激波右侧压力为 ， 试计算激波右试计算激波右侧的速度侧的速度 。（2） 有一扇膨胀波从左向右传播。有一扇膨胀波从左向右传播。 膨胀波左侧物膨胀波左侧物理量为理量为 ； 膨胀波右侧压力为膨胀波右侧压力为 ， 试计算试计算膨胀波右侧的速度膨胀波右侧的速度 。膨胀波要求：要求： 务必写出详细的步骤推导（越详细越好）。切忌照抄书务

76、必写出详细的步骤推导（越详细越好）。切忌照抄书上的公式。上的公式。 2.2 如下如下Sod 激波管问题激波管问题:求出理论解，求出理论解， 并分别画出并分别画出t=0.14时刻时刻 的分布曲线。的分布曲线。Copyright by Li Xinliang663.1 对如下单波方程对如下单波方程 构建的差分格式如下：构建的差分格式如下：试利用试利用Fourier方法，分析其稳定性方法，分析其稳定性4.1 构造高分辨率差分格式，并进行理论分析及数值实验构造高分辨率差分格式，并进行理论分析及数值实验 针对单波方程针对单波方程： 对于空间导数，构造出一种不超过对于空间导数，构造出一种不超过6点格式；并

77、进行点格式；并进行Fourier误差分析，误差分析，画出画出kr,ki的曲线。的曲线。 要求：精度不限；要求：精度不限； 网格基架点数不超过网格基架点数不超过6个；个； 能够分辨的波数范围尽量宽；能够分辨的波数范围尽量宽； （即（即kr,ki曲线近可能接近准确解）曲线近可能接近准确解） 给出差分的具体表达式，给出差分的具体表达式， 画出画出kr,ki的曲线；的曲线； 说明构造格式的阶数，并采用本说明构造格式的阶数，并采用本PPT第第5页的方法给出的精度验证；页的方法给出的精度验证； 形如：另外，进行如下数值验证：另外，进行如下数值验证：空间采用空间采用20个网格点，采用新构造的差分格式离散；时

78、间推进采用个网格点，采用新构造的差分格式离散；时间推进采用3步步Runge-Kutta方法，时间步长可足够小（例如方法，时间步长可足够小（例如0.01）。给出）。给出t=20,50两个时两个时刻的数值解，与精确解比较（画图），并给出数值解的刻的数值解，与精确解比较（画图），并给出数值解的L2模误差。模误差。68Copyright by Li Xinliang提示： 1. 如不使用优化技术，则格式构造方法简单， Taylor展开后解代数方程组即可。 2. 建议尝试使用优化技术 例： 假设格式形式如下如果要求其有5阶精度，则通过Taylor展开可得到6个方程，6个系数可直接解出。我们要求其有4阶精

79、度（当然3阶，2阶也可），于是Taylor展开只能提供5个方程。6个未知数（a1-a6）， 5个方程； 有1个自由参数。 调整这个自由参数，使得kr,ki曲线最为理想。 如何调整？ 1) 可以人工调整，观察kr,ki曲线，选取满意的。 2）可自动调整，设立一个优化目标函数。 例如 调整自由参数，使得该目标函数取最大值。思路：牺牲精度，提高分辨率 69Copyright by Li Xinliang4.2 数值求解数值求解Sod 激波管问题激波管问题 计算其数值解，画出计算其数值解，画出t=0.14时刻密度、速度及压力的分布；并与精确解进行比较时刻密度、速度及压力的分布；并与精确解进行比较（要求

80、画在一张图上）。（要求画在一张图上）。 要求：要求： 1） 空间网格数空间网格数100， 时间推进格式选用时间推进格式选用3阶阶Runge-Kutta,时间步长自选。时间步长自选。 2） 可选用逐点分裂，也可选用特征分裂。可选用逐点分裂，也可选用特征分裂。 3） 建议采用本讲作业题建议采用本讲作业题2（或作业题（或作业题1）自行构造的差分格式计算。）自行构造的差分格式计算。 （作业题（作业题2是激波捕捉格式，效果应当会好些）。是激波捕捉格式，效果应当会好些）。 如果作业题如果作业题1和作业题和作业题2遇到困难，也可采用现有的差分格式遇到困难，也可采用现有的差分格式。70Copyright by

81、 Li Xinliang针对如下针对如下Sod 激波管问题激波管问题 用用5阶阶WENO格式格式计算其数值解，画出计算其数值解，画出t=0.14时刻密度、速度及压力的分布；并与精时刻密度、速度及压力的分布；并与精确解进行比较（要求数值解与精确解画在同一张图上，便于比较）。确解进行比较（要求数值解与精确解画在同一张图上，便于比较）。 要求：要求： 空间网格数空间网格数100， 时间推进格式选用时间推进格式选用3阶阶Runge-Kutta,时间步长自选。时间步长自选。 71Copyright by Li Xinliang作业作业 5.1 针对如下针对如下Shu-Osher 激波激波-密度扰动波干扰

82、问题：密度扰动波干扰问题： 用用5阶阶WENO格式格式计算其数值解，画出计算其数值解，画出t=0.1时刻密度、速度及压力的分布时刻密度、速度及压力的分布 要求：要求： 1） 空间网格数空间网格数200， 时间推进格式选用时间推进格式选用3阶阶Runge-Kutta,时间步长自选。时间步长自选。 2） 使用使用2000个网格点计算，其结果作为个网格点计算，其结果作为“精确解精确解”，与其它结果画在一起，与其它结果画在一起，便于比较。便于比较。72Copyright by Li Xinliang作业作业 5.2初值为:Copyright by Li Xinliang73作业作业 7.1 编制有限体

83、积程序编制有限体积程序, 并计算图示并计算图示钝楔绕流问题钝楔绕流问题r=1几何及流动参数：几何及流动参数： 钝楔的半楔角为钝楔的半楔角为5。来。来流流Mach数数6， Reynold数数10000 （以头半（以头半径及来流参数度量）。径及来流参数度量）。 来流攻角为来流攻角为0，壁面，壁面温度为来流温度的温度为来流温度的4.2倍。倍。计算中，粘性系数可用Sutherland公式计算：本计算中，来流温度设定为要求：要求： 编写程序，给出壁面压力、温度的二维分布云图，并给编写程序，给出壁面压力、温度的二维分布云图，并给出壁面压力及热流。出壁面压力及热流。网格可从“流体中文网”下载 （tecplo

84、t格式）Copyright by Li Xinliang74习题习题 9.1 网格生成网格生成 通过解椭圆型方程生成通过解椭圆型方程生成NACA0012翼型的网格翼型的网格要求：要求： 推荐采用图示的推荐采用图示的C型网格，型网格， 网格点不限；网格点不限； 外边界的位置不限；外边界的位置不限； 可求解无源项的可求解无源项的Laplace方程或有源项的方程或有源项的Poisson方程；方程； 绘制出网格，并给出具体计算说明及公式。绘制出网格，并给出具体计算说明及公式。“翼型数据库大全”http:/www.uiuc.edu/ph/www/m-selig/ads.htmlhttp:/ 0012 翼

85、型（对称翼型）的拟合曲线为翼型（对称翼型）的拟合曲线为 （宋宇宁等（宋宇宁等 “微型飞行器的翼型拟合与模具加工微型飞行器的翼型拟合与模具加工”，电加工与模具，电加工与模具，2002第第5期，期，33-36）习题习题 10.1 求解方腔问题求解方腔问题问题描述：问题描述： 如图示边长为如图示边长为L的方腔，上表面流体以常速度的方腔，上表面流体以常速度U运运动，求解里面的流场（假设流动定常）。动，求解里面的流场（假设流动定常）。 考虑考虑 三种情况三种情况要求：要求： 数值方法不限数值方法不限 （人工压缩性方法、投影法、涡量（人工压缩性方法、投影法、涡量-流函数方法及流函数方法及SIMPLE方方法

86、均可）；法均可）； 空间离散采用差分法，建议采用较高阶精度的方法。空间离散采用差分法，建议采用较高阶精度的方法。 绘制出定常解的流线图。绘制出定常解的流线图。 请详细写明方程及公式的推导过程及计算流程，切勿只上交计算结果。请详细写明方程及公式的推导过程及计算流程，切勿只上交计算结果。 Copyright by Li Xinliang76作业题作业题 11.1以不可压缩槽道湍流为例，推导线性稳定性理论的控制方程及边界条件。要求：以不可压缩槽道湍流为例，推导线性稳定性理论的控制方程及边界条件。要求： 1） 给出扰动量给出扰动量 满足的满足的线性化线性化控制方程及边界条件（必须有推导步骤）控制方程及

87、边界条件（必须有推导步骤） 2） 推导扰动量振幅满足的推导扰动量振幅满足的O-S方程及边界条件方程及边界条件 （必须有详细推导步骤）（必须有详细推导步骤） Copyright by Li Xinliang77作业作业 12.1 试推导不可压缩湍动能试推导不可压缩湍动能 k 及湍能耗散率及湍能耗散率 e e 所满足的控制方程。所满足的控制方程。其中：其中：要求：要求： 必须给出详细的推导过程，切勿只照抄最终公式必须给出详细的推导过程，切勿只照抄最终公式参考文献： 是勋刚 湍流 第三篇提示：提示： step 1) 写出脉动量满足的方程写出脉动量满足的方程 step 2) 两端乘以两端乘以 并平均，

88、即可的并平均，即可的k满足的方程满足的方程 step 3) (1)式两端对式两端对 求导，乘以求导，乘以 后平均，可得后平均，可得e e方程方程Copyright by Li Xinliang78作业作业13.1 熟悉熟悉MPI环境及基本编程环境及基本编程 1） 建立建立MPI运行环境运行环境 （有并行机账户或在微机上安装（有并行机账户或在微机上安装MPI环境）。环境）。 2） 编制如下基本的编制如下基本的MPI程序程序 计算计算S=1+2+3+1000 要求程序可以实现要求程序可以实现N个进程的并行运行且负载尽量均衡。个进程的并行运行且负载尽量均衡。 N可变，程序中可变，程序中使用使用MPI

89、_Comm_Size()函数读入函数读入N。由。由0号进程打印计算结果。号进程打印计算结果。 3）在并行环境上运行，输出结果。）在并行环境上运行，输出结果。 要求：要求： 提交源程序及运行情况的屏幕截图提交源程序及运行情况的屏幕截图Copyright by Li Xinliang7913.2 实现矩阵相乘的并行计算实现矩阵相乘的并行计算矩阵矩阵A, B 均为均为N*N的方阵，试计算矩阵的方阵，试计算矩阵C=AB；使用使用P个进程并行计算（个进程并行计算（N可以被可以被P整除）；整除）；矩阵矩阵A，B及及C均采用分布式存储；均采用分布式存储；A, C按行分割，按行分割， B按列分割存储（见本稿按

90、列分割存储（见本稿 47页）。页）。要求编写计算要求编写计算C矩阵的矩阵的MPI程序，并进行计算。程序，并进行计算。实际计算时，矩阵实际计算时，矩阵A, B请采用如下值，请采用如下值， N设为设为100计算出计算出C矩阵后，请计算矩阵后，请计算 ，并由根节点打印出来。，并由根节点打印出来。将将S值与串行程序的结果进行对比，校验程序的正确性；值与串行程序的结果进行对比，校验程序的正确性；使用使用1,2,4,10个进程进行计算，并利用个进程进行计算，并利用MPI_Wtime( )函数计算程函数计算程序的运行时间；考核加速比及计算效率。序的运行时间；考核加速比及计算效率。要求：要求： 1）提交计算程序；）提交计算程序； 2）使用）使用1，2，4，10个进程计算，提交个进程计算，提交计算结果（计算结果（S值及计算时间）、计算效率及加速比。值及计算时间）、计算效率及加速比。