FiPy

FVM を利用した PDE ソルバーです。

FiPy 3.4
https://www.ctcms.nist.gov/fipy/

良いですね。項に対応する関数（fipy.terms package）を追加するだけで PDE を解くことができます。昔、先輩に紹介いただいたソフトに近いですね。MATLAB も含め有償のソフトはこのように簡単なのでしょうか？というか、これが10年以上前に既に組まれていたのですから驚きです。

https://www.ctcms.nist.gov/fipy/documentation/FAQ.htmlより引用

例えば、拡散のみであれば以下のように書きます。
eq = DiffusionTerm(D)
非定常だと
eq = TransientTerm() == DiffusionTerm(D)
移流項を加えると
eq = TransientTerm() + AdvectionTerm(U) == DiffusionTerm(D)
風上差分と同様の解方は、convection にありました。
eq = TransientTerm() + UpwindConvectionTerm(U) == DiffusionTerm(D)
advection, convection は使い分けられていますね。違いはコチラ↓
https://physics.stackexchange.com/questions/24489/what-exactly-is-the-difference-between-advection-and-convection

これでソース項を除く式の表現ができました。地下水の流速を求めるなら別途式が必要ですが、1文で移流拡散の離散化を任せられるのは手軽ですね。
さらに、∇は以下で表現されていますので、任意の項を作成・追加可能です。
().divergence

コチラの方は日本語で解説されています。
http://penguinitis.g1.xrea.com/study/FiPy/FiPy.html
短いですね。これで２次元や３次元の近似解を求めることができるのですから、手軽としか言いようがありません。


では、試行。
公式のexample と 上記の方のコードを動かしながら、３次元まで拡張です。ここまで短く記述できるのは素晴らしい。
https://www.ctcms.nist.gov/fipy/examples/README.html

from fipy import *

nx = 100
ny = 100
nz = 10
Lx = 1.
Ly = 1.
Lz = 1.
dx = Lx/nx
dy = Ly/ny
dz = Lz/nz

mesh = Grid3D(dx, dy, dz, nx, ny, nz)
phi = CellVariable(name='phi', mesh=mesh, value=.5)

D = 1.
U = (3., 2.,1.)
phi1 = 1.
phi2 = 0.5

phi.constrain(phi1, mesh.facesLeft)
phi.constrain(phi2, mesh.facesRight)

eq = TransientTerm() + UpwindConvectionTerm(coeff=U) == DiffusionTerm(coeff=D)

T=1.
dt=0.1
for step in range(1, int(T/dt)+1):
    print('time step:', dt*step)
    eq.solve(var=phi, dt=dt)
    TSVViewer(vars=phi).plot(filename=str(step)+".tsv")
    VTKViewer(vars=phi).plot(filename=str(step)+".vtk")

結果、2次元まではOKですが、3次元では異常に時間がかかりました（実務では厳しいくらいに）。ラップトップで12分。この程度の計算で12分もかかるのはダメ。ソルバーの変更か並列化が必須でしょう。

Mayavi が動かなかったのと plotly での表示が遅かったので、可視化は別ソフトに頼ることを念頭に VTK 書き出し（plotly での表示は以下）。

import pandas as pd
import plotly.express as px
import plotly.graph_objects as go

df = pd.DataFrame({'x' : mesh.x, 'y': mesh.y, 'z':mesh.z, 'value': phi})

fig = go.Figure(data=go.Volume(
    x=mesh.x,
    y=mesh.y,
    z=mesh.z,
    value=df.value,
    opacity=0.5,
    surface_count=20
    ))

fig.update_layout(margin=dict(l=0, r=0, b=0, t=0))
fig.show()

phi が FiPy のオブジェクトになっていますので、やや扱いづらい。ま、df にはなるので書き換えが可能ならPHREEQCとの連成も容易でしょう。