GCC Code Coverage Report


Directory: src/solver/
Coverage: low: ≥ 0% medium: ≥ 75.0% high: ≥ 90.0%
Coverage Exec / Excl / Total
Lines: 92.3% 131 / 0 / 142
Functions: 100.0% 10 / 0 / 10
Branches: 73.4% 94 / 0 / 128

treatmnt.c
Line Branch Exec Source
1 //-----------------------------------------------------------------------------
2 // treatmnt.c
3 //
4 // Project: EPA SWMM5
5 // Version: 5.2
6 // Date: 11/01/21 (Build 5.2.0)
7 // Author: L. Rossman
8 //
9 // Pollutant treatment functions.
10 //
11 // Update History
12 // ==============
13 // Build 5.1.008:
14 // - A bug in evaluating recursive calls to treatment functions was fixed.
15 // Build 5.2.0:
16 // - Changed enumerated constant used to indicate a math expression error.
17 //-----------------------------------------------------------------------------
18 #define _CRT_SECURE_NO_DEPRECATE
19
20 #include <stdlib.h>
21 #include <string.h>
22 #include "headers.h"
23
24 //-----------------------------------------------------------------------------
25 // Constants
26 //-----------------------------------------------------------------------------
27 static const int PVMAX = 5; // number of process variables
28 enum ProcessVarType {pvHRT, // hydraulic residence time
29 pvDT, // time step duration
30 pvFLOW, // flow rate
31 pvDEPTH, // water height above invert
32 pvAREA}; // storage surface area
33
34 //-----------------------------------------------------------------------------
35 // Shared variables
36 //-----------------------------------------------------------------------------
37 static int ErrCode; // treatment error code
38 static int J; // index of node being analyzed
39 static double Dt; // curent time step (sec)
40 static double Q; // node inflow (cfs)
41 static double V; // node volume (ft3)
42 static double* R; // array of pollut. removals
43 static double* Cin; // node inflow concentrations
44
45 //-----------------------------------------------------------------------------
46 // External functions (declared in funcs.h)
47 //-----------------------------------------------------------------------------
48 // treatmnt_open (called from routing_open)
49 // treatment_close (called from routing_close)
50 // treatmnt_readExpression (called from parseLine in input.c)
51 // treatmnt_delete (called from deleteObjects in project.c)
52 // treatmnt_setInflow (called from qualrout_execute)
53 // treatmnt_treat (called from findNodeQual in qualrout.c)
54
55 //-----------------------------------------------------------------------------
56 // Local functions
57 //-----------------------------------------------------------------------------
58 static int createTreatment(int node);
59 static double getRemoval(int pollut);
60 static int getVariableIndex(char* s);
61 static double getVariableValue(int varCode);
62
63
64 //=============================================================================
65
66 58 int treatmnt_open(void)
67 //
68 // Input: none
69 // Output: returns TRUE if successful, FALSE if not
70 // Purpose: allocates memory for computing pollutant removals by treatment.
71 //
72 {
73 58 R = NULL;
74 58 Cin = NULL;
75
2/2
✓ Branch 0 taken 15 times.
✓ Branch 1 taken 43 times.
58 if ( Nobjects[POLLUT] > 0 )
76 {
77 15 R = (double *) calloc(Nobjects[POLLUT], sizeof(double));
78 15 Cin = (double *) calloc(Nobjects[POLLUT], sizeof(double));
79
2/4
✓ Branch 0 taken 15 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 15 times.
15 if ( R == NULL || Cin == NULL)
80 {
81 report_writeErrorMsg(ERR_MEMORY, "");
82 return FALSE;
83 }
84 }
85 58 return TRUE;
86 }
87
88 //=============================================================================
89
90 58 void treatmnt_close(void)
91 //
92 // Input: none
93 // Output: returns an error code
94 // Purpose: frees memory used for computing pollutant removals by treatment.
95 //
96 {
97
2/2
✓ Branch 0 taken 15 times.
✓ Branch 1 taken 43 times.
58 FREE(R);
98
2/2
✓ Branch 0 taken 15 times.
✓ Branch 1 taken 43 times.
58 FREE(Cin);
99 58 }
100
101 //=============================================================================
102
103 15 int treatmnt_readExpression(char* tok[], int ntoks)
104 //
105 // Input: tok[] = array of string tokens
106 // ntoks = number of tokens
107 // Output: returns an error code
108 // Purpose: reads a treatment expression from a tokenized line of input.
109 //
110 {
111 char s[MAXLINE+1];
112 char* expr;
113 int i, j, k, p;
114 MathExpr* equation; // ptr. to a math. expression
115
116 // --- retrieve node & pollutant
117
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
15 if ( ntoks < 3 ) return error_setInpError(ERR_ITEMS, "");
118 15 j = project_findObject(NODE, tok[0]);
119
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
15 if ( j < 0 ) return error_setInpError(ERR_NAME, tok[0]);
120 15 p = project_findObject(POLLUT, tok[1]);
121
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
15 if ( p < 0 ) return error_setInpError(ERR_NAME, tok[1]);
122
123 // --- concatenate remaining tokens into a single string
124 15 sstrncpy(s, tok[2], MAXLINE);
125
2/2
✓ Branch 0 taken 44 times.
✓ Branch 1 taken 15 times.
59 for ( i=3; i<ntoks; i++)
126 {
127 44 sstrcat(s, " ", MAXLINE);
128 44 sstrcat(s, tok[i], MAXLINE);
129 }
130
131 // --- check treatment type
132
3/6
✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 5 times.
✓ Branch 5 taken 10 times.
15 if ( UCHAR(s[0]) == 'R' ) k = 0;
133
2/6
✗ Branch 0 not taken.
✓ Branch 1 taken 10 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 10 times.
✗ Branch 5 not taken.
10 else if ( UCHAR(s[0]) == 'C' ) k = 1;
134 else return error_setInpError(ERR_KEYWORD, tok[2]);
135
136 // --- start treatment expression after equals sign
137 15 expr = strchr(s, '=');
138
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
15 if ( expr == NULL ) return error_setInpError(ERR_KEYWORD, "");
139 15 else expr++;
140
141 // --- create treatment objects at node j if they don't already exist
142
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 10 times.
15 if ( Node[j].treatment == NULL )
143 {
144
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 5 times.
5 if ( !createTreatment(j) ) return error_setInpError(ERR_MEMORY, "");
145 }
146
147 // --- create a parsed expression tree from the string expr
148 // (getVariableIndex is the function that converts a treatment
149 // variable's name into an index number)
150 15 equation = mathexpr_create(expr, getVariableIndex);
151
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 15 times.
15 if ( equation == NULL )
152 return error_setInpError(ERR_MATH_EXPR, "");
153
154 // --- save the treatment parameters in the node's treatment object
155
1/2
✓ Branch 0 taken 15 times.
✗ Branch 1 not taken.
15 if (Node[j].treatment != NULL)
156 {
157 15 Node[j].treatment[p].treatType = k;
158 15 Node[j].treatment[p].equation = equation;
159 }
160 15 return 0;
161 }
162
163 //=============================================================================
164
165 10183 void treatmnt_delete(int j)
166 //
167 // Input: j = node index
168 // Output: none
169 // Purpose: deletes the treatment objects for each pollutant at a node.
170 //
171 {
172 int p;
173
2/2
✓ Branch 0 taken 5 times.
✓ Branch 1 taken 10178 times.
10183 if ( Node[j].treatment )
174 {
175
2/2
✓ Branch 0 taken 22 times.
✓ Branch 1 taken 5 times.
27 for (p=0; p<Nobjects[POLLUT]; p++)
176 22 mathexpr_delete(Node[j].treatment[p].equation);
177 5 free(Node[j].treatment);
178 }
179 10183 Node[j].treatment = NULL;
180 10183 }
181
182 //=============================================================================
183
184 88564 void treatmnt_setInflow(double qIn, double wIn[])
185 //
186 // Input: j = node index
187 // qIn = flow inflow rate (cfs)
188 // wIn = pollutant mass inflow rate (mass/sec)
189 // Output: none
190 // Purpose: computes and saves array of inflow concentrations to a node.
191 //
192 {
193 int p;
194
2/2
✓ Branch 0 taken 76556 times.
✓ Branch 1 taken 12008 times.
88564 if ( qIn > 0.0 )
195
2/2
✓ Branch 0 taken 159472 times.
✓ Branch 1 taken 76556 times.
236028 for (p = 0; p < Nobjects[POLLUT]; p++) Cin[p] = wIn[p]/qIn;
196 else
197
2/2
✓ Branch 0 taken 24148 times.
✓ Branch 1 taken 12008 times.
36156 for (p = 0; p < Nobjects[POLLUT]; p++) Cin[p] = 0.0;
198 88564 }
199
200 //=============================================================================
201
202 88564 void treatmnt_treat(int j, double q, double v, double tStep)
203 //
204 // Input: j = node index
205 // q = inflow to node (cfs)
206 // v = volume of node (ft3)
207 // tStep = routing time step (sec)
208 // Output: none
209 // Purpose: updates pollutant concentrations at a node after treatment.
210 //
211 {
212 int p; // pollutant index
213 double cOut; // concentration after treatment
214 double massLost; // mass lost by treatment per time step
215 TTreatment* treatment; // pointer to treatment object
216
217 // --- set locally shared variables for node j
218
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 88564 times.
88564 if ( Node[j].treatment == NULL ) return;
219 88564 ErrCode = 0;
220 88564 J = j; // current node
221 88564 Dt = tStep; // current time step
222 88564 Q = q; // current inflow rate
223 88564 V = v; // current node volume
224
225 // --- initialze each removal to indicate no value
226
2/2
✓ Branch 0 taken 183620 times.
✓ Branch 1 taken 88564 times.
272184 for ( p = 0; p < Nobjects[POLLUT]; p++) R[p] = -1.0;
227
228 // --- determine removal of each pollutant
229
2/2
✓ Branch 0 taken 183620 times.
✓ Branch 1 taken 88564 times.
272184 for ( p = 0; p < Nobjects[POLLUT]; p++)
230 {
231 // --- removal is zero if there is no treatment equation
232 183620 treatment = &Node[j].treatment[p];
233
2/2
✓ Branch 0 taken 3787 times.
✓ Branch 1 taken 179833 times.
183620 if ( treatment->equation == NULL ) R[p] = 0.0;
234
235 // --- no removal for removal-type expression when there is no inflow
236
4/4
✓ Branch 0 taken 174423 times.
✓ Branch 1 taken 5410 times.
✓ Branch 2 taken 23961 times.
✓ Branch 3 taken 150462 times.
179833 else if ( treatment->treatType == REMOVAL && q <= ZERO ) R[p] = 0.0;
237
238 // --- otherwise evaluate the treatment expression to find R[p]
239 155872 else getRemoval(p);
240 }
241
242 // --- check for error condition
243
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 88564 times.
88564 if ( ErrCode == ERR_CYCLIC_TREATMENT )
244 {
245 report_writeErrorMsg(ERR_CYCLIC_TREATMENT, Node[J].ID);
246 }
247
248 // --- update nodal concentrations and mass balances
249
2/2
✓ Branch 0 taken 183620 times.
✓ Branch 1 taken 88564 times.
272184 else for ( p = 0; p < Nobjects[POLLUT]; p++ )
250 {
251
2/2
✓ Branch 0 taken 56649 times.
✓ Branch 1 taken 126971 times.
183620 if ( R[p] == 0.0 ) continue;
252 126971 treatment = &Node[j].treatment[p];
253
254 // --- removal-type treatment equations get applied to inflow stream
255
256
2/2
✓ Branch 0 taken 121881 times.
✓ Branch 1 taken 5090 times.
126971 if ( treatment->treatType == REMOVAL )
257 {
258 // --- if no pollutant in inflow then cOut is current nodal concen.
259
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 121881 times.
121881 if ( Cin[p] == 0.0 ) cOut = Node[j].newQual[p];
260
261 // --- otherwise apply removal to influent concen.
262 121881 else cOut = (1.0 - R[p]) * Cin[p];
263
264 // --- cOut can't be greater than mixture concen. at node
265 // (i.e., in case node is a storage unit)
266
2/2
✓ Branch 0 taken 121878 times.
✓ Branch 1 taken 3 times.
121881 cOut = MIN(cOut, Node[j].newQual[p]);
267 }
268
269 // --- concentration-type equations get applied to nodal concentration
270 else
271 {
272 5090 cOut = (1.0 - R[p]) * Node[j].newQual[p];
273 }
274
275 // --- mass lost must account for any initial mass in storage
276 126971 massLost = (Cin[p]*q*tStep + Node[j].oldQual[p]*Node[j].oldVolume -
277 126971 cOut*(q*tStep + Node[j].oldVolume)) / tStep;
278
2/2
✓ Branch 0 taken 3 times.
✓ Branch 1 taken 126968 times.
126971 massLost = MAX(0.0, massLost);
279
280 // --- add mass loss to mass balance totals and revise nodal concentration
281 126971 massbal_addReactedMass(p, massLost);
282 126971 Node[j].newQual[p] = cOut;
283 }
284 }
285
286 //=============================================================================
287
288 5 int createTreatment(int j)
289 //
290 // Input: j = node index
291 // Output: returns TRUE if successful, FALSE if not
292 // Purpose: creates a treatment object for each pollutant at a node.
293 //
294 {
295 int p;
296 5 Node[j].treatment = (TTreatment *) calloc(Nobjects[POLLUT],
297 sizeof(TTreatment));
298
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5 times.
5 if ( Node[j].treatment == NULL )
299 {
300 return FALSE;
301 }
302
2/2
✓ Branch 0 taken 22 times.
✓ Branch 1 taken 5 times.
27 for (p = 0; p < Nobjects[POLLUT]; p++)
303 {
304 22 Node[j].treatment[p].equation = NULL;
305 }
306 5 return TRUE;
307 }
308
309 //=============================================================================
310
311 19 int getVariableIndex(char* s)
312 //
313 // Input: s = name of a process variable or pollutant
314 // Output: returns index of process variable or pollutant
315 // Purpose: finds position of process variable/pollutant in list of names.
316 //
317 {
318 // --- check for a process variable first
319 int k;
320 19 int m = PVMAX; // PVMAX is number of process variables
321
322 19 k = findmatch(s, ProcessVarWords);
323
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 13 times.
19 if ( k >= 0 ) return k;
324
325 // --- then check for a pollutant concentration
326 13 k = project_findObject(POLLUT, s);
327
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 1 time.
13 if ( k >= 0 ) return (k + m);
328
329 // --- finally check for a pollutant removal
330
3/8
✗ Branch 0 not taken.
✓ Branch 1 taken 1 time.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 1 time.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 time.
✗ Branch 7 not taken.
1 if ( UCHAR(s[0]) == 'R' && s[1] == '_')
331 {
332 1 k = project_findObject(POLLUT, s+2);
333
1/2
✓ Branch 0 taken 1 time.
✗ Branch 1 not taken.
1 if ( k >= 0 ) return (Nobjects[POLLUT] + k + m);
334 }
335 return -1;
336 }
337
338 //=============================================================================
339
340 9860 double getVariableValue(int varCode)
341 //
342 // Input: varCode = code number of process variable or pollutant
343 // Output: returns current value of variable
344 // Purpose: finds current value of a process variable or pollutant concen.,
345 // making reference to the node being evaluated which is stored in
346 // shared variable J.
347 //
348 {
349 int p;
350 double a1, a2, y;
351 TTreatment* treatment;
352
353 // --- variable is a process variable
354
2/2
✓ Branch 0 taken 3180 times.
✓ Branch 1 taken 6680 times.
9860 if ( varCode < PVMAX )
355 {
356
5/6
✓ Branch 0 taken 1060 times.
✓ Branch 1 taken 530 times.
✓ Branch 2 taken 530 times.
✓ Branch 3 taken 530 times.
✓ Branch 4 taken 530 times.
✗ Branch 5 not taken.
3180 switch ( varCode )
357 {
358 1060 case pvHRT: // HRT in hours
359
1/2
✓ Branch 0 taken 1060 times.
✗ Branch 1 not taken.
1060 if ( Node[J].type == STORAGE )
360 {
361 1060 return Storage[Node[J].subIndex].hrt / 3600.0;
362 }
363 else return 0.0;
364
365 530 case pvDT:
366 530 return Dt; // time step in seconds
367
368 530 case pvFLOW:
369 530 return Q * UCF(FLOW); // flow in user's units
370
371 530 case pvDEPTH:
372 530 y = (Node[J].oldDepth + Node[J].newDepth) / 2.0;
373 530 return y * UCF(LENGTH); // depth in ft or m
374
375 530 case pvAREA:
376 530 a1 = node_getSurfArea(J, Node[J].oldDepth);
377 530 a2 = node_getSurfArea(J, Node[J].newDepth);
378 530 return (a1 + a2) / 2.0 * UCF(LENGTH) * UCF(LENGTH);
379
380 default: return 0.0;
381 }
382 }
383
384 // --- variable is a pollutant concentration
385
2/2
✓ Branch 0 taken 6150 times.
✓ Branch 1 taken 530 times.
6680 else if ( varCode < PVMAX + Nobjects[POLLUT] )
386 {
387 6150 p = varCode - PVMAX;
388 6150 treatment = &Node[J].treatment[p];
389
2/2
✓ Branch 0 taken 1060 times.
✓ Branch 1 taken 5090 times.
6150 if ( treatment->treatType == REMOVAL ) return Cin[p];
390 5090 return Node[J].newQual[p];
391 }
392
393 // --- variable is a pollutant removal
394 else
395 {
396 530 p = varCode - PVMAX - Nobjects[POLLUT];
397
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 530 times.
530 if ( p >= Nobjects[POLLUT] ) return 0.0;
398 530 return getRemoval(p);
399 }
400 }
401
402 //=============================================================================
403
404 156402 double getRemoval(int p)
405 //
406 // Input: p = pollutant index
407 // Output: returns fractional removal of pollutant
408 // Purpose: computes removal of a specific pollutant
409 //
410 {
411 156402 double c0 = Node[J].newQual[p]; // initial node concentration
412 double r; // removal value
413 TTreatment* treatment;
414
415 // --- case where removal already being computed for another pollutant
416
2/4
✓ Branch 0 taken 156402 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 156402 times.
156402 if ( R[p] > 1.0 || ErrCode )
417 {
418 ErrCode = 1;
419 return 0.0;
420 }
421
422 // --- case where removal already computed
423
3/4
✓ Branch 0 taken 530 times.
✓ Branch 1 taken 155872 times.
✓ Branch 2 taken 530 times.
✗ Branch 3 not taken.
156402 if ( R[p] >= 0.0 && R[p] <= 1.0 ) return R[p];
424
425 // --- set R[p] to value > 1 to show that value is being sought
426 // (prevents infinite recursive calls in case two removals
427 // depend on each other)
428 155872 R[p] = 10.0;
429
430 // --- case where current concen. is zero
431
2/2
✓ Branch 0 taken 28900 times.
✓ Branch 1 taken 126972 times.
155872 if ( c0 == 0.0 )
432 {
433 28900 R[p] = 0.0;
434 28900 return 0.0;
435 }
436
437 // --- apply treatment eqn.
438 126972 treatment = &Node[J].treatment[p];
439 126972 r = mathexpr_eval(treatment->equation, getVariableValue);
440
2/2
✓ Branch 0 taken 1 time.
✓ Branch 1 taken 126971 times.
126972 r = MAX(0.0, r);
441
442 // --- case where treatment eqn. is for removal
443
2/2
✓ Branch 0 taken 121882 times.
✓ Branch 1 taken 5090 times.
126972 if ( treatment->treatType == REMOVAL )
444 {
445
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 121882 times.
121882 r = MIN(1.0, r);
446 121882 R[p] = r;
447 }
448
449 // --- case where treatment eqn. is for effluent concen.
450 else
451 {
452
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 5090 times.
5090 r = MIN(c0, r);
453 5090 R[p] = 1.0 - r/c0;
454 }
455 126972 return R[p];
456 }
457