Application of accelerated fixed-point algorithms to hydrodynamic well-fracture coupling

The coupled simulations of dynamic interactions between the well, hydraulic fractures and reservoir have significant importance in some areas of petroleum reservoir engineering. Several approaches to the problem of coupling between the numerical models of these parts of the full system have been developed in the industry in past years. One of the possible approaches allowing formulation of the problem as a fixed-point problem is studied in the present work. Accelerated Anderson's and Aitken's fixed-point algorithms are applied to the coupling problem. Accelerated algorithms are compared with traditional Picard iterations on the representative set of test cases including ones remarkably problematic for coupling. Relative performance is measured, and the robustness of the algorithms is tested. Accelerated algorithms enable a significant (up to two orders of magnitude) performance boost in some cases and convergent solutions in the cases where simple Picard iterations fail. Based on the analysis, we provide recommendations for the choice of the particular algorithm and tunable relaxation parameter depending on anticipated complexity of the problem.


page 1

page 2

page 3

page 4


An Acceleration of Fixed Point Iterations for M/G/1-type Markov Chains by Means of Relaxation Techniques

We present some accelerated variants of fixed point iterations for compu...

A link between the steepest descent method and fixed-point iterations

We will make a link between the steepest descent method for an unconstra...

Fast K-Means Clustering with Anderson Acceleration

We propose a novel method to accelerate Lloyd's algorithm for K-Means cl...

Effective Response of Heterogeneous Materials using the Recursive Projection Method

This paper applies the Recursive Projection Method (RPM) to the problem ...

A family of fast fixed point iterations for M/G/1-type Markov chains

We consider the problem of computing the minimal nonnegative solution G ...

IFOSMONDI Co-simulation Algorithm with Jacobian-Free Methods in PETSc

IFOSMONDI iterative algorithm for implicit co-simulation of coupled phys...

On the Relationship Between Coupling and Refactoring: An Empirical Viewpoint

[Background] Refactoring has matured over the past twenty years to becom...

1 Introduction

In this paper, we study the robustness and performance of the algorithms for solution of the fixed-point problem in application to pressure-rate coupling between computational models used in petroleum reservoir simulation. Anderson’s AndersonOrig algorithm and a generalization of Aitken’s aitken1927xxv algorithm for multidimensional fixed-point problems, referred to as accelerated algorithms below, are compared to reference Picard iterations picard1890memoire.

Simulation is a common tool to predict and optimize production of hydrocarbons. Often, simulations are focused on either reservoir or well part of the problem, evaluating response of the other part using a simplified steady-state analytical model. This approach may be justified by the huge gap between characteristic times for well and reservoir. However, depending on the nature of production instabilities and transients, temporal and spatial time scales of phenomena occurring in a well and reservoir itself may significantly overlap nennie2007investigation. Accordingly, the coupled simulations taking into account dynamic interaction between well and reservoir and incorporating numerical models for both well and reservoir are needed in certain areas of petroleum reservoir engineering. For the review of particular applications one can refer to nennie2007investigation; sagen2011dynamic.

In dasilva2015review the review of coupled well-reservoir simulations and classification of the possible coupling methods have been presented. Here we follow this classification subdividing the coupling algorithms into fully implicit, implicit and explicit. The fully implicit

approach implies merging of two sets of partial differential equations underlying the numerical simulators into the single set and common discretization. This approach appears to be the most rigorous and is quite commonly used; for example, see 

holmes1998application. However, it may be impractical from the software development perspective, because it would demand re-implementation of major parts of the simulators to be coupled. In the explicit approach, the well and reservoir models are solved sequentially, there are no coupling iterations performed within a timestep and essentially there is a mismatch between pressures and rates in the order of one timestep. This approach seems to be rarely used in well-reservoir coupling and will not be further discussed here. In the present work we adopt the implicit coupling approach that iterates pressure-rate boundary conditions between well and reservoir simulators until they satisfy certain tolerance on each timestep. This approach is described in detail in Section 2. For more examples on applications of this approach, refer to nennie2007investigation and bhat2005coupling.

The coupling problem in the implicit approach may be cast to the form of the fixed-point problem that can be solved by Picard iterations picard1890memoire. The fixed point iteration method (also known as the Picard algorithm) seems to be often selected when it is necessary to quickly implement a simple algorithm to solve a non-linear matrix equation. From our past experience ChertovChaplygin2019, a simple Picard algorithm may converge very slowly. It also depends on an arbitrary parameter (relaxation factor), which is selected via trial and error. An unfortunate choice of the relaxation factor slows down the convergence and may often result in divergence of the numerical solution ChertovChaplygin2019. Manual fine tuning of relaxation parameter is very inconvenient, especially when it is necessary to ensure convergence of all cases in a massive parametric study with many simulator runs.

To our knowledge, application of the accelerated algorithms has been never reported in literature on petroleum well-reservoir simulations. On the contrary, the acceleration methods are used within other domains of multiphysics modeling requiring coupling. For instance, one can mention the broad class of fluid-structure interaction (FSI) problems. Fully implicit, implicit and explicit classes of methods mentioned above correspond to monolithic, partitioned and partitioned staggered classes in FSI terminology degroote2008stability. Aitken relaxation found application within the partitioned approach allowing formulation as a fixed-point problem mok2001accelerated; AitkenUR.

The paper is organized as follows. Following this introduction, Section 2 gives a brief description of the numerical simulation framework for which the coupling algorithm is needed. The coupling algorithm is recast mathematically as the fixed-point problem. In the numerical framework we consider, the fixed-point problem is complicated by the unavailability of analytical derivatives of the target function and high computational cost to evaluate the target function itself, which scales with the size of a numerical problem. As the coupling solver is an integral part of the larger simulation tool, it would also be time-saving if the method did not require tuning of its parameters for each particular case. In Section 3, we describe advanced Anderson’s and Aitken’s iterative algorithms and document some practical aspects of their implementation in our in-house code. These algorithms turn out to be the most suitable methods for coupling in our case given the constraints on the computational cost of the target function and unavailability of the derivatives. We also address the theoretically proven convergence properties of these methods and review some examples of previous experience of their usage in numerical modelling. In section 4, we evaluate in detail the performance of the implemented Anderson’s and Aitken’s iterative algorithms and compare them against the conventional fixed-point (Picard) iterations. The performance of these methods is benchmarked and demonstrated using two numerical problems that are of practical interest for the oil and gas industry and that heavily depend on the well-fracture coupling. The first problem is the cascading failure of hydraulic fractures, where failure of one fracture caused by excessively high flow rates can trigger failure of subsequent fractures. The second problem is the severe slugging during the multiphase flow in a horizontal well with multiple fractures. The cascading failure problem is a particularly hard case for the Picard iterations algorithm and highlights the benefits of the accelerated algorithms. The slugging cases evaluate robustness and performance of the algorithms for a wide range of system parameters and flow regimes. The detailed study of performance provides guidance and recommendations on optimal selection of the control settings of coupling algorithms. Finally, Section 5 gives a summary of key results obtained in this work.

2 Problem overview

2.1 General simulator description

The problem of interest is the transient multiphase flow in the multistage hydraulically fractured well during well startup and early production after fracturing. The system simulated in the numeric experiments presented below is schematically shown in Figure 1. It consists of the surface choke (flow restricting device), used to control the well, the well, conducting fluids’ flow from subsurface to wellhead, the fractures (increased hydraulic conductivity channels created in rock by high pressure pumping of fluid and filled with propping agent) and the reservoir (layer of porous and permeable rock carrying hydrocarbons). In general, the physical processes taking place in different parts of the system have quite distinct spatial and temporal scales. However, our specific interest to describe the early time behaviour of the well implies that characteristic time scales of the different parts of the system overlap. For example, processes of the fractures and well cleanup are mutually dependent. The degree of displacement of fracturing fluid by reservoir hydrocarbons in the well governs hydrostatic head and therefore the bottomhole pressure. The difference between the bottomhole and reservoir pressures in turn defines production rates from the reservoir and the speed of the displacement.

The problem of the coupled simulation in our case is solved by decomposition into several domains. The dimensions of the different domains vary from a 0D choke model, which is treated as pointwise object and described by the set of algebraic relations, 1D cross-section averaged wellbore flow and reservoir inflow models to a 2D averaged fracture flowback model.

Figure 1: Sketch of the simulated system

The principal parts of the numerical simulator are the choke module, wellbore flow module, hydraulic fracture module and reservoir module. Below, we will briefly describe physical models underlying each of these modules. Different approaches are used for coupling of different subdomains. In this work, we are focused on the problem of pressure-rate coupling between the well and fractures; accordingly the coupling between in choke-well and fractures-reservoir systems will be only outlined below.

In practice, the choke is used to create additional pressure drop at the wellhead depending on the flow rates of fluids produced. It allows restricting the flow and controling the well. Because of transient nature of the flow characterized by the gradual displacement of the fracturing fluid by liquid and gaseous hydrocarbons, the choke model should take into account presence of multiple phases and effects related to compressibility of gas. Here we adopt the model presented in alsafran2009predictions. It is formulated as a system of algebraic equations relating composition of the flow, flow rate, choke orifice diameter and pressures upstream and downstream of the choke. From the computational point of view, the choke module in our simulations calculates wellhead upstream choke pressure given the downstream choke pressure, choke opening, specified as boundary conditions, and flow rates through the choke.

The wellbore flow simulator used in this study has been previously described in spesivtsev2017study; spesivtsev2013modeling. The model consists of the 1D cross-section averaged mass and mixture momentum conservation equations. The slip between phases is described according to the drift-flux model shi2005b

. The wellbore module is governed by the wellhead pressure calculated from choke model and inflow from fractures represented as mass source terms. Given the wellhead pressure and inflows, it returns the surface flow rates and distribution of the bottomhole pressure along the well, in particular, values of pressure in points of connection of well with hydraulic fractures and the wellhead flow rates. For a given time moment the problem of pressure-rate coupling in choke-wellbore system may be considered as a problem of simultaneous solution of non-linear algebraic equations. The solution is calculated using a Newton-like algorithm for the fixed values of inflows provided by fractures’ models.

The numerical model of a hydraulic fracture chuprakov2020proppant is the key component of the simulator and incorporates both hydrodynamical and geomechanical aspects of the problem. Fluid flow relative to the proppant pack inside the fracture is governed by the Darcy law. The fluid phase continuity equation solved is two-dimensional and numerically resolves variation of flow rate inside the fracture. The model of fracture closure implemented is local in the sense that it assumes the elastic response of a proppant pack and its resulting width is governed by the local value of effective stress on the pack, mass of proppant per unit area of fracture and compression modulus of the proppant pack. The pack permeability and porosity are also governed barree2018generic by the local effective stress.

The reservoir matrix model provides the inflow of the flowback fluid into the fracture depending on the difference between pressure inside the fracture and the pore pressure. The model used in the present study is relatively simple and is based on the well-known Carter-type solution economides2000. We calculate the inflow to every discretization cell of the fracture according to the local pressure inside the fracture cell, far field reservoir pressure, porosity and compressibility of matrix and viscosity of the fluid. This approximation corresponds to the assumption of bilinear flow in the reservoir. The assumption is justified by the consideration of the early times after flowback start only. Because the Carter-type inflow is a linear function of pressure it may be implicitly approximated then solving for pressure distribution inside the fracture.

2.2 Coupling problem

The need to solve fixed-point nonlinear equations arises from the problem of coupling the rates between well and reservoir numerical models. In the well model, the length of the well is discretized into cells . Let be the mass rate of sources in -th cell. The index denotes a rate of a certain phase, e.g., oil, gas, fracturing fluid, proppant, etc. Depending on these rates and the state of the well at a given moment of time , the wellbore flow simulator returns the pressures in -th cell

. In vector form

The model of reservoir consists of several fractures, each attached to one certain cell of the wellbore. Given pressure in the adjacent cell, the fracture solves flow the problem inside of it and returns the resulting inflow rates.

At each timestep the goal is to match the rates , used as a boundary condition for the well simulator, to the resulting rates . Thus, we face a fixed-point problem:


As solving this equation is an intermediate procedure, done at each time step, the performance of the coupling method is crucial for the simulator’s overall performance. The method must be sufficiently robust to cope with acute changes of the flow regime, which occur in some cases of practical interest. The target function is acquired via numeric simulation, and the analytical form of its Jacobian is inaccessible, making Newton’s method unacceptable. Due to the high computational costs of evaluating , any schemes including line search or backtracking (press2007numerical, Chapter 9.7) are also inappropriate. Methods based on approximation of the Jacobian or its inverse with full rank matrices, such as Broyden methods, are not desirable due to the scale of the problem. Indeed, the dimension of variables vector is . Here is the number of phases of the flow (oil, water, gas, proppant), and , which is the number of cells that have a fracture attached to them, may be as large as , as in one of the cases studied below. This yields a Jacobian with components. Storage of such a large matrix is too demanding for our purposes. Additionally, for such methods, an initial guess for inverse Jacobian is required, and constructing such a guess for each timestep in the absence of an actual analytical Jacobian is a nontrivial problem itself. Thus we require an acceleration algorithm which

  • Stores only a fixed small amount of previous iterates that do not depend on the dimension of

  • Does not evaluate too often, e.g. for backtracking

  • Does not depend on more than one arbitrary parameter (like the relaxation factor in Picard method)

3 Algorithms

From here and further, denotes a sequence of points, generated by some iterative process, is the solution to the fixed-point problem; ; is the residual vector. If not explicitly said otherwise, represents scalar product, is the Euclidean norm in .

3.1 Fixed-point iteration with relaxation (Picard method)


0:   — initial guess,  — relaxation factor
1:  repeat
3:  until checkStopCondition() or

 This is the most straightforward solution. The user-defined parameter determines the length of the algorithm step. The intuition behind its choice is that using close to zero leads to better robustness of the algorithm at the cost of convergence speed. This intuition is supported by strict proof when is a contraction mapping, or, at least non-expansive mapping, see RyuBoydPrimer. In the former case, the convergence is linear with optimal choice of

and convergence speed estimate is

where is the contraction factor of . In the latter case, estimation for the residual is

with optimal value . The latter asymptotic is worse that the former, but the subset of on which the map is non-expansive might be larger than the subset where it is a contraction. Thus, usage of relaxation factor contributes to robustness but decreases speed of convergence.

3.2 Aitken relaxation

The first accelerated method we consider in this paper is Aitken relaxation. The idea of the method is to change the relaxation factor (and, thus, the size of the iteration step) based on the information from the previous iteration. The update method for the vector case was suggested in IronsTuck1969. We used it in a slightly different form, adapted from AitkenUR:  

0:   — initial guess,  — initial relaxation factor
1:  repeat
{Acquire new relaxation factor}
4:  until checkStopCondition() or

 In AitkenLin the authors consider the case where target function can be written as


is a Hermitian matrix. The condition for acceleration of the convergence is that all the eigenvalues of the matrix

have the same sign, where

is identity matrix.

Consider that the last value of relaxation parameter at timestep can be used as initial value at timestep . There is a possibility that such approach can help improve overall performance of the simulator, which is studied below.

3.3 Anderson’s acceleration


0:   — initial guess,  — relaxation factor,
1:  repeat
5:  until checkStopCondition() or

 The rationale under this algorithm is as follows. Consider to be linear: . Then


i.e., Picard step is done from  — the point that minimizes residual over the affine hull of previous iterates. In WalkerPengAnderson the authors prove general equivalence between Anderson’s acceleration and GMRES. In Fang2009TwoCO Anderson acceleration is considered as a member of a broader family of multisecant methods (in particular, Broyden’s family Type-II method).

The existing local convergence theory for general nonlinear case, to the authors’ knowledge, is not promising. Consider, for example, the sufficient conditions of convergence for Anderson acceleration with , proposed in TothKelleyConvergence. When is a contraction, one of the requirements is that contraction factor is sufficiently small, so that . Then the algorithm converges linearly with convergence factor not greater than :

In the same conditions, Picard method would converge with factor . If the value of is small, , so the guaranteed convergence rate for Anderson acceleration is worse. On the other hand, it is reported (see WalkerPengAnderson; LipnikovSvyatskiyVassilevski) that Anderson acceleration actually improves convergence of the coupled iterations. A comprehensive study of application of Anderson’s acceleration to the problem of variably saturated flow is presented in LottWalkerWoodwardYang. The computational experiments, conducted by the authors, demonstrate greater robustness of the algorithm compared to Picard accelerations.

Note that the choice of relaxation parameter is solely up to the researcher. In the works LipnikovSvyatskiyVassilevski; LottWalkerWoodwardYang the authors formulate the Anderson acceleration algorithm without introducing this parameter, in a way which is in our notation equivalent to substituting in formula 3.3.2. In WalkerPengAnderson; TothKelleyConvergence; AndersonOrig both theoretical and numerical findings are provided for the case . In Fang2009TwoCO successful computations with other than are reported, but, unfortunately, no rationale under the particular choice is provided. In the following numeric experiments we will try to address the problem of dependence of convergence speed of Anderson’s method on choice of .

3.4 Implementation details

In our implementation of Anderson’s algorithm, we use only information from two previous points (). The optimization problem (3.3.1) can then be easily solved analytically:


Computational formulas for both  Aitken’s and Anderson’s  algorithms  contain a term in the denominator (see (3.2.1), (3.4.1)). When two consecutive residuals are close to being equal, the mentioned denominator becomes small. In our practice, this resulted in unreasonably large iteration steps, causing oscillation of iterates, which slows down the convergence. Thus, in the numeric implementation, we apply the constraint

In case that constraint is not satisfied at the point , a Picard step with relaxation factor or for Aitken’s and Anderson’s methods, respectively is done . The stop condition in all the algorithms is


where is the required relative tolerance.

4 Numeric experiments

In the series of numerical experiments considered in this section we assess the performance of the implemented accelerated algorithms (Anderson’s and Aitken’s) compared to the widely used conventional fixed-point iterations. For brevity, they are abbreviated as TPA, two-point Anderson’s acceleration; ATK, Aitken’s relaxation and FPI, fixed-point iteration (Picard’s method).

The fixed point problem is formulated in terms of mass rates of each phase normalized by standard densities. Rates were chosen instead of pressure because, although rate error scales with pressure error, the scaling coefficient is dependent on the physical dimensions and properties of multiple system components, which is difficult to estimate beforehand. This scaling may also be highly variable during the simulation. At the same time, formulation in terms of rates ensures that mass balance conservation is directly controlled by the user-selected convergence tolerance of the iterative algorithm.

The algorithms were run with the same relative tolerance . We also impose a constraint on the maximum number of iterations allowed at a physical timestep. If the number of iterations exceeds some predefined value , which is the same for all three solvers, the whole simulation is considered diverged. Under this requirement, all solvers generated the same numerical solution within controlled tolerance consuming different numbers of iterations. The tests are run on the same predefined sequence of timesteps for each of the solvers.

The main goal of the study is to establish recommendations on the usage of these methods and provide guidelines for selection of suitable settings of the numerical algorithms. In particular, it would be useful to establish some rationale for the choice of relaxation factor ( in ATK) in each of the methods. For that, we perform extensive tests with different relaxation factors for each solver.

The performance of the algorithms is evaluated and benchmarked using two numerical problems involving hydrodynamic coupling between the horizontal well and multiple hydraulic fractures, inspired by oilfield applications. The first problem is cascading failure of hydraulic fractures due to aggressive flowback. The second problem is severe slugging during multiphase flow in horizontal well. The cascading failure problem involves rapid changes of flow rates caused by choke adjustment and fracture failure and very well highlights benefits of accelerated algorithm compared to conventional FPI. The slugging case evaluates behavior of iterative algorithms in a wide range of physical system parameters.

In the coupling problem, the time of the individual iteration is defined by the time of evaluation of the function. In our case, it is the total computational time of fracture and wellbore flow simulators. For all three algorithms, this time is much larger than all the additional operations (vector products and summations) needed to generate the next approximation. Thus, the total amount of iterations throughout the simulation is an adequate measure of the overall algorithm performance.

4.1 Cascading failure of fractures

One of the undesired phenomena that sometimes takes place in the oilfield industry is the damage of recently created hydraulic fracture conductivity due to aggressive flowback or excessively high flow rate during initial production. High production rates can cause mobilization of proppant in the near-wellbore zone of the hydraulic fractures experiencing the largest pressure gradients and fluid drag acting on the proppant pack. Washout of proppant from the near-wellbore region can result in drop of fracture conductivity and decrease overall fracture productivity. In the most critical cases, this can lead to formation of a pinch-out zone inside the fracture and disconnect it from the well. During flowback in multistage fractured wells, loss of productivity of one of the fractures can trigger the failure of others. In a system of multiple fractures producing into the same well with a single rate control by the surface choke, the failure of one producing fracture redistributes the full production rate between remaining ones. This happens because as the total flow rate decreases, due to the presence of choke, the bottomhole pressure decreases, and then, in turn, the drawdown and the rate from other fractures increases. As discussed in chertov2020, in case of a large number of producing fractures the total rate produced from the well after redistribution is almost unchanged. The redistribution leads to the increase of production rate from the remaining fractures, a subsequent increase of the drag forces on proppant packs inside them and brings remaining fractures closer to failure. The sequence of hydraulic fracture failures is schematically shown in Figure 2. Accordingly, the failure of hydraulic fractures in multistage wells due to aggressive flowback is a process with positive feedback. If multiple fractures are close to the critical proppant pack stability condition, the failure of the weakest fracture can trigger the cascading failure of many fractures. This phenomenon, named cascading failure of fractures, was introduced in CascadingFailure and studied numerically for different criteria of individual fracture failure in chertov2020.

Figure 2: Sequential failure of hydraulic fractures

Note, that in simulations here we consider an idealized model of fracture failure. It is assumed that fractures immediately and completely lose their conductivity if the fluid velocity in proppant pack exceeds some threshold value. As it was demonstrated in chertov2020, this idealized model of fracture failure tend to generate more steep and dramatic change in overall well productivity compared to explicit modeling of proppant wash-out from fracture and associated gradual degradation of their productivity. The preference on the abrupt fracture failure is justified by the goal of testing and evaluation of coupling algorithms in the most troublesome statement. In terms of the fixed-point problem, it means that that the solution from the previous time step used as the initial guess while solving (2.2.1) can be a poor approximation for the solution at the next time step and behaviour of the function will change significantly form one time step to another.

Figure 3 demonstrates the total production rate from the multistage horizontal well with hydraulic fractures calculated during numerical simulation of the cascading failure case. To give a flavour of the physical properties of the system, the well is composed of m long vertical section and m long horizontal section. There are hydraulic fractures connected to horizontal section with equal spacing of about m between them. Each fracture is connected to the well by a single cell. Once the critical fluid velocity in this cell is reached, the cell is set nonconductive and the entire fracture is disconnected from the well. During the first hours, the choke diameter is kept constant, and wellhead downstream choke pressure is gradually reduced about twice, from  bar to  bar. After that, is kept constant, and choke opening diameter is stepwise increased by every s. After each stepwise increase of the choke, the total flow rate is stepwise increased, as shown in Figure 3. Critical velocity to disconnect each fracture is linearly varied within certain range. Initially, flow rate increase after each choke increase is insufficient to fail any fractures. Then, as the number of producing fractures is still large, each choke increase may fail a few fractures and may also potentially trigger cascading failure of a few more fractures due to rate redistribution. At this stage, the total flow rate is nearly flat between stepwise choke openings. As the number of producing fractures goes down, redistributed increase of the flow rate on remaining fractures goes up, the number of sequentially triggered rate redistributions and the number of fractures dropped after each redistribution goes up. As a result, each stepwise increase of the flow rate is now followed by a decreasing slope of flow rate. In the end, the last choke opening initiates the last cascade of fracture failures, when there is no combination of fractures that can support the flow rate, and all remaining fractures fail in a sequence.

Figure 3: Surface rate for the cascading failure test case

Because FPI is the simplest and most widely applied method, we will estimate the relative acceleration provided by TPA and ATK as , or , where is the total number of iterations for each algorithm. Table 1 gives the summary of the benchmarking on the cascading failure case at different values of relaxation factor parameter. is the number of iterations consumed in the best case, which occurred at relaxation parameter . Relative acceleration in Table 1 is presented for the best run of each algorithm.

0.001 0.01 0.1 0.3 0.5 0.7 1.0
TPA 56,088 18,810 9,672 13,971 2,274 11,910 1,095 1.000 1,095 600
ATK 59,943 59,920 65,892 67,302 75,006 75,074 73,556 0.010 59,920 11
FPI 656,888 Did not converge 0.001 656,888
Table 1: Results of numerical experiments for Cascading Failure case

For all the solvers the required relative tolerance and maximal iterations per timestep . In these conditions, the FPI solver managed to converge only with rather small relaxation factor . The total amount of iterations with FPI was over .  The use of Anderson’s method with allowed acquiring drastic acceleration of around times, though even for the worst choice of the speedup is around . The findings also demonstrate a complicated nonmonotonous dependence of the performance of the TPA solver on the relaxation factor. The general trend that applies for close both to and to is that TPA converges faster for larger relaxation factors. On the other hand, some peculiarities occur in the middle of the range. For example, we can see that for the solver works more than  times faster than for relatively close values and . This suggests that some tuning of may be needed to achieve maximal performance at any given case.

For the ATK solver, acceleration is around times for the best and for the worst. Though not as dramatic as for TPA solver, this also makes the total computational time acceptable for practical studies. It is also interesting that the total amount of iterations is weakly dependent on initial values of , unlike Anderson’s solver. The number of iterations for different relaxation factors varies within approximately of average. This is because in our implementation of Aitken solver, we take on the next timestep equal to  — the last value on previous timestep. With more iterations, less and less depends on the initial choice of . For practical use, that means that ATK provides good acceleration over FPI for all initial values of relaxation parameters and fine tuning of initial is not necessary.

4.2 Severe slugging

Figure 4: Stages of severe slugging

Slugging flow is a special flow regime possible in pipelines and wellbores transporting gas and liquid simultaneously. Slugging flow is characterized by unstable system behavior with significant variations of pressure and flow rates in space and time. These variations can be observed even for steady-state boundary conditions such as constant pressure at the pipeline segment outlet and constant gas and liquid rates at the pipeline inlet. In hilly terrain pipelines, the undulating trajectory may lead to the formation of terrain-induced slugs (for example, taitel1990; dehenau1995; spesivtsev2017study). If the configuration consists of a near-horizontal pipeline or a well connected with a vertical riser (Figure 4), the transient behavior observed in such systems is referred to as severe slugging (for example, fabreetal1990; balino2010). The physical mechanisms leading to the formation of the slug flow are the same in both cases. The effect is caused by the competition between gravity-dominated liquid flow and compressibility of the gas phase. The liquid might accumulate in the lower parts of the system and subsequently block the free gas flow (see Figure 4). The compressibility of gas and continuous inflow to the pipe allows the pressure to build up gradually in the trapped gas volume (Figure 44). Once the gas pressure reaches some threshold, rapid outflow from the system occurs (Figure 4) and the process is repeated (Figure 4). The terrain-induced slugging is generally treated as an undesirable phenomenon, as it may result in the separator overflow and damage of the surface equipment.

Figure 5: Toe-up, horizontal and toe-down wells with bottomhole pressure behavior for different choke diameters

For the parametric study of the process we have designed a simplified model of the well in which slugging can occur. It consists of a vertical and an inclined sections, with lengths km and km respectively. The inclination of the horizontal section is characterized by the angle ; represents the toe-down well and  — toe-up. fractures are located on the inclined section with equal spacing. This simple setup geometry is sufficient to develop a slugging flow regime.

From the computational point of view, this problem was chosen because it can provide a set of test cases which can be defined by a small amount of varying parameters (in particular, choke size and the inclination angle ), but provide qualitatively different flow regimes (oscillatory when slugging occurs, relaxation to steady-state otherwise).

From our previous experience, we anticipated that there were two main factors that may inhibit the convergence in this parametric study. First, we expected fast changes in the solution when the flow regime is oscillatory. When such changes occur, the initial approximation taken from the previous step may be far from the solution on the next step, so the convergence might slow down. Second, we expected that the behavior of the target function can vary with physical parameters. In particular, when the choke diameter is small, the dependence of pressure on rate becomes steeper. Although the analytical form of is not available, this suggests that convergence at small chokes can be hindered.

The goal of the numerical experiments with slugging is to asses the performance of the coupling algorithms on a set of cases with different flow regimes and properties of target function. The dependence of the performance on the choice of relaxation factor is also studied to provide further insights on its role in convergence.

The set of cases is generated by variation of the inclination and the choke opening . The choke opening is kept constant throughout the entire simulation that covers hours. To investigate the dependence of the solvers’ performance on the time step, we performed a parametric study for two time steps: s and s. The relative tolerance was and maximal number of iterations per timestep .

0.1 0.5 1.0 0.1 0.5 1.0 0.1 0.5 1.0
15.0 -5.0 60,123 10,680 10,353 4,465 4,210 4,793 31,373 5,998 4,924
0.0 60,395 9,685 9,393 3,700 4,707 4,246 19,763 4,737 3,741
5.0 60,739 9,990 8,918 3,545 4,496 4,003 29,777 6,075 4,436
40.0 -5.0 62,742 9,746 3,550 3,627 3,680 3,624 30,075 5,979 3,682
0.0 62,591 9,727 3,453 3,584 3,608 3,601 29,749 5,988 3,628
5.0 62,697 9,760 3,455 3,507 3,527 3,531 29,040 5,738 3,611
60.0 -5.0 65,458 10,614 4,352 4,749 4,430 4,668 28,009 6,576 4,392
0.0 63,312 9,815 3,608 3,670 3,722 3,679 27,389 6,002 3,704
5.0 62,300 9,680 3,514 3,620 3,587 3,612 27,494 5,819 3,630
Avg. 62,261 9,966 5,621 3,829 3,996 3,973 28,074 5,879 3,972
Table 2: Results of parametric study for time step s
Figure 6: Maximal acceleration for s

A representative slice of results of the test runs for time step s is compiled in Table 2. The maximum acceleration for TPA and ATK relative to FPI and the map of slugging/non-slugging flow regimes in coordinates is depicted in Figure 6. For the cases with small choke diameter , the slugging did not occur for any inclination angles. For , it appeared only for toe-up and horizontal wells (). For the cases with wide open choke (), oscillations were present at all inclinations. These observations are in qualitative agreement with generally accepted conclusion that choking effectively eliminates or reduces severity of slugging jansen1996elimination.

One notable finding is that the total amount of iterations does not vary from non-oscillatory to oscillatory flow regime with same choke value as much as it varies from smaller choke values to bigger ones. Roughly speaking, all cases with small chokes are harder for all three iterative solvers than all cases with large choke. For a set of cases with fixed choke, oscillatory cases, if present, are only slightly harder than non-oscillatory ones. For the case with the smallest choke and steepest toe-up well all three solvers diverged (marked N/A on Figure 6). The maximum acceleration ( for ATK and for TPA) is achieved in cases with the smallest choke value. For larger chokes, this acceleration is , so it is fair to say that acceleration is achieved only for a limited set of modeling conditions. Based on these data, we can assume that the factor of steeper change of target function is more detrimental for convergence than the factor of generating a bad initial approximation due to oscillations. The former factor only comes to play for the case . All algorithms require more iterations on this case, and accelerated solvers start to demonstrate minor acceleration of about compared to almost no acceleration in nearby cases with large choke.

0.1 0.5 1.0 0.1 0.5 1.0 0.1 0.5 1.0
15.0 -5.0 405,675 62,659 14,179 12,872 12,950 12,858 129,817 26,687 13,208
0.0 410,835 63,855 14,188 13,173 13,243 13,284 169,768 25,840 13,164
5.0 410,930 63,750 14,398 13,217 13,210 13,226 125,375 27,234 13,199
40.0 -5.0 431,409 66,709 15,692 15,830 15,791 15,820 146,306 30,667 16,132
0.0 429,591 66,685 15,648 15,501 15,539 15,535 155,729 30,924 16,137
5.0 429,429 66,969 15,739 15,568 15,637 15,576 151,272 33,076 16,123
60.0 -5.0 462,578 71,354 16,489 16,787 16,730 16,734 136,651 31,470 16,421
0.0 455,892 70,083 16,199 16,488 16,492 16,471 152,419 31,293 16,420
5.0 456,655 70,201 16,180 16,481 16,438 16,421 155,762 32,589 16,270
Avg. 432,554 66,918 15,412 15,101 15,114 15,102 147,011 29,975 15,230
Table 3: Results of parametric study for time step s
Figure 7: Maximal acceleration for s.

The results of computations for the same set of test cases, but with reduced time step s are presented in Table 3 and Figure 7. We can see that the difference in total number of iterations between methods diminishes. This is because as the initial guess at step is the converged solution from previous time step . As becomes smaller, the consequent solutions become closer to each other; the initial approximation becomes better, and all solvers need less steps to converge. Note that, despite the amount of coupling steps per timestep decreased, the number of timesteps and therefore the total amount of coupling iterations for the whole case increased. The maximal acceleration and for ATK and TPA, respectively, was observed on the case with .

15.0 40.0 60.0
-5.0 0.0 5.0 -5.0 0.0 5.0 -5.0 0.0 5.0
0.10 31,373 19,763 29,777 30,075 29,749 29,040 28,009 27,389 27,494
0.25 13,950 9,101 12,351 12,581 12,570 12,616 12,756 12,096 11,816
0.50 5,998 4,737 6,075 5,979 5,988 5,738 6,576 6,002 5,819
0.75 4,357 3,848 3,620 4,007 3,983 3,972 5,332 4,033 3,992
0.95 4,296 3,793 3,763 3,697 3,636 3,605 4,716 3,690 3,596
1.00 4,924 3,741 4,436 3,682 3,628 3,611 4,392 3,704 3,630
Table 4: TPA performance for different values of ,

Table 4 represents the data for the performance of TPA method on the same set of cases for , but at different values of relaxation parameter . One can clearly see that the amount of iterations is decreasing with increasing . In most cases, either , or , but the total amount of iterations for is within of the minimal result. On the other hand, some deviations from that trend are present (underlined in Table 4) for . is for , is for and result for is and worse than optimal, respectively. From these findings we can suggest that is a good default choice for most of the cases, but one could opt for additional fine tuning to achieve the fastest possible convergence rate under specific combination of model parameters.

To analyze the performance of the ATK solver and its dependence on the initial relaxation factor , let us return to the Table 2. A notable finding is that the total number of iterations depends on this parameter much less than it depends on for TPA and FPI solver. To provide a quantitative estimate, we consider the relative difference between , the number of iterations averaged over in given case, and , minimal iterations required in the case. is considered a rough estimate of the expected amount of iterations when is chosen completely arbitrarily. The maximal relative difference of around between and with and inclinations . For three other cases it is less than , and for five other it is less than . The same relative difference is around for TPA, and for FPI method, so usage of a completely arbitrary relaxation factor for these methods without any tuning leads to drastic slowdown of convergence. Overall, in the slugging test case, ATK provides the same or slightly higher performance than TPA with an additional advantage that it does not require fine tuning of relaxation parameter.

Please note once again that all the numerical experiments for Anderson’s acceleration method are done only for a particular case of the algorithm where the next iteration depends only on two previous iterates. Recommendations on the choice of and estimates of performance should be carefully if at all extrapolated on the implementations of Anderson’s algorithm using more than two previous iterates. It was chosen to use implementation with only two iterates because it requires only the most basic vector operations and the least additional memory. Nevertheless, this implementation provided more than sufficient improvement of performance for our goals. Usage of more previous iterates can further improve the performance of the solver, but then the need to efficiently solve the minimization problem 3.3.1 and monitor the conditioning number of the related matrices arises. More details on the multipoint Anderson’s method can be found in WalkerPengAnderson; LottWalkerWoodwardYang.

5 Conclusion

Results of our study indicate that accelerated Anderson’s and Aitken’s algorithms are handy tools for the solution of fixed-point problems arising, for example, in coupling of systems in multiphysics modelling. The widely used and simple fixed-point iterations algorithm — the Picard method — has several downsides. One of them is slow convergence. The other one is dependence on a manually tuned relaxation factor. In some dynamic evolution-type problems, there may be no single value of relaxation factor that ensures convergence for the entire duration of the physical process. Specifics of the coupling problem, when it is necessary to iteratively couple existing numerical modules, often impose specific constraints on the solution. Typical constraints are unavailability of analytical Jacobian and computationally expensive target function. This is the case with the well-fracture coupling, which is one of the important problems in e.g. oilfield application, and which was used here for demonstration. To cope with these constraints, Anderson’s and Aitken’s methods are suggested. Neither of the methods relies on the usage of the Jacobian, its inverse, or their full-rank approximations and has negligible memory and CPU requirements compared to target function calculation. The target function is evaluated only once per iteration. Although it’s not guaranteed mathematically that the global convergence properties for these methods would be better than for the Picard’s method, in practical applications they demonstrate noticeably better performance. In this work, we benchmarked and reconfirmed performance of our implementation of the two-point Anderson’s (TPA) and Aitken’s (ATK) algorithms using practical examples of well-fracture coupling inspired by oilfield applications. We also investigated how performance of these methods depends on algorithm settings and behavior of the simulated problem and generated recommendations for their efficient usage based on numerical experiments.

In the numerical experiments we considered cases with relatively smooth solutions and cases with intermittent, step-like changes in solution. ATK tends to be equal or slightly better than TPA in cases with relatively smooth variations of pressures and rates. At the same, TPA can be superior compared to ATK, by an order of magnitude, in cases with abrupt changes of flow rates, but may require some tuning of relaxation factor.

The data obtained from the numerical experiments generally supports the notion presented in the literature that is a good choice of relaxation factor in Anderson’s algorithm. In the majority of cases, the result was either optimal or sufficiently close to the optimum. On the other hand, there are some notable deviations from this trend. In these cases, fine-tuning of the relaxation parameter may result in faster convergence. The analysis above has shown weak dependence of Aitken’s method performance on the initial choice of relaxation factor.

From these findings, we suggest the following practical strategy of usage of the accelerated solvers. Aitken’s algorithm can be used as a default solver for most of the problems. With an arbitrary , it provides significant acceleration in most of cases. In the worst cases we have seen, it is not slower than the Picard method. For harder problems, e.g., with fast intermittent changes of rate, or small choke, the solution may be further accelerated by application of TPA, which may require fine-tuning of . In that case, a good starting point is a value .


The authors would like to thank Schlumberger for their support and permission to publish this paper. We would also like to acknowledge A. Chaplygin and P. Spesivtsev for various contributions that helped to complete this paper.