dynwave.c
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | //----------------------------------------------------------------------------- | ||
| 2 | // dynwave.c | ||
| 3 | // | ||
| 4 | // Project: EPA SWMM5 | ||
| 5 | // Version: 5.2 | ||
| 6 | // Date: 07/13/23 (Build 5.2.4) | ||
| 7 | // Author: L. Rossman | ||
| 8 | // M. Tryby (EPA) | ||
| 9 | // R. Dickinson (CDM) | ||
| 10 | // | ||
| 11 | // Dynamic wave flow routing functions. | ||
| 12 | // | ||
| 13 | // This module solves the dynamic wave flow routing equations using | ||
| 14 | // Picard Iterations (i.e., a method of successive approximations) | ||
| 15 | // to solve the explicit form of the continuity and momentum equations | ||
| 16 | // for conduits. | ||
| 17 | // | ||
| 18 | // Update History | ||
| 19 | // ============== | ||
| 20 | // Build 5.1.002: | ||
| 21 | // - Only non-ponded nodal surface area is saved for use in | ||
| 22 | // surcharge algorithm. | ||
| 23 | // Build 5.1.007: | ||
| 24 | // - Node losses added to node outflow variable instead of treated | ||
| 25 | // as a separate item when computing change in node flow volume. | ||
| 26 | // Build 5.1.008: | ||
| 27 | // - Module-specific constants moved here from project.c. | ||
| 28 | // - Support added for user-specified minimum variable time step. | ||
| 29 | // - Node crown elevations found here instead of in flowrout.c module. | ||
| 30 | // - OpenMP use to parallelize findLinkFlows() & findNodeDepths(). | ||
| 31 | // - Bug in finding complete list of capacity limited links fixed. | ||
| 32 | // Build 5.1.011: | ||
| 33 | // - Added test for failed memory allocation. | ||
| 34 | // - Fixed illegal array index bug for Ideal Pumps. | ||
| 35 | // Build 5.1.013: | ||
| 36 | // - Include omp.h protected against lack of compiler support for OpenMP. | ||
| 37 | // - SurchargeMethod option used to decide how node surcharging is handled. | ||
| 38 | // - Storage nodes allowed to pressurize if their surcharge depth > 0. | ||
| 39 | // - Minimum flow needed to compute a Courant time step modified. | ||
| 40 | // Build 5.1.014: | ||
| 41 | // - updateNodeFlows() modified to subtract conduit evap. and seepage losses | ||
| 42 | // from downstream node inflow instead of upstream node outflow. | ||
| 43 | // Build 5.1.015: | ||
| 44 | // - Roll back the 5.1.014 change for conduit losses in updateNodeFlows(). | ||
| 45 | // Build 5.2.0: | ||
| 46 | // - Support added for reporting most frequent non-converging links. | ||
| 47 | // Build 5.2.4: | ||
| 48 | // - Conduit evap+seepage outflow split evenly between outflow from | ||
| 49 | // conduit's upstream and non-outfall downstream nodes. | ||
| 50 | //----------------------------------------------------------------------------- | ||
| 51 | #define _CRT_SECURE_NO_DEPRECATE | ||
| 52 | |||
| 53 | #include <stdlib.h> | ||
| 54 | #include <math.h> | ||
| 55 | #include "headers.h" | ||
| 56 | |||
| 57 | //----------------------------------------------------------------------------- | ||
| 58 | // Constants | ||
| 59 | //----------------------------------------------------------------------------- | ||
| 60 | static const double MINTIMESTEP = 0.001; // min. time step (sec) | ||
| 61 | static const double OMEGA = 0.5; // under-relaxation parameter | ||
| 62 | static const double DEFAULT_SURFAREA = 12.566; // Min. nodal surface area (~4 ft diam.) | ||
| 63 | static const double DEFAULT_HEADTOL = 0.005; // Default head tolerance (ft) | ||
| 64 | static const double EXTRAN_CROWN_CUTOFF = 0.96; // crown cutoff for EXTRAN | ||
| 65 | static const double SLOT_CROWN_CUTOFF = 0.985257; // crown cutoff for SLOT | ||
| 66 | static const int DEFAULT_MAXTRIALS = 8; // Max. trials per time step | ||
| 67 | |||
| 68 | |||
| 69 | //----------------------------------------------------------------------------- | ||
| 70 | // Data Structures | ||
| 71 | //----------------------------------------------------------------------------- | ||
| 72 | typedef struct | ||
| 73 | { | ||
| 74 | char converged; // TRUE if iterations for a node done | ||
| 75 | double newSurfArea; // current surface area (ft2) | ||
| 76 | double oldSurfArea; // previous surface area (ft2) | ||
| 77 | double sumdqdh; // sum of dqdh from adjoining links | ||
| 78 | double dYdT; // change in depth w.r.t. time (ft/sec) | ||
| 79 | } TXnode; | ||
| 80 | |||
| 81 | //----------------------------------------------------------------------------- | ||
| 82 | // Shared Variables | ||
| 83 | //----------------------------------------------------------------------------- | ||
| 84 | static double VariableStep; // size of variable time step (sec) | ||
| 85 | static TXnode* Xnode; // extended nodal information | ||
| 86 | |||
| 87 | static double Omega; // actual under-relaxation parameter | ||
| 88 | static int Steps; // number of Picard iterations | ||
| 89 | |||
| 90 | //----------------------------------------------------------------------------- | ||
| 91 | // Function declarations | ||
| 92 | //----------------------------------------------------------------------------- | ||
| 93 | static void initRoutingStep(void); | ||
| 94 | static void initNodeStates(void); | ||
| 95 | static void findBypassedLinks(); | ||
| 96 | static void findLimitedLinks(); | ||
| 97 | |||
| 98 | static void findLinkFlows(double dt); | ||
| 99 | static int isTrueConduit(int link); | ||
| 100 | static void findNonConduitFlow(int link, double dt); | ||
| 101 | static void findNonConduitSurfArea(int link); | ||
| 102 | static double getModPumpFlow(int link, double q, double dt); | ||
| 103 | static void updateNodeFlows(int link); | ||
| 104 | static void updateConvergenceStats(); | ||
| 105 | |||
| 106 | static int findNodeDepths(double dt); | ||
| 107 | static void setNodeDepth(int node, double dt); | ||
| 108 | static double getFloodedDepth(int node, int canPond, double dV, double yNew, | ||
| 109 | double yMax, double dt); | ||
| 110 | |||
| 111 | static double getVariableStep(double maxStep); | ||
| 112 | static double getLinkStep(double tMin, int *minLink); | ||
| 113 | static double getNodeStep(double tMin, int *minNode); | ||
| 114 | |||
| 115 | //============================================================================= | ||
| 116 | |||
| 117 | 41 | void dynwave_init() | |
| 118 | // | ||
| 119 | // Input: none | ||
| 120 | // Output: none | ||
| 121 | // Purpose: initializes dynamic wave routing method. | ||
| 122 | // | ||
| 123 | { | ||
| 124 | int i, j; | ||
| 125 | double z; | ||
| 126 | |||
| 127 | 41 | VariableStep = 0.0; | |
| 128 | 41 | Xnode = (TXnode *) calloc(Nobjects[NODE], sizeof(TXnode)); | |
| 129 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 41 times.
|
41 | if ( Xnode == NULL ) |
| 130 | { | ||
| 131 | ✗ | report_writeErrorMsg(ERR_MEMORY, | |
| 132 | " Not enough memory for dynamic wave routing."); | ||
| 133 | ✗ | return; | |
| 134 | } | ||
| 135 | |||
| 136 | // --- initialize node surface areas & crown elev. | ||
| 137 |
2/2✓ Branch 0 taken 10094 times.
✓ Branch 1 taken 41 times.
|
10135 | for (i = 0; i < Nobjects[NODE]; i++ ) |
| 138 | { | ||
| 139 | 10094 | Xnode[i].newSurfArea = 0.0; | |
| 140 | 10094 | Xnode[i].oldSurfArea = 0.0; | |
| 141 | 10094 | Node[i].crownElev = Node[i].invertElev; | |
| 142 | } | ||
| 143 | |||
| 144 | // --- initialize links & update node crown elevations | ||
| 145 |
2/2✓ Branch 0 taken 10449 times.
✓ Branch 1 taken 41 times.
|
10490 | for (i = 0; i < Nobjects[LINK]; i++) |
| 146 | { | ||
| 147 | 10449 | j = Link[i].node1; | |
| 148 | 10449 | z = Node[j].invertElev + Link[i].offset1 + Link[i].xsect.yFull; | |
| 149 |
2/2✓ Branch 0 taken 4402 times.
✓ Branch 1 taken 6047 times.
|
10449 | Node[j].crownElev = MAX(Node[j].crownElev, z); |
| 150 | |||
| 151 | 10449 | j = Link[i].node2; | |
| 152 | 10449 | z = Node[j].invertElev + Link[i].offset2 + Link[i].xsect.yFull; | |
| 153 |
2/2✓ Branch 0 taken 3396 times.
✓ Branch 1 taken 7053 times.
|
10449 | Node[j].crownElev = MAX(Node[j].crownElev, z); |
| 154 | 10449 | Link[i].flowClass = DRY; | |
| 155 | 10449 | Link[i].dqdh = 0.0; | |
| 156 | } | ||
| 157 | |||
| 158 | // --- set crown cutoff for finding top width of closed conduits | ||
| 159 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 41 times.
|
41 | if ( SurchargeMethod == SLOT ) CrownCutoff = SLOT_CROWN_CUTOFF; |
| 160 | 41 | else CrownCutoff = EXTRAN_CROWN_CUTOFF; | |
| 161 | } | ||
| 162 | |||
| 163 | //============================================================================= | ||
| 164 | |||
| 165 | 41 | void dynwave_close() | |
| 166 | // | ||
| 167 | // Input: none | ||
| 168 | // Output: none | ||
| 169 | // Purpose: frees memory allocated for dynamic wave routing method. | ||
| 170 | // | ||
| 171 | { | ||
| 172 |
1/2✓ Branch 0 taken 41 times.
✗ Branch 1 not taken.
|
41 | FREE(Xnode); |
| 173 | 41 | } | |
| 174 | |||
| 175 | //============================================================================= | ||
| 176 | |||
| 177 | 41 | void dynwave_validate() | |
| 178 | // | ||
| 179 | // Input: none | ||
| 180 | // Output: none | ||
| 181 | // Purpose: adjusts dynamic wave routing options. | ||
| 182 | // | ||
| 183 | { | ||
| 184 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 41 times.
|
41 | if ( MinRouteStep > RouteStep ) MinRouteStep = RouteStep; |
| 185 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 40 times.
|
41 | if ( MinRouteStep < MINTIMESTEP ) MinRouteStep = MINTIMESTEP; |
| 186 |
2/2✓ Branch 0 taken 9 times.
✓ Branch 1 taken 32 times.
|
41 | if ( MinSurfArea == 0.0 ) MinSurfArea = DEFAULT_SURFAREA; |
| 187 | 32 | else MinSurfArea /= UCF(LENGTH) * UCF(LENGTH); | |
| 188 |
2/2✓ Branch 0 taken 9 times.
✓ Branch 1 taken 32 times.
|
41 | if ( HeadTol == 0.0 ) HeadTol = DEFAULT_HEADTOL; |
| 189 | 32 | else HeadTol /= UCF(LENGTH); | |
| 190 |
2/2✓ Branch 0 taken 9 times.
✓ Branch 1 taken 32 times.
|
41 | if ( MaxTrials == 0 ) MaxTrials = DEFAULT_MAXTRIALS; |
| 191 | 41 | } | |
| 192 | |||
| 193 | //============================================================================= | ||
| 194 | |||
| 195 | 367748 | double dynwave_getRoutingStep(double fixedStep) | |
| 196 | // | ||
| 197 | // Input: fixedStep = user-supplied fixed time step (sec) | ||
| 198 | // Output: returns routing time step (sec) | ||
| 199 | // Purpose: computes variable routing time step if applicable. | ||
| 200 | // | ||
| 201 | { | ||
| 202 | // --- use user-supplied fixed step if variable step option turned off | ||
| 203 | // or if its smaller than the min. allowable variable time step | ||
| 204 |
2/2✓ Branch 0 taken 178572 times.
✓ Branch 1 taken 189176 times.
|
367748 | if ( CourantFactor == 0.0 ) return fixedStep; |
| 205 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 189176 times.
|
189176 | if ( fixedStep < MINTIMESTEP ) return fixedStep; |
| 206 | |||
| 207 | // --- at start of simulation (when current variable step is zero) | ||
| 208 | // use the minimum allowable time step | ||
| 209 |
2/2✓ Branch 0 taken 24 times.
✓ Branch 1 taken 189152 times.
|
189176 | if ( VariableStep == 0.0 ) |
| 210 | { | ||
| 211 | 24 | VariableStep = MinRouteStep; | |
| 212 | } | ||
| 213 | |||
| 214 | // --- otherwise compute variable step based on current flow solution | ||
| 215 | 189152 | else VariableStep = getVariableStep(fixedStep); | |
| 216 | |||
| 217 | // --- adjust step to be a multiple of a millisecond | ||
| 218 | 189176 | VariableStep = floor(1000.0 * VariableStep) / 1000.0; | |
| 219 | 189176 | return VariableStep; | |
| 220 | } | ||
| 221 | |||
| 222 | //============================================================================= | ||
| 223 | |||
| 224 | 367734 | int dynwave_execute(double tStep) | |
| 225 | // | ||
| 226 | // Input: links = array of topo sorted links indexes | ||
| 227 | // tStep = time step (sec) | ||
| 228 | // Output: returns number of iterations used | ||
| 229 | // Purpose: routes flows through drainage network over current time step. | ||
| 230 | // | ||
| 231 | { | ||
| 232 | int converged; | ||
| 233 | |||
| 234 | // --- initialize | ||
| 235 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 367734 times.
|
367734 | if ( ErrorCode ) return 0; |
| 236 | 367734 | Steps = 0; | |
| 237 | 367734 | converged = FALSE; | |
| 238 | 367734 | Omega = OMEGA; | |
| 239 | 367734 | initRoutingStep(); | |
| 240 | |||
| 241 | // --- keep iterating until convergence | ||
| 242 |
2/2✓ Branch 0 taken 860935 times.
✓ Branch 1 taken 16762 times.
|
877697 | while ( Steps < MaxTrials ) |
| 243 | { | ||
| 244 | // --- execute a routing step & check for nodal convergence | ||
| 245 | 860935 | initNodeStates(); | |
| 246 | 860935 | findLinkFlows(tStep); | |
| 247 | 860935 | converged = findNodeDepths(tStep); | |
| 248 | 860935 | Steps++; | |
| 249 |
2/2✓ Branch 0 taken 493201 times.
✓ Branch 1 taken 367734 times.
|
860935 | if ( Steps > 1 ) |
| 250 | { | ||
| 251 |
2/2✓ Branch 0 taken 350972 times.
✓ Branch 1 taken 142229 times.
|
493201 | if ( converged ) break; |
| 252 | |||
| 253 | // --- check if link calculations can be skipped in next step | ||
| 254 | 142229 | findBypassedLinks(); | |
| 255 | } | ||
| 256 | } | ||
| 257 |
2/2✓ Branch 0 taken 16762 times.
✓ Branch 1 taken 350972 times.
|
367734 | if ( !converged ) updateConvergenceStats(); |
| 258 | |||
| 259 | // --- identify any capacity-limited conduits | ||
| 260 | 367734 | findLimitedLinks(); | |
| 261 | 367734 | return Steps; | |
| 262 | } | ||
| 263 | |||
| 264 | //============================================================================= | ||
| 265 | |||
| 266 | 16762 | void updateConvergenceStats() | |
| 267 | { | ||
| 268 | int i; | ||
| 269 | 16762 | NonConvergeCount++; | |
| 270 |
2/2✓ Branch 0 taken 1824857 times.
✓ Branch 1 taken 16762 times.
|
1841619 | for (i = 0; i < Nobjects[NODE]; i++) |
| 271 | 1824857 | stats_updateConvergenceStats(i, Xnode[i].converged); | |
| 272 | 16762 | } | |
| 273 | |||
| 274 | //============================================================================= | ||
| 275 | |||
| 276 | 367734 | void initRoutingStep() | |
| 277 | { | ||
| 278 | int i; | ||
| 279 |
2/2✓ Branch 0 taken 25222729 times.
✓ Branch 1 taken 367734 times.
|
25590463 | for (i = 0; i < Nobjects[NODE]; i++) |
| 280 | { | ||
| 281 | 25222729 | Xnode[i].converged = FALSE; | |
| 282 | 25222729 | Xnode[i].dYdT = 0.0; | |
| 283 | } | ||
| 284 |
2/2✓ Branch 0 taken 25213598 times.
✓ Branch 1 taken 367734 times.
|
25581332 | for (i = 0; i < Nobjects[LINK]; i++) |
| 285 | { | ||
| 286 | 25213598 | Link[i].bypassed = FALSE; | |
| 287 | 25213598 | Link[i].surfArea1 = 0.0; | |
| 288 | 25213598 | Link[i].surfArea2 = 0.0; | |
| 289 | } | ||
| 290 | |||
| 291 | // --- a2 preserves conduit area from solution at last time step | ||
| 292 |
2/2✓ Branch 0 taken 24799123 times.
✓ Branch 1 taken 367734 times.
|
25166857 | for ( i = 0; i < Nlinks[CONDUIT]; i++) Conduit[i].a2 = Conduit[i].a1; |
| 293 | 367734 | } | |
| 294 | |||
| 295 | //============================================================================= | ||
| 296 | |||
| 297 | 860935 | void initNodeStates() | |
| 298 | // | ||
| 299 | // Input: none | ||
| 300 | // Output: none | ||
| 301 | // Purpose: initializes node's surface area, inflow & outflow | ||
| 302 | // | ||
| 303 | { | ||
| 304 | int i; | ||
| 305 | |||
| 306 |
2/2✓ Branch 0 taken 67319581 times.
✓ Branch 1 taken 860935 times.
|
68180516 | for (i = 0; i < Nobjects[NODE]; i++) |
| 307 | { | ||
| 308 | // --- initialize nodal surface area | ||
| 309 |
2/2✓ Branch 0 taken 2800436 times.
✓ Branch 1 taken 64519145 times.
|
67319581 | if ( AllowPonding ) |
| 310 | { | ||
| 311 | 2800436 | Xnode[i].newSurfArea = node_getPondedArea(i, Node[i].newDepth); | |
| 312 | } | ||
| 313 | else | ||
| 314 | { | ||
| 315 | 64519145 | Xnode[i].newSurfArea = node_getSurfArea(i, Node[i].newDepth); | |
| 316 | } | ||
| 317 | |||
| 318 | // --- initialize nodal inflow & outflow | ||
| 319 | 67319581 | Node[i].inflow = 0.0; | |
| 320 | 67319581 | Node[i].outflow = Node[i].losses; | |
| 321 |
2/2✓ Branch 0 taken 67277371 times.
✓ Branch 1 taken 42210 times.
|
67319581 | if ( Node[i].newLatFlow >= 0.0 ) |
| 322 | { | ||
| 323 | 67277371 | Node[i].inflow += Node[i].newLatFlow; | |
| 324 | } | ||
| 325 | else | ||
| 326 | { | ||
| 327 | 42210 | Node[i].outflow -= Node[i].newLatFlow; | |
| 328 | } | ||
| 329 | 67319581 | Xnode[i].sumdqdh = 0.0; | |
| 330 | } | ||
| 331 | 860935 | } | |
| 332 | |||
| 333 | //============================================================================= | ||
| 334 | |||
| 335 | 142229 | void findBypassedLinks() | |
| 336 | { | ||
| 337 | int i; | ||
| 338 |
2/2✓ Branch 0 taken 18485769 times.
✓ Branch 1 taken 142229 times.
|
18627998 | for (i = 0; i < Nobjects[LINK]; i++) |
| 339 | { | ||
| 340 |
2/2✓ Branch 0 taken 18227499 times.
✓ Branch 1 taken 258270 times.
|
18485769 | if ( Xnode[Link[i].node1].converged && |
| 341 |
2/2✓ Branch 0 taken 17253336 times.
✓ Branch 1 taken 974163 times.
|
18227499 | Xnode[Link[i].node2].converged ) |
| 342 | 17253336 | Link[i].bypassed = TRUE; | |
| 343 | 1232433 | else Link[i].bypassed = FALSE; | |
| 344 | } | ||
| 345 | 142229 | } | |
| 346 | |||
| 347 | //============================================================================= | ||
| 348 | |||
| 349 | 367734 | void findLimitedLinks() | |
| 350 | // | ||
| 351 | // Input: none | ||
| 352 | // Output: none | ||
| 353 | // Purpose: determines if a conduit link is capacity limited. | ||
| 354 | // | ||
| 355 | { | ||
| 356 | int j, n1, n2, k; | ||
| 357 | double h1, h2; | ||
| 358 | |||
| 359 |
2/2✓ Branch 0 taken 25213598 times.
✓ Branch 1 taken 367734 times.
|
25581332 | for (j = 0; j < Nobjects[LINK]; j++) |
| 360 | { | ||
| 361 | // ---- check only non-dummy conduit links | ||
| 362 |
2/2✓ Branch 1 taken 414475 times.
✓ Branch 2 taken 24799123 times.
|
25213598 | if ( !isTrueConduit(j) ) continue; |
| 363 | |||
| 364 | // --- check that upstream end is full | ||
| 365 | 24799123 | k = Link[j].subIndex; | |
| 366 | 24799123 | Conduit[k].capacityLimited = FALSE; | |
| 367 |
2/2✓ Branch 0 taken 343671 times.
✓ Branch 1 taken 24455452 times.
|
24799123 | if ( Conduit[k].a1 >= Link[j].xsect.aFull ) |
| 368 | { | ||
| 369 | // --- check if HGL slope > conduit slope | ||
| 370 | 343671 | n1 = Link[j].node1; | |
| 371 | 343671 | n2 = Link[j].node2; | |
| 372 | 343671 | h1 = Node[n1].newDepth + Node[n1].invertElev; | |
| 373 | 343671 | h2 = Node[n2].newDepth + Node[n2].invertElev; | |
| 374 |
2/2✓ Branch 0 taken 74187 times.
✓ Branch 1 taken 269484 times.
|
343671 | if ( (h1 - h2) > fabs(Conduit[k].slope) * Conduit[k].length ) |
| 375 | 74187 | Conduit[k].capacityLimited = TRUE; | |
| 376 | } | ||
| 377 | } | ||
| 378 | 367734 | } | |
| 379 | |||
| 380 | //============================================================================= | ||
| 381 | |||
| 382 | 860935 | void findLinkFlows(double dt) | |
| 383 | { | ||
| 384 | int i; | ||
| 385 | |||
| 386 | // --- find new flow in each non-dummy conduit | ||
| 387 | 860935 | #pragma omp parallel num_threads(NumThreads) | |
| 388 | { | ||
| 389 | #pragma omp for | ||
| 390 | for ( i = 0; i < Nobjects[LINK]; i++) | ||
| 391 | { | ||
| 392 | if ( isTrueConduit(i) && !Link[i].bypassed ) | ||
| 393 | dwflow_findConduitFlow(i, Steps, Omega, dt); | ||
| 394 | } | ||
| 395 | } | ||
| 396 | |||
| 397 | // --- update inflow/outflows for nodes attached to non-dummy conduits | ||
| 398 |
2/2✓ Branch 0 taken 67121326 times.
✓ Branch 1 taken 860935 times.
|
67982261 | for ( i = 0; i < Nobjects[LINK]; i++) |
| 399 | { | ||
| 400 |
2/2✓ Branch 1 taken 65741888 times.
✓ Branch 2 taken 1379438 times.
|
67121326 | if ( isTrueConduit(i) ) updateNodeFlows(i); |
| 401 | } | ||
| 402 | |||
| 403 | // --- find new flows for all dummy conduits, pumps & regulators | ||
| 404 |
2/2✓ Branch 0 taken 67121326 times.
✓ Branch 1 taken 860935 times.
|
67982261 | for ( i = 0; i < Nobjects[LINK]; i++) |
| 405 | { | ||
| 406 |
2/2✓ Branch 1 taken 1379438 times.
✓ Branch 2 taken 65741888 times.
|
67121326 | if ( !isTrueConduit(i) ) |
| 407 | { | ||
| 408 |
2/2✓ Branch 0 taken 915642 times.
✓ Branch 1 taken 463796 times.
|
1379438 | if ( !Link[i].bypassed ) findNonConduitFlow(i, dt); |
| 409 | 1379438 | updateNodeFlows(i); | |
| 410 | } | ||
| 411 | } | ||
| 412 | 860935 | } | |
| 413 | |||
| 414 | //============================================================================= | ||
| 415 | |||
| 416 | 226577576 | int isTrueConduit(int j) | |
| 417 | { | ||
| 418 |
3/4✓ Branch 0 taken 222024787 times.
✓ Branch 1 taken 4552789 times.
✓ Branch 2 taken 222024787 times.
✗ Branch 3 not taken.
|
226577576 | return ( Link[j].type == CONDUIT && Link[j].xsect.type != DUMMY ); |
| 419 | } | ||
| 420 | |||
| 421 | //============================================================================= | ||
| 422 | |||
| 423 | 915642 | void findNonConduitFlow(int i, double dt) | |
| 424 | // | ||
| 425 | // Input: i = link index | ||
| 426 | // dt = time step (sec) | ||
| 427 | // Output: none | ||
| 428 | // Purpose: finds new flow in a non-conduit-type link | ||
| 429 | // | ||
| 430 | { | ||
| 431 | double qLast; // previous link flow (cfs) | ||
| 432 | double qNew; // new link flow (cfs) | ||
| 433 | |||
| 434 | // --- get link flow from last iteration | ||
| 435 | 915642 | qLast = Link[i].newFlow; | |
| 436 | 915642 | Link[i].dqdh = 0.0; | |
| 437 | |||
| 438 | // --- get new inflow to link from its upstream node | ||
| 439 | // (link_getInflow returns 0 if flap gate closed or pump is offline) | ||
| 440 | 915642 | qNew = link_getInflow(i); | |
| 441 |
2/2✓ Branch 0 taken 159330 times.
✓ Branch 1 taken 756312 times.
|
915642 | if ( Link[i].type == PUMP ) qNew = getModPumpFlow(i, qNew, dt); |
| 442 | |||
| 443 | // --- find surface area at each end of link | ||
| 444 | 915642 | findNonConduitSurfArea(i); | |
| 445 | |||
| 446 | // --- apply under-relaxation with flow from previous iteration; | ||
| 447 | // --- do not allow flow to change direction without first being 0 | ||
| 448 |
4/4✓ Branch 0 taken 501167 times.
✓ Branch 1 taken 414475 times.
✓ Branch 2 taken 421190 times.
✓ Branch 3 taken 79977 times.
|
915642 | if ( Steps > 0 && Link[i].type != PUMP ) |
| 449 | { | ||
| 450 | 421190 | qNew = (1.0 - Omega) * qLast + Omega * qNew; | |
| 451 |
4/4✓ Branch 0 taken 56930 times.
✓ Branch 1 taken 364260 times.
✓ Branch 2 taken 28011 times.
✓ Branch 3 taken 28919 times.
|
421190 | if ( qNew * qLast < 0.0 ) qNew = 0.001 * SGN(qNew); |
| 452 | } | ||
| 453 | 915642 | Link[i].newFlow = qNew; | |
| 454 | 915642 | } | |
| 455 | |||
| 456 | //============================================================================= | ||
| 457 | |||
| 458 | 159330 | double getModPumpFlow(int i, double q, double dt) | |
| 459 | // | ||
| 460 | // Input: i = link index | ||
| 461 | // q = pump flow from pump curve (cfs) | ||
| 462 | // dt = time step (sec) | ||
| 463 | // Output: returns modified pump flow rate (cfs) | ||
| 464 | // Purpose: modifies pump curve pumping rate depending on amount of water | ||
| 465 | // available at pump's inlet node. | ||
| 466 | // | ||
| 467 | { | ||
| 468 | 159330 | int j = Link[i].node1; // pump's inlet node index | |
| 469 | 159330 | int k = Link[i].subIndex; // pump's index | |
| 470 | double newNetInflow; // inflow - outflow rate (cfs) | ||
| 471 | double netFlowVolume; // inflow - outflow volume (ft3) | ||
| 472 | double y; // node depth (ft) | ||
| 473 | |||
| 474 |
2/2✓ Branch 0 taken 28948 times.
✓ Branch 1 taken 130382 times.
|
159330 | if ( q == 0.0 ) return q; |
| 475 | |||
| 476 | // --- case where inlet node is a storage node: | ||
| 477 | // prevent node volume from going negative | ||
| 478 |
2/2✓ Branch 0 taken 25835 times.
✓ Branch 1 taken 104547 times.
|
130382 | if ( Node[j].type == STORAGE ) return node_getMaxOutflow(j, q, dt); |
| 479 | |||
| 480 | // --- case where inlet is a non-storage node | ||
| 481 |
1/3✗ Branch 0 not taken.
✗ Branch 1 not taken.
✓ Branch 2 taken 104547 times.
|
104547 | switch ( Pump[k].type ) |
| 482 | { | ||
| 483 | // --- for Type1 pump, a volume is computed for inlet node, | ||
| 484 | // so make sure it doesn't go negative | ||
| 485 | ✗ | case TYPE1_PUMP: | |
| 486 | ✗ | return node_getMaxOutflow(j, q, dt); | |
| 487 | |||
| 488 | // --- for other types of pumps, if pumping rate would make depth | ||
| 489 | // at upstream node negative, then set pumping rate = inflow | ||
| 490 | ✗ | case TYPE2_PUMP: | |
| 491 | case TYPE4_PUMP: | ||
| 492 | case TYPE3_PUMP: | ||
| 493 | ✗ | newNetInflow = Node[j].inflow - Node[j].outflow - q; | |
| 494 | ✗ | netFlowVolume = 0.5 * (Node[j].oldNetInflow + newNetInflow ) * dt; | |
| 495 | ✗ | y = Node[j].oldDepth + netFlowVolume / Xnode[j].newSurfArea; | |
| 496 | ✗ | if ( y <= 0.0 ) return Node[j].inflow; | |
| 497 | } | ||
| 498 | 104547 | return q; | |
| 499 | } | ||
| 500 | |||
| 501 | //============================================================================= | ||
| 502 | |||
| 503 | 915642 | void findNonConduitSurfArea(int i) | |
| 504 | // | ||
| 505 | // Input: i = link index | ||
| 506 | // Output: none | ||
| 507 | // Purpose: finds the surface area contributed by a non-conduit | ||
| 508 | // link to its upstream and downstream nodes. | ||
| 509 | // | ||
| 510 | { | ||
| 511 |
2/2✓ Branch 0 taken 156514 times.
✓ Branch 1 taken 759128 times.
|
915642 | if ( Link[i].type == ORIFICE ) |
| 512 | { | ||
| 513 | 156514 | Link[i].surfArea1 = Orifice[Link[i].subIndex].surfArea / 2.; | |
| 514 | } | ||
| 515 | |||
| 516 | // --- no surface area for weirs to maintain SWMM 4 compatibility | ||
| 517 | 759128 | else Link[i].surfArea1 = 0.0; | |
| 518 | |||
| 519 | 915642 | Link[i].surfArea2 = Link[i].surfArea1; | |
| 520 |
2/2✓ Branch 0 taken 907873 times.
✓ Branch 1 taken 7769 times.
|
915642 | if ( Link[i].flowClass == UP_CRITICAL || |
| 521 |
2/2✓ Branch 0 taken 280121 times.
✓ Branch 1 taken 627752 times.
|
915642 | Node[Link[i].node1].type == STORAGE ) Link[i].surfArea1 = 0.0; |
| 522 |
2/2✓ Branch 0 taken 763871 times.
✓ Branch 1 taken 151771 times.
|
915642 | if ( Link[i].flowClass == DN_CRITICAL || |
| 523 |
2/2✓ Branch 0 taken 25 times.
✓ Branch 1 taken 763846 times.
|
915642 | Node[Link[i].node2].type == STORAGE ) Link[i].surfArea2 = 0.0; |
| 524 | 915642 | } | |
| 525 | |||
| 526 | //============================================================================= | ||
| 527 | |||
| 528 | 67121326 | void updateNodeFlows(int i) | |
| 529 | // | ||
| 530 | // Input: i = link index | ||
| 531 | // q = link flow rate (cfs) | ||
| 532 | // Output: none | ||
| 533 | // Purpose: updates cumulative inflow & outflow at link's end nodes. | ||
| 534 | // | ||
| 535 | { | ||
| 536 | int k; | ||
| 537 | 67121326 | int barrels = 1; | |
| 538 | 67121326 | int n1 = Link[i].node1; | |
| 539 | 67121326 | int n2 = Link[i].node2; | |
| 540 | 67121326 | double q = Link[i].newFlow; | |
| 541 | 67121326 | double conduitLossRate = 0.0; | |
| 542 | |||
| 543 | // --- update total inflow & outflow at upstream/downstream nodes | ||
| 544 |
2/2✓ Branch 0 taken 64137490 times.
✓ Branch 1 taken 2983836 times.
|
67121326 | if ( q >= 0.0 ) |
| 545 | { | ||
| 546 | 64137490 | Node[n1].outflow += q; | |
| 547 | 64137490 | Node[n2].inflow += q; | |
| 548 | } | ||
| 549 | else | ||
| 550 | { | ||
| 551 | 2983836 | Node[n1].inflow -= q; | |
| 552 | 2983836 | Node[n2].outflow -= q; | |
| 553 | } | ||
| 554 | |||
| 555 | // --- add any uniform evap & seepage loss from conduit link | ||
| 556 |
2/2✓ Branch 0 taken 65741888 times.
✓ Branch 1 taken 1379438 times.
|
67121326 | if ( Link[i].type == CONDUIT ) |
| 557 | { | ||
| 558 | 65741888 | k = Link[i].subIndex; | |
| 559 | 65741888 | barrels = Conduit[k].barrels; | |
| 560 | 65741888 | conduitLossRate = (Conduit[k].evapLossRate + Conduit[k].seepLossRate) * | |
| 561 | barrels; | ||
| 562 |
2/2✓ Branch 0 taken 1922368 times.
✓ Branch 1 taken 63819520 times.
|
65741888 | if (conduitLossRate > 0.0) |
| 563 | { | ||
| 564 | // --- outfall nodes do not share evap & seepage losses | ||
| 565 |
3/4✓ Branch 0 taken 1922368 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 935582 times.
✓ Branch 3 taken 986786 times.
|
1922368 | if (Node[n1].type != OUTFALL && Node[n2].type != OUTFALL) |
| 566 | 935582 | conduitLossRate /= 2.0; | |
| 567 |
1/2✓ Branch 0 taken 1922368 times.
✗ Branch 1 not taken.
|
1922368 | if (Node[n1].type != OUTFALL) |
| 568 | 1922368 | Node[n1].outflow += conduitLossRate; | |
| 569 |
2/2✓ Branch 0 taken 935582 times.
✓ Branch 1 taken 986786 times.
|
1922368 | if (Node[n2].type != OUTFALL) |
| 570 | 935582 | Node[n2].outflow += conduitLossRate; | |
| 571 | } | ||
| 572 | } | ||
| 573 | |||
| 574 | // --- add surf. area contributions to upstream/downstream nodes | ||
| 575 | 67121326 | Xnode[Link[i].node1].newSurfArea += Link[i].surfArea1 * barrels; | |
| 576 | 67121326 | Xnode[Link[i].node2].newSurfArea += Link[i].surfArea2 * barrels; | |
| 577 | |||
| 578 | // --- update summed value of dqdh at each end node | ||
| 579 | 67121326 | Xnode[Link[i].node1].sumdqdh += Link[i].dqdh; | |
| 580 |
2/2✓ Branch 0 taken 305636 times.
✓ Branch 1 taken 66815690 times.
|
67121326 | if ( Link[i].type == PUMP ) |
| 581 | { | ||
| 582 | 305636 | k = Link[i].subIndex; | |
| 583 |
2/2✓ Branch 0 taken 293260 times.
✓ Branch 1 taken 12376 times.
|
305636 | if ( Pump[k].type != TYPE4_PUMP ) |
| 584 | { | ||
| 585 | 293260 | Xnode[n2].sumdqdh += Link[i].dqdh; | |
| 586 | } | ||
| 587 | } | ||
| 588 | 66815690 | else Xnode[n2].sumdqdh += Link[i].dqdh; | |
| 589 | 67121326 | } | |
| 590 | |||
| 591 | //============================================================================= | ||
| 592 | |||
| 593 | 860935 | int findNodeDepths(double dt) | |
| 594 | // | ||
| 595 | // Input: dt = time step (sec) | ||
| 596 | // Output: returns TRUE if depth change at all non-Outfall nodes is | ||
| 597 | // within the convergence tolerance and FALSE otherwise | ||
| 598 | // Purpose: finds new depth at all nodes and checks if convergence achieved. | ||
| 599 | // | ||
| 600 | { | ||
| 601 | int i; | ||
| 602 | 860935 | double yOld = 0.0; // previous node depth (ft) | |
| 603 | |||
| 604 | // --- compute outfall depths based on flow in connecting link | ||
| 605 |
2/2✓ Branch 1 taken 67121326 times.
✓ Branch 2 taken 860935 times.
|
67982261 | for ( i = 0; i < Nobjects[LINK]; i++ ) link_setOutfallDepth(i); |
| 606 | |||
| 607 | // --- compute new depth for all non-outfall nodes and determine if | ||
| 608 | // depth change from previous iteration is below tolerance | ||
| 609 | 860935 | #pragma omp parallel num_threads(NumThreads) | |
| 610 | { | ||
| 611 | #pragma omp for private(yOld) | ||
| 612 | for ( i = 0; i < Nobjects[NODE]; i++ ) | ||
| 613 | { | ||
| 614 | if ( Node[i].type == OUTFALL ) continue; | ||
| 615 | yOld = Node[i].newDepth; | ||
| 616 | setNodeDepth(i, dt); | ||
| 617 | Xnode[i].converged = TRUE; | ||
| 618 | if ( fabs(yOld - Node[i].newDepth) > HeadTol ) | ||
| 619 | { | ||
| 620 | Xnode[i].converged = FALSE; | ||
| 621 | } | ||
| 622 | } | ||
| 623 | } | ||
| 624 | |||
| 625 | // --- return FALSE if any non-Outfall node failed to converge | ||
| 626 |
2/2✓ Branch 0 taken 41854105 times.
✓ Branch 1 taken 663646 times.
|
42517751 | for (i = 0; i < Nobjects[NODE]; i++) |
| 627 | { | ||
| 628 |
2/2✓ Branch 0 taken 971683 times.
✓ Branch 1 taken 40882422 times.
|
41854105 | if ( Node[i].type == OUTFALL ) continue; |
| 629 |
2/2✓ Branch 0 taken 197289 times.
✓ Branch 1 taken 40685133 times.
|
40882422 | if (Xnode[i].converged == FALSE) return FALSE; |
| 630 | } | ||
| 631 | 663646 | return TRUE; | |
| 632 | } | ||
| 633 | |||
| 634 | //============================================================================= | ||
| 635 | |||
| 636 | 65337949 | void setNodeDepth(int i, double dt) | |
| 637 | // | ||
| 638 | // Input: i = node index | ||
| 639 | // dt = time step (sec) | ||
| 640 | // Output: none | ||
| 641 | // Purpose: sets depth at non-outfall node after current time step. | ||
| 642 | // | ||
| 643 | { | ||
| 644 | int canPond; // TRUE if node can pond overflows | ||
| 645 | int isPonded; // TRUE if node is currently ponded | ||
| 646 | 65337949 | int isSurcharged = FALSE; // TRUE if node is surcharged | |
| 647 | double dQ; // inflow minus outflow at node (cfs) | ||
| 648 | double dV; // change in node volume (ft3) | ||
| 649 | double dy; // change in node depth (ft) | ||
| 650 | double yMax; // max. depth at node (ft) | ||
| 651 | double yOld; // node depth at previous time step (ft) | ||
| 652 | double yLast; // previous node depth (ft) | ||
| 653 | double yNew; // new node depth (ft) | ||
| 654 | double yCrown; // depth to node crown (ft) | ||
| 655 | double surfArea; // node surface area (ft2) | ||
| 656 | double denom; // denominator term | ||
| 657 | double corr; // correction factor | ||
| 658 | double f; // relative surcharge depth | ||
| 659 | |||
| 660 | // --- see if node can pond water above it | ||
| 661 |
4/4✓ Branch 0 taken 2398898 times.
✓ Branch 1 taken 62939051 times.
✓ Branch 2 taken 2027959 times.
✓ Branch 3 taken 370939 times.
|
65337949 | canPond = (AllowPonding && Node[i].pondedArea > 0.0); |
| 662 |
4/4✓ Branch 0 taken 2027959 times.
✓ Branch 1 taken 63309990 times.
✓ Branch 2 taken 177183 times.
✓ Branch 3 taken 1850776 times.
|
65337949 | isPonded = (canPond && Node[i].newDepth > Node[i].fullDepth); |
| 663 | |||
| 664 | // --- initialize values | ||
| 665 | 65337949 | yCrown = Node[i].crownElev - Node[i].invertElev; | |
| 666 | 65337949 | yOld = Node[i].oldDepth; | |
| 667 | 65337949 | yLast = Node[i].newDepth; | |
| 668 | 65337949 | Node[i].overflow = 0.0; | |
| 669 | 65337949 | surfArea = Xnode[i].newSurfArea; | |
| 670 |
2/2✓ Branch 0 taken 62107288 times.
✓ Branch 1 taken 3230661 times.
|
65337949 | surfArea = MAX(surfArea, MinSurfArea); |
| 671 | |||
| 672 | // --- determine average net flow volume into node over the time step | ||
| 673 | 65337949 | dQ = Node[i].inflow - Node[i].outflow; | |
| 674 | 65337949 | dV = 0.5 * (Node[i].oldNetInflow + dQ) * dt; | |
| 675 | |||
| 676 | // --- determine if node is EXTRAN surcharged | ||
| 677 |
1/2✓ Branch 0 taken 65337949 times.
✗ Branch 1 not taken.
|
65337949 | if (SurchargeMethod == EXTRAN) |
| 678 | { | ||
| 679 | // --- ponded nodes don't surcharge | ||
| 680 |
2/2✓ Branch 0 taken 177183 times.
✓ Branch 1 taken 65160766 times.
|
65337949 | if (isPonded) isSurcharged = FALSE; |
| 681 | |||
| 682 | // --- closed storage units that are full are in surcharge | ||
| 683 |
2/2✓ Branch 0 taken 494647 times.
✓ Branch 1 taken 64666119 times.
|
65160766 | else if (Node[i].type == STORAGE) |
| 684 | { | ||
| 685 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 494647 times.
|
494647 | isSurcharged = (Node[i].surDepth > 0.0 && |
| 686 | ✗ | yLast > Node[i].fullDepth); | |
| 687 | } | ||
| 688 | |||
| 689 | // --- surcharge occurs when node depth exceeds top of its highest link | ||
| 690 |
4/4✓ Branch 0 taken 64661995 times.
✓ Branch 1 taken 4124 times.
✓ Branch 2 taken 1323610 times.
✓ Branch 3 taken 63338385 times.
|
64666119 | else isSurcharged = (yCrown > 0.0 && yLast > yCrown); |
| 691 | } | ||
| 692 | |||
| 693 | // --- if node not surcharged, base depth change on surface area | ||
| 694 |
2/2✓ Branch 0 taken 64014339 times.
✓ Branch 1 taken 1323610 times.
|
65337949 | if (!isSurcharged) |
| 695 | { | ||
| 696 | 64014339 | dy = dV / surfArea; | |
| 697 | 64014339 | yNew = yOld + dy; | |
| 698 | |||
| 699 | // --- save non-ponded surface area for use in surcharge algorithm | ||
| 700 |
2/2✓ Branch 0 taken 63837156 times.
✓ Branch 1 taken 177183 times.
|
64014339 | if ( !isPonded ) Xnode[i].oldSurfArea = surfArea; |
| 701 | |||
| 702 | // --- apply under-relaxation to new depth estimate | ||
| 703 |
2/2✓ Branch 0 taken 39912609 times.
✓ Branch 1 taken 24101730 times.
|
64014339 | if ( Steps > 0 ) |
| 704 | { | ||
| 705 | 39912609 | yNew = (1.0 - Omega) * yLast + Omega * yNew; | |
| 706 | } | ||
| 707 | |||
| 708 | // --- don't allow a ponded node to drop much below full depth | ||
| 709 |
4/4✓ Branch 0 taken 177183 times.
✓ Branch 1 taken 63837156 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 177178 times.
|
64014339 | if ( isPonded && yNew < Node[i].fullDepth ) |
| 710 | 5 | yNew = Node[i].fullDepth - FUDGE; | |
| 711 | } | ||
| 712 | |||
| 713 | // --- if node surcharged, base depth change on dqdh | ||
| 714 | // NOTE: depth change is w.r.t depth from previous | ||
| 715 | // iteration; also, do not apply under-relaxation. | ||
| 716 | else | ||
| 717 | { | ||
| 718 | // --- apply correction factor for upstream terminal nodes | ||
| 719 | 1323610 | corr = 1.0; | |
| 720 |
2/2✓ Branch 0 taken 107014 times.
✓ Branch 1 taken 1216596 times.
|
1323610 | if ( Node[i].degree < 0 ) corr = 0.6; |
| 721 | |||
| 722 | // --- allow surface area from last non-surcharged condition | ||
| 723 | // to influence dqdh if depth close to crown depth | ||
| 724 | 1323610 | denom = Xnode[i].sumdqdh; | |
| 725 |
2/2✓ Branch 0 taken 155646 times.
✓ Branch 1 taken 1167964 times.
|
1323610 | if ( yLast < 1.25 * yCrown ) |
| 726 | { | ||
| 727 | 155646 | f = (yLast - yCrown) / yCrown; | |
| 728 | 155646 | denom += (Xnode[i].oldSurfArea/dt - | |
| 729 | 155646 | Xnode[i].sumdqdh) * exp(-15.0 * f); | |
| 730 | } | ||
| 731 | |||
| 732 | // --- compute new estimate of node depth | ||
| 733 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1323610 times.
|
1323610 | if ( denom == 0.0 ) dy = 0.0; |
| 734 | 1323610 | else dy = corr * dQ / denom; | |
| 735 | 1323610 | yNew = yLast + dy; | |
| 736 |
2/2✓ Branch 0 taken 97 times.
✓ Branch 1 taken 1323513 times.
|
1323610 | if ( yNew < yCrown ) yNew = yCrown - FUDGE; |
| 737 | |||
| 738 | // --- don't allow a newly ponded node to rise much above full depth | ||
| 739 |
1/4✗ Branch 0 not taken.
✓ Branch 1 taken 1323610 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
1323610 | if ( canPond && yNew > Node[i].fullDepth ) |
| 740 | ✗ | yNew = Node[i].fullDepth + FUDGE; | |
| 741 | } | ||
| 742 | |||
| 743 | // --- depth cannot be negative | ||
| 744 |
2/2✓ Branch 0 taken 928 times.
✓ Branch 1 taken 65337021 times.
|
65337949 | if ( yNew < 0 ) yNew = 0.0; |
| 745 | |||
| 746 | // --- determine max. non-flooded depth | ||
| 747 | 65337949 | yMax = Node[i].fullDepth; | |
| 748 |
2/2✓ Branch 0 taken 63309990 times.
✓ Branch 1 taken 2027959 times.
|
65337949 | if ( canPond == FALSE ) yMax += Node[i].surDepth; |
| 749 | |||
| 750 | // --- find flooded depth & volume | ||
| 751 |
2/2✓ Branch 0 taken 252011 times.
✓ Branch 1 taken 65085938 times.
|
65337949 | if ( yNew > yMax ) |
| 752 | { | ||
| 753 | 252011 | yNew = getFloodedDepth(i, canPond, dV, yNew, yMax, dt); | |
| 754 | } | ||
| 755 | 65085938 | else Node[i].newVolume = node_getVolume(i, yNew); | |
| 756 | |||
| 757 | // --- compute change in depth w.r.t. time | ||
| 758 | 65337949 | Xnode[i].dYdT = fabs(yNew - yOld) / dt; | |
| 759 | |||
| 760 | // --- save new depth for node | ||
| 761 | 65337949 | Node[i].newDepth = yNew; | |
| 762 | 65337949 | } | |
| 763 | |||
| 764 | //============================================================================= | ||
| 765 | |||
| 766 | 252011 | double getFloodedDepth(int i, int canPond, double dV, double yNew, | |
| 767 | double yMax, double dt) | ||
| 768 | // | ||
| 769 | // Input: i = node index | ||
| 770 | // canPond = TRUE if water can pond over node | ||
| 771 | // isPonded = TRUE if water is currently ponded | ||
| 772 | // dV = change in volume over time step (ft3) | ||
| 773 | // yNew = current depth at node (ft) | ||
| 774 | // yMax = max. depth at node before ponding (ft) | ||
| 775 | // dt = time step (sec) | ||
| 776 | // Output: returns depth at node when flooded (ft) | ||
| 777 | // Purpose: computes depth, volume and overflow for a flooded node. | ||
| 778 | // | ||
| 779 | { | ||
| 780 |
2/2✓ Branch 0 taken 74828 times.
✓ Branch 1 taken 177183 times.
|
252011 | if ( canPond == FALSE ) |
| 781 | { | ||
| 782 | 74828 | Node[i].overflow = dV / dt; | |
| 783 | 74828 | Node[i].newVolume = Node[i].fullVolume; | |
| 784 | 74828 | yNew = yMax; | |
| 785 | } | ||
| 786 | else | ||
| 787 | { | ||
| 788 |
2/2✓ Branch 0 taken 176547 times.
✓ Branch 1 taken 636 times.
|
177183 | Node[i].newVolume = MAX((Node[i].oldVolume+dV), Node[i].fullVolume); |
| 789 | 354366 | Node[i].overflow = (Node[i].newVolume - | |
| 790 |
1/2✓ Branch 0 taken 177183 times.
✗ Branch 1 not taken.
|
177183 | MAX(Node[i].oldVolume, Node[i].fullVolume)) / dt; |
| 791 | } | ||
| 792 |
2/2✓ Branch 0 taken 123145 times.
✓ Branch 1 taken 128866 times.
|
252011 | if ( Node[i].overflow < FUDGE ) Node[i].overflow = 0.0; |
| 793 | 252011 | return yNew; | |
| 794 | |||
| 795 | } | ||
| 796 | |||
| 797 | //============================================================================= | ||
| 798 | |||
| 799 | 189152 | double getVariableStep(double maxStep) | |
| 800 | // | ||
| 801 | // Input: maxStep = user-supplied max. time step (sec) | ||
| 802 | // Output: returns time step (sec) | ||
| 803 | // Purpose: finds time step that satisfies stability criterion but | ||
| 804 | // is no greater than the user-supplied max. time step. | ||
| 805 | // | ||
| 806 | { | ||
| 807 | 189152 | int minLink = -1; // index of link w/ min. time step | |
| 808 | 189152 | int minNode = -1; // index of node w/ min. time step | |
| 809 | double tMin; // allowable time step (sec) | ||
| 810 | double tMinLink; // allowable time step for links (sec) | ||
| 811 | double tMinNode; // allowable time step for nodes (sec) | ||
| 812 | |||
| 813 | // --- find stable time step for links & then nodes | ||
| 814 | 189152 | tMin = maxStep; | |
| 815 | 189152 | tMinLink = getLinkStep(tMin, &minLink); | |
| 816 | 189152 | tMinNode = getNodeStep(tMinLink, &minNode); | |
| 817 | |||
| 818 | // --- use smaller of the link and node time step | ||
| 819 | 189152 | tMin = tMinLink; | |
| 820 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 189152 times.
|
189152 | if ( tMinNode < tMin ) |
| 821 | { | ||
| 822 | ✗ | tMin = tMinNode ; | |
| 823 | ✗ | minLink = -1; | |
| 824 | } | ||
| 825 | |||
| 826 | // --- update count of times the minimum node or link was critical | ||
| 827 | 189152 | stats_updateCriticalTimeCount(minNode, minLink); | |
| 828 | |||
| 829 | // --- don't let time step go below an absolute minimum | ||
| 830 |
2/2✓ Branch 0 taken 35 times.
✓ Branch 1 taken 189117 times.
|
189152 | if ( tMin < MinRouteStep ) tMin = MinRouteStep; |
| 831 | 189152 | return tMin; | |
| 832 | } | ||
| 833 | |||
| 834 | //============================================================================= | ||
| 835 | |||
| 836 | 189152 | double getLinkStep(double tMin, int *minLink) | |
| 837 | // | ||
| 838 | // Input: tMin = critical time step found so far (sec) | ||
| 839 | // Output: minLink = index of link with critical time step; | ||
| 840 | // returns critical time step (sec) | ||
| 841 | // Purpose: finds critical time step for conduits based on Courant criterion. | ||
| 842 | // | ||
| 843 | { | ||
| 844 | int i; // link index | ||
| 845 | int k; // conduit index | ||
| 846 | double q; // conduit flow (cfs) | ||
| 847 | double t; // time step (sec) | ||
| 848 | 189152 | double tLink = tMin; // critical link time step (sec) | |
| 849 | |||
| 850 | // --- examine each conduit link | ||
| 851 |
2/2✓ Branch 0 taken 22944364 times.
✓ Branch 1 taken 189152 times.
|
23133516 | for ( i = 0; i < Nobjects[LINK]; i++ ) |
| 852 | { | ||
| 853 |
2/2✓ Branch 0 taken 22534129 times.
✓ Branch 1 taken 410235 times.
|
22944364 | if ( Link[i].type == CONDUIT ) |
| 854 | { | ||
| 855 | // --- skip conduits with negligible flow, area or Fr | ||
| 856 | 22534129 | k = Link[i].subIndex; | |
| 857 | 22534129 | q = fabs(Link[i].newFlow) / Conduit[k].barrels; | |
| 858 |
2/2✓ Branch 0 taken 21821522 times.
✓ Branch 1 taken 712607 times.
|
22534129 | if ( q <= FUDGE |
| 859 |
1/2✓ Branch 0 taken 21821522 times.
✗ Branch 1 not taken.
|
21821522 | || Conduit[k].a1 <= FUDGE |
| 860 |
2/2✓ Branch 0 taken 647064 times.
✓ Branch 1 taken 21174458 times.
|
21821522 | || Link[i].froude <= 0.01 |
| 861 | 1359671 | ) continue; | |
| 862 | |||
| 863 | // --- compute time step to satisfy Courant condition | ||
| 864 | 21174458 | t = Link[i].newVolume / Conduit[k].barrels / q; | |
| 865 | 21174458 | t = t * Conduit[k].modLength / link_getLength(i); | |
| 866 | 21174458 | t = t * Link[i].froude / (1.0 + Link[i].froude) * CourantFactor; | |
| 867 | |||
| 868 | // --- update critical link time step | ||
| 869 |
2/2✓ Branch 0 taken 47349 times.
✓ Branch 1 taken 21127109 times.
|
21174458 | if ( t < tLink ) |
| 870 | { | ||
| 871 | 47349 | tLink = t; | |
| 872 | 47349 | *minLink = i; | |
| 873 | } | ||
| 874 | } | ||
| 875 | } | ||
| 876 | 189152 | return tLink; | |
| 877 | } | ||
| 878 | |||
| 879 | //============================================================================= | ||
| 880 | |||
| 881 | 189152 | double getNodeStep(double tMin, int *minNode) | |
| 882 | // | ||
| 883 | // Input: tMin = critical time step found so far (sec) | ||
| 884 | // Output: minNode = index of node with critical time step; | ||
| 885 | // returns critical time step (sec) | ||
| 886 | // Purpose: finds critical time step for nodes based on max. allowable | ||
| 887 | // projected change in depth. | ||
| 888 | // | ||
| 889 | { | ||
| 890 | int i; // node index | ||
| 891 | double maxDepth; // max. depth allowed at node (ft) | ||
| 892 | double dYdT; // change in depth per unit time (ft/sec) | ||
| 893 | double t1; // time needed to reach depth limit (sec) | ||
| 894 | 189152 | double tNode = tMin; // critical node time step (sec) | |
| 895 | |||
| 896 | // --- find smallest time so that estimated change in nodal depth | ||
| 897 | // does not exceed safety factor * maxdepth | ||
| 898 |
2/2✓ Branch 0 taken 22774893 times.
✓ Branch 1 taken 189152 times.
|
22964045 | for ( i = 0; i < Nobjects[NODE]; i++ ) |
| 899 | { | ||
| 900 | // --- see if node can be skipped | ||
| 901 |
2/2✓ Branch 0 taken 450242 times.
✓ Branch 1 taken 22324651 times.
|
22774893 | if ( Node[i].type == OUTFALL ) continue; |
| 902 |
2/2✓ Branch 0 taken 502552 times.
✓ Branch 1 taken 21822099 times.
|
22324651 | if ( Node[i].newDepth <= FUDGE) continue; |
| 903 | 22464693 | if ( Node[i].newDepth + FUDGE >= | |
| 904 |
2/2✓ Branch 0 taken 642594 times.
✓ Branch 1 taken 21179505 times.
|
21822099 | Node[i].crownElev - Node[i].invertElev ) continue; |
| 905 | |||
| 906 | // --- define max. allowable depth change using crown elevation | ||
| 907 | 21179505 | maxDepth = (Node[i].crownElev - Node[i].invertElev) * 0.25; | |
| 908 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 21179505 times.
|
21179505 | if ( maxDepth < FUDGE ) continue; |
| 909 | 21179505 | dYdT = Xnode[i].dYdT; | |
| 910 |
2/2✓ Branch 0 taken 19433100 times.
✓ Branch 1 taken 1746405 times.
|
21179505 | if (dYdT < FUDGE ) continue; |
| 911 | |||
| 912 | // --- compute time to reach max. depth & compare with critical time | ||
| 913 | 1746405 | t1 = maxDepth / dYdT; | |
| 914 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1746405 times.
|
1746405 | if ( t1 < tNode ) |
| 915 | { | ||
| 916 | ✗ | tNode = t1; | |
| 917 | ✗ | *minNode = i; | |
| 918 | } | ||
| 919 | } | ||
| 920 | 189152 | return tNode; | |
| 921 | } | ||
| 922 |