diff --git a/src/dwflow.c b/src/dwflow.c index a9382a0db..4ef2464f6 100644 --- a/src/dwflow.c +++ b/src/dwflow.c @@ -5,6 +5,7 @@ // Version: 5.1 // Date: 03/20/14 (Build 5.1.001) // 03/19/15 (Build 5.1.008) +// 03/14/17 (Build 5.1.012) // Author: L. Rossman (EPA) // M. Tryby (EPA) // R. Dickinson (CDM) @@ -15,6 +16,8 @@ // Build 5.1.008: // - Bug in finding if conduit was upstrm/dnstrm full was fixed. // +// Build 5.1.012: +// - Modified uniform loss rate term of conduit momentum equation. //----------------------------------------------------------------------------- #define _CRT_SECURE_NO_DEPRECATE @@ -76,6 +79,8 @@ void dwflow_findConduitFlow(int j, int steps, double omega, double dt) char isFull = FALSE; // TRUE if conduit flowing full char isClosed = FALSE; // TRUE if conduit closed + + // --- adjust isClosed status by any control action if ( Link[j].setting == 0 ) isClosed = TRUE; @@ -208,7 +213,7 @@ void dwflow_findConduitFlow(int j, int steps, double omega, double dt) } // --- 6. term for evap and seepage losses per unit length - dq6 = link_getLossRate(j, qOld, dt) * 1.5 * dt * v / link_getLength(j); //(5.1.008) + dq6 = link_getLossRate(j, qOld, dt) * 2.5 * dt * v / link_getLength(j); //(5.1.012) // --- combine terms to find new conduit flow denom = 1.0 + dq1 + dq5; diff --git a/src/flowrout.c b/src/flowrout.c index d8945e3ae..f1a3b1b10 100644 --- a/src/flowrout.c +++ b/src/flowrout.c @@ -6,6 +6,7 @@ // Date: 03/19/14 (Build 5.1.001) // 09/15/14 (Build 5.1.007) // 03/19/15 (Build 5.1.008) +// 03/14/11 (Build 5.1.012) // Author: L. Rossman (EPA) // M. Tryby (EPA) // @@ -20,6 +21,10 @@ // - Determination of node crown elevations moved to dynwave.c. // - Support added for new way of recording conduit's fullness state. // +// Build 5.1.012: +// - Overflow computed in updateStorageState() must be non-negative. +// - Terminal storage nodes now updated corectly. +// //----------------------------------------------------------------------------- #define _CRT_SECURE_NO_DEPRECATE @@ -570,6 +575,7 @@ void updateStorageState(int i, int j, int links[], double dt) { Node[i].overflow = (v2 - MAX(Node[i].oldVolume, Node[i].fullVolume)) / dt; + if ( Node[i].overflow < FUDGE ) Node[i].overflow = 0.0; //(5.1.012) if ( !AllowPonding || Node[i].pondedArea == 0.0 ) v2 = Node[i].fullVolume; } @@ -632,8 +638,15 @@ void setNewNodeState(int j, double dt) int canPond; // TRUE if ponding can occur at node double newNetInflow; // inflow - outflow at node (cfs) - // --- storage nodes have already been updated - if ( Node[j].type == STORAGE ) return; +//// Following section revised for release 5.1.012. //// //(5.1.012) + // --- update terminal storage nodes + if ( Node[j].type == STORAGE ) + { + if ( Node[j].updated == FALSE ) + updateStorageState(j, Nobjects[LINK], NULL, dt); + return; + } +////////////////////////////////////////////////////////// // --- update stored volume using mid-point integration newNetInflow = Node[j].inflow - Node[j].outflow - Node[j].losses; //(5.1.007) @@ -778,6 +791,11 @@ int steadyflow_execute(int j, double* qin, double* qout, double tStep) } } Conduit[k].a2 = Conduit[k].a1; + + Conduit[k].q1Old = Conduit[k].q1; + Conduit[k].q2Old = Conduit[k].q2; + + Conduit[k].q1 = q; Conduit[k].q2 = q; (*qout) = q * Conduit[k].barrels; diff --git a/src/globals.h b/src/globals.h index a6cc7d6d0..1e7fdbeff 100644 --- a/src/globals.h +++ b/src/globals.h @@ -8,6 +8,7 @@ // 09/15/14 (Build 5.1.007) // 03/19/15 (Build 5.1.008) // 08/01/16 (Build 5.1.011) +// 03/14/17 (Build 5.1.012) // Author: L. Rossman // // Global Variables @@ -27,6 +28,9 @@ // - Added error message text as a variable. // - Added elapsed simulation time (in decimal days) variable. // - Added variables associated with detailed routing events. +// +// Build 5.1.012: +// - InSteadyState variable made local to routing_execute in routing.c. //----------------------------------------------------------------------------- EXTERN TFile @@ -87,8 +91,8 @@ EXTERN int SweepEnd, // Day of year when sweeping ends MaxTrials, // Max. trials for DW routing NumThreads, // Number of parallel threads used //(5.1.008) - NumEvents, // Number of detailed events //(5.1.011) - InSteadyState; // System flows remain constant //(5.1.011) + NumEvents; // Number of detailed events //(5.1.011) + //InSteadyState; // System flows remain constant //(5.1.012) EXTERN double RouteStep, // Routing time step (sec) diff --git a/src/lid.c b/src/lid.c index be1ecddac..a689ad177 100644 --- a/src/lid.c +++ b/src/lid.c @@ -10,6 +10,7 @@ // 04/30/15 (Build 5.1.009) // 08/05/15 (Build 5.1.010) // 08/01/16 (Build 5.1.011) +// 03/14/17 (Build 5.1.012) // Author: L. Rossman (US EPA) // // This module handles all data processing involving LID (Low Impact @@ -60,6 +61,9 @@ // - The top of the storage layer is no longer used as a limit for an // underdrain offset thus allowing upturned drains to be modeled. // - Column headings for the detailed LID report file were modified. +// +// Build 5.1.012: +// - Redefined initialization of wasDry for LID reporting. //----------------------------------------------------------------------------- #define _CRT_SECURE_NO_DEPRECATE @@ -1898,5 +1902,6 @@ void initLidRptFile(char* title, char* lidID, char* subcatchID, TLidUnit* lidUni for ( i = 1; i < colCount; i++) fprintf(f, "\t%s", line9); //... initialize LID dryness state - lidUnit->rptFile->wasDry = TRUE; + lidUnit->rptFile->wasDry = 1; //(5.1.012) + strcpy(lidUnit->rptFile->results, ""); //(5.1.012) } diff --git a/src/lid.h b/src/lid.h index 4b39da3e6..f41fd3976 100644 --- a/src/lid.h +++ b/src/lid.h @@ -6,6 +6,7 @@ // Date: 03/20/14 (Build 5.1.001) // 03/19/15 (Build 5.1.008) // 08/01/16 (Build 5.1.011) +// 03/14/17 (Build 5.1.012) // Author: L. Rossman (US EPA) // // Public interface for LID functions. @@ -19,6 +20,9 @@ // - Water depth replaces moisture content for LID's pavement layer. // - Arguments for lidproc_saveResults() modified. // +// Build 5.1.012: +// - Redefined meaning of wasDry in TLidRptFile structure. +// //----------------------------------------------------------------------------- #ifndef LID_H @@ -141,7 +145,7 @@ typedef struct typedef struct { FILE* file; // file pointer - int wasDry; // true if LID was dry //(5.1.008) + int wasDry; // number of successive dry periods //(5.1.012) char results[256]; // results for current time period //(5.1.008) } TLidRptFile; diff --git a/src/lidproc.c b/src/lidproc.c index 33ca88510..693783447 100644 --- a/src/lidproc.c +++ b/src/lidproc.c @@ -10,6 +10,7 @@ // 04/30/15 (Build 5.1.009) // 08/05/15 (Build 5.1.010) // 08/01/16 (Build 5.1.011) +// 03/14/17 (Build 5.1.012) // Author: L. Rossman (US EPA) // // This module computes the hydrologic performance of an LID (Low Impact @@ -46,6 +47,12 @@ // physically meaningful results. // - Reporting of detailed results re-written. // +// Build 5.1.012: +// - Modified upper limit for soil layer percolation. +// - Modified upper limit on surface infiltration into rain gardens. +// - Modified upper limit on drain flow for LIDs with storage layers. +// - Used re-defined wasDry variable for LID reports to fix duplicate lines. +// //----------------------------------------------------------------------------- #define _CRT_SECURE_NO_DEPRECATE @@ -378,14 +385,15 @@ void lidproc_saveResults(TLidUnit* lidUnit, double ucfRainfall, double ucfRainDe if ( SurfaceInflow < MINFLOW && SurfaceOutflow < MINFLOW && StorageDrain < MINFLOW && - StorageExfil < MINFLOW && - totalEvap < MINFLOW ) isDry = TRUE; + StorageExfil < MINFLOW && + totalEvap < MINFLOW + ) isDry = TRUE; //... update status of HasWetLids if ( !isDry ) HasWetLids = TRUE; - //... write results to LID report file if at a reporting time //(5.1.011) - if ( lidUnit->rptFile && fabs(NewRunoffTime - ReportTime) < 10. ) //(5.1.011) + //... write results to LID report file //(5.1.012) + if ( lidUnit->rptFile ) //(5.1.012) { //... convert rate results to original units (in/hr or mm/hr) ucf = ucfRainfall; @@ -406,11 +414,13 @@ void lidproc_saveResults(TLidUnit* lidUnit, double ucfRainfall, double ucfRainDe rptVars[STOR_DEPTH] = theLidUnit->storageDepth*ucf; //... if the current LID state is wet but the previous state was dry - // then write the saved previous results to the report file thus - // marking the end of a dry period - if ( !isDry && theLidUnit->rptFile->wasDry ) + // for more than one period then write the saved previous results //(5.1.012) + // to the report file thus marking the end of a dry period //(5.10012) + if ( !isDry && theLidUnit->rptFile->wasDry > 1) //(5.1.012) + { fprintf(theLidUnit->rptFile->file, "%s", - theLidUnit->rptFile->results); + theLidUnit->rptFile->results); + } //... write the current results to a string which is saved between // reporting periods @@ -428,20 +438,25 @@ void lidproc_saveResults(TLidUnit* lidUnit, double ucfRainfall, double ucfRainDe { //... if the previous state was wet then write the current // results to file marking the start of a dry period - if ( !theLidUnit->rptFile->wasDry ) + if ( theLidUnit->rptFile->wasDry == 0 ) //(5.1.012) { - fprintf(theLidUnit->rptFile->file, "%s", theLidUnit->rptFile->results); + fprintf(theLidUnit->rptFile->file, "%s", + theLidUnit->rptFile->results); } - theLidUnit->rptFile->wasDry = TRUE; + + //... increment the number of successive dry periods //(5.1.012) + theLidUnit->rptFile->wasDry++; //(5.1.012) } //... if the current LID state is wet else - { //... write the current results to the report file - fprintf(theLidUnit->rptFile->file, "%s", theLidUnit->rptFile->results); - theLidUnit->rptFile->wasDry = FALSE; + fprintf(theLidUnit->rptFile->file, "%s", + theLidUnit->rptFile->results); + + //... re-set the number of successive dry periods to 0 //(5.1.012) + theLidUnit->rptFile->wasDry = 0; //(5.1.012) } } } @@ -518,7 +533,7 @@ void greenRoofFluxRates(double x[], double f[]) //... limit perc rate by available water availVolume = (soilTheta - soilFieldCap) * soilThickness; - maxRate = SurfaceInfil - SoilEvap + MAX(availVolume, 0.0) / Tstep; + maxRate = MAX(availVolume, 0.0) / Tstep - SoilEvap; //(5.1.012) SoilPerc = MIN(SoilPerc, maxRate); SoilPerc = MAX(SoilPerc, 0.0); @@ -542,8 +557,8 @@ void greenRoofFluxRates(double x[], double f[]) else { //... limit drainmat outflow by available storage volume - maxRate = storageDepth * storageVoidFrac / Tstep + - SoilPerc - StorageEvap; + maxRate = storageDepth * storageVoidFrac / Tstep - StorageEvap; //(5.1.012) + if ( storageDepth >= storageThickness ) maxRate += SoilPerc; //(5.1.012) maxRate = MAX(maxRate, 0.0); StorageDrain = MIN(StorageDrain, maxRate); @@ -618,7 +633,7 @@ void biocellFluxRates(double x[], double f[]) //... limit perc rate by available water availVolume = (soilTheta - soilFieldCap) * soilThickness; - maxRate = SurfaceInfil - SoilEvap + MAX(availVolume, 0.0) / Tstep; + maxRate = MAX(availVolume, 0.0) / Tstep - SoilEvap; //(5.1.012) SoilPerc = MIN(SoilPerc, maxRate); SoilPerc = MAX(SoilPerc, 0.0); @@ -640,7 +655,15 @@ void biocellFluxRates(double x[], double f[]) maxRate = MIN(SoilPerc, StorageExfil); SoilPerc = maxRate; StorageExfil = maxRate; - } + +//// Following code segment added to release 5.1.012 //// //(5.1.012) + //... limit surface infil. by unused soil volume + maxRate = (soilPorosity - soilTheta) * soilThickness / Tstep + + SoilPerc + SoilEvap; + SurfaceInfil = MIN(SurfaceInfil, maxRate); +////////////////////////////////////////////////////////// + + } //... storage & soil layers are full else if ( soilTheta >= soilPorosity && storageDepth >= storageThickness ) @@ -668,13 +691,14 @@ void biocellFluxRates(double x[], double f[]) { //... limit storage exfiltration by available storage volume maxRate = SoilPerc - StorageEvap + storageDepth*storageVoidFrac/Tstep; - maxRate = MAX(maxRate, 0.0); StorageExfil = MIN(StorageExfil, maxRate); + StorageExfil = MAX(StorageExfil, 0.0); //... limit underdrain flow by volume above drain offset if ( StorageDrain > 0.0 ) { - maxRate = SoilPerc - StorageExfil - StorageEvap; + maxRate = -StorageExfil - StorageEvap; //(5.1.012) + if ( storageDepth >= storageThickness) maxRate += SoilPerc; //(5.1.012) if ( theLidProc->drain.offset <= storageDepth ) { maxRate += (storageDepth - theLidProc->drain.offset) * @@ -766,14 +790,15 @@ void trenchFluxRates(double x[], double f[]) StorageExfil = MIN(StorageExfil, maxRate); StorageExfil = MAX(StorageExfil, 0.0); - //... limit underdrain flow by available storage volume + //... limit underdrain flow by volume above drain offset if ( StorageDrain > 0.0 ) { - maxRate = StorageInflow - StorageExfil - StorageEvap; + maxRate = -StorageExfil - StorageEvap; //(5.1.012) + if (storageDepth >= storageThickness ) maxRate += StorageInflow; //(5.1.012) if ( theLidProc->drain.offset <= storageDepth ) { maxRate += (storageDepth - theLidProc->drain.offset) * - storageVoidFrac/Tstep; + storageVoidFrac/Tstep; } maxRate = MAX(maxRate, 0.0); StorageDrain = MIN(StorageDrain, maxRate); @@ -869,7 +894,7 @@ void pavementFluxRates(double x[], double f[]) { SoilPerc = getSoilPercRate(soilTheta); availVolume = (soilTheta - soilFieldCap) * soilThickness; - maxRate = PavePerc - SoilEvap + MAX(availVolume, 0.0) / Tstep; + maxRate = MAX(availVolume, 0.0) / Tstep - SoilEvap; //(5.1.012) SoilPerc = MIN(SoilPerc, maxRate); SoilPerc = MAX(SoilPerc, 0.0); } @@ -982,8 +1007,9 @@ void pavementFluxRates(double x[], double f[]) //... limit underdrain flow by volume above drain offset if ( StorageDrain > 0.0 ) { - maxRate = SoilPerc - StorageExfil - StorageEvap; - if ( theLidProc->drain.offset <= storageDepth ) + maxRate = -StorageExfil - StorageEvap; //(5.1.012) + if (storageDepth >= storageThickness ) maxRate += SoilPerc; //(5.1.012) + if ( theLidProc->drain.offset <= storageDepth ) { maxRate += (storageDepth - theLidProc->drain.offset) * storageVoidFrac/Tstep; diff --git a/src/link.c b/src/link.c index 9a058fcf7..f2e207ae4 100644 --- a/src/link.c +++ b/src/link.c @@ -8,6 +8,7 @@ // 03/19/15 (Build 5.1.008) // 08/05/15 (Build 5.1.010) // 08/01/16 (Build 5.1.011) +// 03/14/17 (Build 5.1.012) // Author: L. Rossman (EPA) // M. Tryby (EPA) // @@ -32,6 +33,11 @@ // - Weir shape parameter deprecated. // - Extra geometric parameters ignored for non-conduit open rectangular // cross sections. +// +// Build 5.1.012: +// - Conduit seepage rate now based on flow width, not wetted perimeter. +// - Formula for side flow weir corrected. +// - Crest length contraction adjustments corrected. //----------------------------------------------------------------------------- #define _CRT_SECURE_NO_DEPRECATE @@ -1296,7 +1302,7 @@ double conduit_getLossRate(int j, double q, double tStep) / double depth = 0.5 * (Link[j].oldDepth + Link[j].newDepth); double length; double topWidth; - double wettedPerimeter; + //double wettedPerimeter; //DEPRECATED //(5.1.012) double maxLossRate; double evapLossRate = 0.0, seepLossRate = 0.0, @@ -1319,17 +1325,20 @@ double conduit_getLossRate(int j, double q, double tStep) / { // limit depth to depth at max width if ( depth >= xsect->ywMax ) depth = xsect->ywMax; - + +//// The following section was deprecated for release 5.1.012. //// //(5.1.012) // get wetted perimeter - wettedPerimeter = 0.0; - if ( depth > 0.0 ) - { - wettedPerimeter = xsect_getAofY(xsect, depth) / - xsect_getRofY(xsect, depth); - } +// wettedPerimeter = 0.0; +// if ( depth > 0.0 ) +// { +// wettedPerimeter = xsect_getAofY(xsect, depth) / +// xsect_getRofY(xsect, depth); +// } +///////////////////////////////////////////////////////////////// // compute seepage loss rate across length of conduit - seepLossRate = Link[j].seepRate * wettedPerimeter * length; + seepLossRate = Link[j].seepRate * xsect_getWofY(xsect, depth) * //(5.1.012) + length; seepLossRate *= Adjust.hydconFactor; //(5.1.008) } @@ -2312,9 +2321,11 @@ void weir_getFlow(int j, int k, double head, double dir, int hasFlapGate, length = Link[j].xsect.wMax * UCF(LENGTH); h = head * UCF(LENGTH); +//// Following code segment re-located. //// //(5.1.012) // --- reduce length when end contractions present - length -= 0.1 * Weir[k].endCon * h; - length = MAX(length, 0.0); + //length -= 0.1 * Weir[k].endCon * h; + //length = MAX(length, 0.0); +///////////////////////////////////////////// // --- use appropriate formula for weir flow wType = Weir[k].type; @@ -2323,15 +2334,28 @@ void weir_getFlow(int j, int k, double head, double dir, int hasFlapGate, switch (wType) { case TRANSVERSE_WEIR: + + // --- reduce length when end contractions present //(5.1.012) + length -= 0.1 * Weir[k].endCon * h; //(5.1.012) + length = MAX(length, 0.0); //(5.1.012) *q1 = Weir[k].cDisch1 * length * pow(h, 1.5); break; case SIDEFLOW_WEIR: + + // --- reduce length when end contractions present //(5.1.012) + length -= 0.1 * Weir[k].endCon * h; //(5.1.012) + length = MAX(length, 0.0); //(5.1.012) + // --- weir behaves as a transverse weir under reverse flow if ( dir < 0.0 ) *q1 = Weir[k].cDisch1 * length * pow(h, 1.5); else - *q1 = Weir[k].cDisch1 * length * pow(h, 5./3.); + +//// Corrected formula //// //(5.1.012) +// (see Metcalf & Eddy, Inc., Wastewater Engineering, McGraw-Hill, 1972 p. 164). + *q1 = Weir[k].cDisch1 * pow(length, 0.83) * pow(h, 1.67); + break; case VNOTCH_WEIR: @@ -2341,8 +2365,11 @@ void weir_getFlow(int j, int k, double head, double dir, int hasFlapGate, case TRAPEZOIDAL_WEIR: y = (1.0 - Link[j].setting) * Link[j].xsect.yFull; length = xsect_getWofY(&Link[j].xsect, y) * UCF(LENGTH); - length -= 0.1 * Weir[k].endCon * h; - length = MAX(length, 0.0); + +//// End contractions don't apply to trapezoidal weirs //// //(5.1.012) + //length -= 0.1 * Weir[k].endCon * h; //(5.1.012) + //length = MAX(length, 0.0); //(5.1.012) + *q1 = Weir[k].cDisch1 * length * pow(h, 1.5); *q2 = Weir[k].cDisch2 * Weir[k].slope * pow(h, 2.5); } @@ -2455,7 +2482,7 @@ double weir_getdqdh(int k, double dir, double h, double q1, double q2) case SIDEFLOW_WEIR: // --- weir behaves as a transverse weir under reverse flow if ( dir < 0.0 ) return 1.5 * q1h; - else return 5./3. * q1h; + else return 1.67 * q1h; //(5.1.012) case VNOTCH_WEIR: if ( q2h == 0.0 ) return 2.5 * q1h; // Fully open diff --git a/src/massbal.c b/src/massbal.c index e4a4c8aa8..a5e8c3ccc 100644 --- a/src/massbal.c +++ b/src/massbal.c @@ -8,6 +8,7 @@ // 04/02/15 (Build 5.1.008) // 08/05/15 (Build 5.1.010) // 08/01/16 (Build 5.1.011) +// 03/14/17 (Build 5.1.012) // Author: L. Rossman (EPA) // M. Tryby (EPA) // @@ -26,8 +27,12 @@ // Build 5.1.010: // - Remaining pollutant mass in "dry" elements now added to final storage. // -// Build 5.1.011 +// Build 5.1.011: // - Final stored pollutant mass in links ignored for Steady Flow routing. +// +// Build 5.1.012: +// - Terminal storage nodes no longer treated as non-storage terminal +// nodes are when updating total outflow volume. //----------------------------------------------------------------------------- #define _CRT_SECURE_NO_DEPRECATE @@ -640,7 +645,8 @@ void massbal_updateRoutingTotals(double tStep) for ( j = 0; j < Nobjects[NODE]; j++) { NodeInflow[j] += Node[j].inflow * tStep; - if ( Node[j].type == OUTFALL || Node[j].degree == 0 ) + if ( Node[j].type == OUTFALL || + (Node[j].degree == 0 && Node[j].type != STORAGE) ) //(5.1.012) { NodeOutflow[j] += Node[j].inflow * tStep; } diff --git a/src/project.c b/src/project.c index e061ae759..de6f333fd 100644 --- a/src/project.c +++ b/src/project.c @@ -9,6 +9,7 @@ // 03/19/15 (Build 5.1.008) // 04/30/15 (Build 5.1.009) // 08/01/16 (Build 5.1.011) +// 03/14/17 (Build 5.1.012) // Author: L. Rossman // // Project management functions. @@ -42,6 +43,10 @@ // Build 5.1.011: // - Memory management of hydraulic event dates array added. // +// Build 5.1.012: +// - Minimum conduit slope option initialized to 0 (none). +// - NO/YES no longer accepted as options for NORMAL_FLOW_LIMITED. +// //----------------------------------------------------------------------------- #define _CRT_SECURE_NO_DEPRECATE @@ -571,7 +576,7 @@ int project_readOption(char* s1, char* s2) case NORMAL_FLOW_LTD: m = findmatch(s2, NormalFlowWords); - if ( m < 0 ) m = findmatch(s2, NoYesWords); + //if ( m < 0 ) m = findmatch(s2, NoYesWords); DEPRECATED //(5.1.012) if ( m < 0 ) return error_setInpError(ERR_KEYWORD, s2); NormalFlowLtd = m; break; @@ -789,6 +794,7 @@ void setDefaults() LengtheningStep = 0; // No lengthening of conduits CourantFactor = 0.0; // No variable time step MinSurfArea = 0.0; // Force use of default min. surface area + MinSlope = 0.0; // No user supplied minimum conduit slope //(5.1.012) SkipSteadyState = FALSE; // Do flow routing in steady state periods IgnoreRainfall = FALSE; // Analyze rainfall/runoff IgnoreRDII = FALSE; // Analyze RDII //(5.1.004) diff --git a/src/report.c b/src/report.c index af05f44e8..f6183ea2d 100644 --- a/src/report.c +++ b/src/report.c @@ -8,6 +8,7 @@ // 09/15/14 (Build 5.1.007) // 04/02/15 (Build 5.1.008) // 08/01/16 (Build 5.1.011) +// 03/14/17 (Build 5.1.012) // Author: L. Rossman (EPA) // // Report writing functions. @@ -30,6 +31,9 @@ // - Text of error message saved to global variable ErrorMsg. // - Global variable Warnings incremented after warning message issued. // +// Build 5.1.012: +// - System time step statistics adjusted for time in steady state. +// //----------------------------------------------------------------------------- #define _CRT_SECURE_NO_DEPRECATE @@ -1008,8 +1012,10 @@ void report_writeSysStats(TSysStats* sysStats) // { double x; + double eventStepCount = (double)StepCount - sysStats->steadyStateCount; //(5.1.012) - if ( Nobjects[LINK] == 0 || StepCount == 0 ) return; + if ( Nobjects[LINK] == 0 || StepCount == 0 + || eventStepCount == 0.0 ) return; //(5.1.012) WRITE(""); WRITE("*************************"); WRITE("Routing Time Step Summary"); @@ -1019,19 +1025,19 @@ void report_writeSysStats(TSysStats* sysStats) sysStats->minTimeStep); fprintf(Frpt.file, "\n Average Time Step : %7.2f sec", - sysStats->avgTimeStep / StepCount); + sysStats->avgTimeStep / eventStepCount); //(5.1.012) fprintf(Frpt.file, "\n Maximum Time Step : %7.2f sec", sysStats->maxTimeStep); - x = sysStats->steadyStateCount / StepCount * 100.0; + x = (1.0 - sysStats->avgTimeStep * 1000.0 / NewRoutingTime) * 100.0; //(5.1.012) fprintf(Frpt.file, "\n Percent in Steady State : %7.2f", MIN(x, 100.0)); fprintf(Frpt.file, "\n Average Iterations per Step : %7.2f", - sysStats->avgStepCount / StepCount); + sysStats->avgStepCount / eventStepCount); //(5.1.012) fprintf(Frpt.file, "\n Percent Not Converging : %7.2f", - 100.0 * (double)NonConvergeCount / StepCount); + 100.0 * (double)NonConvergeCount / eventStepCount); //(5.1.012) WRITE(""); } diff --git a/src/roadway.c b/src/roadway.c index 795f9934d..4a1248e01 100644 --- a/src/roadway.c +++ b/src/roadway.c @@ -4,6 +4,7 @@ // Project: EPA SWMM5 // Version: 5.1 // Date: 08/05/15 (Build 5.1.010) +// 03/14/17 (Build 5.1.012) // Author: L. Rossman // // Roadway Weir module for SWMM5 @@ -15,6 +16,8 @@ // weir has the same upstream node but with an offset equal to the height // of the roadway. // +// Build 5.1.012: +// - Entries in discharge coeff. table for gravel roadways corrected. //----------------------------------------------------------------------------- #define _CRT_SECURE_NO_DEPRECATE @@ -40,7 +43,7 @@ static const double Cr_Low_Paved[4][2] = { static const int N_Cr_Low_Gravel = 8; static const double Cr_Low_Gravel[8][2] = { {0.0, 2.5}, {0.5, 2.7}, {1.0, 2.8}, {1.5, 2.9}, {2.0, 2.98}, - {2.5, 3.2}, {3.0, 3.03}, {4.0, 3.05} }; + {2.5, 3.02}, {3.0, 3.03}, {4.0, 3.05} }; //(5.1.012) // Discharge Coefficients for (head / road width) > 0.15 static const int N_Cr_High_Paved = 2; diff --git a/src/routing.c b/src/routing.c index db65622dd..0c70cb912 100644 --- a/src/routing.c +++ b/src/routing.c @@ -8,6 +8,7 @@ // 04/02/15 (Build 5.1.008) // 08/05/15 (Build 5.1.010) // 08/01/16 (Build 5.1.011) +// 03/14/17 (Build 5.1.012) // Author: L. Rossman (EPA) // M. Tryby (EPA) // @@ -32,6 +33,10 @@ // Build 5.1.011: // - Support added for limiting flow routing to specific events. // +// Build 5.1.012: +// - routing_execute() was re-written so that Routing Events and +// Skip Steady Flow options work together correctly. +// //----------------------------------------------------------------------------- #define _CRT_SECURE_NO_DEPRECATE @@ -47,6 +52,7 @@ //----------------------------------------------------------------------------- static int* SortedLinks; static int NextEvent; //(5.1.011) +static int BetweenEvents; //(5.1.012) //----------------------------------------------------------------------------- // External functions (declared in funcs.h) @@ -108,7 +114,8 @@ int routing_open() // --- initialize routing events //(5.1.011) if ( NumEvents > 0 ) sortEvents(); //(5.1.011) NextEvent = 0; //(5.1.011) - InSteadyState = (NumEvents > 0); //(5.1.011) + //InSteadyState = (NumEvents > 0); //(5.1.012) + BetweenEvents = (NumEvents > 0); //(5.1.012) return ErrorCode; } @@ -147,7 +154,7 @@ double routing_getRoutingStep(int routingModel, double fixedStep) if ( Nobjects[LINK] == 0 ) return fixedStep; // --- find largest step possible if between routing events - if ( NumEvents > 0 && InSteadyState ) + if ( NumEvents > 0 && BetweenEvents ) //(5.1.012) { nextTime = MIN(NewRunoffTime, ReportTime); date1 = getDateTime(NewRoutingTime); @@ -169,6 +176,8 @@ double routing_getRoutingStep(int routingModel, double fixedStep) //============================================================================= +//// This function was re-written for release 5.1.012. //// //(5.1.012) + void routing_execute(int routingModel, double routingStep) // // Input: routingModel = routing method code @@ -180,7 +189,7 @@ void routing_execute(int routingModel, double routingStep) int j; int stepCount = 1; int actionCount = 0; -// int inSteadyState = FALSE; //Deleted //(5.1.011) + int inSteadyState = FALSE; DateTime currentDate; double stepFlowError; @@ -189,7 +198,6 @@ void routing_execute(int routingModel, double routingStep) if ( ErrorCode ) return; massbal_updateRoutingTotals(routingStep/2.); -//// Following code block modified for release 5.1.011. //// //(5.1.011) // --- find new link target settings that are not related to // --- control rules (e.g., pump on/off depth limits) for (j=0; j Event[NextEvent].end ) { - InSteadyState = TRUE; + BetweenEvents = TRUE; NextEvent++; } - else if ( currentDate >= Event[NextEvent].start && InSteadyState == TRUE ) + else if ( currentDate >= Event[NextEvent].start && BetweenEvents == TRUE ) { - InSteadyState = FALSE; + BetweenEvents = FALSE; } } - // --- add lateral inflows and evap/seepage losses at nodes //(5.1.007) - if ( InSteadyState == FALSE ) + // --- if not between routing events + if ( BetweenEvents == FALSE ) { + // --- find evap. & seepage losses from storage nodes for (j = 0; j < Nobjects[NODE]; j++) { Node[j].losses = node_getLosses(j, routingStep); } + + // --- add lateral inflows and evap/seepage losses at nodes addExternalInflows(currentDate); addDryWeatherInflows(currentDate); - addWetWeatherInflows(OldRoutingTime); //(5.1.008) - addGroundwaterInflows(OldRoutingTime); //(5.1.008) - addLidDrainInflows(OldRoutingTime); //(5.1.008) + addWetWeatherInflows(OldRoutingTime); + addGroundwaterInflows(OldRoutingTime); + addLidDrainInflows(OldRoutingTime); addRdiiInflows(currentDate); addIfaceInflows(currentDate); - } -//// - // --- check if can skip steady state periods based on flows //(5.1.011) - if ( SkipSteadyState ) - { - if ( OldRoutingTime == 0.0 - || actionCount > 0 - || fabs(stepFlowError) > SysFlowTol - || inflowHasChanged() ) InSteadyState = FALSE; //(5.1.011) - else InSteadyState = TRUE; //(5.1.011) - } + // --- check if can skip steady state periods based on flows + if ( SkipSteadyState ) + { + if ( OldRoutingTime == 0.0 + || actionCount > 0 + || fabs(stepFlowError) > SysFlowTol + || inflowHasChanged() ) inSteadyState = FALSE; + else inSteadyState = TRUE; + } - // --- find new hydraulic state if system has changed - if ( InSteadyState == FALSE ) //(5.1.011) - { - // --- replace old hydraulic state values with current ones - for (j = 0; j < Nobjects[LINK]; j++) link_setOldHydState(j); - for (j = 0; j < Nobjects[NODE]; j++) + // --- find new hydraulic state if system has changed + if ( inSteadyState == FALSE ) { - node_setOldHydState(j); - node_initInflow(j, routingStep); + // --- replace old hydraulic state values with current ones + for (j = 0; j < Nobjects[LINK]; j++) link_setOldHydState(j); + for (j = 0; j < Nobjects[NODE]; j++) + { + node_setOldHydState(j); + node_initInflow(j, routingStep); + } + + // --- route flow through the drainage network + if ( Nobjects[LINK] > 0 ) + { + stepCount = flowrout_execute(SortedLinks, routingModel, routingStep); + } } - // --- route flow through the drainage network - if ( Nobjects[LINK] > 0 ) + // --- route quality through the drainage network + if ( Nobjects[POLLUT] > 0 && !IgnoreQuality ) { - stepCount = flowrout_execute(SortedLinks, routingModel, routingStep); + qualrout_execute(routingStep); } - } - // --- route quality through the drainage network - if ( Nobjects[POLLUT] > 0 && !IgnoreQuality ) - { - qualrout_execute(routingStep); + // --- remove evaporation, infiltration & outflows from system + removeStorageLosses(routingStep); + removeConduitLosses(); + removeOutflows(routingStep); } - - // --- remove evaporation, infiltration & outflows from system - removeStorageLosses(routingStep); - removeConduitLosses(); - removeOutflows(routingStep); //(5.1.008) + else inSteadyState = TRUE; // --- update continuity with new totals // applied over 1/2 of routing step @@ -316,8 +324,8 @@ void routing_execute(int routingModel, double routingStep) // --- update summary statistics if ( RptFlags.flowStats && Nobjects[LINK] > 0 ) { - stats_updateFlowStats(routingStep, getDateTime(NewRoutingTime), //(5.1.008) - stepCount, InSteadyState); //(5.1.011) + stats_updateFlowStats(routingStep, getDateTime(NewRoutingTime), + stepCount, inSteadyState); } } diff --git a/src/runoff.c b/src/runoff.c index 4b4d760cc..fb9401bc0 100644 --- a/src/runoff.c +++ b/src/runoff.c @@ -7,6 +7,7 @@ // 09/15/14 (Build 5.1.007) // 03/19/15 (Build 5.1.008) // 08/01/16 (Build 5.1.011) +// 03/14/17 (Build 5.1.012) // Author: L. Rossman // M. Tryby // @@ -25,6 +26,9 @@ // Build 5.1.011: // - Runoff wet time step kept aligned with reporting times. // - Prior runoff time step used to convert returned outfall volume to flow. +// +// Build 5.1.012: +// - Runoff wet time step no longer kept aligned with reporting times. //----------------------------------------------------------------------------- #define _CRT_SECURE_NO_DEPRECATE @@ -310,14 +314,11 @@ double runoff_getTimeStep(DateTime currentDate) if ( timeStep > 0 && timeStep < maxStep ) maxStep = timeStep; } -//// Following code segment modified for release 5.1.011. //// //(5.1.011) +//// Following code segment modified for release 5.1.012. //// //(5.1.012) // --- determine whether wet or dry time step applies if ( IsRaining || HasSnow || HasRunoff || HasWetLids ) { - // --- when wet stay aligned with report time - timeStep = datetime_timeDiff(getDateTime(ReportTime), currentDate); - if ( timeStep > 0 ) timeStep = MIN(timeStep, WetStep); - else timeStep = WetStep; + timeStep = WetStep; } //// else timeStep = DryStep; diff --git a/src/stats.c b/src/stats.c index 4bd9c2168..6b7ac06a0 100644 --- a/src/stats.c +++ b/src/stats.c @@ -7,6 +7,7 @@ // 09/15/14 (Build 5.1.007) // 03/19/15 (Build 5.1.008) // 08/01/16 (Build 5.1.011) +// 03/14/17 (Build 5.1.012) // Author: L. Rossman (EPA) // R. Dickinson (CDM) // @@ -25,6 +26,10 @@ // - Surcharging is now evaluated only under dynamic wave flow routing and // storage nodes cannot be classified as surcharged. // +// Build 5.1.012: +// - Time step statistics now evaluated only in non-steady state periods. +// - Check for full conduit flow now accounts for number of barrels. +// //----------------------------------------------------------------------------- #define _CRT_SECURE_NO_DEPRECATE @@ -430,21 +435,28 @@ void stats_updateFlowStats(double tStep, DateTime aDate, int stepCount, stats_updateLinkStats(j, tStep, aDate); } - // --- update time step stats - // (skip initial time step for min. value) - if ( OldRoutingTime > 0 ) //(5.1.008) - { - SysStats.minTimeStep = MIN(SysStats.minTimeStep, tStep); - } - SysStats.avgTimeStep += tStep; - SysStats.maxTimeStep = MAX(SysStats.maxTimeStep, tStep); - - // --- update iteration step count stats - SysStats.avgStepCount += stepCount; +//// Following code segment modified for release 5.1.012. //// //(5.1.012) // --- update count of times in steady state SysStats.steadyStateCount += steadyState; + // --- update time step stats if not in steady state + if ( steadyState == FALSE ) + { + // --- skip initial time step for min. value) + if ( OldRoutingTime > 0 ) + { + SysStats.minTimeStep = MIN(SysStats.minTimeStep, tStep); + } + SysStats.avgTimeStep += tStep; + SysStats.maxTimeStep = MAX(SysStats.maxTimeStep, tStep); + + // --- update iteration step count stats + SysStats.avgStepCount += stepCount; + } + +//// + // --- update max. system outfall flow MaxOutfallFlow = MAX(MaxOutfallFlow, SysOutfallFlow); } @@ -645,7 +657,8 @@ void stats_updateLinkStats(int j, double tStep, DateTime aDate) // --- update time conduit is full k = Link[j].subIndex; - if ( q >= Link[j].qFull ) LinkStats[j].timeFullFlow += tStep; + if ( q >= Link[j].qFull * (double)Conduit[k].barrels ) //(5.1.012) + LinkStats[j].timeFullFlow += tStep; if ( Conduit[k].capacityLimited ) LinkStats[j].timeCapacityLimited += tStep; diff --git a/src/statsrpt.c b/src/statsrpt.c index dbb3d5e54..da1626b0d 100644 --- a/src/statsrpt.c +++ b/src/statsrpt.c @@ -21,6 +21,7 @@ // // Build 5.1.011: // - Redundant units conversion on max. reported node depth removed. +// - Node Surcharge table only produced for dynamic wave routing. //----------------------------------------------------------------------------- #define _CRT_SECURE_NO_DEPRECATE @@ -104,7 +105,7 @@ void statsrpt_writeReport() { writeNodeDepths(); writeNodeFlows(); - writeNodeSurcharge(); + if ( RouteModel == DW ) writeNodeSurcharge(); //(5.1.011) writeNodeFlooding(); writeStorageVolumes(); writeOutfallLoads(); diff --git a/src/subcatch.c b/src/subcatch.c index 4867f8101..1a03cd46c 100644 --- a/src/subcatch.c +++ b/src/subcatch.c @@ -9,6 +9,7 @@ // 04/30/15 (Build 5.1.009) // 08/05/15 (Build 5.1.010) // 08/01/16 (Build 5.1.011) +// 03/14/17 (Build 5.1.012) // Author: L. Rossman // // Subcatchment runoff functions. @@ -32,6 +33,10 @@ // Build 5.1.011: // - Subcatchment percent imperviousness not allowed to exceed 100. // +// Build 5.1.012: +// - Subcatchment bottom elevation used instead of aquifer's when +// saving water table value to results file. +// //----------------------------------------------------------------------------- #define _CRT_SECURE_NO_DEPRECATE @@ -874,7 +879,7 @@ void subcatch_getResults(int j, double f, float x[]) { z = (f1 * gw->oldFlow + f * gw->newFlow) * Subcatch[j].area * UCF(FLOW); x[SUBCATCH_GW_FLOW] = (float)z; - z = (Aquifer[gw->aquifer].bottomElev + gw->lowerDepth) * UCF(LENGTH); + z = (gw->bottomElev + gw->lowerDepth) * UCF(LENGTH); //(5.1.012) x[SUBCATCH_GW_ELEV] = (float)z; z = gw->theta; x[SUBCATCH_SOIL_MOIST] = (float)z; diff --git a/src/swmm5.c b/src/swmm5.c index 57791090c..41fa120ee 100644 --- a/src/swmm5.c +++ b/src/swmm5.c @@ -6,6 +6,7 @@ // Date: 03/19/14 (Build 5.1.001) // 03/19/15 (Build 5.1.008) // 08/01/16 (Build 5.1.011) +// 03/14/17 (Build 5.1.012) // Author: L. Rossman // // This is the main module of the computational engine for Version 5 of @@ -31,6 +32,9 @@ // - Changed WarningCode to Warnings (# warnings issued). // - Added swmm_getWarnings() function to retrieve value of Warnings. // - Fixed error code returned on swmm_xxx functions. +// +// Build 5.1.012: +// - #include only used when compiled for Windows. // //----------------------------------------------------------------------------- #define _CRT_SECURE_NO_DEPRECATE @@ -66,13 +70,13 @@ // --- include Windows & exception handling headers #ifdef WINDOWS #include + #include //(5.1.012) #endif #ifdef EXH #include #endif //// -#include #include #include #include diff --git a/src/text.h b/src/text.h index ab74575b1..05075db9d 100644 --- a/src/text.h +++ b/src/text.h @@ -13,6 +13,7 @@ // 04/30/15 (Build 5.1.009) // 08/05/15 (Build 5.1.010) // 08/01/16 (Build 5.1.011) +// 03/14/17 (Build 5.1.012) // Author: L. Rossman // // Text strings @@ -20,7 +21,7 @@ #define FMT01 \ "\n Correct syntax is:\n swmm5 \n" -#define FMT02 "\n... EPA-SWMM 5.1 (Build 5.1.011)\n" //(5.1.011) +#define FMT02 "\n... EPA-SWMM 5.1 (Build 5.1.012)\n" //(5.1.012) #define FMT03 " There are errors.\n" #define FMT04 " There are warnings.\n" @@ -28,7 +29,7 @@ #define FMT06 "\n o Retrieving project data" #define FMT07 "\n o Writing output report" #define FMT08 \ - "\n EPA STORM WATER MANAGEMENT MODEL - VERSION 5.1 (Build 5.1.011)" //(5.1.011) + "\n EPA STORM WATER MANAGEMENT MODEL - VERSION 5.1 (Build 5.1.012)" //(5.1.012) #define FMT09 \ "\n --------------------------------------------------------------" #define FMT10 "\n" diff --git a/src/xsect.c b/src/xsect.c index de103ec91..f09b2a35c 100644 --- a/src/xsect.c +++ b/src/xsect.c @@ -4,6 +4,7 @@ // Project: EPA SWMM5 // Version: 5.1 // Date: 03/20/14 (Build 5.1.001) +// 03/14/17 (Build 5.1.012) // Author: L. Rossman (EPA) // M. Tryby (EPA) // @@ -23,6 +24,9 @@ // A = flow area // R = hyd. radius // S = section factor = A*R^(2/3) +// +// Build 5.1.012: +// - Height at max. width for Modified Baskethandle shape corrected. //----------------------------------------------------------------------------- #define _CRT_SECURE_NO_DEPRECATE @@ -442,7 +446,7 @@ int xsect_setParams(TXsect *xsect, int type, double p[], double ucf) // --- height of circular arc xsect->yBot = xsect->rBot * (1.0 - cos(theta/2.0)); - xsect->ywMax = xsect->yBot; + xsect->ywMax = xsect->yFull - xsect->yBot; //(5.1.012) // --- area of circular arc xsect->aBot = xsect->rBot * xsect->rBot / diff --git a/tools/nrtest-swmm/nrtest_swmm/__init__.py b/tools/nrtest-swmm/nrtest_swmm/__init__.py index 4748a1436..dd8abd1c6 100644 --- a/tools/nrtest-swmm/nrtest_swmm/__init__.py +++ b/tools/nrtest-swmm/nrtest_swmm/__init__.py @@ -1,4 +1,8 @@ # -*- coding: utf-8 -*- +''' +Numerical regression testing (nrtest) plugin for comparing SWMM binary results +files and SWMM text based report files. +''' # system imports import itertools as it @@ -7,10 +11,23 @@ import header_detail_footer as hdf import numpy as np -# project import +# project imports import swmm_reader as sr +__author__ = "Michael Tryby" +__copyright__ = "None" +__credits__ = "Colleen Barr, Maurizio Cingi, Mark Gray, David Hall, Bryant McDonnell" +__license__ = "CC0 1.0 Universal" + +__version__ = "0.2.0" +__date__ = "September 20, 2016" + +__maintainer__ = "Michael Tryby" +__email__ = "tryby.michael@epa.gov" +__status = "Development" + + def swmm_allclose_compare(path_test, path_ref, rtol, atol): ''' Compares results in two SWMM binary files. Using the comparison criteria @@ -25,18 +42,6 @@ def swmm_allclose_compare(path_test, path_ref, rtol, atol): are checked to see if they are equal before being compared using the allclose criteria. This reduces comparison times significantly. - Allclose raises the question, "To what values should the relative and - absolute tolerances be set?" In theory, this turns out to be a difficult - question to answer. - - The choice of the allclose critera itself implies a pragmatic approach to - answering this question. In practice, the tolerances values should be - selected based on the purpose for running the test. - - Values selected to screen for unanticipated changes in results can be more - stringent than those for engine development work that is likely to - significantly change benchmark results. - Arguments: path_test - path to result file being tested path_ref - path to reference result file @@ -51,7 +56,8 @@ def swmm_allclose_compare(path_test, path_ref, rtol, atol): AssertionError() ... ''' - for (test, ref) in it.izip(sr.reader(path_test), sr.reader(path_ref)): + for (test, ref) in it.izip(sr.swmm_output_generator(path_test), + sr.swmm_output_generator(path_ref)): if test.size != ref.size: raise ValueError('Inconsistent lengths') diff --git a/tools/swmm-output/src/outputapi.c b/tools/swmm-output/src/outputapi.c index 7aecab509..fe5f0372f 100644 --- a/tools/swmm-output/src/outputapi.c +++ b/tools/swmm-output/src/outputapi.c @@ -2,7 +2,8 @@ * outputAPI.c * * Author: Colleen Barr -* Modified by: Michael E. Tryby +* Modified by: Michael E. Tryby, +* Bryant McDonnell * * */ @@ -197,14 +198,13 @@ int DLLEXPORT SMO_getUnits(SMOutputAPI* smoapi, SMO_unit code, int* unitFlag) // // Purpose: Returns flow rate units. // -// Note: Concentration units are located after the pollutant ID names and before the object properties start, -// and can differ for each pollutant. They're stored as 4-byte integers with the following codes: -// 0: mg/L -// 1: ug/L -// 2: counts/L -// Probably the best way to do this would not be here -- instead write a function that takes -// NPolluts and ObjPropPos, jump to ObjPropPos, count backward (NPolluts * 4), then read forward -// to get the units for each pollutant +// Returns: +// 0: CFS (cubic feet per second) +// 1: GPM (gallons per minute) +// 2: MGD (million gallons per day) +// 3: CMS (cubic meters per second) +// 4: LPS (liters per second) +// 5: MLD (million liters per day) // { int errorcode = 0; @@ -227,6 +227,33 @@ int DLLEXPORT SMO_getUnits(SMOutputAPI* smoapi, SMO_unit code, int* unitFlag) return errorcode; } + +int DLLEXPORT SMO_getPollutantUnits(SMOutputAPI* smoapi, int pollutantIndex, int* unitFlag) +// +// Purpose: Return integer flag representing the units that the given pollutant is measured in. +// Concentration units are located after the pollutant ID names and before the object properties start, +// and are stored for each pollutant. They're stored as 4-byte integers with the following codes: +// 0: mg/L +// 1: ug/L +// 2: count/L +// +// pollutantIndex: valid values are 0 to Npolluts-1 +{ + int errorcode = 0; + if (smoapi == NULL) errorcode = 410; + else if (smoapi->file == NULL) errorcode = 411; + else if (pollutantIndex < 0 || pollutantIndex >= smoapi->Npolluts) errorcode = 423; + else + { + int offset = smoapi->ObjPropPos - (smoapi->Npolluts - pollutantIndex) * RECORDSIZE; + fseek(smoapi->file, offset, SEEK_SET); + fread(unitFlag, RECORDSIZE, 1, smoapi->file); + } + + return errorcode; +} + + int DLLEXPORT SMO_getStartTime(SMOutputAPI* smoapi, double* time) // // Purpose: Returns start date. diff --git a/tools/swmm-output/src/outputapi.h b/tools/swmm-output/src/outputapi.h index 0e7dfc0c2..2428fa0df 100644 --- a/tools/swmm-output/src/outputapi.h +++ b/tools/swmm-output/src/outputapi.h @@ -2,7 +2,8 @@ * outputAPI.h * * Author: Colleen Barr -* Modified by: Michael E. Tryby +* Modified by: Michael Tryby, +* Bryant McDonnell * * */ @@ -10,7 +11,7 @@ #ifndef OUTPUTAPI_H_ #define OUTPUTAPI_H_ -#define MAXFILENAME 259 // +#define MAXFILENAME 259 // Max characters in file path #define MAXELENAME 31 // Max characters in element name /*------------------- Error Messages --------------------*/ @@ -30,6 +31,7 @@ #define ERR440 "ERROR 440: an unspecified error has occurred" + typedef struct SMOutputAPI SMOutputAPI; // opaque pointer typedef enum { @@ -138,7 +140,8 @@ typedef enum { SMOutputAPI* DLLEXPORT SMO_init(void); int DLLEXPORT SMO_open(SMOutputAPI* smoapi, const char* path); -int DLLEXPORT SMO_getProjectSize(SMOutputAPI* smoapi, SMO_elementCount code, int* count); +int DLLEXPORT SMO_getProjectSize(SMOutputAPI* smoapi, SMO_elementCount code, + int* count); int DLLEXPORT SMO_getUnits(SMOutputAPI* smoapi, SMO_unit code, int* unitFlag); int DLLEXPORT SMO_getStartTime(SMOutputAPI* smoapi, double* time); @@ -150,7 +153,7 @@ int DLLEXPORT SMO_getElementName(SMOutputAPI* smoapi, SMO_elementType type, float* DLLEXPORT SMO_newOutValueSeries(SMOutputAPI* smoapi, long seriesStart, long seriesLength, long* length, int* errcode); float* DLLEXPORT SMO_newOutValueArray(SMOutputAPI* smoapi, SMO_apiFunction func, - SMO_elementType type, long* length, int* errcode); + SMO_elementType type, long* length, int* errcode); int DLLEXPORT SMO_getSubcatchSeries(SMOutputAPI* smoapi, int subcatchIndex, SMO_subcatchAttribute attr, long timeIndex, long length, float* outValueSeries); diff --git a/tools/swmm-output/test/main.c b/tools/swmm-output/test/main.c index b7baa3aa0..6584c237b 100644 --- a/tools/swmm-output/test/main.c +++ b/tools/swmm-output/test/main.c @@ -1,7 +1,8 @@ /* * main.c -* TODO -* Error handling from getID functions (and check other error handling) +* +* Author: Colleen Barr +* Modified by: Michael E. Tryby * */ diff --git a/tools/swmm-reader/main.py b/tools/swmm-reader/main.py index 88fdb6d59..25571bc89 100644 --- a/tools/swmm-reader/main.py +++ b/tools/swmm-reader/main.py @@ -1,9 +1,14 @@ +# -*- coding: utf-8 -*- +''' +Provides entry point main. Useful for development and testing purposes. +''' -import header_detail_footer as hdf +import cStringIO import itertools as it -import numpy as np import time -import cStringIO + +import header_detail_footer as hdf +import numpy as np import swmm_reader as sr @@ -20,8 +25,8 @@ def result_compare(path_test, path_ref, comp_args): start = time.time() - test_reader = sr.reader(path_test) - ref_reader = sr.reader(path_ref) + test_reader = sr.swmm_output_generator(path_test) + ref_reader = sr.swmm_output_generator(path_ref) for test, ref in it.izip(test_reader, ref_reader): total += 1 @@ -37,7 +42,7 @@ def result_compare(path_test, path_ref, comp_args): continue else: try: - np.testing.assert_allclose(test, ref, 1.0e-06, 2*eps) + np.testing.assert_allclose(test, ref, comp_args[0], comp_args[1]) close += 1 except AssertionError as ae: @@ -78,10 +83,11 @@ def report_compare(path_test, path_ref, (comp_args)): import logging -from nrtest.testsuite import TestSuite -from nrtest.compare import compare_testsuite, validate_testsuite from os import listdir from os.path import exists, isfile, isdir, join + +from nrtest.testsuite import TestSuite +from nrtest.compare import compare_testsuite, validate_testsuite from nrtest.execute import execute_testsuite def nrtest_compare(path_test, path_ref, rtol, atol): @@ -151,9 +157,9 @@ def nrtest_execute(app_path, test_path, output_path): # path_test = "C:\\Users\\mtryby\\Workspace\\GitRepo\\Local\\swmm-testsuite\\benchmarks\\v5111\\Example_4\\Example4.out" # path_ref = "C:\\Users\\mtryby\\Workspace\\GitRepo\\Local\\swmm-testsuite\\benchmarks\\v5110\\Example_4\\Example4.out" # result_compare(path_test, path_ref, (0.001, 0.0)) - - path_test = "C:\\Users\\mtryby\\Workspace\\GitRepo\\Local\\swmm-testsuite\\benchmarks\\v5111_x64\\gate_control_3\\gate_control_3.out" - path_ref = "C:\\Users\\mtryby\\Workspace\\GitRepo\\Local\\swmm-testsuite\\benchmarks\\v5111_x86\\gate_control_3\\gate_control_3.out" - print(result_compare(path_test, path_ref, (0.0001, 0.0))) + + path_test = "C:\\Users\\mtryby\\Workspace\\GitRepo\\Local\\swmm-testsuite\\benchmarks\\v5111_bf\\lid_cat\\lid_cat.out" + path_ref = "C:\\Users\\mtryby\\Workspace\\GitRepo\\Local\\swmm-testsuite\\benchmarks\\v5110\\lid_cat\\lid_cat.out" + print(result_compare(path_test, path_ref, (1.0, 0.0))) \ No newline at end of file diff --git a/tools/swmm-reader/swmm_reader/__init__.py b/tools/swmm-reader/swmm_reader/__init__.py index 43bbdc7e6..a85c0a0e1 100644 --- a/tools/swmm-reader/swmm_reader/__init__.py +++ b/tools/swmm-reader/swmm_reader/__init__.py @@ -1,28 +1,43 @@ # -*- coding: utf-8 -*- +''' +The function swmm_output_generator is used to iterate over a swmm binary file. +''' # project import from .swmm_reader import * -def reader(path_ref): + +__author__ = "Michael Tryby" +__copyright = "None" +__credits__ = "Colleen Barr, Maurizio Cingi, Mark Gray, David Hall, Bryant McDonnell" +__license__ = "CC0 1.0 Universal" + +__version__ = "0.2.0" +__maintainer__ = "Michael Tryby" +__email__= "tryby.michael@epa.gov" +__status = "Development" + + +def swmm_output_generator(path_ref): ''' - The function reader() is a generator designed to iterate through a swmm - binary file and yield element results. It is useful for comparing contents - of binary files for numerical regression testing. + swmm_output_generator is designed to iterate over a swmm binary file and + yield element results. It is useful for comparing contents of binary files + for numerical regression testing. - The reader yields a numpy array containing the SWMM element result. + The generator yields a numpy array containing the SWMM element result. Arguments: path_ref - path to result file Raises: - SWMM_Binary_Reader_Error() + SWMM_OutputReaderError() ... ''' - with SWMM_BinaryReader(path_ref) as br: + with SWMM_OutputReader(path_ref) as sor: - for period_index in range(0, br.report_periods()): + for period_index in range(0, sor.report_periods()): for element_type in ElementType: - for element_index in range(0, br.element_count(element_type)): + for element_index in range(0, sor.element_count(element_type)): - yield br.element_result(element_type, period_index, element_index) + yield sor.element_result(element_type, period_index, element_index) \ No newline at end of file diff --git a/tools/swmm-reader/swmm_reader/outputapi.py b/tools/swmm-reader/swmm_reader/outputapi.py index 769a8fabb..cbcf4df49 100644 --- a/tools/swmm-reader/swmm_reader/outputapi.py +++ b/tools/swmm-reader/swmm_reader/outputapi.py @@ -3,11 +3,11 @@ '''Wrapper for outputapi.h Generated with: -C:\Users\mtryby\dev\Anaconda2\Scripts\ctypesgen.py -a -l libswmm-output -L C:\Users\mtryby\Workspace\GitRepo\Local\swmm-output\Release\ -o outputapi.py .\src\outputapi.h +C:\Users\mtryby\dev\Anaconda2\Scripts\ctypesgen.py -a -l libswmm-output -o outputapi.py .\src\outputapi.h Do not modify this file ... -This file HAS been modified! Path to dll was add because dll search is broken on Winblows (See line 592) +THIS FILE HAS BEEN MODIFIED! See WindowsLibraryLoader.genplatformpaths() starting on line 560 ''' __docformat__ = 'restructuredtext' @@ -307,7 +307,7 @@ def __call__(self,*args): # End preamble _libs = {} -_libdirs = ['C:\\Users\\mtryby\\Workspace\\GitRepo\\Local\\swmm-output\\Release\\'] +_libdirs = [] # Begin loader @@ -554,16 +554,31 @@ def load_library(self, libname): def load(self, path): return _WindowsLibrary(path) - + +############################################################################### +###### THIS NEEDS TO BE HAND CODED OR LIBRARY LOAD WILL BREAK ON WINDOWS ###### def getplatformpaths(self, libname): if os.path.sep not in libname: for name in self.name_formats: + # search the directory where this file is executing + dll_in_package_dir = os.path.join(os.path.dirname(os.path.realpath(__file__)), name % libname) + if os.path.exists(dll_in_package_dir): + yield dll_in_package_dir + # search in directories passed as arguments to wrapper generator + for dir in self.other_dirs: + dll_in_other_dirs = os.path.join(dir, name % libname) + if os.path.exists(dll_in_other_dirs): + yield dll_in_other_dirs + # searches current working directory dll_in_current_dir = os.path.abspath(name % libname) if os.path.exists(dll_in_current_dir): - yield dll_in_current_dir + yield dll_in_current_dir + # searches system path path = ctypes.util.find_library(name % libname) if path: yield path +############################################################################### +############################################################################### # Platform switching @@ -590,12 +605,8 @@ def add_library_search_dirs(other_dirs): add_library_search_dirs([]) # Begin libraries -############################################################################### -###### THIS NEEDS TO BE HAND CODED OR LIBRARY LOAD WILL BREAK ON WINDOWS ###### -import site -_libs["libswmm-output"] = load_library(site.getsitepackages()[1] + "\swmm_reader-0.2.0-py2.7.egg\swmm_reader\libswmm-output") -############################################################################### -############################################################################### + +_libs["libswmm-output"] = load_library("libswmm-output") # 1 libraries # End libraries diff --git a/tools/swmm-reader/swmm_reader/swmm_reader.py b/tools/swmm-reader/swmm_reader/swmm_reader.py index cd5582c2b..65579a02d 100644 --- a/tools/swmm-reader/swmm_reader/swmm_reader.py +++ b/tools/swmm-reader/swmm_reader/swmm_reader.py @@ -1,4 +1,8 @@ # -*- coding: utf-8 -*- +''' +The module swmm_reader provides classes used to implement the swmm output +generator. +''' # system imports import ctypes @@ -8,18 +12,29 @@ import numpy as np # project import -import outputapi as oapi +import outputapi + + +__author__ = "Michael Tryby" +__copyright = "None" +__credits__ = "Colleen Barr, Maurizio Cingi, Mark Gray, David Hall, Bryant McDonnell" +__license__ = "CC0 1.0 Universal" + +__version__ = "0.2.0" +__maintainer__ = "Michael Tryby" +__email__= "tryby.michael@epa.gov" +__status = "Development" _err_max_char = 80 class ElementType(enum.Enum): - subcatch = oapi.subcatch - node = oapi.node - link = oapi.link - system = oapi._sys + subcatch = outputapi.subcatch + node = outputapi.node + link = outputapi.link + system = outputapi._sys -class SWMM_BinaryReaderError(Exception): +class SWMM_OutputReaderError(Exception): ''' Custom exception class for SWMM errors. ''' @@ -31,55 +46,55 @@ def __init__(self, error_code, error_message): def __str__(self): return self.message -class SWMM_BinaryReader(): +class SWMM_OutputReader(): ''' - Provides a minimal API used to implement the SWMM result generator. + Provides minimal functionality needed to implement the SWMM output generator. ''' def __init__(self, filename): self.filepath = filename self.ptr_api = ctypes.c_void_p self.ptr_resultbuff = ctypes.c_void_p self.bufflength = ctypes.c_long() - self.getElementResult = {ElementType.subcatch: oapi.SMO_getSubcatchResult, - ElementType.node: oapi.SMO_getNodeResult, - ElementType.link: oapi.SMO_getLinkResult, - ElementType.system: oapi.SMO_getSystemResult} + self.getElementResult = {ElementType.subcatch: outputapi.SMO_getSubcatchResult, + ElementType.node: outputapi.SMO_getNodeResult, + ElementType.link: outputapi.SMO_getLinkResult, + ElementType.system: outputapi.SMO_getSystemResult} def __enter__(self): - self.ptr_api = oapi.SMO_init() - self._error_check(oapi.SMO_open(self.ptr_api, ctypes.c_char_p(self.filepath.encode()))) + self.ptr_api = outputapi.SMO_init() + self._error_check(outputapi.SMO_open(self.ptr_api, ctypes.c_char_p(self.filepath.encode()))) # max system result is vector with 15 elements so should be adequate for result buffer # TODO: What about when there are more than six pollutants defined? error = ctypes.c_long() - self.ptr_resultbuff = oapi.SMO_newOutValueArray(self.ptr_api, ctypes.c_int(oapi.getResult), + self.ptr_resultbuff = outputapi.SMO_newOutValueArray(self.ptr_api, ctypes.c_int(outputapi.getResult), ctypes.c_int(ElementType.system.value), ctypes.byref(self.bufflength), ctypes.byref(error)) self._error_check(error.value) return self def __exit__(self, type, value, traceback): - oapi.SMO_free(self.ptr_resultbuff) - self._error_check(oapi.SMO_close(self.ptr_api)) + outputapi.SMO_free(self.ptr_resultbuff) + self._error_check(outputapi.SMO_close(self.ptr_api)) def _error_message(self, code): error_code = ctypes.c_int(code) - error_message = oapi.String(ctypes.create_string_buffer(_err_max_char)) - oapi.SMO_errMessage(error_code, error_message, _err_max_char) + error_message = outputapi.String(ctypes.create_string_buffer(_err_max_char)) + outputapi.SMO_errMessage(error_code, error_message, _err_max_char) return error_message.data def _error_check(self, err): if err != 0: - raise SWMM_BinaryReaderError(err, self._error_message(err)) + raise SWMM_OutputReaderError(err, self._error_message(err)) def report_periods(self): num_periods = ctypes.c_int() - self._error_check(oapi.SMO_getTimes(self.ptr_api, oapi.numPeriods, ctypes.byref(num_periods))) + self._error_check(outputapi.SMO_getTimes(self.ptr_api, outputapi.numPeriods, ctypes.byref(num_periods))) return num_periods.value def element_count(self, element_type): count = ctypes.c_int() - self._error_check(oapi.SMO_getProjectSize(self.ptr_api, ctypes.c_int(element_type.value), ctypes.byref(count))) + self._error_check(outputapi.SMO_getProjectSize(self.ptr_api, ctypes.c_int(element_type.value), ctypes.byref(count))) return count.value def element_result(self, element_type, time_index, element_index):