diff --git a/src/epanet.c b/src/epanet.c index a7ee920..c9248b9 100755 --- a/src/epanet.c +++ b/src/epanet.c @@ -1668,9 +1668,9 @@ int DLLEXPORT ENgetlinkvalue(int index, int code, float *value) case EN_FLOW: /*** Updated 10/25/00 ***/ - if (S[index] <= CLOSED) v = 0.0; + //if (S[index] <= CLOSED) v = 0.0; - else v = Q[index]*Ucf[FLOW]; + /*else*/ v = Q[index]*Ucf[FLOW]; break; case EN_VELOCITY: diff --git a/src/hydraul.c b/src/hydraul.c index 8a91d6d..6c138f9 100755 --- a/src/hydraul.c +++ b/src/hydraul.c @@ -258,6 +258,9 @@ int nexthyd(long *tstep) else { Htime++; /* Force completion of analysis */ + if (OpenQflag) { + Qtime++; // force completion of wq analysis too + } } *tstep = hydstep; return(errcode); @@ -1049,13 +1052,8 @@ void tanklevels(long tstep) /* Update the tank's volume & water elevation */ n = Tank[i].Node; - - // only adjust tank volume if we're in sequential mode. - // otherwise, the tankmixing function will do it for us. - if (!OpenQflag) { - dv = D[n]*tstep; - Tank[i].V += dv; - } + dv = D[n]*tstep; + Tank[i].V += dv; /*** Updated 6/24/02 ***/ /* Check if tank full/empty within next second */ diff --git a/src/output.c b/src/output.c index ea56057..048a27b 100755 --- a/src/output.c +++ b/src/output.c @@ -161,8 +161,8 @@ int savehyd(long *htime) /* Force flow in closed links to be zero then save flows */ for (i=1; i<=Nlinks; i++) { - if (S[i] <= CLOSED) x[i] = 0.0f; - else x[i] = (REAL4)Q[i]; + //if (S[i] <= CLOSED) x[i] = 0.0f; + /*else*/ x[i] = (REAL4)Q[i]; } fwrite(x+1,sizeof(REAL4),Nlinks,HydFile); @@ -377,7 +377,7 @@ int nodeoutput(int j, REAL4 *x, double ucf) } /* End of nodeoutput */ -int linkoutput(int j, float *x, double ucf) +int linkoutput(int j, REAL4 *x, double ucf) /* **---------------------------------------------------------------- ** Input: j = type of link variable diff --git a/src/quality.c b/src/quality.c index cacbbe9..7311818 100755 --- a/src/quality.c +++ b/src/quality.c @@ -152,7 +152,9 @@ void initqual() for (i=1; i<=Ntanks; i++) Tank[i].V = Tank[i].V0; for (i=1; i<=Nnodes; i++) if (Node[i].S != NULL) Node[i].S->Smass = 0.0; - + + QTankVolumes = calloc(Ntanks, sizeof(double)); // keep track of previous step's tank volumes. + /* Set WQ parameters */ Bucf = 1.0; Tucf = 1.0; @@ -224,8 +226,15 @@ int runqual(long *t) { errcode = gethyd(&hydtime, &hydstep); if (!OpenHflag) { // test for sequential vs stepwise + // sequential Htime = hydtime + hydstep; } + else { + // stepwise + for (int i=1; i<= Ntanks; ++i) { + QTankVolumes[i-1] = Tank[i].V; + } + } } return(errcode); } @@ -249,6 +258,26 @@ int nextqual(long *tstep) *tstep = 0; hydstep = Htime - Qtime; + double *tankVolumes; + // if we're operating in stepwise mode, capture the tank levels so we can restore them later. + if (OpenHflag) { + tankVolumes = calloc(Ntanks, sizeof(double)); + for (int i=1; i<=Ntanks; ++i) { + if (Tank[i].A != 0) { // skip reservoirs + tankVolumes[i-1] = Tank[i].V; + } + } + // restore the previous step's tank volumes + for (int i=1; i<=Ntanks; i++) { + if (Tank[i].A != 0) { // skip reservoirs again + int n = Tank[i].Node; + Tank[i].V = QTankVolumes[i-1]; + H[n] = tankgrade(i,Tank[i].V); + } + } + } + + /* Perform water quality routing over this time step */ if (Qualflag != NONE && hydstep > 0) transport(hydstep); @@ -259,6 +288,19 @@ int nextqual(long *tstep) /* Save final output if no more time steps */ if (!errcode && Saveflag && *tstep == 0) errcode = savefinaloutput(); + + // restore tank levels to post-runH state, if needed. + if (OpenHflag) { + for (int i=1; i<=Ntanks; i++) { + if (Tank[i].A != 0) { // skip reservoirs again + int n = Tank[i].Node; + Tank[i].V = tankVolumes[i-1]; + H[n] = tankgrade(i,Tank[i].V); + } + } + free(tankVolumes); + } + return(errcode); } @@ -326,6 +368,7 @@ int closequal() free(MassIn); free(R); free(XC); + free(QTankVolumes); return(errcode); } @@ -423,6 +466,8 @@ void transport(long tstep) */ { long qtime, dt; + + /* Repeat until elapsed time equals hydraulic time step */ @@ -439,6 +484,8 @@ void transport(long tstep) release(dt); /* Release new nodal flows */ } updatesourcenodes(tstep); /* Update quality at source nodes */ + + } @@ -1023,10 +1070,6 @@ void updatetanks(long dt) default: tankmix1(i,dt); break; } - // if we're operating in stepwise mode, we'll need to update tank head conditions - if (OpenHflag) { - H[n] = tankgrade(i,Tank[i].V); - } } } } diff --git a/src/testLemonTiger.cpp b/src/testLemonTiger.cpp index e521019..7cba39f 100644 --- a/src/testLemonTiger.cpp +++ b/src/testLemonTiger.cpp @@ -1,6 +1,6 @@ #include #include - +#include #include "testLemonTiger.h" #include "toolkit.h" @@ -13,9 +13,19 @@ typedef struct { double head; double demand; double quality; -} singleState_t; +} nodeState_t; -typedef map networkState_t; // nodeIndex, state +typedef struct { + double flow; +} linkState_t; + +typedef map networkNodeState_t; // nodeIndex, state +typedef map networkLinkState_t; // linkIndex, state + +typedef struct { + networkNodeState_t nodeState; + networkLinkState_t linkState; +} networkState_t; typedef map result_t; // time, networkState // access results by, for instance, resultsContainer[time][nodeIndex].head @@ -25,6 +35,7 @@ void checkErr(int err, std::string function); void saveHydResults(networkState_t* networkState); void saveQualResults(networkState_t* networkState); void printResults(result_t* state1, result_t* state2, std::ostream& out); +void compare(result_t* results1, result_t* results2, std::ostream &out); int main(int argc, char * argv[]) { @@ -67,7 +78,7 @@ int main(int argc, char * argv[]) { cout << "Running WQ..." << endl; checkErr( ENopenQ(), "ENopenQ" ); - checkErr( ENinitQ(EN_SAVE), "ENinitQ" ); + checkErr( ENinitQ(EN_NOSAVE), "ENinitQ" ); do { @@ -105,12 +116,14 @@ int main(int argc, char * argv[]) { /* Solve for hydraulics & advance to next time period */ checkErr( ENrunH(&simulationTime), "ENrunH" ); checkErr( ENrunQ(&simulationTime), "ENrunQ" ); - + checkErr( ENnextH(&nextEventH), "ENnextH" ); checkErr( ENnextQ(&nextEventQ), "ENstepQ" ); + saveHydResults(&lemonTigerResults[simulationTime]); saveQualResults(&lemonTigerResults[simulationTime]); + } while (nextEventH > 0); cout << "\t\t\tdone." << endl; @@ -122,8 +135,8 @@ int main(int argc, char * argv[]) { // summarize the results - printResults(&epanetResults, &lemonTigerResults, cout); - + //printResults(&epanetResults, &lemonTigerResults, cout); + compare(&epanetResults, &lemonTigerResults, cout); } catch (int err) { cerr << "exiting with error " << err << endl; @@ -132,15 +145,19 @@ int main(int argc, char * argv[]) { void saveHydResults(networkState_t* networkState) { - int nNodes; - float head, demand; + int nNodes, nLinks; + float head, demand, flow; ENgetcount(EN_NODECOUNT, &nNodes); - - for (int iNode = 1; iNode <= nNodes; iNode++) { + ENgetcount(EN_LINKCOUNT, &nLinks); + for (int iNode = 1; iNode <= nNodes; ++iNode) { ENgetnodevalue(iNode, EN_HEAD, &head); ENgetnodevalue(iNode, EN_DEMAND, &demand); - (*networkState)[iNode].head = head; - (*networkState)[iNode].demand = demand; + (*networkState).nodeState[iNode].head = head; + (*networkState).nodeState[iNode].demand = demand; + } + for (int iLink = 1; iLink <= nLinks; ++iLink) { + ENgetlinkvalue(iLink, EN_FLOW, &flow); + (*networkState).linkState[iLink].flow = flow; } } @@ -152,7 +169,7 @@ void saveQualResults(networkState_t* networkState) { for (int iNode = 1; iNode <= nNodes; iNode++) { ENgetnodevalue(iNode, EN_QUALITY, &quality); - (*networkState)[iNode].quality = quality; + (*networkState).nodeState[iNode].quality = quality; } } @@ -164,7 +181,8 @@ void printResults(result_t* results1, result_t* results2, std::ostream &out) { for (resultIterator = (*results1).begin(); resultIterator != (*results1).end(); ++resultIterator) { // get the current frame const long time = resultIterator->first; - const networkState_t state1 = resultIterator->second; + const networkNodeState_t nodeState1 = resultIterator->second.nodeState; + const networkLinkState_t linkState1 = resultIterator->second.linkState; // see if this time is indexed in the second state container if ((*results2).find(time) == (*results2).end()) { @@ -173,8 +191,8 @@ void printResults(result_t* results1, result_t* results2, std::ostream &out) { } else { // get the second result set's state - const networkState_t state2 = (*results2)[time]; - + const networkNodeState_t networkNodeState2 = (*results2)[time].nodeState; + const networkLinkState_t networkLinkState2 = (*results2)[time].linkState; // print the current simulation time out << left; out << setfill('*') << setw(100) << "*" << endl; @@ -188,34 +206,117 @@ void printResults(result_t* results1, result_t* results2, std::ostream &out) { out << setprecision(OUTPRECISION); // loop through the nodes in the networkState objs, and print out the results for this time period - networkState_t::const_iterator networkIterator; - for (networkIterator = state1.begin(); networkIterator != state1.end(); ++networkIterator) { - int nodeIndex = networkIterator->first; + networkNodeState_t::const_iterator networkNodeIterator; + for (networkNodeIterator = nodeState1.begin(); networkNodeIterator != nodeState1.end(); ++networkNodeIterator) { + int nodeIndex = networkNodeIterator->first; // trusting that all nodes are present... - const singleState_t nodeState1 = networkIterator->second; - const singleState_t nodeState2 = state2.at(nodeIndex); - - // epanet - out << setw(10) << nodeIndex << "|"; - out << setw(COLW) << nodeState1.demand; - out << setw(COLW) << nodeState1.head; - out << setw(COLW) << nodeState1.quality; - - // lemontiger - out << "|"; - out << setw(COLW) << nodeState2.demand; - out << setw(COLW) << nodeState2.head; - out << setw(COLW) << nodeState2.quality; - out << endl; + const nodeState_t nodeState1 = networkNodeIterator->second; + const nodeState_t nodeState2 = networkNodeState2.at(nodeIndex); + if (nodeState1.quality != nodeState2.quality ) { + // epanet + out << setw(10) << nodeIndex << "|"; + out << setw(COLW) << nodeState1.demand; + out << setw(COLW) << nodeState1.head; + out << setw(COLW) << nodeState1.quality; + + // lemontiger + out << "|"; + out << setw(COLW) << nodeState2.demand; + out << setw(COLW) << nodeState2.head; + out << setw(COLW) << nodeState2.quality; + out << endl; + } } + networkLinkState_t::const_iterator networkLinkIterator; + for (networkLinkIterator = linkState1.begin(); networkLinkIterator != linkState1.end(); ++networkLinkIterator) { + int linkIndex = networkLinkIterator->first; + // trusting that all nodes are present... + const linkState_t linkState1 = networkLinkIterator->second; + const linkState_t linkState2 = networkLinkState2.at(linkIndex); + + if ( linkState1.flow != linkState2.flow ) { + // epanet + out << setw(10) << linkIndex << "|"; + out << setw(COLW) << linkState1.flow; + + // lemontiger + out << "|"; + out << setw(COLW) << linkState2.flow; + out << endl; + } + } + + } } } + +void compare(result_t* results1, result_t* results2, std::ostream &out) { + + double sumHeadDiff=0, sumDemandDiff=0, sumQualDiff=0, sumFlowDiff=0; + + result_t::const_iterator resultIterator; + + for (resultIterator = (*results1).begin(); resultIterator != (*results1).end(); ++resultIterator) { + // get the current frame + const long time = resultIterator->first; + const networkNodeState_t nodeState1 = resultIterator->second.nodeState; + const networkLinkState_t linkState1 = resultIterator->second.linkState; + + // see if this time is indexed in the second state container + if ((*results2).find(time) == (*results2).end()) { + // nope. + out << "time " << time << " not found in second result set" << endl; + } + else { + // get the second result set's state + const networkNodeState_t networkNodeState2 = (*results2)[time].nodeState; + const networkLinkState_t networkLinkState2 = (*results2)[time].linkState; + double qualD=0; + + networkNodeState_t::const_iterator networkNodeIterator; + for (networkNodeIterator = nodeState1.begin(); networkNodeIterator != nodeState1.end(); ++networkNodeIterator) { + int nodeIndex = networkNodeIterator->first; + // trusting that all nodes are present... + const nodeState_t nodeState1 = networkNodeIterator->second; + const nodeState_t nodeState2 = networkNodeState2.at(nodeIndex); + + sumHeadDiff += fabs(nodeState1.head - nodeState2.head); + sumDemandDiff += fabs(nodeState1.demand - nodeState2.demand); + + qualD += fabs(nodeState1.quality - nodeState2.quality); + } + //out << "T: " << time << " dq: " << setprecision(20) << qualD << endl; + sumQualDiff += qualD; + + networkLinkState_t::const_iterator networkLinkIterator; + for (networkLinkIterator = linkState1.begin(); networkLinkIterator != linkState1.end(); ++networkLinkIterator) { + int linkIndex = networkLinkIterator->first; + // trusting that all nodes are present... + const linkState_t linkState1 = networkLinkIterator->second; + const linkState_t linkState2 = networkLinkState2.at(linkIndex); + + sumFlowDiff += fabs(linkState1.flow - linkState2.flow); + } + } + } + + int c1 = 18; + int p = 20; + out << setw(c1) << "Head Diff:" << setprecision(p) << sumHeadDiff << endl; + out << setw(c1) << "Demand Diff:" << setprecision(p) << sumDemandDiff << endl; + out << setw(c1) << "Quality Diff:" << setprecision(p) << sumQualDiff << endl; + out << setw(c1) << "Flow Diff:" << setprecision(p) << sumFlowDiff << endl; + + +} + + void checkErr(int err, std::string function) { if (err > 0) { cerr << "Error in " << function << ": " << err << endl; diff --git a/src/types.h b/src/types.h index 3374f8f..343b11f 100755 --- a/src/types.h +++ b/src/types.h @@ -28,7 +28,7 @@ AUTHOR: L. Rossman Definition of 4-byte integers & reals ------------------------------------------- */ -typedef float REAL4; //(2.00.11 - LR) +typedef double REAL4; //(2.00.11 - LR) typedef int INT4; //(2.00.12 - LR) /* diff --git a/src/vars.h b/src/vars.h index b6d9b49..3ea8f49 100755 --- a/src/vars.h +++ b/src/vars.h @@ -153,6 +153,7 @@ AUTHOR: L. Rossman *X, /* General purpose array */ *XC; /* General Purpose array - WQ */ double *H; /* Node heads */ + double *QTankVolumes; STmplist *Patlist; /* Temporary time pattern list */ STmplist *Curvelist; /* Temporary list of curves */ Spattern *Pattern; /* Time patterns */ diff --git a/test/Net3.inp b/test/Net3.inp index e169140..8e39dd1 100644 --- a/test/Net3.inp +++ b/test/Net3.inp @@ -333,7 +333,7 @@ Link 330 OPEN IF Node 1 ABOVE 19.1 [TIMES] Duration 24:00 - Hydraulic Timestep 1:00 + Hydraulic Timestep 1:00 Quality Timestep 0:05 Pattern Timestep 1:00 Pattern Start 0:00