1 Introduction
In this paper, we study the robustness and performance of the algorithms for solution of the fixedpoint problem in application to pressurerate coupling between computational models used in petroleum reservoir simulation. Anderson’s AndersonOrig algorithm and a generalization of Aitken’s aitken1927xxv algorithm for multidimensional fixedpoint 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 steadystate 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 wellreservoir 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 reimplementation 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 wellreservoir coupling and will not be further discussed here. In the present work we adopt the implicit coupling approach that iterates pressurerate 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 fixedpoint 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 nonlinear 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 wellreservoir 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 fluidstructure 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 fixedpoint 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 fixedpoint problem. In the numerical framework we consider, the fixedpoint 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 timesaving 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 inhouse 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 fixedpoint (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 wellfracture 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 crosssection averaged wellbore flow and reservoir inflow models to a 2D averaged fracture flowback model.
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 pressurerate coupling between the well and fractures; accordingly the coupling between in chokewell and fracturesreservoir 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 crosssection averaged mass and mixture momentum conservation equations. The slip between phases is described according to the driftflux 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 pressurerate coupling in chokewellbore system may be considered as a problem of simultaneous solution of nonlinear algebraic equations. The solution is calculated using a Newtonlike 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 twodimensional 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 wellknown Cartertype 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 Cartertype 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 fixedpoint 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 fixedpoint problem:
(2.2.1) 
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 fixedpoint problem; ; is the residual vector. If not explicitly said otherwise, represents scalar product, is the Euclidean norm in .
3.1 Fixedpoint iteration with relaxation (Picard method)
This is the most straightforward solution. The userdefined 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 nonexpansive 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 nonexpansive 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:
(3.2.1) 
In AitkenLin the authors consider the case where target function can be written as
where
is a Hermitian matrix. The condition for acceleration of the convergence is that all the eigenvalues of the matrix
have the same sign, whereis 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
(3.3.1) 
(3.3.2) 
The rationale under this algorithm is as follows. Consider to be linear: . Then
(3.3.3) 
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 TypeII 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:
(3.4.1) 
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
(3.4.2) 
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 fixedpoint iterations. For brevity, they are abbreviated as TPA, twopoint Anderson’s acceleration; ATK, Aitken’s relaxation and FPI, fixedpoint 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 userselected 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 nearwellbore zone of the hydraulic fractures experiencing the largest pressure gradients and fluid drag acting on the proppant pack. Washout of proppant from the nearwellbore 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 pinchout 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.
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 washout 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 fixedpoint 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.
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  
Alg.  
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  — 
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
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 steadystate 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 terraininduced slugs (for example, taitel1990; dehenau1995; spesivtsev2017study). If the configuration consists of a nearhorizontal 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 gravitydominated 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 4, 4). Once the gas pressure reaches some threshold, rapid outflow from the system occurs (Figure 4) and the process is repeated (Figure 4). The terraininduced slugging is generally treated as an undesirable phenomenon, as it may result in the separator overflow and damage of the surface equipment.
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 toedown well and — toeup. 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 steadystate 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 .
Alg.  FPI  ATK  TPA  

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 
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/nonslugging 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 toeup 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 nonoscillatory 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 nonoscillatory ones. For the case with the smallest choke and steepest toeup 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.
Alg.  FPI  ATK  TPA  

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 
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 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 fixedpoint problems arising, for example, in coupling of systems in multiphysics modelling. The widely used and simple fixedpoint 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 evolutiontype 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 wellfracture 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 fullrank 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 twopoint Anderson’s (TPA) and Aitken’s (ATK) algorithms using practical examples of wellfracture 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, steplike 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, finetuning 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 finetuning of . In that case, a good starting point is a value .
Acknowledgements
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.