Non-linear static analysis with the Newton-Raphson method

A non-linear static analysis with Analysis requires four callables, which calculate the external force vector, the internal force vector, the constitutive stiffness matrix and the geometric stiffness matrix. The tangent stiffness matrix \([K_T] = [K_C] + [K_G]\) should be the exact derivative of the internal force vector, which gives quadratic convergence to the Newton-Raphson method, described in newton_raphson._solver_NR().

The class below, used in the structsolve unit tests, wraps analytic problems into these callables:

class Problem(object):
    """Wrap analytic callbacks and record every call to ``calc_fint``

    Each recorded call stores the number of increments already accepted by
    the solver, allowing the iterations of each step to be counted.

    """
    def __init__(self, fext, fint, kC, kG):
        self.fext = fext
        self.fint = fint
        self.kC = kC
        self.kG = kG
        self.an = None
        self.fint_calls = []

    def calc_fext(self, inc=1., silent=True):
        return inc*self.fext

    def calc_fint(self, c, silent=True):
        self.fint_calls.append((len(self.an.increments), c.copy()))
        return self.fint(c)

    def calc_kC(self, c=None, NLgeom=False, silent=True):
        return csc_matrix(self.kC)

    def calc_kG(self, c=None, NLgeom=False, silent=True):
        return csc_matrix(self.kG(c))

    def analysis(self, **kwargs):
        an = Analysis(calc_fext=self.calc_fext, calc_fint=self.calc_fint,
                      calc_kC=self.calc_kC, calc_kG=self.calc_kG)
        for k, v in kwargs.items():
            setattr(an, k, v)
        self.an = an
        return an

    def calls_per_step(self):
        """Number of ``calc_fint`` calls spent to reach each accepted step"""
        steps = [s for s, _ in self.fint_calls]
        return [steps.count(i) for i in range(len(self.an.increments))]

    def solve_reference(self, lbd=1., c0=None, tol=1e-15):
        """Plain dense Newton-Raphson to obtain a reference solution"""
        c = np.linalg.solve(self.kC, lbd*self.fext) if c0 is None else c0
        for _ in range(100):
            R = lbd*self.fext - self.fint(c)
            if np.linalg.norm(R) <= tol*np.linalg.norm(lbd*self.fext):
                return c
            c = c + np.linalg.solve(self.kC + self.kG(c), R)
        raise RuntimeError('reference solution did not converge')

One of these problems is a chain of springs with a cubic hardening, where \(\{F_{int}\} = [K]\{u\} + a \{u\}^3\):

def hardening_springs(P=50.):
    """Chain of 4 springs with cubic hardening: ``fint = K u + a u**3``"""
    n = 4
    K = 100.*(2*np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1))
    K[-1, -1] = 100.
    a = 50.*np.arange(1, n+1)
    fext = np.zeros(n)
    fext[-1] = P
    return Problem(fext, lambda u: K @ u + a*u**3, K,
                   lambda u: np.diag(3*a*u**2))

The full load is applied in a single step, starting from the linear solution, and the order of convergence of the iterations is verified:

def test_nr_full_load_quadratic_convergence():
    p = hardening_springs()
    p.analysis(initialInc=1.)
    p.an.static(NLgeom=True, silent=True)
    check_increments_and_equilibrium(p)
    assert len(p.an.increments) == 1

    ref = p.solve_reference()
    errors = [np.linalg.norm(c - ref)/np.linalg.norm(ref)
              for _, c in p.fint_calls]
    orders = convergence_orders(errors)
    assert len(orders) >= 1, errors
    assert min(orders) > 1.6, (orders, errors)
    # plain Newton from the linear predictor needs 9 residual evaluations to
    # reach machine precision; relTOL requires at most that
    assert len(p.fint_calls) <= 9

The analysis can also be performed in increments, with initialInc defining the size of the first load increment:

def test_nr_incremental_hardening():
    p = hardening_springs()
    p.analysis(initialInc=0.1)
    p.an.static(NLgeom=True, silent=True)
    check_increments_and_equilibrium(p)
    # a good predictor with an exact tangent needs few iterations per step
    assert max(p.calls_per_step()) <= 5, p.calls_per_step()