3.2. Other Critical Time Step Methods
Currently, there are four methods for calculating a critical time step for Presto. First, there is the traditional element-based method. We know that, in general, the element-based time step in Presto can give a fairly conservative estimate for the time step. Second, there is a node-based method for giving a critical time step estimate. Depending on the problem, the node-based method may or may not give an estimate for the critical time step that approaches the theoretical maximum value for a particular model. Although the node-based method can give a larger critical time step estimate than the element-based method, the node-based estimate may still be significantly lower than the maximum theoretical time step for a problem. Finally, there are two methods that use an estimate for the maximum eigenvalue to obtain an estimate for the critical time step. From the maximum eigenvalue, it is possible to derive the theoretical maximum critical time step for a problem via the formula
where \(\Delta t_{crit}\) is the critical time step and \(\theta_{max}\) is the maximum eigenvalue. The two methods employed in Presto to give an estimate for the maximum eigenvalue are the Lanczos method and the power method. The Lanczos method can give an accurate estimate of the maximum eigenvalue for a problem using a very small number of Lanczos vectors compared to the total number of degrees of freedom in a problem. The power method uses a simple iterative process to obtain an estimate for the maximum eigenvalue. The power method is not as powerful a mathematical technique as the Lanczos method for obtaining a maximum eigenvalue estimate. However, it does present another viable option for obtaining a maximum eigenvalue estimate for certain problems and has therefore been included as another method for obtaining the maximum eigenvalue estimate. See [[1], [2]] for more information on the Lanczos method and other related methods in computational linear algebra.
The use of a critical time step from the node-based method or a method based on the maximum eigenvalue estimate is desirable because the larger critical time steps produced by these methods (compared to the element-based method) reduce CPU time. Both methods, however, are not cost effective if they are called every time step to give a critical time step estimate. The cost of doing one node-based estimate or one maximum eigenvalue estimate for the critical time step will not offset the cost benefit of the increase to the critical time step (compared to the element-based time step estimate) over a single time step. Hence, there must be some scheme for
{calling these methods only periodically throughout a calculation and}
{maintaining a larger estimate (than the element-based estimate) for the critical time step in between these calls}
if we are to gain a net benefit from the increase in the critical time step these methods can produce.
If you want to use the maximum eigenvalue estimate for calculating the critical time step and your problem has long periods with a relatively stable time step estimate, the Lanczos method for calculating the maximum eigenvalue will be the preferred method to use. If you have a problem for which the changes in the time step should be monitored more frequently, the power method will be better suited for this problem than the Lanczos method. The preferred use for each of these two different eigenvalue methods will become obvious as you read the background material for these two methods in Section 3.2.1 and Section 3.2.2. Both the Lanczos method and the power method require the use of a variety of control parameters. For some of these control parameters, techniques have been devised to automatically calculate values for these parameters in what should be a robust and reliable method for a wide range of problems. Other parameters are defaulted to values that will optimize the use of a particular eigenvalue calculation method (Lanczos method or power method) with the type of problems that best suit the method that is employed. The default values for these control parameters may change as we gain more experience in using the Lanczos method and the power method.
Detailed discussions of the Lanczos method, the power method, and the node-based method follow. The implementation of the Lanczos method and the power method have many similarities.
3.2.1. Lanczos Method
The Lanczos method, as implemented in Presto, is outlined here so that certain aspects of it can be referenced in subsequent parts of this chapter. In the following set of equations, \(\mathbf{K}_{T}\) is a tangent stiffness matrix, \(\mathbf{M}\) is the mass matrix, and \(\mathbf{r}_{0}\) is an arbitrary starting vector.
The details for this form of the Lanczos method are described in [[3]]. Notice that the Lanczos method is an iterative method. If we use the Lanczos method in a code like Presto to compute the maximum eigenvalue for a particular finite element model, the number of iterations required to give a good estimate for the maximum eigenvalue will depend on the size of the finite element model (the number of nodes and elements), the types of elements in the model, and the material types and material properties used in the model.
The maximum eigenvalue for a particular finite element model gives us the largest estimate we can obtain for the critical time step. Hence, the time step estimate derived from the maximum eigenvalue is our “best” estimate for the critical time step. The estimate for the critical time step based on the maximum eigenvalue can be significantly larger than an element-based time step estimate. However, computing the critical time step (for a given time step) with the Lanczos method is more expensive than computing the critical time step with element-based calculations. Over one time step, it is not possible to recoup the cost of the Lanczos calculations with the increase in the size of the time step over the element-based estimate. Using the Lanczos method for estimating the critical time step in an explicit, transient dynamics code requires a methodology that effectively addresses the computational costs. The following sections outline a cost-effective approach to using the Lanczos method in an explicit, transient dynamics code.
3.2.1.1. Lanczos Method with Constant Time Steps
To explain how to use the Lanczos method in a cost-efficient manner, we must first establish the computational cost of using the Lanczos method. As indicated previously, computing the maximum eigenvalue for a finite element model requires some number of iterations (each iteration in the Lanczos method produces a Lanczos vector) to obtain a good estimate for the maximum eigenvalue. The cost of an iteration (Lanczos vector) is approximately the cost of an internal force calculation. Notice, in preceding equations for the Lanczos method, that the Lanczos method requires the product of the tangent stiffness matrix \(\mathbf{K}_{T}\) with a vector \(\mathbf{p}\). In Presto, we do not construct a tangent stiffness matrix. Instead, we provide the vector \(\mathbf{p}\) for the internal force calculations. The internal force calculations give us the desired matrix \(\times\) vector product of \(\mathbf{K}_{T} \mathbf{p}\).
Over a given time step, the cost of an internal force calculation is the major computational cost. (This assumes no contact. The addition of contact introduces another computationally expensive process into a time step. For our initial discussion, we ignore the cost of contact.) To use the Lanczos method to get a critical time step estimate, one has to call the Lanczos method for some given time step. If the Lanczos method made \(n\) iterations to get a good estimate for the maximum eigenvalue, then the overall cost of the time step would be approximately \(n+1\) times the cost of the internal force calculation. The cost of the time step would be the \(n\) internal force calculations for Lanczos and the actual internal force calculation to advance the time step. If \(n\) is 20 (a minimum for typical problems), the cost of the time step becomes 21 times the cost of the internal force calculation. The critical time step estimate based on the maximum eigenvalue would have to be at least 21 times greater than the element-based critical time step to recoup the cost of the maximum eigenvalue calculation. A typical value for the critical time step estimate based on the maximum value is more in the range of 1.1 to 2.0 times the element-based critical time step estimate. Obviously, the Lanczos method is too expensive to call at every time step for a critical time step estimate.
To explain how the Lanczos method can effectively be used in an explicit, transient dynamics code, we begin with a simple case study. In this case study, we compute the critical time step using the Lanczos method at some time step and then assume that this critical time step value remains constant for a subsequent number of time steps, \(n_{L}\). We only call the Lanczos method once during the \(n_{L}\) time steps. (In reality, the critical time step in an explicit, transient dynamics code like Presto changes with each time step. We address this issue of the changing time step when we present the details for using the Lanczos method in a cost-effective manner.) For this case study, we also assume that the computational cost of an element-based estimate for the critical time step is part of the cost of an internal force calculation. The cost of the element-based time step estimate is a small part of the overall internal force calculations. Finally, for our initial discussion, we assume no contact. We address the issue of contact further in the discussion.
Assume that the Lanczos method computes a global estimate for the critical time step of \(\Delta t_{L}\), which is the value to be used for \(n_{L}\) time steps. At the end of the \(n_{L}\) time steps, the analysis time for the code has been incremented by an amount \(\Delta T\), which is computed as
If the element-based estimate for the time step is \(\Delta t_{e}\) and the number of time steps required to increment the analysis time by \(\Delta T\) is \(n_{e}\), then, for the element-based time step, we have
Because the Lanczos estimate for the critical time step is larger than the element-based estimate, we know that \(n_{e} > n_{L}\). Let us define the ratio \(r\) as
The ratio \(r\) is greater than 1. Now that we have determined the relation between the number of steps required for a Lanczos-based critical time step estimate versus the element-based critical time step estimate to achieve the same analysis time increment, let us examine the computational costs for these two cases in terms of CPU time. Designate the CPU cost for a time step as \(\Delta t_{IF}\). If the number of Lanczos vectors required to obtain the critical time step estimate is \(N_{L}\), then the total computational cost of the Lanczos method and the \(n_{L}\) time steps is
If we use the element-based method, the total computational cost is
Recall that we have chosen \(n_{L}\) and \(n_{e}\) so that we have the same analysis time increment \(\Delta T\) even though we have different critical time steps. Now, we must determine the point at which the computational cost for the Lanczos-based critical time step calculations is the same as the cost for the element-based critical time step calculations. This is the point at which
If we rearrange the above equation to eliminate \(\Delta t_{IF}\) and make use of the ratio \(r\), then we obtain
Consider the case of \(r = 1.25\) and \(N_{L} = 20\). When \(n_{L} = 80\) and \(n_{e} = 100\), the above equations show that the calculations with the Lanczos-based critical time step and the calculations with the element-based time step give the same analysis time for the same computational expense. If we use the Lanczos-based critical time step \(\Delta t_{L}\) for more than eighty iterations, then the Lanczos-based approach becomes cost effective. Our above equations have established the “break-even” point at which it becomes cost effective to use the Lanczos method to reduce computational costs by overcoming the initial cost of the Lanczos calculations with the larger critical time step.
We can build on what we have done thus far to account for contact. Suppose that the computational cost of contact over a time step is some multiple \(m\) of the computational cost of the internal force calculation \(\Delta t_{IF}\). Then the point at which the computational cost for the Lanczos-based calculations is the same as the computational cost for the element-based calculations is
For the case with contact,
Again, consider the case of \(r = 1.25\) and \(N_{L} = 20\). Assume the computational cost of contact calculations is the same as an internal force calculation (\(m = 1\)). For these values, the break-even point is \(n_{L} = 40\) and \(n_{e} = 50\). The added computational cost of the contact calculations results in reaching the break-even point with a smaller number of iterations when compared to the case with no contact.
The above derivations let us calculate a break-even point based on our assumptions of a constant critical time step. Considering that a typical analysis will run for tens of thousands of time steps, something on the order of 100 steps represents a reasonable number of steps to recoup the cost of the Lanczos calculations. Whether the cost of the Lanczos calculations can be recouped in something on the order of 100 calculations depends heavily upon \(N_{L}\). If \(N_{L}\) is sufficiently small, we can recoup the cost of the Lanczos calculations in a reasonable number of steps.
Some computational studies indicate that \(N_{L}\) is in an acceptable range for many problems. The Lanczos method computes a good estimate for the maximum eigenvalue with a small number of Lanczos vectors, \(N_{L}\), compared to the number of degrees of freedom in a problem. Some component studies show that, for a problem with between 250,000 and 350,000 degrees of freedom, one can obtain a good estimate for the maximum eigenvalue with only twenty Lanczos vectors. A large-scale study of a model involving 1.7 million nodes (5.1 million degrees of freedom) showed that only forty-five Lanczos vectors were required to obtain a good estimate of the maximum eigenvalue. These examples demonstrate that the number of Lanczos vectors required for a good maximum eigenvalue estimate is very small when compared to the number of degrees of freedom for a problem. When \(N_{L}\) is in the range of twenty to forty-five, Equations (3.1) and (3.2) show that, with an increase in the time step of about 1.2 to 1.25, we can recoup the cost of the Lanczos method in a reasonable number of time steps.
Now that we have determined we can recoup the cost of the Lanczos calculations in a reasonable number of time steps, let us look at the issue of reusing a Lanczos-based estimate in some manner.
3.2.1.2. Controls for Lanczos Method
As indicated in the above discussion, the Lanczos method can be used in a cost-effective manner in an explicit, transient dynamics code if a Lanczos calculation can be performed and the Lanczos-based estimate for the critical time step can be reused in some way over many subsequent time steps. This section presents an approach for reusing a Lanczos-based estimate over a number of time steps so that we maintain a critical time step estimate that is close to the theoretical maximum value in between the calls to the Lanczos method. The approach discussed here makes use of the element-based critical time step estimate at each time step.
We start our approach with a Lanczos calculation to determine the maximum eigenvalue. The Lanczos method converges to the maximum eigenvalue from below, which means that the method underestimates the maximum eigenvalue. Because the critical time step depends on the inverse of the maximum eigenvalue, we overestimate the critical time step. It is necessary, therefore, to scale back the critical time step estimate from the Lanczos method so that the calculations in the explicit time-stepping scheme do not become unstable. Our approach for determining a scaled-back value for the maximum critical time step makes use of the element-based time step estimate. Again, let \(\Delta t_{L}\) be the critical time step estimate from the Lanczos method and \(\Delta t_{e}\) be the critical time step estimate from the element-based calculations. The scaled-back estimate for the critical time step, \(\Delta t_{s}\), is computed from the equation
where \(f_{s}\) is a scale factor. (A reasonable value for \(f_{s}\) ranges from 0.9 to 0.95 for our problems.) This value of \(f_{s}\) puts \(\Delta t_{s}\) close to and slightly less than a theoretical maximum critical time step.
Once \(\Delta t_{s}\) is determined, the ratio
is computed. This ratio is then used to scale subsequent element-based estimates for the critical time step. If \(\Delta t_{e(n)}\) is the \(n^{th}\) element-based critical time step after the time step where the Lanczos calculations are performed, then the \(n^{th}\) time step after the Lanczos calculations, \(\Delta t_{(n)}\), is
The ratio \(t_{r}\) is used until the next call to the Lanczos method. The next call to the Lanczos method is controlled by one of two mechanisms. With the first mechanism, the user can set the frequency with which the Lanczos method is called. The user can set a parameter so that the Lanczos method is called only once every \(n\) time steps. This number remains fixed throughout an analysis. With the second mechanism, the user can control when the Lanczos method is called based on changes in the element-based critical time step. For this second mechanism, the change in the element-based critical time step estimate is tracked. Suppose the element-based critical time step at the time the Lanczos method was called is \(\Delta t_{e}\). At the \(n^{th}\) step after the call to the Lanczos method, the element-based critical time step is \(\Delta t_{e(n)}\). If the value
is greater than some limit set by the user, then the Lanczos method will be called. If there is a small, monotonic change in the element-based critical time step over a large number of time steps, this second mechanism will result in the Lanczos method being called. Or if there is a large, monotonic change in the element-based critical time step over a few time steps, the Lanczos method will also be called.
These two mechanisms for calling the Lanczos method can be used together. For example, suppose the second mechanism (the mechanism based on a change in the element-based time step) results in a call to the Lanczos method. This resets the counter for the first mechanism (the mechanism using a set number of time steps between calls to the Lanczos method).
This approach for reusing a Lanczos-based time step estimate has been implemented in Presto, and used for a few studies. One of the component studies, as indicated previously, used the same scale factor for \(n_{L} = 1700\) iterations. The break-even point for this problem is \(n_{e} = 45\) time steps (not accounting for contact, which was a part of the component modeling). For this particular problem, the extended use of the Lanczos estimate reduced the computational cost to 56% of what it would have been with the element-based time step.
Not all problems will lend themselves to reusing one Lanczos-based estimate for thousands of time steps. However, if it is possible to use the Lanczos-based estimate for two to three times the number of time steps required to reach the break-even point, we begin to see a noticeable reduction in the total number of time steps required for a problem.
3.2.1.3. Scale Factor for Lanczos Method
When the Lanczos method is called for a given time step, it must appear that the calculations are using the constant tangent stiffness matrix \(\mathbf{K}_{T}\) for all iterations. As indicated previously, we use the internal force calculations to generate the product \(\mathbf{K}_T \mathbf{p}_{j}\) (for the \(j^{th}\) iteration) in the Lanczos calculations. Any vector \(\mathbf{p}_{j}\), as calculated by the Lanczos method, may be such that it represents large-strain behavior and moves the internal force calculations into a nonlinear regime. It is necessary to scale the \(\mathbf{p}_{j}\) vectors so that the internal force calculations are in a small-strain regime, which makes it appear that we are working with a constant tangent stiffness matrix. The vectors \(\mathbf{p}_{j}\) must be scaled so that they represent velocities associated with small strain. When properly scaled vectors are sent to the internal force calculation, the internal force calculation effectively becomes a matrix \(\times\) vector product with a constant tangent stiffness matrix for all iterations during a given call to the Lanczos method.
The scale factor for the \(\mathbf{p}_{j}\) vectors, which we will designate as \(v_{sf}\), must not be too small, as this will create round-off problems and give a bad estimate for the critical time step. If the scale factor is too large, we violate the above restriction of a constant tangent stiffness matrix.
There are two approaches for controlling the scale factor when the Lanczos method is used to compute the maximum eigenvalue. These approaches are discussed in Section 3.2.1.5.
3.2.1.4. Accuracy of Eigenvalue Esimate
Every time a new Lanczos vector is computed, we obtain an additional eigenvalue for our model and, in general, a better estimate for the maximum eigenvalue (and hence a better estimate for the critical time step). The Lanczos method can compute a good value for the maximum eigenvalue with a very small number of total computed eigenvalues compared to the number of degrees of freedom in a problem. There are examples, as previously indicated, of problems with 250,000 to 350,000 degrees of freedom where we have obtained a good estimate of the critical time step with twenty eigenvalues. In one problem with 5.1 million degrees of freedom, we obtained a good estimate of the critical time step with forty-five eigenvalues.
A user could, in theory, determine a reasonable number of eigenvalues necessary for obtaining a good estimate of the maximum eigenvalue based on the above information on model size and the number of eigenvalues required for a good maximum eigenvalue estimate. The user could test the validity of the choice of the number of eigenvalues by increasing the number of eigenvalues slightly and comparing the maximum eigenvalue estimate obtained with the larger number of eigenvalues to the original maximum eigenvalue estimate obtained with the smaller number of eigenvalues. If the change in the two maximum eigenvalue estimates (larger versus smaller total number of eigenvalues) is small, then the original estimate for the number of eigenvalues is reasonably accurate.
As an alternative to directly specifying the number of eigenvalues to be computed, a convergence tolerance could be set on the change in the magnitude of the maximum eigenvalue as additional eigenvalues (Lanczos vectors) are computed. Let \(\theta_{max\_n}\) be the maximum eigenvalue calculated corresponding to \(n\) eigenvalues (Lanczos vectors), and let \(\theta_{max\_n+1}\) be the maximum eigenvalue corresponding to \(n+1\) eigenvalues (Lanczos vectors). The eigenvalues would be computed until
is less than some tolerance. (We will now refer to the value of (3.5) as the convergence measure.) If we calculate the convergence measure for a sequence of maximum eigenvalues computed by the Lanczos method, we will not necessarily see a monotonic decrease in the convergence measure for all time. Typically, the convergence measure will initially show a monotonically decreasing value for an increasing \(n\), and then the convergence measure will show a monotonically increasing value for an increasing \(n\). (The phase in which the convergence measure increases is usually very small, about one to two iterations.) After the phase in which the convergence measure monotonically increases, there is then a phase in which the convergence measure begins a long (over many iterations) monotonic decrease. This behavior (decrease of the convergence measure, slight increase, then long monotonic decrease) is typical of Krylov methods, of which the Lanczos method is one. Using the convergence measure to determine the number of eigenvalues to be computed for a problem is a reasonable option, but some care must be taken in setting the tolerance for the convergence measure.
Both of the techniques discussed above are offered as a way to set the maximum number of eigenvalues required to obtain an accurate eigenvalue estimate for a model. These techniques are discussed further in Section 3.2.1.5.
3.2.1.5. Lanczos Parameters Command Block
BEGIN LANCZOS PARAMETERS <string>lanczos_name
STARTING VECTOR = <string>STRETCH_X|STRETCH_Y|STRETCH_Z|
ISOTHERMAL(ISOTHERMAL)
INCREASE OVER STEPS = <integer>incr_int(5)
NUMBER EIGENVALUES = <integer>num_eig(100)
EIGENVALUE CONVERGENCE TOLERANCE = <real>converge_tol(0.5e-3)
SMALL STRAIN = <real>small_strain(1.0e-3)
VECTOR SCALE = <real>vec_scale(1.0e-5)
SCALE FACTOR = <real>time_scale(1.0)
UPDATE ON TIME STEP CHANGE = <real>tstep_change(0.10)
UPDATE STEP INTERVAL = <integer>step_int(500)
FORCE GLOBAL TIMESTEP = <string>ON|OFF(OFF)
END [LANCZOS PARAMETERS <string>lanczos_name]
If you use the Lanczos method to compute a critical time step, there should be only one LANCZOS PARAMETERS command block, and it should appear in the region. If you have a LANCZOS PARAMETERS command block, you should not specify a NODE BASED TIME STEP PARAMETERS or a POWER METHODS PARAMETERS command block. If you use the Lanczos method to compute the critical time step, the time step increase factor (default or user-specified) will be used to control the increase in the time step estimate; the time step scale factor (default or user-specified) will be used to scale the final time step. These factors are specified in the PARAMETERS FOR PRESTO REGION portion of a TIME STEPPING BLOCK in the TIME CONTROL command block.
The Lanczos method requires some type of starting vector. This is determined by the STARTING VECTOR command line. The various options available for this command line will generate a displacement vector that stretches your model in the \(x\)-, \(y\)-, or \(z\)-direction, or in all three directions (\(x\), \(y\), and \(z\)) at once. (STRETCH_X stretches the model in the \(x\)-direction, STRETCH_Y stretches the model in the \(y\)-direction, STRETCH_Z stretches the model in the \(z\)-direction, and ISOTHERMAL stretches the model in all three directions at once.) The displacement vector then serves as a basis for generating the starting vector \(r_{0}\) in the initialization phase of the Lanczos method. The Lanczos method appears to be fairly insensitive to the choice of a starting vector. The default starting vector is the ISOTHERMAL option. You may encounter cases where the use of one of the other options for a starting vector—STRETCH_X, STRETCH_Y, or STRETCH_Z—may result in a slight increase in accuracy for the critical time step estimate for a given number of eigenvalues. These other options (STRETCH_X, STRETCH_Y, and STRETCH_Z) are offered for these special cases.
The INCREASE OVER STEPS command line determines how many steps, i.e., incr_int, are used to transition from an element-based critical time step estimate to the Lanczos-based estimate at the beginning of an analysis. The user may want to increase from the element-based estimate to the Lanczos-based estimate over a number of time steps if the difference between these two estimates is large. The value of incr_int for this command line defaults to 5.
As indicated in Section 3.2.2.4, some number of eigenvalues must be computed by the Lanczos method to obtain a good estimate for the maximum eigenvalue for a given finite element model. As one option, the user can specify the number of eigenvalues to be computed by the Lanczos method by using the NUMBER EIGENVALUES command line. The default value for num_eig in this command line is 100. As an alternative, the user can specify a tolerance on the convergence measure defined in (3.5). The tolerance is set with the EIGENVALUE CONVERGENCE TOLERANCE command line. The default value for the convergence tolerance, converge_tol, is \(0.5 \times 10^{-3}\). (This default value may change in the future as we gain more experience with the Lanczos method.)
By default, the convergence measure (with the default value of \(0.5 \times 10^{-3}\)) is used to determine the number of eigenvalues to be computed. If neither the EIGENVALUE CONVERGENCE TOLERANCE command line nor the NUMBER EIGENVALUES command line appears in the LANCZOS PARAMETERS command block, then the convergence measure will be used to determine the number of eigenvalues. If the EIGENVALUE CONVERGENCE TOLERANCE command line appears in the command block with a specified tolerance other than the default value, then the number of eigenvalues computed will be determined by the convergence measure and the user-specified tolerance. If the NUMBER EIGENVALUES command line appears in the LANCZOS PARAMETERS command block, then the user-defined number of eigenvalues will be used to set the number of eigenvalues computed by the Lanczos method.
The recommended technique for determining the number of eigenvalues is to use the convergence measure option with the default tolerance. If you use this technique, neither the NUMBER EIGENVALUES command line nor the EIGENVALUE CONVERGENCE TOLERANCE command line should be included in the LANCZOS PARAMETERS command block.
As indicated in Section 3.2.1.3, it is necessary to scale one of the vectors (with a scale factor \(v_{sf}\)) in the Lanczos calculations so that it appears we are always working with a constant tangent stiffness each time we call the Lanczos method. There are two approaches for determining the scale factor. With the first approach, the user can specify the value of a small strain for a particular model by using the SMALL STRAIN command line. (As long as the internal force calculation inside the Lanczos method is working with a vector that approximates small-strain behavior, it will appear that we are working with a constant tangent stiffness inside the Lanczos method.) The value for \(v_{sf}\) is then computed based on your specification for what constitutes a small strain in the model. (The default value for small_strain is \(1.0 \times 10^{-3}\).) With the second approach for determining the scale factor, the user can specify a value for \(v_{sf}\) directly by using the VECTOR SCALE command line. The default value for \(v_{sf}\), the vec_scale parameter in that command line, is \(1.0 \times 10^{-5}\). Many tests have established that this value for the scale factor, \(1.0 \times 10^{-5}\), works well for a range of models encountered at Sandia National Laboratories. If the user wants to use the option of directly specifying the scale factor, it is possible to determine whether a particular scale factor is suitable for a particular problem. Take a scale factor \(v_{sf}\), plus values on either side of it, say, \(0.9 \times v_{sf}\) and \(1.1 \times v_{sf}\). If all three of these scale factor values produce almost the same estimate for the critical time step for a particular model, then the value for \(v_{sf}\) meets the criterion for an acceptable scale factor.
By default, the small-strain approach (with the default value of \(1.0 \times 10^{-3}\)) is used to determine a scale factor. If neither the SMALL STRAIN command line nor the VECTOR SCALE command line appears in the LANCZOS PARAMETERS command block, then the small-strain value will be used to calculate the scale factor. If the SMALL STRAIN command line appears in the command block with any value other than the default value and the VECTOR SCALE command line does not appear in the command block, then the vector scale will be determined by the user-specified value for the small strain. If the VECTOR SCALE command line does appear in the LANCZOS PARAMETERS command block, then the vector scale will be determined by the user-specified value for the vector scale.
The recommended way of determining the scale factor is to use the small-strain approach with the default small-strain value. If you follow this advice, neither the SMALL STRAIN command line nor the VECTOR SCALE command line should be included in the LANCZOS PARAMETERS command block.
As indicated previously, the scale factor must not be too small, as this will create round-off problems and give a bad estimate for the critical time step. If the scale factor is too large, we violate the above restriction of a constant tangent stiffness matrix.
The SCALE FACTOR command line sets the factor \(f_{s}\) in (3.3). The value for \(f_{s}\), i.e., time_scale, is set to 1.0.
The UPDATE ON TIME STEP CHANGE command line sets the value for \(\Delta t_{lim}\) in (3.4). If the change in the element-based critical time step estimate as given by (3.4) exceeds the value for \(\Delta t_{lim}\), then the Lanczos method is called for a new estimate for the critical time step. The default value for \(\Delta t_{lim}\), i.e., tstep_change, is 0.10.
The UPDATE STEP INTERVAL command line sets the number of step intervals at which the Lanczos method is called. If step_int is set to 1000, the Lanczos method will be called every one thousand steps to compute an estimate for the critical time step. (The default value of step_int is 500.) This control mechanism interacts with the control established by the UPDATE ON TIME STEP CHANGE command line. Suppose we have set step_int to 1000, and we have computed 800 steps since the last call to the Lanczos method. If \(\Delta t_{lim}\) has been set to 0.15 and we exceed this value at step 800, then the change in the element-based time step will result in the Lanczos method being called. The counter for keeping track of the number of step intervals since the last Lanczos computation will be reset to zero. The next call to the Lanczos method will occur again in one thousand steps, unless we again exceed the change in the element-based time step.
3.2.2. Power Method
The power method is a simple iterative technique that gives an estimate for the maximum eigenvalue. Consider the following form of the eigenvalue problem, which is useful for our explicit transient dynamics problem.
Reformulate the problem as
Denote the matrix product \(\mathbf{M}^{-1} \mathbf{K}_{T}\) as the matrix operator \(\mathbf{A}\). The eigenvalue problem is now
The maximum eigenvalue for the problem in (3.6) can be computed with the following iterative process:
The value \(\theta_{i}\) is the \(i^{th}\) estimate for the maximum eigenvalue. As \(n\) increases, the power method should yield a more accurate estimate for the maximum eigenvalue. Similar to the Lanczos method, the power method must iterate at some given time step to obtain an estimate for the maximum eigenvalue at that time step. The details for the above form of the power method are described in [[1]].
Like the Lanczos method, the power method requires the use of an internal force calculation. The power method, however, is simpler than the Lanczos method. The power method requires that one vector and one scalar be computed. The Lanczos method requires the computation of five vectors and two scalars. Further, the power method yields a direct estimate for the maximum eigenvalue, whereas the Lanczos method requires that the maximum eigenvalue be extracted from a tridiagonal matrix.
If the sequence of eigenvalues for an analysis problem is
then the rate of convergence to the maximum eigenvalue depends on the ratio
Depending on the magnitudes of \(\mid \lambda_{1} \mid\) and \(\mid \lambda_{2} \mid\), the rate of convergence could be rather slow. The power method, for a given problem, may not converge to as accurate an estimate for the maximum eigenvalue as quickly as the Lanczos method. We should be careful, however, in how we interpret this convergence property for the power method. For many of our problems, we probably have a large number of eigenvalues clustered slightly below the maximum eigenvalue. The power method may give a reasonable estimate for the maximum eigenvalue in a relatively small number of iterations. The power method may converge quickly to some eigenvalue estimate that is slightly below the maximum eigenvalue. Convergence to the maximum eigenvalue, in terms of accuracy to many digits, might then be very slow. However, one would still have an adequate estimate for the maximum eigenvalue to set a value for the critical time step.
One of the advantages of the power method is that we can use the last value for \(\mathbf{x}_{i}\) from a previous call to the power method as a starting vector for our current call to the power method and converge quickly (for some problems) to an accurate estimate for the maximum eigenvalue. Consider the following scenario. We make an initial call to the power method using some arbitrary \(\mathbf{x}_{0}\) starting vector. We obtain a good estimate for the maximum eigenvalue after \(n\) iterations. After some number of time steps, we again call the power method and we use \(\mathbf{x}_{n}\) from the previous call to the power method as our current starting vector. Because we are solving nonlinear problems, the matrix operator, \(\mathbf{A}\) , has most likely changed from the time when we made the initial call to the power method. If the changes in the matrix operator are small, then \(\mathbf{x}_{n}\) should be an excellent choice for the starting vector, and we should converge quickly to a good estimate for the maximum eigenvalue. This behavior, quick convergence to a good eigenvalue estimate using a prior value of \(\mathbf{x}_{n}\) as a starting vector, has been observed in some sample problems. We discuss this issue of restarting the power method using a previous \(\mathbf{x}_{n}\) in more detail in Section 3.2.2.1.
Obtaining a cost-effective and user-friendly implementation of the power method requires addressing the same issues discussed above for an automated implementation of the Lanczos method. First, we must determine some scheme to call the power method on an intermittent basis so that it is cost effective. The control mechanism for calling the Lanczos method on an intermittent basis can also be used for the power method. Second, the \(\mathbf{x}\) vector must be scaled appropriately before it is used in the internal force calculations. (For the Lanczos method, we must scale the vector \(\mathbf{p}\).) Finally, we need to terminate the power method. A tolerance on the value in (3.5) can be set to terminate the method. (Setting a tolerance on (3.5) is the same approach used for the Lanczos method.)
3.2.2.1. Power Method with Constant Time Steps
The cost of making one iteration for the power method is approximately the cost of an internal force calculation. As indicated previously, the cost of one iteration for the Lanczos method (the calculation of a Lanczos vector) is also approximately the cost of an internal force calculation. Therefore, the equations developed in Section 3.2.1.1 are also applicable for determining when the power method can be used in a cost-effective manner. Recall that the equation
in Section 3.2.1.1 gives the number of time steps for the break-even point at which it becomes cost effective for the case of no contact to use the Lanczos method to reduce computational costs by overcoming the initial cost of the Lanczos calculations with the larger critical time step. In the above equation, \(N_{L}\) is the number of Lanczos iterations, and \(r\) is the ratio of the Lanczos-based critical time to the element-based time step estimate. This equation is directly applicable to the power method because the cost of one iteration in the power method is similar to the cost of one iteration in the Lanczos method. (The variable \(N_{L}\) becomes the number of iterations for the power method.)
Now let us return to the scenario in Section 3.2.2 where we look at restarting the power method with the last value for \(\mathbf{x}_{i}\) from a previous call to the power method. We make an initial call to the power method using some arbitrary \(\mathbf{x}_{0}\) starting vector. We obtain a good estimate for the maximum eigenvalue after \(n\) iterations. Our estimate for the critical time step based on the maximum eigenvalue gives us a value of 1.25 for \(r\). After some number of time steps, we again call the power method, and we use \(\mathbf{x}_{n}\) from the previous call to the power method as our current starting vector. In our second call to the power method, the power method converges to a good estimate for the maximum eigenvalue in two iterations. The break-even point where we recoup the cost of the power method (in our second call) is ten time steps. If the power method converges to a good estimate for the maximum eigenvalue in three iterations, then the break-even point where we recoup the cost of the power method (in our second call) is fifteen time steps. There are some problems where we can obtain a good estimate for the maximum eigenvalue in two or three iterations when we use the last value for \(\mathbf{x}_{n}\) from a previous call to the power method as our starting vector for the current call to the power method. For the case where we can converge to a good estimate for the maximum eigenvalue in two or three iterations, we can easily call the power method once every forty or fifty time steps and use the power method to improve the time step estimate in a cost-effective manner. Over a time span of forty or fifty time steps, matrix operator \(\mathbf{A}\), as defined in Equation ((3.6)) may change by only a small amount for a large class of problems. If matrix operator \(\mathbf{A}\) changes by a small amount, then using the last value of \(\mathbf{x}_{n}\) from a previous call to the power method provides an excellent starting vector and gives us quick convergence. This ability of the power method to converge quickly by using previous information makes it a viable scheme for calculating the maximum eigenvalue and producing a critical time step estimate.
For a problem with contact, the equations derived for the Lanczos method (with contact) in Section 3.2.1.1 are also applicable to the power method. For the power method and the Lanczos method, the added computational cost of the contact calculations results in reaching a break-even point with a smaller number of iterations when compared to a case where there is no contact.
3.2.2.2. Controls for Power Method
Using the power method is similar to using the Lanczos method in that we must be able to perform a power method calculation and reuse the power method estimate for the critical time step in some way over a number of subsequent time steps. The approach for reusing a power method estimate for the critical time step over a number of time steps is the same as the approach for reusing a Lanczos-based estimate over a number of time steps as described in Section 3.2.1.2. The approach discussed in Section 3.2.1.2 presents a way to reuse a maximum eigenvalue time step estimate over a number of time steps so that we maintain a critical time step estimate that is close to the theoretical maximum value in between the calls to the maximum eigenvalue calculations. The approach discussed in in Section 3.2.1.2 makes use of the element-based critical time step estimate at each time step.
3.2.2.3. Scale Factor for Power Method
When the power method is called for a given time step, it must appear that the calculations are using a constant matrix operator, \(\mathbf{A}\), for all iterations. The \(\mathbf{A}\) operator is the product of the inverse of the mass matrix \(\mathbf{M}^{-1}\) and the tangent stiffness matrix \(\mathbf{K}_{T}\) for the problem. The mass matrix is a constant. Therefore, we must assure that we are working with what is effectively a constant tangent stiffness matrix for all iterations during a call to the power method. This is the same requirement that we encounter in the Lanczos method (see Section 3.2.1.3). For the power method, we use the internal force calculations to generate the product \(\mathbf{A} \mathbf{x}_{i}\) by first computing the internal forces using \(\mathbf{x}_{i}\), which gives us the product \(\mathbf{K}_{T} \mathbf{x}_{i}\), and then multiplying the internal force results by the inverse of the mass matrix. Any vector \(\mathbf{x}_{i}\), as calculated by the power method, may be such that it represents large-strain behavior and moves the internal force calculations into a nonlinear regime. It is necessary to scale the \(\mathbf{x}_{i}\) vectors so that the internal force calculations are in a small-strain regime, which makes it appear that we are working with a constant tangent stiffness matrix. The vectors \(\mathbf{x}_{i}\) must be scaled so that they represent velocities associated with small-strain calculations. When properly scaled vectors are sent to the internal force calculation, the internal force calculation becomes a matrix \(\times\) vector product with a constant tangent stiffness matrix for all iterations during a given call to the power method.
The scale factor for the \(\mathbf{x}_{i}\) vectors, which we will designate as \(v_{sf}\), must not be too small, as this will create round-off problems and give a bad estimate for the critical time step. If the scale factor is too large, we violate the above restriction of a constant tangent stiffness matrix.
There are two approaches for controlling the scale factor when the power method is used to compute the maximum eigenvalue. These approaches are discussed in Section 3.2.2.5.
3.2.2.4. Accuracy of Eigenvalue Estimate
To use the power method efficiently, it is best not to specify a set number of iterations for every instance when the power method is called. Instead, some tolerance should be set on the convergence measure defined in (3.5). The initial call to the power method may take a significant number of iterations to produce a good estimate for the maximum eigenvalue. On subsequent calls to the power method, the number of iterations to converge to a good estimate for the maximum eigenvalue should be small when we use previous values of \(\mathbf{x}\) as a starting vector.
For the power method, we rely on the convergence measure defined by (3.5) to control the number of iterations. The user can set an upper limit on the number of iterations for the power method. Control of the number of iterations for the power method is discussed in more detail in Section 3.2.2.5.
3.2.2.5. Power Method Parameters Command Block
BEGIN POWER METHOD PARAMETERS <string>powermethod_name
STARTING VECTOR = <string>STRETCH_X|STRETCH_Y|STRETCH_Z|
ISOTHERMAL(ISOTHERMAL)
INCREASE OVER STEPS = <integer>incr_int(5)
NUMBER ITERATIONS = <integer>num_iter(150)
EIGENVALUE CONVERGENCE TOLERANCE =
<real>converge_tol(0.5e-3)
SMALL STRAIN = <real>small_strain(1.0e-3)
VECTOR SCALE = <real>vec_scale(1.0e-5)
SCALE FACTOR = <real>time_scale(0.9)
UPDATE ON TIME STEP CHANGE = <real>tstep_change(0.10)
UPDATE STEP INTERVAL = <integer>step_int(50)
END [POWER METHOD PARAMETERS <string>powermethod_name]
If you use the power method to compute a critical time step, there should be only one POWER METHOD PARAMETERS command block, and it should appear in the region. If you have a POWER METHOD PARAMETERS command block, you should not specify a command block for the node-based method or the Lanczos method. If you use the power method to compute the critical time step, the time step increase factor (default or user-specified) set in the TIME CONTROL command block will be used to control the increase in the time step estimate; the time step scale factor (default or user-specified) set in the TIME CONTROL command block will not be used. These factors are specified in the PARAMETERS FOR PRESTO REGION portion of a TIME STEPPING BLOCK in the TIME CONTROL command block.
The power method requires some type of starting vector. This is determined by the STARTING VECTOR command line. The various options available in this command line will generate a displacement vector that stretches your model in the \(x\)-, \(y\)-, or \(z\)-direction, or in all three directions (\(x\), \(y\), and \(z\)) at once. (STRETCH_X stretches the model in the \(x\)-direction, STRETCH_Y stretches the model in the \(y\)-direction, STRETCH_Z stretches the model in the \(z\)-direction, and ISOTHERMAL stretches the model in all three directions at once.) The displacement vector then serves as a basis for generating the starting vector \(x_{0}\) in the initialization phase of the Lanczos method. The default starting vector is the ISOTHERMAL option. You may encounter cases where use of one of the other options for a starting vector—STRETCH_X, STRETCH_Y, or STRETCH_Z—may result in a slight increase in accuracy for the critical time step estimate for a given number of iterations. These other options (STRETCH_X, STRETCH_Y, and STRETCH_Z) are offered for these special cases.
The INCREASE OVER STEPS command line determines how many steps, incr_int, are used to transition from an element-based critical time step estimate to the power method estimate at the beginning of an analysis. The user may want to increase from the element-based estimate to the power method estimate over a number of time steps if the difference between these two estimates is large. The value of incr_int defaults to 5.
As indicated in Section 3.2.2.5, some number of iterations are required by the power method to obtain a good estimate for the maximum eigenvalue for a given model. The number of iterations is controlled by the tolerance set for the convergence measure given in (3.5). The tolerance is set with the EIGENVALUE CONVERGENCE TOLERANCE command line. The default value for the convergence tolerance, converge_tol is \(0.5 \times 10^{-3}\). (This default value may change in the future as we gain more experience with the power method.)
The user can set an upper limit on the number of iterations for the power method by using the NUMBER ITERATIONS command line. If the number of iterations in the power method exceeds the value set by num_iter, then the power method is terminated. Otherwise, the iterative process in the power method is terminated by the tolerance set on the convergence measure. The default value for num_iter is 150.
As indicated in Section 3.2.2.3, it is necessary to scale the vectors \(\mathbf{x}_{i}\) (with a scale factor \(v_{sf}\)) in the power method calculations so that it appears we are always working with a constant tangent stiffness each time we call the power method. There are two approaches for determining the scale factor. With the first approach, the user can specify the value of a small strain for a particular model by using the SMALL STRAIN command line. (As long as the internal force calculation inside the power method is working with a vector that approximates small-strain behavior, it will appear that we are working with a constant tangent stiffness inside the power method.) The value for \(v_{sf}\) is then computed based on your specification for what constitutes a small strain in the model. (The default value for small_strain is \(1.0 \times 10^{-3}\).) With the second approach for determining the scale factor, the user can specify a value for \(v_{sf}\) directly by using the VECTOR SCALE command line. The default value for \(v_{sf}\), the vec_scale parameter in that command line, is \(1.0 \times 10^{-5}\). A number of tests have established that this value for the scale factor, \(1.0 \times 10^{-5}\), works well for a range of models encountered at Sandia National Laboratories. If the user wants to use the option of directly specifying the scale factor, it is possible to determine whether a particular scale factor is suitable for a particular problem. Take a scale factor \(v_{sf}\), plus values on either side of it, say, \(0.9 \times v_{sf}\) and \(1.1 \times v_{sf}\). If all three of these scale factor values produce almost the same estimate for the critical time step for a particular model, then the value for \(v_{sf}\) meets the criterion for an acceptable scale factor.
By default, the small-strain approach (with the default value of \(1.0 \times 10^{-3}\)) is used to determine a scale factor. If neither the SMALL STRAIN command line nor the VECTOR SCALE command line appears in the POWER METHOD PARAMETERS command block, then the default small-strain value will be used to calculate the scale factor. If the SMALL STRAIN command appears in the command block with any value other than the default value and the VECTOR SCALE command line does not appear in the command block, then the vector scale will be determined by the user-specified value for the small strain. If the VECTOR SCALE command line does appear in the POWER METHOD PARAMETERS command block, then the vector scale will be determined by the user-specified value for the vector scale.
The recommended way of determining the scale factor is to use the small-strain approach with the default small-strain value. Thus, if you follow this advice, neither the SMALL STRAIN command line nor the VECTOR SCALE command line should be included in the POWER METHOD PARAMETERS command block.
As indicated previously, the scale factor must not be too small, as this will create round-off problems and give a bad estimate for the critical time step. If the scale factor is too large, we violate the above restriction of a constant tangent stiffness matrix.
The SCALE FACTOR command line sets the factor \(f_{s}\) in (3.3). The value for \(f_{s}\), i.e., time_scale is set to 0.9. More experience with the power method will help us determine whether a slightly larger value of \(f_{s}\) can be used as a default.
The UPDATE ON TIME STEP CHANGE command line sets the value for \(\Delta t_{lim}\) in (3.4). If the change in the element-based critical time step estimate as given by (3.4) exceeds the value for \(\Delta t_{lim}\), then the power method is called for a new estimate for the critical time step. The default value for \(\Delta t_{lim}\), i.e., tstep_change, is 0.10.
The UPDATE STEP INTERVAL command line sets the number of step intervals at which the power method is called. If step_int is set to 40, the power method will be called every forty steps to compute an estimate for the critical time step. (The default value of step_int is 50.) This control mechanism interacts with the control established by the UPDATE ON TIME STEP CHANGE command line. Suppose we have set step_int to 40, and we have computed thirty-five steps since the last call to the power method. If \(\Delta t_{lim}\) has been set to 0.15 and we exceed this value at step 35, then the change in the element-based time step will result in the power method being called. The counter for keeping track of the number of step intervals since the last power method computation will be reset to zero. The next call to the power method will then occur again in forty steps, unless we again exceed the change in the element-based time step.
3.2.3. Node-Based Method
Now that we have developed schemes to make the Lanczos method and power method cost-effective tools for estimation of the critical time step, let us examine the node-based scheme.
The node-based method in Presto will give an estimate for the critical time step that is greater than or equal to the element-based estimate. It may or may not give an estimate for the time step that is close to the maximum value associated with the maximum eigenvalue for the problem. In general, we assume the node-based method will give us an estimate larger than the element-based estimate, but not significantly larger.
If the node-based scheme is used to determine the critical time step for a block of uniform elements (all the same material), then the estimate from the node-based method will be the same as that for the element-based estimate. The node-based estimate begins to diverge from (become larger than) the element-based estimate as the differences in aspect ratios of the elements attached to a node become larger. We can assume, therefore, that if the stiffest part of our structure has a relatively uniform mesh, the node-based method will not give a significantly larger estimate than the element-based method.
The node-based method costs only a fraction of an internal force calculation. However, the node-based method makes use of element time step estimates. Therefore, to use the node-based method, one must also do the element-based method. The cost of the node-based method is in addition to the element-based method. The modest increase in the critical step estimate from the node-based method is unlikely to offset the added cost of the node-based method.
To make the node-based method cost-effective, we can use the same methodology that is employed for the Lanczos method. Let \(t_{b}\) be the estimate for the critical time step from the node-based method. We define the ratio of the node-based estimate to the element-based estimate as
This ratio is then used to scale subsequent element-based estimates for the critical time step. If \(\Delta t_{e(n)}\) is the \(n^{th}\) element-based critical time step after the time step where the node-based calculations are performed, then the \(n^{th}\) time step after the node-based calculations, \(\Delta t_{(n)}\), is
The ratio \(t_{r}\) is used until the next call to the node-based method. The next call to the node-based method is controlled by one of two mechanisms described in the Lanczos discussion. The node-based method is called after a set number of times or a significant change in the element-based estimate for the critical time step.
3.2.3.1. Node-Based Parameters Command Block
BEGIN NODE BASED TIME STEP PARAMETERS <string>nbased_name
INCREMENT INTERVAL = <integer>incr_int(0)
STEP INTERVAL = <integer>step_int(1)
TIME STEP LIMIT = <real>step_lim(0.10)
END [NODE BASED TIME STEP PARAMETERS <string>nbased_name]
If you use the node-based method to compute a critical time step, there should be only one NODE BASED TIME STEP PARAMETERS command block, and it should appear in the region. If you use the node-based method to compute the critical time step, you should set the time step scale factor to 1.0 and the time step increase factor to a large number (5.0 is an acceptable value for the time step increase factor). These factors are specified in the PARAMETERS FOR PRESTO REGION portion of a TIME STEPPING BLOCK in the TIME CONTROL command block. If you have a NODE BASED TIME STEP PARAMETERS command block, you should not specify a LANCZOS PARAMETERS command block.
The INCREMENT INTERVAL command line determines how many steps, incr_int, are used to transition from an element-based critical time step estimate to the node-based estimate at the beginning of an analysis. The user may want to increase from the element-based estimate to the node-based estimate over a number of time steps if the difference between these two estimates is large. The value of incr_int defaults to 0 so the node based time step is used immediately at the beginning of the analysis.
The STEP INTERVAL command line sets the number of step intervals at which the node-based method is called. If step_int is set to 1000, the node-based method will be called every one thousand steps to compute an estimate for the critical time step. (The default value of step_int is 1.) Increasing the step interval can improve performance, but depending on the problem could run the risk of sending the time step super-critical. This control mechanism interacts with the control established by the TIME STEP LIMIT command line. Suppose we have set step_int to 1000, and we have computed 800 steps since the last call to the node-based method. If \(\Delta t_{lim}\) has been set to 0.15 and we exceed this value at step 800, then the change in the element-based time step will result in the node-based method being called. The counter for keeping track of the number of step intervals since the last node-based computation will be reset to zero. The next call to the node-based method will occur again in one thousand steps, unless we again exceed the change in the element-based time step.
The TIME STEP LIMIT command line sets the value for \(\Delta t_{lim}\),
This is exactly the same as (3.4) presented eleven pages earlier. If the change in the element-based critical time step estimate \(\Delta t_{lim}\), exceeds the value for \(\Delta t_{lim}\), then the node-based method is called for a new estimate for the critical time step. The default value for \(\Delta t_{lim}\), i.e., step_lim, is 0.10.