简介这份资源面向地球物理、水文地质及环境工程领域的学习者与研究人员提供基于开源 Python 库 SimPEG 子模块 SimPEG.FLOW 编写的地下水流方程示例工程用于理解达西定律下的正向模拟与参数反演流程。压缩包共 20 个文件约 45KB以 9 个 py 脚本为核心配合 rst 文档、makefile 构建脚本、yml 持续集成配置及 license、readme 等说明文件目录中可见 Richards 测试与 examples 示例模块结构紧凑、便于按模块阅读源码。已有 492 人学习下载说明其在同类小众方向中具备一定参考价值。读者可从中获取网格构建、边界条件设置、流体力学模型组织方式以及反演优化调用的具体写法并借助文档与示例快速搭建自己的地下水流动模拟脚本适合作为入门 SimPEG.FLOW 的实践起点。1. 从一份“跑不起来”的地下水模型说起前阵子帮一个做水文地质的团队看代码他们手里有个用 Python 写的地下水流数值模型注释里写着“基于 SimPEG 框架”结果换台机器就报依赖冲突网格一加密内存直接爆掉边界条件改了之后结果对不上。这类项目最尴尬的地方在于它不像 Web 项目那样跑起来就能看页面数值模拟的“跑通”和“算对”是两回事。simpegflow 这个资源解决的正是这个痛点——它把地下水流方程稳定流与非稳定流用 SimPEG 的有限体积离散体系完整实现了一遍从网格生成、参数赋值、边界条件组装到线性求解和结果可视化形成一条可复现的链路。适合谁适合已经懂一点地下水动力学、想用 Python 做二次开发的研究生和工程技术人员也适合想拿一个真实 PDE 案例练手 SimPEG 的开发者。下面我按“它是什么 → 怎么跑 → 怎么改”的顺序把这份源码拆开讲。2. 项目定位与技术栈为什么选 SimPEG 而不是自己写求解器2.1 这个项目到底在算什么地下水流方程的核心是达西定律加质量守恒。稳定流情况下写成散度形式就是 ∇·(K∇h) q其中 K 是渗透系数场h 是水头q 是源汇项。非稳定流再多一个储水项 Ss·∂h/∂t。听起来简单但真正落地时有三个麻烦非均质 K 场怎么在网格上表达、定水头边界和通量边界怎么离散、大规模稀疏线性系统怎么高效求解。simpegflow 的价值在于它没有从零造轮子而是把 SimPEG 已经做好的有限体积离散、网格对象和线性求解器接口复用起来只专注写地下水这一层的物理方程和边界条件。从代码组织看它通常会把物理过程封装成一个类似GroundwaterFlow的类继承 SimPEG 的BasePDESimulation或等价的 PDE 基类然后实现make_boundary_conditions、getA、getrhs这些钩子方法。这种写法的好处是网格、插值、矩阵组装这些脏活由框架兜底你改物理模型时只需要动方程本身。我一般会先看这个类的继承链因为它决定了你能白拿哪些能力。2.2 技术栈与依赖关系组件作用备注Python 3.8运行环境建议用 conda 隔离环境SimPEG有限体积离散、网格、求解器核心依赖版本敏感numpy / scipy数组运算与稀疏矩阵scipy.sparse 组装全局矩阵discretize网格与差分算子SimPEG 的底层依赖matplotlib / pyvista结果可视化二维用 matplotlib三维用 pyvistapytest单元测试验证解析解对比这里要强调一点SimPEG 的版本迭代比较快discretize的 API 在不同版本间有过调整所以拿到源码后第一件事是看requirements.txt或environment.yml里锁定的版本别直接pip install simpeg装最新版否则很可能在网格对象的属性名上翻车。血泪经验是数值框架的“最新版”往往意味着接口变动而不是 bug 修复。2.3 目录结构长什么样一个典型的 simpegflow 工程目录大致是这样组织的simpegflow/ ├── simpegflow/ │ ├── __init__.py │ ├── flow.py # 地下水流物理类核心 │ ├── boundary.py # 边界条件封装 │ ├── utils.py # 网格生成、参数场构造 │ └── io.py # 结果导出 ├── examples/ │ ├── steady_2d.py # 二维稳定流算例 │ ├── transient_2d.py # 二维非稳定流算例 │ └── layered_3d.py # 三维层状介质算例 ├── tests/ │ └── test_analytical.py # 与解析解对比 ├── requirements.txt └── README.mdflow.py是你要重点读的文件examples/是你跑通项目的入口tests/是你验证“算对没算对”的依据。很多人上来就改flow.py结果连基准算例都没跑过改完也不知道对错这是大忌。3. 跑通第一个算例从环境搭建到出图3.1 环境准备与依赖安装先建一个干净环境别在系统 Python 里折腾。常见做法是用 conda# 创建独立环境Python 版本按 requirements 来 conda create -n simpegflow python3.9 -y conda activate simpegflow # 先装 SimPEG 及其底层依赖再装项目本身 pip install simpeg discretize pip install -r requirements.txt # 以可编辑模式安装方便后续改源码 pip install -e .pip install -e .这一步很关键它让你改simpegflow/下的源码后不用重新安装就能生效。如果requirements.txt里对 SimPEG 有版本上限务必遵守比如simpeg0.18,0.20这种写法就是在告诉你别越界。装完后跑一句python -c import simpegflow; print(simpegflow.__file__)确认导入路径指向你 clone 下来的目录而不是 site-packages 里的旧副本。3.2 二维稳定流算例逐段拆解打开examples/steady_2d.py核心流程一般分四步建网格、定参数场、设边界、求解。下面这段是我按项目常见写法整理的骨架关键行都加了注释import numpy as np import discretize from simpegflow import GroundwaterFlow # 1. 建二维张量网格x 方向 100 格z 方向 50 格 # extent 控制物理尺寸单位统一用米 mesh discretize.TensorMesh( [np.ones(100) * 1.0, np.ones(50) * 1.0], origin[0.0, 0.0] ) # 2. 构造非均质渗透系数场 K这里用两层介质做演示 K np.ones(mesh.n_cells) * 1e-5 # 背景值单位 m/s K[mesh.cell_centers[:, 1] 25.0] 1e-7 # 下半层低渗 # 3. 实例化物理类把网格和参数场传进去 flow GroundwaterFlow(mesh, KK) # 4. 设置边界左右定水头上下隔水 flow.set_boundary_conditions( head_left10.0, # 左侧水头 10 m head_right5.0, # 右侧水头 5 m no_flow_topTrue, no_flow_bottomTrue ) # 5. 组装并求解线性系统 A h rhs h flow.solve() # 6. 可视化水头场 flow.plot(headh, savefigsteady_2d.png)逐段说明网格用TensorMesh是因为矩形区域用张量网格最省事np.ones(100)*1.0表示每个格子 1 米宽总长 100 米。K 场按单元中心坐标分层赋值这是模拟层状介质的标准做法注意单位是 m/s量级搞错结果会差好几个数量级。set_boundary_conditions里定水头边界会进入右端项隔水边界对应零通量在有限体积里就是该面通量置零。solve()内部会调用getA()组装稀疏矩阵再交给 scipy 的求解器。最后plot只是薄封装真正出图靠 matplotlib 的pcolormesh。3.3 非稳定流算例的时间步设置非稳定流算例transient_2d.py多了一个时间循环。常见写法是显式或隐式时间离散项目里一般用后向欧拉保证稳定# 时间参数总时长 100 天步长 1 天 t_total 100.0 * 86400.0 # 换算成秒 dt 1.0 * 86400.0 n_steps int(t_total / dt) # 储水系数场单位 1/m Ss np.ones(mesh.n_cells) * 1e-4 flow GroundwaterFlow(mesh, KK, SsSs) flow.set_boundary_conditions(head_left10.0, head_right5.0) h flow.initial_head(10.0) # 初始水头全场 10 m for step in range(n_steps): # 每个时间步组装 (A Ss*V/dt) h^{n1} rhs Ss*V/dt * h^n h flow.solve_transient(h, dtdt)这里dt的单位必须和 K、Ss 的单位自洽K 用 m/s、Ss 用 1/m时间就得用秒所以天数要乘 86400。后向欧拉的矩阵每个时间步都要重新组装如果 Ss 和 dt 不变其实可以预分解一次矩阵反复用这是后面二次开发可以优化的点。跑完这个算例你会看到水头从初始的均匀 10 米逐渐向左右边界驱动的稳态过渡这个过程图最能说明模型是否合理。4. 核心模块拆解矩阵组装与边界条件实现4.1 有限体积离散的矩阵从哪来地下水流方程的有限体积离散本质是对每个单元做质量守恒流入减流出等于源汇。SimPEG 的discretize提供了面通量算子项目里getA方法通常这样写def getA(self): # K 在面上取调和平均保证层状介质界面通量连续 K_face self.mesh.aveF2CC_1D(self.K) # 简化示意实际按面方向取 # 面通量系数矩阵diag(K_face * 面积 / 距离) Mf self.mesh.get_face_inner_product(self.K) # 散度算子把面通量映射到单元净流出 D self.mesh.face_divergence # 组装 A D * Mf * GG 是梯度算子 A D Mf self.mesh.cell_gradient return A.tocsr() # 转 CSR 格式求解更快逻辑说明get_face_inner_product把渗透系数和几何因子面积除以格心距揉成一个对角矩阵作用在梯度上得到面通量face_divergence再把面通量累加成单元的净流出。两者相乘就是离散后的 ∇·(K∇)。参数上K 的赋值位置要和网格面匹配如果 K 定义在单元中心框架会自动做平均但层状介质界面处用算术平均会高估通量合格做法是调和平均这也是很多人结果偏大的原因。tocsr()不能省COO 格式在迭代求解时慢得离谱。4.2 边界条件的三种类型怎么落地边界条件是这个项目里最容易改错的地方。定水头边界Dirichlet在离散时是把对应单元的水头固定矩阵行替换为单位行、右端项设为给定水头通量边界Neumann是修改右端项把给定通量加到对应面的散度里隔水边界就是通量为零什么都不加。项目里boundary.py一般会提供三个函数def apply_dirichlet(A, rhs, mesh, cells, head_value): # 把指定单元的行清零对角线置 1右端项设为水头值 for c in cells: A[c, :] 0.0 A[c, c] 1.0 rhs[c] head_value return A, rhs def apply_neumann(rhs, mesh, faces, flux_value): # 通量边界直接加到右端项注意符号流入为正 rhs mesh.face_divergence.T flux_value return rhs注意apply_dirichlet里逐行清零在单元数大时很慢常见优化是用 scipy 的lil_matrix先改再转 CSR或者用布尔索引批量操作。符号约定也要盯紧达西定律里通量方向和水头梯度方向相反边界通量流入计算域时右端项取正搞反了结果会整体偏移。我见过有人把定水头边界加在错误的单元索引上模型照样能解但物理上完全不对所以加完边界一定要打印边界单元的水头核对。4.3 线性求解器的选择与预条件A组装出来后是个大型稀疏非对称矩阵非稳定流时更明显直接scipy.sparse.linalg.spsolve在小规模下没问题网格一过十万单元就吃力。项目里通常会暴露求解器选项from scipy.sparse.linalg import cg, spilu, LinearOperator # 对称正定时用共轭梯度配不完全 LU 预条件 ilu spilu(A.tocsc()) M LinearOperator(A.shape, ilu.solve) h, info cg(A, rhs, MM, rtol1e-8, maxiter1000) if info ! 0: raise RuntimeError(f求解未收敛info{info})参数说明rtol是相对残差容差地下水问题一般 1e-8 够用maxiter给太小会误报不收敛。info非零就是黑匣子报警别忽略先检查矩阵是否奇异——最常见原因是整个区域没有定水头边界纯隔水边界下方程只有差一个常数的解矩阵必然奇异。这个坑我在早期项目里踩过模型跑完水头全是同一个值排查半天才发现是边界漏设。5. 避坑与排查数值模型最容易翻车的五个地方5.1 现象结果量级离谱水头差了几百倍原因通常是单位不统一。K 用 m/s、长度用 km、时间用天混着用必然出错。解决在项目入口统一单位制我一般强制全部转 SI长度米、时间秒、K 米每秒在utils.py里写一个assert_units函数做量纲检查。改完先跑一个一维解析解对比误差在 1% 以内再往下做。5.2 现象网格加密后内存爆掉原因是全局矩阵按单元数平方增长三维网格几十万单元时 CSR 也扛不住。解决优先用迭代求解器加预条件别用直接法网格只在关心区域加密外围用粗网格必要时用float32存 K 场省内存。如果还不行考虑区域分解但那是另一个量级的改造。5.3 现象非稳定流结果随时间震荡原因是时间步太大或用了显式格式。后向欧拉无条件稳定但如果dt相对储水项太小矩阵条件数会变差。解决先估算特征时间Ss·L²/Kdt取它的十分之一量级同时检查 Ss 是否全为零Ss 为零时非稳定项消失方程退化成稳定流时间循环等于白跑。5.4 现象边界附近水头出现非物理尖峰原因是定水头边界加在了单元中心而不是面上或者边界单元同时被加了两种边界。解决明确每个边界只施加一种条件加完后打印边界单元的A行和rhs值核对用mesh.cell_centers定位边界单元时注意坐标是格心不是格点。5.5 现象换机器后 import 报错或结果不一致原因是 SimPEG 版本不同导致网格属性名或默认插值方式变化。解决锁死requirements.txt版本用 conda 导出environment.yml在tests/里保留解析解对比测试换环境后先跑测试再跑算例测试不过就别信结果。6. 二次开发把模型改成你自己的水文地质问题拿到这份源码后真正有价值的是把它改成你关心的场景。改造路径我一般按三步走先换参数场再换边界最后换方程。换参数场最直接K和Ss都是按单元赋值的数组你可以从钻孔数据插值生成也可以读栅格文件。常见做法是写一个build_K_from_layers(mesh, layers)函数把地层界面坐标转成单元索引。换边界时如果项目只提供了定水头和隔水你需要补一个通量边界或降水入渗边界入渗本质是面通量加到顶部单元的右端项即可注意单位是 m/s 乘面积。换方程是进阶操作。比如加一个非线性项潜水面问题里 K 依赖饱和厚度就需要在getA里让 K 随水头更新外层套 Newton 迭代。下面是一个改造骨架def solve_nonlinear(self, h0, tol1e-6, max_iter20): h h0.copy() for it in range(max_iter): # 用当前水头更新饱和厚度重算 K_eff K_eff self.K * self.saturated_thickness(h) A self.assemble_A(K_eff) rhs self.assemble_rhs(h) h_new self.linear_solve(A, rhs) # 水头变化小于容差即收敛 if np.linalg.norm(h_new - h) / np.linalg.norm(h_new) tol: return h_new h h_new raise RuntimeError(非线性迭代未收敛)验证方法上别只看图好不好看。我的习惯是每个改造都配一个解析解或已发表基准算例做对比比如一维稳定流的线性水头分布、二维承压含水层的 Theis 解退化情形。tests/test_analytical.py里通常已经有模板照着加一个test_your_case就行。改造完先跑测试再跑自己的算例两步都过才算数。最后说个我自己的教训早期我改边界条件时图省事直接在solve()返回后手动改几个单元的水头看起来结果对了但质量守恒彻底破坏后续做粒子追踪时通量场完全对不上。从那以后我每次改边界都强制走一遍“组装后打印 A 行 核对 rhs 跑守恒检查”的流程再也不敢在结果上打补丁。希望这份拆解能帮你少走几步弯路把 simpegflow 真正跑成自己手里的工具。本文还有配套的精品资源点击获取