指南：從稀疏矩陣裝配到并行求解的有限元開發(fā))
簡介PETSc-FEM是一套基于PETSc并行科學(xué)計算庫的有限元求解代碼面向需要處理大規(guī)模偏微分方程問題的科研人員與工程師。它支持線性、多項式及高階有限元空間可靈活應(yīng)對不同復(fù)雜度的幾何結(jié)構(gòu)與物理過程集成多重網(wǎng)格、AMG等預(yù)處理技術(shù)以及GMRES、BiCGStab等高效Krylov迭代求解器保證并行計算性能。求解結(jié)果支持VTK格式輸出可直接配合ParaView、VisIt進行后處理分析代碼提供C、C和Fortran接口并配有示例與文檔便于用戶集成和快速上手。壓縮包為tgz格式大小約12.83MB內(nèi)含源代碼、編譯腳本、示例問題、測試案例與配套文檔方便使用者快速了解代碼結(jié)構(gòu)、運行流程并復(fù)用核心模塊。已有198人學(xué)習(xí)適合流體力學(xué)、固體力學(xué)、地球物理等領(lǐng)域的數(shù)值模擬開發(fā)者借助該代碼可快速搭建并驗證自己的并行有限元求解流程。 搞有限元仿真的人十有八九都經(jīng)歷過這種階段網(wǎng)格畫好了單元剛度矩陣推導(dǎo)出來了結(jié)果一到“裝配、求解、并行”這一層就卡住了。自己寫線性求解器要么是收斂慢得離譜要么是并行一上進程就開始互相踩內(nèi)存。后來我接觸到PETScFEM這一套基于PETSc的有限元開源代碼才真正把精力從“怎么求解”里解放出來專注回“怎么建?！薄Ｈ绻阋舱跒橛邢拊绦虻牡讓忧蠼夂筒⑿袛U展發(fā)愁這篇文章就是奔著你的痛點去的。PETScFEM不是一個孤立的軟件它底層依賴PETSc這個科學(xué)計算界的老牌開源庫。PETSc提供了Vec、Mat、DM、SNES、TS這一整套抽象層剛好對應(yīng)有限元里的向量、稀疏矩陣、網(wǎng)格管理、非線性迭代和時間積分。換句話說你只需要把“單元層面的計算”寫好剩下的線性代數(shù)、并行通信、預(yù)條件處理全部交給底層這比自己從零開始寫一套求解器要靠譜得多。這篇文章我會從選型思路、編譯部署、核心代碼拆解、非線性與瞬態(tài)擴展再到我實際踩過的坑完整走一遍。1. PETScFEM到底解決什么問題1.1 有限元編程最痛苦的那三件事如果你自己從零寫過有限元程序一定有過這種感受單元剛度矩陣的推導(dǎo)是體力活但真正讓人崩潰的是后面三件事。第一件是稀疏矩陣的存儲和組裝。有限元裝配出來的矩陣是高度稀疏的但到底用CSR還是COO行索引怎么排非零元預(yù)分配多少這些細節(jié)直接決定一個十萬自由度的問題是要幾秒還是要幾分鐘。很多人第一次寫有限元都栽在矩陣裝配的效率上。第二件是線性方程組的求解。直接法在幾百萬元度下還能湊合一旦網(wǎng)格細化到千萬級稀疏直接法內(nèi)存直接爆炸。換迭代法的話KSP選什么、預(yù)處理器用什么完全是一個經(jīng)驗活沒有高人指點很難調(diào)通。第三件是并行化。OpenMP一開線程競爭就把性能打回原形MPI一上邊界節(jié)點的通信邏輯寫到你懷疑人生。而PETScFEM把這三大痛點全部封裝好了你要做的只是描述“單元是什么”剩下的全局矩陣裝配、分布式存儲、通信和求解全部由庫替你完成。1.2 PETSc生態(tài)里的“有限元四件套”PETSc本身不是一個有限元庫它更像是一個科學(xué)計算工具箱。針對有限元它提供了四個核心組件剛好對應(yīng)完整求解流程的每一步。第一個是Vec和Mat分別是并行向量和并行稀疏矩陣。你在單元里算出來的局部剛度矩陣通過MatSetValues加到全局矩陣上PETSc自動處理跨進程的索引映射不需要你自己管節(jié)點編號的全局約定。第二個是DM也就是網(wǎng)格管理對象。DMAbstract是基類DMPlex是處理非結(jié)構(gòu)化網(wǎng)格的主力。DMPlex可以用DMPlexCreateFromFile直接讀入Gmsh或Exodus格式的網(wǎng)格自動完成網(wǎng)格分塊、節(jié)點重排序、邊界標記等處理。第三個是SNES非線性求解器。它把牛頓法的整體框架幫你搭好了你只需要提供殘量函數(shù)和Jacobian組裝函數(shù)甚至Jacobian可以用有限差分近似。對于材料非線性、幾何非線性這類問題SNES配合線搜索和信賴域策略比你自己寫迭代穩(wěn)太多。第四個是TS時間積分器。從最簡單的歐拉法到BDF、RK、Crank-NicolsonTS都內(nèi)置了而且支持自適應(yīng)步長。你只需要給出質(zhì)量矩陣和右端項時間步進的事情它全包。PETScFEM就是把這些組件串起來形成一個面向有限元的開發(fā)框架讓非數(shù)值計算背景的人也能寫出可擴展的并行有限元程序。2. 從源碼編譯到跑通第一個例程2.1 環(huán)境準備與configure要點我建議你第一次接觸PETScFEM時千萬別自己手動從官網(wǎng)逐個下載依賴直接用PETSc的下載腳本一把梭。PETSc極大簡化了構(gòu)建流程它會自動拉取BLAS、LAPACK、MPI、HDF5、Metis等底層庫。git clone -b release https://gitlab.com/petsc/petsc.git petsc cd petsc ./configure --with-debugging0 --with-ccmpicc --with-cxxmpicxx --with-fcmpif90 \ --download-mpich --download-fblaslapack --download-hdf5 --download-metis --download-parmetis \ --download-superlu_dist --download-mumps make PETSC_DIR$PWD PETSC_ARCHarch-linux-c-opt all幾個參數(shù)需要你特別留意--with-debugging0編譯的是優(yōu)化版速度快幾倍調(diào)試階段也可以用debug版便于定位問題--download-superlu_dist和--download-mumps這兩個稀疏直接法求解器對中小規(guī)模問題很友好當?shù)ú皇諗繒r它們是堅強的后盾。編譯過程大概需要十五到三十分鐘根據(jù)機器性能不同有差異。我看到很多人在這一步卡住多半是網(wǎng)絡(luò)問題導(dǎo)致依賴下載失敗建議提前配好代理或鏡像源。此外如果機器上沒有MPI環(huán)境直接讓PETSc自動下載MPICH是最省事的方式避免系統(tǒng)自帶MPI和PETSc編譯選項沖突。2.2 最小可運行示例穩(wěn)態(tài)熱傳導(dǎo)跑通一個最簡單的穩(wěn)態(tài)熱傳導(dǎo)問題是理解PETScFEM工作流的最佳路徑。這里我給出一個核心代碼骨架基于PETSc的C API完整的可編譯代碼建議直接參考PETSc自帶的示例。#include petscdmplex.h #include petscsnes.h #include petscds.h int main(int argc, char **argv) { PetscInitialize(argc, argv, NULL, NULL); DM dm; Vec x; Mat J; SNES snes; // 1. 創(chuàng)建并讀入網(wǎng)格 DMPlexCreateFromFile(PETSC_COMM_WORLD, mesh.msh, PETSC_TRUE, dm); DMPlexDistribute(dm, 0, NULL, dm); DMSetFromOptions(dm); DMSetUp(dm); // 2. 創(chuàng)建有限元離散對象 PetscDS ds; DMSetField(dm, 0, NULL, (PetscObject)NULL, 1); DMCreateDS(dm, ds); // 3. 創(chuàng)建非線性求解器和向量 DMCreateGlobalVector(dm, x); DMSetSolution(dm, x); DMCreateMatrix(dm, J); SNESCreate(PETSC_COMM_WORLD, snes); SNESSetDM(snes, dm); SNESSetApplicationContext(snes, NULL); // 4. 求解 SNESSolve(snes, NULL, x); // 5. 輸出結(jié)果 PetscViewer v; PetscViewerASCIIOpen(PETSC_COMM_WORLD, solution.vtu, v); VecView(x, v); PetscViewerDestroy(v); SNESDestroy(snes); DMDestroy(dm); PetscFinalize(); return 0; }這里你不需要理解每一行的細節(jié)但要注意整個流程的模式創(chuàng)建DM讀網(wǎng)格、創(chuàng)建DS管理離散化、創(chuàng)建SNES迭代、求解、輸出。后續(xù)做任何問題基本都是在這個骨架上替換“弱形式定義”的部分。實際運行時只要寫好網(wǎng)格文件用一行命令就能啟動mpiexec -n 4 ./heat -dm_plex_filename mesh.msh -petscspace_degree 1 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps看到輸出里SNES迭代收斂并且VTU文件生成恭喜你第一個并行有限元程序已經(jīng)跑起來了。3. 核心計算流程的代碼拆解3.1 網(wǎng)格讀入DMPlex工作模式DMPlex是PETSc管理非結(jié)構(gòu)化網(wǎng)格的主力。相比于傳統(tǒng)的數(shù)據(jù)結(jié)構(gòu)它對有限元最友好的地方在于所有實體頂點、邊、面、體統(tǒng)一抽象成“點”并通過錐和支撐的關(guān)系描述拓撲連接。這種設(shè)計帶來幾個直接好處。第一網(wǎng)格加密和自適應(yīng)時新生成的單元不需要重建整個數(shù)據(jù)結(jié)構(gòu)只修改局部拓撲關(guān)系。第二并行分區(qū)時DMPlex利用Metis或Parmetis自動完成圖劃分把網(wǎng)格分布到各進程同時記錄交界面的ghost單元。第三邊界條件通過DMPlexGetLabel獲取邊界編號你不用手工遍歷所有面判斷邊界類型。讀入網(wǎng)格時我推薦優(yōu)先使用Gmsh的msh格式PETSc對它的支持最成熟。需要注意Gmsh保存時盡量選擇ASCII格式版本兼容性最好。網(wǎng)格文件里的物理分組Physical Tag會直接轉(zhuǎn)換為DMPlex的Label后續(xù)設(shè)置Dirichlet邊界或Neumann邊界時直接用標簽名引用即可。如果你的網(wǎng)格里面含有高階單元比如二階三角形在DMPlex里也可以通過-petscspace_degree 2直接控制基函數(shù)階次不需要額外處理幾何節(jié)點的重排。這個東西對非結(jié)構(gòu)化網(wǎng)格來說很省心。3.2 弱形式定義與單元裝配PETScFEM里最核心也最容易勸退新手的部分是PetscDS。它的作用是把你的偏微分方程弱形式告訴PETSc然后由PETSc在你提供的單元上自動計算單元剛度矩陣和殘量。穩(wěn)態(tài)泊松方程是理解這個過程的最佳例子。弱形式為∫ Ω ?v · ?u dΩ ∫ Ω v f dΩ在PetscDS里你通過PetscDSAddBoundary和殘量函數(shù)定義來實現(xiàn)。殘量函數(shù)接收每個積分點的場值、梯度以及測試函數(shù)基函數(shù)返回局部殘量向量。對線性問題你也可以直接提供Jacobian組裝函數(shù)對非線性問題則必須提供。static PetscErrorCode f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[]) { f0[0] x[0] x[1]; /* 源項 f x y */ return PETSC_SUCCESS; }你不需要管單元循環(huán)、積分點坐標、映射到全局自由度這些事情。DMPlex和PetscDS自動遍歷所有單元調(diào)用你定義的函數(shù)將局部矩陣累加到全局矩陣。這個過程對分布式內(nèi)存是完全透明的每個進程只處理自己分到的單元交叉邊界處由PETSc的Communicator機制自動做全局組裝。3.3 線性求解KSP與預(yù)處理器選型關(guān)鍵點線性求解器是整個有限元計算最容易出問題、也最值得調(diào)整的地方。很多線性問題裝配矩陣很快但KSP迭代不收斂程序一直卡在Linear solve did not converge due to DIVERGED_INDEFINITE_PC這類錯誤上。我的經(jīng)驗是遵循一個簡單的選型規(guī)則二維中小規(guī)模問題自由度50萬直接用直接法最省心-ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps。MUMPS的內(nèi)存和魯棒性表現(xiàn)優(yōu)異。大規(guī)模橢圓型問題用共軛梯度法配合超節(jié)點分解預(yù)處理器-ksp_type cg -pc_type gamg。GAMG對橢圓算子效果極其優(yōu)秀接近最優(yōu)。大規(guī)模非對稱問題如對流占優(yōu)的對流擴散首選GMRES加ILU-ksp_type gmres -pc_type ilu。ILU的填充級別可以調(diào)-pc_factor_levels 2是常見的起點。這里最關(guān)鍵的一條經(jīng)驗是永遠不要用默認的KSP和PC配置去跑復(fù)雜問題默認的KSPGMRES PCJACOBI在稍大一點的問題上幾乎必然失敗。每次求解前先花五分鐘測試不同KSP/PC組合的收斂表現(xiàn)這個時間遠比你后面排查發(fā)散問題省得多。4. 向真實問題擴展非線性與時間相關(guān)計算4.1 非線性問題用SNES真實工程問題極少是線性的。材料非線性、大變形、接觸等都會讓控制方程變成F(u)0的形式。PETSc的SNES你用三個函數(shù)就能套住殘量函數(shù)、Jacobian函數(shù)、單調(diào)性初值。殘量函數(shù)的形式和線性問題中定義弱形式類似但f0返回的是非線性殘量。Jacobian函數(shù)需要額外輸出矩陣J?F/?u。好消息是Jacobian可以不精確提供通過SNESSetJacobian加上-snes_fd_color可以自動用有限差分近似對于單元數(shù)量少的問題非常好用。SNESSetFunction(snes, r, FormResidual, ctx); SNESSetJacobian(snes, J, J, FormJacobian, ctx);這里我對你的建議是不要一上來就追求解析Jacobian先用顏色有限差分跑通流程確認離散化沒有錯誤后再逐步替換成解析式。解析Jacobian能顯著提升收斂速度但手寫推導(dǎo)容易出錯得不償失。SNES的收斂控制推薦顯式設(shè)置-snes_type newtonls -snes_linesearch_type bt -snes_max_it 50 -snes_atol 1e-10 -snes_rtol 1e-8遇到不收斂時優(yōu)先檢查殘量初值是否合理。很多時候不是求解器的問題而是初始猜測給得太差。4.2 瞬態(tài)問題用TS對含時間的偏微分方程PETSc的TS模塊會讓你的生活輕松非常多。你不需要手動寫出時間步進循環(huán)只需要提供空間離散后的右端項函數(shù)和質(zhì)量矩陣。對熱傳導(dǎo)這類拋物線問題典型的TS配置如下-ts_type beuler -ts_dt 0.01 -ts_max_time 1.0 -ts_monitor對應(yīng)代碼里只需要TSCreate(PETSC_COMM_WORLD, ts); TSSetDM(ts, dm); TSSetProblemType(ts, TS_LINEAR); TSSetRHSFunction(ts, NULL, FormRHS, ctx); TSSetRHSJacobian(ts, J, J, FormJacobian, ctx); TSSetTimeStep(ts, 0.01); TSSetMaxTime(ts, 1.0);TS最大的優(yōu)勢是自適應(yīng)時間步。對剛性問題-ts_type rosw或-ts_type bdf配合-ts_adapt_type basic會自動調(diào)節(jié)步長在保證穩(wěn)定性的前提下盡量走大步長這個特性在反應(yīng)擴散、流固耦合這類多時間尺度問題中尤為重要。需要特別提醒的是TS對質(zhì)量矩陣是半隱式處理還是顯式處理取決于你提供的RHS Jacobian里是否包含質(zhì)量項。這類細節(jié)最容易讓人困惑建議先從最簡單的BEuler開始跑通后再換高階方法。4.3 并行策略與性能調(diào)優(yōu)PETScFEM的并行能力是它最大的賣點之一但并行效率并非天然好需要你注意幾個使用層面的事項。首先網(wǎng)格文件分區(qū)。如果你用Gmsh生成了一個大網(wǎng)格首次加載時PETSc會做一次分區(qū)平衡但這個過程默認只做一次。建議用DMPlexDistribute顯式調(diào)用分區(qū)并通過-dm_plex_partition_overlap 1增加ghost層的重疊避免邊界通信成為瓶頸。其次矩陣預(yù)分配。PETSc做矩陣裝配時如果非零元預(yù)分配不夠會觸發(fā)動態(tài)增加這一步會顯著拖慢性能并且?guī)韮?nèi)存碎片。通過DMCreateMatrix配合-info查看裝配階段日志確認MatAssemblyBegin沒有大量reallocation告警。如果有使用MatSetOption和MatMPIAIJSetPreallocation手動預(yù)分配。第三負載均衡。非結(jié)構(gòu)化網(wǎng)格很容易出現(xiàn)某些進程分到的單元計算量大、某些進程很閑的情況。Parmetis的分區(qū)效果通常比Metis好因為它會考慮進程間的通信開銷。在configure時我就建議你下載parmetis原因就在這里。實測下來一個百萬自由度的線彈性問題4進程并行時加速比能達到3.5左右16進程能到12以上前提是預(yù)處理器的設(shè)置跟得上。如果發(fā)現(xiàn)擴展性上不去優(yōu)先檢查是否是KSP內(nèi)部全局通信太多換用-ksp_type cg和-pc_type gamg這類局部預(yù)處理組合能有效改善。5. 踩坑記錄與排查技巧5.1 編譯期高頻錯誤我見過太多人在編譯階段就卡住這里把幾個高頻問題列出來方便你對照查找。依賴下載失敗configure時--download-*聯(lián)網(wǎng)失敗多數(shù)是網(wǎng)絡(luò)策略導(dǎo)致建議手動下載源碼包放進petsc/目錄下的packages文件夾PETSc會自動識別并跳過下載。MPI版本沖突系統(tǒng)MPI和PETSc自帶MPICH混用會導(dǎo)致頭文件不匹配。最穩(wěn)妥的做法是全程只用PETSc自動下載的MPICH不要動系統(tǒng)MPI。32位索引溢出問題規(guī)模超過20億自由度才會遇到但如果你要跑超大算例需要在configure時加上--with-64-bit-indices否則索引溢出會出現(xiàn)莫名其妙的內(nèi)存崩潰。編譯期間遇到報錯第一件事不是去查錯誤碼而是打開$PETSC_DIR/$PETSC_ARCH/lib/petsc/conf/petscvariables確認編譯選項是否符合預(yù)期。很多時候是路徑或環(huán)境變量沒對齊。5.2 求解不收斂的排查清單求解器不收斂是有限元最消耗精力的環(huán)節(jié)我整理了一個排查順序按這個順序查通常能定位問題。第一檢查矩陣是否對稱正定。對彈性力學(xué)、熱傳導(dǎo)這類問題理論上矩陣一定對稱正定。如果你裝配結(jié)束后用-ksp_monitor_true_residual觀察殘差曲線發(fā)現(xiàn)殘差震蕩或發(fā)散先確認是不是有約束沒處理干凈導(dǎo)致矩陣奇異。第二檢查邊界條件是否齊全。結(jié)構(gòu)問題如果沒有施加足夠的位移約束矩陣就會奇異CG法在第一次迭代就會因為DIVERGED_INDEFINITE_MAT退出。解決方法是先檢查網(wǎng)格Label里邊界標記是否正確傳給了PetscDS。第三檢查預(yù)處理器是否需要調(diào)換。ILU在矩陣條件數(shù)很差的情況下不穩(wěn)定換成GAMG或直接法先驗證“離散本身是否沒問題”如果直接法能收斂而迭代法不行問題基本出在PC選型而不是模型定義。第四檢查網(wǎng)格質(zhì)量?；儐卧獣?dǎo)致單元Jacobian為負或接近零直接破壞全局矩陣性質(zhì)。對這類問題先跑Gmsh的Mesh-Check確認最小單元質(zhì)量及時修復(fù)網(wǎng)格。5.3 并行調(diào)試的幾個實用技巧并行程序調(diào)試比串行麻煩一個量級有幾個技巧能讓你的頭發(fā)少掉一些。調(diào)試時先把進程數(shù)設(shè)為1確認串行無誤后再上并行。用-start_in_debugger可以在每個進程啟動時掛起調(diào)試器但更實用的方式是在關(guān)鍵函數(shù)里加PetscPrintf打印當前進程的局部信息再通過-log_view查看各進程的耗時分布。另一個常見問題是結(jié)果在串行正確、并行卻出錯這幾乎總是邊界通信或全局索引映射的問題。在PETSc里用DMPlexGetFullData檢查每個進程的本地網(wǎng)格和ghost數(shù)據(jù)或者用VecView輸出分片向量看邊界節(jié)點的值是否在相鄰進程間一致。還有一個隱藏問題隨機數(shù)或全局變量。有限元程序里如果用了全局累積的計數(shù)器或者非確定性隨機數(shù)并行時不同進程的計算順序不同結(jié)果就會有微小差異這類bug極難復(fù)現(xiàn)建議從一開始就避免在單元計算中使用全局狀態(tài)。6. 個人實操中的幾點額外體會最后再說幾個我自己的習(xí)慣。第一個是全套遵循PETScFEM的命名和管理方式不要把Fortran和C的接口混著用。PETSc的C接口和Fortran接口雖然底層一樣但錯誤處理方式不同混用很容易在出錯時拿不到完整調(diào)用棧。第二個是善用PETSc的-info和-log_view。很多人只在出問題時才開日志平時關(guān)掉節(jié)省輸出。我的建議是任何新算例第一次跑至少開一次-log_view看看矩陣裝配耗時、求解器迭代次數(shù)、每一步的耗時占比這些信息能直接告訴你優(yōu)化方向。第三個是保留調(diào)試版的PETSc構(gòu)建。調(diào)試版速度慢但能捕獲內(nèi)存越界、未初始化變量這類隱蔽錯誤。我一般會同時保留arch-linux-c-opt和arch-linux-c-debug兩套構(gòu)建開發(fā)調(diào)試用debug正式跑算例用opt。還有一點是關(guān)于社區(qū)資源。PETSc官方文檔和示例非常豐富$PETSC_DIR/src/snes/tutorials和$PETSC_DIR/src/ts/tutorials里面幾乎有你能想到的所有經(jīng)典算例從線性彈性到不可壓縮流體都有。遇到問題先翻tutorials再上GitLab提issue比你在搜索引擎里瞎找效率高得多。PETScFEM這條路前期確實需要一點耐心主要是概念多、接口抽象層次高但一旦你把第一個問題完整跑通之后再換模型、換求解器、上并行就會順暢非常多。這篇分享里的所有配置和調(diào)試經(jīng)驗都是我在實際項目中一條一條總結(jié)出來的希望它幫你少走幾個彎路。本文還有配套的精品資源點擊獲取 PETScFEM實戰(zhàn)指南：從稀疏矩陣裝配到并行求解的有限元開發(fā)
簡介PETSc-FEM是一套基于PETSc并行科學(xué)計算庫的有限元求解代碼面向需要處理大規(guī)模偏微分方程問題的科研人員與工程師。它支持線性、多項式及高階有限元空間可靈活應(yīng)對不同復(fù)雜度的幾何結(jié)構(gòu)與物理過程集成多重網(wǎng)格、AMG等預(yù)處理技術(shù)以及GMRES、BiCGStab等高效Krylov迭代求解器保證并行計算性能。求解結(jié)果支持VTK格式輸出可直接配合ParaView、VisIt進行后處理分析代碼提供C、C和Fortran接口并配有示例與文檔便于用戶集成和快速上手。壓縮包為tgz格式大小約12.83MB內(nèi)含源代碼、編譯腳本、示例問題、測試案例與配套文檔方便使用者快速了解代碼結(jié)構(gòu)、運行流程并復(fù)用核心模塊。已有198人學(xué)習(xí)適合流體力學(xué)、固體力學(xué)、地球物理等領(lǐng)域的數(shù)值模擬開發(fā)者借助該代碼可快速搭建并驗證自己的并行有限元求解流程。 搞有限元仿真的人十有八九都經(jīng)歷過這種階段網(wǎng)格畫好了單元剛度矩陣推導(dǎo)出來了結(jié)果一到“裝配、求解、并行”這一層就卡住了。自己寫線性求解器要么是收斂慢得離譜要么是并行一上進程就開始互相踩內(nèi)存。后來我接觸到PETScFEM這一套基于PETSc的有限元開源代碼才真正把精力從“怎么求解”里解放出來專注回“怎么建?！薄Ｈ绻阋舱跒橛邢拊绦虻牡讓忧蠼夂筒⑿袛U展發(fā)愁這篇文章就是奔著你的痛點去的。PETScFEM不是一個孤立的軟件它底層依賴PETSc這個科學(xué)計算界的老牌開源庫。PETSc提供了Vec、Mat、DM、SNES、TS這一整套抽象層剛好對應(yīng)有限元里的向量、稀疏矩陣、網(wǎng)格管理、非線性迭代和時間積分。換句話說你只需要把“單元層面的計算”寫好剩下的線性代數(shù)、并行通信、預(yù)條件處理全部交給底層這比自己從零開始寫一套求解器要靠譜得多。這篇文章我會從選型思路、編譯部署、核心代碼拆解、非線性與瞬態(tài)擴展再到我實際踩過的坑完整走一遍。1. PETScFEM到底解決什么問題1.1 有限元編程最痛苦的那三件事如果你自己從零寫過有限元程序一定有過這種感受單元剛度矩陣的推導(dǎo)是體力活但真正讓人崩潰的是后面三件事。第一件是稀疏矩陣的存儲和組裝。有限元裝配出來的矩陣是高度稀疏的但到底用CSR還是COO行索引怎么排非零元預(yù)分配多少這些細節(jié)直接決定一個十萬自由度的問題是要幾秒還是要幾分鐘。很多人第一次寫有限元都栽在矩陣裝配的效率上。第二件是線性方程組的求解。直接法在幾百萬元度下還能湊合一旦網(wǎng)格細化到千萬級稀疏直接法內(nèi)存直接爆炸。換迭代法的話KSP選什么、預(yù)處理器用什么完全是一個經(jīng)驗活沒有高人指點很難調(diào)通。第三件是并行化。OpenMP一開線程競爭就把性能打回原形MPI一上邊界節(jié)點的通信邏輯寫到你懷疑人生。而PETScFEM把這三大痛點全部封裝好了你要做的只是描述“單元是什么”剩下的全局矩陣裝配、分布式存儲、通信和求解全部由庫替你完成。1.2 PETSc生態(tài)里的“有限元四件套”PETSc本身不是一個有限元庫它更像是一個科學(xué)計算工具箱。針對有限元它提供了四個核心組件剛好對應(yīng)完整求解流程的每一步。第一個是Vec和Mat分別是并行向量和并行稀疏矩陣。你在單元里算出來的局部剛度矩陣通過MatSetValues加到全局矩陣上PETSc自動處理跨進程的索引映射不需要你自己管節(jié)點編號的全局約定。第二個是DM也就是網(wǎng)格管理對象。DMAbstract是基類DMPlex是處理非結(jié)構(gòu)化網(wǎng)格的主力。DMPlex可以用DMPlexCreateFromFile直接讀入Gmsh或Exodus格式的網(wǎng)格自動完成網(wǎng)格分塊、節(jié)點重排序、邊界標記等處理。第三個是SNES非線性求解器。它把牛頓法的整體框架幫你搭好了你只需要提供殘量函數(shù)和Jacobian組裝函數(shù)甚至Jacobian可以用有限差分近似。對于材料非線性、幾何非線性這類問題SNES配合線搜索和信賴域策略比你自己寫迭代穩(wěn)太多。第四個是TS時間積分器。從最簡單的歐拉法到BDF、RK、Crank-NicolsonTS都內(nèi)置了而且支持自適應(yīng)步長。你只需要給出質(zhì)量矩陣和右端項時間步進的事情它全包。PETScFEM就是把這些組件串起來形成一個面向有限元的開發(fā)框架讓非數(shù)值計算背景的人也能寫出可擴展的并行有限元程序。2. 從源碼編譯到跑通第一個例程2.1 環(huán)境準備與configure要點我建議你第一次接觸PETScFEM時千萬別自己手動從官網(wǎng)逐個下載依賴直接用PETSc的下載腳本一把梭。PETSc極大簡化了構(gòu)建流程它會自動拉取BLAS、LAPACK、MPI、HDF5、Metis等底層庫。git clone -b release https://gitlab.com/petsc/petsc.git petsc cd petsc ./configure --with-debugging0 --with-ccmpicc --with-cxxmpicxx --with-fcmpif90 \ --download-mpich --download-fblaslapack --download-hdf5 --download-metis --download-parmetis \ --download-superlu_dist --download-mumps make PETSC_DIR$PWD PETSC_ARCHarch-linux-c-opt all幾個參數(shù)需要你特別留意--with-debugging0編譯的是優(yōu)化版速度快幾倍調(diào)試階段也可以用debug版便于定位問題--download-superlu_dist和--download-mumps這兩個稀疏直接法求解器對中小規(guī)模問題很友好當?shù)ú皇諗繒r它們是堅強的后盾。編譯過程大概需要十五到三十分鐘根據(jù)機器性能不同有差異。我看到很多人在這一步卡住多半是網(wǎng)絡(luò)問題導(dǎo)致依賴下載失敗建議提前配好代理或鏡像源。此外如果機器上沒有MPI環(huán)境直接讓PETSc自動下載MPICH是最省事的方式避免系統(tǒng)自帶MPI和PETSc編譯選項沖突。2.2 最小可運行示例穩(wěn)態(tài)熱傳導(dǎo)跑通一個最簡單的穩(wěn)態(tài)熱傳導(dǎo)問題是理解PETScFEM工作流的最佳路徑。這里我給出一個核心代碼骨架基于PETSc的C API完整的可編譯代碼建議直接參考PETSc自帶的示例。#include petscdmplex.h #include petscsnes.h #include petscds.h int main(int argc, char **argv) { PetscInitialize(argc, argv, NULL, NULL); DM dm; Vec x; Mat J; SNES snes; // 1. 創(chuàng)建并讀入網(wǎng)格 DMPlexCreateFromFile(PETSC_COMM_WORLD, mesh.msh, PETSC_TRUE, dm); DMPlexDistribute(dm, 0, NULL, dm); DMSetFromOptions(dm); DMSetUp(dm); // 2. 創(chuàng)建有限元離散對象 PetscDS ds; DMSetField(dm, 0, NULL, (PetscObject)NULL, 1); DMCreateDS(dm, ds); // 3. 創(chuàng)建非線性求解器和向量 DMCreateGlobalVector(dm, x); DMSetSolution(dm, x); DMCreateMatrix(dm, J); SNESCreate(PETSC_COMM_WORLD, snes); SNESSetDM(snes, dm); SNESSetApplicationContext(snes, NULL); // 4. 求解 SNESSolve(snes, NULL, x); // 5. 輸出結(jié)果 PetscViewer v; PetscViewerASCIIOpen(PETSC_COMM_WORLD, solution.vtu, v); VecView(x, v); PetscViewerDestroy(v); SNESDestroy(snes); DMDestroy(dm); PetscFinalize(); return 0; }這里你不需要理解每一行的細節(jié)但要注意整個流程的模式創(chuàng)建DM讀網(wǎng)格、創(chuàng)建DS管理離散化、創(chuàng)建SNES迭代、求解、輸出。后續(xù)做任何問題基本都是在這個骨架上替換“弱形式定義”的部分。實際運行時只要寫好網(wǎng)格文件用一行命令就能啟動mpiexec -n 4 ./heat -dm_plex_filename mesh.msh -petscspace_degree 1 -ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps看到輸出里SNES迭代收斂并且VTU文件生成恭喜你第一個并行有限元程序已經(jīng)跑起來了。3. 核心計算流程的代碼拆解3.1 網(wǎng)格讀入DMPlex工作模式DMPlex是PETSc管理非結(jié)構(gòu)化網(wǎng)格的主力。相比于傳統(tǒng)的數(shù)據(jù)結(jié)構(gòu)它對有限元最友好的地方在于所有實體頂點、邊、面、體統(tǒng)一抽象成“點”并通過錐和支撐的關(guān)系描述拓撲連接。這種設(shè)計帶來幾個直接好處。第一網(wǎng)格加密和自適應(yīng)時新生成的單元不需要重建整個數(shù)據(jù)結(jié)構(gòu)只修改局部拓撲關(guān)系。第二并行分區(qū)時DMPlex利用Metis或Parmetis自動完成圖劃分把網(wǎng)格分布到各進程同時記錄交界面的ghost單元。第三邊界條件通過DMPlexGetLabel獲取邊界編號你不用手工遍歷所有面判斷邊界類型。讀入網(wǎng)格時我推薦優(yōu)先使用Gmsh的msh格式PETSc對它的支持最成熟。需要注意Gmsh保存時盡量選擇ASCII格式版本兼容性最好。網(wǎng)格文件里的物理分組Physical Tag會直接轉(zhuǎn)換為DMPlex的Label后續(xù)設(shè)置Dirichlet邊界或Neumann邊界時直接用標簽名引用即可。如果你的網(wǎng)格里面含有高階單元比如二階三角形在DMPlex里也可以通過-petscspace_degree 2直接控制基函數(shù)階次不需要額外處理幾何節(jié)點的重排。這個東西對非結(jié)構(gòu)化網(wǎng)格來說很省心。3.2 弱形式定義與單元裝配PETScFEM里最核心也最容易勸退新手的部分是PetscDS。它的作用是把你的偏微分方程弱形式告訴PETSc然后由PETSc在你提供的單元上自動計算單元剛度矩陣和殘量。穩(wěn)態(tài)泊松方程是理解這個過程的最佳例子。弱形式為∫ Ω ?v · ?u dΩ ∫ Ω v f dΩ在PetscDS里你通過PetscDSAddBoundary和殘量函數(shù)定義來實現(xiàn)。殘量函數(shù)接收每個積分點的場值、梯度以及測試函數(shù)基函數(shù)返回局部殘量向量。對線性問題你也可以直接提供Jacobian組裝函數(shù)對非線性問題則必須提供。static PetscErrorCode f0_u(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], const PetscReal x[], PetscScalar f0[]) { f0[0] x[0] x[1]; /* 源項 f x y */ return PETSC_SUCCESS; }你不需要管單元循環(huán)、積分點坐標、映射到全局自由度這些事情。DMPlex和PetscDS自動遍歷所有單元調(diào)用你定義的函數(shù)將局部矩陣累加到全局矩陣。這個過程對分布式內(nèi)存是完全透明的每個進程只處理自己分到的單元交叉邊界處由PETSc的Communicator機制自動做全局組裝。3.3 線性求解KSP與預(yù)處理器選型關(guān)鍵點線性求解器是整個有限元計算最容易出問題、也最值得調(diào)整的地方。很多線性問題裝配矩陣很快但KSP迭代不收斂程序一直卡在Linear solve did not converge due to DIVERGED_INDEFINITE_PC這類錯誤上。我的經(jīng)驗是遵循一個簡單的選型規(guī)則二維中小規(guī)模問題自由度50萬直接用直接法最省心-ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps。MUMPS的內(nèi)存和魯棒性表現(xiàn)優(yōu)異。大規(guī)模橢圓型問題用共軛梯度法配合超節(jié)點分解預(yù)處理器-ksp_type cg -pc_type gamg。GAMG對橢圓算子效果極其優(yōu)秀接近最優(yōu)。大規(guī)模非對稱問題如對流占優(yōu)的對流擴散首選GMRES加ILU-ksp_type gmres -pc_type ilu。ILU的填充級別可以調(diào)-pc_factor_levels 2是常見的起點。這里最關(guān)鍵的一條經(jīng)驗是永遠不要用默認的KSP和PC配置去跑復(fù)雜問題默認的KSPGMRES PCJACOBI在稍大一點的問題上幾乎必然失敗。每次求解前先花五分鐘測試不同KSP/PC組合的收斂表現(xiàn)這個時間遠比你后面排查發(fā)散問題省得多。4. 向真實問題擴展非線性與時間相關(guān)計算4.1 非線性問題用SNES真實工程問題極少是線性的。材料非線性、大變形、接觸等都會讓控制方程變成F(u)0的形式。PETSc的SNES你用三個函數(shù)就能套住殘量函數(shù)、Jacobian函數(shù)、單調(diào)性初值。殘量函數(shù)的形式和線性問題中定義弱形式類似但f0返回的是非線性殘量。Jacobian函數(shù)需要額外輸出矩陣J?F/?u。好消息是Jacobian可以不精確提供通過SNESSetJacobian加上-snes_fd_color可以自動用有限差分近似對于單元數(shù)量少的問題非常好用。SNESSetFunction(snes, r, FormResidual, ctx); SNESSetJacobian(snes, J, J, FormJacobian, ctx);這里我對你的建議是不要一上來就追求解析Jacobian先用顏色有限差分跑通流程確認離散化沒有錯誤后再逐步替換成解析式。解析Jacobian能顯著提升收斂速度但手寫推導(dǎo)容易出錯得不償失。SNES的收斂控制推薦顯式設(shè)置-snes_type newtonls -snes_linesearch_type bt -snes_max_it 50 -snes_atol 1e-10 -snes_rtol 1e-8遇到不收斂時優(yōu)先檢查殘量初值是否合理。很多時候不是求解器的問題而是初始猜測給得太差。4.2 瞬態(tài)問題用TS對含時間的偏微分方程PETSc的TS模塊會讓你的生活輕松非常多。你不需要手動寫出時間步進循環(huán)只需要提供空間離散后的右端項函數(shù)和質(zhì)量矩陣。對熱傳導(dǎo)這類拋物線問題典型的TS配置如下-ts_type beuler -ts_dt 0.01 -ts_max_time 1.0 -ts_monitor對應(yīng)代碼里只需要TSCreate(PETSC_COMM_WORLD, ts); TSSetDM(ts, dm); TSSetProblemType(ts, TS_LINEAR); TSSetRHSFunction(ts, NULL, FormRHS, ctx); TSSetRHSJacobian(ts, J, J, FormJacobian, ctx); TSSetTimeStep(ts, 0.01); TSSetMaxTime(ts, 1.0);TS最大的優(yōu)勢是自適應(yīng)時間步。對剛性問題-ts_type rosw或-ts_type bdf配合-ts_adapt_type basic會自動調(diào)節(jié)步長在保證穩(wěn)定性的前提下盡量走大步長這個特性在反應(yīng)擴散、流固耦合這類多時間尺度問題中尤為重要。需要特別提醒的是TS對質(zhì)量矩陣是半隱式處理還是顯式處理取決于你提供的RHS Jacobian里是否包含質(zhì)量項。這類細節(jié)最容易讓人困惑建議先從最簡單的BEuler開始跑通后再換高階方法。4.3 并行策略與性能調(diào)優(yōu)PETScFEM的并行能力是它最大的賣點之一但并行效率并非天然好需要你注意幾個使用層面的事項。首先網(wǎng)格文件分區(qū)。如果你用Gmsh生成了一個大網(wǎng)格首次加載時PETSc會做一次分區(qū)平衡但這個過程默認只做一次。建議用DMPlexDistribute顯式調(diào)用分區(qū)并通過-dm_plex_partition_overlap 1增加ghost層的重疊避免邊界通信成為瓶頸。其次矩陣預(yù)分配。PETSc做矩陣裝配時如果非零元預(yù)分配不夠會觸發(fā)動態(tài)增加這一步會顯著拖慢性能并且?guī)韮?nèi)存碎片。通過DMCreateMatrix配合-info查看裝配階段日志確認MatAssemblyBegin沒有大量reallocation告警。如果有使用MatSetOption和MatMPIAIJSetPreallocation手動預(yù)分配。第三負載均衡。非結(jié)構(gòu)化網(wǎng)格很容易出現(xiàn)某些進程分到的單元計算量大、某些進程很閑的情況。Parmetis的分區(qū)效果通常比Metis好因為它會考慮進程間的通信開銷。在configure時我就建議你下載parmetis原因就在這里。實測下來一個百萬自由度的線彈性問題4進程并行時加速比能達到3.5左右16進程能到12以上前提是預(yù)處理器的設(shè)置跟得上。如果發(fā)現(xiàn)擴展性上不去優(yōu)先檢查是否是KSP內(nèi)部全局通信太多換用-ksp_type cg和-pc_type gamg這類局部預(yù)處理組合能有效改善。5. 踩坑記錄與排查技巧5.1 編譯期高頻錯誤我見過太多人在編譯階段就卡住這里把幾個高頻問題列出來方便你對照查找。依賴下載失敗configure時--download-*聯(lián)網(wǎng)失敗多數(shù)是網(wǎng)絡(luò)策略導(dǎo)致建議手動下載源碼包放進petsc/目錄下的packages文件夾PETSc會自動識別并跳過下載。MPI版本沖突系統(tǒng)MPI和PETSc自帶MPICH混用會導(dǎo)致頭文件不匹配。最穩(wěn)妥的做法是全程只用PETSc自動下載的MPICH不要動系統(tǒng)MPI。32位索引溢出問題規(guī)模超過20億自由度才會遇到但如果你要跑超大算例需要在configure時加上--with-64-bit-indices否則索引溢出會出現(xiàn)莫名其妙的內(nèi)存崩潰。編譯期間遇到報錯第一件事不是去查錯誤碼而是打開$PETSC_DIR/$PETSC_ARCH/lib/petsc/conf/petscvariables確認編譯選項是否符合預(yù)期。很多時候是路徑或環(huán)境變量沒對齊。5.2 求解不收斂的排查清單求解器不收斂是有限元最消耗精力的環(huán)節(jié)我整理了一個排查順序按這個順序查通常能定位問題。第一檢查矩陣是否對稱正定。對彈性力學(xué)、熱傳導(dǎo)這類問題理論上矩陣一定對稱正定。如果你裝配結(jié)束后用-ksp_monitor_true_residual觀察殘差曲線發(fā)現(xiàn)殘差震蕩或發(fā)散先確認是不是有約束沒處理干凈導(dǎo)致矩陣奇異。第二檢查邊界條件是否齊全。結(jié)構(gòu)問題如果沒有施加足夠的位移約束矩陣就會奇異CG法在第一次迭代就會因為DIVERGED_INDEFINITE_MAT退出。解決方法是先檢查網(wǎng)格Label里邊界標記是否正確傳給了PetscDS。第三檢查預(yù)處理器是否需要調(diào)換。ILU在矩陣條件數(shù)很差的情況下不穩(wěn)定換成GAMG或直接法先驗證“離散本身是否沒問題”如果直接法能收斂而迭代法不行問題基本出在PC選型而不是模型定義。第四檢查網(wǎng)格質(zhì)量?；儐卧獣?dǎo)致單元Jacobian為負或接近零直接破壞全局矩陣性質(zhì)。對這類問題先跑Gmsh的Mesh-Check確認最小單元質(zhì)量及時修復(fù)網(wǎng)格。5.3 并行調(diào)試的幾個實用技巧并行程序調(diào)試比串行麻煩一個量級有幾個技巧能讓你的頭發(fā)少掉一些。調(diào)試時先把進程數(shù)設(shè)為1確認串行無誤后再上并行。用-start_in_debugger可以在每個進程啟動時掛起調(diào)試器但更實用的方式是在關(guān)鍵函數(shù)里加PetscPrintf打印當前進程的局部信息再通過-log_view查看各進程的耗時分布。另一個常見問題是結(jié)果在串行正確、并行卻出錯這幾乎總是邊界通信或全局索引映射的問題。在PETSc里用DMPlexGetFullData檢查每個進程的本地網(wǎng)格和ghost數(shù)據(jù)或者用VecView輸出分片向量看邊界節(jié)點的值是否在相鄰進程間一致。還有一個隱藏問題隨機數(shù)或全局變量。有限元程序里如果用了全局累積的計數(shù)器或者非確定性隨機數(shù)并行時不同進程的計算順序不同結(jié)果就會有微小差異這類bug極難復(fù)現(xiàn)建議從一開始就避免在單元計算中使用全局狀態(tài)。6. 個人實操中的幾點額外體會最后再說幾個我自己的習(xí)慣。第一個是全套遵循PETScFEM的命名和管理方式不要把Fortran和C的接口混著用。PETSc的C接口和Fortran接口雖然底層一樣但錯誤處理方式不同混用很容易在出錯時拿不到完整調(diào)用棧。第二個是善用PETSc的-info和-log_view。很多人只在出問題時才開日志平時關(guān)掉節(jié)省輸出。我的建議是任何新算例第一次跑至少開一次-log_view看看矩陣裝配耗時、求解器迭代次數(shù)、每一步的耗時占比這些信息能直接告訴你優(yōu)化方向。第三個是保留調(diào)試版的PETSc構(gòu)建。調(diào)試版速度慢但能捕獲內(nèi)存越界、未初始化變量這類隱蔽錯誤。我一般會同時保留arch-linux-c-opt和arch-linux-c-debug兩套構(gòu)建開發(fā)調(diào)試用debug正式跑算例用opt。還有一點是關(guān)于社區(qū)資源。PETSc官方文檔和示例非常豐富$PETSC_DIR/src/snes/tutorials和$PETSC_DIR/src/ts/tutorials里面幾乎有你能想到的所有經(jīng)典算例從線性彈性到不可壓縮流體都有。遇到問題先翻tutorials再上GitLab提issue比你在搜索引擎里瞎找效率高得多。PETScFEM這條路前期確實需要一點耐心主要是概念多、接口抽象層次高但一旦你把第一個問題完整跑通之后再換模型、換求解器、上并行就會順暢非常多。這篇分享里的所有配置和調(diào)試經(jīng)驗都是我在實際項目中一條一條總結(jié)出來的希望它幫你少走幾個彎路。本文還有配套的精品資源點擊獲取