简介PETSc-FEM是一套基于PETSc并行科学计算库的有限元求解代码面向需要处理大规模偏微分方程问题的科研人员与工程师。它支持线性、多项式及高阶有限元空间可灵活应对不同复杂度的几何结构与物理过程集成多重网格、AMG等预处理技术以及GMRES、BiCGStab等高效Krylov迭代求解器保证并行计算性能。求解结果支持VTK格式输出可直接配合ParaView、VisIt进行后处理分析代码提供C、C和Fortran接口并配有示例与文档便于用户集成和快速上手。压缩包为tgz格式大小约12.83MB内含源代码、编译脚本、示例问题、测试案例与配套文档方便使用者快速了解代码结构、运行流程并复用核心模块。已有198人学习适合流体力学、固体力学、地球物理等领域的数值模拟开发者借助该代码可快速搭建并验证自己的并行有限元求解流程。 搞有限元仿真的人十有八九都经历过这种阶段网格画好了单元刚度矩阵推导出来了结果一到“装配、求解、并行”这一层就卡住了。自己写线性求解器要么是收敛慢得离谱要么是并行一上进程就开始互相踩内存。后来我接触到PETScFEM这一套基于PETSc的有限元开源代码才真正把精力从“怎么求解”里解放出来专注回“怎么建模”。如果你也正在为有限元程序的底层求解和并行扩展发愁这篇文章就是奔着你的痛点去的。PETScFEM不是一个孤立的软件它底层依赖PETSc这个科学计算界的老牌开源库。PETSc提供了Vec、Mat、DM、SNES、TS这一整套抽象层刚好对应有限元里的向量、稀疏矩阵、网格管理、非线性迭代和时间积分。换句话说你只需要把“单元层面的计算”写好剩下的线性代数、并行通信、预条件处理全部交给底层这比自己从零开始写一套求解器要靠谱得多。这篇文章我会从选型思路、编译部署、核心代码拆解、非线性与瞬态扩展再到我实际踩过的坑完整走一遍。1. PETScFEM到底解决什么问题1.1 有限元编程最痛苦的那三件事如果你自己从零写过有限元程序一定有过这种感受单元刚度矩阵的推导是体力活但真正让人崩溃的是后面三件事。第一件是稀疏矩阵的存储和组装。有限元装配出来的矩阵是高度稀疏的但到底用CSR还是COO行索引怎么排非零元预分配多少这些细节直接决定一个十万自由度的问题是要几秒还是要几分钟。很多人第一次写有限元都栽在矩阵装配的效率上。第二件是线性方程组的求解。直接法在几百万元度下还能凑合一旦网格细化到千万级稀疏直接法内存直接爆炸。换迭代法的话KSP选什么、预处理器用什么完全是一个经验活没有高人指点很难调通。第三件是并行化。OpenMP一开线程竞争就把性能打回原形MPI一上边界节点的通信逻辑写到你怀疑人生。而PETScFEM把这三大痛点全部封装好了你要做的只是描述“单元是什么”剩下的全局矩阵装配、分布式存储、通信和求解全部由库替你完成。1.2 PETSc生态里的“有限元四件套”PETSc本身不是一个有限元库它更像是一个科学计算工具箱。针对有限元它提供了四个核心组件刚好对应完整求解流程的每一步。第一个是Vec和Mat分别是并行向量和并行稀疏矩阵。你在单元里算出来的局部刚度矩阵通过MatSetValues加到全局矩阵上PETSc自动处理跨进程的索引映射不需要你自己管节点编号的全局约定。第二个是DM也就是网格管理对象。DMAbstract是基类DMPlex是处理非结构化网格的主力。DMPlex可以用DMPlexCreateFromFile直接读入Gmsh或Exodus格式的网格自动完成网格分块、节点重排序、边界标记等处理。第三个是SNES非线性求解器。它把牛顿法的整体框架帮你搭好了你只需要提供残量函数和Jacobian组装函数甚至Jacobian可以用有限差分近似。对于材料非线性、几何非线性这类问题SNES配合线搜索和信赖域策略比你自己写迭代稳太多。第四个是TS时间积分器。从最简单的欧拉法到BDF、RK、Crank-NicolsonTS都内置了而且支持自适应步长。你只需要给出质量矩阵和右端项时间步进的事情它全包。PETScFEM就是把这些组件串起来形成一个面向有限元的开发框架让非数值计算背景的人也能写出可扩展的并行有限元程序。2. 从源码编译到跑通第一个例程2.1 环境准备与configure要点我建议你第一次接触PETScFEM时千万别自己手动从官网逐个下载依赖直接用PETSc的下载脚本一把梭。PETSc极大简化了构建流程它会自动拉取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几个参数需要你特别留意--with-debugging0编译的是优化版速度快几倍调试阶段也可以用debug版便于定位问题--download-superlu_dist和--download-mumps这两个稀疏直接法求解器对中小规模问题很友好当迭代法不收敛时它们是坚强的后盾。编译过程大概需要十五到三十分钟根据机器性能不同有差异。我看到很多人在这一步卡住多半是网络问题导致依赖下载失败建议提前配好代理或镜像源。此外如果机器上没有MPI环境直接让PETSc自动下载MPICH是最省事的方式避免系统自带MPI和PETSc编译选项冲突。2.2 最小可运行示例稳态热传导跑通一个最简单的稳态热传导问题是理解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. 创建并读入网格 DMPlexCreateFromFile(PETSC_COMM_WORLD, mesh.msh, PETSC_TRUE, dm); DMPlexDistribute(dm, 0, NULL, dm); DMSetFromOptions(dm); DMSetUp(dm); // 2. 创建有限元离散对象 PetscDS ds; DMSetField(dm, 0, NULL, (PetscObject)NULL, 1); DMCreateDS(dm, ds); // 3. 创建非线性求解器和向量 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. 输出结果 PetscViewer v; PetscViewerASCIIOpen(PETSC_COMM_WORLD, solution.vtu, v); VecView(x, v); PetscViewerDestroy(v); SNESDestroy(snes); DMDestroy(dm); PetscFinalize(); return 0; }这里你不需要理解每一行的细节但要注意整个流程的模式创建DM读网格、创建DS管理离散化、创建SNES迭代、求解、输出。后续做任何问题基本都是在这个骨架上替换“弱形式定义”的部分。实际运行时只要写好网格文件用一行命令就能启动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文件生成恭喜你第一个并行有限元程序已经跑起来了。3. 核心计算流程的代码拆解3.1 网格读入DMPlex工作模式DMPlex是PETSc管理非结构化网格的主力。相比于传统的数据结构它对有限元最友好的地方在于所有实体顶点、边、面、体统一抽象成“点”并通过锥和支撑的关系描述拓扑连接。这种设计带来几个直接好处。第一网格加密和自适应时新生成的单元不需要重建整个数据结构只修改局部拓扑关系。第二并行分区时DMPlex利用Metis或Parmetis自动完成图划分把网格分布到各进程同时记录交界面的ghost单元。第三边界条件通过DMPlexGetLabel获取边界编号你不用手工遍历所有面判断边界类型。读入网格时我推荐优先使用Gmsh的msh格式PETSc对它的支持最成熟。需要注意Gmsh保存时尽量选择ASCII格式版本兼容性最好。网格文件里的物理分组Physical Tag会直接转换为DMPlex的Label后续设置Dirichlet边界或Neumann边界时直接用标签名引用即可。如果你的网格里面含有高阶单元比如二阶三角形在DMPlex里也可以通过-petscspace_degree 2直接控制基函数阶次不需要额外处理几何节点的重排。这个东西对非结构化网格来说很省心。3.2 弱形式定义与单元装配PETScFEM里最核心也最容易劝退新手的部分是PetscDS。它的作用是把你的偏微分方程弱形式告诉PETSc然后由PETSc在你提供的单元上自动计算单元刚度矩阵和残量。稳态泊松方程是理解这个过程的最佳例子。弱形式为∫ Ω ∇v · ∇u dΩ ∫ Ω v f dΩ在PetscDS里你通过PetscDSAddBoundary和残量函数定义来实现。残量函数接收每个积分点的场值、梯度以及测试函数基函数返回局部残量向量。对线性问题你也可以直接提供Jacobian组装函数对非线性问题则必须提供。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; }你不需要管单元循环、积分点坐标、映射到全局自由度这些事情。DMPlex和PetscDS自动遍历所有单元调用你定义的函数将局部矩阵累加到全局矩阵。这个过程对分布式内存是完全透明的每个进程只处理自己分到的单元交叉边界处由PETSc的Communicator机制自动做全局组装。3.3 线性求解KSP与预处理器选型关键点线性求解器是整个有限元计算最容易出问题、也最值得调整的地方。很多线性问题装配矩阵很快但KSP迭代不收敛程序一直卡在Linear solve did not converge due to DIVERGED_INDEFINITE_PC这类错误上。我的经验是遵循一个简单的选型规则二维中小规模问题自由度50万直接用直接法最省心-ksp_type preonly -pc_type lu -pc_factor_mat_solver_type mumps。MUMPS的内存和鲁棒性表现优异。大规模椭圆型问题用共轭梯度法配合超节点分解预处理器-ksp_type cg -pc_type gamg。GAMG对椭圆算子效果极其优秀接近最优。大规模非对称问题如对流占优的对流扩散首选GMRES加ILU-ksp_type gmres -pc_type ilu。ILU的填充级别可以调-pc_factor_levels 2是常见的起点。这里最关键的一条经验是永远不要用默认的KSP和PC配置去跑复杂问题默认的KSPGMRES PCJACOBI在稍大一点的问题上几乎必然失败。每次求解前先花五分钟测试不同KSP/PC组合的收敛表现这个时间远比你后面排查发散问题省得多。4. 向真实问题扩展非线性与时间相关计算4.1 非线性问题用SNES真实工程问题极少是线性的。材料非线性、大变形、接触等都会让控制方程变成F(u)0的形式。PETSc的SNES你用三个函数就能套住残量函数、Jacobian函数、单调性初值。残量函数的形式和线性问题中定义弱形式类似但f0返回的是非线性残量。Jacobian函数需要额外输出矩阵J∂F/∂u。好消息是Jacobian可以不精确提供通过SNESSetJacobian加上-snes_fd_color可以自动用有限差分近似对于单元数量少的问题非常好用。SNESSetFunction(snes, r, FormResidual, ctx); SNESSetJacobian(snes, J, J, FormJacobian, ctx);这里我对你的建议是不要一上来就追求解析Jacobian先用颜色有限差分跑通流程确认离散化没有错误后再逐步替换成解析式。解析Jacobian能显著提升收敛速度但手写推导容易出错得不偿失。SNES的收敛控制推荐显式设置-snes_type newtonls -snes_linesearch_type bt -snes_max_it 50 -snes_atol 1e-10 -snes_rtol 1e-8遇到不收敛时优先检查残量初值是否合理。很多时候不是求解器的问题而是初始猜测给得太差。4.2 瞬态问题用TS对含时间的偏微分方程PETSc的TS模块会让你的生活轻松非常多。你不需要手动写出时间步进循环只需要提供空间离散后的右端项函数和质量矩阵。对热传导这类抛物线问题典型的TS配置如下-ts_type beuler -ts_dt 0.01 -ts_max_time 1.0 -ts_monitor对应代码里只需要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最大的优势是自适应时间步。对刚性问题-ts_type rosw或-ts_type bdf配合-ts_adapt_type basic会自动调节步长在保证稳定性的前提下尽量走大步长这个特性在反应扩散、流固耦合这类多时间尺度问题中尤为重要。需要特别提醒的是TS对质量矩阵是半隐式处理还是显式处理取决于你提供的RHS Jacobian里是否包含质量项。这类细节最容易让人困惑建议先从最简单的BEuler开始跑通后再换高阶方法。4.3 并行策略与性能调优PETScFEM的并行能力是它最大的卖点之一但并行效率并非天然好需要你注意几个使用层面的事项。首先网格文件分区。如果你用Gmsh生成了一个大网格首次加载时PETSc会做一次分区平衡但这个过程默认只做一次。建议用DMPlexDistribute显式调用分区并通过-dm_plex_partition_overlap 1增加ghost层的重叠避免边界通信成为瓶颈。其次矩阵预分配。PETSc做矩阵装配时如果非零元预分配不够会触发动态增加这一步会显著拖慢性能并且带来内存碎片。通过DMCreateMatrix配合-info查看装配阶段日志确认MatAssemblyBegin没有大量reallocation告警。如果有使用MatSetOption和MatMPIAIJSetPreallocation手动预分配。第三负载均衡。非结构化网格很容易出现某些进程分到的单元计算量大、某些进程很闲的情况。Parmetis的分区效果通常比Metis好因为它会考虑进程间的通信开销。在configure时我就建议你下载parmetis原因就在这里。实测下来一个百万自由度的线弹性问题4进程并行时加速比能达到3.5左右16进程能到12以上前提是预处理器的设置跟得上。如果发现扩展性上不去优先检查是否是KSP内部全局通信太多换用-ksp_type cg和-pc_type gamg这类局部预处理组合能有效改善。5. 踩坑记录与排查技巧5.1 编译期高频错误我见过太多人在编译阶段就卡住这里把几个高频问题列出来方便你对照查找。依赖下载失败configure时--download-*联网失败多数是网络策略导致建议手动下载源码包放进petsc/目录下的packages文件夹PETSc会自动识别并跳过下载。MPI版本冲突系统MPI和PETSc自带MPICH混用会导致头文件不匹配。最稳妥的做法是全程只用PETSc自动下载的MPICH不要动系统MPI。32位索引溢出问题规模超过20亿自由度才会遇到但如果你要跑超大算例需要在configure时加上--with-64-bit-indices否则索引溢出会出现莫名其妙的内存崩溃。编译期间遇到报错第一件事不是去查错误码而是打开$PETSC_DIR/$PETSC_ARCH/lib/petsc/conf/petscvariables确认编译选项是否符合预期。很多时候是路径或环境变量没对齐。5.2 求解不收敛的排查清单求解器不收敛是有限元最消耗精力的环节我整理了一个排查顺序按这个顺序查通常能定位问题。第一检查矩阵是否对称正定。对弹性力学、热传导这类问题理论上矩阵一定对称正定。如果你装配结束后用-ksp_monitor_true_residual观察残差曲线发现残差震荡或发散先确认是不是有约束没处理干净导致矩阵奇异。第二检查边界条件是否齐全。结构问题如果没有施加足够的位移约束矩阵就会奇异CG法在第一次迭代就会因为DIVERGED_INDEFINITE_MAT退出。解决方法是先检查网格Label里边界标记是否正确传给了PetscDS。第三检查预处理器是否需要调换。ILU在矩阵条件数很差的情况下不稳定换成GAMG或直接法先验证“离散本身是否没问题”如果直接法能收敛而迭代法不行问题基本出在PC选型而不是模型定义。第四检查网格质量。畸变单元会导致单元Jacobian为负或接近零直接破坏全局矩阵性质。对这类问题先跑Gmsh的Mesh-Check确认最小单元质量及时修复网格。5.3 并行调试的几个实用技巧并行程序调试比串行麻烦一个量级有几个技巧能让你的头发少掉一些。调试时先把进程数设为1确认串行无误后再上并行。用-start_in_debugger可以在每个进程启动时挂起调试器但更实用的方式是在关键函数里加PetscPrintf打印当前进程的局部信息再通过-log_view查看各进程的耗时分布。另一个常见问题是结果在串行正确、并行却出错这几乎总是边界通信或全局索引映射的问题。在PETSc里用DMPlexGetFullData检查每个进程的本地网格和ghost数据或者用VecView输出分片向量看边界节点的值是否在相邻进程间一致。还有一个隐藏问题随机数或全局变量。有限元程序里如果用了全局累积的计数器或者非确定性随机数并行时不同进程的计算顺序不同结果就会有微小差异这类bug极难复现建议从一开始就避免在单元计算中使用全局状态。6. 个人实操中的几点额外体会最后再说几个我自己的习惯。第一个是全套遵循PETScFEM的命名和管理方式不要把Fortran和C的接口混着用。PETSc的C接口和Fortran接口虽然底层一样但错误处理方式不同混用很容易在出错时拿不到完整调用栈。第二个是善用PETSc的-info和-log_view。很多人只在出问题时才开日志平时关掉节省输出。我的建议是任何新算例第一次跑至少开一次-log_view看看矩阵装配耗时、求解器迭代次数、每一步的耗时占比这些信息能直接告诉你优化方向。第三个是保留调试版的PETSc构建。调试版速度慢但能捕获内存越界、未初始化变量这类隐蔽错误。我一般会同时保留arch-linux-c-opt和arch-linux-c-debug两套构建开发调试用debug正式跑算例用opt。还有一点是关于社区资源。PETSc官方文档和示例非常丰富$PETSC_DIR/src/snes/tutorials和$PETSC_DIR/src/ts/tutorials里面几乎有你能想到的所有经典算例从线性弹性到不可压缩流体都有。遇到问题先翻tutorials再上GitLab提issue比你在搜索引擎里瞎找效率高得多。PETScFEM这条路前期确实需要一点耐心主要是概念多、接口抽象层次高但一旦你把第一个问题完整跑通之后再换模型、换求解器、上并行就会顺畅非常多。这篇分享里的所有配置和调试经验都是我在实际项目中一条一条总结出来的希望它帮你少走几个弯路。本文还有配套的精品资源点击获取