JIANG/EPANET
1
0
Fork 0
Code Issues Pull Requests Actions 4 Packages Projects Releases Wiki Activity

Merge pull request #7 from samhatchett/next

Browse Source 
Next
This commit is contained in:
Sam Hatchett
2013-09-30 13:30:18 -07:00
parent 9413dfe69f 10a5ae4cf6
commit c3adaa79e9
17 changed files with 451 additions and 392 deletions
Download Patch File Download Diff File Expand all files Collapse all files

0
.gitignore vendored Normal file → Executable file
View File

0
README.md Normal file → Executable file
View File

0
build/MSVS/LemonTigerJ.sln Normal file → Executable file
View File

0
build/MSVS/LemonTigerJ.vcxproj Normal file → Executable file
View File

20
build/Xcode/epanet/epanet.xcodeproj/project.pbxproj Normal file → Executable file
View File

@@ -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;

0
build/Xcode/epanet/epanet.xcodeproj/project.xcworkspace/contents.xcworkspacedata generated Normal file → Executable file
View File

143
include/epanet2.h Normal file → Executable file
View File

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

116
src/epanet.c
View File

@@ -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; n<demandIdx && d != NULL; d=d->next) n++;
if(n!=demandIdx) return(253);
d->Base = baseDemand/Ucf[FLOW];
}
return 0;
}
int DLLEXPORT ENgetdemandpattern(int nodeIndex, int demandIdx, int *pattIdx)
{
Pdemand d;

234
src/hydraul.c
View File

@@ -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];

8
src/input3.c
View File

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

55
src/output.c
View File

@@ -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];

182
src/quality.c
View File

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

36
src/report.c
View File

@@ -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++;

32
src/rules.c
View File

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

17
src/vars.h
View File

@@ -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

0
test/Net3.inp Normal file → Executable file
View File

0
test/simplenet.inp Normal file → Executable file
View File