Newton-Raphson solver (structsolve.newton_raphson)

Newton-Raphson solver for non-linear static analyses

The module-level constants control the step size and the divergence checks, and they are shared with the arc-length solvers of structsolve.arc_length.

The solver is used through Analysis.static() with Analysis.NL_method = 'NR'.

structsolve.newton_raphson._solver_NR(an, 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. See _check_convergence() for the convergence criteria. 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. The divergence and too-slow checks are only activated after NUM_ITER_BEFORE_DIVERGENCE_CHECK iterations, 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 an.final_load_factor, 1.0 by default. The analysis stops earlier, keeping the converged increments, when the load increment becomes smaller than an.minInc or when an.stop_criterion(an) returns True after a converged step, see stopping.

The callables of an that declare the optional keyword arguments inc, c or lbd receive the current load factor and configuration, see callbacks.

Parameters:

anAnalysis

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

silentbool, optional

A boolean to tell whether the log messages should be printed.

initialIncfloat or None, optional

Initial load increment, an.initialInc is used if 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().

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, see callbacks, 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 \|q\|\) is an additional reference of the relative 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 kT = kC + kG is rebuilt at every iteration, which gives quadratic convergence when kT is the exact Jacobian of fint. With an.modified_NR = True the tangent stiffness matrix 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, the linear constitutive stiffness matrix \(\boldsymbol{K}_{C0}\) is used instead of the tangent at the first iteration.

With an.line_search = True, the full Newton step (eta = 1) is tried first and it is accepted if it satisfies a sufficient-decrease test on the residual norm. Otherwise eta is reduced by backtracking, see _NR_iterations().

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.

With an.constraint_values, the values \(\boldsymbol{d}\) of the constraints \(\boldsymbol{C}^\top \boldsymbol{c} = \lambda \boldsymbol{d}\) of an.constraints, the linear solution of the first predictor is that of the bordered system with the values \(\boldsymbol{d}\), see static._factorize(), such that \(\boldsymbol{C}^\top \Delta \boldsymbol{c} = \Delta\lambda \boldsymbol{d}\) for every predictor, scaled from the previous converged increment, and the corrections keep it, their border being null. A prescribed displacement then needs no external load: the reactions of the constraints, part of the internal force vector, are the reference of the convergence criterion.

References

[1]

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

Convergence and divergence checks

structsolve.newton_raphson._check_convergence(an, R, fint, fext_total, fext_ref=None, q_ref=None, w=None, Kc_abs=None, full_output=False)

Check the convergence of the residual force vector

With an.convergence = 'crisfield', the default, the convergence is achieved when Crisfield’s \(\zeta \le\) an.relTOL, the residual and the reference force vectors being measured with the scaled norm of the weights w over the free degrees of freedom. With an.convergence = 'relative' it is achieved when either:

• ||R|| <= relTOL*max(||fext_total||, ||fint||, ||fext_ref||, ||q_ref||), using the Euclidean norm, if an.relTOL is not None

• max(|R|) <= absTOL, if an.absTOL is not None

an.absTOL is also accepted as an additional criterion with 'crisfield', and, with both criteria, the round-off floor of convergence.roundoff_floor() when Kc_abs is given. See convergence.check_convergence().

Parameters:

anAnalysis

Analysis object with the attributes convergence, relTOL and absTOL.

Rarray-like

Residual force vector.

fintarray-like

Internal force vector, of the same evaluation as R.

fext_totalarray-like

External force vector at the current load factor.

fext_refarray-like or None, optional

Additional reference force vector, e.g. the external force vector of the current increment of the arc-length solvers.

q_refarray-like or None, optional

Additional reference force vector, the load factor times \(-\partial \boldsymbol{R}/\partial \lambda\) when calc_fext accepts lbd, see callbacks.

warray-like or None, optional

Weights of Crisfield’s criterion, see convergence.crisfield_weights(), required with an.convergence = 'crisfield'.

Kc_absarray-like or None, optional

\(|K_{C0}| |c|\) of the round-off floor, see convergence.check_convergence(); None disables it.

full_outputbool, optional

Whether to also return roundoff.

Returns:

convergedbool

Whether the convergence criterion is satisfied.

Rnormfloat

Norm of the residual of the criterion, \(\|\boldsymbol{R}\|_{\boldsymbol{d}}\) with 'crisfield', the Euclidean norm with 'relative', used by the divergence checks and the line search.

Rrelfloat

Relative residual, Crisfield’s \(\zeta\) with 'crisfield', ||R||/max(||fext_total||, ||fint||, ||fext_ref||, ||q_ref||) with 'relative'.

roundoffbool

Only with full_output=True: whether the convergence is achieved by the round-off floor only.

Raises:

ValueError

When an.convergence is not 'crisfield' or 'relative', when an.relTOL and an.absTOL are both None, or when w is missing with 'crisfield'.

Notes

Crisfield’s criterion, see convergence, does not depend on the units of each degree of freedom, e.g. translations and rotations, whereas the Euclidean norm of the relative criterion adds the squares of forces and of moments. The reference fext_ref is useful when the load factor is close to zero, where fext_total and fint are small. The reference q_ref is needed under a prescribed displacement alone, where the external and the internal force vectors of the reduced problem both vanish at equilibrium and cannot be the reference.

structsolve.newton_raphson._check_divergence(an, Rnorm, Rnorms, iteration, silent=False, Rrel=None)

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 performed, and, after NUM_ITER_BEFORE_DIVERGENCE_CHECK iterations, when the residual diverges or when the convergence is too slow. With an.early_divergence = True a diverging step is also stopped earlier.

Parameters:

anAnalysis

Analysis object with the convergence criteria.

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

A boolean to tell whether the log messages should be printed.

Rrelfloat or None, optional

Relative residual of the current iteration, see _check_convergence(), used by the early divergence check.

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.

With an.early_divergence = True, a diverging step is stopped without waiting for NUM_ITER_BEFORE_DIVERGENCE_CHECK iterations, such that its load or arc-length increment is cut before the iterations are spent, when:

• the relative residual Rrel is larger than EARLY_DIVERGENCE_MAX_RREL, at any iteration;

• \(j >\) EARLY_DIVERGENCE_NUM_ITER and \(r_j >\) EARLY_DIVERGENCE_FACTOR \(\min(r_1, ..., r_{j-1})\).

These are the criteria with which the Riks analyses of a tow-steered cylinder of the DCSE website (2026) gave up the increments that would not converge beyond the limit point.

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

Newton-Raphson iterations at a fixed load factor

Corrects the predictor dc at the load factor total until the convergence criteria of _check_convergence() are 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, the early divergence check by an.early_divergence.

Parameters:

anAnalysis

Analysis object.

carray-like

Solution at the last converged state.

dcarray-like

Predictor of the solution increment.

totalfloat

Load factor.

fextarray-like

External force vector for a unit load factor.

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.

silentbool, optional

A boolean to tell whether the log messages should be printed.

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, see _check_convergence(). \(\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.

qarray-like 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.

warray-like or None, optional

Weights of Crisfield’s criterion, see convergence.crisfield_weights(), computed from calc_kC() when None and an.convergence = 'crisfield'.

K0_abssparse matrix or None, optional

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

projcallable or None, optional

Projection of the residual that removes the reactions of an.constraints, see _residual_projector(), computed when None and an.constraints is given.

Returns:

convergedbool

Whether the iterations converged.

dcarray-like

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

qarray-like 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) = \|R(c + \Delta c + \eta \, \delta c)\|^2\), 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, 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., 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

Constants

structsolve.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

structsolve.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

structsolve.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

structsolve.newton_raphson.LINE_SEARCH_ALPHA = 0.0001

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

structsolve.newton_raphson.INC_GROWTH_FACTOR = 1.1111

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

structsolve.newton_raphson.INC_CUT_FACTOR = 0.5

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

Inertia of the tangent stiffness matrix

structsolve.newton_raphson._converged_state(an, c, lbd, silent=False)

Factorization of the tangent stiffness matrix at a converged state, recording its inertia

Used with an.inertia = True. The tangent stiffness matrix at c and the load factor lbd is factorized with an.linear_solver, see static._factorize(), and the number of its negative eigenvalues, the sign of its determinant and the number of perturbed pivots are appended to an.num_negative_eigvals, an.det_signs and an.num_perturbed_pivots; for a singular matrix None, 0 and None.

Parameters:

anAnalysis

Analysis object.

carray-like

Converged solution.

lbdfloat

Converged load factor.

silentbool, optional

A boolean to tell whether the log messages should be printed.

Returns:

solve_kTcallable or None

The factorization, see static._factorize(), None when the matrix is singular.

structsolve.newton_raphson._eigvals_inertia(solve_kT, k, silent=False, s=None)

The number of negative real eigenvalues among the k of smallest magnitude of the SuperLU factorization solve_kT, set as its num_negative_eigvals, see Analysis.inertia_num_eigvals

Nothing is done when solve_kT already has a count, the exact one of PARDISO. The eigenvalues are computed by static._smallest_eigvals() with the LU factorization of solve_kT, on the degrees of freedom whose row or column is not null and on the space of the constraints, scaled by s, of all the degrees of freedom, when given: those of \(\boldsymbol{K}_T \boldsymbol{v} = \mu \boldsymbol{S}^{-2} \boldsymbol{v}\), \(\boldsymbol{S} = \mathrm{diag}(s)\). The analyses pass the scaling of _inertia_scaling(), \(\boldsymbol{S}^{-2} = |\mathrm{diag}(\boldsymbol{K}_0)|\), the diagonal of the linear stiffness matrix, fixed along the path: a kinetic criterion with a lumped stiffness-proportional mass, whose eigenvalues vary continuously along the path and cross zero where \(\boldsymbol{K}_T\) is singular, whatever the units of the unknowns. In the unknowns themselves, the eigenvalues of smallest magnitude are those of the unknowns of smallest stiffness, e.g. the high-order terms of a Ritz model, not those that cross zero: the cantilever cylinder of Schweizerhof and Ramm (1984) in panels (Sanders, 14 x 8 terms, 260 unknowns), beyond its bifurcation under follower pressure, had its negative eigenvalue 130th by magnitude without the scaling. The attributes set on solve_kT are:

• num_negative_eigvals: the number of real eigenvalues that are negative, None when ARPACK did not converge;

• num_complex_negative: the number of complex eigenvalues with a negative real part, among the k;

• inertia_num_eigvals: k;

• inertia_complete: False when the parity of the count differs from the sign of the determinant, i.e. an odd number of negative real eigenvalues lies beyond the k of smallest magnitude, with a warning.

Returns:

nnegint or None

solve_kT.num_negative_eigvals.