OpenSees Cloud

OpenSees AMI

Finite Differences

Original Post - 27 Nov 2021 - Michael H. Scott

Visit Structural Analysis Is Simple on Substack.

A previous post showed how to compute response sensitivity by the DDM, or direct differentiation method. Comparisons with finite difference calculations verified that the DDM results were correct. In this post, I’ll dig a little deeper into finite differences.

The advantage of the finite difference method (FDM) is it will work for any model parameter–you do not have to implement any special methods within the core of OpenSees. But the FDM’s major disadvantage is it requires at least one full re-analysis for every parameter.

Why are finite difference computations important? First, some of the OpenSees reliability modules use finite differences internally to compute search directions for gradient-based algorithms. Second, if a third party application uses OpenSees as a black box, e.g., for optimization, finite differences are the most foolproof way to get derivatives of the response.

Computing finite differences in OpenSees is straightforward, as will be shown for pushover analysis of the steel moment-resisting frame below. It’s the same frame from the DDM post, where you can find more information on the model.

The yield stress, Fy, and elastic modulus, E, of the columns will be considered as parameters, although you can use as many parameters as you like. The sensitivity of the roof displacement and base shear with respect to these parameters will be computed by finite differences.

The code shown below uses the reset and updateParameter functions within a loop over the model parameters. If the setParameter and updateParameter methods are not implemented for an element or constitutive model, you can put the model definition inside the loop and manage the parameters yourself.

#
# Define model
#

ops.parameter(1)
ops.parameter(2)
for ele in [1,2,3,4]: # Column elements
   ops.addToParameter(1,'element',ele,'E')
   ops.addToParameter(2,'element',ele,'Fy')

#
# Define regular analysis options
#

# The mean response
U0 = [0]*Nsteps
for i in range(Nsteps):
   ops.analyze(1)
   U0[i] = ops.nodeDisp(6,1)

# Do finite difference for each parameter
dU = {}
eps = 1.0e-3 # Or whatever
for param in ops.getParamTags():
   ops.reset()
   dU[param] = [0]*Nsteps

   h0 = ops.getParamValue(param)
   dh = eps*h0

   # Forward perturbation
   ops.updateParameter(param,h0+dh)
   for i in range(Nsteps):
      ops.analyze(1)
      U = ops.nodeDisp(6,1)
      dU[param][i] = (U-U0[i])/dh

   # Reset to original value
   ops.updateParameter(param,h0)

The above code implements forward finite difference (FFD) where each parameter is perturbed forward and the difference is taken with respect to the mean response, i.e., \((U(\theta+\Delta\theta)-U(\theta))/\Delta\theta\). Note that it’s important to set the parameter back to its original value at the end of the loop; otherwise, you’ll be taking differences about successively different points.

You can also do central finite difference (CFD–not to be confused with computational fluid dynamics), where you perturb forward and backward, then take the difference about the center, i.e., \((U(\theta+\Delta\theta/2)-U(\theta-\Delta\theta/2))/\Delta\theta\). CFD is more accurate than FFD, but requires two re-analyses for each parameter.

for param in ops.getParamTags():
   ops.reset()
   dU[param] = [0]*Nsteps

   h0 = ops.getParamValue(param)
   dh = eps*h0

   # Forward perturbation
   ops.updateParameter(param,h0+0.5*dh)
   for i in range(Nsteps):
      ops.analyze(1)
      U = ops.nodeDisp(6,1)
      dU[param][i] = U

   # Backward perturbation
   ops.updateParameter(param,h0-0.5*dh)
   for i in range(Nsteps):
      ops.analyze(1)
      U = ops.nodeDisp(6,1)
      dU[param][i] = (dU[param][i]-U)/dh

   # Reset to original value
   ops.updateParameter(param,h0)

The sensitivity of the roof displacement computed by DDM and FFD and CFD with parameter perturbations of \(\Delta\theta=0.1\theta\) are shown below.

For the large perturbation, the finite difference approaches are not close to the DDM result, and the CFD does better than the FFD, but at twice the computational cost. Decreasing the perturbation to \(\Delta\theta=0.01\theta\), the finite difference results improve.

Further reductions in the parameter perturbation, e.g., less than \(0.001\theta\), lead to indiscernible differences in the finite difference and DDM results for the pushover analysis of this frame.