GCC Code Coverage Report


Directory: src/solver/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 90.6% 310 / 0 / 342
Functions: 100.0% 23 / 0 / 23
Branches: 83.2% 218 / 0 / 262

routing.c
Line Branch Exec Source
1 //-----------------------------------------------------------------------------
2 // routing.c
3 //
4 // Project: EPA SWMM5
5 // Version: 5.2
6 // Date: 11/01/21 (Build 5.2.0)
7 // Author: L. Rossman
8 // M. Tryby (EPA)
9 //
10 // Conveyance system routing functions.
11 //
12 // Update History
13 // ==============
14 // Build 5.1.007:
15 // - Nodal evap/seepage losses computed using conditions at start of time step.
16 // - DWF pollutant concentrations ignored if DWF is negative.
17 // - Separate mass balance accounting made for storage evap. & seepage.
18 // - Nodal mass balance accounting for negative lateral inflows corrected.
19 // Build 5.1.008:
20 // - Initialization of flow and quality routing systems moved here from swmm5.c.
21 // - Lateral inflows now evaluated at start (not end) of time step.
22 // - Flows from LID drains included in lateral inflows.
23 // - Conduit evap/seepage losses multiplied by number of barrels before
24 // being added into mass balances.
25 // Build 5.1.010:
26 // - Time when a link's setting is changed is recorded.
27 // Build 5.1.011:
28 // - Support added for limiting flow routing to specific events.
29 // Build 5.1.012:
30 // - routing_execute() was re-written so that Routing Events and
31 // Skip Steady Flow options work together correctly.
32 // Build 5.1.013:
33 // - Support added for evaluating controls rules at RuleStep time interval.
34 // - Back flow through Outfall nodes now treated as External Inflows for
35 // mass balance purposes.
36 // - Global infiltration factor for storage seepage set in routing_execute.
37 // Build 5.2.0:
38 // - Support added for street flow capture and sewer backflow thru inlets.
39 // - Shell sort replaces insertion sort for sorting Event array.
40 //-----------------------------------------------------------------------------
41 #define _CRT_SECURE_NO_DEPRECATE
42
43 #include <stdio.h>
44 #include <stdlib.h>
45 #include <string.h>
46 #include <math.h>
47 #include "headers.h"
48 #include "lid.h"
49 //-----------------------------------------------------------------------------
50 // Shared variables
51 //-----------------------------------------------------------------------------
52 static int* SortedLinks;
53 static int NextEvent;
54 static int BetweenEvents;
55 static double NewRuleTime;
56
57 //-----------------------------------------------------------------------------
58 // External functions (declared in funcs.h)
59 //-----------------------------------------------------------------------------
60 // routing_open (called by swmm_start in swmm5.c)
61 // routing_getRoutingStep (called by swmm_step in swmm5.c)
62 // routing_execute (called by swmm_step in swmm5.c)
63 // routing_close (called by swmm_end in swmm5.c)
64
65 //-----------------------------------------------------------------------------
66 // Function declarations
67 //-----------------------------------------------------------------------------
68 static int evaluateControlRules(DateTime currentDate, double routingStep);
69 static void sortEvents(void);
70 static int isBetweenEvents(DateTime currentDate);
71 static int isInSteadyState(int actionCount, double stepFlowError);
72 static int inflowHasChanged(void);
73
74 static void initSystemInflows();
75 static void addSystemInflows(DateTime currentDate, double routingStep);
76 static void addExternalInflows(DateTime currentDate);
77 static void addDryWeatherInflows(DateTime currentDate);
78 static void addWetWeatherInflows(double routingTime);
79 static void addGroundwaterInflows(double routingTime);
80 static void addRdiiInflows(DateTime currentDate);
81 static void addIfaceInflows(DateTime currentDate);
82 static void addLidDrainInflows(double routingTime);
83
84 static int routeFlow(int routingModel, double routingStep);
85 static void removeSystemOutflows(double routingStep);
86 static void removeStorageLosses(double tStep);
87 static void removeConduitLosses(void);
88 static void removeOutflows(double tStep);
89
90
91 //=============================================================================
92
93 58 int routing_open()
94 //
95 // Input: none
96 // Output: returns an error code
97 // Purpose: initializes the routing analyzer.
98 //
99 {
100 // --- open treatment system
101
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 58 times.
58 if ( !treatmnt_open() ) return ErrorCode;
102
103 // --- topologically sort the links
104 58 SortedLinks = NULL;
105
2/2
✓ Branch 0 taken 48 times.
✓ Branch 1 taken 10 times.
58 if ( Nobjects[LINK] > 0 )
106 {
107 48 SortedLinks = (int *) calloc(Nobjects[LINK], sizeof(int));
108
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 48 times.
48 if ( !SortedLinks )
109 {
110 report_writeErrorMsg(ERR_MEMORY, "");
111 return ErrorCode;
112 }
113 48 toposort_sortLinks(SortedLinks);
114
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 48 times.
48 if ( ErrorCode ) return ErrorCode;
115 }
116
117 // --- open any routing interface files
118 58 iface_openRoutingFiles();
119
120 // --- initialize flow and quality routing systems
121 58 flowrout_init(RouteModel);
122
2/2
✓ Branch 0 taken 57 times.
✓ Branch 1 taken 1 time.
58 if ( Fhotstart1.mode == NO_FILE ) qualrout_init();
123
124 // --- initialize routing events
125
2/2
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 49 times.
58 if ( NumEvents > 0 ) sortEvents();
126 58 NextEvent = 0;
127 58 BetweenEvents = (NumEvents > 0);
128 58 NewRuleTime = 0.0;
129 58 return ErrorCode;
130 }
131
132 //=============================================================================
133
134 58 void routing_close(int routingModel)
135 //
136 // Input: routingModel = routing method code
137 // Output: none
138 // Purpose: closes down the routing analyzer.
139 //
140 {
141 // --- close any routing interface files
142 58 iface_closeRoutingFiles();
143
144 // --- free allocated memory
145 58 flowrout_close(routingModel);
146 58 treatmnt_close();
147
2/2
✓ Branch 0 taken 48 times.
✓ Branch 1 taken 10 times.
58 FREE(SortedLinks);
148 58 }
149
150 //=============================================================================
151
152 1066421 double routing_getRoutingStep(int routingModel, double fixedStep)
153 //
154 // Input: routingModel = routing method code
155 // fixedStep = user-supplied time step (sec)
156 // Output: returns a routing time step (sec)
157 // Purpose: determines time step used for flow routing at current time period.
158 //
159 {
160 double date1, date2, nextTime;
161 1066421 double routingStep = 0.0, nextRuleTime, nextRoutingTime;
162
163
2/2
✓ Branch 0 taken 42480 times.
✓ Branch 1 taken 1023941 times.
1066421 if ( Nobjects[LINK] == 0 ) return fixedStep;
164
165 // --- find largest step possible if between routing events
166
4/4
✓ Branch 0 taken 209150 times.
✓ Branch 1 taken 814791 times.
✓ Branch 2 taken 35929 times.
✓ Branch 3 taken 173221 times.
1023941 if ( NumEvents > 0 && BetweenEvents )
167 {
168
2/2
✓ Branch 0 taken 35162 times.
✓ Branch 1 taken 767 times.
35929 nextTime = MIN(NewRunoffTime, ReportTime);
169 35929 date1 = getDateTime(NewRoutingTime);
170 35929 date2 = getDateTime(nextTime);
171
4/4
✓ Branch 0 taken 18622 times.
✓ Branch 1 taken 17307 times.
✓ Branch 2 taken 18067 times.
✓ Branch 3 taken 555 times.
35929 if ( date2 > date1 && date2 < Event[NextEvent].start )
172 {
173 18067 routingStep = (nextTime - NewRoutingTime) / 1000.0;
174 }
175 else
176 {
177 17862 date1 = getDateTime(NewRoutingTime + 1000.0 * fixedStep);
178
2/2
✓ Branch 0 taken 17847 times.
✓ Branch 1 taken 15 times.
17862 if ( date1 < Event[NextEvent].start ) return fixedStep;
179 }
180 }
181
182 // --- otherwise use a regular flow-routing based time step
183
2/2
✓ Branch 0 taken 988027 times.
✓ Branch 1 taken 18067 times.
1006094 if (routingStep == 0.0)
184 {
185 988027 routingStep = flowrout_getRoutingStep(routingModel, fixedStep);
186 }
187
188 // --- determine if control rule time interval reached
189
2/2
✓ Branch 0 taken 1921 times.
✓ Branch 1 taken 1004173 times.
1006094 if (RuleStep > 0)
190 {
191 1921 nextRuleTime = NewRuleTime + 1000. * RuleStep;
192 1921 nextRoutingTime = NewRoutingTime + 1000. * routingStep;
193
2/2
✓ Branch 0 taken 480 times.
✓ Branch 1 taken 1441 times.
1921 if (nextRoutingTime >= nextRuleTime)
194 {
195 480 routingStep = (nextRuleTime - NewRoutingTime) / 1000.0;
196 }
197 }
198 1006094 return routingStep;
199 }
200
201 //=============================================================================
202
203 1066421 void routing_execute(int routingModel, double routingStep)
204 //
205 // Input: routingModel = routing method code
206 // routingStep = routing time step (sec)
207 // Output: none
208 // Purpose: executes the routing process at the current time period.
209 //
210 {
211 1066421 int trialsCount = 1; // trials required to solve flow routing
212 1066421 int actionCount = 0; // number of control actions taken
213 1066421 int inSteadyState = TRUE; // system is in steady state
214 DateTime currentDate; // date at start of routing step
215 double stepFlowError; // 1 - (system outflow) / (system inflow)
216
217
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1066421 times.
1066421 if ( ErrorCode ) return;
218
219 // --- update mass balance totals over previous half time step
220 1066421 massbal_updateRoutingTotals(routingStep/2.);
221
222 // --- take any applicable control rule actions
223 1066421 currentDate = getDateTime(NewRoutingTime);
224 1066421 actionCount = evaluateControlRules(currentDate, routingStep);
225
226 // --- initialize mass balance and system inflow variables
227 1066421 stepFlowError = massbal_getStepFlowError();
228 1066421 massbal_initTimeStepTotals();
229 1066421 initSystemInflows();
230
231 // --- check that current date falls within a user-speficied event period
232 1066421 BetweenEvents = isBetweenEvents(currentDate);
233
2/2
✓ Branch 0 taken 1030493 times.
✓ Branch 1 taken 35928 times.
1066421 if (BetweenEvents == FALSE)
234 {
235 // --- apply current inflows to conveyance system
236 1030493 addSystemInflows(currentDate, routingStep);
237 1030493 inlet_findCapturedFlows(routingStep);
238
239 // --- route flows if system is not in steady state
240 1030493 inSteadyState = isInSteadyState(actionCount, stepFlowError);
241
1/2
✓ Branch 0 taken 1030493 times.
✗ Branch 1 not taken.
1030493 if (inSteadyState == FALSE)
242 1030493 trialsCount = routeFlow(routingModel, routingStep);
243
244 // --- route water quality constituents
245
3/4
✓ Branch 0 taken 201935 times.
✓ Branch 1 taken 828558 times.
✓ Branch 2 taken 201935 times.
✗ Branch 3 not taken.
1030493 if (Nobjects[POLLUT] > 0 && !IgnoreQuality)
246 {
247 201935 inlet_adjustQualInflows();
248 201935 qualrout_execute(routingStep);
249 }
250
251 // --- update mass balance totals for flows leaving the system
252 1030493 removeSystemOutflows(routingStep);
253 1030493 inlet_adjustQualOutflows();
254
255 // --- update time step & flow routing statistics
256
2/2
✓ Branch 0 taken 988013 times.
✓ Branch 1 taken 42480 times.
1030493 if (Nobjects[LINK] > 0)
257 {
258 988013 stats_updateFlowStats(routingStep, getDateTime(NewRoutingTime));
259 988013 stats_updateTimeStepStats(routingStep, trialsCount, inSteadyState);
260 }
261 }
262
263 // --- update mass balance totals over the current half time step
264 1066421 massbal_updateRoutingTotals(routingStep / 2.);
265 }
266
267 //=============================================================================
268
269 1066421 int evaluateControlRules(DateTime currentDate, double routingStep)
270 {
271 int j;
272 1066421 int actionCount = 0;
273
274 // --- find new link target settings that are not related to
275 // --- control rules (e.g., pump on/off depth limits)
276
2/2
✓ Branch 1 taken 56221721 times.
✓ Branch 2 taken 1066421 times.
57288142 for (j=0; j<Nobjects[LINK]; j++) link_setTargetSetting(j);
277
278 // --- evaluate control rules if next evaluation time reached
279
4/4
✓ Branch 0 taken 1921 times.
✓ Branch 1 taken 1064500 times.
✓ Branch 2 taken 480 times.
✓ Branch 3 taken 1441 times.
1066421 if (RuleStep == 0 || fabs(NewRoutingTime - NewRuleTime) < 1.0)
280 {
281 1064980 controls_evaluate(currentDate, currentDate - StartDateTime,
282 routingStep / SECperDAY);
283 }
284
285 // --- change each link's actual setting if it differs from its target
286
2/2
✓ Branch 0 taken 56221721 times.
✓ Branch 1 taken 1066421 times.
57288142 for (j=0; j<Nobjects[LINK]; j++)
287 {
288
2/2
✓ Branch 0 taken 102 times.
✓ Branch 1 taken 56221619 times.
56221721 if ( Link[j].targetSetting != Link[j].setting )
289 {
290 // --- update time when link was switched between open & closed
291
2/2
✓ Branch 0 taken 85 times.
✓ Branch 1 taken 17 times.
102 if ( Link[j].targetSetting * Link[j].setting == 0.0 )
292 85 Link[j].timeLastSet = currentDate;
293
294 // --- implement the change in the link's setting
295 102 link_setSetting(j, routingStep);
296 102 actionCount++;
297 }
298 }
299
300 // --- update value of elapsed routing time (in milliseconds)
301 1066421 OldRoutingTime = NewRoutingTime;
302 1066421 NewRoutingTime = NewRoutingTime + 1000.0 * routingStep;
303
304 // --- see if control rule evaluation time should be advanced
305
2/2
✓ Branch 0 taken 478 times.
✓ Branch 1 taken 1065943 times.
1066421 if (fabs(NewRoutingTime - (NewRuleTime + 1000.0*RuleStep)) < 1)
306 478 NewRuleTime += 1000.0 * RuleStep;
307 1066421 return actionCount;
308 }
309
310 //=============================================================================
311
312 1066421 void initSystemInflows()
313 {
314 int j;
315
316 // --- replace old water quality state with new state
317
2/2
✓ Branch 0 taken 237863 times.
✓ Branch 1 taken 828558 times.
1066421 if ( Nobjects[POLLUT] > 0 )
318 {
319
2/2
✓ Branch 1 taken 20940478 times.
✓ Branch 2 taken 237863 times.
21178341 for (j=0; j<Nobjects[NODE]; j++) node_setOldQualState(j);
320
2/2
✓ Branch 1 taken 21771441 times.
✓ Branch 2 taken 237863 times.
22009304 for (j=0; j<Nobjects[LINK]; j++) link_setOldQualState(j);
321 }
322
323 // --- set infiltration factor for storage unit seepage
324 // (-1 argument indicates global factor is used)
325 1066421 infil_setInfilFactor(-1);
326
327 // --- initialize lateral inflows at nodes
328
2/2
✓ Branch 0 taken 55932339 times.
✓ Branch 1 taken 1066421 times.
56998760 for (j = 0; j < Nobjects[NODE]; j++)
329 {
330 55932339 Node[j].oldLatFlow = Node[j].newLatFlow;
331 55932339 Node[j].newLatFlow = 0.0;
332 }
333 1066421 }
334
335 //=============================================================================
336
337 1066421 int isBetweenEvents(DateTime currentDate)
338 {
339 // --- if no events defined then result is always false
340
2/2
✓ Branch 0 taken 857271 times.
✓ Branch 1 taken 209150 times.
1066421 if ( NumEvents == 0 ) return FALSE;
341
342 // --- currrent event period has ended so result is true
343
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 209144 times.
209150 if ( currentDate > Event[NextEvent].end )
344 {
345 6 NextEvent++;
346 6 return TRUE;
347 }
348
349 // --- we've entered the next event period so result is false
350
2/2
✓ Branch 0 taken 173222 times.
✓ Branch 1 taken 35922 times.
209144 else if ( currentDate >= Event[NextEvent].start )
351 {
352 173222 return FALSE;
353 }
354 35922 return TRUE;
355 }
356
357 //=============================================================================
358
359 1030493 void addSystemInflows(DateTime currentDate, double routingStep)
360 {
361 int j;
362
363 // --- find evap. & seepage losses from storage nodes
364
2/2
✓ Branch 0 taken 37894467 times.
✓ Branch 1 taken 1030493 times.
38924960 for (j = 0; j < Nobjects[NODE]; j++)
365 37894467 Node[j].losses = node_getLosses(j, routingStep);
366
367 // --- add lateral inflows at nodes
368 1030493 addExternalInflows(currentDate);
369 1030493 addDryWeatherInflows(currentDate);
370 1030493 addWetWeatherInflows(OldRoutingTime);
371 1030493 addGroundwaterInflows(OldRoutingTime);
372 1030493 addLidDrainInflows(OldRoutingTime);
373 1030493 addRdiiInflows(currentDate);
374 1030493 addIfaceInflows(currentDate);
375
376 // --- initialize node inflow for quality routing
377
2/2
✓ Branch 0 taken 37894467 times.
✓ Branch 1 taken 1030493 times.
38924960 for (j = 0; j < Nobjects[NODE]; j++)
378
2/2
✓ Branch 0 taken 18599877 times.
✓ Branch 1 taken 19294590 times.
37894467 Node[j].qualInflow = MAX(0.0, Node[j].newLatFlow);
379 1030493 }
380
381 //=============================================================================
382
383 1030493 int isInSteadyState(int actionCount, double stepFlowError)
384 {
385 // --- check if can skip steady state periods based on flows
386
2/2
✓ Branch 0 taken 422 times.
✓ Branch 1 taken 1030071 times.
1030493 if ( SkipSteadyState )
387 {
388
1/2
✓ Branch 0 taken 422 times.
✗ Branch 1 not taken.
422 if ( OldRoutingTime == 0.0
389
2/2
✓ Branch 0 taken 412 times.
✓ Branch 1 taken 10 times.
422 || actionCount > 0
390
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 411 times.
412 || fabs(stepFlowError) > SysFlowTol
391
1/2
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
422 || inflowHasChanged() ) return FALSE;
392 else return TRUE;
393 }
394 1030071 return FALSE;
395 }
396
397 //=============================================================================
398
399 1030493 int routeFlow(int routingModel, double routingStep)
400 {
401 int j;
402 1030493 int stepCount = 1;
403
404 // --- replace old hydraulic state values with current ones
405
2/2
✓ Branch 1 taken 37208897 times.
✓ Branch 2 taken 1030493 times.
38239390 for (j = 0; j < Nobjects[LINK]; j++) link_setOldHydState(j);
406
2/2
✓ Branch 1 taken 37894467 times.
✓ Branch 2 taken 1030493 times.
38924960 for (j = 0; j < Nobjects[NODE]; j++) node_setOldHydState(j);
407
408 // --- initialize node inflows to lateral flows, outflows to evap +
409 // seepage losses, & overflows to excess stored volume
410
2/2
✓ Branch 0 taken 37894467 times.
✓ Branch 1 taken 1030493 times.
38924960 for (j = 0; j < Nobjects[NODE]; j++)
411 37894467 node_initFlows(j, routingStep);
412
413 // --- route flow through the drainage network
414
2/2
✓ Branch 0 taken 988013 times.
✓ Branch 1 taken 42480 times.
1030493 if ( Nobjects[LINK] > 0 )
415 {
416 988013 stepCount = flowrout_execute(SortedLinks, routingModel, routingStep);
417 }
418
419 // --- save overflows at inlet capture nodes as inlet backflow
420 1030493 return stepCount;
421 }
422
423 //=============================================================================
424
425 1030493 void removeSystemOutflows(double routingStep)
426 {
427 // --- remove evaporation, infiltration & outflows from system
428 1030493 removeStorageLosses(routingStep);
429 1030493 removeConduitLosses();
430 1030493 removeOutflows(routingStep);
431 1030493 }
432
433 //=============================================================================
434
435 1030493 void addExternalInflows(DateTime currentDate)
436 //
437 // Input: currentDate = current date/time
438 // Output: none
439 // Purpose: adds direct external inflows to nodes at current date.
440 //
441 {
442 int j, p;
443 double q, w;
444 TExtInflow* inflow;
445
446 // --- for each node with a defined external inflow
447
2/2
✓ Branch 0 taken 37894467 times.
✓ Branch 1 taken 1030493 times.
38924960 for (j = 0; j < Nobjects[NODE]; j++)
448 {
449 // --- get flow inflow
450 37894467 q = Node[j].apiExtInflow;
451 37894467 inflow = Node[j].extInflow;
452
2/2
✓ Branch 0 taken 1191123 times.
✓ Branch 1 taken 36726409 times.
37917532 while ( inflow )
453 {
454
2/2
✓ Branch 0 taken 1168058 times.
✓ Branch 1 taken 23065 times.
1191123 if ( inflow->type == FLOW_INFLOW )
455 {
456 1168058 q += inflow_getExtInflow(inflow, currentDate);
457 1168058 break;
458 }
459 23065 else inflow = inflow->next;
460 }
461
2/2
✓ Branch 0 taken 37426464 times.
✓ Branch 1 taken 468003 times.
37894467 if ( fabs(q) < FLOW_TOL ) q = 0.0;
462
463 // --- add flow inflow to node's lateral inflow
464 37894467 Node[j].newLatFlow += q;
465
1/2
✓ Branch 0 taken 37894467 times.
✗ Branch 1 not taken.
37894467 if (q >= 0.0)
466 37894467 massbal_addInflowFlow(EXTERNAL_INFLOW, q);
467 else
468 {
469 massbal_addOutflowFlow(-q, FALSE);
470 continue;
471 }
472
473 // --- add on any inflow (i.e., reverse flow) through an outfall
474
3/4
✓ Branch 0 taken 1309094 times.
✓ Branch 1 taken 36585373 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1309094 times.
37894467 if ( Node[j].type == OUTFALL && Node[j].oldNetInflow < 0.0 )
475 {
476 q = q - Node[j].oldNetInflow;
477 }
478
479 // --- get pollutant mass inflows
480 37894467 inflow = Node[j].extInflow;
481
2/2
✓ Branch 0 taken 1191123 times.
✓ Branch 1 taken 37894467 times.
39085590 while ( inflow )
482 {
483
2/2
✓ Branch 0 taken 23065 times.
✓ Branch 1 taken 1168058 times.
1191123 if ( inflow->type != FLOW_INFLOW )
484 {
485 23065 p = inflow->param;
486 23065 w = inflow_getExtInflow(inflow, currentDate);
487
1/2
✓ Branch 0 taken 23065 times.
✗ Branch 1 not taken.
23065 if ( inflow->type == CONCEN_INFLOW ) w *= q;
488 23065 Node[j].newQual[p] += w;
489 23065 massbal_addInflowQual(EXTERNAL_INFLOW, p, w);
490 }
491 1191123 inflow = inflow->next;
492 }
493 }
494 1030493 }
495
496 //=============================================================================
497
498 1030493 void addDryWeatherInflows(DateTime currentDate)
499 //
500 // Input: currentDate = current date/time
501 // Output: none
502 // Purpose: adds dry weather inflows to nodes at current date.
503 //
504 {
505 int j, p;
506 int month, day, hour;
507 double q, w;
508 TDwfInflow* inflow;
509
510 // --- get month (zero-based), day-of-week (zero-based),
511 // & hour-of-day for routing date/time
512 1030493 month = datetime_monthOfYear(currentDate) - 1;
513 1030493 day = datetime_dayOfWeek(currentDate) - 1;
514 1030493 hour = datetime_hourOfDay(currentDate);
515
516 // --- for each node with a defined dry weather inflow
517
2/2
✓ Branch 0 taken 37894467 times.
✓ Branch 1 taken 1030493 times.
38924960 for (j = 0; j < Nobjects[NODE]; j++)
518 {
519 37894467 inflow = Node[j].dwfInflow;
520
2/2
✓ Branch 0 taken 20877604 times.
✓ Branch 1 taken 17016863 times.
37894467 if ( !inflow ) continue;
521
522 // --- get flow inflow (i.e., the inflow whose param code is -1)
523 17016863 q = 0.0;
524
2/2
✓ Branch 0 taken 17016863 times.
✓ Branch 1 taken 422 times.
17017285 while ( inflow )
525 {
526
2/2
✓ Branch 0 taken 17016441 times.
✓ Branch 1 taken 422 times.
17016863 if ( inflow->param < 0 )
527 {
528 17016441 q = inflow_getDwfInflow(inflow, month, day, hour);
529 17016441 break;
530 }
531 422 inflow = inflow->next;
532 }
533
2/2
✓ Branch 0 taken 12709 times.
✓ Branch 1 taken 17004154 times.
17016863 if ( fabs(q) < FLOW_TOL ) q = 0.0;
534
535 // --- add flow inflow to node's lateral inflow
536 17016863 Node[j].newLatFlow += q;
537 17016863 massbal_addInflowFlow(DRY_WEATHER_INFLOW, q);
538
539 // --- stop if inflow is non-positive
540
2/2
✓ Branch 0 taken 12709 times.
✓ Branch 1 taken 17004154 times.
17016863 if ( q <= 0.0 ) continue;
541
542 // --- add default DWF pollutant inflows
543
2/2
✓ Branch 0 taken 246448 times.
✓ Branch 1 taken 17004154 times.
17250602 for ( p = 0; p < Nobjects[POLLUT]; p++)
544 {
545
2/2
✓ Branch 0 taken 61612 times.
✓ Branch 1 taken 184836 times.
246448 if ( Pollut[p].dwfConcen > 0.0 )
546 {
547 61612 w = q * Pollut[p].dwfConcen;
548 61612 Node[j].newQual[p] += w;
549 61612 massbal_addInflowQual(DRY_WEATHER_INFLOW, p, w);
550 }
551 }
552
553 // --- get pollutant mass inflows
554 17004154 inflow = Node[j].dwfInflow;
555
2/2
✓ Branch 0 taken 17004154 times.
✓ Branch 1 taken 17004154 times.
34008308 while ( inflow )
556 {
557
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 17004154 times.
17004154 if ( inflow->param >= 0 )
558 {
559 p = inflow->param;
560 w = q * inflow_getDwfInflow(inflow, month, day, hour);
561 Node[j].newQual[p] += w;
562 massbal_addInflowQual(DRY_WEATHER_INFLOW, p, w);
563
564 // --- subtract off any default inflow
565 if ( Pollut[p].dwfConcen > 0.0 )
566 {
567 w = q * Pollut[p].dwfConcen;
568 Node[j].newQual[p] -= w;
569 massbal_addInflowQual(DRY_WEATHER_INFLOW, p, -w);
570 }
571 }
572 17004154 inflow = inflow->next;
573 }
574 }
575 1030493 }
576
577 //=============================================================================
578
579 1030493 void addWetWeatherInflows(double routingTime)
580 //
581 // Input: routingTime = elasped time (millisec)
582 // Output: none
583 // Purpose: adds runoff inflows to nodes at current elapsed time.
584 //
585 {
586 int i, j, p;
587 double q, w;
588 double f;
589
590 // --- find where current routing time lies between latest runoff times
591
2/2
✓ Branch 0 taken 653142 times.
✓ Branch 1 taken 377351 times.
1030493 if ( Nobjects[SUBCATCH] == 0 ) return;
592 377351 f = (routingTime - OldRunoffTime) / (NewRunoffTime - OldRunoffTime);
593
2/2
✓ Branch 0 taken 7101 times.
✓ Branch 1 taken 370250 times.
377351 if ( f < 0.0 ) f = 0.0;
594
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 377351 times.
377351 if ( f > 1.0 ) f = 1.0;
595
596 // for each subcatchment outlet node,
597 // add interpolated runoff flow & pollutant load to node's inflow
598
2/2
✓ Branch 0 taken 2097083 times.
✓ Branch 1 taken 377351 times.
2474434 for (i = 0; i < Nobjects[SUBCATCH]; i++)
599 {
600 2097083 j = Subcatch[i].outNode;
601
2/2
✓ Branch 0 taken 2083801 times.
✓ Branch 1 taken 13282 times.
2097083 if ( j >= 0)
602 {
603 // add runoff flow to lateral inflow
604 2083801 q = subcatch_getWtdOutflow(i, f); // current runoff flow
605 2083801 Node[j].newLatFlow += q;
606 2083801 massbal_addInflowFlow(WET_WEATHER_INFLOW, q);
607
608 // add pollutant load
609
2/2
✓ Branch 0 taken 3354004 times.
✓ Branch 1 taken 2083801 times.
5437805 for (p = 0; p < Nobjects[POLLUT]; p++)
610 {
611 3354004 w = surfqual_getWtdWashoff(i, p, f);
612 3354004 Node[j].newQual[p] += w;
613 3354004 massbal_addInflowQual(WET_WEATHER_INFLOW, p, w);
614 }
615 }
616 }
617 }
618
619 //=============================================================================
620
621 1030493 void addGroundwaterInflows(double routingTime)
622 //
623 // Input: routingTime = elasped time (millisec)
624 // Output: none
625 // Purpose: adds groundwater inflows to nodes at current elapsed time.
626 //
627 {
628 int i, j, p;
629 double q, w;
630 double f;
631 TGroundwater* gw;
632
633 // --- find where current routing time lies between latest runoff times
634
2/2
✓ Branch 0 taken 653142 times.
✓ Branch 1 taken 377351 times.
1030493 if ( Nobjects[SUBCATCH] == 0 ) return;
635 377351 f = (routingTime - OldRunoffTime) / (NewRunoffTime - OldRunoffTime);
636
2/2
✓ Branch 0 taken 7101 times.
✓ Branch 1 taken 370250 times.
377351 if ( f < 0.0 ) f = 0.0;
637
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 377351 times.
377351 if ( f > 1.0 ) f = 1.0;
638
639 // --- for each subcatchment
640
2/2
✓ Branch 0 taken 2097083 times.
✓ Branch 1 taken 377351 times.
2474434 for (i = 0; i < Nobjects[SUBCATCH]; i++)
641 {
642 // --- see if subcatch contains groundwater
643 2097083 gw = Subcatch[i].groundwater;
644
2/2
✓ Branch 0 taken 183766 times.
✓ Branch 1 taken 1913317 times.
2097083 if ( gw )
645 {
646 // --- identify node receiving groundwater flow
647 183766 j = gw->node;
648
1/2
✓ Branch 0 taken 183766 times.
✗ Branch 1 not taken.
183766 if ( j >= 0 )
649 {
650 // add groundwater flow to lateral inflow
651 183766 q = ( (1.0 - f)*(gw->oldFlow) + f*(gw->newFlow) )
652 183766 * Subcatch[i].area;
653
2/2
✓ Branch 0 taken 13635 times.
✓ Branch 1 taken 170131 times.
183766 if ( fabs(q) < FLOW_TOL ) continue;
654 170131 Node[j].newLatFlow += q;
655 170131 massbal_addInflowFlow(GROUNDWATER_INFLOW, q);
656
657 // add pollutant load (for positive inflow)
658
2/2
✓ Branch 0 taken 136182 times.
✓ Branch 1 taken 33949 times.
170131 if ( q > 0.0 )
659 {
660
2/2
✓ Branch 0 taken 533408 times.
✓ Branch 1 taken 136182 times.
669590 for (p = 0; p < Nobjects[POLLUT]; p++)
661 {
662 533408 w = q * Pollut[p].gwConcen;
663 533408 Node[j].newQual[p] += w;
664 533408 massbal_addInflowQual(GROUNDWATER_INFLOW, p, w);
665 }
666 }
667 }
668 }
669 }
670 }
671
672 //=============================================================================
673
674 1030493 void addLidDrainInflows(double routingTime)
675 //
676 // Input: routingTime = elasped time (millisec)
677 // Output: none
678 // Purpose: adds inflows to nodes receiving LID drain flow.
679 //
680 {
681 int j;
682 double f;
683
684 // for each subcatchment
685
2/2
✓ Branch 0 taken 653142 times.
✓ Branch 1 taken 377351 times.
1030493 if ( Nobjects[SUBCATCH] == 0 ) return;
686 377351 f = (routingTime - OldRunoffTime) / (NewRunoffTime - OldRunoffTime);
687
2/2
✓ Branch 0 taken 7101 times.
✓ Branch 1 taken 370250 times.
377351 if ( f < 0.0 ) f = 0.0;
688
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 377351 times.
377351 if ( f > 1.0 ) f = 1.0;
689
2/2
✓ Branch 0 taken 2097083 times.
✓ Branch 1 taken 377351 times.
2474434 for (j = 0; j < Nobjects[SUBCATCH]; j++)
690 {
691
3/4
✓ Branch 0 taken 2097083 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 48241 times.
✓ Branch 3 taken 2048842 times.
2097083 if ( Subcatch[j].area > 0.0 && Subcatch[j].lidArea > 0.0 )
692 48241 lid_addDrainInflow(j, f);
693 }
694 }
695
696 //=============================================================================
697
698 1030493 void addRdiiInflows(DateTime currentDate)
699 //
700 // Input: currentDate = current date/time
701 // Output: none
702 // Purpose: adds RDII inflows to nodes at current date.
703 //
704 {
705 int i, j, p;
706 double q, w;
707 int numRdiiNodes;
708
709 // --- see if any nodes have RDII at current date
710 1030493 numRdiiNodes = rdii_getNumRdiiFlows(currentDate);
711
712 // --- add RDII flow to each node's lateral inflow
713
2/2
✓ Branch 0 taken 10971860 times.
✓ Branch 1 taken 1030493 times.
12002353 for (i=0; i<numRdiiNodes; i++)
714 {
715 10971860 rdii_getRdiiFlow(i, &j, &q);
716
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 10971860 times.
10971860 if ( j < 0 ) continue;
717
2/2
✓ Branch 0 taken 7223500 times.
✓ Branch 1 taken 3748360 times.
10971860 if ( fabs(q) < FLOW_TOL ) continue;
718 3748360 Node[j].newLatFlow += q;
719 3748360 massbal_addInflowFlow(RDII_INFLOW, q);
720
721 // add pollutant load (for positive inflow)
722
1/2
✓ Branch 0 taken 3748360 times.
✗ Branch 1 not taken.
3748360 if ( q > 0.0 )
723 {
724
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3748360 times.
3748360 for (p = 0; p < Nobjects[POLLUT]; p++)
725 {
726 w = q * Pollut[p].rdiiConcen;
727 Node[j].newQual[p] += w;
728 massbal_addInflowQual(RDII_INFLOW, p, w);
729 }
730 }
731 }
732 1030493 }
733
734 //=============================================================================
735
736 1030493 void addIfaceInflows(DateTime currentDate)
737 //
738 // Input: currentDate = current date/time
739 // Output: none
740 // Purpose: adds inflows from routing interface file to nodes at current date.
741 //
742 {
743 int i, j, p;
744 double q, w;
745 int numIfaceNodes;
746
747 // --- see if any nodes have interface inflows at current date
748
2/2
✓ Branch 0 taken 1026172 times.
✓ Branch 1 taken 4321 times.
1030493 if ( Finflows.mode != USE_FILE ) return;
749 4321 numIfaceNodes = iface_getNumIfaceNodes(currentDate);
750
751 // --- add interface flow to each node's lateral inflow
752
2/2
✓ Branch 0 taken 4321 times.
✓ Branch 1 taken 4321 times.
8642 for (i=0; i<numIfaceNodes; i++)
753 {
754 4321 j = iface_getIfaceNode(i);
755
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4321 times.
4321 if ( j < 0 ) continue;
756 4321 q = iface_getIfaceFlow(i);
757
2/2
✓ Branch 0 taken 845 times.
✓ Branch 1 taken 3476 times.
4321 if ( fabs(q) < FLOW_TOL ) continue;
758 3476 Node[j].newLatFlow += q;
759 3476 massbal_addInflowFlow(EXTERNAL_INFLOW, q);
760
761 // add pollutant load (for positive inflow)
762
1/2
✓ Branch 0 taken 3476 times.
✗ Branch 1 not taken.
3476 if ( q > 0.0 )
763 {
764
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 3476 times.
3476 for (p = 0; p < Nobjects[POLLUT]; p++)
765 {
766 w = q * iface_getIfaceQual(i, p);
767 Node[j].newQual[p] += w;
768 massbal_addInflowQual(EXTERNAL_INFLOW, p, w);
769 }
770 }
771 }
772 }
773
774 //=============================================================================
775
776 1 int inflowHasChanged()
777 //
778 // Input: none
779 // Output: returns TRUE if external inflows or outfall flows have changed
780 // from the previous time step
781 // Purpose: checks if the hydraulic state of the system has changed from
782 // the previous time step.
783 //
784 {
785 int j;
786 double diff, qOld, qNew;
787
788 // --- check if external inflows or outfall flows have changed
789
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
1 for (j = 0; j < Nobjects[NODE]; j++)
790 {
791 1 qOld = Node[j].oldLatFlow;
792 1 qNew = Node[j].newLatFlow;
793
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
1 if ( fabs(qOld) > TINY ) diff = (qNew / qOld) - 1.0;
794
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
1 else if ( fabs(qNew) > TINY ) diff = 1.0;
795 else diff = 0.0;
796
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
1 if ( fabs(diff) > LatFlowTol ) return TRUE;
797 if ( Node[j].type == OUTFALL || Node[j].degree == 0 )
798 {
799 qOld = Node[j].oldFlowInflow;
800 qNew = Node[j].inflow;
801 if ( fabs(qOld) > TINY ) diff = (qNew / qOld) - 1.0;
802 else if ( fabs(qNew) > TINY ) diff = 1.0;
803 else diff = 0.0;
804 if ( fabs(diff) > LatFlowTol ) return TRUE;
805 }
806 }
807 return FALSE;
808 }
809
810 //=============================================================================
811
812 1030493 void removeStorageLosses(double tStep)
813 //
814 // Input: tStep = routing time step (sec)
815 // Output: none
816 // Purpose: adds flow rate lost from all storage nodes due to evaporation
817 // & seepage in current time step to overall mass balance totals.
818 //
819 {
820 int i;
821 1030493 double evapLoss = 0.0,
822 1030493 exfilLoss = 0.0;
823
824 // --- check each storage node
825
2/2
✓ Branch 0 taken 37894467 times.
✓ Branch 1 taken 1030493 times.
38924960 for ( i = 0; i < Nobjects[NODE]; i++ )
826 {
827
2/2
✓ Branch 0 taken 143130 times.
✓ Branch 1 taken 37751337 times.
37894467 if (Node[i].type == STORAGE)
828 {
829 // --- update total system storage losses
830 143130 evapLoss += Storage[Node[i].subIndex].evapLoss;
831 143130 exfilLoss += Storage[Node[i].subIndex].exfilLoss;
832 }
833 }
834
835 // --- add loss rates (ft3/sec) to time step's mass balance
836 1030493 massbal_addNodeLosses(evapLoss/tStep, exfilLoss/tStep);
837 1030493 }
838
839 //=============================================================================
840
841 1030493 void removeConduitLosses()
842 //
843 // Input: none
844 // Output: none
845 // Purpose: adds flow rate lost from all conduits due to evaporation
846 // & seepage over current time step to overall mass balance.
847 //
848 {
849 int i, k;
850 double barrels,
851 1030493 evapLoss = 0.0,
852 1030493 seepLoss = 0.0;
853
854
2/2
✓ Branch 0 taken 37208897 times.
✓ Branch 1 taken 1030493 times.
38239390 for ( i = 0; i < Nobjects[LINK]; i++ )
855 {
856
2/2
✓ Branch 0 taken 36789022 times.
✓ Branch 1 taken 419875 times.
37208897 if (Link[i].type == CONDUIT)
857 {
858 // --- retrieve number of barrels
859 36789022 k = Link[i].subIndex;
860 36789022 barrels = Conduit[k].barrels;
861
862 // --- update total conduit losses
863 36789022 evapLoss += Conduit[k].evapLossRate * barrels;
864 36789022 seepLoss += Conduit[k].seepLossRate * barrels;
865 }
866 }
867 1030493 massbal_addLinkLosses(evapLoss, seepLoss);
868 1030493 }
869
870 //=============================================================================
871
872 1030493 void removeOutflows(double tStep)
873 //
874 // Input: none
875 // Output: none
876 // Purpose: finds flows that leave the system and adds these to mass
877 // balance totals.
878 //
879 {
880 int i, p, k;
881 int isFlooded;
882 double q, w, v;
883
884
2/2
✓ Branch 0 taken 37894467 times.
✓ Branch 1 taken 1030493 times.
38924960 for ( i = 0; i < Nobjects[NODE]; i++ )
885 {
886 // --- accumulate inflow volume & pollut. load at outfalls
887
4/4
✓ Branch 0 taken 1309094 times.
✓ Branch 1 taken 36585373 times.
✓ Branch 2 taken 655306 times.
✓ Branch 3 taken 653788 times.
37894467 if ( Node[i].type == OUTFALL && Node[i].inflow > 0.0 )
888 {
889 655306 k = Node[i].subIndex;
890
2/2
✓ Branch 0 taken 1405 times.
✓ Branch 1 taken 653901 times.
655306 if ( Outfall[k].routeTo >= 0 )
891 {
892 1405 v = Node[i].inflow * tStep;
893 1405 Outfall[k].vRouted += v;
894
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1405 times.
1405 for (p = 0; p < Nobjects[POLLUT]; p++)
895 Outfall[k].wRouted[p] += Node[i].newQual[p] * v;
896 }
897 }
898
899 // --- update mass balance with flow and mass leaving the system
900 // through outfalls and flooded interior nodes
901 37894467 q = node_getSystemOutflow(i, &isFlooded);
902
2/2
✓ Branch 0 taken 693070 times.
✓ Branch 1 taken 37201397 times.
37894467 if ( q > 0.0 )
903 {
904 693070 massbal_addOutflowFlow(q, isFlooded);
905
2/2
✓ Branch 0 taken 399453 times.
✓ Branch 1 taken 693070 times.
1092523 for ( p = 0; p < Nobjects[POLLUT]; p++ )
906 {
907 399453 w = q * Node[i].newQual[p];
908 399453 massbal_addOutflowQual(p, w, isFlooded);
909 }
910 }
911 37201397 else massbal_addInflowFlow(EXTERNAL_INFLOW, -q);
912
913 // --- update mass balance with mass leaving system through negative
914 // lateral inflows (lateral flow was previously accounted for)
915 37894467 q = Node[i].newLatFlow;
916
2/2
✓ Branch 0 taken 20452 times.
✓ Branch 1 taken 37874015 times.
37894467 if ( q < 0.0 )
917 {
918
2/2
✓ Branch 0 taken 1851 times.
✓ Branch 1 taken 20452 times.
22303 for ( p = 0; p < Nobjects[POLLUT]; p++ )
919 {
920 1851 w = -q * Node[i].newQual[p];
921 1851 massbal_addOutflowQual(p, w, FALSE);
922 }
923 }
924 }
925 1030493 }
926
927 //=============================================================================
928
929 9 void sortEvents()
930 //
931 // Input: none
932 // Output: none
933 // Purpose: sorts the entries of the Event array in chronological order.
934 //
935 {
936 int i, j, gap;
937 TEvent temp;
938
939 // Apply shell sort to event list
940
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 9 times.
11 for (gap = NumEvents/2; gap >= 1; gap /= 2)
941 {
942
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 2 times.
6 for (i = gap; i < NumEvents; i += gap)
943 {
944 4 temp = Event[i];
945 4 j = i - gap;
946
2/4
✓ Branch 0 taken 4 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 4 times.
4 while (j >= 0 && Event[j].start > temp.start)
947 {
948 Event[j+gap] = Event[j];
949 j -= gap;
950 }
951
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
4 if (j != i-gap) Event[j+gap] = temp;
952 }
953 }
954
955 // Adjust for overlapping events
956
2/2
✓ Branch 0 taken 4 times.
✓ Branch 1 taken 9 times.
13 for (i = 0; i < NumEvents-1; i++)
957 {
958
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 4 times.
4 if ( Event[i].end > Event[i+1].start ) Event[i].end = Event[i+1].start;
959 }
960 9 }
961
962 //=============================================================================
963