From ad48088eada1c952572a653609d191fd779cdd4f Mon Sep 17 00:00:00 2001 From: MehmetMelihKosucu Date: Thu, 15 Jul 2021 12:40:47 +0300 Subject: [PATCH] EPANET-RWCGGA - EPANET 3 hydraulic solver is enhanced with Rigid Water Column Global Gradient Algorithm i.e. unsteady compressible flow analysis algorithm of WDNs. - Valve representation which produces resistance according to valve opening is defined on EPANET 3. - In order to provide convergence globally Convergence Tracking Control Method is added to EPANET 3. - The developments and additions are performed on three WDNs varying complexity. --- CMakeLists.txt | 2 + src/CLI/main.cpp | 8 +- src/Core/constants.h | 2 +- src/Core/datamanager.cpp | 36 +- src/Core/datamanager.h | 3 +- src/Core/diagnostics.cpp | 2 +- src/Core/diagnostics.h | 2 +- src/Core/epanet3.cpp | 89 +++- src/Core/error.cpp | 2 +- src/Core/error.h | 2 +- src/Core/hydbalance.cpp | 57 +- src/Core/hydbalance.h | 14 +- src/Core/hydengine.cpp | 83 ++- src/Core/hydengine.h | 9 +- src/Core/network.cpp | 2 +- src/Core/network.h | 2 +- src/Core/options.cpp | 28 +- src/Core/options.h | 4 +- src/Core/project.cpp | 2 +- src/Core/project.h | 3 +- src/Core/qualbalance.cpp | 2 +- src/Core/qualbalance.h | 2 +- src/Core/qualengine.cpp | 2 +- src/Core/qualengine.h | 2 +- src/Core/units.cpp | 2 +- src/Core/units.h | 2 +- src/Elements/control.cpp | 2 +- src/Elements/control.h | 2 +- src/Elements/curve.cpp | 2 +- src/Elements/curve.h | 2 +- src/Elements/demand.cpp | 2 +- src/Elements/demand.h | 2 +- src/Elements/element.cpp | 2 +- src/Elements/element.h | 2 +- src/Elements/emitter.cpp | 2 +- src/Elements/emitter.h | 2 +- src/Elements/junction.cpp | 18 +- src/Elements/junction.h | 4 +- src/Elements/link.cpp | 9 +- src/Elements/link.h | 8 +- src/Elements/node.cpp | 10 +- src/Elements/node.h | 8 +- src/Elements/pattern.cpp | 2 +- src/Elements/pattern.h | 2 +- src/Elements/pipe.cpp | 12 +- src/Elements/pipe.h | 3 +- src/Elements/pump.cpp | 7 +- src/Elements/pump.h | 3 +- src/Elements/pumpcurve.cpp | 2 +- src/Elements/pumpcurve.h | 2 +- src/Elements/qualsource.cpp | 2 +- src/Elements/qualsource.h | 2 +- src/Elements/reservoir.cpp | 4 +- src/Elements/reservoir.h | 4 +- src/Elements/tank.cpp | 3 +- src/Elements/tank.h | 3 +- src/Elements/valve.cpp | 171 +++++- src/Elements/valve.h | 12 +- src/Input/controlparser.cpp | 2 +- src/Input/controlparser.h | 2 +- src/Input/curveparser.cpp | 2 +- src/Input/curveparser.h | 2 +- src/Input/inputparser.cpp | 2 +- src/Input/inputparser.h | 2 +- src/Input/inputreader.cpp | 2 +- src/Input/inputreader.h | 2 +- src/Input/linkparser.cpp | 130 ++--- src/Input/linkparser.h | 2 +- src/Input/nodeparser.cpp | 2 +- src/Input/nodeparser.h | 2 +- src/Input/optionparser.cpp | 8 +- src/Input/optionparser.h | 2 +- src/Input/patternparser.cpp | 2 +- src/Input/patternparser.h | 2 +- src/Models/demandmodel.cpp | 2 +- src/Models/demandmodel.h | 2 +- src/Models/headlossmodel.cpp | 2 +- src/Models/headlossmodel.h | 2 +- src/Models/leakagemodel.cpp | 2 +- src/Models/leakagemodel.h | 2 +- src/Models/pumpenergy.cpp | 2 +- src/Models/pumpenergy.h | 2 +- src/Models/qualmodel.cpp | 2 +- src/Models/qualmodel.h | 2 +- src/Models/tankmixmodel.cpp | 2 +- src/Models/tankmixmodel.h | 2 +- src/Output/outputfile.cpp | 2 +- src/Output/outputfile.h | 2 +- src/Output/projectwriter.cpp | 8 +- src/Output/projectwriter.h | 2 +- src/Output/reportfields.cpp | 2 +- src/Output/reportfields.h | 2 +- src/Output/reportwriter.cpp | 8 +- src/Output/reportwriter.h | 2 +- src/Solvers/ggasolver.cpp | 231 +++++++-- src/Solvers/ggasolver.h | 12 +- src/Solvers/hydsolver.cpp | 9 +- src/Solvers/hydsolver.h | 4 +- src/Solvers/ltdsolver.cpp | 2 +- src/Solvers/ltdsolver.h | 2 +- src/Solvers/matrixsolver.cpp | 2 +- src/Solvers/matrixsolver.h | 2 +- src/Solvers/qualsolver.cpp | 2 +- src/Solvers/qualsolver.h | 2 +- src/Solvers/rwcggasolver.cpp | 915 +++++++++++++++++++++++++++++++++ src/Solvers/rwcggasolver.h | 89 ++++ src/Solvers/sparspaksolver.cpp | 2 +- src/Solvers/sparspaksolver.h | 2 +- src/Utilities/graph.cpp | 2 +- src/Utilities/graph.h | 2 +- src/Utilities/mempool.cpp | 2 +- src/Utilities/mempool.h | 2 +- src/Utilities/segpool.cpp | 2 +- src/Utilities/segpool.h | 2 +- src/Utilities/utilities.cpp | 2 +- src/Utilities/utilities.h | 2 +- src/epanet3.h | 22 +- 117 files changed, 1917 insertions(+), 290 deletions(-) create mode 100644 src/Solvers/rwcggasolver.cpp create mode 100644 src/Solvers/rwcggasolver.h diff --git a/CMakeLists.txt b/CMakeLists.txt index 84fe302..de9897e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -64,6 +64,7 @@ src/Output/projectwriter.cpp src/Output/reportfields.cpp src/Output/reportwriter.cpp src/Solvers/ggasolver.cpp +src/Solvers/rwcggasolver.cpp src/Solvers/hydsolver.cpp src/Solvers/ltdsolver.cpp src/Solvers/matrixsolver.cpp @@ -124,6 +125,7 @@ src/Output/projectwriter.h src/Output/reportfields.h src/Output/reportwriter.h src/Solvers/ggasolver.h +src/Solvers/rwcggasolver.h src/Solvers/hydsolver.h src/Solvers/ltdsolver.h src/Solvers/matrixsolver.h diff --git a/src/CLI/main.cpp b/src/CLI/main.cpp index 99e8482..d674eef 100644 --- a/src/CLI/main.cpp +++ b/src/CLI/main.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). @@ -11,6 +11,9 @@ #include "epanet3.h" #include +#include + +//namespace plt = matplotlibcpp; int main(int argc, char* argv[]) { @@ -30,5 +33,8 @@ int main(int argc, char* argv[]) // ... run a full EPANET analysis EN_runEpanet(f1, f2, f3); //system("PAUSE"); + return 0; } + + diff --git a/src/Core/constants.h b/src/Core/constants.h index 212b6ef..92e0e81 100644 --- a/src/Core/constants.h +++ b/src/Core/constants.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Core/datamanager.cpp b/src/Core/datamanager.cpp index d2a2051..017849d 100644 --- a/src/Core/datamanager.cpp +++ b/src/Core/datamanager.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -272,6 +272,40 @@ int DataManager::getLinkValue(int index, int param, double* value, Network* nw) //----------------------------------------------------------------------------- +int DataManager::setLinkValue(int index, int param, double value, Network* nw) +{ + if (index < 0 || index >= nw->count(Element::LINK)) return 205; + Link* link = nw->link(index); + switch (param) + { + case EN_DIAMETER: + link->diameter = value / nw->ucf(Units::DIAMETER); + link->setLossFactor(); + link->setResistance(nw); + break; + case EN_MINORLOSS: + link->lossCoeff = value; + link->setLossFactor(); + break; + case EN_INITSTATUS: + link->initStatus = value; + break; + case EN_INITSETTING: + link->initSetting = value; + break; + case EN_FLOW: + link->flow = value / nw->ucf(Units::FLOW); + break; + case EN_STATUS: + link->status = value; + break; + case EN_ENERGY: + break; // TO BE ADDED + } + return 0; +} + +//----------------------------------------------------------------------------- int getTankValue(int param, Node* node, double* value, Network* nw) { double lcf = nw->ucf(Units::LENGTH); diff --git a/src/Core/datamanager.h b/src/Core/datamanager.h index d2e07bc..f280f97 100644 --- a/src/Core/datamanager.h +++ b/src/Core/datamanager.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -28,6 +28,7 @@ struct DataManager static int getLinkType(int index, int* type, Network* nw); static int getLinkNodes(int index, int* fromNode, int* toNode, Network* nw); static int getLinkValue(int index, int param, double* value, Network* nw); + static int setLinkValue(int index, int param, double v, Network* nw); }; #endif // DATAMANAGER_H_ diff --git a/src/Core/diagnostics.cpp b/src/Core/diagnostics.cpp index f27bc8d..b25747d 100644 --- a/src/Core/diagnostics.cpp +++ b/src/Core/diagnostics.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Core/diagnostics.h b/src/Core/diagnostics.h index ec98469..32158f7 100644 --- a/src/Core/diagnostics.h +++ b/src/Core/diagnostics.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Core/epanet3.cpp b/src/Core/epanet3.cpp index 8776032..b1f76d2 100644 --- a/src/Core/epanet3.cpp +++ b/src/Core/epanet3.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). @@ -6,7 +6,7 @@ */ /////////////////////////////////////////////// -// Implementation of EPANET 3's API library // +// Implementation of EPANET 3.1's API library // ///////////////'''''/////////////////////////// // TO DO: @@ -26,6 +26,7 @@ #include using namespace Epanet; +using namespace std; #define project(p) ((Project *)p) @@ -49,9 +50,26 @@ int EN_runEpanet(const char* inpFile, const char* rptFile, const char* outFile) Project p; int err = 0; + // Output files in text format, and variables which are existing in output file + + //std::ofstream myfile("D:\\EPANET_3\\Networks\\RWC\\Onizuka\\Onizuka1986CCV.txt"); + //std::ofstream myfile("D:\\EPANET_3\\Networks\\RWC\\Nault2016\\Nault2016CCV.txt"); + std::ofstream myfile("D:\\EPANET_3\\Networks\\RWC\\Nault2018\\Nault2018.txt"); + //std::ofstream myfile("D:\\EPANET_3\\Networks\\RWC\\Nault2018EPA3-EPS.txt"); + + int IndexP5, IndexP8, IndexT1, IndexJ2; + double f5, f8, h1, h2; //Onizuka (1986) + + int IndexJ7, IndexJ18, IndexJ25, IndexP34, IndexPS2P3, IndexP2, IndexP25, IndexR1, IndexR2, IndexR3; + double h7, h18, h25, f34, fPS2P3, f2, f25; // Nault and Karney (2016) + + int IndexJ70, IndexJ105, IndexP447, IndexP1364; + double h70, h105, f447, f1364; // Nault and Karney(2018) + // ... initialize execution time clock clock_t start_t = clock(); + for (;;) { // ... open the command line files and load network data @@ -66,6 +84,42 @@ int EN_runEpanet(const char* inpFile, const char* rptFile, const char* outFile) if ( (err = p.initSolver(false)) ) break; std::cout << "\n "; + /*int ErrorP5 = EN_getLinkIndex("P5", &IndexP5, p.getNetwork()); + int ErrorP8 = EN_getLinkIndex("P8", &IndexP8, p.getNetwork()); + int ErrorT1 = EN_getNodeIndex("T1", &IndexT1, p.getNetwork()); + int ErrorJ2 = EN_getNodeIndex("J2", &IndexJ2, p.getNetwork()); + + double ErrorValP5 = EN_getLinkValue(IndexP5, EN_FLOW, &f5, p.getNetwork()); + double ErrorValP8 = EN_getLinkValue(IndexP8, EN_FLOW, &f8, p.getNetwork()); + double ErrorValT1 = EN_getNodeValue(IndexT1, EN_HEAD, &h1, p.getNetwork()); + double ErrorValJ2 = EN_getNodeValue(IndexJ2, EN_HEAD, &h2, p.getNetwork()); // Onizuka (1986) */ + + /*int ErrorJ7 = EN_getNodeIndex("J7", &IndexJ7, p.getNetwork()); + int ErrorJ18 = EN_getNodeIndex("J18", &IndexJ18, p.getNetwork()); + int ErrorJ25 = EN_getNodeIndex("J25", &IndexJ25, p.getNetwork()); + int ErrorP34 = EN_getLinkIndex("P34", &IndexP34, p.getNetwork()); + int ErrorPS2P3 = EN_getLinkIndex("PS2-P3", &IndexPS2P3, p.getNetwork()); + int ErrorP2 = EN_getLinkIndex("P2", &IndexP2, p.getNetwork()); + int ErrorP25 = EN_getLinkIndex("P25", &IndexP25, p.getNetwork()); + + double ErrorValJ7 = EN_getNodeValue(IndexJ7, EN_HEAD, &h7, p.getNetwork()); + double ErrorValJ18 = EN_getNodeValue(IndexJ18, EN_HEAD, &h18, p.getNetwork()); + double ErrorValJ25 = EN_getNodeValue(IndexJ25, EN_HEAD, &h25, p.getNetwork()); + double ErrorValP34 = EN_getLinkValue(IndexP34, EN_FLOW, &f34, p.getNetwork()); + double ErrorValPS2P3 = EN_getLinkValue(IndexPS2P3, EN_FLOW, &fPS2P3, p.getNetwork()); + double ErrorValP2 = EN_getLinkValue(IndexP2, EN_FLOW, &f2, p.getNetwork()); + double ErrorValP25 = EN_getLinkValue(IndexP25, EN_FLOW, &f25, p.getNetwork()); // Nault and Karney (2016) */ + + int ErrorJ70 = EN_getNodeIndex("J70", &IndexJ70, p.getNetwork()); + int ErrorJ105 = EN_getNodeIndex("J105", &IndexJ105, p.getNetwork()); + int ErrorP447 = EN_getLinkIndex("P447", &IndexP447, p.getNetwork()); + int ErrorP1364 = EN_getLinkIndex("P1364", &IndexP1364, p.getNetwork()); + + double ErrorValJ70 = EN_getNodeValue(IndexJ70, EN_HEAD, &h70, p.getNetwork()); + double ErrorValJ105 = EN_getNodeValue(IndexJ105, EN_HEAD, &h105, p.getNetwork()); + double ErrorValP447 = EN_getLinkValue(IndexP447, EN_FLOW, &f447, p.getNetwork()); + double ErrorValP1364 = EN_getLinkValue(IndexP1364, EN_FLOW, &f1364, p.getNetwork()); // Nault and Karney (2018) */ + // ... step through each time period int t = 0; int tstep = 0; @@ -73,13 +127,37 @@ int EN_runEpanet(const char* inpFile, const char* rptFile, const char* outFile) { std::cout << "\r Solving network at " //r << Utilities::getTime(t+tstep) << " hrs ... "; - + // ... run solver to compute hydraulics err = p.runSolver(&t); p.writeMsgLog(); // ... advance solver to next period in time while solving for water quality if ( !err ) err = p.advanceSolver(&tstep); + + /*double ErrorValP5 = EN_getLinkValue(IndexP5, EN_FLOW, &f5, p.getNetwork()); + double ErrorValP8 = EN_getLinkValue(IndexP8, EN_FLOW, &f8, p.getNetwork()); + double ErrorValT1 = EN_getNodeValue(IndexT1, EN_HEAD, &h1, p.getNetwork()); + double ErrorValJ2 = EN_getNodeValue(IndexJ2, EN_HEAD, &h2, p.getNetwork()); + + myfile << t << " " << f5 << " " << f8 << " " << h1 << " " << h2 << "\n"; // */ + + /*double ErrorValJ7 = EN_getNodeValue(IndexJ7, EN_HEAD, &h7, p.getNetwork()); + double ErrorValJ18 = EN_getNodeValue(IndexJ18, EN_HEAD, &h18, p.getNetwork()); + double ErrorValJ25 = EN_getNodeValue(IndexJ25, EN_HEAD, &h25, p.getNetwork()); + double ErrorValP34 = EN_getLinkValue(IndexP34, EN_FLOW, &f34, p.getNetwork()); + double ErrorValPS2P3 = EN_getLinkValue(IndexPS2P3, EN_FLOW, &fPS2P3, p.getNetwork()); + double ErrorValP2 = EN_getLinkValue(IndexP2, EN_FLOW, &f2, p.getNetwork()); + double ErrorValP25 = EN_getLinkValue(IndexP25, EN_FLOW, &f25, p.getNetwork()); + myfile << t << " " << h7 << " " << h18 << " " << h25 << " " << f34 << " " << fPS2P3 << " " << f2 << " " << f25 << "\n"; // */ + + double ErrorValJ70 = EN_getNodeValue(IndexJ70, EN_HEAD, &h70, p.getNetwork()); + double ErrorValJ105 = EN_getNodeValue(IndexJ105, EN_HEAD, &h105, p.getNetwork()); + double ErrorValP447 = EN_getLinkValue(IndexP447, EN_FLOW, &f447, p.getNetwork()); + double ErrorValP1364 = EN_getLinkValue(IndexP1364, EN_FLOW, &f1364, p.getNetwork()); + + myfile << t << " " << h70 << " " << h105 << " " << f447 << " " << f1364 << "\n"; // */ + } while (tstep > 0 && !err ); break; } @@ -339,5 +417,10 @@ int EN_getLinkValue(int index, int param, double* value, EN_Project p) return DataManager::getLinkValue(index, param, value, project(p)->getNetwork()); } +int EN_setLinkValue(int index, int param, double value, EN_Project p) +{ + return DataManager::setLinkValue(index, param, value, project(p)->setNetwork()); +} + } // end of namespace diff --git a/src/Core/error.cpp b/src/Core/error.cpp index 7541eb9..0a4886c 100644 --- a/src/Core/error.cpp +++ b/src/Core/error.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Core/error.h b/src/Core/error.h index c9348d6..f10cf3e 100644 --- a/src/Core/error.h +++ b/src/Core/error.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Core/hydbalance.cpp b/src/Core/hydbalance.cpp index 01ab46c..e86476a 100644 --- a/src/Core/hydbalance.cpp +++ b/src/Core/hydbalance.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). @@ -16,9 +16,14 @@ #include "network.h" #include "Elements/node.h" #include "Elements/link.h" +#include "Elements/valve.h" +#include "Elements/node.h" +#include "Core/hydengine.h" +#include "Solvers/rwcggasolver.h" #include #include +#include using namespace std; void findNodeOutflows(double lamda, double dH[], double xQ[], Network* nw); @@ -36,7 +41,9 @@ double HydBalance::evaluate( double dH[], // change in nodal heads double dQ[], // change in link flows double xQ[], // nodal inflow minus outflow - Network* nw) // network being analyzed + Network* nw, // network being analyzed + int currentTime, + double tstep) { // ... initialize which elements have the maximum errors maxFlowErr = 0.0; @@ -54,7 +61,7 @@ double HydBalance::evaluate( // ... find the error norm in satisfying conservation of energy // (updating xQ with internal link flows) - double norm = findHeadErrorNorm(lamda, dH, dQ, xQ, nw); + double norm = findHeadErrorNorm(lamda, dH, dQ, xQ, nw, currentTime, tstep); // ... update xQ with external outflows @@ -78,7 +85,7 @@ double HydBalance::evaluate( // Find the error norm in satisfying the head loss equation across each link. double HydBalance::findHeadErrorNorm( - double lamda, double dH[], double dQ[], double xQ[], Network* nw) + double lamda, double dH[], double dQ[], double xQ[], Network* nw, int currentTime, double tstep) { double norm = 0.0; double count = 0.0; @@ -89,12 +96,18 @@ double HydBalance::findHeadErrorNorm( int linkCount = nw->count(Element::LINK); for (int i = 0; i < linkCount; i++) { + // ... identify link's end nodes Link* link = nw->link(i); int n1 = link->fromNode->index; int n2 = link->toNode->index; + if (link->status == Link::LINK_CLOSED) + { + int mmk = 32; + } + // ... apply updated flow to end node flow balances double flowChange = lamda * dQ[i]; @@ -104,26 +117,43 @@ double HydBalance::findHeadErrorNorm( // ... update network's max. flow change + previousMaxFlowChange = maxFlowChange; + double err = abs(flowChange); if ( err > maxFlowChange ) { maxFlowChange = err; maxFlowChangeLink = i; } - + // ... compute head loss and its gradient (head loss is saved // ... to link->hLoss and its gradient to link->hGrad) //******************************************************************* link->findHeadLoss(nw, flow); //******************************************************************* + + double unsteadyTerm = 0; - // ... evaluate head loss error + // ... evaluate head loss error according to Steady and Unsteady Flow Conditions + - double h1 = link->fromNode->head + lamda * dH[n1]; - double h2 = link->toNode->head + lamda * dH[n2]; + if (currentTime == 0 || nw->option(Options::HYD_SOLVER) == "GGA" ) + { + unsteadyTerm = 0; + } + + else + { + unsteadyTerm = (link->inertialTerm) * (link->flow - link->pastFlow) / tstep; + } + + h1 = link->fromNode->head + lamda * dH[n1]; + h2 = link->toNode->head + lamda * dH[n2]; if ( link->hGrad == 0.0 ) link->hLoss = h1 - h2; - err = h1 - h2 - link->hLoss; - if ( abs(err) > maxHeadErr ) + //err = h1 - h2 - link->hLoss; + err = unsteadyTerm - h1 + h2 + link->hLoss; + + if ( abs(err) > maxHeadErr ) { maxHeadErr = abs(err); maxHeadErrLink = i; @@ -319,7 +349,14 @@ double findTotalFlowChange(double lamda, double dQ[], Network* nw) for ( int i = 0; i < nw->count(Element::LINK); i++ ) { + Link* link = nw->link(i); + + if (link->status == Link::LINK_CLOSED) + { + int mmk = 32; + } + dq = lamda * dQ[i]; dqSum += abs(dq); qSum += abs(link->flow + dq); diff --git a/src/Core/hydbalance.h b/src/Core/hydbalance.h index 3bda602..09b6cb1 100644 --- a/src/Core/hydbalance.h +++ b/src/Core/hydbalance.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -24,7 +24,15 @@ struct HydBalance { double maxFlowErr; //!< max. flow error (cfs) double maxHeadErr; //!< max. head loss error (ft) + double h1ini; + double h1; + double h2ini; + double h2; + double phloss; + double gHn; + double previousMaxHeadErr; //!< previous max. head loss error (ft) double maxFlowChange; //!< max. flow change (cfs) + double previousMaxFlowChange; //!< previous max. flow change (cfs) double totalFlowChange; //!< (summed flow changes) / (summed flows) int maxHeadErrLink; //!< link with max. head loss error @@ -32,9 +40,9 @@ struct HydBalance int maxFlowChangeLink; //!< link with max. flow change double evaluate( - double lamda, double dH[], double dQ[], double xQ[], Network* nw); + double lamda, double dH[], double dQ[], double xQ[], Network* nw, int currentTime, double tstep); double findHeadErrorNorm( - double lamda, double dH[], double dQ[], double xQ[], Network* nw); + double lamda, double dH[], double dQ[], double xQ[], Network* nw, int currentTime, double tstep); double findFlowErrorNorm(double xQ[], Network* nw); }; diff --git a/src/Core/hydengine.cpp b/src/Core/hydengine.cpp index 2320dac..7748df6 100644 --- a/src/Core/hydengine.cpp +++ b/src/Core/hydengine.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). @@ -173,12 +173,18 @@ int HydEngine::solve(int* t) //if ( network->option(Options::REPORT_TRIALS) ) network->msgLog << endl; int trials = 0; - int statusCode = hydSolver->solve(hydStep, trials); + int statusCode = hydSolver->solve(hydStep, trials, currentTime); if ( statusCode == HydSolver::SUCCESSFUL && isPressureDeficient() ) { statusCode = resolvePressureDeficiency(trials); } + + for (Link* link : network->links) + { + link->previousStatus = link->status; + } + reportDiagnostics(statusCode, trials); if ( halted ) throw SystemError(SystemError::HYDRAULICS_SOLVER_FAILURE); return statusCode; @@ -213,6 +219,9 @@ void HydEngine::advance(int* tstep) updateEnergyUsage(); updateTanks(); + pastJunction(); + pastLink(); + // ... advance time counters currentTime += hydStep; @@ -289,7 +298,7 @@ void HydEngine::updateCurrentConditions() int p = network->option(Options::DEMAND_PATTERN); if ( p >= 0 ) patternFactor = network->pattern(p)->currentFactor(); - + // ... update node conditions for (Node* node : network->nodes) @@ -299,6 +308,26 @@ void HydEngine::updateCurrentConditions() // ... set its fixed grade state (for tanks & reservoirs) node->setFixedGrade(); + + if (node->type() == Node::JUNCTION) + { + double ph; + node->pastHead = node->head; + node->ph = node->head; + + } + else if (node->type() == Node::RESERVOIR) + { + node->pastHead = node->head; + node->ph = node->head; + } + + else if (node->type() == Node::TANK) + { + node->pastHead = node->head; + node->ph = node->head; + } + } // ... update link conditions @@ -308,6 +337,13 @@ void HydEngine::updateCurrentConditions() // ... open a temporarily closed link //if ( link->status >= Link::TEMP_CLOSED ) link->status = Link::LINK_OPEN; + if (link->type() == Link::PIPE || Link::PUMP || Link::VALVE) + { + link->pastFlow = link->flow; + link->pastHloss = link->hLoss; + link->pastSetting = link->setting; + } + // ... apply pattern-based pump or valve setting link->applyControlPattern(network->msgLog); } @@ -352,7 +388,7 @@ int HydEngine::resolvePressureDeficiency(int& trials) // set to fixed grade (which occurred in isPressureDeficient()) if ( reportTrials ) network->msgLog << s_ReSolve1; - int statusCode = hydSolver->solve(hydStep, trials2); + int statusCode = hydSolver->solve(hydStep, trials2, currentTime); if ( statusCode == HydSolver::FAILED_ILL_CONDITIONED ) return statusCode; // ... adjust actual demands for the pressure deficient junctions @@ -375,7 +411,7 @@ int HydEngine::resolvePressureDeficiency(int& trials) network->msgLog << "\n\n " << count1 << s_Reductions1; network->msgLog << s_ReSolve2; } - statusCode = hydSolver->solve(hydStep, trials3); + statusCode = hydSolver->solve(hydStep, trials3, currentTime); // ... check once more for any remaining pressure deficiencies @@ -401,7 +437,7 @@ int HydEngine::resolvePressureDeficiency(int& trials) network->msgLog << "\n " << count2 << s_Reductions2; network->msgLog << s_ReSolve2 << "\n"; } - statusCode = hydSolver->solve(hydStep, trials4); + statusCode = hydSolver->solve(hydStep, trials4, currentTime); } trials += trials2 + trials3 + trials4; @@ -592,7 +628,9 @@ void HydEngine::updateTanks() { Tank* tank = static_cast(node); tank->pastHead = tank->head; + tank->ph = tank->head; tank->pastVolume = tank->volume; + tank->pastArea = tank->area; tank->pastOutflow = tank->outflow; node->fixedGrade = true; tank->updateVolume(hydStep); @@ -601,6 +639,39 @@ void HydEngine::updateTanks() } } +void HydEngine::pastJunction() +{ + for (Node* node : network->nodes) + { + if (node->type() == Node::JUNCTION) + { + double ph; + node->pastHead = node->head; + node->ph = node->head; + + } + else if (node->type() == Node::RESERVOIR) + { + node->pastHead = node->head; + node->ph = node->head; + } + } +} + + +void HydEngine::pastLink() +{ + for (Link* link : network->links) + { + if (link->type() == Link::PIPE || link->type() == Link::PUMP || link->type() == Link::VALVE) + { + link->pastFlow = link->flow; + link->pastHloss = link->hLoss; + link->pastSetting = link->setting; + } + } +} + //----------------------------------------------------------------------------- // Advances all time patterns. diff --git a/src/Core/hydengine.h b/src/Core/hydengine.h index 236f676..2046a5b 100644 --- a/src/Core/hydengine.h +++ b/src/Core/hydengine.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). @@ -43,6 +43,7 @@ class HydEngine int getElapsedTime() { return currentTime; } double getPeakKwatts() { return peakKwatts; } + int currentTime; //!< current simulation time (sec) private: @@ -54,7 +55,7 @@ class HydEngine // Engine components Network* network; //!< network being analyzed - HydSolver* hydSolver; //!< steady state hydraulic solver + HydSolver* hydSolver; //!< steady state or rwc unsteady hydraulic solver MatrixSolver* matrixSolver; //!< sparse matrix solver // HydFile* hydFile; //!< hydraulics file accessor @@ -65,7 +66,7 @@ class HydEngine int startTime; //!< starting time of day (sec) int rptTime; //!< current reporting time (sec) int hydStep; //!< hydraulic time step (sec) - int currentTime; //!< current simulation time (sec) + int timeOfDay; //!< current time of day (sec) double peakKwatts; //!< peak energy usage (kwatts) std::string timeStepReason; //!< reason for taking next time step @@ -75,6 +76,8 @@ class HydEngine void initMatrixSolver(); int getTimeStep(); + void pastJunction(); + void pastLink(); int timeToPatternChange(int tstep); int timeToActivateControl(int tstep); int timeToCloseTank(int tstep); diff --git a/src/Core/network.cpp b/src/Core/network.cpp index 5606110..54eed6a 100644 --- a/src/Core/network.cpp +++ b/src/Core/network.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Core/network.h b/src/Core/network.h index 8009af5..c1f29d8 100644 --- a/src/Core/network.h +++ b/src/Core/network.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Core/options.cpp b/src/Core/options.cpp index f39f0a1..a473797 100644 --- a/src/Core/options.cpp +++ b/src/Core/options.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -33,11 +33,17 @@ static const char* flowUnitsWords[] = // Keywords for PressureUnits enumeration in options.h static const char* pressureUnitsWords[] = {"PSI", "METERS", "PKA", 0}; +// Keywords for Hyd_Solver enumeration in options.h +static const char* hydSolverWords[] = { "GGA", "RWCGGA", 0 }; + // Headloss formula keywords static const char* headlossModelWords[] = {"H-W", "D-W", "C-M", 0}; // Hydraulic Newton solver step size method names -static const char* stepSizingWords[] = {"FULL", "RELAXATION", "LINESEARCH", 0}; +static const char* stepSizingWords[] = {"FULL", "RELAXATION", "LINESEARCH", "BRF", "ARF", 0}; + +// Valve representation types names +static const char* valveRepWords[] = { "Toe", "Cd", 0 }; static const char* ifUnbalancedWords[] = {"STOP", "CONTINUE", 0}; @@ -80,6 +86,7 @@ void Options::setDefaults() stringOptions[LEAKAGE_MODEL] = "NONE"; stringOptions[HYD_SOLVER] = "GGA"; stringOptions[STEP_SIZING] = "FULL"; + stringOptions[VALVE_REP_TYPE] = "Cd"; stringOptions[MATRIX_SOLVER] = "SPARSPAK"; stringOptions[DEMAND_PATTERN_NAME] = ""; stringOptions[QUAL_MODEL] = "NONE"; @@ -109,7 +116,7 @@ void Options::setDefaults() valueOptions[MINIMUM_PRESSURE] = 0.0; valueOptions[SERVICE_PRESSURE] = 0.0; valueOptions[PRESSURE_EXPONENT] = 0.5; - valueOptions[EMITTER_EXPONENT] = 0.5; + valueOptions[EMITTER_EXPONENT] = 1.18; valueOptions[DEMAND_MULTIPLIER] = 1.0; valueOptions[RELATIVE_ACCURACY] = 0.0; @@ -117,6 +124,7 @@ void Options::setDefaults() valueOptions[FLOW_TOLERANCE] = 0.0; valueOptions[FLOW_CHANGE_LIMIT] = 0.0; valueOptions[TIME_WEIGHT] = 0.0; + valueOptions[TEMP_DISC_PARA] = 0.0; valueOptions[ENERGY_PRICE] = 0.0; valueOptions[PEAKING_CHARGE] = 0.0; @@ -162,12 +170,24 @@ int Options::setOption(StringOption option, const string& value) stringOptions[HEADLOSS_MODEL] = headlossModelWords[i]; break; + case HYD_SOLVER: + i = Utilities::findFullMatch(value, hydSolverWords); + if (i < 0) return InputError::INVALID_KEYWORD; + stringOptions[HYD_SOLVER] = hydSolverWords[i]; + break; + case STEP_SIZING: i = Utilities::findFullMatch(value, stepSizingWords); if (i < 0) return InputError::INVALID_KEYWORD; stringOptions[STEP_SIZING] = stepSizingWords[i]; break; + case VALVE_REP_TYPE: + i = Utilities::findFullMatch(value, valveRepWords); + if (i < 0) return InputError::INVALID_KEYWORD; + stringOptions[VALVE_REP_TYPE] = valveRepWords[i]; + break; + case DEMAND_MODEL: i = Utilities::findFullMatch(value, demandModelWords); if (i < 0) return InputError::INVALID_KEYWORD; @@ -364,6 +384,8 @@ string Options::hydOptionsToStr() } s << setw(w) << "TIME_WEIGHT"; s << valueOptions[TIME_WEIGHT] << "\n"; + s << setw(w) << "TEMP_DISC_PARA"; + s << valueOptions[TEMP_DISC_PARA] << "\n"; s << setw(w) << "STEP_SIZING"; s << stringOptions[STEP_SIZING] << "\n"; s << setw(w) << "IF_UNBALANCED"; diff --git a/src/Core/options.h b/src/Core/options.h index 8b1edfd..e68293a 100644 --- a/src/Core/options.h +++ b/src/Core/options.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). @@ -48,6 +48,7 @@ class Options LEAKAGE_MODEL, //!< Name of pipe leakage model used HYD_SOLVER, //!< Name of hydraulic solver method STEP_SIZING, //!< Name of Newton step size method + VALVE_REP_TYPE, //!< Name of Valve representation type MATRIX_SOLVER, //!< Name of sparse matrix eqn. solver DEMAND_PATTERN_NAME, //!< Name of global demand pattern @@ -106,6 +107,7 @@ class Options FLOW_TOLERANCE, //!< Convergence tolerance for flow balance FLOW_CHANGE_LIMIT, //!< Max. flow change for convergence TIME_WEIGHT, //!< Time weighting for variable head tanks + TEMP_DISC_PARA, //!< Temporal Discretization Parameter // Water quality options MOLEC_DIFFUSIVITY, //!< Chemical's molecular diffusivity (ft2/sec) diff --git a/src/Core/project.cpp b/src/Core/project.cpp index 81dc3a2..dbf8ca8 100644 --- a/src/Core/project.cpp +++ b/src/Core/project.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). diff --git a/src/Core/project.h b/src/Core/project.h index f7c3e9f..20e01a4 100644 --- a/src/Core/project.h +++ b/src/Core/project.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -67,6 +67,7 @@ namespace Epanet void writeMsgLog(std::ostream& out); void writeMsgLog(); Network* getNetwork() { return &network; } + Network* setNetwork() { return &network; } private: diff --git a/src/Core/qualbalance.cpp b/src/Core/qualbalance.cpp index f20073f..8627807 100644 --- a/src/Core/qualbalance.cpp +++ b/src/Core/qualbalance.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). diff --git a/src/Core/qualbalance.h b/src/Core/qualbalance.h index 4468428..367d48d 100644 --- a/src/Core/qualbalance.h +++ b/src/Core/qualbalance.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Core/qualengine.cpp b/src/Core/qualengine.cpp index fb591d8..1e04d57 100644 --- a/src/Core/qualengine.cpp +++ b/src/Core/qualengine.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Core/qualengine.h b/src/Core/qualengine.h index f5b6a17..30dc25e 100644 --- a/src/Core/qualengine.h +++ b/src/Core/qualengine.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Core/units.cpp b/src/Core/units.cpp index 2ae99ef..d8d21f5 100644 --- a/src/Core/units.cpp +++ b/src/Core/units.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Core/units.h b/src/Core/units.h index 26b1f1a..af31799 100644 --- a/src/Core/units.h +++ b/src/Core/units.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Elements/control.cpp b/src/Elements/control.cpp index eae7cbe..b02bc5b 100644 --- a/src/Elements/control.cpp +++ b/src/Elements/control.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Elements/control.h b/src/Elements/control.h index 97a995b..6c7fca7 100644 --- a/src/Elements/control.h +++ b/src/Elements/control.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Elements/curve.cpp b/src/Elements/curve.cpp index 9ad6fd1..8c8e156 100644 --- a/src/Elements/curve.cpp +++ b/src/Elements/curve.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Elements/curve.h b/src/Elements/curve.h index 4da4ec6..4024721 100644 --- a/src/Elements/curve.h +++ b/src/Elements/curve.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Elements/demand.cpp b/src/Elements/demand.cpp index 9d9cbe7..a77b6eb 100644 --- a/src/Elements/demand.cpp +++ b/src/Elements/demand.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Elements/demand.h b/src/Elements/demand.h index 4d7dd5a..8a52c62 100644 --- a/src/Elements/demand.h +++ b/src/Elements/demand.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Elements/element.cpp b/src/Elements/element.cpp index 94c3330..30b91c0 100644 --- a/src/Elements/element.cpp +++ b/src/Elements/element.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Elements/element.h b/src/Elements/element.h index 97adc9e..2cace80 100644 --- a/src/Elements/element.h +++ b/src/Elements/element.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Elements/emitter.cpp b/src/Elements/emitter.cpp index 68c3732..dc933e5 100644 --- a/src/Elements/emitter.cpp +++ b/src/Elements/emitter.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). diff --git a/src/Elements/emitter.h b/src/Elements/emitter.h index bfc6012..c7c6de8 100644 --- a/src/Elements/emitter.h +++ b/src/Elements/emitter.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). diff --git a/src/Elements/junction.cpp b/src/Elements/junction.cpp index 15ba55f..c0ca6a1 100644 --- a/src/Elements/junction.cpp +++ b/src/Elements/junction.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). @@ -16,15 +16,19 @@ using namespace std; //----------------------------------------------------------------------------- // Junction Constructor //----------------------------------------------------------------------------- -Junction::Junction(string name_): - Node(name_), - pMin(MISSING), - pFull(MISSING), - emitter(nullptr) +Junction::Junction(string name_) : +Node(name_), +pMin(MISSING), +pFull(MISSING), +emitter(nullptr) +//pastHead(0.0), +//ph(0.0) { + } + //----------------------------------------------------------------------------- // Junction Destructor //----------------------------------------------------------------------------- @@ -81,6 +85,8 @@ void Junction::convertUnits(Network* nw) void Junction::initialize(Network* nw) { head = elev + (pFull - pMin) / 2.0;; + //pastHead = elev + (pFull - pMin) / 2.0; // Head value on the previous timestep + ph = elev + (pFull - pMin) / 2.0; // synonym of past head quality = initQual; actualDemand = 0.0; outflow = 0.0; diff --git a/src/Elements/junction.h b/src/Elements/junction.h index 2555090..6df3b2f 100644 --- a/src/Elements/junction.h +++ b/src/Elements/junction.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). @@ -44,6 +44,8 @@ class Junction: public Node double pMin; //!< minimum pressure head to have demand (ft) double pFull; //!< pressure head required for full demand (ft) Emitter* emitter; //!< emitter object + double pastHead; //!< Head on the previous time step + double ph; //!< synonym of past head }; #endif diff --git a/src/Elements/link.cpp b/src/Elements/link.cpp index 46ec7e5..91df7df 100644 --- a/src/Elements/link.cpp +++ b/src/Elements/link.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -37,12 +37,16 @@ Link::Link(string name_) : lossCoeff(0.0), initSetting(1.0), status(0), + previousStatus(0), flow(0.0), + pastFlow(0.0), leakage(0.0), hLoss(0.0), + pastHloss(0.0), hGrad(0.0), setting(0.0), - quality(0.0) + quality(0.0), + inertialTerm(0.0) {} /// Destructor @@ -84,6 +88,7 @@ void Link::initialize(bool reInitFlow) else setInitFlow(); } leakage = 0.0; + inertialTerm = 0; } //----------------------------------------------------------------------------- diff --git a/src/Elements/link.h b/src/Elements/link.h index 6d51ea6..652bd10 100644 --- a/src/Elements/link.h +++ b/src/Elements/link.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -53,6 +53,7 @@ class Link: public Element virtual void setInitStatus(int s) {} virtual void setInitSetting(double s) {} virtual void setResistance(Network* nw) {} + virtual void setLossFactor() {} // Retrieves hydraulic variables virtual double getVelocity() {return 0.0;} @@ -103,12 +104,17 @@ class Link: public Element // Computed Variables int status; //!< current status + int previousStatus; //!< Status on the previous time step double flow; //!< flow rate (cfs) + double pastFlow; // Yeni double leakage; //!< leakage rate (cfs) double hLoss; //!< head loss (ft) + double pastHloss; // Yeni double hGrad; //!< head loss gradient (ft/cfs) double setting; //!< current setting + double pastSetting; double quality; //!< avg. quality concen. (mass/ft3) + double inertialTerm; // Yeni }; #endif diff --git a/src/Elements/node.cpp b/src/Elements/node.cpp index dc0e6d2..e264cbb 100644 --- a/src/Elements/node.cpp +++ b/src/Elements/node.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). @@ -28,6 +28,10 @@ Node::Node(string name_) : qualSource(nullptr), fixedGrade(false), head(0.0), + h1ini(0.0), + h2ini(0.0), + pastHead(0.0), + ph(0.0), qGrad(0.0), fullDemand(0.0), actualDemand(0.0), @@ -69,6 +73,10 @@ Node* Node::factory(int type_, string name_, MemPool* memPool) void Node::initialize(Network* nw) { head = elev; + pastHead = elev; // head on the previous timestep + ph = elev; // synonym of the past head + h1ini = elev; + h2ini = elev; quality = initQual; if ( qualSource ) qualSource->quality = quality; actualDemand = 0.0; diff --git a/src/Elements/node.h b/src/Elements/node.h index 751f72d..3c86366 100644 --- a/src/Elements/node.h +++ b/src/Elements/node.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). @@ -69,6 +69,12 @@ class Node: public Element // Computed Variables bool fixedGrade; //!< fixed grade status double head; //!< hydraulic head (ft) + double h1ini; //!< hydraulic head of upstream node in the initial moment of iteration + double h2ini; //!< hydraulic head of downstream node in the initial moment of iteration + //double h1; + //double h2; + double pastHead; //!< Hydraulic Head in previous timestep + double ph; //!< synonym of past head double qGrad; //!< gradient of outflow w.r.t. head (cfs/ft) double fullDemand; //!< full demand required (cfs) double actualDemand; //!< actual demand delivered (cfs) diff --git a/src/Elements/pattern.cpp b/src/Elements/pattern.cpp index c4bdec6..5d95705 100644 --- a/src/Elements/pattern.cpp +++ b/src/Elements/pattern.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Elements/pattern.h b/src/Elements/pattern.h index 436218e..6f9bbe7 100644 --- a/src/Elements/pattern.h +++ b/src/Elements/pattern.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Elements/pipe.cpp b/src/Elements/pipe.cpp index f7fbc92..420397f 100644 --- a/src/Elements/pipe.cpp +++ b/src/Elements/pipe.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -99,10 +99,17 @@ void Pipe::setInitFlow() { // ... flow at velocity of 1 ft/s flow = PI * diameter * diameter / 4.0; + pastFlow = 0; + pastHloss = 0; } //----------------------------------------------------------------------------- +void Pipe::setLossFactor() +{ + lossFactor = 0.02517 * lossCoeff / pow(diameter, 4); +} + double Pipe::getVelocity() { double area = PI * diameter * diameter / 4.0; @@ -131,14 +138,17 @@ void Pipe::findHeadLoss(Network* nw, double q) if ( status == LINK_CLOSED || status == TEMP_CLOSED ) { HeadLossModel::findClosedHeadLoss(q, hLoss, hGrad); + inertialTerm = MIN_GRADIENT; } else { nw->headLossModel->findHeadLoss(this, q, hLoss, hGrad); if ( hasCheckValve ) HeadLossModel::addCVHeadLoss(q, hLoss, hGrad); + inertialTerm = (length * 4) / (32.174 * PI * diameter * diameter); // Pipe Inertial Term } } + //----------------------------------------------------------------------------- double Pipe::findLeakage(Network* nw, double h, double& dqdh) diff --git a/src/Elements/pipe.h b/src/Elements/pipe.h index 35f367d..506e780 100644 --- a/src/Elements/pipe.h +++ b/src/Elements/pipe.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -37,6 +37,7 @@ class Pipe: public Link void setInitStatus(int s); void setInitSetting(double s); void setResistance(Network* nw); + void setLossFactor(); double getRe(const double q, const double viscos); double getResistance() {return resistance;} diff --git a/src/Elements/pump.cpp b/src/Elements/pump.cpp index 15715a3..294bdad 100644 --- a/src/Elements/pump.cpp +++ b/src/Elements/pump.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -72,6 +72,9 @@ void Pump::setInitFlow() { // ... initial flow is design point of pump curve flow = pumpCurve.qInit * initSetting; + pastFlow = 0; + pastHloss = 0; + pastSetting = 0; } //----------------------------------------------------------------------------- @@ -82,6 +85,7 @@ void Pump::findHeadLoss(Network* nw, double q) if ( speed == 0.0 || status == LINK_CLOSED || status == TEMP_CLOSED ) { HeadLossModel::findClosedHeadLoss(q, hLoss, hGrad); + inertialTerm = MIN_GRADIENT; } // --- get head loss from pump curve and add a check valve @@ -90,6 +94,7 @@ void Pump::findHeadLoss(Network* nw, double q) { pumpCurve.findHeadLoss(speed, q, hLoss, hGrad); if ( !isHpPump() ) HeadLossModel::addCVHeadLoss(q, hLoss, hGrad); + inertialTerm = 10.76; // Approximate value for pump inertial term } } diff --git a/src/Elements/pump.h b/src/Elements/pump.h index acb5afe..d5b594f 100644 --- a/src/Elements/pump.h +++ b/src/Elements/pump.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -43,6 +43,7 @@ class Pump: public Link void setInitFlow(); void setInitStatus(int s); void setInitSetting(double s); + double getSetting(Network* nw) { return speed; } bool isHpPump() { return pumpCurve.isConstHP(); } diff --git a/src/Elements/pumpcurve.cpp b/src/Elements/pumpcurve.cpp index 83131ce..e6f03ce 100644 --- a/src/Elements/pumpcurve.cpp +++ b/src/Elements/pumpcurve.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Elements/pumpcurve.h b/src/Elements/pumpcurve.h index 7d46a34..ba8b565 100644 --- a/src/Elements/pumpcurve.h +++ b/src/Elements/pumpcurve.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Elements/qualsource.cpp b/src/Elements/qualsource.cpp index 9f465ee..2a870ee 100644 --- a/src/Elements/qualsource.cpp +++ b/src/Elements/qualsource.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Elements/qualsource.h b/src/Elements/qualsource.h index 710d96e..e3691d4 100644 --- a/src/Elements/qualsource.h +++ b/src/Elements/qualsource.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Elements/reservoir.cpp b/src/Elements/reservoir.cpp index 97bb4e4..a5b7a10 100644 --- a/src/Elements/reservoir.cpp +++ b/src/Elements/reservoir.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -21,6 +21,8 @@ Reservoir::Reservoir(string name_): { fullDemand = 0.0; fixedGrade = true; + pastHead = head; + ph = head; } Reservoir::~Reservoir() {} diff --git a/src/Elements/reservoir.h b/src/Elements/reservoir.h index b28718e..80cbe8c 100644 --- a/src/Elements/reservoir.h +++ b/src/Elements/reservoir.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -36,6 +36,8 @@ class Reservoir: public Node int type() { return Node::RESERVOIR; } void convertUnits(Network* nw); void setFixedGrade(); + double pastHead; //!< water elev. in previous time period (ft) + double ph; //!< synonym of past head // Properties Pattern* headPattern; //!< time pattern for reservoir's head diff --git a/src/Elements/tank.cpp b/src/Elements/tank.cpp index 83a3041..8c4436f 100644 --- a/src/Elements/tank.cpp +++ b/src/Elements/tank.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -33,6 +33,7 @@ Tank::Tank(string name_) : ucfLength(1.0), pastHead(0.0), pastVolume(0.0), + pastArea(0.0), pastOutflow(0.0) { fullDemand = 0.0; diff --git a/src/Elements/tank.h b/src/Elements/tank.h index e9772ba..8406950 100644 --- a/src/Elements/tank.h +++ b/src/Elements/tank.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -67,6 +67,7 @@ class Tank: public Node double ucfLength; //!< units conversion factor for length double pastHead; //!< water elev. in previous time period (ft) double pastVolume; //!< volume in previous time period (ft3) + double pastArea; //!< area in previous time period (ft2) double pastOutflow; //!< outflow in previous time period (cfs) }; diff --git a/src/Elements/valve.cpp b/src/Elements/valve.cpp index e4db7ed..185c9ba 100644 --- a/src/Elements/valve.cpp +++ b/src/Elements/valve.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -6,6 +6,7 @@ */ #include "valve.h" +#include "pattern.h" #include "node.h" #include "curve.h" #include "Core/network.h" @@ -16,7 +17,9 @@ #include using namespace std; -const char* Valve::ValveTypeWords[] = {"PRV", "PSV", "FCV", "TCV", "PBV", "GPV"}; +// A new valve type is introduced to EPANET. The name of the Valve is Closure Control Valve (CCV) + +const char* Valve::ValveTypeWords[] = {"PRV", "PSV", "FCV", "TCV", "PBV", "GPV", "CCV"}; const double MIN_LOSS_COEFF = 0.1; // default minor loss coefficient @@ -28,6 +31,7 @@ Valve::Valve(string name_) : Link(name_), valveType(TCV), lossFactor(0.0), + settingPattern(nullptr), hasFixedStatus(false), elev(0.0) { @@ -110,6 +114,11 @@ void Valve::setInitSetting(double s) hasFixedStatus = false; } +void Valve::setLossFactor() +{ + lossFactor = 0.02517 * lossCoeff / pow(diameter, 4); +} + //----------------------------------------------------------------------------- // Initialize a valve's state. @@ -134,6 +143,14 @@ void Valve::setInitFlow() { flow = setting; } + else if (valveType == CCV) + { + if (setting == 0) flow = ZERO_FLOW; + else flow = PI * diameter * diameter / 4.0; + } // */ + pastFlow = 0; + pastHloss = 0; + pastSetting = 0; } //----------------------------------------------------------------------------- @@ -178,6 +195,7 @@ double Valve::getSetting(Network* nw) void Valve::findHeadLoss(Network* nw, double q) { hLoss = 0.0; + hGrad = 0.0; // ... valve is temporarily closed (e.g., tries to drain an empty tank) @@ -185,6 +203,8 @@ void Valve::findHeadLoss(Network* nw, double q) if ( status == TEMP_CLOSED) { HeadLossModel::findClosedHeadLoss(q, hLoss, hGrad); + + inertialTerm = MIN_GRADIENT; } // ... valve has fixed status (OPEN or CLOSED) @@ -194,18 +214,37 @@ void Valve::findHeadLoss(Network* nw, double q) if (status == LINK_CLOSED) { HeadLossModel::findClosedHeadLoss(q, hLoss, hGrad); + inertialTerm = MIN_GRADIENT; } - else if (status == LINK_OPEN) findOpenHeadLoss(q); + else if (status == LINK_OPEN) + { + findOpenHeadLoss(q); + inertialTerm = MIN_GRADIENT; + } } // ... head loss for active valves depends on valve type else switch (valveType) { - case PBV: findPbvHeadLoss(q); break; - case TCV: findTcvHeadLoss(q); break; - case GPV: findGpvHeadLoss(nw, q); break; - case FCV: findFcvHeadLoss(q); break; + case PBV: findPbvHeadLoss(q); inertialTerm = MIN_GRADIENT; break; + case TCV: findTcvHeadLoss(q); inertialTerm = MIN_GRADIENT; break; + case CCV: + if(setting == 0) + { + status = LINK_CLOSED; + HeadLossModel::findClosedHeadLoss(q, hLoss, hGrad); + inertialTerm = MIN_GRADIENT; + } + else if (setting != 0) + { + status = LINK_OPEN; + findCcvHeadLoss(nw, q); + inertialTerm = 10.765 / (32.174 * PI * diameter * diameter); // Approximate value for valve inertial term + } + break; + case GPV: findGpvHeadLoss(nw, q); inertialTerm = MIN_GRADIENT; break; + case FCV: findFcvHeadLoss(q); inertialTerm = MIN_GRADIENT; break; // ... PRVs & PSVs without fixed status can be either // OPEN, CLOSED, or ACTIVE. @@ -214,8 +253,10 @@ void Valve::findHeadLoss(Network* nw, double q) if ( status == LINK_CLOSED ) { HeadLossModel::findClosedHeadLoss(q, hLoss, hGrad); + inertialTerm = MIN_GRADIENT; } else if ( status == LINK_OPEN ) findOpenHeadLoss(q); + inertialTerm = MIN_GRADIENT; break; } } @@ -285,6 +326,44 @@ void Valve::findTcvHeadLoss(double q) //----------------------------------------------------------------------------- +// Find the head loss and its gradient for a closure control valve. + +void Valve::findCcvHeadLoss(Network* nw, double q) +{ + //... save open valve loss factor + + double f = lossFactor; + + if (nw->option(Options::VALVE_REP_TYPE) == "Toe") + { + double valve_conductance = 16.96; // ft^2.5 / s = 0.87 m^(2.5) / s + double vc_square = valve_conductance * valve_conductance; + double Toe = setting; // Globe Valve + double Toe_Square = Toe * Toe; + lossFactor = 1 / (vc_square * Toe_Square); // (Nault and Karney, 2016) */ + } + + else if (nw->option(Options::VALVE_REP_TYPE) == "Cd") + { + double d2 = diameter * diameter; + double FullArea = (d2 * PI) / 4; // Area of valve + //double Area = FullArea * (-0.3575 * pow(setting, 4) + 0.2766 * pow(setting, 3) - 0.2233 * pow(setting, 2) + 1.3056 * setting - 0.0005); + double FullArea2 = FullArea * FullArea; + double Cd = -1.1293 * pow(setting, 6) + 3.3823 * pow(setting, 5) - 3.443 * pow(setting, 4) + 0.5671 * pow(setting, 3) + 1.0371 * pow(setting, 2) - 0.0037 * setting; // Globe Valve according to Tullis (1989) + double Cd2 = Cd * Cd; + double g = 32.174; // gravitational acceleration + lossFactor = (1 / Cd2 - 1) * (1 / (2 * g * FullArea2)); // */ + } + + // setting value is between 0 and 1 + + findOpenHeadLoss(q); + + // ... restore open valve loss factor + + //lossFactor = f; +} + // Find the head loss and its derivative for a general purpose valve. void Valve::findGpvHeadLoss(Network* nw, double q) @@ -431,21 +510,61 @@ int Valve::updatePsvStatus(double q, double h1, double h2) // Change the setting of a valve. -bool Valve::changeSetting( - double newSetting, bool makeChange, const string reason, ostream& msgLog) +bool Valve::changeSetting(double newSetting, bool makeChange, const string reason, ostream& msgLog) { - if ( newSetting != setting ) - { - if ( makeChange ) - { - hasFixedStatus = false; - status = Link::LINK_OPEN; - msgLog << "\n " << reason; - setting = newSetting; - } - return true; - } - return false; + if (valveType == CCV) + { + if (setting != newSetting) + { + if (status == Link::LINK_CLOSED && newSetting == 0.0) + { + setting = newSetting; + return false; + } + + if (makeChange) + { + if (newSetting == 0.0) + { + status = Link::LINK_CLOSED; + flow = ZERO_FLOW; + } + else status = Link::LINK_OPEN; + msgLog << "\n " << reason; + setting = newSetting; + } + return true; + } + return false; + } + + else + { + if (setting != newSetting) + { + if (status == Link::LINK_CLOSED) // && newSetting == 0.0) + { + setting = newSetting; + return false; + } + + if (makeChange) + { + if (newSetting == 0.0) + { + status = Link::LINK_CLOSED; + flow = ZERO_FLOW; + } + else status = Link::LINK_OPEN; + msgLog << "\n " << reason; + setting = newSetting; + } + return true; + } + return false; + } + + } //----------------------------------------------------------------------------- @@ -488,3 +607,13 @@ void Valve::validateStatus(Network* nw, double qTol) default: break; } } + +void Valve::applyControlPattern(ostream& msgLog) +{ + if (settingPattern) + { + changeSetting(settingPattern->currentFactor(), true, "setting pattern", msgLog); + + + } +} diff --git a/src/Elements/valve.h b/src/Elements/valve.h index 593a97d..6cb1b3c 100644 --- a/src/Elements/valve.h +++ b/src/Elements/valve.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -16,6 +16,7 @@ #include class Network; +class Pattern; //! \class Valve //! \brief A Link that controls flow or pressure. @@ -32,7 +33,8 @@ class Valve: public Link FCV, //!< flow control valve TCV, //!< throttle control valve PBV, //!< pressure breaker valve - GPV //!< general purpose valve + GPV, //!< general purpose valve + CCV //!< closure control valve }; static const char* ValveTypeWords[]; @@ -49,6 +51,7 @@ class Valve: public Link void setInitFlow(); void setInitStatus(int s); void setInitSetting(double s); + void setLossFactor(); void initialize(bool initFlow); bool isPRV(); @@ -65,6 +68,7 @@ class Valve: public Link const std::string reason, std::ostream& msgLog); void validateStatus(Network* nw, double qTol); + bool makeChange; double getVelocity(); double getRe(const double q, const double viscos); @@ -73,6 +77,9 @@ class Valve: public Link // Properties ValveType valveType; //!< valve type double lossFactor; //!< minor loss factor + Pattern* settingPattern; //!< Setting Pattern + void applyControlPattern(std::ostream& msgLog); + protected: void findOpenHeadLoss(double q); @@ -80,6 +87,7 @@ class Valve: public Link void findTcvHeadLoss(double q); void findGpvHeadLoss(Network* nw, double q); void findFcvHeadLoss(double q); + void findCcvHeadLoss(Network* nw, double q); int updatePrvStatus(double q, double h1, double h2); int updatePsvStatus(double q, double h1, double h2); diff --git a/src/Input/controlparser.cpp b/src/Input/controlparser.cpp index 716e650..8de2d0d 100644 --- a/src/Input/controlparser.cpp +++ b/src/Input/controlparser.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Input/controlparser.h b/src/Input/controlparser.h index 9496620..06f00ab 100644 --- a/src/Input/controlparser.h +++ b/src/Input/controlparser.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Input/curveparser.cpp b/src/Input/curveparser.cpp index 29d6f81..777bcc3 100644 --- a/src/Input/curveparser.cpp +++ b/src/Input/curveparser.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Input/curveparser.h b/src/Input/curveparser.h index 58f51ca..6b13906 100644 --- a/src/Input/curveparser.h +++ b/src/Input/curveparser.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Input/inputparser.cpp b/src/Input/inputparser.cpp index 30c5c97..03681a2 100644 --- a/src/Input/inputparser.cpp +++ b/src/Input/inputparser.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). diff --git a/src/Input/inputparser.h b/src/Input/inputparser.h index b477df7..63abd14 100644 --- a/src/Input/inputparser.h +++ b/src/Input/inputparser.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Input/inputreader.cpp b/src/Input/inputreader.cpp index aeff2b1..ca43937 100644 --- a/src/Input/inputreader.cpp +++ b/src/Input/inputreader.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). diff --git a/src/Input/inputreader.h b/src/Input/inputreader.h index 8207a7c..5c97359 100644 --- a/src/Input/inputreader.h +++ b/src/Input/inputreader.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Input/linkparser.cpp b/src/Input/linkparser.cpp index 8e94104..a065ad6 100644 --- a/src/Input/linkparser.cpp +++ b/src/Input/linkparser.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -29,7 +29,7 @@ static const char* w_Effic = "EFFIC"; static const char* w_OPEN = "OPEN"; static const char* w_CLOSED = "CLOSED"; static const char* w_CV = "CV"; -static const char* valveTypeWords[] = {"PRV", "PSV", "FCV", "TCV", "PBV", "GPV", 0}; +static const char* valveTypeWords[] = {"PRV", "PSV", "FCV", "TCV", "PBV", "GPV", "CCV", 0}; //----------------------------------------------------------------------------- // Local functions @@ -360,60 +360,74 @@ void parsePumpData(Pump* pump, Network* nw, vector& tokenList) void parseValveData(Valve* valve, Network* network, vector& tokenList) { - // Contents of tokenList are: - // 0 - valve ID - // 1 - upstream node ID - // 2 - downstream node ID - // 3 - diameter - // 4 - valve type - // 5 - valve setting - // 6 - minor loss coeff. (optional) - - // ... check for enough input tokens - - if ( tokenList.size() < 6 ) throw InputError(InputError::TOO_FEW_ITEMS, ""); - string* tokens = &tokenList[0]; - - // ... read diameter - - if ( !Utilities::parseNumber(tokens[3], valve->diameter)|| - valve->diameter <= 0.0 ) - { - throw InputError(InputError::INVALID_NUMBER, tokens[3]); - } - - // ... read valve type - - int vType = Utilities::findMatch(tokens[4], valveTypeWords); - if ( vType < 0 ) throw InputError(InputError::INVALID_KEYWORD, tokens[4]); - valve->valveType = (Valve::ValveType)vType; - - // ... read index of head loss curve for General Purpose Valve - - if ( valve->valveType == Valve::GPV ) - { - int c = network->indexOf(Element::CURVE, tokens[5]); - if ( c < 0 ) throw InputError(InputError::UNDEFINED_OBJECT, tokens[5]); - valve->initSetting = c; - } - - // ... read numerical setting for other types of valves - else - { - if ( !Utilities::parseNumber(tokens[5], valve->initSetting) ) - { - throw InputError(InputError::INVALID_NUMBER, tokens[5]); - } - } - - // ... read optional minor loss coeff. - - if ( tokenList.size() > 6 ) - { - if ( !Utilities::parseNumber(tokens[6], valve->lossCoeff) || - valve->lossCoeff < 0.0 ) - { - throw InputError(InputError::INVALID_NUMBER, tokens[6]); - } - } + // Contents of tokenList are: + // 0 - valve ID + // 1 - upstream node ID + // 2 - downstream node ID + // 3 - diameter + // 4 - valve type + // 5 - valve setting + // 6 - minor loss coeff. (optional) + // 7 - valve setting pattern (optional) + + // ... check for enough input tokens + + if (tokenList.size() < 6) throw InputError(InputError::TOO_FEW_ITEMS, ""); + string* tokens = &tokenList[0]; + + // ... read diameter + + if (!Utilities::parseNumber(tokens[3], valve->diameter) || + valve->diameter <= 0.0) + { + throw InputError(InputError::INVALID_NUMBER, tokens[3]); + } + + // ... read valve type + + int vType = Utilities::findMatch(tokens[4], valveTypeWords); + if (vType < 0) throw InputError(InputError::INVALID_KEYWORD, tokens[4]); + valve->valveType = (Valve::ValveType)vType; + + // ... read index of head loss curve for General Purpose Valve + + if (valve->valveType == Valve::GPV) + { + int c = network->indexOf(Element::CURVE, tokens[5]); + if (c < 0) throw InputError(InputError::UNDEFINED_OBJECT, tokens[5]); + valve->initSetting = c; + } + + // ... read numerical setting for other types of valves + else + { + if (!Utilities::parseNumber(tokens[5], valve->initSetting)) + { + throw InputError(InputError::INVALID_NUMBER, tokens[5]); + } + } + + // ... read optional minor loss coeff. + + if (tokenList.size() > 6) + { + if (!Utilities::parseNumber(tokens[6], valve->lossCoeff) || + valve->lossCoeff < 0.0) + { + throw InputError(InputError::INVALID_NUMBER, tokens[6]); + } + } + + + // read optional setting pattern + + if (tokenList.size() > 7 && tokens[7] != "*") + { + valve->settingPattern = network->pattern(tokens[7]); + if (!valve->settingPattern) + { + throw InputError(InputError::UNDEFINED_OBJECT, tokens[2]); + } + } } + diff --git a/src/Input/linkparser.h b/src/Input/linkparser.h index 16af078..0564857 100644 --- a/src/Input/linkparser.h +++ b/src/Input/linkparser.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Input/nodeparser.cpp b/src/Input/nodeparser.cpp index 5a40b2b..2a629cd 100644 --- a/src/Input/nodeparser.cpp +++ b/src/Input/nodeparser.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Input/nodeparser.h b/src/Input/nodeparser.h index c478cb7..ad6137f 100644 --- a/src/Input/nodeparser.h +++ b/src/Input/nodeparser.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Input/optionparser.cpp b/src/Input/optionparser.cpp index 822be2a..4abaf36 100644 --- a/src/Input/optionparser.cpp +++ b/src/Input/optionparser.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -22,7 +22,7 @@ static const char* stringOptionKeywords[] = {"HYDRAULICS_FILE", "", "", // placeholders for file names "MAP_FILE", "HEADLOSS_MODEL", "DEMAND_MODEL", "LEAKAGE_MODEL", - "HYDRAULIC_SOLVER", "STEP_SIZING", "MATRIX_SOLVER", "", + "HYD_SOLVER", "STEP_SIZING", "VALVE_REP_TYPE", "MATRIX_SOLVER", "", "QUALITY_MODEL", "QUALITY_NAME", "QUALITY_UNITS", 0}; // ... Keywords for IndexOption enumeration in options.h @@ -46,8 +46,8 @@ static const char* valueOptionKeywords[] = "MINIMUM_PRESSURE", "SERVICE_PRESSURE", "PRESSURE_EXPONENT", "EMITTER_EXPONENT", "LEAKAGE_COEFF1", "LEAKAGE_COEFF2", "RELATIVE_ACCURACY", "HEAD_TOLERANCE", "FLOW_TOLERANCE", - "FLOW_CHANGE_LIMIT", "TIME_WEIGHT", "SPECIFIC_DIFFUSIVITY", - "QUALITY_TOLERANCE", 0}; + "FLOW_CHANGE_LIMIT", "TIME_WEIGHT", "TEMP_DISC_PARA", "SPECIFIC_DIFFUSIVITY", + "QUALITY_TOLERANCE", 0}; // ... Keywords for TimeOption enumeration in options.h static const char* timeOptionKeywords[] = diff --git a/src/Input/optionparser.h b/src/Input/optionparser.h index 478add5..184d151 100644 --- a/src/Input/optionparser.h +++ b/src/Input/optionparser.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Input/patternparser.cpp b/src/Input/patternparser.cpp index 7747795..970f618 100644 --- a/src/Input/patternparser.cpp +++ b/src/Input/patternparser.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Input/patternparser.h b/src/Input/patternparser.h index 5867cda..9b17887 100644 --- a/src/Input/patternparser.h +++ b/src/Input/patternparser.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Models/demandmodel.cpp b/src/Models/demandmodel.cpp index 1251684..d8dd522 100644 --- a/src/Models/demandmodel.cpp +++ b/src/Models/demandmodel.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). diff --git a/src/Models/demandmodel.h b/src/Models/demandmodel.h index 6f4ee64..7ff8f67 100644 --- a/src/Models/demandmodel.h +++ b/src/Models/demandmodel.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Models/headlossmodel.cpp b/src/Models/headlossmodel.cpp index 897642c..2be350b 100644 --- a/src/Models/headlossmodel.cpp +++ b/src/Models/headlossmodel.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Models/headlossmodel.h b/src/Models/headlossmodel.h index 9987a47..ff70d55 100644 --- a/src/Models/headlossmodel.h +++ b/src/Models/headlossmodel.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Models/leakagemodel.cpp b/src/Models/leakagemodel.cpp index 8b71b6b..8de04eb 100644 --- a/src/Models/leakagemodel.cpp +++ b/src/Models/leakagemodel.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Models/leakagemodel.h b/src/Models/leakagemodel.h index 75dc4b2..ff9e255 100644 --- a/src/Models/leakagemodel.h +++ b/src/Models/leakagemodel.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). diff --git a/src/Models/pumpenergy.cpp b/src/Models/pumpenergy.cpp index 9279a44..9a27572 100644 --- a/src/Models/pumpenergy.cpp +++ b/src/Models/pumpenergy.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Models/pumpenergy.h b/src/Models/pumpenergy.h index c76228f..227b732 100644 --- a/src/Models/pumpenergy.h +++ b/src/Models/pumpenergy.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Models/qualmodel.cpp b/src/Models/qualmodel.cpp index a9a3816..4910f8f 100644 --- a/src/Models/qualmodel.cpp +++ b/src/Models/qualmodel.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Models/qualmodel.h b/src/Models/qualmodel.h index f1ad0fe..3f50974 100644 --- a/src/Models/qualmodel.h +++ b/src/Models/qualmodel.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Models/tankmixmodel.cpp b/src/Models/tankmixmodel.cpp index dfd28c1..68ebd7a 100644 --- a/src/Models/tankmixmodel.cpp +++ b/src/Models/tankmixmodel.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Models/tankmixmodel.h b/src/Models/tankmixmodel.h index 26458a0..1926811 100644 --- a/src/Models/tankmixmodel.h +++ b/src/Models/tankmixmodel.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Output/outputfile.cpp b/src/Output/outputfile.cpp index 32b9821..a77c985 100644 --- a/src/Output/outputfile.cpp +++ b/src/Output/outputfile.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Output/outputfile.h b/src/Output/outputfile.h index e41d072..d75f7dc 100644 --- a/src/Output/outputfile.h +++ b/src/Output/outputfile.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Output/projectwriter.cpp b/src/Output/projectwriter.cpp index f29c58c..3ba3095 100644 --- a/src/Output/projectwriter.cpp +++ b/src/Output/projectwriter.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -257,6 +257,12 @@ void ProjectWriter::writeValves() link->convertSetting(network, link->initSetting); fout << setw(12) << cf * link->initSetting << "\n"; } + if (valve->settingPattern) + { + fout << setw(8) << "PATTERN"; + fout << setw(16) << valve->settingPattern->name; + } + fout << "\n"; } } } diff --git a/src/Output/projectwriter.h b/src/Output/projectwriter.h index 5b73b74..87c5fb6 100644 --- a/src/Output/projectwriter.h +++ b/src/Output/projectwriter.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Output/reportfields.cpp b/src/Output/reportfields.cpp index 03a0f6f..ac01889 100644 --- a/src/Output/reportfields.cpp +++ b/src/Output/reportfields.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Output/reportfields.h b/src/Output/reportfields.h index 7db4571..fdd8adf 100644 --- a/src/Output/reportfields.h +++ b/src/Output/reportfields.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Output/reportwriter.cpp b/src/Output/reportwriter.cpp index 82c9ea6..9be4d9f 100644 --- a/src/Output/reportwriter.cpp +++ b/src/Output/reportwriter.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -43,7 +43,7 @@ void ReportWriter::writeHeading() sout << "\n * E P A N E T *"; sout << "\n * Hydraulic and Water Quality *"; sout << "\n * Analysis for Pipe Networks *"; - sout << "\n * Version 3.0.000 *"; + sout << "\n * Version 3.1.000 *"; sout << "\n ******************************************************************\n"; } @@ -96,6 +96,7 @@ void ReportWriter::writeSummary(string inpFileName) sout << "\n Number of Pumps ............... " << nPumps; sout << "\n Number of Valves .............. " << nValves; + sout << "\n Hydraulic Solver .............. " << network->option(Options::HYD_SOLVER); // Two different solver could be selected by this option (GGA and RWC-GGA) sout << "\n Head Loss Model ............... " << network->option(Options::HEADLOSS_MODEL); sout << "\n Leakage Model ................. " << network->option(Options::LEAKAGE_MODEL); sout << "\n Demand Model .................. " << network->option(Options::DEMAND_MODEL); @@ -104,6 +105,9 @@ void ReportWriter::writeSummary(string inpFileName) sout << "\n Head Tolerance ................ " << network->option(Options::HEAD_TOLERANCE); sout << "\n Flow Tolerance ................ " << network->option(Options::FLOW_TOLERANCE); sout << "\n Flow Change Limit ............. " << network->option(Options::FLOW_CHANGE_LIMIT); + sout << "\n Temporal Disc. Para. .......... " << network->option(Options::TEMP_DISC_PARA); // A Discretization Parameter which is in between 0.5 and 1. This option is valid for only RWC-GGA Hydraulic Solver. + sout << "\n Step Sizing Method ............ " << network->option(Options::STEP_SIZING); // An option to provide convergence in case of pressure dependent demand. Full, ARF, and BRF are preferences. + sout << "\n CCV Representation Type ....... " << network->option(Options::VALVE_REP_TYPE); // CCV has two types: Toe and Cd Valves. sout << "\n Quality Model ................. " << network->option(Options::QUAL_MODEL); if ( network->option(Options::QUAL_TYPE) == Options::TRACE ) diff --git a/src/Output/reportwriter.h b/src/Output/reportwriter.h index 984e68e..39f8504 100644 --- a/src/Output/reportwriter.h +++ b/src/Output/reportwriter.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Solvers/ggasolver.cpp b/src/Solvers/ggasolver.cpp index c60839c..bb6210c 100644 --- a/src/Solvers/ggasolver.cpp +++ b/src/Solvers/ggasolver.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). @@ -6,12 +6,9 @@ */ //////////////////////////////////////////////////////////////////////// - // Implementation of the Global Gradient Algorithm hydraulic solver // + // Implementation of the Global Gradient Algorithm hydraulic solver and Convergence Tracking Control Method // //////////////////////////////////////////////////////////////////////// - // TO DO: - // - implement the line search procedure and consider moving it to - // its own module so it can be used by other solvers #include "ggasolver.h" #include "matrixsolver.h" @@ -21,6 +18,7 @@ #include "Elements/junction.h" #include "Elements/tank.h" #include "Elements/link.h" +#include "Elements/valve.h" #include #include @@ -51,7 +49,7 @@ static const double ErrorThreshold = 1.0; static const double Huge = numeric_limits::max(); // step sizing enumeration -enum StepSizing {FULL, RELAXATION, LINESEARCH}; +enum StepSizing {FULL, RELAXATION, LINESEARCH, BRF, ARF}; //----------------------------------------------------------------------------- @@ -81,6 +79,10 @@ GGASolver::GGASolver(Network* nw, MatrixSolver* ms) : HydSolver(nw, ms) stepSizing = RELAXATION; else if (network->option(Options::STEP_SIZING) == "LINESEARCH" ) stepSizing = LINESEARCH; + else if (network->option(Options::STEP_SIZING) == "BRF") + stepSizing = BRF; + else if (network->option(Options::STEP_SIZING) == "ARF") + stepSizing = ARF; else stepSizing = FULL; errorNorm = 0.0; @@ -102,7 +104,7 @@ GGASolver::~GGASolver() // Solve network for heads and flows -int GGASolver::solve(double tstep_, int& trials) +int GGASolver::solve(double tstep_, int& trials, int currentTime) { // ... initialize variables @@ -121,6 +123,9 @@ int GGASolver::solve(double tstep_, int& trials) theta = min(theta, 1.0); if ( theta > 0.0 ) theta = max(theta, 0.5); + minErrorNorm = 1000000000; + dl = 1.0; + // ... set values for convergence limits setConvergenceLimits(); @@ -141,44 +146,129 @@ int GGASolver::solve(double tstep_, int& trials) if ( statusChanged ) { - oldErrorNorm = findErrorNorm(0.0); + oldErrorNorm = findErrorNorm(0.0, currentTime, tstep); lamda = 1.0; } statusChanged = false; - - // ... find changes in heads and flows - - int errorCode = findHeadChanges(); - if ( errorCode >= 0 ) - { - Node* node = network->node(errorCode); - network->msgLog << endl << s_IllConditioned << node->name; - return HydSolver::FAILED_ILL_CONDITIONED; - } - findFlowChanges(); - - // ... find step size to take for head/flow changes - // (which evaluates new gradients for next trial) - - lamda = findStepSize(trials); - updateSolution(lamda); - - // ... check for convergence - - if ( reportTrials ) reportTrial(trials, lamda); - converged = hasConverged(); - - // ... if close to convergence then check for any link status changes - - if ( converged ) //|| errorNorm < ErrorThreshold ) - { - statusChanged = linksChangedStatus(); - } - - // ... check if the current solution can be accepted - - if ( converged && !statusChanged ) break; - trials++; + dl = 1.000000000000000; + + int errorCode = findHeadChanges(); + if (errorCode >= 0) + { + Node* node = network->node(errorCode); + network->msgLog << endl << s_IllConditioned << node->name; + return HydSolver::FAILED_ILL_CONDITIONED; + } + findFlowChanges(); + + if (stepSizing == ARF) + { + double lamda = 1.0; + errorNorm = findErrorNorm(lamda, currentTime, tstep); + updateSolution(lamda); + + if (errorNorm < oldErrorNorm) + { + if (reportTrials) reportTrial(trials, lamda); + converged = hasConverged(); + + // ... if close to convergence then check for any link status changes + + if (converged) //|| errorNorm < ErrorThreshold ) + { + statusChanged = linksChangedStatus(); + } + + // ... check if the current solution can be accepted + + if (converged && !statusChanged) break; + trials++; + } + else + { + // ... find changes in heads and flows + int errorCode = findHeadChanges(); + if (errorCode >= 0) + { + Node* node = network->node(errorCode); + network->msgLog << endl << s_IllConditioned << node->name; + return HydSolver::FAILED_ILL_CONDITIONED; + } + findFlowChanges(); // */ + int counter = 0; + while (minErrorNorm >= oldErrorNorm) + { + dl *= 0.25; + if (dl < 0.001) break; + lamda = findStepSize(trials, currentTime); + updateSolution(lamda); + errorNorm = minErrorNorm; + } + if (reportTrials) reportTrial(trials, lamda); + converged = hasConverged(); + + // ... if close to convergence then check for any link status changes + + if (converged) //|| errorNorm < ErrorThreshold ) + { + statusChanged = linksChangedStatus(); + } + + // ... check if the current solution can be accepted + + if (converged && !statusChanged) break; + trials++; + } + } + + else if (stepSizing == BRF) + { + int counter = 0; + while (minErrorNorm >= oldErrorNorm) + { + dl *= 0.25; + if (dl < 0.001) break; + lamda = findStepSize(trials, currentTime); + updateSolution(lamda); + errorNorm = minErrorNorm; // */ + } + if (reportTrials) reportTrial(trials, lamda); + converged = hasConverged(); + + // ... if close to convergence then check for any link status changes + + if (converged) // || errorNorm < ErrorThreshold ) + { + statusChanged = linksChangedStatus(); + } + + // ... check if the current solution can be accepted + + if (converged && !statusChanged) break; + trials++; + } + else + { + lamda = findStepSize(trials, currentTime); + updateSolution(lamda); + + // ... check for convergence + + if (reportTrials) reportTrial(trials, lamda); + converged = hasConverged(); + + // ... if close to convergence then check for any link status changes + + if (converged) //|| errorNorm < ErrorThreshold ) + { + statusChanged = linksChangedStatus(); + } + + // ... check if the current solution can be accepted + + if (converged && !statusChanged) break; + trials++; + } } //if ( reportTrials ) network->msgLog << s_HlossEvals << hLossEvalCount; if ( trials > trialsLimit ) return HydSolver::FAILED_NO_CONVERGENCE; @@ -341,44 +431,79 @@ void GGASolver::findFlowChanges() // Find how much of the head and flow changes to apply to a new solution. -double GGASolver::findStepSize(int trials) +double GGASolver::findStepSize(int trials, int currentTime) { // ... find the new error norm at full step size double lamda = 1.0; - errorNorm = findErrorNorm(lamda); + errorNorm = findErrorNorm(lamda, currentTime, tstep); if ( stepSizing == RELAXATION && oldErrorNorm < ErrorThreshold ) { lamda = 0.5; - double errorNorm2 = findErrorNorm(lamda); + double errorNorm2 = findErrorNorm(lamda, currentTime, tstep); if ( errorNorm2 < errorNorm ) errorNorm = errorNorm2; else { lamda = 1.0; - errorNorm = findErrorNorm(lamda); + errorNorm = findErrorNorm(lamda, currentTime, tstep); } } // ... if called for, implement a line search procedure // to find the best step size lamda to take - if ( stepSizing == LINESEARCH && trials > 1 ) - { - //////////// TO BE IMPLEMENTED /////////////// - } - return lamda; + if (stepSizing == BRF || stepSizing == ARF) + { + { + minErrorNorm = 0; + double testError = 1000000; + lambdaNumber = 1 / dl; + Lambda.resize(lambdaNumber, 0); + + + memset(&Lambda[0], 0, lambdaNumber*sizeof(double)); + + for (int i = 0; i < lambdaNumber; i++) + { + Lambda[i] += (i + 1)*dl; + + int errorCode = findHeadChanges(); + if (errorCode >= 0) + { + Node* node = network->node(errorCode); + network->msgLog << endl << s_IllConditioned << node->name; + return HydSolver::FAILED_ILL_CONDITIONED; + } + findFlowChanges(); //*/ + + updateSolution(Lambda[i]); + + errorNorm = findErrorNorm(Lambda[i], currentTime, tstep); + + if (errorNorm < testError) + { + testError = errorNorm; + lamda = Lambda[i]; + updateSolution(Lambda[i]); + } + } + minErrorNorm = testError; + } + return lamda; + } + return lamda; } //----------------------------------------------------------------------------- // Compute the error norm associated with a given step size. -double GGASolver::findErrorNorm(double lamda) +double GGASolver::findErrorNorm(double lamda, int currentTime, double tstep) { hLossEvalCount++; return hydBalance.evaluate(lamda, (double*)&dH[0], (double*)&dQ[0], - (double*)&xQ[0], network); + (double*)&xQ[0], network, currentTime, tstep); } //----------------------------------------------------------------------------- diff --git a/src/Solvers/ggasolver.h b/src/Solvers/ggasolver.h index cc36045..baaada9 100644 --- a/src/Solvers/ggasolver.h +++ b/src/Solvers/ggasolver.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -27,7 +27,7 @@ class GGASolver : public HydSolver GGASolver(Network* nw, MatrixSolver* ms); ~GGASolver(); - int solve(double tstep, int& trials); + int solve(double tstep, int& trials, int currentTime); private: @@ -35,6 +35,7 @@ class GGASolver : public HydSolver int linkCount; // number of network links int hLossEvalCount; // number of head loss evaluations int stepSizing; // Newton step sizing method + int lambdaNumber; int trialsLimit; // limit on number of trials bool reportTrials; // report summary of each trial @@ -44,6 +45,8 @@ class GGASolver : public HydSolver double flowRatioLimit; // allowable total flow change / total flow double tstep; // time step (sec) double theta; // time weighting constant + double minErrorNorm; + double dl; double errorNorm; // solution error norm double oldErrorNorm; // previous error norm @@ -52,6 +55,7 @@ class GGASolver : public HydSolver std::vector dH; // head change at each node (ft) std::vector dQ; // flow change in each link (cfs) std::vector xQ; // node flow imbalances (cfs) + std::vector Lambda; // Functions that assemble linear equation coefficients void setFixedGradeNodes(); @@ -63,12 +67,12 @@ class GGASolver : public HydSolver // Functions that update the hydraulic solution int findHeadChanges(); void findFlowChanges(); - double findStepSize(int trials); + double findStepSize(int trials, int currentTime); void updateSolution(double lamda); // Functions that check for convergence void setConvergenceLimits(); - double findErrorNorm(double lamda); + double findErrorNorm(double lamda, int currentTime, double tstep); bool hasConverged(); bool linksChangedStatus(); void reportTrial(int trials, double lamda); diff --git a/src/Solvers/hydsolver.cpp b/src/Solvers/hydsolver.cpp index b130c34..ccf09b0 100644 --- a/src/Solvers/hydsolver.cpp +++ b/src/Solvers/hydsolver.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -9,6 +9,7 @@ // Include header files for the different hydraulic solvers here. #include "ggasolver.h" +#include "rwcggasolver.h" using namespace std; @@ -21,5 +22,9 @@ HydSolver::~HydSolver() {} HydSolver* HydSolver::factory(const string name, Network* nw, MatrixSolver* ms) { if (name == "GGA") return new GGASolver(nw, ms); - return nullptr; + // return nullptr; + + else if (name == "RWCGGA") return new RWCGGASolver(nw, ms); + return nullptr; + } diff --git a/src/Solvers/hydsolver.h b/src/Solvers/hydsolver.h index 7cbff35..a88f4e1 100644 --- a/src/Solvers/hydsolver.h +++ b/src/Solvers/hydsolver.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -36,7 +36,7 @@ class HydSolver HydSolver(Network* nw, MatrixSolver* ms); virtual ~HydSolver(); static HydSolver* factory(const std::string name, Network* nw, MatrixSolver* ms); - virtual int solve(double tstep, int& trials) = 0; + virtual int solve(double tstep, int& trials, int currentTime) = 0; protected: diff --git a/src/Solvers/ltdsolver.cpp b/src/Solvers/ltdsolver.cpp index a8f0083..8e1a963 100644 --- a/src/Solvers/ltdsolver.cpp +++ b/src/Solvers/ltdsolver.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Solvers/ltdsolver.h b/src/Solvers/ltdsolver.h index b7cf192..72878e8 100644 --- a/src/Solvers/ltdsolver.h +++ b/src/Solvers/ltdsolver.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Solvers/matrixsolver.cpp b/src/Solvers/matrixsolver.cpp index f010a26..7c3e0fd 100644 --- a/src/Solvers/matrixsolver.cpp +++ b/src/Solvers/matrixsolver.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Solvers/matrixsolver.h b/src/Solvers/matrixsolver.h index a5b7a23..04fbed4 100644 --- a/src/Solvers/matrixsolver.h +++ b/src/Solvers/matrixsolver.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Solvers/qualsolver.cpp b/src/Solvers/qualsolver.cpp index 747ae69..240b9ee 100644 --- a/src/Solvers/qualsolver.cpp +++ b/src/Solvers/qualsolver.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Solvers/qualsolver.h b/src/Solvers/qualsolver.h index a336c60..0968c84 100644 --- a/src/Solvers/qualsolver.h +++ b/src/Solvers/qualsolver.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Solvers/rwcggasolver.cpp b/src/Solvers/rwcggasolver.cpp new file mode 100644 index 0000000..32ae2d3 --- /dev/null +++ b/src/Solvers/rwcggasolver.cpp @@ -0,0 +1,915 @@ +/* EPANET 3.1 + * + * Copyright (c) 2016 Open Water Analytics + * Distributed under the MIT License (see the LICENSE file for details). + * + */ + + //////////////////////////////////////////////////////////////////////// + // Implementation of the Rigid Water Column Global Gradient Algorithm hydraulic solver and Convergence Tracking Control Method // + //////////////////////////////////////////////////////////////////////// + + + + +#include "epanet3.h" +#include "rwcggasolver.h" +#include "matrixsolver.h" +#include "Core/network.h" +#include "Core/constants.h" +#include "Elements/control.h" +#include "Elements/junction.h" +#include "Core/project.h" +#include "Core/datamanager.h" +#include "Core/constants.h" +#include "Core/error.h" +#include "Utilities/utilities.h" + +#include "Elements/tank.h" +#include "Elements/link.h" +#include "Core/hydengine.h" +#include "Core/hydbalance.h" +#include "Elements/valve.h" + + +#include +#include +#include +#include //for debugging +#include +#include +#include + + +using namespace std; + +static const string s_Trial = " Trial "; +static const string s_IllConditioned = " Hydraulic matrix ill-conditioned at node "; +static const string s_StepSize = " Step Size = "; +static const string s_TotalError = " Error Norm = "; +static const string s_HlossEvals = " Head Loss Evaluations = "; +static const string s_HeadError = " Head Error = "; +static const string s_ForLink = " for Link "; +static const string s_FlowError = " Flow Error = "; +static const string s_AtNode = " at Node "; +static const string s_FlowChange = " Flow Change = "; +static const string s_TotFlowChange = " Total Flow Change Ratio = "; +static const string s_NodeLabel = " Node "; +static const string s_FGChange = " Fixed Grade Status changed to "; + +//----------------------------------------------------------------------------- + +// error norm threshold for status checking +static const double ErrorThreshold = 1.0; +static const double Huge = numeric_limits::max(); + +// step sizing enumeration +enum StepSizing {FULL, RELAXATION, LINESEARCH, BRF, ARF}; + +//----------------------------------------------------------------------------- + +// Constructor + +RWCGGASolver::RWCGGASolver(Network* nw, MatrixSolver* ms) : HydSolver(nw, ms) +{ + nodeCount = network->count(Element::NODE); + linkCount = network->count(Element::LINK); + + dH.resize(nodeCount, 0); // nodal head changes + dQ.resize(linkCount, 0); // link flow changes + xQ.resize(nodeCount, 0); // nodal excess flow (inflow - outflow) + + hLossEvalCount = 0; + trialsLimit = 0; + reportTrials = network->option(Options::REPORT_TRIALS); + + headErrLimit = 0.0; + flowErrLimit = 0.0; + flowChangeLimit = 0.0; + flowRatioLimit = 0.0; + tstep = 0.0; + theta = 0.0; + kappa = 0.0; + dhmaxpast = 0.0; + maxHeadErrPast = 0.0; + + if (network->option(Options::STEP_SIZING) == "RELAXATION") + stepSizing = RELAXATION; + else if (network->option(Options::STEP_SIZING) == "LINESEARCH") + stepSizing = LINESEARCH; + else if (network->option(Options::STEP_SIZING) == "BRF") + stepSizing = BRF; + else if (network->option(Options::STEP_SIZING) == "ARF") + stepSizing = ARF; + else stepSizing = FULL; + + errorNorm = 0.0; + oldErrorNorm = 0.0; +} + +//----------------------------------------------------------------------------- + +// Destructor + +RWCGGASolver::~RWCGGASolver() +{ + dH.clear(); + dQ.clear(); + xQ.clear(); +} + +std::ofstream dosyam("D:\\EPANET_3\\Networks\\RWC\\kappa.txt"); + +//----------------------------------------------------------------------------- + +// Solve network for heads and flows + +int RWCGGASolver::solve(double tstep_, int& trials, int currentTime) + +{ + // ... initialize variables + + double lamda = 1.0; + bool statusChanged = true; + bool converged = false; + + errorNorm = Huge; + hLossEvalCount = 0; + tstep = tstep_; + trials = 1; + + // ... get time weighting option for tank updating + + theta = network->option(Options::TIME_WEIGHT); + theta = min(theta, 1.0); + if ( theta > 0.0 ) theta = max(theta, 0.5); + + dosyam << kappa << "\n"; // */ + + kappa = network->option(Options::TEMP_DISC_PARA); + kappa = min(kappa, 1.0); + if (kappa < 0) kappa = 0; // */ + + minErrorNorm = 1000000000; + dl = 1.0; + + // ... set values for convergence limits + setConvergenceLimits(); + + // ... perform Newton iterations + + while ( trials <= trialsLimit ) + { + // ... save current error norm + if (currentTime == 161) + { + dhmaxpast = 0; + maxHeadErrPast = 0; + } + + oldErrorNorm = errorNorm; + + // ... determine which nodes have fixed heads (e.g., PRVs) + + setFixedGradeNodes(); + + // ... re-compute error norm if links change status + + if ( statusChanged ) + { + oldErrorNorm = findErrorNorm(0.0, currentTime, tstep); + lamda = 1.0; + } + statusChanged = false; + + dl = 1.000000000000000; + + // ... find changes in heads and flows + int errorCode = findHeadChanges(); + if (errorCode >= 0) + { + Node* node = network->node(errorCode); + network->msgLog << endl << s_IllConditioned << node->name; + return HydSolver::FAILED_ILL_CONDITIONED; + } + findFlowChanges(); + + if (stepSizing == ARF) + { + double lamda = 1.0; + errorNorm = findErrorNorm(lamda, currentTime, tstep); + updateSolution(lamda); + + if (errorNorm < oldErrorNorm) + { + if (reportTrials) reportTrial(trials, lamda); + converged = hasConverged(); + + // ... if close to convergence then check for any link status changes + + if (converged) //|| errorNorm < ErrorThreshold ) + { + statusChanged = linksChangedStatus(); + } + + // ... check if the current solution can be accepted + + if (converged && !statusChanged) break; + trials++; + } + else + { + // ... find changes in heads and flows + int errorCode = findHeadChanges(); + if (errorCode >= 0) + { + Node* node = network->node(errorCode); + network->msgLog << endl << s_IllConditioned << node->name; + return HydSolver::FAILED_ILL_CONDITIONED; + } + findFlowChanges(); // */ + int counter = 0; + + // ... a while loop which is used to minimum error norm among various error norms + while (minErrorNorm >= oldErrorNorm) + { + dl *= 0.25; + if (dl < 0.001) break; + lamda = findStepSize(trials, currentTime); + updateSolution(lamda); + errorNorm = minErrorNorm; + } + if (reportTrials) reportTrial(trials, lamda); + converged = hasConverged(); + + // ... if close to convergence then check for any link status changes + + if (converged) //|| errorNorm < ErrorThreshold ) + { + statusChanged = linksChangedStatus(); + } + + // ... check if the current solution can be accepted + + if (converged && !statusChanged) break; + trials++; + } + } + + else if (stepSizing == BRF) + { + int counter = 0; + + // ... a while loop which is used to minimum error norm among various error norms + while (minErrorNorm >= oldErrorNorm) + { + dl *= 0.25; + if (dl < 0.001) break; + lamda = findStepSize(trials, currentTime); + updateSolution(lamda); + errorNorm = minErrorNorm; // */ + } + if (reportTrials) reportTrial(trials, lamda); + converged = hasConverged(); + + // ... if close to convergence then check for any link status changes + + if (converged) // || errorNorm < ErrorThreshold ) + { + statusChanged = linksChangedStatus(); + } + + // ... check if the current solution can be accepted + + if (converged && !statusChanged) break; + trials++; + } + else + { + lamda = findStepSize(trials, currentTime); + updateSolution(lamda); + + // ... check for convergence + + if (reportTrials) reportTrial(trials, lamda); + converged = hasConverged(); + + // ... if close to convergence then check for any link status changes + + if (converged) //|| errorNorm < ErrorThreshold ) + { + statusChanged = linksChangedStatus(); + } + + // ... check if the current solution can be accepted + + if (converged && !statusChanged) break; + trials++; + } + } + //if ( reportTrials ) network->msgLog << s_HlossEvals << hLossEvalCount; + if ( trials > trialsLimit ) return HydSolver::FAILED_NO_CONVERGENCE; + return HydSolver::SUCCESSFUL; +} + +//----------------------------------------------------------------------------- + +// Establish error limits for convergence of heads and flows + +void RWCGGASolver::setConvergenceLimits() +{ + // ... maximum trials + trialsLimit = network->option(Options::MAX_TRIALS); + + // ... limit on relative flow change ratio (old EPANET2 criterion) + flowRatioLimit = network->option(Options::RELATIVE_ACCURACY); + + // ... tolerance in satisfying hydraulic balances + headErrLimit = network->option(Options::HEAD_TOLERANCE) / + network->ucf(Units::LENGTH); + flowErrLimit = network->option(Options::FLOW_TOLERANCE) / + network->ucf(Units::FLOW); + + // ... limit on largest flow change + flowChangeLimit = network->option(Options::FLOW_CHANGE_LIMIT) / + network->ucf(Units::FLOW); + + // ... use a default head error limit if need be + if ( flowRatioLimit == 0.0 && headErrLimit == 0.0 && + flowErrLimit == 0.0 && flowChangeLimit == 0.0 ) + { + headErrLimit = 0.005; + } + + // ... convert missing limits to a huge number + if ( flowRatioLimit == 0.0 ) flowRatioLimit = Huge; + if ( headErrLimit == 0.0 ) headErrLimit = Huge; + if ( flowErrLimit == 0.0 ) flowErrLimit = Huge; + if ( flowChangeLimit == 0.0 ) flowChangeLimit = Huge;} + +//----------------------------------------------------------------------------- + +// Adjust fixed grade status of specific nodes. + +void RWCGGASolver::setFixedGradeNodes() +{ + Node* node; + + // ... change fixed grade status for PRV/PSV nodes + + for (Link* link : network->links) + { + // ... check if link is a PRV or PSV + + if ( link->isPRV() ) node = link->toNode; + else if ( link->isPSV() ) node = link->fromNode; + else continue; + + // ... set the fixed grade status of the valve's control node + + if ( link->status == Link::VALVE_ACTIVE ) + { + node->fixedGrade = true; + node->head = link->setting + node->elev; + } + else node->fixedGrade = false; + } + + // ... after time 0, tstep will be non-zero and tank levels + // will be non-fixed if time weighting (theta) is non-zero + + if ( theta > 0.0 && tstep > 0.0 ) + { + for (Node* tankNode : network->nodes) + { + if ( tankNode->type() == Node::TANK ) tankNode->fixedGrade = false; + } + } + +} + +//----------------------------------------------------------------------------- + +// Find changes in nodal heads by solving a linearized system of equations. + +int RWCGGASolver::findHeadChanges() +{ + // ... setup the coeff. matrix of the RWCGGA linearized system + + setMatrixCoeffs(); + + // ... temporarily use the head change array dH[] to store new heads + + double *h = &dH[0]; + + // ... solve the linearized RWCGGA system for new nodal heads + // (matrixSolver returns a negative integer if it runs successfully; + // otherwise it returns the index of the row that caused it to fail.) + + int errorCode = matrixSolver->solve(nodeCount, h); + if ( errorCode >= 0 ) return errorCode; + + // ... save new heads as head changes + + for (int i = 0; i < nodeCount; i++) + { + dH[i] = h[i] - network->node(i)->head; + } + + // ... return a negative number indicating that + // the matrix solver ran successfully + + return -1; +} + +//----------------------------------------------------------------------------- + +// Find the changes in link flows resulting from a set of nodal head changes. + +void RWCGGASolver::findFlowChanges() +{ + for (int i = 0; i < linkCount; i++) + { + // ... get link object and its end node indexes + + dQ[i] = 0.0; + Link* link = network->link(i); + int n1 = link->fromNode->index; + int n2 = link->toNode->index; + + // ... flow change for pressure regulating valves + + if (link->hGrad == 0.0) + { + if (link->isPRV()) dQ[i] = -xQ[n2] - link->flow; + if (link->isPSV()) dQ[i] = xQ[n1] - link->flow; + continue; + } + + if (tstep == 0) // || network->link->type() == Link::VALVE || network->link->type() == Link::PUMP) + { + + // ... apply GGA flow change formula: + + double dh = (link->fromNode->head + dH[n1]) - + (link->toNode->head + dH[n2]); + double dq = (link->hLoss - dh) / link->hGrad; + + // ... special case to prevent negative flow in constant HP pumps + + if (link->isHpPump() && + link->status == Link::LINK_OPEN && + dq > link->flow) dq = link->flow / 2.0; + + // ... save flow change + + dQ[i] = -dq; + } + else + { + // ... apply RWCGGA flow change formula: + + double dh = (link->fromNode->head + dH[n1]) - (link->toNode->head + dH[n2]); + double dhpast = (link->fromNode->pastHead) - (link->toNode->pastHead); + double pastTerms = ((1 - kappa) / kappa) * (link->pastHloss - dhpast); + //double flows = (link->inertialTerm / (kappa * tstep)) * (link->flow - link->pastFlow); + double flows = (link->inertialTerm / (kappa * tstep)) * (link->pastFlow) + link->hGrad * link->flow; + double dq = -(dh - link->hLoss + flows - pastTerms) / (link->hGrad + (link->inertialTerm / (kappa * tstep)))+link->flow; + + // ... special case to prevent negative flow in constant HP pumps + + if (link->isHpPump() && + link->status == Link::LINK_OPEN && + dq > link->flow) dq = link->flow / 2.0; + + // ... save flow change + + dQ[i] = -dq; + } + } +} + +//----------------------------------------------------------------------------- + +// Find how much of the head and flow changes to apply to a new solution. + +double RWCGGASolver::findStepSize(int trials, int currentTime) +{ + // ... find the new error norm at full step size + + double lamda = 1.0; + errorNorm = findErrorNorm(lamda, currentTime, tstep); + + if (stepSizing == RELAXATION && oldErrorNorm < ErrorThreshold) + { + lamda = 0.5; + double errorNorm2 = findErrorNorm(lamda, currentTime, tstep); + if (errorNorm2 < errorNorm) errorNorm = errorNorm2; + else + { + lamda = 1.0; + errorNorm = findErrorNorm(lamda, currentTime, tstep); + } + } + + // ... if called for, implement a lamda search procedure + // to find the best step size lamda to take + + if (stepSizing == ARF || stepSizing == BRF) + { + { + minErrorNorm = 0; + double testError = 1000000; + lambdaNumber = 1 / dl; + Lambda.resize(lambdaNumber, 0); + + memset(&Lambda[0], 0, lambdaNumber*sizeof(double)); + + for (int i = 0; i < lambdaNumber; i++) + { + Lambda[i] += (i + 1)*dl; + + int errorCode = findHeadChanges(); + if (errorCode >= 0) + { + Node* node = network->node(errorCode); + network->msgLog << endl << s_IllConditioned << node->name; + return HydSolver::FAILED_ILL_CONDITIONED; + } + findFlowChanges(); //*/ + errorNorm = findErrorNorm(Lambda[i], currentTime, tstep); + + if (errorNorm < testError) + { + testError = errorNorm; + lamda = Lambda[i]; + updateSolution(Lambda[i]); + } + } + minErrorNorm = testError; + } + return lamda; + } + return lamda; +} + +//----------------------------------------------------------------------------- + +// Compute the error norm associated with a given step size. + +double RWCGGASolver::findErrorNorm(double lamda, int currentTime, double tstep) +{ + + hLossEvalCount++; + return hydBalance.evaluate(lamda, (double*)&dH[0], (double*)&dQ[0], + (double*)&xQ[0], network, currentTime, tstep); +} + +//----------------------------------------------------------------------------- + +// Update heads and flows for a given step size. + +void RWCGGASolver::updateSolution(double lamda) +{ + for (int i = 0; i < nodeCount; i++) network->node(i)->head += lamda * dH[i]; + for (int i = 0; i < linkCount; i++) network->link(i)->flow += lamda * dQ[i]; +} + +//----------------------------------------------------------------------------- + +// Check if current solution meets convergence limits + +bool RWCGGASolver::hasConverged() +{ + return ( hydBalance.maxHeadErr < headErrLimit ) && + ( hydBalance.maxFlowErr < flowErrLimit ) && + ( hydBalance.maxFlowChange < flowChangeLimit ) && + ( hydBalance.totalFlowChange < flowRatioLimit ); +} + +//----------------------------------------------------------------------------- + +void RWCGGASolver::reportTrial(int trials, double lamda) +{ + network->msgLog << endl << endl << s_Trial << trials << ":"; + network->msgLog << endl << s_StepSize << lamda; + network->msgLog << endl << s_TotalError << errorNorm; + + // ... report link with maximum head loss error + + network->msgLog << endl << s_HeadError << + hydBalance.maxHeadErr * network->ucf(Units::LENGTH) << " " << + network->getUnits(Units::LENGTH); + if (hydBalance.maxHeadErrLink >= 0) + { + network->msgLog << s_ForLink << network->link(hydBalance.maxHeadErrLink)->name; + } + + // ... report node with maximum flow balance error + + network->msgLog << endl << s_FlowError << + hydBalance.maxFlowErr * network->ucf(Units::FLOW) << + " " << network->getUnits(Units::FLOW); + if (hydBalance.maxFlowErrNode >= 0) + { + network->msgLog << s_AtNode << network->node(hydBalance.maxFlowErrNode)->name; + } + + // ... report link with largest change in flow rate + + network->msgLog << endl << s_FlowChange << hydBalance.maxFlowChange * + network->ucf(Units::FLOW) << " " << network->getUnits(Units::FLOW); + if ( hydBalance.maxFlowChangeLink >= 0 ) + { + network->msgLog << s_ForLink << network->link(hydBalance.maxFlowChangeLink)->name; + } + + // ... report total link flow change relative to total link flow + + network->msgLog << endl << s_TotFlowChange << hydBalance.totalFlowChange; +} + +//----------------------------------------------------------------------------- + +// Compute the coefficient matrix of the linearized set of equations for heads. + +void RWCGGASolver::setMatrixCoeffs() +{ + + memset(&xQ[0], 0, nodeCount*sizeof(double)); + matrixSolver->reset(); + setLinkCoeffs(); + setNodeCoeffs(); + setValveCoeffs(); +} + +//----------------------------------------------------------------------------- + + +// Compute matrix coefficients for link head loss gradients. + +void RWCGGASolver::setLinkCoeffs() +{ + for (int j = 0; j < linkCount; j++) + { + // ... skip links with zero head gradient + // (e.g. active pressure regulating valves) + + Link* link = network->link(j); + if ( link->hGrad == 0.0 ) continue; + + // ... identify end nodes of link + + Node* node1 = link->fromNode; + Node* node2 = link->toNode; + int n1 = node1->index; + int n2 = node2->index; + + // ... update node flow balances + + xQ[n1] -= link->flow; + xQ[n2] += link->flow; + + // ... a is contribution to coefficient matrix + // b is contribution to right hand side + + if (tstep == 0) + { + double a = 1.0 / link->hGrad; + double b = a * link->hLoss; + + // ... update off-diagonal coeff. of matrix if both start and + // end nodes are not fixed grade + + if (!node1->fixedGrade && !node2->fixedGrade) + { + matrixSolver->addToOffDiag(j, -a); + } + + // ... if start node has fixed grade, then apply a to r.h.s. + // of that node's row; + + if (node1->fixedGrade) + { + matrixSolver->addToRhs(n2, a * node1->head); + } + + // ... otherwise add a to row's diagonal coeff. and + // add b to its r.h.s. + + else + { + matrixSolver->addToDiag(n1, a); + matrixSolver->addToRhs(n1, b); + } + + // ... do the same for the end node, except subtract b from r.h.s + + if (node2->fixedGrade) + { + matrixSolver->addToRhs(n1, a * node2->head); + } + else + { + matrixSolver->addToDiag(n2, a); + matrixSolver->addToRhs(n2, -b); + } + } + + else + { + // a and b are upgraded according to RWC-GGA + + double a = 1.0 / ((link->hGrad) + (link->inertialTerm / (kappa * tstep))); + double b = a * ((link->hLoss) + ((1 - kappa) / kappa) * (link->pastHloss - (node1->pastHead - node2->pastHead)) - (link->inertialTerm / (kappa * tstep)) * (link->pastFlow) - link->hGrad * link->flow) + link->flow; + + // ... update off-diagonal coeff. of matrix if both start and + // end nodes are not fixed grade + + if (!node1->fixedGrade && !node2->fixedGrade) + { + matrixSolver->addToOffDiag(j, -a); + } + + // ... if start node has fixed grade, then apply a to r.h.s. + // of that node's row; + + if (node1->fixedGrade) + { + matrixSolver->addToRhs(n2, a * node1->head); + } + + // ... otherwise add a to row's diagonal coeff. and + // add b to its r.h.s. + + else + { + matrixSolver->addToDiag(n1, a); + matrixSolver->addToRhs(n1, b); + } + + // ... do the same for the end node, except subtract b from r.h.s + + if (node2->fixedGrade) + { + matrixSolver->addToRhs(n1, a * node2->head); + } + else + { + matrixSolver->addToDiag(n2, a); + matrixSolver->addToRhs(n2, -b); + } + } + } +} + +//----------------------------------------------------------------------------- + +// Compute matrix coefficients for dynamic tanks and external node outflows. + +void RWCGGASolver::setNodeCoeffs() +{ + for (int i = 0; i < nodeCount; i++) + { + // ... if node's head not fixed + + Node* node = network->node(i); + if ( !node->fixedGrade ) + { + // ... for dynamic tanks, add area terms to row i + // of the head solution matrix & r.h.s. vector + + if ( node->type() == Node::TANK && theta != 0.0 ) + { + Tank* tank = static_cast(node); + + if (tank->head == tank->pastHead) + { + double a = tank->area / (theta * tstep); + matrixSolver->addToDiag(i, a); + + a = a * tank->pastHead + (1.0 - theta) * tank->pastOutflow / theta; + matrixSolver->addToRhs(i, a); // */ + } + + else + { + double a = (tank->area + ((tank->area - tank->pastArea) /(tank->head - tank->pastHead)) * (tank->head)) / (theta * tstep); + matrixSolver->addToDiag(i, a); + + double b = (tank->pastArea) * tank->pastHead / (theta * tstep) + (1.0 - theta) * tank->pastOutflow / theta; + double c = b + ((tank->area - tank->pastArea) / (tank->head - tank->pastHead)) * (tank->head) / (theta * tstep); + matrixSolver->addToRhs(i, c); // + } + } + + // ... for junctions, add effect of external outflows + + else if ( node->type() == Node::JUNCTION ) + { + // ... update junction's net inflow + xQ[i] -= node->outflow; + matrixSolver->addToDiag(i, node->qGrad); + matrixSolver->addToRhs(i, node->qGrad * node->head); + } + + // ... add node's net inflow to r.h.s. row + matrixSolver->addToRhs(i, (double)xQ[i]); + } + + // ... if node has fixed head, force solution to produce it + + else + { + matrixSolver->setDiag(i, 1.0); + matrixSolver->setRhs(i, node->head); + } + } +} + +//----------------------------------------------------------------------------- + +// Compute matrix coefficients for pressure regulating valves. + +void RWCGGASolver::setValveCoeffs() +{ + for (Link* link : network->links) + { + // ... skip links that are not active pressure regulating valves + + if ( link->hGrad > 0.0 ) continue; + + // ... determine end node indexes of link + + int n1 = link->fromNode->index; + int n2 = link->toNode->index; + + // ... add net inflow of downstream node of a PRV to the + // r.h.s. row of its upstream node + + if ( link->isPRV() ) + { + matrixSolver->addToRhs(n1, (double)xQ[n2]); + } + + // ... add net inflow of upstream node of a PSV to the + // r.h.s. row of its downstream node + + if ( link->isPSV() ) + { + matrixSolver->addToRhs(n2, (double)xQ[n1]); + } + } +} + +//----------------------------------------------------------------------------- + +// Check if any links change status at the current trial solution. + +bool RWCGGASolver::linksChangedStatus() +{ + bool result = false; + for (Link* link : network->links) + { + // ... get head at each end of link + + double h1 = link->fromNode->head; + double h2 = link->toNode->head; + double q = link->flow; + + // ... update link's status + + int oldStatus = link->status; + if ( link->status >= Link::TEMP_CLOSED ) link->status = Link::LINK_OPEN; + link->updateStatus(q, h1, h2); + + // ... check for flow into full or out of empty tanks + + if ( link->status > Link::LINK_CLOSED ) + { + if ( link->fromNode->isClosed(q) || link->toNode->isClosed(-q) ) + { + link->status = Link::TEMP_CLOSED; + link->flow = ZERO_FLOW; + } + } + + //... write status change to message log + + if (oldStatus != link->status) + { + if ( reportTrials ) + { + network->msgLog << endl << link->writeStatusChange(oldStatus); + } + result = true; + } + } + //if ( result && reportTrials ) network->msgLog << endl; + + // --- look for status changes caused by pressure switch controls + if (Control::applyPressureControls(network)) + result = true; + + return result; +} diff --git a/src/Solvers/rwcggasolver.h b/src/Solvers/rwcggasolver.h new file mode 100644 index 0000000..ed76779 --- /dev/null +++ b/src/Solvers/rwcggasolver.h @@ -0,0 +1,89 @@ +/* EPANET 3.1 + * + * Copyright (c) 2016 Open Water Analytics + * Licensed under the terms of the MIT License (see the LICENSE file for details). + * + */ + +//! \file rwcggasolver.h +//! \brief Describes the RWCGGASolver class. + +#ifndef RWCGGASOLVER_H_ +#define RWCGGASOLVER_H_ + +#include "Solvers/hydsolver.h" +#include "Core/hydbalance.h" + +#include + +class HydSolver; + +//! \class RWCGGASolver +//! \brief A hydraulic solver based on RWC Global Gradient Algorithm. + +class RWCGGASolver : public HydSolver +{ + public: + + RWCGGASolver(Network* nw, MatrixSolver* ms); + ~RWCGGASolver(); + int solve(double tstep, int& trials, int currentTime); + double tstep; // time step (sec) + + private: + + int nodeCount; // number of network nodes + int linkCount; // number of network links + int hLossEvalCount; // number of head loss evaluations + int stepSizing; // Newton step sizing method + int lambdaNumber; + + int trialsLimit; // limit on number of trials + bool reportTrials; // report summary of each trial + double headErrLimit; // allowable head error (ft) + double flowErrLimit; // allowable flow error (cfs) + double flowChangeLimit; // allowable flow change (cfs) + double flowRatioLimit; // allowable total flow change / total flow + + double theta; // time weighting constant + double kappa; // temporal discretization parameter + double epsilon; // convergence parameter + double dhmaxpast; + double maxHeadErrPast; + double maxHeadError; + double minErrorNorm; + double dl; + + double errorNorm; // solution error norm + double oldErrorNorm; // previous error norm + HydBalance hydBalance; // hydraulic balance results + + std::vector dH; // head change at each node (ft) + std::vector dQ; // flow change in each link (cfs) + std::vector xQ; // node flow imbalances (cfs) + std::vector Lambda; + + // Functions that assemble linear equation coefficients + void setFixedGradeNodes(); + void setMatrixCoeffs(); + void setLinkCoeffs(); + void setNodeCoeffs(); + void setValveCoeffs(); + void setInertialTerm(); + + + // Functions that update the hydraulic solution + int findHeadChanges(); + void findFlowChanges(); + double findStepSize(int trials, int currentTime); + void updateSolution(double lamda); + + // Functions that check for convergence + void setConvergenceLimits(); + double findErrorNorm(double lamda, int currentTime, double tstep); + bool hasConverged(); + bool linksChangedStatus(); + void reportTrial(int trials, double lamda); +}; + +#endif diff --git a/src/Solvers/sparspaksolver.cpp b/src/Solvers/sparspaksolver.cpp index 8b658a7..4e1be7f 100644 --- a/src/Solvers/sparspaksolver.cpp +++ b/src/Solvers/sparspaksolver.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Solvers/sparspaksolver.h b/src/Solvers/sparspaksolver.h index 6900544..c9f868b 100644 --- a/src/Solvers/sparspaksolver.h +++ b/src/Solvers/sparspaksolver.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Utilities/graph.cpp b/src/Utilities/graph.cpp index 6487f94..303e120 100644 --- a/src/Utilities/graph.cpp +++ b/src/Utilities/graph.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). diff --git a/src/Utilities/graph.h b/src/Utilities/graph.h index 0d7a44f..f91d6ca 100644 --- a/src/Utilities/graph.h +++ b/src/Utilities/graph.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Utilities/mempool.cpp b/src/Utilities/mempool.cpp index 227d5ce..2e46e02 100644 --- a/src/Utilities/mempool.cpp +++ b/src/Utilities/mempool.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Utilities/mempool.h b/src/Utilities/mempool.h index 53a6aad..f542c48 100644 --- a/src/Utilities/mempool.h +++ b/src/Utilities/mempool.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Utilities/segpool.cpp b/src/Utilities/segpool.cpp index b138ca2..88a538d 100644 --- a/src/Utilities/segpool.cpp +++ b/src/Utilities/segpool.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Utilities/segpool.h b/src/Utilities/segpool.h index a8f62dc..d1ff86a 100644 --- a/src/Utilities/segpool.h +++ b/src/Utilities/segpool.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Utilities/utilities.cpp b/src/Utilities/utilities.cpp index 0687d66..f6de68a 100644 --- a/src/Utilities/utilities.cpp +++ b/src/Utilities/utilities.cpp @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). diff --git a/src/Utilities/utilities.h b/src/Utilities/utilities.h index 0d02e20..de0525c 100644 --- a/src/Utilities/utilities.h +++ b/src/Utilities/utilities.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Distributed under the MIT License (see the LICENSE file for details). diff --git a/src/epanet3.h b/src/epanet3.h index 797cef2..1faff17 100644 --- a/src/epanet3.h +++ b/src/epanet3.h @@ -1,4 +1,4 @@ -/* EPANET 3 +/* EPANET 3.1 * * Copyright (c) 2016 Open Water Analytics * Licensed under the terms of the MIT License (see the LICENSE file for details). @@ -96,15 +96,16 @@ enum NodeTypes { EN_TANK}; //2 enum LinkTypes { - EN_CVPIPE, //0 - EN_PIPE, //1 - EN_PUMP, //2 - EN_PRV, //3 - EN_PSV, //4 - EN_PBV, //5 - EN_FCV, //6 - EN_TCV, //7 - EN_GPV}; //8 + EN_CVPIPE, //0 + EN_PIPE, //1 + EN_PUMP, //2 + EN_PRV, //3 + EN_PSV, //4 + EN_PBV, //5 + EN_FCV, //6 + EN_TCV, //7 + EN_GPV, //8 + EN_CCV}; //9 enum QualModelTypes { EN_NONE, //0 @@ -206,6 +207,7 @@ int EN_getLinkId(int, char *, EN_Project); int EN_getLinkType(int, int *, EN_Project); int EN_getLinkNodes(int, int *, int *, EN_Project); int EN_getLinkValue(int, int, double *, EN_Project); +int EN_setLinkValue(int, int, double, EN_Project); //================================================================================== /* TO BE ADDED