diff --git a/.gitignore b/.gitignore old mode 100644 new mode 100755 diff --git a/README.md b/README.md old mode 100644 new mode 100755 diff --git a/build/MSVS/LemonTigerJ.sln b/build/MSVS/LemonTigerJ.sln old mode 100644 new mode 100755 diff --git a/build/MSVS/LemonTigerJ.vcxproj b/build/MSVS/LemonTigerJ.vcxproj old mode 100644 new mode 100755 diff --git a/build/Xcode/epanet/epanet.xcodeproj/project.pbxproj b/build/Xcode/epanet/epanet.xcodeproj/project.pbxproj old mode 100644 new mode 100755 index 30d1018..2f904ac --- a/build/Xcode/epanet/epanet.xcodeproj/project.pbxproj +++ b/build/Xcode/epanet/epanet.xcodeproj/project.pbxproj @@ -40,7 +40,7 @@ 22322FA41068369500641384 /* report.c in Sources */ = {isa = PBXBuildFile; fileRef = 22322F7E1068369500641384 /* report.c */; }; 22322FA51068369500641384 /* rules.c in Sources */ = {isa = PBXBuildFile; fileRef = 22322F7F1068369500641384 /* rules.c */; }; 22322FA61068369500641384 /* smatrix.c in Sources */ = {isa = PBXBuildFile; fileRef = 22322F801068369500641384 /* smatrix.c */; }; - 22322FAA106836BC00641384 /* epanet2.h in Headers */ = {isa = PBXBuildFile; fileRef = 22322FA9106836B000641384 /* epanet2.h */; }; + 22322FAA106836BC00641384 /* epanet2.h in Headers */ = {isa = PBXBuildFile; fileRef = 22322FA9106836B000641384 /* epanet2.h */; settings = {ATTRIBUTES = (Public, ); }; }; 2255753F17551234009946B1 /* epanet.c in Sources */ = {isa = PBXBuildFile; fileRef = 22322F711068369500641384 /* epanet.c */; }; 2255754017551234009946B1 /* hash.c in Sources */ = {isa = PBXBuildFile; fileRef = 22322F731068369500641384 /* hash.c */; }; 2255754117551234009946B1 /* hydraul.c in Sources */ = {isa = PBXBuildFile; fileRef = 22322F751068369500641384 /* hydraul.c */; }; @@ -54,6 +54,7 @@ 2255754917551234009946B1 /* report.c in Sources */ = {isa = PBXBuildFile; fileRef = 22322F7E1068369500641384 /* report.c */; }; 2255754A17551234009946B1 /* rules.c in Sources */ = {isa = PBXBuildFile; fileRef = 22322F7F1068369500641384 /* rules.c */; }; 2255754B17551234009946B1 /* smatrix.c in Sources */ = {isa = PBXBuildFile; fileRef = 22322F801068369500641384 /* smatrix.c */; }; + 226537E0179EDEEB00258C60 /* epanet2.h in Headers */ = {isa = PBXBuildFile; fileRef = 22322FA9106836B000641384 /* epanet2.h */; settings = {ATTRIBUTES = (Public, ); }; }; /* End PBXBuildFile section */ /* Begin PBXContainerItemProxy section */ @@ -180,6 +181,7 @@ isa = PBXHeadersBuildPhase; buildActionMask = 2147483647; files = ( + 226537E0179EDEEB00258C60 /* epanet2.h in Headers */, ); runOnlyForDeploymentPostprocessing = 0; }; @@ -258,7 +260,7 @@ 08FB7793FE84155DC02AAC07 /* Project object */ = { isa = PBXProject; attributes = { - LastUpgradeCheck = 0460; + LastUpgradeCheck = 0500; }; buildConfigurationList = 1DEB914E08733D8E0010E9CD /* Build configuration list for PBXProject "epanet" */; compatibilityVersion = "Xcode 3.2"; @@ -390,16 +392,20 @@ 1DEB914F08733D8E0010E9CD /* Debug */ = { isa = XCBuildConfiguration; buildSettings = { - ARCHS = "$(ARCHS_STANDARD_32_64_BIT)"; + CLANG_WARN_BOOL_CONVERSION = YES; CLANG_WARN_CONSTANT_CONVERSION = YES; + CLANG_WARN_EMPTY_BODY = YES; CLANG_WARN_ENUM_CONVERSION = YES; CLANG_WARN_INT_CONVERSION = YES; CLANG_WARN__DUPLICATE_METHOD_MATCH = YES; GCC_C_LANGUAGE_STANDARD = gnu99; GCC_OPTIMIZATION_LEVEL = 0; GCC_VERSION = com.apple.compilers.llvm.clang.1_0; + GCC_WARN_64_TO_32_BIT_CONVERSION = YES; GCC_WARN_ABOUT_RETURN_TYPE = YES; + GCC_WARN_UNDECLARED_SELECTOR = YES; GCC_WARN_UNINITIALIZED_AUTOS = YES; + GCC_WARN_UNUSED_FUNCTION = YES; GCC_WARN_UNUSED_VARIABLE = YES; HEADER_SEARCH_PATHS = macinclude; ONLY_ACTIVE_ARCH = YES; @@ -410,15 +416,19 @@ 1DEB915008733D8E0010E9CD /* Release */ = { isa = XCBuildConfiguration; buildSettings = { - ARCHS = "$(ARCHS_STANDARD_32_64_BIT)"; + CLANG_WARN_BOOL_CONVERSION = YES; CLANG_WARN_CONSTANT_CONVERSION = YES; + CLANG_WARN_EMPTY_BODY = YES; CLANG_WARN_ENUM_CONVERSION = YES; CLANG_WARN_INT_CONVERSION = YES; CLANG_WARN__DUPLICATE_METHOD_MATCH = YES; GCC_C_LANGUAGE_STANDARD = gnu99; GCC_VERSION = com.apple.compilers.llvm.clang.1_0; + GCC_WARN_64_TO_32_BIT_CONVERSION = YES; GCC_WARN_ABOUT_RETURN_TYPE = YES; + GCC_WARN_UNDECLARED_SELECTOR = YES; GCC_WARN_UNINITIALIZED_AUTOS = YES; + GCC_WARN_UNUSED_FUNCTION = YES; GCC_WARN_UNUSED_VARIABLE = YES; HEADER_SEARCH_PATHS = macinclude; SDKROOT = ""; @@ -461,7 +471,6 @@ isa = XCBuildConfiguration; buildSettings = { ALWAYS_SEARCH_USER_PATHS = NO; - ARCHS = "$(ARCHS_STANDARD_64_BIT)"; CLANG_CXX_LANGUAGE_STANDARD = "gnu++0x"; CLANG_CXX_LIBRARY = "libc++"; CLANG_WARN_EMPTY_BODY = YES; @@ -482,7 +491,6 @@ isa = XCBuildConfiguration; buildSettings = { ALWAYS_SEARCH_USER_PATHS = NO; - ARCHS = "$(ARCHS_STANDARD_64_BIT)"; CLANG_CXX_LANGUAGE_STANDARD = "gnu++0x"; CLANG_CXX_LIBRARY = "libc++"; CLANG_WARN_EMPTY_BODY = YES; diff --git a/build/Xcode/epanet/epanet.xcodeproj/project.xcworkspace/contents.xcworkspacedata b/build/Xcode/epanet/epanet.xcodeproj/project.xcworkspace/contents.xcworkspacedata old mode 100644 new mode 100755 diff --git a/include/epanet2.h b/include/epanet2.h old mode 100644 new mode 100755 index d851a43..0543a15 --- a/include/epanet2.h +++ b/include/epanet2.h @@ -1,7 +1,7 @@ /* ******************************************************************* - TOOLKIT.H - Prototypes for EPANET Functions Exported to DLL Toolkit + EPANET2.H - Prototypes for EPANET Functions Exported to DLL Toolkit VERSION: 2.00 DATE: 5/8/00 @@ -9,8 +9,8 @@ 3/1/01 8/15/07 (2.00.11) 2/14/08 (2.00.12) - AUTHOR: L. Rossman - US EPA - NRMRL + AUTHORS: L. Rossman - US EPA - NRMRL + OpenWaterAnalytics members: see git stats for contributors ******************************************************************* */ @@ -177,96 +177,81 @@ #if defined(__cplusplus) extern "C" { #endif - int DLLEXPORT ENepanet(char *, char *, char *, void (*) (char *)); + int DLLEXPORT ENepanet(char *inpFile, char *rptFile, char *binOutFile, void (*callback) (char *)); - int DLLEXPORT ENopen(char *, char *, char *); - int DLLEXPORT ENsaveinpfile(char *); - int DLLEXPORT ENclose(void); + int DLLEXPORT ENopen(char *inpFile, char *rptFile, char *binOutFile); + int DLLEXPORT ENsaveinpfile(char *filename); + int DLLEXPORT ENclose(); - int DLLEXPORT ENsolveH(void); - int DLLEXPORT ENsaveH(void); - int DLLEXPORT ENopenH(void); - int DLLEXPORT ENinitH(int); - int DLLEXPORT ENrunH(long *); - int DLLEXPORT ENnextH(long *tstep); - int DLLEXPORT ENcloseH(void); - int DLLEXPORT ENsavehydfile(char *); - int DLLEXPORT ENusehydfile(char *); + int DLLEXPORT ENsolveH(); + int DLLEXPORT ENsaveH(); + int DLLEXPORT ENopenH(); + int DLLEXPORT ENinitH(int initFlag); + int DLLEXPORT ENrunH(long *currentTime); + int DLLEXPORT ENnextH(long *tStep); + int DLLEXPORT ENcloseH(); + int DLLEXPORT ENsavehydfile(char *filename); + int DLLEXPORT ENusehydfile(char *filename); - int DLLEXPORT ENsolveQ(void); - int DLLEXPORT ENopenQ(void); - int DLLEXPORT ENinitQ(int); - int DLLEXPORT ENrunQ(long *); - int DLLEXPORT ENnextQ(long *); - int DLLEXPORT ENstepQ(long *); - int DLLEXPORT ENcloseQ(void); + int DLLEXPORT ENsolveQ(); + int DLLEXPORT ENopenQ(); + int DLLEXPORT ENinitQ(int saveFlag); + int DLLEXPORT ENrunQ(long *currentTime); + int DLLEXPORT ENnextQ(long *tStep); + int DLLEXPORT ENstepQ(long *timeLeft); + int DLLEXPORT ENcloseQ(); - int DLLEXPORT ENwriteline(char *); - int DLLEXPORT ENreport(void); - int DLLEXPORT ENresetreport(void); - int DLLEXPORT ENsetreport(char *); + int DLLEXPORT ENwriteline(char *line); + int DLLEXPORT ENreport(); + int DLLEXPORT ENresetreport(); + int DLLEXPORT ENsetreport(char *reportFormat); - int DLLEXPORT ENgetcontrol(int, int *, int *, EN_API_FLOAT_TYPE *, - int *, EN_API_FLOAT_TYPE *); - int DLLEXPORT ENgetcount(int, int *); - int DLLEXPORT ENgetoption(int, EN_API_FLOAT_TYPE *); - int DLLEXPORT ENgettimeparam(int, long *); - int DLLEXPORT ENgetflowunits(int *); - int DLLEXPORT ENgetpatternindex(char *, int *); - int DLLEXPORT ENgetpatternid(int, char *); - int DLLEXPORT ENgetpatternlen(int, int *); - int DLLEXPORT ENgetpatternvalue(int, int, EN_API_FLOAT_TYPE *); + int DLLEXPORT ENgetcontrol(int controlIndex, int *controlType, int *linkIdx, EN_API_FLOAT_TYPE *setting, int *nodeIdx, EN_API_FLOAT_TYPE *level); + int DLLEXPORT ENgetcount(int code, int *count); + int DLLEXPORT ENgetoption(int code, EN_API_FLOAT_TYPE *value); + int DLLEXPORT ENgettimeparam(int code, long *value); + int DLLEXPORT ENgetflowunits(int *code); + int DLLEXPORT ENgetpatternindex(char *id, int *index); + int DLLEXPORT ENgetpatternid(int index, char *id); + int DLLEXPORT ENgetpatternlen(int index, int *len); + int DLLEXPORT ENgetpatternvalue(int index, int period, EN_API_FLOAT_TYPE *value); int DLLEXPORT ENgetaveragepatternvalue(int index, EN_API_FLOAT_TYPE *value); - int DLLEXPORT ENgetqualtype(int *, int *); - int DLLEXPORT ENgeterror(int, char *, int); - int DLLEXPORT ENgetstatistic(int code, int* value); + int DLLEXPORT ENgetqualtype(int *qualcode, int *tracenode); + int DLLEXPORT ENgeterror(int errcode, char *errmsg, int maxLen); + int DLLEXPORT ENgetstatistic(int code, EN_API_FLOAT_TYPE* value); - int DLLEXPORT ENgetnodeindex(char *, int *); - int DLLEXPORT ENgetnodeid(int, char *); - int DLLEXPORT ENgetnodetype(int, int *); - int DLLEXPORT ENgetnodevalue(int, int, EN_API_FLOAT_TYPE *); - int DLLEXPORT ENgetcoord(int , EN_API_FLOAT_TYPE *, EN_API_FLOAT_TYPE *); + int DLLEXPORT ENgetnodeindex(char *id, int *index); + int DLLEXPORT ENgetnodeid(int index, char *id); + int DLLEXPORT ENgetnodetype(int index, int *code); + int DLLEXPORT ENgetnodevalue(int index, int code, EN_API_FLOAT_TYPE *value); + int DLLEXPORT ENgetcoord(int index, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y); - int DLLEXPORT ENgetnumdemands(int, int *); - int DLLEXPORT ENgetbasedemand(int, int, EN_API_FLOAT_TYPE *); - int DLLEXPORT ENgetdemandpattern(int, int, int *); + int DLLEXPORT ENgetnumdemands(int nodeIndex, int *numDemands); + int DLLEXPORT ENgetbasedemand(int nodeIndex, int demandIdx, EN_API_FLOAT_TYPE *baseDemand); + int DLLEXPORT ENgetdemandpattern(int nodeIndex, int demandIdx, int *pattIdx); - int DLLEXPORT ENgetlinkindex(char *, int *); - int DLLEXPORT ENgetlinkid(int, char *); - int DLLEXPORT ENgetlinktype(int, int *); - int DLLEXPORT ENgetlinknodes(int, int *, int *); - int DLLEXPORT ENgetlinkvalue(int, int, EN_API_FLOAT_TYPE *); + int DLLEXPORT ENgetlinkindex(char *id, int *index); + int DLLEXPORT ENgetlinkid(int index, char *id); + int DLLEXPORT ENgetlinktype(int index, int *code); + int DLLEXPORT ENgetlinknodes(int index, int *node1, int *node2); + int DLLEXPORT ENgetlinkvalue(int index, int code, EN_API_FLOAT_TYPE *value); - int DLLEXPORT ENgetcurve(int curveIndex, int *nValues, EN_API_FLOAT_TYPE **xValues, EN_API_FLOAT_TYPE **yValues); + int DLLEXPORT ENgetcurve(int curveIndex, char* id, int *nValues, EN_API_FLOAT_TYPE **xValues, EN_API_FLOAT_TYPE **yValues); - int DLLEXPORT ENgetversion(int *); + int DLLEXPORT ENgetversion(int *version); - int DLLEXPORT ENsetcontrol(int, int, int, EN_API_FLOAT_TYPE, int, EN_API_FLOAT_TYPE); - int DLLEXPORT ENsetnodevalue(int, int, EN_API_FLOAT_TYPE); - int DLLEXPORT ENsetlinkvalue(int, int, EN_API_FLOAT_TYPE); - int DLLEXPORT ENaddpattern(char *); - int DLLEXPORT ENsetpattern(int, EN_API_FLOAT_TYPE *, int); - int DLLEXPORT ENsetpatternvalue(int, int, EN_API_FLOAT_TYPE); - int DLLEXPORT ENsettimeparam(int, long); - int DLLEXPORT ENsetoption(int, EN_API_FLOAT_TYPE); - int DLLEXPORT ENsetstatusreport(int); + int DLLEXPORT ENsetcontrol(int cindex, int ctype, int lindex, EN_API_FLOAT_TYPE setting, int nindex, EN_API_FLOAT_TYPE level); + int DLLEXPORT ENsetnodevalue(int index, int code, EN_API_FLOAT_TYPE v); + int DLLEXPORT ENsetlinkvalue(int index, int code, EN_API_FLOAT_TYPE v); + int DLLEXPORT ENaddpattern(char *id); + int DLLEXPORT ENsetpattern(int index, EN_API_FLOAT_TYPE *f, int len); + int DLLEXPORT ENsetpatternvalue(int index, int period, EN_API_FLOAT_TYPE value); + int DLLEXPORT ENsettimeparam(int code, long value); + int DLLEXPORT ENsetoption(int code, EN_API_FLOAT_TYPE v); + int DLLEXPORT ENsetstatusreport(int code); int DLLEXPORT ENsetqualtype(int qualcode, char *chemname, char *chemunits, char *tracenode); - - //LemonTiger functions - /* See testLT.c for a LemonTiger test */ - - //LT equivalent to ENopenH() + ENopenQ() + ENinitH() + ENinitQ() - int DLLEXPORT ENopeninitHQ(); - - //LT equivalent to ENrunQ() + ENnextQ(); - int DLLEXPORT ENrunnextHQ(long*, long*); - - //LT equivalent to ENrunQ() + ENstepQ(); - int DLLEXPORT ENrunstepHQ(long*, long*); - - //LT equivalent to ENcloseH() + ENcloseQ(); - int DLLEXPORT ENcloseHQ(); + int DLLEXPORT ENsetbasedemand(int nodeIndex, int demandIdx, EN_API_FLOAT_TYPE baseDemand); #if defined(__cplusplus) } diff --git a/src/epanet.c b/src/epanet.c index 1a10546..2c4aa09 100755 --- a/src/epanet.c +++ b/src/epanet.c @@ -1248,7 +1248,7 @@ int DLLEXPORT ENgeterror(int errcode, char *errmsg, int n) else return(0); } -int DLLEXPORT ENgetstatistic(int code, int* value) +int DLLEXPORT ENgetstatistic(int code, EN_API_FLOAT_TYPE* value) /*---------------------------------------------------------------- ** Input: code = type of simulation statistic to retrieve ** Output: value = value of requested statistic @@ -1259,10 +1259,10 @@ int DLLEXPORT ENgetstatistic(int code, int* value) { switch (code) { case EN_ITERATIONS: - *value = _iterations; + *value = (EN_API_FLOAT_TYPE)_iterations; break; case EN_RELATIVEERROR: - *value = _relativeError; + *value = (EN_API_FLOAT_TYPE)_relativeError; break; default: break; @@ -1443,19 +1443,19 @@ int DLLEXPORT ENgetnodevalue(int index, int code, EN_API_FLOAT_TYPE *value) break; //(2.00.11 - LR) case EN_DEMAND: - v = D[index]*Ucf[FLOW]; + v = NodeDemand[index]*Ucf[FLOW]; break; case EN_HEAD: - v = H[index]*Ucf[HEAD]; + v = NodeHead[index]*Ucf[HEAD]; break; case EN_PRESSURE: - v = (H[index] - Node[index].El)*Ucf[PRESSURE]; + v = (NodeHead[index] - Node[index].El)*Ucf[PRESSURE]; break; case EN_QUALITY: - v = C[index]*Ucf[QUALITY]; + v = NodeQual[index]*Ucf[QUALITY]; break; /*** New parameters added for retrieval begins here ***/ //(2.00.12 - LR) @@ -1492,7 +1492,7 @@ int DLLEXPORT ENgetnodevalue(int index, int code, EN_API_FLOAT_TYPE *value) v = (Tank[index-Njuncs].Hmin - Node[index].El) * Ucf[ELEV]; } break; - + case EN_MAXLEVEL: v = 0.0; if ( index > Njuncs ) @@ -1518,7 +1518,7 @@ int DLLEXPORT ENgetnodevalue(int index, int code, EN_API_FLOAT_TYPE *value) case EN_TANKVOLUME: if (index <= Njuncs) return(251); - v = tankvolume(index-Njuncs, H[index])*Ucf[VOLUME]; + v = tankvolume(index-Njuncs, NodeHead[index])*Ucf[VOLUME]; break; default: return(251); @@ -1685,7 +1685,7 @@ int DLLEXPORT ENgetlinkvalue(int index, int code, EN_API_FLOAT_TYPE *value) case EN_FLOW: /*** Updated 10/25/00 ***/ - if (S[index] <= CLOSED) v = 0.0; + if (LinkStatus[index] <= CLOSED) v = 0.0; else v = Q[index]*Ucf[FLOW]; break; @@ -1693,7 +1693,7 @@ int DLLEXPORT ENgetlinkvalue(int index, int code, EN_API_FLOAT_TYPE *value) if (Link[index].Type == PUMP) v = 0.0; /*** Updated 11/19/01 ***/ - else if (S[index] <= CLOSED) v = 0.0; + else if (LinkStatus[index] <= CLOSED) v = 0.0; else { @@ -1706,26 +1706,31 @@ int DLLEXPORT ENgetlinkvalue(int index, int code, EN_API_FLOAT_TYPE *value) case EN_HEADLOSS: /*** Updated 11/19/01 ***/ - if (S[index] <= CLOSED) v = 0.0; + if (LinkStatus[index] <= CLOSED) v = 0.0; else { - h = H[Link[index].N1] - H[Link[index].N2]; + h = NodeHead[Link[index].N1] - NodeHead[Link[index].N2]; if (Link[index].Type != PUMP) h = ABS(h); v = h*Ucf[HEADLOSS]; } break; case EN_STATUS: - if (S[index] <= CLOSED) v = 0.0; + if (LinkStatus[index] <= CLOSED) v = 0.0; else v = 1.0; break; case EN_SETTING: - if (Link[index].Type == PIPE || Link[index].Type == CV) + if (Link[index].Type == PIPE || Link[index].Type == CV) { return(ENgetlinkvalue(index, EN_ROUGHNESS, value)); - if (K[index] == MISSING) v = 0.0; - else v = K[index]; + } + if (LinkSetting[index] == MISSING) { + v = 0.0; + } + else { + v = LinkSetting[index]; + } switch (Link[index].Type) { case PRV: @@ -1750,7 +1755,7 @@ int DLLEXPORT ENgetlinkvalue(int index, int code, EN_API_FLOAT_TYPE *value) } -int DLLEXPORT ENgetcurve(int curveIndex, int *nValues, EN_API_FLOAT_TYPE **xValues, EN_API_FLOAT_TYPE **yValues) +int DLLEXPORT ENgetcurve(int curveIndex, char* id, int *nValues, EN_API_FLOAT_TYPE **xValues, EN_API_FLOAT_TYPE **yValues) /*---------------------------------------------------------------- ** Input: curveIndex = curve index ** Output: *nValues = number of points on curve @@ -1776,6 +1781,7 @@ int DLLEXPORT ENgetcurve(int curveIndex, int *nValues, EN_API_FLOAT_TYPE **xVal pointY[iPoint] = (EN_API_FLOAT_TYPE)y; } + strncpy(id, curve.ID, MAXID); *nValues = nPoints; *xValues = pointX; *yValues = pointY; @@ -1910,7 +1916,7 @@ int DLLEXPORT ENsetnodevalue(int index, int code, EN_API_FLOAT_TYPE v) Tank[j].Hmin += value; Tank[j].Hmax += value; Node[index].El += value; - H[index] += value; + NodeHead[index] += value; } break; @@ -1967,7 +1973,9 @@ int DLLEXPORT ENsetnodevalue(int index, int code, EN_API_FLOAT_TYPE v) source->Pat = 0; Node[index].S = source; } - if (code == EN_SOURCEQUAL) source->C0 = value; + if (code == EN_SOURCEQUAL) { + source->C0 = value; + } else if (code == EN_SOURCEPAT) { j = ROUND(value); @@ -1991,7 +1999,7 @@ int DLLEXPORT ENsetnodevalue(int index, int code, EN_API_FLOAT_TYPE v) Tank[j].Hmin = Tank[j].H0; Tank[j].Hmax = Tank[j].H0; Node[index].El = Tank[j].H0; - H[index] = Tank[j].H0; + NodeHead[index] = Tank[j].H0; } else { @@ -2000,7 +2008,9 @@ int DLLEXPORT ENsetnodevalue(int index, int code, EN_API_FLOAT_TYPE v) || value < Tank[j].Hmin) return(202); Tank[j].H0 = value; Tank[j].V0 = tankvolume(j, Tank[j].H0); - H[index] = Tank[j].H0; + // Resetting Volume in addition to initial volume + Tank[j].V = Tank[j].V0; + NodeHead[index] = Tank[j].H0; } break; @@ -2155,7 +2165,7 @@ int DLLEXPORT ENsetlinkvalue(int index, int code, EN_API_FLOAT_TYPE v) if (code == EN_INITSTATUS) setlinkstatus(index, s, &Link[index].Stat, &Link[index].Kc); else - setlinkstatus(index, s, &S[index], &K[index]); + setlinkstatus(index, s, &LinkStatus[index], &LinkSetting[index]); break; case EN_INITSETTING: @@ -2182,7 +2192,7 @@ int DLLEXPORT ENsetlinkvalue(int index, int code, EN_API_FLOAT_TYPE v) if (code == EN_INITSETTING) setlinksetting(index, value, &Link[index].Stat, &Link[index].Kc); else - setlinksetting(index, value, &S[index], &K[index]); + setlinksetting(index, value, &LinkStatus[index], &LinkSetting[index]); } break; @@ -2750,13 +2760,13 @@ void initpointers() **---------------------------------------------------------------- */ { - D = NULL; - C = NULL; - H = NULL; + NodeDemand = NULL; + NodeQual = NULL; + NodeHead = NULL; Q = NULL; - R = NULL; - S = NULL; - K = NULL; + PipeRateCoeff = NULL; + LinkStatus = NULL; + LinkSetting = NULL; OldStat = NULL; Node = NULL; @@ -2819,13 +2829,13 @@ int allocdata() { n = MaxNodes + 1; Node = (Snode *) calloc(n, sizeof(Snode)); - D = (double *) calloc(n, sizeof(double)); - C = (double *) calloc(n, sizeof(double)); - H = (double *) calloc(n, sizeof(double)); + NodeDemand = (double *) calloc(n, sizeof(double)); + NodeQual = (double *) calloc(n, sizeof(double)); + NodeHead = (double *) calloc(n, sizeof(double)); ERRCODE(MEMCHECK(Node)); - ERRCODE(MEMCHECK(D)); - ERRCODE(MEMCHECK(C)); - ERRCODE(MEMCHECK(H)); + ERRCODE(MEMCHECK(NodeDemand)); + ERRCODE(MEMCHECK(NodeQual)); + ERRCODE(MEMCHECK(NodeHead)); } /* Allocate memory for network links */ @@ -2834,12 +2844,12 @@ int allocdata() n = MaxLinks + 1; Link = (Slink *) calloc(n, sizeof(Slink)); Q = (double *) calloc(n, sizeof(double)); - K = (double *) calloc(n, sizeof(double)); - S = (char *) calloc(n, sizeof(char)); + LinkSetting = (double *) calloc(n, sizeof(double)); + LinkStatus = (char *) calloc(n, sizeof(char)); ERRCODE(MEMCHECK(Link)); ERRCODE(MEMCHECK(Q)); - ERRCODE(MEMCHECK(K)); - ERRCODE(MEMCHECK(S)); + ERRCODE(MEMCHECK(LinkSetting)); + ERRCODE(MEMCHECK(LinkStatus)); } /* Allocate memory for tanks, sources, pumps, valves, */ @@ -2950,12 +2960,12 @@ void freedata() Psource source; /* Free memory for computed results */ - free(D); - free(C); - free(H); + free(NodeDemand); + free(NodeQual); + free(NodeHead); free(Q); - free(K); - free(S); + free(LinkSetting); + free(LinkStatus); /* Free memory for node data */ if (Node != NULL) @@ -3285,6 +3295,22 @@ int DLLEXPORT ENgetbasedemand(int nodeIndex, int demandIdx, EN_API_FLOAT_TYPE * } return 0; } + +int DLLEXPORT ENsetbasedemand(int nodeIndex, int demandIdx, EN_API_FLOAT_TYPE baseDemand) +{ + Pdemand d; + int n=1; + /* Check for valid arguments */ + if (!Openflag) return(102); + if (nodeIndex <= 0 || nodeIndex > Nnodes) return(203); + if (nodeIndex <= Njuncs) { + for(d=Node[nodeIndex].D; nnext) n++; + if(n!=demandIdx) return(253); + d->Base = baseDemand/Ucf[FLOW]; + } + return 0; +} + int DLLEXPORT ENgetdemandpattern(int nodeIndex, int demandIdx, int *pattIdx) { Pdemand d; diff --git a/src/hydraul.c b/src/hydraul.c index bdf0213..e9bcbb7 100755 --- a/src/hydraul.c +++ b/src/hydraul.c @@ -123,10 +123,10 @@ void inithyd(int initflag) for (i=1; i<=Ntanks; i++) { Tank[i].V = Tank[i].V0; - H[Tank[i].Node] = Tank[i].H0; + NodeHead[Tank[i].Node] = Tank[i].H0; /*** Updated 10/25/00 ***/ - D[Tank[i].Node] = 0.0; + NodeDemand[Tank[i].Node] = 0.0; OldStat[Nlinks+i] = TEMPCLOSED; } @@ -140,24 +140,24 @@ void inithyd(int initflag) for (i=1; i<=Nlinks; i++) { /* Initialize status and setting */ - S[i] = Link[i].Stat; - K[i] = Link[i].Kc; + LinkStatus[i] = Link[i].Stat; + LinkSetting[i] = Link[i].Kc; /* Start active control valves in ACTIVE position */ //(2.00.11 - LR) if ( (Link[i].Type == PRV || Link[i].Type == PSV || Link[i].Type == FCV) //(2.00.11 - LR) && (Link[i].Kc != MISSING) - ) S[i] = ACTIVE; //(2.00.11 - LR) + ) LinkStatus[i] = ACTIVE; //(2.00.11 - LR) /*** Updated 3/1/01 ***/ /* Initialize flows if necessary */ - if (S[i] <= CLOSED) Q[i] = QZERO; + if (LinkStatus[i] <= CLOSED) Q[i] = QZERO; else if (ABS(Q[i]) <= QZERO || initflag > 0) - initlinkflow(i, S[i], K[i]); + initlinkflow(i, LinkStatus[i], LinkSetting[i]); /* Save initial status */ - OldStat[i] = S[i]; + OldStat[i] = LinkStatus[i]; } /* Reset pump energy usage */ @@ -205,7 +205,7 @@ int runhyd(long *t) if (Statflag) writehydstat(iter,relerr); /* solution info */ - _relativeError = (int)relerr; + _relativeError = relerr; _iterations = iter; /*** Updated 3/1/01 ***/ @@ -375,7 +375,7 @@ void setlinkflow(int k, double dh) /* use approx. inverse of formula. */ if (Formflag == DW) { - x = -log(K[k]/3.7/Link[k].Diam); + x = -log(LinkSetting[k]/3.7/Link[k].Diam); y = sqrt(ABS(dh)/Link[k].R/1.32547); Q[k] = x*y; } @@ -402,17 +402,17 @@ void setlinkflow(int k, double dh) /* For custom pump curve, interpolate from curve */ if (Pump[p].Ptype == CUSTOM) { - dh = -dh*Ucf[HEAD]/SQR(K[k]); + dh = -dh*Ucf[HEAD]/SQR(LinkSetting[k]); i = Pump[p].Hcurve; Q[k] = interp(Curve[i].Npts,Curve[i].Y,Curve[i].X, - dh)*K[k]/Ucf[FLOW]; + dh)*LinkSetting[k]/Ucf[FLOW]; } /* Otherwise use inverse of power curve */ else { - h0 = -SQR(K[k])*Pump[p].H0; - x = pow(K[k],2.0-Pump[p].N); + h0 = -SQR(LinkSetting[k])*Pump[p].H0; + x = pow(LinkSetting[k],2.0-Pump[p].N); x = ABS(h0-dh)/(Pump[p].R*x), y = 1.0/Pump[p].N; Q[k] = pow(x,y); @@ -606,7 +606,7 @@ void demands() if (djunc > 0.0) Dsystem += djunc; sum += djunc; } - D[i] = sum; + NodeDemand[i] = sum; } /* Update head at fixed grade nodes with time patterns. */ @@ -619,7 +619,7 @@ void demands() { k = p % (long) Pattern[j].Length; i = Tank[n].Node; - H[i] = Node[i].El*Pattern[j].F[k]; + NodeHead[i] = Node[i].El*Pattern[j].F[k]; } } } @@ -632,7 +632,7 @@ void demands() { i = Pump[n].Link; k = p % (long) Pattern[j].Length; - setlinksetting(i, Pattern[j].F[k], &S[i], &K[i]); + setlinksetting(i, Pattern[j].F[k], &LinkStatus[i], &LinkSetting[i]); } } } /* End of demands */ @@ -664,8 +664,8 @@ int controls() /* Link is controlled by tank level */ if ((n = Control[i].Node) > 0 && n > Njuncs) { - h = H[n]; - vplus = ABS(D[n]); + h = NodeHead[n]; + vplus = ABS(NodeDemand[n]); v1 = tankvolume(n-Njuncs,h); v2 = tankvolume(n-Njuncs,Control[i].Grade); if (Control[i].Type == LOWLEVEL && v1 <= v2 + vplus) @@ -689,16 +689,16 @@ int controls() /* Update link status & pump speed or valve setting */ if (reset == 1) { - if (S[k] <= CLOSED) s1 = CLOSED; + if (LinkStatus[k] <= CLOSED) s1 = CLOSED; else s1 = OPEN; s2 = Control[i].Status; - k1 = K[k]; + k1 = LinkSetting[k]; k2 = k1; if (Link[k].Type > PIPE) k2 = Control[i].Setting; if (s1 != s2 || k1 != k2) { - S[k] = s2; - K[k] = k2; + LinkStatus[k] = s2; + LinkSetting[k] = k2; if (Statflag) writecontrolaction(k,i); // if (s1 != s2) initlinkflow(k, S[k], K[k]); setsum++; @@ -764,8 +764,8 @@ void tanktimestep(long *tstep) { if (Tank[i].A == 0.0) continue; /* Skip reservoirs */ n = Tank[i].Node; - h = H[n]; /* Current tank grade */ - q = D[n]; /* Flow into tank */ + h = NodeHead[n]; /* Current tank grade */ + q = NodeDemand[n]; /* Flow into tank */ if (ABS(q) <= QZERO) continue; if (q > 0.0 && h < Tank[i].Hmax) { @@ -802,8 +802,8 @@ void controltimestep(long *tstep) if ( (n = Control[i].Node) > 0) /* Node control: */ { if ((j = n-Njuncs) <= 0) continue; /* Node is a tank */ - h = H[n]; /* Current tank grade */ - q = D[n]; /* Flow into tank */ + h = NodeHead[n]; /* Current tank grade */ + q = NodeDemand[n]; /* Flow into tank */ if (ABS(q) <= QZERO) continue; if ( (h < Control[i].Grade && @@ -838,8 +838,8 @@ void controltimestep(long *tstep) /* Check if rule actually changes link status or setting */ k = Control[i].Link; if ( - (Link[k].Type > PIPE && K[k] != Control[i].Setting) || - (S[k] != Control[i].Status) + (Link[k].Type > PIPE && LinkSetting[k] != Control[i].Setting) || + (LinkStatus[k] != Control[i].Status) ) *tstep = t; } @@ -954,7 +954,7 @@ void addenergy(long hstep) { /* Skip closed pumps */ k = Pump[j].Link; - if (S[k] <= CLOSED) continue; + if (LinkStatus[k] <= CLOSED) continue; q = MAX(QZERO, ABS(Q[k])); /* Find pump-specific energy cost */ @@ -1001,7 +1001,7 @@ void getenergy(int k, double *kw, double *eff) /*** Updated 6/24/02 ***/ /* No energy if link is closed */ - if (S[k] <= CLOSED) + if (LinkStatus[k] <= CLOSED) { *kw = 0.0; *eff = 0.0; @@ -1011,7 +1011,7 @@ void getenergy(int k, double *kw, double *eff) /* Determine flow and head difference */ q = ABS(Q[k]); - dh = ABS(H[Link[k].N1] - H[Link[k].N2]); + dh = ABS(NodeHead[Link[k].N1] - NodeHead[Link[k].N2]); /* For pumps, find effic. at current flow */ if (Link[k].Type == PUMP) @@ -1053,18 +1053,18 @@ void tanklevels(long tstep) /* Update the tank's volume & water elevation */ n = Tank[i].Node; - dv = D[n]*tstep; + dv = NodeDemand[n]*tstep; Tank[i].V += dv; /*** Updated 6/24/02 ***/ /* Check if tank full/empty within next second */ - if (Tank[i].V + D[n] >= Tank[i].Vmax) { + if (Tank[i].V + NodeDemand[n] >= Tank[i].Vmax) { Tank[i].V = Tank[i].Vmax; } - else if (Tank[i].V - D[n] <= Tank[i].Vmin) { + else if (Tank[i].V - NodeDemand[n] <= Tank[i].Vmin) { Tank[i].V = Tank[i].Vmin; } - H[n] = tankgrade(i,Tank[i].V); + NodeHead[n] = tankgrade(i,Tank[i].V); } } /* End of tanklevels */ @@ -1188,7 +1188,7 @@ int netsolve(int *iter, double *relerr) /* Update current solution. */ /* (Row[i] = row of solution matrix corresponding to node i). */ - for (i=1; i<=Njuncs; i++) H[i] = F[Row[i]]; /* Update heads */ + for (i=1; i<=Njuncs; i++) NodeHead[i] = F[Row[i]]; /* Update heads */ newerr = newflows(); /* Update flows */ *relerr = newerr; @@ -1244,7 +1244,7 @@ int netsolve(int *iter, double *relerr) } /* Add any emitter flows to junction demands */ - for (i=1; i<=Njuncs; i++) D[i] += E[i]; + for (i=1; i<=Njuncs; i++) NodeDemand[i] += E[i]; return(errcode); } /* End of netsolve */ @@ -1274,15 +1274,15 @@ int badvalve(int n) Link[k].Type == PSV || Link[k].Type == FCV) { - if (S[k] == ACTIVE) + if (LinkStatus[k] == ACTIVE) { if (Statflag == FULL) { sprintf(Msg,FMT61,clocktime(Atime,Htime),Link[k].ID); writeline(Msg); } - if (Link[k].Type == FCV) S[k] = XFCV; - else S[k] = XPRESSURE; + if (Link[k].Type == FCV) LinkStatus[k] = XFCV; + else LinkStatus[k] = XPRESSURE; return(1); } } @@ -1313,25 +1313,25 @@ int valvestatus() for (i=1; i<=Nvalves; i++) /* Examine each valve */ { k = Valve[i].Link; /* Link index of valve */ - if (K[k] == MISSING) continue; /* Valve status fixed */ + if (LinkSetting[k] == MISSING) continue; /* Valve status fixed */ n1 = Link[k].N1; /* Start & end nodes */ n2 = Link[k].N2; - s = S[k]; /* Save current status */ + s = LinkStatus[k]; /* Save current status */ // if (s != CLOSED /* No change if flow is */ //(2.00.11 - LR) // && ABS(Q[k]) < Qtol) continue; /* negligible. */ //(2.00.11 - LR) switch (Link[k].Type) /* Evaluate new status: */ { - case PRV: hset = Node[n2].El + K[k]; - S[k] = prvstatus(k,s,hset,H[n1],H[n2]); + case PRV: hset = Node[n2].El + LinkSetting[k]; + LinkStatus[k] = prvstatus(k,s,hset,NodeHead[n1],NodeHead[n2]); break; - case PSV: hset = Node[n1].El + K[k]; - S[k] = psvstatus(k,s,hset,H[n1],H[n2]); + case PSV: hset = Node[n1].El + LinkSetting[k]; + LinkStatus[k] = psvstatus(k,s,hset,NodeHead[n1],NodeHead[n2]); break; //// FCV status checks moved back into the linkstatus() function //// //(2.00.12 - LR) -// case FCV: S[k] = fcvstatus(k,s,H[n1],H[n2]); //(2.00.12 - LR) +// case FCV: S[k] = fcvstatus(k,s,NodeHead[n1],NodeHead[n2]); //(2.00.12 - LR) // break; //(2.00.12 - LR) default: continue; @@ -1342,9 +1342,9 @@ int valvestatus() /* This strategy improves convergence. */ /* Check for status change */ - if (s != S[k]) + if (s != LinkStatus[k]) { - if (Statflag == FULL) writestatchange(k,s,S[k]); + if (Statflag == FULL) writestatchange(k,s,LinkStatus[k]); change = TRUE; } } @@ -1375,29 +1375,29 @@ int linkstatus() { n1 = Link[k].N1; n2 = Link[k].N2; - dh = H[n1] - H[n2]; + dh = NodeHead[n1] - NodeHead[n2]; /* Re-open temporarily closed links (status = XHEAD or TEMPCLOSED) */ - status = S[k]; - if (status == XHEAD || status == TEMPCLOSED) S[k] = OPEN; + status = LinkStatus[k]; + if (status == XHEAD || status == TEMPCLOSED) LinkStatus[k] = OPEN; /* Check for status changes in CVs and pumps */ - if (Link[k].Type == CV) S[k] = cvstatus(S[k],dh,Q[k]); - if (Link[k].Type == PUMP && S[k] >= OPEN && K[k] > 0.0) //(2.00.11 - LR) - S[k] = pumpstatus(k,-dh); + if (Link[k].Type == CV) LinkStatus[k] = cvstatus(LinkStatus[k],dh,Q[k]); + if (Link[k].Type == PUMP && LinkStatus[k] >= OPEN && LinkSetting[k] > 0.0) //(2.00.11 - LR) + LinkStatus[k] = pumpstatus(k,-dh); /* Check for status changes in non-fixed FCVs */ - if (Link[k].Type == FCV && K[k] != MISSING) //(2.00.12 - LR)// - S[k] = fcvstatus(k,status,H[n1],H[n2]); //(2.00.12 - LR)// + if (Link[k].Type == FCV && LinkSetting[k] != MISSING) //(2.00.12 - LR)// + LinkStatus[k] = fcvstatus(k,status,NodeHead[n1],NodeHead[n2]); //(2.00.12 - LR)// /* Check for flow into (out of) full (empty) tanks */ if (n1 > Njuncs || n2 > Njuncs) tankstatus(k,n1,n2); /* Note change in link status; do not revise link flow */ //(2.00.11 - LR) - if (status != S[k]) + if (status != LinkStatus[k]) { change = TRUE; - if (Statflag == FULL) writestatchange(k,status,S[k]); + if (Statflag == FULL) writestatchange(k,status,LinkStatus[k]); //if (S[k] <= CLOSED) Q[k] = QZERO; //(2.00.11 - LR) //else setlinkflow(k, dh); //(2.00.11 - LR) @@ -1449,7 +1449,7 @@ char pumpstatus(int k, double dh) /* Prevent reverse flow through pump */ p = PUMPINDEX(k); if (Pump[p].Ptype == CONST_HP) hmax = BIG; - else hmax = SQR(K[k])*Pump[p].Hmax; + else hmax = SQR(LinkSetting[k])*Pump[p].Hmax; if (dh > hmax + Htol) return(XHEAD); /*** Flow higher than pump curve no longer results in a status change ***/ //(2.00.11 - LR) @@ -1477,7 +1477,7 @@ char prvstatus(int k, char s, double hset, double h1, double h2) double htol = Htol; status = s; - if (K[k] == MISSING) return(status); /* Status fixed by user */ + if (LinkSetting[k] == MISSING) return(status); /* Status fixed by user */ hml = Link[k].Km*SQR(Q[k]); /* Head loss when open */ /*** Status rules below have changed. ***/ //(2.00.11 - LR) @@ -1527,7 +1527,7 @@ char psvstatus(int k, char s, double hset, double h1, double h2) double htol = Htol; status = s; - if (K[k] == MISSING) return(status); /* Status fixed by user */ + if (LinkSetting[k] == MISSING) return(status); /* Status fixed by user */ hml = Link[k].Km*SQR(Q[k]); /* Head loss when open */ /*** Status rules below have changed. ***/ //(2.00.11 - LR) @@ -1580,9 +1580,15 @@ char fcvstatus(int k, char s, double h1, double h2) { char status; /* New valve status */ status = s; - if (h1 - h2 < -Htol) status = XFCV; - else if ( Q[k] < -Qtol ) status = XFCV; //(2.00.11 - LR) - else if (s == XFCV && Q[k] >= K[k]) status = ACTIVE; + if (h1 - h2 < -Htol) { + status = XFCV; + } + else if ( Q[k] < -Qtol ) { + status = XFCV; //(2.00.11 - LR) + } + else if (s == XFCV && Q[k] >= LinkSetting[k]) { + status = ACTIVE; + } return(status); } @@ -1615,39 +1621,39 @@ void tankstatus(int k, int n1, int n2) n2 = n; q = -q; } - h = H[n1] - H[n2]; + h = NodeHead[n1] - NodeHead[n2]; /* Skip reservoirs & closed links */ - if (Tank[i].A == 0.0 || S[k] <= CLOSED) return; + if (Tank[i].A == 0.0 || LinkStatus[k] <= CLOSED) return; /* If tank full, then prevent flow into it */ - if (H[n1] >= Tank[i].Hmax - Htol) + if (NodeHead[n1] >= Tank[i].Hmax - Htol) { /* Case 1: Link is a pump discharging into tank */ if ( Link[k].Type == PUMP ) { - if (Link[k].N2 == n1) S[k] = TEMPCLOSED; + if (Link[k].N2 == n1) LinkStatus[k] = TEMPCLOSED; } /* Case 2: Downstream head > tank head */ /* (i.e., an open outflow check valve would close) */ - else if (cvstatus(OPEN, h, q) == CLOSED) S[k] = TEMPCLOSED; + else if (cvstatus(OPEN, h, q) == CLOSED) LinkStatus[k] = TEMPCLOSED; } /* If tank empty, then prevent flow out of it */ - if (H[n1] <= Tank[i].Hmin + Htol) + if (NodeHead[n1] <= Tank[i].Hmin + Htol) { /* Case 1: Link is a pump discharging from tank */ if ( Link[k].Type == PUMP) { - if (Link[k].N1 == n1) S[k] = TEMPCLOSED; + if (Link[k].N1 == n1) LinkStatus[k] = TEMPCLOSED; } /* Case 2: Tank head > downstream head */ /* (i.e., a closed outflow check valve would open) */ - else if (cvstatus(CLOSED, h, q) == OPEN) S[k] = TEMPCLOSED; + else if (cvstatus(CLOSED, h, q) == OPEN) LinkStatus[k] = TEMPCLOSED; } } /* End of tankstatus */ @@ -1681,10 +1687,10 @@ int pswitch() { /* Determine if control conditions are satisfied */ if (Control[i].Type == LOWLEVEL - && H[n] <= Control[i].Grade + Htol ) + && NodeHead[n] <= Control[i].Grade + Htol ) reset = 1; if (Control[i].Type == HILEVEL - && H[n] >= Control[i].Grade - Htol ) + && NodeHead[n] >= Control[i].Grade - Htol ) reset = 1; } @@ -1692,28 +1698,28 @@ int pswitch() if (reset == 1) { change = 0; - s = S[k]; + s = LinkStatus[k]; if (Link[k].Type == PIPE) { if (s != Control[i].Status) change = 1; } if (Link[k].Type == PUMP) { - if (K[k] != Control[i].Setting) change = 1; + if (LinkSetting[k] != Control[i].Setting) change = 1; } if (Link[k].Type >= PRV) { - if (K[k] != Control[i].Setting) change = 1; - else if (K[k] == MISSING && + if (LinkSetting[k] != Control[i].Setting) change = 1; + else if (LinkSetting[k] == MISSING && s != Control[i].Status) change = 1; } /* If a change occurs, update status & setting */ if (change) { - S[k] = Control[i].Status; - if (Link[k].Type > PIPE) K[k] = Control[i].Setting; - if (Statflag == FULL) writestatchange(k,s,S[k]); + LinkStatus[k] = Control[i].Status; + if (Link[k].Type > PIPE) LinkSetting[k] = Control[i].Setting; + if (Statflag == FULL) writestatchange(k,s,LinkStatus[k]); /* Re-set flow if status has changed */ // if (S[k] != s) initlinkflow(k, S[k], K[k]); @@ -1741,7 +1747,7 @@ double newflows() int k, n, n1, n2; /* Initialize net inflows (i.e., demands) at tanks */ - for (n=Njuncs+1; n <= Nnodes; n++) D[n] = 0.0; + for (n=Njuncs+1; n <= Nnodes; n++) NodeDemand[n] = 0.0; /* Initialize sum of flows & corrections */ qsum = 0.0; @@ -1761,7 +1767,7 @@ double newflows() n1 = Link[k].N1; n2 = Link[k].N2; - dh = H[n1] - H[n2]; + dh = NodeHead[n1] - NodeHead[n2]; dq = Y[k] - P[k]*dh; /* Adjust flow change by the relaxation factor */ //(2.00.11 - LR) @@ -1780,10 +1786,10 @@ double newflows() dqsum += ABS(dq); /* Update net flows to tanks */ - if ( S[k] > CLOSED ) //(2.00.12 - LR) + if ( LinkStatus[k] > CLOSED ) //(2.00.12 - LR) { - if (n1 > Njuncs) D[n1] -= Q[k]; - if (n2 > Njuncs) D[n2] += Q[k]; + if (n1 > Njuncs) NodeDemand[n1] -= Q[k]; + if (n2 > Njuncs) NodeDemand[n2] += Q[k]; } } @@ -1861,7 +1867,7 @@ void linkcoeffs() case PRV: case PSV: /* If valve status fixed then treat as pipe */ /* otherwise ignore the valve for now. */ - if (K[k] == MISSING) valvecoeff(k); //pipecoeff(k); //(2.00.11 - LR) + if (LinkSetting[k] == MISSING) valvecoeff(k); //pipecoeff(k); //(2.00.11 - LR) else continue; break; default: continue; @@ -1877,13 +1883,13 @@ void linkcoeffs() Aii[Row[n1]] += P[k]; /* Diagonal coeff. */ F[Row[n1]] += Y[k]; /* RHS coeff. */ } - else F[Row[n2]] += (P[k]*H[n1]); /* Node n1 is a tank */ + else F[Row[n2]] += (P[k]*NodeHead[n1]); /* Node n1 is a tank */ if (n2 <= Njuncs) /* Node n2 is junction */ { Aii[Row[n2]] += P[k]; /* Diagonal coeff. */ F[Row[n2]] -= Y[k]; /* RHS coeff. */ } - else F[Row[n1]] += (P[k]*H[n2]); /* Node n2 is a tank */ + else F[Row[n1]] += (P[k]*NodeHead[n2]); /* Node n2 is a tank */ } } /* End of linkcoeffs */ @@ -1904,7 +1910,7 @@ void nodecoeffs() /* flow imbalance & add imbalance to RHS array F. */ for (i=1; i<=Njuncs; i++) { - X[i] -= D[i]; + X[i] -= NodeDemand[i]; F[Row[i]] += X[i]; } } /* End of nodecoeffs */ @@ -1925,7 +1931,7 @@ void valvecoeffs() for (i=1; i<=Nvalves; i++) /* Examine each valve */ { k = Valve[i].Link; /* Link index of valve */ - if (K[k] == MISSING) continue; /* Valve status fixed */ + if (LinkSetting[k] == MISSING) continue; /* Valve status fixed */ n1 = Link[k].N1; /* Start & end nodes */ n2 = Link[k].N2; switch (Link[k].Type) /* Call valve-specific */ @@ -1992,7 +1998,7 @@ double emitflowchange(int i) p = 1/RQtol; else p = 1.0/p; - return(E[i]/Qexp - p*(H[i] - Node[i].El)); + return(E[i]/Qexp - p*(NodeHead[i] - Node[i].El)); } @@ -2019,7 +2025,7 @@ void pipecoeff(int k) dfdq; /* Derivative of fric. fact. */ /* For closed pipe use headloss formula: h = CBIG*q */ - if (S[k] <= CLOSED) + if (LinkStatus[k] <= CLOSED) { P[k] = 1.0/CBIG; Y[k] = Q[k]; @@ -2148,8 +2154,10 @@ void pumpcoeff(int k) r, /* Flow resistance coeff. */ n; /* Flow exponent coeff. */ + double setting = LinkSetting[k]; + /* Use high resistance pipe if pump closed or cannot deliver head */ - if (S[k] <= CLOSED || K[k] == 0.0) + if (LinkStatus[k] <= CLOSED || setting == 0.0) { //pipecoeff(k); //(2.00.11 - LR) P[k] = 1.0/CBIG; //(2.00.11 - LR) @@ -2166,7 +2174,7 @@ void pumpcoeff(int k) { /* Find intercept (h0) & slope (r) of pump curve */ /* line segment which contains speed-adjusted flow. */ - curvecoeff(Pump[p].Hcurve,q/K[k],&h0,&r); + curvecoeff(Pump[p].Hcurve, q/setting, &h0, &r); /* Determine head loss coefficients. */ Pump[p].H0 = -h0; @@ -2175,9 +2183,9 @@ void pumpcoeff(int k) } /* Adjust head loss coefficients for pump speed. */ - h0 = SQR(K[k])*Pump[p].H0; + h0 = SQR(setting)*Pump[p].H0; n = Pump[p].N; - r = Pump[p].R*pow(K[k],2.0-n); + r = Pump[p].R*pow(setting,2.0-n); if (n != 1.0) r = n*r*pow(q,n-1.0); /* Compute inverse headloss gradient (P) and flow correction factor (Y) */ @@ -2241,7 +2249,7 @@ void gpvcoeff(int k) /*** Updated 9/7/00 ***/ /* Treat as a pipe if valve closed */ - if (S[k] == CLOSED) valvecoeff(k); //pipecoeff(k); //(2.00.11 - LR) + if (LinkStatus[k] == CLOSED) valvecoeff(k); //pipecoeff(k); //(2.00.11 - LR) /* Otherwise utilize headloss curve */ /* whose index is stored in K */ @@ -2253,7 +2261,7 @@ void gpvcoeff(int k) /*** Updated 10/25/00 ***/ /*** Updated 12/29/00 ***/ - curvecoeff((int)ROUND(K[k]),q,&h0,&r); + curvecoeff((int)ROUND(LinkSetting[k]),q,&h0,&r); /* Compute inverse headloss gradient (P) */ /* and flow correction factor (Y). */ @@ -2273,19 +2281,19 @@ void pbvcoeff(int k) */ { /* If valve fixed OPEN or CLOSED then treat as a pipe */ - if (K[k] == MISSING || K[k] == 0.0) valvecoeff(k); //pipecoeff(k); //(2.00.11 - LR) + if (LinkSetting[k] == MISSING || LinkSetting[k] == 0.0) valvecoeff(k); //pipecoeff(k); //(2.00.11 - LR) /* If valve is active */ else { /* Treat as a pipe if minor loss > valve setting */ - if (Link[k].Km*SQR(Q[k]) > K[k]) valvecoeff(k); //pipecoeff(k); //(2.00.11 - LR) + if (Link[k].Km*SQR(Q[k]) > LinkSetting[k]) valvecoeff(k); //pipecoeff(k); //(2.00.11 - LR) /* Otherwise force headloss across valve to be equal to setting */ else { P[k] = CBIG; - Y[k] = K[k]*CBIG; + Y[k] = LinkSetting[k]*CBIG; } } } /* End of pbvcoeff */ @@ -2306,8 +2314,8 @@ void tcvcoeff(int k) km = Link[k].Km; /* If valve not fixed OPEN or CLOSED, compute its loss coeff. */ - if (K[k] != MISSING) - Link[k].Km = 0.02517*K[k]/(SQR(Link[k].Diam)*SQR(Link[k].Diam)); + if (LinkSetting[k] != MISSING) + Link[k].Km = 0.02517*LinkSetting[k]/(SQR(Link[k].Diam)*SQR(Link[k].Diam)); /* Then apply usual pipe formulas */ valvecoeff(k); //pipecoeff(k); //(2.00.11 - LR) @@ -2333,9 +2341,9 @@ void prvcoeff(int k, int n1, int n2) double hset; /* Valve head setting */ i = Row[n1]; /* Matrix rows of nodes */ j = Row[n2]; - hset = Node[n2].El + K[k]; /* Valve setting */ + hset = Node[n2].El + LinkSetting[k]; /* Valve setting */ - if (S[k] == ACTIVE) + if (LinkStatus[k] == ACTIVE) { /* Set coeffs. to force head at downstream @@ -2379,9 +2387,9 @@ void psvcoeff(int k, int n1, int n2) double hset; /* Valve head setting */ i = Row[n1]; /* Matrix rows of nodes */ j = Row[n2]; - hset = Node[n1].El + K[k]; /* Valve setting */ + hset = Node[n1].El + LinkSetting[k]; /* Valve setting */ - if (S[k] == ACTIVE) + if (LinkStatus[k] == ACTIVE) { /* Set coeffs. to force head at upstream @@ -2423,7 +2431,7 @@ void fcvcoeff(int k, int n1, int n2) { int i,j; /* Rows in solution matrix */ double q; /* Valve flow setting */ - q = K[k]; + q = LinkSetting[k]; i = Row[n1]; j = Row[n2]; @@ -2432,7 +2440,7 @@ void fcvcoeff(int k, int n1, int n2) flow setting as external demand at upstream node and external supply at downstream node. */ - if (S[k] == ACTIVE) + if (LinkStatus[k] == ACTIVE) { X[n1] -= q; F[i] -= q; @@ -2474,7 +2482,7 @@ void valvecoeff(int k) double p; // Valve is closed. Use a very small matrix coeff. - if (S[k] <= CLOSED) + if (LinkStatus[k] <= CLOSED) { P[k] = 1.0/CBIG; Y[k] = Q[k]; diff --git a/src/input3.c b/src/input3.c index a19bfdc..d0922f6 100755 --- a/src/input3.c +++ b/src/input3.c @@ -99,9 +99,9 @@ int juncdata() demand->Pat = p; demand->next = Node[Njuncs].D; Node[Njuncs].D = demand; - D[Njuncs] = y; + NodeDemand[Njuncs] = y; } - else D[Njuncs] = MISSING; + else NodeDemand[Njuncs] = MISSING; /*** end of update ***/ return(0); } /* end of juncdata */ @@ -682,11 +682,11 @@ int demanddata() /*** Updated 6/24/02 ***/ demand = Node[j].D; - if (demand && D[j] != MISSING) + if (demand && NodeDemand[j] != MISSING) { demand->Base = y; demand->Pat = p; - D[j] = MISSING; + NodeDemand[j] = MISSING; } /*** End of update ***/ diff --git a/src/output.c b/src/output.c index 4237736..0e18bb8 100755 --- a/src/output.c +++ b/src/output.c @@ -152,30 +152,30 @@ int savehyd(long *htime) fwrite(&t,sizeof(INT4),1,HydFile); /* Save current nodal demands (D) */ - for (i=1; i<=Nnodes; i++) x[i] = (REAL4)D[i]; + for (i=1; i<=Nnodes; i++) x[i] = (REAL4)NodeDemand[i]; fwrite(x+1,sizeof(REAL4),Nnodes,HydFile); /* Copy heads (H) to buffer of floats (x) and save buffer */ - for (i=1; i<=Nnodes; i++) x[i] = (REAL4)H[i]; + for (i=1; i<=Nnodes; i++) x[i] = (REAL4)NodeHead[i]; fwrite(x+1,sizeof(REAL4),Nnodes,HydFile); /* Force flow in closed links to be zero then save flows */ for (i=1; i<=Nlinks; i++) { - if (S[i] <= CLOSED) x[i] = 0.0f; + if (LinkStatus[i] <= CLOSED) x[i] = 0.0f; else x[i] = (REAL4)Q[i]; } fwrite(x+1,sizeof(REAL4),Nlinks,HydFile); /* Copy link status to buffer of floats (x) & write buffer */ - for (i=1; i<=Nlinks; i++) x[i] = (REAL4)S[i]; + for (i=1; i<=Nlinks; i++) x[i] = (REAL4)LinkStatus[i]; fwrite(x+1,sizeof(REAL4),Nlinks,HydFile); /* Save link settings & check for successful write-to-disk */ /* (We assume that if any of the previous fwrites failed, */ /* then this one will also fail.) */ - for (i=1; i<=Nlinks; i++) x[i] = (REAL4)K[i]; + for (i=1; i<=Nlinks; i++) x[i] = (REAL4)LinkSetting[i]; if (fwrite(x+1,sizeof(REAL4),Nlinks,HydFile) < (unsigned)Nlinks) errcode = 308; free(x); @@ -286,19 +286,19 @@ int readhyd(long *hydtime) *hydtime = t; if (fread(x+1,sizeof(REAL4),Nnodes,HydFile) < (unsigned)Nnodes) result = 0; - else for (i=1; i<=Nnodes; i++) D[i] = x[i]; + else for (i=1; i<=Nnodes; i++) NodeDemand[i] = x[i]; if (fread(x+1,sizeof(REAL4),Nnodes,HydFile) < (unsigned)Nnodes) result = 0; - else for (i=1; i<=Nnodes; i++) H[i] = x[i]; + else for (i=1; i<=Nnodes; i++) NodeHead[i] = x[i]; if (fread(x+1,sizeof(REAL4),Nlinks,HydFile) < (unsigned)Nlinks) result = 0; else for (i=1; i<=Nlinks; i++) Q[i] = x[i]; if (fread(x+1,sizeof(REAL4),Nlinks,HydFile) < (unsigned)Nlinks) result = 0; - else for (i=1; i<=Nlinks; i++) S[i] = (char) x[i]; + else for (i=1; i<=Nlinks; i++) LinkStatus[i] = (char) x[i]; if (fread(x+1,sizeof(REAL4),Nlinks,HydFile) < (unsigned)Nlinks) result = 0; - else for (i=1; i<=Nlinks; i++) K[i] = x[i]; + else for (i=1; i<=Nlinks; i++) LinkSetting[i] = x[i]; free(x); return result; @@ -361,16 +361,16 @@ int nodeoutput(int j, REAL4 *x, double ucf) switch(j) { case DEMAND: for (i=1; i<=Nnodes; i++) - x[i] = (REAL4)(D[i]*ucf); + x[i] = (REAL4)(NodeDemand[i]*ucf); break; case HEAD: for (i=1; i<=Nnodes; i++) - x[i] = (REAL4)(H[i]*ucf); + x[i] = (REAL4)(NodeHead[i]*ucf); break; case PRESSURE: for (i=1; i<=Nnodes; i++) - x[i] = (REAL4)((H[i] - Node[i].El)*ucf); + x[i] = (REAL4)((NodeHead[i] - Node[i].El)*ucf); break; case QUALITY: for (i=1; i<=Nnodes; i++) - x[i] = (REAL4)(C[i]*ucf); + x[i] = (REAL4)(NodeQual[i]*ucf); } /* Write x[1] to x[Nnodes] to output file */ @@ -413,10 +413,10 @@ int linkoutput(int j, REAL4 *x, double ucf) break; case HEADLOSS: for (i=1; i<=Nlinks; i++) { - if (S[i] <= CLOSED) x[i] = 0.0f; + if (LinkStatus[i] <= CLOSED) x[i] = 0.0f; else { - h = H[Link[i].N1] - H[Link[i].N2]; + h = NodeHead[Link[i].N1] - NodeHead[Link[i].N2]; if (Link[i].Type != PUMP) h = ABS(h); if (Link[i].Type <= PIPE) x[i] = (REAL4)(1000.0*h/Link[i].Len); @@ -428,25 +428,26 @@ int linkoutput(int j, REAL4 *x, double ucf) x[i] = (REAL4)(avgqual(i)*ucf); break; case STATUS: for (i=1; i<=Nlinks; i++) - x[i] = (REAL4)S[i]; + x[i] = (REAL4)LinkStatus[i]; break; case SETTING: for (i=1; i<=Nlinks; i++) { - if (K[i] != MISSING) + double setting = LinkSetting[i]; + if (setting != MISSING) switch (Link[i].Type) { case CV: - case PIPE: x[i] = (REAL4)K[i]; + case PIPE: x[i] = (REAL4)setting; break; - case PUMP: x[i] = (REAL4)K[i]; + case PUMP: x[i] = (REAL4)setting; break; case PRV: case PSV: - case PBV: x[i] = (REAL4)(K[i]*Ucf[PRESSURE]); + case PBV: x[i] = (REAL4)(setting*Ucf[PRESSURE]); break; - case FCV: x[i] = (REAL4)(K[i]*Ucf[FLOW]); + case FCV: x[i] = (REAL4)(setting*Ucf[FLOW]); break; - case TCV: x[i] = (REAL4)K[i]; + case TCV: x[i] = (REAL4)setting; break; default: x[i] = 0.0f; } @@ -455,7 +456,7 @@ int linkoutput(int j, REAL4 *x, double ucf) break; case REACTRATE: /* Overall reaction rate in mass/L/day */ if (Qualflag == NONE) memset(x,0,(Nlinks+1 )*sizeof(REAL4)); - else for (i=1; i<=Nlinks; i++) x[i] = (REAL4)(R[i]*ucf); + else for (i=1; i<=Nlinks; i++) x[i] = (REAL4)(PipeRateCoeff[i]*ucf); break; case FRICTION: /* f = 2ghd/(Lu^2) where f = friction factor */ /* u = velocity, g = grav. accel., h = head */ @@ -464,7 +465,7 @@ int linkoutput(int j, REAL4 *x, double ucf) { if (Link[i].Type <= PIPE && ABS(Q[i]) > TINY) { - h = ABS(H[Link[i].N1] - H[Link[i].N2]); + h = ABS(NodeHead[Link[i].N1] - NodeHead[Link[i].N2]); f = 39.725*h*pow(Link[i].Diam,5)/Link[i].Len/SQR(Q[i]); x[i] = (REAL4)f; } @@ -642,11 +643,11 @@ int savetimestat(REAL4 *x, char objtype) /* Update internal output variables where applicable */ if (objtype == NODEHDR) switch (j) { - case DEMAND: for (i=1; i<=n; i++) D[i] = x[i]/Ucf[DEMAND]; + case DEMAND: for (i=1; i<=n; i++) NodeDemand[i] = x[i]/Ucf[DEMAND]; break; - case HEAD: for (i=1; i<=n; i++) H[i] = x[i]/Ucf[HEAD]; + case HEAD: for (i=1; i<=n; i++) NodeHead[i] = x[i]/Ucf[HEAD]; break; - case QUALITY: for (i=1; i<=n; i++) C[i] = x[i]/Ucf[QUALITY]; + case QUALITY: for (i=1; i<=n; i++) NodeQual[i] = x[i]/Ucf[QUALITY]; break; } else if (j == FLOW) for (i=1; i<=n; i++) Q[i] = x[i]/Ucf[FLOW]; diff --git a/src/quality.c b/src/quality.c index d17af7c..be9193d 100755 --- a/src/quality.c +++ b/src/quality.c @@ -106,10 +106,10 @@ int openqual() if (SegPool == NULL) errcode = 101; //(2.00.11 - LR) /* Allocate scratch array & reaction rate array*/ - XC = (double *) calloc(MAX((Nnodes+1),(Nlinks+1)),sizeof(double)); - R = (double *) calloc((Nlinks+1), sizeof(double)); - ERRCODE(MEMCHECK(XC)); - ERRCODE(MEMCHECK(R)); + TempQual = (double *) calloc(MAX((Nnodes+1),(Nlinks+1)),sizeof(double)); + PipeRateCoeff = (double *) calloc((Nlinks+1), sizeof(double)); + ERRCODE(MEMCHECK(TempQual)); + ERRCODE(MEMCHECK(PipeRateCoeff)); /* Allocate memory for WQ solver */ n = Nlinks+Ntanks+1; @@ -148,7 +148,7 @@ void initqual() int i; /* Initialize quality, tank volumes, & source mass flows */ - for (i=1; i<=Nnodes; i++) C[i] = Node[i].C0; + for (i=1; i<=Nnodes; i++) NodeQual[i] = Node[i].C0; for (i=1; i<=Ntanks; i++) Tank[i].C = Node[Tank[i].Node].C0; for (i=1; i<=Ntanks; i++) Tank[i].V = Tank[i].V0; for (i=1; i<=Nnodes; i++) { @@ -165,7 +165,7 @@ void initqual() if (Qualflag != NONE) { /* Initialize WQ at trace node (if applicable) */ - if (Qualflag == TRACE) C[TraceNode] = 100.0; + if (Qualflag == TRACE) NodeQual[TraceNode] = 100.0; /* Compute Schmidt number */ if (Diffus > 0.0) @@ -203,6 +203,8 @@ void initqual() Qtime = 0; Rtime = Rstart; Nperiods = 0; + + initsegs(); } @@ -240,7 +242,7 @@ int runqual(long *t) for (int i=1; i<= Nlinks; ++i) { - if (S[i] <= CLOSED) { + if (LinkStatus[i] <= CLOSED) { QLinkFlow[i-1] = Q[i]; } } @@ -255,7 +257,7 @@ int runqual(long *t) for (int i=1; i<= Nlinks; ++i) { - if (S[i] <= CLOSED) { + if (LinkStatus[i] <= CLOSED) { QLinkFlow[i-1] = Q[i]; } } @@ -304,13 +306,13 @@ int nextqual(long *tstep) if (Tank[i].A != 0) { // skip reservoirs again int n = Tank[i].Node; Tank[i].V = QTankVolumes[i-1]; - H[n] = tankgrade(i,Tank[i].V); + NodeHead[n] = tankgrade(i,Tank[i].V); } } // restore the previous step's pipe link flows for (int i=1; i<=Nlinks; i++) { - if (S[i] <= CLOSED) { + if (LinkStatus[i] <= CLOSED) { Q[i] = 0.0; } } @@ -334,12 +336,12 @@ int nextqual(long *tstep) if (Tank[i].A != 0) { // skip reservoirs again int n = Tank[i].Node; Tank[i].V = tankVolumes[i-1]; - H[n] = tankgrade(i,Tank[i].V); + NodeHead[n] = tankgrade(i,Tank[i].V); } } for (int i=1; i<=Nlinks; ++i) { - if (S[i] <= CLOSED) { + if (LinkStatus[i] <= CLOSED) { Q[i] = QLinkFlow[i-1]; } } @@ -412,8 +414,8 @@ int closequal() free(FlowDir); free(VolIn); free(MassIn); - free(R); - free(XC); + free(PipeRateCoeff); + free(TempQual); free(QTankVolumes); free(QLinkFlow); return(errcode); @@ -463,12 +465,20 @@ int gethyd(long *hydtime, long *hydstep) { /* Compute reaction rate coeffs. */ - if (Reactflag && Qualflag != AGE) ratecoeffs(); - + if (Reactflag && Qualflag != AGE) { + ratecoeffs(); + } + /* Initialize pipe segments (at time 0) or */ /* else re-orient segments if flow reverses.*/ - if (Qtime == 0) initsegs(); - else reorientsegs(); + //if (Qtime == 0) + // initsegs(); + //else + // if hydraulics are open, or if we're in sequential mode (where qtime can increase) + if (OpenHflag || Qtime != 0) { + reorientsegs(); + } + } return(errcode); } @@ -550,8 +560,10 @@ void initsegs() { /* Establish flow direction */ - FlowDir[k] = '+'; - if (Q[k] < 0.) FlowDir[k] = '-'; + FlowDir[k] = '+'; + if (Q[k] < 0.) { + FlowDir[k] = '-'; + } /* Set segs to zero */ LastSeg[k] = NULL; @@ -559,7 +571,7 @@ void initsegs() /* Find quality of downstream node */ j = DOWN_NODE(k); - if (j <= Njuncs) c = C[j]; + if (j <= Njuncs) c = NodeQual[j]; else c = Tank[j-Njuncs].C; /* Fill link with single segment with this quality */ @@ -617,9 +629,13 @@ void reorientsegs() { /* Find new flow direction */ - newdir = '+'; - if (Q[k] == 0.0) newdir = FlowDir[k]; - else if (Q[k] < 0.0) newdir = '-'; + newdir = '+'; + if (Q[k] == 0.0) { + newdir = FlowDir[k]; + } + else if (Q[k] < 0.0) { + newdir = '-'; + } /* If direction changes, then reverse order of segments */ /* (first to last) and save new direction */ @@ -683,8 +699,8 @@ void updatesegs(long dt) } /* Normalize volume-weighted reaction rate */ - if (vsum > 0.0) R[k] = rsum/vsum/dt*SECperDAY; - else R[k] = 0.0; + if (vsum > 0.0) PipeRateCoeff[k] = rsum/vsum/dt*SECperDAY; + else PipeRateCoeff[k] = 0.0; } } @@ -765,7 +781,7 @@ void accumulate(long dt) /* Re-set memory used to accumulate mass & volume */ memset(VolIn,0,(Nnodes+1)*sizeof(double)); memset(MassIn,0,(Nnodes+1)*sizeof(double)); - memset(XC,0,(Nnodes+1)*sizeof(double)); + memset(TempQual,0,(Nnodes+1)*sizeof(double)); /* Compute average conc. of segments adjacent to each node */ /* (For use if there is no transport through the node) */ @@ -784,9 +800,13 @@ void accumulate(long dt) VolIn[j]++; } } - for (k=1; k<=Nnodes; k++) - if (VolIn[k] > 0.0) XC[k] = MassIn[k]/VolIn[k]; - + + for (k=1; k<=Nnodes; k++) { + if (VolIn[k] > 0.0) { + TempQual[k] = MassIn[k]/VolIn[k]; + } + } + /* Move mass from first segment of each pipe into downstream node */ memset(VolIn,0,(Nnodes+1)*sizeof(double)); memset(MassIn,0,(Nnodes+1)*sizeof(double)); @@ -804,7 +824,7 @@ void accumulate(long dt) { VolIn[j] += v; seg = FirstSeg[k]; - cseg = C[i]; + cseg = NodeQuali]; if (seg != NULL) cseg = seg->c; MassIn[j] += v*cseg; removesegs(k); @@ -866,27 +886,33 @@ void updatenodes(long dt) ** Purpose: updates concentration at all nodes to mixture of accumulated ** inflow from connecting pipes. ** -** Note: Does not account for source flow effects. XC[i] contains +** Note: Does not account for source flow effects. TempQual[i] contains ** average concen. of segments adjacent to node i, used in case ** there was no inflow into i. **--------------------------------------------------------------------------- */ { - int i; - - /* Update junction quality */ - for (i=1; i<=Njuncs; i++) - { - if (D[i] < 0.0) VolIn[i] -= D[i]*dt; - if (VolIn[i] > 0.0) C[i] = MassIn[i]/VolIn[i]; - else C[i] = XC[i]; - } - - /* Update tank quality */ - updatetanks(dt); - - /* For flow tracing, set source node concen. to 100. */ - if (Qualflag == TRACE) C[TraceNode] = 100.0; + int i; + + /* Update junction quality */ + for (i=1; i<=Njuncs; i++) + { + if (NodeDemand[i] < 0.0) { + VolIn[i] -= NodeDemand[i]*dt; + } + if (VolIn[i] > 0.0) { + NodeQual[i] = MassIn[i]/VolIn[i]; + } + else { + NodeQual[i] = TempQual[i]; + } + } + + /* Update tank quality */ + updatetanks(dt); + + /* For flow tracing, set source node concen. to 100. */ + if (Qualflag == TRACE) NodeQual[TraceNode] = 100.0; } @@ -908,14 +934,14 @@ void sourceinput(long dt) /* Establish a flow cutoff which indicates no outflow from a node */ qcutoff = 10.0*TINY; - /* Zero-out the work array XC */ - memset(XC,0,(Nnodes+1)*sizeof(double)); + /* Zero-out the work array TempQual */ + memset(TempQual,0,(Nnodes+1)*sizeof(double)); if (Qualflag != CHEM) return; /* Consider each node */ for (n=1; n<=Nnodes; n++) { - + double thisDemand = NodeDemand[n]; /* Skip node if no WQ source */ source = Node[n].S; if (source == NULL) continue; @@ -923,7 +949,7 @@ void sourceinput(long dt) /* Find total flow volume leaving node */ if (n <= Njuncs) volout = VolIn[n]; /* Junctions */ - else volout = VolIn[n] - D[n]*dt; /* Tanks */ + else volout = VolIn[n] - (thisDemand * dt); /* Tanks */ qout = volout / (double) dt; /* Evaluate source input only if node outflow > cutoff flow */ @@ -939,13 +965,13 @@ void sourceinput(long dt) case CONCEN: /* Only add source mass if demand is negative */ - if (D[n] < 0.0) + if (thisDemand < 0.0) { - massadded = -s*D[n]*dt; + massadded = -s*thisDemand*dt; /* If node is a tank then set concen. to 0. */ /* (It will be re-set to true value in updatesourcenodes()) */ - if (n > Njuncs) C[n] = 0.0; + if (n > Njuncs) NodeQual[n] = 0.0; } else massadded = 0.0; break; @@ -959,9 +985,13 @@ void sourceinput(long dt) /* Mass added is difference between source */ /* & node concen. times outflow volume */ case SETPOINT: - if (s > C[n]) massadded = (s-C[n])*volout; - else massadded = 0.0; - break; + if (s > NodeQual[n]) { + massadded = (s-NodeQual[n])*volout; + } + else { + massadded = 0.0; + } + break; /* Flow-Paced Booster Source: */ /* Mass added = source concen. times outflow volume */ @@ -971,7 +1001,7 @@ void sourceinput(long dt) } /* Source concen. contribution = (mass added / outflow volume) */ - XC[n] = massadded/volout; + TempQual[n] = massadded/volout; /* Update total mass added for time period & simulation */ source->Smass += massadded; @@ -987,8 +1017,8 @@ void sourceinput(long dt) if (Tank[j].A == 0.0) { n = Njuncs + j; - volout = VolIn[n] - D[n]*dt; - if (volout > 0.0) Wsource += volout*C[n]; + volout = VolIn[n] - NodeDemand[n]*dt; + if (volout > 0.0) Wsource += volout*NodeQual[n]; } } } @@ -1022,7 +1052,7 @@ void release(long dt) v = q*dt; /* Include source contribution in quality released from node. */ - c = C[n] + XC[n]; + c = NodeQual[n] + TempQual[n]; /* If link has a last seg, check if its quality */ /* differs from that of the flow released from node.*/ @@ -1051,7 +1081,7 @@ void updatesourcenodes(long dt) ** Input: dt = current WQ time step ** Output: none ** Purpose: updates quality at source nodes. -** (XC[n] = concen. added by source at node n) +** (TempQual[n] = concen. added by source at node n) **--------------------------------------------------- */ { @@ -1067,13 +1097,13 @@ void updatesourcenodes(long dt) if (source == NULL) continue; /* Add source to current node concen. */ - C[n] += XC[n]; + NodeQual[n] += TempQual[n]; /* For tanks, node concen. = internal concen. */ if (n > Njuncs) { i = n - Njuncs; - if (Tank[i].A > 0.0) C[n] = Tank[i].C; + if (Tank[i].A > 0.0) NodeQual[n] = Tank[i].C; } /* Normalize mass added at source to time step */ @@ -1100,7 +1130,7 @@ void updatetanks(long dt) /* Use initial quality for reservoirs */ if (Tank[i].A == 0.0) { - C[n] = Node[n].C0; + NodeQual[n] = Node[n].C0; } /* Update tank WQ based on mixing model */ else { @@ -1143,7 +1173,7 @@ void updatetanks(long dt) // /* Update tank volume & nodal quality */ // Tank[i].V += D[n]*dt; -// C[n] = Tank[i].C; +// NodeQual[n] = Tank[i].C; //} @@ -1168,7 +1198,7 @@ void tankmix1(int i, long dt) /* Determine tank & volumes */ vold = Tank[i].V; n = Tank[i].Node; - Tank[i].V += D[n]*dt; + Tank[i].V += NodeDemand[n]*dt; vin = VolIn[n]; /* Compute inflow concen. */ @@ -1181,7 +1211,7 @@ void tankmix1(int i, long dt) c = MIN(c, cmax); c = MAX(c, 0.0); Tank[i].C = c; - C[n] = Tank[i].C; + NodeQual[n] = Tank[i].C; } /*** Updated 10/25/00 ***/ @@ -1218,7 +1248,7 @@ void tankmix2(int i, long dt) /* Find inflows & outflows */ n = Tank[i].Node; - vnet = D[n]*dt; + vnet = NodeDemand[n]*dt; vin = VolIn[n]; if (vin > 0.0) cin = MassIn[n]/vin; else cin = 0.0; @@ -1274,7 +1304,7 @@ void tankmix2(int i, long dt) /* represent quality of tank since this is where */ /* outflow begins to flow from */ Tank[i].C = seg1->c; - C[n] = Tank[i].C; + NodeQual[n] = Tank[i].C; } @@ -1309,7 +1339,7 @@ void tankmix3(int i, long dt) /* Find inflows & outflows */ n = Tank[i].Node; - vnet = D[n]*dt; + vnet = NodeDemand[n]*dt; vin = VolIn[n]; vout = vin - vnet; if (vin > 0.0) cin = MassIn[n]/VolIn[n]; @@ -1350,7 +1380,7 @@ void tankmix3(int i, long dt) /* to represent overall quality of tank */ if (vsum > 0.0) Tank[i].C = csum/vsum; else Tank[i].C = FirstSeg[k]->c; - C[n] = Tank[i].C; + NodeQual[n] = Tank[i].C; /* Add new last segment for new flow entering tank */ if (vin > 0.0) @@ -1400,7 +1430,7 @@ void tankmix4(int i, long dt) /* Find inflows & outflows */ n = Tank[i].Node; - vnet = D[n]*dt; + vnet = NodeDemand[n]*dt; vin = VolIn[n]; if (vin > 0.0) cin = MassIn[n]/VolIn[n]; else cin = 0.0; @@ -1468,7 +1498,7 @@ void tankmix4(int i, long dt) /* Reported tank quality is mixture of flow released and any inflow */ Tank[i].C = (csum + MassIn[n])/(vsum + vin); } - C[n] = Tank[i].C; + NodeQual[n] = Tank[i].C; } @@ -1523,7 +1553,7 @@ double avgqual(int k) seg = seg->prev; } if (vsum > 0.0) return(msum/vsum); - else return( (C[Link[k].N1] + C[Link[k].N2])/2. ); + else return( (NodeQual[Link[k].N1] + NodeQual[Link[k].N2])/2. ); } @@ -1544,7 +1574,7 @@ void ratecoeffs() kw = Link[k].Kw; if (kw != 0.0) kw = piperate(k); Link[k].Rc = kw; - R[k] = 0.0; + PipeRateCoeff[k] = 0.0; } } /* End of ratecoeffs */ diff --git a/src/report.c b/src/report.c index b021461..deffe03 100755 --- a/src/report.c +++ b/src/report.c @@ -296,15 +296,15 @@ void writehydstat(int iter, double relerr) for (i=1; i<=Ntanks; i++) { n = Tank[i].Node; - if (ABS(D[n]) < 0.001) newstat = CLOSED; - else if (D[n] > 0.0) newstat = FILLING; - else if (D[n] < 0.0) newstat = EMPTYING; + if (ABS(NodeDemand[n]) < 0.001) newstat = CLOSED; + else if (NodeDemand[n] > 0.0) newstat = FILLING; + else if (NodeDemand[n] < 0.0) newstat = EMPTYING; else newstat = OldStat[Nlinks+i]; if (newstat != OldStat[Nlinks+i]) { if (Tank[i].A > 0.0) sprintf(s1,FMT50,atime,Node[n].ID,StatTxt[newstat], - (H[n]-Node[n].El)*Ucf[HEAD],Field[HEAD].Units); + (NodeHead[n]-Node[n].El)*Ucf[HEAD],Field[HEAD].Units); else sprintf(s1,FMT51,atime,Node[n].ID,StatTxt[newstat]); writeline(s1); OldStat[Nlinks+i] = newstat; @@ -314,15 +314,15 @@ void writehydstat(int iter, double relerr) /* Display status changes for links */ for (i=1; i<=Nlinks; i++) { - if (S[i] != OldStat[i]) + if (LinkStatus[i] != OldStat[i]) { if (Htime == 0) sprintf(s1,FMT52,atime,LinkTxt[Link[i].Type],Link[i].ID, - StatTxt[S[i]]); + StatTxt[LinkStatus[i]]); else sprintf(s1,FMT53,atime,LinkTxt[Link[i].Type],Link[i].ID, - StatTxt[OldStat[i]],StatTxt[S[i]]); + StatTxt[OldStat[i]],StatTxt[LinkStatus[i]]); writeline(s1); - OldStat[i] = S[i]; + OldStat[i] = LinkStatus[i]; } } writeline(" "); @@ -763,7 +763,7 @@ void writestatchange(int k, char s1, char s2) { /*** Updated 10/25/00 ***/ - setting = K[k]; //Link[k].Kc; + setting = LinkSetting[k]; //Link[k].Kc; switch (Link[k].Type) { @@ -871,7 +871,7 @@ int writehydwarn(int iter, double relerr) /* Check for negative pressures */ for (i=1; i<=Njuncs; i++) { - if (H[i] < Node[i].El && D[i] > 0.0) + if (NodeHead[i] < Node[i].El && NodeDemand[i] > 0.0) { sprintf(Msg,WARN06,clocktime(Atime,Htime)); if (Messageflag) writeline(Msg); @@ -884,10 +884,10 @@ int writehydwarn(int iter, double relerr) for (i=1; i<=Nvalves; i++) { j = Valve[i].Link; - if (S[j] >= XFCV) + if (LinkStatus[j] >= XFCV) { sprintf(Msg,WARN05,LinkTxt[Link[j].Type],Link[j].ID, - StatTxt[S[j]],clocktime(Atime,Htime)); + StatTxt[LinkStatus[j]],clocktime(Atime,Htime)); if (Messageflag) writeline(Msg); flag = 5; } @@ -897,10 +897,10 @@ int writehydwarn(int iter, double relerr) for (i=1; i<=Npumps; i++) { j = Pump[i].Link; - s = S[j]; //(2.00.11 - LR) - if (S[j] >= OPEN) //(2.00.11 - LR) + s = LinkStatus[j]; //(2.00.11 - LR) + if (LinkStatus[j] >= OPEN) //(2.00.11 - LR) { //(2.00.11 - LR) - if (Q[j] > K[j]*Pump[i].Qmax) s = XFLOW; //(2.00.11 - LR) + if (Q[j] > LinkSetting[j]*Pump[i].Qmax) s = XFLOW; //(2.00.11 - LR) if (Q[j] < 0.0) s = XHEAD; //(2.00.11 - LR) } //(2.00.11 - LR) if (s == XHEAD || s == XFLOW) //(2.00.11 - LR) @@ -984,7 +984,7 @@ int disconnected() mcount = Ntanks; for (i=1; i<=Njuncs; i++) { - if (D[i] < 0.0) + if (NodeDemand[i] < 0.0) { mcount++; nodelist[mcount] = i; @@ -999,7 +999,7 @@ int disconnected() count = 0; for (i=1; i<=Njuncs; i++) { - if (!marked[i] && D[i] != 0.0) /* Skip if no demand */ + if (!marked[i] && NodeDemand[i] != 0.0) /* Skip if no demand */ { count++; if (count <= MAXCOUNT && Messageflag) @@ -1068,7 +1068,7 @@ void marknodes(int m, int *nodelist, char *marked) } /* Mark connection node if link not closed */ - if (S[k] > CLOSED) + if (LinkStatus[k] > CLOSED) { marked[j] = 1; m++; diff --git a/src/rules.c b/src/rules.c index aa2e003..7ad038f 100755 --- a/src/rules.c +++ b/src/rules.c @@ -704,7 +704,7 @@ int checkstatus(struct Premise *p) case IS_OPEN: case IS_CLOSED: case IS_ACTIVE: - i = S[p->index]; + i = LinkStatus[p->index]; if (i <= CLOSED) j = IS_CLOSED; else if (i == ACTIVE) j = IS_ACTIVE; else j = IS_OPEN; @@ -736,20 +736,20 @@ int checkvalue(struct Premise *p) /*** Updated 10/25/00 ***/ case r_DEMAND: if (p->object == r_SYSTEM) x = Dsystem*Ucf[DEMAND]; - else x = D[i]*Ucf[DEMAND]; + else x = NodeDemand[i]*Ucf[DEMAND]; break; case r_HEAD: - case r_GRADE: x = H[i]*Ucf[HEAD]; + case r_GRADE: x = NodeHead[i]*Ucf[HEAD]; break; - case r_PRESSURE: x = (H[i]-Node[i].El)*Ucf[PRESSURE]; + case r_PRESSURE: x = (NodeHead[i]-Node[i].El)*Ucf[PRESSURE]; break; - case r_LEVEL: x = (H[i]-Node[i].El)*Ucf[HEAD]; + case r_LEVEL: x = (NodeHead[i]-Node[i].El)*Ucf[HEAD]; break; case r_FLOW: x = ABS(Q[i])*Ucf[FLOW]; break; - case r_SETTING: if (K[i] == MISSING) return(0); - x = K[i]; + case r_SETTING: if (LinkSetting[i] == MISSING) return(0); + x = LinkSetting[i]; switch (Link[i].Type) { case PRV: @@ -761,14 +761,14 @@ int checkvalue(struct Premise *p) case r_FILLTIME: if (i <= Njuncs) return(0); j = i-Njuncs; if (Tank[j].A == 0.0) return(0); - if (D[i] <= TINY) return(0); - x = (Tank[j].Vmax - Tank[j].V)/D[i]; + if (NodeDemand[i] <= TINY) return(0); + x = (Tank[j].Vmax - Tank[j].V)/NodeDemand[i]; break; case r_DRAINTIME: if (i <= Njuncs) return(0); j = i-Njuncs; if (Tank[j].A == 0.0) return(0); - if (D[i] >= -TINY) return(0); - x = (Tank[j].Vmin - Tank[j].V)/D[i]; + if (NodeDemand[i] >= -TINY) return(0); + x = (Tank[j].Vmin - Tank[j].V)/NodeDemand[i]; break; default: return(0); } @@ -875,21 +875,21 @@ int takeactions() flag = FALSE; a = item->action; k = a->link; - s = S[k]; - v = K[k]; + s = LinkStatus[k]; + v = LinkSetting[k]; x = a->setting; /* Switch link from closed to open */ if (a->status == IS_OPEN && s <= CLOSED) { - setlinkstatus(k, 1, &S[k], &K[k]); + setlinkstatus(k, 1, &LinkStatus[k], &LinkSetting[k]); flag = TRUE; } /* Switch link from not closed to closed */ else if (a->status == IS_CLOSED && s > CLOSED) { - setlinkstatus(k, 0, &S[k], &K[k]); + setlinkstatus(k, 0, &LinkStatus[k], &LinkSetting[k]); flag = TRUE; } @@ -905,7 +905,7 @@ int takeactions() } if (ABS(x-v) > tol) { - setlinksetting(k, x, &S[k], &K[k]); + setlinksetting(k, x, &LinkStatus[k], &LinkSetting[k]); flag = TRUE; } } diff --git a/src/vars.h b/src/vars.h index 4378123..59b92ad 100755 --- a/src/vars.h +++ b/src/vars.h @@ -144,17 +144,17 @@ AUTHOR: L. Rossman SField Field[MAXVAR]; /* Output reporting fields */ /* Array pointers not allocated and freed in same routine */ - char *S, /* Link status */ + char *LinkStatus, /* Link status */ *OldStat; /* Previous link/tank status */ - double *D, /* Node actual demand */ - *C, /* Node actual quality */ + double *NodeDemand, /* Node actual demand */ + *NodeQual, /* Node actual quality */ *E, /* Emitter flows */ - *K, /* Link settings */ + *LinkSetting, /* Link settings */ *Q, /* Link flows */ - *R, /* Pipe reaction rate */ + *PipeRateCoeff, /* Pipe reaction rate */ *X, /* General purpose array */ - *XC; /* General purpose array */ -EXTERN double *H; /* Node heads */ + *TempQual; /* General purpose array for water quality */ +EXTERN double *NodeHead; /* Node heads */ EXTERN double *QTankVolumes; EXTERN double *QLinkFlow; EXTERN STmplist *Patlist; /* Temporary time pattern list */ @@ -171,7 +171,8 @@ EXTERN Svalve *Valve; /* Valve data */ EXTERN Scontrol *Control; /* Control data */ EXTERN HTtable *Nht, *Lht; /* Hash tables for ID labels */ EXTERN Padjlist *Adjlist; /* Node adjacency lists */ -EXTERN int _relativeError, _iterations; /* Info about hydraulic solution */ +EXTERN double _relativeError; +EXTERN int _iterations; /* Info about hydraulic solution */ /* ** NOTE: Hydraulic analysis of the pipe network at a given point in time diff --git a/test/Net3.inp b/test/Net3.inp old mode 100644 new mode 100755 diff --git a/test/simplenet.inp b/test/simplenet.inp old mode 100644 new mode 100755