基于SimPEG的地下水数值模拟实战：从跑通到二次开发

简介：这份资源面向地球物理、水文地质及环境工程领域的学习者与研究人员，提供基于开源 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 技术栈与依赖关系

这里要强调一点：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.md

flow.py是你要重点读的文件，examples/是你跑通项目的入口，tests/是你验证“算对没算对”的依据。很多人上来就改flow.py，结果连基准算例都没跑过，改完也不知道对错，这是大忌。

3. 跑通第一个算例：从环境搭建到出图

3.1 环境准备与依赖安装

先建一个干净环境，别在系统 Python 里折腾。常见做法是用 conda：

# 创建独立环境，Python 版本按 requirements 来 conda create -n simpegflow python=3.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 有版本上限，务必遵守，比如simpeg>=0.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, K=K) # 4. 设置边界：左右定水头，上下隔水 flow.set_boundary_conditions( head_left=10.0, # 左侧水头 10 m head_right=5.0, # 右侧水头 5 m no_flow_top=True, no_flow_bottom=True ) # 5. 组装并求解线性系统 A h = rhs h = flow.solve() # 6. 可视化水头场 flow.plot(head=h, savefig="steady_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, K=K, Ss=Ss) flow.set_boundary_conditions(head_left=10.0, head_right=5.0) h = flow.initial_head(10.0) # 初始水头全场 10 m for step in range(n_steps): # 每个时间步组装 (A + Ss*V/dt) h^{n+1} = rhs + Ss*V/dt * h^n h = flow.solve_transient(h, dt=dt)

这里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 * G，G 是梯度算子 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, M=M, rtol=1e-8, maxiter=1000) 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²/K，dt取它的十分之一量级；同时检查 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, tol=1e-6, max_iter=20): 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 真正跑成自己手里的工具。

本文还有配套的精品资源，点击获取