Newton-Raphson solver (pyfe3d.solver.newton_raphson)

Load-controlled Newton-Raphson solver of the non-linear static analyses of Analysis, selected with Analysis.NL_method = 'NR'. The module-level constants control the step size and the divergence checks, and they are shared with the arc-length solvers of arc_length.

pyfe3d.solver.newton_raphson.DIVERGENCE_FACTOR = 1000.0

A step diverged when its residual norm becomes larger than this factor times the residual norm at the start of the step

pyfe3d.solver.newton_raphson.INC_CUT_FACTOR = 0.5

Factor used to reduce the load increment when a step fails to converge

pyfe3d.solver.newton_raphson.INC_GROWTH_FACTOR = 1.1111

Factor used to increase the load increment after a step that converged without being cut

pyfe3d.solver.newton_raphson.LINE_SEARCH_ALPHA = 0.0001

Constant of the sufficient-decrease (Armijo) test of the line search

pyfe3d.solver.newton_raphson.NUM_ITER_BEFORE_DIVERGENCE_CHECK = 5

Number of iterations of a step before the divergence and too-slow checks are activated. A non-monotonic residual in the first iterations is normal, especially in postbuckling

pyfe3d.solver.newton_raphson.TOO_SLOW_WINDOW = 5

The convergence is too slow when the smallest residual norm of the last TOO_SLOW_WINDOW iterations is not at least too_slow_TOL (relative) smaller than the smallest residual norm of the previous iterations

pyfe3d.solver.newton_raphson._NR_iterations(an, pb, c, dc, total, fext, w, kT0=None, silent=False, lbd_dependent=False, q=None, K0_abs=None)

Newton-Raphson iterations at a fixed load factor

Corrects the predictor dc at the load factor total until the convergence criterion of check_convergence() is satisfied, or until the iterations are stopped by _check_divergence() or by a singular tangent stiffness matrix. The tangent stiffness matrix is calc_kC + calc_kG at the current iterate, see callbacks.calc_kT(); its updates and the line search are controlled by an.modified_NR, an.compute_every_n, an.line_search and an.max_iter_line_search.

Parameters:

anAnalysis

Analysis object with the parameters of the solver.

pbobject

Object with the callables, an itself or a FreeDofs.

cnp.ndarray

Solution at the last converged state.

dcnp.ndarray

Predictor of the solution increment.

totalfloat

Load factor, dimensionless.

fextnp.ndarray

External force vector for a unit load factor.

wnp.ndarray or None

Weights of Crisfield’s criterion, see crisfield_weights(), None with the relative criterion.

kT0sparse_matrix or None, optional

Tangent stiffness matrix used in the first iteration when an.modified_NR=True. If None it is calculated at c + dc. Default None.

silentbool, optional

Whether the log messages are suppressed. Default False.

lbd_dependentbool, optional

Whether calc_fext accepts lbd, see callbacks; total*q, with \(\boldsymbol{q} = -\partial \boldsymbol{R}/\partial \lambda\), is then an additional reference of the convergence criterion. \(\boldsymbol{q}\) is evaluated together with the tangent stiffness matrix, at the same iterate, and the last one evaluated is the reference until the next update of the tangent. Default False.

qnp.ndarray or None, optional

The last \(\boldsymbol{q}\) evaluated before these iterations, the reference until the first update of the tangent stiffness matrix, used only when lbd_dependent=True. Default None.

K0_abssparse_matrix or None, optional

\(|\boldsymbol{K}_{C0}|\) of the round-off floor of the convergence criteria, see check_convergence(), computed from pb.calc_kC() when None. Default None.

Returns:

convergedbool

Whether the iterations converged.

dcnp.ndarray

Solution increment, such that c + dc is the solution.

qnp.ndarray or None

The last \(\boldsymbol{q}\) evaluated, None when lbd_dependent=False.

Notes

At each iteration, with the load factor \(\lambda\) = total fixed,

\[\boldsymbol{R} = \boldsymbol{F}_{int}(\boldsymbol{c} + \Delta \boldsymbol{c}, \lambda) - \lambda \boldsymbol{F}_{ext}, \qquad \boldsymbol{K}_T \delta \boldsymbol{c} = -\boldsymbol{R}, \qquad \Delta \boldsymbol{c} \leftarrow \Delta \boldsymbol{c} + \eta \delta \boldsymbol{c}\]

With full Newton-Raphson, \(\boldsymbol{K}_T\) is computed and factorized at every iteration at \(c + \Delta c\), and the convergence is quadratic when it is the exact Jacobian \(\partial \boldsymbol{R}/\partial \boldsymbol{c}\) ([1], Section 1.2.2, p. 9). With an.modified_NR = True, \(\boldsymbol{K}_T\) is computed at the first iteration, or given as kT0, and then every an.compute_every_n iterations, its LU factorization being reused in the other iterations.

\(\eta = 1\) without line search. With an.line_search = True, with \(\phi(\eta) = \|\boldsymbol{R}(\boldsymbol{c} + \Delta \boldsymbol{c} + \eta \, \delta \boldsymbol{c})\|^2\) in the norm of the convergence criterion, the scaled norm \(\|\cdot\|_{\boldsymbol{d}}\) of Crisfield’s criterion or the Euclidean norm, the step \(\eta\) is accepted when it satisfies the sufficient-decrease (Armijo) test

\[\phi(\eta) \le (1 - 2 \alpha \eta) \phi(0)\]

with \(\alpha\) = LINE_SEARCH_ALPHA, which assumes \(\phi'(0) = -2\phi(0)\), exact for a Newton direction in any norm whose weights do not change along the line, i.e. when \(\boldsymbol{K}_T\) is the exact Jacobian at \(c + \Delta c\); with an.modified_NR = True it is an approximation. The full step \(\eta = 1\) is tried first. Otherwise \(\eta\) is replaced by the minimizer \(\phi(0)/A\) of the quadratic interpolation of \(\phi(0)\), \(\phi'(0)\) and \(\phi(\eta)\), with \(A = (\phi(\eta) - \phi(0) + 2\phi(0)\eta)/\eta^2\), or by \(\eta/2\) when \(A\) is not finite or not positive, safeguarded in \([0.1 \eta, 0.5 \eta]\). After an.max_iter_line_search reductions the last \(\eta\) is accepted. The internal force vector of the accepted point is reused in the residual of the next iteration.

The line search is adapted from the backtracking line search of [2] (Section 6.3.2, Algorithm 6.3.5, pp. 126-129, and Algorithm A6.3.1, pp. 325-327), applied to nonlinear equations with the merit function \(f = \frac{1}{2}\|D_F F\|^2\) ([2], Section 6.5, Eq. (6.5.3), p. 147, and Algorithm NEFN, p. 303), whose directional derivative along the Newton step is \(-F^\top F = -2f\) (p. 148, with \(D_F\) the identity). With \(F = R\), whose Jacobian is \(\boldsymbol{K}_T\), and \(\phi = 2f\), \(D_F\) being the identity with the Euclidean norm or \(\mathrm{diag}(\sqrt{w_j})\) with Crisfield’s criterion, the test above is their condition (6.3.3) (p. 118) with their \(\alpha = 10^{-4}\) (p. 126), and the first reduction is their quadratic backtrack, Eq. (6.3.17) (p. 127). It differs from Algorithm A6.3.1 in that: the quadratic interpolation through the last trial \(\eta\) is used at every reduction, whereas [2] uses, from the second backtrack on, the cubic interpolation of \(f(0)\), \(f'(0)\) and the last two trial values, Eq. (6.3.18) (p. 128); the bound \(0.5 \eta\) is also imposed at the first reduction, where [2] imposes only the lower bound \(0.1\), the quadratic backtrack being then smaller than \(1/(2(1 - \alpha))\) (p. 127; Algorithm A6.3.1, steps 10.3c.1E.4 and 10.3c.4); and [2] has no maximum number of reductions, but limits the length of the step to maxstep and stops when the step length becomes smaller than a minimum, returning the current point.

References

[1]

Crisfield, M. A., Non-linear Finite Element Analysis of Solids and Structures, Volume 1: Essentials, John Wiley & Sons, Chichester, 1991.

[2] (1,2,3,4,5)

Dennis, J. E., Jr. and Schnabel, R. B., Numerical Methods for Unconstrained Optimization and Nonlinear Equations, SIAM, Philadelphia, 1996 (first published by Prentice-Hall, Englewood Cliffs, NJ, 1983), pages of the SIAM edition. https://doi.org/10.1137/1.9781611971200

pyfe3d.solver.newton_raphson._abs_matrix(kC0)

Elementwise absolute value \(|\boldsymbol{K}_{C0}|\) of the linear stiffness matrix kC0, sparse, for the round-off floor of the convergence criteria, see roundoff_floor()

pyfe3d.solver.newton_raphson._check_divergence(an, Rnorm, Rnorms, iteration, silent=False)

Check if the iterations of a step must be stopped without convergence

The iterations are stopped when the residual is not finite, when an.maxNumIter corrections were done, and, after NUM_ITER_BEFORE_DIVERGENCE_CHECK iterations, when the residual diverges or when the convergence is too slow.

Parameters:

anAnalysis

Analysis object with the attributes maxNumIter and too_slow_TOL.

Rnormfloat

Norm of the residual at the current iteration, the norm of the convergence criterion, see check_convergence().

Rnormslist

Norms of the residual at the previous iterations of the step; the current norm is appended to it.

iterationint

Current iteration, starting at 1.

silentbool, optional

Whether the log messages are suppressed. Default False.

Returns:

stopbool

Whether the iterations must be stopped.

Notes

With \(r_i\) the residual norm at the iteration \(i\) of the step, \(r_1\) being the residual of the predictor, the iterations are stopped at the iteration \(j\) when:

• \(r_j\) is not finite;

• \(j >\) an.maxNumIter, i.e. an.maxNumIter corrections were done;

• \(j >\) NUM_ITER_BEFORE_DIVERGENCE_CHECK and \(r_j >\) DIVERGENCE_FACTOR \(r_1\) (divergence);

• \(j >\) NUM_ITER_BEFORE_DIVERGENCE_CHECK, \(j > W\) and \(\min(r_{j-W+1}, ..., r_j) > (1 -\) an.too_slow_TOL \() \min(r_1, ..., r_{j-W})\), with \(W\) = TOO_SLOW_WINDOW (too slow convergence).

The divergence and too-slow checks are delayed because a non-monotonic residual in the first iterations is normal, especially in postbuckling. The checks only compare norms of the same step with each other, and with the scaled norm of Crisfield’s criterion they do not depend on the units of the degrees-of-freedom.

pyfe3d.solver.newton_raphson._log_iteration(an, iteration, R, zeta, silent, lbd=None, roundoff=False)

Message of one equilibrium iteration, with \(\zeta\) of Crisfield’s criterion or the relative residual of the relative criterion, marked when the iteration converged by the round-off floor only, see check_convergence()

pyfe3d.solver.newton_raphson._norm(v, w)

Norm of the convergence criterion: \(\|v\|_{\boldsymbol{d}}\) with the weights w of Crisfield’s criterion, the Euclidean norm when w is None

pyfe3d.solver.newton_raphson._solver_NR(an, pb, silent=False, initialInc=None)

Newton-Raphson solver with load control

Used by Analysis.static() when NL_method='NR'. The load factor goes from zero to one in increments; at each increment the equilibrium of the internal force vector with the external force vector scaled by the load factor is solved with the iterations of _NR_iterations(), starting from a predictor extrapolated from the previous converged increment, and the converged load factors and solutions are appended to an.increments and an.cs.

The variant of the iterations is selected by an.modified_NR, an.compute_every_n, an.kT_initial_state, an.line_search and an.max_iter_line_search, and the convergence criterion by an.convergence, see check_convergence(). A step that does not converge within an.maxNumIter iterations, that diverges or that converges too slowly is repeated with the load increment multiplied by INC_CUT_FACTOR, see _check_divergence(). After a step that converged without being cut, the load increment is multiplied by INC_GROWTH_FACTOR, limited by an.maxInc. The last step is adjusted such that the analysis finishes at a load factor of exactly 1.0. The analysis stops earlier, keeping the converged increments, when the load increment becomes smaller than an.minInc.

Parameters:

anAnalysis

Analysis object; the converged load factors and solutions are appended to an.increments and an.cs.

pbobject

Object with the callables, an itself or a FreeDofs.

silentbool, optional

Whether the log messages are suppressed. Default False.

initialIncfloat or None, optional

Initial load increment, dimensionless; an.initialInc is used if None. Default None.

Notes

The load factor \(\lambda\) goes from zero to one in increments, solving at each increment

\[\boldsymbol{R} = \boldsymbol{F}_{int}(\boldsymbol{c}) - \lambda \boldsymbol{F}_{ext} = \boldsymbol{0}\]

Each increment starts from a predictor that extrapolates the previous converged increment of the solution, scaled by the ratio between the new and the previous load increments. The first increment uses the linear solution as predictor: \(\Delta\lambda \boldsymbol{K}_{C0}^{-1} \boldsymbol{F}_{ext}\), with \(\boldsymbol{K}_{C0}\) the linear stiffness matrix calc_kC(), which also gives the scaling of Crisfield’s criterion, computed once for the whole analysis, see convergence.

The callables that accept the keyword argument inc receive the load factor \(\lambda\) of the step, see callbacks, which allows configuration-dependent loads such as follower pressures.

When calc_fext accepts the keyword argument lbd, e.g. for prescribed displacements imposed exactly, the load vector \(\boldsymbol{q} = -\partial \boldsymbol{R}/\partial \lambda\) at \(\boldsymbol{c} = \boldsymbol{0}\), \(\lambda = 0\) gives the predictor of the first increment, its linear solution \(\boldsymbol{K}_{C0}^{-1} \boldsymbol{q}\), instead of the linear solution of the reference vector, which is null under a prescribed displacement alone, and \(\lambda \boldsymbol{q}\) is an additional reference of the convergence criterion. \(\boldsymbol{q}\) is evaluated together with the tangent stiffness matrix, at the same iterate, i.e. calc_fext is called once per update of the tangent stiffness matrix, and the last \(\boldsymbol{q}\) evaluated is the reference of the next convergence checks, also in the next step.

By default full Newton-Raphson is used, i.e. the tangent stiffness matrix is rebuilt at every iteration, which gives quadratic convergence when it is the exact Jacobian of the internal force vector. With an.modified_NR = True it is rebuilt only at the first iteration of each step and every an.compute_every_n iterations, and its LU factorization is reused in the other iterations, at the cost of a linear convergence rate; at the first step with an.kT_initial_state = False, \(\boldsymbol{K}_{C0}\) is used instead of the tangent at the first iteration. The full and modified Newton-Raphson methods, combined with incremental predictors, are described in [1] (Sections 1.2.2 and 1.2.3, pp. 8-13), with the quadratic and the linear convergence rates.

The step that would end within max(an.minInc, 1e-12) of the load factor 1.0, or beyond it, is replaced by the step to exactly 1.0.

References

[1]

Crisfield, M. A., Non-linear Finite Element Analysis of Solids and Structures, Volume 1: Essentials, John Wiley & Sons, Chichester, 1991.

pyfe3d.solver.newton_raphson._weights(an, kC0)

Weights of Crisfield’s criterion from the linear stiffness matrix kC0, see crisfield_weights(), or None with an.convergence = 'relative'