GCC Code Coverage Report

Directory: src/solver/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 88.6% 954 / 0 / 1077
Functions: 98.1% 53 / 0 / 54
Branches: 66.5% 525 / 0 / 790
link.c
Line Branch Exec Source
1 //-----------------------------------------------------------------------------
2 // link.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 // Conveyance system link functions
11 //
12 // Update History
13 // ==============
14 // Build 5.1.007:
15 // - Optional surcharging of weirs introduced.
16 // Build 5.1.008:
17 // - Bug in finding flow through surcharged weir fixed.
18 // - Bug in finding if conduit is upstrm/dnstrm full fixed.
19 // - Monthly conductivity adjustment applied to conduit seepage.
20 // - Conduit seepage limited by conduit's flow rate.
21 // Build 5.1.010:
22 // - Support added for new ROADWAY_WEIR object.
23 // - Time of last setting change initialized for links.
24 // Build 5.1.011:
25 // - Crest elevation of regulator links raised to downstream invert.
26 // - Fixed converting roadWidth weir parameter to internal units.
27 // - Weir shape parameter deprecated.
28 // - Extra geometric parameters ignored for non-conduit open rectangular
29 // cross sections.
30 // Build 5.1.012:
31 // - Conduit seepage rate now based on flow width, not wetted perimeter.
32 // - Formula for side flow weir corrected.
33 // - Crest length contraction adjustments corrected.
34 // Build 5.1.013:
35 // - Maximum depth adjustments made for storage units that can surcharge.
36 // - Support added for head-dependent weir coefficient curves.
37 // - Adjustment of regulator link crest offset to match downstream node invert
38 // now only done for Dynamic Wave flow routing.
39 // Build 5.1.014:
40 // - Conduit evap. and seepage losses initialized to 0 in conduit_initState()
41 // and not allowed to exceed current flow rate in conduit_getLossRate().
42 // Build 5.2.0:
43 // - Support added for Streets and Inlets.
44 // - Support added for variable speed pumps.
45 // Build 5.2.1
46 // - Warning no longer issued when conduit elevation drop < MIN_DELTA_Z.
47 // Build 5.2.2:
48 // - Warning for conduit elevation drop < MIN_DELTA_Z restored.
49 // Build 5.2.4:
50 // - Conduit evap+seepage loss under DW routing limited by conduit volume.
51 //-----------------------------------------------------------------------------
52 #define _CRT_SECURE_NO_DEPRECATE
53
54 #include <string.h>
55 #include <stdlib.h>
56 #include <math.h>
57 #include "headers.h"
58 #include "inlet.h"
59
60 //-----------------------------------------------------------------------------
61 // Constants
62 //-----------------------------------------------------------------------------
63 static const double MIN_DELTA_Z = 0.001; // minimum elevation change for conduit
64 // slopes (ft)
65
66 //-----------------------------------------------------------------------------
67 // External functions (declared in funcs.h)
68 //-----------------------------------------------------------------------------
69 // link_readParams (called by parseLine in input.c)
70 // link_readXsectParams (called by parseLine in input.c)
71 // link_readLossParams (called by parseLine in input.c)
72 // link_validate (called by project_validate in project.c)
73 // link_initState (called by initObjects in swmm5.c)
74 // link_setOldHydState (called by routing_execute in routing.c)
75 // link_setOldQualState (called by routing_execute in routing.c)
76 // link_setTargetSetting (called by routing_execute in routing.c)
77 // link_setSetting (called by routing_execute in routing.c)
78 // link_getResults (called by output_saveLinkResults)
79 // link_getLength (called in dwflow.c, kinwave.c & flowrout.c)
80 // link_getFroude (called in dwflow.c)
81 // link_getInflow (called in flowrout.c & dynwave.c)
82 // link_setOutfallDepth (called in flowrout.c & dynwave.c)
83 // link_getYcrit (called by link_setOutfallDepth & in dwflow.c)
84 // link_getYnorm (called by conduit_initState, link_setOutfallDepth & in dwflow.c)
85 // link_getVelocity (called by link_getResults & stats_updateLinkStats)
86 // link_getPower (called by stats_updateLinkStats in stats.c)
87 // link_getLossRate (called in dwflow.c, kinwave.c & flowrout.c)
88
89 //-----------------------------------------------------------------------------
90 // Local functions
91 //-----------------------------------------------------------------------------
92 static void link_setParams(int j, int type, int n1, int n2, int k, double x[]);
93 static void link_convertOffsets(int j);
94 static double link_getOffsetHeight(int j, double offset, double elev);
95
96 static int conduit_readParams(int j, int k, char* tok[], int ntoks);
97 static void conduit_validate(int j, int k);
98 static void conduit_initState(int j, int k);
99 static void conduit_reverse(int j, int k);
100 static double conduit_getLength(int j);
101 static double conduit_getLengthFactor(int j, int k, double roughness);
102 static double conduit_getSlope(int j);
103 static double conduit_getInflow(int j);
104 static double conduit_getLossRate(int j, int routeModel, double q,
105 double tstep);
106
107 static int pump_readParams(int j, int k, char* tok[], int ntoks);
108 static void pump_validate(int j, int k);
109 static void pump_initState(int j, int k);
110 static double pump_getInflow(int j);
111
112 static int orifice_readParams(int j, int k, char* tok[], int ntoks);
113 static void orifice_validate(int j, int k);
114 static void orifice_setSetting(int j, double tstep);
115 static double orifice_getWeirCoeff(int j, int k, double h);
116 static double orifice_getInflow(int j);
117 static double orifice_getFlow(int j, int k, double head, double f,
118 int hasFlapGate);
119
120 static int weir_readParams(int j, int k, char* tok[], int ntoks);
121 static void weir_validate(int j, int k);
122 static void weir_setSetting(int j);
123 static double weir_getInflow(int j);
124 static double weir_getOpenArea(int j, double y);
125 static void weir_getFlow(int j, int k, double head, double dir,
126 int hasFlapGate, double* q1, double* q2);
127 static double weir_getOrificeFlow(int j, double head, double y, double cOrif);
128 static double weir_getdqdh(int k, double dir, double h, double q1, double q2);
129
130 static int outlet_readParams(int j, int k, char* tok[], int ntoks);
131 static double outlet_getFlow(int k, double head);
132 static double outlet_getInflow(int j);
133
134
135 //=============================================================================
136
137 10508 int link_readParams(int j, int type, int k, char* tok[], int ntoks)
138 //
139 // Input: j = link index
140 // type = link type code
141 // k = link type index
142 // tok[] = array of string tokens
143 // ntoks = number of tokens
144 // Output: returns an error code
145 // Purpose: reads parameters for a specific type of link from a
146 // tokenized line of input data.
147 //
148 {
149
5/6
✓ Branch 0 taken 9772 times.
✓ Branch 1 taken 106 times.
✓ Branch 2 taken 159 times.
✓ Branch 3 taken 469 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 10508 switch ( type )
150 {
151 9772 case CONDUIT: return conduit_readParams(j, k, tok, ntoks);
152 106 case PUMP: return pump_readParams(j, k, tok, ntoks);
153 159 case ORIFICE: return orifice_readParams(j, k, tok, ntoks);
154 469 case WEIR: return weir_readParams(j, k, tok, ntoks);
155 2 case OUTLET: return outlet_readParams(j, k, tok, ntoks);
156 ✗ default: return 0;
157 }
158 }
159
160 //=============================================================================
161
162 10400 int link_readXsectParams(char* tok[], int ntoks)
163 //
164 // Input: tok[] = array of string tokens
165 // ntoks = number of tokens
166 // Output: returns an error code
167 // Purpose: reads a link's cross section parameters from a tokenized
168 // line of input data.
169 // Formats:
170 // Link Shape Geom1 Geom2 Geom3 Geom4 (Barrels Culvert)
171 // Link IRREGULAR TransectID
172 // Link STREET StreetID
173 //
174 {
175 int i, j, k;
176 double x[4];
177
178 // --- check for minimum number of tokens
179
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10400 times.
 10400 if (ntoks < 3) return error_setInpError(ERR_ITEMS, "");
180
181 // --- get index of link
182 10400 j = project_findObject(LINK, tok[0]);
183
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10400 times.
 10400 if ( j < 0 ) return error_setInpError(ERR_NAME, tok[0]);
184
185 // --- get code of xsection shape
186 10400 k = findmatch(tok[1], XsectTypeWords);
187
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10400 times.
 10400 if ( k < 0 ) return error_setInpError(ERR_KEYWORD, tok[1]);
188
189 // --- assign default number of barrels to conduit
190
2/2
✓ Branch 0 taken 9772 times.
✓ Branch 1 taken 628 times.
 10400 if ( Link[j].type == CONDUIT ) Conduit[Link[j].subIndex].barrels = 1;
191
192 // --- assume link is not a culvert
193 10400 Link[j].xsect.culvertCode = 0;
194
195 // --- for irregular shape, find index of transect object
196
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 10399 times.
 10400 if ( k == IRREGULAR )
197 {
198 1 i = project_findObject(TRANSECT, tok[2]);
199
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( i < 0 ) return error_setInpError(ERR_NAME, tok[2]);
200 1 Link[j].xsect.type = k;
201 1 Link[j].xsect.transect = i;
202 1 return 0;
203 }
204
205 // --- for street cross section, find index of Street object
206
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 10387 times.
 10399 else if (k == STREET_XSECT)
207 {
208 12 i = project_findObject(STREET, tok[2]);
209
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 12 times.
 12 if (i < 0) return error_setInpError(ERR_NAME, tok[2]);
210 12 Link[j].xsect.type = k;
211 12 Link[j].xsect.transect = i;
212 12 return 0;
213 }
214
215 else
216 {
217 // --- check that geometric parameters are present
218
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10387 times.
 10387 if (ntoks < 6) return error_setInpError(ERR_ITEMS, "");
219
220 // --- parse max. depth & shape curve for a custom shape
221
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 10377 times.
 10387 if ( k == CUSTOM )
222 {
223
2/4
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 10 times.
 10 if ( !getDouble(tok[2], &x[0]) || x[0] <= 0.0 )
224 ✗ return error_setInpError(ERR_NUMBER, tok[2]);
225 10 i = project_findObject(CURVE, tok[3]);
226
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 10 if ( i < 0 ) return error_setInpError(ERR_NAME, tok[3]);
227 10 Link[j].xsect.type = k;
228 10 Link[j].xsect.transect = i;
229 10 Link[j].xsect.yFull = x[0] / UCF(LENGTH);
230 }
231
232 // --- parse and save geometric parameters
233
2/2
✓ Branch 0 taken 41508 times.
✓ Branch 1 taken 10377 times.
 51885 else for (i = 2; i <= 5; i++)
234 {
235
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 41508 times.
 41508 if ( !getDouble(tok[i], &x[i-2]) )
236 ✗ return error_setInpError(ERR_NUMBER, tok[i]);
237 }
238
239 // --- ignore extra parameters for non-conduit open rectangular shapes
240
4/4
✓ Branch 0 taken 628 times.
✓ Branch 1 taken 9759 times.
✓ Branch 2 taken 466 times.
✓ Branch 3 taken 162 times.
 10387 if ( Link[j].type != CONDUIT && k == RECT_OPEN )
241 {
242 466 x[2] = 0.0;
243 466 x[3] = 0.0;
244 }
245
1/2
✗ Branch 2 not taken.
✓ Branch 3 taken 10387 times.
 10387 if ( !xsect_setParams(&Link[j].xsect, k, x, UCF(LENGTH)) )
246 {
247 ✗ return error_setInpError(ERR_NUMBER, "");
248 }
249
250 // --- parse number of barrels if present
251
3/4
✓ Branch 0 taken 9759 times.
✓ Branch 1 taken 628 times.
✓ Branch 2 taken 9759 times.
✗ Branch 3 not taken.
 10387 if ( Link[j].type == CONDUIT && ntoks >= 7 )
252 {
253 9759 i = atoi(tok[6]);
254
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9759 times.
 9759 if ( i <= 0 ) return error_setInpError(ERR_NUMBER, tok[6]);
255 9759 else Conduit[Link[j].subIndex].barrels = (char)i;
256 }
257
258 // --- parse culvert code if present
259
4/4
✓ Branch 0 taken 9759 times.
✓ Branch 1 taken 628 times.
✓ Branch 2 taken 5943 times.
✓ Branch 3 taken 3816 times.
 10387 if ( Link[j].type == CONDUIT && ntoks >= 8 )
260 {
261 5943 i = atoi(tok[7]);
262
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5943 times.
 5943 if ( i < 0 ) return error_setInpError(ERR_NUMBER, tok[7]);
263 5943 else Link[j].xsect.culvertCode = i;
264 }
265 }
266 10387 return 0;
267 }
268
269 //=============================================================================
270
271 184 int link_readLossParams(char* tok[], int ntoks)
272 //
273 // Input: tok[] = array of string tokens
274 // ntoks = number of tokens
275 // Output: returns an error code
276 // Purpose: reads local loss parameters for a link from a tokenized
277 // line of input data.
278 //
279 // Format: LinkID cInlet cOutlet cAvg FlapGate(YES/NO) SeepRate
280 //
281 {
282 int i, j, k;
283 double x[3];
284 184 double seepRate = 0.0;
285
286
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 184 times.
 184 if ( ntoks < 4 ) return error_setInpError(ERR_ITEMS, "");
287 184 j = project_findObject(LINK, tok[0]);
288
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 184 times.
 184 if ( j < 0 ) return error_setInpError(ERR_NAME, tok[0]);
289
2/2
✓ Branch 0 taken 552 times.
✓ Branch 1 taken 184 times.
 736 for (i=1; i<=3; i++)
290 {
291
2/4
✓ Branch 1 taken 552 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 552 times.
 552 if ( ! getDouble(tok[i], &x[i-1]) || x[i-1] < 0.0 )
292 ✗ return error_setInpError(ERR_NUMBER, tok[i]);
293 }
294 184 k = 0;
295
1/2
✓ Branch 0 taken 184 times.
✗ Branch 1 not taken.
 184 if ( ntoks >= 5 )
296 {
297 184 k = findmatch(tok[4], NoYesWords);
298
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 184 times.
 184 if ( k < 0 ) return error_setInpError(ERR_KEYWORD, tok[4]);
299 }
300
1/2
✓ Branch 0 taken 184 times.
✗ Branch 1 not taken.
 184 if ( ntoks >= 6 )
301 {
302
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 184 times.
 184 if ( ! getDouble(tok[5], &seepRate) )
303 ✗ return error_setInpError(ERR_NUMBER, tok[5]);
304 }
305 184 Link[j].cLossInlet = x[0];
306 184 Link[j].cLossOutlet = x[1];
307 184 Link[j].cLossAvg = x[2];
308 184 Link[j].hasFlapGate = k;
309 184 Link[j].seepRate = seepRate / UCF(RAINFALL);
310 184 return 0;
311 }
312
313 //=============================================================================
314
315 10508 void link_setParams(int j, int type, int n1, int n2, int k, double x[])
316 //
317 // Input: j = link index
318 // type = link type code
319 // n1 = index of upstream node
320 // n2 = index of downstream node
321 // k = index of link's sub-type
322 // x = array of parameter values
323 // Output: none
324 // Purpose: sets parameters for a link.
325 //
326 {
327 10508 Link[j].node1 = n1;
328 10508 Link[j].node2 = n2;
329 10508 Link[j].type = type;
330 10508 Link[j].subIndex = k;
331 10508 Link[j].offset1 = 0.0;
332 10508 Link[j].offset2 = 0.0;
333 10508 Link[j].q0 = 0.0;
334 10508 Link[j].qFull = 0.0;
335 10508 Link[j].setting = 1.0;
336 10508 Link[j].targetSetting = 1.0;
337 10508 Link[j].hasFlapGate = 0;
338 10508 Link[j].qLimit = 0.0; // 0 means that no limit is defined
339 10508 Link[j].direction = 1;
340
341
5/6
✓ Branch 0 taken 9772 times.
✓ Branch 1 taken 106 times.
✓ Branch 2 taken 159 times.
✓ Branch 3 taken 469 times.
✓ Branch 4 taken 2 times.
✗ Branch 5 not taken.
 10508 switch (type)
342 {
343 9772 case CONDUIT:
344 9772 Conduit[k].length = x[0] / UCF(LENGTH);
345 9772 Conduit[k].modLength = Conduit[k].length;
346 9772 Conduit[k].roughness = x[1];
347 9772 Link[j].offset1 = x[2] / UCF(LENGTH);
348 9772 Link[j].offset2 = x[3] / UCF(LENGTH);
349 9772 Link[j].q0 = x[4] / UCF(FLOW);
350 9772 Link[j].qLimit = x[5] / UCF(FLOW);
351 9772 break;
352
353 106 case PUMP:
354 106 Pump[k].pumpCurve = (int)x[0];
355 106 Link[j].hasFlapGate = FALSE;
356 106 Pump[k].initSetting = x[1];
357 106 Pump[k].yOn = x[2] / UCF(LENGTH);
358 106 Pump[k].yOff = x[3] / UCF(LENGTH);
359 106 Pump[k].xMin = 0.0;
360 106 Pump[k].xMax = 0.0;
361 106 break;
362
363 159 case ORIFICE:
364 159 Orifice[k].type = (int)x[0];
365 159 Link[j].offset1 = x[1] / UCF(LENGTH);
366 159 Link[j].offset2 = Link[j].offset1;
367 159 Orifice[k].cDisch = x[2];
368 159 Link[j].hasFlapGate = (x[3] > 0.0) ? 1 : 0;
369 159 Orifice[k].orate = x[4] * 3600.0;
370 159 break;
371
372 469 case WEIR:
373 469 Weir[k].type = (int)x[0];
374 469 Link[j].offset1 = x[1] / UCF(LENGTH);
375 469 Link[j].offset2 = Link[j].offset1;
376 469 Weir[k].cDisch1 = x[2];
377 469 Link[j].hasFlapGate = (x[3] > 0.0) ? 1 : 0;
378 469 Weir[k].endCon = x[4];
379 469 Weir[k].cDisch2 = x[5];
380 469 Weir[k].canSurcharge = (int)x[6];
381 469 Weir[k].roadWidth = x[7] / UCF(LENGTH);
382 469 Weir[k].roadSurface = (int)x[8];
383 469 Weir[k].cdCurve = (int)x[9];
384 469 break;
385
386 2 case OUTLET:
387 2 Link[j].offset1 = x[0] / UCF(LENGTH);
388 2 Link[j].offset2 = Link[j].offset1;
389 2 Outlet[k].qCoeff = x[1];
390 2 Outlet[k].qExpon = x[2];
391 2 Outlet[k].qCurve = (int)x[3];
392 2 Link[j].hasFlapGate = (x[4] > 0.0) ? 1 : 0;
393 2 Outlet[k].curveType = (int)x[5];
394
395 2 xsect_setParams(&Link[j].xsect, DUMMY, NULL, 0.0);
396 2 break;
397
398 }
399 10508 }
400
401 //=============================================================================
402
403 10508 void link_validate(int j)
404 //
405 // Input: j = link index
406 // Output: none
407 // Purpose: validates a link's properties.
408 //
409 {
410 int n;
411
412
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 10503 times.
 10508 if ( LinkOffsets == ELEV_OFFSET ) link_convertOffsets(j);
413
5/5
✓ Branch 0 taken 9772 times.
✓ Branch 1 taken 106 times.
✓ Branch 2 taken 159 times.
✓ Branch 3 taken 469 times.
✓ Branch 4 taken 2 times.
 10508 switch ( Link[j].type )
414 {
415 9772 case CONDUIT: conduit_validate(j, Link[j].subIndex); break;
416 106 case PUMP: pump_validate(j, Link[j].subIndex); break;
417 159 case ORIFICE: orifice_validate(j, Link[j].subIndex); break;
418 469 case WEIR: weir_validate(j, Link[j].subIndex); break;
419 }
420
421 // --- check if crest of regulator opening < invert of downstream node
422
2/2
✓ Branch 0 taken 630 times.
✓ Branch 1 taken 9878 times.
 10508 switch ( Link[j].type )
423 {
424 630 case ORIFICE:
425 case WEIR:
426 case OUTLET:
427 630 if ( Node[Link[j].node1].invertElev + Link[j].offset1 <
428
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 630 times.
 630 Node[Link[j].node2].invertElev )
429 {
430 ✗ if (RouteModel == DW)
431 {
432 ✗ Link[j].offset1 = Node[Link[j].node2].invertElev -
433 ✗ Node[Link[j].node1].invertElev;
434 ✗ report_writeWarningMsg(WARN10b, Link[j].ID);
435 }
436 ✗ else report_writeWarningMsg(WARN10a, Link[j].ID);
437 }
438 }
439
440 // --- force max. depth of end nodes to be >= link crown height
441 // at non-storage nodes
442
443 // --- skip pumps and bottom orifices
444
2/2
✓ Branch 0 taken 10402 times.
✓ Branch 1 taken 106 times.
 10508 if ( Link[j].type == PUMP ||
445
2/2
✓ Branch 0 taken 159 times.
✓ Branch 1 taken 10243 times.
 10402 (Link[j].type == ORIFICE &&
446
2/2
✓ Branch 0 taken 14 times.
✓ Branch 1 taken 145 times.
 265 Orifice[Link[j].subIndex].type == BOTTOM_ORIFICE) ) return;
447
448 // --- extend upstream node's full depth to link's crown elevation
449 10388 n = Link[j].node1;
450
3/4
✓ Branch 0 taken 61 times.
✓ Branch 1 taken 10327 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 61 times.
 10388 if ( Node[n].type != STORAGE || Node[n].surDepth > 0.0 )
451 {
452
2/2
✓ Branch 0 taken 10292 times.
✓ Branch 1 taken 35 times.
 10327 Node[n].fullDepth = MAX(Node[n].fullDepth,
453 Link[j].offset1 + Link[j].xsect.yFull);
454 }
455
456 // --- do same for downstream node only for conduit links
457 10388 n = Link[j].node2;
458
3/4
✓ Branch 0 taken 82 times.
✓ Branch 1 taken 10306 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 82 times.
 10388 if ( (Node[n].type != STORAGE || Node[n].surDepth > 0.0) &&
459
2/2
✓ Branch 0 taken 9718 times.
✓ Branch 1 taken 588 times.
 10306 Link[j].type == CONDUIT )
460 {
461
2/2
✓ Branch 0 taken 9312 times.
✓ Branch 1 taken 406 times.
 9718 Node[n].fullDepth = MAX(Node[n].fullDepth,
462 Link[j].offset2 + Link[j].xsect.yFull);
463 }
464 }
465
466 //=============================================================================
467
468 5 void link_convertOffsets(int j)
469 //
470 // Input: j = link index
471 // Output: none
472 // Purpose: converts offset elevations to offset heights for a link.
473 //
474 {
475 double elev;
476
477 5 elev = Node[Link[j].node1].invertElev;
478 5 Link[j].offset1 = link_getOffsetHeight(j, Link[j].offset1, elev);
479
1/2
✓ Branch 0 taken 5 times.
✗ Branch 1 not taken.
 5 if ( Link[j].type == CONDUIT )
480 {
481 5 elev = Node[Link[j].node2].invertElev;
482 5 Link[j].offset2 = link_getOffsetHeight(j, Link[j].offset2, elev);
483 }
484 ✗ else Link[j].offset2 = Link[j].offset1;
485 5 }
486
487 //=============================================================================
488
489 10 double link_getOffsetHeight(int j, double offset, double elev)
490 //
491 // Input: j = link index
492 // offset = link elevation offset (ft)
493 // elev = node invert elevation (ft)
494 // Output: returns offset distance above node invert (ft)
495 // Purpose: finds offset height for one end of a link.
496 //
497 {
498
2/4
✓ Branch 0 taken 10 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 10 times.
 10 if ( offset <= MISSING || Link[j].type == PUMP) return 0.0;
499 10 offset -= elev;
500
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 10 if ( offset >= 0.0 ) return offset;
501
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
 10 if ( offset >= -MIN_DELTA_Z ) return 0.0;
502 10 report_writeWarningMsg(WARN03, Link[j].ID);
503 10 return 0.0;
504 }
505
506 //=============================================================================
507
508 10508 void link_initState(int j)
509 //
510 // Input: j = link index
511 // Output: none
512 // Purpose: initializes a link's state variables at start of simulation.
513 //
514 {
515 int p;
516
517 // --- initialize hydraulic state
518 10508 Link[j].oldFlow = Link[j].q0;
519 10508 Link[j].newFlow = Link[j].q0;
520 10508 Link[j].oldDepth = 0.0;
521 10508 Link[j].newDepth = 0.0;
522 10508 Link[j].oldVolume = 0.0;
523 10508 Link[j].newVolume = 0.0;
524 10508 Link[j].setting = 1.0;
525 10508 Link[j].targetSetting = 1.0;
526 10508 Link[j].timeLastSet = StartDate;
527 10508 Link[j].inletControl = FALSE;
528 10508 Link[j].normalFlow = FALSE;
529
2/2
✓ Branch 0 taken 9772 times.
✓ Branch 1 taken 736 times.
 10508 if ( Link[j].type == CONDUIT ) conduit_initState(j, Link[j].subIndex);
530
2/2
✓ Branch 0 taken 106 times.
✓ Branch 1 taken 10402 times.
 10508 if ( Link[j].type == PUMP ) pump_initState(j, Link[j].subIndex);
531
532 // --- initialize water quality state
533
2/2
✓ Branch 0 taken 27182 times.
✓ Branch 1 taken 10508 times.
 37690 for (p = 0; p < Nobjects[POLLUT]; p++)
534 {
535 27182 Link[j].oldQual[p] = 0.0;
536 27182 Link[j].newQual[p] = 0.0;
537 27182 Link[j].totalLoad[p] = 0.0;
538 }
539 10508 }
540
541 //=============================================================================
542
543 12925100 double link_getInflow(int j)
544 //
545 // Input: j = link index
546 // Output: returns link flow rate (cfs)
547 // Purpose: finds total flow entering a link during current time step.
548 //
549 {
550
2/2
✓ Branch 0 taken 31329 times.
✓ Branch 1 taken 12893771 times.
 12925100 if ( Link[j].setting == 0 ) return 0.0;
551
5/6
✓ Branch 0 taken 11998612 times.
✓ Branch 1 taken 141238 times.
✓ Branch 2 taken 154967 times.
✓ Branch 3 taken 541350 times.
✓ Branch 4 taken 57604 times.
✗ Branch 5 not taken.
 12893771 switch ( Link[j].type )
552 {
553 11998612 case CONDUIT: return conduit_getInflow(j);
554 141238 case PUMP: return pump_getInflow(j);
555 154967 case ORIFICE: return orifice_getInflow(j);
556 541350 case WEIR: return weir_getInflow(j);
557 57604 case OUTLET: return outlet_getInflow(j);
558 ✗ default: return node_getOutflow(Link[j].node1, j);
559 }
560 }
561
562 //=============================================================================
563
564 37208897 void link_setOldHydState(int j)
565 //
566 // Input: j = link index
567 // Output: none
568 // Purpose: replaces link's old hydraulic state values with current ones.
569 //
570 {
571 int k;
572
573 37208897 Link[j].oldDepth = Link[j].newDepth;
574 37208897 Link[j].oldFlow = Link[j].newFlow;
575 37208897 Link[j].oldVolume = Link[j].newVolume;
576
577
2/2
✓ Branch 0 taken 36789022 times.
✓ Branch 1 taken 419875 times.
 37208897 if ( Link[j].type == CONDUIT )
578 {
579 36789022 k = Link[j].subIndex;
580 36789022 Conduit[k].q1Old = Conduit[k].q1;
581 36789022 Conduit[k].q2Old = Conduit[k].q2;
582 }
583 37208897 }
584
585 //=============================================================================
586
587 21771441 void link_setOldQualState(int j)
588 //
589 // Input: j = link index
590 // Output: none
591 // Purpose: replaces link's old water quality state values with current ones.
592 //
593 {
594 int p;
595
2/2
✓ Branch 0 taken 81872039 times.
✓ Branch 1 taken 21771441 times.
 103643480 for (p = 0; p < Nobjects[POLLUT]; p++)
596 {
597 81872039 Link[j].oldQual[p] = Link[j].newQual[p];
598 81872039 Link[j].newQual[p] = 0.0;
599 }
600 21771441 }
601
602 //=============================================================================
603
604 56222688 void link_setTargetSetting(int j)
605 //
606 // Input: j = link index
607 // Output: none
608 // Purpose: updates a link's target setting.
609 //
610 {
611 int k, n1;
612
2/2
✓ Branch 0 taken 337486 times.
✓ Branch 1 taken 55885202 times.
 56222688 if ( Link[j].type == PUMP )
613 {
614 337486 k = Link[j].subIndex;
615 337486 n1 = Link[j].node1;
616 337486 Link[j].targetSetting = Link[j].setting;
617
2/2
✓ Branch 0 taken 259957 times.
✓ Branch 1 taken 77529 times.
 337486 if ( Pump[k].yOff > 0.0 &&
618
2/2
✓ Branch 0 taken 12300 times.
✓ Branch 1 taken 247657 times.
 259957 Link[j].setting > 0.0 &&
619
2/2
✓ Branch 0 taken 60 times.
✓ Branch 1 taken 12240 times.
 12300 Node[n1].newDepth < Pump[k].yOff ) Link[j].targetSetting = 0.0;
620
2/2
✓ Branch 0 taken 279820 times.
✓ Branch 1 taken 57666 times.
 337486 if ( Pump[k].yOn > 0.0 &&
621
2/2
✓ Branch 0 taken 247657 times.
✓ Branch 1 taken 32163 times.
 279820 Link[j].setting == 0.0 &&
622
2/2
✓ Branch 0 taken 11 times.
✓ Branch 1 taken 247646 times.
 247657 Node[n1].newDepth > Pump[k].yOn ) Link[j].targetSetting = 1.0;
623 }
624 56222688 }
625
626 //=============================================================================
627
628 1069 void link_setSetting(int j, double tstep)
629 //
630 // Input: j = link index
631 // tstep = time step over which setting is adjusted
632 // Output: none
633 // Purpose: updates a link's setting as a result of a control action.
634 //
635 {
636
2/2
✓ Branch 0 taken 39 times.
✓ Branch 1 taken 1030 times.
 1069 if ( Link[j].type == ORIFICE ) orifice_setSetting(j, tstep);
637
2/2
✓ Branch 0 taken 78 times.
✓ Branch 1 taken 952 times.
 1030 else if ( Link[j].type == WEIR ) weir_setSetting(j);
638 952 else Link[j].setting = Link[j].targetSetting;
639 1069 }
640
641 //=============================================================================
642
643 48680957 int link_setFlapGate(int j, int n1, int n2, double q)
644 //
645 // Input: j = link index
646 // n1 = index of node on upstream end of link
647 // n2 = index of node on downstream end of link
648 // q = signed flow value (value and units don't matter)
649 // Output: returns TRUE if there is reverse flow through a flap gate
650 // associated with the link.
651 // Purpose: based on the sign of the flow, determines if a flap gate
652 // associated with the link should close or not.
653 //
654 {
655 48680957 int n = -1;
656
657 // --- check for reverse flow through link's flap gate
658
2/2
✓ Branch 0 taken 168894 times.
✓ Branch 1 taken 48512063 times.
 48680957 if ( Link[j].hasFlapGate )
659 {
660
2/2
✓ Branch 0 taken 15989 times.
✓ Branch 1 taken 152905 times.
 168894 if ( q * (double)Link[j].direction < 0.0 ) return TRUE;
661 }
662
663 // --- check for Outfall with flap gate node on inflow end of link
664
2/2
✓ Branch 0 taken 2446160 times.
✓ Branch 1 taken 46218808 times.
 48664968 if ( q < 0.0 ) n = n2;
665
2/2
✓ Branch 0 taken 46218802 times.
✓ Branch 1 taken 2446166 times.
 48664968 if ( q > 0.0 ) n = n1;
666
2/2
✓ Branch 0 taken 48664962 times.
✓ Branch 1 taken 6 times.
 48664968 if ( n >= 0 &&
667
2/2
✓ Branch 0 taken 17862 times.
✓ Branch 1 taken 48647100 times.
 48664962 Node[n].type == OUTFALL &&
668
2/2
✓ Branch 0 taken 15661 times.
✓ Branch 1 taken 2201 times.
 17862 Outfall[Node[n].subIndex].hasFlapGate ) return TRUE;
669 48649307 return FALSE;
670 }
671
672 //=============================================================================
673
674 1498302 void link_getResults(int j, double f, float x[])
675 //
676 // Input: j = link index
677 // f = time weighting factor
678 // Output: x = array of weighted results
679 // Purpose: retrieves time-weighted average of old and new results for a link.
680 //
681 {
682 int p; // pollutant index
683 double y, // depth
684 q, // flow
685 u, // velocity
686 v, // volume
687 c; // capacity, setting or concentration
688 1498302 double f1 = 1.0 - f;
689
690 1498302 y = f1*Link[j].oldDepth + f*Link[j].newDepth;
691 1498302 q = f1*Link[j].oldFlow + f*Link[j].newFlow;
692 1498302 v = f1*Link[j].oldVolume + f*Link[j].newVolume;
693 1498302 u = link_getVelocity(j, q, y);
694 1498302 c = 0.0;
695
2/2
✓ Branch 0 taken 1478366 times.
✓ Branch 1 taken 19936 times.
 1498302 if (Link[j].type == CONDUIT)
696 {
697
1/2
✓ Branch 0 taken 1478366 times.
✗ Branch 1 not taken.
 1478366 if (Link[j].xsect.type != DUMMY)
698 1478366 c = xsect_getAofY(&Link[j].xsect, y) / Link[j].xsect.aFull;
699 }
700 19936 else c = Link[j].setting;
701
702 // --- override time weighting for pump flow between on/off states
703
4/4
✓ Branch 0 taken 2340 times.
✓ Branch 1 taken 1495962 times.
✓ Branch 2 taken 1021 times.
✓ Branch 3 taken 1319 times.
 1498302 if (Link[j].type == PUMP && Link[j].oldFlow*Link[j].newFlow == 0.0)
704 {
705
1/2
✓ Branch 0 taken 1021 times.
✗ Branch 1 not taken.
 1021 if ( f >= f1 ) q = Link[j].newFlow;
706 ✗ else q = Link[j].oldFlow;
707 }
708
709 1498302 y *= UCF(LENGTH);
710 1498302 v *= UCF(VOLUME);
711 1498302 q *= UCF(FLOW) * (double)Link[j].direction;
712 1498302 u *= UCF(LENGTH) * (double)Link[j].direction;
713 1498302 x[LINK_DEPTH] = (float)y;
714 1498302 x[LINK_FLOW] = (float)q;
715 1498302 x[LINK_VELOCITY] = (float)u;
716 1498302 x[LINK_VOLUME] = (float)v;
717 1498302 x[LINK_CAPACITY] = (float)c;
718
719
4/4
✓ Branch 0 taken 1494846 times.
✓ Branch 1 taken 3456 times.
✓ Branch 2 taken 336445 times.
✓ Branch 3 taken 1494846 times.
 1834747 if ( !IgnoreQuality ) for (p = 0; p < Nobjects[POLLUT]; p++)
720 {
721 336445 c = f1*Link[j].oldQual[p] + f*Link[j].newQual[p];
722 336445 x[LINK_QUAL+p] = (float)c;
723 }
724 1498302 }
725
726 //=============================================================================
727
728 67130808 void link_setOutfallDepth(int j)
729 //
730 // Input: j = link index
731 // Output: none
732 // Purpose: sets depth at outfall node connected to link j.
733 //
734 {
735 int k; // conduit index
736 int n; // outfall node index
737 double z; // invert offset height (ft)
738 double q; // flow rate (cfs)
739 67130808 double yCrit = 0.0; // critical flow depth (ft)
740 67130808 double yNorm = 0.0; // normal flow depth (ft)
741
742 // --- find which end node of link is an outfall
743
2/2
✓ Branch 0 taken 1977695 times.
✓ Branch 1 taken 65153113 times.
 67130808 if ( Node[Link[j].node2].type == OUTFALL )
744 {
745 1977695 n = Link[j].node2;
746 1977695 z = Link[j].offset2;
747 }
748
2/2
✓ Branch 0 taken 1553 times.
✓ Branch 1 taken 65151560 times.
 65153113 else if ( Node[Link[j].node1].type == OUTFALL )
749 {
750 1553 n = Link[j].node1;
751 1553 z = Link[j].offset1;
752 }
753 65151560 else return;
754
755 // --- find both normal & critical depth for current flow
756
2/2
✓ Branch 0 taken 1952133 times.
✓ Branch 1 taken 27115 times.
 1979248 if ( Link[j].type == CONDUIT )
757 {
758 1952133 k = Link[j].subIndex;
759 1952133 q = fabs(Link[j].newFlow / Conduit[k].barrels);
760 1952133 yNorm = link_getYnorm(j, q);
761 1952133 yCrit = link_getYcrit(j, q);
762 }
763
764 // --- set new depth at node
765 1979248 node_setOutletDepth(n, yNorm, yCrit, z);
766 }
767
768 //=============================================================================
769
770 3763349 double link_getYcrit(int j, double q)
771 //
772 // Input: j = link index
773 // q = link flow rate (cfs)
774 // Output: returns critical depth (ft)
775 // Purpose: computes critical depth for given flow rate.
776 //
777 {
778 3763349 return xsect_getYcrit(&Link[j].xsect, q);
779 }
780
781 //=============================================================================
782
783 3773121 double link_getYnorm(int j, double q)
784 //
785 // Input: j = link index
786 // q = link flow rate (cfs)
787 // Output: returns normal depth (ft)
788 // Purpose: computes normal depth for given flow rate.
789 //
790 {
791 int k;
792 double s, a, y;
793
794
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3773121 times.
 3773121 if ( Link[j].type != CONDUIT ) return 0.0;
795
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3773121 times.
 3773121 if ( Link[j].xsect.type == DUMMY ) return 0.0;
796 3773121 q = fabs(q);
797 3773121 k = Link[j].subIndex;
798
2/2
✓ Branch 0 taken 178889 times.
✓ Branch 1 taken 3594232 times.
 3773121 if ( q > Conduit[k].qMax ) q = Conduit[k].qMax;
799
2/2
✓ Branch 0 taken 384342 times.
✓ Branch 1 taken 3388779 times.
 3773121 if ( q <= 0.0 ) return 0.0;
800 3388779 s = q / Conduit[k].beta;
801 3388779 a = xsect_getAofS(&Link[j].xsect, s);
802 3388779 y = xsect_getYofA(&Link[j].xsect, a);
803 3388779 return y;
804 }
805
806 //=============================================================================
807
808 132087011 double link_getLength(int j)
809 //
810 // Input: j = link index
811 // Output: returns length (ft)
812 // Purpose: finds true length of a link.
813 //
814 {
815
1/2
✓ Branch 0 taken 132087011 times.
✗ Branch 1 not taken.
 132087011 if ( Link[j].type == CONDUIT ) return conduit_getLength(j);
816 ✗ return 0.0;
817 }
818
819 //=============================================================================
820
821 38707199 double link_getVelocity(int j, double flow, double depth)
822 //
823 // Input: j = link index
824 // flow = link flow rate (cfs)
825 // depth = link flow depth (ft)
826 // Output: returns flow velocity (fps)
827 // Purpose: finds flow velocity given flow and depth.
828 //
829 {
830 double area;
831 38707199 double veloc = 0.0;
832 int k;
833
834
2/2
✓ Branch 0 taken 14456668 times.
✓ Branch 1 taken 24250531 times.
 38707199 if ( depth <= 0.01 ) return 0.0;
835
2/2
✓ Branch 0 taken 24075065 times.
✓ Branch 1 taken 175466 times.
 24250531 if ( Link[j].type == CONDUIT )
836 {
837 24075065 k = Link[j].subIndex;
838 24075065 flow /= Conduit[k].barrels;
839 24075065 area = xsect_getAofY(&Link[j].xsect, depth);
840
2/2
✓ Branch 0 taken 24075056 times.
✓ Branch 1 taken 9 times.
 24075065 if (area > FUDGE ) veloc = flow / area;
841 }
842 24250531 return veloc;
843 }
844
845 //=============================================================================
846
847 63452045 double link_getFroude(int j, double v, double y)
848 //
849 // Input: j = link index
850 // v = flow velocity (fps)
851 // y = flow depth (ft)
852 // Output: returns Froude Number
853 // Purpose: computes Froude Number for given velocity and flow depth
854 //
855 {
856 63452045 TXsect* xsect = &Link[j].xsect;
857
858 // --- return 0 if link is not a conduit
859
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 63452045 times.
 63452045 if ( Link[j].type != CONDUIT ) return 0.0;
860
861 // --- return 0 if link empty or closed conduit is full
862
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 63452045 times.
 63452045 if ( y <= FUDGE ) return 0.0;
863
2/2
✓ Branch 1 taken 61417136 times.
✓ Branch 2 taken 2034909 times.
 63452045 if ( !xsect_isOpen(xsect->type) &&
864
2/2
✓ Branch 0 taken 738333 times.
✓ Branch 1 taken 60678803 times.
 61417136 xsect->yFull - y <= FUDGE ) return 0.0;
865
866 // --- compute hydraulic depth
867 62713712 y = xsect_getAofY(xsect, y) / xsect_getWofY(xsect, y);
868
869 // --- compute Froude No.
870 62713712 return fabs(v) / sqrt(GRAVITY * y);
871 }
872
873 //=============================================================================
874
875 70325 double link_getPower(int j)
876 //
877 // Input: j = link index
878 // Output: returns power consumed by link in kwatts
879 // Purpose: computes power consumed by head loss (or head gain) of
880 // water flowing through a link
881 //
882 {
883 70325 int n1 = Link[j].node1;
884 70325 int n2 = Link[j].node2;
885 70325 double dh = (Node[n1].invertElev + Node[n1].newDepth) -
886 70325 (Node[n2].invertElev + Node[n2].newDepth);
887 70325 double q = fabs(Link[j].newFlow);
888 70325 return fabs(dh) * q / 8.814 * KWperHP;
889 }
890
891 //=============================================================================
892
893 60241010 double link_getLossRate(int j, int routeModel, double q, double tstep)
894 //
895 // Input: j = link index
896 // routeModel = flow routing model type
897 // q = flow rate (ft3/sec)
898 // tstep = time step (sec)
899 // Output: returns uniform loss rate in link (ft3/sec)
900 // Purpose: computes rate at which flow volume is lost in a link due to
901 // evaporation and seepage.
902 //
903 {
904
1/2
✓ Branch 0 taken 60241010 times.
✗ Branch 1 not taken.
 60241010 if ( Link[j].type == CONDUIT )
905 60241010 return conduit_getLossRate(j, routeModel, q, tstep);
906 ✗ else return 0.0;
907 }
908
909 //=============================================================================
910
911 48346631 char link_getFullState(double a1, double a2, double aFull)
912 //
913 // Input: a1 = upstream link area (ft2)
914 // a2 = downstream link area (ft2)
915 // aFull = area of full conduit
916 // Output: returns fullness state of a link
917 // Purpose: determines if a link is upstream, downstream or completely full.
918 //
919 {
920
2/2
✓ Branch 0 taken 1100919 times.
✓ Branch 1 taken 47245712 times.
 48346631 if ( a1 >= aFull )
921 {
922
2/2
✓ Branch 0 taken 860924 times.
✓ Branch 1 taken 239995 times.
 1100919 if ( a2 >= aFull ) return ALL_FULL;
923 239995 else return UP_FULL;
924 }
925
2/2
✓ Branch 0 taken 1737857 times.
✓ Branch 1 taken 45507855 times.
 47245712 if ( a2 >= aFull ) return DN_FULL;
926 45507855 return 0;
927 }
928
929 //=============================================================================
930 // C O N D U I T M E T H O D S
931 //=============================================================================
932
933 9772 int conduit_readParams(int j, int k, char* tok[], int ntoks)
934 //
935 // Input: j = link index
936 // k = conduit index
937 // tok[] = array of string tokens
938 // ntoks = number of tokens
939 // Output: returns an error code
940 // Purpose: reads conduit parameters from a tokenzed line of input.
941 //
942 {
943 int n1, n2;
944 double x[6];
945 char* id;
946
947 // --- check for valid ID and end node IDs
948
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9772 times.
 9772 if ( ntoks < 7 ) return error_setInpError(ERR_ITEMS, "");
949 9772 id = project_findID(LINK, tok[0]); // link ID
950
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9772 times.
 9772 if ( id == NULL ) return error_setInpError(ERR_NAME, tok[0]);
951 9772 n1 = project_findObject(NODE, tok[1]); // upstrm. node
952
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9772 times.
 9772 if ( n1 < 0 ) return error_setInpError(ERR_NAME, tok[1]);
953 9772 n2 = project_findObject(NODE, tok[2]); // dwnstrm. node
954
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9772 times.
 9772 if ( n2 < 0 ) return error_setInpError(ERR_NAME, tok[2]);
955
956 // --- parse length & Mannings N
957
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 9772 times.
 9772 if ( !getDouble(tok[3], &x[0]) )
958 ✗ return error_setInpError(ERR_NUMBER, tok[3]);
959
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 9772 times.
 9772 if ( !getDouble(tok[4], &x[1]) )
960 ✗ return error_setInpError(ERR_NUMBER, tok[4]);
961
962 // --- parse offsets
963
3/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 9767 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5 times.
 9772 if ( LinkOffsets == ELEV_OFFSET && *tok[5] == '*' ) x[2] = MISSING;
964
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 9772 times.
 9772 else if ( !getDouble(tok[5], &x[2]) )
965 ✗ return error_setInpError(ERR_NUMBER, tok[5]);
966
3/4
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 9767 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5 times.
 9772 if ( LinkOffsets == ELEV_OFFSET && *tok[6] == '*' ) x[3] = MISSING;
967
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 9772 times.
 9772 else if ( !getDouble(tok[6], &x[3]) )
968 ✗ return error_setInpError(ERR_NUMBER, tok[6]);
969
970 // --- parse optional parameters
971 9772 x[4] = 0.0; // init. flow
972
1/2
✓ Branch 0 taken 9772 times.
✗ Branch 1 not taken.
 9772 if ( ntoks >= 8 )
973 {
974
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 9772 times.
 9772 if ( !getDouble(tok[7], &x[4]) )
975 ✗ return error_setInpError(ERR_NUMBER, tok[7]);
976 }
977 9772 x[5] = 0.0;
978
1/2
✓ Branch 0 taken 9772 times.
✗ Branch 1 not taken.
 9772 if ( ntoks >= 9 )
979 {
980
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 9772 times.
 9772 if ( !getDouble(tok[8], &x[5]) )
981 ✗ return error_setInpError(ERR_NUMBER, tok[8]);
982 }
983
984 // --- add parameters to data base
985 9772 Link[j].ID = id;
986 9772 link_setParams(j, CONDUIT, n1, n2, k, x);
987 9772 return 0;
988 }
989
990 //=============================================================================
991
992 9772 void conduit_validate(int j, int k)
993 //
994 // Input: j = link index
995 // k = conduit index
996 // Output: none
997 // Purpose: validates a conduit's properties.
998 //
999 {
1000 double aa;
1001 double lengthFactor, roughness, slope;
1002
1003 // --- a storage node cannot have a dummy outflow link
1004
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 9772 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 9772 if ( Link[j].xsect.type == DUMMY && RouteModel == DW )
1005 {
1006 ✗ if ( Node[Link[j].node1].type == STORAGE )
1007 {
1008 ✗ report_writeErrorMsg(ERR_DUMMY_LINK, Node[Link[j].node1].ID);
1009 ✗ return;
1010 }
1011 }
1012
1013 // --- if custom xsection, then set its parameters
1014
2/2
✓ Branch 0 taken 10 times.
✓ Branch 1 taken 9762 times.
 9772 if ( Link[j].xsect.type == CUSTOM )
1015 10 xsect_setCustomXsectParams(&Link[j].xsect);
1016
1017 // --- if irreg. xsection, assign transect roughness to conduit
1018
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 9771 times.
 9772 if ( Link[j].xsect.type == IRREGULAR )
1019 {
1020 1 xsect_setIrregXsectParams(&Link[j].xsect);
1021 1 Conduit[k].roughness = Transect[Link[j].xsect.transect].roughness;
1022 }
1023
1024 // --- if street xsection, then set its parameters
1025
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 9760 times.
 9772 if (Link[j].xsect.type == STREET_XSECT)
1026 {
1027 12 xsect_setStreetXsectParams(&Link[j].xsect);
1028 12 Conduit[k].roughness = Street[Link[j].xsect.transect].roughness;
1029 }
1030
1031 // --- if force main xsection, adjust units on D-W roughness height
1032
2/2
✓ Branch 0 taken 134 times.
✓ Branch 1 taken 9638 times.
 9772 if ( Link[j].xsect.type == FORCE_MAIN )
1033 {
1034
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 126 times.
 134 if ( ForceMainEqn == D_W ) Link[j].xsect.rBot /= UCF(RAINDEPTH);
1035
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 134 times.
 134 if ( Link[j].xsect.rBot <= 0.0 )
1036 ✗ report_writeErrorMsg(ERR_XSECT, Link[j].ID);
1037 }
1038
1039 // --- check for valid length & roughness
1040
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9772 times.
 9772 if ( Conduit[k].length <= 0.0 )
1041 ✗ report_writeErrorMsg(ERR_LENGTH, Link[j].ID);
1042
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9772 times.
 9772 if ( Conduit[k].roughness <= 0.0 )
1043 ✗ report_writeErrorMsg(ERR_ROUGHNESS, Link[j].ID);
1044
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9772 times.
 9772 if ( Conduit[k].barrels <= 0 )
1045 ✗ report_writeErrorMsg(ERR_BARRELS, Link[j].ID);
1046
1047 // --- check for valid xsection
1048
1/2
✓ Branch 0 taken 9772 times.
✗ Branch 1 not taken.
 9772 if ( Link[j].xsect.type != DUMMY )
1049 {
1050
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9772 times.
 9772 if ( Link[j].xsect.type < 0 )
1051 ✗ report_writeErrorMsg(ERR_NO_XSECT, Link[j].ID);
1052
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9772 times.
 9772 else if ( Link[j].xsect.aFull <= 0.0 )
1053 ✗ report_writeErrorMsg(ERR_XSECT, Link[j].ID);
1054 }
1055
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9772 times.
 9772 if ( ErrorCode ) return;
1056
1057 // --- check for negative offsets
1058
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 9771 times.
 9772 if ( Link[j].offset1 < 0.0 )
1059 {
1060 1 report_writeWarningMsg(WARN03, Link[j].ID);
1061 1 Link[j].offset1 = 0.0;
1062 }
1063
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 9771 times.
 9772 if ( Link[j].offset2 < 0.0 )
1064 {
1065 1 report_writeWarningMsg(WARN03, Link[j].ID);
1066 1 Link[j].offset2 = 0.0;
1067 }
1068
1069 // --- adjust conduit offsets for partly filled circular xsection
1070
2/2
✓ Branch 0 taken 175 times.
✓ Branch 1 taken 9597 times.
 9772 if ( Link[j].xsect.type == FILLED_CIRCULAR )
1071 {
1072 175 Link[j].offset1 += Link[j].xsect.yBot;
1073 175 Link[j].offset2 += Link[j].xsect.yBot;
1074 }
1075
1076 // --- compute conduit slope
1077 9772 slope = conduit_getSlope(j);
1078 9772 Conduit[k].slope = slope;
1079
1080 // --- reverse orientation of conduit if using dynamic wave routing
1081 // and slope is negative
1082
4/4
✓ Branch 0 taken 9718 times.
✓ Branch 1 taken 54 times.
✓ Branch 2 taken 467 times.
✓ Branch 3 taken 9251 times.
 9772 if ( RouteModel == DW &&
1083 467 slope < 0.0 &&
1084
1/2
✓ Branch 0 taken 467 times.
✗ Branch 1 not taken.
 467 Link[j].xsect.type != DUMMY )
1085 {
1086 467 conduit_reverse(j, k);
1087 }
1088
1089 // --- get equivalent Manning roughness for Force Mains
1090 // for use when pipe is partly full
1091 9772 roughness = Conduit[k].roughness;
1092
4/4
✓ Branch 0 taken 9718 times.
✓ Branch 1 taken 54 times.
✓ Branch 2 taken 134 times.
✓ Branch 3 taken 9584 times.
 9772 if ( RouteModel == DW && Link[j].xsect.type == FORCE_MAIN )
1093 {
1094 134 roughness = forcemain_getEquivN(j, k);
1095 }
1096
1097 // --- adjust roughness for meandering natural channels
1098
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 9771 times.
 9772 if ( Link[j].xsect.type == IRREGULAR )
1099 {
1100 1 lengthFactor = Transect[Link[j].xsect.transect].lengthFactor;
1101 1 roughness *= sqrt(lengthFactor);
1102 }
1103
1104 // --- lengthen conduit if lengthening option is in effect
1105 9772 lengthFactor = 1.0;
1106
2/2
✓ Branch 0 taken 9718 times.
✓ Branch 1 taken 54 times.
 9772 if ( RouteModel == DW &&
1107
2/2
✓ Branch 0 taken 9608 times.
✓ Branch 1 taken 110 times.
 9718 LengtheningStep > 0.0 &&
1108
1/2
✓ Branch 0 taken 9608 times.
✗ Branch 1 not taken.
 9608 Link[j].xsect.type != DUMMY )
1109 {
1110 9608 lengthFactor = conduit_getLengthFactor(j, k, roughness);
1111 }
1112
1113
2/2
✓ Branch 0 taken 2559 times.
✓ Branch 1 taken 7213 times.
 9772 if ( lengthFactor != 1.0 )
1114 {
1115 2559 Conduit[k].modLength = lengthFactor * conduit_getLength(j);
1116 2559 slope /= lengthFactor;
1117 2559 roughness = roughness / sqrt(lengthFactor);
1118 }
1119
1120 // --- compute roughness factor used when computing friction
1121 // slope term in Dynamic Wave flow routing
1122
1123 // --- special case for non-Manning Force Mains
1124 // (roughness factor for full flow is saved in xsect.sBot)
1125
4/4
✓ Branch 0 taken 9718 times.
✓ Branch 1 taken 54 times.
✓ Branch 2 taken 134 times.
✓ Branch 3 taken 9584 times.
 9772 if ( RouteModel == DW && Link[j].xsect.type == FORCE_MAIN )
1126 {
1127 134 Link[j].xsect.sBot =
1128 134 forcemain_getRoughFactor(j, lengthFactor);
1129 }
1130 9772 Conduit[k].roughFactor = GRAVITY * SQR(roughness/PHI);
1131
1132 // --- compute full flow through cross section
1133
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9772 times.
 9772 if ( Link[j].xsect.type == DUMMY ) Conduit[k].beta = 0.0;
1134 9772 else Conduit[k].beta = PHI * sqrt(fabs(slope)) / roughness;
1135 9772 Link[j].qFull = Link[j].xsect.sFull * Conduit[k].beta;
1136 9772 Conduit[k].qMax = Link[j].xsect.sMax * Conduit[k].beta;
1137
1138 // --- see if flow is supercritical most of time
1139 // by comparing normal & critical velocities.
1140 // (factor of 0.3 is for circular pipe 95% full)
1141 // NOTE: this factor was used in the past for a modified version of
1142 // Kinematic Wave routing but is now deprecated.
1143 9772 aa = Conduit[k].beta / sqrt(32.2) *
1144 9772 pow(Link[j].xsect.yFull, 0.1666667) * 0.3;
1145
2/2
✓ Branch 0 taken 1300 times.
✓ Branch 1 taken 8472 times.
 9772 if ( aa >= 1.0 ) Conduit[k].superCritical = TRUE;
1146 8472 else Conduit[k].superCritical = FALSE;
1147
1148 // --- set value of hasLosses flag
1149
2/2
✓ Branch 0 taken 9743 times.
✓ Branch 1 taken 29 times.
 9772 if ( Link[j].cLossInlet == 0.0 &&
1150
1/2
✓ Branch 0 taken 9743 times.
✗ Branch 1 not taken.
 9743 Link[j].cLossOutlet == 0.0 &&
1151
1/2
✓ Branch 0 taken 9743 times.
✗ Branch 1 not taken.
 9743 Link[j].cLossAvg == 0.0
1152 9743 ) Conduit[k].hasLosses = FALSE;
1153 29 else Conduit[k].hasLosses = TRUE;
1154 }
1155
1156 //=============================================================================
1157
1158 467 void conduit_reverse(int j, int k)
1159 //
1160 // Input: j = link index
1161 // k = conduit index
1162 // Output: none
1163 // Purpose: reverses direction of a conduit
1164 //
1165 {
1166 int i;
1167 double z;
1168 double cLoss;
1169
1170 // --- reverse end nodes
1171 467 i = Link[j].node1;
1172 467 Link[j].node1 = Link[j].node2;
1173 467 Link[j].node2 = i;
1174
1175 // --- reverse node offsets
1176 467 z = Link[j].offset1;
1177 467 Link[j].offset1 = Link[j].offset2;
1178 467 Link[j].offset2 = z;
1179
1180 // --- reverse loss coeffs.
1181 467 cLoss = Link[j].cLossInlet;
1182 467 Link[j].cLossInlet = Link[j].cLossOutlet;
1183 467 Link[j].cLossOutlet = cLoss;
1184
1185 // --- reverse direction & slope
1186 467 Conduit[k].slope = -Conduit[k].slope;
1187 467 Link[j].direction *= (signed char)-1;
1188
1189 // --- reverse initial flow value
1190 467 Link[j].q0 = -Link[j].q0;
1191 467 }
1192
1193 //=============================================================================
1194
1195 181034358 double conduit_getLength(int j)
1196 //
1197 // Input: j = link index
1198 // Output: returns conduit's length (ft)
1199 // Purpose: finds true length of a conduit.
1200 //
1201 // Note: for irregular natural channels, user inputs length of main
1202 // channel (for FEMA purposes) but program should use length
1203 // associated with entire flood plain. Transect.lengthFactor
1204 // is the ratio of these two lengths.
1205 //
1206 {
1207 181034358 int k = Link[j].subIndex;
1208 int t;
1209
2/2
✓ Branch 0 taken 181025603 times.
✓ Branch 1 taken 8755 times.
 181034358 if ( Link[j].xsect.type != IRREGULAR ) return Conduit[k].length;
1210 8755 t = Link[j].xsect.transect;
1211
2/4
✓ Branch 0 taken 8755 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 8755 times.
 8755 if ( t < 0 || t >= Nobjects[TRANSECT] ) return Conduit[k].length;
1212 8755 return Conduit[k].length / Transect[t].lengthFactor;
1213 }
1214
1215 //=============================================================================
1216
1217 9608 double conduit_getLengthFactor(int j, int k, double roughness)
1218 //
1219 // Input: j = link index
1220 // k = conduit index
1221 // roughness = conduit Manning's n
1222 // Output: returns factor by which a conduit should be lengthened
1223 // Purpose: computes amount of conduit lengthing to improve numerical stability.
1224 //
1225 // The following form of the Courant criterion is used:
1226 // L = t * v * (1 + Fr) / Fr
1227 // where L = conduit length, t = time step, v = velocity, & Fr = Froude No.
1228 // After substituting Fr = v / sqrt(gy), where y = flow depth, we get:
1229 // L = t * ( sqrt(gy) + v )
1230 //
1231 {
1232 double ratio;
1233 double yFull;
1234 double vFull;
1235 double tStep;
1236
1237 // --- evaluate flow depth and velocity at full normal flow condition
1238 9608 yFull = Link[j].xsect.yFull;
1239
2/2
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 9596 times.
 9608 if ( xsect_isOpen(Link[j].xsect.type) )
1240 {
1241 12 yFull = Link[j].xsect.aFull / xsect_getWofY(&Link[j].xsect, yFull);
1242 }
1243 9608 vFull = PHI / roughness * Link[j].xsect.sFull *
1244 9608 sqrt(fabs(Conduit[k].slope)) / Link[j].xsect.aFull;
1245
1246 // --- determine ratio of Courant length to actual length
1247
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9608 times.
 9608 if ( LengtheningStep == 0.0 ) tStep = RouteStep;
1248
2/2
✓ Branch 0 taken 3546 times.
✓ Branch 1 taken 6062 times.
 9608 else tStep = MIN(RouteStep, LengtheningStep);
1249 9608 ratio = (sqrt(GRAVITY*yFull) + vFull) * tStep / conduit_getLength(j);
1250
1251 // --- return max. of 1.0 and ratio
1252
2/2
✓ Branch 0 taken 2559 times.
✓ Branch 1 taken 7049 times.
 9608 if ( ratio > 1.0 ) return ratio;
1253 7049 else return 1.0;
1254 }
1255
1256 //=============================================================================
1257
1258 9772 double conduit_getSlope(int j)
1259 //
1260 // Input: j = link index
1261 // Output: returns conduit slope
1262 // Purpose: computes conduit slope.
1263 //
1264 {
1265 double elev1, elev2, delta, slope;
1266 9772 double length = conduit_getLength(j);
1267
1268 // --- check that elevation drop > minimum allowable drop
1269 9772 elev1 = Link[j].offset1 + Node[Link[j].node1].invertElev;
1270 9772 elev2 = Link[j].offset2 + Node[Link[j].node2].invertElev;
1271 9772 delta = fabs(elev1 - elev2);
1272
2/2
✓ Branch 0 taken 23 times.
✓ Branch 1 taken 9749 times.
 9772 if ( delta < MIN_DELTA_Z )
1273 {
1274 23 report_writeWarningMsg(WARN04, Link[j].ID);
1275 23 delta = MIN_DELTA_Z;
1276 }
1277
1278 // --- elevation drop cannot exceed conduit length
1279
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 9771 times.
 9772 if ( delta >= length )
1280 {
1281 1 report_writeWarningMsg(WARN08, Link[j].ID);
1282 1 slope = delta / length;
1283 }
1284
1285 // --- slope = elev. drop / horizontal distance
1286 9771 else slope = delta / sqrt(SQR(length) - SQR(delta));
1287
1288 // -- check that slope exceeds minimum allowable slope
1289
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 9772 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 9772 if ( MinSlope > 0.0 && slope < MinSlope )
1290 {
1291 ✗ report_writeWarningMsg(WARN05, Link[j].ID);
1292 ✗ slope = MinSlope;
1293 // keep min. slope positive for SF or KW routing
1294 ✗ if (RouteModel == SF || RouteModel == KW) return slope;
1295 }
1296
1297 // --- change sign for adverse slope
1298
2/2
✓ Branch 0 taken 467 times.
✓ Branch 1 taken 9305 times.
 9772 if ( elev1 < elev2 ) slope = -slope;
1299 9772 return slope;
1300 }
1301
1302 //=============================================================================
1303
1304 9772 void conduit_initState(int j, int k)
1305 //
1306 // Input: j = link index
1307 // k = conduit index
1308 // Output: none
1309 // Purpose: sets initial conduit depth to normal depth of initial flow
1310 //
1311 {
1312 9772 Link[j].newDepth = link_getYnorm(j, Link[j].q0 / Conduit[k].barrels);
1313 9772 Link[j].oldDepth = Link[j].newDepth;
1314 9772 Conduit[k].evapLossRate = 0.0;
1315 9772 Conduit[k].seepLossRate = 0.0;
1316 9772 }
1317
1318 //=============================================================================
1319
1320 11998612 double conduit_getInflow(int j)
1321 //
1322 // Input: j = link index
1323 // Output: returns flow in link (cfs)
1324 // Purpose: finds inflow to conduit from upstream node.
1325 //
1326 {
1327 11998612 double qIn = node_getOutflow(Link[j].node1, j);
1328
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 11998612 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 11998612 if ( Link[j].qLimit > 0.0 ) qIn = MIN(qIn, Link[j].qLimit);
1329 11998612 return qIn;
1330 }
1331
1332 //=============================================================================
1333
1334 60241010 double conduit_getLossRate(int j, int routeModel, double q, double tstep)
1335 //
1336 // Input: j = link index
1337 // routeModel = type of flow routing model
1338 // q = current link flow rate (cfs)
1339 // tstep = current routing time step (sec)
1340 // Output: returns rate of evaporation & seepage losses (ft3/sec)
1341 // Purpose: computes volumetric rate of water evaporation & seepage
1342 // from a conduit (per barrel).
1343 //
1344 {
1345 TXsect *xsect;
1346 60241010 double depth = 0.5 * (Link[j].oldDepth + Link[j].newDepth);
1347 double length;
1348 double width, topWidth;
1349 60241010 double evapLossRate = 0.0,
1350 60241010 seepLossRate = 0.0,
1351 60241010 totalLossRate = 0.0;
1352
1353
2/2
✓ Branch 0 taken 48925408 times.
✓ Branch 1 taken 11315602 times.
 60241010 if ( depth > FUDGE )
1354 {
1355 48925408 xsect = &Link[j].xsect;
1356 48925408 length = conduit_getLength(j);
1357
1358 // --- find evaporation rate for open conduits
1359
4/4
✓ Branch 1 taken 1980976 times.
✓ Branch 2 taken 46944432 times.
✓ Branch 3 taken 1230897 times.
✓ Branch 4 taken 750079 times.
 48925408 if ( xsect_isOpen(xsect->type) && Evap.rate > 0.0 )
1360 {
1361 1230897 topWidth = xsect_getWofY(xsect, depth);
1362 1230897 evapLossRate = topWidth * length * Evap.rate;
1363 }
1364
1365 // --- compute seepage loss rate
1366
2/2
✓ Branch 0 taken 9763 times.
✓ Branch 1 taken 48915645 times.
 48925408 if ( Link[j].seepRate > 0.0 )
1367 {
1368 // limit depth to depth at max width
1369
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9763 times.
 9763 if (xsect->type == RECT_CLOSED) width = xsect->wMax;
1370 else
1371 {
1372
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 9763 times.
 9763 if ( depth >= xsect->ywMax ) depth = xsect->ywMax;
1373 9763 width = xsect_getWofY(xsect, depth);
1374 }
1375
1376 // compute seepage loss rate across length of conduit
1377 9763 seepLossRate = Link[j].seepRate * width * length;
1378 9763 seepLossRate *= Adjust.hydconFactor;
1379 }
1380
1381 // --- compute total loss rate
1382 48925408 totalLossRate = evapLossRate + seepLossRate;
1383
1384 // --- limit total loss rate to current volume for DW routing
1385 // or current link flow rate otherwise
1386
2/2
✓ Branch 0 taken 48249618 times.
✓ Branch 1 taken 675790 times.
 48925408 if (routeModel == DW) q = Link[j].newVolume / tstep;
1387
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 675790 times.
 675790 else q = ABS(q);
1388
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 48925408 times.
 48925408 if (totalLossRate > q)
1389 {
1390 ✗ evapLossRate = evapLossRate * q / totalLossRate;
1391 ✗ seepLossRate = seepLossRate * q / totalLossRate;
1392 ✗ totalLossRate = q;
1393 }
1394 }
1395
1396 60241010 Conduit[Link[j].subIndex].evapLossRate = evapLossRate;
1397 60241010 Conduit[Link[j].subIndex].seepLossRate = seepLossRate;
1398 60241010 return totalLossRate;
1399 }
1400
1401
1402 //=============================================================================
1403 // P U M P M E T H O D S
1404 //=============================================================================
1405
1406 106 int pump_readParams(int j, int k, char* tok[], int ntoks)
1407 //
1408 // Input: j = link index
1409 // k = pump index
1410 // tok[] = array of string tokens
1411 // ntoks = number of tokens
1412 // Output: returns an error code
1413 // Purpose: reads pump parameters from a tokenized line of input.
1414 //
1415 {
1416 int m;
1417 int n1, n2;
1418 double x[4];
1419 char* id;
1420
1421 // --- check for valid ID and end node IDs
1422
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 106 times.
 106 if ( ntoks < 3 ) return error_setInpError(ERR_ITEMS, "");
1423 106 id = project_findID(LINK, tok[0]);
1424
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 106 times.
 106 if ( id == NULL ) return error_setInpError(ERR_NAME, tok[0]);
1425 106 n1 = project_findObject(NODE, tok[1]);
1426
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 106 times.
 106 if ( n1 < 0 ) return error_setInpError(ERR_NAME, tok[1]);
1427 106 n2 = project_findObject(NODE, tok[2]);
1428
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 106 times.
 106 if ( n2 < 0 ) return error_setInpError(ERR_NAME, tok[2]);
1429
1430 // --- parse curve name
1431 106 x[0] = -1.;
1432
1/2
✓ Branch 0 taken 106 times.
✗ Branch 1 not taken.
 106 if ( ntoks >= 4 )
1433 {
1434
2/2
✓ Branch 1 taken 83 times.
✓ Branch 2 taken 23 times.
 106 if ( !strcomp(tok[3],"*") )
1435 {
1436 83 m = project_findObject(CURVE, tok[3]);
1437
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 83 times.
 83 if ( m < 0 ) return error_setInpError(ERR_NAME, tok[3]);
1438 83 x[0] = m;
1439 }
1440 }
1441
1442 // --- parse init. status if present
1443 106 x[1] = 1.0;
1444
1/2
✓ Branch 0 taken 106 times.
✗ Branch 1 not taken.
 106 if ( ntoks >= 5 )
1445 {
1446 106 m = findmatch(tok[4], OffOnWords);
1447
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 106 times.
 106 if ( m < 0 ) return error_setInpError(ERR_KEYWORD, tok[4]);
1448 106 x[1] = m;
1449 }
1450
1451 // --- parse startup/shutoff depths if present
1452 106 x[2] = 0.0;
1453
1/2
✓ Branch 0 taken 106 times.
✗ Branch 1 not taken.
 106 if ( ntoks >= 6 )
1454 {
1455
2/4
✓ Branch 1 taken 106 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 106 times.
 106 if ( !getDouble(tok[5], &x[2]) || x[2] < 0.0)
1456 ✗ return error_setInpError(ERR_NUMBER, tok[5]);
1457 }
1458 106 x[3] = 0.0;
1459
1/2
✓ Branch 0 taken 106 times.
✗ Branch 1 not taken.
 106 if ( ntoks >= 7 )
1460 {
1461
2/4
✓ Branch 1 taken 106 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 106 times.
 106 if ( !getDouble(tok[6], &x[3]) || x[3] < 0.0 )
1462 ✗ return error_setInpError(ERR_NUMBER, tok[6]);
1463 }
1464
1465 // --- add parameters to pump object
1466 106 Link[j].ID = id;
1467 106 link_setParams(j, PUMP, n1, n2, k, x);
1468 106 return 0;
1469 }
1470
1471 //=============================================================================
1472
1473 106 void pump_validate(int j, int k)
1474 //
1475 // Input: j = link index
1476 // k = pump index
1477 // Output: none
1478 // Purpose: validates a pump's properties
1479 //
1480 {
1481 int m, n1;
1482 double x, y;
1483
1484 106 Link[j].xsect.yFull = 0.0;
1485
1486 // --- check for valid curve type
1487 106 m = Pump[k].pumpCurve;
1488
2/2
✓ Branch 0 taken 23 times.
✓ Branch 1 taken 83 times.
 106 if ( m < 0 )
1489 {
1490 23 Pump[k].type = IDEAL_PUMP;
1491 }
1492 else
1493 {
1494
1/2
✓ Branch 0 taken 83 times.
✗ Branch 1 not taken.
 83 if ( Curve[m].curveType < PUMP1_CURVE ||
1495
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 83 times.
 83 Curve[m].curveType > PUMP5_CURVE )
1496 ✗ report_writeErrorMsg(ERR_NO_CURVE, Link[j].ID);
1497
1498 // --- store pump curve type with pump's parameters
1499 else
1500 {
1501 83 Pump[k].type = Curve[m].curveType - PUMP1_CURVE;
1502
1/2
✓ Branch 1 taken 83 times.
✗ Branch 2 not taken.
 83 if ( table_getFirstEntry(&Curve[m], &x, &y) )
1503 {
1504 83 Link[j].qFull = y;
1505 83 Pump[k].xMin = x;
1506 83 Pump[k].xMax = x;
1507
2/2
✓ Branch 1 taken 101 times.
✓ Branch 2 taken 83 times.
 184 while ( table_getNextEntry(&Curve[m], &x, &y) )
1508 {
1509
2/2
✓ Branch 0 taken 93 times.
✓ Branch 1 taken 8 times.
 101 Link[j].qFull = MAX(y, Link[j].qFull);
1510 101 Pump[k].xMax = x;
1511 }
1512 }
1513 83 Link[j].qFull /= UCF(FLOW);
1514 }
1515 }
1516
1517 // --- check that shutoff depth < startup depth
1518
3/4
✓ Branch 0 taken 92 times.
✓ Branch 1 taken 14 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 92 times.
 106 if ( Pump[k].yOn > 0.0 && Pump[k].yOn <= Pump[k].yOff )
1519 ✗ report_writeErrorMsg(ERR_PUMP_LIMITS, Link[j].ID);
1520
1521 // --- assign wet well volume to inlet node of Type 1 pump
1522
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 105 times.
 106 if ( Pump[k].type == TYPE1_PUMP )
1523 {
1524 1 n1 = Link[j].node1;
1525
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( Node[n1].type != STORAGE )
1526 ✗ Node[n1].fullVolume = MAX(Node[n1].fullVolume,
1527 Pump[k].xMax / UCF(VOLUME));
1528 }
1529
1530 106 }
1531
1532 //=============================================================================
1533
1534 106 void pump_initState(int j, int k)
1535 //
1536 // Input: j = link index
1537 // k = pump index
1538 // Output: none
1539 // Purpose: initializes pump conditions at start of a simulation
1540 //
1541 {
1542 106 Link[j].setting = Pump[k].initSetting;
1543 106 Link[j].targetSetting = Pump[k].initSetting;
1544 106 }
1545
1546 //=============================================================================
1547
1548 141238 double pump_getInflow(int j)
1549 //
1550 // Input: j = link index
1551 // Output: returns pump flow (cfs)
1552 // Purpose: finds flow produced by a pump.
1553 //
1554 {
1555 int k, m;
1556 int n1, n2;
1557 double vol, depth, head;
1558 141238 double qIn, qIn1, dh = 0.001;
1559 141238 double s = 1.0; // speed setting
1560
1561 141238 k = Link[j].subIndex;
1562 141238 m = Pump[k].pumpCurve;
1563 141238 n1 = Link[j].node1;
1564 141238 n2 = Link[j].node2;
1565
1566 // --- no flow if setting is closed
1567 141238 Link[j].flowClass = NO;
1568 141238 Link[j].setting = Link[j].targetSetting;
1569
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 141238 times.
 141238 if ( Link[j].setting == 0.0 ) return 0.0;
1570
1571 // --- pump flow = node inflow for IDEAL_PUMP
1572
2/2
✓ Branch 0 taken 104547 times.
✓ Branch 1 taken 36691 times.
 141238 if ( Pump[k].type == IDEAL_PUMP )
1573 104547 qIn = Node[n1].inflow + Node[n1].overflow;
1574
1575 // --- pumping rate depends on pump curve type
1576
4/5
✓ Branch 0 taken 2160 times.
✓ Branch 1 taken 2160 times.
✓ Branch 2 taken 27122 times.
✓ Branch 3 taken 5249 times.
✗ Branch 4 not taken.
 36691 else switch(Curve[m].curveType)
1577 {
1578 2160 case PUMP1_CURVE:
1579 2160 vol = Node[n1].newVolume * UCF(VOLUME);
1580 2160 qIn = table_intervalLookup(&Curve[m], vol) / UCF(FLOW);
1581
1582 // --- check if off of pump curve
1583
2/4
✓ Branch 0 taken 2160 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2160 times.
 2160 if ( vol < Pump[k].xMin || vol > Pump[k].xMax )
1584 ✗ Link[j].flowClass = YES;
1585 2160 break;
1586
1587 2160 case PUMP2_CURVE:
1588 2160 depth = Node[n1].newDepth * UCF(LENGTH);
1589 2160 qIn = table_intervalLookup(&Curve[m], depth) / UCF(FLOW);
1590
1591 // --- check if off of pump curve
1592
2/4
✓ Branch 0 taken 2160 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2160 times.
 2160 if ( depth < Pump[k].xMin || depth > Pump[k].xMax )
1593 ✗ Link[j].flowClass = YES;
1594 2160 break;
1595
1596 27122 case PUMP3_CURVE:
1597 case PUMP5_CURVE:
1598
2/2
✓ Branch 0 taken 2200 times.
✓ Branch 1 taken 24922 times.
 27122 if (Curve[m].curveType == PUMP5_CURVE) s = Link[j].setting;
1599 27122 head = ((Node[n2].newDepth + Node[n2].invertElev) -
1600 27122 (Node[n1].newDepth + Node[n1].invertElev)) / s / s;
1601
2/2
✓ Branch 0 taken 2458 times.
✓ Branch 1 taken 24664 times.
 27122 head = MAX(head, 0.0) * UCF(LENGTH);
1602 27122 qIn = table_lookup(&Curve[m], head) / UCF(FLOW);
1603
1604 // --- compute dQ/dh (slope of pump curve) and
1605 // reverse sign since flow decreases with increasing head
1606 27122 Link[j].dqdh = -table_getSlope(&Curve[m], head) *
1607 27122 UCF(LENGTH) / UCF(FLOW) / s;
1608
1609 // --- check if off of pump curve
1610
3/4
✓ Branch 0 taken 4366 times.
✓ Branch 1 taken 22756 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4366 times.
 27122 if (head < Pump[k].xMin || head > Pump[k].xMax)
1611 22756 Link[j].flowClass = YES;
1612 27122 break;
1613
1614 5249 case PUMP4_CURVE:
1615 5249 depth = Node[n1].newDepth;
1616 5249 qIn = table_lookup(&Curve[m], depth*UCF(LENGTH)) / UCF(FLOW);
1617
1618 // --- compute dQ/dh (slope of pump curve)
1619 5249 qIn1 = table_lookup(&Curve[m], (depth+dh)*UCF(LENGTH)) / UCF(FLOW);
1620 5249 Link[j].dqdh = (qIn1 - qIn) / dh;
1621
1622 // --- check if off of pump curve
1623 5249 depth *= UCF(LENGTH);
1624
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5249 times.
 5249 if ( depth < Pump[k].xMin ) Link[j].flowClass = DN_DRY;
1625
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5249 times.
 5249 if ( depth > Pump[k].xMax ) Link[j].flowClass = UP_DRY;
1626 5249 break;
1627
1628 ✗ default: qIn = 0.0;
1629 }
1630
1631 // --- do not allow reverse flow through pump
1632
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 141238 times.
 141238 if ( qIn < 0.0 ) qIn = 0.0;
1633 141238 return qIn * Link[j].setting;
1634 }
1635
1636
1637 //=============================================================================
1638 // O R I F I C E M E T H O D S
1639 //=============================================================================
1640
1641 159 int orifice_readParams(int j, int k, char* tok[], int ntoks)
1642 //
1643 // Input: j = link index
1644 // k = orifice index
1645 // tok[] = array of string tokens
1646 // ntoks = number of tokens
1647 // Output: returns an error code
1648 // Purpose: reads orifice parameters from a tokenized line of input.
1649 //
1650 {
1651 int m;
1652 int n1, n2;
1653 double x[5];
1654 char* id;
1655
1656 // --- check for valid ID and end node IDs
1657
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 159 times.
 159 if ( ntoks < 6 ) return error_setInpError(ERR_ITEMS, "");
1658 159 id = project_findID(LINK, tok[0]);
1659
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 159 times.
 159 if ( id == NULL ) return error_setInpError(ERR_NAME, tok[0]);
1660 159 n1 = project_findObject(NODE, tok[1]);
1661
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 159 times.
 159 if ( n1 < 0 ) return error_setInpError(ERR_NAME, tok[1]);
1662 159 n2 = project_findObject(NODE, tok[2]);
1663
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 159 times.
 159 if ( n2 < 0 ) return error_setInpError(ERR_NAME, tok[2]);
1664
1665 // --- parse orifice parameters
1666 159 m = findmatch(tok[3], OrificeTypeWords);
1667
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 159 times.
 159 if ( m < 0 ) return error_setInpError(ERR_KEYWORD, tok[3]);
1668 159 x[0] = m; // type
1669
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 159 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 159 if ( LinkOffsets == ELEV_OFFSET && *tok[4] == '*' ) x[1] = MISSING;
1670
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 159 times.
 159 else if ( ! getDouble(tok[4], &x[1]) ) // crest height
1671 ✗ return error_setInpError(ERR_NUMBER, tok[4]);
1672
2/4
✓ Branch 1 taken 159 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 159 times.
 159 if ( ! getDouble(tok[5], &x[2]) || x[2] < 0.0 ) // cDisch
1673 ✗ return error_setInpError(ERR_NUMBER, tok[5]);
1674 159 x[3] = 0.0;
1675
1/2
✓ Branch 0 taken 159 times.
✗ Branch 1 not taken.
 159 if ( ntoks >= 7 )
1676 {
1677 159 m = findmatch(tok[6], NoYesWords);
1678
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 159 times.
 159 if ( m < 0 ) return error_setInpError(ERR_KEYWORD, tok[6]);
1679 159 x[3] = m; // flap gate
1680 }
1681 159 x[4] = 0.0;
1682
2/2
✓ Branch 0 taken 158 times.
✓ Branch 1 taken 1 time.
 159 if ( ntoks >= 8 )
1683 {
1684
2/4
✓ Branch 1 taken 158 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 158 times.
 158 if ( ! getDouble(tok[7], &x[4]) || x[4] < 0.0 ) // orate
1685 ✗ return error_setInpError(ERR_NUMBER, tok[7]);
1686 }
1687
1688 // --- add parameters to orifice object
1689 159 Link[j].ID = id;
1690 159 link_setParams(j, ORIFICE, n1, n2, k, x);
1691 159 return 0;
1692 }
1693
1694 //=============================================================================
1695
1696 159 void orifice_validate(int j, int k)
1697 //
1698 // Input: j = link index
1699 // k = orifice index
1700 // Output: none
1701 // Purpose: validates an orifice's properties
1702 //
1703 {
1704 159 int err = 0;
1705
1706 // --- check for valid xsection
1707
2/2
✓ Branch 0 taken 117 times.
✓ Branch 1 taken 42 times.
 159 if ( Link[j].xsect.type != RECT_CLOSED
1708
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 117 times.
 117 && Link[j].xsect.type != CIRCULAR ) err = ERR_REGULATOR_SHAPE;
1709
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 159 times.
 159 if ( err > 0 )
1710 {
1711 ✗ report_writeErrorMsg(err, Link[j].ID);
1712 ✗ return;
1713 }
1714
1715 // --- check for negative offset
1716
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 159 times.
 159 if ( Link[j].offset1 < 0.0 ) Link[j].offset1 = 0.0;
1717
1718 // --- compute partial flow adjustment
1719 159 orifice_setSetting(j, 0.0);
1720
1721 // --- compute an equivalent length
1722 159 Orifice[k].length = 2.0 * RouteStep * sqrt(GRAVITY * Link[j].xsect.yFull);
1723
2/2
✓ Branch 0 taken 102 times.
✓ Branch 1 taken 57 times.
 159 Orifice[k].length = MAX(200.0, Orifice[k].length);
1724 159 Orifice[k].surfArea = 0.0;
1725 }
1726
1727 //=============================================================================
1728
1729 198 void orifice_setSetting(int j, double tstep)
1730 //
1731 // Input: j = link index
1732 // tstep = time step over which setting is adjusted (sec)
1733 // Output: none
1734 // Purpose: updates an orifice's setting as a result of a control action.
1735 //
1736 {
1737 198 int k = Link[j].subIndex;
1738 double delta, step;
1739 double h, f;
1740
1741 // --- case where adjustment rate is instantaneous
1742
3/4
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 182 times.
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 198 if ( Orifice[k].orate == 0.0 || tstep == 0.0)
1743 198 Link[j].setting = Link[j].targetSetting;
1744
1745 // --- case where orifice setting depends on time step
1746 else
1747 {
1748 ✗ delta = Link[j].targetSetting - Link[j].setting;
1749 ✗ step = tstep / Orifice[k].orate;
1750 ✗ if ( step + 0.001 >= fabs(delta) )
1751 ✗ Link[j].setting = Link[j].targetSetting;
1752 ✗ else Link[j].setting += SGN(delta) * step;
1753 }
1754
1755 // --- find effective orifice discharge coeff.
1756 198 h = Link[j].setting * Link[j].xsect.yFull;
1757 198 f = xsect_getAofY(&Link[j].xsect, h) * sqrt(2.0 * GRAVITY);
1758 198 Orifice[k].cOrif = Orifice[k].cDisch * f;
1759
1760 // --- find equiv. discharge coeff. for when weir flow occurs
1761 198 Orifice[k].cWeir = orifice_getWeirCoeff(j, k, h) * f;
1762 198 }
1763
1764 //=============================================================================
1765
1766 198 double orifice_getWeirCoeff(int j, int k, double h)
1767 //
1768 // Input: j = link index
1769 // k = orifice index
1770 // h = height of orifice opening (ft)
1771 // Output: returns a discharge coefficient (ft^1/2)
1772 // Purpose: computes the discharge coefficient for an orifice
1773 // at the critical depth where weir flow begins.
1774 //
1775 {
1776 double w, aOverL;
1777
1778 // --- this is for bottom orifices
1779
2/2
✓ Branch 0 taken 16 times.
✓ Branch 1 taken 182 times.
 198 if ( Orifice[k].type == BOTTOM_ORIFICE )
1780 {
1781 // --- find critical height above opening where orifice flow
1782 // turns into weir flow. It equals (Co/Cw)*(Area/Length)
1783 // where Co is the orifice coeff., Cw is the weir coeff/sqrt(2g),
1784 // Area is the area of the opening, and Length = circumference
1785 // of the opening. For a basic sharp crested weir, Cw = 0.414.
1786
2/2
✓ Branch 0 taken 8 times.
✓ Branch 1 taken 8 times.
 16 if (Link[j].xsect.type == CIRCULAR) aOverL = h / 4.0;
1787 else
1788 {
1789 8 w = Link[j].xsect.wMax;
1790 8 aOverL = (h*w) / (2.0*(h+w));
1791 }
1792 16 h = Orifice[k].cDisch / 0.414 * aOverL;
1793 16 Orifice[k].hCrit = h;
1794 }
1795
1796 // --- this is for side orifices
1797 else
1798 {
1799 // --- critical height is simply height of opening
1800 182 Orifice[k].hCrit = h;
1801
1802 // --- head on orifice is distance to center line
1803 182 h = h / 2.0;
1804 }
1805
1806 // --- return a coefficient for the critical depth
1807 198 return Orifice[k].cDisch * sqrt(h);
1808 }
1809
1810 //=============================================================================
1811
1812 154967 double orifice_getInflow(int j)
1813 //
1814 // Input: j = link index
1815 // Output: returns orifice flow rate (cfs)
1816 // Purpose: finds the flow through an orifice.
1817 //
1818 {
1819 int k, n1, n2;
1820 double head, h1, h2, y1, dir;
1821 double f;
1822 154967 double hcrest = 0.0;
1823 154967 double hcrown = 0.0;
1824 double hmidpt;
1825 double q, ratio;
1826
1827 // --- get indexes of end nodes and link's orifice
1828 154967 n1 = Link[j].node1;
1829 154967 n2 = Link[j].node2;
1830 154967 k = Link[j].subIndex;
1831
1832 // --- find heads at upstream & downstream nodes
1833
1/2
✓ Branch 0 taken 154967 times.
✗ Branch 1 not taken.
 154967 if ( RouteModel == DW )
1834 {
1835 154967 h1 = Node[n1].newDepth + Node[n1].invertElev;
1836 154967 h2 = Node[n2].newDepth + Node[n2].invertElev;
1837 }
1838 else
1839 {
1840 ✗ h1 = Node[n1].newDepth + Node[n1].invertElev;
1841 ✗ h2 = Node[n1].invertElev;
1842 }
1843
2/2
✓ Branch 0 taken 89783 times.
✓ Branch 1 taken 65184 times.
 154967 dir = (h1 >= h2) ? +1.0 : -1.0;
1844
1845 // --- exchange h1 and h2 for reverse flow
1846 154967 y1 = Node[n1].newDepth;
1847
2/2
✓ Branch 0 taken 65184 times.
✓ Branch 1 taken 89783 times.
 154967 if ( dir < 0.0 )
1848 {
1849 65184 head = h1;
1850 65184 h1 = h2;
1851 65184 h2 = head;
1852 65184 y1 = Node[n2].newDepth;
1853 }
1854
1855 // --- orifice is a bottom orifice (oriented in horizontal plane)
1856
2/2
✓ Branch 0 taken 1688 times.
✓ Branch 1 taken 153279 times.
 154967 if ( Orifice[k].type == BOTTOM_ORIFICE )
1857 {
1858 // --- compute crest elevation
1859 1688 hcrest = Node[n1].invertElev + Link[j].offset1;
1860
1861 // --- compute head on orifice
1862
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1688 times.
 1688 if (h1 < hcrest) head = 0.0;
1863
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1688 times.
 1688 else if (h2 > hcrest) head = h1 - h2;
1864 1688 else head = h1 - hcrest;
1865
1866 // --- find fraction of critical height for which weir flow occurs
1867 1688 f = head / Orifice[k].hCrit;
1868
1/2
✓ Branch 0 taken 1688 times.
✗ Branch 1 not taken.
 1688 f = MIN(f, 1.0);
1869 }
1870
1871 // --- otherwise orifice is a side orifice (oriented in vertical plane)
1872 else
1873 {
1874 // --- compute elevations of orifice crest and crown
1875 153279 hcrest = Node[n1].invertElev + Link[j].offset1;
1876 153279 hcrown = hcrest + Link[j].xsect.yFull * Link[j].setting;
1877 153279 hmidpt = (hcrest + hcrown) / 2.0;
1878
1879 // --- compute degree of inlet submergence
1880
3/4
✓ Branch 0 taken 144639 times.
✓ Branch 1 taken 8640 times.
✓ Branch 2 taken 144639 times.
✗ Branch 3 not taken.
 153279 if ( h1 < hcrown && hcrown > hcrest )
1881 144639 f = (h1 - hcrest) / (hcrown - hcrest);
1882 8640 else f = 1.0;
1883
1884 // --- compute head on orifice
1885
2/2
✓ Branch 0 taken 144639 times.
✓ Branch 1 taken 8640 times.
 153279 if ( f < 1.0 ) head = h1 - hcrest;
1886
1/2
✓ Branch 0 taken 8640 times.
✗ Branch 1 not taken.
 8640 else if ( h2 < hmidpt ) head = h1 - hmidpt;
1887 ✗ else head = h1 - h2;
1888 }
1889
1890 // --- return if head is negligible or flap gate closed
1891
4/6
✓ Branch 0 taken 152085 times.
✓ Branch 1 taken 2882 times.
✓ Branch 2 taken 152085 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 152085 times.
 307052 if ( head <= FUDGE || y1 <= FUDGE ||
1892 152085 link_setFlapGate(j, n1, n2, dir) )
1893 {
1894 2882 Link[j].newDepth = 0.0;
1895 2882 Link[j].flowClass = DRY;
1896 2882 Orifice[k].surfArea = FUDGE * Orifice[k].length;
1897 2882 Link[j].dqdh = 0.0;
1898 2882 return 0.0;
1899 }
1900
1901 // --- determine flow class
1902 152085 Link[j].flowClass = SUBCRITICAL;
1903
2/2
✓ Branch 0 taken 13264 times.
✓ Branch 1 taken 138821 times.
 152085 if ( hcrest > h2 )
1904 {
1905
2/2
✓ Branch 0 taken 12477 times.
✓ Branch 1 taken 787 times.
 13264 if ( dir == 1.0 ) Link[j].flowClass = DN_CRITICAL;
1906 787 else Link[j].flowClass = UP_CRITICAL;
1907 }
1908
1909 // --- compute flow depth and surface area
1910 152085 y1 = Link[j].xsect.yFull * Link[j].setting;
1911
2/2
✓ Branch 0 taken 151298 times.
✓ Branch 1 taken 787 times.
 152085 if ( Orifice[k].type == SIDE_ORIFICE )
1912 {
1913 151298 Link[j].newDepth = y1 * f;
1914 151298 Orifice[k].surfArea =
1915 151298 xsect_getWofY(&Link[j].xsect, Link[j].newDepth) *
1916 151298 Orifice[k].length;
1917 }
1918 else
1919 {
1920 787 Link[j].newDepth = y1;
1921 787 Orifice[k].surfArea = xsect_getAofY(&Link[j].xsect, y1);
1922 }
1923
1924 // --- find flow through the orifice
1925 152085 q = dir * orifice_getFlow(j, k, head, f, Link[j].hasFlapGate);
1926
1927 // --- apply Villemonte eqn. to correct for submergence
1928
4/4
✓ Branch 0 taken 143445 times.
✓ Branch 1 taken 8640 times.
✓ Branch 2 taken 129344 times.
✓ Branch 3 taken 14101 times.
 152085 if ( f < 1.0 && h2 > hcrest )
1929 {
1930 129344 ratio = (h2 - hcrest) / (h1 - hcrest);
1931 129344 q *= pow( (1.0 - pow(ratio, 1.5)), 0.385);
1932 }
1933 152085 return q;
1934 }
1935
1936 //=============================================================================
1937
1938 152085 double orifice_getFlow(int j, int k, double head, double f, int hasFlapGate)
1939 //
1940 // Input: j = link index
1941 // k = orifice index
1942 // head = head across orifice
1943 // f = fraction of critical depth filled
1944 // hasFlapGate = flap gate indicator
1945 // Output: returns flow through an orifice
1946 // Purpose: computes flow through an orifice as a function of head.
1947 //
1948 {
1949 double area, q;
1950 double veloc, hLoss;
1951
1952 // --- case where orifice is closed
1953
2/4
✓ Branch 0 taken 152085 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 152085 times.
 152085 if ( head == 0.0 || f <= 0.0 )
1954 {
1955 ✗ Link[j].dqdh = 0.0;
1956 ✗ return 0.0;
1957 }
1958
1959 // --- case where inlet depth is below critical depth;
1960 // orifice behaves as a weir
1961
2/2
✓ Branch 0 taken 143445 times.
✓ Branch 1 taken 8640 times.
 152085 else if ( f < 1.0 )
1962 {
1963 143445 q = Orifice[k].cWeir * pow(f, 1.5);
1964 143445 Link[j].dqdh = 1.5 * q / (f * Orifice[k].hCrit);
1965 }
1966
1967 // --- case where normal orifice flow applies
1968 else
1969 {
1970 8640 q = Orifice[k].cOrif * sqrt(head);
1971 8640 Link[j].dqdh = q / (2.0 * head);
1972 }
1973
1974 // --- apply ARMCO adjustment for headloss from flap gate
1975
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 152085 times.
 152085 if ( hasFlapGate )
1976 {
1977 // --- compute velocity for current orifice flow
1978 ✗ area = xsect_getAofY(&Link[j].xsect,
1979 ✗ Link[j].setting * Link[j].xsect.yFull);
1980 ✗ veloc = q / area;
1981
1982 // --- compute head loss from gate
1983 ✗ hLoss = (4.0 / GRAVITY) * veloc * veloc *
1984 ✗ exp(-1.15 * veloc / sqrt(head) );
1985
1986 // --- update head (for orifice flow)
1987 // or critical depth fraction (for weir flow)
1988 ✗ if ( f < 1.0 )
1989 {
1990 ✗ f = f - hLoss/Orifice[k].hCrit;
1991 ✗ if ( f < 0.0 ) f = 0.0;
1992 }
1993 else
1994 {
1995 ✗ head = head - hLoss;
1996 ✗ if ( head < 0.0 ) head = 0.0;
1997 }
1998
1999 // --- make recursive call to this function, with hasFlapGate
2000 // set to false, to find flow values at adjusted head value
2001 ✗ q = orifice_getFlow(j, k, head, f, FALSE);
2002 }
2003 152085 return q;
2004 }
2005
2006 //=============================================================================
2007 // W E I R M E T H O D S
2008 //=============================================================================
2009
2010 469 int weir_readParams(int j, int k, char* tok[], int ntoks)
2011 //
2012 // Input: j = link index
2013 // k = weir index
2014 // tok[] = array of string tokens
2015 // ntoks = number of tokens
2016 // Output: returns an error code
2017 // Purpose: reads weir parameters from a tokenized line of input.
2018 //
2019 {
2020 int m;
2021 int n1, n2;
2022 double x[10];
2023 char* id;
2024
2025 // --- check for valid ID and end node IDs
2026
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 469 times.
 469 if ( ntoks < 6 ) return error_setInpError(ERR_ITEMS, "");
2027 469 id = project_findID(LINK, tok[0]);
2028
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 469 times.
 469 if ( id == NULL ) return error_setInpError(ERR_NAME, tok[0]);
2029 469 n1 = project_findObject(NODE, tok[1]);
2030
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 469 times.
 469 if ( n1 < 0 ) return error_setInpError(ERR_NAME, tok[1]);
2031 469 n2 = project_findObject(NODE, tok[2]);
2032
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 469 times.
 469 if ( n2 < 0 ) return error_setInpError(ERR_NAME, tok[2]);
2033
2034 // --- parse weir parameters
2035 469 m = findmatch(tok[3], WeirTypeWords);
2036
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 469 times.
 469 if ( m < 0 ) return error_setInpError(ERR_KEYWORD, tok[3]);
2037 469 x[0] = m; // type
2038
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 469 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 469 if ( LinkOffsets == ELEV_OFFSET && *tok[4] == '*' ) x[1] = MISSING;
2039
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 469 times.
 469 else if ( ! getDouble(tok[4], &x[1]) ) // height
2040 ✗ return error_setInpError(ERR_NUMBER, tok[4]);
2041
2/4
✓ Branch 1 taken 469 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 469 times.
 469 if ( ! getDouble(tok[5], &x[2]) || x[2] < 0.0 ) // cDisch1
2042 ✗ return error_setInpError(ERR_NUMBER, tok[5]);
2043 469 x[3] = 0.0;
2044 469 x[4] = 0.0;
2045 469 x[5] = 0.0;
2046 469 x[6] = 1.0;
2047 469 x[7] = 0.0;
2048 469 x[8] = 0.0;
2049 469 x[9] = -1.0;
2050
2/4
✓ Branch 0 taken 469 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 469 times.
✗ Branch 3 not taken.
 469 if ( ntoks >= 7 && *tok[6] != '*' )
2051 {
2052 469 m = findmatch(tok[6], NoYesWords);
2053
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 469 times.
 469 if ( m < 0 ) return error_setInpError(ERR_KEYWORD, tok[6]);
2054 469 x[3] = m; // flap gate
2055 }
2056
2/4
✓ Branch 0 taken 469 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 469 times.
✗ Branch 3 not taken.
 469 if ( ntoks >= 8 && *tok[7] != '*' )
2057 {
2058
2/4
✓ Branch 1 taken 469 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 469 times.
 469 if ( ! getDouble(tok[7], &x[4]) || x[4] < 0.0 ) // endCon
2059 ✗ return error_setInpError(ERR_NUMBER, tok[7]);
2060 }
2061
2/4
✓ Branch 0 taken 469 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 469 times.
✗ Branch 3 not taken.
 469 if ( ntoks >= 9 && *tok[8] != '*' )
2062 {
2063
2/4
✓ Branch 1 taken 469 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 469 times.
 469 if ( ! getDouble(tok[8], &x[5]) || x[5] < 0.0 ) // cDisch2
2064 ✗ return error_setInpError(ERR_NUMBER, tok[8]);
2065 }
2066
2067
2/4
✓ Branch 0 taken 469 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 469 times.
✗ Branch 3 not taken.
 469 if ( ntoks >= 10 && *tok[9] != '*' )
2068 {
2069 469 m = findmatch(tok[9], NoYesWords);
2070
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 469 times.
 469 if ( m < 0 ) return error_setInpError(ERR_KEYWORD, tok[9]);
2071 469 x[6] = m; // canSurcharge
2072 }
2073
2074
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 467 times.
 469 if ( (m = (int)x[0]) == ROADWAY_WEIR )
2075 {
2076
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if ( ntoks >= 11 ) // road width
2077 {
2078
2/4
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 2 times.
 2 if ( ! getDouble(tok[10], &x[7]) || x[7] < 0.0 )
2079 ✗ return error_setInpError(ERR_NUMBER, tok[10]);
2080 }
2081
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if ( ntoks >= 12 ) // road surface
2082 {
2083
2/2
✓ Branch 1 taken 1 time.
✓ Branch 2 taken 1 time.
 2 if ( strcomp(tok[11], "PAVED") ) x[8] = 1.0;
2084
1/2
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
 1 else if ( strcomp(tok[11], "GRAVEL") ) x[8] = 2.0;
2085 }
2086 }
2087
2088
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 469 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 469 if (ntoks >= 13 && *tok[12] != '*')
2089 {
2090 ✗ m = project_findObject(CURVE, tok[12]); // coeff. curve
2091 ✗ if (m < 0) return error_setInpError(ERR_NAME, tok[12]);
2092 ✗ x[9] = m;
2093 }
2094
2095 // --- add parameters to weir object
2096 469 Link[j].ID = id;
2097 469 link_setParams(j, WEIR, n1, n2, k, x);
2098 469 return 0;
2099 }
2100
2101 //=============================================================================
2102
2103 469 void weir_validate(int j, int k)
2104 //
2105 // Input: j = link index
2106 // k = weir index
2107 // Output: none
2108 // Purpose: validates a weir's properties
2109 //
2110 {
2111 469 int err = 0;
2112 double q, q1, q2, head;
2113
2114 // --- check for valid cross section
2115
3/4
✓ Branch 0 taken 466 times.
✓ Branch 1 taken 2 times.
✓ Branch 2 taken 1 time.
✗ Branch 3 not taken.
 469 switch ( Weir[k].type)
2116 {
2117 466 case TRANSVERSE_WEIR:
2118 case SIDEFLOW_WEIR:
2119 case ROADWAY_WEIR:
2120
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 466 times.
 466 if ( Link[j].xsect.type != RECT_OPEN ) err = ERR_REGULATOR_SHAPE;
2121 466 Weir[k].slope = 0.0;
2122 466 break;
2123
2124 2 case VNOTCH_WEIR:
2125
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( Link[j].xsect.type != TRIANGULAR ) err = ERR_REGULATOR_SHAPE;
2126 else
2127 {
2128 2 Weir[k].slope = Link[j].xsect.sBot;
2129 }
2130 2 break;
2131
2132 1 case TRAPEZOIDAL_WEIR:
2133
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
 1 if ( Link[j].xsect.type != TRAPEZOIDAL ) err = ERR_REGULATOR_SHAPE;
2134 else
2135 {
2136 1 Weir[k].slope = Link[j].xsect.sBot;
2137 }
2138 1 break;
2139 }
2140
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 469 times.
 469 if ( err > 0 )
2141 {
2142 ✗ report_writeErrorMsg(err, Link[j].ID);
2143 ✗ return;
2144 }
2145
2146 // --- check for negative offset
2147
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 469 times.
 469 if ( Link[j].offset1 < 0.0 ) Link[j].offset1 = 0.0;
2148
2149 // --- compute an equivalent length
2150 469 Weir[k].length = 2.0 * RouteStep * sqrt(GRAVITY * Link[j].xsect.yFull);
2151
2/2
✓ Branch 0 taken 19 times.
✓ Branch 1 taken 450 times.
 469 Weir[k].length = MAX(200.0, Weir[k].length);
2152 469 Weir[k].surfArea = 0.0;
2153
2154 // --- find flow through weir when water level equals weir height
2155 469 head = Link[j].xsect.yFull;
2156 469 weir_getFlow(j, k, head, 1.0, FALSE, &q1, &q2);
2157 469 q = q1 + q2;
2158
2159 // --- compute equivalent orifice coeff. (for CFS flow units)
2160 469 head = head / 2.0; // head seen by equivalent orifice
2161 469 Weir[k].cSurcharge = q / sqrt(head);
2162 }
2163
2164 //=============================================================================
2165
2166 78 void weir_setSetting(int j)
2167 //
2168 // Input: j = link index
2169 // Output: none
2170 // Purpose: updates a weir's setting as a result of a control action.
2171 //
2172 {
2173 78 int k = Link[j].subIndex;
2174 double h, q, q1, q2;
2175
2176 // --- adjust weir setting
2177 78 Link[j].setting = Link[j].targetSetting;
2178
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
 78 if ( !Weir[k].canSurcharge ) return;
2179
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 78 times.
 78 if ( Weir[k].type == ROADWAY_WEIR ) return;
2180
2181 // --- find orifice coeff. for surcharged flow
2182
2/2
✓ Branch 0 taken 7 times.
✓ Branch 1 taken 71 times.
 78 if ( Link[j].setting == 0.0 ) Weir[k].cSurcharge = 0.0;
2183 else
2184 {
2185 // --- find flow through weir when water level equals weir height
2186 71 h = Link[j].setting * Link[j].xsect.yFull;
2187 71 weir_getFlow(j, k, h, 1.0, FALSE, &q1, &q2);
2188 71 q = q1 + q2;
2189
2190 // --- compute equivalent orifice coeff. (for CFS flow units)
2191 71 h = h / 2.0; // head seen by equivalent orifice
2192 71 Weir[k].cSurcharge = q / sqrt(h);
2193 }
2194 }
2195
2196 //=============================================================================
2197
2198 541350 double weir_getInflow(int j)
2199 //
2200 // Input: j = link index
2201 // Output: returns weir flow rate (cfs)
2202 // Purpose: finds the flow over a weir.
2203 //
2204 {
2205 int n1; // index of upstream node
2206 int n2; // index of downstream node
2207 int k; // index of weir
2208 double q1; // flow through central part of weir (cfs)
2209 double q2; // flow through end sections of weir (cfs)
2210 double head; // head on weir (ft)
2211 double h1; // upstrm nodal head (ft)
2212 double h2; // downstrm nodal head (ft)
2213 double hcrest; // head at weir crest (ft)
2214 double hcrown; // head at weir crown (ft)
2215 double y; // water depth in weir (ft)
2216 double dir; // direction multiplier
2217 double ratio;
2218 541350 double weirPower[] = {1.5, // transverse weir
2219 5./3., // side flow weir
2220 2.5, // v-notch weir
2221 1.5}; // trapezoidal weir
2222
2223 541350 n1 = Link[j].node1;
2224 541350 n2 = Link[j].node2;
2225 541350 k = Link[j].subIndex;
2226
1/2
✓ Branch 0 taken 541350 times.
✗ Branch 1 not taken.
 541350 if ( RouteModel == DW )
2227 {
2228 541350 h1 = Node[n1].newDepth + Node[n1].invertElev;
2229 541350 h2 = Node[n2].newDepth + Node[n2].invertElev;
2230 }
2231 else
2232 {
2233 ✗ h1 = Node[n1].newDepth + Node[n1].invertElev;
2234 ✗ h2 = Node[n1].invertElev;
2235 }
2236
2/2
✓ Branch 0 taken 493033 times.
✓ Branch 1 taken 48317 times.
 541350 dir = (h1 > h2) ? +1.0 : -1.0;
2237
2238 // --- exchange h1 and h2 for reverse flow
2239
2/2
✓ Branch 0 taken 48317 times.
✓ Branch 1 taken 493033 times.
 541350 if ( dir < 0.0 )
2240 {
2241 48317 head = h1;
2242 48317 h1 = h2;
2243 48317 h2 = head;
2244 }
2245
2246 // --- find head of weir's crest and crown
2247 541350 hcrest = Node[n1].invertElev + Link[j].offset1;
2248 541350 hcrown = hcrest + Link[j].xsect.yFull;
2249
2250 // --- treat a roadway weir as a special case
2251
2/2
✓ Branch 0 taken 72008 times.
✓ Branch 1 taken 469342 times.
 541350 if ( Weir[k].type == ROADWAY_WEIR )
2252 72008 return roadway_getInflow(j, dir, hcrest, h1, h2);
2253
2254 // --- adjust crest ht. for partially open weir
2255 469342 hcrest += (1.0 - Link[j].setting) * Link[j].xsect.yFull;
2256
2257 // --- compute head relative to weir crest
2258 469342 head = h1 - hcrest;
2259
2260 // --- return if head is negligible or flap gate closed
2261 469342 Link[j].dqdh = 0.0;
2262
5/6
✓ Branch 0 taken 263359 times.
✓ Branch 1 taken 205983 times.
✓ Branch 2 taken 263359 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 594 times.
✓ Branch 5 taken 262765 times.
 732701 if ( head <= FUDGE || hcrest >= hcrown ||
2263 263359 link_setFlapGate(j, n1, n2, dir) )
2264 {
2265 206577 Link[j].newDepth = 0.0;
2266 206577 Link[j].flowClass = DRY;
2267 206577 return 0.0;
2268 }
2269
2270 // --- determine flow class
2271 262765 Link[j].flowClass = SUBCRITICAL;
2272
2/2
✓ Branch 0 taken 74268 times.
✓ Branch 1 taken 188497 times.
 262765 if ( hcrest > h2 )
2273 {
2274
2/2
✓ Branch 0 taken 67286 times.
✓ Branch 1 taken 6982 times.
 74268 if ( dir == 1.0 ) Link[j].flowClass = DN_CRITICAL;
2275 6982 else Link[j].flowClass = UP_CRITICAL;
2276 }
2277
2278 // --- compute new equivalent surface area
2279
2/2
✓ Branch 0 taken 220361 times.
✓ Branch 1 taken 42404 times.
 262765 y = Link[j].xsect.yFull - (hcrown - MIN(h1, hcrown));
2280 262765 Weir[k].surfArea = xsect_getWofY(&Link[j].xsect, y) * Weir[k].length;
2281
2282 // --- head is above crown
2283
2/2
✓ Branch 0 taken 42404 times.
✓ Branch 1 taken 220361 times.
 262765 if ( h1 >= hcrown )
2284 {
2285 // --- use equivalent orifice if weir can surcharge
2286
1/2
✓ Branch 0 taken 42404 times.
✗ Branch 1 not taken.
 42404 if ( Weir[k].canSurcharge )
2287 {
2288 42404 y = (hcrest + hcrown) / 2.0;
2289
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 42404 times.
 42404 if ( h2 < y ) head = h1 - y;
2290 42404 else head = h1 - h2;
2291 42404 y = hcrown - hcrest;
2292 42404 q1 = weir_getOrificeFlow(j, head, y, Weir[k].cSurcharge);
2293 42404 Link[j].newDepth = y;
2294 42404 return dir * q1;
2295 }
2296
2297 // --- otherwise limit head to height of weir opening
2298 ✗ else head = hcrown - hcrest;
2299 }
2300
2301 // --- use weir eqn. to find flows through central (q1)
2302 // and end sections (q2) of weir
2303 220361 weir_getFlow(j, k, head, dir, Link[j].hasFlapGate, &q1, &q2);
2304
2305 // --- apply Villemonte eqn. to correct for submergence
2306
2/2
✓ Branch 0 taken 146093 times.
✓ Branch 1 taken 74268 times.
 220361 if ( h2 > hcrest )
2307 {
2308 146093 ratio = (h2 - hcrest) / (h1 - hcrest);
2309 146093 q1 *= pow( (1.0 - pow(ratio, weirPower[Weir[k].type])), 0.385);
2310
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 146093 times.
 146093 if ( q2 > 0.0 )
2311 ✗ q2 *= pow( (1.0 - pow(ratio, weirPower[VNOTCH_WEIR])), 0.385);
2312 }
2313
2314 // --- return total flow through weir
2315
1/2
✓ Branch 0 taken 220361 times.
✗ Branch 1 not taken.
 220361 Link[j].newDepth = MIN((h1 - hcrest), Link[j].xsect.yFull);
2316 220361 return dir * (q1 + q2);
2317 }
2318
2319 //=============================================================================
2320
2321 220901 void weir_getFlow(int j, int k, double head, double dir, int hasFlapGate,
2322 double* q1, double* q2)
2323 //
2324 // Input: j = link index
2325 // k = weir index
2326 // head = head across weir (ft)
2327 // dir = flow direction indicator
2328 // hasFlapGate = flap gate indicator
2329 // Output: q1 = flow through central portion of weir (cfs)
2330 // q2 = flow through end sections of weir (cfs)
2331 // Purpose: computes flow over weir given head.
2332 //
2333 {
2334 double length;
2335 double h;
2336 double y;
2337 double hLoss;
2338 double area;
2339 double veloc;
2340 int wType;
2341 220901 int cdCurve = Weir[k].cdCurve;
2342 220901 double cDisch1 = Weir[k].cDisch1;
2343
2344 // --- q1 = flow through central portion of weir,
2345 // q2 = flow through end sections of trapezoidal weir
2346 220901 *q1 = 0.0;
2347 220901 *q2 = 0.0;
2348 220901 Link[j].dqdh = 0.0;
2349
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 220901 times.
 220901 if ( head <= 0.0 ) return;
2350
2351 // --- convert weir length & head to original units
2352 220901 length = Link[j].xsect.wMax * UCF(LENGTH);
2353 220901 h = head * UCF(LENGTH);
2354
2355 // --- lookup tabulated discharge coeff.
2356
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 220901 times.
 220901 if ( cdCurve >= 0 ) cDisch1 = table_lookup(&Curve[cdCurve], h);
2357
2358 // --- use appropriate formula for weir flow
2359 220901 wType = Weir[k].type;
2360
2/2
✓ Branch 0 taken 536 times.
✓ Branch 1 taken 220365 times.
 220901 if ( wType == VNOTCH_WEIR &&
2361
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 536 times.
 536 Link[j].setting < 1.0 ) wType = TRAPEZOIDAL_WEIR;
2362
5/5
✓ Branch 0 taken 119000 times.
✓ Branch 1 taken 200 times.
✓ Branch 2 taken 536 times.
✓ Branch 3 taken 101163 times.
✓ Branch 4 taken 2 times.
 220901 switch (wType)
2363 {
2364 119000 case TRANSVERSE_WEIR:
2365
2366 // --- reduce length when end contractions present
2367 119000 length -= 0.1 * Weir[k].endCon * h;
2368
1/2
✓ Branch 0 taken 119000 times.
✗ Branch 1 not taken.
 119000 length = MAX(length, 0.0);
2369 119000 *q1 = cDisch1 * length * pow(h, 1.5);
2370 119000 break;
2371
2372 200 case SIDEFLOW_WEIR:
2373
2374 // --- reduce length when end contractions present
2375 200 length -= 0.1 * Weir[k].endCon * h;
2376
1/2
✓ Branch 0 taken 200 times.
✗ Branch 1 not taken.
 200 length = MAX(length, 0.0);
2377
2378 // --- weir behaves as a transverse weir under reverse flow
2379
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 200 times.
 200 if ( dir < 0.0 )
2380 ✗ *q1 = cDisch1 * length * pow(h, 1.5);
2381 else
2382
2383 // Corrected formula (see Metcalf & Eddy, Inc.,
2384 // Wastewater Engineering, McGraw-Hill, 1972 p. 164).
2385 200 *q1 = cDisch1 * pow(length, 0.83) * pow(h, 1.67);
2386
2387 200 break;
2388
2389 536 case VNOTCH_WEIR:
2390 536 *q1 = cDisch1 * Weir[k].slope * pow(h, 2.5);
2391 536 break;
2392
2393 101163 case TRAPEZOIDAL_WEIR:
2394 101163 y = (1.0 - Link[j].setting) * Link[j].xsect.yFull;
2395 101163 length = xsect_getWofY(&Link[j].xsect, y) * UCF(LENGTH);
2396 101163 *q1 = cDisch1 * length * pow(h, 1.5);
2397 101163 *q2 = Weir[k].cDisch2 * Weir[k].slope * pow(h, 2.5);
2398 }
2399
2400 // --- convert CMS flows to CFS
2401
2/2
✓ Branch 0 taken 7519 times.
✓ Branch 1 taken 213382 times.
 220901 if ( UnitSystem == SI )
2402 {
2403 7519 *q1 /= M3perFT3;
2404 7519 *q2 /= M3perFT3;
2405 }
2406
2407 // --- apply ARMCO adjustment for headloss from flap gate
2408
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 220901 times.
 220901 if ( hasFlapGate )
2409 {
2410 // --- compute flow area & velocity for current weir flow
2411 ✗ area = weir_getOpenArea(j, head);
2412 ✗ if ( area > TINY )
2413 {
2414 ✗ veloc = (*q1 + *q2) / area;
2415
2416 // --- compute headloss and subtract from original head
2417 ✗ hLoss = (4.0 / GRAVITY) * veloc * veloc *
2418 ✗ exp(-1.15 * veloc / sqrt(head) );
2419 ✗ head = head - hLoss;
2420 ✗ if ( head < 0.0 ) head = 0.0;
2421
2422 // --- make recursive call to this function, with hasFlapGate
2423 // set to false, to find flow values at adjusted head value
2424 ✗ weir_getFlow(j, k, head, dir, FALSE, q1, q2);
2425 }
2426 }
2427 220901 Link[j].dqdh = weir_getdqdh(k, dir, head, *q1, *q2);
2428 }
2429
2430 //=============================================================================
2431
2432 42404 double weir_getOrificeFlow(int j, double head, double y, double cOrif)
2433 //
2434 // Input: j = link index
2435 // head = head across weir (ft)
2436 // y = height of upstream water level above weir crest (ft)
2437 // cOrif = orifice flow coefficient
2438 // Output: returns flow through weir
2439 // Purpose: finds flow through a surcharged weir using the orifice equation.
2440 //
2441 {
2442 double a, q, v, hloss;
2443
2444 // --- evaluate the orifice flow equation
2445 42404 q = cOrif * sqrt(head);
2446
2447 // --- apply Armco adjustment if weir has a flap gate
2448
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 42404 times.
 42404 if ( Link[j].hasFlapGate )
2449 {
2450 ✗ a = weir_getOpenArea(j, y);
2451 ✗ if ( a > 0.0 )
2452 {
2453 ✗ v = q / a;
2454 ✗ hloss = (4.0 / GRAVITY) * v * v * exp(-1.15 * v / sqrt(y) );
2455 ✗ head -= hloss;
2456 ✗ head = MAX(head, 0.0);
2457 ✗ q = cOrif * sqrt(head);
2458 }
2459 }
2460
1/2
✓ Branch 0 taken 42404 times.
✗ Branch 1 not taken.
 42404 if ( head > 0.0 ) Link[j].dqdh = q / (2.0 * head);
2461 ✗ else Link[j].dqdh = 0.0;
2462 42404 return q;
2463 }
2464
2465 //=============================================================================
2466
2467 ✗ double weir_getOpenArea(int j, double y)
2468 //
2469 // Input: j = link index
2470 // y = depth of water above weir crest (ft)
2471 // Output: returns area between weir crest and y (ft2)
2472 // Purpose: finds flow area through a weir.
2473 //
2474 {
2475 double z, zy;
2476
2477 // --- find offset of weir crest due to control setting
2478 ✗ z = (1.0 - Link[j].setting) * Link[j].xsect.yFull;
2479
2480 // --- ht. of crest + ht of water above crest
2481 ✗ zy = z + y;
2482 ✗ zy = MIN(zy, Link[j].xsect.yFull);
2483
2484 // --- return difference between area of offset + water depth
2485 // and area of just the offset
2486 ✗ return xsect_getAofY(&Link[j].xsect, zy) -
2487 ✗ xsect_getAofY(&Link[j].xsect, z);
2488 }
2489
2490 //=============================================================================
2491
2492 220901 double weir_getdqdh(int k, double dir, double h, double q1, double q2)
2493 {
2494 double q1h;
2495 double q2h;
2496
2497
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 220901 times.
 220901 if ( fabs(h) < FUDGE ) return 0.0;
2498 220901 q1h = fabs(q1/h);
2499 220901 q2h = fabs(q2/h);
2500
2501
5/5
✓ Branch 0 taken 119000 times.
✓ Branch 1 taken 200 times.
✓ Branch 2 taken 536 times.
✓ Branch 3 taken 101163 times.
✓ Branch 4 taken 2 times.
 220901 switch (Weir[k].type)
2502 {
2503 119000 case TRANSVERSE_WEIR: return 1.5 * q1h;
2504
2505 200 case SIDEFLOW_WEIR:
2506 // --- weir behaves as a transverse weir under reverse flow
2507
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 200 times.
 200 if ( dir < 0.0 ) return 1.5 * q1h;
2508 200 else return 1.67 * q1h;
2509
2510 536 case VNOTCH_WEIR:
2511
1/2
✓ Branch 0 taken 536 times.
✗ Branch 1 not taken.
 536 if ( q2h == 0.0 ) return 2.5 * q1h; // Fully open
2512 ✗ else return 1.5 * q1h + 2.5 * q2h; // Partly open
2513
2514 101163 case TRAPEZOIDAL_WEIR: return 1.5 * q1h + 2.5 * q2h;
2515 }
2516 2 return 0.0;
2517 }
2518
2519
2520 //=============================================================================
2521 // O U T L E T D E V I C E M E T H O D S
2522 //=============================================================================
2523
2524 2 int outlet_readParams(int j, int k, char* tok[], int ntoks)
2525 //
2526 // Input: j = link index
2527 // k = outlet index
2528 // tok[] = array of string tokens
2529 // ntoks = number of tokens
2530 // Output: returns an error code
2531 // Purpose: reads outlet parameters from a tokenized line of input.
2532 //
2533 {
2534 int i, m, n;
2535 int n1, n2;
2536 double x[6];
2537 char* id;
2538 char* s;
2539
2540 // --- check for valid ID and end node IDs
2541
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( ntoks < 6 ) return error_setInpError(ERR_ITEMS, "");
2542 2 id = project_findID(LINK, tok[0]);
2543
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( id == NULL ) return error_setInpError(ERR_NAME, tok[0]);
2544 2 n1 = project_findObject(NODE, tok[1]);
2545
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( n1 < 0 ) return error_setInpError(ERR_NAME, tok[1]);
2546 2 n2 = project_findObject(NODE, tok[2]);
2547
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( n2 < 0 ) return error_setInpError(ERR_NAME, tok[2]);
2548
2549 // --- get height above invert
2550
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 2 if ( LinkOffsets == ELEV_OFFSET && *tok[3] == '*' ) x[0] = MISSING;
2551 else
2552 {
2553
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 2 if ( ! getDouble(tok[3], &x[0]) )
2554 ✗ return error_setInpError(ERR_NUMBER, tok[3]);
2555
2/4
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 2 times.
 2 if ( LinkOffsets == DEPTH_OFFSET && x[0] < 0.0 ) x[0] = 0.0;
2556 }
2557
2558 // --- see if outlet flow relation is tabular or functional
2559 2 m = findmatch(tok[4], RelationWords);
2560
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( m < 0 ) return error_setInpError(ERR_KEYWORD, tok[4]);
2561 2 x[1] = 0.0;
2562 2 x[2] = 0.0;
2563 2 x[3] = -1.0;
2564 2 x[4] = 0.0;
2565
2566 // --- see if rating curve is head or depth based
2567 2 x[5] = NODE_DEPTH; //default is depth-based
2568 2 s = strtok(tok[4], "/"); //parse token for
2569 2 s = strtok(NULL, "/"); // qualifier term
2570
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 2 if ( strcomp(s, w_HEAD) ) x[5] = NODE_HEAD; //check if its "HEAD"
2571
2572 // --- get params. for functional outlet device
2573
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if ( m == FUNCTIONAL )
2574 {
2575
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( ntoks < 7 ) return error_setInpError(ERR_ITEMS, "");
2576
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 2 if ( ! getDouble(tok[5], &x[1]) )
2577 ✗ return error_setInpError(ERR_NUMBER, tok[5]);
2578
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 2 times.
 2 if ( ! getDouble(tok[6], &x[2]) )
2579 ✗ return error_setInpError(ERR_NUMBER, tok[6]);
2580 2 n = 7;
2581 }
2582
2583 // --- get name of outlet rating curve
2584 else
2585 {
2586 ✗ i = project_findObject(CURVE, tok[5]);
2587 ✗ if ( i < 0 ) return error_setInpError(ERR_NAME, tok[5]);
2588 ✗ x[3] = i;
2589 ✗ n = 6;
2590 }
2591
2592 // --- check if flap gate specified
2593
1/2
✓ Branch 0 taken 2 times.
✗ Branch 1 not taken.
 2 if ( ntoks > n)
2594 {
2595 2 i = findmatch(tok[n], NoYesWords);
2596
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 2 times.
 2 if ( i < 0 ) return error_setInpError(ERR_KEYWORD, tok[n]);
2597 2 x[4] = i;
2598 }
2599
2600 // --- add parameters to outlet object
2601 2 Link[j].ID = id;
2602 2 link_setParams(j, OUTLET, n1, n2, k, x);
2603 2 return 0;
2604 }
2605
2606 //=============================================================================
2607
2608 57604 double outlet_getInflow(int j)
2609 //
2610 // Input: j = link index
2611 // Output: outlet flow rate (cfs)
2612 // Purpose: finds the flow through an outlet.
2613 //
2614 {
2615 int k, n1, n2;
2616 double head, hcrest, h1, h2, y1, dir;
2617
2618 // --- get indexes of end nodes
2619 57604 n1 = Link[j].node1;
2620 57604 n2 = Link[j].node2;
2621 57604 k = Link[j].subIndex;
2622
2623 // --- find heads at upstream & downstream nodes
2624
1/2
✓ Branch 0 taken 57604 times.
✗ Branch 1 not taken.
 57604 if ( RouteModel == DW )
2625 {
2626 57604 h1 = Node[n1].newDepth + Node[n1].invertElev;
2627 57604 h2 = Node[n2].newDepth + Node[n2].invertElev;
2628 }
2629 else
2630 {
2631 ✗ h1 = Node[n1].newDepth + Node[n1].invertElev;
2632 ✗ h2 = Node[n1].invertElev;
2633 }
2634
1/2
✓ Branch 0 taken 57604 times.
✗ Branch 1 not taken.
 57604 dir = (h1 >= h2) ? +1.0 : -1.0;
2635
2636 // --- exchange h1 and h2 for reverse flow
2637 57604 y1 = Node[n1].newDepth;
2638
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 57604 times.
 57604 if ( dir < 0.0 )
2639 {
2640 ✗ y1 = h1;
2641 ✗ h1 = h2;
2642 ✗ h2 = y1;
2643 ✗ y1 = Node[n2].newDepth;
2644 }
2645
2646 // --- for a NODE_DEPTH rating curve the effective head across the
2647 // outlet is the depth above the crest elev. while for a NODE_HEAD
2648 // curve it is the difference between upstream & downstream heads
2649 57604 hcrest = Node[n1].invertElev + Link[j].offset1;
2650
1/4
✗ Branch 0 not taken.
✓ Branch 1 taken 57604 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
 57604 if ( Outlet[k].curveType == NODE_HEAD && RouteModel == DW )
2651 ✗ head = h1 - MAX(h2, hcrest);
2652 57604 else head = h1 - hcrest;
2653
2654 // --- no flow if either no effective head difference,
2655 // no upstream water available, or closed flap gate
2656
4/6
✓ Branch 0 taken 14402 times.
✓ Branch 1 taken 43202 times.
✓ Branch 2 taken 14402 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 14402 times.
 72006 if ( head <= FUDGE || y1 <= FUDGE ||
2657 14402 link_setFlapGate(j, n1, n2, dir) )
2658 {
2659 43202 Link[j].newDepth = 0.0;
2660 43202 Link[j].flowClass = DRY;
2661 43202 return 0.0;
2662 }
2663
2664 // --- otherwise use rating curve to compute flow
2665 14402 Link[j].newDepth = head;
2666 14402 Link[j].flowClass = SUBCRITICAL;
2667 14402 return dir * Link[j].setting * outlet_getFlow(k, head);
2668 }
2669
2670 //=============================================================================
2671
2672 14402 double outlet_getFlow(int k, double head)
2673 //
2674 // Input: k = outlet index
2675 // head = head across outlet (ft)
2676 // Output: returns outlet flow rate (cfs)
2677 // Purpose: computes flow rate through an outlet given head.
2678 //
2679 {
2680 int m;
2681 double h;
2682
2683 // --- convert head to original units
2684 14402 h = head * UCF(LENGTH);
2685
2686 // --- look-up flow in rating curve table if provided
2687 14402 m = Outlet[k].qCurve;
2688
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 14402 times.
 14402 if ( m >= 0 ) return table_lookup(&Curve[m], h) / UCF(FLOW);
2689
2690 // --- otherwise use function to find flow
2691 14402 else return Outlet[k].qCoeff * pow(h, Outlet[k].qExpon) / UCF(FLOW);
2692 }
2693