SMARTSOL〜The AIM Implementation
第二回目は順応陰解法(AIM)の実装編をお届けします。AIM法の陽･陰解法への切り替えについて閾値法と局所行列安定法が提唱されている。This second thread presents an implementation of the Adaptive Implicit Method (AIM), in which the Threshold Method as well as the Local Matrix Stability Analysis Method are proposed as the switching criterion between the explicit and implicit AIM method.
Threshold Method(THR, 閾値法)
理論は至極単純である。要は、ある時刻に於ける独立変数の変化量が閾値を超える場合、陽解法を陰解法へ切り替える手法である。通常、坑井及びその隣接グリッドは陰変数として初期設定するが、それ以外のグリッドに対しては陽変数設定とする。シミュレーションの計算過程の独立変数の変化量から陽変数を陰変数へ切り替える事。欠点としては陰変数から陽変数への切り替え設定が出来ないことである。これに対応するため、下記局所行列安定法が考案される。
Local Matrix Stability Method(LMS, 局所行列安定法)
実装方法に関する具体的な理論は引用文献(References１)に基づく。
ファンら(1988年Fung et al)によって開発された方法で、局所的なヤコビ行列の振れ(Local Amplification Matrix)の数値安定性を基に切り替えを判断する方法である。陽変数から陰変数のみならず陰変数から陽変数へ双方向に切り替える事が可能となる。難点は閾値法に比べて演算負荷が多い事である。
多相流体流動方程式の線形化に伴うニュートン反復法は、多相流体流動方程式\(f_i^{n+1}\)、独立変数\(x_i^{n+1}\)から次式で表せる。
$$x_i^{n+1,k+1} = x_i^{n+1,k} – \left( \frac{\partial f_i^{n+1,k}}{\partial x_j^{n+1.k}} \right)^{-1} f_j^{n+1.k} \tag{7.5}$$
ここで、\(x_i^{n+1,k+1}\)は時刻\(n+1\)における変数\(x_i\)の第\(k+1\)近似解を示す。\(f_i^{n+1,k}\)は時刻\(n+1\)における第\(i\)方程式の第\(k\)近似解を示す。\(i\)は多相流体流動方程式若しくは独立変数を示す。
ニュートン法が収束した時点における離散化多相流体流動方程式は次式となる。
$$df_j^{n+1,k+1} = \left( \frac{\partial f_j^{n+1,k+1}}{\partial x_i^{n,m}} \right) dx_i^{n,m} + \left( \frac{\partial f_j^{n+1,k+1}}{\partial x_i^{n+1,k+1}} \right) dx_i^{n+1,k+1}=0 \tag{7.6}$$
ここで、添字\(m\)は時間\(n\)に於けるニュートン法の最終反復回数を示す。
時刻\(n+1\)と時刻\(n\)における演算誤差を其々\(\epsilon_i^{n}\)、\(\epsilon_i^{n+1}\)とすると、上式は
$$\left( \frac{\partial f_j^{n+1,k+1}}{\partial x_i^{n+1.k+1}} \right) \epsilon_i^{n+1} = \left( \frac{\partial f_j^{n+1,k+1}}{\partial x_i^{n,m}} \right) \epsilon_i^{n} \tag{7.7}$$
となる。
変換行列\(\mathbf{T} \equiv \left( \frac{\partial f_j^{n+1,k+1}}{\partial x_i^{n+1.k+1}} \right)^{-1} \left( \frac{\partial f_j^{n+1,k+1}}{\partial x_i^{n,m}} \right)\)を導入して次式の演算誤差を得る。
$$\epsilon_i^{n+1} = \mathbf{T} \epsilon_i^{n} \tag{7.8}$$
任意の時間幅における数値安定性を確保する為の充分条件は次式となる。
$$\|\epsilon_i^{n+1}\| \leq \|\epsilon_i^{n}\| \tag{7.9}$$
依って、変換行列\(\mathbf{T}\)が正規であれば、同充分条件の必要条件は次式で与えられる。
$$\rho(T) := {\max_{1≦i≦n_e} |\alpha_i(\mathbf{T})|} \leq 1 \tag{7.10}$$
ここで、\(\rho\)は変換行列\(\mathbf{T}\)のスペクトラム半径、\(\alpha_i\)は変換行列\(\mathbf{T}\)の固有値である。スペクトル半径とは固有値の絶対値の最小上界である。つまり「任意の\(a \in \mathbf A\)に対して\(a \le m\)を満たす\(m\in \mathbb R\)のうち最小のもの」を示す。又、独立変数\(x_i\)の総数\(n_e\)は三相ブラックオイル貯留層シミュレータにおいて３変数である。
任意の格子における数値安定性について限定すると、式(7.5)は次式に書き換えられる。
$$\left( \frac{\partial f_l^{n+1,k+1}}{\partial x_i^{n+1.k+1}} \right) \epsilon_i^{n+1} = \left( \frac{\partial f_l^{n+1,k+1}}{\partial x_i^{n,m}} \right) \epsilon_i^{n} \tag{7.11}$$
ここで、\(\ell\)は相(油、ガス、水)を示す。
即ち、任意の時間幅、格子における数値安定性を確保するには、先の充分条件(7.9)並びに(7.10)から次式が成立しなければならない。
$$\begin{align} \rho(T) &:= {\max_{1≦i≦n_e} |\alpha_i(\mathbf{T})|} \leq 1 \\
&\|\epsilon_\ell^{n+1}\| \leq \|\epsilon_\ell^{n}\| \end{align} \tag{7.12}$$
同条件を実際に適用するに当って、多相流体流動方程式\(f_i^{n+1}\)に於ける時刻\(n\)の蓄積項を無視する必要がある。
具体的には、時間幅\(\Delta t\)にて時間差分近似によって得られた対象格子における次式(7.11-1)蓄積項に関して、時刻\(n\)における蓄積項を無視し、
$$\begin{align} \int_{V_{i,j,k}} \frac{\partial}{\partial t} \Big (\varphi \frac{S_\ell}{B_\ell} \Big )\partial V &\approx \frac{V_{i,j,k}}{\Delta t}\Delta t \Big (\varphi \frac{S_\ell}{B_\ell} \Big )_{i,j,k} \\
&\approx \frac{V_{i,j,k}}{\Delta t} \Big [ \Big (\varphi \frac{S_\ell}{B_\ell} \Big )^{n+1} - \Big (\varphi \frac{S_\ell}{B_\ell} \Big )^{n} \Big ]_{i,j,k} \end{align} \tag{7.13-1}$$
次式の時刻\(n+1\)のみの蓄積項に帰結する。
$$\int_{V_{i,j,k}} \frac{\partial}{\partial t}\Big (\varphi \frac{S_\ell}{B_\ell} \Big )\partial V \approx \frac{V_{i,j,k}}{\Delta t}\Delta t \Big (\varphi \frac{S_\ell}{B_\ell} \Big )_{i,j,k} \approx \frac{V_{i,j,k}}{\Delta t} \Big (\varphi \frac{S_\ell}{B_\ell} \Big )^{n+1}_{i,j,k} \tag{7.13-2}$$
ここで、\(V\)対象格子の容積、孔隙率\(\varphi\)、相\(\ell\)飽和率\(S_\ell\)、相\(\ell\)容積係数\(B_\ell\)とする。
詳しくは、引用文献１に於けるAppendix (A-10)及び(A-11)式を参照されたし。
Implementation
自家製シミュレータ”SMART”にはAIM閾値法のみ組み込んでいるが、”iSMART”にはAIM閾値法と局所行列安定法の双方を実装している。これは”SMART”のヤコビ行列は解析微分によって導出している為、数値微分法を用いる”iSMART”に比べて収束性は格段に良く完全陰解法でも充分に実用的であり敢えてAIM法を組み込む必要がないものの、完全陰解法に加えてIMPES法も標準で組み込んでいる為、AIM閾値法を興味本位に組み込んでいる。
一方、”iSMART”は初歩的なシミュレータで機能はPrimitiveである。ヤコビ行列は有限差分法による数値微分法にて導出している関係で収束性は悪い。そのため演算速度を高める手法として、１変数の独立変数（坑底圧力）を基づく坑井モデルの簡素化に加えて、AIM法による演算負荷の軽減並びに速度の向上を追求している。
実装に当たっての留意点を下記に纏める。SPE Comparative Solution Projectの検証結果についても付記する。
Threshold Method
独立変数(飽和率、飽和圧力、ガス由比Rs/Rv)の時刻変化量に対する閾値を設け、同閾値を超える変化量のグリッドに対しては陽から陰解法に転換する。
- IMEXの閾値設定は独立変数の時刻変化量(NORM)に対応させる。IMEX NORM Default設定は油層圧力、飽和率=435.0psia、0.1である。一方、iSMART及びSMARTでは収束性と統合性の見地からECLIPSE100のIMPES設定(200.0psia、0.1)に対応させる。
- IMEX Default設定の閾値は飽和圧力、飽和率=0.0125、0.0125である。iSMARTの独立変数はIMEX同等、閾値は同じく0.125/0.125に設定する。一方、SMARTの独立変数はECLIPSE100と同じく飽和率、Rs/Rvであるが、飽和率のみを閾値対象として0.01と若干低減する。
Local Matrix Stability Method
- iSMARTの具体的な実装では、時刻\(n+1\)、\(n\)に置けるヤコビ行列を検証するに当たり、IMPES法にて変換行列\(\mathbf{T}\)を構築する。その上でその数値安定性から陽⇔陰解法の転換をチェックする。この前提が正しいかどうかは文献からは同定できないものの、演算過程に置ける局所行列安定性の挙動から妥当と判断している。しかし、引用文献１のような計算結果が得られていない。
- 固有値の演算負荷は高い事もあって、時間ステップ10回毎にLocal Matrix Stabilityの算定を陰解法グリッドのFrontのみチェックする事として演算負荷の低減に努めている。
- 陽転換係数は0.8として陰転換係数1.0に比べて軽減して発振誤差を低減する。
Result
| SPE Case | Implicit | CPU Time Intel®Core™ i3 M370 | Total No. of Timesteps | % of Implicit Grids |
| SPE 1 | FIM | 0.5 sec | 84 | 100 |
| – | AIM-LMS | 1.1 sec | 48 | 100 |
| – | AIM-THR | 1.0 sec | 52 | 71 |
| SPE 2 | FIM | 0.5 sec | 34 | 100 |
| – | AIM-LMS | 1.0 sec | 92 | 97 |
| – | AIM-THR | 0.7 sec | 65 | 85 |
| SPE 9 | FIM | 47.0 sec | 66 | 100 |
| – | AIM-LMS | 70.9 sec | 68 | 100 |
| – | AIM-THR | 525 sec | 147 | 96 |
| Case | Implicit | CPU Time Intel®Core™ i3 M370 | Total No. of Timesteps | % of Implicit Grids |
| SPE 1 | FIM | 0.67 sec | 40 | 100 |
| – | AIM-THR | 0.74 sec | 51 | 62 |
| SPE 2 | FIM | 0.27 sec | 34 | 100 |
| – | AIM-THR | 0.38 sec | 51 | 53 |
| SPE 9 | FIM(ILU) | 40.2 sec | 47 | 100 |
| – | FIM(NF) | 8.5 sec | 42 | 100 |
| – | AIM-THR | 52.0 sec | 57 | 41 |
| SPE Case | Implicit | CPU Time Intel®Core™ i3 M370 | Total No. of Timesteps | % of Implicit Grids |
| SPE 9 | FIM | 11.7 sec | 61 | 100 |
| – | AIM-LMS | 11.3 sec | 41 | 58 |
| – | AIM-THR | 10.4 sec | 42 | 35 |
Conclusion
自家製シミュレータの実装AIM法の演算結果を踏まえ以下、所見を纏める。
算法からの推論
- AIM法は完全陰解法に比べて時刻ステップが増える為、ニュートン法の総逐次反復法は多くなろう
- 陰解法格子のFrontのみに対して固有値の算定対象としているが、固有値算定の演算負荷は高いので、AIM閾値法の方がAIM局所行列安定法に比べて、より実用的な手法と考えられる
- 更に、AIM法に関わる行列計算処理の負荷を考慮するとAIM法の効果は限定的とも推察
計算結果の総括
- iSMART: 全ケースに於いてFIM法が優る。更に、不均質性の顕著な大規模な課題(SPE9)に対してはAIM閾値法は局所行列安定法に比べて演算負荷は高くなり速度は低下する。特に、数値微分を用いている事もあって、ILU前処理法はNF前処理法に比べて大幅にPerformanceが低下する。対応策としてIncidence MatrixのScaling法に付いて検討したが顕著な効果なし。CPR法が有効かもしれない
- SMART: FIM法はIMEX、ECLIPSEと比べても遜色ないPerformanceである。
演算速度はFIM法がAIM法を優る。ILU前処理法に比べてNF前処理法の方が収束性がRobustである。iSMART同様、不均質性の著しい大規模な課題(SPE9)に対してはILU前処理法のPerformanceはNF前処理法に比べて劣る - IMEX: FIM法、AIM法双方とも演算速度は変わりなし。
AIM-THR法はAIM-LSM法と遜色ない演算速度で納得できる。しかし、FIM法のTimestepがAIM法に比べて大幅に上回るのは納得できぬ
考察
FIM法はAIM法に比べて演算速度に変わりなし。行列解法をFIM法に特化することでFIM法の演算速度が更に向上する可能性がある。詰まり、ECLIPSEがFIM法に拘る証と認識する。
iSMARTに自動微分を実装後に再度レビューしたい。
以上
Seize the day!
References
- Larry S-K. Fung, David A. Collins, and Long H. Nghiem. An adaptive-implicit switching criterion based on numerical stability analysis. SPE16003, SPE Reservoir Engineering, 4(1):45-51, 1989.
- P.A. Forsyth Jr. and P.H. Sammon. Practical considerations for adaptive implicit methods in reservoir simulation. Journal of Computational Physics, 62:265–281, 1986.
- G.W. Thomas and D.H. Thurnau. Reservoir simulation using an adaptive implicit method. SPE-AIME, 23:759–768, 1983.
〆