kinwave.c
| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | //----------------------------------------------------------------------------- | ||
| 2 | // kinwave.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 | // Kinematic wave flow routing functions. | ||
| 11 | // | ||
| 12 | // Update History | ||
| 13 | // ============== | ||
| 14 | // Build 5.1.008: | ||
| 15 | // - Conduit inflow passed to function that computes conduit losses. | ||
| 16 | // Build 5.1.014: | ||
| 17 | // - Arguments to function link_getLossRate changed. | ||
| 18 | // Build 5.2.4: | ||
| 19 | // - Arguments to function link_getLossRate changed again. | ||
| 20 | //----------------------------------------------------------------------------- | ||
| 21 | #define _CRT_SECURE_NO_DEPRECATE | ||
| 22 | |||
| 23 | #include <math.h> | ||
| 24 | #include "headers.h" | ||
| 25 | #include "findroot.h" | ||
| 26 | |||
| 27 | //----------------------------------------------------------------------------- | ||
| 28 | // Constants | ||
| 29 | //----------------------------------------------------------------------------- | ||
| 30 | static const double WX = 0.6; // distance weighting | ||
| 31 | static const double WT = 0.6; // time weighting | ||
| 32 | static const double EPSIL = 0.001; // convergence criterion | ||
| 33 | |||
| 34 | //----------------------------------------------------------------------------- | ||
| 35 | // Shared variables | ||
| 36 | //----------------------------------------------------------------------------- | ||
| 37 | static double Beta1; | ||
| 38 | static double C1; | ||
| 39 | static double C2; | ||
| 40 | static double Afull; | ||
| 41 | static double Qfull; | ||
| 42 | static TXsect* pXsect; | ||
| 43 | |||
| 44 | //----------------------------------------------------------------------------- | ||
| 45 | // External functions (declared in funcs.h) | ||
| 46 | //----------------------------------------------------------------------------- | ||
| 47 | // kinwave_execute (called by flowrout_execute) | ||
| 48 | |||
| 49 | //----------------------------------------------------------------------------- | ||
| 50 | // Local functions | ||
| 51 | //----------------------------------------------------------------------------- | ||
| 52 | static int solveContinuity(double qin, double ain, double* aout); | ||
| 53 | static void evalContinuity(double a, double* f, double* df, void* p); | ||
| 54 | |||
| 55 | //============================================================================= | ||
| 56 | |||
| 57 | 100920 | int kinwave_execute(int j, double* qinflow, double* qoutflow, double tStep) | |
| 58 | // | ||
| 59 | // Input: j = link index | ||
| 60 | // qinflow = inflow at current time (cfs) | ||
| 61 | // tStep = time step (sec) | ||
| 62 | // Output: qoutflow = outflow at current time (cfs), | ||
| 63 | // returns number of iterations used | ||
| 64 | // Purpose: finds outflow over time step tStep given flow entering a | ||
| 65 | // conduit using Kinematic Wave flow routing. | ||
| 66 | // | ||
| 67 | // | ||
| 68 | // ^ q3 | ||
| 69 | // t | | ||
| 70 | // | qin, ain |-------------------| qout, aout | ||
| 71 | // | | Flow ---> | | ||
| 72 | // |----> x q1, a1 |-------------------| q2, a2 | ||
| 73 | // | ||
| 74 | // | ||
| 75 | { | ||
| 76 | int k; | ||
| 77 | 100920 | int result = 1; | |
| 78 | double dxdt, dq; | ||
| 79 | double ain, aout; | ||
| 80 | double qin, qout; | ||
| 81 | double a1, a2, q1, q2, q3; | ||
| 82 | |||
| 83 | // --- no routing for non-conduit link | ||
| 84 | 100920 | (*qoutflow) = (*qinflow); | |
| 85 |
2/2✓ Branch 0 taken 5400 times.
✓ Branch 1 taken 95520 times.
|
100920 | if ( Link[j].type != CONDUIT ) return result; |
| 86 | |||
| 87 | // --- no routing for dummy xsection | ||
| 88 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 95520 times.
|
95520 | if ( Link[j].xsect.type == DUMMY ) return result; |
| 89 | |||
| 90 | // --- assign module-level variables | ||
| 91 | 95520 | pXsect = &Link[j].xsect; | |
| 92 | 95520 | Qfull = Link[j].qFull; | |
| 93 | 95520 | Afull = Link[j].xsect.aFull; | |
| 94 | 95520 | k = Link[j].subIndex; | |
| 95 | 95520 | Beta1 = Conduit[k].beta / Qfull; | |
| 96 | |||
| 97 | // --- normalize previous flows | ||
| 98 | 95520 | q1 = Conduit[k].q1 / Qfull; | |
| 99 | 95520 | q2 = Conduit[k].q2 / Qfull; | |
| 100 | |||
| 101 | // --- normalize inflow | ||
| 102 | 95520 | qin = (*qinflow) / Conduit[k].barrels / Qfull; | |
| 103 | |||
| 104 | // --- compute evaporation and infiltration loss rate | ||
| 105 | 95520 | q3 = link_getLossRate(j, KW, qin*Qfull, tStep) / Qfull; | |
| 106 | |||
| 107 | // --- normalize previous areas | ||
| 108 | 95520 | a1 = Conduit[k].a1 / Afull; | |
| 109 | 95520 | a2 = Conduit[k].a2 / Afull; | |
| 110 | |||
| 111 | // --- use full area when inlet flow >= full flow | ||
| 112 |
2/2✓ Branch 0 taken 2317 times.
✓ Branch 1 taken 93203 times.
|
95520 | if ( qin >= 1.0 ) ain = 1.0; |
| 113 | |||
| 114 | // --- get normalized inlet area corresponding to inlet flow | ||
| 115 | 93203 | else ain = xsect_getAofS(pXsect, qin/Beta1) / Afull; | |
| 116 | |||
| 117 | // --- check for no flow | ||
| 118 |
4/4✓ Branch 0 taken 44944 times.
✓ Branch 1 taken 50576 times.
✓ Branch 2 taken 43314 times.
✓ Branch 3 taken 1630 times.
|
95520 | if ( qin <= TINY && q2 <= TINY ) |
| 119 | { | ||
| 120 | 43314 | qout = 0.0; | |
| 121 | 43314 | aout = 0.0; | |
| 122 | } | ||
| 123 | |||
| 124 | // --- otherwise solve finite difference form of continuity eqn. | ||
| 125 | else | ||
| 126 | { | ||
| 127 | // --- compute constant factors | ||
| 128 | 52206 | dxdt = link_getLength(j) / tStep * Afull / Qfull; | |
| 129 | 52206 | dq = q2 - q1; | |
| 130 | 52206 | C1 = dxdt * WT / WX; | |
| 131 | 52206 | C2 = (1.0 - WT) * (ain - a1); | |
| 132 | 52206 | C2 = C2 - WT * a2; | |
| 133 | 52206 | C2 = C2 * dxdt / WX; | |
| 134 | 52206 | C2 = C2 + (1.0 - WX) / WX * dq - qin; | |
| 135 | 52206 | C2 = C2 + q3 / WX; | |
| 136 | |||
| 137 | // --- starting guess for aout is value from previous time step | ||
| 138 | 52206 | aout = a2; | |
| 139 | |||
| 140 | // --- solve continuity equation for aout | ||
| 141 | 52206 | result = solveContinuity(qin, ain, &aout); | |
| 142 | |||
| 143 | // --- report error if continuity eqn. not solved | ||
| 144 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 52206 times.
|
52206 | if ( result == -1 ) |
| 145 | { | ||
| 146 | ✗ | report_writeErrorMsg(ERR_KINWAVE, Link[j].ID); | |
| 147 | ✗ | return 1; | |
| 148 | } | ||
| 149 |
2/2✓ Branch 0 taken 2376 times.
✓ Branch 1 taken 49830 times.
|
52206 | if ( result <= 0 ) result = 1; |
| 150 | |||
| 151 | // --- compute normalized outlet flow from outlet area | ||
| 152 | 52206 | qout = Beta1 * xsect_getSofA(pXsect, aout*Afull); | |
| 153 |
2/2✓ Branch 0 taken 2317 times.
✓ Branch 1 taken 49889 times.
|
52206 | if ( qin > 1.0 ) qin = 1.0; |
| 154 | } | ||
| 155 | |||
| 156 | // --- save new flows and areas | ||
| 157 | 95520 | Conduit[k].q1 = qin * Qfull; | |
| 158 | 95520 | Conduit[k].a1 = ain * Afull; | |
| 159 | 95520 | Conduit[k].q2 = qout * Qfull; | |
| 160 | 95520 | Conduit[k].a2 = aout * Afull; | |
| 161 | 191040 | Conduit[k].fullState = | |
| 162 | 95520 | link_getFullState(Conduit[k].a1, Conduit[k].a2, Afull); | |
| 163 | 95520 | (*qinflow) = Conduit[k].q1 * Conduit[k].barrels; | |
| 164 | 95520 | (*qoutflow) = Conduit[k].q2 * Conduit[k].barrels; | |
| 165 | 95520 | return result; | |
| 166 | } | ||
| 167 | |||
| 168 | //============================================================================= | ||
| 169 | |||
| 170 | 52206 | int solveContinuity(double qin, double ain, double* aout) | |
| 171 | // | ||
| 172 | // Input: qin = upstream normalized flow | ||
| 173 | // ain = upstream normalized area | ||
| 174 | // aout = downstream normalized area | ||
| 175 | // Output: new value for aout; returns an error code | ||
| 176 | // Purpose: solves continuity equation f(a) = Beta1*S(a) + C1*a + C2 = 0 | ||
| 177 | // for 'a' using the Newton-Raphson root finder function. | ||
| 178 | // Return code has the following meanings: | ||
| 179 | // >= 0 number of function evaluations used | ||
| 180 | // -1 Newton function failed | ||
| 181 | // -2 flow always above max. flow | ||
| 182 | // -3 flow always below zero | ||
| 183 | // | ||
| 184 | // Note: pXsect (pointer to conduit's cross-section), and constants Beta1, | ||
| 185 | // C1, and C2 are module-level shared variables assigned values | ||
| 186 | // in kinwave_execute(). | ||
| 187 | // | ||
| 188 | { | ||
| 189 | int n; // # evaluations or error code | ||
| 190 | double aLo, aHi, aTmp; // lower/upper bounds on a | ||
| 191 | double fLo, fHi; // lower/upper bounds on f | ||
| 192 | 52206 | double tol = EPSIL; // absolute convergence tol. | |
| 193 | |||
| 194 | // --- first determine bounds on 'a' so that f(a) passes through 0. | ||
| 195 | |||
| 196 | // --- set upper bound to area at full flow | ||
| 197 | 52206 | aHi = 1.0; | |
| 198 | 52206 | fHi = 1.0 + C1 + C2; | |
| 199 | |||
| 200 | // --- try setting lower bound to area where section factor is maximum | ||
| 201 | 52206 | aLo = xsect_getAmax(pXsect) / Afull; | |
| 202 |
1/2✓ Branch 0 taken 52206 times.
✗ Branch 1 not taken.
|
52206 | if ( aLo < aHi ) |
| 203 | { | ||
| 204 | 52206 | fLo = ( Beta1 * pXsect->sMax ) + (C1 * aLo) + C2; | |
| 205 | } | ||
| 206 | ✗ | else fLo = fHi; | |
| 207 | |||
| 208 | // --- if fLo and fHi have same sign then set lower bound to 0 | ||
| 209 |
2/2✓ Branch 0 taken 52204 times.
✓ Branch 1 taken 2 times.
|
52206 | if ( fHi*fLo > 0.0 ) |
| 210 | { | ||
| 211 | 52204 | aHi = aLo; | |
| 212 | 52204 | fHi = fLo; | |
| 213 | 52204 | aLo = 0.0; | |
| 214 | 52204 | fLo = C2; | |
| 215 | } | ||
| 216 | |||
| 217 | // --- proceed with search for root if fLo and fHi have different signs | ||
| 218 |
2/2✓ Branch 0 taken 49830 times.
✓ Branch 1 taken 2376 times.
|
52206 | if ( fHi*fLo <= 0.0 ) |
| 219 | { | ||
| 220 | // --- start search at midpoint of lower/upper bounds | ||
| 221 | // if initial value outside of these bounds | ||
| 222 |
3/4✓ Branch 0 taken 49830 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1 time.
✓ Branch 3 taken 49829 times.
|
49830 | if ( *aout < aLo || *aout > aHi ) *aout = 0.5*(aLo + aHi); |
| 223 | |||
| 224 | // --- if fLo > fHi then switch aLo and aHi | ||
| 225 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 49830 times.
|
49830 | if ( fLo > fHi ) |
| 226 | { | ||
| 227 | ✗ | aTmp = aLo; | |
| 228 | ✗ | aLo = aHi; | |
| 229 | ✗ | aHi = aTmp; | |
| 230 | } | ||
| 231 | |||
| 232 | // --- call the Newton root finder method passing it the | ||
| 233 | // evalContinuity function to evaluate the function | ||
| 234 | // and its derivatives | ||
| 235 | 49830 | n = findroot_Newton(aLo, aHi, aout, tol, evalContinuity, NULL); | |
| 236 | |||
| 237 | // --- check if root finder succeeded | ||
| 238 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 49830 times.
|
49830 | if ( n <= 0 ) n = -1; |
| 239 | } | ||
| 240 | |||
| 241 | // --- if lower/upper bound functions both negative then use full flow | ||
| 242 |
2/2✓ Branch 0 taken 2308 times.
✓ Branch 1 taken 68 times.
|
2376 | else if ( fLo < 0.0 ) |
| 243 | { | ||
| 244 |
2/2✓ Branch 0 taken 2298 times.
✓ Branch 1 taken 10 times.
|
2308 | if ( qin > 1.0 ) *aout = ain; |
| 245 | 10 | else *aout = 1.0; | |
| 246 | 2308 | n = -2; | |
| 247 | } | ||
| 248 | |||
| 249 | // --- if lower/upper bound functions both positive then use no flow | ||
| 250 |
1/2✓ Branch 0 taken 68 times.
✗ Branch 1 not taken.
|
68 | else if ( fLo > 0 ) |
| 251 | { | ||
| 252 | 68 | *aout = 0.0; | |
| 253 | 68 | n = -3; | |
| 254 | } | ||
| 255 | ✗ | else n = -1; | |
| 256 | 52206 | return n; | |
| 257 | } | ||
| 258 | |||
| 259 | //============================================================================= | ||
| 260 | |||
| 261 | 75865 | void evalContinuity(double a, double* f, double* df, void* p) | |
| 262 | // | ||
| 263 | // Input: a = outlet normalized area | ||
| 264 | // Output: f = value of continuity eqn. | ||
| 265 | // df = derivative of continuity eqn. | ||
| 266 | // Purpose: computes value of continuity equation (f) and its derivative (df) | ||
| 267 | // w.r.t. normalized area for link with normalized outlet area 'a'. | ||
| 268 | // | ||
| 269 | { | ||
| 270 | 75865 | *f = (Beta1 * xsect_getSofA(pXsect, a*Afull)) + (C1 * a) + C2; | |
| 271 | 75865 | *df = (Beta1 * Afull * xsect_getdSdA(pXsect, a*Afull)) + C1; | |
| 272 | 75865 | } | |
| 273 | |||
| 274 | //============================================================================= | ||
| 275 |