先般紹介した”SMARTSOL” Matrix Solverのコーディング編の開始です。Matrix Solverとは、油層内媒体を移動する流体の物質収支を偏微分方程式(Mass Conservation Equation)を線形化し得られる疎な係数行列をもつ連立一次方程式を解くサブルーチンでリザーバーシミュレーションモデルの中核である。当然座学に付きあって頂ける程興味津々なテーマでもないので、先ずは実際にその係数行列を構築して連立一次方程式の解法手順を外観しましょう。実学主体の方が理解しやすいかも。Let’s commence the practical coding for the ”SMARSOL” Matrix Solver introduced earlier. Matrix Solver is the subroutine librares that solve the simultaneous linear equations with a sparse incidence matrix incurred after linearizing the Mass Conservation Equations by the Newton-Rapson method which define the material balance of the fluid movement in the reservoir porous media. It forms the core of the reservoir simulation model. The subject matter is not so exciting enough to study through lectures. It would be a practically amicable approach to understand as well as be much easier to learn in practice by actually building up test matrices incurred in the reservoir simulation model.

Preface

リザーバーシミュレーションモデルに於ける偏微分方程式は連続的な事象で一般解をもたない為、コンピュータで扱えるように有限差分法等による空間離散化及びAIM法による時間離散化で求めた非線形方程式をニュートン法にて線形近似した大規模な非対称疎なヤコビ行列を係数行列(Incidence Matrix)にもつ連立1次方程式を数値解法によって効率的に近似解を求めるものである。

大規模で非対称な疎行列を取り扱う場合、限られた記憶容量で如何に効率よく速く解くかが重要となる。直接法は解くべき問題のサイズが比較的小さい場合には有効な手法であるが、問題が大きくなると演算時間や記憶領域が急激に増大し実用性を失う為、一般に反復法が用いられる。

リザーバーシミュレーションモデルの分野に於ける大規模非対称疎行列の反復解法として､アイセンシュタットら(1983年Eisenstat, Elman, and Schultz)が共役勾配法の非対称問題への拡張として考案した一般共役残差法(Generalized Conjugate Residual Method、GCR法)や一般最小残差法(Generalized Minimum Residual Method、GMRES法)がよく知られている｡中でも最も簡便で有効な手法のORTHOMIN法(ORTHOgonalization MlNimization Method)は一般共役残差法の実用形で、ビンサム(1976年Vinsome)によって考案された。

なお、大規模非対称疎係数行列に対して前処理(Preconditioning)を行い反復解の高速化を図る手段も併せて取られる。疎な係数行列は固有値の分布が狭く密集し単位行列に近い程反復法の収束性が向上する事が知られており、前処理は固有値分布を改善し収束性を向上することを目的とする。リザーバーシミュレーションモデルに於ける前処理方法として、

• ILU前処理(Incomplete LU Factorization)
• NF前処理(Nested Factorisation)

が主流である｡ILU前処理法はガウス消去法を完全に行わず不完全に行い係数行列を上下の三角行列に近似分解する方法で、CMGやVIPの商用リザーバーシミュレータに用いられているごく一般的な定石技法である。一方、Nested Factorization法はILU前処理法のようにLU分解するのではなく､係数行列の元の構造を保持した状態で係数行列をブロック三重対角構造に分解する技法でECLIPSE商用シミュレータで開発された。特異な技法だが対角成分(Diagonal)に誤差行列(Error Matrix)に集約させることで卓越した収束性向上性能を示す。その収束性はロジックの性格上、格子順序付け方向(Axis Direction)に左右される特性をもつが、本法に出会った時の衝撃・感動は今でも忘れられない。

本シリーズでは､先ず､直接解法の基本解法であるガウス消去法による三角化分解の一般原理を復習した後、非零成分の構造を予め決めず任意に設定できるよう汎用性を持たせたILU前処理法を説明する。そのタ収束性の向上に資する重要な取組みについても補足する。個人的思入れが深いNF前処理法については別立てで詳説したい。また､現在のシミュレータの重要概念である時間離散化手法(AIM)による係数行列の構築・解法、加えて、強結合坑井モデル(Strongly Coupled Well Model)の処理についても事例を通して説明する予定である。最後に、反復解法の代表的な手法である一般最小残差法並びにORTHOMIN法にてシリーズを完結する計画である。

具体的なロジックは該当Paperを渉猟して頂くのが最も包括的に理解する早道であろうが、次回以降俺なりのポイントを講義することとし、本稿では手っ取り早く係数行列の構築に関する実学を紹介するに留める。

Test Incidence Matrix

三次元領域を想定した最小限の有限差分格子システムを想定したテスト係数行列(Test Incidence Matrices)を作成する。本稿のポイントは格子の順番付けを行う際の方向性によって係数行列の違いを認識することを学ぶ。

Gridding

具体的な例題として､対象格子領域が４＊３＊２=２４格子システムで２坑井(生産井・圧入井各々)を有する図１のような三次元問題を考える。

Null格子を10格子設定する。Active格子は全部で14格子となる｡内、AIM(Adaptive Implicit)格子は3格子、残り全ての格子はFIM(Fully Implicit)格子とする。又、坑井のLayer数は坑井１が３枚､坑井２が２枚の多層仕上げ坑井とする｡Internal SystemとはNull格子を除いたActive格子のみを対象とする格子システムを示す。

格子システムは６通りの組み合わせがあるが、格子システムによる収束性の差異を比較検証するために、小規模な格子システム故その差異は明白ではないものの、収束性が最も優れた格子システムから最悪システムであろう３種類の格子システムを採用する。具体的には、３軸(\(x,y,z)\)の格子方向に応じた(２-３-４)、(３-２-４)、及び(４-３-２)の３ケースを検討する。因みに、”SMARTSOL”は最適な格子方向に補正する機能を有する。

Incidence Matrix

三次元システムを対象とするので､各格子間の流れは７点有限差分展開(7-Point Operator)を施すことで図２に示す係数行列を得る｡例えば､(４-３-２)システムに於いては、格子1は上流格子2、5、及び13に対して流れが生じる｡このように格子全体について全ての流れを調べて係数行列を構築する。

ILU前処理

ILU前処理に於ける係数行列を下図２−１に示す。４-３-２格子システムは対角成分に非零要素が集中しており他のシステムに比べて収束性向は秀でていると思料。

何れ説明するが、ガウス消去法の掃き出しの過程で坑井格子は重要な要素である故、係数行列に明示的に表示している。

NF前処理

NF前処理に於ける係数行列を下図２−２に示す。図にはブロック三重対角行列を形成するブロック分割を明示的に表示している。２-３-４若しくは３-２-４格子システムは他のシステムに比べブロック分割が多く演算に於いて収束性向は優位と思料。

リザーバーシミュレータに於いて、一般的に垂直方向(Z軸)のTransmissibility(流体の易動度)はHarmonic Averageになる為、平面(X、Y軸)方向に比べて低下する傾向がある。本稿ではその影響について考慮していないが(”SMARTSOL”は格子方向最適化補正機能を有する)、この垂直軸のTransmissibilityがNF前処理の収束性に大きな影響を与えるファクターである。

係数行列分解上坑井格子の存在は無関係であるが、ILU前処理法に比べてこの点は、前処理演算速度の観点からは重要な要素である。

今回は格子の順番付けによる係数行列の仕上がり具合の違いの紹介で実学を終わりとします。次回からしばらく講義主体の取組みになる見込み。

以上
つづく