GCC Code Coverage Report


Directory: src/solver/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 89.8% 237 / 0 / 264
Functions: 100.0% 17 / 0 / 17
Branches: 78.5% 164 / 0 / 209

flowrout.c
Line Branch Exec Source
1 //-----------------------------------------------------------------------------
2 // flowrout.c
3 //
4 // Project: EPA SWMM5
5 // Version: 5.2
6 // Date: 06/12/23 (Build 5.2.4)
7 // Author: L. Rossman
8 // M. Tryby (EPA)
9 //
10 // Flow routing functions.
11 //
12 // Update History
13 // ==============
14 // Build 5.1.007:
15 // - updateStorageState() modified in response to node outflow being
16 // initialized with current evap & seepage losses in routing_execute().
17 // Build 5.1.008:
18 // - Determination of node crown elevations moved to dynwave.c.
19 // - Support added for new way of recording conduit's fullness state.
20 // Build 5.1.012:
21 // - Overflow computed in updateStorageState() must be non-negative.
22 // - Terminal storage nodes now updated corectly.
23 // Build 5.1.014:
24 // - Arguments to function link_getLossRate changed.
25 // Build 5.2.0:
26 // - Correction made to updating state of terminal storage nodes.
27 // Build 5.2.1:
28 // - For storage routing, after convergence the reported depth is now
29 // based on the last volume found rather than the next trial depth.
30 // Build 5.2.4:
31 // - Arguments to link_getLossRate changed.
32 //-----------------------------------------------------------------------------
33 #define _CRT_SECURE_NO_DEPRECATE
34
35 #include <stdlib.h>
36 #include <math.h>
37 #include "headers.h"
38
39 //-----------------------------------------------------------------------------
40 // Constants
41 //-----------------------------------------------------------------------------
42 static const double OMEGA = 0.55; // under-relaxation parameter
43 static const int MAXITER = 10; // max. iterations for storage updating
44 static const double STOPTOL = 0.005; // storage updating stopping tolerance
45
46 //-----------------------------------------------------------------------------
47 // External functions (declared in funcs.h)
48 //-----------------------------------------------------------------------------
49 // flowrout_init (called by routing_open)
50 // flowrout_close (called by routing_close)
51 // flowrout_getRoutingStep (called routing_getRoutingStep)
52 // flowrout_execute (called routing_execute)
53
54 //-----------------------------------------------------------------------------
55 // Local functions
56 //-----------------------------------------------------------------------------
57 static void initLinkDepths(void);
58 static void initNodeDepths(void);
59 static void initNodes(void);
60 static void initLinks(int routingModel);
61 static void validateTreeLayout(void);
62 static void validateGeneralLayout(void);
63 static void updateStorageState(int i, int j, int links[], double dt);
64 static double getStorageOutflow(int node, int j, int links[], double dt);
65 static double getLinkInflow(int link, double dt);
66 static void setNewNodeState(int node, double dt);
67 static void setNewLinkState(int link);
68 static void updateNodeDepth(int node, double y);
69 static int steadyflow_execute(int link, double* qin, double* qout,
70 double tStep);
71
72
73 //=============================================================================
74
75 58 void flowrout_init(int routingModel)
76 //
77 // Input: routingModel = routing model code
78 // Output: none
79 // Purpose: initializes flow routing system.
80 //
81 {
82 // --- initialize for dynamic wave routing
83
2/2
✓ Branch 0 taken 41 times.
✓ Branch 1 taken 17 times.
58 if ( routingModel == DW )
84 {
85 // --- check for valid conveyance network layout
86 41 validateGeneralLayout();
87 41 dynwave_init();
88
89 // --- initialize node & link depths if not using a hotstart file
90
2/2
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 1 time.
41 if ( Fhotstart1.mode == NO_FILE )
91 {
92 40 initNodeDepths();
93 40 initLinkDepths();
94 }
95 }
96
97 // --- validate network layout for kinematic wave routing
98 17 else validateTreeLayout();
99
100 // --- initialize node & link volumes
101 58 initNodes();
102 58 initLinks(routingModel);
103 58 }
104
105 //=============================================================================
106
107 58 void flowrout_close(int routingModel)
108 //
109 // Input: routingModel = routing method code
110 // Output: none
111 // Purpose: closes down routing method used.
112 //
113 {
114
2/2
✓ Branch 0 taken 41 times.
✓ Branch 1 taken 17 times.
58 if ( routingModel == DW ) dynwave_close();
115 58 }
116
117 //=============================================================================
118
119 988027 double flowrout_getRoutingStep(int routingModel, double fixedStep)
120 //
121 // Input: routingModel = type of routing method used
122 // fixedStep = user-assigned max. routing step (sec)
123 // Output: returns adjusted value of routing time step (sec)
124 // Purpose: finds variable time step for dynamic wave routing.
125 //
126 {
127
2/2
✓ Branch 0 taken 367748 times.
✓ Branch 1 taken 620279 times.
988027 if ( routingModel == DW )
128 {
129 367748 return dynwave_getRoutingStep(fixedStep);
130 }
131 620279 return fixedStep;
132 }
133
134 //=============================================================================
135
136 988013 int flowrout_execute(int links[], int routingModel, double tStep)
137 //
138 // Input: links = array of link indexes in topo-sorted order (per routing model)
139 // routingModel = type of routing method used
140 // tStep = routing time step (sec)
141 // Output: returns number of computational steps taken
142 // Purpose: routes flow through conveyance network over current time step.
143 //
144 {
145 int i, j;
146 int n1; // upstream node of link
147 double qin; // link inflow (cfs)
148 double qout; // link outflow (cfs)
149 double steps; // computational step count
150
151 // --- set overflows to drain any ponded water
152
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 988013 times.
988013 if ( ErrorCode ) return 0;
153
2/2
✓ Branch 0 taken 37851987 times.
✓ Branch 1 taken 988013 times.
38840000 for (j = 0; j < Nobjects[NODE]; j++)
154 {
155 37851987 Node[j].updated = FALSE;
156 37851987 Node[j].overflow = 0.0;
157
2/2
✓ Branch 0 taken 37708857 times.
✓ Branch 1 taken 143130 times.
37851987 if ( Node[j].type != STORAGE
158
2/2
✓ Branch 0 taken 88273 times.
✓ Branch 1 taken 37620584 times.
37708857 && Node[j].newVolume > Node[j].fullVolume )
159 {
160 88273 Node[j].overflow = (Node[j].newVolume - Node[j].fullVolume)/tStep;
161 }
162 }
163
164 // --- execute dynamic wave routing if called for
165
2/2
✓ Branch 0 taken 367734 times.
✓ Branch 1 taken 620279 times.
988013 if ( routingModel == DW )
166 {
167 367734 return dynwave_execute(tStep);
168 }
169
170 // --- otherwise examine each link, moving from upstream to downstream
171 620279 steps = 0.0;
172
2/2
✓ Branch 0 taken 11995299 times.
✓ Branch 1 taken 620279 times.
12615578 for (i = 0; i < Nobjects[LINK]; i++)
173 {
174 // --- see if upstream node is a storage unit whose state needs updating
175 11995299 j = links[i];
176 11995299 n1 = Link[j].node1;
177
2/2
✓ Branch 0 taken 14040 times.
✓ Branch 1 taken 11981259 times.
11995299 if ( Node[n1].type == STORAGE ) updateStorageState(n1, i, links, tStep);
178
179 // --- retrieve inflow at upstream end of link
180 11995299 qin = getLinkInflow(j, tStep);
181
182 // --- route flow through link
183
2/2
✓ Branch 0 taken 11894379 times.
✓ Branch 1 taken 100920 times.
11995299 if ( routingModel == SF )
184 11894379 steps += steadyflow_execute(j, &qin, &qout, tStep);
185 else
186 100920 steps += kinwave_execute(j, &qin, &qout, tStep);
187 11995299 Link[j].newFlow = qout;
188
189 // adjust outflow at upstream node and inflow at downstream node
190 11995299 Node[ Link[j].node1 ].outflow += qin;
191 11995299 Node[ Link[j].node2 ].inflow += qout;
192 }
193
1/2
✓ Branch 0 taken 620279 times.
✗ Branch 1 not taken.
620279 if ( Nobjects[LINK] > 0 ) steps /= Nobjects[LINK];
194
195 // --- update state of each non-updated node and link
196
2/2
✓ Branch 1 taken 12629258 times.
✓ Branch 2 taken 620279 times.
13249537 for ( j=0; j<Nobjects[NODE]; j++) setNewNodeState(j, tStep);
197
2/2
✓ Branch 1 taken 11995299 times.
✓ Branch 2 taken 620279 times.
12615578 for ( j=0; j<Nobjects[LINK]; j++) setNewLinkState(j);
198 620279 return (int)(steps+0.5);
199 }
200
201 //=============================================================================
202
203 17 void validateTreeLayout()
204 //
205 // Input: none
206 // Output: none
207 // Purpose: validates tree-like conveyance system layout used for Steady
208 // and Kinematic Wave flow routing
209 //
210 {
211 int j;
212
213 // --- check nodes
214
2/2
✓ Branch 0 taken 89 times.
✓ Branch 1 taken 17 times.
106 for ( j = 0; j < Nobjects[NODE]; j++ )
215 {
216
4/4
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 34 times.
✓ Branch 2 taken 11 times.
✓ Branch 3 taken 40 times.
89 switch ( Node[j].type )
217 {
218 // --- dividers must have only 2 outlet links
219 4 case DIVIDER:
220
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
4 if ( Node[j].degree > 2 )
221 {
222 report_writeErrorMsg(ERR_DIVIDER, Node[j].ID);
223 }
224 4 break;
225
226 // --- outfalls cannot have any outlet links
227 34 case OUTFALL:
228
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 34 times.
34 if ( Node[j].degree > 0 )
229 {
230 report_writeErrorMsg(ERR_OUTFALL, Node[j].ID);
231 }
232 34 break;
233
234 // --- storage nodes can have multiple outlets
235 11 case STORAGE: break;
236
237 // --- all other nodes allowed only one outlet link
238 40 default:
239
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 40 times.
40 if ( Node[j].degree > 1 )
240 {
241 report_writeErrorMsg(ERR_MULTI_OUTLET, Node[j].ID);
242 }
243 }
244 }
245
246 // --- check links
247
2/2
✓ Branch 0 taken 59 times.
✓ Branch 1 taken 17 times.
76 for (j=0; j<Nobjects[LINK]; j++)
248 {
249
2/3
✓ Branch 0 taken 54 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
59 switch ( Link[j].type )
250 {
251 // --- non-dummy conduits cannot have adverse slope
252 54 case CONDUIT:
253
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 54 times.
54 if ( Conduit[Link[j].subIndex].slope < 0.0 &&
254 Link[j].xsect.type != DUMMY )
255 {
256 report_writeErrorMsg(ERR_SLOPE, Link[j].ID);
257 }
258 54 break;
259
260 // --- regulator links must be outlets of storage nodes
261 case ORIFICE:
262 case WEIR:
263 case OUTLET:
264 if ( Node[Link[j].node1].type != STORAGE )
265 {
266 report_writeErrorMsg(ERR_REGULATOR, Link[j].ID);
267 }
268 }
269 }
270 17 }
271
272 //=============================================================================
273
274 41 void validateGeneralLayout()
275 //
276 // Input: none
277 // Output: nonw
278 // Purpose: validates general conveyance system layout.
279 //
280 {
281 int i, j;
282 41 int outletCount = 0;
283
284 // --- use node inflow attribute to count inflow connections
285
2/2
✓ Branch 0 taken 10094 times.
✓ Branch 1 taken 41 times.
10135 for ( i=0; i<Nobjects[NODE]; i++ ) Node[i].inflow = 0.0;
286
287 // --- examine each link
288
2/2
✓ Branch 0 taken 10449 times.
✓ Branch 1 taken 41 times.
10490 for ( j = 0; j < Nobjects[LINK]; j++ )
289 {
290 // --- update inflow link count of downstream node
291 10449 i = Link[j].node1;
292
2/2
✓ Branch 0 taken 10442 times.
✓ Branch 1 taken 7 times.
10449 if ( Node[i].type != OUTFALL ) i = Link[j].node2;
293 10449 Node[i].inflow += 1.0;
294
295 // --- if link is dummy link or ideal pump then it must
296 // be the only link exiting the upstream node
297
3/4
✓ Branch 0 taken 9718 times.
✓ Branch 1 taken 731 times.
✓ Branch 2 taken 9718 times.
✗ Branch 3 not taken.
10449 if ( (Link[j].type == CONDUIT && Link[j].xsect.type == DUMMY) ||
298
2/2
✓ Branch 0 taken 101 times.
✓ Branch 1 taken 10348 times.
10449 (Link[j].type == PUMP &&
299
2/2
✓ Branch 0 taken 23 times.
✓ Branch 1 taken 78 times.
101 Pump[Link[j].subIndex].type == IDEAL_PUMP) )
300 {
301 23 i = Link[j].node1;
302
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 23 times.
23 if ( Link[j].direction < 0 ) i = Link[j].node2;
303
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 23 times.
23 if ( Node[i].degree > 1 )
304 {
305 report_writeErrorMsg(ERR_DUMMY_LINK, Node[i].ID);
306 }
307 }
308 }
309
310 // --- check each node to see if it qualifies as an outlet node
311 // (meaning that degree = 0)
312
2/2
✓ Branch 0 taken 10094 times.
✓ Branch 1 taken 41 times.
10135 for ( i = 0; i < Nobjects[NODE]; i++ )
313 {
314 // --- if node is of type Outfall, check that it has only 1
315 // connecting link (which can either be an outflow or inflow link)
316
2/2
✓ Branch 0 taken 314 times.
✓ Branch 1 taken 9780 times.
10094 if ( Node[i].type == OUTFALL )
317 {
318
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 314 times.
314 if ( Node[i].degree + (int)Node[i].inflow > 1 )
319 {
320 report_writeErrorMsg(ERR_OUTFALL, Node[i].ID);
321 }
322 314 else outletCount++;
323 }
324 }
325
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 41 times.
41 if ( outletCount == 0 ) report_writeErrorMsg(ERR_NO_OUTLETS, "");
326
327 // --- reset node inflows back to zero
328
2/2
✓ Branch 0 taken 10094 times.
✓ Branch 1 taken 41 times.
10135 for ( i = 0; i < Nobjects[NODE]; i++ )
329 {
330
2/2
✓ Branch 0 taken 1334 times.
✓ Branch 1 taken 8760 times.
10094 if ( Node[i].inflow == 0.0 ) Node[i].degree = -Node[i].degree;
331 10094 Node[i].inflow = 0.0;
332 }
333 41 }
334
335 //=============================================================================
336
337 40 void initNodeDepths(void)
338 //
339 // Input: none
340 // Output: none
341 // Purpose: sets initial depth at nodes for Dynamic Wave flow routing.
342 //
343 {
344 int i; // link or node index
345 int n; // node index
346 double y; // node water depth (ft)
347
348 // --- use Node[].inflow as a temporary accumulator for depth in
349 // connecting links and Node[].outflow as a temporary counter
350 // for the number of connecting links
351
2/2
✓ Branch 0 taken 9178 times.
✓ Branch 1 taken 40 times.
9218 for (i = 0; i < Nobjects[NODE]; i++)
352 {
353 9178 Node[i].inflow = 0.0;
354 9178 Node[i].outflow = 0.0;
355 }
356
357 // --- total up flow depths in all connecting links into nodes
358
2/2
✓ Branch 0 taken 9482 times.
✓ Branch 1 taken 40 times.
9522 for (i = 0; i < Nobjects[LINK]; i++)
359 {
360
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 9474 times.
9482 if ( Link[i].newDepth > FUDGE ) y = Link[i].newDepth + Link[i].offset1;
361 9474 else y = 0.0;
362 9482 n = Link[i].node1;
363 9482 Node[n].inflow += y;
364 9482 Node[n].outflow += 1.0;
365 9482 n = Link[i].node2;
366 9482 Node[n].inflow += y;
367 9482 Node[n].outflow += 1.0;
368 }
369
370 // --- if no user-supplied depth then set initial depth at non-storage/
371 // non-outfall nodes to average of depths in connecting links
372
2/2
✓ Branch 0 taken 9178 times.
✓ Branch 1 taken 40 times.
9218 for ( i = 0; i < Nobjects[NODE]; i++ )
373 {
374
2/2
✓ Branch 0 taken 278 times.
✓ Branch 1 taken 8900 times.
9178 if ( Node[i].type == OUTFALL ) continue;
375
2/2
✓ Branch 0 taken 94 times.
✓ Branch 1 taken 8806 times.
8900 if ( Node[i].type == STORAGE ) continue;
376
2/2
✓ Branch 0 taken 51 times.
✓ Branch 1 taken 8755 times.
8806 if ( Node[i].initDepth > 0.0 ) continue;
377
2/2
✓ Branch 0 taken 8751 times.
✓ Branch 1 taken 4 times.
8755 if ( Node[i].outflow > 0.0 )
378 {
379 8751 Node[i].newDepth = Node[i].inflow / Node[i].outflow;
380 }
381 }
382
383 // --- compute initial depths at all outfall nodes
384
2/2
✓ Branch 1 taken 9482 times.
✓ Branch 2 taken 40 times.
9522 for ( i = 0; i < Nobjects[LINK]; i++ ) link_setOutfallDepth(i);
385 40 }
386
387 //=============================================================================
388
389 40 void initLinkDepths()
390 //
391 // Input: none
392 // Output: none
393 // Purpose: sets initial flow depths in conduits under Dyn. Wave routing.
394 //
395 {
396 int i; // link index
397 double y, y1, y2; // depths (ft)
398
399 // --- examine each link
400
2/2
✓ Branch 0 taken 9482 times.
✓ Branch 1 taken 40 times.
9522 for (i = 0; i < Nobjects[LINK]; i++)
401 {
402 // --- examine each conduit
403
2/2
✓ Branch 0 taken 8852 times.
✓ Branch 1 taken 630 times.
9482 if ( Link[i].type == CONDUIT )
404 {
405 // --- skip conduits with user-assigned initial flows
406 // (their depths have already been set to normal depth)
407
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8844 times.
8852 if ( Link[i].q0 != 0.0 ) continue;
408
409 // --- set depth to average of depths at end nodes
410 8844 y1 = Node[Link[i].node1].newDepth - Link[i].offset1;
411
2/2
✓ Branch 0 taken 8446 times.
✓ Branch 1 taken 398 times.
8844 y1 = MAX(y1, 0.0);
412
2/2
✓ Branch 0 taken 8827 times.
✓ Branch 1 taken 17 times.
8844 y1 = MIN(y1, Link[i].xsect.yFull);
413 8844 y2 = Node[Link[i].node2].newDepth - Link[i].offset2;
414
2/2
✓ Branch 0 taken 7516 times.
✓ Branch 1 taken 1328 times.
8844 y2 = MAX(y2, 0.0);
415
2/2
✓ Branch 0 taken 8709 times.
✓ Branch 1 taken 135 times.
8844 y2 = MIN(y2, Link[i].xsect.yFull);
416 8844 y = 0.5 * (y1 + y2);
417
2/2
✓ Branch 0 taken 258 times.
✓ Branch 1 taken 8586 times.
8844 y = MAX(y, FUDGE);
418 8844 Link[i].newDepth = y;
419 }
420 }
421 40 }
422
423 //=============================================================================
424
425 58 void initNodes()
426 //
427 // Input: none
428 // Output: none
429 // Purpose: sets initial inflow/outflow and volume for each node
430 //
431 {
432 int i;
433
434
2/2
✓ Branch 0 taken 10183 times.
✓ Branch 1 taken 58 times.
10241 for ( i = 0; i < Nobjects[NODE]; i++ )
435 {
436 // --- initialize node inflow and outflow
437 10183 Node[i].inflow = Node[i].newLatFlow;
438 10183 Node[i].outflow = 0.0;
439
440 // --- initialize node volume
441 10183 Node[i].newVolume = 0.0;
442
2/2
✓ Branch 0 taken 6420 times.
✓ Branch 1 taken 3763 times.
10183 if ( AllowPonding &&
443
2/2
✓ Branch 0 taken 5269 times.
✓ Branch 1 taken 1151 times.
6420 Node[i].pondedArea > 0.0 &&
444
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5269 times.
5269 Node[i].newDepth > Node[i].fullDepth )
445 {
446 Node[i].newVolume = Node[i].fullVolume +
447 (Node[i].newDepth - Node[i].fullDepth) *
448 Node[i].pondedArea;
449 }
450 10183 else Node[i].newVolume = node_getVolume(i, Node[i].newDepth);
451 }
452
453 // --- update nodal inflow/outflow at ends of each link
454 // (needed for Steady Flow & Kin. Wave routing)
455
2/2
✓ Branch 0 taken 10508 times.
✓ Branch 1 taken 58 times.
10566 for ( i = 0; i < Nobjects[LINK]; i++ )
456 {
457
1/2
✓ Branch 0 taken 10508 times.
✗ Branch 1 not taken.
10508 if ( Link[i].newFlow >= 0.0 )
458 {
459 10508 Node[Link[i].node1].outflow += Link[i].newFlow;
460 10508 Node[Link[i].node2].inflow += Link[i].newFlow;
461 }
462 else
463 {
464 Node[Link[i].node1].inflow -= Link[i].newFlow;
465 Node[Link[i].node2].outflow -= Link[i].newFlow;
466 }
467 }
468 58 }
469
470 //=============================================================================
471
472 58 void initLinks(int routingModel)
473 //
474 // Input: none
475 // Output: none
476 // Purpose: sets initial upstream/downstream conditions in links.
477 //
478 {
479 int i; // link index
480 int k; // conduit or pump index
481
482 // --- examine each link
483
2/2
✓ Branch 0 taken 10508 times.
✓ Branch 1 taken 58 times.
10566 for ( i = 0; i < Nobjects[LINK]; i++ )
484 {
485
2/2
✓ Branch 0 taken 21 times.
✓ Branch 1 taken 10487 times.
10508 if ( routingModel == SF) Link[i].newFlow = 0.0;
486
487 // --- otherwise if link is a conduit
488
2/2
✓ Branch 0 taken 9751 times.
✓ Branch 1 taken 736 times.
10487 else if ( Link[i].type == CONDUIT )
489 {
490 // --- assign initial flow to both ends of conduit
491 9751 k = Link[i].subIndex;
492 9751 Conduit[k].q1 = Link[i].newFlow / Conduit[k].barrels;
493 9751 Conduit[k].q2 = Conduit[k].q1;
494
495 // --- find areas based on initial flow depth
496 9751 Conduit[k].a1 = xsect_getAofY(&Link[i].xsect, Link[i].newDepth);
497 9751 Conduit[k].a2 = Conduit[k].a1;
498
499 // --- compute initial volume from area
500 {
501 9751 Link[i].newVolume = Conduit[k].a1 * link_getLength(i) *
502 9751 Conduit[k].barrels;
503 }
504 9751 Link[i].oldVolume = Link[i].newVolume;
505 }
506 }
507 58 }
508
509 //=============================================================================
510
511 12009458 double getLinkInflow(int j, double dt)
512 //
513 // Input: j = link index
514 // dt = routing time step (sec)
515 // Output: returns link inflow (cfs)
516 // Purpose: finds flow into upstream end of link at current time step under
517 // Steady or Kin. Wave routing.
518 //
519 {
520 12009458 int n1 = Link[j].node1;
521 double q;
522
2/2
✓ Branch 0 taken 10846 times.
✓ Branch 1 taken 11998612 times.
12009458 if ( Link[j].type == CONDUIT ||
523
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10846 times.
10846 Link[j].type == PUMP ||
524
0/2
✗ Branch 0 not taken.
✗ Branch 1 not taken.
12009458 Node[n1].type == STORAGE ) q = link_getInflow(j);
525 else q = 0.0;
526 12009458 return node_getMaxOutflow(n1, q, dt);
527 }
528
529 //=============================================================================
530
531 14040 void updateStorageState(int i, int j, int links[], double dt)
532 //
533 // Input: i = index of storage node
534 // j = current position in links array
535 // links = array of topo-sorted link indexes
536 // dt = routing time step (sec)
537 // Output: none
538 // Purpose: updates depth and volume of a storage node using successive
539 // approximation with under-relaxation for Steady or Kin. Wave
540 // routing.
541 //
542 {
543 int iter; // iteration counter
544 int stopped; // TRUE when iterations stop
545 double vFixed; // fixed terms of flow balance eqn.
546 double v2; // new volume estimate (ft3)
547 double d1; // initial value of storage depth (ft)
548 double d2; // updated value of storage depth (ft)
549
550 // --- see if storage node needs updating
551
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 14040 times.
14040 if ( Node[i].type != STORAGE ) return;
552
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 14040 times.
14040 if ( Node[i].updated ) return;
553
554 // --- compute terms of flow balance eqn.
555 // v2 = v1 + (inflow - outflow)*dt
556 // that do not depend on storage depth at end of time step
557 14040 vFixed = Node[i].oldVolume +
558 14040 0.5 * (Node[i].oldNetInflow + Node[i].inflow -
559 14040 Node[i].outflow) * dt;
560 14040 d1 = Node[i].newDepth;
561
562 // --- iterate finding outflow (which depends on depth) and subsequent
563 // new volume and depth until negligible depth change occurs
564 14040 iter = 1;
565 14040 stopped = FALSE;
566
3/4
✓ Branch 0 taken 28199 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 14159 times.
✓ Branch 3 taken 14040 times.
28199 while ( iter < MAXITER && !stopped )
567 {
568 // --- find new volume from flow balance eqn.
569 14159 v2 = vFixed - 0.5 * getStorageOutflow(i, j, links, dt) * dt;
570
571 // --- limit volume to full volume if no ponding
572 // and compute overflow rate
573
2/2
✓ Branch 0 taken 1230 times.
✓ Branch 1 taken 12929 times.
14159 v2 = MAX(0.0, v2);
574 14159 Node[i].overflow = 0.0;
575
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 14159 times.
14159 if ( v2 > Node[i].fullVolume )
576 {
577 Node[i].overflow = (v2 - MAX(Node[i].oldVolume,
578 Node[i].fullVolume)) / dt;
579 if ( Node[i].overflow < FUDGE ) Node[i].overflow = 0.0;
580 if ( !AllowPonding || Node[i].pondedArea == 0.0 )
581 v2 = Node[i].fullVolume;
582 }
583
584 // --- update node's volume & depth
585 14159 Node[i].newVolume = v2;
586 14159 d2 = node_getDepth(i, v2);
587 14159 Node[i].newDepth = d2;
588
589 // --- use under-relaxation to estimate new depth value
590 // and stop if close enough to previous value
591 14159 d2 = (1.0 - OMEGA)*d1 + OMEGA*d2;
592
2/2
✓ Branch 0 taken 14040 times.
✓ Branch 1 taken 119 times.
14159 if ( fabs(d2 - d1) <= STOPTOL ) stopped = TRUE;
593
594 // --- update old depth with new value and continue to iterate
595 14159 d1 = d2;
596 14159 iter++;
597 }
598
599 // --- mark node as being updated
600 14040 Node[i].updated = TRUE;
601 }
602
603 //=============================================================================
604
605 14159 double getStorageOutflow(int i, int j, int links[], double dt)
606 //
607 // Input: i = index of storage node
608 // j = current position in links array
609 // links = array of topo-sorted link indexes
610 // dt = routing time step (sec)
611 // Output: returns total outflow from storage node (cfs)
612 // Purpose: computes total flow released from a storage node.
613 //
614 {
615 int k, m;
616 14159 double outflow = 0.0;
617
618
2/2
✓ Branch 0 taken 25745 times.
✓ Branch 1 taken 2573 times.
28318 for (k = j; k < Nobjects[LINK]; k++)
619 {
620 25745 m = links[k];
621
2/2
✓ Branch 0 taken 11586 times.
✓ Branch 1 taken 14159 times.
25745 if ( Link[m].node1 != i ) break;
622 14159 outflow += getLinkInflow(m, dt);
623 }
624 14159 return outflow;
625 }
626
627 //=============================================================================
628
629 12629258 void setNewNodeState(int j, double dt)
630 //
631 // Input: j = node index
632 // dt = time step (sec)
633 // Output: none
634 // Purpose: updates state of node after current time step
635 // for Steady Flow or Kinematic Wave flow routing.
636 //
637 {
638 int canPond; // TRUE if ponding can occur at node
639 double newNetInflow; // inflow - outflow at node (cfs)
640
641 // --- update terminal storage nodes
642
2/2
✓ Branch 0 taken 14040 times.
✓ Branch 1 taken 12615218 times.
12629258 if ( Node[j].type == STORAGE )
643 {
644
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 14040 times.
14040 if ( Node[j].updated == FALSE )
645 updateStorageState(j, Nobjects[LINK], NULL, dt);
646 14040 return;
647 }
648
649 // --- update stored volume
650 12615218 newNetInflow = Node[j].inflow - Node[j].outflow - Node[j].losses;
651 12615218 Node[j].newVolume = Node[j].oldVolume + newNetInflow * dt;
652
2/2
✓ Branch 0 taken 12549053 times.
✓ Branch 1 taken 66165 times.
12615218 if ( Node[j].newVolume < FUDGE ) Node[j].newVolume = 0.0;
653
654 // --- determine any overflow lost from system
655 12615218 Node[j].overflow = 0.0;
656
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 12615218 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
12615218 canPond = (AllowPonding && Node[j].pondedArea > 0.0);
657
2/2
✓ Branch 0 taken 66165 times.
✓ Branch 1 taken 12549053 times.
12615218 if ( Node[j].newVolume > Node[j].fullVolume )
658 {
659
1/2
✓ Branch 0 taken 66165 times.
✗ Branch 1 not taken.
66165 Node[j].overflow = (Node[j].newVolume - MAX(Node[j].oldVolume,
660 66165 Node[j].fullVolume)) / dt;
661
2/2
✓ Branch 0 taken 2190 times.
✓ Branch 1 taken 63975 times.
66165 if ( Node[j].overflow < FUDGE ) Node[j].overflow = 0.0;
662
1/2
✓ Branch 0 taken 66165 times.
✗ Branch 1 not taken.
66165 if ( !canPond ) Node[j].newVolume = Node[j].fullVolume;
663 }
664
665 // --- compute a depth from volume
666 // (depths at upstream nodes are subsequently adjusted in
667 // setNewLinkState to reflect depths in connected conduit)
668 12615218 Node[j].newDepth = node_getDepth(j, Node[j].newVolume);
669 }
670
671 //=============================================================================
672
673 11995299 void setNewLinkState(int j)
674 //
675 // Input: j = link index
676 // Output: none
677 // Purpose: updates state of link after current time step under
678 // Steady Flow or Kinematic Wave flow routing
679 //
680 {
681 int k;
682 double a, y1, y2;
683
684 11995299 Link[j].newDepth = 0.0;
685 11995299 Link[j].newVolume = 0.0;
686
687
2/2
✓ Branch 0 taken 11989899 times.
✓ Branch 1 taken 5400 times.
11995299 if ( Link[j].type == CONDUIT )
688 {
689 // --- find avg. depth from entry/exit conditions
690 11989899 k = Link[j].subIndex;
691 11989899 a = 0.5 * (Conduit[k].a1 + Conduit[k].a2);
692 11989899 Link[j].newVolume = a * link_getLength(j) * Conduit[k].barrels;
693 11989899 y1 = xsect_getYofA(&Link[j].xsect, Conduit[k].a1);
694 11989899 y2 = xsect_getYofA(&Link[j].xsect, Conduit[k].a2);
695 11989899 Link[j].newDepth = 0.5 * (y1 + y2);
696
697 // --- update depths at end nodes
698 11989899 updateNodeDepth(Link[j].node1, y1 + Link[j].offset1);
699 11989899 updateNodeDepth(Link[j].node2, y2 + Link[j].offset2);
700
701 // --- check if capacity limited
702
2/2
✓ Branch 0 taken 2317 times.
✓ Branch 1 taken 11987582 times.
11989899 if ( Conduit[k].a1 >= Link[j].xsect.aFull )
703 {
704 2317 Conduit[k].capacityLimited = TRUE;
705 2317 Conduit[k].fullState = ALL_FULL;
706 }
707 else
708 {
709 11987582 Conduit[k].capacityLimited = FALSE;
710 11987582 Conduit[k].fullState = 0;
711 }
712 }
713 11995299 }
714
715 //=============================================================================
716
717 23979798 void updateNodeDepth(int i, double y)
718 //
719 // Input: i = node index
720 // y = flow depth (ft)
721 // Output: none
722 // Purpose: updates water depth at a node with a possibly higher value.
723 //
724 {
725 // --- storage nodes were updated elsewhere
726
2/2
✓ Branch 0 taken 8640 times.
✓ Branch 1 taken 23971158 times.
23979798 if ( Node[i].type == STORAGE ) return;
727
728 // --- if non-outfall node is flooded, then use full depth
729
3/4
✓ Branch 0 taken 23341639 times.
✓ Branch 1 taken 629519 times.
✓ Branch 2 taken 23341639 times.
✗ Branch 3 not taken.
23971158 if ( Node[i].type != OUTFALL && Node[i].degree > 0 &&
730
2/2
✓ Branch 0 taken 4634 times.
✓ Branch 1 taken 23337005 times.
23341639 Node[i].overflow > 0.0 ) y = Node[i].fullDepth;
731
732 // --- if current new depth below y
733
2/2
✓ Branch 0 taken 728982 times.
✓ Branch 1 taken 23242176 times.
23971158 if ( Node[i].newDepth < y )
734 {
735 // --- update new depth
736 728982 Node[i].newDepth = y;
737
738 // --- depth cannot exceed full depth (if value exists)
739
2/4
✓ Branch 0 taken 728982 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 728982 times.
728982 if ( Node[i].fullDepth > 0.0 && y > Node[i].fullDepth )
740 {
741 Node[i].newDepth = Node[i].fullDepth;
742 }
743 }
744 }
745
746 //=============================================================================
747
748 11894379 int steadyflow_execute(int j, double* qin, double* qout, double tStep)
749 //
750 // Input: j = link index
751 // qin = inflow to link (cfs)
752 // tStep = time step (sec)
753 // Output: qin = adjusted inflow to link (limited by flow capacity) (cfs)
754 // qout = link's outflow (cfs)
755 // returns 1 if successful
756 // Purpose: performs steady flow routing through a single link.
757 //
758 {
759 int k;
760 double s;
761 double q;
762
763 // --- use Manning eqn. to compute flow area for conduits
764
1/2
✓ Branch 0 taken 11894379 times.
✗ Branch 1 not taken.
11894379 if ( Link[j].type == CONDUIT )
765 {
766 11894379 k = Link[j].subIndex;
767 11894379 q = (*qin) / Conduit[k].barrels;
768
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11894379 times.
11894379 if ( Link[j].xsect.type == DUMMY ) Conduit[k].a1 = 0.0;
769 else
770 {
771 // --- adjust flow for evap and infil losses
772 11894379 q -= link_getLossRate(j, SF, q, tStep);
773
774 // --- flow can't exceed full flow
775
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11894379 times.
11894379 if ( q > Link[j].qFull )
776 {
777 q = Link[j].qFull;
778 Conduit[k].a1 = Link[j].xsect.aFull;
779 (*qin) = q * Conduit[k].barrels;
780 }
781
782 // --- infer flow area from flow rate
783 else
784 {
785 11894379 s = q / Conduit[k].beta;
786 11894379 Conduit[k].a1 = xsect_getAofS(&Link[j].xsect, s);
787 }
788 }
789 11894379 Conduit[k].a2 = Conduit[k].a1;
790
791 11894379 Conduit[k].q1Old = Conduit[k].q1;
792 11894379 Conduit[k].q2Old = Conduit[k].q2;
793
794 11894379 Conduit[k].q1 = q;
795 11894379 Conduit[k].q2 = q;
796 11894379 (*qout) = q * Conduit[k].barrels;
797 }
798 else (*qout) = (*qin);
799 11894379 return 1;
800 }
801
802 //=============================================================================
803