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 | 44 | 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 | 44 | VariableStep = 0.0; | |
| 128 | 44 | Xnode = (TXnode *) calloc(Nobjects[NODE], sizeof(TXnode)); | |
| 129 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 44 times.
|
44 | 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 334 times.
✓ Branch 1 taken 44 times.
|
378 | for (i = 0; i < Nobjects[NODE]; i++ ) |
| 138 | { | ||
| 139 | 334 | Xnode[i].newSurfArea = 0.0; | |
| 140 | 334 | Xnode[i].oldSurfArea = 0.0; | |
| 141 | 334 | Node[i].crownElev = Node[i].invertElev; | |
| 142 | } | ||
| 143 | |||
| 144 | // --- initialize links & update node crown elevations | ||
| 145 |
2/2✓ Branch 0 taken 273 times.
✓ Branch 1 taken 44 times.
|
317 | for (i = 0; i < Nobjects[LINK]; i++) |
| 146 | { | ||
| 147 | 273 | j = Link[i].node1; | |
| 148 | 273 | z = Node[j].invertElev + Link[i].offset1 + Link[i].xsect.yFull; | |
| 149 |
2/2✓ Branch 0 taken 71 times.
✓ Branch 1 taken 202 times.
|
273 | Node[j].crownElev = MAX(Node[j].crownElev, z); |
| 150 | |||
| 151 | 273 | j = Link[i].node2; | |
| 152 | 273 | z = Node[j].invertElev + Link[i].offset2 + Link[i].xsect.yFull; | |
| 153 |
2/2✓ Branch 0 taken 88 times.
✓ Branch 1 taken 185 times.
|
273 | Node[j].crownElev = MAX(Node[j].crownElev, z); |
| 154 | 273 | Link[i].flowClass = DRY; | |
| 155 | 273 | Link[i].dqdh = 0.0; | |
| 156 | } | ||
| 157 | |||
| 158 | // --- set crown cutoff for finding top width of closed conduits | ||
| 159 |
2/2✓ Branch 0 taken 1 time.
✓ Branch 1 taken 43 times.
|
44 | if ( SurchargeMethod == SLOT ) CrownCutoff = SLOT_CROWN_CUTOFF; |
| 160 | 43 | else CrownCutoff = EXTRAN_CROWN_CUTOFF; | |
| 161 | } | ||
| 162 | |||
| 163 | //============================================================================= | ||
| 164 | |||
| 165 | 44 | 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 44 times.
✗ Branch 1 not taken.
|
44 | FREE(Xnode); |
| 173 | 44 | } | |
| 174 | |||
| 175 | //============================================================================= | ||
| 176 | |||
| 177 | 44 | 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 44 times.
|
44 | if ( MinRouteStep > RouteStep ) MinRouteStep = RouteStep; |
| 185 |
2/2✓ Branch 0 taken 2 times.
✓ Branch 1 taken 42 times.
|
44 | if ( MinRouteStep < MINTIMESTEP ) MinRouteStep = MINTIMESTEP; |
| 186 |
2/2✓ Branch 0 taken 17 times.
✓ Branch 1 taken 27 times.
|
44 | if ( MinSurfArea == 0.0 ) MinSurfArea = DEFAULT_SURFAREA; |
| 187 | 27 | else MinSurfArea /= UCF(LENGTH) * UCF(LENGTH); | |
| 188 |
2/2✓ Branch 0 taken 18 times.
✓ Branch 1 taken 26 times.
|
44 | if ( HeadTol == 0.0 ) HeadTol = DEFAULT_HEADTOL; |
| 189 | 26 | else HeadTol /= UCF(LENGTH); | |
| 190 |
2/2✓ Branch 0 taken 18 times.
✓ Branch 1 taken 26 times.
|
44 | if ( MaxTrials == 0 ) MaxTrials = DEFAULT_MAXTRIALS; |
| 191 | 44 | } | |
| 192 | |||
| 193 | //============================================================================= | ||
| 194 | |||
| 195 | 431267 | 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 180732 times.
✓ Branch 1 taken 250535 times.
|
431267 | if ( CourantFactor == 0.0 ) return fixedStep; |
| 205 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 250535 times.
|
250535 | 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 31 times.
✓ Branch 1 taken 250504 times.
|
250535 | if ( VariableStep == 0.0 ) |
| 210 | { | ||
| 211 | 31 | VariableStep = MinRouteStep; | |
| 212 | } | ||
| 213 | |||
| 214 | // --- otherwise compute variable step based on current flow solution | ||
| 215 | 250504 | else VariableStep = getVariableStep(fixedStep); | |
| 216 | |||
| 217 | // --- adjust step to be a multiple of a millisecond | ||
| 218 | 250535 | VariableStep = floor(1000.0 * VariableStep) / 1000.0; | |
| 219 | 250535 | return VariableStep; | |
| 220 | } | ||
| 221 | |||
| 222 | //============================================================================= | ||
| 223 | |||
| 224 | 431255 | 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 431255 times.
|
431255 | if ( ErrorCode ) return 0; |
| 236 | 431255 | Steps = 0; | |
| 237 | 431255 | converged = FALSE; | |
| 238 | 431255 | Omega = OMEGA; | |
| 239 | 431255 | initRoutingStep(); | |
| 240 | |||
| 241 | // --- keep iterating until convergence | ||
| 242 |
2/2✓ Branch 0 taken 986324 times.
✓ Branch 1 taken 16738 times.
|
1003062 | while ( Steps < MaxTrials ) |
| 243 | { | ||
| 244 | // --- execute a routing step & check for nodal convergence | ||
| 245 | 986324 | initNodeStates(); | |
| 246 | 986324 | findLinkFlows(tStep); | |
| 247 | 986324 | converged = findNodeDepths(tStep); | |
| 248 | 986324 | Steps++; | |
| 249 |
2/2✓ Branch 0 taken 555069 times.
✓ Branch 1 taken 431255 times.
|
986324 | if ( Steps > 1 ) |
| 250 | { | ||
| 251 |
2/2✓ Branch 0 taken 414517 times.
✓ Branch 1 taken 140552 times.
|
555069 | if ( converged ) break; |
| 252 | |||
| 253 | // --- check if link calculations can be skipped in next step | ||
| 254 | 140552 | findBypassedLinks(); | |
| 255 | } | ||
| 256 | } | ||
| 257 |
2/2✓ Branch 0 taken 16738 times.
✓ Branch 1 taken 414517 times.
|
431255 | if ( !converged ) updateConvergenceStats(); |
| 258 | |||
| 259 | // --- identify any capacity-limited conduits | ||
| 260 | 431255 | findLimitedLinks(); | |
| 261 | 431255 | return Steps; | |
| 262 | } | ||
| 263 | |||
| 264 | //============================================================================= | ||
| 265 | |||
| 266 | 16738 | void updateConvergenceStats() | |
| 267 | { | ||
| 268 | int i; | ||
| 269 | 16738 | NonConvergeCount++; | |
| 270 |
2/2✓ Branch 0 taken 970889 times.
✓ Branch 1 taken 16738 times.
|
987627 | for (i = 0; i < Nobjects[NODE]; i++) |
| 271 | 970889 | stats_updateConvergenceStats(i, Xnode[i].converged); | |
| 272 | 16738 | } | |
| 273 | |||
| 274 | //============================================================================= | ||
| 275 | |||
| 276 | 431255 | void initRoutingStep() | |
| 277 | { | ||
| 278 | int i; | ||
| 279 |
2/2✓ Branch 0 taken 5136258 times.
✓ Branch 1 taken 431255 times.
|
5567513 | for (i = 0; i < Nobjects[NODE]; i++) |
| 280 | { | ||
| 281 | 5136258 | Xnode[i].converged = FALSE; | |
| 282 | 5136258 | Xnode[i].dYdT = 0.0; | |
| 283 | } | ||
| 284 |
2/2✓ Branch 0 taken 4749083 times.
✓ Branch 1 taken 431255 times.
|
5180338 | for (i = 0; i < Nobjects[LINK]; i++) |
| 285 | { | ||
| 286 | 4749083 | Link[i].bypassed = FALSE; | |
| 287 | 4749083 | Link[i].surfArea1 = 0.0; | |
| 288 | 4749083 | Link[i].surfArea2 = 0.0; | |
| 289 | } | ||
| 290 | |||
| 291 | // --- a2 preserves conduit area from solution at last time step | ||
| 292 |
2/2✓ Branch 0 taken 4379092 times.
✓ Branch 1 taken 431255 times.
|
4810347 | for ( i = 0; i < Nlinks[CONDUIT]; i++) Conduit[i].a2 = Conduit[i].a1; |
| 293 | 431255 | } | |
| 294 | |||
| 295 | //============================================================================= | ||
| 296 | |||
| 297 | 986324 | 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 17405312 times.
✓ Branch 1 taken 986324 times.
|
18391636 | for (i = 0; i < Nobjects[NODE]; i++) |
| 307 | { | ||
| 308 | // --- initialize nodal surface area | ||
| 309 |
2/2✓ Branch 0 taken 1619881 times.
✓ Branch 1 taken 15785431 times.
|
17405312 | if ( AllowPonding ) |
| 310 | { | ||
| 311 | 1619881 | Xnode[i].newSurfArea = node_getPondedArea(i, Node[i].newDepth); | |
| 312 | } | ||
| 313 | else | ||
| 314 | { | ||
| 315 | 15785431 | Xnode[i].newSurfArea = node_getSurfArea(i, Node[i].newDepth); | |
| 316 | } | ||
| 317 | |||
| 318 | // --- initialize nodal inflow & outflow | ||
| 319 | 17405312 | Node[i].inflow = 0.0; | |
| 320 | 17405312 | Node[i].outflow = Node[i].losses; | |
| 321 |
2/2✓ Branch 0 taken 17363102 times.
✓ Branch 1 taken 42210 times.
|
17405312 | if ( Node[i].newLatFlow >= 0.0 ) |
| 322 | { | ||
| 323 | 17363102 | Node[i].inflow += Node[i].newLatFlow; | |
| 324 | } | ||
| 325 | else | ||
| 326 | { | ||
| 327 | 42210 | Node[i].outflow -= Node[i].newLatFlow; | |
| 328 | } | ||
| 329 | 17405312 | Xnode[i].sumdqdh = 0.0; | |
| 330 | } | ||
| 331 | 986324 | } | |
| 332 | |||
| 333 | //============================================================================= | ||
| 334 | |||
| 335 | 140552 | void findBypassedLinks() | |
| 336 | { | ||
| 337 | int i; | ||
| 338 |
2/2✓ Branch 0 taken 7707334 times.
✓ Branch 1 taken 140552 times.
|
7847886 | for (i = 0; i < Nobjects[LINK]; i++) |
| 339 | { | ||
| 340 |
2/2✓ Branch 0 taken 7452706 times.
✓ Branch 1 taken 254628 times.
|
7707334 | if ( Xnode[Link[i].node1].converged && |
| 341 |
2/2✓ Branch 0 taken 6517890 times.
✓ Branch 1 taken 934816 times.
|
7452706 | Xnode[Link[i].node2].converged ) |
| 342 | 6517890 | Link[i].bypassed = TRUE; | |
| 343 | 1189444 | else Link[i].bypassed = FALSE; | |
| 344 | } | ||
| 345 | 140552 | } | |
| 346 | |||
| 347 | //============================================================================= | ||
| 348 | |||
| 349 | 431255 | 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 4749083 times.
✓ Branch 1 taken 431255 times.
|
5180338 | for (j = 0; j < Nobjects[LINK]; j++) |
| 360 | { | ||
| 361 | // ---- check only non-dummy conduit links | ||
| 362 |
2/2✓ Branch 1 taken 369991 times.
✓ Branch 2 taken 4379092 times.
|
4749083 | if ( !isTrueConduit(j) ) continue; |
| 363 | |||
| 364 | // --- check that upstream end is full | ||
| 365 | 4379092 | k = Link[j].subIndex; | |
| 366 | 4379092 | Conduit[k].capacityLimited = FALSE; | |
| 367 |
2/2✓ Branch 0 taken 196779 times.
✓ Branch 1 taken 4182313 times.
|
4379092 | if ( Conduit[k].a1 >= Link[j].xsect.aFull ) |
| 368 | { | ||
| 369 | // --- check if HGL slope > conduit slope | ||
| 370 | 196779 | n1 = Link[j].node1; | |
| 371 | 196779 | n2 = Link[j].node2; | |
| 372 | 196779 | h1 = Node[n1].newDepth + Node[n1].invertElev; | |
| 373 | 196779 | h2 = Node[n2].newDepth + Node[n2].invertElev; | |
| 374 |
2/2✓ Branch 0 taken 82414 times.
✓ Branch 1 taken 114365 times.
|
196779 | if ( (h1 - h2) > fabs(Conduit[k].slope) * Conduit[k].length ) |
| 375 | 82414 | Conduit[k].capacityLimited = TRUE; | |
| 376 | } | ||
| 377 | } | ||
| 378 | 431255 | } | |
| 379 | |||
| 380 | //============================================================================= | ||
| 381 | |||
| 382 | 986324 | void findLinkFlows(double dt) | |
| 383 | { | ||
| 384 | int i; | ||
| 385 | |||
| 386 | // --- find new flow in each non-dummy conduit | ||
| 387 | 986324 | #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 16282039 times.
✓ Branch 1 taken 986324 times.
|
17268363 | for ( i = 0; i < Nobjects[LINK]; i++) |
| 399 | { | ||
| 400 |
2/2✓ Branch 1 taken 15067777 times.
✓ Branch 2 taken 1214262 times.
|
16282039 | if ( isTrueConduit(i) ) updateNodeFlows(i); |
| 401 | } | ||
| 402 | |||
| 403 | // --- find new flows for all dummy conduits, pumps & regulators | ||
| 404 |
2/2✓ Branch 0 taken 16282039 times.
✓ Branch 1 taken 986324 times.
|
17268363 | for ( i = 0; i < Nobjects[LINK]; i++) |
| 405 | { | ||
| 406 |
2/2✓ Branch 1 taken 1214262 times.
✓ Branch 2 taken 15067777 times.
|
16282039 | if ( !isTrueConduit(i) ) |
| 407 | { | ||
| 408 |
2/2✓ Branch 0 taken 829186 times.
✓ Branch 1 taken 385076 times.
|
1214262 | if ( !Link[i].bypassed ) findNonConduitFlow(i, dt); |
| 409 | 1214262 | updateNodeFlows(i); | |
| 410 | } | ||
| 411 | } | ||
| 412 | 986324 | } | |
| 413 | |||
| 414 | //============================================================================= | ||
| 415 | |||
| 416 | 53595200 | int isTrueConduit(int j) | |
| 417 | { | ||
| 418 |
3/4✓ Branch 0 taken 49582423 times.
✓ Branch 1 taken 4012777 times.
✓ Branch 2 taken 49582423 times.
✗ Branch 3 not taken.
|
53595200 | return ( Link[j].type == CONDUIT && Link[j].xsect.type != DUMMY ); |
| 419 | } | ||
| 420 | |||
| 421 | //============================================================================= | ||
| 422 | |||
| 423 | 829186 | 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 | 829186 | qLast = Link[i].newFlow; | |
| 436 | 829186 | 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 | 829186 | qNew = link_getInflow(i); | |
| 441 |
2/2✓ Branch 0 taken 47655 times.
✓ Branch 1 taken 781531 times.
|
829186 | if ( Link[i].type == PUMP ) qNew = getModPumpFlow(i, qNew, dt); |
| 442 | |||
| 443 | // --- find surface area at each end of link | ||
| 444 | 829186 | 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 459195 times.
✓ Branch 1 taken 369991 times.
✓ Branch 2 taken 435301 times.
✓ Branch 3 taken 23894 times.
|
829186 | if ( Steps > 0 && Link[i].type != PUMP ) |
| 449 | { | ||
| 450 | 435301 | qNew = (1.0 - Omega) * qLast + Omega * qNew; | |
| 451 |
4/4✓ Branch 0 taken 57670 times.
✓ Branch 1 taken 377631 times.
✓ Branch 2 taken 28222 times.
✓ Branch 3 taken 29448 times.
|
435301 | if ( qNew * qLast < 0.0 ) qNew = 0.001 * SGN(qNew); |
| 452 | } | ||
| 453 | 829186 | Link[i].newFlow = qNew; | |
| 454 | 829186 | } | |
| 455 | |||
| 456 | //============================================================================= | ||
| 457 | |||
| 458 | 47655 | 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 | 47655 | int j = Link[i].node1; // pump's inlet node index | |
| 469 | 47655 | 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 21887 times.
✓ Branch 1 taken 25768 times.
|
47655 | if ( q == 0.0 ) return q; |
| 475 | |||
| 476 | // --- case where inlet node is a storage node: | ||
| 477 | // prevent node volume from going negative | ||
| 478 |
1/2✓ Branch 0 taken 25768 times.
✗ Branch 1 not taken.
|
25768 | if ( Node[j].type == STORAGE ) return node_getMaxOutflow(j, q, dt); |
| 479 | |||
| 480 | // --- case where inlet is a non-storage node | ||
| 481 | ✗ | 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 | ✗ | return q; | |
| 499 | } | ||
| 500 | |||
| 501 | //============================================================================= | ||
| 502 | |||
| 503 | 829186 | 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 160663 times.
✓ Branch 1 taken 668523 times.
|
829186 | if ( Link[i].type == ORIFICE ) |
| 512 | { | ||
| 513 | 160663 | Link[i].surfArea1 = Orifice[Link[i].subIndex].surfArea / 2.; | |
| 514 | } | ||
| 515 | |||
| 516 | // --- no surface area for weirs to maintain SWMM 4 compatibility | ||
| 517 | 668523 | else Link[i].surfArea1 = 0.0; | |
| 518 | |||
| 519 | 829186 | Link[i].surfArea2 = Link[i].surfArea1; | |
| 520 |
2/2✓ Branch 0 taken 822204 times.
✓ Branch 1 taken 6982 times.
|
829186 | if ( Link[i].flowClass == UP_CRITICAL || |
| 521 |
2/2✓ Branch 0 taken 365240 times.
✓ Branch 1 taken 456964 times.
|
829186 | Node[Link[i].node1].type == STORAGE ) Link[i].surfArea1 = 0.0; |
| 522 |
2/2✓ Branch 0 taken 676077 times.
✓ Branch 1 taken 153109 times.
|
829186 | if ( Link[i].flowClass == DN_CRITICAL || |
| 523 |
2/2✓ Branch 0 taken 2881 times.
✓ Branch 1 taken 673196 times.
|
829186 | Node[Link[i].node2].type == STORAGE ) Link[i].surfArea2 = 0.0; |
| 524 | 829186 | } | |
| 525 | |||
| 526 | //============================================================================= | ||
| 527 | |||
| 528 | 16282039 | 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 | 16282039 | int barrels = 1; | |
| 538 | 16282039 | int n1 = Link[i].node1; | |
| 539 | 16282039 | int n2 = Link[i].node2; | |
| 540 | 16282039 | double q = Link[i].newFlow; | |
| 541 | 16282039 | double conduitLossRate = 0.0; | |
| 542 | |||
| 543 | // --- update total inflow & outflow at upstream/downstream nodes | ||
| 544 |
2/2✓ Branch 0 taken 16119701 times.
✓ Branch 1 taken 162338 times.
|
16282039 | if ( q >= 0.0 ) |
| 545 | { | ||
| 546 | 16119701 | Node[n1].outflow += q; | |
| 547 | 16119701 | Node[n2].inflow += q; | |
| 548 | } | ||
| 549 | else | ||
| 550 | { | ||
| 551 | 162338 | Node[n1].inflow -= q; | |
| 552 | 162338 | Node[n2].outflow -= q; | |
| 553 | } | ||
| 554 | |||
| 555 | // --- add any uniform evap & seepage loss from conduit link | ||
| 556 |
2/2✓ Branch 0 taken 15067777 times.
✓ Branch 1 taken 1214262 times.
|
16282039 | if ( Link[i].type == CONDUIT ) |
| 557 | { | ||
| 558 | 15067777 | k = Link[i].subIndex; | |
| 559 | 15067777 | barrels = Conduit[k].barrels; | |
| 560 | 15067777 | conduitLossRate = (Conduit[k].evapLossRate + Conduit[k].seepLossRate) * | |
| 561 | barrels; | ||
| 562 |
2/2✓ Branch 0 taken 1993944 times.
✓ Branch 1 taken 13073833 times.
|
15067777 | if (conduitLossRate > 0.0) |
| 563 | { | ||
| 564 | // --- outfall nodes do not share evap & seepage losses | ||
| 565 |
3/4✓ Branch 0 taken 1993944 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 978291 times.
✓ Branch 3 taken 1015653 times.
|
1993944 | if (Node[n1].type != OUTFALL && Node[n2].type != OUTFALL) |
| 566 | 978291 | conduitLossRate /= 2.0; | |
| 567 |
1/2✓ Branch 0 taken 1993944 times.
✗ Branch 1 not taken.
|
1993944 | if (Node[n1].type != OUTFALL) |
| 568 | 1993944 | Node[n1].outflow += conduitLossRate; | |
| 569 |
2/2✓ Branch 0 taken 978291 times.
✓ Branch 1 taken 1015653 times.
|
1993944 | if (Node[n2].type != OUTFALL) |
| 570 | 978291 | Node[n2].outflow += conduitLossRate; | |
| 571 | } | ||
| 572 | } | ||
| 573 | |||
| 574 | // --- add surf. area contributions to upstream/downstream nodes | ||
| 575 | 16282039 | Xnode[Link[i].node1].newSurfArea += Link[i].surfArea1 * barrels; | |
| 576 | 16282039 | Xnode[Link[i].node2].newSurfArea += Link[i].surfArea2 * barrels; | |
| 577 | |||
| 578 | // --- update summed value of dqdh at each end node | ||
| 579 | 16282039 | Xnode[Link[i].node1].sumdqdh += Link[i].dqdh; | |
| 580 |
2/2✓ Branch 0 taken 168465 times.
✓ Branch 1 taken 16113574 times.
|
16282039 | if ( Link[i].type == PUMP ) |
| 581 | { | ||
| 582 | 168465 | k = Link[i].subIndex; | |
| 583 |
1/2✓ Branch 0 taken 168465 times.
✗ Branch 1 not taken.
|
168465 | if ( Pump[k].type != TYPE4_PUMP ) |
| 584 | { | ||
| 585 | 168465 | Xnode[n2].sumdqdh += Link[i].dqdh; | |
| 586 | } | ||
| 587 | } | ||
| 588 | 16113574 | else Xnode[n2].sumdqdh += Link[i].dqdh; | |
| 589 | 16282039 | } | |
| 590 | |||
| 591 | //============================================================================= | ||
| 592 | |||
| 593 | 986324 | 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 | 986324 | 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 16282039 times.
✓ Branch 2 taken 986324 times.
|
17268363 | 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 | 986324 | #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 8325766 times.
✓ Branch 1 taken 802748 times.
|
9128514 | for (i = 0; i < Nobjects[NODE]; i++) |
| 627 | { | ||
| 628 |
2/2✓ Branch 0 taken 1107440 times.
✓ Branch 1 taken 7218326 times.
|
8325766 | if ( Node[i].type == OUTFALL ) continue; |
| 629 |
2/2✓ Branch 0 taken 183576 times.
✓ Branch 1 taken 7034750 times.
|
7218326 | if (Xnode[i].converged == FALSE) return FALSE; |
| 630 | } | ||
| 631 | 802748 | return TRUE; | |
| 632 | } | ||
| 633 | |||
| 634 | //============================================================================= | ||
| 635 | |||
| 636 | 15324514 | 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 | 15324514 | 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 1250773 times.
✓ Branch 1 taken 14073741 times.
✓ Branch 2 taken 864615 times.
✓ Branch 3 taken 386158 times.
|
15324514 | canPond = (AllowPonding && Node[i].pondedArea > 0.0); |
| 662 |
4/4✓ Branch 0 taken 864615 times.
✓ Branch 1 taken 14459899 times.
✓ Branch 2 taken 177183 times.
✓ Branch 3 taken 687432 times.
|
15324514 | isPonded = (canPond && Node[i].newDepth > Node[i].fullDepth); |
| 663 | |||
| 664 | // --- initialize values | ||
| 665 | 15324514 | yCrown = Node[i].crownElev - Node[i].invertElev; | |
| 666 | 15324514 | yOld = Node[i].oldDepth; | |
| 667 | 15324514 | yLast = Node[i].newDepth; | |
| 668 | 15324514 | Node[i].overflow = 0.0; | |
| 669 | 15324514 | surfArea = Xnode[i].newSurfArea; | |
| 670 |
2/2✓ Branch 0 taken 13010540 times.
✓ Branch 1 taken 2313974 times.
|
15324514 | surfArea = MAX(surfArea, MinSurfArea); |
| 671 | |||
| 672 | // --- determine average net flow volume into node over the time step | ||
| 673 | 15324514 | dQ = Node[i].inflow - Node[i].outflow; | |
| 674 | 15324514 | dV = 0.5 * (Node[i].oldNetInflow + dQ) * dt; | |
| 675 | |||
| 676 | // --- determine if node is EXTRAN surcharged | ||
| 677 |
2/2✓ Branch 0 taken 15307530 times.
✓ Branch 1 taken 16984 times.
|
15324514 | if (SurchargeMethod == EXTRAN) |
| 678 | { | ||
| 679 | // --- ponded nodes don't surcharge | ||
| 680 |
2/2✓ Branch 0 taken 177183 times.
✓ Branch 1 taken 15130347 times.
|
15307530 | if (isPonded) isSurcharged = FALSE; |
| 681 | |||
| 682 | // --- closed storage units that are full are in surcharge | ||
| 683 |
2/2✓ Branch 0 taken 559889 times.
✓ Branch 1 taken 14570458 times.
|
15130347 | else if (Node[i].type == STORAGE) |
| 684 | { | ||
| 685 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 559889 times.
|
559889 | 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 14566334 times.
✓ Branch 1 taken 4124 times.
✓ Branch 2 taken 361084 times.
✓ Branch 3 taken 14205250 times.
|
14570458 | 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 14963430 times.
✓ Branch 1 taken 361084 times.
|
15324514 | if (!isSurcharged) |
| 695 | { | ||
| 696 | 14963430 | dy = dV / surfArea; | |
| 697 | 14963430 | yNew = yOld + dy; | |
| 698 | |||
| 699 | // --- save non-ponded surface area for use in surcharge algorithm | ||
| 700 |
2/2✓ Branch 0 taken 14786247 times.
✓ Branch 1 taken 177183 times.
|
14963430 | if ( !isPonded ) Xnode[i].oldSurfArea = surfArea; |
| 701 | |||
| 702 | // --- apply under-relaxation to new depth estimate | ||
| 703 |
2/2✓ Branch 0 taken 10604278 times.
✓ Branch 1 taken 4359152 times.
|
14963430 | if ( Steps > 0 ) |
| 704 | { | ||
| 705 | 10604278 | 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 14786247 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 177178 times.
|
14963430 | 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 | 361084 | corr = 1.0; | |
| 720 |
2/2✓ Branch 0 taken 144370 times.
✓ Branch 1 taken 216714 times.
|
361084 | 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 | 361084 | denom = Xnode[i].sumdqdh; | |
| 725 |
2/2✓ Branch 0 taken 23283 times.
✓ Branch 1 taken 337801 times.
|
361084 | if ( yLast < 1.25 * yCrown ) |
| 726 | { | ||
| 727 | 23283 | f = (yLast - yCrown) / yCrown; | |
| 728 | 23283 | denom += (Xnode[i].oldSurfArea/dt - | |
| 729 | 23283 | 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 361084 times.
|
361084 | if ( denom == 0.0 ) dy = 0.0; |
| 734 | 361084 | else dy = corr * dQ / denom; | |
| 735 | 361084 | yNew = yLast + dy; | |
| 736 |
2/2✓ Branch 0 taken 119 times.
✓ Branch 1 taken 360965 times.
|
361084 | 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 361084 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
361084 | 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 984 times.
✓ Branch 1 taken 15323530 times.
|
15324514 | if ( yNew < 0 ) yNew = 0.0; |
| 745 | |||
| 746 | // --- determine max. non-flooded depth | ||
| 747 | 15324514 | yMax = Node[i].fullDepth; | |
| 748 |
2/2✓ Branch 0 taken 14459899 times.
✓ Branch 1 taken 864615 times.
|
15324514 | if ( canPond == FALSE ) yMax += Node[i].surDepth; |
| 749 | |||
| 750 | // --- find flooded depth & volume | ||
| 751 |
2/2✓ Branch 0 taken 264586 times.
✓ Branch 1 taken 15059928 times.
|
15324514 | if ( yNew > yMax ) |
| 752 | { | ||
| 753 | 264586 | yNew = getFloodedDepth(i, canPond, dV, yNew, yMax, dt); | |
| 754 | } | ||
| 755 | 15059928 | else Node[i].newVolume = node_getVolume(i, yNew); | |
| 756 | |||
| 757 | // --- compute change in depth w.r.t. time | ||
| 758 | 15324514 | Xnode[i].dYdT = fabs(yNew - yOld) / dt; | |
| 759 | |||
| 760 | // --- save new depth for node | ||
| 761 | 15324514 | Node[i].newDepth = yNew; | |
| 762 | 15324514 | } | |
| 763 | |||
| 764 | //============================================================================= | ||
| 765 | |||
| 766 | 264586 | 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 87403 times.
✓ Branch 1 taken 177183 times.
|
264586 | if ( canPond == FALSE ) |
| 781 | { | ||
| 782 | 87403 | Node[i].overflow = dV / dt; | |
| 783 | 87403 | Node[i].newVolume = Node[i].fullVolume; | |
| 784 | 87403 | 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 141441 times.
|
264586 | if ( Node[i].overflow < FUDGE ) Node[i].overflow = 0.0; |
| 793 | 264586 | return yNew; | |
| 794 | |||
| 795 | } | ||
| 796 | |||
| 797 | //============================================================================= | ||
| 798 | |||
| 799 | 250504 | 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 | 250504 | int minLink = -1; // index of link w/ min. time step | |
| 808 | 250504 | 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 | 250504 | tMin = maxStep; | |
| 815 | 250504 | tMinLink = getLinkStep(tMin, &minLink); | |
| 816 | 250504 | tMinNode = getNodeStep(tMinLink, &minNode); | |
| 817 | |||
| 818 | // --- use smaller of the link and node time step | ||
| 819 | 250504 | tMin = tMinLink; | |
| 820 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 250504 times.
|
250504 | 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 | 250504 | 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 250469 times.
|
250504 | if ( tMin < MinRouteStep ) tMin = MinRouteStep; |
| 831 | 250504 | return tMin; | |
| 832 | } | ||
| 833 | |||
| 834 | //============================================================================= | ||
| 835 | |||
| 836 | 250504 | 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 | 250504 | double tLink = tMin; // critical link time step (sec) | |
| 849 | |||
| 850 | // --- examine each conduit link | ||
| 851 |
2/2✓ Branch 0 taken 2423320 times.
✓ Branch 1 taken 250504 times.
|
2673824 | for ( i = 0; i < Nobjects[LINK]; i++ ) |
| 852 | { | ||
| 853 |
2/2✓ Branch 0 taken 2066305 times.
✓ Branch 1 taken 357015 times.
|
2423320 | if ( Link[i].type == CONDUIT ) |
| 854 | { | ||
| 855 | // --- skip conduits with negligible flow, area or Fr | ||
| 856 | 2066305 | k = Link[i].subIndex; | |
| 857 | 2066305 | q = fabs(Link[i].newFlow) / Conduit[k].barrels; | |
| 858 |
2/2✓ Branch 0 taken 1862433 times.
✓ Branch 1 taken 203872 times.
|
2066305 | if ( q <= FUDGE |
| 859 |
1/2✓ Branch 0 taken 1862433 times.
✗ Branch 1 not taken.
|
1862433 | || Conduit[k].a1 <= FUDGE |
| 860 |
2/2✓ Branch 0 taken 172410 times.
✓ Branch 1 taken 1690023 times.
|
1862433 | || Link[i].froude <= 0.01 |
| 861 | 376282 | ) continue; | |
| 862 | |||
| 863 | // --- compute time step to satisfy Courant condition | ||
| 864 | 1690023 | t = Link[i].newVolume / Conduit[k].barrels / q; | |
| 865 | 1690023 | t = t * Conduit[k].modLength / link_getLength(i); | |
| 866 | 1690023 | t = t * Link[i].froude / (1.0 + Link[i].froude) * CourantFactor; | |
| 867 | |||
| 868 | // --- update critical link time step | ||
| 869 |
2/2✓ Branch 0 taken 6285 times.
✓ Branch 1 taken 1683738 times.
|
1690023 | if ( t < tLink ) |
| 870 | { | ||
| 871 | 6285 | tLink = t; | |
| 872 | 6285 | *minLink = i; | |
| 873 | } | ||
| 874 | } | ||
| 875 | } | ||
| 876 | 250504 | return tLink; | |
| 877 | } | ||
| 878 | |||
| 879 | //============================================================================= | ||
| 880 | |||
| 881 | 250504 | 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 | 250504 | 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 2631166 times.
✓ Branch 1 taken 250504 times.
|
2881670 | for ( i = 0; i < Nobjects[NODE]; i++ ) |
| 899 | { | ||
| 900 | // --- see if node can be skipped | ||
| 901 |
2/2✓ Branch 0 taken 500389 times.
✓ Branch 1 taken 2130777 times.
|
2631166 | if ( Node[i].type == OUTFALL ) continue; |
| 902 |
2/2✓ Branch 0 taken 129963 times.
✓ Branch 1 taken 2000814 times.
|
2130777 | if ( Node[i].newDepth <= FUDGE) continue; |
| 903 | 2239373 | if ( Node[i].newDepth + FUDGE >= | |
| 904 |
2/2✓ Branch 0 taken 238559 times.
✓ Branch 1 taken 1762255 times.
|
2000814 | Node[i].crownElev - Node[i].invertElev ) continue; |
| 905 | |||
| 906 | // --- define max. allowable depth change using crown elevation | ||
| 907 | 1762255 | maxDepth = (Node[i].crownElev - Node[i].invertElev) * 0.25; | |
| 908 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 1762255 times.
|
1762255 | if ( maxDepth < FUDGE ) continue; |
| 909 | 1762255 | dYdT = Xnode[i].dYdT; | |
| 910 |
2/2✓ Branch 0 taken 1297079 times.
✓ Branch 1 taken 465176 times.
|
1762255 | if (dYdT < FUDGE ) continue; |
| 911 | |||
| 912 | // --- compute time to reach max. depth & compare with critical time | ||
| 913 | 465176 | t1 = maxDepth / dYdT; | |
| 914 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 465176 times.
|
465176 | if ( t1 < tNode ) |
| 915 | { | ||
| 916 | ✗ | tNode = t1; | |
| 917 | ✗ | *minNode = i; | |
| 918 | } | ||
| 919 | } | ||
| 920 | 250504 | return tNode; | |
| 921 | } | ||
| 922 |