fake data was being written to the binary file. if a link element is nominally closed, even if there's actually flow calculated there, the flow is disregarded upon writing results. so when the flow is read back in, it's zero -- which is different than it would have been if it had stayed in memory. i wonder what that does to mass balance?

this version works, at least with Net3, simplenet, and SampleTown -- but the results are non-canonical, as the current workaround for the issue above writes tiny flow values to the binary file for closed elements. presumably we should check for closed status in the transport() function.

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:

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;
-
+      dv = D[n]*tstep;
-      // only adjust tank volume if we're in sequential mode.
+      Tank[i].V += dv;
-      // otherwise, the tankmixing function will do it for us.
-      if (!OpenQflag) {
-        dv = D[n]*tstep;
-        Tank[i].V += dv;
-      }
 
       /*** Updated 6/24/02 ***/
       /* Check if tank full/empty within next second */

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;
+     //if (S[i] <= CLOSED) x[i] = 0.0f;
-      else x[i] = (REAL4)Q[i];
+      /*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                         

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);
-        }
       }
    }
 }

src/testLemonTiger.cpp

@@ -1,6 +1,6 @@
 #include <map>
 #include <iomanip>
-
+#include <math.h>
 #include "testLemonTiger.h"
 #include "toolkit.h"
 
@@ -13,9 +13,19 @@ typedef struct {
   double head;
   double demand;
   double quality;
-} singleState_t;
+} nodeState_t;
 
-typedef map<int, singleState_t> networkState_t; // nodeIndex, state
+typedef struct {
+  double flow;
+} linkState_t;
+
+typedef map<int, nodeState_t> networkNodeState_t; // nodeIndex, state
+typedef map<int, linkState_t> networkLinkState_t; // linkIndex, state
+
+typedef struct {
+  networkNodeState_t nodeState;
+  networkLinkState_t linkState;
+} networkState_t;
 typedef map<long, networkState_t> 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;
+  int nNodes, nLinks;
-  float head, demand;
+  float head, demand, flow;
   ENgetcount(EN_NODECOUNT, &nNodes);
-  
+  ENgetcount(EN_LINKCOUNT, &nLinks);
-  for (int iNode = 1; iNode <= nNodes; iNode++) {
+  for (int iNode = 1; iNode <= nNodes; ++iNode) {
     ENgetnodevalue(iNode, EN_HEAD, &head);
     ENgetnodevalue(iNode, EN_DEMAND, &demand);
-    (*networkState)[iNode].head = head;
+    (*networkState).nodeState[iNode].head = head;
-    (*networkState)[iNode].demand = demand;
+    (*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;
+      networkNodeState_t::const_iterator networkNodeIterator;
-      for (networkIterator = state1.begin(); networkIterator != state1.end(); ++networkIterator) {
+      for (networkNodeIterator = nodeState1.begin(); networkNodeIterator != nodeState1.end(); ++networkNodeIterator) {
-        int nodeIndex = networkIterator->first;
+        int nodeIndex = networkNodeIterator->first;
         // trusting that all nodes are present...
-        const singleState_t nodeState1 = networkIterator->second;
+        const nodeState_t nodeState1 = networkNodeIterator->second;
-        const singleState_t nodeState2 = state2.at(nodeIndex);
+        const nodeState_t nodeState2 = networkNodeState2.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;
         
+        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;

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)
 
 /*

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                */

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