前回は有限体積法の線形方程式を解く共役勾配法という解法について説明しました。
今回はシリーズ最終回です。共役勾配法をPythonでプログラミングして有限体積法のプログラムを完成させたいと思います。
目次
共役勾配法のプログラム
前回説明したとおり共役勾配法のアルゴリズムは以下のとおりです。
- 適当な初期値 ${\bf T}_0$ を決める。
- ${\bf r}_0 = {\bf b} - {\bf A} {\bf T}_0$、 ${\bf d}_0 = {\bf r}_0$ とする。
- 収束するまで以下の計算を繰り返す。
① $\alpha_k=({\bf d}_k, {\bf r}_k)/({\bf d}_k, {\bf A}{\bf d}_k) $
② ${\bf T}_{k+1} = {\bf T}_k + \alpha_k {\bf d}_k$
③ ${\bf r}_{k+1} = {\bf r}_k - \alpha_k{\bf A} {\bf d}_k$
④ 残差 ${\bf r}_{k+1}$ で収束判定。収束なら終了。
⑤ $\beta_k= ({\bf r}_{k+1}, {\bf r}_{k+1})/({\bf r}_k, {\bf r}_k)$
⑥ ${\bf d}_{k+1} = {\bf r}_{k+1} + \beta_k {\bf d}_k $
今回の有限体積法でのPythonプログラムは以下のようになります。
d = np.zeros_like(temp)
r = np.zeros_like(temp)
ad = np.zeros_like(temp)まず配列を用意します。d は方向ベクトル ${\bf d}$、r は残差ベクトル ${\bf r}$、ad は行列 ${\bf A}$ と ${\bf d}$ の積 ${\bf A}{\bf d}$ です。
#-- initial residual
for icell in range(1, ncell+1):
r[icell] = bc[icell] - ac[icell] * temp[icell]
for iface in range(4):
inbor = cnbor[icell][iface]
r[icell] -= af[icell][iface] * temp[inbor] ここは、残差ベクトルの初期値 ${\bf r}_0 = {\bf b} - {\bf A} {\bf T}_0$ を計算しています。
d = np.copy(r)
rr = np.dot(r, r)
bnrm = np.linalg.norm(bc)d = np.copy(r) は ${\bf d}_0 = {\bf r}_0$ で、r を d にコピーしています。
rr = np.dot(r, r) は、r と r の内積を計算しています。あとで $\beta$ を計算するときに使います。
bnrm = np.linalg.norm(bc) は、bc のノルム $||{\bf b}||$ を計算しています。相対残差を計算するときに使います。
#-- iteration
for iter in range(1, itermax+1):共役勾配法の反復計算のforループです。iter は反復回数で 1 から itermax まで回ります。
#-- A*d
for icell in range(1, ncell+1):
ad[icell] = ac[icell] * d[icell]
for iface in range(4):
inbor = cnbor[icell][iface]
ad[icell] += af[icell][iface] * d[inbor] ${\bf A}{\bf d}$ を計算しています。
#-- alpha
alpha = np.dot(d, r) / np.dot(d, ad)$\alpha$ を計算します。
#-- temperature
temp += alpha * d温度 ${\bf T}$ を更新します。
#-- residual
r -= alpha * ad
resd = np.linalg.norm(r) / bnrm残差ベクトル ${\bf r}$ を更新して、相対残差 $||{\bf r}||/||{\bf b}||$ を計算します。
#-- convergence check
if resd <= eps:
break 収束判定をチェックします。残差 resd が収束判定値 eps 以下になると、break で反復のforループを抜けます。
eps は最初に定義している小さな数です。今回は 1.0e-5 としています。残差がどのくらいになれば収束とみなすかは、実際の問題では厳密な指標を示すのは難しいです。収束判定値 eps が大きいと収束が早いですが、定常状態になっていない場合があります。逆に eps が小さいと真の解により近づきますが、収束までの反復回数が多くなり計算時間がかかります。計算結果やあとで紹介するモニター点の履歴などを見ながら収束を総合的に判断する必要があります。
#-- beta
beta = np.dot(r, r) / rr
rr = np.dot(r, r)
#-- d
d = r + beta * d最後に $\beta$ を計算し、${\bf d}$ を更新します。rr は次のステップ用に r と r の内積を格納しています。
共役勾配法のアルゴリズムはシンプルなので、プログラムもアルゴリズムの通りでとてもシンプルに書けます。
２次元熱伝導問題の有限体積法プログラム
以上でプログラムの骨格は出来たので、今回の２次元熱伝導問題の有限体積法プログラムの完成版を載せておきます。
以下はプログラムの補足です。
モニター点の表示
今回の問題では板の中央の温度を知りたいとありました。そこで、指定した位置の値を画面に表示することにします。有限体積法はセルごとに値を持っているため、求めたい点が含まれるセルの値を表示することにしましょう。
#-- monitor
pmon = np.array([0.5, 0.5])
#-- find monitor cell
imon = 0
flag = False
for i in range(1, nx+1):
for j in range(1, ny+1):
icell = i + (j - 1) * nx
x = (i - 1) * dx
y = (j - 1) * dy
if pmon[0] >= x and pmon[0] <= (x+dx) and pmon[1] >= y and pmon[1] <= (y+dy):
imon = icell
flag = True
break
if flag:
break まず pmon でモニターしたい位置の座標を指定します。そして、その点がどのセルに含まれるか見つけていきます。
forループで i、j でループします。x、y はそのセルの左下角の座標です。ここでは単純にモニター点のx座標、y座標とセルの点を比較して、セルに含まれるかどうかを判定しています。imon がモニター点のセル番号になります。セルが見つかれば break で j 側のループを抜けます。ただし、forループは２重になっているので、flag を用意し i 側のループも抜けるようにしています。
結果のグラフと可視化
結果表示として、残差とモニター点の温度の履歴、温度のコンター図を表示してみましょう。これまでの講座でグラフや図の表示はよく出てきました。その応用です。
import matplotlib.pyplot as plt可視化用に matplotlib をインポートします。
iter_plt = []
resd_plt = []
temp_plt = []まず、共役勾配法の計算の前にグラフ表示のデータ格納用に配列を用意しておきます。iter_plt は反復回数、resd_plt は残差、temp_plt はモニター点の温度を格納する配列です。
#-- plot data
iter_plt.append(iter)
resd_plt.append(resd)
temp_plt.append(temp[imon])
print('iter {:6d} , resd {:12.4e} , temp {:12.4e}'.format(iter, resd, temp[imon]))共役勾配法の反復ループの中でこれらの配列に、それぞれの値を格納しています。同時に、print() で画面にも出力します。
# post processing ====================================================
#-- residual and monitor graph
fig1 = plt.figure()
ax1 = fig1.add_subplot(2, 1, 1)
ax1.plot(iter_plt, resd_plt)
ax1.set(ylabel='Residual')
ax1.set_yscale('log')
ax2 = fig1.add_subplot(2, 1, 2)
ax2.plot(iter_plt, temp_plt)
ax2.set(xlabel='Iteration', ylabel='Monitor') 計算が終わったあとに、グラフを出力します。ここでは、一つの figure に残差とモニター温度の二つのグラフ（axes）を表示します。残差は通常、対数表示で見るので set_yscale('log') で対数グラフにしています。
#-- contour plot
fig2 = plt.figure()
ax3 = fig2.add_subplot(1, 1, 1)
gx = np.arange(0.0, lx+0.5*dx, dx)
gy = np.arange(0.0, ly+0.5*dy, dy)
X, Y = np.meshgrid(gx, gy)
val = np.zeros((ny, nx))
for i in range(1, nx+1):
for j in range(1, ny+1):
icell = i + (j - 1) * nx
val[j-1][i-1] = temp[icell]
im = ax3.pcolormesh(X, Y, val, cmap='jet')
ax3.set_xlabel('x')
ax3.set_ylabel('y')
ax3.set_aspect('equal')
fig2.colorbar(im)
plt.show()次の figure には温度コンター図を表示します。meshgrid で表示用のグリッドを作成しています。配列 val に表示用の温度データを格納します。表示用の配列は、i と j が入れ替わった順番になっていることに注意してください。pcolormesh() でコンター図を作成しています。
計算実行
これで、２次元熱伝導を解くための有限体積法のプログラムが完成しました。Pythonで実行してみてください。
iter 1 , resd 1.4015e-01 , temp 3.0000e+02
iter 2 , resd 1.0976e-01 , temp 3.0000e+02
iter 3 , resd 7.8111e-02 , temp 3.0000e+02
iter 4 , resd 6.4608e-02 , temp 3.0000e+02
iter 5 , resd 5.3684e-02 , temp 3.0000e+02
...
iter 252 , resd 1.1651e-05 , temp 4.3572e+02
iter 253 , resd 1.0807e-05 , temp 4.3572e+02
iter 254 , resd 1.0389e-05 , temp 4.3572e+02
iter 255 , resd 1.0008e-05 , temp 4.3572e+02
iter 256 , resd 9.0604e-06 , temp 4.3572e+02最初は大きかった残差が徐々に小さくなっていく様子がわかります。モニターしている温度は初期値の300Kから変化していき、435.72Kに収束していきました。共役勾配法は256回の反復（イタレーション）で収束しました。
残差とモニター温度のグラフを見ると収束状況がよくわかります。
温度は150イタレーション以降はほぼ一定となり収束しているようです。このように、残差だけでなく目的の量の推移も合わせて見ていくとよいでしょう。
結果の温度コンター図は以下のようになりました。
右下の500Kのところから左上の300Kの箇所へ徐々に温度が変化している様子がわかります。今回の条件では $y \ge 0.8$ で熱伝導率が小さいため、温度の低い箇所はその領域に抑えられています。
このプログラムを少し変更すれば、境界条件を変えたり、違う熱伝導率の組み合わせで計算することができます。また、発熱量 qc を与えると発熱を伴う計算も行うことができます。以下の結果は中央付近に２箇所発熱を与えた場合の例です。
他の条件でもいろいろ計算してみてください。
まとめ
今回のシリーズでは、２次元熱伝導方程式をシミュレーションするプログラムを作りながら、有限体積法について学びました。有限体積法がどのような手法なのか理解できたでしょうか。
流体解析に適用したり、多様な形状のセルを使ったり、非定常問題を解いたりするためには、さらに複雑な定式化やスキーム、解法などが必要になります。また、計算速度を高速化したり、大規模な問題に適用したりするためには、データの格納方法や参照方法、他のアルゴリズムやプログラム言語を考える必要も出てきます。ただ、有限体積法をざっくり把握するには、ひとつの取っ掛かりになったのではないかと思います。
汎用の解析ソフトを使うにしても、その中でどのような手法で解かれているかを知っておくのは役に立ちます。「残差」とか「メッシュ」とか「解法」とか、言葉をブラックボックスにしておくより、きちんとイメージできていた方が解析やシミュレーションの仕事の助けになるかもしれません。
流体解析もまた機会があれば書いてみたいと思います。