dwflow.c
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | //----------------------------------------------------------------------------- | ||
| 2 | // dwflow.c | ||
| 3 | // | ||
| 4 | // Project: EPA SWMM5 | ||
| 5 | // Version: 5.2 | ||
| 6 | // Date: 06/12/23 (Build 5.2.4) | ||
| 7 | // Author: L. Rossman | ||
| 8 | // M. Tryby (EPA) | ||
| 9 | // R. Dickinson (CDM) | ||
| 10 | // | ||
| 11 | // Solves the momentum equation for flow in a conduit under dynamic wave | ||
| 12 | // flow routing. | ||
| 13 | // | ||
| 14 | // Update History | ||
| 15 | // ============== | ||
| 16 | // Build 5.1.008: | ||
| 17 | // - Bug in finding if conduit was upstrm/dnstrm full was fixed. | ||
| 18 | // Build 5.1.012: | ||
| 19 | // - Modified uniform loss rate term of conduit momentum equation. | ||
| 20 | // Build 5.1.013: | ||
| 21 | // - Preissmann slot surcharge option implemented. | ||
| 22 | // - Changed sign of uniform loss rate term (dq6) in flow updating equation. | ||
| 23 | // Build 5.1.014: | ||
| 24 | // - Conduit evap. and seepage loss initialized to 0 in dwflow_findConduitFlow. | ||
| 25 | // - Most current flow (qLast) used instead of previous time period flow | ||
| 26 | // (qOld) in call to link_getLossRate. | ||
| 27 | // Build 5.2.1: | ||
| 28 | // - Implements the new option to skip checking for normal flow limitations. | ||
| 29 | // Build 5.2.4: | ||
| 30 | // - Arguments to function link_getLossRate changed. | ||
| 31 | //----------------------------------------------------------------------------- | ||
| 32 | #define _CRT_SECURE_NO_DEPRECATE | ||
| 33 | |||
| 34 | #include <math.h> | ||
| 35 | #include "headers.h" | ||
| 36 | |||
| 37 | static const double MAXVELOCITY = 50.; // max. allowable velocity (ft/sec) | ||
| 38 | |||
| 39 | static int getFlowClass(int link, double q, double h1, double h2, | ||
| 40 | double y1, double y2, double* criticalDepth, double* normalDepth, | ||
| 41 | double* fasnh); | ||
| 42 | static void findSurfArea(int link, double q, double length, double* h1, | ||
| 43 | double* h2, double* y1, double* y2); | ||
| 44 | static double findLocalLosses(int link, double a1, double a2, double aMid, | ||
| 45 | double q); | ||
| 46 | |||
| 47 | static double getWidth(TXsect* xsect, double y); | ||
| 48 | static double getSlotWidth(TXsect* xsect, double y); | ||
| 49 | static double getArea(TXsect* xsect, double y, double wSlot); | ||
| 50 | static double getHydRad(TXsect* xsect, double y); | ||
| 51 | |||
| 52 | static double checkNormalFlow(int j, double q, double y1, double y2, | ||
| 53 | double a1, double r1); | ||
| 54 | |||
| 55 | //============================================================================= | ||
| 56 | |||
| 57 | 50609586 | void dwflow_findConduitFlow(int j, int steps, double omega, double dt) | |
| 58 | // | ||
| 59 | // Input: j = link index | ||
| 60 | // steps = number of iteration steps taken | ||
| 61 | // omega = under-relaxation parameter | ||
| 62 | // dt = time step (sec) | ||
| 63 | // Output: returns new flow value (cfs) | ||
| 64 | // Purpose: updates flow in conduit link by solving finite difference | ||
| 65 | // form of continuity and momentum equations. | ||
| 66 | // | ||
| 67 | { | ||
| 68 | int k; // index of conduit | ||
| 69 | int n1, n2; // indexes of end nodes | ||
| 70 | double z1, z2; // upstream/downstream invert elev. (ft) | ||
| 71 | double h1, h2; // upstream/dounstream flow heads (ft) | ||
| 72 | double y1, y2; // upstream/downstream flow depths (ft) | ||
| 73 | double a1, a2; // upstream/downstream flow areas (ft2) | ||
| 74 | double r1; // upstream hyd. radius (ft) | ||
| 75 | double yMid, rMid, aMid; // mid-stream or avg. values of y, r, & a | ||
| 76 | double aWtd, rWtd; // upstream weighted area & hyd. radius | ||
| 77 | double qLast; // flow from previous iteration (cfs) | ||
| 78 | double qOld; // flow from previous time step (cfs) | ||
| 79 | double aOld; // area from previous time step (ft2) | ||
| 80 | double v; // velocity (ft/sec) | ||
| 81 | double rho; // upstream weighting factor | ||
| 82 | double sigma; // inertial damping factor | ||
| 83 | double length; // effective conduit length (ft) | ||
| 84 | double wSlot; // Preissmann slot width (ft) | ||
| 85 | double dq1, dq2, dq3, dq4, dq5, // terms in momentum eqn. | ||
| 86 | dq6; // term for evap and infil losses | ||
| 87 | double denom; // denominator of flow update formula | ||
| 88 | double q; // new flow value (cfs) | ||
| 89 | double barrels; // number of barrels in conduit | ||
| 90 | 50609586 | TXsect* xsect = &Link[j].xsect; // ptr. to conduit's cross section data | |
| 91 | 50609586 | char isFull = FALSE; // TRUE if conduit flowing full | |
| 92 | 50609586 | char isClosed = FALSE; // TRUE if conduit closed | |
| 93 | |||
| 94 | |||
| 95 | |||
| 96 | // --- adjust isClosed status by any control action | ||
| 97 |
2/2✓ Branch 0 taken 24574 times.
✓ Branch 1 taken 50585012 times.
|
50609586 | if ( Link[j].setting == 0 ) isClosed = TRUE; |
| 98 | |||
| 99 | // --- get flow from last time step & previous iteration | ||
| 100 | 50609586 | k = Link[j].subIndex; | |
| 101 | 50609586 | barrels = Conduit[k].barrels; | |
| 102 | 50609586 | qOld = Link[j].oldFlow / barrels; | |
| 103 | 50609586 | qLast = Conduit[k].q1; | |
| 104 | 50609586 | Conduit[k].evapLossRate = 0.0; | |
| 105 | 50609586 | Conduit[k].seepLossRate = 0.0; | |
| 106 | |||
| 107 | // --- get most current heads at upstream and downstream ends of conduit | ||
| 108 | 50609586 | n1 = Link[j].node1; | |
| 109 | 50609586 | n2 = Link[j].node2; | |
| 110 | 50609586 | z1 = Node[n1].invertElev + Link[j].offset1; | |
| 111 | 50609586 | z2 = Node[n2].invertElev + Link[j].offset2; | |
| 112 | 50609586 | h1 = Node[n1].newDepth + Node[n1].invertElev; | |
| 113 | 50609586 | h2 = Node[n2].newDepth + Node[n2].invertElev; | |
| 114 |
2/2✓ Branch 0 taken 50380795 times.
✓ Branch 1 taken 228791 times.
|
50609586 | h1 = MAX(h1, z1); |
| 115 |
2/2✓ Branch 0 taken 49642828 times.
✓ Branch 1 taken 966758 times.
|
50609586 | h2 = MAX(h2, z2); |
| 116 | |||
| 117 | // --- get unadjusted upstream and downstream flow depths in conduit | ||
| 118 | // (flow depth = head in conduit - elev. of conduit invert) | ||
| 119 | 50609586 | y1 = h1 - z1; | |
| 120 | 50609586 | y2 = h2 - z2; | |
| 121 |
2/2✓ Branch 0 taken 49181042 times.
✓ Branch 1 taken 1428544 times.
|
50609586 | y1 = MAX(y1, FUDGE); |
| 122 |
2/2✓ Branch 0 taken 48462287 times.
✓ Branch 1 taken 2147299 times.
|
50609586 | y2 = MAX(y2, FUDGE); |
| 123 | |||
| 124 | // --- flow depths can't exceed full depth of conduit if slot not used | ||
| 125 |
1/2✓ Branch 0 taken 50609586 times.
✗ Branch 1 not taken.
|
50609586 | if ( SurchargeMethod != SLOT ) |
| 126 | { | ||
| 127 |
2/2✓ Branch 0 taken 49519735 times.
✓ Branch 1 taken 1089851 times.
|
50609586 | y1 = MIN(y1, xsect->yFull); |
| 128 |
2/2✓ Branch 0 taken 47851831 times.
✓ Branch 1 taken 2757755 times.
|
50609586 | y2 = MIN(y2, xsect->yFull); |
| 129 | } | ||
| 130 | |||
| 131 | // -- get area from solution at previous time step | ||
| 132 | 50609586 | aOld = Conduit[k].a2; | |
| 133 |
2/2✓ Branch 0 taken 49097627 times.
✓ Branch 1 taken 1511959 times.
|
50609586 | aOld = MAX(aOld, FUDGE); |
| 134 | |||
| 135 | // --- use Courant-modified length instead of conduit's actual length | ||
| 136 | 50609586 | length = Conduit[k].modLength; | |
| 137 | |||
| 138 | // --- find surface area contributions to upstream and downstream nodes | ||
| 139 | // based on previous iteration's flow estimate | ||
| 140 | 50609586 | findSurfArea(j, qLast, length, &h1, &h2, &y1, &y2); | |
| 141 | |||
| 142 | // --- compute area at each end of conduit & hyd. radius at upstream end | ||
| 143 | 50609586 | wSlot = getSlotWidth(xsect, y1); | |
| 144 | 50609586 | a1 = getArea(xsect, y1, wSlot); | |
| 145 | 50609586 | r1 = getHydRad(xsect, y1); | |
| 146 | 50609586 | wSlot = getSlotWidth(xsect, y2); | |
| 147 | 50609586 | a2 = getArea(xsect, y2, wSlot); | |
| 148 | |||
| 149 | // --- compute area & hyd. radius at midpoint | ||
| 150 | 50609586 | yMid = 0.5 * (y1 + y2); | |
| 151 | 50609586 | wSlot = getSlotWidth(xsect, yMid); | |
| 152 | 50609586 | aMid = getArea(xsect, yMid, wSlot); | |
| 153 | 50609586 | rMid = getHydRad(xsect, yMid); | |
| 154 | |||
| 155 | // --- alternate approach not currently used, but might produce better | ||
| 156 | // Bernoulli energy balance for steady flows | ||
| 157 | //aMid = (a1+a2)/2.0; | ||
| 158 | //rMid = (r1+getHydRad(xsect,y2))/2.0; | ||
| 159 | |||
| 160 | // --- check if conduit is flowing full | ||
| 161 |
2/2✓ Branch 0 taken 1119494 times.
✓ Branch 1 taken 49490092 times.
|
50609586 | if ( y1 >= xsect->yFull && |
| 162 |
2/2✓ Branch 0 taken 858626 times.
✓ Branch 1 taken 260868 times.
|
1119494 | y2 >= xsect->yFull) isFull = TRUE; |
| 163 | |||
| 164 | // --- set new flow to zero if conduit is dry or if flap gate is closed | ||
| 165 |
2/2✓ Branch 0 taken 49689375 times.
✓ Branch 1 taken 920211 times.
|
50609586 | if ( Link[j].flowClass == DRY || |
| 166 |
2/2✓ Branch 0 taken 49199017 times.
✓ Branch 1 taken 490358 times.
|
49689375 | Link[j].flowClass == UP_DRY || |
| 167 |
3/4✓ Branch 0 taken 49199017 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 49174453 times.
✓ Branch 3 taken 24564 times.
|
49199017 | Link[j].flowClass == DN_DRY || |
| 168 |
2/2✓ Branch 0 taken 923342 times.
✓ Branch 1 taken 48251111 times.
|
49174453 | isClosed || |
| 169 | aMid <= FUDGE ) | ||
| 170 | { | ||
| 171 | 2358475 | Conduit[k].a1 = 0.5 * (a1 + a2); | |
| 172 | 2358475 | Conduit[k].q1 = 0.0;; | |
| 173 | 2358475 | Conduit[k].q2 = 0.0; | |
| 174 | 2358475 | Link[j].dqdh = GRAVITY * dt * aMid / length * barrels; | |
| 175 | 2358475 | Link[j].froude = 0.0; | |
| 176 |
1/2✓ Branch 0 taken 2358475 times.
✗ Branch 1 not taken.
|
2358475 | Link[j].newDepth = MIN(yMid, Link[j].xsect.yFull); |
| 177 | 2358475 | Link[j].newVolume = Conduit[k].a1 * link_getLength(j) * barrels; | |
| 178 | 2358475 | Link[j].newFlow = 0.0; | |
| 179 | 2358475 | return; | |
| 180 | } | ||
| 181 | |||
| 182 | // --- compute velocity from last flow estimate | ||
| 183 | 48251111 | v = qLast / aMid; | |
| 184 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 48251111 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
48251111 | if ( fabs(v) > MAXVELOCITY ) v = MAXVELOCITY * SGN(qLast); |
| 185 | |||
| 186 | // --- compute Froude No. | ||
| 187 | 48251111 | Link[j].froude = link_getFroude(j, v, yMid); | |
| 188 |
2/2✓ Branch 0 taken 47103639 times.
✓ Branch 1 taken 1147472 times.
|
48251111 | if ( Link[j].flowClass == SUBCRITICAL && |
| 189 |
2/2✓ Branch 0 taken 11350865 times.
✓ Branch 1 taken 35752774 times.
|
47103639 | Link[j].froude > 1.0 ) Link[j].flowClass = SUPCRITICAL; |
| 190 | |||
| 191 | // --- find inertial damping factor (sigma) | ||
| 192 |
2/2✓ Branch 0 taken 17006584 times.
✓ Branch 1 taken 31244527 times.
|
48251111 | if ( Link[j].froude <= 0.5 ) sigma = 1.0; |
| 193 |
2/2✓ Branch 0 taken 12097612 times.
✓ Branch 1 taken 19146915 times.
|
31244527 | else if ( Link[j].froude >= 1.0 ) sigma = 0.0; |
| 194 | 19146915 | else sigma = 2.0 * (1.0 - Link[j].froude); | |
| 195 | |||
| 196 | // --- get upstream-weighted area & hyd. radius based on damping factor | ||
| 197 | // (modified version of R. Dickinson's slope weighting) | ||
| 198 | 48251111 | rho = 1.0; | |
| 199 |
6/6✓ Branch 0 taken 47392485 times.
✓ Branch 1 taken 858626 times.
✓ Branch 2 taken 45030591 times.
✓ Branch 3 taken 2361894 times.
✓ Branch 4 taken 44494289 times.
✓ Branch 5 taken 536302 times.
|
48251111 | if ( !isFull && qLast > 0.0 && h1 >= h2 ) rho = sigma; |
| 200 | 48251111 | aWtd = a1 + (aMid - a1) * rho; | |
| 201 | 48251111 | rWtd = r1 + (rMid - r1) * rho; | |
| 202 | |||
| 203 | // --- determine how much inertial damping to apply | ||
| 204 |
2/2✓ Branch 0 taken 4821254 times.
✓ Branch 1 taken 43429857 times.
|
48251111 | if ( InertDamping == NO_DAMPING ) sigma = 1.0; |
| 205 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 43429857 times.
|
43429857 | else if ( InertDamping == FULL_DAMPING ) sigma = 0.0; |
| 206 | |||
| 207 | // --- use full inertial damping if closed conduit is surcharged | ||
| 208 |
4/4✓ Branch 0 taken 858626 times.
✓ Branch 1 taken 47392485 times.
✓ Branch 3 taken 737415 times.
✓ Branch 4 taken 121211 times.
|
48251111 | if ( isFull && !xsect_isOpen(xsect->type) ) sigma = 0.0; |
| 209 | |||
| 210 | // --- compute terms of momentum eqn.: | ||
| 211 | // --- 1. friction slope term | ||
| 212 |
4/4✓ Branch 0 taken 28557 times.
✓ Branch 1 taken 48222554 times.
✓ Branch 2 taken 23626 times.
✓ Branch 3 taken 4931 times.
|
48251111 | if ( xsect->type == FORCE_MAIN && isFull ) |
| 213 | 23626 | dq1 = dt * forcemain_getFricSlope(j, fabs(v), rMid); | |
| 214 | 48227485 | else dq1 = dt * Conduit[k].roughFactor / pow(rWtd, 1.33333) * fabs(v); | |
| 215 | |||
| 216 | // --- 2. energy slope term | ||
| 217 | 48251111 | dq2 = dt * GRAVITY * aWtd * (h2 - h1) / length; | |
| 218 | |||
| 219 | // --- 3 & 4. inertial terms | ||
| 220 | 48251111 | dq3 = 0.0; | |
| 221 | 48251111 | dq4 = 0.0; | |
| 222 |
2/2✓ Branch 0 taken 36739187 times.
✓ Branch 1 taken 11511924 times.
|
48251111 | if ( sigma > 0.0 ) |
| 223 | { | ||
| 224 | 36739187 | dq3 = 2.0 * v * (aMid - aOld) * sigma; | |
| 225 | 36739187 | dq4 = dt * v * v * (a2 - a1) / length * sigma; | |
| 226 | } | ||
| 227 | |||
| 228 | // --- 5. local losses term | ||
| 229 | 48251111 | dq5 = 0.0; | |
| 230 |
2/2✓ Branch 0 taken 160998 times.
✓ Branch 1 taken 48090113 times.
|
48251111 | if ( Conduit[k].hasLosses ) |
| 231 | { | ||
| 232 | 160998 | dq5 = findLocalLosses(j, a1, a2, aMid, qLast) / 2.0 / length * dt; | |
| 233 | } | ||
| 234 | |||
| 235 | // --- 6. term for evap and seepage losses per unit length | ||
| 236 | 48251111 | dq6 = link_getLossRate(j, DW, qLast, dt) * 2.5 * dt * v / link_getLength(j); | |
| 237 | |||
| 238 | // --- combine terms to find new conduit flow | ||
| 239 | 48251111 | denom = 1.0 + dq1 + dq5; | |
| 240 | 48251111 | q = (qOld - dq2 + dq3 + dq4 + dq6) / denom; | |
| 241 | |||
| 242 | // --- compute derivative of flow w.r.t. head | ||
| 243 | 48251111 | Link[j].dqdh = 1.0 / denom * GRAVITY * dt * aWtd / length * barrels; | |
| 244 | |||
| 245 | // --- check if any flow limitation applies | ||
| 246 | 48251111 | Link[j].inletControl = FALSE; | |
| 247 | 48251111 | Link[j].normalFlow = FALSE; | |
| 248 |
2/2✓ Branch 0 taken 45865268 times.
✓ Branch 1 taken 2385843 times.
|
48251111 | if ( q > 0.0 ) |
| 249 | { | ||
| 250 | // --- check for inlet controlled culvert flow | ||
| 251 |
4/4✓ Branch 0 taken 35999 times.
✓ Branch 1 taken 45829269 times.
✓ Branch 2 taken 34053 times.
✓ Branch 3 taken 1946 times.
|
45865268 | if ( xsect->culvertCode > 0 && !isFull ) |
| 252 | 34053 | q = culvert_getInflow(j, q, h1); | |
| 253 | |||
| 254 | // --- check for normal flow limitation based on surface slope & Fr | ||
| 255 |
3/4✓ Branch 0 taken 45831215 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 44804860 times.
✓ Branch 3 taken 1026355 times.
|
45831215 | else if (NormalFlowLtd != NEITHER && y1 < Link[j].xsect.yFull && |
| 256 |
2/2✓ Branch 0 taken 12276272 times.
✓ Branch 1 taken 32528588 times.
|
44804860 | ( Link[j].flowClass == SUBCRITICAL || |
| 257 |
2/2✓ Branch 0 taken 11239859 times.
✓ Branch 1 taken 1036413 times.
|
12276272 | Link[j].flowClass == SUPCRITICAL )) |
| 258 | 43768447 | q = checkNormalFlow(j, q, y1, y2, a1, r1); | |
| 259 | } | ||
| 260 | |||
| 261 | // --- apply under-relaxation weighting between new & old flows; | ||
| 262 | // --- do not allow change in flow direction without first being zero | ||
| 263 |
2/2✓ Branch 0 taken 24628898 times.
✓ Branch 1 taken 23622213 times.
|
48251111 | if ( steps > 0 ) |
| 264 | { | ||
| 265 | 24628898 | q = (1.0 - omega) * qLast + omega * q; | |
| 266 |
4/4✓ Branch 0 taken 2671 times.
✓ Branch 1 taken 24626227 times.
✓ Branch 2 taken 870 times.
✓ Branch 3 taken 1801 times.
|
24628898 | if ( q * qLast < 0.0 ) q = 0.001 * SGN(q); |
| 267 | } | ||
| 268 | |||
| 269 | // --- check if user-supplied flow limit applies | ||
| 270 |
2/2✓ Branch 0 taken 86000 times.
✓ Branch 1 taken 48165111 times.
|
48251111 | if ( Link[j].qLimit > 0.0 ) |
| 271 | { | ||
| 272 |
3/4✓ Branch 0 taken 41237 times.
✓ Branch 1 taken 44763 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 41237 times.
|
86000 | if ( fabs(q) > Link[j].qLimit ) q = SGN(q) * Link[j].qLimit; |
| 273 | } | ||
| 274 | |||
| 275 | // --- check for reverse flow with closed flap gate | ||
| 276 |
2/2✓ Branch 1 taken 31056 times.
✓ Branch 2 taken 48220055 times.
|
48251111 | if ( link_setFlapGate(j, n1, n2, q) ) q = 0.0; |
| 277 | |||
| 278 | // --- do not allow flow out of a dry node | ||
| 279 | // (as suggested by R. Dickinson) | ||
| 280 |
3/4✓ Branch 0 taken 45169757 times.
✓ Branch 1 taken 3081354 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 45169757 times.
|
48251111 | if( q > FUDGE && Node[n1].newDepth <= FUDGE ) q = FUDGE; |
| 281 |
4/4✓ Branch 0 taken 2337241 times.
✓ Branch 1 taken 45913870 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 2337236 times.
|
48251111 | if( q < -FUDGE && Node[n2].newDepth <= FUDGE ) q = -FUDGE; |
| 282 | |||
| 283 | // --- save new values of area, flow, depth, & volume | ||
| 284 | 48251111 | Conduit[k].a1 = aMid; | |
| 285 | 48251111 | Conduit[k].q1 = q; | |
| 286 | 48251111 | Conduit[k].q2 = q; | |
| 287 |
1/2✓ Branch 0 taken 48251111 times.
✗ Branch 1 not taken.
|
48251111 | Link[j].newDepth = MIN(yMid, xsect->yFull); |
| 288 | 48251111 | aMid = (a1 + a2) / 2.0; | |
| 289 | // aMid = MIN(aMid, xsect->aFull); //Slot can have aMid > aFull | ||
| 290 | 48251111 | Conduit[k].fullState = link_getFullState(a1, a2, xsect->aFull); | |
| 291 | 48251111 | Link[j].newVolume = aMid * link_getLength(j) * barrels; | |
| 292 | 48251111 | Link[j].newFlow = q * barrels; | |
| 293 | } | ||
| 294 | |||
| 295 | //============================================================================= | ||
| 296 | |||
| 297 | 49750960 | int getFlowClass(int j, double q, double h1, double h2, double y1, double y2, | |
| 298 | double *yC, double *yN, double* fasnh) | ||
| 299 | // | ||
| 300 | // Input: j = conduit link index | ||
| 301 | // q = current conduit flow (cfs) | ||
| 302 | // h1 = head at upstream end of conduit (ft) | ||
| 303 | // h2 = head at downstream end of conduit (ft) | ||
| 304 | // y1 = upstream flow depth in conduit (ft) | ||
| 305 | // y2 = downstream flow depth in conduit (ft) | ||
| 306 | // yC = critical flow depth (ft) | ||
| 307 | // yN = normal flow depth (ft) | ||
| 308 | // fasnh = fraction between norm. & crit. depth | ||
| 309 | // Output: returns flow classification code | ||
| 310 | // Purpose: determines flow class for a conduit based on depths at each end. | ||
| 311 | // | ||
| 312 | { | ||
| 313 | int n1, n2; // indexes of upstrm/downstrm nodes | ||
| 314 | int flowClass; // flow classification code | ||
| 315 | double ycMin, ycMax; // min/max critical depths (ft) | ||
| 316 | double z1, z2; // offsets of conduit inverts (ft) | ||
| 317 | |||
| 318 | // --- get upstream & downstream node indexes | ||
| 319 | 49750960 | n1 = Link[j].node1; | |
| 320 | 49750960 | n2 = Link[j].node2; | |
| 321 | |||
| 322 | // --- get upstream & downstream conduit invert offsets | ||
| 323 | 49750960 | z1 = Link[j].offset1; | |
| 324 | 49750960 | z2 = Link[j].offset2; | |
| 325 | |||
| 326 | // --- base offset of an outfall conduit on outfall's depth | ||
| 327 |
3/4✓ Branch 0 taken 1547 times.
✓ Branch 1 taken 49749413 times.
✓ Branch 2 taken 1547 times.
✗ Branch 3 not taken.
|
49750960 | if ( Node[n1].type == OUTFALL ) z1 = MAX(0.0, (z1 - Node[n1].newDepth)); |
| 328 |
3/4✓ Branch 0 taken 1917017 times.
✓ Branch 1 taken 47833943 times.
✓ Branch 2 taken 1917017 times.
✗ Branch 3 not taken.
|
49750960 | if ( Node[n2].type == OUTFALL ) z2 = MAX(0.0, (z2 - Node[n2].newDepth)); |
| 329 | |||
| 330 | // --- default class is SUBCRITICAL | ||
| 331 | 49750960 | flowClass = SUBCRITICAL; | |
| 332 | 49750960 | *fasnh = 1.0; | |
| 333 | |||
| 334 | // --- case where both ends of conduit are wet | ||
| 335 |
4/4✓ Branch 0 taken 48322416 times.
✓ Branch 1 taken 1428544 times.
✓ Branch 2 taken 47095328 times.
✓ Branch 3 taken 1227088 times.
|
49750960 | if ( y1 > FUDGE && y2 > FUDGE ) |
| 336 | { | ||
| 337 |
2/2✓ Branch 0 taken 2283819 times.
✓ Branch 1 taken 44811509 times.
|
47095328 | if ( q < 0.0 ) |
| 338 | { | ||
| 339 | // --- upstream end at critical depth if flow depth is | ||
| 340 | // below conduit's critical depth and an upstream | ||
| 341 | // conduit offset exists | ||
| 342 |
2/2✓ Branch 0 taken 33729 times.
✓ Branch 1 taken 2250090 times.
|
2283819 | if ( z1 > 0.0 ) |
| 343 | { | ||
| 344 | 33729 | *yN = link_getYnorm(j, fabs(q)); | |
| 345 | 33729 | *yC = link_getYcrit(j, fabs(q)); | |
| 346 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 33729 times.
|
33729 | ycMin = MIN(*yN, *yC); |
| 347 |
2/2✓ Branch 0 taken 30122 times.
✓ Branch 1 taken 3607 times.
|
33729 | if ( y1 < ycMin ) flowClass = UP_CRITICAL; |
| 348 | } | ||
| 349 | } | ||
| 350 | |||
| 351 | // --- case of normal direction flow | ||
| 352 | else | ||
| 353 | { | ||
| 354 | // --- downstream end at smaller of critical and normal depth | ||
| 355 | // if downstream flow depth below this and a downstream | ||
| 356 | // conduit offset exists | ||
| 357 |
2/2✓ Branch 0 taken 948782 times.
✓ Branch 1 taken 43862727 times.
|
44811509 | if ( z2 > 0.0 ) |
| 358 | { | ||
| 359 | 948782 | *yN = link_getYnorm(j, fabs(q)); | |
| 360 | 948782 | *yC = link_getYcrit(j, fabs(q)); | |
| 361 |
2/2✓ Branch 0 taken 310297 times.
✓ Branch 1 taken 638485 times.
|
948782 | ycMin = MIN(*yN, *yC); |
| 362 |
2/2✓ Branch 0 taken 693950 times.
✓ Branch 1 taken 254832 times.
|
948782 | ycMax = MAX(*yN, *yC); |
| 363 |
2/2✓ Branch 0 taken 373139 times.
✓ Branch 1 taken 575643 times.
|
948782 | if ( y2 < ycMin ) flowClass = DN_CRITICAL; |
| 364 |
2/2✓ Branch 0 taken 287004 times.
✓ Branch 1 taken 288639 times.
|
575643 | else if ( y2 < ycMax ) |
| 365 | { | ||
| 366 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 287004 times.
|
287004 | if ( ycMax - ycMin < FUDGE ) *fasnh = 0.0; |
| 367 | 287004 | else *fasnh = (ycMax - y2) / (ycMax - ycMin); | |
| 368 | } | ||
| 369 | } | ||
| 370 | } | ||
| 371 | } | ||
| 372 | |||
| 373 | // --- case where no flow at either end of conduit | ||
| 374 |
4/4✓ Branch 0 taken 1428544 times.
✓ Branch 1 taken 1227088 times.
✓ Branch 2 taken 920211 times.
✓ Branch 3 taken 508333 times.
|
2655632 | else if ( y1 <= FUDGE && y2 <= FUDGE ) flowClass = DRY; |
| 375 | |||
| 376 | // --- case where downstream end of pipe is wet, upstream dry | ||
| 377 |
2/2✓ Branch 0 taken 508333 times.
✓ Branch 1 taken 1227088 times.
|
1735421 | else if ( y2 > FUDGE ) |
| 378 | { | ||
| 379 | // --- flow classification is UP_DRY if downstream head < | ||
| 380 | // invert of upstream end of conduit | ||
| 381 |
2/2✓ Branch 0 taken 490358 times.
✓ Branch 1 taken 17975 times.
|
508333 | if ( h2 < Node[n1].invertElev + Link[j].offset1 ) flowClass = UP_DRY; |
| 382 | |||
| 383 | // --- otherwise, the downstream head will be >= upstream | ||
| 384 | // conduit invert creating a flow reversal and upstream end | ||
| 385 | // should be at critical depth, providing that an upstream | ||
| 386 | // offset exists (otherwise subcritical condition is maintained) | ||
| 387 |
2/2✓ Branch 0 taken 8574 times.
✓ Branch 1 taken 9401 times.
|
17975 | else if ( z1 > 0.0 ) |
| 388 | { | ||
| 389 | 8574 | *yN = link_getYnorm(j, fabs(q)); | |
| 390 | 8574 | *yC = link_getYcrit(j, fabs(q)); | |
| 391 | 8574 | flowClass = UP_CRITICAL; | |
| 392 | } | ||
| 393 | } | ||
| 394 | |||
| 395 | // --- case where upstream end of pipe is wet, downstream dry | ||
| 396 | else | ||
| 397 | { | ||
| 398 | // --- flow classification is DN_DRY if upstream head < | ||
| 399 | // invert of downstream end of conduit | ||
| 400 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1227088 times.
|
1227088 | if ( h1 < Node[n2].invertElev + Link[j].offset2 ) flowClass = DN_DRY; |
| 401 | |||
| 402 | // --- otherwise flow at downstream end should be at critical depth | ||
| 403 | // providing that a downstream offset exists (otherwise | ||
| 404 | // subcritical condition is maintained) | ||
| 405 |
2/2✓ Branch 0 taken 820131 times.
✓ Branch 1 taken 406957 times.
|
1227088 | else if ( z2 > 0.0 ) |
| 406 | { | ||
| 407 | 820131 | *yN = link_getYnorm(j, fabs(q)); | |
| 408 | 820131 | *yC = link_getYcrit(j, fabs(q)); | |
| 409 | 820131 | flowClass = DN_CRITICAL; | |
| 410 | } | ||
| 411 | } | ||
| 412 | 49750960 | return flowClass; | |
| 413 | } | ||
| 414 | |||
| 415 | //============================================================================= | ||
| 416 | |||
| 417 | 50609586 | void findSurfArea(int j, double q, double length, double* h1, double* h2, | |
| 418 | double* y1, double* y2) | ||
| 419 | // | ||
| 420 | // Input: j = conduit link index | ||
| 421 | // q = current conduit flow (cfs) | ||
| 422 | // length = conduit length (ft) | ||
| 423 | // h1 = head at upstream end of conduit (ft) | ||
| 424 | // h2 = head at downstream end of conduit (ft) | ||
| 425 | // y1 = upstream flow depth (ft) | ||
| 426 | // y2 = downstream flow depth (ft) | ||
| 427 | // Output: updated values of h1, h2, y1, & y2; | ||
| 428 | // Purpose: assigns surface area of conduit to its up and downstream nodes. | ||
| 429 | // | ||
| 430 | { | ||
| 431 | int n1, n2; // indexes of upstrm/downstrm nodes | ||
| 432 | double flowDepth1; // flow depth at upstrm end (ft) | ||
| 433 | double flowDepth2; // flow depth at downstrm end (ft) | ||
| 434 | double flowDepthMid; // flow depth at midpt. (ft) | ||
| 435 | double width1; // top width at upstrm end (ft) | ||
| 436 | double width2; // top width at downstrm end (ft) | ||
| 437 | double widthMid; // top width at midpt. (ft) | ||
| 438 | 50609586 | double surfArea1 = 0.0; // surface area at upstream node (ft2) | |
| 439 | 50609586 | double surfArea2 = 0.0; // surface area st downstrm node (ft2) | |
| 440 | double criticalDepth; // critical flow depth (ft) | ||
| 441 | double normalDepth; // normal flow depth (ft) | ||
| 442 | double fullDepth; // full depth (ft) | ||
| 443 | 50609586 | double fasnh = 1.0; // fraction between norm. & crit. depth | |
| 444 | 50609586 | TXsect* xsect = &Link[j].xsect; // pointer to cross-section data | |
| 445 | |||
| 446 | // --- get node indexes & current flow depths | ||
| 447 | 50609586 | n1 = Link[j].node1; | |
| 448 | 50609586 | n2 = Link[j].node2; | |
| 449 | 50609586 | flowDepth1 = *y1; | |
| 450 | 50609586 | flowDepth2 = *y2; | |
| 451 | |||
| 452 | 50609586 | normalDepth = (flowDepth1 + flowDepth2) / 2.0; | |
| 453 | 50609586 | criticalDepth = normalDepth; | |
| 454 | |||
| 455 | // --- find conduit's flow classification | ||
| 456 | 50609586 | fullDepth = xsect->yFull; | |
| 457 |
4/4✓ Branch 0 taken 1119494 times.
✓ Branch 1 taken 49490092 times.
✓ Branch 2 taken 858626 times.
✓ Branch 3 taken 260868 times.
|
50609586 | if (flowDepth1 >= fullDepth && flowDepth2 >= fullDepth) |
| 458 | { | ||
| 459 | 858626 | Link[j].flowClass = SUBCRITICAL; | |
| 460 | } | ||
| 461 | 49750960 | else Link[j].flowClass = getFlowClass(j, q, *h1, *h2, *y1, *y2, | |
| 462 | &criticalDepth, &normalDepth, &fasnh); | ||
| 463 | |||
| 464 | // --- add conduit's surface area to its end nodes depending on flow class | ||
| 465 |
5/7✓ Branch 0 taken 47967051 times.
✓ Branch 1 taken 38696 times.
✓ Branch 2 taken 1193270 times.
✓ Branch 3 taken 490358 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 920211 times.
✗ Branch 6 not taken.
|
50609586 | switch ( Link[j].flowClass ) |
| 466 | { | ||
| 467 | 47967051 | case SUBCRITICAL: | |
| 468 | 47967051 | flowDepthMid = 0.5 * (flowDepth1 + flowDepth2); | |
| 469 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 47967051 times.
|
47967051 | if ( flowDepthMid < FUDGE ) flowDepthMid = FUDGE; |
| 470 | 47967051 | width1 = getWidth(xsect, flowDepth1); | |
| 471 | 47967051 | width2 = getWidth(xsect, flowDepth2); | |
| 472 | 47967051 | widthMid = getWidth(xsect, flowDepthMid); | |
| 473 | 47967051 | surfArea1 = (width1 + widthMid) * length / 4.; | |
| 474 | 47967051 | surfArea2 = (widthMid + width2) * length / 4. * fasnh; | |
| 475 | 47967051 | break; | |
| 476 | |||
| 477 | 38696 | case UP_CRITICAL: | |
| 478 | 38696 | flowDepth1 = criticalDepth; | |
| 479 |
2/2✓ Branch 0 taken 734 times.
✓ Branch 1 taken 37962 times.
|
38696 | if ( normalDepth < criticalDepth ) flowDepth1 = normalDepth; |
| 480 |
2/2✓ Branch 0 taken 38412 times.
✓ Branch 1 taken 284 times.
|
38696 | flowDepth1 = MAX(flowDepth1, FUDGE); |
| 481 | 38696 | *h1 = Node[n1].invertElev + Link[j].offset1 + flowDepth1; | |
| 482 | 38696 | flowDepthMid = 0.5 * (flowDepth1 + flowDepth2); | |
| 483 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 38696 times.
|
38696 | if ( flowDepthMid < FUDGE ) flowDepthMid = FUDGE; |
| 484 | 38696 | width2 = getWidth(xsect, flowDepth2); | |
| 485 | 38696 | widthMid = getWidth(xsect, flowDepthMid); | |
| 486 | 38696 | surfArea2 = (widthMid + width2) * length * 0.5; | |
| 487 | 38696 | break; | |
| 488 | |||
| 489 | 1193270 | case DN_CRITICAL: | |
| 490 | 1193270 | flowDepth2 = criticalDepth; | |
| 491 |
2/2✓ Branch 0 taken 576200 times.
✓ Branch 1 taken 617070 times.
|
1193270 | if ( normalDepth < criticalDepth ) flowDepth2 = normalDepth; |
| 492 |
2/2✓ Branch 0 taken 1097837 times.
✓ Branch 1 taken 95433 times.
|
1193270 | flowDepth2 = MAX(flowDepth2, FUDGE); |
| 493 | 1193270 | *h2 = Node[n2].invertElev + Link[j].offset2 + flowDepth2; | |
| 494 | 1193270 | width1 = getWidth(xsect, flowDepth1); | |
| 495 | 1193270 | flowDepthMid = 0.5 * (flowDepth1 + flowDepth2); | |
| 496 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1193270 times.
|
1193270 | if ( flowDepthMid < FUDGE ) flowDepthMid = FUDGE; |
| 497 | 1193270 | widthMid = getWidth(xsect, flowDepthMid); | |
| 498 | 1193270 | surfArea1 = (width1 + widthMid) * length * 0.5; | |
| 499 | 1193270 | break; | |
| 500 | |||
| 501 | 490358 | case UP_DRY: | |
| 502 | 490358 | flowDepth1 = FUDGE; | |
| 503 | 490358 | flowDepthMid = 0.5 * (flowDepth1 + flowDepth2); | |
| 504 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 490358 times.
|
490358 | if ( flowDepthMid < FUDGE ) flowDepthMid = FUDGE; |
| 505 | 490358 | width1 = getWidth(xsect, flowDepth1); | |
| 506 | 490358 | width2 = getWidth(xsect, flowDepth2); | |
| 507 | 490358 | widthMid = getWidth(xsect, flowDepthMid); | |
| 508 | |||
| 509 | // --- assign avg. surface area of downstream half of conduit | ||
| 510 | // to the downstream node | ||
| 511 | 490358 | surfArea2 = (widthMid + width2) * length / 4.; | |
| 512 | |||
| 513 | // --- if there is no free-fall at upstream end, assign the | ||
| 514 | // upstream node the avg. surface area of the upstream half | ||
| 515 |
2/2✓ Branch 0 taken 472710 times.
✓ Branch 1 taken 17648 times.
|
490358 | if ( Link[j].offset1 <= 0.0 ) |
| 516 | { | ||
| 517 | 472710 | surfArea1 = (width1 + widthMid) * length / 4.; | |
| 518 | } | ||
| 519 | 490358 | break; | |
| 520 | |||
| 521 | ✗ | case DN_DRY: | |
| 522 | ✗ | flowDepth2 = FUDGE; | |
| 523 | ✗ | flowDepthMid = 0.5 * (flowDepth1 + flowDepth2); | |
| 524 | ✗ | if ( flowDepthMid < FUDGE ) flowDepthMid = FUDGE; | |
| 525 | ✗ | width1 = getWidth(xsect, flowDepth1); | |
| 526 | ✗ | width2 = getWidth(xsect, flowDepth2); | |
| 527 | ✗ | widthMid = getWidth(xsect, flowDepthMid); | |
| 528 | |||
| 529 | // --- assign avg. surface area of upstream half of conduit | ||
| 530 | // to the upstream node | ||
| 531 | ✗ | surfArea1 = (widthMid + width1) * length / 4.; | |
| 532 | |||
| 533 | // --- if there is no free-fall at downstream end, assign the | ||
| 534 | // downstream node the avg. surface area of the downstream half | ||
| 535 | ✗ | if ( Link[j].offset2 <= 0.0 ) | |
| 536 | { | ||
| 537 | ✗ | surfArea2 = (width2 + widthMid) * length / 4.; | |
| 538 | } | ||
| 539 | ✗ | break; | |
| 540 | |||
| 541 | 920211 | case DRY: | |
| 542 | 920211 | surfArea1 = FUDGE * length / 2.0; | |
| 543 | 920211 | surfArea2 = surfArea1; | |
| 544 | 920211 | break; | |
| 545 | } | ||
| 546 | 50609586 | Link[j].surfArea1 = surfArea1; | |
| 547 | 50609586 | Link[j].surfArea2 = surfArea2; | |
| 548 | 50609586 | *y1 = flowDepth1; | |
| 549 | 50609586 | *y2 = flowDepth2; | |
| 550 | 50609586 | } | |
| 551 | |||
| 552 | //============================================================================= | ||
| 553 | |||
| 554 | 160998 | double findLocalLosses(int j, double a1, double a2, double aMid, double q) | |
| 555 | // | ||
| 556 | // Input: j = link index | ||
| 557 | // a1 = upstream area (ft2) | ||
| 558 | // a2 = downstream area (ft2) | ||
| 559 | // aMid = midpoint area (ft2) | ||
| 560 | // q = flow rate (cfs) | ||
| 561 | // Output: returns local losses (ft/sec) | ||
| 562 | // Purpose: computes local losses term of momentum equation. | ||
| 563 | // | ||
| 564 | { | ||
| 565 | 160998 | double losses = 0.0; | |
| 566 | 160998 | q = fabs(q); | |
| 567 |
1/2✓ Branch 0 taken 160998 times.
✗ Branch 1 not taken.
|
160998 | if ( a1 > FUDGE ) losses += Link[j].cLossInlet * (q/a1); |
| 568 |
2/2✓ Branch 0 taken 160299 times.
✓ Branch 1 taken 699 times.
|
160998 | if ( a2 > FUDGE ) losses += Link[j].cLossOutlet * (q/a2); |
| 569 |
1/2✓ Branch 0 taken 160998 times.
✗ Branch 1 not taken.
|
160998 | if ( aMid > FUDGE ) losses += Link[j].cLossAvg * (q/aMid); |
| 570 | 160998 | return losses; | |
| 571 | } | ||
| 572 | |||
| 573 | //============================================================================= | ||
| 574 | |||
| 575 | 299664917 | double getSlotWidth(TXsect* xsect, double y) | |
| 576 | { | ||
| 577 | 299664917 | double yNorm = y / xsect->yFull; | |
| 578 | |||
| 579 | // --- return 0.0 if slot surcharge method not used | ||
| 580 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 299664917 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
299664917 | if (SurchargeMethod != SLOT || xsect_isOpen(xsect->type) || |
| 581 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
299664917 | yNorm < CrownCutoff) return 0.0; |
| 582 | |||
| 583 | // --- for depth > 1.78 * pipe depth, slot width = 1% of max. width | ||
| 584 | ✗ | if (yNorm > 1.78) return 0.01 * xsect->wMax; | |
| 585 | |||
| 586 | // --- otherwise use the Sjoberg formula | ||
| 587 | ✗ | return xsect->wMax * 0.5423 * exp(-pow(yNorm, 2.4)); | |
| 588 | } | ||
| 589 | |||
| 590 | //============================================================================= | ||
| 591 | |||
| 592 | 147836159 | double getWidth(TXsect* xsect, double y) | |
| 593 | // | ||
| 594 | // Input: xsect = ptr. to conduit cross section | ||
| 595 | // y = flow depth (ft) | ||
| 596 | // Output: returns top width (ft) | ||
| 597 | // Purpose: computes top width of flow surface in conduit. | ||
| 598 | // | ||
| 599 | { | ||
| 600 | 147836159 | double wSlot = getSlotWidth(xsect, y); | |
| 601 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 147836159 times.
|
147836159 | if (wSlot > 0.0) return wSlot; |
| 602 |
4/4✓ Branch 0 taken 5121705 times.
✓ Branch 1 taken 142714454 times.
✓ Branch 3 taken 4317269 times.
✓ Branch 4 taken 804436 times.
|
147836159 | if (y / xsect->yFull >= CrownCutoff && !xsect_isOpen(xsect->type)) |
| 603 | 4317269 | y = CrownCutoff * xsect->yFull; | |
| 604 | 147836159 | return xsect_getWofY(xsect, y); | |
| 605 | } | ||
| 606 | |||
| 607 | //============================================================================= | ||
| 608 | |||
| 609 | 151828758 | double getArea(TXsect* xsect, double y, double wSlot) | |
| 610 | // | ||
| 611 | // Input: xsect = ptr. to conduit cross section | ||
| 612 | // y = flow depth (ft) | ||
| 613 | // Output: returns flow area (ft2) | ||
| 614 | // Purpose: computes area of flow cross-section in a conduit. | ||
| 615 | // | ||
| 616 | { | ||
| 617 |
2/2✓ Branch 0 taken 4748085 times.
✓ Branch 1 taken 147080673 times.
|
151828758 | if ( y >= xsect->yFull ) return xsect->aFull + (y - xsect->yFull) * wSlot; |
| 618 | 147080673 | return xsect_getAofY(xsect, y); | |
| 619 | } | ||
| 620 | |||
| 621 | //============================================================================= | ||
| 622 | |||
| 623 | 101219172 | double getHydRad(TXsect* xsect, double y) | |
| 624 | // | ||
| 625 | // Input: xsect = ptr. to conduit cross section | ||
| 626 | // y = flow depth (ft) | ||
| 627 | // Output: returns hydraulic radius (ft) | ||
| 628 | // Purpose: computes hydraulic radius of flow cross-section in a conduit. | ||
| 629 | // | ||
| 630 | { | ||
| 631 |
2/2✓ Branch 0 taken 1978120 times.
✓ Branch 1 taken 99241052 times.
|
101219172 | if (y >= xsect->yFull) return xsect->rFull; |
| 632 | 99241052 | return xsect_getRofY(xsect, y); | |
| 633 | } | ||
| 634 | |||
| 635 | //============================================================================= | ||
| 636 | |||
| 637 | 43768447 | double checkNormalFlow(int j, double q, double y1, double y2, double a1, | |
| 638 | double r1) | ||
| 639 | // | ||
| 640 | // Input: j = link index | ||
| 641 | // q = link flow found from dynamic wave equations (cfs) | ||
| 642 | // y1 = flow depth at upstream end (ft) | ||
| 643 | // y2 = flow depth at downstream end (ft) | ||
| 644 | // a1 = flow area at upstream end (ft2) | ||
| 645 | // r1 = hyd. radius at upstream end (ft) | ||
| 646 | // Output: returns modifed flow in link (cfs) | ||
| 647 | // Purpose: checks if flow in link should be replaced by normal flow. | ||
| 648 | // | ||
| 649 | { | ||
| 650 | 43768447 | int check = FALSE; | |
| 651 | 43768447 | int k = Link[j].subIndex; | |
| 652 | 43768447 | int n1 = Link[j].node1; | |
| 653 | 43768447 | int n2 = Link[j].node2; | |
| 654 |
4/4✓ Branch 0 taken 43768030 times.
✓ Branch 1 taken 417 times.
✓ Branch 2 taken 1612198 times.
✓ Branch 3 taken 42155832 times.
|
43768447 | int hasOutfall = (Node[n1].type == OUTFALL || Node[n2].type == OUTFALL); |
| 655 | double qNorm; | ||
| 656 | double f1; | ||
| 657 | |||
| 658 | // --- check if water surface slope < conduit slope | ||
| 659 |
2/6✓ Branch 0 taken 43768447 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 43768447 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
43768447 | if ( NormalFlowLtd == SLOPE || NormalFlowLtd == BOTH || hasOutfall ) |
| 660 | { | ||
| 661 |
2/2✓ Branch 0 taken 27170206 times.
✓ Branch 1 taken 16598241 times.
|
43768447 | if ( y1 < y2) check = TRUE; |
| 662 | } | ||
| 663 | |||
| 664 | // --- check if Fr >= 1.0 at upstream end of conduit | ||
| 665 |
6/8✓ Branch 0 taken 16598241 times.
✓ Branch 1 taken 27170206 times.
✓ Branch 2 taken 16598241 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 16598241 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 15284739 times.
✓ Branch 7 taken 1313502 times.
|
43768447 | if ( !check && (NormalFlowLtd == FROUDE || NormalFlowLtd == BOTH) && |
| 666 | !hasOutfall ) | ||
| 667 | { | ||
| 668 |
3/4✓ Branch 0 taken 15284739 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 15200934 times.
✓ Branch 3 taken 83805 times.
|
15284739 | if ( y1 > FUDGE && y2 > FUDGE ) |
| 669 | { | ||
| 670 | 15200934 | f1 = link_getFroude(j, q/a1, y1); | |
| 671 |
2/2✓ Branch 0 taken 2917754 times.
✓ Branch 1 taken 12283180 times.
|
15200934 | if ( f1 >= 1.0 ) check = TRUE; |
| 672 | } | ||
| 673 | } | ||
| 674 | |||
| 675 | // --- check if normal flow < dynamic flow | ||
| 676 |
2/2✓ Branch 0 taken 30087960 times.
✓ Branch 1 taken 13680487 times.
|
43768447 | if ( check ) |
| 677 | { | ||
| 678 | 30087960 | qNorm = Conduit[k].beta * a1 * pow(r1, 2./3.); | |
| 679 |
2/2✓ Branch 0 taken 21691156 times.
✓ Branch 1 taken 8396804 times.
|
30087960 | if ( qNorm < q ) |
| 680 | { | ||
| 681 | 21691156 | Link[j].normalFlow = TRUE; | |
| 682 | 21691156 | return qNorm; | |
| 683 | } | ||
| 684 | } | ||
| 685 | 22077291 | return q; | |
| 686 | } | ||
| 687 |