2.2 SOR法による連立一次方程式の解法

2.2.1 連立一次方程式の解法

前節の離散式(2-1-6),(2-1-9),(2-1-10)をまとめると次式の連立一次方程式になります。

(2-2-1)
ここでi,jは行列の行番号と列番号であり前節と異なることに注意してください。 行列の大きさはN=(Nx+1)*(Ny+1)*(Nz+1)です。

連立一次方程式の解法には直接法と反復法があります。
直接法は主に密行列(行列要素がほとんど非0である行列)で用いられます。 LU分解法が代表的なものです。密行列は境界要素法などで現れます。
反復法は主に疎行列(各行の非0の行列要素が少数個である行列)で用いられます。 共役勾配法とSOR法(successive over-relaxation method, 逐次的過剰緩和法)[3][4][5]が代表的なものです。疎行列は差分法などで現れます。
本シミュレーターは差分法を用いているために反復法を採用します。 共役勾配法はたくさんの作業配列が必要となり必要メモリーが増えますが、 SOR法は作業配列が不要であるために最小限のメモリーで計算することができ、 従って大規模問題を解くことができます。 以上の理由から本シミュレーターではSOR法を採用します。

2.2.2 反復法

ここでは各種の反復法を説明します。
反復法では適当な初期値xi(0)を与えて解が収束するまで計算を繰り返します。 ここで右上の数字は反復回数を表します。

Jacobi法
Jacobi法のアルゴリズムは次式の通りです。 更新した解を別途保存しておき最後にすべての変数を更新します。 この方法はメモリーが2倍になる欠点があります。

(2-2-2)

Gauss-Seidel法
Gauss-Seidel法のアルゴリズムは次式の通りです。 変数をその都度更新します。このために余分な配列が不要になりプログラミングも簡単になります。

(2-2-3)

SOR法
SOR法のアルゴリズムは次式の通りです。いったんGauss-Seidel法と同じ式で計算し、 さらに緩和係数ωで外挿します。ωを加速パラメーターと呼ぶこともあります。 1≦ω<2の範囲で使用します。 ωが2に近いとき最も収束が速くなります。ω=1のときはGauss-Seidel法と一致します。 SOR法はGauss-Seidel法と同じ演算量で収束を速くする優れたアルゴリズムですが、 最適なωを予め知ることができないという欠点があります。

(2-2-4)

red-black SOR法
この方法は並列計算のために行列の性質を利用したものです[4]。odd-even SOR法と呼ぶこともあります。 式(2-2-4)からわかるように変数の更新を並列に実行すると変数への参照が競合し結果に再現性がありません。 ところが式(2-1-6)を見ると、自分以外はi,j,kのいずれかが1だけ異なる変数を参照しています。 これから各反復計算の第1ステップではi+j+kが偶数となる変数のみを更新し、 続いて第2ステップではi+j+kが奇数となる変数のみを更新します。 これによって変数の競合がなくなり、スレッド数によらず並列計算で同じ結果が得られます。 反復回数当たりの計算時間も変わりません。
本シミュレーターではred-black SOR法を採用します。

2.2.3 行列の特徴

本シミュレーターの行列については以下の特徴があります。

これにより、式(2-2-2),(2-2-3),(2-2-4)の右辺のbiは0であり、 左辺のiは電極上を除きます。

2.2.4 反復計算の初期値

方程式は電圧について線形なので電圧を一次変換しても行列は変わりません。 電極電圧の最大値をVmax、最小値をVminとします。 反復計算の初期値を次式で与えます。 これによって最大電圧=+1、最小電圧=-1と正規化されて無次元量になりますので、 電極の電圧によらず収束状況が一定になります。

(2-2-5)
反復計算が収束した後、電極を含む全領域の電圧を次式で戻します。
(2-2-6)

2.2.5 収束判定条件

収束判定条件は次式とします。ここでεは十分小さい数です。 左辺のr(k)を残差(residual)と呼びます。

(2-2-7)

2.2.6 反復法の収束の比較

2.2.2で述べた各種の反復法の収束状況を比較します。
問題図は図2-2-1の通りです。セル数=100*100*100です。


図2-2-1 ベンチマーク問題

図2-2-2に各手法の収束状況を示します。
これからJacobi法とGauss-seidel法は収束が遅く、SOR法により収束が大幅が速くなることがわかります。

図2-2-2 反復法の各手法の収束状況(Nx=Ny=Nz=100)

2.2.7 必要メモリー

本シミュレーターではGPUで計算を行うことも考慮して単精度で計算します。 図2-2-3に単精度と倍精度の収束状況を示します。 実用上は収束判定条件は1e-5で十分なので単精度で問題ないことがわかります。
計算に必要な配列は各節点における電圧と比誘電率と電極情報です。 このうち比誘電率と電極情報は255種類のテーブルへの間接的な参照とします。 これによりそれぞれ1バイトで済みます。 結局計算に必要なメモリーは6*Nx*Ny*Nzバイトです。 例えばNx=Ny=Nz=100のとき6MB、Nx=Ny=Nz=1000のとき6GBです。 ただし上で述べたように誘電体の比誘電率の種類と電極の電圧の種類は255通り以下です。


図2-2-3 単精度と倍精度の収束状況(red-black SOR法, ω=1.95, Nx=Ny=Nz=100)

2.2.8 緩和係数

図2.2.4に緩和係数ωを変えたときの収束状況を示します。 ωが2に近づくほど収束が速いことがわかります。


図2-2-4 緩和係数と収束状況の関係(red-black SOR法, Nx=Ny=Nz=100)

収束付近(残差=1e-5付近)では反復回数と残差の対数が直線関係にありますので、 収束の速さは指定した残差(例えば1e-5)に達する反復回数で評価することができます。

2.2.9 セル数と最適緩和係数の関係

ここでは収束が最も速くなる緩和係数を最適緩和係数と呼びωoptと表します。 上で述べたように最適緩和係数は予め知ることができません。 最適な緩和係数を知るために複数回の計算を行うことは無意味ですので、 最適緩和係数のおおよその値を求める式があれば便利です。
図2-2-5はセル数(Nx X Ny X Nz)をいろいろ変えたときの、緩和係数ωと収束するまでの反復回数の関係です。 収束反復回数が最小になる緩和係数が最適緩和係数ωoptです。 図からセル数が大きくなると最適緩和係数ωoptも大きくなることがわかります。


図2-2-5 セル数と最適緩和係数の関係(red-black SOR法, 収束判定条件=1e-5, 矢印は式(2-2-8))

図2-2-5に次式で評価される最適緩和係数ωoptを矢印で示しています。 これから次式は最適緩和係数ωoptのよい近似値であることがわかります。

(2-2-8)
実際は最適緩和係数ωoptはセルサイズのアスペクト比やセルサイズの不均一度で変わる可能性がありますので、 計算モデルを変えながら多数の計算を行うときは予め緩和係数を変えて計算し最適緩和係数ωoptを求めておくことは有効です。