Failure to Solve

Solving a system of simultaneous linear equations, canonically referred to as solving Ax=b in math speak, is at the heart of every equilibrium solution algorithm for nonlinear analysis. In the context of OpenSees, A is the effective tangent stiffness matrix, x is the vector of displacement increments, and b is the residual force vector.

However, there’s several reasons solving Ax=b can fail and they all boil down to rigid body motion, or mechanisms. In most cases, user error is to blame and some software will automatically correct the issue.

OpenSees is not one of those softwares. When the solution of Ax=b fails, you’ll see something like the following message among the small firestorm of output.

The first two lines of output tell you the solver failed. The first line will vary depending on which solver you use, but they all say essentially the same thing.

Now, here are common scenarios that lead to a failed solution of Ax=b.

Suppose we would like to model a cantilever column with two nodes. The solution will fail if there is a disconnected, or floating, node or if there is not sufficient boundary conditions to prevent rigid body motion, e.g., a pin instead of fixed support for the column.

Similarly, for simple beams in 3D, you will have to prevent rigid body spinning about the beam’s longitudinal axis.

Note that you could get a solution out of the two scenarios shown above by adding mass to all DOFs and performing a dynamic analysis. But rigid body dynamics is outside the scope of most earthquake engineering applications.

Also note that checking your model for floating nodes is straightforward. Build a list of connected node tags from ops.getEleTags() and ops.eleNodes(), then see if any nodes defined in the model (from ops.getNodeTags()) are not in the connected nodes list. The code below uses the “exclusive or” (XOR) operator, but there’s other ways to perform this check.

def floatingNodes():
    
    connectedNodes = []
    for ele in ops.getEleTags():
        for nd in ops.eleNodes(ele):
            connectedNodes.append(nd)
   
    definedNodes = ops.getNodeTags()

    # Use XOR operator, ^
    return list(set(connectedNodes) ^ set(definedNodes))

You can define an equivalent function in Tcl.

Now consider the same column but with a rotational spring at its base. The model has three nodes with two nodes collocated at the base for the zero length spring element.

With the two nodes at the base connected only through their rotational DOFs, the translational DOFs of these nodes must be constrained to have the same motion–otherwise, the solver will fail. This is simply another case of preventing rigid body motion.

# Constrain DOFs 1 and 2 of nodes I and J to be the same
ops.equalDOF(ndI,ndJ,1,2) 

Sure, for this particular case, we could just pin the node on top of the spring. However, the equalDOF is the way to go for models with rotational springs at beam-column joints where there are no external support conditions.

Even with a successful model definition, you can run into problems as the simulation advances through time. Consider the case of a displacement-controlled pushover analysis where the post-yield moment-rotation response of the spring is perfectly plastic. The tangent stiffness of the spring is zero and you will encounter an error with solvers that do not pivot. Use UmfPack or SparseGeneral -piv instead of ProfileSPD.

While the examples in this post were of simple cantilever columns, you’ll see similar situations in larger models of structural systems where you can have global or local mechanisms.

Although not addressed here, there are analogies for geotechnical models, most commonly p-y springs with zero stiffness or lack of equalDOF between spring nodes. You can also face issues with the equation solver based on which constraint handler you use for multi-point constraints.