Convergence and troubleshooting of non-linear analyses#
The non-linear analyses of Analysis solve the equilibrium
\(\boldsymbol{R} = \boldsymbol{F}_{int}(\boldsymbol{c}) - \lambda \boldsymbol{F}_{ext} = \boldsymbol{0}\) in increments of the
load factor (Newton-Raphson) or of the arc length (Riks and Crisfield). The
attributes of Analysis listed below control the iterations of each
increment and the size of the increments.
Convergence criteria#
The criterion is selected with Analysis.convergence, see
convergence and newton_raphson._check_convergence():
'crisfield', the default: Crisfield’s criterion. Each component of a force-like vector is scaled with the diagonal of the linear stiffness matrixcalc_kC(), \(d_j = K_{{C0}_{jj}}\), the entries of the diagonal scaling vector \(\boldsymbol{d}\), over the free degrees of freedom, i.e. those whose row or column ofcalc_kC()is not null,\[\|v\|_{\boldsymbol{d}} = \sqrt{\sum_j v_j^2/|d_j|}, \qquad \zeta = \frac{\|\boldsymbol{R}\|_{\boldsymbol{d}}}{\max\left( \|\lambda \boldsymbol{F}_{ext}\|_{\boldsymbol{d}}, \|\boldsymbol{F}_{int}\|_{\boldsymbol{d}}, \|\Delta\lambda \boldsymbol{F}_{ext}\|_{\boldsymbol{d}}, \|\lambda \, \partial \boldsymbol{R}/\partial \lambda\|_{\boldsymbol{d}}\right)}\]and the iteration converged when \(\zeta \le\)
relTOL, default1e-9. \(\Delta\lambda \boldsymbol{F}_{ext}\) is a reference of the arc-length methods only, and \(-\lambda \, \partial \boldsymbol{R}/\partial \lambda\) of load vectors that depend on the load factor, seecallbacks. The criterion does not depend on the units of each degree of freedom, e.g. translations and rotations, and the residual at the constrained degrees of freedom, e.g. a support reaction, does not enter it.'relative': the Euclidean norm of \(\boldsymbol{R}\) is not larger thanrelTOLtimes the largest Euclidean norm of the same references, over all the degrees of freedom, the only criterion up to version 0.8.0, with the defaultrelTOL=1e-6then. It does not depend on the units of the model as a whole, but it adds the squares of forces and of moments.
With either criterion, absTOL, default None, also accepts an
iteration whose largest absolute value of \(\boldsymbol{R}\) is not larger than
absTOL, a force in the units of the model. None disables a
tolerance, and at least one must be given. The divergence and too-slow
checks and the line search use the norm of the selected criterion.
A tolerance below the round-off error of \(\boldsymbol{R}\) cannot be met, and every step would then fail, e.g. under large self-equilibrated axial forces near a null load. With either criterion, an iteration has therefore also converged when its residual is at the round-off floor,
with the Euclidean norm with 'relative', \(\varepsilon\) the machine
precision, \(|\boldsymbol{K}_{C0}|\) and \(|\boldsymbol{c}|\) the elementwise absolute values of
calc_kC() and of the solution, and \(C\) =
convergence.ROUNDOFF_FACTOR \(= 3\), see convergence.roundoff_floor();
\(\zeta\) may then be larger than relTOL, and the log messages mark these
iterations. ROUNDOFF_FACTOR = 0 disables the floor.
Increments and iterations#
initialInc: first load increment; with the arc-length methods, the first arc length, corresponding to a load factor increment ofinitialIncalong the initial tangent.maxIncandminInc: largest and smallest increment. The analysis stops when a step fails with an increment smaller thanminInc.maxNumIter: corrections per step before the step is repeated with a smaller increment;too_slow_TOLand the module constants ofnewton_raphsondetect diverging and too slow iterations.maxArcLength: cumulative arc length after which the arc-length methods stop, witharc_length.MAX_NUM_STEPSconverged steps at most.modified_NR,compute_every_n,kT_initial_stateandline_search: variants of the iterations, seeAnalysis.
Detecting a run that did not finish#
The analyses finish at the load factor exactly 1.0 unless they stop
earlier, and an.increments and an.cs keep only the converged
increments. A run that did not finish is detected with:
an.static(NLgeom=True)
if not an.increments or an.increments[-1] != 1.:
... # stopped at an.increments[-1], or no increment converged
The log messages of the solvers give the reason, e.g. WARNING: Minimum
step size ... achieved!; they are not printed with silent=True, see
Notation and sign conventions. A finished run is not necessarily a continuous path:
beyond a limit point, the load-controlled Newton-Raphson method may converge
to another branch of the equilibrium path, which a jump of the solution
between increments reveals. The arc-length methods trace the path through
the limit point.
Checking the tangent stiffness matrix#
The tangent stiffness matrix calc_kC + calc_kG must be the derivative
of calc_fint, \(\boldsymbol{K}_T = \partial \boldsymbol{F}_{int}/\partial \boldsymbol{c}\). An inconsistent
tangent stiffness matrix, e.g. a missing geometric stiffness matrix, does
not change the converged solution, which is controlled by the residual, but
the convergence becomes linear and may fail, such that the steps are cut
and the analysis may stop. Two checks:
the order of convergence of the iterations, quadratic with an exact tangent stiffness matrix and full Newton-Raphson, in the log messages of the iterations;
central finite differences of
calc_fintcompared withcalc_kC + calc_kGat a few configurations.
Examples#
The analytic problems of tests/analytic_problems.py, see
Non-linear static analysis with the Newton-Raphson method. A finite-difference check of the tangent
stiffness matrix:
def tangent_error(p, c, eps=1e-6):
"""Relative difference between the tangent stiffness matrix of the
callables and the central finite differences of ``calc_fint`` at ``c``"""
an = p.analysis()
kT = (an.calc_kC(c, NLgeom=True) + an.calc_kG(c, NLgeom=True)).toarray()
kT_fd = np.zeros_like(kT)
for j in range(c.size):
dc = np.zeros_like(c)
dc[j] = eps*max(1., abs(c[j]))
kT_fd[:, j] = (p.fint(c + dc) - p.fint(c - dc))/(2*dc[j])
return np.linalg.norm(kT - kT_fd)/np.linalg.norm(kT_fd)
def test_tangent_consistency_finite_differences():
"""The tangent stiffness matrix must be the derivative of the internal
force vector, which a finite-difference check verifies at any
configuration"""
c = np.array([0.1, -0.2, 0.3, 0.5])
assert tangent_error(hardening_springs(), c) < 1e-8
assert tangent_error(von_mises_truss(), c[:2]) < 1e-8
# without kG the tangent is inconsistent
p = hardening_springs()
p.kG = lambda u: np.zeros((4, 4))
assert tangent_error(p, c) > 1e-2
An inconsistent tangent stiffness matrix:
def test_inconsistent_tangent_converges_linearly():
"""An inconsistent tangent stiffness matrix converges linearly, with many
more iterations per step, and here not at all at the higher load
factors: the analysis stops before the load factor 1.0"""
p = hardening_springs()
p.analysis(initialInc=0.1)
p.an.static(NLgeom=True, silent=True)
check_increments_and_equilibrium(p)
assert max(p.calls_per_step()) <= 5
q = hardening_springs()
ref = q.solve_reference(lbd=0.1)
q.kG = lambda u: np.zeros((4, 4))
q.analysis(initialInc=0.1)
q.an.static(NLgeom=True, silent=True)
assert q.an.increments[-1] < 1.
assert max(q.calls_per_step()) > 20
# the first attempt, at the load factor 0.1, stops after maxNumIter
# corrections without reaching relTOL, the step is then cut
assert q.an.increments[0] < 0.1
errors = [np.linalg.norm(c - ref)/np.linalg.norm(ref)
for _, c in q.fint_calls[:q.an.maxNumIter]]
orders = convergence_orders(errors)
assert len(orders) >= 1
assert max(orders) < 1.2, orders
A tolerance that cannot be met, with the round-off floor disabled, and a
run stopped by maxArcLength:
def test_nr_unreachable_tolerance_stops_without_increments(monkeypatch):
"""A tolerance below the round-off cannot be met without the round-off
floor of the convergence criteria, disabled here: every step is cut
until the increment is smaller than minInc, and no increment is
stored"""
from structsolve import convergence
monkeypatch.setattr(convergence, 'ROUNDOFF_FACTOR', 0.)
p = hardening_springs()
p.analysis(relTOL=1e-20)
p.an.static(NLgeom=True, silent=True)
assert p.an.increments == []
assert p.an.cs == []
@pytest.mark.parametrize('method', METHODS)
def test_arc_length_maxArcLength(method):
p = von_mises_truss(frac=1.5)
p.analysis(NL_method=method, initialInc=0.05, maxArcLength=0.3)
p.an.static(NLgeom=True, silent=True)
check_equilibrium(p)
inc = np.asarray(p.an.increments)
assert 0 < inc[-1] < 1.
Load control beyond a limit point:
def test_nr_jumps_over_limit_point():
"""Beyond the limit load, load control converges to another branch of
the equilibrium path: the run finishes at the load factor 1.0, with a
jump of the solution that the arc-length methods do not have"""
p = von_mises_truss(frac=1.5)
p.analysis(initialInc=0.05)
p.an.static(NLgeom=True, silent=True)
check_increments_and_equilibrium(p)
w = np.array([c[0] for c in p.an.cs])
assert np.diff(w).max() > 1.
# the same final state as the arc-length methods, on the inverted
# configuration
assert w[-1] > 2.
q = von_mises_truss(frac=1.5)
q.analysis(NL_method='arc_length_riks', initialInc=0.05)
q.an.static(NLgeom=True, silent=True)
wq = np.array([c[0] for c in q.an.cs])
assert np.diff(wq).max() < 0.5
np.testing.assert_allclose(p.an.cs[-1], q.an.cs[-1], rtol=1e-5)
Both criteria give the same increments and solutions when the forces and the displacements of a model are scaled:
@pytest.mark.parametrize('method', METHODS + ['NR'])
def test_units_independence(method):
frac = 1.5 if method != 'NR' else 0.99
results = []
for force, disp in [(1., 1.), (1.e6, 1.), (1.e-4, 1.e-3)]:
p = scaled_problem(von_mises_truss(frac=frac), force, disp)
p.analysis(NL_method=method, initialInc=0.05)
p.an.static(NLgeom=True, silent=True)
check_equilibrium(p)
results.append((np.asarray(p.an.increments),
disp*np.asarray(p.an.cs), len(p.fint_calls)))
inc0, cs0, calls0 = results[0]
assert inc0[-1] == 1.0
for inc, cs, calls in results[1:]:
assert calls == calls0
np.testing.assert_allclose(inc, inc0, rtol=1e-6, atol=1e-9)
np.testing.assert_allclose(cs, cs0, rtol=1e-6, atol=1e-9)
Only Crisfield’s criterion is independent of the units of each degree of freedom, and only it ignores the residual of a support:
def rescaled_dof_problem(p, s):
"""Same problem with the degree of freedom 1 measured in other units,
``c1' = s c1``, e.g. a rotation: its force is ``F1' = F1/s`` and the
stiffness matrices ``K' = S^-1 K S^-1``, ``S = diag(1, s)``"""
Sinv = np.array([1., 1./s])
return Problem(Sinv*p.fext,
lambda c: Sinv*p.fint(Sinv*c),
Sinv[:, None]*p.kC*Sinv[None, :],
lambda c: Sinv[:, None]*p.kG(Sinv*c)*Sinv[None, :])
def test_crisfield_criterion_units_invariance():
r"""Measuring one degree of freedom in other units does not change
Crisfield's `\zeta`, hence neither the number of iterations nor the
increments of the Newton-Raphson method, whereas it changes the relative
criterion, whose Euclidean norm adds the squares of forces of different
units. The arc-length constraint, with the Euclidean norm of the
displacements, is not invariant, see :func:`.arc_length._solver_arc_length`"""
method, frac = 'NR', 0.99
scales = [1., 1.e3, 1.e-3]
results = {}
for convergence in ('crisfield', 'relative'):
for s in scales:
p = rescaled_dof_problem(von_mises_truss(frac=frac), s)
p.analysis(NL_method=method, initialInc=0.05,
convergence=convergence, relTOL=1e-9)
p.an.static(NLgeom=True, silent=True)
cs = np.asarray(p.an.cs)/np.array([1., s])
results[convergence, s] = (np.asarray(p.an.increments), cs,
len(p.fint_calls))
inc0, cs0, calls0 = results['crisfield', 1.]
assert inc0[-1] == 1.
for s in scales[1:]:
inc, cs, calls = results['crisfield', s]
assert calls == calls0
np.testing.assert_allclose(inc, inc0, rtol=1e-12)
np.testing.assert_allclose(cs, cs0, rtol=1e-9, atol=1e-12)
calls = [results['relative', s][2] for s in scales]
assert len(set(calls)) > 1, calls
@pytest.mark.parametrize('method', METHODS + ['NR'])
def test_crisfield_criterion_free_dofs_only(method):
r"""A support as a null row and column of the stiffness matrices, with
the internal force vector of all the degrees of freedom, including the
support reaction: Crisfield's criterion measures the free degrees of
freedom only and gives the solution of the model without the support
degree of freedom, whereas the reaction keeps the relative criterion
from converging"""
p = hardening_springs()
n = p.kC.shape[0]
# dof 0 is the support, connected to the dof 1 by the first spring
K = np.zeros((n + 1, n + 1))
K[1:, 1:] = p.kC
def fint(c):
f = np.zeros(n + 1)
f[1:] = p.fint(c[1:])
f[0] = -100.*c[1]
return f
def kG(c):
G = np.zeros((n + 1, n + 1))
G[1:, 1:] = p.kG(c[1:])
return G
q = Problem(np.concatenate(([0.], p.fext)), fint, K, kG)
q.analysis(NL_method=method, initialInc=0.1, minInc=1e-3)
q.an.static(NLgeom=True, silent=True)
assert q.an.increments[-1] == 1.
np.testing.assert_allclose(q.an.cs[-1][1:], p.solve_reference(),
rtol=1e-8)
assert q.an.cs[-1][0] == 0.
r = Problem(q.fext, fint, K, kG)
r.analysis(NL_method=method, initialInc=0.1, minInc=1e-3,
convergence='relative')
r.an.static(NLgeom=True, silent=True)
assert not r.an.increments or r.an.increments[-1] < 1.