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 | 9724625 | 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 | 9724625 | TXsect* xsect = &Link[j].xsect; // ptr. to conduit's cross section data | |
| 91 | 9724625 | char isFull = FALSE; // TRUE if conduit flowing full | |
| 92 | 9724625 | char isClosed = FALSE; // TRUE if conduit closed | |
| 93 | |||
| 94 | |||
| 95 | |||
| 96 | // --- adjust isClosed status by any control action | ||
| 97 |
2/2✓ Branch 0 taken 26014 times.
✓ Branch 1 taken 9698611 times.
|
9724625 | if ( Link[j].setting == 0 ) isClosed = TRUE; |
| 98 | |||
| 99 | // --- get flow from last time step & previous iteration | ||
| 100 | 9724625 | k = Link[j].subIndex; | |
| 101 | 9724625 | barrels = Conduit[k].barrels; | |
| 102 | 9724625 | qOld = Link[j].oldFlow / barrels; | |
| 103 | 9724625 | qLast = Conduit[k].q1; | |
| 104 | 9724625 | Conduit[k].evapLossRate = 0.0; | |
| 105 | 9724625 | Conduit[k].seepLossRate = 0.0; | |
| 106 | |||
| 107 | // --- get most current heads at upstream and downstream ends of conduit | ||
| 108 | 9724625 | n1 = Link[j].node1; | |
| 109 | 9724625 | n2 = Link[j].node2; | |
| 110 | 9724625 | z1 = Node[n1].invertElev + Link[j].offset1; | |
| 111 | 9724625 | z2 = Node[n2].invertElev + Link[j].offset2; | |
| 112 | 9724625 | h1 = Node[n1].newDepth + Node[n1].invertElev; | |
| 113 | 9724625 | h2 = Node[n2].newDepth + Node[n2].invertElev; | |
| 114 |
2/2✓ Branch 0 taken 9605974 times.
✓ Branch 1 taken 118651 times.
|
9724625 | h1 = MAX(h1, z1); |
| 115 |
2/2✓ Branch 0 taken 9135908 times.
✓ Branch 1 taken 588717 times.
|
9724625 | 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 | 9724625 | y1 = h1 - z1; | |
| 120 | 9724625 | y2 = h2 - z2; | |
| 121 |
2/2✓ Branch 0 taken 9230709 times.
✓ Branch 1 taken 493916 times.
|
9724625 | y1 = MAX(y1, FUDGE); |
| 122 |
2/2✓ Branch 0 taken 8605676 times.
✓ Branch 1 taken 1118949 times.
|
9724625 | y2 = MAX(y2, FUDGE); |
| 123 | |||
| 124 | // --- flow depths can't exceed full depth of conduit if slot not used | ||
| 125 |
2/2✓ Branch 0 taken 9716133 times.
✓ Branch 1 taken 8492 times.
|
9724625 | if ( SurchargeMethod != SLOT ) |
| 126 | { | ||
| 127 |
2/2✓ Branch 0 taken 8887375 times.
✓ Branch 1 taken 828758 times.
|
9716133 | y1 = MIN(y1, xsect->yFull); |
| 128 |
2/2✓ Branch 0 taken 8473656 times.
✓ Branch 1 taken 1242477 times.
|
9716133 | y2 = MIN(y2, xsect->yFull); |
| 129 | } | ||
| 130 | |||
| 131 | // -- get area from solution at previous time step | ||
| 132 | 9724625 | aOld = Conduit[k].a2; | |
| 133 |
2/2✓ Branch 0 taken 8824949 times.
✓ Branch 1 taken 899676 times.
|
9724625 | aOld = MAX(aOld, FUDGE); |
| 134 | |||
| 135 | // --- use Courant-modified length instead of conduit's actual length | ||
| 136 | 9724625 | length = Conduit[k].modLength; | |
| 137 | |||
| 138 | // --- find surface area contributions to upstream and downstream nodes | ||
| 139 | // based on previous iteration's flow estimate | ||
| 140 | 9724625 | findSurfArea(j, qLast, length, &h1, &h2, &y1, &y2); | |
| 141 | |||
| 142 | // --- compute area at each end of conduit & hyd. radius at upstream end | ||
| 143 | 9724625 | wSlot = getSlotWidth(xsect, y1); | |
| 144 | 9724625 | a1 = getArea(xsect, y1, wSlot); | |
| 145 | 9724625 | r1 = getHydRad(xsect, y1); | |
| 146 | 9724625 | wSlot = getSlotWidth(xsect, y2); | |
| 147 | 9724625 | a2 = getArea(xsect, y2, wSlot); | |
| 148 | |||
| 149 | // --- compute area & hyd. radius at midpoint | ||
| 150 | 9724625 | yMid = 0.5 * (y1 + y2); | |
| 151 | 9724625 | wSlot = getSlotWidth(xsect, yMid); | |
| 152 | 9724625 | aMid = getArea(xsect, yMid, wSlot); | |
| 153 | 9724625 | 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 862171 times.
✓ Branch 1 taken 8862454 times.
|
9724625 | if ( y1 >= xsect->yFull && |
| 162 |
2/2✓ Branch 0 taken 554967 times.
✓ Branch 1 taken 307204 times.
|
862171 | 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 9480847 times.
✓ Branch 1 taken 243778 times.
|
9724625 | if ( Link[j].flowClass == DRY || |
| 166 |
2/2✓ Branch 0 taken 9231939 times.
✓ Branch 1 taken 248908 times.
|
9480847 | Link[j].flowClass == UP_DRY || |
| 167 |
3/4✓ Branch 0 taken 9231939 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 9205936 times.
✓ Branch 3 taken 26003 times.
|
9231939 | Link[j].flowClass == DN_DRY || |
| 168 |
2/2✓ Branch 0 taken 923447 times.
✓ Branch 1 taken 8282489 times.
|
9205936 | isClosed || |
| 169 | aMid <= FUDGE ) | ||
| 170 | { | ||
| 171 | 1442136 | Conduit[k].a1 = 0.5 * (a1 + a2); | |
| 172 | 1442136 | Conduit[k].q1 = 0.0;; | |
| 173 | 1442136 | Conduit[k].q2 = 0.0; | |
| 174 | 1442136 | Link[j].dqdh = GRAVITY * dt * aMid / length * barrels; | |
| 175 | 1442136 | Link[j].froude = 0.0; | |
| 176 |
1/2✓ Branch 0 taken 1442136 times.
✗ Branch 1 not taken.
|
1442136 | Link[j].newDepth = MIN(yMid, Link[j].xsect.yFull); |
| 177 | 1442136 | Link[j].newVolume = Conduit[k].a1 * link_getLength(j) * barrels; | |
| 178 | 1442136 | Link[j].newFlow = 0.0; | |
| 179 | 1442136 | return; | |
| 180 | } | ||
| 181 | |||
| 182 | // --- compute velocity from last flow estimate | ||
| 183 | 8282489 | v = qLast / aMid; | |
| 184 |
3/4✓ Branch 0 taken 4 times.
✓ Branch 1 taken 8282485 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
|
8282489 | if ( fabs(v) > MAXVELOCITY ) v = MAXVELOCITY * SGN(qLast); |
| 185 | |||
| 186 | // --- compute Froude No. | ||
| 187 | 8282489 | Link[j].froude = link_getFroude(j, v, yMid); | |
| 188 |
2/2✓ Branch 0 taken 7521947 times.
✓ Branch 1 taken 760542 times.
|
8282489 | if ( Link[j].flowClass == SUBCRITICAL && |
| 189 |
2/2✓ Branch 0 taken 2545910 times.
✓ Branch 1 taken 4976037 times.
|
7521947 | Link[j].froude > 1.0 ) Link[j].flowClass = SUPCRITICAL; |
| 190 | |||
| 191 | // --- find inertial damping factor (sigma) | ||
| 192 |
2/2✓ Branch 0 taken 2511867 times.
✓ Branch 1 taken 5770622 times.
|
8282489 | if ( Link[j].froude <= 0.5 ) sigma = 1.0; |
| 193 |
2/2✓ Branch 0 taken 3129815 times.
✓ Branch 1 taken 2640807 times.
|
5770622 | else if ( Link[j].froude >= 1.0 ) sigma = 0.0; |
| 194 | 2640807 | 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 | 8282489 | rho = 1.0; | |
| 199 |
6/6✓ Branch 0 taken 7727522 times.
✓ Branch 1 taken 554967 times.
✓ Branch 2 taken 7656540 times.
✓ Branch 3 taken 70982 times.
✓ Branch 4 taken 7603366 times.
✓ Branch 5 taken 53174 times.
|
8282489 | if ( !isFull && qLast > 0.0 && h1 >= h2 ) rho = sigma; |
| 200 | 8282489 | aWtd = a1 + (aMid - a1) * rho; | |
| 201 | 8282489 | rWtd = r1 + (rMid - r1) * rho; | |
| 202 | |||
| 203 | // --- determine how much inertial damping to apply | ||
| 204 |
2/2✓ Branch 0 taken 4409847 times.
✓ Branch 1 taken 3872642 times.
|
8282489 | if ( InertDamping == NO_DAMPING ) sigma = 1.0; |
| 205 |
2/2✓ Branch 0 taken 216046 times.
✓ Branch 1 taken 3656596 times.
|
3872642 | 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 554967 times.
✓ Branch 1 taken 7727522 times.
✓ Branch 3 taken 431264 times.
✓ Branch 4 taken 123703 times.
|
8282489 | if ( isFull && !xsect_isOpen(xsect->type) ) sigma = 0.0; |
| 209 | |||
| 210 | // --- compute terms of momentum eqn.: | ||
| 211 | // --- 1. friction slope term | ||
| 212 |
3/4✓ Branch 0 taken 23056 times.
✓ Branch 1 taken 8259433 times.
✓ Branch 2 taken 23056 times.
✗ Branch 3 not taken.
|
8282489 | if ( xsect->type == FORCE_MAIN && isFull ) |
| 213 | 23056 | dq1 = dt * forcemain_getFricSlope(j, fabs(v), rMid); | |
| 214 | 8259433 | else dq1 = dt * Conduit[k].roughFactor / pow(rWtd, 1.33333) * fabs(v); | |
| 215 | |||
| 216 | // --- 2. energy slope term | ||
| 217 | 8282489 | dq2 = dt * GRAVITY * aWtd * (h2 - h1) / length; | |
| 218 | |||
| 219 | // --- 3 & 4. inertial terms | ||
| 220 | 8282489 | dq3 = 0.0; | |
| 221 | 8282489 | dq4 = 0.0; | |
| 222 |
2/2✓ Branch 0 taken 5716640 times.
✓ Branch 1 taken 2565849 times.
|
8282489 | if ( sigma > 0.0 ) |
| 223 | { | ||
| 224 | 5716640 | dq3 = 2.0 * v * (aMid - aOld) * sigma; | |
| 225 | 5716640 | dq4 = dt * v * v * (a2 - a1) / length * sigma; | |
| 226 | } | ||
| 227 | |||
| 228 | // --- 5. local losses term | ||
| 229 | 8282489 | dq5 = 0.0; | |
| 230 |
2/2✓ Branch 0 taken 158231 times.
✓ Branch 1 taken 8124258 times.
|
8282489 | if ( Conduit[k].hasLosses ) |
| 231 | { | ||
| 232 | 158231 | 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 | 8282489 | dq6 = link_getLossRate(j, DW, qLast, dt) * 2.5 * dt * v / link_getLength(j); | |
| 237 | |||
| 238 | // --- combine terms to find new conduit flow | ||
| 239 | 8282489 | denom = 1.0 + dq1 + dq5; | |
| 240 | 8282489 | q = (qOld - dq2 + dq3 + dq4 + dq6) / denom; | |
| 241 | |||
| 242 | // --- compute derivative of flow w.r.t. head | ||
| 243 | 8282489 | Link[j].dqdh = 1.0 / denom * GRAVITY * dt * aWtd / length * barrels; | |
| 244 | |||
| 245 | // --- check if any flow limitation applies | ||
| 246 | 8282489 | Link[j].inletControl = FALSE; | |
| 247 | 8282489 | Link[j].normalFlow = FALSE; | |
| 248 |
2/2✓ Branch 0 taken 8229253 times.
✓ Branch 1 taken 53236 times.
|
8282489 | if ( q > 0.0 ) |
| 249 | { | ||
| 250 | // --- check for inlet controlled culvert flow | ||
| 251 |
4/4✓ Branch 0 taken 35999 times.
✓ Branch 1 taken 8193254 times.
✓ Branch 2 taken 34053 times.
✓ Branch 3 taken 1946 times.
|
8229253 | 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 8195200 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 7380013 times.
✓ Branch 3 taken 815187 times.
|
8195200 | else if (NormalFlowLtd != NEITHER && y1 < Link[j].xsect.yFull && |
| 256 |
2/2✓ Branch 0 taken 3123168 times.
✓ Branch 1 taken 4256845 times.
|
7380013 | ( Link[j].flowClass == SUBCRITICAL || |
| 257 |
2/2✓ Branch 0 taken 2434956 times.
✓ Branch 1 taken 688212 times.
|
3123168 | Link[j].flowClass == SUPCRITICAL )) |
| 258 | 6691801 | 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 4624540 times.
✓ Branch 1 taken 3657949 times.
|
8282489 | if ( steps > 0 ) |
| 264 | { | ||
| 265 | 4624540 | q = (1.0 - omega) * qLast + omega * q; | |
| 266 |
4/4✓ Branch 0 taken 531 times.
✓ Branch 1 taken 4624009 times.
✓ Branch 2 taken 229 times.
✓ Branch 3 taken 302 times.
|
4624540 | 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 8196489 times.
|
8282489 | 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 14074 times.
✓ Branch 2 taken 8268415 times.
|
8282489 | 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 7614859 times.
✓ Branch 1 taken 667630 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 7614859 times.
|
8282489 | if( q > FUDGE && Node[n1].newDepth <= FUDGE ) q = FUDGE; |
| 281 |
4/4✓ Branch 0 taken 38300 times.
✓ Branch 1 taken 8244189 times.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 38299 times.
|
8282489 | if( q < -FUDGE && Node[n2].newDepth <= FUDGE ) q = -FUDGE; |
| 282 | |||
| 283 | // --- save new values of area, flow, depth, & volume | ||
| 284 | 8282489 | Conduit[k].a1 = aMid; | |
| 285 | 8282489 | Conduit[k].q1 = q; | |
| 286 | 8282489 | Conduit[k].q2 = q; | |
| 287 |
2/2✓ Branch 0 taken 8278624 times.
✓ Branch 1 taken 3865 times.
|
8282489 | Link[j].newDepth = MIN(yMid, xsect->yFull); |
| 288 | 8282489 | aMid = (a1 + a2) / 2.0; | |
| 289 | // aMid = MIN(aMid, xsect->aFull); //Slot can have aMid > aFull | ||
| 290 | 8282489 | Conduit[k].fullState = link_getFullState(a1, a2, xsect->aFull); | |
| 291 | 8282489 | Link[j].newVolume = aMid * link_getLength(j) * barrels; | |
| 292 | 8282489 | Link[j].newFlow = q * barrels; | |
| 293 | } | ||
| 294 | |||
| 295 | //============================================================================= | ||
| 296 | |||
| 297 | 9169658 | 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 | 9169658 | n1 = Link[j].node1; | |
| 320 | 9169658 | n2 = Link[j].node2; | |
| 321 | |||
| 322 | // --- get upstream & downstream conduit invert offsets | ||
| 323 | 9169658 | z1 = Link[j].offset1; | |
| 324 | 9169658 | z2 = Link[j].offset2; | |
| 325 | |||
| 326 | // --- base offset of an outfall conduit on outfall's depth | ||
| 327 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 9169658 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
9169658 | if ( Node[n1].type == OUTFALL ) z1 = MAX(0.0, (z1 - Node[n1].newDepth)); |
| 328 |
3/4✓ Branch 0 taken 1931261 times.
✓ Branch 1 taken 7238397 times.
✓ Branch 2 taken 1931261 times.
✗ Branch 3 not taken.
|
9169658 | if ( Node[n2].type == OUTFALL ) z2 = MAX(0.0, (z2 - Node[n2].newDepth)); |
| 329 | |||
| 330 | // --- default class is SUBCRITICAL | ||
| 331 | 9169658 | flowClass = SUBCRITICAL; | |
| 332 | 9169658 | *fasnh = 1.0; | |
| 333 | |||
| 334 | // --- case where both ends of conduit are wet | ||
| 335 |
4/4✓ Branch 0 taken 8675742 times.
✓ Branch 1 taken 493916 times.
✓ Branch 2 taken 7800571 times.
✓ Branch 3 taken 875171 times.
|
9169658 | if ( y1 > FUDGE && y2 > FUDGE ) |
| 336 | { | ||
| 337 |
2/2✓ Branch 0 taken 23666 times.
✓ Branch 1 taken 7776905 times.
|
7800571 | 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 1 time.
✓ Branch 1 taken 23665 times.
|
23666 | if ( z1 > 0.0 ) |
| 343 | { | ||
| 344 | 1 | *yN = link_getYnorm(j, fabs(q)); | |
| 345 | 1 | *yC = link_getYcrit(j, fabs(q)); | |
| 346 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | ycMin = MIN(*yN, *yC); |
| 347 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
|
1 | 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 515816 times.
✓ Branch 1 taken 7261089 times.
|
7776905 | if ( z2 > 0.0 ) |
| 358 | { | ||
| 359 | 515816 | *yN = link_getYnorm(j, fabs(q)); | |
| 360 | 515816 | *yC = link_getYcrit(j, fabs(q)); | |
| 361 |
2/2✓ Branch 0 taken 221688 times.
✓ Branch 1 taken 294128 times.
|
515816 | ycMin = MIN(*yN, *yC); |
| 362 |
2/2✓ Branch 0 taken 349568 times.
✓ Branch 1 taken 166248 times.
|
515816 | ycMax = MAX(*yN, *yC); |
| 363 |
2/2✓ Branch 0 taken 303529 times.
✓ Branch 1 taken 212287 times.
|
515816 | if ( y2 < ycMin ) flowClass = DN_CRITICAL; |
| 364 |
2/2✓ Branch 0 taken 153328 times.
✓ Branch 1 taken 58959 times.
|
212287 | else if ( y2 < ycMax ) |
| 365 | { | ||
| 366 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 153328 times.
|
153328 | if ( ycMax - ycMin < FUDGE ) *fasnh = 0.0; |
| 367 | 153328 | 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 493916 times.
✓ Branch 1 taken 875171 times.
✓ Branch 2 taken 243778 times.
✓ Branch 3 taken 250138 times.
|
1369087 | 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 250138 times.
✓ Branch 1 taken 875171 times.
|
1125309 | 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 248908 times.
✓ Branch 1 taken 1230 times.
|
250138 | 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 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1230 times.
|
1230 | else if ( z1 > 0.0 ) |
| 388 | { | ||
| 389 | ✗ | *yN = link_getYnorm(j, fabs(q)); | |
| 390 | ✗ | *yC = link_getYcrit(j, fabs(q)); | |
| 391 | ✗ | 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 875171 times.
|
875171 | 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 540301 times.
✓ Branch 1 taken 334870 times.
|
875171 | else if ( z2 > 0.0 ) |
| 406 | { | ||
| 407 | 540301 | *yN = link_getYnorm(j, fabs(q)); | |
| 408 | 540301 | *yC = link_getYcrit(j, fabs(q)); | |
| 409 | 540301 | flowClass = DN_CRITICAL; | |
| 410 | } | ||
| 411 | } | ||
| 412 | 9169658 | return flowClass; | |
| 413 | } | ||
| 414 | |||
| 415 | //============================================================================= | ||
| 416 | |||
| 417 | 9724625 | 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 | 9724625 | double surfArea1 = 0.0; // surface area at upstream node (ft2) | |
| 439 | 9724625 | 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 | 9724625 | double fasnh = 1.0; // fraction between norm. & crit. depth | |
| 444 | 9724625 | TXsect* xsect = &Link[j].xsect; // pointer to cross-section data | |
| 445 | |||
| 446 | // --- get node indexes & current flow depths | ||
| 447 | 9724625 | n1 = Link[j].node1; | |
| 448 | 9724625 | n2 = Link[j].node2; | |
| 449 | 9724625 | flowDepth1 = *y1; | |
| 450 | 9724625 | flowDepth2 = *y2; | |
| 451 | |||
| 452 | 9724625 | normalDepth = (flowDepth1 + flowDepth2) / 2.0; | |
| 453 | 9724625 | criticalDepth = normalDepth; | |
| 454 | |||
| 455 | // --- find conduit's flow classification | ||
| 456 | 9724625 | fullDepth = xsect->yFull; | |
| 457 |
4/4✓ Branch 0 taken 862171 times.
✓ Branch 1 taken 8862454 times.
✓ Branch 2 taken 554967 times.
✓ Branch 3 taken 307204 times.
|
9724625 | if (flowDepth1 >= fullDepth && flowDepth2 >= fullDepth) |
| 458 | { | ||
| 459 | 554967 | Link[j].flowClass = SUBCRITICAL; | |
| 460 | } | ||
| 461 | 9169658 | 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 |
4/7✓ Branch 0 taken 8388109 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 843830 times.
✓ Branch 3 taken 248908 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 243778 times.
✗ Branch 6 not taken.
|
9724625 | switch ( Link[j].flowClass ) |
| 466 | { | ||
| 467 | 8388109 | case SUBCRITICAL: | |
| 468 | 8388109 | flowDepthMid = 0.5 * (flowDepth1 + flowDepth2); | |
| 469 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 8388109 times.
|
8388109 | if ( flowDepthMid < FUDGE ) flowDepthMid = FUDGE; |
| 470 | 8388109 | width1 = getWidth(xsect, flowDepth1); | |
| 471 | 8388109 | width2 = getWidth(xsect, flowDepth2); | |
| 472 | 8388109 | widthMid = getWidth(xsect, flowDepthMid); | |
| 473 | 8388109 | surfArea1 = (width1 + widthMid) * length / 4.; | |
| 474 | 8388109 | surfArea2 = (widthMid + width2) * length / 4. * fasnh; | |
| 475 | 8388109 | break; | |
| 476 | |||
| 477 | ✗ | case UP_CRITICAL: | |
| 478 | ✗ | flowDepth1 = criticalDepth; | |
| 479 | ✗ | if ( normalDepth < criticalDepth ) flowDepth1 = normalDepth; | |
| 480 | ✗ | flowDepth1 = MAX(flowDepth1, FUDGE); | |
| 481 | ✗ | *h1 = Node[n1].invertElev + Link[j].offset1 + flowDepth1; | |
| 482 | ✗ | flowDepthMid = 0.5 * (flowDepth1 + flowDepth2); | |
| 483 | ✗ | if ( flowDepthMid < FUDGE ) flowDepthMid = FUDGE; | |
| 484 | ✗ | width2 = getWidth(xsect, flowDepth2); | |
| 485 | ✗ | widthMid = getWidth(xsect, flowDepthMid); | |
| 486 | ✗ | surfArea2 = (widthMid + width2) * length * 0.5; | |
| 487 | ✗ | break; | |
| 488 | |||
| 489 | 843830 | case DN_CRITICAL: | |
| 490 | 843830 | flowDepth2 = criticalDepth; | |
| 491 |
2/2✓ Branch 0 taken 441265 times.
✓ Branch 1 taken 402565 times.
|
843830 | if ( normalDepth < criticalDepth ) flowDepth2 = normalDepth; |
| 492 |
2/2✓ Branch 0 taken 753628 times.
✓ Branch 1 taken 90202 times.
|
843830 | flowDepth2 = MAX(flowDepth2, FUDGE); |
| 493 | 843830 | *h2 = Node[n2].invertElev + Link[j].offset2 + flowDepth2; | |
| 494 | 843830 | width1 = getWidth(xsect, flowDepth1); | |
| 495 | 843830 | flowDepthMid = 0.5 * (flowDepth1 + flowDepth2); | |
| 496 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 843830 times.
|
843830 | if ( flowDepthMid < FUDGE ) flowDepthMid = FUDGE; |
| 497 | 843830 | widthMid = getWidth(xsect, flowDepthMid); | |
| 498 | 843830 | surfArea1 = (width1 + widthMid) * length * 0.5; | |
| 499 | 843830 | break; | |
| 500 | |||
| 501 | 248908 | case UP_DRY: | |
| 502 | 248908 | flowDepth1 = FUDGE; | |
| 503 | 248908 | flowDepthMid = 0.5 * (flowDepth1 + flowDepth2); | |
| 504 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 248908 times.
|
248908 | if ( flowDepthMid < FUDGE ) flowDepthMid = FUDGE; |
| 505 | 248908 | width1 = getWidth(xsect, flowDepth1); | |
| 506 | 248908 | width2 = getWidth(xsect, flowDepth2); | |
| 507 | 248908 | widthMid = getWidth(xsect, flowDepthMid); | |
| 508 | |||
| 509 | // --- assign avg. surface area of downstream half of conduit | ||
| 510 | // to the downstream node | ||
| 511 | 248908 | 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 248905 times.
✓ Branch 1 taken 3 times.
|
248908 | if ( Link[j].offset1 <= 0.0 ) |
| 516 | { | ||
| 517 | 248905 | surfArea1 = (width1 + widthMid) * length / 4.; | |
| 518 | } | ||
| 519 | 248908 | 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 | 243778 | case DRY: | |
| 542 | 243778 | surfArea1 = FUDGE * length / 2.0; | |
| 543 | 243778 | surfArea2 = surfArea1; | |
| 544 | 243778 | break; | |
| 545 | } | ||
| 546 | 9724625 | Link[j].surfArea1 = surfArea1; | |
| 547 | 9724625 | Link[j].surfArea2 = surfArea2; | |
| 548 | 9724625 | *y1 = flowDepth1; | |
| 549 | 9724625 | *y2 = flowDepth2; | |
| 550 | 9724625 | } | |
| 551 | |||
| 552 | //============================================================================= | ||
| 553 | |||
| 554 | 158231 | 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 | 158231 | double losses = 0.0; | |
| 566 | 158231 | q = fabs(q); | |
| 567 |
1/2✓ Branch 0 taken 158231 times.
✗ Branch 1 not taken.
|
158231 | if ( a1 > FUDGE ) losses += Link[j].cLossInlet * (q/a1); |
| 568 |
2/2✓ Branch 0 taken 157532 times.
✓ Branch 1 taken 699 times.
|
158231 | if ( a2 > FUDGE ) losses += Link[j].cLossOutlet * (q/a2); |
| 569 |
1/2✓ Branch 0 taken 158231 times.
✗ Branch 1 not taken.
|
158231 | if ( aMid > FUDGE ) losses += Link[j].cLossAvg * (q/aMid); |
| 570 | 158231 | return losses; | |
| 571 | } | ||
| 572 | |||
| 573 | //============================================================================= | ||
| 574 | |||
| 575 | 56772586 | double getSlotWidth(TXsect* xsect, double y) | |
| 576 | { | ||
| 577 | 56772586 | double yNorm = y / xsect->yFull; | |
| 578 | |||
| 579 | // --- return 0.0 if slot surcharge method not used | ||
| 580 |
3/4✓ Branch 0 taken 50949 times.
✓ Branch 1 taken 56721637 times.
✓ Branch 3 taken 50949 times.
✗ Branch 4 not taken.
|
56772586 | if (SurchargeMethod != SLOT || xsect_isOpen(xsect->type) || |
| 581 |
2/2✓ Branch 0 taken 27839 times.
✓ Branch 1 taken 23110 times.
|
56772586 | yNorm < CrownCutoff) return 0.0; |
| 582 | |||
| 583 | // --- for depth > 1.78 * pipe depth, slot width = 1% of max. width | ||
| 584 |
2/2✓ Branch 0 taken 8272 times.
✓ Branch 1 taken 14838 times.
|
23110 | if (yNorm > 1.78) return 0.01 * xsect->wMax; |
| 585 | |||
| 586 | // --- otherwise use the Sjoberg formula | ||
| 587 | 14838 | return xsect->wMax * 0.5423 * exp(-pow(yNorm, 2.4)); | |
| 588 | } | ||
| 589 | |||
| 590 | //============================================================================= | ||
| 591 | |||
| 592 | 27598711 | 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 | 27598711 | double wSlot = getSlotWidth(xsect, y); | |
| 601 |
2/2✓ Branch 0 taken 11555 times.
✓ Branch 1 taken 27587156 times.
|
27598711 | if (wSlot > 0.0) return wSlot; |
| 602 |
4/4✓ Branch 0 taken 2893937 times.
✓ Branch 1 taken 24693219 times.
✓ Branch 3 taken 2070867 times.
✓ Branch 4 taken 823070 times.
|
27587156 | if (y / xsect->yFull >= CrownCutoff && !xsect_isOpen(xsect->type)) |
| 603 | 2070867 | y = CrownCutoff * xsect->yFull; | |
| 604 | 27587156 | return xsect_getWofY(xsect, y); | |
| 605 | } | ||
| 606 | |||
| 607 | //============================================================================= | ||
| 608 | |||
| 609 | 29173875 | 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 2675830 times.
✓ Branch 1 taken 26498045 times.
|
29173875 | if ( y >= xsect->yFull ) return xsect->aFull + (y - xsect->yFull) * wSlot; |
| 618 | 26498045 | return xsect_getAofY(xsect, y); | |
| 619 | } | ||
| 620 | |||
| 621 | //============================================================================= | ||
| 622 | |||
| 623 | 19449250 | 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 1417252 times.
✓ Branch 1 taken 18031998 times.
|
19449250 | if (y >= xsect->yFull) return xsect->rFull; |
| 632 | 18031998 | return xsect_getRofY(xsect, y); | |
| 633 | } | ||
| 634 | |||
| 635 | //============================================================================= | ||
| 636 | |||
| 637 | 6691801 | 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 | 6691801 | int check = FALSE; | |
| 651 | 6691801 | int k = Link[j].subIndex; | |
| 652 | 6691801 | int n1 = Link[j].node1; | |
| 653 | 6691801 | int n2 = Link[j].node2; | |
| 654 |
3/4✓ Branch 0 taken 6691801 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1646912 times.
✓ Branch 3 taken 5044889 times.
|
6691801 | 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 |
3/6✓ Branch 0 taken 6584010 times.
✓ Branch 1 taken 107791 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6584010 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
6691801 | if ( NormalFlowLtd == SLOPE || NormalFlowLtd == BOTH || hasOutfall ) |
| 660 | { | ||
| 661 |
2/2✓ Branch 0 taken 4216347 times.
✓ Branch 1 taken 2475454 times.
|
6691801 | if ( y1 < y2) check = TRUE; |
| 662 | } | ||
| 663 | |||
| 664 | // --- check if Fr >= 1.0 at upstream end of conduit | ||
| 665 |
7/8✓ Branch 0 taken 2475454 times.
✓ Branch 1 taken 4216347 times.
✓ Branch 2 taken 2475454 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2421024 times.
✓ Branch 5 taken 54430 times.
✓ Branch 6 taken 1096231 times.
✓ Branch 7 taken 1324793 times.
|
6691801 | if ( !check && (NormalFlowLtd == FROUDE || NormalFlowLtd == BOTH) && |
| 666 | !hasOutfall ) | ||
| 667 | { | ||
| 668 |
3/4✓ Branch 0 taken 1096231 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1091019 times.
✓ Branch 3 taken 5212 times.
|
1096231 | if ( y1 > FUDGE && y2 > FUDGE ) |
| 669 | { | ||
| 670 | 1091019 | f1 = link_getFroude(j, q/a1, y1); | |
| 671 |
2/2✓ Branch 0 taken 537911 times.
✓ Branch 1 taken 553108 times.
|
1091019 | if ( f1 >= 1.0 ) check = TRUE; |
| 672 | } | ||
| 673 | } | ||
| 674 | |||
| 675 | // --- check if normal flow < dynamic flow | ||
| 676 |
2/2✓ Branch 0 taken 4754258 times.
✓ Branch 1 taken 1937543 times.
|
6691801 | if ( check ) |
| 677 | { | ||
| 678 | 4754258 | qNorm = Conduit[k].beta * a1 * pow(r1, 2./3.); | |
| 679 |
2/2✓ Branch 0 taken 3760531 times.
✓ Branch 1 taken 993727 times.
|
4754258 | if ( qNorm < q ) |
| 680 | { | ||
| 681 | 3760531 | Link[j].normalFlow = TRUE; | |
| 682 | 3760531 | return qNorm; | |
| 683 | } | ||
| 684 | } | ||
| 685 | 2931270 | return q; | |
| 686 | } | ||
| 687 |