APPENDIX A: DIAGNOSING COMPUTATIONAL PROBLEMS IN FEQ

SECTION FOUR

4. Failure during Unsteady-Flow Computations

Once the steady-flow computations are complete and the final checks are made, FEQ simulations may fail during unsteady-flow computations. Some failures are abrupt; that is, the computations appear to be going smoothly, and suddenly, with no prior evidence of problems, FEQ issues an error message and stops. Other failures are more subtle, generally signaled by the time step becoming smaller and smaller as the computations continue. This reaches a point where progress in the simulation is so slow that the program either terminates because the user-given minimum time step has been reached or the user terminates the program manually (with a control-break on most DOS systems).

The most common cause of an abrupt failure is one or more time-series tables containing an inadequate range of values for all arguments encountered during the computations. Every time-series table referenced in the unsteady-flow model must have a range of arguments that starts at or before the start time (STIME) and ends at or after the end time (ETIME) for the simulation. Remember that values for decimal numbers cannot be specified exactly by the binary representation used in digital computers. Consequently, small round-off and truncation errors can cause the internal value for time used in FEQ and the value stored in the table for time to not agree exactly. To avoid this problem make sure that all time-series tables extend beyond any span of simulation by small amounts. Then small errors will not cause a table overflow or underflow.

4.1 Frozen-time failures

The steady-flow computations in FEQ do not include options for all the possible control structures and conditions that may exist in the unsteady-flow model. Therefore, the steady-flow results will be close to but not exactly the same as the steady-flow results that would be obtained from FEQ if the boundary conditions are held fixed for an extended time period. If the differences are too large, the unsteady-flow computations may fail to converge. The same problem can occur if the convergence tolerance for steady-flow is too large. The primary cause of the problems is large differences in water-surface elevation relative to the depth.If they are too large, the rapid transients generated in the hydraulic system can produce unrealistic flow patterns that will stop the computational process.

Frozen-time computations are unsteady-flow computations done with the time held fixed (frozen) so that the boundary conditions do not change. If all water-surface elevations and flows used and computed in the backwater computations were exactly consistent with the unsteady-flow computations then the frozen-time computations would be steady-flow but they would be computed as unsteady-flow. However, in all but simple models there can be small differences in the estimated initial conditions. These small differences will cause unsteady flow during the frozen time computations. The model is simulating small dam breaks at each water-surface discontinuity. During these computations these transients should become much smaller so that the computations after time is unfrozen will begin at a defined initial condition and not one that depends on the existence of discrepancies in the initial conditions. The example output file ends with the following message:

TIME STEP TOO SMALL. RUN TERMINATED. 

Thus if the computations fail during frozen-time simulation, first look for a large difference in water-surface elevation relative to the depth. Next, check the boundary conditions that involve flow and time. If there is a large difference in flow in the branch and the flow imposed by the boundary condition, especially if the flow at the boundary is larger than the flow in the branch, the same computational problems can occur. If the boundary condition was forced to zero flow or an error in the model specification drains water from the branch, then the problem may be that the branch is becoming dry and the depths are too shallow. The problem can often be solved by gradually increasing the amount of discharge estimated for the steady-flow initial conditions in the branches near the error.

4.2 Zero pivots

Another cause of abrupt failure is much more subtle and fortunately infrequent. Extensive checks are done in FEQ on the Network Matrix specification to make sure the various pieces of the network are connected properly. The messages issued are usually quite specific about the problem so that it can be easily corrected. However, on rare occasions, a problem arises that can be difficult to trace. This results in the "zero pivot" error message. In the solution process, FEQ must solve a linear system of equations while solving the non-linear system of equations. In order to preserve the structure of the linear system, FEQ does not use any row or column interchanges. Such interchanges, designed to reduce round-off and truncation error, would radically alter the structure of the non-zero coefficients in the linear system and result in a large increase in the computations needed. In the process of linear solution, FEQ must divide by the current diagonal element in the linear system, called the pivot element. If this number is zero, the linear solution cannot continue. The linear system is still valid and has a solution but it cannot be found without interchanging rows in the linear system of equations.

It is not possible to guarantee that the linear system will be such that it can be solved without row interchanges. The value of the current diagonal element changes as the computations progress so that there is no simple way of avoiding the problem. FEQ expends considerable computational time to create the linear system of equations so that a zero pivot is avoided most of the time. However, some guidelines are presented below on what to do if a zero pivot is encountered.

The error message gives a row and column number in the coefficient matrix of the linear system. This is sometimes helpful in finding the part of the model that causes the problem. As the structure of the linear system is developed, FEQ creates a table that gives the row numbers that each entry in the Network Matrix creates. The output from addpnd2a.o is shown in Figure 1. The input order is the order of input for each Network Matrix item, and is given as the input items are processed, as shown in line 399. The equation number given for each Network Matrix Input entry in lines 428 and following is the LAST equation number created by the entry. Line 428 shows a boundary condition that creates only one equation, the first equation in the system. Branch 10 creates equation numbers 2 through 45. Equation 45 is the last equation for branch 10 and is generated by the last computational element, that is, the most downstream element on branch 10. Branches create two or more equations; dummy branches and level-pool reservoirs create two equations; and all other options in FEQ create only one equation.

Line 468 shows a negative branch number. However, no negative branch numbers are allowed in the input. This is the means used by FEQ to signal that the equations for branch 30 are written in reverse order, from downstream to upstream, instead of the normal order, from upstream to downstream. Therefore, the last equation for branch 30 (equation 235) was generated by the first element, that is, the most upstream element on branch 30.

FIGURE 1. Structure of the Network Matrix

   Line           Line contents
 Number
 ------ ----------------------------------------------------------------
    399     Code=   5 Input order#=   40
    400     Type=   6 Ups head node= D30 Dwns head node=  U3 Q node=  U3
    401     Tab# for pos flow=  506 Tab# for neg flow=  507
    402     Tab# for adjustment factor=    0 Base elev.=    677.90
    403     Tab# for variable base elevation=    0
    ...
    412  
    413   BNODE=
    414  
    415 0Number of elements in matrix.
    416  BDY NODE  LENGTH
    417       U10    1192
    418        D2    1827
    419       U30    1263
    420       F50    1306
    421 0BOUNDARY NODE USED= U10
    422  
    423  
    424  Order of equations in Network Matrix.
    425  
    426  Equation Input Summary of Input
    427    Number Order
    428         1     1  Code=  6 U10
    429        45     2  Code=  1 Bra=  10
    430        46     4  Code= 13  D10  U11
    431        47     3  Code=  2  D10  U11   F9
    432        49     7  Code=  1 Bra=  11
    433        50     9  Code= 13  D11   U1
      ...
    441        81    11  Code= 15  F11  F12    1
    442        82    14  Code=  3  F49  F12
    443        83    12  Code=  2  F49  F10  F12
      ...
    457       190    36  Code=  3   D3   F3
    458       191    29  Code=  5 Unn=  F2 Dnn=  U4 Qnn=  F2
    459       193    32  Code= 15   F1   F2   -1
    460       194    35  Code=  3   D3   F1

    461       195    28  Code=  5 Unn=  F6 Dnn=  U4 Qnn=  F6
    462       197    31  Code= 15   F5   F6   -1
    463       198    34  Code=  3   D3   F5
    464       199    37  Code=  2   D3   F3   F1   F5
    465       217    38  Code=  1 Bra=  -3
    466       218    40  Code=  5 Unn= D30 Dnn=  U3 Qnn=  U3
    467       219    39  Code=  2  D30   U3
    468       235    41  Code=  1 Bra= -30
    469       236    42  Code=  6 U30

This output table is useful in locating the zero pivot in the matrix. The problem occurred above the row given in the message, but how far above this row is unknown. The following actions should eliminate the error:

1. Rearrange the specification of connections so that the same hydraulic system is represented but the connections are given in a different order. For example, if two or more Code 3 values are at one junction, select a different node as the common elevation node.

2. Explicitly give a beginning node value, BNODE. If BNODE is left blank, FEQ will try all valid beginning nodes and will then select the node that gives the smallest element count for the matrix.

3. Sometimes a zero pivot results from a problem in the specification of the network that FEQ does not detect. Inexperienced users often inadvertently create zero-pivot problems, which suggests that there is a distinct pattern of development which minimizes such errors. The pattern that reduces zero pivots, however, has not been rigorously defined. Here are some general rules that may help reduce zero-pivot errors:

3.1 Never put dummy branches in series. It may work but is poor practice, because the null branches that have been replaced by dummy branches would not work if placed in series. If needed, a very short (about 1 foot) branch could be placed between dummy branches. The storage and friction introduced by the real branch are negligible.

3.2 If a tributary stream enters immediately above or immediately below a culvert or a bridge or other structure, insert a short (up to ten feet long), real branch so that the tributary enters in its own junction distinct from the culvert or bridge.

3.3 If two culverts or bridges or similar structures are right at a junction, insert the needed number of short real branches so that the culvert or bridge would have a distinct junction.

4. For level-pool reservoirs represented by type 3 capacity tables, if the first elevation input is the initial condition of the reservoir, and the volume is zero, then the area should have a non-zero value.

5. Many times, the addition of an interpolated cross section in the branches upstream and downstream of the branch identified in the zero-pivot error message row number will solve the problem.

6. There is often more than one way to represent the same hydraulic system using FEQ. If the zero-pivot error message has a row number that is at or close to a component of the system that can be represented in a different way, then try doing so if the previous actions prove unsuccessful.

4.3 Table overflows

A potential problem with one- or two-dimensional function tables is insufficient extent of the tabulated function. Check for warning messages about the maximum level being exceeded in one or more tables. When a level in a function table is not high enough, a warning message is issued, and the maximum value in the table is used to continue the computations. The maximum value in the table might have been exceeded during the iterations leading to a converged result that was below the top of the table. If this is the case, there will probably not be any major computational problem. However, if the flows are much larger than expected or if some other error has distorted the model, then there may be hundreds or thousands of warnings about overtopped function tables. In this case, the run may fail as the time step declines. Clearly, the reasons for the overtopped function tables must be determined. If the run is completed, FEQ will output a list of the one- and two-dimensional tables that were exceeded at least once by converged results. This list can be helpful in locating the tables that may need to be computed to higher levels.

Table below range

The cause of this error is discussed in section 3.4. The error can usually be avoided by extending the capacity table for the level-pool reservoirs to an elevation of 0.0 or invoking the elevation-to-depth argument translation option.