Full Equations (FEQ) Model for the Solution of the Full, Dynamic Equations of Motion for One-Dimensional Unsteady Flow in Open Channels and Through Control Structures

U.S. GEOLOGICAL SURVEY WATER-RESOURCES INVESTIGATIONS REPORT 96-4240

9.3 Stopping Criteria for Newton's Method

The Jacobian matrix is solved by use of a direct method based on a variant of Gaussian elimination, as described in the previous section, to compute the successive corrections to the estimated unknowns. Criteria are needed to determine when to stop the iteration. Selection of stopping criteria is difficult, as indicated by Hamming (1973, p. 68-70) and Dennis and Schnabel (1983, p. 159-161), because determining an acceptable difference between iterations is complicated. Under good conditions, the corrections from Newton's method should decrease rapidly. This usually results when the conditions for convergence of the method are met. In other cases, however, the corrections may not decrease or may decrease slowly, so alternative actions must be taken in applying FEQ to successfully complete the computations.

Two forms of convergence criteria, relative and absolute, can be used. A relative criterion involves the size of the correction or the size of the residual function relative to some other quantity. An absolute criterion directly involves the size of the correction or the residual function. No one criterion will work for all flow conditions. A combination of relative and absolute criteria is applied in FEQ simulation.

A relative criterion for closeness works well in FEQ simulation if the quantity tested is not small. If the quantity tested becomes too small, then the relative convergence criterion is supplanted by absolute criteria in FEQ simulation so that computational convergence is obtained. For example, the relative change in flow rate is meaningful so long as the flow does not become too small. Obviously, a flow of zero cannot be used to define a relative correction. The user must define what flow is too small because FEQ could be applied to simulate the lower Mississippi River or it could be applied to simulate a 5-ft wide brook. Thus, the user specifies a value called QSMALL that is added to the absolute value of the flow to yield the quantity defining the relative change in flow given the correction to the flow rate from Newton's method. If QSMALL is too small, and the user-specified relative change criterion, EPSSYS, is too small, then the criteria for stopping will not be met, and the time step will be reduced in FEQ simulation as described later in this section. QSMALL is usually some small fraction of the flow range of interest. For example, if the flows of interest are greater than 100 ft ³ /s, then a value of QSMALL of 0.1 ft ³ /s would not be appropriate because the flows are only known at best to within 5 percent. Use of a QSMALL value of 5 to 10 percent of the maximum flows of interest has worked well in typical applications of FEQ. This means that the relative change accepted for the low flows might be considerably larger than EPSSYS.

The stopping criterion for water-surface-height values at locations having a cross-sectional area is defined in terms of the relative change in area. This criterion works well if the depth and area are not too small. The stopping criteria for small flow areas are considered differently than for small flows. If the correction to water-surface height is less than an input value, ABSTOL, then the computations have converged at that point and the relative change is taken as zero. Otherwise, the relative change in the area is computed and must be less than the value of EPSSYS before convergence is achieved. Again, this means that the relative change in area accepted for small areas is larger than for large areas. The elevations at the flow-path end nodes on dummy branches are subject to the relative criterion of change in elevation divided by the current depth at the node.

Relative and absolute criteria also are used for level-pool reservoirs. The relative criterion is , where S_R is the storage volume in the level-pool reservoir. The absolute criterion is fixed internally at 0.001 ft. The change in storage is given by the product of the correction in water-surface elevation and the current surface area of the reservoir. The storage volume in the reservoir can become quite small. Thus, a lower limit on the value of storage used in calculating the relative convergence is set in FEQ. This lower limit is the larger of 1,000 ft ³ or the volume of 1 ft on the current surface area of the reservoir.

The criteria are applied to each unknown in the equation system. The rule for convergence for the equation system is that all unknowns must satisfy convergence criteria simultaneously. Thus, the maximum value of all relative changes must be less than EPSSYS before convergence is achieved. Experience has shown that, in some cases, only one or two variables will prevent convergence, whereas relative changes in all other variables are only a small fraction of EPSSYS. This is often the case when convergence is slow and the time step is being reduced to solve some computational problem. To increase the robustness of the solution scheme, another convergence criterion is added. Convergence may be declared with this criterion if no more than a user-input number, NUMGT, of unknowns does not meet the current convergence criteria but do meet a less restrictive, secondary convergence criterion. Thus, EPSSYS and ABSTOL together become the primary convergence criteria. NUMGT and the secondary relative change criterion become the secondary criteria for convergence. The number of unknowns not meeting the primary convergence criteria is output at each iteration in model simulation.

In most cases, all unknowns meet the convergence criteria, but disregarding local convergence problems greatly increases the robustness of the unsteady-flow analysis when flows and depths are small. Because the affected flows and depths are small relative to the flows and depths of primary interest, the effect of these additional convergence criteria on the final results is small. Results from several tests have shown that the maximum flows and stages obtained are usually the same or are within the uncertainty implied by the convergence criteria.

If the convergence criteria are not met within a user-supplied maximum number of attempts, the time step is reduced with a user-supplied factor, the integration weight, W_T, is incremented, again with a user-supplied factor, and a solution is computed again. This process continues until the time step becomes too small to continue or until the convergence criteria are met. The minimum time step is specified by the user. A weighted average of the number of iterations to convergence is maintained, and the average is used to increase or decrease the time step. A maximum time step, a minimum time step, and several weights to control the time step are specified by the user. Thus, the time step is under model control during the computations and will vary as the rate of convergence to a solution varies throughout a simulation.

The value used for EPSSYS must reflect the estimation of the stopping criteria based on the results of the most recent correction. All computed corrections are used. Because Newton's method converges rapidly when the root is approached, a 5-percent relative correction made to an unknown value often implies a relative correction of one-tenth or less of that amount in the next correction. Experience with FEQ simulation has supported this conclusion. A reduction of EPSSYS by a factor of 5, from 0.05 to 0.01, generally makes little difference in the results because most of the time steps have converged to a relative correction of 0.01 or less. Run time increases of about 30 percent have been typical with reduction in the system relative tolerance by a factor of 5 or more. Moreover, the stopping criteria are applied to the maximum value anywhere in the model. Commonly, most of the unknowns have relative corrections of a factor of 10 or more smaller than the maximum relative correction in the system.