第十七講-并行計算

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

常見的大型地球科學并行計算軟件系統(tǒng)GCNfileSourceCodeGeneratorPDEfileGIOfilePreprocessinginterfaceDatastructureMeshgenerationAllthesourcecodeforfiniteelementcomputinginterfacetosolversPostprocessingVisualizationtoolsinterfaceMeshgenerationtoolsComputingdomainconstructiontoolsMeshgenerationlibMeshpartitionlibDatapartitionlibSequentialandParallelsolverslibFinishedmodelslibAuxiliaryroutineslibAmira,GID,OpenDX,ParallelView,RocksclustersVizArchitectureprofileofparallelfiniteelementprogramgenerator/Computingenvironment張懷，2010pFEPG并行有限元程序自動生成系統(tǒng)的組織架構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計算地球動力學的學科交叉特色-機遇與挑戰(zhàn)并存并行有限元方法基礎理論研究及其自動生成實現(xiàn)1997-2003：Lagrange乘子法，區(qū)域分解算法(DomainDecompositionMethod)，非連續(xù)變形分析算法(LagrangeMultiplierDiscontinuousDeformationAnalysis)的算法與并行自動生成實現(xiàn);2002-2009：預條件Krylov子空間算法求解器的優(yōu)化改進(Preconditionered

KrylovSubspaceIterativeSolvers)與并行有限元自動生成技術的理論與算法實現(xiàn)；多場強耦合系統(tǒng)并行有限元求解算法；國內空白，國際上目前只有類似的串行系統(tǒng)。國際上有類似有限元方法的自動求解軟件包(Femlab,DiffPack,PDE2D等)，但是沒有自動生成并行有限元源程序技術；已經(jīng)形成商業(yè)產(chǎn)品pFEPG(ParallelFiniteElementProgramGenerator)；計算地球動力學若干模型的并行有限元計算算法研究與實現(xiàn)(2003-)板塊與活動構造（ActiveTectonics）-粘彈性并行有限元模型；區(qū)域性和全球地震波數(shù)值模擬-大規(guī)模并行顯式有限元算法與并行實現(xiàn);地幔對流與大陸演化過程-多場強耦合系統(tǒng)求解算法;基礎模型–線彈性模型與對流擴散溫度場方程;計算地球動力學基礎架構中的若干輔助系統(tǒng)–跟蹤與科研需要(2005-)海量數(shù)據(jù)并行可視化設備(ParallelVisualization);

非結構化網(wǎng)格與并行CVT/SCVT網(wǎng)格并行生成算法；Preliminaryresults通過高效的分區(qū)算法實現(xiàn)，用戶能自由指定分區(qū)，并自動完成各個節(jié)點的負載均衡；1,200,000quadrilateralstructuredmeshpartition,sub-domainswithdifferentcolors(left)andsub-domainboundaries(right)50,000triangleunstructuredmeshpartition,(partitionedby10)研究探討一些地球科學問題中的基礎科學問題：1區(qū)域性活動構造-Maxwell粘彈性模型;2斷層系統(tǒng)的應力演化與地震活動性關系-粘-彈-塑性模型;3區(qū)域性板塊構造演化歷史–強非線性的非牛頓流體熱演化問題、大變形（幾何非線性）；4地球自轉快速變化后的全球應力場變化分析；5區(qū)域性應力演化歷史與地震活動性分析；關系到國際民生問題的研究-大的工程應用方面的基礎研究工作6二維、三維地震波強地面運動數(shù)值模擬；7海嘯數(shù)值模擬，危險性分析；8港口安全性評估；基礎設置建設其它94D海量數(shù)據(jù)體的并行可視化和體繪制的并行顯示；

10GPU計算與GPU+CPU層次結構計算的算法研究；

10建立計算地球動力學網(wǎng)格計算基礎架構（國際合作，跟蹤研究）；地球動力學方面的科學研究內容

能不能做一些關系國計民生的問題？地震海嘯港口安全滑坡泥石流颶風……地震災害預測與救助指導E-Science系統(tǒng)現(xiàn)有工作基礎：全球地震波傳播數(shù)值模擬；區(qū)域地震波數(shù)值模擬-采用有限元方法，具有差分法和譜元法不具備的優(yōu)勢，同時采用合理的并行顯格式算法，同樣可以達到差分法的效率和譜元法的高分辨率；突出的優(yōu)點是可以考慮各種類型的單元；PeakGroundAcceleration(PGA)分析；功率譜密度（PowerSpectralDensity）分析，確定區(qū)域本征值分布特征；非結構化網(wǎng)格生成技術-可以充分考慮多尺度問題；數(shù)據(jù)庫累積；并行可視化系統(tǒng)；網(wǎng)格國際合作，整合即可完成本研究項目；PFEPG并行有限元程序自動生成系統(tǒng)并行程序與串行程序的區(qū)別在于：1.是否分區(qū)，2.是否存在主從節(jié)點，3.是否存在節(jié)點之間的數(shù)據(jù)交換。PFEPG工作原理圖客戶端——PFEPG服務器——并行機三者之間的關系應用pFEPG進行并行計算：

客戶機不是并行機的節(jié)點pFEPG并行程序結構解析算法基礎：區(qū)域分解求解問題范圍：動態(tài)問題穩(wěn)態(tài)問題線性問題非線性問題元件程序組成（六個模塊）前處理分區(qū)模塊（partition）START、BFT、Solve、E和U元件程序注：其中E和U元件程序根據(jù)用戶給出的偏微分方程表達式和算法表達式由系統(tǒng)自動生成，其余幾個元件程序由系統(tǒng)給出，不隨表達式的變動而變動。本系統(tǒng)的基本思想是元件化程序設計和程序自動生成技術相結合，根據(jù)區(qū)域分解法的特點，將有限元程序分成：前處理分區(qū)、START、BFT、Solve、E和U等六個模塊，其中E和U元件程序程序根據(jù)用戶給出的偏微分方程表達式和算法表達式由系統(tǒng)自動生成，其余幾個元件程序由系統(tǒng)給出，不隨表達式的變動而變動。PFEPG并行有限元程序自動生成系統(tǒng)前處理分區(qū)模塊：該模塊又分成四個元件程序：partition、getpart、spart、sgetpart元件；該模塊的功能事根據(jù)用戶提供的前處理數(shù)據(jù)將整個求解區(qū)域根據(jù)用戶的要求合理劃分成數(shù)個子區(qū)域。前處理分區(qū)模塊主進程：Partition.fGetpart.fPartition.dat從進程：Spart.fSgetpart.f該功能模塊是根據(jù)用戶提供的前處理數(shù)據(jù)將整個求解區(qū)域根據(jù)用戶的要求合理劃分成數(shù)個子區(qū)域，partition給出各個子區(qū)域所包含的整體節(jié)點號，局部節(jié)點號和整體節(jié)點編號的對應關系，并將各子區(qū)域的一些前處理信息，包括節(jié)點數(shù)、時間信息等，以及各區(qū)域所含節(jié)點的坐標值分塊發(fā)送給與各子區(qū)域對應的從進程。Getpart根據(jù)partition提供的各子區(qū)域所含節(jié)點信息得到各自的局部單元信息、初邊值信息、節(jié)點ID信息以及各子區(qū)域節(jié)點的各自由度對應的整體方程號等，并發(fā)送各相應從進程。Spart元件接收來自主進程partition元件的坐標值；Sgetpart接收來自主進程getpart元件發(fā)送的消息。前處理分區(qū)模塊有限元網(wǎng)格分區(qū)示例分區(qū)圖START元件：功能：在各子區(qū)域內給出每個節(jié)點的各個自由度與將要形成的代數(shù)方程組的變量（方程號）的對應關系。得到整體剛度矩陣的存儲格式，該元件程序只在從進程執(zhí)行。BFT元件：功能：計算邊值，更新時間。該元件由兩個程序：bftm和bfts組成，分別由主進程和從進程執(zhí)行，bftm由主進程執(zhí)行，主要是給出每一時刻解的邊值，對時間進行更新，進行結束與否的判斷，并將判斷結果和新的時間和各區(qū)域解的邊值發(fā)送給從進程

；bftm.fbfts由從進程執(zhí)行，接收由主進程發(fā)送的新的時間和解的邊值。bfts.fE元件：功能：該元件程序用于單元計算，包括單剛，單質和單元載荷計算，并把它們由節(jié)點自由度表示轉換成代數(shù)方程組的變量表示，同時處理邊界約束條件，并形成代數(shù)方程組的系數(shù)矩陣和右端項。SOLV元件：功能：組裝和求解整體線性代數(shù)方程組。由兩個程序組成：solv和msolvsolv由從進程執(zhí)行Msolv由主進程執(zhí)行U元件：功能：方程組解轉換回節(jié)點自由度值，并進行其它后處理計算包括u元件程序和mu元件程序U元件由從進程執(zhí)行，進行解的轉換和后處理，并將結果傳輸?shù)街鬟M程Mu元件程序由主進程執(zhí)行，接收來自于各子區(qū)域的解，整理得到整體區(qū)域上的解，更新子區(qū)域公共邊界解，結果存盤處理。為了使并行程序運行結果與串行程序運行結果保持一致，同時避免由于并行I/O所帶來的程序設計上的復雜性，我們規(guī)定從進程不參與磁盤文件的I/O，所有磁盤文件的I/O操作均由主進程完成。需要保存的數(shù)據(jù)寫入預先開辟的內存段中，實現(xiàn)數(shù)據(jù)文件與特定內存段之間的對應關系。這些數(shù)據(jù)文件保存都將是一些數(shù)組，單個的數(shù)據(jù)被作為公共變量在各元件程序之間進行傳遞。數(shù)據(jù)結構眾所周知，并行程序的調試非常困難，PFEPG系統(tǒng)的這種工作模式最大程度上降低了程序的調試難度。用戶上傳的有限元表達式和由服務器生成的計算源程序均是文本文件的形式，其大小一百K以內，在Internet上傳輸速度很快，一般來說只需十幾秒即可生成用戶所需的全部源程序，并傳輸?shù)接脩舻谋镜貦C。下面說明PFEPG系統(tǒng)的使用步驟：1．填寫Partition文件Partition.dat文件中是一些前處理分區(qū)程序需要的常數(shù)，內容一般如下41221其中，第一行的第一個整數(shù)4表示總的分塊數(shù)，第二個整數(shù)1表示采用的分區(qū)算法（1表示采用圖論的分區(qū)算法，2表示利用坐標分區(qū)的算法）；第二行的第一個2表示x方向剖分數(shù)，第二個2表示y方向剖分數(shù)，1表示z方向剖分數(shù)。第一行的第一個數(shù)是第二行三個數(shù)的乘積，在做大規(guī)模計算時，需要修改partition.dat文件來增加分塊數(shù)，但這個關系需要保持。如果用戶不填寫此文件，系統(tǒng)會產(chǎn)生如上的缺省文件。2．用戶填寫完公式文件（PDE,GCN,GIO等）并通過界面?zhèn)魉偷椒掌骱?，運行mpi命令如下：mpigcn文件的前綴名即可產(chǎn)生并行有限元源程序，并傳回客戶機。數(shù)據(jù)文件說明partition.dat功能：指定和記錄計算時分區(qū)數(shù)目和方法格式：41221M：總的分區(qū)數(shù)(>=2)N：分區(qū)方法，只可取1,2，分別表示按圖論的分塊方法和按坐標的分塊方法MxMyMz:三個坐標軸方向的分區(qū)數(shù)M=Mx*My*Mz…MN…MxMyMz3．把以下文件傳送到并行機全部并行fortran源程序（*.f文件），fegen.b、runmpi、time0、nbefile、partition.dat、*.io。如果使用的客戶端是Window系統(tǒng)下的，對于上傳的文件，一般需要在并行機上運行如下命令dos2unix*.*將上傳的文件轉成Unix文件格式。對于由多臺PC機或工作站組成的并行機，PFEPG客戶端安裝在其中一臺PC或工作站上的情況，不需要傳送任何文件。4．在并行機上運行如下命令csh–xfegen.b通過編譯器產(chǎn)生并行的可執(zhí)行程序。5．運行前處理，mpigid，cshpreproc.mpi。6．在并行機上運行如下命令csh–xrunmpi進行并行計算。7．運行后處理，cshpostproc.mpi。pFEPG程序代碼淺析pFEPG封裝的幾個基本通信函數(shù)sendint(des,src,idata)：從源進程src發(fā)送一個整型數(shù)idata到目的進程dessendr(des,src,rdata)：從源進程src發(fā)送一個實型數(shù)rdata到目的進程dessendai(des,src,idata,len)：從源進程src發(fā)送一個長度為len的整型數(shù)組idata到目的進程dessendar(des,src,rdata,len)：從源進程src發(fā)送一個長度為len的實型數(shù)組rdata到目的進程despFEPG程序代碼淺析pFEPG封裝的幾個基本通信函數(shù)recvint(des,src,idata)：從源進程src接收一個整型數(shù)idata到目的進程desrecvr(des,src,rdata)：從源進程src接收一個實型數(shù)rdata到目的進程desrecvai(des,src,idata,len)：從源進程src接收一個長度為len的整型數(shù)組idata到目的進程desrecvar(des,src,rdata,len)：從源進程src接收一個長度為len的整型數(shù)組rdata到目的進程despFEPG程序代碼淺析代碼中的幾個固定變量和數(shù)組說明幾個常用的固定變量：numblk：分塊數(shù)/從進程個數(shù)nsource：=0，主進程進程號iblk：塊號/從進程進程號主要用于消息通信pFEPG程序代碼淺析代碼中的幾個固定變量和數(shù)組說明幾個常用的固定數(shù)組：knode_iblk(numblk)：子區(qū)域節(jié)點總數(shù)n_update(numblk)：子區(qū)域內部節(jié)點數(shù)mglnod(knode)：整體節(jié)點編號對應的子區(qū)域節(jié)點編號mlgnod(numblk,max_lknode)：子區(qū)域節(jié)點編號對應的整體節(jié)點編號iorder(numblk,max_lknode)：子區(qū)域節(jié)點標識，1為內部節(jié)點，0為擴展節(jié)點主要用于解的轉換和整理算例：三維熱固耦合問題

熱固耦合問題模型控制方程熱傳導方程彈性體方程其中為彈性體位移，為作用于彈性體上的體積力。分別表示應力和應變向量。A是一個矩陣微分算子D為彈性矩陣Q為溫度場產(chǎn)生的熱應力向量，由下式確定：原型文件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