Follower loads and unsymmetric matrices

Loads that depend on the configuration, such as follower (hydrostatic) pressures, are supported through optional keyword arguments of the callables of Analysis, see callbacks: calc_fint, calc_kC and calc_kG receive the load factor inc when they declare it, and calc_fext receives the current configuration c, returning the load vector \(-\partial \boldsymbol{R}/\partial \lambda\) used by the arc-length methods. Their load stiffness makes the tangent stiffness matrix unsymmetric in general, which the solvers support: solve() and the non-linear solvers use LU factorizations, lb() selects the solvers of unsymmetric matrices when needed, symmetric=None, and freq() uses general eigenvalue solvers by default.

Non-linear analyses

Hardening springs under a load \(\lambda(\boldsymbol{q}_0 - \boldsymbol{K}_{F_{ref}} \boldsymbol{c})\) with an unsymmetric load stiffness matrix \(\boldsymbol{K}_F = -\partial \boldsymbol{F}/\partial \boldsymbol{c} = \lambda \boldsymbol{K}_{F_{ref}}\): the callables receive the load factor and the configuration, and Newton-Raphson and the arc-length methods converge to the equilibrium path. The Riks method converges quadratically only with the load vector of the current configuration:

class FollowerProblem(object):
    def __init__(self, K, a, q0, Q, with_c=True):
        self.K = K
        self.a = a
        self.q0 = q0
        self.Q = Q
        self.with_c = with_c
        self.an = None
        # (inc, ||R||/||lbd*q||, number of converged increments) of each
        # calc_fint call
        self.calls = []

    def fint_el(self, c):
        return self.K @ c + self.a*c**3

    def residual(self, c, lbd):
        return self.fint_el(c) - lbd*(self.q0 + self.Q @ c)

    def calc_fext(self, inc=1., silent=True, c=None):
        if c is None:
            return inc*self.q0
        return inc*(self.q0 + self.Q @ c)

    def calc_fint(self, c, inc=1., silent=True):
        R = self.residual(c, inc)
        ref = max(np.linalg.norm(inc*(self.q0 + self.Q @ c)), 1e-300)
        step = len(self.an.increments) if self.an is not None else 0
        self.calls.append((inc, np.linalg.norm(R)/ref, step))
        return self.fint_el(c) - inc*(self.Q @ c)

    def calc_kC(self, c=None, NLgeom=False, inc=1., silent=True):
        if not NLgeom:
            return csr_matrix(self.K)
        return csr_matrix(self.K + np.diag(3*self.a*c**2) - inc*self.Q)

    def calc_kG(self, c=None, NLgeom=False, silent=True):
        return csr_matrix(self.K.shape)

    def analysis(self, **kwargs):
        if self.with_c:
            calc_fext = self.calc_fext
        else:
            calc_fext = lambda inc=1., silent=True: self.calc_fext(inc)
        an = Analysis(calc_fext, self.calc_fint, self.calc_kC, self.calc_kG)
        for k, v in kwargs.items():
            setattr(an, k, v)
        self.an = an
        return an

def test_riks_quadratic_only_with_configuration_load_vector():
    r"""The bordered system of Riks is the Newton-Raphson method of the
    augmented system only when -dR/dlbd = q0 + Q c, i.e. with ``calc_fext``
    accepting ``c``; with the constant reference vector the corrections
    converge linearly"""
    worst = {}
    for with_c in (True, False):
        pb = make_problem(with_c=with_c)
        an = pb.analysis(NL_method='arc_length_riks', initialInc=0.5,
                         relTOL=1e-13, maxNumIter=60)
        an.static(NLgeom=True, silent=True)
        assert np.isclose(an.increments[-1], 1.)
        # residuals of the iterations of the first step, whose first value
        # is the one of the tangent predictor
        first = [e for inc, e, step in pb.calls if step == 0]
        worst[with_c] = (min(orders(first)) if orders(first) else None,
                         first)
    p_c, errors_c = worst[True]
    assert p_c is not None and p_c > 1.7, errors_c
    p_const, errors_const = worst[False]
    assert len(errors_const) > len(errors_c), (errors_c, errors_const)

Inertia of unsymmetric tangent stiffness matrices

The symmetric indefinite factorization of PARDISO counts the negative eigenvalues of a symmetric tangent stiffness matrix (Analysis.inertia with linear_solver='pardiso'); an unsymmetric one is factorized by SuperLU, whose LU factorization gives only the sign of the determinant, the parity of the number of negative real eigenvalues. With Analysis.inertia_num_eigvals = k, the number of negative real eigenvalues among the \(k\) of smallest magnitude is computed at every converged increment by shift-invert at \(\sigma = 0\) with ARPACK, with the LU factorization of the converged state, and recorded in Analysis.num_negative_eigvals, such that max_inertia_jump applies, see newton_raphson._eigvals_inertia().

It detects the divergence points, where a real eigenvalue crosses zero, the static critical points of the path: limit points and bifurcations. It does not detect flutter: a complex pair \(a \pm i b\) whose real part becomes negative is not a static critical point, the tangent staying regular, and it is not counted (the number of such pairs among the \(k\) is logged), the kinetic criterion of freq() being needed, see below. Negative real eigenvalues beyond the \(k\) of smallest magnitude are missed; an odd number of them is reported, since the parity of the count then differs from the sign of the determinant.

A problem whose tangent is \(\mathrm{diag}(k_0, \boldsymbol{S} (\boldsymbol{D}_0 + \lambda \boldsymbol{D}_1) \boldsymbol{S}^{-1})\), with two real eigenvalues crossing zero at \(\lambda = 0.27\) and \(0.63\) and a complex pair at \(0.47\): the count goes 0, 1, 2, the sign of the determinant +1, -1, +1, and the count agrees with numpy.linalg.eigvals() at every increment:

@pytest.mark.parametrize('linear_solver',
                         ['splu'] + (['pardiso'] if has_pardiso() else []))
@pytest.mark.parametrize('method', METHODS)
def test_count_along_path(method, linear_solver):
    """Two real eigenvalues cross zero, 0 -> 1 -> 2, the sign of the
    determinant returning to +1, and a complex pair crosses the imaginary
    axis, which is not counted"""
    p = LinearUnsymmetric([0.27, 0.63], complex_crossings=[0.47])
    an = p.analysis(NL_method=method, inertia=True, inertia_num_eigvals=6,
                    linear_solver=linear_solver, initialInc=0.1, maxInc=0.1)
    with pytest.warns(RuntimeWarning) if linear_solver == 'pardiso' \
            else _nowarn():
        an.static(NLgeom=True, silent=True)
    assert an.increments[-1] == 1.
    counts = an.num_negative_eigvals
    assert len(counts) == len(an.increments)
    for lbd, c, nneg, det in zip(an.increments, an.cs, counts, an.det_signs):
        np.testing.assert_allclose(c[0], lbd/p.k0, rtol=1e-12)
        assert nneg == dense_count(p.K(lbd), 6) == dense_count(p.K(lbd))
        assert nneg == int(lbd > 0.27) + int(lbd > 0.63)
        assert det == (-1)**nneg
    assert sorted(set(counts)) == [0, 1, 2]
    assert an.det_signs[-1] == 1
    # the complex pair has a negative real part at the end
    lam = np.linalg.eigvals(p.K(1.))
    assert np.any((lam.real < 0) & (np.abs(lam.imag) > 1.))
    assert an.num_perturbed_pivots == [None]*len(counts)

Linear static analysis

With a configuration-dependent load, Analysis.static(NLgeom=False) solves the geometrically linear problem \(\boldsymbol{K}_{C0} \boldsymbol{c} = \boldsymbol{F}(\boldsymbol{c})\):

def test_analysis_linear_static_with_configuration_dependent_load():
    pb = linear_follower_problem()
    an = pb.analysis()
    an.relTOL = 1e-12
    increments, cs = an.static(NLgeom=False, silent=True)
    ref = np.linalg.solve(pb.K0 - pb.Q, pb.F0)
    assert np.allclose(cs[0], ref, rtol=1e-9)
    # without c the load is the reference one, as before
    pb_old = linear_follower_problem(with_c=False)
    increments, cs = pb_old.analysis().static(NLgeom=False, silent=True)
    assert np.allclose(cs[0], np.linalg.solve(pb.K0, pb.F0))

Linear buckling with unsymmetric matrices

@pytest.mark.parametrize('sparse_solver,max_dense_size',
                         [(True, 2000), (True, 0), (False, 2000)])
def test_lb_unsymmetric_real(sparse_solver, max_dense_size):
    r"""Unsymmetric KG with real eigenvalues, against scipy.linalg.eig"""
    rng = np.random.default_rng(3)
    n = 40
    K = spd(n, rng)
    # KG = S + small skew part keeps the eigenvalues real
    S = -spd(n, rng)/n
    W = rng.standard_normal((n, n))
    KG = S + 1e-3*(W - W.T)
    assert not is_symmetric(KG)
    eigvals, eigvecs = lb(csr_matrix(K), csr_matrix(KG), silent=True,
                          sparse_solver=sparse_solver, num_eigvalues=6,
                          max_dense_size=max_dense_size)
    assert np.isrealobj(eigvals) and np.isrealobj(eigvecs)
    mu = eig(KG, K, right=False)
    lam = np.sort((-1/mu).real)
    lam = lam[lam > 0]
    assert np.allclose(eigvals[:5], lam[:5], rtol=1e-8)
    for i in range(5):
        u = eigvecs[:, i]
        r = K @ u + eigvals[i]*(KG @ u)
        assert np.linalg.norm(r) <= 1e-8*np.linalg.norm(K @ u)

Under follower loads, many negative or complex load multipliers may have a smaller modulus than the critical one, e.g. when the pressure is reversed. The sparse solver of lb() searches the lowest positive real load multiplier along the positive real axis, see Eigenvalue solvers. The cantilevers of Schweizerhof and Ramm (1984), with the matrices of panels, original and reversed:

@pytest.mark.parametrize('name, sign, approx', FOLLOWER)
def test_lb_follower_shift_search(name, sign, approx):
    """Regression test: the sparse fallback of version 0.5.3, eigs() on
    K^-1 KG, silently returned only negative or complex load multipliers for
    the reversed problems. The shift search returns the lowest positive real
    load multiplier of the full dense eig(), first, without the dense
    fallback"""
    import warnings
    from structsolve.linear_buckling import DenseFallbackWarning

    K, KG = load_follower(name, sign)
    ref = lowest_positive_real(lb(K, KG, silent=True, sparse_solver=False)[0])
    np.testing.assert_allclose(ref, approx, rtol=1e-5)
    with warnings.catch_warnings():
        warnings.simplefilter('error', DenseFallbackWarning)
        eigvals, eigvecs = lb(K, KG, silent=True, max_dense_size=0)
    assert eigvals.shape == (25,) and eigvecs.shape == (K.shape[0], 25)
    assert np.isreal(eigvals[0]) and eigvals[0].real > 0
    np.testing.assert_allclose(eigvals[0].real, ref, rtol=1e-8)
    assert lowest_positive_real(eigvals) == eigvals[0].real
    res = relative_residuals(K, KG, eigvals, eigvecs)
    assert res[0] < 1e-8 and np.all(res < 1e-5)

Flutter: the kinetic criterion

Ziegler’s double pendulum under a follower force never diverges, its stiffness matrix is never singular, and lb() finds no critical load. It loses stability by flutter at \(P_f = (7/2 - \sqrt{2}) k/l\), where two natural frequencies of freq() coalesce and become complex:

def ziegler(P, extra=0):
    r"""``K(P)`` and ``M``, with ``extra`` decoupled unit oscillators added
    such that the sparse solvers, limited to ``k < size - 1``, can be used"""
    k, l, m = K_SPRING, L_BAR, MASS
    K = np.array([[2*k - P*l, -k + P*l], [-k, k]])
    M = m*l**2*np.array([[3., 1.], [1., 1.]])
    if extra:
        rng = np.random.default_rng(0)
        K = block_diag(K, np.diag(1e3 + rng.random(extra)))
        M = block_diag(M, np.eye(extra))
    return K, M

@pytest.mark.parametrize('sparse_solver,extra', [(False, 0), (True, 10)])
def test_freq_ziegler_flutter(sparse_solver, extra):
    r"""The general solver of :func:`freq` gives real frequencies below the
    flutter load and complex ones above it; the flutter load found by
    bisection is the analytical one"""
    kw = dict(silent=True, sparse_solver=sparse_solver, num_eigvalues=4,
              sort=False)
    below = freq(*map(csr_matrix, ziegler(0.9*P_FLUTTER, extra)), **kw)[0]
    above = freq(*map(csr_matrix, ziegler(1.1*P_FLUTTER, extra)), **kw)[0]
    assert max_imag(below) <= 1e-8*np.abs(below).max()
    assert max_imag(above) > 1e-3*np.abs(above).max()
    a, b = 0.9*P_FLUTTER, 1.1*P_FLUTTER
    for _ in range(40):
        c = 0.5*(a + b)
        lam2 = freq(*map(csr_matrix, ziegler(c, extra)), **kw)[0]
        if max_imag(lam2) > 1e-6*np.abs(lam2).max():
            b = c
        else:
            a = c
    assert np.isclose(0.5*(a + b), P_FLUTTER, rtol=1e-4)

The follower pressure loads of panels use all these features, see Shell.add_pressure_load(..., follower=True).