Linear and nonlinear solvers for simulating high-temperature multiphase flow within large-scale engineered subsurface systems

Multiphase flow simulation is well-known to be computationally demanding, and modeling large-scale engineered subsurface systems entails significant additional numerical challenges. These challenges arise from: (a) the presence of small-scale discrete features like shafts, tunnels, waste packages, and barriers; (b) the need to accurately represent both the waste form processes at the small spatial scale of the repository and the large-scale transport processes throughout heterogeneous geological formations; (c) the strong contrast in material properties such as porosity and permeability, as well as the nonlinear constitutive relations for multiphase flow; (d) high-temperature heat sources underground (e.g., due to decay of high level nuclear wastes) cause nearby water to boil off into steam, leading to a dry-out condition in porous media, with subsequent re-saturation of the heat source area. Numerical solution is based on discretization of the coupled system of nonlinear governing equations and solving a linear system of equations at each Newton-Raphson iteration. Practical problems require a very large number of unknowns that must be solved efficiently using iterative methods in parallel on high-performance computers. The unique challenges noted above can lead to an ill-conditioned Jacobian matrix and non-convergence with Newton’s method due to discontinuous nonlinearity in constitutive models. Moreover, practical applications such as nuclear repositories, carbon sequestration sites and geothermal reservoirs can require numerous Monte-Carlo simulations to explore uncertainly in material properties, geological heterogeneity, failure scenarios, or other factors; governmental regulatory agencies can mandate these as part of performance and safety assessments. Finally, some applications like nuclear waste repository require simulations over a million years. Hence there is a need for flexible, robust, and computationally efficient methods for multiphase flow in large-scale engineered subsurface systems. We apply the open-source simulator PFLOTRAN which employs a finite volume discretization and uses the PETSc parallel framework. We evaluate the performance of several preconditioners for the iterative solution of the linearized Jacobian system; these range from stabilized-biconjugate-gradient with block-Jacobi preconditioning (BCGS) to methods adopted from reservoir modeling, such as the constrained pressure residual (CPR) two-stage preconditioner and flexible generalized residual solver (FGMRES). We also implement within PETSc the general-purpose nonlinear solver, Newton trust-region dogleg Cauchy (NTRDC), which truncates the Newton update or modifies the update with a Cauchy solution that is within the quadratic model trust-region of the objective function and Newton trust-region (NTR). We demonstrate the effectiveness with two large-scale simulations for a series of test problems with increasing difficulty. In one numerical experiment, we find that the NTRDC and FGMRES-CPR-ABF (FCA) preconditioners generally perform best for the test problem having the most extreme nonlinear processes, achieving a 50x speed-up compared with BCGS. The most ill-conditioned and extreme nonlinear simulations do not converge with BCGS, but they do complete with NTRDC and FCA. In the other test problem, simulations with high-temperature heat sources causing extreme nonlinear processes with many state changes in the domain do not converge with conventional NR, but they do complete with the trust-region variants. We also investigate the strong scalability of each method and demonstrate the impact of node-packing upon parallel performance on modern processor architectures