前処理法(Preconditioner)としてIncomplete LU FactorizationとNested Factorizationを紹介してきた。自分の現役時代のシミュレーション環境ではこれにて充分であったが、最近より複雑な地層を対象とする大規模な油層シミュレーションの環境下ではより効果的な解法としてCPR前処理(Constrained Pressure Residual Preconditioner)なる手法が用いられているようだ。まだin-house貯留層シミュレーターSMARTのMatrix Solution Package (SMARTSOL)には実装していないが、その概要及び実装の必要性等につき考察する。

Preface

上流石油業界に置ける貯留層数値モデルに於いて近年の大規模非対称行列の連立一次方程式の解法としては元の方程式をヤコビ係数行列(ヤコビアン、ヤコビ行列)と残差ベクトルから成るクリロフ部分空間に射影して解く反復法が多く利用されている。その中でORTHOMIN法とGMRES法が貯留層流体の解法として有効であることが示されている。

これら反復法の収束性は係数行列の固有値分布に依存する。固有値分布が少なく且つ単位行列に近い程収束が早い。元の係数行列に良く似た前処理行列\(M\)を適用することで固有値分布を改善することができる。この前処理行列\(M\)による固有値分布の改善プロセスを前処理法と呼ぶ。

CPR前処理法に於いて、先に述べた圧力と飽和率離散化モデルの特性が線形解法の収束性を左右することから、ヤコビ係数行列\(J\)を圧力ブロックと飽和率ブロックに分離して各々を反復法にて近似解を求めるプロセスを考える。

先ず、ヤコビ係数行列を圧力ブロック行列\(J_{pp}\)、飽和率ブロック行列\(J_{ss}\)、及び圧力･飽和率結合行列\(J_{ps}\)と\(J_{sp}\)に分割する。同様に、解ベクトルと残差ベクトルも圧力と飽和率に夫々分割する。圧力及び飽和率変数を各々\(x_{p}\)と\(x_{s}\)、圧力と飽和率に関連する残差ベクトルを夫々\(r_{p}\)と\(r_{s}\)と定める。詰まり、ヤコビ係数行列を次のように分割する。

$$J\Delta x=\begin{bmatrix} J_{pp} & J_{ps} \\ J_{sp} & J_{ss} \end{bmatrix}
\begin{bmatrix} \Delta x_{p} \\ \Delta x_{s} \end{bmatrix} = r= \begin{bmatrix} r_{p} \\ r_{s} \end{bmatrix} \tag{1}$$

ヤコビ係数行列は総格子数\(Nc\)に対して各格子毎に\(Ne\)組の方程式で構成される。

今、\((Nc⋅Ne) \times Nc\)からなる行列\(C\)とその転置行列\(C^{T}\)を導入し、圧力ブロックと同じ次元の\(Nc\)から成る単位行列\(I_{p}\)にて、次式のように定義する。

$$C=\begin{bmatrix} I_{p} \\ 0 \end{bmatrix}, ~~C^\mathrm{T}=\begin{bmatrix} I_{p}~0 \end{bmatrix} \tag{2}$$

変換行列\(C\)にてヤコビ係数行列を圧力係数行列に縮小する。

$$\begin{align}J_{p}&=C^\mathrm{T}MJC \\ r_{p}&=C^\mathrm{T}MR\end{align}\tag{3}$$

ここで、シューア補行列(Shur compliment)を盛り込んだ前処理行列\(M\)を導入することでヤコビ係数行列全体から圧力ブロックを分離する。

$$M = \begin{bmatrix} I_{p} & – J_{ps}\hat{J}^{-1}_{ss} \\ 0 & I_{s} \end{bmatrix} \tag{4}$$

ここで、\(\hat{J}_{ss}\)行列は飽和率ブロック係数行列である。一般的に\(J_{ss}\)行列は密行列である為逆行列への変換は効率的ではない。そこで、対角行列\(\text{diag}(J_{ss})\)に単純近似し演算負荷を軽減する。

例えば、前処理行列\(M\)の設定方法としてIMPES離散化モデルを睨んだ設定も考えうる。次式のように、\(J_{ps}\)行列と\(J_{ss}\)行列の各列の和を要素とするベクトルで構成される\(\text{colsum}\)演算子を用いることで構築する。

$$M=\begin{bmatrix} I_{p} & – \text{colsum}(J_{ps}) \text{colsum}(J_{ss})^{-1} \\ 0 & I_{s} \end{bmatrix} \tag{4-1}$$

因みに、先の(3)式の圧力係数行列と残差ベクトルは(2)式から次式らに帰結する。

$$\begin{align}J_{p} &= \begin{bmatrix} I_{p} & 0 \end{bmatrix}
\begin{bmatrix} I_{p} & -J_{ps}\hat{J}^{-1}_{ss} \\ 0 & I_{s} \end{bmatrix}
\begin{bmatrix} J_{pp} & J_{ps} \\ J_{sp} & J_{ss} \end{bmatrix}
\begin{bmatrix} I_{p} \\ 0 \end{bmatrix} = J_{pp} – J_{ps}\hat{J}^{-1}_{ss}J_{sp} \\
r_{p} &= \begin{bmatrix} I_{p} & 0 \end{bmatrix}
\begin{bmatrix} I_{p} & -J_{ps}\hat{J}^{-1}_{ss} \\ 0 & I_{s} \end{bmatrix}
\begin{bmatrix} R_{p} \\ R_{s} \end{bmatrix} = R_{p} – J_{ps}\hat{J}^{-1}_{ss}R_{s}\end{align}\tag{5}$$

Practical CPR Algorithm

CPR前処理法を用いた反復法による線形方程式の解法は次の２段階プロセスにて実行する。CPR前処理法は順応一般化最小残差法(Flexible GMRES, FGMRES)の反復法と共に用いるのが一般的である。下記に具体的な算定手順を説明する。

１）圧力ブロックの残差ベクトル\(r_{p}\)を算定する
$$r_{p}=C^\mathrm{T}MR\tag{6-1}$$

２）圧力ブロック係数行列と残差ベクトルを満足する連立一次方程式から暫定圧力解\(x_{p}\)を求める
$$J_{p}x_{p}=r_{p}\tag{6-2}$$
圧力ブロックに置ける暫定圧力解はFGMRES反復法にて求める。この際、前処理行列(\(M_{p}\))はILU(0)やNF前処理法にて求める。この暫定圧力解を求めるプロセスを第１段の内部ループによる求解プロセスと云う。

３）先に得られた暫定圧力解にてヤコビ係数行列全体\(J\)を再計算する。加えて元の残差ベクトル\(r\)も再計算する

$$\tilde{r}=r-JCx_{p}\tag{6-3}$$

４）先にUpdateしたヤコビ係数行列と残差ベクトルから成る連立一次方程式から反復解\(\tilde{x}\)を求める。この暫定飽和率解を求めるプロセスを第２段の外部ループによる求解プロセスと云う。

$$\tilde{x}=\tilde{M}^{-1}\tilde{r}\tag{6-4}$$

５）得られた反復解\(\tilde{x}\)から解\(x\)を更新する

$$x=\tilde{x}+Cx_{p}\tag{6-5}$$

上記ステップ１〜５にて、圧力・飽和率解が許容範囲内に収束するまで、ILU(0)やNF前処理法にてFGMRES反復を繰り返すことになる。

一点、CPR前処理法は完全陰解法に対して有効な手法であるがIMPES解法に対しては適用できないことに留意する必要がある。

Summary

総括すると、CPR前処理法のアルゴリズムは強結合坑井モデル乃至順応陰解法モデルに於いて全体システムから坑井乃至圧力を其々を分離する手法と全く同手法である。異なる点は変換行列\(C\)を新たに組込む事、２段階ループによる反復計算が必要になる事である。

当然ながら、１段のNested FactorizationやILU Preconditioner前処理反復法よりも２段階のCPR前処理反復法の方が計算負荷は高いが、CPR前処理反復法は少ない反復回にて求解できる為、通常の１段反復法で収束性の劣る大規模なシミュレーションモデルに対しては有効な手法との検証が確認されている。尤も、NF前処理法との比較検証例は少なくその優位性は不明であるが、肌感として検討する価値は高いと思料。

しかしながら、俺の限定的で貧弱なTest Matrix環境下ではCPR前処理に基づく反復法の優位性を確認することは難しい。依って、現段階ではSMARTSOL Matrix Solution Packageには実装していない。腐るほど暇になったらSMARTSOLに盛り込む予定である。

実際のコーディングの具体例を記載できなかったことで今回の投稿は期待はずれの感は否めないが、以上にて終わりとする。

そろそろ山行に向けた鍛錬を始めようかな。

以上
Seize the day!

References

1. ECLIPSE Reservoir Simulator, Technical Description
2. Wallis, J. R., Kendall, R. P., and Little, T. E. 1985. Constrained residual acceleration of conjugate residual methods. In: SPE Reservoir Simulation Symposium. doi:10.2118/13536-MS.
3. M. Cusini, A. A. Lukyanov, J. R. Natvig, and H. Ha-jibeygi. Constrained pressure residual multiscale (CPR-MS) method for fully implicit simulation of multiphase flow in porous media. J. Comput. Phys., 299:472–486, doi: 10.1016/j.jcp.2015.07.019.
4. R. Piazza. 2019. Performance Assessment of Linear Solvers for Fully Implicit Reservoir Simulation. Dissertação de Mestrado – Departamento de Engenharia Mecânica, Pontifícia Universidade Católica do Rio de Janeiro.

〆