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