第十七讲-并行计算

八、并行计算CIG-ComputationalInfrastructureforGeodynamics/GeoFEMProjecthttp://geofem.tokyo.rist.or.jp/QuakeSimProject/TeraShake–TypicalFramework/projects/scec/terashake/GEONGRIDProject:

常见的大型地球科学并行计算软件系统GCNfileSourceCodeGeneratorPDEfileGIOfilePreprocessinginterfaceDatastructureMeshgenerationAllthesourcecodeforfiniteelementcomputinginterfacetosolversPostprocessingVisualizationtoolsinterfaceMeshgenerationtoolsComputingdomainconstructiontoolsMeshgenerationlibMeshpartitionlibDatapartitionlibSequentialandParallelsolverslibFinishedmodelslibAuxiliaryroutineslibAmira,GID,OpenDX,ParallelView,RocksclustersVizArchitectureprofileofparallelfiniteelementprogramgenerator/Computingenvironment张怀，2010pFEPG并行有限元程序自动生成系统的组织架构dispuvcoorxyfuncfunafunbfuncshap%1%2gaus%3mass%1load=fufv$c6pe=prmt(1)$c6pv=prmt(2)$c6fu=prmt(3)$c6fv=prmt(4)$c6fact=pe/(1.+pv)/(1.-2.*pv)funcfuna=+[u/x]funb=+[v/y]func=+[u/y]+[v/x]stifdist=+[funa;funa]*fact*(1.-pv)+[funa;funb]*fact*(pv)+[funb;funa]*fact*(pv)+[funb;funb]*fact*(1.-pv)+[func;func]*fact*(0.5-pv)*es,em,ef,Estifn,Estifv,*es(k,k),em(k),ef(k),Estifn(k,k),Estifv(kk),

goto(1,2),ityp1callseuq4g2(r,coef,prmt,es,em,ec,ef,ne)goto32callseugl2g2(r,coef,prmt,es,em,ec,ef,ne)goto33continueDOJ=1,NMATEPRMT(J)=EMATE((IMATE-1)*NMATE+J)EnddoPRMT(NMATE+1)=TIMEPRMT(NMATE+2)=DTprmt(nmate+3)=imateprmt(nmate+4)=numOtherelementmatrixcomputingSubsPDEexpressionContainsinformationofthephysicalmodel,suchasvariablesandequationsforgeneratingelementstiffnessmatrix.FortranSegmentscodesthatrealizethephysicalmodelatelementlevel.variablesequationAutomatedCodeGeneratorStep1:FromPDEexpressiontoFortransegmentsSegment1Segment2Segment3Segment4Step2:FromalgorithmexpressiontoFortransegments

doi=1,kdoj=1,kestifn(i,j)=0.0enddoenddodoi=1,kestifn(i,i)=estifn(i,i)doj=1,kestifn(i,j)=estifn(i,j)+es(i,j)enddoenddo

U(IDGF,NODI)=U(IDGF,NODI)*+ef(i)defistifSmassMloadFtypeemdtylstep0equationmatrix=[S]FORC=[F]SOLUTIONUwrite(s,unod)UendAlgorithmExpressionContainsinformationforformingglobalstiffnessmatrixforthemodel.FortranSegmentscodesthatrealizethephysicalmodelatgloballevel.StiffnessmatrixSegment5Segment6SUBROUTINEETSUB(KNODE,KDGOF,IT,KCOOR,KELEM,K,KK,*NUMEL,ITYP,NCOOR,NUM,TIME,DT,NODVAR,COOR,NODE,#SUBET.sub*U)implicitdoubleprecision(a-h,o-z)DIMENSIONNODVAR(KDGOF,KNODE),COOR(KCOOR,KNODE),*U(KDGOF,KNODE),EMATE(300),#SUBDIM.sub*R(500),PRMT(500),COEF(500),LM(500)#SUBFORT.sub#ELEM.subCWRITE(*,*)'ESEMEF='CWRITE(*,18)(EF(I),I=1,K)#MATRIX.subL=0M=0I=0DO700INOD=1,NNE………U(IDGF,NODI)=U(IDGF,NODI)#LVL.subDO500JNOD=1,NNE………500CONTINUE700CONTINUE………returnendProgramStencilFortranSegmentsgeneratedStep3:PlugFortransegmentsintoastencil,formingfinalFEprogramSegment1Segment2Segment4Segment3Segment5Segment6…………..Fromphysicalmodeltonumericalcomputation:

AdifferentapproachPartialDifferenceEquations(PDEs)PDE2FEMserverDataGrid(GEONandothers)

High-performancecomputingVianetworktouploadVianetworktouploadNumericalresultsVianetworktodownloadpost-processingAnalyzeandvisualitionAutomatedsourcecodegeneratorAutomaticsourcecodegeneratorfuncfuna=+[u/x]

………funf=+[u/y]+[v/x]

………dist=+[funa;funa]*d(1,1)+[funa;funb]*d(1,2)+[funa;func]*d(1,3)+[funb;funa]*d(2,1)+[funb;funb]*d(2,2)+[funb;func]*d(2,3)+[func;funa]*d(3,1)+[func;funb]*d(3,2)+[func;func]*d(3,3)+[fund;fund]*d(4,4)+[fune;fune]*d(5,5)+[funf;funf]*d(6,6)

load=+[u]*fu+[v]*fv+[w]*fw-[funa]*f(1)-[funb]*f(2)-[func]*f(3)-[fund]*f(4)-[fune]*f(5)-[funf]*f(6)PDEsCompletesourcecodeFEMModelingLanguageDataGrid(GEONandothers)PhysicalmodelModelresultsHPCCData=>???SWFSWF计算地球动力学的学科交叉特色-机遇与挑战并存并行有限元方法基础理论研究及其自动生成实现1997-2003：Lagrange乘子法，区域分解算法(DomainDecompositionMethod)，非连续变形分析算法(LagrangeMultiplierDiscontinuousDeformationAnalysis)的算法与并行自动生成实现;2002-2009：预条件Krylov子空间算法求解器的优化改进(Preconditionered

KrylovSubspaceIterativeSolvers)与并行有限元自动生成技术的理论与算法实现；多场强耦合系统并行有限元求解算法；国内空白，国际上目前只有类似的串行系统。国际上有类似有限元方法的自动求解软件包(Femlab,DiffPack,PDE2D等)，但是没有自动生成并行有限元源程序技术；已经形成商业产品pFEPG(ParallelFiniteElementProgramGenerator)；计算地球动力学若干模型的并行有限元计算算法研究与实现(2003-)板块与活动构造（ActiveTectonics）-粘弹性并行有限元模型；区域性和全球地震波数值模拟-大规模并行显式有限元算法与并行实现;地幔对流与大陆演化过程-多场强耦合系统求解算法;基础模型–线弹性模型与对流扩散温度场方程;计算地球动力学基础架构中的若干辅助系统–跟踪与科研需要(2005-)海量数据并行可视化设备(ParallelVisualization);

非结构化网格与并行CVT/SCVT网格并行生成算法；Preliminaryresults通过高效的分区算法实现，用户能自由指定分区，并自动完成各个节点的负载均衡；1,200,000quadrilateralstructuredmeshpartition,sub-domainswithdifferentcolors(left)andsub-domainboundaries(right)50,000triangleunstructuredmeshpartition,(partitionedby10)研究探讨一些地球科学问题中的基础科学问题：1区域性活动构造-Maxwell粘弹性模型;2断层系统的应力演化与地震活动性关系-粘-弹-塑性模型;3区域性板块构造演化历史–强非线性的非牛顿流体热演化问题、大变形（几何非线性）；4地球自转快速变化后的全球应力场变化分析；5区域性应力演化历史与地震活动性分析；关系到国际民生问题的研究-大的工程应用方面的基础研究工作6二维、三维地震波强地面运动数值模拟；7海啸数值模拟，危险性分析；8港口安全性评估；基础设置建设其它94D海量数据体的并行可视化和体绘制的并行显示；

10GPU计算与GPU+CPU层次结构计算的算法研究；

10建立计算地球动力学网格计算基础架构（国际合作，跟踪研究）；地球动力学方面的科学研究内容

能不能做一些关系国计民生的问题？地震海啸港口安全滑坡泥石流飓风……地震灾害预测与救助指导E-Science系统现有工作基础：全球地震波传播数值模拟；区域地震波数值模拟-采用有限元方法，具有差分法和谱元法不具备的优势，同时采用合理的并行显格式算法，同样可以达到差分法的效率和谱元法的高分辨率；突出的优点是可以考虑各种类型的单元；PeakGroundAcceleration(PGA)分析；功率谱密度（PowerSpectralDensity）分析，确定区域本征值分布特征；非结构化网格生成技术-可以充分考虑多尺度问题；数据库累积；并行可视化系统；网格国际合作，整合即可完成本研究项目；PFEPG并行有限元程序自动生成系统并行程序与串行程序的区别在于：1.是否分区，2.是否存在主从节点，3.是否存在节点之间的数据交换。PFEPG工作原理图客户端——PFEPG服务器——并行机三者之间的关系应用pFEPG进行并行计算：

客户机不是并行机的节点pFEPG并行程序结构解析算法基础：区域分解求解问题范围：动态问题稳态问题线性问题非线性问题元件程序组成（六个模块）前处理分区模块（partition）START、BFT、Solve、E和U元件程序注：其中E和U元件程序根据用户给出的偏微分方程表达式和算法表达式由系统自动生成，其余几个元件程序由系统给出，不随表达式的变动而变动。本系统的基本思想是元件化程序设计和程序自动生成技术相结合，根据区域分解法的特点，将有限元程序分成：前处理分区、START、BFT、Solve、E和U等六个模块，其中E和U元件程序程序根据用户给出的偏微分方程表达式和算法表达式由系统自动生成，其余几个元件程序由系统给出，不随表达式的变动而变动。PFEPG并行有限元程序自动生成系统前处理分区模块：该模块又分成四个元件程序：partition、getpart、spart、sgetpart元件；该模块的功能事根据用户提供的前处理数据将整个求解区域根据用户的要求合理划分成数个子区域。前处理分区模块主进程：Partition.fGetpart.fPartition.dat从进程：Spart.fSgetpart.f该功能模块是根据用户提供的前处理数据将整个求解区域根据用户的要求合理划分成数个子区域，partition给出各个子区域所包含的整体节点号，局部节点号和整体节点编号的对应关系，并将各子区域的一些前处理信息，包括节点数、时间信息等，以及各区域所含节点的坐标值分块发送给与各子区域对应的从进程。Getpart根据partition提供的各子区域所含节点信息得到各自的局部单元信息、初边值信息、节点ID信息以及各子区域节点的各自由度对应的整体方程号等，并发送各相应从进程。Spart元件接收来自主进程partition元件的坐标值；Sgetpart接收来自主进程getpart元件发送的消息。前处理分区模块有限元网格分区示例分区图START元件：功能：在各子区域内给出每个节点的各个自由度与将要形成的代数方程组的变量（方程号）的对应关系。得到整体刚度矩阵的存储格式，该元件程序只在从进程执行。BFT元件：功能：计算边值，更新时间。该元件由两个程序：bftm和bfts组成，分别由主进程和从进程执行，bftm由主进程执行，主要是给出每一时刻解的边值，对时间进行更新，进行结束与否的判断，并将判断结果和新的时间和各区域解的边值发送给从进程

；bftm.fbfts由从进程执行，接收由主进程发送的新的时间和解的边值。bfts.fE元件：功能：该元件程序用于单元计算，包括单刚，单质和单元载荷计算，并把它们由节点自由度表示转换成代数方程组的变量表示，同时处理边界约束条件，并形成代数方程组的系数矩阵和右端项。SOLV元件：功能：组装和求解整体线性代数方程组。由两个程序组成：solv和msolvsolv由从进程执行Msolv由主进程执行U元件：功能：方程组解转换回节点自由度值，并进行其它后处理计算包括u元件程序和mu元件程序U元件由从进程执行，进行解的转换和后处理，并将结果传输到主进程Mu元件程序由主进程执行，接收来自于各子区域的解，整理得到整体区域上的解，更新子区域公共边界解，结果存盘处理。为了使并行程序运行结果与串行程序运行结果保持一致，同时避免由于并行I/O所带来的程序设计上的复杂性，我们规定从进程不参与磁盘文件的I/O，所有磁盘文件的I/O操作均由主进程完成。需要保存的数据写入预先开辟的内存段中，实现数据文件与特定内存段之间的对应关系。这些数据文件保存都将是一些数组，单个的数据被作为公共变量在各元件程序之间进行传递。数据结构众所周知，并行程序的调试非常困难，PFEPG系统的这种工作模式最大程度上降低了程序的调试难度。用户上传的有限元表达式和由服务器生成的计算源程序均是文本文件的形式，其大小一百K以内，在Internet上传输速度很快，一般来说只需十几秒即可生成用户所需的全部源程序，并传输到用户的本地机。下面说明PFEPG系统的使用步骤：1．填写Partition文件Partition.dat文件中是一些前处理分区程序需要的常数，内容一般如下41221其中，第一行的第一个整数4表示总的分块数，第二个整数1表示采用的分区算法（1表示采用图论的分区算法，2表示利用坐标分区的算法）；第二行的第一个2表示x方向剖分数，第二个2表示y方向剖分数，1表示z方向剖分数。第一行的第一个数是第二行三个数的乘积，在做大规模计算时，需要修改partition.dat文件来增加分块数，但这个关系需要保持。如果用户不填写此文件，系统会产生如上的缺省文件。2．用户填写完公式文件（PDE,GCN,GIO等）并通过界面传送到服务器后，运行mpi命令如下：mpigcn文件的前缀名即可产生并行有限元源程序，并传回客户机。数据文件说明partition.dat功能：指定和记录计算时分区数目和方法格式：41221M：总的分区数(>=2)N：分区方法，只可取1,2，分别表示按图论的分块方法和按坐标的分块方法MxMyMz:三个坐标轴方向的分区数M=Mx*My*Mz…MN…MxMyMz3．把以下文件传送到并行机全部并行fortran源程序（*.f文件），fegen.b、runmpi、time0、nbefile、partition.dat、*.io。如果使用的客户端是Window系统下的，对于上传的文件，一般需要在并行机上运行如下命令dos2unix*.*将上传的文件转成Unix文件格式。对于由多台PC机或工作站组成的并行机，PFEPG客户端安装在其中一台PC或工作站上的情况，不需要传送任何文件。4．在并行机上运行如下命令csh–xfegen.b通过编译器产生并行的可执行程序。5．运行前处理，mpigid，cshpreproc.mpi。6．在并行机上运行如下命令csh–xrunmpi进行并行计算。7．运行后处理，cshpostproc.mpi。pFEPG程序代码浅析pFEPG封装的几个基本通信函数sendint(des,src,idata)：从源进程src发送一个整型数idata到目的进程dessendr(des,src,rdata)：从源进程src发送一个实型数rdata到目的进程dessendai(des,src,idata,len)：从源进程src发送一个长度为len的整型数组idata到目的进程dessendar(des,src,rdata,len)：从源进程src发送一个长度为len的实型数组rdata到目的进程despFEPG程序代码浅析pFEPG封装的几个基本通信函数recvint(des,src,idata)：从源进程src接收一个整型数idata到目的进程desrecvr(des,src,rdata)：从源进程src接收一个实型数rdata到目的进程desrecvai(des,src,idata,len)：从源进程src接收一个长度为len的整型数组idata到目的进程desrecvar(des,src,rdata,len)：从源进程src接收一个长度为len的整型数组rdata到目的进程despFEPG程序代码浅析代码中的几个固定变量和数组说明几个常用的固定变量：numblk：分块数/从进程个数nsource：=0，主进程进程号iblk：块号/从进程进程号主要用于消息通信pFEPG程序代码浅析代码中的几个固定变量和数组说明几个常用的固定数组：knode_iblk(numblk)：子区域节点总数n_update(numblk)：子区域内部节点数mglnod(knode)：整体节点编号对应的子区域节点编号mlgnod(numblk,max_lknode)：子区域节点编号对应的整体节点编号iorder(numblk,max_lknode)：子区域节点标识，1为内部节点，0为扩展节点主要用于解的转换和整理算例：三维热固耦合问题

热固耦合问题模型控制方程热传导方程弹性体方程其中为弹性体位移，为作用于弹性体上的体积力。分别表示应力和应变向量。A是一个矩阵微分算子D为弹性矩阵Q为温度场产生的热应力向量，由下式确定：原型文件et.gcnet.gioheat.vdesdisp.vdessdisp.vde客户端需要准备的文件heat.vdedispucoorxyzshap%1%2gaus%3mateekq1.0d01;0.0d0;vectxxyzstifdist=+[u/x_i;u/x_i]*ekload=+[u]*qendsdisp.vdeFUNC$c6fact=pe/(1.+pv)/(1.-2.*pv)$c6ft=(1.+pv)*fact*alfaex=+[u/x]

ey=+[v/y]

ez=+[w/z]eyz=+[w/y]+[v/z]exz=+[w/x]+[u/z]exy=+[u/y]+[v/x]

stifdist=+[em_i;em_j]*de_i_j*fact+[en_i;en_i]*fact*(0.5-pv)

load=+[u_i]*fd_i+[em_i]*fe_i

enddispu,v,wcoorx,y,zfuncexeyezeyzexzexyshap%1%2gaus%3coeftnmatepepvalfafxfyfz1.0e9;0.3;1.0e-5;0;0;7e3;vectemexeyezvecteneyzexzexyvectuuvwvectfdfxfyfzvectfeft*tnft*tnft*tnmatrde33(1.-pv)pvpvpv(1.-pv)pvpvpv(1.-pv)ssdisp.vdedispsx,sy,sz,syz,sxz,sxycoorx,y,zcoefu,v,w,tnshap%1%2gaus%3mass%11.0matepepvalfafxfyfz1.0e9;0.3;1.0e-5;0;0;7e3;vectfssxsyszsyzsxzsxyvectfffsxfsyfszfsyzfsxzfsxyvectemexeyezvecteneyzexzexyvectfmfsxfsyfszvectfnfsyzfsxzfsxyvectfefttfttfttmatrde33(1.-pv)pvpvpv(1.-pv)pvpvpv(1.-pv)stif$c6fact=pe/(1.+pv)/(1.-2.*pv)$cvftt=(1.+pv)*alfa*tn*fact$cvex=+{u/x}$cvey=+{v/y}$cvez=+{w/z}$cveyz=+{w/y}+{v/z}$cvexz=+{w/x}+{u/z}$cvexy=+{u/y}+{v/x}$cvfm_i=+fact*de_i_j