簡介這份資源面向地球物理、水文地質(zhì)及環(huán)境工程領域的學習者與研究人員提供基于開源 Python 庫 SimPEG 子模塊 SimPEG.FLOW 編寫的地下水流方程示例工程用于理解達西定律下的正向模擬與參數(shù)反演流程。壓縮包共 20 個文件約 45KB以 9 個 py 腳本為核心配合 rst 文檔、makefile 構(gòu)建腳本、yml 持續(xù)集成配置及 license、readme 等說明文件目錄中可見 Richards 測試與 examples 示例模塊結(jié)構(gòu)緊湊、便于按模塊閱讀源碼。已有 492 人學習下載說明其在同類小眾方向中具備一定參考價值。讀者可從中獲取網(wǎng)格構(gòu)建、邊界條件設置、流體力學模型組織方式以及反演優(yōu)化調(diào)用的具體寫法并借助文檔與示例快速搭建自己的地下水流動模擬腳本適合作為入門 SimPEG.FLOW 的實踐起點。1. 從一份“跑不起來”的地下水模型說起前陣子幫一個做水文地質(zhì)的團隊看代碼他們手里有個用 Python 寫的地下水流數(shù)值模型注釋里寫著“基于 SimPEG 框架”結(jié)果換臺機器就報依賴沖突網(wǎng)格一加密內(nèi)存直接爆掉邊界條件改了之后結(jié)果對不上。這類項目最尷尬的地方在于它不像 Web 項目那樣跑起來就能看頁面數(shù)值模擬的“跑通”和“算對”是兩回事。simpegflow 這個資源解決的正是這個痛點——它把地下水流方程穩(wěn)定流與非穩(wěn)定流用 SimPEG 的有限體積離散體系完整實現(xiàn)了一遍從網(wǎng)格生成、參數(shù)賦值、邊界條件組裝到線性求解和結(jié)果可視化形成一條可復現(xiàn)的鏈路。適合誰適合已經(jīng)懂一點地下水動力學、想用 Python 做二次開發(fā)的研究生和工程技術人員也適合想拿一個真實 PDE 案例練手 SimPEG 的開發(fā)者。下面我按“它是什么 → 怎么跑 → 怎么改”的順序把這份源碼拆開講。2. 項目定位與技術棧為什么選 SimPEG 而不是自己寫求解器2.1 這個項目到底在算什么地下水流方程的核心是達西定律加質(zhì)量守恒。穩(wěn)定流情況下寫成散度形式就是 ?·(K?h) q其中 K 是滲透系數(shù)場h 是水頭q 是源匯項。非穩(wěn)定流再多一個儲水項 Ss·?h/?t。聽起來簡單但真正落地時有三個麻煩非均質(zhì) K 場怎么在網(wǎng)格上表達、定水頭邊界和通量邊界怎么離散、大規(guī)模稀疏線性系統(tǒng)怎么高效求解。simpegflow 的價值在于它沒有從零造輪子而是把 SimPEG 已經(jīng)做好的有限體積離散、網(wǎng)格對象和線性求解器接口復用起來只專注寫地下水這一層的物理方程和邊界條件。從代碼組織看它通常會把物理過程封裝成一個類似GroundwaterFlow的類繼承 SimPEG 的BasePDESimulation或等價的 PDE 基類然后實現(xiàn)make_boundary_conditions、getA、getrhs這些鉤子方法。這種寫法的好處是網(wǎng)格、插值、矩陣組裝這些臟活由框架兜底你改物理模型時只需要動方程本身。我一般會先看這個類的繼承鏈因為它決定了你能白拿哪些能力。2.2 技術棧與依賴關系組件作用備注Python 3.8運行環(huán)境建議用 conda 隔離環(huán)境SimPEG有限體積離散、網(wǎng)格、求解器核心依賴版本敏感numpy / scipy數(shù)組運算與稀疏矩陣scipy.sparse 組裝全局矩陣discretize網(wǎng)格與差分算子SimPEG 的底層依賴matplotlib / pyvista結(jié)果可視化二維用 matplotlib三維用 pyvistapytest單元測試驗證解析解對比這里要強調(diào)一點SimPEG 的版本迭代比較快discretize的 API 在不同版本間有過調(diào)整所以拿到源碼后第一件事是看requirements.txt或environment.yml里鎖定的版本別直接pip install simpeg裝最新版否則很可能在網(wǎng)格對象的屬性名上翻車。血淚經(jīng)驗是數(shù)值框架的“最新版”往往意味著接口變動而不是 bug 修復。2.3 目錄結(jié)構(gòu)長什么樣一個典型的 simpegflow 工程目錄大致是這樣組織的simpegflow/ ├── simpegflow/ │ ├── __init__.py │ ├── flow.py # 地下水流物理類核心 │ ├── boundary.py # 邊界條件封裝 │ ├── utils.py # 網(wǎng)格生成、參數(shù)場構(gòu)造 │ └── io.py # 結(jié)果導出 ├── examples/ │ ├── steady_2d.py # 二維穩(wěn)定流算例 │ ├── transient_2d.py # 二維非穩(wěn)定流算例 │ └── layered_3d.py # 三維層狀介質(zhì)算例 ├── tests/ │ └── test_analytical.py # 與解析解對比 ├── requirements.txt └── README.mdflow.py是你要重點讀的文件examples/是你跑通項目的入口tests/是你驗證“算對沒算對”的依據(jù)。很多人上來就改flow.py結(jié)果連基準算例都沒跑過改完也不知道對錯這是大忌。3. 跑通第一個算例從環(huán)境搭建到出圖3.1 環(huán)境準備與依賴安裝先建一個干凈環(huán)境別在系統(tǒng) Python 里折騰。常見做法是用 conda# 創(chuàng)建獨立環(huán)境Python 版本按 requirements 來 conda create -n simpegflow python3.9 -y conda activate simpegflow # 先裝 SimPEG 及其底層依賴再裝項目本身 pip install simpeg discretize pip install -r requirements.txt # 以可編輯模式安裝方便后續(xù)改源碼 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 二維穩(wěn)定流算例逐段拆解打開examples/steady_2d.py核心流程一般分四步建網(wǎng)格、定參數(shù)場、設邊界、求解。下面這段是我按項目常見寫法整理的骨架關鍵行都加了注釋import numpy as np import discretize from simpegflow import GroundwaterFlow # 1. 建二維張量網(wǎng)格x 方向 100 格z 方向 50 格 # extent 控制物理尺寸單位統(tǒng)一用米 mesh discretize.TensorMesh( [np.ones(100) * 1.0, np.ones(50) * 1.0], origin[0.0, 0.0] ) # 2. 構(gòu)造非均質(zhì)滲透系數(shù)場 K這里用兩層介質(zhì)做演示 K np.ones(mesh.n_cells) * 1e-5 # 背景值單位 m/s K[mesh.cell_centers[:, 1] 25.0] 1e-7 # 下半層低滲 # 3. 實例化物理類把網(wǎng)格和參數(shù)場傳進去 flow GroundwaterFlow(mesh, KK) # 4. 設置邊界左右定水頭上下隔水 flow.set_boundary_conditions( head_left10.0, # 左側(cè)水頭 10 m head_right5.0, # 右側(cè)水頭 5 m no_flow_topTrue, no_flow_bottomTrue ) # 5. 組裝并求解線性系統(tǒng) A h rhs h flow.solve() # 6. 可視化水頭場 flow.plot(headh, savefigsteady_2d.png)逐段說明網(wǎng)格用TensorMesh是因為矩形區(qū)域用張量網(wǎng)格最省事np.ones(100)*1.0表示每個格子 1 米寬總長 100 米。K 場按單元中心坐標分層賦值這是模擬層狀介質(zhì)的標準做法注意單位是 m/s量級搞錯結(jié)果會差好幾個數(shù)量級。set_boundary_conditions里定水頭邊界會進入右端項隔水邊界對應零通量在有限體積里就是該面通量置零。solve()內(nèi)部會調(diào)用getA()組裝稀疏矩陣再交給 scipy 的求解器。最后plot只是薄封裝真正出圖靠 matplotlib 的pcolormesh。3.3 非穩(wěn)定流算例的時間步設置非穩(wěn)定流算例transient_2d.py多了一個時間循環(huán)。常見寫法是顯式或隱式時間離散項目里一般用后向歐拉保證穩(wěn)定# 時間參數(shù)總時長 100 天步長 1 天 t_total 100.0 * 86400.0 # 換算成秒 dt 1.0 * 86400.0 n_steps int(t_total / dt) # 儲水系數(shù)場單位 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時間就得用秒所以天數(shù)要乘 86400。后向歐拉的矩陣每個時間步都要重新組裝如果 Ss 和 dt 不變其實可以預分解一次矩陣反復用這是后面二次開發(fā)可以優(yōu)化的點。跑完這個算例你會看到水頭從初始的均勻 10 米逐漸向左右邊界驅(qū)動的穩(wěn)態(tài)過渡這個過程圖最能說明模型是否合理。4. 核心模塊拆解矩陣組裝與邊界條件實現(xiàn)4.1 有限體積離散的矩陣從哪來地下水流方程的有限體積離散本質(zhì)是對每個單元做質(zhì)量守恒流入減流出等于源匯。SimPEG 的discretize提供了面通量算子項目里getA方法通常這樣寫def getA(self): # K 在面上取調(diào)和平均保證層狀介質(zhì)界面通量連續(xù) K_face self.mesh.aveF2CC_1D(self.K) # 簡化示意實際按面方向取 # 面通量系數(shù)矩陣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() # 轉(zhuǎn) CSR 格式求解更快邏輯說明get_face_inner_product把滲透系數(shù)和幾何因子面積除以格心距揉成一個對角矩陣作用在梯度上得到面通量face_divergence再把面通量累加成單元的凈流出。兩者相乘就是離散后的 ?·(K?)。參數(shù)上K 的賦值位置要和網(wǎng)格面匹配如果 K 定義在單元中心框架會自動做平均但層狀介質(zhì)界面處用算術平均會高估通量合格做法是調(diào)和平均這也是很多人結(jié)果偏大的原因。tocsr()不能省COO 格式在迭代求解時慢得離譜。4.2 邊界條件的三種類型怎么落地邊界條件是這個項目里最容易改錯的地方。定水頭邊界Dirichlet在離散時是把對應單元的水頭固定矩陣行替換為單位行、右端項設為給定水頭通量邊界Neumann是修改右端項把給定通量加到對應面的散度里隔水邊界就是通量為零什么都不加。項目里boundary.py一般會提供三個函數(shù)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里逐行清零在單元數(shù)大時很慢常見優(yōu)化是用 scipy 的lil_matrix先改再轉(zhuǎn) CSR或者用布爾索引批量操作。符號約定也要盯緊達西定律里通量方向和水頭梯度方向相反邊界通量流入計算域時右端項取正搞反了結(jié)果會整體偏移。我見過有人把定水頭邊界加在錯誤的單元索引上模型照樣能解但物理上完全不對所以加完邊界一定要打印邊界單元的水頭核對。4.3 線性求解器的選擇與預條件A組裝出來后是個大型稀疏非對稱矩陣非穩(wěn)定流時更明顯直接scipy.sparse.linalg.spsolve在小規(guī)模下沒問題網(wǎng)格一過十萬單元就吃力。項目里通常會暴露求解器選項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})參數(shù)說明rtol是相對殘差容差地下水問題一般 1e-8 夠用maxiter給太小會誤報不收斂。info非零就是黑匣子報警別忽略先檢查矩陣是否奇異——最常見原因是整個區(qū)域沒有定水頭邊界純隔水邊界下方程只有差一個常數(shù)的解矩陣必然奇異。這個坑我在早期項目里踩過模型跑完水頭全是同一個值排查半天才發(fā)現(xiàn)是邊界漏設。5. 避坑與排查數(shù)值模型最容易翻車的五個地方5.1 現(xiàn)象結(jié)果量級離譜水頭差了幾百倍原因通常是單位不統(tǒng)一。K 用 m/s、長度用 km、時間用天混著用必然出錯。解決在項目入口統(tǒng)一單位制我一般強制全部轉(zhuǎn) SI長度米、時間秒、K 米每秒在utils.py里寫一個assert_units函數(shù)做量綱檢查。改完先跑一個一維解析解對比誤差在 1% 以內(nèi)再往下做。5.2 現(xiàn)象網(wǎng)格加密后內(nèi)存爆掉原因是全局矩陣按單元數(shù)平方增長三維網(wǎng)格幾十萬單元時 CSR 也扛不住。解決優(yōu)先用迭代求解器加預條件別用直接法網(wǎng)格只在關心區(qū)域加密外圍用粗網(wǎng)格必要時用float32存 K 場省內(nèi)存。如果還不行考慮區(qū)域分解但那是另一個量級的改造。5.3 現(xiàn)象非穩(wěn)定流結(jié)果隨時間震蕩原因是時間步太大或用了顯式格式。后向歐拉無條件穩(wěn)定但如果dt相對儲水項太小矩陣條件數(shù)會變差。解決先估算特征時間Ss·L2/Kdt取它的十分之一量級同時檢查 Ss 是否全為零Ss 為零時非穩(wěn)定項消失方程退化成穩(wěn)定流時間循環(huán)等于白跑。5.4 現(xiàn)象邊界附近水頭出現(xiàn)非物理尖峰原因是定水頭邊界加在了單元中心而不是面上或者邊界單元同時被加了兩種邊界。解決明確每個邊界只施加一種條件加完后打印邊界單元的A行和rhs值核對用mesh.cell_centers定位邊界單元時注意坐標是格心不是格點。5.5 現(xiàn)象換機器后 import 報錯或結(jié)果不一致原因是 SimPEG 版本不同導致網(wǎng)格屬性名或默認插值方式變化。解決鎖死requirements.txt版本用 conda 導出environment.yml在tests/里保留解析解對比測試換環(huán)境后先跑測試再跑算例測試不過就別信結(jié)果。6. 二次開發(fā)把模型改成你自己的水文地質(zhì)問題拿到這份源碼后真正有價值的是把它改成你關心的場景。改造路徑我一般按三步走先換參數(shù)場再換邊界最后換方程。換參數(shù)場最直接K和Ss都是按單元賦值的數(shù)組你可以從鉆孔數(shù)據(jù)插值生成也可以讀柵格文件。常見做法是寫一個build_K_from_layers(mesh, layers)函數(shù)把地層界面坐標轉(zhuǎn)成單元索引。換邊界時如果項目只提供了定水頭和隔水你需要補一個通量邊界或降水入滲邊界入滲本質(zhì)是面通量加到頂部單元的右端項即可注意單位是 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(非線性迭代未收斂)驗證方法上別只看圖好不好看。我的習慣是每個改造都配一個解析解或已發(fā)表基準算例做對比比如一維穩(wěn)定流的線性水頭分布、二維承壓含水層的 Theis 解退化情形。tests/test_analytical.py里通常已經(jīng)有模板照著加一個test_your_case就行。改造完先跑測試再跑自己的算例兩步都過才算數(shù)。最后說個我自己的教訓早期我改邊界條件時圖省事直接在solve()返回后手動改幾個單元的水頭看起來結(jié)果對了但質(zhì)量守恒徹底破壞后續(xù)做粒子追蹤時通量場完全對不上。從那以后我每次改邊界都強制走一遍“組裝后打印 A 行 核對 rhs 跑守恒檢查”的流程再也不敢在結(jié)果上打補丁。希望這份拆解能幫你少走幾步彎路把 simpegflow 真正跑成自己手里的工具。本文還有配套的精品資源點擊獲取