GCC Code Coverage Report

Directory: src/solver/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 92.4% 549 / 0 / 594
Functions: 100.0% 30 / 0 / 30
Branches: 71.1% 249 / 0 / 350
node.c
Line Branch Exec Source
1 //-----------------------------------------------------------------------------
2 // node.c
3 //
4 // Project: EPA SWMM5
5 // Version: 5.2
6 // Date: 10/29/22 (Build 5.2.2)
7 // Author: L. Rossman
8 //
9 // Conveyance system node functions.
10 //
11 // Update History
12 // ==============
13 // Build 5.1.007:
14 // - Ponded area property for storage nodes deprecated.
15 // - Support for Green-Ampt seepage from bottom and sides of storage node added.
16 // - Storage node evap. & seepage losses now computed at start of each routing
17 // time step.
18 // Build 5.1.008:
19 // - Support added for sending outfall discharge to a subcatchment.
20 // Build 5.1.010:
21 // - Storage losses now based on node's new volume instead of old volume.
22 // Build 5.1.013:
23 // - A surcharge depth can now be applied to storage nodes.
24 // - A negative inflow is now assigned to an Outfall node with backflow.
25 // Build 5.1.014:
26 // - Fixed bug in storage_losses() that affected storage exfiltration.
27 // Build 5.1.015:
28 // - Fatal error issued if a storage node's area curve produces a negative
29 // volume when extrapolated to the node's full depth.
30 // Build 5.2.0:
31 // - Support added Streets and Inlets.
32 // - Support added for analytical storage shapes.
33 // Build 5.2.1:
34 // - Warning no longer issued when node full depth is increased to match
35 // crown of highest connecting link.
36 // - a2 term for paraboloid shaped storage units was corrected
37 // Build 5.2.2:
38 // - Warning restored for node full depth being increased to crown of highest
39 // connecting link.
40 //-----------------------------------------------------------------------------
41 #define _CRT_SECURE_NO_DEPRECATE
42
43 #include <stdlib.h>
44 #include <string.h>
45 #include <math.h>
46 #include "headers.h"
47 #include "findroot.h"
48
49 //-----------------------------------------------------------------------------
50 // Local Declarations
51 //-----------------------------------------------------------------------------
52 typedef struct
53 {
54 int k; // storage unit index
55 double v; // storage unit volume (ft3)
56 } TStorageVol;
57
58 //-----------------------------------------------------------------------------
59 // External functions (declared in funcs.h)
60 //-----------------------------------------------------------------------------
61 // node_readParams (called from readNode in input.c)
62 // node_validate (called from project_validate)
63 // node_initState (called from project_init)
64 // node_setOldHydState (called from routing_execute)
65 // node_setOldQualState (called from routing_execute)
66 // node_initFlows (called from routing_execute)
67 // node_setOutletDepth (called from routing_execute)
68 // node_getLosses (called from routing_execute)
69 // node_getSystemOutflow (called from removeOutflows in routing.c)
70 // node_getResults (called from output_saveNodeResults)
71 // node_getPondedArea (called from initNodeStates in dynwave.c)
72 // node_getOutflow (called from link_getInflow & conduit_getInflow)
73 // node_getMaxOutflow (called from flowrout.c and dynwave.c)
74 // node_getSurfArea
75 // node_getDepth
76 // node_getVolume
77
78 //-----------------------------------------------------------------------------
79 // Local functions
80 //-----------------------------------------------------------------------------
81 static void node_setParams(int j, int type, int k, double x[]);
82 static int junc_readParams(int j, int k, char* tok[], int ntoks);
83
84 static int outfall_readParams(int j, int k, char* tok[], int ntoks);
85 static void outfall_setOutletDepth(int j, double yNorm, double yCrit, double z);
86
87 static int storage_readParams(int j, int k, char* tok[], int ntoks);
88 static double storage_getDepth(int j, double v);
89 static double storage_getVolume(int j, double d);
90 static double storage_getSurfArea(int j, double d);
91 static void storage_getVolDiff(double y, double* f, double* df, void* p);
92 static double storage_getOutflow(int j, int i);
93 static double storage_getLosses(int j, double tStep);
94
95 static int divider_readParams(int j, int k, char* tok[], int ntoks);
96 static void divider_validate(int j);
97 static double divider_getOutflow(int j, int link);
98
99
100 //=============================================================================
101
102 10183 int node_readParams(int j, int type, int k, char* tok[], int ntoks)
103 //
104 // Input: j = node index
105 // type = node type code
106 // k = index of node type
107 // tok[] = array of string tokens
108 // ntoks = number of tokens
109 // Output: returns an error code
110 // Purpose: reads node properties from a tokenized line of input.
111 //
112 {
113
4/5
✓ Branch 0 taken 9714 times.
✓ Branch 1 taken 348 times.
✓ Branch 2 taken 117 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 10183 switch ( type )
114 {
115 9714 case JUNCTION: return junc_readParams(j, k, tok, ntoks);
116 348 case OUTFALL: return outfall_readParams(j, k, tok, ntoks);
117 117 case STORAGE: return storage_readParams(j, k, tok, ntoks);
118 4 case DIVIDER: return divider_readParams(j, k, tok, ntoks);
119 ✗ default: return 0;
120 }
121 }
122
123 //=============================================================================
124
125 10183 void node_setParams(int j, int type, int k, double x[])
126 //
127 // Input: j = node index
128 // type = node type code
129 // k = index of node type
130 // x[] = array of property values
131 // Output: none
132 // Purpose: assigns property values to a node.
133 //
134 {
135 10183 Node[j].type = type;
136 10183 Node[j].subIndex = k;
137 10183 Node[j].invertElev = x[0] / UCF(LENGTH);
138 10183 Node[j].crownElev = Node[j].invertElev;
139 10183 Node[j].initDepth = 0.0;
140 10183 Node[j].newVolume = 0.0;
141 10183 Node[j].fullVolume = 0.0;
142 10183 Node[j].fullDepth = 0.0;
143 10183 Node[j].surDepth = 0.0;
144 10183 Node[j].pondedArea = 0.0;
145 10183 Node[j].degree = 0;
146 10183 Node[j].inlet = NO_INLET;
147
4/5
✓ Branch 0 taken 9714 times.
✓ Branch 1 taken 348 times.
✓ Branch 2 taken 117 times.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
 10183 switch (type)
148 {
149 9714 case JUNCTION:
150 9714 Node[j].fullDepth = x[1] / UCF(LENGTH);
151 9714 Node[j].initDepth = x[2] / UCF(LENGTH);
152 9714 Node[j].surDepth = x[3] / UCF(LENGTH);
153 9714 Node[j].pondedArea = x[4] / (UCF(LENGTH)*UCF(LENGTH));
154 9714 break;
155
156 348 case OUTFALL:
157 348 Outfall[k].type = (int)x[1];
158 348 Outfall[k].fixedStage = x[2] / UCF(LENGTH);
159 348 Outfall[k].tideCurve = (int)x[3];
160 348 Outfall[k].stageSeries = (int)x[4];
161 348 Outfall[k].hasFlapGate = (char)x[5];
162 348 Outfall[k].routeTo = (int)x[6];
163 348 Outfall[k].wRouted = NULL;
164
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 347 times.
 348 if ( Outfall[k].routeTo >= 0 )
165 {
166 1 Outfall[k].wRouted =
167 1 (double *) calloc(Nobjects[POLLUT], sizeof(double));
168 }
169 348 break;
170
171 117 case STORAGE:
172 117 Node[j].fullDepth = x[1] / UCF(LENGTH);
173 117 Node[j].initDepth = x[2] / UCF(LENGTH);
174 117 Storage[k].shape = (int)x[3];
175 117 Storage[k].a1 = x[4];
176 117 Storage[k].a2 = x[5];
177 117 Storage[k].a0 = x[6];
178 117 Storage[k].aCurve = (int)x[7];
179 117 Node[j].surDepth = x[8] / UCF(LENGTH);
180 117 Storage[k].fEvap = x[9];
181 117 break;
182
183 4 case DIVIDER:
184 4 Divider[k].link = (int)x[1];
185 4 Divider[k].type = (int)x[2];
186 4 Divider[k].flowCurve = (int)x[3];
187 4 Divider[k].qMin = x[4] / UCF(FLOW);
188 4 Divider[k].dhMax = x[5];
189 4 Divider[k].cWeir = x[6];
190 4 Node[j].fullDepth = x[7] / UCF(LENGTH);
191 4 Node[j].initDepth = x[8] / UCF(LENGTH);
192 4 Node[j].surDepth = x[9] / UCF(LENGTH);
193 4 Node[j].pondedArea = x[10] / (UCF(LENGTH)*UCF(LENGTH));
194 4 break;
195 }
196 10183 }
197
198 //=============================================================================
199
200 10183 void node_validate(int j)
201 //
202 // Input: j = node index
203 // Output: none
204 // Purpose: validates a node's properties.
205 //
206 {
207 TDwfInflow* inflow;
208
209 // --- see if full depth was increased to accommodate conduit crown
210
4/4
✓ Branch 0 taken 434 times.
✓ Branch 1 taken 9749 times.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 433 times.
 10183 if ( Node[j].fullDepth > Node[j].oldDepth && Node[j].oldDepth > 0.0 )
211 {
212 1 report_writeWarningMsg(WARN02, Node[j].ID);
213 }
214
215 // --- check that initial depth does not exceed max. depth
216
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10183 times.
 10183 if ( Node[j].initDepth > Node[j].fullDepth + Node[j].surDepth )
217 ✗ report_writeErrorMsg(ERR_NODE_DEPTH, Node[j].ID);
218
219 // --- check for negative volume for storage node at full depth
220
2/2
✓ Branch 0 taken 117 times.
✓ Branch 1 taken 10066 times.
 10183 if (Node[j].type == STORAGE)
221
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 117 times.
 117 if (node_getVolume(j, Node[j].fullDepth) < 0.0)
222 ✗ report_writeErrorMsg(ERR_STORAGE_VOLUME, Node[j].ID);
223
224
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 10179 times.
 10183 if ( Node[j].type == DIVIDER ) divider_validate(j);
225
226 // --- initialize dry weather inflows
227 10183 inflow = Node[j].dwfInflow;
228
2/2
✓ Branch 0 taken 3943 times.
✓ Branch 1 taken 10183 times.
 14126 while (inflow)
229 {
230 3943 inflow_initDwfInflow(inflow);
231 3943 inflow = inflow->next;
232 }
233 10183 }
234
235 //=============================================================================
236
237 10183 void node_initState(int j)
238 //
239 // Input: j = node index
240 // Output: none
241 // Purpose: initializes a node's state variables at start of simulation.
242 //
243 {
244 int p, k;
245
246 // --- initialize depth
247 10183 Node[j].oldDepth = Node[j].initDepth;
248 10183 Node[j].newDepth = Node[j].oldDepth;
249 10183 Node[j].crownElev = Node[j].invertElev;
250
251 10183 Node[j].fullVolume = node_getVolume(j, Node[j].fullDepth);
252 10183 Node[j].oldVolume = node_getVolume(j, Node[j].oldDepth);
253 10183 Node[j].newVolume = Node[j].oldVolume;
254
255 // --- initialize water quality state
256
2/2
✓ Branch 0 taken 25771 times.
✓ Branch 1 taken 10183 times.
 35954 for (p = 0; p < Nobjects[POLLUT]; p++)
257 {
258 25771 Node[j].oldQual[p] = 0.0;
259 25771 Node[j].newQual[p] = 0.0;
260 }
261
262 // --- initialize any inflow
263 10183 Node[j].oldLatFlow = 0.0;
264 10183 Node[j].newLatFlow = 0.0;
265 10183 Node[j].apiExtInflow = 0.0;
266 10183 Node[j].losses = 0.0;
267
268 // --- initialize storage nodes
269
2/2
✓ Branch 0 taken 117 times.
✓ Branch 1 taken 10066 times.
 10183 if ( Node[j].type == STORAGE )
270 {
271 // --- set hydraulic residence time to 0
272 117 k = Node[j].subIndex;
273 117 Storage[k].hrt = 0.0;
274
275 // --- initialize exfiltration properties
276
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 114 times.
 117 if ( Storage[k].exfil ) exfil_initState(k);
277 }
278
279 // --- initialize flow stream routed from outfall onto a subcatchment
280
2/2
✓ Branch 0 taken 348 times.
✓ Branch 1 taken 9835 times.
 10183 if ( Node[j].type == OUTFALL )
281 {
282 348 k = Node[j].subIndex;
283
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 347 times.
 348 if ( Outfall[k].routeTo >= 0 )
284 {
285 1 Outfall[k].vRouted = 0.0;
286
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 for (p = 0; p < Nobjects[POLLUT]; p++) Outfall[k].wRouted[p] = 0.0;
287 }
288 }
289 10183 }
290
291 //=============================================================================
292
293 37894467 void node_setOldHydState(int j)
294 //
295 // Input: j = node index
296 // Output: none
297 // Purpose: replaces a node's old hydraulic state values with new ones.
298 //
299 {
300 37894467 Node[j].oldDepth = Node[j].newDepth;
301 37894467 Node[j].oldVolume = Node[j].newVolume;
302 37894467 Node[j].oldFlowInflow = Node[j].inflow;
303 37894467 Node[j].oldNetInflow = Node[j].inflow - Node[j].outflow;
304 37894467 }
305
306 //=============================================================================
307
308 20940478 void node_setOldQualState(int j)
309 //
310 // Input: j = node index
311 // Output: none
312 // Purpose: replaces a node's old water quality state values with new ones.
313 //
314 {
315 int p;
316
2/2
✓ Branch 0 taken 78198686 times.
✓ Branch 1 taken 20940478 times.
 99139164 for (p = 0; p < Nobjects[POLLUT]; p++)
317 {
318 78198686 Node[j].oldQual[p] = Node[j].newQual[p];
319 78198686 Node[j].newQual[p] = 0.0;
320 }
321 20940478 }
322
323 //=============================================================================
324
325 37894467 void node_initFlows(int j, double tStep)
326 //
327 // Input: j = node index
328 // tStep = time step (sec)
329 // Output: none
330 // Purpose: initializes a node's inflow/outflow/overflow at start of time step.
331 //
332 {
333 // --- initialize inflow & outflow
334 37894467 Node[j].inflow = Node[j].newLatFlow;
335 37894467 Node[j].outflow = Node[j].losses;
336
337 // --- set overflow to any excess stored volume
338
2/2
✓ Branch 0 taken 88273 times.
✓ Branch 1 taken 37806194 times.
 37894467 if ( Node[j].newVolume > Node[j].fullVolume )
339 88273 Node[j].overflow = (Node[j].newVolume - Node[j].fullVolume) / tStep;
340 37806194 else Node[j].overflow = 0.0;
341 37894467 }
342
343 //=============================================================================
344
345 12629377 double node_getDepth(int j, double v)
346 //
347 // Input: j = node index
348 // v = volume (ft3)
349 // Output: returns depth of water at a node (ft)
350 // Purpose: computes a node's water depth from its volume.
351 //
352 {
353
2/2
✓ Branch 0 taken 14159 times.
✓ Branch 1 taken 12615218 times.
 12629377 switch ( Node[j].type )
354 {
355 14159 case STORAGE: return storage_getDepth(j, v);
356 12615218 default: return 0.0;
357 }
358 }
359
360 //=============================================================================
361
362 65116604 double node_getVolume(int j, double d)
363 //
364 // Input: j = node index
365 // d = water depth (ft)
366 // Output: returns volume of water at a node (ft3)
367 // Purpose: computes volume stored at a node from its water depth.
368 //
369 {
370
2/2
✓ Branch 0 taken 485671 times.
✓ Branch 1 taken 64630933 times.
 65116604 switch ( Node[j].type )
371 {
372 485671 case STORAGE: return storage_getVolume(j, d);
373
374 64630933 default:
375
2/2
✓ Branch 0 taken 64630825 times.
✓ Branch 1 taken 108 times.
 64630933 if ( Node[j].fullDepth > 0.0 )
376 64630825 return Node[j].fullVolume * (d / Node[j].fullDepth);
377 108 else return 0.0;
378 }
379 }
380
381 //=============================================================================
382
383 67143458 double node_getSurfArea(int j, double d)
384 //
385 // Input: j = node index
386 // d = water depth (ft)
387 // Output: returns surface area of water at a node (ft2)
388 // Purpose: computes surface area of water stored at a node from water depth.
389 //
390 {
391
2/2
✓ Branch 0 taken 495707 times.
✓ Branch 1 taken 66647751 times.
 67143458 switch (Node[j].type)
392 {
393 495707 case STORAGE: return storage_getSurfArea(j, d);
394 66647751 default: return 0.0;
395 }
396 }
397
398 //=============================================================================
399
400 11998612 double node_getOutflow(int j, int k)
401 //
402 // Input: j = node index
403 // k = link index
404 // Output: returns flow rate (cfs)
405 // Purpose: computes outflow from node available for inflow into a link.
406 //
407 {
408
3/3
✓ Branch 0 taken 1920 times.
✓ Branch 1 taken 17353 times.
✓ Branch 2 taken 11979339 times.
 11998612 switch ( Node[j].type )
409 {
410 1920 case DIVIDER: return divider_getOutflow(j, k);
411 17353 case STORAGE: return storage_getOutflow(j, k);
412 11979339 default: return Node[j].inflow + Node[j].overflow;
413 }
414 }
415
416 //=============================================================================
417
418 12035293 double node_getMaxOutflow(int j, double q, double tStep)
419 //
420 // Input: j = node index
421 // q = original outflow rate (cfs)
422 // tStep = time step (sec)
423 // Output: returns modified flow rate (cfs)
424 // Purpose: limits outflow rate from a node with storage volume.
425 //
426 {
427 double qMax;
428
2/2
✓ Branch 0 taken 54034 times.
✓ Branch 1 taken 11981259 times.
 12035293 if ( Node[j].fullVolume > 0.0 )
429 {
430 54034 qMax = Node[j].inflow + Node[j].oldVolume / tStep;
431
2/2
✓ Branch 0 taken 2731 times.
✓ Branch 1 taken 51303 times.
 54034 if ( q > qMax ) q = qMax;
432 }
433
2/2
✓ Branch 0 taken 11315709 times.
✓ Branch 1 taken 719584 times.
 12035293 return MAX(0.0, q);
434 }
435
436 //=============================================================================
437
438 37894467 double node_getSystemOutflow(int j, int *isFlooded)
439 //
440 // Input: j = node index
441 // isFlooded = TRUE if node becomes flooded
442 // Output: returns flow rate lost from system (cfs)
443 // Purpose: computes flow rate at outfalls and flooded nodes.
444 //
445 {
446 37894467 double outflow = 0.0;;
447
448 // --- assume there is no flooding
449 37894467 *isFlooded = FALSE;
450
451 // --- if node is an outfall
452
2/2
✓ Branch 0 taken 1309094 times.
✓ Branch 1 taken 36585373 times.
 37894467 if ( Node[j].type == OUTFALL )
453 {
454 // --- node receives inflow from outfall conduit
455
2/2
✓ Branch 0 taken 1308539 times.
✓ Branch 1 taken 555 times.
 1309094 if ( Node[j].outflow == 0.0 ) outflow = Node[j].inflow;
456
457 // --- node sends flow into outfall conduit
458 // (therefore it has a negative outflow)
459 else
460 {
461
1/2
✓ Branch 0 taken 555 times.
✗ Branch 1 not taken.
 555 if ( Node[j].inflow == 0.0 )
462 {
463 555 outflow = -Node[j].outflow;
464 555 Node[j].inflow = fabs(outflow);
465 }
466 }
467
468 // --- set overflow and volume to 0
469 1309094 Node[j].overflow = 0.0;
470 1309094 Node[j].newVolume = 0.0;
471 }
472
473 // --- node is a terminal node under Steady or Kin. Wave routing
474
2/2
✓ Branch 0 taken 11994339 times.
✓ Branch 1 taken 24591034 times.
 36585373 else if ( RouteModel != DW &&
475
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 11994339 times.
 11994339 Node[j].degree == 0 &&
476 ✗ Node[j].type != STORAGE
477 )
478 {
479 ✗ if ( Node[j].outflow == 0.0 ) outflow = Node[j].inflow;
480 ✗ Node[j].overflow = 0.0;
481 ✗ Node[j].newVolume = 0.0;
482 }
483
484 // --- otherwise node is an interior node and any
485 // overflow is considered as system outflow and flooding
486 else
487 {
488
2/2
✓ Branch 0 taken 36497100 times.
✓ Branch 1 taken 88273 times.
 36585373 if ( Node[j].newVolume <= Node[j].fullVolume)
489 36497100 outflow = Node[j].overflow;
490
2/2
✓ Branch 0 taken 37764 times.
✓ Branch 1 taken 36547609 times.
 36585373 if ( outflow > 0.0 ) *isFlooded = TRUE;
491 }
492 37894467 return outflow;
493 }
494
495 //=============================================================================
496
497 1931254 void node_getResults(int j, double f, float x[])
498 //
499 // Input: j = node index
500 // f = weighting factor
501 // x[] = array of nodal reporting variables
502 // Output: none
503 // Purpose: computes weighted average of old and new results at a node.
504 //
505 {
506 int p;
507 double z;
508 1931254 double f1 = 1.0 - f;
509
510 1931254 z = (f1 * Node[j].oldDepth + f * Node[j].newDepth) * UCF(LENGTH);
511 1931254 x[NODE_DEPTH] = (float)z;
512 1931254 z = Node[j].invertElev * UCF(LENGTH);
513 1931254 x[NODE_HEAD] = x[NODE_DEPTH] + (float)z;
514 1931254 z = (f1*Node[j].oldVolume + f*Node[j].newVolume) * UCF(VOLUME);
515 1931254 x[NODE_VOLUME] = (float)z;
516 1931254 z = (f1*Node[j].oldLatFlow + f*Node[j].newLatFlow) * UCF(FLOW);
517 1931254 x[NODE_LATFLOW] = (float)z;
518 1931254 z = (f1*Node[j].oldFlowInflow + f*Node[j].inflow) * UCF(FLOW);
519 1931254 x[NODE_INFLOW] = (float)z;
520 1931254 z = Node[j].overflow * UCF(FLOW);
521 1931254 x[NODE_OVERFLOW] = (float)z;
522
523
4/4
✓ Branch 0 taken 1927798 times.
✓ Branch 1 taken 3456 times.
✓ Branch 2 taken 827538 times.
✓ Branch 3 taken 1927798 times.
 2758792 if ( !IgnoreQuality ) for (p = 0; p < Nobjects[POLLUT]; p++)
524 {
525 827538 z = f1*Node[j].oldQual[p] + f*Node[j].newQual[p];
526 827538 x[NODE_QUAL+p] = (float)z;
527 }
528 1931254 }
529
530 //=============================================================================
531
532 1979248 void node_setOutletDepth(int j, double yNorm, double yCrit, double z)
533 //
534 // Input: j = node index
535 // yNorm = normal flow depth (ft)
536 // yCrit = critical flow depth (ft)
537 // z = offset of connecting outfall link from node invert (ft)
538 // Output: none
539 // Purpose: sets water depth at a node that serves as an outlet point.
540 //
541 {
542
1/3
✗ Branch 0 not taken.
✓ Branch 1 taken 1979248 times.
✗ Branch 2 not taken.
 1979248 switch (Node[j].type)
543 {
544 // --- do nothing if outlet is a storage unit
545 ✗ case STORAGE:
546 ✗ return;
547
548 // --- if outlet is a designated outfall then use outfall's specs
549 1979248 case OUTFALL:
550 1979248 outfall_setOutletDepth(j, yNorm, yCrit, z);
551 1979248 break;
552
553 // --- for all other nodes, use min. of critical & normal depths
554 ✗ default:
555 ✗ if ( z > 0.0 ) Node[j].newDepth = 0.0;
556 ✗ else Node[j].newDepth = MIN(yNorm, yCrit);
557 }
558 }
559
560 //=============================================================================
561
562 2800436 double node_getPondedArea(int j, double d)
563 //
564 // Input: j = node index
565 // d = water depth (ft)
566 // Output: returns surface area of water at a node (ft2)
567 // Purpose: computes surface area of water at a node based on depth.
568 //
569 {
570 double a;
571
572 // --- use regular getSurfArea function if node not flooded
573
4/4
✓ Branch 0 taken 188112 times.
✓ Branch 1 taken 2612324 times.
✓ Branch 2 taken 10929 times.
✓ Branch 3 taken 177183 times.
 2800436 if ( d <= Node[j].fullDepth || Node[j].pondedArea == 0.0 )
574 {
575 2623253 return node_getSurfArea(j, d);
576 }
577
578 // --- compute ponded depth
579 177183 d = d - Node[j].fullDepth;
580
581 // --- use ponded area for flooded node
582 177183 a = Node[j].pondedArea;
583
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 177183 times.
 177183 if ( a <= 0.0 ) a = node_getSurfArea(j, Node[j].fullDepth);
584 177183 return a;
585 }
586
587 //=============================================================================
588
589 37894467 double node_getLosses(int j, double tStep)
590 //
591 // Input: j = node index
592 // tStep = time step (sec)
593 // Output: returns water loss rate at node (ft3)
594 // Purpose: computes the rates of evaporation and infiltration over a given
595 // time step for a node.
596 //
597 {
598
2/2
✓ Branch 0 taken 143130 times.
✓ Branch 1 taken 37751337 times.
 37894467 if ( Node[j].type == STORAGE ) return storage_getLosses(j, tStep);
599 37751337 else return 0.0;
600 }
601
602 //=============================================================================
603 // J U N C T I O N M E T H O D S
604 //=============================================================================
605
606 9714 int junc_readParams(int j, int k, char* tok[], int ntoks)
607 //
608 // Input: j = node index
609 // k = junction index
610 // tok[] = array of string tokens
611 // ntoks = number of tokens
612 // Output: returns an error message
613 // Purpose: reads a junction's properties from a tokenized line of input.
614 //
615 // Format of input line is:
616 // nodeID elev maxDepth initDepth surDepth aPond
617 {
618 int i;
619 double x[6];
620 char* id;
621
622
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9714 times.
 9714 if ( ntoks < 2 ) return error_setInpError(ERR_ITEMS, "");
623 9714 id = project_findID(NODE, tok[0]);
624
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9714 times.
 9714 if ( id == NULL ) return error_setInpError(ERR_NAME, tok[0]);
625
626 // --- parse invert elev., max. depth, init. depth, surcharged depth,
627 // & ponded area values
628
2/2
✓ Branch 0 taken 48570 times.
✓ Branch 1 taken 9714 times.
 58284 for ( i = 1; i <= 5; i++ )
629 {
630 48570 x[i-1] = 0.0;
631
1/2
✓ Branch 0 taken 48570 times.
✗ Branch 1 not taken.
 48570 if ( i < ntoks )
632 {
633
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 48570 times.
 48570 if ( ! getDouble(tok[i], &x[i-1]) )
634 ✗ return error_setInpError(ERR_NUMBER, tok[i]);
635 }
636 }
637
638 // --- check for non-negative values (except for invert elev.)
639
2/2
✓ Branch 0 taken 38856 times.
✓ Branch 1 taken 9714 times.
 48570 for ( i = 1; i <= 4; i++ )
640 {
641
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 38856 times.
 38856 if ( x[i] < 0.0 ) return error_setInpError(ERR_NUMBER, tok[i+1]);
642 }
643
644 // --- add parameters to junction object
645 9714 Node[j].ID = id;
646 9714 node_setParams(j, JUNCTION, k, x);
647 9714 return 0;
648 }
649
650 //=============================================================================
651 // S T O R A G E M E T H O D S
652 //=============================================================================
653
654 117 int storage_readParams(int j, int k, char* tok[], int ntoks)
655 //
656 // Input: j = node index
657 // k = storage unit index
658 // tok[] = array of string tokens
659 // ntoks = number of tokens
660 // Output: returns an error message
661 // Purpose: reads a storage unit's properties from a tokenized line of input.
662 //
663 // Format of input line is:
664 // nodeID elev maxDepth initDepth curveType a1 a2 a0 surDepth fEvap (infil)
665 // nodeID elev maxDepth initDepth TABULAR curveID surDepth fEvap (infil)
666 // x[0] x[1] x[2] x[3] x[4..7] x[8] x[9]
667 {
668 int i, m, n;
669 double x[10], y[3];
670 double A, B; //base semi-axis length & width for conical shape
671 double L, W; //base length & width for pyramidal shape
672 double Z; //run over rise for conical & pyramidal sides
673 char* id;
674
675 // --- get ID name
676
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 117 times.
 117 if ( ntoks < 6 ) return error_setInpError(ERR_ITEMS, "");
677 117 id = project_findID(NODE, tok[0]);
678
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 117 times.
 117 if ( id == NULL ) return error_setInpError(ERR_NAME, tok[0]);
679
680 // --- get invert elev, max. depth, & init. depth
681
2/2
✓ Branch 0 taken 351 times.
✓ Branch 1 taken 117 times.
 468 for ( i = 1; i <= 3; i++ )
682 {
683
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 351 times.
 351 if ( ! getDouble(tok[i], &x[i-1]) )
684 ✗ return error_setInpError(ERR_NUMBER, tok[i]);
685 }
686
687 // --- get surf. area relation type
688 117 m = findmatch(tok[4], RelationWords);
689
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 117 times.
 117 if ( m < 0 ) return error_setInpError(ERR_KEYWORD, tok[4]);
690 117 x[3] = m;
691 117 x[4] = 0.0; // a1
692 117 x[5] = 0.0; // a2
693 117 x[6] = 0.0; // a0
694 117 x[7] = -1.0; // curveID
695 117 x[8] = 0.0; // surDepth
696 117 x[9] = 0.0; // fEvap
697
698 // --- get surf. area curve name
699
2/2
✓ Branch 0 taken 56 times.
✓ Branch 1 taken 61 times.
 117 if (m == TABULAR)
700 {
701 56 i = project_findObject(CURVE, tok[5]);
702
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 56 times.
 56 if (i < 0) return error_setInpError(ERR_NAME, tok[5]);
703 56 x[7] = i;
704 56 n = 6;
705 }
706
707 // --- get surf. area function coeffs.
708 else
709 {
710
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 61 times.
 61 if (ntoks < 8) return error_setInpError(ERR_ITEMS, "");
711
2/2
✓ Branch 0 taken 183 times.
✓ Branch 1 taken 61 times.
 244 for (i=5; i<=7; i++)
712 {
713
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 183 times.
 183 if ( ! getDouble(tok[i], &y[i-5]) )
714 ✗ return error_setInpError(ERR_NUMBER, tok[i]);
715 }
716 61 n = 8;
717 }
718
719 // --- check for valid data
720
3/3
✓ Branch 0 taken 53 times.
✓ Branch 1 taken 8 times.
✓ Branch 2 taken 56 times.
 117 switch (m)
721 {
722 53 case FUNCTIONAL:
723 // area at 0 depth can't be negative
724
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 53 times.
 53 if (y[2] < 0.0) return error_setInpError(ERR_NUMBER, tok[7]);
725 53 break;
726
727 8 case CYLINDRICAL:
728 case CONICAL:
729 case PARABOLOID:
730 case PYRAMIDAL:
731 // length or width can't be <= 0, slope can't be < 0
732
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
 8 if (y[0] <= 0.0 ) return error_setInpError(ERR_NUMBER, tok[5]);
733
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
 8 if (y[1] <= 0.0) return error_setInpError(ERR_NUMBER, tok[6]);
734
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 8 times.
 8 if (y[2] < 0.0) return error_setInpError(ERR_NUMBER, tok[7]);
735 8 break;
736 }
737 // --- height of paraboloid shape can't be 0
738
3/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 115 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
 117 if (m == PARABOLOID && y[2] == 0.0)
739 ✗ return error_setInpError(ERR_NUMBER, tok[7]);
740
741 // --- convert supplied parameters to coeffs. in surface area equation
742
6/6
✓ Branch 0 taken 53 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 2 times.
✓ Branch 3 taken 2 times.
✓ Branch 4 taken 2 times.
✓ Branch 5 taken 56 times.
 117 switch (m)
743 {
744 53 case FUNCTIONAL:
745 53 x[4] = y[0];
746 53 x[5] = y[1];
747 53 x[6] = y[2];
748 53 break;
749
750 2 case CYLINDRICAL:
751 2 A = y[0] / 2.; // base semi-axis length
752 2 B = y[1] / 2.; // base semi axis width
753 2 x[4] = 0.0; // a1 term
754 2 x[5] = 0.0; // a2 term
755 2 x[6] = PI * A * B; // a0 term
756 2 break;
757
758 2 case CONICAL:
759 2 A = y[0] / 2.; // base semi-axis length
760 2 B = y[1] / 2.; // base semi axis width
761 2 Z = y[2]; // side slope
762 2 x[4] = 2.0 * PI * B * Z; // linear coeff.
763 2 x[5] = PI * B / A * Z * Z; // quadratic coeff.
764 2 x[6] = PI * A * B; // constant
765 2 break;
766
767 2 case PARABOLOID:
768 2 A = y[0] / 2.; // top semi-axis length
769 2 B = y[1] / 2.; // top semi-axis width
770 2 Z = y[2]; // top height
771 2 x[4] = PI * A * B / Z; // a1 term
772 2 x[5] = 0.0; // a2 term
773 2 x[6] = 0.0; // a0 term
774 2 break;
775
776 2 case PYRAMIDAL:
777 2 L = y[0];
778 2 W = y[1];
779 2 Z = y[2];
780 2 x[4] = 2.0 * (L + W) * Z; // linear coeff.
781 2 x[5] = 4.0 * Z * Z; // quadratic coeff.
782 2 x[6] = L * W; // constant
783 2 break;
784 }
785
786 // --- get surcharge depth if present
787
1/2
✓ Branch 0 taken 117 times.
✗ Branch 1 not taken.
 117 if ( ntoks > n)
788 {
789
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 117 times.
 117 if ( ! getDouble(tok[n], &x[8]) )
790 ✗ return error_setInpError(ERR_NUMBER, tok[n]);
791 117 n++;
792 }
793
794 // --- get evaporation fraction if present
795
1/2
✓ Branch 0 taken 117 times.
✗ Branch 1 not taken.
 117 if ( ntoks > n )
796 {
797
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 117 times.
 117 if ( ! getDouble(tok[n], &x[9]) )
798 ✗ return error_setInpError(ERR_NUMBER, tok[n]);
799 117 n++;
800 }
801
802 // --- add parameters to storage unit object
803 117 Node[j].ID = id;
804 117 node_setParams(j, STORAGE, k, x);
805
806 // --- read exfiltration parameters if present
807
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 114 times.
 117 if ( ntoks > n ) return exfil_readStorageParams(k, tok, ntoks, n);
808 114 return 0;
809 }
810
811 //=============================================================================
812
813 14159 double storage_getDepth(int j, double v)
814 //
815 // Input: j = node index
816 // v = volume (ft3)
817 // Output: returns depth of water at a storage node (ft)
818 // Purpose: computes a storage node's water depth from its volume.
819 //
820 {
821 14159 int k = Node[j].subIndex;
822 14159 int i = Storage[k].aCurve;
823 14159 int shape = Storage[k].shape;
824 14159 double a0 = Storage[k].a0;
825 14159 double a1 = Storage[k].a1;
826 14159 double a2 = Storage[k].a2;
827 double d, e;
828 TStorageVol storageVol;
829
830 // --- return max depth if a max. volume has been computed
831 // and volume is > max. volume
832
1/2
✓ Branch 0 taken 14159 times.
✗ Branch 1 not taken.
 14159 if ( Node[j].fullVolume > 0.0
833
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 14159 times.
 14159 && v >= Node[j].fullVolume ) return Node[j].fullDepth;
834
2/2
✓ Branch 0 taken 1230 times.
✓ Branch 1 taken 12929 times.
 14159 if ( v == 0.0 ) return 0.0;
835
836 // --- convert volume to user's units
837 12929 v *= UCF(VOLUME);
838 12929 storageVol.k = k;
839 12929 storageVol.v = v;
840
841
5/6
✓ Branch 0 taken 5669 times.
✓ Branch 1 taken 1452 times.
✓ Branch 2 taken 1453 times.
✓ Branch 3 taken 1453 times.
✓ Branch 4 taken 2902 times.
✗ Branch 5 not taken.
 12929 switch (shape)
842 {
843 5669 case TABULAR:
844 5669 i = Storage[k].aCurve;
845
1/2
✓ Branch 0 taken 5669 times.
✗ Branch 1 not taken.
 5669 if (i >= 0)
846 5669 d = table_getStorageDepth(&Curve[i], v);
847 ✗ else d = 0.0;
848 5669 break;
849
850 1452 case CYLINDRICAL:
851 // area = a0; v = a0*d;
852 1452 d = v / a0;
853 1452 break;
854
855 1453 case PARABOLOID:
856 // area = a1*d; v = (a1/2)*d^2
857 1453 d = sqrt(2.0 * v / a1);
858 1453 break;
859
860 1453 case FUNCTIONAL:
861 // area = a0 + a1; v = (a0 + a1) * d
862
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1453 times.
 1453 if (a2 == 0.0)
863 {
864 ✗ d = v / (a0 + a1);
865 }
866 // area = a1*d^a2; v = a1/(a2+1)*d^(a2+1)
867
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1453 times.
 1453 else if (a0 == 0.0)
868 {
869 ✗ e = 1.0 / (a2 + 1.0);
870 ✗ d = pow(v / (a1 * e), e);
871 }
872 // area = a0 + a1*d; v = a0*d + (a1/2)*d^2
873
2/4
✓ Branch 0 taken 1453 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1453 times.
✗ Branch 3 not taken.
 1453 else if (a2 == 1.0 && a1 > 0.0)
874 {
875 1453 d = (sqrt(a0*a0 + 2.*a1*v) - a0) / a1;
876 }
877 else
878 // area = a0 + a1*d^a2
879 {
880 ✗ d = v / (a0 + a1);
881 ✗ findroot_Newton(0.0, Node[j].fullDepth*UCF(LENGTH), &d,
882 0.001, storage_getVolDiff, &storageVol);
883 }
884 1453 break;
885
886 2902 case CONICAL:
887 case PYRAMIDAL:
888 // area = a0 + a1*d + a2*d^2; v = a0*d + (a1/2)*d^2 + (a2/3)*d^3
889 2902 d = v / a0;
890 2902 findroot_Newton(0.0, Node[j].fullDepth*UCF(LENGTH), &d,
891 0.001, storage_getVolDiff, &storageVol);
892 2902 break;
893
894 ✗ default:
895 ✗ d = 0.0;
896 }
897
898 12929 d /= UCF(LENGTH);
899
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 12929 times.
 12929 if ( d > Node[j].fullDepth )
900 ✗ d = Node[j].fullDepth;
901 12929 return d;
902 }
903
904 //=============================================================================
905
906 5804 void storage_getVolDiff(double y, double* f, double* df, void* p)
907 //
908 // Input: y = depth of water (ft)
909 // p = pointer to a TStorageVol object
910 // Output: f = volume of water (ft3)
911 // df = dVolume/dDepth ( = surface area)(ft2)
912 // Purpose: computes volume difference and its derivative at a storage node
913 // using the node's area versus depth function.
914 //
915 {
916 int k;
917 TStorageVol* storageVol;
918
919 // --- cast void pointer p to a TStorageVol object
920 5804 storageVol = (TStorageVol *)p;
921 5804 k = storageVol->k;
922
923 // --- compute volume & surface area at depth y
924 5804 *f = storage_getVolume(k, y) - storageVol->v;
925 5804 *df = storage_getSurfArea(k, y);
926 5804 }
927
928 //=============================================================================
929
930 491475 double storage_getVolume(int j, double d)
931 //
932 // Input: j = node index
933 // d = depth (ft)
934 // Output: returns volume of stored water (ft3)
935 // Purpose: computes a storage node's water volume from its depth.
936 //
937 {
938 491475 int k = Node[j].subIndex;
939 int i;
940 double n, v;
941
942 // --- return full volume if depth >= max. depth
943
2/2
✓ Branch 0 taken 3990 times.
✓ Branch 1 taken 487485 times.
 491475 if ( d == 0.0 ) return 0.0;
944
2/2
✓ Branch 0 taken 234 times.
✓ Branch 1 taken 487251 times.
 487485 if ( d >= Node[j].fullDepth
945
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 234 times.
 234 && Node[j].fullVolume > 0.0 ) return Node[j].fullVolume;
946
947
3/4
✓ Branch 0 taken 44927 times.
✓ Branch 1 taken 379114 times.
✓ Branch 2 taken 63444 times.
✗ Branch 3 not taken.
 487485 switch (Storage[k].shape)
948 {
949 // --- for tabular shape function, use end area method
950 44927 case TABULAR:
951 44927 i = Storage[k].aCurve;
952
1/2
✓ Branch 0 taken 44927 times.
✗ Branch 1 not taken.
 44927 if (i >= 0)
953 89854 return table_getStorageVolume(&Curve[i], d*UCF(LENGTH)) /
954 44927 UCF(VOLUME);
955 ✗ else return 0.0;
956
957 // --- for FUNCTIONAL relation, integrate a0 + a1*d^a2
958 379114 case FUNCTIONAL:
959 379114 d *= UCF(LENGTH);
960 379114 n = Storage[k].a2 + 1.0;
961 379114 v = (Storage[k].a0 * d) + Storage[k].a1 / n * pow(d, n);
962 379114 return v / UCF(VOLUME);
963
964 // --- for other shapes evaluate cubic eqn. a0*d + (a1/2)*d^2 + (a2/3)*d^3
965 63444 case CYLINDRICAL:
966 case CONICAL:
967 case PARABOLOID:
968 case PYRAMIDAL:
969 63444 d *= UCF(LENGTH);
970 63444 v = d * (Storage[k].a0 + d * (Storage[k].a1 / 2.0 + d * Storage[k].a2 / 3.0));
971 63444 return v / UCF(VOLUME);
972
973 ✗ default: return 0.0;
974 }
975 }
976
977 //=============================================================================
978
979 577948 double storage_getSurfArea(int j, double d)
980 //
981 // Input: j = node index
982 // d = depth (ft)
983 // Output: returns surface area (ft2)
984 // Purpose: computes a storage node's surface area from its water depth.
985 //
986 {
987 577948 double area = 0.0;
988 577948 int k = Node[j].subIndex;
989 int i;
990
991
3/4
✓ Branch 0 taken 66696 times.
✓ Branch 1 taken 413276 times.
✓ Branch 2 taken 97976 times.
✗ Branch 3 not taken.
 577948 switch (Storage[k].shape)
992 {
993 // --- for tabular shape function, use table look-up
994 66696 case TABULAR:
995 66696 i = Storage[k].aCurve;
996
1/2
✓ Branch 0 taken 66696 times.
✗ Branch 1 not taken.
 66696 if (i >= 0)
997 66696 area = table_lookupEx(&Curve[i], d*UCF(LENGTH));
998 66696 break;
999
1000 // --- for FUNCTIONAL relation, evaluate a0 + a1*d^a2
1001 413276 case FUNCTIONAL:
1002 826552 area = Storage[k].a0 + Storage[k].a1 *
1003 413276 pow(d*UCF(LENGTH), Storage[k].a2);
1004 413276 break;
1005
1006 // --- for other shapes, evaluate quadratic a0 + a1*d + a2*d^2
1007 97976 case CYLINDRICAL:
1008 case CONICAL:
1009 case PARABOLOID:
1010 case PYRAMIDAL:
1011 97976 d *= UCF(LENGTH);
1012 97976 area = Storage[k].a0 + d * (Storage[k].a1 + d * Storage[k].a2);
1013 97976 break;
1014
1015 ✗ default: return 0.0;
1016 }
1017 577948 return area / UCF(LENGTH) / UCF(LENGTH);
1018 }
1019
1020 //=============================================================================
1021
1022 17353 double storage_getOutflow(int j, int i)
1023 //
1024 // Input: j = node index
1025 // i = link index
1026 // Output: returns flow from storage node into conduit link (cfs)
1027 // Purpose: finds outflow from a storage node into its connecting conduit link
1028 // ( non-conduit links have their own getInflow functions).
1029 //
1030 {
1031 int k;
1032 double a, y;
1033
1034 // --- link must be a conduit
1035
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 17353 times.
 17353 if ( Link[i].type != CONDUIT ) return 0.0;
1036
1037 // --- find depth of water in conduit
1038 17353 y = Node[j].newDepth - Link[i].offset1;
1039
1040 // --- return 0 if conduit empty or full flow if full
1041
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 17353 times.
 17353 if ( y <= 0.0 ) return 0.0;
1042
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 17353 times.
 17353 if ( y >= Link[i].xsect.yFull ) return Link[i].qFull;
1043
1044 // --- if partially full, return normal flow
1045 17353 k = Link[i].subIndex;
1046 17353 a = xsect_getAofY(&Link[i].xsect, y);
1047 17353 return Conduit[k].beta * xsect_getSofA(&Link[i].xsect, a);
1048 }
1049
1050 //=============================================================================
1051
1052 143130 double storage_getLosses(int j, double tStep)
1053 //
1054 // Input: j = node index
1055 // tStep = time step (sec)
1056 // Output: returns evaporation + infiltration rate (cfs)
1057 // Purpose: computes combined rate of water evaporated & infiltrated from
1058 // a storage node.
1059 //
1060 {
1061 int k;
1062 double depth;
1063 double area;
1064 143130 double evapRate = 0.0;
1065 143130 double exfilRate = 0.0;
1066 143130 double totalLoss = 0.0;
1067 double lossRatio;
1068 TExfil* exfil;
1069
1070 // --- get node's evap. rate (ft/s) & exfiltration object
1071 143130 k = Node[j].subIndex;
1072 143130 evapRate = Evap.rate * Storage[k].fEvap;
1073 143130 exfil = Storage[k].exfil;
1074
1075 // --- if either of these apply
1076
4/4
✓ Branch 0 taken 85776 times.
✓ Branch 1 taken 57354 times.
✓ Branch 2 taken 19083 times.
✓ Branch 3 taken 66693 times.
 143130 if ( evapRate > 0.0 || exfil != NULL)
1077 {
1078 // --- obtain storage depth & surface area
1079 76437 depth = Node[j].newDepth;
1080 76437 area = storage_getSurfArea(j, depth);
1081
1082 // --- compute evap rate over this area (cfs)
1083
2/2
✓ Branch 0 taken 75985 times.
✓ Branch 1 taken 452 times.
 76437 if (Node[j].newVolume > FUDGE)
1084 75985 evapRate = area * evapRate;
1085
1086 // --- find exfiltration rate (cfs) through bottom and side banks
1087
2/2
✓ Branch 0 taken 19083 times.
✓ Branch 1 taken 57354 times.
 76437 if ( exfil != NULL )
1088 {
1089 19083 exfilRate = exfil_getLoss(exfil, tStep, depth, area);
1090 }
1091
1092 // --- total loss over time step cannot exceed stored volume
1093 76437 totalLoss = (evapRate + exfilRate) * tStep;
1094
2/2
✓ Branch 0 taken 450 times.
✓ Branch 1 taken 75987 times.
 76437 if ( totalLoss > Node[j].newVolume )
1095 {
1096 450 lossRatio = Node[j].newVolume / totalLoss;
1097 450 evapRate *= lossRatio;
1098 450 exfilRate *= lossRatio;
1099 }
1100 }
1101
1102 // --- save evap & infil losses at the node
1103 143130 Storage[Node[j].subIndex].evapLoss = evapRate * tStep;
1104 143130 Storage[Node[j].subIndex].exfilLoss = exfilRate * tStep;
1105 143130 return evapRate + exfilRate;
1106 }
1107
1108
1109 //=============================================================================
1110 // D I V I D E R M E T H O D S
1111 //=============================================================================
1112
1113 4 int divider_readParams(int j, int k, char* tok[], int ntoks)
1114 //
1115 // Input: j = node index
1116 // k = divider index
1117 // tok[] = array of string tokens
1118 // ntoks = number of tokens
1119 // Output: returns an error message
1120 // Purpose: reads a flow divider's properties from a tokenized line of input.
1121 //
1122 // Format of input line is:
1123 // nodeID elev divLink TABULAR curveID (optional params)
1124 // nodeID elev divLink OVERFLOW (optional params)
1125 // nodeID elev divLink CUTOFF qCutoff (optional params)
1126 // nodeID elev divLink WEIR qMin dhMax cWeir (optional params)
1127 // where optional params are:
1128 // maxDepth initDepth surDepth aPond
1129 //
1130 {
1131 int i, m, m1, m2, n;
1132 double x[11];
1133 char* id;
1134
1135 // --- get ID name
1136
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if ( ntoks < 4 ) return error_setInpError(ERR_ITEMS, "");
1137 4 id = project_findID(NODE, tok[0]);
1138
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if ( id == NULL ) return error_setInpError(ERR_NAME, tok[0]);
1139
1140 // --- get invert elev.
1141
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 4 times.
 4 if ( ! getDouble(tok[1], &x[0]) ) return error_setInpError(ERR_NUMBER, tok[1]);
1142
1143 // --- initialize parameter values
1144
2/2
✓ Branch 0 taken 40 times.
✓ Branch 1 taken 4 times.
 44 for ( i=1; i<11; i++) x[i] = 0.0;
1145
1146 // --- check if no diverted link supplied
1147
2/4
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
 4 if ( strlen(tok[2]) == 0 || strcmp(tok[2], "*") == 0 ) x[1] = -1.0;
1148
1149 // --- otherwise get index of diverted link
1150 else
1151 {
1152 4 m1 = project_findObject(LINK, tok[2]);
1153
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if ( m1 < 0 ) return error_setInpError(ERR_NAME, tok[2]);
1154 4 x[1] = m1;
1155 }
1156
1157 // --- get divider type
1158 4 n = 4;
1159 4 m1 = findmatch(tok[3], DividerTypeWords);
1160
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
 4 if ( m1 < 0 ) return error_setInpError(ERR_KEYWORD, tok[3]);
1161 4 x[2] = m1;
1162
1163 // --- get index of flow diversion curve for Tabular divider
1164 4 x[3] = -1;
1165
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3 times.
 4 if ( m1 == TABULAR_DIVIDER )
1166 {
1167
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( ntoks < 5 ) return error_setInpError(ERR_ITEMS, "");
1168 1 m2 = project_findObject(CURVE, tok[4]);
1169
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( m2 < 0 ) return error_setInpError(ERR_NAME, tok[4]);
1170 1 x[3] = m2;
1171 1 n = 5;
1172 }
1173
1174 // --- get cutoff flow for Cutoff divider
1175
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3 times.
 4 if ( m1 == CUTOFF_DIVIDER )
1176 {
1177
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( ntoks < 5 ) return error_setInpError(ERR_ITEMS, "");
1178
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
 1 if ( ! getDouble(tok[4], &x[4]) )
1179 ✗ return error_setInpError(ERR_NUMBER, tok[4]);
1180 1 n = 5;
1181 }
1182
1183 // --- get qmin, dhMax, & cWeir for Weir divider
1184
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3 times.
 4 if ( m1 == WEIR_DIVIDER )
1185 {
1186
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( ntoks < 7 ) return error_setInpError(ERR_ITEMS, "");
1187
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 1 time.
 4 for (i=4; i<7; i++)
1188
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
 3 if ( ! getDouble(tok[i], &x[i]) )
1189 ✗ return error_setInpError(ERR_NUMBER, tok[i]);
1190 1 n = 7;
1191 }
1192
1193 // --- no parameters needed for Overflow divider
1194
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3 times.
 4 if ( m1 == OVERFLOW_DIVIDER ) n = 4;
1195
1196 // --- retrieve optional full depth, init. depth, surcharged depth
1197 // & ponded area
1198 4 m = 7;
1199
3/4
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 4 times.
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 20 for (i=n; i<ntoks && m<11; i++)
1200 {
1201
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 16 times.
 16 if ( ! getDouble(tok[i], &x[m]) )
1202 {
1203 ✗ return error_setInpError(ERR_NUMBER, tok[i]);
1204 }
1205 16 m++;
1206 }
1207
1208 // --- add parameters to data base
1209 4 Node[j].ID = id;
1210 4 node_setParams(j, DIVIDER, k, x);
1211 4 return 0;
1212 }
1213
1214 //=============================================================================
1215
1216 4 void divider_validate(int j)
1217 //
1218 // Input: j = node index
1219 // Output: none
1220 // Purpose: validates a flow divider's properties.
1221 //
1222 {
1223 int i, k;
1224
1225 // --- check that diverted link is attached to divider
1226 4 k = Node[j].subIndex;
1227 4 i = Divider[k].link;
1228
2/6
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
 4 if ( i < 0 || ( Link[i].node1 != j && Link[i].node2 != j) )
1229 {
1230 ✗ report_writeErrorMsg(ERR_DIVIDER_LINK, Node[j].ID);
1231 }
1232
1233 // --- validate parameters supplied for weir-type divider
1234
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 3 times.
 4 if ( Divider[k].type == WEIR_DIVIDER )
1235 {
1236
2/4
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 time.
 1 if ( Divider[k].dhMax <= 0.0 || Divider[k].cWeir <= 0.0 )
1237 ✗ report_writeErrorMsg(ERR_WEIR_DIVIDER, Node[j].ID);
1238 else
1239 {
1240 // --- find flow when weir is full
1241 2 Divider[k].qMax = Divider[k].cWeir * pow(Divider[k].dhMax, 1.5)
1242 1 / UCF(FLOW);
1243
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( Divider[k].qMin > Divider[k].qMax )
1244 ✗ report_writeErrorMsg(ERR_WEIR_DIVIDER, Node[j].ID);
1245 }
1246 }
1247 4 }
1248
1249 //=============================================================================
1250
1251 1920 double divider_getOutflow(int j, int k)
1252 //
1253 // Input: j = node index
1254 // k = index of diversion link
1255 // Output: returns diverted flow rate (cfs)
1256 // Purpose: computes flow sent through divider node into its diversion link.
1257 //
1258 // NOTE: requires that links be previously sorted so that the non-diversion
1259 // link always gets evaluated before the diversion link
1260 {
1261 int i; // index of divider node
1262 int m; // index of diverted flow table
1263 double qIn; // inflow to divider
1264 double qOut; // diverted outflow
1265 double f; // fraction of weir divider full
1266
1267 1920 qIn = Node[j].inflow + Node[j].overflow;
1268 1920 i = Node[j].subIndex;
1269
4/5
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 480 times.
✓ Branch 2 taken 480 times.
✓ Branch 3 taken 480 times.
✗ Branch 4 not taken.
 1920 switch ( Divider[i].type )
1270 {
1271 480 case CUTOFF_DIVIDER:
1272
2/2
✓ Branch 0 taken 200 times.
✓ Branch 1 taken 280 times.
 480 if ( qIn <= Divider[i].qMin ) qOut = 0.0;
1273 280 else qOut = qIn - Divider[i].qMin;
1274 480 break;
1275
1276 480 case OVERFLOW_DIVIDER:
1277 // --- outflow sent into non-diversion link is simply node's inflow
1278
2/2
✓ Branch 0 taken 240 times.
✓ Branch 1 taken 240 times.
 480 if ( k != Divider[i].link ) qOut = qIn;
1279
1280 // --- diversion link receives any excess of node's inflow and
1281 // outflow sent previously into non-diversion link
1282 240 else qOut = qIn - Node[j].outflow;
1283
2/2
✓ Branch 0 taken 241 times.
✓ Branch 1 taken 239 times.
 480 if ( qOut < FLOW_TOL ) qOut = 0.0;
1284 480 return qOut;
1285
1286 480 case WEIR_DIVIDER:
1287 // --- no flow if inflow < qMin
1288
2/2
✓ Branch 0 taken 120 times.
✓ Branch 1 taken 360 times.
 480 if ( qIn <= Divider[i].qMin ) qOut = 0.0;
1289
1290 // --- otherwise use weir eqn.
1291 else
1292 {
1293 // --- find fractional depth of flow over weir
1294 360 f = (qIn - Divider[i].qMin) /
1295 360 (Divider[i].qMax - Divider[i].qMin);
1296
1297 // --- if weir surcharged, use orifice eqn.
1298
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 360 times.
 360 if ( f > 1.0 ) qOut = Divider[i].qMax * sqrt(f);
1299
1300 // --- otherwise use weir eqn.
1301 360 else qOut = Divider[i].cWeir *
1302 360 pow(f*Divider[i].dhMax, 1.5) / UCF(FLOW);
1303 }
1304 480 break;
1305
1306 480 case TABULAR_DIVIDER:
1307 480 m = Divider[i].flowCurve;
1308
1/2
✓ Branch 0 taken 480 times.
✗ Branch 1 not taken.
 480 if ( m >= 0 )
1309 480 qOut = table_lookup(&Curve[m], qIn * UCF(FLOW)) / UCF(FLOW);
1310 ✗ else qOut = 0.0;
1311 480 break;
1312
1313 ✗ default: qOut = 0.0;
1314 }
1315
1316 // --- make sure outflow doesn't exceed inflow
1317
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1440 times.
 1440 if ( qOut > qIn ) qOut = qIn;
1318
1319 // --- if link k not the diversion link, then re-define qOut as
1320 // the undiverted flow
1321
2/2
✓ Branch 0 taken 720 times.
✓ Branch 1 taken 720 times.
 1440 if ( k != Divider[i].link )
1322 {
1323 720 qOut = qIn - qOut;
1324 }
1325 1440 return qOut;
1326 }
1327
1328
1329 //=============================================================================
1330 // O U T F A L L M E T H O D S
1331 //=============================================================================
1332
1333 348 int outfall_readParams(int j, int k, char* tok[], int ntoks)
1334 //
1335 // Input: j = node index
1336 // k = outfall index
1337 // tok[] = array of string tokens
1338 // ntoks = number of tokens
1339 // Output: returns an error message
1340 // Purpose: reads an outfall's properties from a tokenized line of input.
1341 //
1342 // Format of input line is:
1343 // nodeID elev FIXED fixedStage (flapGate) (routeTo)
1344 // nodeID elev TIDAL curveID (flapGate) (routeTo)
1345 // nodeID elev TIMESERIES tseriesID (flapGate) (routTo)
1346 // nodeID elev FREE (flapGate) (routeTo)
1347 // nodeID elev NORMAL (flapGate) (routeTo)
1348 //
1349 {
1350 int i, m, n;
1351 double x[7];
1352 char* id;
1353
1354
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 348 times.
 348 if ( ntoks < 3 ) return error_setInpError(ERR_ITEMS, "");
1355 348 id = project_findID(NODE, tok[0]); // node ID
1356
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 348 times.
 348 if ( id == NULL )
1357 ✗ return error_setInpError(ERR_NAME, tok[0]);
1358
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 348 times.
 348 if ( ! getDouble(tok[1], &x[0]) ) // invert elev.
1359 ✗ return error_setInpError(ERR_NUMBER, tok[1]);
1360 348 i = findmatch(tok[2], OutfallTypeWords); // outfall type
1361
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 348 times.
 348 if ( i < 0 ) return error_setInpError(ERR_KEYWORD, tok[2]);
1362 348 x[1] = i; // outfall type
1363 348 x[2] = 0.0; // fixed stage
1364 348 x[3] = -1.; // tidal curve
1365 348 x[4] = -1.; // tide series
1366 348 x[5] = 0.; // flap gate
1367 348 x[6] = -1.; // route to subcatch
1368
1369 348 n = 4;
1370
2/2
✓ Branch 0 taken 247 times.
✓ Branch 1 taken 101 times.
 348 if ( i >= FIXED_OUTFALL )
1371 {
1372
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 247 times.
 247 if ( ntoks < 4 ) return error_setInpError(ERR_ITEMS, "");
1373 247 n = 5;
1374
3/4
✓ Branch 0 taken 19 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 226 times.
✗ Branch 3 not taken.
 247 switch ( i )
1375 {
1376 19 case FIXED_OUTFALL: // fixed stage
1377
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 19 times.
 19 if ( ! getDouble(tok[3], &x[2]) )
1378 ✗ return error_setInpError(ERR_NUMBER, tok[3]);
1379 19 break;
1380 2 case TIDAL_OUTFALL: // tidal curve
1381 2 m = project_findObject(CURVE, tok[3]);
1382
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( m < 0 ) return error_setInpError(ERR_NAME, tok[3]);
1383 2 x[3] = m;
1384 2 break;
1385 226 case TIMESERIES_OUTFALL: // stage time series
1386 226 m = project_findObject(TSERIES, tok[3]);
1387
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 226 times.
 226 if ( m < 0 ) return error_setInpError(ERR_NAME, tok[3]);
1388 226 x[4] = m;
1389 226 Tseries[m].refersTo = TIMESERIES_OUTFALL;
1390 }
1391 }
1392
2/2
✓ Branch 0 taken 347 times.
✓ Branch 1 taken 1 time.
 348 if ( ntoks == n )
1393 {
1394 347 m = findmatch(tok[n-1], NoYesWords); // flap gate
1395
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 347 times.
 347 if ( m < 0 ) return error_setInpError(ERR_KEYWORD, tok[n-1]);
1396 347 x[5] = m;
1397 }
1398
1399
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 347 times.
 348 if ( ntoks == n+1)
1400 {
1401 1 m = project_findObject(SUBCATCH, tok[n]);
1402
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( m < 0 ) return error_setInpError(ERR_NAME, tok[n]);
1403 1 x[6] = m;
1404 }
1405
1406 348 Node[j].ID = id;
1407 348 node_setParams(j, OUTFALL, k, x);
1408 348 return 0;
1409 }
1410
1411 //=============================================================================
1412
1413 1979248 void outfall_setOutletDepth(int j, double yNorm, double yCrit, double z)
1414 //
1415 // Input: j = node index
1416 // yNorm = normal flow depth in outfall conduit (ft)
1417 // yCrit = critical flow depth in outfall conduit (ft)
1418 // z = height to outfall conduit invert (ft)
1419 // Output: none
1420 // Purpose: sets water depth at an outfall node.
1421 //
1422 {
1423 double x, y; // x,y values in table
1424 double yNew; // new depth above invert elev. (ft)
1425 double stage; // water elevation at outfall (ft)
1426 int k; // table index
1427 1979248 int i = Node[j].subIndex; // outfall index
1428 DateTime currentDate; // current date/time in days
1429
1430
5/6
✓ Branch 0 taken 1651815 times.
✓ Branch 1 taken 48901 times.
✓ Branch 2 taken 223052 times.
✓ Branch 3 taken 2892 times.
✓ Branch 4 taken 52588 times.
✗ Branch 5 not taken.
 1979248 switch ( Outfall[i].type )
1431 {
1432 1651815 case FREE_OUTFALL:
1433
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1651815 times.
 1651815 if ( z > 0.0 ) Node[j].newDepth = 0.0;
1434
2/2
✓ Branch 0 taken 508221 times.
✓ Branch 1 taken 1143594 times.
 1651815 else Node[j].newDepth = MIN(yNorm, yCrit);
1435 1700716 return;
1436
1437 48901 case NORMAL_OUTFALL:
1438
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 48901 times.
 48901 if ( z > 0.0 ) Node[j].newDepth = 0.0;
1439 48901 else Node[j].newDepth = yNorm;
1440 48901 return;
1441
1442 223052 case FIXED_OUTFALL:
1443 223052 stage = Outfall[i].fixedStage;
1444 223052 break;
1445
1446 2892 case TIDAL_OUTFALL:
1447 2892 k = Outfall[i].tideCurve;
1448 2892 table_getFirstEntry(&Curve[k], &x, &y);
1449 2892 currentDate = NewRoutingTime / MSECperDAY;
1450 2892 x += ( currentDate - floor(currentDate) ) * 24.0;
1451 2892 stage = table_lookup(&Curve[k], x) / UCF(LENGTH);
1452 2892 break;
1453
1454 52588 case TIMESERIES_OUTFALL:
1455 52588 k = Outfall[i].stageSeries;
1456 52588 currentDate = StartDateTime + NewRoutingTime / MSECperDAY;
1457 52588 stage = table_tseriesLookup(&Tseries[k], currentDate, TRUE) /
1458 52588 UCF(LENGTH);
1459 52588 break;
1460 ✗ default: stage = Node[j].invertElev;
1461 }
1462
1463 // --- now determine depth at node given outfall stage elev.
1464
1465 // --- let critical flow depth be min. of critical & normal depth
1466
2/2
✓ Branch 0 taken 242939 times.
✓ Branch 1 taken 35593 times.
 278532 yCrit = MIN(yCrit, yNorm);
1467
1468 // --- if elev. of critical depth is below outfall stage elev. then
1469 // the outfall stage determines node depth
1470
2/2
✓ Branch 0 taken 261290 times.
✓ Branch 1 taken 17242 times.
 278532 if ( yCrit + z + Node[j].invertElev < stage )
1471 {
1472 261290 yNew = stage - Node[j].invertElev;
1473 }
1474
1475 // --- otherwise if the outfall conduit lies above the outfall invert
1476
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 17242 times.
 17242 else if ( z > 0.0 )
1477 {
1478 // --- if the outfall stage lies below the bottom of the outfall
1479 // conduit then the result is distance from node invert to stage
1480 ✗ if ( stage < Node[j].invertElev + z )
1481 ✗ yNew = MAX(0.0, (stage - Node[j].invertElev));
1482
1483 // --- otherwise stage lies between bottom of conduit and critical
1484 // depth in conduit so result is elev. of critical depth
1485 ✗ else yNew = z + yCrit;
1486 }
1487
1488 // --- and for case where there is no conduit offset and outfall stage
1489 // lies below critical depth, then node depth = critical depth
1490 17242 else yNew = yCrit;
1491 278532 Node[j].newDepth = yNew;
1492 }
1493
1494 //=============================================================================
1495