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

Merge remote-tracking branch 'origin/owa_next'

Browse Source 
# Conflicts:
#	build/Linux/Makefile
#	build/Linux/Makefile_Rev
#	src/epanet.c
#	src/output.c
#	src/toolkit.h
This commit is contained in:
Sam Hatchett
2015-06-09 13:34:13 -04:00
parent ee65efefc1 71b017e9ff
commit 2563ae62c6
47 changed files with 2945 additions and 2521 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/enumstxt.h
View File

@@ -11,6 +11,9 @@ AUTHOR: L. Rossman
**********************************************************************
*/
#ifndef ENUMSTXT_H
#define ENUMSTXT_H
char *NodeTxt[] = {t_JUNCTION,
t_RESERVOIR,
t_TANK};
@@ -133,3 +136,4 @@ char *Fldname[] = {t_ELEV, t_DEMAND, t_HEAD,
char *LogoTxt[] = {LOGO1,LOGO2,LOGO3,LOGO4,LOGO5,LOGO6,NULL};
#endif

432
src/epanet.c
View File

@@ -96,9 +96,9 @@ This module calls the following functions that reside in other modules:
writelogo()
writereport()
HASH.C
HTcreate()
HTfind()
HTfree()
ENHashTablecreate()
ENHashTableFind()
ENHashTableFree()
The macro ERRCODE(x) is defined in TYPES.H. It says if the current
value of the error code variable (errcode) is not fatal (< 100) then
@@ -107,17 +107,6 @@ execute function x and set the error code equal to its return value.
*******************************************************************************
*/
/*** New compile directives ***/ //(2.00.11 - LR)
//#define CLE /* Compile as a command line executable */
//#define SOL /* Compile as a shared object library */
//#define DLL /* Compile as a Windows DLL */
/*** Following lines are deprecated ***/ //(2.00.11 - LR)
//#ifdef DLL
//#include <windows.h>
//#include <float.h>
//#endif
/*** Need to define WINDOWS to use the getTmpName function ***/ //(2.00.12 - LR)
// --- define WINDOWS
#undef WINDOWS
@@ -137,14 +126,14 @@ execute function x and set the error code equal to its return value.
#endif
#include <math.h>
#include <float.h> //(2.00.12 - LR)
#include "hash.h"
#include "text.h"
#include "types.h"
#include "enumstxt.h"
#include "funcs.h"
#define EXTERN
#include "vars.h"
#include "toolkit.h"
#include "epanet2.h"
void (* viewprog) (char *); /* Pointer to progress viewing function */
@@ -318,6 +307,7 @@ int DLLEXPORT ENopen(char *f1, char *f2, char *f3)
/* Free temporary linked lists used for Patterns & Curves */
freeTmplist(Patlist);
freeTmplist(Curvelist);
freeTmplist(Coordlist);
/* If using previously saved hydraulics then open its file */
if (Hydflag == USE) ERRCODE(openhydfile());
@@ -754,7 +744,8 @@ int DLLEXPORT ENopenQ()
OpenQflag = FALSE;
SaveQflag = FALSE;
if (!Openflag) return(102);
if (!SaveHflag) return(104);
// !LT! todo - check for SaveHflag / set sequential/step mode
//if (!SaveHflag) return(104);
/* Open WQ solver */
ERRCODE(openqual());
@@ -974,7 +965,7 @@ int DLLEXPORT ENgetversion(int *v)
int DLLEXPORT ENgetcontrol(int cindex, int *ctype, int *lindex,
float *setting, int *nindex, float *level)
EN_API_FLOAT_TYPE *setting, int *nindex, EN_API_FLOAT_TYPE *level)
/*----------------------------------------------------------------
** Input: cindex = control index (position of control statement
** in the input file, starting from 1)
@@ -1020,9 +1011,9 @@ int DLLEXPORT ENgetcontrol(int cindex, int *ctype, int *lindex,
else if (*nindex > 0)
lvl = (Control[cindex].Grade - Node[*nindex].El)*Ucf[PRESSURE];
else
lvl = (float)Control[cindex].Time;
*setting = (float)s;
*level = (float)lvl;
lvl = (EN_API_FLOAT_TYPE)Control[cindex].Time;
*setting = (EN_API_FLOAT_TYPE)s;
*level = (EN_API_FLOAT_TYPE)lvl;
return(0);
}
@@ -1053,7 +1044,7 @@ int DLLEXPORT ENgetcount(int code, int *count)
}
int DLLEXPORT ENgetoption(int code, float *value)
int DLLEXPORT ENgetoption(int code, EN_API_FLOAT_TYPE *value)
/*----------------------------------------------------------------
** Input: code = option code (see TOOLKIT.H)
** Output: *value = option value
@@ -1063,7 +1054,7 @@ int DLLEXPORT ENgetoption(int code, float *value)
*/
{
double v = 0.0;
*value = 0.0f;
*value = 0.0;
if (!Openflag) return(102);
switch (code)
{
@@ -1079,7 +1070,7 @@ int DLLEXPORT ENgetoption(int code, float *value)
break;
default: return(251);
}
*value = (float)v;
*value = (EN_API_FLOAT_TYPE)v;
return(0);
}
@@ -1196,7 +1187,7 @@ int DLLEXPORT ENgetpatternlen(int index, int *len)
}
int DLLEXPORT ENgetpatternvalue(int index, int period, float *value)
int DLLEXPORT ENgetpatternvalue(int index, int period, EN_API_FLOAT_TYPE *value)
/*----------------------------------------------------------------
** Input: index = index of time pattern
** period = pattern time period
@@ -1206,11 +1197,11 @@ int DLLEXPORT ENgetpatternvalue(int index, int period, float *value)
** and pattern
**----------------------------------------------------------------
*/
{ *value = 0.0f;
{ *value = 0.0;
if (!Openflag) return(102);
if (index < 1 || index > Npats) return(205);
if (period < 1 || period > Pattern[index].Length) return(251);
*value = (float)Pattern[index].F[period-1];
*value = (EN_API_FLOAT_TYPE)Pattern[index].F[period-1];
return(0);
}
@@ -1233,6 +1224,19 @@ int DLLEXPORT ENgetqualtype(int *qualcode, int *tracenode)
return(0);
}
int DLLEXPORT ENgetqualinfo(int *qualcode, char *chemname, char *chemunits, int *tracenode)
{
ENgetqualtype(qualcode, tracenode);
if (Qualflag == TRACE) {
strncpy(chemname, "", MAXID);
strncpy(chemunits, "dimensionless", MAXID);
}
else {
strncpy(chemname,ChemName,MAXID);
strncpy(chemunits,ChemUnits,MAXID);
}
return 0;
}
int DLLEXPORT ENgeterror(int errcode, char *errmsg, int n)
/*----------------------------------------------------------------
@@ -1258,7 +1262,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
@@ -1269,10 +1273,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;
@@ -1345,7 +1349,22 @@ int DLLEXPORT ENgetnodetype(int index, int *code)
}
int DLLEXPORT ENgetnodevalue(int index, int code, float *value)
int DLLEXPORT ENgetcoord(int index, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y)
/*----------------------------------------------------------------
** Input: index = node index
** Output: *x = value of node's coordinate
** *x = value of node's coordinate
** Returns: error code
** Purpose: retrieves coordinate x, y for a node
**----------------------------------------------------------------
*/
{
*x = Coord[index].X[0];
*y = Coord[index].Y[0];
return 0;
}
int DLLEXPORT ENgetnodevalue(int index, int code, EN_API_FLOAT_TYPE *value)
/*----------------------------------------------------------------
** Input: index = node index
** code = node parameter code (see TOOLKIT.H)
@@ -1360,7 +1379,7 @@ int DLLEXPORT ENgetnodevalue(int index, int code, float *value)
Psource source;
/* Check for valid arguments */
*value = 0.0f;
*value = 0.0;
if (!Openflag) return(102);
if (index <= 0 || index > Nnodes) return(203);
@@ -1438,19 +1457,19 @@ int DLLEXPORT ENgetnodevalue(int index, int code, float *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)
@@ -1487,7 +1506,7 @@ int DLLEXPORT ENgetnodevalue(int index, int code, float *value)
v = (Tank[index-Njuncs].Hmin - Node[index].El) * Ucf[ELEV];
}
break;
case EN_MAXLEVEL:
v = 0.0;
if ( index > Njuncs )
@@ -1513,12 +1532,12 @@ int DLLEXPORT ENgetnodevalue(int index, int code, float *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);
}
*value = (float)v;
*value = (EN_API_FLOAT_TYPE)v;
return(0);
}
@@ -1603,7 +1622,7 @@ int DLLEXPORT ENgetlinknodes(int index, int *node1, int *node2)
}
int DLLEXPORT ENgetlinkvalue(int index, int code, float *value)
int DLLEXPORT ENgetlinkvalue(int index, int code, EN_API_FLOAT_TYPE *value)
/*------------------------------------------------------------------
** Input: index = link index
** code = link parameter code (see TOOLKIT.H)
@@ -1616,7 +1635,7 @@ int DLLEXPORT ENgetlinkvalue(int index, int code, float *value)
double a,h,q, v = 0.0;
/* Check for valid arguments */
*value = 0.0f;
*value = 0.0;
if (!Openflag) return(102);
if (index <= 0 || index > Nlinks) return(204);
@@ -1643,9 +1662,12 @@ int DLLEXPORT ENgetlinkvalue(int index, int code, float *value)
break;
case EN_MINORLOSS:
v = Link[index].Km;
if (Link[index].Type != PUMP)
{
v = Link[index].Km;
v *= (SQR(Link[index].Diam)*SQR(Link[index].Diam)/0.02517);
}
else v = 0.0;
break;
case EN_INITSTATUS:
@@ -1677,8 +1699,7 @@ int DLLEXPORT ENgetlinkvalue(int index, int code, float *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;
@@ -1686,7 +1707,7 @@ int DLLEXPORT ENgetlinkvalue(int index, int code, float *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
{
@@ -1699,26 +1720,31 @@ int DLLEXPORT ENgetlinkvalue(int index, int code, float *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:
@@ -1743,12 +1769,12 @@ int DLLEXPORT ENgetlinkvalue(int index, int code, float *value)
default: return(251);
}
*value = (float)v;
*value = (EN_API_FLOAT_TYPE)v;
return(0);
}
int DLLEXPORT ENgetcurve(int curveIndex, int *nValues, float **xValues, float **yValues) // !sph
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
@@ -1764,17 +1790,17 @@ int DLLEXPORT ENgetcurve(int curveIndex, int *nValues, float **xValues, float *
Scurve curve = Curve[curveIndex];
int nPoints = curve.Npts;
float *pointX = calloc(nPoints, sizeof(float));
float *pointY = calloc(nPoints, sizeof(float));
EN_API_FLOAT_TYPE *pointX = calloc(nPoints, sizeof(EN_API_FLOAT_TYPE));
EN_API_FLOAT_TYPE *pointY = calloc(nPoints, sizeof(EN_API_FLOAT_TYPE));
int iPoint;
for (iPoint = 0; iPoint < nPoints; iPoint++) {
double x = curve.X[iPoint] * Ucf[LENGTH];
double y = curve.Y[iPoint] * Ucf[VOLUME];
pointX[iPoint] = (float)x;
pointY[iPoint] = (float)y;
pointX[iPoint] = (EN_API_FLOAT_TYPE)x;
pointY[iPoint] = (EN_API_FLOAT_TYPE)y;
}
strncpy(id, curve.ID, MAXID);
*nValues = nPoints;
*xValues = pointX;
*yValues = pointY;
@@ -1782,6 +1808,7 @@ int DLLEXPORT ENgetcurve(int curveIndex, int *nValues, float **xValues, float *
return err;
}
/*
----------------------------------------------------------------
Functions for changing network data
@@ -1790,7 +1817,7 @@ int DLLEXPORT ENgetcurve(int curveIndex, int *nValues, float **xValues, float *
int DLLEXPORT ENsetcontrol(int cindex, int ctype, int lindex,
float setting, int nindex, float level)
EN_API_FLOAT_TYPE setting, int nindex, EN_API_FLOAT_TYPE level)
/*----------------------------------------------------------------
** Input: cindex = control index (position of control statement
** in the input file, starting from 1)
@@ -1878,7 +1905,7 @@ int DLLEXPORT ENsetcontrol(int cindex, int ctype, int lindex,
}
int DLLEXPORT ENsetnodevalue(int index, int code, float v)
int DLLEXPORT ENsetnodevalue(int index, int code, EN_API_FLOAT_TYPE v)
/*----------------------------------------------------------------
** Input: index = node index
** code = node parameter code (see TOOLKIT.H)
@@ -1908,7 +1935,7 @@ int DLLEXPORT ENsetnodevalue(int index, int code, float v)
Tank[j].Hmin += value;
Tank[j].Hmax += value;
Node[index].El += value;
H[index] += value;
NodeHead[index] += value;
}
break;
@@ -1965,14 +1992,16 @@ int DLLEXPORT ENsetnodevalue(int index, int code, float 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);
if (j < 0 || j > Npats) return(205);
source->Pat = j;
}
else
else // code == EN_SOURCETYPE
{
j = ROUND(value);
if ( j < CONCEN || j > FLOWPACED) return(251);
@@ -1989,7 +2018,7 @@ int DLLEXPORT ENsetnodevalue(int index, int code, float 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
{
@@ -1998,7 +2027,9 @@ int DLLEXPORT ENsetnodevalue(int index, int code, float 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;
@@ -2087,7 +2118,7 @@ int DLLEXPORT ENsetnodevalue(int index, int code, float v)
}
int DLLEXPORT ENsetlinkvalue(int index, int code, float v)
int DLLEXPORT ENsetlinkvalue(int index, int code, EN_API_FLOAT_TYPE v)
/*----------------------------------------------------------------
** Input: index = link index
** code = link parameter code (see TOOLKIT.H)
@@ -2153,7 +2184,7 @@ int DLLEXPORT ENsetlinkvalue(int index, int code, float 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:
@@ -2180,7 +2211,7 @@ int DLLEXPORT ENsetlinkvalue(int index, int code, float 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;
@@ -2274,7 +2305,7 @@ int DLLEXPORT ENaddpattern(char *id)
}
int DLLEXPORT ENsetpattern(int index, float *f, int n)
int DLLEXPORT ENsetpattern(int index, EN_API_FLOAT_TYPE *f, int n)
/*----------------------------------------------------------------
** Input: index = time pattern index
** *f = array of pattern multipliers
@@ -2303,7 +2334,7 @@ int DLLEXPORT ENsetpattern(int index, float *f, int n)
}
int DLLEXPORT ENsetpatternvalue(int index, int period, float value)
int DLLEXPORT ENsetpatternvalue(int index, int period, EN_API_FLOAT_TYPE value)
/*----------------------------------------------------------------
** Input: index = time pattern index
** period = time pattern period
@@ -2334,62 +2365,93 @@ int DLLEXPORT ENsettimeparam(int code, long value)
{
if (!Openflag) return(102);
if (OpenHflag || OpenQflag) {
// --> there's nothing wrong with changing certain time parameters during a simulation run
if (code != EN_DURATION) {
// --> there's nothing wrong with changing certain time parameters during a simulation run, or before the run has started.
// todo -- how to tell?
/*
if (code == EN_DURATION || code == EN_HTIME || code == EN_REPORTSTEP || code == EN_DURATION || Htime == 0) {
// it's ok
}
else {
return(109);
}
*/
}
if (value < 0) return(202);
switch(code)
{
case EN_DURATION: Dur = value;
if (Rstart > Dur) Rstart = 0;
break;
case EN_HYDSTEP: if (value == 0) return(202);
Hstep = value;
Hstep = MIN(Pstep, Hstep);
Hstep = MIN(Rstep, Hstep);
Qstep = MIN(Qstep, Hstep);
break;
case EN_QUALSTEP: if (value == 0) return(202);
Qstep = value;
Qstep = MIN(Qstep, Hstep);
break;
case EN_PATTERNSTEP: if (value == 0) return(202);
Pstep = value;
if (Hstep > Pstep) Hstep = Pstep;
break;
case EN_PATTERNSTART: Pstart = value;
break;
case EN_REPORTSTEP: if (value == 0) return(202);
Rstep = value;
if (Hstep > Rstep) Hstep = Rstep;
break;
case EN_REPORTSTART: if (Rstart > Dur) return(202);
Rstart = value;
break;
case EN_RULESTEP: if (value == 0) return(202);
Rulestep = value;
Rulestep = MIN(Rulestep, Hstep);
break;
case EN_STATISTIC: if (value > RANGE) return(202);
Tstatflag = (char)value;
break;
case EN_HTIME: Htime = value;
break;
default: return(251);
switch(code)
{
case EN_DURATION:
Dur = value;
if (Rstart > Dur) Rstart = 0;
break;
case EN_HYDSTEP:
if (value == 0) return(202);
Hstep = value;
Hstep = MIN(Pstep, Hstep);
Hstep = MIN(Rstep, Hstep);
Qstep = MIN(Qstep, Hstep);
break;
case EN_QUALSTEP:
if (value == 0) return(202);
Qstep = value;
Qstep = MIN(Qstep, Hstep);
break;
case EN_PATTERNSTEP:
if (value == 0) return(202);
Pstep = value;
if (Hstep > Pstep) Hstep = Pstep;
break;
case EN_PATTERNSTART:
Pstart = value;
break;
case EN_REPORTSTEP:
if (value == 0) return(202);
Rstep = value;
if (Hstep > Rstep) Hstep = Rstep;
break;
case EN_REPORTSTART:
if (Rstart > Dur) return(202);
Rstart = value;
break;
case EN_RULESTEP:
if (value == 0) return(202);
Rulestep = value;
Rulestep = MIN(Rulestep, Hstep);
break;
case EN_STATISTIC:
if (value > RANGE) return(202);
Tstatflag = (char)value;
break;
case EN_HTIME:
Htime = value;
break;
case EN_QTIME:
Qtime = value;
break;
default:
return(251);
}
return(0);
}
int DLLEXPORT ENsetoption(int code, float v)
int DLLEXPORT ENsetoption(int code, EN_API_FLOAT_TYPE v)
/*----------------------------------------------------------------
** Input: code = option code (see TOOLKIT.H)
** v = option value
** Output: none
** Returns: error code
** Purpose: sets value for an analysis option
** Returns: error code
** Purpose: sets value for an analysis option
**----------------------------------------------------------------
*/
{
@@ -2751,13 +2813,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;
@@ -2768,10 +2830,12 @@ void initpointers()
Pattern = NULL;
Curve = NULL;
Control = NULL;
Coord = NULL;
X = NULL;
Patlist = NULL;
Curvelist = NULL;
Coordlist = NULL;
Adjlist = NULL;
Aii = NULL;
Aij = NULL;
@@ -2784,8 +2848,8 @@ void initpointers()
XLNZ = NULL;
NZSUB = NULL;
LNZ = NULL;
Nht = NULL;
Lht = NULL;
NodeHashTable = NULL;
LinkHashTable = NULL;
initrules();
}
@@ -2803,10 +2867,10 @@ int allocdata()
int errcode = 0;
/* Allocate node & link ID hash tables */
Nht = HTcreate();
Lht = HTcreate();
ERRCODE(MEMCHECK(Nht));
ERRCODE(MEMCHECK(Lht));
NodeHashTable = ENHashTableCreate();
LinkHashTable = ENHashTableCreate();
ERRCODE(MEMCHECK(NodeHashTable));
ERRCODE(MEMCHECK(LinkHashTable));
/* Allocate memory for network nodes */
/*************************************************************
@@ -2818,13 +2882,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 */
@@ -2833,12 +2897,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, */
@@ -2851,12 +2915,14 @@ int allocdata()
Control = (Scontrol *) calloc(MaxControls+1,sizeof(Scontrol));
Pattern = (Spattern *) calloc(MaxPats+1, sizeof(Spattern));
Curve = (Scurve *) calloc(MaxCurves+1, sizeof(Scurve));
Coord = (Scoord *) calloc(MaxNodes+1, sizeof(Scoord));
ERRCODE(MEMCHECK(Tank));
ERRCODE(MEMCHECK(Pump));
ERRCODE(MEMCHECK(Valve));
ERRCODE(MEMCHECK(Control));
ERRCODE(MEMCHECK(Pattern));
ERRCODE(MEMCHECK(Curve));
ERRCODE(MEMCHECK(Coord));
}
/* Initialize pointers used in patterns, curves, and demand category lists */
@@ -2874,7 +2940,19 @@ int allocdata()
Curve[n].X = NULL;
Curve[n].Y = NULL;
}
for (n=0; n<=MaxNodes; n++) Node[n].D = NULL;
for (n=0; n<=MaxNodes; n++)
{
// node demand
Node[n].D = NULL;
/* Allocate memory for coord data */
Coord[n].X = (double *) calloc(1, sizeof(double));
Coord[n].Y = (double *) calloc(1, sizeof(double));
if (Coord[n].X == NULL || Coord[n].Y == NULL) return(101);
Coord[n].X[0] = 0;
Coord[n].Y[0] = 0;
}
}
/* Allocate memory for rule base (see RULES.C) */
@@ -2935,12 +3013,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)
@@ -2991,8 +3069,8 @@ void freedata()
freerules();
/* Free hash table memory */
if (Nht != NULL) HTfree(Nht);
if (Lht != NULL) HTfree(Lht);
if (NodeHashTable != NULL) ENHashTableFree(NodeHashTable);
if (LinkHashTable != NULL) ENHashTableFree(LinkHashTable);
}
@@ -3106,7 +3184,7 @@ int findnode(char *id)
**----------------------------------------------------------------
*/
{
return(HTfind(Nht,id));
return(ENHashTableFind(NodeHashTable,id));
}
@@ -3119,7 +3197,7 @@ int findlink(char *id)
**----------------------------------------------------------------
*/
{
return(HTfind(Lht,id));
return(ENHashTableFind(LinkHashTable,id));
}
@@ -3181,6 +3259,8 @@ char *geterrmsg(int errcode)
case 307: strcpy(Msg,ERR307); break;
case 308: strcpy(Msg,ERR308); break;
case 309: strcpy(Msg,ERR309); break;
case 401: strcpy(Msg,ERR401); break;
default: strcpy(Msg,"");
}
return(Msg);
@@ -3251,22 +3331,43 @@ int DLLEXPORT ENgetnumdemands(int nodeIndex, int *numDemands)
*numDemands=n;
return 0;
}
int DLLEXPORT ENgetbasedemand(int nodeIndex, int demandIdx, float *baseDemand)
int DLLEXPORT ENgetbasedemand(int nodeIndex, int demandIdx, EN_API_FLOAT_TYPE *baseDemand)
{
Pdemand d;
int n=0;
/* Check for valid arguments */
if (!Openflag) return(102);
if (nodeIndex <= 0 || nodeIndex > Nnodes) return(203);
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);
*baseDemand=d->Base*Ucf[FLOW];
return 0;
*baseDemand=(EN_API_FLOAT_TYPE)(d->Base*Ucf[FLOW]);
}
else {
*baseDemand=(EN_API_FLOAT_TYPE)(0.0);
}
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;
int n=0;
int n=1;
/* Check for valid arguments */
if (!Openflag) return(102);
if (nodeIndex <= 0 || nodeIndex > Nnodes) return(203);
@@ -3276,5 +3377,28 @@ int DLLEXPORT ENgetdemandpattern(int nodeIndex, int demandIdx, int *pattIdx)
return 0;
}
int DLLEXPORT ENgetaveragepatternvalue(int index, EN_API_FLOAT_TYPE *value)
/*----------------------------------------------------------------
** Input: index = index of time pattern
** period = pattern time period
** Output: *value = pattern multiplier
** Returns: error code
** Purpose: retrieves multiplier for a specific time period
** and pattern
**----------------------------------------------------------------
*/
{ *value = 0.0;
if (!Openflag) return(102);
if (index < 1 || index > Npats) return(205);
//if (period < 1 || period > Pattern[index].Length) return(251);
int i;
for (i=0; i<Pattern[index].Length; i++) {
*value+=Pattern[index].F[i];
}
*value/=(EN_API_FLOAT_TYPE)Pattern[index].Length;
return(0);
}
/*************************** END OF EPANET.C ***************************/

9
src/funcs.h
View File

@@ -25,6 +25,10 @@ AUTHOR: L. Rossman
** NOTE: The exportable functions that can be called
** via the DLL are prototyped in TOOLKIT.H.
*/
#ifndef FUNCS_H
#define FUNCS_H
void initpointers(void); /* Initializes pointers */
int allocdata(void); /* Allocates memory */
void freeTmplist(STmplist *); /* Frees items in linked list */
@@ -61,11 +65,13 @@ int addnodeID(int, char *); /* Adds node ID to data base */
int addlinkID(int, char *); /* Adds link ID to data base */
int addpattern(char *); /* Adds pattern to data base */
int addcurve(char *); /* Adds curve to data base */
int addcoord(char *); /* Adds coord to data base */
STmplist *findID(char *, STmplist *); /* Locates ID on linked list */
int unlinked(void); /* Checks for unlinked nodes */
int getpumpparams(void); /* Computes pump curve coeffs.*/
int getpatterns(void); /* Gets pattern data from list*/
int getcurves(void); /* Gets curve data from list */
int getcoords(void); /* Gets coordinate data from list */
int findmatch(char *,char *[]); /* Finds keyword in line */
int match(char *, char *); /* Checks for word match */
int gettokens(char *); /* Tokenizes input line */
@@ -82,6 +88,7 @@ int pumpdata(void); /* Processes pump data */
int valvedata(void); /* Processes valve data */
int patterndata(void); /* Processes pattern data */
int curvedata(void); /* Processes curve data */
int coordata(void); /* Processes coordinate data */
int demanddata(void); /* Processes demand data */
int controldata(void); /* Processes simple controls */
int energydata(void); /* Processes energy data */
@@ -279,3 +286,5 @@ int saveepilog(void); /* Saves output file epilog */
/* ------------ INPFILE.C --------------*/
int saveinpfile(char *); /* Saves network to text file */
#endif

186
src/hash.c
View File

@@ -1,22 +1,22 @@
/*-----------------------------------------------------------------------------
** hash.c
**
** Implementation of a simple Hash Table for string storage & retrieval
**
** Written by L. Rossman
** Last Updated on 6/19/03
**
** The hash table data structure (HTable) is defined in "hash.h".
** Interface Functions:
** HTcreate() - creates a hash table
** HTinsert() - inserts a string & its index value into a hash table
** HTfind() - retrieves the index value of a string from a table
** HTfree() - frees a hash table
**
*********************************************************************
** NOTE: This is a modified version of the original HASH.C module.
*********************************************************************
*/
** hash.c
**
** Implementation of a simple Hash Table for string storage & retrieval
**
** Written by L. Rossman
** Last Updated on 6/19/03
**
** The hash table data structure (HTable) is defined in "hash.h".
** Interface Functions:
** HTcreate() - creates a hash table
** HTinsert() - inserts a string & its index value into a hash table
** HTfind() - retrieves the index value of a string from a table
** HTfree() - frees a hash table
**
*********************************************************************
** NOTE: This is a modified version of the original HASH.C module.
*********************************************************************
*/
#ifndef __APPLE__
#include <malloc.h>
@@ -26,89 +26,97 @@
#include <string.h>
#include "hash.h"
/* Use Fletcher's checksum to compute 2-byte hash of string */
unsigned int hash(char *str)
unsigned int _enHash(char *str);
unsigned int _enHash(char *str)
{
unsigned int sum1= 0, check1;
unsigned long sum2= 0L;
while( '\0' != *str )
{
sum1 += (*str);
str++;
if ( 255 <= sum1 ) sum1 -= 255;
sum2 += sum1;
unsigned int hash = 5381;
int c;
while ((c = *str++)) {
hash = ((hash << 5) + hash) + c; /* hash * 33 + c */
}
unsigned int retHash = hash % ENHASHTABLEMAXSIZE;
return retHash;
}
ENHashTable *ENHashTableCreate()
{
int i;
ENHashTable *ht = (ENHashTable *) calloc(ENHASHTABLEMAXSIZE, sizeof(ENHashTable));
if (ht != NULL) {
for (i=0; i<ENHASHTABLEMAXSIZE; i++) {
ht[i] = NULL;
}
check1= sum2;
check1 %= 255;
check1= 255 - (sum1+check1) % 255;
sum1= 255 - (sum1+check1) % 255;
return( ( ( check1 << 8 ) | sum1 ) % HTMAXSIZE);
}
return(ht);
}
HTtable *HTcreate()
int ENHashTableInsert(ENHashTable *ht, char *key, int data)
{
int i;
HTtable *ht = (HTtable *) calloc(HTMAXSIZE, sizeof(HTtable));
if (ht != NULL) for (i=0; i<HTMAXSIZE; i++) ht[i] = NULL;
return(ht);
unsigned int i = _enHash(key);
ENHashEntry *entry;
if ( i >= ENHASHTABLEMAXSIZE ) {
return(0);
}
entry = (ENHashEntry *) malloc(sizeof(ENHashEntry));
if (entry == NULL) {
return(0);
}
entry->key = key;
entry->data = data;
entry->next = ht[i];
ht[i] = entry;
return(1);
}
int HTinsert(HTtable *ht, char *key, int data)
int ENHashTableFind(ENHashTable *ht, char *key)
{
unsigned int i = hash(key);
struct HTentry *entry;
if ( i >= HTMAXSIZE ) return(0);
entry = (struct HTentry *) malloc(sizeof(struct HTentry));
if (entry == NULL) return(0);
entry->key = key;
entry->data = data;
entry->next = ht[i];
ht[i] = entry;
return(1);
unsigned int i = _enHash(key);
ENHashEntry *entry;
if ( i >= ENHASHTABLEMAXSIZE ) {
return(NOTFOUND);
}
entry = ht[i];
while (entry != NULL)
{
if ( strcmp(entry->key,key) == 0 ) {
return(entry->data);
}
entry = entry->next;
}
return(NOTFOUND);
}
int HTfind(HTtable *ht, char *key)
char *ENHashTableFindKey(ENHashTable *ht, char *key)
{
unsigned int i = hash(key);
struct HTentry *entry;
if ( i >= HTMAXSIZE ) return(NOTFOUND);
entry = ht[i];
while (entry != NULL)
{
if ( strcmp(entry->key,key) == 0 ) return(entry->data);
entry = entry->next;
}
return(NOTFOUND);
unsigned int i = _enHash(key);
ENHashEntry *entry;
if ( i >= ENHASHTABLEMAXSIZE ) {
return(NULL);
}
entry = ht[i];
while (entry != NULL)
{
if ( strcmp(entry->key,key) == 0 ) {
return(entry->key);
}
entry = entry->next;
}
return(NULL);
}
char *HTfindKey(HTtable *ht, char *key)
void ENHashTableFree(ENHashTable *ht)
{
unsigned int i = hash(key);
struct HTentry *entry;
if ( i >= HTMAXSIZE ) return(NULL);
entry = ht[i];
while (entry != NULL)
{
if ( strcmp(entry->key,key) == 0 ) return(entry->key);
entry = entry->next;
}
return(NULL);
}
void HTfree(HTtable *ht)
{
struct HTentry *entry,
*nextentry;
int i;
for (i=0; i<HTMAXSIZE; i++)
{
entry = ht[i];
while (entry != NULL)
{
nextentry = entry->next;
free(entry);
entry = nextentry;
}
}
free(ht);
ENHashEntry *entry, *nextentry;
int i;
for (i=0; i<ENHASHTABLEMAXSIZE; i++)
{
entry = ht[i];
while (entry != NULL)
{
nextentry = entry->next;
free(entry);
entry = nextentry;
}
}
free(ht);
}

24
src/hash.h
View File

@@ -4,21 +4,25 @@
**
*/
#define HTMAXSIZE 1999
#ifndef HASH_H
#define HASH_H
#define ENHASHTABLEMAXSIZE 128000
#define NOTFOUND 0
struct HTentry
typedef struct HTentryStruct
{
char *key;
int data;
struct HTentry *next;
};
struct HTentryStruct *next;
} ENHashEntry;
typedef struct HTentry *HTtable;
typedef ENHashEntry *ENHashTable;
HTtable *HTcreate(void);
int HTinsert(HTtable *, char *, int);
int HTfind(HTtable *, char *);
char *HTfindKey(HTtable *, char *);
void HTfree(HTtable *);
ENHashTable *ENHashTableCreate(void);
int ENHashTableInsert(ENHashTable *, char *, int);
int ENHashTableFind(ENHashTable *, char *);
char *ENHashTableFindKey(ENHashTable *, char *);
void ENHashTableFree(ENHashTable *);
#endif

240
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 */
@@ -259,6 +259,9 @@ int nexthyd(long *tstep)
else
{
Htime++; /* Force completion of analysis */
if (OpenQflag) {
Qtime++; // force completion of wq analysis too
}
}
*tstep = hydstep;
return(errcode);
@@ -372,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;
}
@@ -399,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);
@@ -603,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. */
@@ -616,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];
}
}
}
@@ -629,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 */
@@ -661,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)
@@ -686,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++;
@@ -761,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)
{
@@ -799,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 &&
@@ -835,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;
}
@@ -951,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 */
@@ -998,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;
@@ -1008,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)
@@ -1050,15 +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) Tank[i].V = Tank[i].Vmax;
if (Tank[i].V - D[n] <= Tank[i].Vmin) Tank[i].V = Tank[i].Vmin;
H[n] = tankgrade(i,Tank[i].V);
if (Tank[i].V + NodeDemand[n] >= Tank[i].Vmax) {
Tank[i].V = Tank[i].Vmax;
}
else if (Tank[i].V - NodeDemand[n] <= Tank[i].Vmin) {
Tank[i].V = Tank[i].Vmin;
}
NodeHead[n] = tankgrade(i,Tank[i].V);
}
} /* End of tanklevels */
@@ -1182,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;
@@ -1238,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 */
@@ -1268,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);
}
}
@@ -1307,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;
@@ -1336,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;
}
}
@@ -1369,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)
@@ -1443,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)
@@ -1471,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)
@@ -1521,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)
@@ -1574,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);
}
@@ -1609,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 */
@@ -1675,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;
}
@@ -1686,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]);
@@ -1735,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;
@@ -1755,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)
@@ -1774,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];
}
}
@@ -1855,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;
@@ -1871,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 */
@@ -1898,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 */
@@ -1919,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 */
@@ -1986,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));
}
@@ -2013,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];
@@ -2142,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)
@@ -2160,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;
@@ -2169,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) */
@@ -2235,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 */
@@ -2247,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). */
@@ -2267,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 */
@@ -2300,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)
@@ -2327,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
@@ -2373,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
@@ -2417,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];
@@ -2426,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;
@@ -2468,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];

139
src/input2.c
View File

@@ -48,6 +48,7 @@ char *Tok[MAXTOKS]; /* Array of token strings */
/* Used in INPUT3.C: */
STmplist *PrevPat; /* Pointer to pattern list element */
STmplist *PrevCurve; /* Pointer to curve list element */
STmplist *PrevCoord; /* Pointer to coordinate list element */
/* Defined in enumstxt.h in EPANET.C */
extern char *SectTxt[]; /* Input section keywords */
@@ -78,7 +79,8 @@ int netsize()
MaxRules = 0;
MaxCurves = 0;
sect = -1;
MaxCoords = 0;
/* Add a default pattern 0 */
MaxPats = -1;
addpattern("");
@@ -106,20 +108,22 @@ int netsize()
/* Add to count of current component */
switch(sect)
{
case _JUNCTIONS: MaxJuncs++; break;
case _RESERVOIRS:
case _TANKS: MaxTanks++; break;
case _PIPES: MaxPipes++; break;
case _PUMPS: MaxPumps++; break;
case _VALVES: MaxValves++; break;
case _CONTROLS: MaxControls++; break;
case _RULES: addrule(tok); break; /* See RULES.C */
case _PATTERNS: errcode = addpattern(tok);
break;
case _CURVES: errcode = addcurve(tok);
break;
}
{
case _JUNCTIONS: MaxJuncs++; break;
case _RESERVOIRS:
case _TANKS: MaxTanks++; break;
case _PIPES: MaxPipes++; break;
case _PUMPS: MaxPumps++; break;
case _VALVES: MaxValves++; break;
case _CONTROLS: MaxControls++; break;
case _RULES: addrule(tok); break; /* See RULES.C */
case _PATTERNS: errcode = addpattern(tok);
break;
case _CURVES: errcode = addcurve(tok);
break;
// case _COORDS: errcode = addcoord(tok); //06.02.2010-woohn
// break;
}
if (errcode) break;
}
@@ -172,6 +176,8 @@ int readdata()
Npats = MaxPats;
PrevPat = NULL;
PrevCurve = NULL;
PrevCoord = NULL;
sect = -1;
errsum = 0;
@@ -239,6 +245,7 @@ int readdata()
/* Get pattern & curve data from temp. lists */
if (!errcode) errcode = getpatterns();
if (!errcode) errcode = getcurves();
//if (!errcode) errcode = getcoords();
if (!errcode) errcode = getpumpparams();
/* Free input buffer */
@@ -293,7 +300,7 @@ int newline(int sect, char *line)
case _OPTIONS: return(optiondata());
/* Data in these sections are not used for any computations */
case _COORDS: return(0);
case _COORDS: return (0); //return(coordata());
case _LABELS: return(0);
case _TAGS: return(0);
case _VERTICES: return(0);
@@ -416,7 +423,7 @@ int addnodeID(int n, char *id)
{
if (findnode(id)) return(0); /* see EPANET.C */
strncpy(Node[n].ID, id, MAXID);
HTinsert(Nht, Node[n].ID, n); /* see HASH.C */
ENHashTableInsert(NodeHashTable, Node[n].ID, n); /* see HASH.C */
return(1);
}
@@ -433,7 +440,7 @@ int addlinkID(int n, char *id)
{
if (findlink(id)) return(0); /* see EPANET.C */
strncpy(Link[n].ID, id, MAXID);
HTinsert(Lht, Link[n].ID, n); /* see HASH.C */
ENHashTableInsert(LinkHashTable, Link[n].ID, n); /* see HASH.C */
return(1);
}
@@ -513,6 +520,43 @@ int addcurve(char *id)
return(0);
}
int addcoord(char *id)
/*
**-------------------------------------------------------------
** Input: id = curve ID label
** Output: returns error code
** Purpose: adds a new curve to the database
**--------------------------------------------------------------
*/
{
STmplist *c;
/* Check if ID is same as last one processed */
if (Coordlist != NULL && strcmp(id,Coordlist->ID) == 0) return(0);
/* Check that coordinate was not already created */
if (findID(id,Coordlist) == NULL)
{
/* Update coordinate count & create new list element */
(MaxCoords)++;
c = (STmplist *) malloc(sizeof(STmplist));
if (c == NULL) {
return(101);
}
else {
/* Initialize list element properties */
// c->i = MaxCoords; // bug! if coordinates are not in the same order as junctions, then this is a BAD assumption
// do this later: c->i = findnode(id);
strncpy(c->ID,id,MAXID);
c->x = NULL;
c->y = NULL;
c->next = Coordlist;
Coordlist = c;
}
}
return(0);
}
STmplist *findID(char *id, STmplist *list)
/*
@@ -705,6 +749,65 @@ int getcurves(void)
return(0);
}
int getcoords(void)
/*
**-----------------------------------------------------------
** Input: none
** Output: returns error code
** Purpose: retrieves curve data from temporary linked list
**-----------------------------------------------------------
*/
{
int i,j,n;
double x;
SFloatlist *xFloatList, *yFloatList;
STmplist *coordinateList;
/* Start at head of coordinate list */
coordinateList = Coordlist;
/* Traverse list of coordinates */
while (coordinateList != NULL)
{
// BAD! ---> i = coordinateList->i;
i = findnode(coordinateList->ID);
if (i >= 1 && i <= MaxNodes)
{
/* Save coordinate ID */
strcpy(Coord[i].ID, coordinateList->ID);
n = 1; //Coord[i].Npts
/* Traverse list of x,y data */
x = BIG;
xFloatList = coordinateList->x;
yFloatList = coordinateList->y;
j = n - 1;
while (xFloatList != NULL && yFloatList != NULL && j >= 0)
{
/* Check that x data is in ascending order */
if (xFloatList->value >= x)
{
sprintf(Msg,ERR230,coordinateList->ID);
writeline(Msg);
return(200);
}
x = xFloatList->value;
/* Save x,y data in Curve structure */
Coord[i].X[j] = xFloatList->value;
xFloatList = xFloatList->next;
Coord[i].Y[j] = yFloatList->value;
yFloatList = yFloatList->next;
j--;
}
}
coordinateList = coordinateList->next;
}
return(0);
}
int findmatch(char *line, char *keyword[])
/*

63
src/input3.c
View File

@@ -42,6 +42,8 @@ extern char *Fldname[];
extern char *Tok[MAXTOKS];
extern STmplist *PrevPat;
extern STmplist *PrevCurve;
extern STmplist *PrevCoord;
extern int Ntokens;
@@ -97,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 */
@@ -577,6 +579,59 @@ int curvedata()
return(0);
}
int coordata()
/*
**--------------------------------------------------------------
** Input: none
** Output: returns error code
** Purpose: processes coordinate data
** Format:
** [COORD]
** id x y
**--------------------------------------------------------------
*/
{
double x,y;
SFloatlist *fx, *fy;
STmplist *c;
/* Check for valid curve ID */
if (Ntokens < 3) return(201);
if (
PrevCoord != NULL &&
strcmp(Tok[0],PrevCoord->ID) == 0
) c = PrevCoord;
else c = findID(Tok[0],Coordlist);
// c = findID(Tok[0],Coordlist);
if (c == NULL) return(205);
/* Check for valid data */
if (!getfloat(Tok[1],&x)) return(202);
if (!getfloat(Tok[2],&y)) return(202);
/* Add new data point to curve's linked list */
fx = (SFloatlist *) malloc(sizeof(SFloatlist));
fy = (SFloatlist *) malloc(sizeof(SFloatlist));
if (fx == NULL || fy == NULL) return(101);
fx->value = x;
fx->next = c->x;
c->x = fx;
fy->value = y;
fy->next = c->y;
c->y = fy;
//Curve[c->i].Npts++;
/* Save the pointer to this curve */
PrevCoord = c;
return(0);
/* Save coordn data */
//Coord[Njuncs].X = x;
//Coord[Njuncs].Y = y;
} /* end of coordata */
int demanddata()
/*
@@ -627,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 ***/

7
src/mempool.h
View File

@@ -6,6 +6,10 @@
** The type alloc_handle_t provides an opaque reference to the
** alloc pool - only the alloc routines know its structure.
*/
#ifndef MEMPOOL_H
#define MEMPOOL_H
#ifndef DLLEXPORT
#ifdef DLL
#ifdef __cplusplus
@@ -24,6 +28,7 @@
#endif
#endif
typedef struct
{
long dummy;
@@ -34,3 +39,5 @@ DLLEXPORT char *Alloc(long);
DLLEXPORT alloc_handle_t *AllocSetPool(alloc_handle_t *);
DLLEXPORT void AllocReset(void);
DLLEXPORT void AllocFreePool(void);
#endif

61
src/output.c
View File

@@ -152,29 +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;
else x[i] = (REAL4)Q[i];
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);
@@ -285,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;
@@ -338,7 +339,6 @@ int saveoutput()
/* Write out node results, then link results */
for (j=DEMAND; j<=QUALITY; j++) ERRCODE(nodeoutput(j,x,Ucf[j]));
for (j=FLOW; j<=FRICTION; j++) ERRCODE(linkoutput(j,x,Ucf[j]));
free(x);
return(errcode);
} /* End of saveoutput */
@@ -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 */
@@ -380,7 +380,7 @@ int nodeoutput(int j, REAL4 *x, double ucf)
} /* End of nodeoutput */
int linkoutput(int j, float *x, double ucf)
int linkoutput(int j, REAL4 *x, double ucf)
/*
**----------------------------------------------------------------
** Input: j = type of link variable
@@ -413,10 +413,10 @@ int linkoutput(int j, float *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, float *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, float *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, float *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];

312
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*/
X = (double *) calloc(MAX((Nnodes+1),(Nlinks+1)),sizeof(double));
R = (double *) calloc((Nlinks+1), sizeof(double));
ERRCODE(MEMCHECK(X));
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,12 +148,16 @@ 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++)
if (Node[i].S != NULL) Node[i].S->Smass = 0.0;
for (i=1; i<=Nnodes; i++) {
if (Node[i].S != NULL) Node[i].S->Smass = 0.0;
}
QTankVolumes = calloc(Ntanks, sizeof(double)); // keep track of previous step's tank volumes.
QLinkFlow = calloc(Nlinks, sizeof(double)); // keep track of previous step's link flows.
/* Set WQ parameters */
Bucf = 1.0;
Tucf = 1.0;
@@ -161,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)
@@ -189,13 +193,18 @@ void initqual()
Wsource = 0.0;
/* Re-position hydraulics file */
if (!OpenHflag) {
fseek(HydFile,HydOffset,SEEK_SET);
}
/* Set elapsed times to zero */
Htime = 0;
Qtime = 0;
Rtime = Rstart;
Nperiods = 0;
initsegs();
}
@@ -221,8 +230,40 @@ int runqual(long *t)
if (Qtime == Htime)
{
errcode = gethyd(&hydtime, &hydstep);
Htime = hydtime + hydstep;
if (!OpenHflag) { // test for sequential vs stepwise
// sequential
Htime = hydtime + hydstep;
}
else {
// stepwise calculation - hydraulic results are already in memory
for (int i=1; i<= Ntanks; ++i) {
QTankVolumes[i-1] = Tank[i].V;
}
for (int i=1; i<= Nlinks; ++i)
{
if (LinkStatus[i] <= CLOSED) {
QLinkFlow[i-1] = Q[i];
}
}
}
}
else {
// stepwise calculation
for (int i=1; i<= Ntanks; ++i) {
QTankVolumes[i-1] = Tank[i].V;
}
for (int i=1; i<= Nlinks; ++i)
{
if (LinkStatus[i] <= CLOSED) {
QLinkFlow[i-1] = Q[i];
}
}
}
return(errcode);
}
@@ -243,8 +284,41 @@ int nextqual(long *tstep)
/* Determine time step */
*tstep = 0;
hydstep = Htime - Qtime;
// hydstep = Htime - Qtime;
if (Htime <= Dur) hydstep = Htime - Qtime;
else hydstep = 0;
double *tankVolumes;
// if we're operating in stepwise mode, capture the tank levels so we can restore them later.
if (OpenHflag) {
tankVolumes = calloc(Ntanks, sizeof(double));
for (int i=1; i<=Ntanks; ++i) {
if (Tank[i].A != 0) { // skip reservoirs
tankVolumes[i-1] = Tank[i].V;
}
}
// restore the previous step's tank volumes
for (int i=1; i<=Ntanks; i++) {
if (Tank[i].A != 0) { // skip reservoirs again
int n = Tank[i].Node;
Tank[i].V = QTankVolumes[i-1];
NodeHead[n] = tankgrade(i,Tank[i].V);
}
}
// restore the previous step's pipe link flows
for (int i=1; i<=Nlinks; i++) {
if (LinkStatus[i] <= CLOSED) {
Q[i] = 0.0;
}
}
}
/* Perform water quality routing over this time step */
if (Qualflag != NONE && hydstep > 0) transport(hydstep);
@@ -255,6 +329,26 @@ int nextqual(long *tstep)
/* Save final output if no more time steps */
if (!errcode && Saveflag && *tstep == 0) errcode = savefinaloutput();
// restore tank levels to post-runH state, if needed.
if (OpenHflag) {
for (int i=1; i<=Ntanks; i++) {
if (Tank[i].A != 0) { // skip reservoirs again
int n = Tank[i].Node;
Tank[i].V = tankVolumes[i-1];
NodeHead[n] = tankgrade(i,Tank[i].V);
}
}
for (int i=1; i<=Nlinks; ++i) {
if (LinkStatus[i] <= CLOSED) {
Q[i] = QLinkFlow[i-1];
}
}
free(tankVolumes);
}
return(errcode);
}
@@ -320,8 +414,10 @@ int closequal()
free(FlowDir);
free(VolIn);
free(MassIn);
free(R);
free(X);
free(PipeRateCoeff);
free(TempQual);
free(QTankVolumes);
free(QLinkFlow);
return(errcode);
}
@@ -344,10 +440,14 @@ int gethyd(long *hydtime, long *hydstep)
{
int errcode = 0;
/* Read hydraulic results from file */
if (!readhyd(hydtime)) return(307);
if (!readhydstep(hydstep)) return(307);
Htime = *hydtime;
// if hydraulics are not open, then we're operating in sequential mode.
// else hydraulics are open, so use the hydraulic results in memory rather than reading from the temp file.
if (!OpenHflag) {
/* Read hydraulic results from file */
if (!readhyd(hydtime)) return(307);
if (!readhydstep(hydstep)) return(307);
Htime = *hydtime;
}
/* Save current results to output file */
if (Htime >= Rtime)
@@ -365,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);
}
@@ -415,7 +523,7 @@ void transport(long tstep)
*/
{
long qtime, dt;
/* Repeat until elapsed time equals hydraulic time step */
AllocSetPool(SegPool); //(2.00.11 - LR)
@@ -431,6 +539,7 @@ void transport(long tstep)
release(dt); /* Release new nodal flows */
}
updatesourcenodes(tstep); /* Update quality at source nodes */
}
@@ -451,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;
@@ -460,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 */
@@ -518,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 */
@@ -584,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;
}
}
@@ -666,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(X,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) */
@@ -685,9 +800,13 @@ void accumulate(long dt)
VolIn[j]++;
}
}
for (k=1; k<=Nnodes; k++)
if (VolIn[k] > 0.0) X[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));
@@ -705,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);
@@ -767,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. X[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] = X[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;
}
@@ -809,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 X */
memset(X,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;
@@ -824,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 */
@@ -840,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;
@@ -860,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 */
@@ -872,7 +1001,7 @@ void sourceinput(long dt)
}
/* Source concen. contribution = (mass added / outflow volume) */
X[n] = massadded/volout;
TempQual[n] = massadded/volout;
/* Update total mass added for time period & simulation */
source->Smass += massadded;
@@ -888,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];
}
}
}
@@ -923,7 +1052,7 @@ void release(long dt)
v = q*dt;
/* Include source contribution in quality released from node. */
c = C[n] + X[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.*/
@@ -952,7 +1081,7 @@ void updatesourcenodes(long dt)
** Input: dt = current WQ time step
** Output: none
** Purpose: updates quality at source nodes.
** (X[n] = concen. added by source at node n)
** (TempQual[n] = concen. added by source at node n)
**---------------------------------------------------
*/
{
@@ -968,13 +1097,13 @@ void updatesourcenodes(long dt)
if (source == NULL) continue;
/* Add source to current node concen. */
C[n] += X[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 */
@@ -997,21 +1126,22 @@ void updatetanks(long dt)
/* Examine each reservoir & tank */
for (i=1; i<=Ntanks; i++)
{
n = Tank[i].Node;
/* Use initial quality for reservoirs */
if (Tank[i].A == 0.0)
{
n = Tank[i].Node;
C[n] = Node[n].C0;
NodeQual[n] = Node[n].C0;
}
/* Update tank WQ based on mixing model */
else switch(Tank[i].MixModel)
{
case MIX2: tankmix2(i,dt); break;
case FIFO: tankmix3(i,dt); break;
case LIFO: tankmix4(i,dt); break;
default: tankmix1(i,dt); break;
else {
switch(Tank[i].MixModel)
{
case MIX2: tankmix2(i,dt); break;
case FIFO: tankmix3(i,dt); break;
case LIFO: tankmix4(i,dt); break;
default: tankmix1(i,dt); break;
}
}
}
}
@@ -1043,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;
//}
@@ -1068,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. */
@@ -1081,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 ***/
@@ -1118,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;
@@ -1174,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;
}
@@ -1209,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];
@@ -1250,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)
@@ -1300,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;
@@ -1368,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;
}
@@ -1423,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. );
}
@@ -1443,8 +1573,8 @@ void ratecoeffs()
{
kw = Link[k].Kw;
if (kw != 0.0) kw = piperate(k);
Link[k].R = kw;
R[k] = 0.0;
Link[k].Rc = kw;
PipeRateCoeff[k] = 0.0;
}
} /* End of ratecoeffs */
@@ -1526,7 +1656,7 @@ double pipereact(int k, double c, double v, long dt)
/* Otherwise find bulk & wall reaction rates */
rbulk = bulkrate(c,Link[k].Kb,BulkOrder)*Bucf;
rwall = wallrate(c,Link[k].Diam,Link[k].Kw,Link[k].R);
rwall = wallrate(c,Link[k].Diam,Link[k].Kw,Link[k].Rc);
/* Find change in concentration over timestep */
dcbulk = rbulk*(double)dt;

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

332
src/testLemonTiger.cpp Executable file
View File

@@ -0,0 +1,332 @@
#include <map>
#include <iomanip>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "testLemonTiger.h"
#include "toolkit.h"
#define COLW 15
#define OUTPRECISION 6
using namespace std;
typedef struct {
double head;
double demand;
double quality;
} nodeState_t;
typedef struct {
double flow;
} linkState_t;
typedef map<int, nodeState_t> networkNodeState_t; // nodeIndex, state
typedef map<int, linkState_t> networkLinkState_t; // linkIndex, state
typedef struct {
networkNodeState_t nodeState;
networkLinkState_t linkState;
} networkState_t;
typedef map<long, networkState_t> result_t; // time, networkState
// access results by, for instance, resultsContainer[time][nodeIndex].head
void checkErr(int err, std::string function);
void saveHydResults(networkState_t* networkState);
void saveQualResults(networkState_t* networkState);
void printResults(result_t* state1, result_t* state2, std::ostream& out);
void compare(result_t* results1, result_t* results2, std::ostream &out);
int main(int argc, char * argv[]) {
// create storage structures for results.
result_t epanetResults, lemonTigerResults;
cout << "Lemon Tiger TEST" << endl
<< "________________" << endl;
long simulationTime = 0;
long nextEventH = 0, nextEventQ = 0;
long simTimeRemaining = 0;
try {
/* Batch solver (old epanet) */
cout << "*****Original EPANET results******" << endl;
checkErr( ENopen(argv[1], argv[2], (char*)""), "ENopen" );
checkErr( ENopenH(), "ENopenH" );
checkErr( ENinitH(EN_SAVE), "ENinitH" );
cout << "Running hydraulics..." << endl;
do {
/* Solve for hydraulics & advance to next time period */
checkErr( ENrunH(&simulationTime), "ENrunH" );
checkErr( ENnextH(&nextEventH), "ENnextH" );
// gather hydraulic results
saveHydResults(&epanetResults[simulationTime]);
} while (nextEventH > 0);
// hydraulics are done
checkErr( ENcloseH(), "ENcloseH" );
cout << "\t\t\tdone." << endl;
cout << "Running WQ..." << endl;
checkErr( ENopenQ(), "ENopenQ" );
checkErr( ENinitQ(EN_NOSAVE), "ENinitQ" );
do {
checkErr( ENrunQ(&simulationTime), "ENrunQ" );
checkErr( ENnextQ(&nextEventH), "ENstepQ" );
// gather quality results
saveQualResults(&epanetResults[simulationTime]);
} while (nextEventH > 0);
// water quality is done
checkErr( ENcloseQ(), "ENcloseQ" );
cout << "\t\t\tdone." << endl;
// everything is done
checkErr( ENclose(), "ENclose" );
nextEventH = 0;
simTimeRemaining = 0;
simulationTime = 0;
/* stepwise solver (LemonTiger) */
cout << "*****LemonTiger results******" << endl;
checkErr( ENopen(argv[1], argv[2], (char*)""), "ENopen" );
checkErr( ENopenH(), "ENopenH" );
checkErr( ENinitH(EN_NOSAVE), "ENinitH" );
checkErr( ENopenQ(), "ENopenQ" );
checkErr( ENinitQ(EN_NOSAVE), "ENinitQ" );
cout << "Running stepwise hydraulics and water quality..." << endl;
do {
/* Solve for hydraulics & advance to next time period */
checkErr( ENrunH(&simulationTime), "ENrunH" );
checkErr( ENrunQ(&simulationTime), "ENrunQ" );
checkErr( ENnextH(&nextEventH), "ENnextH" );
checkErr( ENnextQ(&nextEventQ), "ENstepQ" );
saveHydResults(&lemonTigerResults[simulationTime]);
saveQualResults(&lemonTigerResults[simulationTime]);
} while (nextEventH > 0);
cout << "\t\t\tdone." << endl;
// all done
checkErr( ENcloseH(), "ENcloseH" );
checkErr( ENcloseQ(), "ENcloseQ" );
checkErr( ENclose(), "ENclose" );
// summarize the results
printResults(&epanetResults, &lemonTigerResults, cout);
compare(&epanetResults, &lemonTigerResults, cout);
} catch (int err) {
cerr << "exiting with error " << err << endl;
}
}
void saveHydResults(networkState_t* networkState) {
int nNodes, nLinks;
float head, demand, flow;
ENgetcount(EN_NODECOUNT, &nNodes);
ENgetcount(EN_LINKCOUNT, &nLinks);
for (int iNode = 1; iNode <= nNodes; ++iNode) {
ENgetnodevalue(iNode, EN_HEAD, &head);
ENgetnodevalue(iNode, EN_DEMAND, &demand);
(*networkState).nodeState[iNode].head = head;
(*networkState).nodeState[iNode].demand = demand;
}
for (int iLink = 1; iLink <= nLinks; ++iLink) {
ENgetlinkvalue(iLink, EN_FLOW, &flow);
(*networkState).linkState[iLink].flow = flow;
}
}
void saveQualResults(networkState_t* networkState) {
int nNodes;
float quality;
ENgetcount(EN_NODECOUNT, &nNodes);
for (int iNode = 1; iNode <= nNodes; iNode++) {
ENgetnodevalue(iNode, EN_QUALITY, &quality);
(*networkState).nodeState[iNode].quality = quality;
}
}
void printResults(result_t* results1, result_t* results2, std::ostream &out) {
result_t::const_iterator resultIterator;
for (resultIterator = (*results1).begin(); resultIterator != (*results1).end(); ++resultIterator) {
// get the current frame
const long time = resultIterator->first;
const networkNodeState_t networkNodeState1= resultIterator->second.nodeState;
//nodeState1 = resultIterator->second.nodeState;
const networkLinkState_t networkLinkState1 = resultIterator->second.linkState;
//linkState1 = resultIterator->second.linkState;
// see if this time is indexed in the second state container
if ((*results2).find(time) == (*results2).end()) {
// nope.
out << "time " << time << " not found in second result set" << endl;
}
else {
// get the second result set's state
const networkNodeState_t networkNodeState2 = (*results2)[time].nodeState;
const networkLinkState_t networkLinkState2 = (*results2)[time].linkState;
// print the current simulation time
out << left;
out << setfill('*') << setw(100) << "*" << endl;
out << setfill(' ');
out << setw(4) << "T = " << setw(6) << time;
out << "|" << setw(3*COLW) << "EPANET";
out << "|" << setw(3*COLW) << "LemonTiger" << endl;
out << setw(10) << "Index" << "|";
out << setw(COLW) << "Demand" << setw(COLW) << "Head" << setw(COLW) << "Quality" << "|";
out << setw(COLW) << "Demand" << setw(COLW) << "Head" << setw(COLW) << "Quality" << endl;
out << setprecision(OUTPRECISION);
// loop through the nodes in the networkState objs, and print out the results for this time period
networkNodeState_t::const_iterator networkNodeIterator;
for (networkNodeIterator = networkNodeState1.begin(); networkNodeIterator != networkNodeState1.end(); ++networkNodeIterator) {
int nodeIndex = networkNodeIterator->first;
// trusting that all nodes are present...
const nodeState_t nodeState1 = networkNodeIterator->second;
const nodeState_t nodeState2 = networkNodeState2.at(nodeIndex);
if (nodeState1.quality != nodeState2.quality ) {
// epanet
out << setw(10) << nodeIndex << "|";
out << setw(COLW) << nodeState1.demand;
out << setw(COLW) << nodeState1.head;
out << setw(COLW) << nodeState1.quality;
// lemontiger
out << "|";
out << setw(COLW) << nodeState2.demand;
out << setw(COLW) << nodeState2.head;
out << setw(COLW) << nodeState2.quality;
out << endl;
}
}
networkLinkState_t::const_iterator networkLinkIterator;
for (networkLinkIterator = networkLinkState1.begin(); networkLinkIterator != networkLinkState1.end(); ++networkLinkIterator) {
int linkIndex = networkLinkIterator->first;
// trusting that all nodes are present...
const linkState_t linkState1 = networkLinkIterator->second;
const linkState_t linkState2 = networkLinkState2.at(linkIndex);
if ( linkState1.flow != linkState2.flow ) {
// epanet
out << setw(10) << linkIndex << "|";
out << setw(COLW) << linkState1.flow;
// lemontiger
out << "|";
out << setw(COLW) << linkState2.flow;
out << endl;
}
}
}
}
}
void compare(result_t* results1, result_t* results2, std::ostream &out) {
double sumHeadDiff=0, sumDemandDiff=0, sumQualDiff=0, sumFlowDiff=0;
result_t::const_iterator resultIterator;
for (resultIterator = (*results1).begin(); resultIterator != (*results1).end(); ++resultIterator) {
// get the current frame
const long time = resultIterator->first;
const networkNodeState_t nodeState1 = resultIterator->second.nodeState;
const networkLinkState_t linkState1 = resultIterator->second.linkState;
// see if this time is indexed in the second state container
if ((*results2).find(time) == (*results2).end()) {
// nope.
out << "time " << time << " not found in second result set" << endl;
}
else {
// get the second result set's state
const networkNodeState_t networkNodeState2 = (*results2)[time].nodeState;
const networkLinkState_t networkLinkState2 = (*results2)[time].linkState;
double qualD=0;
networkNodeState_t::const_iterator networkNodeIterator;
for (networkNodeIterator = nodeState1.begin(); networkNodeIterator != nodeState1.end(); ++networkNodeIterator) {
int nodeIndex = networkNodeIterator->first;
// trusting that all nodes are present...
const nodeState_t nodeState1 = networkNodeIterator->second;
const nodeState_t nodeState2 = networkNodeState2.at(nodeIndex);
sumHeadDiff += fabs(nodeState1.head - nodeState2.head);
sumDemandDiff += fabs(nodeState1.demand - nodeState2.demand);
qualD += fabs(nodeState1.quality - nodeState2.quality);
}
//out << "T: " << time << " dq: " << setprecision(20) << qualD << endl;
sumQualDiff += qualD;
networkLinkState_t::const_iterator networkLinkIterator;
for (networkLinkIterator = linkState1.begin(); networkLinkIterator != linkState1.end(); ++networkLinkIterator) {
int linkIndex = networkLinkIterator->first;
// trusting that all nodes are present...
const linkState_t linkState1 = networkLinkIterator->second;
const linkState_t linkState2 = networkLinkState2.at(linkIndex);
sumFlowDiff += fabs(linkState1.flow - linkState2.flow);
}
}
}
int c1 = 18;
int p = 20;
out << setw(c1) << "Head Diff:" << setprecision(p) << sumHeadDiff << endl;
out << setw(c1) << "Demand Diff:" << setprecision(p) << sumDemandDiff << endl;
out << setw(c1) << "Quality Diff:" << setprecision(p) << sumQualDiff << endl;
out << setw(c1) << "Flow Diff:" << setprecision(p) << sumFlowDiff << endl;
}
void checkErr(int err, std::string function) {
if (err > 0) {
cerr << "Error in " << function << ": " << err << endl;
char errmsg[1024];
ENgeterror(err, errmsg, 1024);
cerr << errmsg << endl;
throw err;
}
}

15
src/testLemonTiger.h Executable file
View File

@@ -0,0 +1,15 @@
//
// testLemonTiger.h
// epanet
//
// Created by Sam Hatchett on 2/1/13.
//
//
#ifndef __epanet__testLemonTiger__
#define __epanet__testLemonTiger__
#include <iostream>
#include <vector>
#endif /* defined(__epanet__testLemonTiger__) */

5
src/text.h
View File

@@ -14,6 +14,8 @@ AUTHOR: L. Rossman
****************************************************
*/
/* ------------ Keyword Dictionary ---------- */
#ifndef TEXT_H
#define TEXT_H
#define w_USE "USE"
#define w_SAVE "SAVE"
@@ -501,6 +503,8 @@ AUTHOR: L. Rossman
#define ERR308 "File Error 308: cannot save results to file."
#define ERR309 "File Error 309: cannot save results to report file."
#define ERR401 "Sync Error 401: Qstep is not dividable by Hstep. Can't sync."
#define R_ERR201 "Input Error 201: syntax error in following line of "
#define R_ERR202 "Input Error 202: illegal numeric value in following line of "
#define R_ERR203 "Input Error 203: undefined node in following line of "
@@ -528,3 +532,4 @@ AUTHOR: L. Rossman
#define WARN5 "WARNING: Valves cannot deliver enough flow."
#define WARN6 "WARNING: System has negative pressures."
#endif

508
src/toolkit.h
View File

@@ -1,255 +1,255 @@
/*
*******************************************************************
TOOLKIT.H - Prototypes for EPANET Functions Exported to DLL Toolkit
VERSION: 2.00
DATE: 5/8/00
10/25/00
3/1/01
8/15/07 (2.00.11)
2/14/08 (2.00.12)
AUTHOR: L. Rossman
US EPA - NRMRL
*******************************************************************
*/
#ifndef DLLEXPORT
#ifdef DLL
#ifdef __cplusplus
#define DLLEXPORT extern "C" __declspec(dllexport)
#else
/*
*******************************************************************
TOOLKIT.H - Prototypes for EPANET Functions Exported to DLL Toolkit
VERSION: 2.00
DATE: 5/8/00
10/25/00
3/1/01
8/15/07 (2.00.11)
2/14/08 (2.00.12)
AUTHOR: L. Rossman
US EPA - NRMRL
*******************************************************************
*/
#ifndef DLLEXPORT
#ifdef DLL
#ifdef __cplusplus
#define DLLEXPORT extern "C" __declspec(dllexport)
#else
#define DLLEXPORT __declspec(dllexport) __stdcall
#endif
#elif defined(CYGWIN)
#define DLLEXPORT __stdcall
#else
#ifdef __cplusplus
#define DLLEXPORT
#else
#define DLLEXPORT
#endif
#endif
#endif
// --- Define the EPANET toolkit constants
#define EN_ELEVATION 0 /* Node parameters */
#define EN_BASEDEMAND 1
#define EN_PATTERN 2
#define EN_EMITTER 3
#define EN_INITQUAL 4
#define EN_SOURCEQUAL 5
#define EN_SOURCEPAT 6
#define EN_SOURCETYPE 7
#define EN_TANKLEVEL 8
#define EN_DEMAND 9
#define EN_HEAD 10
#define EN_PRESSURE 11
#define EN_QUALITY 12
#define EN_SOURCEMASS 13
#define EN_INITVOLUME 14
#define EN_MIXMODEL 15
#define EN_MIXZONEVOL 16
#define EN_TANKDIAM 17
#define EN_MINVOLUME 18
#define EN_VOLCURVE 19
#define EN_MINLEVEL 20
#define EN_MAXLEVEL 21
#define EN_MIXFRACTION 22
#define EN_TANK_KBULK 23
#define EN_TANKVOLUME 24
#define EN_MAXVOLUME 25
#define EN_DIAMETER 0 /* Link parameters */
#define EN_LENGTH 1
#define EN_ROUGHNESS 2
#define EN_MINORLOSS 3
#define EN_INITSTATUS 4
#define EN_INITSETTING 5
#define EN_KBULK 6
#define EN_KWALL 7
#define EN_FLOW 8
#define EN_VELOCITY 9
#define EN_HEADLOSS 10
#define EN_STATUS 11
#define EN_SETTING 12
#define EN_ENERGY 13
#define EN_LINKQUAL 14 /* TNT */
#define EN_LINKPATTERN 15
#define EN_DURATION 0 /* Time parameters */
#define EN_HYDSTEP 1
#define EN_QUALSTEP 2
#define EN_PATTERNSTEP 3
#define EN_PATTERNSTART 4
#define EN_REPORTSTEP 5
#define EN_REPORTSTART 6
#define EN_RULESTEP 7
#define EN_STATISTIC 8
#define EN_PERIODS 9
#define EN_STARTTIME 10 /* Added TNT 10/2/2009 */
#define EN_HTIME 11
#define EN_HALTFLAG 12
#define EN_NEXTEVENT 13
#define EN_ITERATIONS 0
#define EN_RELATIVEERROR 1
#define EN_NODECOUNT 0 /* Component counts */
#define EN_TANKCOUNT 1
#define EN_LINKCOUNT 2
#define EN_PATCOUNT 3
#define EN_CURVECOUNT 4
#define EN_CONTROLCOUNT 5
#define EN_JUNCTION 0 /* Node types */
#define EN_RESERVOIR 1
#define EN_TANK 2
#define EN_CVPIPE 0 /* Link types. */
#define EN_PIPE 1 /* See LinkType in TYPES.H */
#define EN_PUMP 2
#define EN_PRV 3
#define EN_PSV 4
#define EN_PBV 5
#define EN_FCV 6
#define EN_TCV 7
#define EN_GPV 8
#define EN_NONE 0 /* Quality analysis types. */
#define EN_CHEM 1 /* See QualType in TYPES.H */
#define EN_AGE 2
#define EN_TRACE 3
#define EN_CONCEN 0 /* Source quality types. */
#define EN_MASS 1 /* See SourceType in TYPES.H. */
#define EN_SETPOINT 2
#define EN_FLOWPACED 3
#define EN_CFS 0 /* Flow units types. */
#define EN_GPM 1 /* See FlowUnitsType */
#define EN_MGD 2 /* in TYPES.H. */
#define EN_IMGD 3
#define EN_AFD 4
#define EN_LPS 5
#define EN_LPM 6
#define EN_MLD 7
#define EN_CMH 8
#define EN_CMD 9
#define EN_TRIALS 0 /* Misc. options */
#define EN_ACCURACY 1
#define EN_TOLERANCE 2
#define EN_EMITEXPON 3
#define EN_DEMANDMULT 4
#define EN_LOWLEVEL 0 /* Control types. */
#define EN_HILEVEL 1 /* See ControlType */
#define EN_TIMER 2 /* in TYPES.H. */
#define EN_TIMEOFDAY 3
#define EN_AVERAGE 1 /* Time statistic types. */
#define EN_MINIMUM 2 /* See TstatType in TYPES.H */
#define EN_MAXIMUM 3
#define EN_RANGE 4
#define EN_MIX1 0 /* Tank mixing models */
#define EN_MIX2 1
#define EN_FIFO 2
#define EN_LIFO 3
#define EN_NOSAVE 0 /* Save-results-to-file flag */
#define EN_SAVE 1
#define EN_INITFLOW 10 /* Re-initialize flows flag */
#define EN_CONST_HP 0 /* constant horsepower */
#define EN_POWER_FUNC 1 /* power function */
#define EN_CUSTOM 2 /* user-defined custom curve */
// --- Declare the EPANET toolkit functions
#if defined(__cplusplus)
extern "C" {
#endif
int DLLEXPORT ENepanet(char *, char *, char *, void (*) (char *));
int DLLEXPORT ENopen(char *, char *, char *);
int DLLEXPORT ENsaveinpfile(char *);
int DLLEXPORT ENclose(void);
int DLLEXPORT ENsolveH(void);
int DLLEXPORT ENsaveH(void);
int DLLEXPORT ENopenH(void);
int DLLEXPORT ENinitH(int);
int DLLEXPORT ENrunH(long *);
int DLLEXPORT ENnextH(long *);
int DLLEXPORT ENcloseH(void);
int DLLEXPORT ENsavehydfile(char *);
int DLLEXPORT ENusehydfile(char *);
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 ENwriteline(char *);
int DLLEXPORT ENreport(void);
int DLLEXPORT ENresetreport(void);
int DLLEXPORT ENsetreport(char *);
int DLLEXPORT ENgetcontrol(int, int *, int *, float *,
int *, float *);
int DLLEXPORT ENgetcount(int, int *);
int DLLEXPORT ENgetoption(int, float *);
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, float *);
int DLLEXPORT ENgetqualtype(int *, int *);
int DLLEXPORT ENgeterror(int, char *, int);
int DLLEXPORT ENgetstatistic(int code, int* value);
int DLLEXPORT ENgetnodeindex(char *, int *);
int DLLEXPORT ENgetnodeid(int, char *);
int DLLEXPORT ENgetnodetype(int, int *);
int DLLEXPORT ENgetnodevalue(int, int, float *);
int DLLEXPORT ENgetnumdemands(int, int *);
int DLLEXPORT ENgetbasedemand(int, int, float *);
int DLLEXPORT ENgetdemandpattern(int, int, int *);
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, float *);
int DLLEXPORT ENgetcurve(int curveIndex, int *nValues, float **xValues, float **yValues);
int DLLEXPORT ENgetheadcurve(int, char *);
int DLLEXPORT ENgetpumptype(int, int *);
int DLLEXPORT ENgetversion(int *);
int DLLEXPORT ENsetcontrol(int, int, int, float, int, float);
int DLLEXPORT ENsetnodevalue(int, int, float);
int DLLEXPORT ENsetlinkvalue(int, int, float);
int DLLEXPORT ENaddpattern(char *);
int DLLEXPORT ENsetpattern(int, float *, int);
int DLLEXPORT ENsetpatternvalue(int, int, float);
int DLLEXPORT ENsettimeparam(int, long);
int DLLEXPORT ENsetoption(int, float);
int DLLEXPORT ENsetstatusreport(int);
int DLLEXPORT ENsetqualtype(int, char *, char *, char *);
#if defined(__cplusplus)
}
#endif
#endif
#elif defined(CYGWIN)
#define DLLEXPORT __stdcall
#else
#ifdef __cplusplus
#define DLLEXPORT
#else
#define DLLEXPORT
#endif
#endif
#endif
// --- Define the EPANET toolkit constants
#define EN_ELEVATION 0 /* Node parameters */
#define EN_BASEDEMAND 1
#define EN_PATTERN 2
#define EN_EMITTER 3
#define EN_INITQUAL 4
#define EN_SOURCEQUAL 5
#define EN_SOURCEPAT 6
#define EN_SOURCETYPE 7
#define EN_TANKLEVEL 8
#define EN_DEMAND 9
#define EN_HEAD 10
#define EN_PRESSURE 11
#define EN_QUALITY 12
#define EN_SOURCEMASS 13
#define EN_INITVOLUME 14
#define EN_MIXMODEL 15
#define EN_MIXZONEVOL 16
#define EN_TANKDIAM 17
#define EN_MINVOLUME 18
#define EN_VOLCURVE 19
#define EN_MINLEVEL 20
#define EN_MAXLEVEL 21
#define EN_MIXFRACTION 22
#define EN_TANK_KBULK 23
#define EN_TANKVOLUME 24
#define EN_MAXVOLUME 25
#define EN_DIAMETER 0 /* Link parameters */
#define EN_LENGTH 1
#define EN_ROUGHNESS 2
#define EN_MINORLOSS 3
#define EN_INITSTATUS 4
#define EN_INITSETTING 5
#define EN_KBULK 6
#define EN_KWALL 7
#define EN_FLOW 8
#define EN_VELOCITY 9
#define EN_HEADLOSS 10
#define EN_STATUS 11
#define EN_SETTING 12
#define EN_ENERGY 13
#define EN_LINKQUAL 14 /* TNT */
#define EN_LINKPATTERN 15
#define EN_DURATION 0 /* Time parameters */
#define EN_HYDSTEP 1
#define EN_QUALSTEP 2
#define EN_PATTERNSTEP 3
#define EN_PATTERNSTART 4
#define EN_REPORTSTEP 5
#define EN_REPORTSTART 6
#define EN_RULESTEP 7
#define EN_STATISTIC 8
#define EN_PERIODS 9
#define EN_STARTTIME 10 /* Added TNT 10/2/2009 */
#define EN_HTIME 11
#define EN_HALTFLAG 12
#define EN_NEXTEVENT 13
#define EN_ITERATIONS 0
#define EN_RELATIVEERROR 1
#define EN_NODECOUNT 0 /* Component counts */
#define EN_TANKCOUNT 1
#define EN_LINKCOUNT 2
#define EN_PATCOUNT 3
#define EN_CURVECOUNT 4
#define EN_CONTROLCOUNT 5
#define EN_JUNCTION 0 /* Node types */
#define EN_RESERVOIR 1
#define EN_TANK 2
#define EN_CVPIPE 0 /* Link types. */
#define EN_PIPE 1 /* See LinkType in TYPES.H */
#define EN_PUMP 2
#define EN_PRV 3
#define EN_PSV 4
#define EN_PBV 5
#define EN_FCV 6
#define EN_TCV 7
#define EN_GPV 8
#define EN_NONE 0 /* Quality analysis types. */
#define EN_CHEM 1 /* See QualType in TYPES.H */
#define EN_AGE 2
#define EN_TRACE 3
#define EN_CONCEN 0 /* Source quality types. */
#define EN_MASS 1 /* See SourceType in TYPES.H. */
#define EN_SETPOINT 2
#define EN_FLOWPACED 3
#define EN_CFS 0 /* Flow units types. */
#define EN_GPM 1 /* See FlowUnitsType */
#define EN_MGD 2 /* in TYPES.H. */
#define EN_IMGD 3
#define EN_AFD 4
#define EN_LPS 5
#define EN_LPM 6
#define EN_MLD 7
#define EN_CMH 8
#define EN_CMD 9
#define EN_TRIALS 0 /* Misc. options */
#define EN_ACCURACY 1
#define EN_TOLERANCE 2
#define EN_EMITEXPON 3
#define EN_DEMANDMULT 4
#define EN_LOWLEVEL 0 /* Control types. */
#define EN_HILEVEL 1 /* See ControlType */
#define EN_TIMER 2 /* in TYPES.H. */
#define EN_TIMEOFDAY 3
#define EN_AVERAGE 1 /* Time statistic types. */
#define EN_MINIMUM 2 /* See TstatType in TYPES.H */
#define EN_MAXIMUM 3
#define EN_RANGE 4
#define EN_MIX1 0 /* Tank mixing models */
#define EN_MIX2 1
#define EN_FIFO 2
#define EN_LIFO 3
#define EN_NOSAVE 0 /* Save-results-to-file flag */
#define EN_SAVE 1
#define EN_INITFLOW 10 /* Re-initialize flows flag */
#define EN_CONST_HP 0 /* constant horsepower */
#define EN_POWER_FUNC 1 /* power function */
#define EN_CUSTOM 2 /* user-defined custom curve */
// --- Declare the EPANET toolkit functions
#if defined(__cplusplus)
extern "C" {
#endif
int DLLEXPORT ENepanet(char *, char *, char *, void (*) (char *));
int DLLEXPORT ENopen(char *, char *, char *);
int DLLEXPORT ENsaveinpfile(char *);
int DLLEXPORT ENclose(void);
int DLLEXPORT ENsolveH(void);
int DLLEXPORT ENsaveH(void);
int DLLEXPORT ENopenH(void);
int DLLEXPORT ENinitH(int);
int DLLEXPORT ENrunH(long *);
int DLLEXPORT ENnextH(long *);
int DLLEXPORT ENcloseH(void);
int DLLEXPORT ENsavehydfile(char *);
int DLLEXPORT ENusehydfile(char *);
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 ENwriteline(char *);
int DLLEXPORT ENreport(void);
int DLLEXPORT ENresetreport(void);
int DLLEXPORT ENsetreport(char *);
int DLLEXPORT ENgetcontrol(int, int *, int *, float *,
int *, float *);
int DLLEXPORT ENgetcount(int, int *);
int DLLEXPORT ENgetoption(int, float *);
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, float *);
int DLLEXPORT ENgetqualtype(int *, int *);
int DLLEXPORT ENgeterror(int, char *, int);
int DLLEXPORT ENgetstatistic(int code, int* value);
int DLLEXPORT ENgetnodeindex(char *, int *);
int DLLEXPORT ENgetnodeid(int, char *);
int DLLEXPORT ENgetnodetype(int, int *);
int DLLEXPORT ENgetnodevalue(int, int, float *);
int DLLEXPORT ENgetnumdemands(int, int *);
int DLLEXPORT ENgetbasedemand(int, int, float *);
int DLLEXPORT ENgetdemandpattern(int, int, int *);
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, float *);
int DLLEXPORT ENgetcurve(int curveIndex, int *nValues, float **xValues, float **yValues);
int DLLEXPORT ENgetheadcurve(int, char *);
int DLLEXPORT ENgetpumptype(int, int *);
int DLLEXPORT ENgetversion(int *);
int DLLEXPORT ENsetcontrol(int, int, int, float, int, float);
int DLLEXPORT ENsetnodevalue(int, int, float);
int DLLEXPORT ENsetlinkvalue(int, int, float);
int DLLEXPORT ENaddpattern(char *);
int DLLEXPORT ENsetpattern(int, float *, int);
int DLLEXPORT ENsetpatternvalue(int, int, float);
int DLLEXPORT ENsettimeparam(int, long);
int DLLEXPORT ENsetoption(int, float);
int DLLEXPORT ENsetstatusreport(int);
int DLLEXPORT ENsetqualtype(int, char *, char *, char *);
#if defined(__cplusplus)
}
#endif

15
src/types.h
View File

@@ -17,6 +17,8 @@ AUTHOR: L. Rossman
**********************************************************************
*/
#ifndef TYPES_H
#define TYPES_H
/*********************************************************/
/* All floats have been re-declared as doubles (7/3/07). */
@@ -26,7 +28,7 @@ AUTHOR: L. Rossman
Definition of 4-byte integers & reals
-------------------------------------------
*/
typedef float REAL4; //(2.00.11 - LR)
typedef double REAL4; //(2.00.11 - LR)
typedef int INT4; //(2.00.12 - LR)
/*
@@ -35,7 +37,7 @@ typedef int INT4; /
-----------------------------
*/
/*** Updated ***/
#define CODEVERSION 20012 //(2.00.12 - LR)
#define CODEVERSION 20100
#define MAGICNUMBER 516114521
#define VERSION 200
#define EOFMARK 0x1A /* Use 0x04 for UNIX systems */
@@ -165,6 +167,13 @@ typedef struct /* CURVE OBJECT */
double *Y; /* Y-values */
} Scurve;
typedef struct /* Coord OBJECT */
{
char ID[MAXID+1]; /* Coord ID */
double *X; /* X-values */
double *Y; /* Y-values */
} Scoord;
struct Sdemand /* DEMAND CATEGORY OBJECT */
{
double Base; /* Baseline demand */
@@ -206,6 +215,7 @@ typedef struct /* LINK OBJECT */
double Kb; /* Bulk react. coeff */
double Kw; /* Wall react. coeff */
double R; /* Flow resistance */
double Rc; /* Reaction cal */
char Type; /* Link type */
char Stat; /* Initial status */
char Rpt; /* Reporting flag */
@@ -451,3 +461,4 @@ enum HdrType /* Type of table heading */
NODEHDR, /* Node Results */
LINKHDR}; /* Link Results */
#endif

56
src/vars.h
View File

@@ -11,15 +11,21 @@ AUTHOR: L. Rossman
************************************************************************
*/
EXTERN FILE *InFile, /* Input file pointer */
#ifndef VARS_H
#define VARS_H
#include <stdio.h>
#include "hash.h"
FILE *InFile, /* Input file pointer */
*OutFile, /* Output file pointer */
*RptFile, /* Report file pointer */
*HydFile, /* Hydraulics file pointer */
*TmpOutFile; /* Temporary file handle */
EXTERN long HydOffset, /* Hydraulics file byte offset */
long HydOffset, /* Hydraulics file byte offset */
OutOffset1, /* 1st output file byte offset */
OutOffset2; /* 2nd output file byte offset */
EXTERN char Msg[MAXMSG+1], /* Text of output message */
char Msg[MAXMSG+1], /* Text of output message */
InpFname[MAXFNAME+1], /* Input file name */
Rpt1Fname[MAXFNAME+1], /* Primary report file name */
Rpt2Fname[MAXFNAME+1], /* Secondary report file name */
@@ -59,7 +65,7 @@ EXTERN char Msg[MAXMSG+1], /* Text of output message */
OpenQflag, /* Quality system opened flag */
SaveQflag, /* Quality results saved flag */
Saveflag; /* General purpose save flag */
EXTERN int MaxNodes, /* Node count from input file */
int MaxNodes, /* Node count from input file */
MaxLinks, /* Link count from input file */
MaxJuncs, /* Junction count */
MaxPipes, /* Pipe count */
@@ -70,6 +76,7 @@ EXTERN int MaxNodes, /* Node count from input file */
MaxRules, /* Rule count */
MaxPats, /* Pattern count */
MaxCurves, /* Curve count */
MaxCoords, /* Coords count */
Nnodes, /* Number of network nodes */
Ntanks, /* Number of tanks */
Njuncs, /* Number of junction nodes */
@@ -81,6 +88,7 @@ EXTERN int MaxNodes, /* Node count from input file */
Nrules, /* Number of control rules */
Npats, /* Number of time patterns */
Ncurves, /* Number of data curves */
Ncoords, /* Number of Coords */
Nperiods, /* Number of reporting periods */
Ncoeffs, /* Number of non-0 matrix coeffs*/
DefPat, /* Default demand pattern */
@@ -91,7 +99,7 @@ EXTERN int MaxNodes, /* Node count from input file */
PageSize, /* Lines/page in output report */
CheckFreq, /* Hydraulics solver parameter */
MaxCheck; /* Hydraulics solver parameter */
EXTERN double Ucf[MAXVAR], /* Unit conversion factors */
double Ucf[MAXVAR], /* Unit conversion factors */
Ctol, /* Water quality tolerance */
Htol, /* Hydraulic head tolerance */
Qtol, /* Flow rate tolerance */
@@ -120,7 +128,7 @@ EXTERN double Ucf[MAXVAR], /* Unit conversion factors */
Wwall, /* Avg. wall reaction rate */
Wtank, /* Avg. tank reaction rate */
Wsource; /* Avg. mass inflow */
EXTERN long Tstart, /* Starting time of day (sec) */
long Tstart, /* Starting time of day (sec) */
Hstep, /* Nominal hyd. time step (sec) */
Qstep, /* Quality time step (sec) */
Pstep, /* Time pattern time step (sec) */
@@ -133,32 +141,38 @@ EXTERN long Tstart, /* Starting time of day (sec) */
Hydstep, /* Actual hydraulic time step */
Rulestep, /* Rule evaluation time step */
Dur; /* Duration of simulation (sec) */
EXTERN SField Field[MAXVAR]; /* Output reporting fields */
SField Field[MAXVAR]; /* Output reporting fields */
/* Array pointers not allocated and freed in same routine */
EXTERN char *S, /* Link status */
char *LinkStatus, /* Link status */
*OldStat; /* Previous link/tank status */
EXTERN 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 */
*X; /* General purpose array */
EXTERN double *H; /* Node heads */
*PipeRateCoeff, /* Pipe reaction rate */
*X, /* General purpose array */
*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 */
EXTERN STmplist *Curvelist; /* Temporary list of curves */
EXTERN STmplist *Coordlist; /* Temporary list of coordinates*/
EXTERN Spattern *Pattern; /* Time patterns */
EXTERN Scurve *Curve; /* Curve data */
EXTERN Scoord *Coord; /* Coordinate data */
EXTERN Snode *Node; /* Node data */
EXTERN Slink *Link; /* Link data */
EXTERN Stank *Tank; /* Tank data */
EXTERN Spump *Pump; /* Pump data */
EXTERN Svalve *Valve; /* Valve data */
EXTERN Scontrol *Control; /* Control data */
EXTERN HTtable *Nht, *Lht; /* Hash tables for ID labels */
EXTERN ENHashTable *NodeHashTable, *LinkHashTable; /* 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
@@ -180,18 +194,20 @@ EXTERN int _relativeError, _iterations; /* Info about hydraulic solution */
** The following arrays are used to efficiently manage this sparsity:
*/
EXTERN double *Aii, /* Diagonal coeffs. of A */
double *Aii, /* Diagonal coeffs. of A */
*Aij, /* Non-zero, off-diagonal coeffs. of A */
*F; /* Right hand side coeffs. */
EXTERN double *P, /* Inverse headloss derivatives */
double *P, /* Inverse headloss derivatives */
*Y; /* Flow correction factors */
EXTERN int *Order, /* Node-to-row of A */
int *Order, /* Node-to-row of A */
*Row, /* Row-to-node of A */
*Ndx; /* Index of link's coeff. in Aij */
/*
** The following arrays store the positions of the non-zero coeffs.
** of the lower triangular portion of A whose values are stored in Aij:
*/
EXTERN int *XLNZ, /* Start position of each column in NZSUB */
int *XLNZ, /* Start position of each column in NZSUB */
*NZSUB, /* Row index of each coeff. in each column */
*LNZ; /* Position of each coeff. in Aij array */
#endif