中科院计算流体力学最新讲义CFD2011-第14讲-MPI并行程序设计初步2.ppt

计算流体力学讲义第六讲MPI并行程序设计 2 李新亮lixl 力学所主楼219 82543801 知识点 阻塞通信与非阻塞通信非连续数据的发送与接收OpenMP并行程序设计初步 1 CopyrightbyLiXinliang 讲义 课件上传至 流体中文网 流体论坛 CFD基础理论 也可到如下网址下载 http cid MPI程序的运行原理 服务器 前端机 编译可执行代码复制N份 每个节点运行一份调用MPI库函数得到每个节点号my id根据my id不同 程序执行情况不同调用MPI库函数进行通讯 MPI编程的基本思想 主从式 对等式 2 CopyrightbyLiXinliang 重点 对等式程序设计 知识回顾 CopyrightbyLiXinliang 3 a exe 对等式设计 对等式 程序设计思想 如果我是其中一个进程 我应当做 完成我需要完成的任务 站在其中一个进程的角度思考 基本的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 消息接收MPI Recv buf count datatype source tag comm status ierr 4 CopyrightbyLiXinliang MPI的消息发送机制 两步进行MPI Send A 发送MPI Recv B 接收 发送变量A 接收到变量B 配合使用 5 CopyrightbyLiXinliang 一 阻塞式通信与非阻塞式通信 阻塞式发送与接收 MPI Send A MPI Recv B 6 CopyrightbyLiXinliang MPI Send 返回后缓冲区可释放sum callMPI Send sum sum 变量可重复利用MPI Recv 返回后缓冲区数据可使用CallMPI Recv sum1 Sum sum0 sum1 7 CopyrightbyLiXinliang 非阻塞发送 启动发送 立即返回 计算 通信完成 释放发送缓冲区 发送消息 非阻塞接收 启动接收 立即返回 计算 通信完成 引用接收数据 接收消息 计算与通信重叠 非阻塞消息发送与接收 8 CopyrightbyLiXinliang 非阻塞消息发送MPI ISend buf count datatype dest tag comm request ierr Inbuf count datatype dest tag commOutrequest ierrRequest 返回的非阻塞通信对象 整数 非阻塞消息接收MPI IRecv buf count datatype source tag comm request ierr Inbuf count datatype source tag commOutrequest ierr非阻塞通信的完成MPI Wait request status ierr 等待消息收发完成MPI Test request flag stutus ierr MPI Waitall const request array status ierr 等待多个消息完成InrequestOutstatus flag logical型 9 CopyrightbyLiXinliang 非阻塞通信调用后立即返回 缓冲区不能立即使用Sum 计算某变量MPI Isend sum 发送该变量sum 不能给变量重新赋值 发送可能尚未完成 MPI Irecv sum1 sum sum0 sum1数据不能立即使用 接收可能未完成 MPI Isend sum request CallMPI Wait request status ierr Sum MPI Irecv sum1 request CallMPI Wait request status ierr Sum sum0 sum1 10 CopyrightbyLiXinliang 利用通信与计算重叠技术提高效率 例 计算差分串行程序realA N N B N N h Doi 1 NB I 1 A I 2 A I 1 hB I N A I N A I N 1 henddoDoj 2 N 1Doi 1 NB i j A i j 1 A i j 1 2 h EnddoEnddo 0 J 1 2 3 N 1 N i 1i 2i N 11 CopyrightbyLiXinliang 并行程序 以两个进程并行为例realA N N 2 B N N 2 A1 N h If myid eq 0 thencallMPI send A 1 N 2 N MPI real 1 99 MPI Comm world ierr callMPI recv A1 N MPI real 1 99 MPI Comm World status ierr ElsecallMPI recv A1 N MPI real 0 99 MPI Comm World status ierr callMPI send A 1 1 N MPI real 0 99 MPI Comm world ierr endif 0 1 J 1 2 N 2 A 1 N 2 A 2 N 2 A 3 N 2 A N N 2 12 CopyrightbyLiXinliang If myid eq 0 thenDoi 1 NB i 1 A i 2 A i 1 hB i N A1 i A i N 1 2 h EnddoElseDoi 1 NB i 1 A i 2 A1 i 2 h B i N A i N A i N 1 hEnddoendifDoj 2 N 1Doi 1 NB i j A i j 1 A i j 1 2 h EnddoEnddo 0 1 J 1 2 N 2 特点 先收发边界信息再进行计算缺点 通信过程中CPU空闲 13 CopyrightbyLiXinliang 内边界 通信与计算重叠realA N N 2 B N N 2 A1 N hintegermyid ierr req1 req2 status If myid eq 0 thencallMPI ISend A 1 N 2 N MPI real 1 99 MPI Comm world req1 ierr callMPI Irecv A1 N MPI real 1 99 MPI Comm World req2 ierr ElsecallMPI Irecv A1 N MPI real 0 99 MPI Comm World req2 ierr callMPI Isend A 1 1 N MPI real 0 99 MPI Comm world req1 ierr endif 0 1 J 1 2 N 2 14 CopyrightbyLiXinliang Doj 2 N 1Doi 1 NB i j A i j 1 A i j 1 2 h EnddoEnddoCallMPI wait req2 statue ierr If myid eq 0 thenDoi 1 NB I 1 A I 2 A I 1 hB I N A1 i A I N 1 2 h EnddoElseDoi 1 NB I 1 A I 2 A1 i 2 h B I N A1 i A I N 1 hEnddoendif 0 1 J 1 2 N 2 特点 传递边界信息同时进行计算 内点 读取系统时间doubleprecisiontimetime MPI Wtime 15 CopyrightbyLiXinliang 二 如何收发非连续数据例如 发送数组的一行A 100 50 发送A 1 1 A 1 2 A 1 3 A 1 1 A 1 2 A 1 3 方法1 多次发送通信开销大 效率低 16 CopyrightbyLiXinliang 方法2 将发送的数据拷贝到连续的数组中dimensionA 100 50 B 50 If myid eq 0 thenDoi 1 50B i A 1 i EnddocallMPI Send B 50 MPI REAL 1 99 MPI COMM WORLD ierr ElsecallMPI Recv B 50 MPI Real 0 99 Doi 1 50A 1 i B i Enddoendif 不足 额外的内存占用额外的拷贝操作通信不复杂的情况 内存拷贝工作量不大 该方法也可以采用 效果还可以 17 CopyrightbyLiXinliang 方法3 构建新的数据结构 Count 块的数量 blocklength 每块的元素个数Stride 跨度 各块起始元素之间的距离 Oldtype 旧数据类型 Newtype 新数据类型 整数 例 integerMY TYPECallMPI TYPE VECTOR 4 1 3 MPI REAL MY TYPE ierr CallMPI TYPE Commit MY TYPE ierr 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数组的一行RealA 3 4 A 1 在内存中的排列次序 18 CopyrightbyLiXinliang 例 发送三维数组中的一个面 Fortran 数组 realA 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 MPI Type Vector P M M N MPI Real My Type ierr 通信3 连续分布 无需构造新类型 19 CopyrightbyLiXinliang MPI TYPE INDEXED count array of blocklengths array of displacements oldtype newtype ierr 构造数据类型更灵活的函数 直接指定每块的元素个数及偏移量 块的数量 整数 每块元素的个数 整形数组 每块的偏移量 整形数组 例 数组realA 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 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 ierrcount Ndok 1 Nblocklengthes k kdisplacements k k 1 NenddocallMPI TYPE INDEXED count blocklengths displacements MPI REAL newtype ierr CallMPI TYPE Commit Newtype ierr callMPI Send A 1 1 1 Newtype 20 CopyrightbyLiXinliang N 三 MPI的通信域和组预定义通讯域MPI Comm World 包含所有进程的组通讯域的分割MPI Comm Split comm color key New Comm 0 2 1 4 3 5 7 6 8 9 10 11 Color相同的进程在同一组根据key的大小排序 key相同时按原ID排序 例如 12个进程 分成3行4列Integermyid 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 column CallMPI Comm rank Comm Raw myid raw ierr CallMPI Comm rank Comm line myid line ierr MPI Comm World RAW Column Color 分组标准 Key 排序依据如相同 按原ID排 提交新定义的组 否则新组无效 不要忘记 计算行号 列号 21 CopyrightbyLiXinliang 例 计算差分三维分割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得到组内编号建立三个方向通讯的数据结构4 通信 计算内点差分5 计算边界差分 0 2 1 4 3 5 7 6 8 9 10 11 MPI Comm World 22 CopyrightbyLiXinliang Parameter M1 M NM N1 N NN P1 P NP RealA 0 M1 1 0 N1 1 0 P1 1 Integermyid Comm X Comm Y Comm Z id X id Y id Z request 12 CallMPI Comm Rank MPI Comm World myid ierr CallMPI Comm Split MPI Comm World mod myid NM 0 Comm X ierr CallMPI Comm Split MPI Comm World mod myid NM NN NM 0 Comm Y ierr CallMPI Comm Split MPI Comm World myid NM NN 0 Comm Z ierr CallMPI Comm Rank Comm X id x ierr CallMPI Comm Rank Comm Y id y ierr CallMPI Comm Rank Comm Z id z ierr 定义三个方向的通信域 23 CopyrightbyLiXinliang CallMPI Type Vector N1 2 P1 2 1 M1 2 MPI real Type X ierr CallMPI Type Vector P1 2 N1 2 M1 2 N1 2 MPI real Type Y ierr CallMPI Type Commit Type X ierr CallMPI 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 NMCallMPI Isend A 1 0 0 1 TYPE X id X Pre 99 Comm X request 1 ierr CallMPI Isend A M1 0 0 1 TYPE X id X next 99 Comm X request 2 ierr CallMPI Irecv A 0 0 0 1 TYPE X id X next 99 Comm X request 3 ierr CallMPI Irecv A M1 1 0 0 1 TYPE X id X Pre 99 Comm X request 4 ierr 定义新的数据结构 24 CopyrightbyLiXinliang Dok 2 P1 1Doj 2 N1 1Doi 2 M1 1Ax 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 EnddoEnddoEnddocallMPI Wait All 12 request status ierr dok 1 P1doj 1 N1Ax 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 enddoEnddo 内点 边界点 25 CopyrightbyLiXinliang 四 分布数组的文件存储分布数组realA M m1 N n1 存储方式1 每个进程存储到独立的文件realA M m1 N n1 character len 50 filename write filename file I4 4 dat myidopen 55 file filename form unformatted write 55 Aclose 55 file 0000 datfile 0001 datfile 0002 dat 优点 程序简单缺点 数据文件多 不易处理 改变处理器数目时需特殊处理 0 1 2 3 26 CopyrightbyLiXinliang 分布数组realA M m1 N n1 存储方式2 收集到0节点存储存储到一个文件缺点 改变处理器规模时 需要处理存储方式3 收集到0节点 重新装配成大数组收集A M m1 N n1 组成A0 M N realA0 M N A M m1 N n1 A1 M m1 N n1 if myid eq 0 thendok 0 m1 n1callMPI recv A1 M m1 N n1 MPI real k A0 i global j global A1 i j 把A1装配到A0enddoWrite 33 A0elsecallMPI Send A endif 0 1 2 3 0 1 2 3 0 27 CopyrightbyLiXinliang 存储方式4 按列搜集后存储 RealAj M If myid eq 0 thenopen 33 file A dat form binary doj 1 N收集矩阵A0的第j列存储到Aj write 33 AjenddoElse endif 第1列第2列第3列 优点 存储的数据形式与内存中A0的存放格式一致 存储的文件串行程序可直接读取realA M N open 55 file A dat form binary read 55 Aclose 55 28 CopyrightbyLiXinliang 存储方式5并行IO MPI2 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 offset buff const datatype status ierr offset偏移 buff缓冲区 const数目 29 CopyrightbyLiXinliang Part3实例教学 CFD程序的MPI实现实例 1 用拟谱方法求解不可压N S方程实例 2 用流水线方法计算紧致差分常用的优化方法 30 CopyrightbyLiXinliang 回顾基本的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 消息接收MPI Recv buf count datatype source tag comm status ierr 31 CopyrightbyLiXinliang 非阻塞消息发送MPI ISend buf count datatype dest tag comm request ierr Inbuf count datatype dest tag commOutrequest ierrRequest 返回的非阻塞通信对象 整数 非阻塞消息接收MPI IRecv buf count datatype source tag comm request ierr Inbuf count datatype source tag commOutrequest ierr非阻塞通信的完成MPI Wait request status ierr 等待消息收发完成MPI Test request flag stutus ierr MPI Waitall const request array status ierr 等待多个消息完成InrequestOutstatus flag logical型 32 CopyrightbyLiXinliang 发送非连续数据 构建新的数据结构 MPI TYPE VECTOR count blocklength stride oldtype newtype ierr Count 块的数量 blocklength 每块的元素个数Stride 跨度 各块起始元素之间的距离 Oldtype 旧数据类型 Newtype 新数据类型 整数 例 integerMY TYPECallMPI TYPE VECTOR 50 1 100 MPI REAL MY TYPE ierr CallMPI TYPE Commit MY TYPE ierr 33 CopyrightbyLiXinliang 通讯域的分割MPI Comm Split comm color key New Comm 0 2 1 4 3 5 7 6 8 9 10 11 Color相同的进程在同一组根据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 CallMPI Comm rank Comm Raw myid raw ierr CallMPI Comm rank Comm line myid line ierr MPI Comm World 34 CopyrightbyLiXinliang 实例1 用 拟 谱方法求解二维不可压N S方程 2p 物理模型 周期性边界条件 按照给定能谱布置初始流动研究流动的演化规律 35 CopyrightbyLiXinliang Fourier变换 1D Fourier变换的特点 求导数 乘积 困难 非线性项 卷积 计算量巨大 在物理空间计算 Fourier变换的快速算法FFT 36 CopyrightbyLiXinliang 二维Fourier变换 两次一维Fourier变换 37 CopyrightbyLiXinliang 求解步骤 1 读入初值2 调用FFT得到3 计算调用FFT得到4 计算调用FFT得到5 计算6 积分求出下一时间步的值7 调用FFT得到8 循环3 7 直到给定的时间 38 CopyrightbyLiXinliang 实际计算中 要采用抑制混淆误差的措施 程序的并行化 二维FFT 二维FFT 调用两次一维FFT一维FFT算法复杂 并行化难度大二维FFT的并行 重新分布 SubroutineFFT2d nx ny u integernx nyComplexu nx ny Fu nx ny u1 ny u2 nx doi 1 nxu1 u i callFFT1d ny u1 Fu i u1 enddodoj 1 nyu2 Fu j callFFT1d nx u1 u j u1 enddoend 39 CopyrightbyLiXinliang 数据重分布的实现 A1 M P N A2 M N P 1 2 3 4 a b c d 对等式编程思想 我 需要完成的工作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 成A2 40 CopyrightbyLiXinliang 程序 SubroutineRedistibute ItoJ A1 A2 M N P IntegerM N P k ierr status MPI Status Size realA1 M P N A2 M N P B1 M P N P P B2 M P N P P dok 1 PB1 P A1 k 1 N P 1 k N P callMPI Send B1 M N P P MPI Real k 1 Enddodok 1 PcallMPI Recv B2 M N P P MPI Real k 1 A2 k 1 M P 1 k M P B2 P Enddoend问题 全部发送 发送成功后再启动接收 容易死锁 按行分布 按列分布 41 CopyrightbyLiXinliang SubroutineRedistibute ItoJ A1 A2 M N P IntegerM N P k ierr status MPI Status Size realA1 M P N A2 M N P B1 M P N P P B2 M P N P P dok 1 PB1 P A1 k 1 N P 1 k N P id send myid kmodPid recv myid kmodPcallMPI Send B1 M N P P MPI Real id send callMPI Recv B2 M N P P MPI Real id recv A2 k 1 M P 1 k M P B2 P Enddoend问题 按顺序发送 接收 不易死锁 42 CopyrightbyLiXinliang 数据全交换 MPI AlltoAll sendbuf sendcount sendtype recvbuf recvcount recvtype comm ierr sendbuf发送缓冲区 首地址 recvbuf接收缓冲区 首地址 sendcount发送数目recvcount接收数目sendtype发送类型recvtype接收类型Comm通信域ierr整数 返回错误值 0为成功 To0 To1 To2 To3 Sendbuf的数据格式 sendcount From0 From1 From2 From3 Recvbuf的数据格式 recvcount 43 CopyrightbyLiXinliang 程序 SubroutineRedistibute ItoJ A1 A2 M N P IntegerM N P k ierr status MPI Status Size realA1 M P N A2 M N P B1 M P N P P B2 M P N P P dok 1 PB1 P A1 k 1 N P 1 k N P enddocallMPI AlltoAll B1 M N P P MPI Real B2 M N P P MPI Real MPI Comm World ierr dok 1 PA2 k 1 M P 1 k M P B2 P Enddoend问题 无法做到计算与通信重叠 44 CopyrightbyLiXinliang 二维并行FFT的实现 输入数据 输出数据均为按列分布 1 调用一维FFT实现i 方向的变换u u12 重新分布数据 按列 按行 u1 u2调用一维FFT实现j 方向的变换u2 Fu2重新分布数据 按行 按列 Fu2 Fu 45 CopyrightbyLiXinliang 实例 2 利用流水线实现紧致差分的并行化 紧致型差分格式 相同网格点上引入更多信息 性能更优化 是的差分逼近 普通差分格式 显式给出Fi的表达式 紧致型差分格式 隐式给出Fi的表达式 6阶中心 6阶对称紧致 Lele 5阶迎风紧致 Fu j 2j 1jj 1j 2 46 CopyrightbyLiXinliang 普通差分格式 直接计算导数 并行容易 紧致格式的计算 递推 递推公式 计算出 由边界条件或边界格式给出 2 由递推计算出全部导数 后面的数据必须等待前一步计算完成 无法并行 47 CopyrightbyLiXinliang 二维问题 流水线法求解 流水线示意图 步骤 1 计算d 2 fork 1 M 如果myid 0 计算F k 0 否则从myid 1接收F k 0 fori 1 N1 N1 N P 计算F k i 如果myid P 1向myid 1发送F k N1 缺点 通信次数过多 48 CopyrightbyLiXinliang 通信次数过于频繁 解决方法 分块流水线 步骤 1 计算d 2 forkp 1 MP 如果myid 0 计算F kp 0 否则从myid 1接收F kp 0 forj 1 N1 N1 N P 计算F kp j 如果myid P 1向myid 1发送F kp N1 F kp i 表示第kp块 49 CopyrightbyLiXinliang 对称紧致格式 追赶法 令 则 代入 1 得 对比 2 得 边界处导数可由边界条件或边界格式给出 则 步骤 1 2 由 3 式递推 得到3 4 由 2 式递推 得到 特点 两次递推 并行方法与前文类似 50 CopyrightbyLiXinliang 常用的并行优化方法1 通信与计算重叠采用非阻塞通信Isend Irecv2 用重复计算代替通信3 拆分长消息 合并短消息4 优化通信方式 51 CopyrightbyLiXinliang 用重复计算代替通信例如 计算差分u分布存储 f u 为u的函数 方法1 计算出v f u 通信得到uN 1 vN 1计算差分 方法2 计算出v f u 通信得到uN 1 边界外 计算出vN 1 f uN 1 计算差分 方法2 计算量大 通信量小当函数f u 不复杂时 可提高效率 1 2NN 1 52 CopyrightbyLiXinliang 长消息切割成多个短消息发送 接收callMPI Send A 1 100000 MPI Real 1 改为 dom 1 10callMPI Send A m 1 10000 1 10000 MPI real 1 enddo长消息 非缓冲 短消息 缓冲 缓冲区 MPI Send 缓冲区 MPI Send MPI Recv MPI Recv 53 CopyrightbyLiXinliang 合并短消息 dom 1 100callMPI Send A 1 m 1 MPI real 1 enddo改为dom 1 100B m A 1 m enddocallMPI Send B 1 100 MPI Real 1 54 CopyrightbyLiXinliang 优化通信方式 例 数据散发0号进程 数据A 100 散发给0 99方式1 0进程执行100次MPI Send其他进程执行MPI RecvMPI Scatter 采用该算法方式2 0进程把A 100 切割成10份 发送给10个进程10个进程接收A1 10 后再散发 55 CopyrightbyLiXinliang OpenMP并行编程入门 一 特点1 针对共享内存计算机结构全部CPU 线程均可访问内存2 程序改动量小 实现方便 以编译指示符为主 3 适用于小规模并行或与MPI配合进行大规模并行 内存 CPU 核心 CPU 核心 CPU 核心 1台PC机 1个计算节点 共享内存构架 CPU 内存 CPU 内存 CPU 内存 外部网络 节点1 节点2 Cluster结构 分布内存构架 print code1 OMPPARALLELprint code2 OMPENDPARALLELprint code3 end 例1 test1 f90 编译 在深腾7000 运行结果 屏幕截图 iforttest1 f90 openmp 添加 openmp选项 运行 1 设置线程数 并行执行的数目 exportOMP NUM THREADS 4 例如 4个 2 执行 a out 显示结果 code1code2code2code2code2code3 并行域中的代码执行了4次 Test2 f90 print code1 OMPPARALLELprint code2 OMPPARALLELprint code3 OMPENDPARALLEL OMPENDPARALLELprint code4 end DO循环分解 openMP最常用的并行方法 OMPPARALLEL OMPDOdok 1 12print kenddo OMPENDDO OMPENDPARALLELend 示例 线程0k 1 2 3 线程1k 4 5 6 线程2k 7 8 9 线程2k 10 11 12 OMPPARALLEL OMPDO OMPPARALLELDO 简写 运行结果 屏幕截图 运行结果 123789456101112 线程0 线程2 线程1 线程3 implicitnoneinteger parameter N 100000000integer kreal 8 dimension allocatable x y zreal 8 time1 time2 OMP get wtimeallocate x N y N z N time1 OMP get wtime OMPPARALLELDOSHARED x y z PRIVATE k dok 1 Nx k k 1 d0 N 1 d0 y k k 1 d0 N 1 d0 z k x k y k enddo OMPENDPARALLELDO time2 OMP get wtime deallocate x y z print TotalWallTimeis time2 time1end 例 test4 屏幕截图 采用单线程执行 耗时2 15秒采用2线程执行 耗时1 43秒采用4线程执行 耗时1 28秒 三 OpenMP的数据结构 共享与私有 OMPPARALLELDOdok 1 6print kenddo OMPENDPARALLELDOend 线程0k 线程1k 循环变量k在两个线程中的值是不同的 K是一个进程私有变量 PRIVATE 共享变量 全体进程均可访问的公共变量私有变量 各个进程私有的变量 x 8 0 y x 2 0 OMPPARALLELDOSHARED x y PRIVATE k z dok 1 6z k x yprint x y zenddo OMPENDPARALLELDOend 线程0k z 线程1k z x y 私有变量 公共变量 例 将下面代码并行化 Integer parameter N 1024Real dimension N x y zRealr 给x y赋值 Dok 1 Nr sqrt x k x k y k y k z k 1 0 1 0 r Enddo 关键 分析哪些是共享变量 哪些是私有变量 显然 r k是私有变量 其他均为共享变量 OMPPARALLELDOSHARED DEFAULT PREATE r k Dok 1 Nr sqrt x k x k y k y k z k 1 0 1 0 r Enddo OMPENDPARALLELDO 四 OpenCFD EC3D的OpenMP并行化举例 OMPPARALLELDOdok 1 nz 1doj 1 ny 1doi 1 nx 1d 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 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 enddoenddoenddo OMPENDPARALLELDO OMPPARALLELDODEFAULT 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 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 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 5d0zj 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 5d0yk 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 5d0zk 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 5d0Jac 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 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 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 64 CopyrightbyLiXinliang 计算流体力学 课程2011习题汇总 1推导无量纲的Navier Stokes方程1 2对于一维Euler方程组推导Jocabian矩阵以及中的表达式 要求 给出具体推导过程 切忌从书上抄录公式 越详细越好 CopyrightbyLiXinliang 65 2 1公式推导 1 一激波从左向右传播 激波左侧物理量为 激波右侧压力为 试计算激波右侧的速度 2 有一扇膨胀波从左向右传播 膨胀波左侧物理量为 膨胀波右侧压力为 试计算膨胀波右侧的速度 膨胀波 要求 务必写出详细的步骤推导 越详细越好 切忌照抄书上的公式 2 2如下Sod激波管问题 求出理论解 并分别画出t 0 14时刻的分布曲线 CopyrightbyLiXinliang 66 3 1对如下单波方程 构建的差分格式如下 试利用Fourier方法 分析其稳定性 4 1构造高分辨率差分格式 并进行理论分析及数值实验 针对单波方程 对于空间导数 构造出一种不超过6点格式 并进行Fourier误差分析 画出kr ki的曲线 要求 精度不限 网格基架点数不超过6个 能够分辨的波数范围尽量宽 即kr ki曲线近可能接近准确解 给出差分的具体表达式 画出kr ki的曲线 说明构造格式的阶数 并采用本PPT第5页的方法给出的精度验证 形如 另外 进行如下数值验证 空间采用20个网格点 采用新构造的差分格式离散 时间推进采用3步Runge Kutta方法 时间步长可足够小 例如0 01 给出t 20 50两个时刻的数值解 与精确解比较 画图 并给出数值解的L2模误差 68 CopyrightbyLiXinliang 提示 1 如不使用优化技术 则格式构造方法简单 Taylor展开后解代数方程组即可 2 建议尝试使用优化技术 例 假设格式形式如下 如果要求其有5阶精度 则通过Taylor展开可得到6个方程 6个系数可直接解出 我们要求其有4阶精度 当然3阶 2阶也可 于是Taylor展开只能提供5个方程 6个未知数 a1 a6 5个方程 有1个自由参数 调整这个自由参数 使得kr ki曲线最为理想 如何调整 1 可以人工调整 观察kr ki曲线 选取满意的 2 可自动调整 设立一个优化目标函数 例如调整自由参数 使得该目标函数取最大值 思路 牺牲精度 提高分辨率 69 CopyrightbyLiXinliang 4 2数值求解Sod激波管问题 计算其数值解