“SMART” is an inhouse three-dimensional, three-phase black oil reservoir simulator in the petroleum upstream industry, in which its function and capability are at minimum compatible with the industry standard simulator “ECLIPSE”. “iSMART” is the initially primitive version of “SMART”, which incorporates the essential functions with even limited, especially it features the adaptive implicit capability and the strongly coupled well model. For purpose of studying a general approaches of fundamental formulations of reservoir equations together with the strongly coupled well modeling, and as an ultimate purpose, how to develop a reservoir simulator, therefore, ”iSMART” would be the best candidate as the practical material. Here, I wish to address an overview of ”SMART/iSMART” prior to sharing with the ”SMARTSOL” large sparse matrix solver in series. “SMART”は独自開発した石油上流業界に於ける三次元三相ブラックオイル油層シミュレーターである。その機能と用途は業界標準シミュレーター”ECLIPSE”との互換性を最低限有する。”iSMART”は”SMART”の暫定初期バージョンであるが、最低限なるも根幹となる機能、特に当時”State of Art”であった順応陰解法と油層との強結合坑井モデルを組込んでいる。依って、油層流動方程式の基本的な定式化の一般的なアプローチ、油層との強結合坑井モデルの解法、及び最終的な目的として油層シミュレーターを開発する手法を学ぶ上で、”iSMART”は最も適した教材と云えよう。”SMARTSOL”大規模疎行列解法をシリーズで紹介するに当って、先ず”SMART/iSMART”の概要を紹介したい。

大学の専門は理工系であったがどちらかと云うと文系の教材に興味を抱き全くコンピュータと無縁な学校生活を送ったが、会社に入ってからは独自の世界を構築できるコンピュータの世界に興味を持った。要は一種の会社実世界からの逃避だったかも知れない。俺の新入社員の頃はパワハラ極地の時代、酷かったからなあ。

俺がパソコンに関わり始めたのはパソコンの黎明期に当たる1980年代初頭、丁度CP/MとMSDOSの覇権争いが始まった頃でもある。当時日本ではPC-9801が最も人気があったが、性格が捩れた俺は富士通のFM-11を買って遊んでいた。ハードよりもソフト開発に興味があり、8087数値演算コプロセッサを別途購入して数値計算のソフトの高速化を図っていたものだ。8087は東京エレクトロン本社へ出向いて購入。チョンガーだった俺は当時のお金で11万円払ったと記憶している。婚約指輪と同程度の値段だったような。

結果的にCP/MをパクったMSDOSがOSの覇権を一気に獲得できたのはMicrosoft社が商業的に上手だったからであろうが、この頃同社の株を買っていれば大金持ちになれたかも？

石油上流業界に於いてもコンピュータ応用分野で著しい進展があった。最も衝撃的だったのはラプラス変換を用いた坑井圧力テスト解析法の出現である。こうした進展は油層流動モデル(以降、Reservoir Simulation、シミュレータ)の開発もパソコンにも適用された。一方、国内に於いても海外の潮流に遅れまいと石油公団でも海外の大学・ソフト開発業界と連携した強化プロジェクトを推進中であった。このような前向きの環境の中で、俺は業務の合間をぬって見様見真似でソフト開発に追従していたものである。

Reservoir Simulation

1980年代初頭から始まったシミュレータの開発は一気に進み、1980年代終わりには主要なロジックは完成している。それ以降はモデルのResolutionの精微化は進んだが新たな技術進展は乏しい。今後の課題はAIとの融合した最適な自動モデルの構築並びに自動ヒストリーマッチング機能の分野ではなかろうか。

一般的にシミュレータの開発に於いて、最も苦労するのは大規模疎行列解法(以降、Matrix Solver、マトリックスソルバー)であると言われるが、俺的には坑井モデルの構築・解法が最もTime Consumingで頭を使ったテーマである。

大規模疎行列の解法は既に確立されており数学的な算法に沿って単純にコーディングするだけである。只、同解法に対して順応陰解法(以下、AIM、Adaptive Implicit)を導入する際には大規模疎行列の解法手順のあらゆる段階で変更が必要となる。同様に、坑井モデルは油層との強結合(Strongly Coupled Well Model)の概念は確立されているもののその解法はやや煩雑で少なからぬ処理が必要となる。一方、坑井グループ集団制御モデル(Collective Control Model)は定番のロジックはなく(参考資料あり)、プログラマーの器量が最も発揮できる分野であろう。

“iSMART”の坑井モデルは坑井単体に限定しているが、”SMART”ではグループ設定を可能ならしめる坑井グループ集団制御モデルを組込んだ包括的な坑井モデルへと深化している。一般的なシミュレータは坑底圧力(Bottom Hole Flowing Pressure)のみを独立変数とした坑井モデルであるが、”SMART”はWellbore Friction Modelの導入を睨んでマルチセグメント坑井(Multisegment Well)を視野に多独立変数に対応可能な坑井モデルを組込んでいる。更に、順応陰解法にも対応し、油層グリッド毎の圧力・飽和率との同時解法を達成している。因みに、Matrix Solver ”SMARTSOL”は”SMART”、”iSMART”で同一である。

“iSMART”と”SMART”の機能の比較表を下表に示す。

Dual Porosity、Dual Permeability、Psuedo Miscible、API Tracking、Local Grid Refinement、Wellbore Friction Model等のオプション機能は特異な環境に対するシミュレーションとなる為、”SMART”では扱っていない。尚、”成分型数値モデル”を必要とする条件設定もありうるが、実業では成分型シミュレーションモデルの必要性を感じなかったので、PVT Behaviour Packageの開発のみに留まり成分型モデルの開発には到らなかった。

Reservoir Simulationに関する概要についてはPetroWikiで包括的に纏められている。ご参考迄。

尚、些細なことであるが、”iSMART”は当時の標準であったFORTRAN77に沿ってCodingしている。FORTRAN90に準拠した”SMART”に比べてCodingの柔軟性には劣る。

Performance

過去シミュレータの開発時代、石油上流技術協会(SPE, Society of Petroleum Engineers)でシミュレータの性能評価を競うプロジェクト(CSP, Comparative Solution Projects)が企画されている。これは当時の最新シミュレータの開発動向を探る上で不可欠なテーマを取扱っている。ブラックオイルモデルの基本的な性能比較を睨んだCSPプロジェクトは１,２,７,及び９で特にCSP９は非均質な油層を対象とした算法難易度の高いケースである。CSPプロジェクトの概要はペンシルバニア州立情報科学技術大学が公開している情報で整理されている。

ここではCSP９を事例として”SMART/iSMART”の性能を紹介する。
対象ケースはグリッド数は24x25x15、総数9000個、坑井数は26坑(内、生産井25坑、圧入井1坑)である。坑井配置は下図参照の事。

Comparison by Commercial Models

本事例に於ける”SMART”の性能を商業シミュレータと比較した結果を下表に示す。併せて、”iSMART”の結果も示す。Numerical Control、Time Step Sizeは”ECLIPSE/IMEX”を含めてDefault設定のままでCPU TIMEの短縮、収束性を向上させる試みは一切行っていない。
Gistは以下。

1. “SMART”
   商業シミュレータに比し遜色のない性能を示す。
   前処理”NF”ベストケースは商業シミュレータ”ECLIPSE”に比べて格段に優位な性能を示す。収束性及び演算時間共略互角である。下表1-1対1-4参照。更に、Newton Methodの収束判定条件を変更した場合、演算速度は格段に遅くなるが収束性は保持し互角である。下表対1-2対1-4参照。但、Time Step2箇所で反復計算が発生している(赤太字箇所)。一方、前処理”ILU”ベストケースは、”IMEX”に対して秀でた収束性能を示すが演算速度は劣る。下表2-1対2-4-1参照。
   • Matrix Preconditioning (前処理手法)
     俺の実態感としては、”ECLIPSE”で提唱された”NF”(Nested Factorization)の方が総合的にRobustな前処理法であると判断している。”IMEX”等に導入されている一般的な前処理手法である”ILU”(Incomplete LU Factorization)の方が優位、とBehieらは指摘しているが、原理的には”NF”前処理の方が収束性が優る。特に、大規模且つ独立変数が多い坑井モデルに対してはその傾向が顕著に発揮すると思料する。
     前処理は後述する線形化によって得られる疎係数行列の固有値分布を改善することで収束性の向上を図る技法である。”NF”前処理法は、疎係数行列を上下三角行列に近似分解する”ILU”前処理法と異なり、疎係数行列の元構造を保持しつつ対角成分(Diagonal)に誤差行列(Error Matrix)を集約するアルゴリズムで、高位な固有値分布を期待できる。加えて、些細な事であるが坑井の開閉の都度要求される係数行列(Incidence Matrix)の”Fill-in”の位置を算定する”Symbolic Factorization”の必要性も不要である。後述するように対象課題、坑井独立変数の差異等モデルによっては前処理手法の優位性が如実に反映される結果を発揮する。
   • Newton Method
     Newton Methodは油層内を移動する流体の物質収支を規定する偏微分方程式の求解の為にて線形化する技法である。”E”、”I”、或いは”C”は収束判定基準の指標で、”E”は”ECLIPSE”、”I”は”IMEX”で採用されている収束判定基準を示す。これらの理論・実証については各社User’s Manualを参照して頂くのが宜しいかと思うが、”E”判定は”Material Balance Error”と”Residual Error”双方を判定対象とする。一方、”I”判定は”Material Balance Error (L1 Norm)”のみを判定対象とする緩和条件設定(Relaxation)である。特に、”ILU”前処理ケースに於いて”E”設定を課する場合、収束性・演算速度共IMPES並に低下する非実用的な性能となる。加えて、Time Stepで多数反復計算が発生している。下表2-2及びNewton Methodに於けるErrorとSpeedに関する比較表を参照の事。これは前述したPreconditioning法の理論的背景に密接な関係があり、ILU前処理では”E”条件下では非実用的な設定となってしまう。又、”C”判定はErrorとSpeedを両立させるCompromiseな判定基準を意味し”iSMART”項で触れる。実際として偏微分方程式の非線形性(Linearity)は然程大きくなく、”ECLIPSE/IMEX”同様、判定基準値(Tolerance)を1.0^-3に設定している。
   • 独立変数
     ”SMART”及び”ECLIPSE”に於ける坑井独立変数は３個であるが、”IMEX”等一般シミュレータは１変数である。詳細については次項”iSMART”で後述する。
   • Automatic Time Step Control
     ”ECLIPSE”の判定基準は前回のTime Stepに於ける飽和率の変化量から算出したTime Truncation Errorに基づく。一方、”IMEX”は具体的な判定ロジックは明らかでないが、当該Time Stepに於ける圧力・飽和率双方の変化量から判定していると推察。”SMART”は”ECLIPSE”と同じ手法を採用している。
   • Linear Solver Convergence Tolerance
     DefaultのToleranceは”ECLIPSE/IMEX””同様1.0^-4である。”SMARTSOL”では1.0^-6と厳しい設定をDefaultとしている。同Toleranceを変更するとNewton反復回数は変化するがTime Stepに対する改善は見られないようだ。
   • Matrix Iterative Method(反復法解法)
     “SMARTSOL”に組込んであるIterative Methodは”ORTHOMIN”と”GMRES”であるが、双方とも性能は互角。
   • ”AIM”
     ”AIM”に特化した商業シミュレータも存在するが、”FIM”(以下、陰解法、Fully Implicit)に対する優位性は確認できない。
   • CPU Time
     相対的に商業シミュレータは迅速である。行列計算の最適化、配列の次元配置・処理方法等にに関与しているとも想定されるが仔細は不明。尚、”SMARTSOL”は、”CRS(Compressed Row Storage)”形式をベースに見やすいように行列要素を独自の手法に則り格納しているが、処理スピードの観点から改善の余地有りと認識している。
2. “iSMART”
   負荷の軽いモデルなこともあって、総じて市販シミュレータと同等のパフォーマンスを示す。しかし、Time Step反復計算が数回発生し、演算速度も劣っている。Primitiveなモデルである故”SMART”との差を論じることは合理的ではなく、特記事項のみ下記する。
   • Newton Method
     独立変数(圧力と飽和率)の変化量を判定対象に加えて、先に述べた”E”判定の収束性の高速化を図る”C”判定は、特にILU前処理法に於いて”E”判定基準下で収束性が緩慢となる過程で有効に働く。
   • Time Stepの反復計算は解くべき独立変数の算定値が現実的な範囲から逸脱する為に発生する反復計算である。一義的な原因は、”iSMART”ではニュートン法の微分項に数値微分(前進差分近似)値を用いている事による。”SMART”では全て解析微分値を求めて収束性の向上を図っているように、解析微分に変更することで改善が期待できよう。
   • 坑井独立変数
     ”iSMART”は”IMEX”同様、坑底圧力のみを解くべき独立変数とした坑井モデルを組込んである。一方、”SMART”は”ECLIPSE”と同じく独立変数は３変数とし、発展性のある尤度の高い設定である。理論的に独立変数が少ないほど処理時間が迅速となり、加えて収束性に対する耐性は高いと考えられる。確かに、坑底圧力のみを独立変数とする坑井モデルの場合でも実用的な精度でCross Flowモデルを構築できる為、初期的なモデルで多用されてきた経緯がある。

参考迄に、代表的な入力ファイルを添付する。

唯ー保有するTabletのLenovo TB-8740FによるConsole出力のDumpファイルも添付する。

3D View

”IMEX”の三次元可視化ソフトを用いた結果を表示する。下図では最終Time Stepに於ける油層圧力、油、ガス、水飽和率を示している。これらは以前”IMEX”の試用時に作成したものである。

Wayforward

以上、独自開発した”iSMART/SMART”シミュレータの概要・性能はお分かり頂けたと思う。今後、Large Sparse Matrix Solver “SMARTSOL”のソースコードを参照しながら具体的な大規模疎行列の数値解法手順について共有・説明していく予定。

途中で息切れするかも知れないが、Life Workとして取り組む所存。
乞うご期待！

Appendix

1. Intel oneAPI Fortran Compiling Option

専らifortを利用している。

Release Mode
    FFLAGS = -fast
Debug Mode
    FFLAGS = -O0 -CB -check uninit -init=snan -fpe0 -ftrapuv -traceback -g

gfortranでは以下の不具合が発生する。理由は不明

• -Ofast最適化ではmachine epsilon演算が上手く計算されない
• print表記がScreen上で反映されない

Release Mode
   FFLAGS = -ffree-line-length-none -fbackslash -O

2. References

1. Akand,W.I.,Kamy S.:”A Review on SPE’s Comparative Solution Projects (CSPs),” Journal of Petroleum Science Research (JPSR) Volume 2 Issue 4, October 2013.
   1. Odeh,A.S.:”Comparison of Solutions to a Three-Dimensional Black-Oil Reservoir Simulation Problem,” JPT 33 (1): 13–25. SPE9723, 1981.
   2. Weinstein,H.G., Chappelear,J.E., and Nolen,J.S.:”Second Comparative Solution Project: A Three-Phase Coning Study,” JPT 38 (3): 345-353, SPE10489, 1986.
   3. Kenyon,D.:”Third SPE Comparative Solution Project: Gas Cycling of Retrograde Condensate Reservoirs,” JPT 39 (8): 981-997, SPE12278, 1987.
   4. Aziz,K., Ramesh,A.B., and Woo,P.T.:”Fourth SPE Comparative Solution Project: Comparison of Steam Injection Simulators,” JPT 39 (12): 1576–1584, SPE13510, 1987.
   5. Killough,J.E. and Kossack,C.A.:”Fifth Comparative Solution Project: Evaluation of Miscible Flood Simulators,” Presented at the SPE Symposium on Reservoir Simulation, San Antonio, SPE16000, 1987.
   6. Firoozabadi,A. and Thomas,L.K.:”Sixth SPE Comparative Solution Project: Dual-Porosity Simulators,” JPT 42 (6): 710-715, 762-763, SPE18741, 1990.
   7. Nghiem,L., Collins,D.A., and Sharma,R.:”Seventh SPE Comparative Solution Project: Modelling of Horizontal Wells in Reservoir Simulation,” Presented at the SPE Symposium on Reservoir Simulation, Anaheim, SPE21221, 1991.
   8. Quandalle,P.:”Eighth SPE Comparative Solution Project: Gridding Techniques in Reservoir Simulation,” Presented at the SPE Symposium on Reservoir Simulation, New Orleans, SPE25263,1993.
   9. Killough,J.E.:”Ninth SPE Comparative Solution Project: A Reexamination of Black-Oil Simulation,” Presented at the SPE Reservoir Simulation Symposium, San Antonio, SPE29110, 1995.
   10. Christie,M.A. and Blunt,M.J.:”Tenth SPE Comparative Solution Project: A Comparison of Upscaling Techniques,” SPEREE 4 (4): 308–317, SPE72469, 2001.
2. Reservoir Simulation Basics
   1. Van Everdingen,A.F. and Hurst,W.:“The Application of the Laplace Transformation to Flow Problems in Reservoirs,” Transaction of AIME, 1949.
   2. Thomas,L.K., Lumpkin,W.B. and Reheis,C.M.:“Reservoir Simulation of Variable Bubble-Point Problems,” SPEJ, Feburay 1976.
   3. Aziz,K., and Settari,A.:”Petroleum Reservoir Simulation,” Applied Science Publishers Ltd., London, 1979.
   4. Abou-Kassem,J.H. and Aziz,K.:“Analytical Well Models for Reservoir Simulation,” SPEJ, August 1985.
   5. Forsyth,P.A. and Sammon,P.H.:“Gas Phase Appearance and Disapperance in Fully Implicit Black Oil Simulation,” SPEJ, October 1984.
   6. Nghiem,L.:“An Integral Approach to Discretizing the Reservoir Fiow Equations,” Transactions of AIME, 285, 1988.
3. FIM
   1. Bansal,P.P., Harper,J.L., McDonald,A.E., Moreland,E.E., and Odeh,A.S.:”A Strongly Coupled, Fully Implicit,Three Dimensional,Three Phase Reservoir Simulator,” SPE8329, 1979.
   2. Au,A.D.K., Behie,A., Rubin,B., and Vinsome,K.:“Techniques for Fully Implicit Reservoir Simulation,” SPE9302, 1980.
   3. Behie,A. and Vinsome,P.K.W.:“Block Iterative Methods for Fully Implicit Reservoir Simulation,” SPEJ, October 1982.
   4. Lim.K.T., Schiozer,D.J., and Aziz,K.:“A New Approach for Residual and Jacobian Array Construction in Reservoir Simulators,” SPE Computer Applications, August 1995.
4. Well Model
   1. Trimble,R.H. and McDonald,A.B.:“A Strongly Coupled, Fully Implicit Three-Dimensional, Three-Phase Well Coning Model,” SPEJ, August 1981.
   2. Behie,A., Collins,D., Forsyth,P.A.Jr. and Sammon,P.H.:“Fully Coupled Multiblock Wells in Oil Simulation,” SPEJ August 1985.
   3. Holmes,J.A.:“Enhancement to the Strongly Coupled, Fully Implicit Well Model: Wellbore Crossflow Modeling and Collective Well Control,” SPE12259, 1983.
   4. Modline,A.D., Coats,K.H., and Wells,M.W.:”A Super Position Method for Representing Wellbore Crossflow in Reservoir Simulation” SPERE, August 1992.
   5. Coats,K.H., Thomas,L.K., and Pierson,R.G.:“Compositional and Black Oil Reservoir Simulation,” SPE29111, 1995.
   6. Holmes,J.A., Barkve,T., and Lund,O.:”Application of a multisegment well model to simulate flow in advanced well,” SPE50646, 1998.
5. AIM
   1. Thomas,G.W. and Thurnau, D.H.: “Reservoir Simulation Using an Adaptive Implicit Method,” SPE10120, 1981.
   2. Thomas,G.W. and Thurnau, D.H.: “The Mathematical Basis of the Adaptive Implicit Method,” SPE10495, 1982.
   3. Forsyth,P.A. and Sammon, P.H.: “Practical Considerations for Adaptive Implicit Methods in Reservoir Simulations,” J. Comp. Phys., 1986.
   4. Fung,L.S., Collins,D.A. and Nghiem,L.X.:“An Adaptive-Implicit Switching Criterion Based on Numerical Stability Analysis,” SPE16003, 1989.
6. Matrix Solver
   1. T.A.Davis,:”Direct methods for sparse linear systems,” SIAM, Philadelphia, PA, 2006.
   2. Y.Saad,:”Iterative Methods for Sparse Linear Systems,” SIAM, Philadelphia, PA, 2003.
   3. Price,H.S. and Coats,K.H.:“Direct Methods in Reservoir Simulation,” SPEJ, June 1974.
   4. Appleyard,J.R.:”Nested Factorization,” Presented at the SPE Reservoir Simulation Symposium, San Francisco. SPE12264, 1983.
   5. Behie,A.:“Comparison of Nested Factorization, Constrained Pressure Residual, and Incomplete Factorization Preconditionings,” SPE13331, 1985.
   6. Vinsome,P.K.W.:“ORTHOMIN, An Iterative Method for Solving Sparse Sets of Simultaneous Linear Equations,” SPE5729, 1976.
   7. Saad,Y. and Schultz, M.H.:”GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems,” SIAMJ. Sci. Stat. Comput., 7, 1986.

以上