Accelerators
In nonlinear finite element analysis, iterative solution schemes such as the Newton-Raphson method may converge slowly or fail to converge altogether when the system exhibits strong nonlinearities. Solution accelerators improve the convergence behaviour of the iterative process without requiring a change to the underlying Jacobian. numgeo provides two fundamentally different acceleration strategies:
- Line search, which scales the magnitude of the iterative correction along its original direction to minimise the energy potential, and
- Anderson acceleration, which combines information from successive iterations to modify both the magnitude and the direction of the update.
Both approaches are selected via the type parameter on the *Accelerator keyword (previsously *Line search).
Line search
The rationale behind Line Search is that the direction found by the Newton-Raphson method is often a good direction, but the step size (magnitude of the solution) is not. Furthermore, it is cheaper to compute the residual for several points along \(\boldsymbol{c}\) rather than form and factor a new system Jacobian.
The Line Search algorithm uses a prediction of the iterative solution increment \(\boldsymbol{c}^i\) as obtained by the Newton-Raphson algorithm and scales this vector by a value to minimize the energy potential1. \(\Pi = \boldsymbol{r}^T \boldsymbol{c}^{(i)}\). Note that only the magnitude is scaled, the direction of the prediction remains unchanged. While the local minimum of the energy potential represents the equilibrium, the minimum in the line search direction can be regarded as the best solution in the predicted direction. The scaled iterative increment reads:
where \(\lambda\) is a scalar scaling variable. For \(\lambda > 1.0\) an extrapolation is performed. The scaling factor is bound to \(\lambda_{min} \leq \lambda \leq \lambda_{max}\), per default numgeo uses \(\lambda_{min}=0.25\) and \(\lambda_{max}=1.0\) (no extrapolation). A minimum of \(\Pi\) in the line search direction requires that the derivative of \(\Pi\) to \(\eta\) must be zero:
The above can be interpreted as follows: at the minimum, the residual \(\boldsymbol{r}\) is orthogonal to the direction \(\boldsymbol{c}\). Equation (2) can be solved by iterative refinement of \(\lambda\). In numgeo, the purpose of the Line Search is to accelerate the Newton-Raphson method or to ''help'' finding convergence where none would be achieved otherwise.
Relaxation
The simplest form of step-size control applies a constant, user-defined scaling factor \(\lambda\) to every correction:
Setting \(\lambda < 1\) produces under-relaxation, which can stabilise iterations in strongly nonlinear problems at the cost of slower convergence.
Usage
Relaxation is selected with Type = Relaxation. It is the cheapest option and requires no additional residual evaluations. It is most useful when divergence is caused by overshoot rather than a poor search direction.
Linear Line Search
In its simplest form a linear relation between the potential at the beginning and the end of the present increment is assumed. The residual for \(\lambda=0\) is known from the previous increment, and the residual for \(\lambda=1\) is known from the present increment. Assuming a linear relation in between yields the value of \(\lambda\) without extra calculations:
where \(\boldsymbol{r}_0\) is the residual before the correction was applied and \(\boldsymbol{r}\) is the residual after a full correction \(\lambda = 1\). Evaluating \(\boldsymbol{r}\) at \(\lambda = 1\) requires one additional residual assembly, which makes this method more expensive per iteration than the other variants.
Usage
Linear Line Search is selected with Type = Linear. Because it requires one extra residual assembly per iteration, it should only be used when the cost of an additional assembly is small relative to the expected reduction in iteration count.
Multi-field simulations
In simulations involving multiple active physical fields, such as a consolidation analysis using coupled two-phase elements, the scaling factor \(\lambda\) is calculated separately for each active field.
Back-looking Line Search
The back-looking line search avoids the extra residual assembly of the linear method by re-using information that is already available from the previous iteration. At iteration \(k\), the residuals \(\boldsymbol{r}_{k-1}\) and \(\boldsymbol{r}_k\) are both known, as is the previous correction \(\boldsymbol{c}_{k-1}\). Assuming the same linear model as in the linear line search but applied to quantities from the previous step, the scaling factor is computed as:
This \(\lambda\) is mathematically the optimal scaling for the previous correction \(\boldsymbol{c}_{k-1}\). It is then applied to the current correction \(\boldsymbol{c}_k\) under the assumption that successive corrections are approximately parallel, which is a reasonable approximation in modified Newton-Raphson iterations where the Jacobian is held constant.
Usage
Back-looking Line Search is selected with Type = Back-Looking. It requires no additional residual evaluations and is therefore as cheap as relaxation, while adapting the scaling factor to the current state of the iteration.
Multi-field simulations
As with the linear line search, the scaling factor \(\lambda\) is calculated separately for each active degree-of-freedom field (displacements, pore water pressure, volumetric Jacobian \(J\)).
Limitations
The back-looking approach is most effective when successive corrections point in similar directions. Near the formation of a failure mechanism, where the plastic zone evolves between iterations, the direction of \(\boldsymbol{c}\) can change substantially and the transferred scaling factor becomes less reliable. In such situations, Anderson acceleration may be more effective.
Anderson acceleration
When the global Jacobian used in the Newton-Raphson iteration is not updated at every iteration or changes only mildly (for instance when an elastic tangent is used in place of the consistent elasto-plastic tangent) the iteration reduces to a fixed-point iteration of the form:
where \(\mathbf{K}_e\) is the (constant) elastic stiffness matrix and \(\boldsymbol{r}\) is the out-of-balance force vector. Such a scheme converges at best linearly with rate
where \(\mathbf{K}_\mathrm{ep}\) is the (unknown) elasto-plastic tangent. As plasticity spreads and \(\mathbf{K}_\mathrm{ep}\) softens, \(\rho\) approaches unity and convergence becomes arbitrarily slow. This is particularly problematic in strength reduction analyses (safety factor analyses), where the global failure mechanism causes \(\rho \to 1\) at the critical load level.
Anderson acceleration2 (also known as Anderson mixing) is a general technique for accelerating the convergence of fixed-point iterations. It achieves this by forming the next iterate not from a single correction, but from a linear combination of the most recent corrections, chosen to minimize the residual in a least-squares sense. Unlike line search, which only scales the magnitude of the correction, Anderson acceleration modifies both the magnitude and the direction of the update.
Algorithm (depth-1)
numgeo implements Anderson acceleration at depth \(m=1\), which requires only the current and one previous iteration. At iteration \(k \geq 1\) the following quantities are available:
| Symbol | Description |
|---|---|
| \(\boldsymbol{u}^{(k)}\) | current cumulative solution increment |
| \(\boldsymbol{c}^{(k)}\) | current correction from the linear solver |
| \(\boldsymbol{u}^{(k-1)}\) | cumulative solution increment from the previous iteration |
| \(\boldsymbol{c}^{(k-1)}\) | correction from the previous iteration |
Define the differences
The mixing parameter \(\alpha\) is determined by minimising \(\|\boldsymbol{c}^{(k)} - \alpha\,\Delta\boldsymbol{c}\|^2\), which yields:
The accelerated update then reads:
For \(\alpha = 0\) this reduces to the standard modified Newton-Raphson step. When the fixed-point map is affine (i.e. when the Jacobian is exactly constant), the least-squares minimisation is exact and depth-1 Anderson acceleration is equivalent to a secant update that projects the iterate to the root of the plane spanned by the last two residuals3.
First iteration
At iteration \(k=0\) of each increment, no previous history is available. numgeo performs a standard unaccelerated step and stores \(\boldsymbol{u}^{(0)}\) and \(\boldsymbol{c}^{(0)}\) for use at \(k=1\).
Safeguards
The mixing parameter \(\alpha\) is clamped to \(|\alpha| \leq \alpha_\mathrm{max}\), where \(\alpha_\mathrm{max}\) is controlled by the LAMBDA_MAX parameter on the *Acceleration keyword. This prevents over-extrapolation when the least-squares problem is poorly conditioned.
If \(\|\Delta\boldsymbol{c}\|^2\) falls below a numerical tolerance, indicating that the corrections have not changed between iterations, acceleration is skipped and a plain modified Newton-Raphson step is taken.
The stored history is cleared at the start of every new increment. In strength reduction analyses this means that the acceleration restarts for each new factor-of-safety level, which is necessary because the constitutive parameters (and therefore the fixed-point map) change between increments.
Not combinable with line search
Anderson acceleration and line search are mutually exclusive. The METHOD parameter on *LINE_SEARCH selects one or the other. Running both simultaneously is not supported, because Anderson acceleration assumes that the correction \(\boldsymbol{c}^{(k)}\) has not been rescaled by a line search factor.
Applicability
Anderson acceleration is most effective in problems where the Jacobian is held constant (or nearly constant) across iterations, such as modified Newton-Raphson with an elastic tangent. It is particularly well-suited for strength reduction (safety factor) analyses, where the iteration count near the critical load level can be reduced substantially. It can also be used in standard implicit analyses, though the benefit is smaller when the Jacobian already approximates the consistent tangent.
-
Reference Manual__
-
Strictly speaking ''energy potential'' is the correct terminology for the physical behavior, e.g. in case of plastic deformation. However, this poses no problem for the algorithmic implementation within an increment. ↩
-
Anderson, D.G. (1965). Iterative procedures for nonlinear integral equations. Journal of the ACM, 12(4), 547–560. ↩
-
Walker, H.F. and Ni, P. (2011). Anderson acceleration for fixed-point iterations. SIAM Journal on Numerical Analysis, 49(4), 1715–1735. ↩