Shutting Off the Containment Unit

Original Post - 01 Mar 2023 - Michael H. Scott

Visit Structural Analysis Is Simple on Substack.

If you’ve ever used OpenSees–even if you’re a geotech–you’ve used the force-based element.

When Remo implemented the force-based element, it was the only material nonlinear frame element available in OpenSees (G3 at the time); thus, the original name nonlinearBeamColumn. Only after a standard displacement-based frame element (dispBeamColumn) was added did we change the name from nonlinearBeamColumn to forceBeamColumn. Other material nonlinear formulations have come since.

That the force-based element in OpenSees is based on the non-iterative formulation developed by Neuenhofer and Filippou (1997) is technically true. However, the non-iterative formulation is not the default behavior. By default, the maximum number of internal iterations, maxIter, is set to 10, not 1. You have to use the -iter option to make maxIter=1.

However, there are some strange boundary cases where using maxIter=1 can grind the force-based element to a halt, even for a simple model.

Run this example and watch the convergence test churn on the last dozen or so time steps.

wipe
model basic -ndm 2 -ndf 3

node 1 0 0; fix 1 1 1 0
node 2 100 0; fix 2 1 1 0

uniaxialMaterial Steel02 1 50 29000 0.005
section Fiber 1 {
   patch rect 1 20 1 -10.0 -5.0 10.0 5.0
}

geomTransf Linear 1
element forceBeamColumn 1 1 2 4 1 1 -iter 1 1e-12

integrator LoadControl 500.0
timeSeries Linear 1
pattern Plain 1 1 {
   load 2 0 0 1.0
}

test NormDispIncr 1e-6 10 1
analysis Static

analyze 100

I chose Tcl because it’s easier to use diagnostic tools like gprof and valgrind with a stand-alone executable. You’ll see the same churn on the last dozen time steps in an equivalent OpenSeesPy model.

Anyway, I set out to determine why maxIter=1 causes huge problems for this simple example.

But first, some background on the OpenSees force-based element state determination is necessary. Given basic deformations, coming from the nodal displacements by way of the geometric transformation, the element attempts three internal Newton algorithms:

Regular Newton-Raphson, always using the current element flexibility.
Newton-Raphson with the initial element flexibility on the first iteration, then current flexibility on subsequent iterations.
Newton-Raphson with the initial element flexibility on every iteration.

A few things to note about these algorithms within the force-based element:

Despite the code comments, the sequence of algorithms is programmed as 1, 3, then 2, i.e., the element will try iteration with the initial flexibility (3) before trying initial then current flexibility (2).
To give algorithm 3 a fighting chance, the number of internal iterations is set to 10*maxIter for this algorithm only.
If maxIter=1, algorithms 2 and 3 are the same, i.e., algorithm 2 will never make it to “subsequent iterations”.

So, with maxIter=1 and the 1-3-2 order of algorithms, algorithm 2, the best backup plan for failure of algorithm 1, never gets enacted. If 10 iterations with constant flexibility fails, one iteration with constant flexibility will also fail.

The failure of the 1-3-2 algorithm sequence in the force-based element state determination sets off subdivisions of those basic deformations coming from the geometric transformation. Four subdivisions are allowed, each cutting the basic deformations by an order of magnitude: 1/10, 1/100, 1/1000, then 1/10000.

If the 1-3-2 algorithm continues failing, the element will chug on the 1/10000 subdivision and the run-time for simple models goes nuts when maxIter=1.

Of course none of this is an issue with the default maxIter=10. But we want the element to work when maxIter=1, right?

The fix is simple:

Change the order of algorithms from 1-3-2 to 1-2-3 putting the initial flexibility iteration last.
Allow algorithm 2 a few more iterations, e.g., maxIter+5, to do its thing.

I’m not really sure you can call it a fix though. It’s just a more sensible backup plan for failure of algorithm 1, which is inevitable when there is material yielding and maxIter=1. So, the element is still not really non-iterative, but I’m not ready to make a go at Excalibur and all the responsibility that comes with it.

Anyway, before Frank dismisses PR #1168 and calls me a “tosser”, pull my branch and see how the updates affect your IDAs of 2D and 3D reinforced concrete frames using force-based elements with fiber sections of Concrete23 and Steel08.

To get the branch, issue these two commands in your favorite git client:

% git checkout -b mhscott-fbc-update master
% git pull git@github.com:mhscott/OpenSees.git fbc-update

Then compile OpenSees (Tcl or Py) and run your IDAs. Do the analyses run faster? Do the results change?

If the results do change, your model was fickle to begin with–consider retracting published papers and writing errata. Or, for works in progress, change the narrative. No, seriously, if your results change, I’d love to hear about it.

It’s never too late to fix a problem. And it’s never too late to revert a pull request.

If the maxIter=10 default was an IDA ghost containment unit, I’m Walter Peck, the EPA inspector played by William Atherton (who also played a prickish role in Die Hard) shutting down the Ghostbusters. However, unlike Peck, I promise not to blame you if all IDA-hell breaks loose with this change in the force-based e

OpenSees Cloud

OpenSees AMI

Shutting Off the Containment Unit