Abstract
Two efficient and scalable numerical solution methods will be compared using exact Jacobians to solve the fully coupled Newton systems arising during fully implicit discretization of the equations for two-phase flow in porous media. These methods use algebraic multigrid (AMG) to solve the linear systems in every Newton step. The algebraic multigrid methods rely on (i) a Schur Complement Reduction (SCR-AMG) and (ii) a Constrained Pressure Residual method (CPR-AMG) to decouple elliptic and hyperbolic contributions. Both methods employ automatic differentiation (AD) to calculate exact Jacobians and are compared with finite difference (FD) approximations of the Jacobian. The superiority of AD is shown by several numerical test cases from the field of CO2 geo-sequestration comprising two- and three-dimensional examples. A weak scaling test on JUQUEEN, a BlueGene/Q supercomputer, demonstrates the efficiency and scalability of both methods. To achieve maximum comparability and reproducibility, the Portable Extensible Toolkit for Scientific Computation (PETSc) is used as framework for the comparison of all solvers.
Article PDF
Similar content being viewed by others
References
Aziz, K., Settari, A.: Petroleum Reservoir Simulation. Applied Science Publishers, London (1979)
Balay, S., Gropp, W.D., McInnes, L.C., Smith, B.F.: Modern Software Tools in Scientific Computing, Birkhȧuser Arge, E., Bruaset, A.M., Langtangen, H.P. (eds.) MA, Boston (1997)
Balay, S., Abhyankar, S., Adams, M.F., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Dener, A., Eijkhout, V., Gropp, W.D., Karpeyev, D., Kaushik, D., Knepley, M.G., May, D.A., McInnes, L.C., Mills, R.T., Munson, T., Rupp, K., Sanan, P., Smith, B.F., Zampini, S., Zhang, H., Zhang, H.: PETSc users manual Tech. Rep. ANL-95/11 - Revision 3.13, Argonne National Laboratory, Lemont, IL, USA, https://www.mcs.anl.gov/petsc (2020)
Balay, S., Abhyankar, S., Adams, M.F., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Dener, A., Eijkhout, V., Gropp, W.D., Karpeyev, D., Kaushik, D., Knepley, M.G., May, D.A., McInnes, L.C., Mills, R.T., Munson, T., Rupp, K., Sanan, P., Smith, B.F., Zampini, S., Zhang, H., Zhang, H.: PETSc Web page http://www.mcs.anl.gov/petsc (2020)
Bastian, P.: Numerical Computation of Multiphase Flows in Porous Media. Habilitation thesis, Technische Fakultät, Christian-Albrechts-Universität Kiel https://conan.iwr.uni-heidelberg.de/data/people/peter/pdf/Bastian_habilitationthesis.pdf (1999)
Bastian, P., Helmig, R.: Efficient fully-coupled solution techniques for two-phase flow in porous media: Parallel multigrid solution and large scale computations. Adv. Water Resour. 23(3), 199–216 (1999). https://doi.org/10.1016/S0309-1708(99)00014-7
Bear, J.: Hydraulics of Groundwater McGraw-Hill New York (1979)
Brooks, R.J., Corey, A.T.: Hydraulic Properties of Porous Media, vol 3 Colorado State University Hydrology Paper Fort Collins CO (1964)
Buckley, S.E., Leverett, M.C.: Mechanism of fluid displacement in sands. Transactions of the AIME 146(01), 107–116 (1942)
Bui, Q.M., Elman, H.C., Moulton, J.D.: Algebraic Multigrid Preconditioners for Multiphase Flow in Porous Media. SIAM Journal on Scientific Computing 39(5), S662–S680 (2017). https://doi.org/10.1137/16M1082652
Burdine, N.T.: Relative Permeability Calculations from Pore-Size Distribution Data. Petroleum Transactions AIME 198, 71–77 (1953). https://doi.org/10.2118/225-G
Büsing, H.: Efficient Solution Techniques for Multi-phase Flow in Porous Media. In: Lirkov I., Margenov, S., Large-Scale Scientific Computing, Springer International Publishing, Cham, Switzerland. https://doi.org/10.1007/978-3-319-73441-5_63, pp. 572–579 (2018)
Chavent, G., Jaffrè, J.: Mathematical Models and Finite Elementes for Reservoir Simulation North-Holland Amsterdam (1986)
Clauser, C.: Numerical Simulation of Reactive Flow in Hot Aquifers: SHEMAT and processing SHEMAT Springer Heidelberg-Berlin (2003)
Darcy, H.: Les Fontaines Publiques de la Ville de Dijon Dalmont Paris (1856)
Dawson, C.N., Klíe, H., Wheeler, M.F., Woodward, C.S.: A parallel, implicit, cell-centered method for two-phase flow with a preconditioned Newton–K,rylov solver. Computational Geosciences 1(3-4), 215–249 (1997). https://doi.org/10.1023/A:1011521413158
Douglas, J. Jr., Peaceman, D., Rachford, H. Jr.: A Method for Calculating Multi-Dimensional Immiscible Displacement. Transactions of the Metallurgical Society of AIME 216, 297–308 (1959)
Falgout, R.D., Yang, U.M.: hypre: A Library of High Performance Preconditioners Sloot PMA Hoekstra AG Tan CJK Dongarra JJ eds Computational Science — ICCS 2002 Springer Heidelberg-Berlin (2002)
Falgout, R.D., Jones, J.E., Yang, U.M. Bruaset, A.M., Tveito, A. (eds.): The Design and Implementation of hypre, a Library of Parallel High Performance Preconditioners. Springer, Heidelberg-Berlin (2006). https://doi.org/10.1007/3-540-31619-1_8
van Genuchten, M.T.: A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils. Soil Science Society of America Journal 44, 892–898 (1980). https://doi.org/10.2136/sssaj1980.03615995004400050002x
Gries, S.: On the Convergence of System-AMG in Reservoir Simulation. SPE J. 23(02), 589–597 (2018). https://doi.org/10.2118/182630-PA
Gries, S., Stüben, K., Brown, G.L., Chen, D., Collins, D.A., et al.: Preconditioning for Efficiently Applying Algebraic Multigrid in Fully Implicit Reservoir Simulations. SPE J. 19(04), 726–736 (2014). https://doi.org/10.2118/163608-PA
Griewank, A.: Evaluating Derivatives: Principles and Techniques of Algorithmic Differentiation Society for Industrial and Applied Mathematics SIAM Philadelphia PA (2000)
Hackbusch, W.: Multi-Grid Methods and Applications Springer Heidelberg-Berlin (1985)
Hascoet, L., Pascual, V.: The Tapenade Automatic Differentiation tool: principles, model, and specification. ACM Transactions on Mathematical Software TOMS 39(3), 20 (2013). https://doi.org/10.1145/2450153.2450158
Henson, V., Yang, U.: BoomerAMG: A parallel algebraic multigrid solver and preconditioner. Appl. Numer. Math. 41(1), 155–177 (2002). https://doi.org/10.1016/S0168-9274(01)00115-5
Jenny, P., Lee, S., Tchelepi, H.: Adaptive fully implicit multi-scale finite-volume method for multi-phase flow and transport in heterogeneous porous media. J. Comput. Phys. 217(2), 627–641 (2006). https://doi.org/10.1016/j.jcp.2006.01.028
Kayum, S., Cancelliere, M., Rogowski, M., Al-Zawawi, A.: Application of Algebraic Multigrid in Fully Implicit Massive Reservoir Simulations. In: SPE Europec featured at 81st EAGE Conference and Exhibition, London, England, UK, Society of Petroleum Engineers Richardson. https://doi.org/10.2118/195472-MS. TX, USA (2019)
Klíe, H.M., Ramè, M., Wheeler, M.: Two-stage Preconditioners for inexact Newton Methods in multi-phase reservoir simulation. http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.30.7551&rep=rep1&type=pdf (1996)
Klíe, H.M., Monteagudo, J., Hoteit, H., Rodriguez, A.A.: Towards a New Generation of Physics Driven Solvers for Black Oil and Compositional Flow Simulation. In: SPE Reservoir Simulation Symposium, The Woodlands, TX, USA, Society of Petroleum Engineers, Richardson, TX, USA. https://doi.org/10.2118/118752-MS (2009)
Lacroix, S., Vassilevski, Y.V., Wheeler, M.F.: Decoupling preconditioners in the implicit parallel accurate reservoir simulator IPARS. Numerical Linear Algebra with Applications 8(8), 537–549 (2001). https://doi.org/10.1002/nla.264
Lacroix, S., Vassilevski, Y., Wheeler, J., Wheeler, M.: Iterative Solution Methods for Modeling Multiphase Flow in Porous Media Fully Implicitly. SIAM J. Sci. Comput. 25(3), 905–926 (2003). https://doi.org/10.1137/S106482750240443X
Mary, T.: Block low-rank multifrontal solvers : complexity, performance, and scalability. PhD thesis, Université Paul Sabatier - Toulouse III https://tel.archives-ouvertes.fr/tel-01929478 (2017)
Mishev, I.D., Shaw, J.S., Lu, P.: Numerical Experiments with AMG Solver in Reservoir Simulation. In: SPE Reservoir Simulation Symposium, The Woodlands, TX, USA, Society of Petroleum Engineers, Richardson, TX, USA. https://doi.org/10.2118/141765-MS (2011)
Mualem, Y.: A New Model for Predicting the Hydraulic Conductivity of Unsaturated Porous Media. Water Resour. Res. 12, 513–522 (1976). https://doi.org/10.1029/WR012i003p00513
Pape, H., Clauser, C., Iffland, J.: Permeability prediction based on fractal pore-space geometry. Geophysics 64(5), 1447–1460 (1999). https://doi.org/10.1190/1.1444649
Rall, LB: Automatic Differentiation: Techniques and Applications Springer New York (1981)
Rath, V., Wolf, A., Bücker, M.: Joint three-dimensional inversion of coupled groundwater flow and heat transfer based on automatic differentiation: sensitivity calculation, verification and synthetic examples. Geophys. J. Int. 167, 453–466 (2006). https://doi.org/10.1111/j.1365-246X.2006.03074.x
Saad, Y.: A Flexible Inner-Outer Preconditioned GMRES Algorithm. SIAM J. Sci. Comput. 14 (2), 461–469 (1993). https://doi.org/10.1137/0914028
Saad, Y.: Iterative Methods for Sparse Linear Systems 2nd edn Society for Industrial and Applied Mathematics Philadelphia, PA (2003)
Saad, Y., Schultz, M.H.: GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems. SIAM Journal on Scientific and Statistical Computing 7(3), 856–869 (1986). https://doi.org/10.1137/0907058
Scheichl, R., Masson, R., Wendebourg, J.: Decoupling and block preconditioning for sedimentary basin simulations. Computational Geosciences 7 (4), 295–318 (2003). https://doi.org/10.1023/B:COMG.0000005244.61636.4e
Singh, V.P., Cavanagh, A., Hansen, H., Nazarian, B., Iding, M., Ringrose, P.S., et al.: Reservoir Modeling of CO2 Plume Behavior Calibrated Against Monitoring Data From Sleipner, Norway. In: SPE Annual Technical Conference and Exhibition, Florence, Italy, Society of Petroleum Engineers Richardson TX USA. https://doi.org/10.2118/134891-MS (2010)
SPE: SPE 10th Comparative Solution Project Description of Model 2: SPE Project Website: http://www.spe.org/web/csp/datasets/set02.htm Richardson TX (2001)
Stüben, K., Clees, T., Klie, H., Lu, B., Wheeler, M.F.: Algebraic Multigrid Methods (AMG) for the Efficient Solution of Fully Implicit Formulations in Reservoir Simulation. In: SPE Reservoir Simulation Symposium, Houston, TX, USA, Society of Petroleum Engineers, Richardson, TX, USA. https://doi.org/10.2118/105832-MS (2007)
Van der Vorst, H.A.: Bi-CGStab: A Fast and Smoothly Converging Variant of Bi-CG for the Solution of Nonsymmetric Linear Systems. SIAM J. Sci. Stat. Comput. 13(2), 631–644 (1992). https://doi.org/10.1137/0913035
Wallis, J., et al.: Incomplete Gaussian Elimination as a Preconditioning for Generalized Conjugate Gradient Acceleration. In: SPE Reservoir Simulation Symposium, San Francisco, CA, Society of Petroleum Engineers, Richardson, TX, USA. https://doi.org/10.2118/12265-MS (1983)
Wallis, J.R., Kendall, R., Little, T., et al.: Constrained Residual Acceleration of Conjugate Residual Methods. https://doi.org/10.2118/13536-MS (1985)
Funding
Open Access funding provided by Projekt DEAL.
Author information
Authors and Affiliations
Corresponding author
Additional information
Publisher’s note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The research leading to these results has received funding from the European Union’s Horizon2020 Research and Innovation Program under grant agreement No. 640573 (Project DESCRAMBLE) and No. 676629 (Project EoCoE)
Appendices
Appendix 1: PETSc command line options for the different solvers
This appendix shows the different command line options used to invoke the linear solvers in PETSc.
The following options for the weak scaling test are shared by all the solvers.
-ksp_atol 1e-50 -ksp_rtol 1e-6 -ksp_max_it 100 -ksp_type fgmres
1.1 SCR-AMG
-pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition a11 -fieldsplit_0_ksp_type preonly -fieldsplit_0_pc_type hypre -fieldsplit_0_pc_hypre_type boomeramg -fieldsplit_1_ksp_type gmres -fieldsplit_1_pc_type hypre -fieldsplit_1_pc_hypre_type boomeramg -fieldsplit_1_ksp_max_it 10 -fieldsplit_1_ksp_rtol 1e-2
1.2 CPR-AMG
-pc_type composite -pc_composite_type multiplicative -pc_composite_pcs fieldsplit,bjacobi -sub_0_ksp_type fgmres -sub_0_pc_fieldsplit_type additive -sub_0_fieldsplit_0_ksp_type gmres -sub_0_fieldsplit_0_pc_type hypre -sub_0_fieldsplit_0_pc_hypre_type boomeramg -sub_0_fieldsplit_1_ksp_type preonly -sub_0_fieldsplit_1_pc_type none -sub_1_sub_pc_type ilu
Appendix 2: Code verification using the Buckley-Leverett problem
The Buckley-Leverett problem describes the immiscible displacement of oil by water in a porous medium (cf. [9]). It is widely used for the verification of numerical models. Figure 13 shows the analytical solution for the problem and compares it with numerical solutions with different cell sizes. It is clear that the numerical solution approaches the analytical for increasingly smaller cell sizes.
Appendix 3: Discussion of gravity numbers
The gravitational number is defined as:
and relates gravitational and viscous forces. Here, \({v}_{\text {cr}}=\frac {\phi l_{\text {cr}}}{t_{\text {cr}}}\) is the characteristic velocity with a characteristic length lcr and a characteristic time tcr. Permeability is a scalar, which means that it is assumed to be isotropic.
For our test case 1, we have a reservoir thickness of 100 m, thus lcr = 100m. Simulation time is 67 days and consequently tcr = 5, 788, 800s. Densities are chosen as ρn = 450 kg m− 3 and ρw = 1000 kg m− 3. Non-wetting viscosity is μn = 3.0 × 10− 5 Pa s and permeability is \(\mathbb {K}=10^{-12}\) m2. All in all, we have:
In addition, we consider two test cases: (i) injection into a hot aquifer with ρn = 250kg m− 3 and μn = 3.0 × 10− 5 Pa s and (ii) injection into a deep aquifer with ρn = 650kg m− 3 and μn = 8.0 × 10− 5 Pa s. This leads to gravitational numbers of (i) Gr ≈ 70.99 and (ii) Gr ≈ 12.42.
Comparing Table 1 with Table 6, we have a higher number of Newton and linear iterations for the hot aquifer. This is in line with the higher gravitational number of Gr ≈ 70.99 compared with Gr ≈ 52.06 for the original test case.
Table 7 shows results for the deep aquifer with a low gravitational number of Gr ≈ 12.42. Here, we have a lower number of Newton and linear iterations. Again, this is in line with the lower gravitational number. Analogous results hold for the diffusion-dominated case.
Gravity numbers for the Sleipner and SPE10B case are Gr ≈ 88.96 and Gr ≈ 555.20, respectively. To calculate a single gravitational number, we used mean permeabilities and porosities. The high gravitational number of the SPE10B case reflects its high degree of difficulty.
Appendix 4: Discussion of CFL numbers
The Courant–Friedrichs–Lewy (CFL) condition can be derived using the saturation equation from the fractional flow formulation:
Here, vtot = vw + vn is the total velocity and \(f_{w}=\frac {\lambda _{w}}{\lambda _{w}+\lambda _{n}}\) is the fractional flow function. Then, the CFL number C is defined as:
Figure 14 shows the CFL number for every time step for the advection- and diffusion-dominated two-dimensional test case for the original, hot, and deep aquifer. The figure shows results from the finite difference approximation of the Jacobian and SCR-AMG.
For the advection-dominated case, we have CFL numbers of over 100 in the 12th time step. For the hot aquifer, CFL numbers reach values of approximately 125 also in the 12th time step and for the deep aquifer values of over 50 in the 13th time step. Similarly, for the diffusion-dominated case, we have CFL numbers of over 140 in the 13th time step. For the hot aquifer, CFL numbers reach a value of over 105 in the 15th time step and for the deep reservoir over 110 in the 13th time step. The high CFL numbers indicate that the fully implicit fully coupled approach is indeed well justified.
Appendix 5: Influence of anisotropy
In this section, the influence of anisotropy on solver performance is examined. The diffusion-dominated two-dimensional test case from Section 4 is the basis for the comparison. Results for anisotropy factors of 1:100 and 1:10,000 are shown in Tables 8 and 9, respectively.
Comparing Table 8 with Table 2, performance of CPR-AMG decreases. An increase in the number of time steps (factor of 4.21 for FD and 4.69 for AD), Newton iterations, and linear iterations is observed. The performance of SCR-AMG is affected less by anisotropy with a slight increase in the number of time steps (factor of 1.41 for FD and 1.50 for AD), Newton iterations, and linear iterations.
Looking at Table 9, performance of CPR-AMG drastically decreases, both for the FD and the AD case. The number of time steps increases by a factor of 31.90 (FD) and 34.19 (AD) compared with the original values with no anisotropy. The performance of SCR-AMG is significantly more stable. The number of time steps increases only by a factor of 1.59 (FD) and 1.77 (AD) compared with the values with no anisotropy.
In summary, CPR-AMG is affected much more by anisotropy compared with SCR-AMG. Interestingly, this is also true for the advection-dominated case, where CPR-AMG originally outperformed SCR-AMG.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
About this article
Cite this article
Büsing, H. Efficient solution techniques for two-phase flow in heterogeneous porous media using exact Jacobians. Comput Geosci 25, 163–177 (2021). https://doi.org/10.1007/s10596-020-09995-w
Received:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s10596-020-09995-w
Keywords
- Algebraic multigrid (AMG)
- Schur complement reduction (SCR-AMG)
- Constrained pressure residual (CPR-AMG)
- Multiphase flow in porous media
- Automatic differentiation (AD)
- CO2 geo-sequestration