4.2. Conjugate Gradient Solver

The core nonlinear preconditioned conjugate gradient solver is controlled through the CG command block. This must be nested within a SOLVER command block, regardless of whether multilevel controls are used or not.

BEGIN CG
  #
  # convergence commands
  # Default is not defined.
  TARGET RESIDUAL = <real>target_resid
    [DURING <string list>period_names]
  TARGET RELATIVE RESIDUAL = <real>target_rel_resid(1.0e-4)
    [DURING <string list>period_names]
  # Default is 10 times target residual
  # or 1.0e+20 if a multilevel solver is defined
  # and control failure is not defined.
  ACCEPTABLE RESIDUAL = <real>accept_resid
    [DURING <string list>period_names]
  # Default is 10 times target relative residual
  # or 1.0e+20 if a multilevel solver is defined
  # and control failure is not defined.
  ACCEPTABLE RELATIVE RESIDUAL = <real>accept_rel_resid
    [DURING <string list>period_names]
  # Default is not defined.
  REFERENCE = EXTERNAL|INTERNAL|BELYTSCHKO|RESIDUAL|ENERGY
    (EXTERNAL) [DURING <string list>period_names]
  MINIMUM REFERENCE = NONE|FRACTION_OF_INITIAL(NONE) <real>min_ref(0.0)
    [DURING <string list>period_names]
  RESIDUAL NORM TYPE = ALL|TRANSLATION|SCALE_RB_ROTATIONS
    (ALL) [DURING <string list>period_names]
  MINIMUM RESIDUAL IMPROVEMENT = <real>resid_improvement
    [DURING <string list>period_names]
  MINIMUM ITERATIONS = <integer>min_iter(0)
    [DURING <string list>period_names]
  RESIDUAL ROUNDOFF TOLERANCE = <real>(1.0e-15)
  #
  # default = 100 for tangent preconditioner.
  # default = max(NumNodes,1000) for all other preconditioners.
  MAXIMUM ITERATIONS = <integer>max_iter
    [DURING <string list>period_names]
  #
  # preconditioner commands
  PRECONDITIONER = BLOCK|BLOCK_INITIAL|DIAGONAL|
    ELASTIC|IDENTITY|PROBE|TANGENT(ELASTIC) [<real>scaling_factor]
  PRECONDITIONER ITERATION UPDATE = <integer>iter_update
  BALANCE PROBE = <integer>balance_probe(1)
    [DURING <string list>period_names]
  NODAL PROBE FACTOR = <real>probe_factor(1.0e-6)
    [DURING <string list>period_names]
  NODAL DIAGONAL SCALE = <real>nodal_diag_scale(0.0)
    [DURING <string list>period_names]
  NODAL DIAGONAL SHIFT = <real>nodal_diag_shift(0.0)
    [DURING <string list>period_names]
  BEGIN FULL TANGENT PRECONDITIONER
    # ... see Section 4.3
  END [FULL TANGENT PRECONDITIONER]
  #
  # line search command, default is secant
  LINE SEARCH ACTUAL|TANGENT [DURING <string list>period_names]
  LINE SEARCH SECANT [<real>scale_factor]
  LINE SEARCH BACKTRACK [<real>min_step(1.0e-3)
                         <real>reduction(0.95)
                         <real>init_step(1.0)
                         <enum>NONE|SECANT(NONE)
                         <real>dispinc(infinity)]
  LINE SEARCH TYPE MULTI = <real>num_domain(2)
    [DURING <string list>period_names]
  #
  # diagnostic output commands
  ITERATION PRINT = <integer>iprint(25)
    [DURING <string list>period_names]
  # off by default.
  ITERATION PLOT = <integer>iplot
    [DURING <string list>period_names]
  ITERATION PLOT OUTPUT BLOCKS = <string list>plot_blocks
  #
  # cg algorithm commands
  ITERATION RESET = <integer>iter_reset(infinity)
    [DURING <string list>period_names]
  ORTHOGONALITY MEASURE FOR RESET = <real>ortho_reset(0.5)
    [DURING <string list>period_names]
  RESET LIMITS <integer>iter_start <integer>iter_reset
    <real>reset_growth <real>reset_orthogonality
    [DURING <string list>period_names]
  BETA METHOD = FletcherReeves|PolakRibiere|Newton|
    PolakRibierePlus(PolakRibiere)
    [DURING <string list>period_names]
  # line search step length bounds
  MINIMUM STEP LENGTH = <real>min_step(-infinity)
    [DURING <string list>period_names]
  MAXIMUM STEP LENGTH = <real>max_step(infinity)
    [DURING <string list>period_names]
END [CG]

Section 4.2.1 through Section 4.2.4 describe the components of the CG command block.

4.2.1. Convergence Commands

# Default is not defined.
TARGET RESIDUAL = <real>target_resid
  [DURING <string list>period_names]
TARGET RELATIVE RESIDUAL = <real>target_rel_resid(1.0e-4)
  [DURING <string list>period_names]
ACCEPTABLE RESIDUAL = <real>accept_resid(1.0e-3)
  [DURING <string list>period_names]
# Default is 10 times target relative residual
# or 1e20 if a multilevel solver is defined
# and control failure is not defined.
ACCEPTABLE RELATIVE RESIDUAL = <real>accept_rel_resid
  [DURING <string list>period_names]
# Default is not defined.
REFERENCE = EXTERNAL|INTERNAL|BELYTSCHKO|RESIDUAL|ENERGY(EXTERNAL)
  [DURING <string list>period_names]
MINIMUM REFERENCE = NONE|FRACTION_OF_INITIAL(NONE) <real>min_ref(0.0)
  [DURING <string list>period_names]
RESIDUAL NORM TYPE = ALL|TRANSLATION|SCALE_RB_ROTATIONS(ALL)
  [DURING <string list>period_names]
MINIMUM RESIDUAL IMPROVEMENT = <real>resid_improvement
  [DURING <string list>period_names]
MINIMUM ITERATIONS = <integer>min_iter(0)
  [DURING <string list>period_names]
# default = 100 for tangent preconditioner
# default = max(NumNodes,1000) for all other preconditioners.
MAXIMUM ITERATIONS = <integer>max_iter
  [DURING <string list>period_names]
RESIDUAL ROUNDOFF TOLERANCE = <real>1.0e-15)
  [DURING <string list>period_names]

The nonlinear preconditioned CG solver iterates to decrease the residual until the residual is deemed to be sufficiently small, based on user-specified convergence criteria. The command lines listed above are placed in the CG command block and used to control these convergence criteria.

Solver convergence is measured by computing the \(L^2\) norm of the residual and comparing that residual norm with target convergence criteria specified by the user. Convergence can be monitored either directly in terms of the residual norm, or in terms of a relative residual, which is the residual norm divided by a reference quantity that is indicative of the current loading conditions on a model. Basing convergence on the relative residual is often helpful because doing so ensures that the convergence target is meaningful for the model. There are some situations, however, when it is better to check for convergence based on the absolute residual rather than on the relative residual. Checking the absolute residual is especially helpful when the model passes through a stage during which there are no loads. The actual residual and the relative residual can be used for checking convergence at the same time.

All the convergence command lines in the CG command block can have different values for different time periods by using the DURING specification. This specification consists of the key word DURING followed by a list of time periods and is appended to the end of the command line. Each of the time periods included in the list must correspond to the name of a TIME STEPPING command block (i.e., time_block_name) specified in the TIME CONTROL command block. If the DURING specification is omitted from the command line, the value of the parameter in the command line is applicable for the entire analysis. This value can be overwritten for specific periods by repeating the same command line and using the DURING specification. In other words, multiple occurrences of the same command line can be included in the command block, where one without a DURING specification applies to many time periods and others with unique DURING specifications apply to specific time periods.

The TARGET RESIDUAL command line specifies the target convergence criterion in terms of the actual residual norm. The TARGET RELATIVE RESIDUAL command line specifies a convergence criterion in terms of the relative residual. If both command lines are included in the CG command block, the solver will accept the solution as converged if either the target residual or the target relative residual is below the specified values. Note that use of the term “target” in this discussion means “desired.”

The solver also allows the user to define acceptable convergence criteria. If the solution has not converged to the specified targets before the maximum number of iterations is reached, the residual is checked against the acceptable convergence criteria. These criteria are specified via the ACCEPTABLE RESIDUAL and ACCEPTABLE RELATIVE RESIDUAL command lines. The concepts of acceptable residual and acceptable relative residual are the same as those used for the target limits, i.e., in the TARGET RESIDUAL and TARGET RELATIVE RESIDUAL command lines. If the solution has not met the target criteria but has met the acceptable criteria, the solution is allowed to proceed. The default value for each acceptable criterion command line is 10 times the value of its corresponding target criterion command line, e.g., if the value of target_resid was 1.0e-6, the value of accept_resid would be 1.0e-5.

When a multilevel solver is being used, the acceptable residual and acceptable relative residual are set to \(1e^{20}\) as long as control failure is not being used. This is because the multilevel controls ultimately check and determine convergence, and allowing the multilevel controls to continue iterating even when a model problem fails to converge can allow the overall problem to converge.

Solutions that meet only the acceptable criteria are noted in the log file. If the acceptable criteria are not met, the solution at that load step has failed to converge, and the solver exits. If adaptive time stepping, as discussed in Section 4.10.2, is active, a solution of the load step may be attempted with a smaller time step. Otherwise, the code will exit at this point with an error.

If relative residuals are given in the convergence criteria, the REFERENCE command line can be used to select the method for computing the reference load. This command line has several options: EXTERNAL, INTERNAL, BELYTSCHKO, RESIDUAL and ENERGY. When using the default option (EXTERNAL), the reference value is taken as the maximum of the \(L^2\) norm of either the current externally applied or reaction loads. If the model has force boundary conditions, the external force vector used for this norm is composed of the nodal forces resulting from those boundary conditions. If there are no prescribed forces, the reaction forces at all prescribed kinematic boundary conditions are used instead. The INTERNAL option uses the norm of the internal forces as the reference load. This option is helpful in situations where a structure has no externally applied boundary conditions, such as in the case of a thermally loaded structure. The BELYTSCHKO option uses the maximum of the norm of the external load, the norm of the reactions, the norm of the internal forces, and the norm of the inertial forces. The RESIDUAL option denotes that the residual from the initial residual should be used as the reference quantity. The initial residual is computed from the first iteration of the solver.

The ENERGY option to the REFERENCE command line is a modification of the EXTERNAL option where energy conjugate norms of both the residual and reference are used. The energy conjugate norm of the residual is given by \(R \cdot K^{-1} \cdot R\) where \(K\) is the diagonal of the stiffness matrix. For the energy conjugate norm of the external forces, \(R\) is replaced by the external nodal forces. The primary advantage of the energy norm is that forces and moments are converted to consistent and compatible energy units.

The MINIMUM REFERENCE offers the ability to set a value that the reference cannot drop below. Using the option NONE will use the value of min_ref for the reference if the computed reference is smaller. The option FRACTION_OF_INITIAL will compute a minimum reference value that is a fraction min_ref of the reference at the first iteration of the loadstep. A M will appear by the reference value in the log file if the minimum reference was used. If the line command is not used the minimum reference is zero. Having a minimum reference can make the use of relative tolerances still possible in the cases that the reference value disappears as the problem is solved.

Convergence can be difficult in some problems with rigid bodies without any applied moments. In these cases, it has been found advantageous to scale the rotational residuals on the rigid bodies by an approximation of the moment arm of the rigid bodies about the reference point. This can be accomplished through use of the RESIDUAL NORM TYPE command with the keyword SCALE_RB_ROTATIONS. The other options are ALL and TRANSLATION. A value of ALL is the default and treats all values the same in both the residual norm and the reference norm. A value of TRANSLATION uses only the translational components of residual and forces in the norm computations. Though this command is also available for other levels in a multilevel solver, this command is only appropriate for the core CG solver of a Sierra/SM region and should not be used in either RVE or Control Modes regions.

The default ALL option means that the residual includes all translational and rotational degrees of freedom. Note, translational forces and rotational moments have different units. Thus the ALL option may artificially increase or decrease the contribution of rotational degrees of freedom to the total residual depending on the units used. For example a model with a length unit in millimeters may tend to increase the contribution of rotational degrees of freedom to the convergence check (due to the numerically large moment arm) whereas a model with a length unit of meters may tend to decrease the contribution of rotational degrees of freedom to the convergence check (due to a numerically small moment arm.) The ALL option is mathematically inconsistent, but remains the default for historical compatibility. Also in practice the inconsistency is generally irrelevant as translational and rotational degrees of freedom tend to converge at the same rate so that convergence of both quantities is obtained simultaneously.

The TRANSLATION option to the RESIDUAL NORM TYPE command computes the residual norm over only translational degrees of freedom and ignores all rotational components. This option may be useful if the model is failing to converge rotational degrees of freedom as may occur on problems with very thick shells. The SCALE_RB_ROTATIONS option computes the residual over translational degrees of freedom plus scaled rigid body rotational degrees of freedom. The rigid body rotational degrees of freedom are divided by the rigid body characteristic size to produce consistent units for a unified residual norm. The SCALE_RB_ROTATIONS option excludes all other rotational degrees of freedom, such as those on shells, from the residual norm. The SCALE_RB_ROTATIONS option is recommended if a model contains large rigid bodies with applied external moment boundary conditions. Also note another option to avoid any issues encountered by residual unit incompatibility is to use the ENERGY option for the solver REFERENCE command.

Warning

The RESIDUAL NORM TYPE command should not be used in RVE or Control Modes regions.

The MINIMUM RESIDUAL IMPROVEMENT command line stipulates that the CG solver must reduce the initial residual by a specified amount to obtain convergence. If this command line is included in the CG command block, the improvement condition must be met in addition to the standard target criteria. The parameter resid_improvement is a real number between 0.0 and 1.0 that specifies the amount by which the residual must be improved as a fraction of the initial residual in the current model problem. Thus, to stipulate that the residual must be smaller than 10% of the initial residual, one would set the value of resid_improvement equal to 0.9. The MINIMUM RESIDUAL IMPROVEMENT command line is primarily useful in the context of the multilevel solver. The model problems presented to the core solver sometimes begin with low residuals and require few iterations to converge. This command line forces the core solver to further improve the residual, which can accelerate the convergence of the other solver levels.

The MAXIMUM ITERATIONS and MINIMUM ITERATIONS command lines are used to specify the maximum and minimum number of core nonlinear solver iterations, respectively. These limits apply for a load step if the CG solver is used in the stand-alone mode or for a model problem if the CG solver is used as the core solver in the multilevel solver. The default maximum number of iterations is set to \(max(N_{nodes},1000)\) where \(N_{nodes}\) is the global number of nodes in the mesh. The linear CG algorithm is guaranteed to converge to the solution in a number of iterations equal to the number of equations in the matrix. However, since we are solving a nonlinear system of equations using a nonlinear CG algorithm, the number of iterations to convergence can be larger than the number of equations . The default value for the minimum number of iterations, min_iter, is zero. If a number greater than zero is specified, the solver will iterate at least that many times, regardless of whether the convergence criteria are met.

The RESIDUAL ROUNDOFF TOLERANCE command line specifies a tolerance used to compute an approximately zero residual. Due to finite precision arithmetic, the solver can only drop the residual so far before additional corrections to the nodal coordinate field cannot be numerically represented. The approximately zero residual is the residual that would be produced by a perturbation of the displacement by the smallest machine representable perturbation. The default value for this command is the accuracy of double precision arithmetic (1.0e-15). This tolerance number can be set to zero to disable this feature or set higher to converge to larger, but still near machine limit, residuals.

The log file annotates which convergence criterion was met in each CG iteration using markings shown in Table Table 4.1.

Table 4.1 Log File Convergence Markings

Mark

Convergence Property

>M

Convergence is reached, but still have not iterated the minimum number of specified iterations

>A

Reached maximum iterations and still not converged

<T

Residual converged to target tolerance of either relative or absolute criteria

<A

Reached maximum iterations and converged to acceptable tolerance of either relative or absolute criteria

0

Residual is converged due to being approximately zero to within machine precision

4.2.2. Preconditioner Commands

PRECONDITIONER = BLOCK|BLOCK_INITIAL|DIAGONAL|
  ELASTIC|IDENTITY|PROBE|TANGENT|(ELASTIC) [<real>scaling_factor]
PRECONDITIONER ITERATION UPDATE = <integer>iter_update
  [DURING <string list>period_names]
BALANCE PROBE = <integer>balance_probe(1)
  [DURING <string list>period_names]
NODAL PROBE FACTOR = <real>probe_factor(1.0e-6)
  [DURING <string list>period_names]
NODAL DIAGONAL SCALE = <real>nodal_diag_scale(0.0)
  [DURING <string list>period_names]
NODAL DIAGONAL SHIFT = <real>nodal_diag_shift(0.0)
  [DURING <string list>period_names]
BEGIN FULL TANGENT PRECONDITIONER
  # ... see :numref:`implicit-preconditioner`
END [FULL TANGENT PRECONDITIONER]

These command lines are used to select the type of preconditioner used by the CG algorithm, as well as the method used to form the preconditioner for certain types of preconditioners. The preconditioner is applied to the residual vector to obtain a gradient direction, which in turn is used to find a search direction. The most effective preconditioner is one that closely approximates the effect of multiplying by the inverse of the tangent stiffness matrix.

There are two basic types of preconditioners available: nodal and full tangent. The nodal preconditioners provide only the diagonal terms of the full tangent stiffness matrix or 3-by-3-block diagonal matrices along the diagonal of the stiffness matrix. The diagonal nature of a nodal preconditioner allows it to be inverted inexpensively and also use very little memory. Nodal preconditioners are so termed because they only account for the stiffness at each node in isolation and ignore the coupling between nodes that is accounted for in the off-diagonal terms of a full stiffness matrix. The result of this approximation is that with a nodal preconditioner, in a single CG iteration, the equilibration of a residual at a node can only cause movement in nodes that are directly connected by an element to that node. Thus, as the model size increases, the number of iterations also increases. Nodal preconditioners require many iterations, but they are often efficient because the iterations are inexpensive, especially if the problem is “blocky” in nature.

The PRECONDITIONER command line allows the user to select a nodal preconditioner or the full tangent preconditioner with most defaults. As defined below, this command line has several options. Some options are synonyms for other options.

  • The BLOCK and ELASTIC options, which are synonyms, specify the default preconditioner. These options provide a nodal preconditioner with 3-by-3-block matrices that are computed based on the elastic material properties. These matrices are updated at every model problem to account for the current geometry.

  • The BLOCK_INITIAL option, which is a variant of the BLOCK preconditioner, forms the preconditioner at the beginning of the analysis but never recomputes it for efficiency.

  • The DIAGONAL preconditioner is formed in the same way as the BLOCK preconditioner, but only uses the terms on the diagonal.

  • The PROBE option forms a 3-by-3-block diagonal preconditioner by probing the stiffness of nodes rather than by using the elastic stiffness.

  • The IDENTITY option uses an identity matrix for the preconditioner, meaning that the residual is used as the gradient direction. This option is only of academic interest and is included for completeness.

  • The TANGENT option is a quick way of specifying the FULL TANGENT PRECONDITIONER with most defaults used. A small amount (1.0e-10) of tangentdiagonal scaling is assigned by default when using PRECONDITIONER = TANGENT.See Section 4.3 for a further description of the full tangent preconditioner.

The PRECONDITIONER ITERATION UPDATE command line controls the frequency at which the full tangent and nodal probe preconditioner are updated. By default, the preconditioner is only updated at the first iteration of each model problem. If this command line is used, updates will occur at the specified frequency iter_update. If a model experiences highly nonlinear behavior during a model problem, it may be useful to use this command to cause more frequent preconditioner updates. Such problems may converge better when iter_update is set to about 10 for the full tangent case. Realize that frequent updates may reduce the number of CG iterations but increase the computational cost.

The BALANCE PROBE and NODAL PROBE FACTOR command lines control the behavior of the probing algorithm that is used to obtain the stiffness for the probe preconditioner (selected via the PROBE option in the PRECONDITIONER command line). For the probe preconditioner, the tangent stiffness \(K\),

\[K_{ij} = \frac{\partial F^{int}_i}{\partial x_j},\]

where \(F^{int}\) is the internal force and \(x\) is the displacement, is approximated using finite differences. The nodal probe preconditioner only computes the 3-by-3-block diagonal version of the stiffness matrix. The NODAL PROBE FACTOR command line controls the probing distance, with the value of probe_factor defaulting to 1.0e-6. The probe factor controls the probing distance relative to element size. Smaller values may give a better tangent approximation, but these values could result in the generation of round-off errors.

The BALANCE PROBE command line selects the type of finite differencing scheme used to probe for the stiffness. The value of balance_probe can be set to 0, 1, or 2, as explained below.

  • Setting balance_probe to 0, the default for the nodal probe preconditioner, causes the preconditioner to be formed by forward finite differencing. Probing occurs in only one direction, resulting in first-order error in the preconditioner. The stiffness computed when balance_probe is set to 0 appears as

    \[K_{ij} = \frac{F^{int}_i(x+\delta e_j) - F^{int}_i(x)}{\delta},\]

    where \(\delta\) is the probing distance, controlled by the NODAL PROBE FACTOR command line, and \(e_j\) is a unit vector in the \(j\)th equation direction.

  • Setting balance_probe to 1, the default for the full tangent preconditioner, causes the preconditioner to be formed with central differencing. Probing at each element degree of freedom is performed in both positive and negative directions, leading to an estimate of the tangent stiffness with second-order error. This approach can be necessary in problems with material nonlinearity, where a probe that stretches rather than compresses the element can result in orders of magnitude differences in estimated stiffness. The stiffness computed when balance_probe is set to 1 appears as

    \[K_{ij} = \frac{F^{int}_i(x+\delta e_j) - F^{int}_i(x-\delta e_j)}{2\delta}.\]

  • The value of balance_probe may also be set to 2, which allows the full tangent preconditioner to obtain central finite differencing with fourth-order error. This approach requires twice the work as central differencing but is more accurate. The stiffness computed when balance_probe is set to 2 appears as

    \[K_{ij} = \frac{-F^{int}_i(x+2\delta e_j) + 8F^{int}(x+\delta e_j) -8F^{int}(x-\delta e_j) + F^{int}_i(x-2\delta e_j)}{12\delta}.\]

The NODAL DIAGONAL SCALE command line specifies a nodal diagonal scale factor, which is used to scale the diagonal entries of the nodal preconditioner. This is useful in some situations where a model may be softening due to material effects. The diagonal terms of the nodal preconditioner are scaled by \((1.0 + \texttt{nodal\_diagonal\_scale})\). The default value of nodal_diagonal_scale is 0.0, which means the diagonals are not scaled. A value of 0.1 scales up all the diagonal terms by 10 percent.

The NODAL DIAGONAL SHIFT command line specifies a nodal diagonal shift, which is summed into the diagonal entries of the nodal preconditioner. This is useful in some situations where a model may be softening due to material effects. The diagonal terms of the nodal preconditioner are shifted by \(( \texttt{num\_nodes} * \texttt{tangent\_diagonal\_shift})\). The default value of nodal_diagonal_shift is 0.0, which means the diagonals are not shifted.

The FULL TANGENT PRECONDITIONER command block enables the full tangent preconditioner and specifies options that apply specifically to that type of preconditioner. See Section 4.3 for a description of the options available. The FULL TANGENT PRECONDITIONER command block cannot be used with the PRECONDITIONER command line in a CG command block. The user must select only one type of preconditioner for an analysis.

4.2.3. Line Search Command

LINE SEARCH ACTUAL|TANGENT [DURING <string list>period_names]
LINE SEARCH SECANT [<real>scale_factor]
LINE SEARCH BACKTRACK [<real>min_step(1.0e-3)
                       <real>reduction(0.95)
                       <real>init_step(1.0)
                       <enum>NONE|SECANT(NONE)
                       <real>dispinc(infinity)]
LINE SEARCH TYPE MULTI = <real>num_domain(2)
  [DURING <string list>period_names]

During each CG iteration, after the search direction is computed, a line search is used to find an optimal scaling factor that when applied to the search vector will result in a minimized residual. The line search algorithm is controlled with the LINE SEARCH command line. The SECANT option is the default, and is used if the LINE SEARCH command line is not present.

  • The ACTUAL line search is a single-step quadratic line search that is equivalent to the corresponding JAS3D line search option. The line search requires an additional evaluation for each iteration of the element internal forces and thus the material response. For a given search direction \(s\), residual \(r\), and velocity \(v\), the step length \(\alpha\) is computed as

    \[\alpha = - \frac{s^T r(v)}{s^T(F_{int}(v+s) - F_{int}(v))},\]

    where the residual is computed based on the internal force \(F_{int}\) and on the external force \(F_{ext}\) as \(r = F_{int}(v) - F_{ext}(v)\).

  • The default SECANT line search is a slight variation of the ACTUAL line search described above. Like the ACTUAL line search, this line search requires an additional internal force evaluation for each iteration. The step length \(\alpha\) is computed as

    \[\alpha = - \frac{s^T r(v)}{s^T(r(v+s) - r(v))}.\]

    The optional scale_factor can be used with the secant line search to scale the search vector \(s\) in the equation above. If scale_factor is omitted, no scaling is done. If scale_factor is provided on the command line, \(s\) is scaled by the product of scale_factor and the dimension of the smallest element. Setting the scale factor to a small value can be helpful to avoid inverting elements during the evaluation of \(r(v+s)\). The SECANT line search is recommended for problems in which the external force is highly nonlinear (e.g., a beam-bending problem with a pressure distribution).

  • The TANGENT line search avoids the additional internal force calculation that is inherent in the other line search types by using an approximate tangent modulus computed in the last call to the material subroutine for each element. The step length \(\alpha\) is computed as

    \[\alpha = - \frac{s^T r(v)}{s^TKs},\]

where \(K\) is the approximate tangent stiffness.

Warning

The TANGENT line search capability is not supported for any element formulation besides uniform gradient hex.

  • The BACKTRACK line search uses a cutback search to find the optimal step size. The solver will first try the init_step step size. If that step size does not meet the desired reduction it will continue searching by halving the step size until the desired reduction in residual is met, 20 iterations, or the min_stepsize is reached. The optimal step size is selected from all those attempted. The BACKTRACK line search will tend to be more expensive than other line search methods but has potential to be more robust in some circumstances such as problems with very sharp nonlinearity. The optional SECANT_INIT will add the secant step size to the step sizes attempted and will prefer it if successful. The dispinc will further limit the step size so that the displacement at any node is less than the value specified.

  • The MULTI line search first subdivides the line search space into multiple domains and then applies the secant line search method on each domain. The best line search scaling factor found in any sub domain is ultimately used. The multi line search option may in some circumstances provide a more optimal or robust line search value than the secant, though at a significantly increased cost per iteration. The number of subdomains used is defined by the sub_domain option which defaults to two.

4.2.4. CG Algorithm Commands

ITERATION RESET = <integer>iter_reset(infinity)
  [DURING <string list>period_names]
ORTHOGONALITY MEASURE FOR RESET = <real>ortho_reset(0.5)
  [DURING <string list>period_names]
RESET LIMITS <integer>iter_start <integer>iter_reset
  <real>reset_growth <real>reset_orthogonality
  [DURING <string list>period_names]
BETA METHOD = FletcherReeves|PolakRibiere|Newton|
  PolakRibierePlus(PolakRibiere)
  [DURING <string list>period_names]
MINIMUM STEP LENGTH = <real>min_step(-infinity)
  [DURING <string list>period_names]
MAXIMUM STEP LENGTH = <real>max_step(infinity)
  [DURING <string list>period_names]

The behavior of the CG algorithm can be changed using the command lines listed above. The CG algorithm orthogonalizes the current search direction against the previous search direction at each iteration. During an iteration where a reset occurs, the orthogonalization is skipped, resulting in a search in the steepest descent direction. These command lines can be appended with the DURING specification, as discussed in Section 4.2.1.

The ITERATION RESET command line specifies that the CG algorithm be reset after every iter_reset iterations. The default value is set equal to the maximum integer, which defaults the conjugate gradient algorithm to not reset.

For the CG algorithm to work well, each new search direction is constructed orthogonal to the previous search directions. If the current search direction has a significant component in one of the previous search directions, error can be introduced that will not be corrected in future iterations, resulting in a lack of convergence. At every iteration, the relative lack of orthogonality between the current gradient direction and the previous search direction is computed. When these become far from orthogonal, a reset can be performed to ignore the previous search directions and use the current gradient direction as the next search direction, which is equivalent to setting the orthogonalization parameter \(\beta\) to zero for that iteration (see (4.1)). The ORTHOGONALITY MEASURE FOR RESET command line specifies that a reset will occur if the lack of orthogonality exceeds the value of ortho_reset. The default value of 0.5 results in few resets. Setting ortho_reset to a value as tight as 0.01 can result in many resets and potentially improved convergence on some problems.

The RESET LIMITS command line is also used to control CG resets. This command line is provided to achieve compatibility with JAS3D, and it behaves exactly the same as its JAS3D counterpart. The command line, however, should not be used when either the ITERATION RESET command line or the ORTHOGONALITY MEASURE FOR RESET command line is included in the CG command block. Up to four parameters can optionally be specified in the RESET LIMITS command line. The first, iter_start, sets the number of iterations to wait before looking for a minimum residual. The second, iter_reset, sets the number of iterations to allow between finding a minimum and restarting the iterative algorithm. The third, reset_growth, sets the amount of growth in the residual norm that would indicate divergence and thus trigger a reset. Finally, reset_orthogonality sets the relative lack of orthogonality between the current gradient direction and the previous search direction that would trigger a reset.

The BETA METHOD command line allows different formulas to be used in computing the orthogonalization parameter \(\beta\) in the CG algorithm,

(4.1)\[s_k = g_k + \beta s_{k-1},\]

where \(s_k\) is the current iteration’s search direction, \(s_{k-1}\) is the previous iteration’s search direction, and \(g_k\) is the current iteration’s gradient direction.

In the below equations, \(r_k\) is the current iteration’s residual vector, \(r_{k-1}\) is the previous iteration’s residual vector, and \(g_{k-1}\) is the previous iteration’s gradient vector. The PolakRibierePlus formula forces beta to be positive. A steepest descent direction is taken when beta becomes negative.

  • For the FletcherReeves option the \(\beta\) scalar is computed as

    \[\beta_k = \frac{r^T_k g_k}{r^T_{k-1} g_{k-1}}.\]

  • For the PolakRibiere option (the default) \(\beta\) is computed as

    \[\beta_k = \frac{r^T_k (g_k - g_{k-1})}{r^T_{k-1} g_{k-1}}.\]

  • For the PolakRibierePlus option \(\beta\) is computed as

    \[\beta_k = {\rm max}(\beta^{PR}_k, 0).\]

  • For the Newton option the nonlinear conjugate gradient algorithm is transformed into a Quasi-Newton method without a line search. In order to best approximate a Newton solve, this method should be used with a full tangent preconditioner that is updated on every iteration.

The MINIMUM STEP LENGTH and MAXIMUM STEP LENGTH commands set the minimum step length and maximum step length that the CG algorithm can take along a search direction. The default values are set so that any positive or negative step length will be accepted. This option can be used to limit motion of the model along any search direction and can help improve convergence sometimes.

4.2.5. Diagnostic Output Commands

ITERATION PRINT = <integer>iter_print
  [DURING <string list>period_names]
ITERATION PLOT = <integer>iter_plot
  [DURING <string list>period_names]
ITERATION PLOT OUTPUT BLOCKS  = <string list>plot_blocks

The command lines listed above can be used to control the output of diagnostic information from the CG solver. The ITERATION PRINT and ITERATION PLOT command lines can be appended with the DURING specification, as discussed in Section 4.2.1.

The ITERATION PRINT command line controls the frequency at which the convergence information is printed to the log file. The default value of iter_print is 25 if a nodal preconditioner is used, and it is 1 if a full tangent preconditioner is used.

For a description of ITERATION PLOT and ITERATION PLOT OUTPUT BLOCKS, see Section 4.1.1.

4.2.6. CG Input Example

begin solver
  begin cg
    target relative residual = 1.0e-07
    Maximum Iterations = 5000
    preconditioner iteration update = 100
    preconditioner = elastic
  end
end