diff --git a/include/epanet2.bas b/include/epanet2.bas index 789a441..0826a3b 100644 --- a/include/epanet2.bas +++ b/include/epanet2.bas @@ -128,6 +128,8 @@ Public Const EN_ACCURACY = 1 Public Const EN_TOLERANCE = 2 Public Const EN_EMITEXPON = 3 Public Const EN_DEMANDMULT = 4 +Public Const EN_HEADERROR = 5 +Public Const EN_FLOWCHANGE = 6 Public Const EN_LOWLEVEL = 0 ' Control types Public Const EN_HILEVEL = 1 diff --git a/include/epanet2.h b/include/epanet2.h index ebd367b..81b469e 100644 --- a/include/epanet2.h +++ b/include/epanet2.h @@ -212,7 +212,9 @@ typedef enum { EN_ACCURACY = 1, EN_TOLERANCE = 2, EN_EMITEXPON = 3, - EN_DEMANDMULT = 4 + EN_DEMANDMULT = 4, + EN_HEADERROR = 5, + EN_FLOWCHANGE = 6 } EN_Option; typedef enum { diff --git a/include/epanet2.vb b/include/epanet2.vb index 44a47a1..fcb71af 100644 --- a/include/epanet2.vb +++ b/include/epanet2.vb @@ -121,6 +121,8 @@ Public Const EN_ACCURACY = 1 Public Const EN_TOLERANCE = 2 Public Const EN_EMITEXPON = 3 Public Const EN_DEMANDMULT = 4 +Public Const EN_HEADERROR = 5 +Public Const EN_FLOWCHANGE = 6 Public Const EN_LOWLEVEL = 0 ' Control types Public Const EN_HILEVEL = 1 diff --git a/src/epanet.c b/src/epanet.c index 602e0ed..f0f0ff8 100644 --- a/src/epanet.c +++ b/src/epanet.c @@ -71,11 +71,12 @@ This module calls the following functions that reside in other modules: runhyd() nexthyd() closehyd() - resistance() tankvolume() getenergy() setlinkstatus() setlinksetting() + HYDCOEFFS + resistcoeff() QUALITY.C openqual() initqual() @@ -129,6 +130,7 @@ execute function x and set the error code equal to its return value. #include "text.h" #include "types.h" #define EXTERN +////////////////////////////////////////////#include "epanet2.h" #include "vars.h" /**************************************************************** @@ -1341,6 +1343,14 @@ int DLLEXPORT EN_getoption(EN_Project *pr, EN_Option code, case EN_DEMANDMULT: v = hyd->Dmult; break; + + case EN_HEADERROR: + v = hyd->HeadErrorLimit * Ucf[HEAD]; + break; + case EN_FLOWCHANGE: + v = hyd->FlowChangeLimit * Ucf[FLOW]; + break; + default: return (251); } @@ -2648,7 +2658,7 @@ int DLLEXPORT EN_setlinkvalue(EN_Project *p, int index, int code, r = Link[index].Diam / value; /* Ratio of old to new diam */ Link[index].Km *= SQR(r) * SQR(r); /* Adjust minor loss factor */ Link[index].Diam = value; /* Update diameter */ - resistance(p, index); /* Update resistance factor */ + resistcoeff(p, index); /* Update resistance coeff. */ } break; @@ -2657,7 +2667,7 @@ int DLLEXPORT EN_setlinkvalue(EN_Project *p, int index, int code, if (value <= 0.0) return (202); Link[index].Len = value / Ucf[ELEV]; - resistance(p, index); + resistcoeff(p, index); } break; @@ -2668,7 +2678,7 @@ int DLLEXPORT EN_setlinkvalue(EN_Project *p, int index, int code, Link[index].Kc = value; if (hyd->Formflag == DW) Link[index].Kc /= (1000.0 * Ucf[ELEV]); - resistance(p, index); + resistcoeff(p, index); } break; @@ -3171,6 +3181,18 @@ int DLLEXPORT EN_setoption(EN_Project *p, int code, EN_API_FLOAT_TYPE v) return (202); hyd->Dmult = value; break; + + case EN_HEADERROR: + if (value < 0.0) + return (202); + hyd->HeadErrorLimit = value / Ucf[HEAD]; + break; + case EN_FLOWCHANGE: + if (value < 0.0) + return (202); + hyd->FlowChangeLimit = value / Ucf[FLOW]; + break; + default: return (251); } @@ -3419,7 +3441,7 @@ int openhydfile(EN_Project *p) out->HydFile = NULL; switch (out->Hydflag) { case SCRATCH: - getTmpName(out->HydFname); + getTmpName(p, out->HydFname); out->HydFile = fopen(out->HydFname, "w+b"); break; case SAVE: @@ -3512,7 +3534,7 @@ int openoutfile(EN_Project *p) // else if ( (OutFile = tmpfile()) == NULL) else { - getTmpName(out->OutFname); + getTmpName(p, out->OutFname); if ((out->OutFile = fopen(out->OutFname, "w+b")) == NULL) { writecon(FMT08); @@ -3530,7 +3552,7 @@ int openoutfile(EN_Project *p) if (!errcode) { if (rep->Tstatflag != SERIES) { // if ( (TmpOutFile = tmpfile()) == NULL) errcode = 304; - getTmpName(out->TmpFname); + getTmpName(p, out->TmpFname); out->TmpOutFile = fopen(out->TmpFname, "w+b"); if (out->TmpOutFile == NULL) errcode = 304; @@ -3837,7 +3859,7 @@ void freedata(EN_Project *p) ---------------------------------------------------------------- */ -char *getTmpName(char *fname) +char *getTmpName(EN_Project *p, char *fname) // // Input: fname = file name string // Output: returns pointer to file name @@ -3847,25 +3869,24 @@ char *getTmpName(char *fname) #ifdef _WIN32 - char* name = NULL; + char* name = NULL; - // --- use Windows _tempnam function to get a pointer to an - // unused file name that begins with "en" - name = _tempnam(NULL, "en"); - if (name == NULL) return NULL; + // --- use Windows _tempnam function to get a pointer to an + // unused file name that begins with "en" + name = _tempnam(NULL, "en"); + if (name == NULL) return NULL; - // --- copy the file name to fname - if (strlen(name) < MAXFNAME) strncpy(fname, name, MAXFNAME); - else fname = NULL; + // --- copy the file name to fname + if (strlen(name) < MAXFNAME) strncpy(fname, name, MAXFNAME); + else fname = NULL; - // --- free the pointer returned by _tempnam - if (name) free(name); + // --- free the pointer returned by _tempnam + if (name) free(name); /* /////////////////// DEPRECATED ///////////////////////////////////// // --- for Windows systems: #ifdef WINDOWS - out_file_t *out = &p->out_files; // --- use system function tmpnam() to create a temporary file name char name[MAXFNAME + 1]; int n; diff --git a/src/funcs.h b/src/funcs.h index 79490c6..b65f194 100755 --- a/src/funcs.h +++ b/src/funcs.h @@ -31,27 +31,28 @@ AUTHOR: L. Rossman #include "types.h" -void initpointers(EN_Project *p); /* Initializes pointers */ -int allocdata(EN_Project *p); /* Allocates memory */ -void freeTmplist(STmplist *); /* Frees items in linked list */ -void freeFloatlist(SFloatlist *); /* Frees list of floats */ -void freedata(EN_Project *p); /* Frees allocated memory */ -int openfiles(EN_Project *p, char *,char *,char *); /* Opens input & report files */ -int openhydfile(EN_Project *p); /* Opens hydraulics file */ -int openoutfile(EN_Project *p); /* Opens binary output file */ -int strcomp(char *, char *); /* Compares two strings */ -char* getTmpName(char* fname); /* Gets temporary file name */ -double interp(int, double *,double *, double); /* Interpolates a data curve */ +void initpointers(EN_Project *pr); /* Initializes pointers */ +int allocdata(EN_Project *pr); /* Allocates memory */ +void freeTmplist(STmplist *); /* Frees items in linked list */ +void freeFloatlist(SFloatlist *); /* Frees list of floats */ +void freedata(EN_Project *pr); /* Frees allocated memory */ +int openfiles(EN_Project *pr, char *,char *,char *); /* Opens input & report files */ +int openhydfile(EN_Project *pr); /* Opens hydraulics file */ +int openoutfile(EN_Project *pr); /* Opens binary output file */ +int strcomp(char *, char *); /* Compares two strings */ +char* getTmpName(EN_Project *p, char* fname); /* Gets temporary file name */ +double interp(int n, double x[], double y[], + double xx); /* Interpolates a data curve */ -int findnode(EN_Network *n, char *); /* Finds node's index from ID */ -int findlink(EN_Network *n, char *); /* Finds link's index from ID */ -int findtank(EN_Network *n, int); /* Find tank index from node index */ // (AH) -int findvalve(EN_Network *n, int); /* Find valve index from node index */ // (AH) -int findpump(EN_Network *n, int); /* Find pump index from node index */ // (AH) -char *geterrmsg(int errcode, char *msg); /* Gets text of error message */ -void errmsg(EN_Project *p, int); /* Reports program error */ -void writecon(char *); /* Writes text to console */ -void writewin(void (*vp)(char *), char *); /* Passes text to calling app */ +int findnode(EN_Network *n, char *); /* Finds node's index from ID */ +int findlink(EN_Network *n, char *); /* Finds link's index from ID */ +int findtank(EN_Network *n, int); /* Find tank index from node index */ // (AH) +int findvalve(EN_Network *n, int); /* Find valve index from node index */ // (AH) +int findpump(EN_Network *n, int); /* Find pump index from node index */ // (AH) +char *geterrmsg(int errcode, char *msg); /* Gets text of error message */ +void errmsg(EN_Project *p, int); /* Reports program error */ +void writecon(char *); /* Writes text to console */ +void writewin(void (*vp)(char *), char *); /* Passes text to calling app */ /* ------- INPUT1.C --------------------*/ int getdata(EN_Project *pr); /* Gets network data */ @@ -77,8 +78,9 @@ int getpatterns(EN_Project *pr); /* Gets pattern data from li int getcurves(EN_Project *pr); /* Gets curve data from list */ int findmatch(char *, char *[]); /* Finds keyword in line */ int match(const char *, const char *); /* Checks for word match */ -int gettokens(char *s, char** Tok, int maxToks, char *comment); /* Tokenizes input line */ -int getfloat(char *, double *); /* Converts string to double */ +int gettokens(char *s, char** Tok, int maxToks, + char *comment); /* Tokenizes input line */ +int getfloat(char *, double *); /* Converts string to double */ double hour(char *, char *); /* Converts time to hours */ int setreport(EN_Project *pr, char *); /* Processes reporting command*/ void inperrmsg(EN_Project *pr, int,int,char *); /* Input error message */ @@ -111,7 +113,8 @@ int powercurve(double, double, double, /* Coeffs. of power pump cur double, double, double *, double *, double *); int valvecheck(EN_Project *pr, int, int, int); /* Checks valve placement */ -void changestatus(EN_Network *net, int, StatType, double); /* Changes status of a link */ +void changestatus(EN_Network *net, int, StatType, + double); /* Changes status of a link */ /* -------------- RULES.C --------------*/ void initrules(rules_t *rules); /* Initializes rule base */ @@ -136,13 +139,14 @@ void writerelerr(EN_Project *pr, int, double); /* Writes convergence error void writestatchange(EN_Project *pr, int,char,char); /* Writes link status change */ void writecontrolaction(EN_Project *pr, int, int); /* Writes control action taken*/ void writeruleaction(EN_Project *pr, int, char *); /* Writes rule action taken */ -int writehydwarn(EN_Project *pr, int,double); /* Writes hydraulic warnings */ +int writehydwarn(EN_Project *pr, int,double); /* Writes hydraulic warnings */ void writehyderr(EN_Project *pr, int); /* Writes hydraulic error msg.*/ int disconnected(EN_Project *pr); /* Checks for disconnections */ void marknodes(EN_Project *pr, int, int *, char *); /* Identifies connected nodes */ void getclosedlink(EN_Project *pr, int, char *); /* Finds a disconnecting link */ void writelimits(EN_Project *pr, int,int); /* Writes reporting limits */ -int checklimits(report_options_t *rep, double *,int,int); /* Checks variable limits */ +int checklimits(report_options_t *rep, double *, + int,int); /* Checks variable limits */ void writetime(EN_Project *pr, char *); /* Writes current clock time */ char *clocktime(char *, long); /* Converts time to hrs:min */ char *fillstr(char *, char, int); /* Fills string with character*/ @@ -150,10 +154,7 @@ int getnodetype(EN_Network *net, int); /* Determines node typ /* --------- HYDRAUL.C -----------------*/ int openhyd(EN_Project *pr); /* Opens hydraulics solver */ - -/*** Updated 3/1/01 ***/ void inithyd(EN_Project *pr, int initFlags); /* Re-sets initial conditions */ - int runhyd(EN_Project *pr, long *); /* Solves 1-period hydraulics */ int nexthyd(EN_Project *pr, long *); /* Moves to next time period */ void closehyd(EN_Project *pr); /* Closes hydraulics solver */ @@ -164,64 +165,31 @@ void initlinkflow(EN_Project *pr, int, char, void setlinkflow(EN_Project *pr, int, double); /* Sets link flow via headloss*/ void setlinkstatus(EN_Project *pr, int, char, StatType *, double *); /* Sets link status */ - void setlinksetting(EN_Project *pr, int, double, StatType *, double *); /* Sets pump/valve setting */ -void resistance(EN_Project *pr, int); /* Computes resistance coeff. */ void demands(EN_Project *pr); /* Computes current demands */ int controls(EN_Project *pr); /* Controls link settings */ long timestep(EN_Project *pr); /* Computes new time step */ int tanktimestep(EN_Project *pr, long *); /* Time till tanks fill/drain */ void controltimestep(EN_Project *pr, long *); /* Time till control action */ void ruletimestep(EN_Project *pr, long *); /* Time till rule action */ + void addenergy(EN_Project *pr, long); /* Accumulates energy usage */ void getenergy(EN_Project *pr, int, double *, double *); /* Computes link energy use */ + void tanklevels(EN_Project *pr, long); /* Computes new tank levels */ double tankvolume(EN_Project *pr, int,double); /* Finds tank vol. from grade */ double tankgrade(EN_Project *pr, int,double); /* Finds tank grade from vol. */ -int netsolve(EN_Project *pr, int *,double *); /* Solves network equations */ -int badvalve(EN_Project *pr, int); /* Checks for bad valve */ -int valvestatus(EN_Project *pr); /* Updates valve status */ -int linkstatus(EN_Project *pr); /* Updates link status */ -StatType cvstatus(EN_Project *pr, StatType, - double,double); /* Updates CV status */ -StatType pumpstatus(EN_Project *pr, int,double); /* Updates pump status */ -StatType prvstatus(EN_Project *pr, int,StatType, - double,double,double); /* Updates PRV status */ - -StatType psvstatus(EN_Project *pr, int,StatType, - double,double,double); /* Updates PSV status */ - -StatType fcvstatus(EN_Project *pr, int,StatType, - double,double); /* Updates FCV status */ - -void tankstatus(EN_Project *pr, int,int,int); /* Checks if tank full/empty */ -int pswitch(EN_Project *pr); /* Pressure switch controls */ -double newflows(EN_Project *pr); /* Updates link flows */ -void newcoeffs(EN_Project *pr); /* Computes matrix coeffs. */ -void linkcoeffs(EN_Project *pr); /* Computes link coeffs. */ -void nodecoeffs(EN_Project *pr); /* Computes node coeffs. */ -void valvecoeffs(EN_Project *pr); /* Computes valve coeffs. */ -void pipecoeff(EN_Project *pr, int); /* Computes pipe coeff. */ -double DWcoeff(EN_Project *pr, int, double *); /* Computes D-W coeff. */ -void pumpcoeff(EN_Project *pr, int); /* Computes pump coeff. */ -/*** Updated 10/25/00 ***/ -/*** Updated 12/29/00 ***/ -void curvecoeff(EN_Project *pr, int,double, - double *,double *); /* Computes curve coeffs. */ - +/* ----------- HYDSOLVER.C - ----------*/ +int hydsolve(EN_Project *pr, int *,double *); /* Solves network equations */ -void gpvcoeff(EN_Project *pr, int iLink); /* Computes GPV coeff. */ -void pbvcoeff(EN_Project *pr, int iLink); /* Computes PBV coeff. */ -void tcvcoeff(EN_Project *pr, int iLink); /* Computes TCV coeff. */ -void prvcoeff(EN_Project *pr, int iLink, int n1, int n2); /* Computes PRV coeff. */ -void psvcoeff(EN_Project *pr, int iLink, int n1, int n2); /* Computes PSV coeff. */ -void fcvcoeff(EN_Project *pr, int iLink, int n1, int n2); /* Computes FCV coeff. */ -void emittercoeffs(EN_Project *pr); /* Computes emitter coeffs. */ -double emitflowchange(EN_Project *pr, int); /* Computes new emitter flow */ +/* ----------- HYDCOEFFS.C --------------*/ +void resistcoeff(EN_Project *pr, int k); /* Finds pipe flow resistance */ +void hlosscoeff(EN_Project *pr, int k); /* Finds link head loss coeff */ +void matrixcoeffs(EN_Project *pr); /* Finds hyd. matrix coeffs. */ /* ----------- SMATRIX.C ---------------*/ int createsparse(EN_Project *pr); /* Creates sparse matrix */ @@ -242,7 +210,7 @@ int storesparse(EN_Project *pr, int); /* Stores sparse matrix int ordersparse(hydraulics_t *h, int); /* Orders matrix storage */ void transpose(int,int *,int *, /* Transposes sparse matrix */ int *,int *,int *,int *,int *); -int linsolve(solver_t *s, int); /* Solves set of linear eqns. */ +int linsolve(solver_t *s, int); /* Solves set of inear eqns. */ /* ----------- QUALITY.C ---------------*/ int openqual(EN_Project *pr); /* Opens WQ solver system */ @@ -278,11 +246,10 @@ double pipereact(EN_Project *pr, int,double, double tankreact(EN_Project *pr, double,double, double,long); /* Reacts water in a tank */ double bulkrate(EN_Project *pr, double,double, - double); /* Finds bulk reaction rate */ + double); /* Finds bulk reaction rate */ double wallrate(EN_Project *pr, double,double, double,double); /* Finds wall reaction rate */ - /* ------------ OUTPUT.C ---------------*/ int savenetdata(EN_Project *pr); /* Saves basic data to file */ int savehyd(EN_Project *pr, long *); /* Saves hydraulic solution */ diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c new file mode 100644 index 0000000..9b830d9 --- /dev/null +++ b/src/hydcoeffs.c @@ -0,0 +1,933 @@ +/* +********************************************************************* + +HYDCOEFFS.C -- hydraulic coefficients for the EPANET Program + +******************************************************************* +*/ + +#include +#include +#ifndef __APPLE__ +#include +#else +#include +#endif +#include +#include "types.h" +#include "funcs.h" + +// Constants used for computing Darcy-Weisbach friction factor +const double A1 = 0.314159265359e04; // 1000*PI +const double A2 = 0.157079632679e04; // 500*PI +const double A3 = 0.502654824574e02; // 16*PI +const double A4 = 6.283185307; // 2*PI +const double A8 = 4.61841319859; // 5.74*(PI/4)^.9 +const double A9 = -8.685889638e-01; // -2/ln(10) +const double AA = -1.5634601348; // -2*.9*2/ln(10) +const double AB = 3.28895476345e-03; // 5.74/(4000^.9) +const double AC = -5.14214965799e-03; // AA*AB + +// External functions +//void resistcoeff(EN_Project *pr, int k); +//void hlosscoeff(EN_Project *pr, int k); +//void matrixcoeffs(EN_Project *pr); + +// Local functions +static void linkcoeffs(EN_Project *pr); +static void nodecoeffs(EN_Project *pr); +static void valvecoeffs(EN_Project *pr); +static void emittercoeffs(EN_Project *pr); + +static void pipecoeff(EN_Project *pr, int k); +static void DWpipecoeff(EN_Project *pr, int k); +static double frictionFactor(EN_Project *pr, int k, double *dfdq); + +static void pumpcoeff(EN_Project *pr, int k); +static void curvecoeff(EN_Project *pr, int i, double q, double *h0, double *r); + +static void valvecoeff(EN_Project *pr, int k); +static void gpvcoeff(EN_Project *pr, int k); +static void pbvcoeff(EN_Project *pr, int k); +static void tcvcoeff(EN_Project *pr, int k); +static void prvcoeff(EN_Project *pr, int k, int n1, int n2); +static void psvcoeff(EN_Project *pr, int k, int n1, int n2); +static void fcvcoeff(EN_Project *pr, int k, int n1, int n2); + + + +void resistcoeff(EN_Project *pr, int k) +/* +**-------------------------------------------------------------------- +** Input: k = link index +** Output: none +** Purpose: computes link flow resistance coefficient +**-------------------------------------------------------------------- +*/ +{ + double e, d, L; + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + Slink *link = &net->Link[k]; + + switch (link->Type) { + + // ... Link is a pipe. Compute resistance based on headloss formula. + // Friction factor for D-W formula gets included during head loss + // calculation. + case EN_CVPIPE: + case EN_PIPE: + e = link->Kc; // Roughness coeff. + d = link->Diam; // Diameter + L = link->Len; // Length + switch (hyd->Formflag) + { + case HW: + link->R = 4.727*L / pow(e, hyd->Hexp) / pow(d, 4.871); + break; + case DW: + link->R = L / 2.0 / 32.2 / d / SQR(PI*SQR(d) / 4.0); + break; + case CM: + link->R = SQR(4.0*e / (1.49*PI*d*d)) * pow((d / 4.0), -1.333)*L; + } + break; + + // ... Link is a pump. Use huge resistance. + case EN_PUMP: + link->R = CBIG; + break; + + // ... For all other links (e.g. valves) use a small resistance + default: + link->R = CSMALL; + break; + } +} + + +void hlosscoeff(EN_Project *pr, int k) +/* +**-------------------------------------------------------------- +** Input: k = link index +** Output: none +** Purpose: computes P and Y coefficients for a link +**-------------------------------------------------------------- +*/ +{ + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + Slink *link = &net->Link[k]; + + switch (link->Type) + { + case EN_CVPIPE: + case EN_PIPE: + pipecoeff(pr, k); + break; + case EN_PUMP: + pumpcoeff(pr, k); + break; + case EN_PBV: + pbvcoeff(pr, k); + break; + case EN_TCV: + tcvcoeff(pr, k); + break; + case EN_GPV: + gpvcoeff(pr, k); + break; + case EN_FCV: + case EN_PRV: + case EN_PSV: + if (hyd->LinkSetting[k] == MISSING) { + valvecoeff(pr, k); + } + else sol->P[k] = 0.0; + } +} + + +void matrixcoeffs(EN_Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: computes coefficients of linearized network eqns. +**-------------------------------------------------------------- +*/ +{ + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + EN_Network *net = &pr->network; + + memset(sol->Aii, 0, (net->Nnodes + 1) * sizeof(double)); + memset(sol->Aij, 0, (hyd->Ncoeffs + 1) * sizeof(double)); + memset(sol->F, 0, (net->Nnodes + 1) * sizeof(double)); + memset(hyd->X_tmp, 0, (net->Nnodes + 1) * sizeof(double)); + linkcoeffs(pr); + emittercoeffs(pr); + nodecoeffs(pr); + valvecoeffs(pr); +} + + +void linkcoeffs(EN_Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: computes matrix coefficients for links +**-------------------------------------------------------------- +*/ +{ + int k, n1, n2; + + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + Slink *link; + + // Examine each link of network */ + for (k = 1; k <= net->Nlinks; k++) + { + if (sol->P[k] == 0.0) continue; + link = &net->Link[k]; + n1 = link->N1; // Start node of link + n2 = link->N2; // End node of link + + // Update net nodal inflows (X), solution matrix (A) and RHS array (F) + // (Use covention that flow out of node is (-), flow into node is (+)) + hyd->X_tmp[n1] -= hyd->LinkFlows[k]; + hyd->X_tmp[n2] += hyd->LinkFlows[k]; + + // Off-diagonal coeff. + sol->Aij[sol->Ndx[k]] -= sol->P[k]; + + // Node n1 is junction + if (n1 <= net->Njuncs) + { + sol->Aii[sol->Row[n1]] += sol->P[k]; // Diagonal coeff. + sol->F[sol->Row[n1]] += sol->Y[k]; // RHS coeff. + } + + // Node n1 is a tank + else { + sol->F[sol->Row[n2]] += (sol->P[k] * hyd->NodeHead[n1]); + } + + // Node n2 is junction + if (n2 <= net->Njuncs) { + sol->Aii[sol->Row[n2]] += sol->P[k]; // Diagonal coeff. + sol->F[sol->Row[n2]] -= sol->Y[k]; // RHS coeff. + } + + // Node n2 is a tank + else { + sol->F[sol->Row[n1]] += (sol->P[k] * hyd->NodeHead[n2]); + } + } +} + + +void nodecoeffs(EN_Project *pr) +/* +**---------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: completes calculation of nodal flow imbalance (X) +** & flow correction (F) arrays +**---------------------------------------------------------------- +*/ +{ + int i; + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + EN_Network *net = &pr->network; + + // For junction nodes, subtract demand flow from net + // flow imbalance & add imbalance to RHS array F. + for (i = 1; i <= net->Njuncs; i++) + { + hyd->X_tmp[i] -= hyd->NodeDemand[i]; + sol->F[sol->Row[i]] += hyd->X_tmp[i]; + } +} + + +void valvecoeffs(EN_Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: computes matrix coeffs. for PRVs, PSVs & FCVs +** whose status is not fixed to OPEN/CLOSED +**-------------------------------------------------------------- +*/ +{ + int i, k, n1, n2; + + hydraulics_t *hyd = &pr->hydraulics; + EN_Network *net = &pr->network; + Slink *link; + Svalve *valve; + + // Examine each valve + for (i = 1; i <= net->Nvalves; i++) + { + valve = &net->Valve[i]; + k = valve->Link; // Link index of valve + link = &net->Link[k]; + if (hyd->LinkSetting[k] == MISSING) { + continue; // Valve status fixed + } + n1 = link->N1; // Start & end nodes + n2 = link->N2; + switch (link->Type) // Call valve-specific function + { + case EN_PRV: + prvcoeff(pr, k, n1, n2); + break; + case EN_PSV: + psvcoeff(pr, k, n1, n2); + break; + case EN_FCV: + fcvcoeff(pr, k, n1, n2); + break; + default: continue; + } + } +} + + +void emittercoeffs(EN_Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: computes matrix coeffs. for emitters +** +** Note: Emitters consist of a fictitious pipe connected to +** a fictitious reservoir whose elevation equals that +** of the junction. The headloss through this pipe is +** Ke*(Flow)^hyd->Qexp, where Ke = emitter headloss coeff. +**-------------------------------------------------------------- +*/ +{ + int i; + double ke; + double p; + double q; + double y; + double z; + + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + EN_Network *net = &pr->network; + Snode *node; + + for (i = 1; i <= net->Njuncs; i++) + { + node = &net->Node[i]; + if (node->Ke == 0.0) { + continue; + } + ke = MAX(CSMALL, node->Ke); // emitter coeff. + q = hyd->EmitterFlows[i]; // emitter flow + z = ke * pow(ABS(q), hyd->Qexp); // emitter head loss + p = hyd->Qexp * z / ABS(q); // head loss gradient + if (p < hyd->RQtol) { + p = 1.0 / hyd->RQtol; + } + else { + p = 1.0 / p; // inverse head loss gradient + } + y = SGN(q)*z*p; // head loss / gradient + sol->Aii[sol->Row[i]] += p; // addition to main diagonal + sol->F[sol->Row[i]] += y + p * node->El; // addition to r.h.s. + hyd->X_tmp[i] -= q; // addition to net node inflow + } +} + + +void pipecoeff(EN_Project *pr, int k) +/* +**-------------------------------------------------------------- +** Input: k = link index +** Output: none +** Purpose: computes P & Y coefficients for pipe k +** +** P = inverse head loss gradient = 1/(dh/dQ) +** Y = flow correction term = h*P +**-------------------------------------------------------------- +*/ +{ + double hpipe, // Normal head loss + hml, // Minor head loss + ml, // Minor loss coeff. + p, // q*(dh/dq) + q, // Abs. value of flow + r; // Resistance coeff. + + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + Slink *link = &pr->network.Link[k]; + + // For closed pipe use headloss formula: h = CBIG*q + if (hyd->LinkStatus[k] <= CLOSED) + { + sol->P[k] = 1.0 / CBIG; + sol->Y[k] = hyd->LinkFlows[k]; + return; + } + + // ... head loss formula is Darcy-Weisbach + if (hyd->Formflag == DW) { + DWpipecoeff(pr, k); + return; + } + + // Evaluate headloss coefficients + q = ABS(hyd->LinkFlows[k]); // Absolute flow + ml = link->Km; // Minor loss coeff. + r = link->R; // Resistance coeff. + + // Use large P coefficient for small flow resistance product + if ( (r+ml)*q < hyd->RQtol) + { + sol->P[k] = 1.0 / hyd->RQtol; + sol->Y[k] = hyd->LinkFlows[k] / hyd->Hexp; + return; + } + + // Compute P and Y coefficients + hpipe = r*pow(q, hyd->Hexp); // Friction head loss + p = hyd->Hexp*hpipe; // Q*dh(friction)/dQ + if (ml > 0.0) + { + hml = ml*q*q; // Minor head loss + p += 2.0*hml; // Q*dh(Total)/dQ + } + else hml = 0.0; + p = hyd->LinkFlows[k] / p; // 1 / (dh/dQ) + sol->P[k] = ABS(p); + sol->Y[k] = p*(hpipe + hml); +} + + +void DWpipecoeff(EN_Project *pr, int k) +/* +**-------------------------------------------------------------- +** Input: k = link index +** Output: none +** Purpose: computes pipe head loss coeffs. for Darcy-Weisbach +** formula. +**-------------------------------------------------------------- +*/ +{ + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + Slink *link = &pr->network.Link[k]; + + double q = ABS(hyd->LinkFlows[k]); + double dfdq = 0.0; + double r, r1, f, ml, p, hloss; + + ml = link->Km; // Minor loss coeff. + r = link->R; // Resistance coeff. + f = frictionFactor(pr,k,&dfdq); // D-W friction factor + r1 = f*r+ml; + + // Use large P coefficient for small flow resistance product + if (r1*q < hyd->RQtol) + { + sol->P[k] = 1.0/hyd->RQtol; + sol->Y[k] = hyd->LinkFlows[k]/hyd->Hexp; + return; + } + + // Compute P and Y coefficients + hloss = r1*SQR(q); // Total head loss + p = 2.0*r1*q; // |dHloss/dQ| + // + dfdq*r*q*q; // Ignore df/dQ term + p = 1.0 / p; + sol->P[k] = p; + sol->Y[k] = SGN(hyd->LinkFlows[k]) * hloss * p; +} + + +double frictionFactor(EN_Project *pr, int k, double *dfdq) +/* +**-------------------------------------------------------------- +** Input: k = link index +** Output: returns friction factor and +** replaces dfdq (derivative of f w.r.t. flow) +** Purpose: computes Darcy-Weisbach friction factor and its +** derivative as a function of Reynolds Number (Re). +** +** Note: Current formulas for dfdq need to be corrected +** so dfdq returned as 0. +**-------------------------------------------------------------- +*/ +{ + double q, // Abs. value of flow + f; // Friction factor + double x1,x2,x3,x4, + y1,y2,y3, + fa,fb,r; + double s,w; + + hydraulics_t *hyd = &pr->hydraulics; + Slink *link = &pr->network.Link[k]; + + *dfdq = 0.0; + if (link->Type > EN_PIPE) + return(1.0); // Only apply to pipes + q = ABS(hyd->LinkFlows[k]); + s = hyd->Viscos * link->Diam; + w = q/s; // w = Re(Pi/4) + + // For Re >= 4000 use Colebrook Formula + if (w >= A1) + { + y1 = A8/pow(w,0.9); + y2 = link->Kc/(3.7*link->Diam) + y1; + y3 = A9*log(y2); + f = 1.0/SQR(y3); + /* *dfdq = (2.0+AA*y1/(y2*y3))*f; */ /* df/dq */ + } + + // For Re > 2000 use Interpolation Formula + else if (w > A2) + { + y2 = link->Kc/(3.7*link->Diam) + AB; + y3 = A9*log(y2); + fa = 1.0/SQR(y3); + fb = (2.0+AC/(y2*y3))*fa; + r = w/A2; + x1 = 7.0*fa - fb; + x2 = 0.128 - 17.0*fa + 2.5*fb; + x3 = -0.128 + 13.0*fa - (fb+fb); + x4 = r*(0.032 - 3.0*fa + 0.5*fb); + f = x1 + r*(x2 + r*(x3+x4)); + /* *dfdq = (x1 + x1 + r*(3.0*x2 + r*(4.0*x3 + 5.0*x4))); */ + } + + // For Re > 8 (Laminar Flow) use Hagen-Poiseuille Formula + else if (w > A4) + { + f = A3*s/q; // 16 * PI * Viscos * Diam / Flow + /* *dfdq = A3*s; */ + } + else + { + f = 8.0; + *dfdq = 0.0; + } + return(f); + +} + + +void pumpcoeff(EN_Project *pr, int k) +/* +**-------------------------------------------------------------- +** Input: k = link index +** Output: none +** Purpose: computes P & Y coeffs. for pump in link k +**-------------------------------------------------------------- +*/ +{ + int p; // Pump index + double h0, // Shutoff head + q, // Abs. value of flow + r, // Flow resistance coeff. + n; // Flow exponent coeff. + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + double setting = hyd->LinkSetting[k]; + Spump *pump; + + // Use high resistance pipe if pump closed or cannot deliver head + if (hyd->LinkStatus[k] <= CLOSED || setting == 0.0) { + sol->P[k] = 1.0 / CBIG; + sol->Y[k] = hyd->LinkFlows[k]; + return; + } + + q = ABS(hyd->LinkFlows[k]); + q = MAX(q, TINY); + + // Obtain reference to pump object + p = findpump(&pr->network, k); + pump = &pr->network.Pump[p]; + + // Get pump curve coefficients for custom pump curve. + if (pump->Ptype == CUSTOM) + { + // Find intercept (h0) & slope (r) of pump curve + // line segment which contains speed-adjusted flow. + curvecoeff(pr, pump->Hcurve, q / setting, &h0, &r); + + // Determine head loss coefficients. + pump->H0 = -h0; + pump->R = -r; + pump->N = 1.0; + } + + // Adjust head loss coefficients for pump speed. + h0 = SQR(setting) * pump->H0; + n = pump->N; + r = pump->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) + sol->P[k] = 1.0 / MAX(r, hyd->RQtol); + sol->Y[k] = hyd->LinkFlows[k] / n + sol->P[k] * h0; +} + + +void curvecoeff(EN_Project *pr, int i, double q, double *h0, double *r) +/* +**------------------------------------------------------------------- +** Input: i = curve index +** q = flow rate +** Output: *h0 = head at zero flow (y-intercept) +** *r = dHead/dFlow (slope) +** Purpose: computes intercept and slope of head v. flow curve +** at current flow. +**------------------------------------------------------------------- +*/ +{ + int k1, k2, npts; + double *x, *y; + Scurve *curve; + + // Remember that curve is stored in untransformed units + q *= pr->Ucf[FLOW]; + curve = &pr->network.Curve[i]; + x = curve->X; // x = flow + y = curve->Y; // y = head + npts = curve->Npts; + + // Find linear segment of curve that brackets flow q + k2 = 0; + while (k2 < npts && x[k2] < q) + k2++; + + if (k2 == 0) + k2++; + + else if (k2 == npts) + k2--; + + k1 = k2 - 1; + + // Compute slope and intercept of this segment + *r = (y[k2] - y[k1]) / (x[k2] - x[k1]); + *h0 = y[k1] - (*r)*x[k1]; + + // Convert units + *h0 = (*h0) / pr->Ucf[HEAD]; + *r = (*r) * pr->Ucf[FLOW] / pr->Ucf[HEAD]; +} + + +void gpvcoeff(EN_Project *pr, int k) +/* +**-------------------------------------------------------------- +** Input: k = link index +** Output: none +** Purpose: computes P & Y coeffs. for general purpose valve +**-------------------------------------------------------------- +*/ +{ + double h0, // Headloss curve intercept + q, // Abs. value of flow + r; // Flow resistance coeff. + + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + + // Treat as a pipe if valve closed + if (hyd->LinkStatus[k] == CLOSED) { + valvecoeff(pr, k); + } + + // Otherwise utilize headloss curve + // whose index is stored in K + else { + // Find slope & intercept of headloss curve. + q = ABS(hyd->LinkFlows[k]); + q = MAX(q, TINY); + curvecoeff(pr, (int)ROUND(hyd->LinkSetting[k]), q, &h0, &r); + + // Compute inverse headloss gradient (P) + // and flow correction factor (Y). + sol->P[k] = 1.0 / MAX(r, hyd->RQtol); + sol->Y[k] = sol->P[k] * (h0 + r*q) * SGN(hyd->LinkFlows[k]); + } +} + + +void pbvcoeff(EN_Project *pr, int k) +/* +**-------------------------------------------------------------- +** Input: k = link index +** Output: none +** Purpose: computes P & Y coeffs. for pressure breaker valve +**-------------------------------------------------------------- +*/ +{ + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + Slink *link = &pr->network.Link[k]; + + // If valve fixed OPEN or CLOSED then treat as a pipe + if (hyd->LinkSetting[k] == MISSING || hyd->LinkSetting[k] == 0.0) { + valvecoeff(pr, k); + } + + // If valve is active + else { + // Treat as a pipe if minor loss > valve setting + if (link->Km * SQR(hyd->LinkFlows[k]) > hyd->LinkSetting[k]) { + valvecoeff(pr, k); + } + // Otherwise force headloss across valve to be equal to setting + else { + sol->P[k] = CBIG; + sol->Y[k] = hyd->LinkSetting[k] * CBIG; + } + } +} + + +void tcvcoeff(EN_Project *pr, int k) +/* +**-------------------------------------------------------------- +** Input: k = link index +** Output: none +** Purpose: computes P & Y coeffs. for throttle control valve +**-------------------------------------------------------------- +*/ +{ + double km; + hydraulics_t *hyd = &pr->hydraulics; + Slink *link = &pr->network.Link[k]; + + // Save original loss coeff. for open valve + km = link->Km; + + // If valve not fixed OPEN or CLOSED, compute its loss coeff. + if (hyd->LinkSetting[k] != MISSING) { + link->Km = 0.02517 * hyd->LinkSetting[k] / (SQR(link->Diam)*SQR(link->Diam)); + } + + // Then apply usual valve formula + valvecoeff(pr, k); + + // Restore original loss coeff. + link->Km = km; +} + + +void prvcoeff(EN_Project *pr, int k, int n1, int n2) +/* +**-------------------------------------------------------------- +** Input: k = link index +** n1 = upstream node of valve +** n2 = downstream node of valve +** Output: none +** Purpose: computes solution matrix coeffs. for pressure +** reducing valves +**-------------------------------------------------------------- +*/ +{ + int i, j; // Rows of solution matrix + double hset; // Valve head setting + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + + i = sol->Row[n1]; // Matrix rows of nodes + j = sol->Row[n2]; + hset = pr->network.Node[n2].El + + hyd->LinkSetting[k]; // Valve setting + + if (hyd->LinkStatus[k] == ACTIVE) + { + + // Set coeffs. to force head at downstream + // node equal to valve setting & force flow + // to equal to flow imbalance at downstream node. + + sol->P[k] = 0.0; + sol->Y[k] = hyd->LinkFlows[k] + hyd->X_tmp[n2]; // Force flow balance + sol->F[j] += (hset * CBIG); // Force head = hset + sol->Aii[j] += CBIG; // at downstream node + if (hyd->X_tmp[n2] < 0.0) { + sol->F[i] += hyd->X_tmp[n2]; + } + return; + } + + // For OPEN, CLOSED, or XPRESSURE valve + // compute matrix coeffs. using the valvecoeff() function. + + valvecoeff(pr, k); + sol->Aij[sol->Ndx[k]] -= sol->P[k]; + sol->Aii[i] += sol->P[k]; + sol->Aii[j] += sol->P[k]; + sol->F[i] += (sol->Y[k] - hyd->LinkFlows[k]); + sol->F[j] -= (sol->Y[k] - hyd->LinkFlows[k]); +} + + +void psvcoeff(EN_Project *pr, int k, int n1, int n2) +/* +**-------------------------------------------------------------- +** Input: k = link index +** n1 = upstream node of valve +** n2 = downstream node of valve +** Output: none +** Purpose: computes solution matrix coeffs. for pressure +** sustaining valve +**-------------------------------------------------------------- +*/ +{ + int i, j; // Rows of solution matrix + double hset; // Valve head setting + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + + i = sol->Row[n1]; // Matrix rows of nodes + j = sol->Row[n2]; + hset = pr->network.Node[n1].El + + hyd->LinkSetting[k]; // Valve setting + + if (hyd->LinkStatus[k] == ACTIVE) + { + // Set coeffs. to force head at upstream + // node equal to valve setting & force flow + // equal to flow imbalance at upstream node. + + sol->P[k] = 0.0; + sol->Y[k] = hyd->LinkFlows[k] - hyd->X_tmp[n1]; // Force flow balance + sol->F[i] += (hset * CBIG); // Force head = hset + sol->Aii[i] += CBIG; // at upstream node + if (hyd->X_tmp[n1] > 0.0) { + sol->F[j] += hyd->X_tmp[n1]; + } + return; + } + + // For OPEN, CLOSED, or XPRESSURE valve + // compute matrix coeffs. using the valvecoeff() function. + + valvecoeff(pr, k); + sol->Aij[sol->Ndx[k]] -= sol->P[k]; + sol->Aii[i] += sol->P[k]; + sol->Aii[j] += sol->P[k]; + sol->F[i] += (sol->Y[k] - hyd->LinkFlows[k]); + sol->F[j] -= (sol->Y[k] - hyd->LinkFlows[k]); +} + + +void fcvcoeff(EN_Project *pr, int k, int n1, int n2) +/* +**-------------------------------------------------------------- +** Input: k = link index +** n1 = upstream node of valve +** n2 = downstream node of valve +** Output: none +** Purpose: computes solution matrix coeffs. for flow control +** valve +**-------------------------------------------------------------- +*/ +{ + int i, j; // Rows in solution matrix + double q; // Valve flow setting + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + + q = hyd->LinkSetting[k]; + i = hyd->solver.Row[n1]; + j = hyd->solver.Row[n2]; + + // If valve active, break network at valve and treat + // flow setting as external demand at upstream node + // and external supply at downstream node. + + if (hyd->LinkStatus[k] == ACTIVE) + { + hyd->X_tmp[n1] -= q; + sol->F[i] -= q; + hyd->X_tmp[n2] += q; + sol->F[j] += q; + sol->P[k] = 1.0 / CBIG; + sol->Aij[sol->Ndx[k]] -= sol->P[k]; + sol->Aii[i] += sol->P[k]; + sol->Aii[j] += sol->P[k]; + sol->Y[k] = hyd->LinkFlows[k] - q; + } + + // Otherwise treat valve as an open pipe + + else + { + valvecoeff(pr, k); + sol->Aij[sol->Ndx[k]] -= sol->P[k]; + sol->Aii[i] += sol->P[k]; + sol->Aii[j] += sol->P[k]; + sol->F[i] += (sol->Y[k] - hyd->LinkFlows[k]); + sol->F[j] -= (sol->Y[k] - hyd->LinkFlows[k]); + } +} + + +void valvecoeff(EN_Project *pr, int k) +/* +**-------------------------------------------------------------- +** Input: k = link index +** Output: none +** Purpose: computes solution matrix coeffs. for a completely +** open, closed, or throttled control valve. +**-------------------------------------------------------------- +*/ +{ + double p; + + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + Slink *link = &net->Link[k]; + + double flow = hyd->LinkFlows[k]; + + // Valve is closed. Use a very small matrix coeff. + if (hyd->LinkStatus[k] <= CLOSED) + { + sol->P[k] = 1.0 / CBIG; + sol->Y[k] = flow; + return; + } + + // Account for any minor headloss through the valve + if (link->Km > 0.0) + { + p = 2.0 * link->Km * fabs(flow); + if (p < hyd->RQtol) { + p = hyd->RQtol; + } + sol->P[k] = 1.0 / p; + sol->Y[k] = flow / 2.0; + } + else + { + sol->P[k] = 1.0 / hyd->RQtol; + sol->Y[k] = flow; + } +} diff --git a/src/hydraul.c b/src/hydraul.c index f4eb45c..545ed37 100755 --- a/src/hydraul.c +++ b/src/hydraul.c @@ -29,13 +29,13 @@ AUTHOR: L. Rossman closehyd() -- called from ENcloseH() in EPANET.C tankvolume() -- called from ENsetnodevalue() in EPANET.C setlinkstatus(), - setlinksetting(), - resistance()-- all called from ENsetlinkvalue() in EPANET.C + setlinksetting() -- all called from ENsetlinkvalue() in EPANET.C External functions called by this module are: createsparse() -- see SMATRIX.C freesparse() -- see SMATRIX.C - linsolve() -- see SMATRIX.C + resistcoeff() -- see HYDCOEFFS.C + hydsolve() -- see HYDSOLVER.C checkrules() -- see RULES.C interp() -- see EPANET.C savehyd() -- see OUTPUT.C @@ -54,34 +54,11 @@ AUTHOR: L. Rossman #include #endif #include -#include "hash.h" #include "text.h" #include "types.h" -#include "epanet2.h" #include "funcs.h" -#define EXTERN extern -#include "vars.h" #define QZERO 1.e-6 /* Equivalent to zero flow */ -#define CBIG 1.e8 /* Big coefficient */ -#define CSMALL 1.e-6 /* Small coefficient */ - -/* Constants used for computing Darcy-Weisbach friction factor */ -#define A1 0.314159265359e04 /* 1000*PI */ -#define A2 0.157079632679e04 /* 500*PI */ -#define A3 0.502654824574e02 /* 16*PI */ -#define A4 6.283185307 /* 2*PI */ -#define A8 4.61841319859 /* 5.74*(PI/4)^.9 */ -#define A9 -8.685889638e-01 /* -2/ln(10) */ -#define AA -1.5634601348 /* -2*.9*2/ln(10) */ -#define AB 3.28895476345e-03 /* 5.74/(4000^.9) */ -#define AC -5.14214965799e-03 /* AA*AB */ - -/*** Updated 3/1/01 ***/ - -/* Function to find flow coeffs. through open/closed valves */ -void valvecoeff(EN_Project *pr, int k); - int openhyd(EN_Project *pr) /* @@ -92,15 +69,15 @@ int openhyd(EN_Project *pr) *-------------------------------------------------------------- */ { - int i; - int errcode = 0; - ERRCODE(createsparse(pr)); /* See SMATRIX.C */ - ERRCODE(allocmatrix(pr)); /* Allocate solution matrices */ - for (i=1; i <= pr->network.Nlinks; i++) { /* Initialize flows */ - Slink *link = &pr->network.Link[i]; - initlinkflow(pr, i, link->Stat, link->Kc); - } - return(errcode); + int i; + int errcode = 0; + ERRCODE(createsparse(pr)); /* See SMATRIX.C */ + ERRCODE(allocmatrix(pr)); /* Allocate solution matrices */ + for (i=1; i <= pr->network.Nlinks; i++) { /* Initialize flows */ + Slink *link = &pr->network.Link[i]; + initlinkflow(pr, i, link->Stat, link->Kc); + } + return(errcode); } @@ -115,72 +92,80 @@ void inithyd(EN_Project *pr, int initflag) **-------------------------------------------------------------- */ { - int i,j; + int i,j; - time_options_t *time = &pr->time_options; - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - out_file_t *out = &pr->out_files; + time_options_t *time = &pr->time_options; + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + out_file_t *out = &pr->out_files; + Stank *tank; + Slink *link; + Spump *pump; - /* Initialize tanks */ - for (i=1; i <= net->Ntanks; i++) { - Stank *tank = &net->Tank[i]; - tank->V = tank->V0; - hyd->NodeHead[tank->Node] = tank->H0; - hyd->NodeDemand[tank->Node] = 0.0; - hyd->OldStat[net->Nlinks+i] = TEMPCLOSED; - } - - /* Initialize emitter flows */ - memset(hyd->EmitterFlows,0,(net->Nnodes+1)*sizeof(double)); - for (i=1; i <= net->Njuncs; i++) { - if (net->Node[i].Ke > 0.0) { - hyd->EmitterFlows[i] = 1.0; + /* Initialize tanks */ + for (i=1; i <= net->Ntanks; i++) { + tank = &net->Tank[i]; + tank->V = tank->V0; + hyd->NodeHead[tank->Node] = tank->H0; + hyd->NodeDemand[tank->Node] = 0.0; + hyd->OldStat[net->Nlinks+i] = TEMPCLOSED; } - } - /* Initialize links */ - for (i=1; i <= net->Nlinks; i++) { - Slink *link = &net->Link[i]; - /* Initialize status and setting */ - hyd->LinkStatus[i] = link->Stat; - hyd->LinkSetting[i] = link->Kc; - /* Start active control valves in ACTIVE position */ - if ( - (link->Type == EN_PRV || link->Type == EN_PSV - || link->Type == EN_FCV) - && (link->Kc != MISSING) - ) hyd->LinkStatus[i] = ACTIVE; + /* Initialize emitter flows */ + memset(hyd->EmitterFlows,0,(net->Nnodes+1)*sizeof(double)); + for (i=1; i <= net->Njuncs; i++) { + if (net->Node[i].Ke > 0.0) { + hyd->EmitterFlows[i] = 1.0; + } + } + + /* Initialize links */ + for (i=1; i <= net->Nlinks; i++) { + link = &net->Link[i]; + + /* Initialize status and setting */ + hyd->LinkStatus[i] = link->Stat; + hyd->LinkSetting[i] = link->Kc; + + /* Compute flow resistance */ + resistcoeff(pr, i); + + /* Start active control valves in ACTIVE position */ + if ( + (link->Type == EN_PRV || link->Type == EN_PSV + || link->Type == EN_FCV) && (link->Kc != MISSING) + ) hyd->LinkStatus[i] = ACTIVE; /*** Updated 3/1/01 ***/ - /* Initialize flows if necessary */ - if (hyd->LinkStatus[i] <= CLOSED) hyd->LinkFlows[i] = QZERO; - else if (ABS(hyd->LinkFlows[i]) <= QZERO || initflag > 0) - initlinkflow(pr, i, hyd->LinkStatus[i], hyd->LinkSetting[i]); + /* Initialize flows if necessary */ + if (hyd->LinkStatus[i] <= CLOSED) hyd->LinkFlows[i] = QZERO; + else if (ABS(hyd->LinkFlows[i]) <= QZERO || initflag > 0) + initlinkflow(pr, i, hyd->LinkStatus[i], hyd->LinkSetting[i]); - /* Save initial status */ - hyd->OldStat[i] = hyd->LinkStatus[i]; - } + /* Save initial status */ + hyd->OldStat[i] = hyd->LinkStatus[i]; + } - /* Reset pump energy usage */ - for (i=1; i <= net->Npumps; i++) - { - for (j=0; j<6; j++) { - net->Pump[i].Energy[j] = 0.0; - } - } + /* Initialize pump energy usage */ + for (i=1; i <= net->Npumps; i++) + { + pump = &net->Pump[i]; + for (j = 0; j < MAX_ENERGY_STATS; j++) { + pump->Energy[j] = 0.0; + } + } - /* Re-position hydraulics file */ - if (pr->save_options.Saveflag) { - fseek(out->HydFile,out->HydOffset,SEEK_SET); - } + /* Re-position hydraulics file */ + if (pr->save_options.Saveflag) { + fseek(out->HydFile,out->HydOffset,SEEK_SET); + } /*** Updated 3/1/01 ***/ - /* Initialize current time */ - hyd->Haltflag = 0; - time->Htime = 0; - time->Hydstep = 0; - time->Rtime = time->Rstep; + /* Initialize current time */ + hyd->Haltflag = 0; + time->Htime = 0; + time->Hydstep = 0; + time->Rtime = time->Rstep; } @@ -194,42 +179,42 @@ int runhyd(EN_Project *pr, long *t) **-------------------------------------------------------------- */ { - int iter; /* Iteration count */ - int errcode; /* Error code */ - double relerr; /* Solution accuracy */ + int iter; /* Iteration count */ + int errcode; /* Error code */ + double relerr; /* Solution accuracy */ - hydraulics_t *hyd = &pr->hydraulics; - time_options_t *time = &pr->time_options; - report_options_t *rep = &pr->report; + hydraulics_t *hyd = &pr->hydraulics; + time_options_t *time = &pr->time_options; + report_options_t *rep = &pr->report; - /* Find new demands & control actions */ - *t = time->Htime; - demands(pr); - controls(pr); + /* Find new demands & control actions */ + *t = time->Htime; + demands(pr); + controls(pr); - /* Solve network hydraulic equations */ - errcode = netsolve(pr,&iter,&relerr); - if (!errcode) { - /* Report new status & save results */ - if (rep->Statflag) { - writehydstat(pr,iter,relerr); - } + /* Solve network hydraulic equations */ + errcode = hydsolve(pr,&iter,&relerr); + if (!errcode) { + /* Report new status & save results */ + if (rep->Statflag) { + writehydstat(pr,iter,relerr); + } - /* solution info */ - hyd->relativeError = relerr; - hyd->iterations = iter; + /* solution info */ + hyd->relativeError = relerr; + hyd->iterations = iter; /*** Updated 3/1/01 ***/ - /* If system unbalanced and no extra trials */ - /* allowed, then activate the Haltflag. */ - if (relerr > hyd->Hacc && hyd->ExtraIter == -1) { - hyd->Haltflag = 1; - } + /* If system unbalanced and no extra trials */ + /* allowed, then activate the Haltflag. */ + if (relerr > hyd->Hacc && hyd->ExtraIter == -1) { + hyd->Haltflag = 1; + } - /* Report any warning conditions */ - if (!errcode) { - errcode = writehydwarn(pr,iter,relerr); - } + /* Report any warning conditions */ + if (!errcode) { + errcode = writehydwarn(pr,iter,relerr); + } } return(errcode); } /* end of runhyd */ @@ -459,7 +444,8 @@ void setlinkflow(EN_Project *pr, int k, double dh) dh = -dh * pr->Ucf[HEAD] / SQR(hyd->LinkSetting[k]); i = net->Pump[p].Hcurve; curve = &net->Curve[i]; - hyd->LinkFlows[k] = interp(curve->Npts,curve->Y,curve->X,dh) * hyd->LinkSetting[k] / pr->Ucf[FLOW]; + hyd->LinkFlows[k] = interp(curve->Npts,curve->Y,curve->X,dh) * + hyd->LinkSetting[k] / pr->Ucf[FLOW]; } /* Otherwise use inverse of power curve */ @@ -565,55 +551,6 @@ void setlinksetting(EN_Project *pr, int index, double value, StatType *s, doubl } -void resistance(EN_Project *pr, int k) -/* -**-------------------------------------------------------------------- -** Input: k = link index -** Output: none -** Purpose: computes link flow resistance -**-------------------------------------------------------------------- -*/ -{ - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - - double e,d,L; - - Slink *link = &net->Link[k]; - link->R = CSMALL; - - - switch (link->Type) { - /* Link is a pipe. Compute resistance based on headloss formula. */ - /* Friction factor for D-W formula gets included during solution */ - /* process in pipecoeff() function. */ - case EN_CVPIPE: - case EN_PIPE: - e = link->Kc; /* Roughness coeff. */ - d = link->Diam; /* Diameter */ - L = link->Len; /* Length */ - switch(hyd->Formflag) - { - case HW: - link->R = 4.727*L/pow(e,hyd->Hexp)/pow(d,4.871); - break; - case DW: - link->R = L/2.0/32.2/d/SQR(PI*SQR(d)/4.0); - break; - case CM: - link->R = SQR(4.0*e/(1.49*PI*d*d)) * pow((d/4.0),-1.333)*L; - } - break; - - /* Link is a pump. Use negligible resistance. */ - case EN_PUMP: - link->R = CBIG; //CSMALL; - break; - default: - break; - } -} - void demands(EN_Project *pr) /* @@ -699,7 +636,8 @@ int controls(EN_Project *pr) double k1, k2; char s1, s2; Slink *link; - + Scontrol *control; + EN_Network *net = &pr->network; time_options_t *top = &pr->time_options; hydraulics_t *hyd = &pr->hydraulics; @@ -708,13 +646,13 @@ int controls(EN_Project *pr) setsum = 0; for (i=1; i <= net->Ncontrols; i++) { - Scontrol *control = &net->Control[i]; + control = &net->Control[i]; /* Make sure that link is defined */ reset = 0; - if ( (k = control->Link) <= 0) { - continue; - } - link = &net->Link[k]; + if ( (k = control->Link) <= 0) { + continue; + } + link = &net->Link[k]; /* Link is controlled by tank level */ if ((n = control->Node) > 0 && n > net->Njuncs) { @@ -1028,7 +966,7 @@ void addenergy(EN_Project *pr, long hstep) { int i,j,k; long m,n; - double c0,c, /* Energy cost (cost/kwh) */ + double c0,c, /* Energy cost (cost/kwh) */ f0, /* Energy cost factor */ dt, /* Time interval (hr) */ e, /* Pump efficiency (fraction) */ @@ -1038,6 +976,7 @@ void addenergy(EN_Project *pr, long hstep) EN_Network *net = &pr->network; hydraulics_t *hyd = &pr->hydraulics; + Spump *pump; /* Determine current time interval in hours */ if (pr->time_options.Dur == 0) { @@ -1053,7 +992,7 @@ void addenergy(EN_Project *pr, long hstep) if (dt == 0.0) { return; } - n = (pr->time_options.Htime + pr->time_options.Pstart) / pr->time_options.Pstep; + n = (pr->time_options.Htime + pr->time_options.Pstart) / pr->time_options.Pstep; /* Compute default energy cost at current time */ c0 = hyd->Ecost; @@ -1067,17 +1006,17 @@ void addenergy(EN_Project *pr, long hstep) /* Examine each pump */ for (j=1; j <= net->Npumps; j++) { - Spump *pump = &net->Pump[j]; /* Skip closed pumps */ + pump = &net->Pump[j]; k = pump->Link; - if (hyd->LinkStatus[k] <= CLOSED) { + if (hyd->LinkStatus[k] <= CLOSED) { continue; - } + } q = MAX(QZERO, ABS(hyd->LinkFlows[k])); /* Find pump-specific energy cost */ - if (pump->Ecost > 0.0) { - c = pump->Ecost; + if (pump->Ecost > 0.0) { + c = pump->Ecost; } else { c = c0; @@ -1096,12 +1035,12 @@ void addenergy(EN_Project *pr, long hstep) psum += p; /* Update pump's cumulative statistics */ - pump->Energy[0] += dt; /* Time on-line */ - pump->Energy[1] += e*dt; /* Effic.-hrs */ - pump->Energy[2] += p/q*dt; /* kw/cfs-hrs */ - pump->Energy[3] += p*dt; /* kw-hrs */ - pump->Energy[4] = MAX(pump->Energy[4],p); - pump->Energy[5] += c*p*dt; /* cost-hrs. */ + pump->Energy[PCNT_ONLINE] += dt; + pump->Energy[PCNT_EFFIC] += e*dt; + pump->Energy[KWH_PER_FLOW] += p/q*dt; + pump->Energy[TOTAL_KWH] += p*dt; + pump->Energy[MAX_KW] = MAX(pump->Energy[MAX_KW],p); + pump->Energy[TOTAL_COST] += c*p*dt; } /* Update maximum kw value */ @@ -1119,12 +1058,16 @@ void getenergy(EN_Project *pr, int k, double *kw, double *eff) **---------------------------------------------------------------- */ { - int i,j; - double dh, q, e; - double q4eff; //q4eff=flow at nominal speed to compute efficiency + int i, // efficiency curve index + j; // pump index + double dh, // head across pump (ft) + q, // flow through pump (cfs) + e; // pump efficiency + double q4eff; // flow at nominal pump speed of 1.0 + double speed; // current speed setting Scurve *curve; - EN_Network *net = &pr->network; + EN_Network *net = &pr->network; hydraulics_t *hyd = &pr->hydraulics; Slink *link = &net->Link[k]; @@ -1147,12 +1090,15 @@ void getenergy(EN_Project *pr, int k, double *kw, double *eff) { j = findpump(net,k); e = hyd->Epump; + speed = hyd->LinkSetting[k]; if ( (i = net->Pump[j].Ecurve) > 0) - { - q4eff = q / hyd->LinkSetting[k]; + { + q4eff = q / speed * pr->Ucf[FLOW]; curve = &net->Curve[i]; - e = interp(curve->Npts,curve->X, curve->Y, q4eff * pr->Ucf[FLOW]); - } + e = interp(curve->Npts,curve->X, curve->Y, q4eff); + /* Sarbu and Borza pump speed adjustment */ + e = 100.0 - ((100.0-e) * pow(1.0/speed, 0.1)); + } e = MIN(e, 100.0); e = MAX(e, 1.0); e /= 100.0; @@ -1219,6 +1165,8 @@ double tankvolume(EN_Project *pr, int i, double h) int j; EN_Network *net = &pr->network; Stank *tank = &net->Tank[i]; + Scurve *curve; + /* Use level*area if no volume curve */ j = tank->Vcurve; if (j == 0) { @@ -1228,8 +1176,10 @@ double tankvolume(EN_Project *pr, int i, double h) /* If curve exists, interpolate on h to find volume v */ /* remembering that volume curve is in original units.*/ else { - Scurve *curve = &net->Curve[j]; - return(interp(curve->Npts, curve->X, curve->Y, (h - net->Node[tank->Node].El) * pr->Ucf[HEAD]) / pr->Ucf[VOLUME]); + curve = &net->Curve[j]; + return(interp(curve->Npts, curve->X, curve->Y, + (h - net->Node[tank->Node].El) * + pr->Ucf[HEAD]) / pr->Ucf[VOLUME]); } } /* End of tankvolume */ @@ -1265,1688 +1215,4 @@ double tankgrade(EN_Project *pr, int i, double v) } /* End of tankgrade */ -int netsolve(EN_Project *pr, int *iter, double *relerr) -/* -**------------------------------------------------------------------- -** Input: none -** Output: *iter = # of iterations to reach solution -** *relerr = convergence error in solution -** returns error code -** Purpose: solves network nodal equations for heads and flows -** using Todini's Gradient algorithm -** -*** Updated 9/7/00 *** -*** Updated 2.00.11 *** -*** Updated 2.00.12 *** -** Notes: Status checks on CVs, pumps and pipes to tanks are made -** every hyd->CheckFreq iteration, up until MaxCheck iterations -** are reached. Status checks on control valves are made -** every iteration if DampLimit = 0 or only when the -** convergence error is at or below DampLimit. If DampLimit -** is > 0 then future computed flow changes are only 60% of -** their full value. A complete status check on all links -** is made when convergence is achieved. If convergence is -** not achieved in hyd->MaxIter trials and hyd->ExtraIter > 0 then -** another hyd->ExtraIter trials are made with no status changes -** made to any links and a warning message is generated. -** -** This procedure calls linsolve() which appears in SMATRIX.C. -**------------------------------------------------------------------- -*/ -{ - int i; /* Node index */ - int errcode = 0; /* Node causing solution error */ - int nextcheck; /* Next status check trial */ - int maxtrials; /* Max. trials for convergence */ - double newerr; /* New convergence error */ - int valveChange; /* Valve status change flag */ - int statChange; - - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - solver_t *sol = &hyd->solver; - report_options_t *rep = &pr->report; - - /* Initialize status checking & relaxation factor */ - nextcheck = hyd->CheckFreq; - hyd->RelaxFactor = 1.0; - - /* Repeat iterations until convergence or trial limit is exceeded. */ - /* (hyd->ExtraIter used to increase trials in case of status cycling.) */ - if (pr->report.Statflag == FULL) { - writerelerr(pr,0,0); - } - maxtrials = hyd->MaxIter; - if (hyd->ExtraIter > 0) { - maxtrials += hyd->ExtraIter; - } - *iter = 1; - while (*iter <= maxtrials) { - /* - ** Compute coefficient matrices A & F and solve A*H = F - ** where H = heads, A = Jacobian coeffs. derived from - ** head loss gradients, & F = flow correction terms. - ** Solution for H is returned in F from call to linsolve(). - */ - newcoeffs(pr); - errcode = linsolve(&hyd->solver,net->Njuncs); - - /* Take action depending on error code */ - if (errcode < 0) { - break; /* Memory allocation problem */ - } - if (errcode > 0) { /* Ill-conditioning problem */ - /* If control valve causing problem, fix its status & continue, */ - /* otherwise end the iterations with no solution. */ - if (badvalve(pr,sol->Order[errcode])) { - continue; - } - else break; - } - - /* Update current solution. */ - /* (Row[i] = row of solution matrix corresponding to node i). */ - for (i=1; i<=net->Njuncs; i++) { - hyd->NodeHead[i] = sol->F[sol->Row[i]]; /* Update heads */ - } - newerr = newflows(pr); /* Update flows */ - *relerr = newerr; - - /* Write convergence error to status report if called for */ - if (rep->Statflag == FULL) { - writerelerr(pr, *iter,*relerr); - } - - /* Apply solution damping & check for change in valve status */ - hyd->RelaxFactor = 1.0; - valveChange = FALSE; - if ( hyd->DampLimit > 0.0 ) { - if( *relerr <= hyd->DampLimit ) { - hyd->RelaxFactor = 0.6; - valveChange = valvestatus(pr); - } - } - else { - valveChange = valvestatus(pr); - } - - /* Check for convergence */ - if (*relerr <= hyd->Hacc) { - /* We have convergence. Quit if we are into extra iterations. */ - if (*iter > hyd->MaxIter) { - break; - } - - /* Quit if no status changes occur. */ - statChange = FALSE; - if (valveChange) { - statChange = TRUE; - } - if (linkstatus(pr)) { - statChange = TRUE; - } - if (pswitch(pr)) { - statChange = TRUE; - } - if (!statChange) { - break; - } - - /* We have a status change so continue the iterations */ - nextcheck = *iter + hyd->CheckFreq; - } - - /* No convergence yet. See if its time for a periodic status */ - /* check on pumps, CV's, and pipes connected to tanks. */ - else if (*iter <= hyd->MaxCheck && *iter == nextcheck) - { - linkstatus(pr); - nextcheck += hyd->CheckFreq; - } - (*iter)++; - } - - /* Iterations ended. Report any errors. */ - if (errcode < 0) errcode = 101; /* Memory allocation error */ - else if (errcode > 0) - { - writehyderr(pr, sol->Order[errcode]); /* Ill-conditioned eqns. error */ - errcode = 110; - } - - /* Add any emitter flows to junction demands */ - for (i=1; i<=net->Njuncs; i++) { - hyd->NodeDemand[i] += hyd->EmitterFlows[i]; - } - return(errcode); -} /* End of netsolve */ - - -int badvalve(EN_Project *pr, int n) -/* - **----------------------------------------------------------------- -** Input: n = node index -** Output: returns 1 if node n belongs to an active control valve, -** 0 otherwise -** Purpose: determines if a node belongs to an active control valve -** whose setting causes an inconsistent set of eqns. If so, -** the valve status is fixed open and a warning condition -** is generated. -**----------------------------------------------------------------- -*/ -{ - int i,k,n1,n2; - Slink *link; - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - report_options_t *rep = &pr->report; - time_options_t *time = &pr->time_options; - - for (i=1; i <= net->Nvalves; i++) { - k = net->Valve[i].Link; - link = &net->Link[k]; - n1 = link->N1; - n2 = link->N2; - if (n == n1 || n == n2) { - EN_LinkType t = link->Type; - if (t == EN_PRV || t == EN_PSV || t == EN_FCV) { - if (hyd->LinkStatus[k] == ACTIVE) { - if (rep->Statflag == FULL) { - sprintf(pr->Msg,FMT61,clocktime(rep->Atime,time->Htime),link->ID); - writeline(pr, pr->Msg); - } - if (link->Type == EN_FCV) { - hyd->LinkStatus[k] = XFCV; - } - else { - hyd->LinkStatus[k] = XPRESSURE; - } - return(1); - } - } - return(0); - } - } - return(0); -} - - -int valvestatus(EN_Project *pr) -/* -**----------------------------------------------------------------- -** Input: none -** Output: returns 1 if any pressure or flow control valve -** changes status, 0 otherwise -** Purpose: updates status for PRVs & PSVs whose status -** is not fixed to OPEN/CLOSED -**----------------------------------------------------------------- -*/ -{ - int change = FALSE, /* Status change flag */ - i,k, /* Valve & link indexes */ - n1,n2; /* Start & end nodes */ - StatType status; /* Valve status settings */ - double hset; /* Valve head setting */ - Slink *link; - - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - report_options_t *rep = &pr->report; - - for (i=1; i <= net->Nvalves; i++) /* Examine each valve */ - { - k = net->Valve[i].Link; /* Link index of valve */ - link = &net->Link[k]; - if (hyd->LinkSetting[k] == MISSING) { - continue; /* Valve status fixed */ - } - n1 = link->N1; /* Start & end nodes */ - n2 = link->N2; - status = hyd->LinkStatus[k]; /* Save current status */ - -// if (s != CLOSED /* No change if flow is */ -// && ABS(hyd->LinkFlows[k]) < hyd->Qtol) continue; /* negligible. */ - - switch (link->Type) /* Evaluate new status: */ - { - case EN_PRV: - hset = net->Node[n2].El + hyd->LinkSetting[k]; - hyd->LinkStatus[k] = prvstatus(pr,k,status,hset,hyd->NodeHead[n1],hyd->NodeHead[n2]); - break; - case EN_PSV: - hset = net->Node[n1].El + hyd->LinkSetting[k]; - hyd->LinkStatus[k] = psvstatus(pr,k,status,hset,hyd->NodeHead[n1],hyd->NodeHead[n2]); - break; - -//// FCV status checks moved back into the linkstatus() function //// -// case FCV: S[k] = fcvstatus(k,s,hyd->NodeHead[n1],hyd->NodeHead[n2]); -// break; - - default: - continue; - } - -/*** Updated 9/7/00 ***/ - /* Do not reset flow in valve if its status changes. */ - /* This strategy improves convergence. */ - - /* Check for status change */ - if (status != hyd->LinkStatus[k]) - { - if (rep->Statflag == FULL) { - writestatchange(pr, k,status,hyd->LinkStatus[k]); - } - change = TRUE; - } - } - return(change); -} /* End of valvestatus() */ - - -/*** Updated 9/7/00 ***/ -int linkstatus(EN_Project *pr) -/* -**-------------------------------------------------------------- -** Input: none -** Output: returns 1 if any link changes status, 0 otherwise -** Purpose: determines new status for pumps, CVs, FCVs & pipes -** to tanks. -**-------------------------------------------------------------- -*/ -{ - int change = FALSE, /* Status change flag */ - k, /* Link index */ - n1, /* Start node index */ - n2; /* End node index */ - double dh; /* Head difference */ - StatType status; /* Current status */ - - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - report_options_t *rep = &pr->report; - - /* Examine each Slink */ - for (k=1; k <= net->Nlinks; k++) - { - Slink *link = &net->Link[k]; - n1 = link->N1; - n2 = link->N2; - dh = hyd->NodeHead[n1] - hyd->NodeHead[n2]; - - /* Re-open temporarily closed links (status = XHEAD or TEMPCLOSED) */ - status = hyd->LinkStatus[k]; - if (status == XHEAD || status == TEMPCLOSED) { - hyd->LinkStatus[k] = OPEN; - } - - /* Check for status changes in CVs and pumps */ - if (link->Type == EN_CVPIPE) { - hyd->LinkStatus[k] = cvstatus(pr,hyd->LinkStatus[k],dh,hyd->LinkFlows[k]); - } - if (link->Type == EN_PUMP && hyd->LinkStatus[k] >= OPEN && hyd->LinkSetting[k] > 0.0) { - hyd->LinkStatus[k] = pumpstatus(pr,k,-dh); - } - - /* Check for status changes in non-fixed FCVs */ - if (link->Type == EN_FCV && hyd->LinkSetting[k] != MISSING) { // - hyd->LinkStatus[k] = fcvstatus(pr,k,status,hyd->NodeHead[n1],hyd->NodeHead[n2]); // - } - - /* Check for flow into (out of) full (empty) tanks */ - if (n1 > net->Njuncs || n2 > net->Njuncs) { - tankstatus(pr,k,n1,n2); - } - - /* Note change in link status; do not revise link flow */ - if (status != hyd->LinkStatus[k]) { - change = TRUE; - if (rep->Statflag == FULL) { - writestatchange(pr,k,status,hyd->LinkStatus[k]); - } - - //if (S[k] <= CLOSED) hyd->LinkFlows[k] = QZERO; - //else setlinkflow(k, dh); - } - } - return(change); -} /* End of linkstatus */ - - -StatType cvstatus(EN_Project *pr, StatType s, double dh, double q) -/* -**-------------------------------------------------- -** Input: s = current status -** dh = headloss -** q = flow -** Output: returns new link status -** Purpose: updates status of a check valve. -**-------------------------------------------------- -*/ -{ - hydraulics_t *hyd = &pr->hydraulics; - - /* Prevent reverse flow through CVs */ - if (ABS(dh) > hyd->Htol) - { - if (dh < -hyd->Htol) { - return(CLOSED); - } - else if (q < -hyd->Qtol) { - return(CLOSED); - } - else { - return(OPEN); - } - } - else - { - if (q < -hyd->Qtol) { - return(CLOSED); - } - else { - return(s); - } - } -} - - -StatType pumpstatus(EN_Project *pr, int k, double dh) -/* -**-------------------------------------------------- -** Input: k = link index -** dh = head gain -** Output: returns new pump status -** Purpose: updates status of an open pump. -**-------------------------------------------------- -*/ -{ - int p; - double hmax; - hydraulics_t *hyd = &pr->hydraulics; - EN_Network *net = &pr->network; - - /* Prevent reverse flow through pump */ - p = findpump(net,k); - if (net->Pump[p].Ptype == CONST_HP) { - hmax = BIG; - } - else { - hmax = SQR(hyd->LinkSetting[k]) * net->Pump[p].Hmax; - } - if (dh > hmax + hyd->Htol) { - return(XHEAD); - } - -/*** Flow higher than pump curve no longer results in a status change ***/ - /* Check if pump cannot deliver flow */ - //if (hyd->LinkFlows[k] > K[k]*Pump[p].Qmax + hyd->Qtol) return(XFLOW); - return(OPEN); -} - - -StatType prvstatus(EN_Project *pr, int k, StatType s, double hset, double h1, double h2) -/* -**----------------------------------------------------------- -** Input: k = link index -** s = current status -** hset = valve head setting -** h1 = head at upstream node -** h2 = head at downstream node -** Output: returns new valve status -** Purpose: updates status of a pressure reducing valve. -**----------------------------------------------------------- -*/ -{ - StatType status; /* New valve status */ - double hml; /* Minor headloss */ - hydraulics_t *hyd = &pr->hydraulics; - - double htol = hyd->Htol; - Slink *link = &pr->network.Link[k]; - - status = s; - if (hyd->LinkSetting[k] == MISSING) - return(status); /* Status fixed by user */ - hml = link->Km*SQR(hyd->LinkFlows[k]); /* Head loss when open */ - - /*** Status rules below have changed. ***/ - - switch (s) - { - case ACTIVE: - if (hyd->LinkFlows[k] < -hyd->Qtol) { - status = CLOSED; - } - else if (h1-hml < hset-htol) { - status = OPEN; - } - else { - status = ACTIVE; - } - break; - case OPEN: - if (hyd->LinkFlows[k] < -hyd->Qtol) { - status = CLOSED; - } - else if (h2 >= hset+htol) { - status = ACTIVE; - } - else { - status = OPEN; - } - break; - case CLOSED: - if ( h1 >= hset+htol && h2 < hset-htol) { - status = ACTIVE; - } - else if (h1 < hset-htol && h1 > h2+htol) { - status = OPEN; - } - else { - status = CLOSED; - } - break; - case XPRESSURE: - if (hyd->LinkFlows[k] < -hyd->Qtol) { - status = CLOSED; - } - break; - default: - break; - } - return(status); -} - - -StatType psvstatus(EN_Project *pr, int k, StatType s, double hset, double h1, double h2) -/* -**----------------------------------------------------------- -** Input: k = link index -** s = current status -** hset = valve head setting -** h1 = head at upstream node -** h2 = head at downstream node -** Output: returns new valve status -** Purpose: updates status of a pressure sustaining valve. -**----------------------------------------------------------- -*/ -{ - StatType status; /* New valve status */ - double hml; /* Minor headloss */ - hydraulics_t *hyd = &pr->hydraulics; - - double htol = hyd->Htol; - - Slink *link = &pr->network.Link[k]; - - status = s; - if (hyd->LinkSetting[k] == MISSING) { - return(status); /* Status fixed by user */ - } - hml = link->Km*SQR(hyd->LinkFlows[k]); /* Head loss when open */ - - /*** Status rules below have changed. ***/ - - switch (s) - { - case ACTIVE: - if (hyd->LinkFlows[k] < -hyd->Qtol) { - status = CLOSED; - } - else if (h2+hml > hset+htol) { - status = OPEN; - } - else { - status = ACTIVE; - } - break; - case OPEN: - if (hyd->LinkFlows[k] < -hyd->Qtol) { - status = CLOSED; - } - else if (h1 < hset-htol) { - status = ACTIVE; - } - else { - status = OPEN; - } - break; - case CLOSED: - if (h2 > hset+htol && h1 > h2+htol) { - status = OPEN; - } - else if (h1 >= hset+htol && h1 > h2+htol) { - status = ACTIVE; - } - else { - status = CLOSED; - } - break; - case XPRESSURE: - if (hyd->LinkFlows[k] < -hyd->Qtol) { - status = CLOSED; - } - break; - default: - break; - } - return(status); -} - - -StatType fcvstatus(EN_Project *pr, int k, StatType s, double h1, double h2) -/* -**----------------------------------------------------------- -** Input: k = link index -** s = current status -** h1 = head at upstream node -** h2 = head at downstream node -** Output: returns new valve status -** Purpose: updates status of a flow control valve. -** -** Valve status changes to XFCV if flow reversal. -** If current status is XFCV and current flow is -** above setting, then valve becomes active. -** If current status is XFCV, and current flow -** positive but still below valve setting, then -** status remains same. -**----------------------------------------------------------- -*/ -{ - StatType status; /* New valve status */ - hydraulics_t *hyd = &pr->hydraulics; - - status = s; - if (h1 - h2 < -hyd->Htol) { - status = XFCV; - } - else if ( hyd->LinkFlows[k] < -hyd->Qtol ) { - status = XFCV; - } - else if (s == XFCV && hyd->LinkFlows[k] >= hyd->LinkSetting[k]) { - status = ACTIVE; - } - return(status); -} - - -/*** Updated 9/7/00 ***/ -/*** Updated 11/19/01 ***/ -void tankstatus(EN_Project *pr, int k, int n1, int n2) -/* -**---------------------------------------------------------------- -** Input: k = link index -** n1 = start node of link -** n2 = end node of link -** Output: none -** Purpose: closes link flowing into full or out of empty tank -**---------------------------------------------------------------- -*/ -{ - int i,n; - double h,q; - Stank *tank; - - hydraulics_t *hyd = &pr->hydraulics; - EN_Network *net = &pr->network; - Slink *link = &net->Link[k]; - - /* Make node n1 be the tank */ - q = hyd->LinkFlows[k]; - i = n1 - net->Njuncs; - if (i <= 0) { - i = n2 - net->Njuncs; - if (i <= 0) { - return; - } - n = n1; - n1 = n2; - n2 = n; - q = -q; - } - h = hyd->NodeHead[n1] - hyd->NodeHead[n2]; - tank = &net->Tank[i]; - /* Skip reservoirs & closed links */ - if (tank->A == 0.0 || hyd->LinkStatus[k] <= CLOSED) { - return; - } - - /* If tank full, then prevent flow into it */ - if (hyd->NodeHead[n1] >= tank->Hmax - hyd->Htol) - { - - /* Case 1: Link is a pump discharging into tank */ - if ( link->Type == EN_PUMP ) { - if (link->N2 == n1) - hyd->LinkStatus[k] = TEMPCLOSED; - } - - /* Case 2: Downstream head > tank head */ - /* (i.e., an open outflow check valve would close) */ - else if (cvstatus(pr, OPEN, h, q) == CLOSED) { - hyd->LinkStatus[k] = TEMPCLOSED; - } - } - - /* If tank empty, then prevent flow out of it */ - if (hyd->NodeHead[n1] <= tank->Hmin + hyd->Htol) { - - /* Case 1: Link is a pump discharging from tank */ - if ( link->Type == EN_PUMP) { - if (link->N1 == n1) { - hyd->LinkStatus[k] = TEMPCLOSED; - } - } - - /* Case 2: Tank head > downstream head */ - /* (i.e., a closed outflow check valve would open) */ - else if (cvstatus(pr, CLOSED, h, q) == OPEN) { - hyd->LinkStatus[k] = TEMPCLOSED; - } - } -} /* End of tankstatus */ - - -int pswitch(EN_Project *pr) -/* -**-------------------------------------------------------------- -** Input: none -** Output: returns 1 if status of any link changes, 0 if not -** Purpose: adjusts settings of links controlled by junction -** pressures after a hydraulic solution is found -**-------------------------------------------------------------- -*/ -{ - int i, /* Control statement index */ - k, /* Link being controlled */ - n, /* Node controlling Slink */ - reset, /* Flag on control conditions */ - change, /* Flag for status or setting change */ - anychange = 0; /* Flag for 1 or more changes */ - char s; /* Current link status */ - - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - report_options_t *rep = &pr->report; - - /* Check each control statement */ - for (i=1; i <= net->Ncontrols; i++) - { - reset = 0; - if ( (k = net->Control[i].Link) <= 0) { - continue; - } - - /* Determine if control based on a junction, not a tank */ - if ( (n = net->Control[i].Node) > 0 && n <= net->Njuncs) { - /* Determine if control conditions are satisfied */ - if (net->Control[i].Type == LOWLEVEL - && hyd->NodeHead[n] <= net->Control[i].Grade + hyd->Htol ) { - reset = 1; - } - if (net->Control[i].Type == HILEVEL - && hyd->NodeHead[n] >= net->Control[i].Grade - hyd->Htol ) { - reset = 1; - } - } - - /* Determine if control forces a status or setting change */ - if (reset == 1) - { - Slink *link = &net->Link[k]; - change = 0; - s = hyd->LinkStatus[k]; - if (link->Type == EN_PIPE) { - if (s != net->Control[i].Status) { - change = 1; - } - } - if (link->Type == EN_PUMP) { - if (hyd->LinkSetting[k] != net->Control[i].Setting) { - change = 1; - } - } - if (link->Type >= EN_PRV) { - if (hyd->LinkSetting[k] != net->Control[i].Setting) { - change = 1; - } - else if (hyd->LinkSetting[k] == MISSING && s != net->Control[i].Status) { - change = 1; - } - } - - /* If a change occurs, update status & setting */ - if (change) { - hyd->LinkStatus[k] = net->Control[i].Status; - if (link->Type > EN_PIPE) { - hyd->LinkSetting[k] = net->Control[i].Setting; - } - if (rep->Statflag == FULL) { - writestatchange(pr, k,s,hyd->LinkStatus[k]); - } - - /* Re-set flow if status has changed */ -// if (S[k] != s) initlinkflow(k, S[k], K[k]); - anychange = 1; - } - } - } - return(anychange); -} /* End of pswitch */ - - -double newflows(EN_Project *pr) -/* -**---------------------------------------------------------------- -** Input: none -** Output: returns solution convergence error -** Purpose: updates link flows after new nodal heads computed -**---------------------------------------------------------------- -*/ -{ - double dh, /* Link head loss */ - dq; /* Link flow change */ - double dqsum, /* Network flow change */ - qsum; /* Network total flow */ - int k, n, n1, n2; - - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - solver_t *s = &hyd->solver; - - /* Initialize net inflows (i.e., demands) at tanks */ - for (n = net->Njuncs + 1; n <= net->Nnodes; n++) { - hyd->NodeDemand[n] = 0.0; - } - - /* Initialize sum of flows & corrections */ - qsum = 0.0; - dqsum = 0.0; - - /* Update flows in all links */ - for (k=1; k<=net->Nlinks; k++) - { - Slink *link = &net->Link[k]; - /* - ** Apply flow update formula: - ** dq = Y - P*(new head loss) - ** P = 1/(dh/dq) - ** Y = P*(head loss based on current flow) - ** where P & Y were computed in newcoeffs(). - */ - - n1 = link->N1; - n2 = link->N2; - dh = hyd->NodeHead[n1] - hyd->NodeHead[n2]; - dq = s->Y[k] - s->P[k] * dh; - - /* Adjust flow change by the relaxation factor */ - dq *= hyd->RelaxFactor; - - /* Prevent flow in constant HP pumps from going negative */ - if (link->Type == EN_PUMP) { - n = findpump(net,k); - if (net->Pump[n].Ptype == CONST_HP && dq > hyd->LinkFlows[k]) { - dq = hyd->LinkFlows[k]/2.0; - } - } - hyd->LinkFlows[k] -= dq; - - /* Update sum of absolute flows & flow corrections */ - qsum += ABS(hyd->LinkFlows[k]); - dqsum += ABS(dq); - - /* Update net flows to tanks */ - if ( hyd->LinkStatus[k] > CLOSED ) - { - if (n1 > net->Njuncs) { - hyd->NodeDemand[n1] -= hyd->LinkFlows[k]; - } - if (n2 > net->Njuncs) { - hyd->NodeDemand[n2] += hyd->LinkFlows[k]; - } - } - - } - - /* Update emitter flows */ - for (k=1; k<=net->Njuncs; k++) - { - if (net->Node[k].Ke == 0.0) { - continue; - } - dq = emitflowchange(pr,k); - hyd->EmitterFlows[k] -= dq; - qsum += ABS(hyd->EmitterFlows[k]); - dqsum += ABS(dq); - } - - /* Return ratio of total flow corrections to total flow */ - if (qsum > hyd->Hacc) { - return(dqsum/qsum); - } - else { - return(dqsum); - } -} /* End of newflows */ - - -void newcoeffs(EN_Project *pr) -/* -**-------------------------------------------------------------- -** Input: none -** Output: none -** Purpose: computes coefficients of linearized network eqns. -**-------------------------------------------------------------- -*/ -{ - hydraulics_t *hyd = &pr->hydraulics; - solver_t *s = &hyd->solver; - EN_Network *n = &pr->network; - - memset(s->Aii,0,(n->Nnodes+1)*sizeof(double)); /* Reset coeffs. to 0 */ - memset(s->Aij,0,(hyd->Ncoeffs+1)*sizeof(double)); - memset(s->F,0,(n->Nnodes+1)*sizeof(double)); - memset(hyd->X_tmp,0,(n->Nnodes+1)*sizeof(double)); - memset(s->P,0,(n->Nlinks+1)*sizeof(double)); - memset(s->Y,0,(n->Nlinks+1)*sizeof(double)); - linkcoeffs(pr); /* Compute link coeffs. */ - emittercoeffs(pr); /* Compute emitter coeffs.*/ - nodecoeffs(pr); /* Compute node coeffs. */ - valvecoeffs(pr); /* Compute valve coeffs. */ -} /* End of newcoeffs */ - - -void linkcoeffs(EN_Project *pr) -/* -**-------------------------------------------------------------- -** Input: none -** Output: none -** Purpose: computes solution matrix coefficients for links -**-------------------------------------------------------------- -*/ -{ - int k,n1,n2; - - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - solver_t *s = &hyd->solver; - EN_Network *n = &pr->network; - - /* Examine each link of network */ - for (k=1; k<=net->Nlinks; k++) - { - Slink *link = &net->Link[k]; - n1 = link->N1; /* Start node of Slink */ - n2 = link->N2; /* End node of link */ - - /* Compute P[k] = 1 / (dh/dQ) and Y[k] = h * P[k] */ - /* for each link k (where h = link head loss). */ - /* FCVs, PRVs, and PSVs with non-fixed status */ - /* are analyzed later. */ - - switch (link->Type) - { - case EN_CVPIPE: - case EN_PIPE: - pipecoeff(pr, k); - break; - case EN_PUMP: - pumpcoeff(pr, k); - break; - case EN_PBV: - pbvcoeff(pr, k); - break; - case EN_TCV: - tcvcoeff(pr, k); - break; - case EN_GPV: - gpvcoeff(pr, k); - break; - case EN_FCV: - case EN_PRV: - case EN_PSV: /* If valve status fixed then treat as pipe, otherwise ignore the valve for now. */ - if (hyd->LinkSetting[k] == MISSING) { - valvecoeff(pr, k); //pipecoeff(k); - } - else { - continue; - } - break; - default: - continue; - } - - /* Update net nodal inflows (X), solution matrix (A) and RHS array (F) */ - /* (Use covention that flow out of node is (-), flow into node is (+)) */ - hyd->X_tmp[n1] -= hyd->LinkFlows[k]; - hyd->X_tmp[n2] += hyd->LinkFlows[k]; - s->Aij[s->Ndx[k]] -= s->P[k]; /* Off-diagonal coeff. */ - if (n1 <= n->Njuncs) /* Node n1 is junction */ - { - s->Aii[s->Row[n1]] += s->P[k]; /* Diagonal coeff. */ - s->F[s->Row[n1]] += s->Y[k]; /* RHS coeff. */ - } - else { - s->F[s->Row[n2]] += (s->P[k]*hyd->NodeHead[n1]); /* Node n1 is a tank */ - } - if (n2 <= n->Njuncs) { /* Node n2 is junction */ - s->Aii[s->Row[n2]] += s->P[k]; /* Diagonal coeff. */ - s->F[s->Row[n2]] -= s->Y[k]; /* RHS coeff. */ - } - else { - s->F[s->Row[n1]] += (s->P[k] * hyd->NodeHead[n2]); /* Node n2 is a tank */ - } - } -} /* End of linkcoeffs */ - - -void nodecoeffs(EN_Project *pr) -/* -**---------------------------------------------------------------- -** Input: none -** Output: none -** Purpose: completes calculation of nodal flow imbalance (X) -** & flow correction (F) arrays -**---------------------------------------------------------------- -*/ -{ - int i; - hydraulics_t *hyd = &pr->hydraulics; - solver_t *s = &hyd->solver; - EN_Network *n = &pr->network; - - /* For junction nodes, subtract demand flow from net */ - /* flow imbalance & add imbalance to RHS array F. */ - for (i=1; i <= n->Njuncs; i++) - { - hyd->X_tmp[i] -= hyd->NodeDemand[i]; - s->F[s->Row[i]] += hyd->X_tmp[i]; - } -} /* End of nodecoeffs */ - - -void valvecoeffs(EN_Project *pr) -/* -**-------------------------------------------------------------- -** Input: none -** Output: none -** Purpose: computes matrix coeffs. for PRVs, PSVs & FCVs -** whose status is not fixed to OPEN/CLOSED -**-------------------------------------------------------------- -*/ -{ - int i,k,n1,n2; - - hydraulics_t *hyd = &pr->hydraulics; - EN_Network *n = &pr->network; - Slink *link; - Svalve *valve; - - for (i=1; i<=n->Nvalves; i++) /* Examine each valve */ - { - valve = &n->Valve[i]; - k = valve->Link; /* Link index of valve */ - link = &n->Link[k]; - if (hyd->LinkSetting[k] == MISSING) { - continue; /* Valve status fixed */ - } - n1 = link->N1; /* Start & end nodes */ - n2 = link->N2; - switch (link->Type) /* Call valve-specific */ - { /* function */ - case EN_PRV: - prvcoeff(pr, k,n1,n2); - break; - case EN_PSV: - psvcoeff(pr, k,n1,n2); - break; - case EN_FCV: - fcvcoeff(pr, k,n1,n2); - break; - default: continue; - } - } -} /* End of valvecoeffs */ - - -void emittercoeffs(EN_Project *pr) -/* -**-------------------------------------------------------------- -** Input: none -** Output: none -** Purpose: computes matrix coeffs. for emitters -** -** Note: Emitters consist of a fictitious pipe connected to -** a fictitious reservoir whose elevation equals that -** of the junction. The headloss through this pipe is -** Ke*(Flow)^hyd->Qexp, where Ke = emitter headloss coeff. -**-------------------------------------------------------------- -*/ -{ - int i; - double ke; - double p; - double q; - double y; - double z; - - hydraulics_t *hyd = &pr->hydraulics; - solver_t *s = &hyd->solver; - EN_Network *n = &pr->network; - - for (i=1; i <= n->Njuncs; i++) - { - Snode *node = &n->Node[i]; - if (node->Ke == 0.0) { - continue; - } - ke = MAX(CSMALL, node->Ke); - q = hyd->EmitterFlows[i]; - z = ke * pow(ABS(q),hyd->Qexp); - p = hyd->Qexp * z / ABS(q); - if (p < hyd->RQtol) { - p = 1.0 / hyd->RQtol; - } - else { - p = 1.0 / p; - } - y = SGN(q)*z*p; - s->Aii[s->Row[i]] += p; - s->F[s->Row[i]] += y + p * node->El; - hyd->X_tmp[i] -= q; - } -} - - -double emitflowchange(EN_Project *pr, int i) -/* -**-------------------------------------------------------------- -** Input: i = node index -** Output: returns change in flow at an emitter node -** Purpose: computes flow change at an emitter node -**-------------------------------------------------------------- -*/ -{ - double ke, p; - hydraulics_t *hyd = &pr->hydraulics; - EN_Network *n = &pr->network; - Snode *node = &n->Node[i]; - - ke = MAX(CSMALL, node->Ke); - p = hyd->Qexp * ke * pow(ABS(hyd->EmitterFlows[i]),(hyd->Qexp - 1.0)); - if (p < hyd->RQtol) { - p = 1/hyd->RQtol; - } - else { - p = 1.0/p; - } - return(hyd->EmitterFlows[i] / hyd->Qexp - p * (hyd->NodeHead[i] - node->El)); -} - - -void pipecoeff(EN_Project *pr, int k) -/* -**-------------------------------------------------------------- -** Input: k = link index -** Output: none -** Purpose: computes P & Y coefficients for pipe k -** -** P = inverse head loss gradient = 1/(dh/dQ) -** Y = flow correction term = h*P -**-------------------------------------------------------------- -*/ -{ - double hpipe, /* Normal head loss */ - hml, /* Minor head loss */ - ml, /* Minor loss coeff. */ - p, /* q*(dh/dq) */ - q, /* Abs. value of flow */ - r, /* Resistance coeff. */ - r1, /* Total resistance factor */ - f, /* D-W friction factor */ - dfdq; /* Derivative of fric. fact. */ - - hydraulics_t *hyd = &pr->hydraulics; - solver_t *s = &hyd->solver; - Slink *link = &pr->network.Link[k]; - - /* For closed pipe use headloss formula: h = CBIG*q */ - if (hyd->LinkStatus[k] <= CLOSED) - { - s->P[k] = 1.0/CBIG; - s->Y[k] = hyd->LinkFlows[k]; - return; - } - - /* Evaluate headloss coefficients */ - q = ABS(hyd->LinkFlows[k]); /* Absolute flow */ - ml = link->Km; /* Minor loss coeff. */ - r = link->R; /* Resistance coeff. */ - f = 1.0; /* D-W friction factor */ - if (hyd->Formflag == DW) f = DWcoeff(pr,k,&dfdq); - r1 = f*r+ml; - - /* Use large P coefficient for small flow resistance product */ - if (r1*q < hyd->RQtol) - { - s->P[k] = 1.0/hyd->RQtol; - s->Y[k] = hyd->LinkFlows[k]/hyd->Hexp; - return; - } - - /* Compute P and Y coefficients */ - if (hyd->Formflag == DW) /* D-W eqn. */ - { - hpipe = r1*SQR(q); /* Total head loss */ - p = 2.0*r1*q; /* |dh/dQ| */ - /* + dfdq*r*q*q;*/ /* Ignore df/dQ term */ - p = 1.0/p; - s->P[k] = p; - s->Y[k] = SGN(hyd->LinkFlows[k])*hpipe*p; - } - else /* H-W or C-M eqn. */ - { - hpipe = r*pow(q,hyd->Hexp); /* Friction head loss */ - p = hyd->Hexp*hpipe; /* Q*dh(friction)/dQ */ - if (ml > 0.0) - { - hml = ml*q*q; /* Minor head loss */ - p += 2.0*hml; /* Q*dh(Total)/dQ */ - } - else hml = 0.0; - p = hyd->LinkFlows[k]/p; /* 1 / (dh/dQ) */ - s->P[k] = ABS(p); - s->Y[k] = p*(hpipe + hml); - } -} /* End of pipecoeff */ - - -double DWcoeff(EN_Project *pr, int k, double *dfdq) -/* -**-------------------------------------------------------------- -** Input: k = link index -** Output: returns Darcy-Weisbach friction factor -** Purpose: computes Darcy-Weisbach friction factor -** -** Uses interpolating polynomials developed by -** E. Dunlop for transition flow from 2000 < Re < 4000. -** -** df/dq term is ignored as it slows convergence rate. -**-------------------------------------------------------------- -*/ -{ - double q, /* Abs. value of flow */ - f; /* Friction factor */ - double x1,x2,x3,x4, - y1,y2,y3, - fa,fb,r; - double s,w; - - hydraulics_t *hyd = &pr->hydraulics; - Slink *link = &pr->network.Link[k]; - - *dfdq = 0.0; - if (link->Type > EN_PIPE) - return(1.0); /* Only apply to pipes */ - q = ABS(hyd->LinkFlows[k]); - s = hyd->Viscos * link->Diam; - w = q/s; /* w = Re(Pi/4) */ - if (w >= A1) /* Re >= 4000; Colebrook Formula */ - { - y1 = A8/pow(w,0.9); - y2 = link->Kc/(3.7*link->Diam) + y1; - y3 = A9*log(y2); - f = 1.0/SQR(y3); - /* *dfdq = (2.0+AA*y1/(y2*y3))*f; */ /* df/dq */ - } - else if (w > A2) /* Re > 2000; Interpolation formula */ - { - y2 = link->Kc/(3.7*link->Diam) + AB; - y3 = A9*log(y2); - fa = 1.0/SQR(y3); - fb = (2.0+AC/(y2*y3))*fa; - r = w/A2; - x1 = 7.0*fa - fb; - x2 = 0.128 - 17.0*fa + 2.5*fb; - x3 = -0.128 + 13.0*fa - (fb+fb); - x4 = r*(0.032 - 3.0*fa + 0.5*fb); - f = x1 + r*(x2 + r*(x3+x4)); - /* *dfdq = (x1 + x1 + r*(3.0*x2 + r*(4.0*x3 + 5.0*x4))); */ - } - else if (w > A4) /* Laminar flow: Hagen-Poiseuille Formula */ - { - f = A3*s/q; - /* *dfdq = A3*s; */ - } - else - { - f = 8.0; - *dfdq = 0.0; - } - return(f); -} /* End of DWcoeff */ - - -/*** Updated 10/25/00 ***/ -/*** Updated 12/29/00 ***/ -void pumpcoeff(EN_Project *pr, int k) -/* -**-------------------------------------------------------------- -** Input: k = link index -** Output: none -** Purpose: computes P & Y coeffs. for pump in link k -**-------------------------------------------------------------- -*/ -{ - int p; /* Pump index */ - double h0, /* Shutoff head */ - q, /* Abs. value of flow */ - r, /* Flow resistance coeff. */ - n; /* Flow exponent coeff. */ - hydraulics_t *hyd; - solver_t *s; - double setting; - Spump *pump; - - hyd = &pr->hydraulics; - s = &hyd->solver; - - setting = hyd->LinkSetting[k]; - - /* Use high resistance pipe if pump closed or cannot deliver head */ - if (hyd->LinkStatus[k] <= CLOSED || setting == 0.0) { - s->P[k] = 1.0/CBIG; - s->Y[k] = hyd->LinkFlows[k]; - return; - } - - q = ABS(hyd->LinkFlows[k]); - q = MAX(q,TINY); - p = findpump(&pr->network,k); - - pump = &pr->network.Pump[p]; - /* Get pump curve coefficients for custom pump curve. */ - if (pump->Ptype == CUSTOM) - { - /* Find intercept (h0) & slope (r) of pump curve */ - /* line segment which contains speed-adjusted flow. */ - curvecoeff(pr,pump->Hcurve, q/setting, &h0, &r); - - /* Determine head loss coefficients. */ - pump->H0 = -h0; - pump->R = -r; - pump->N = 1.0; - } - - /* Adjust head loss coefficients for pump speed. */ - h0 = SQR(setting) * pump->H0; - n = pump->N; - r = pump->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) */ - s->P[k] = 1.0/MAX(r,hyd->RQtol); - s->Y[k] = hyd->LinkFlows[k]/n + s->P[k]*h0; -} /* End of pumpcoeff */ - - -/*** Updated 10/25/00 ***/ -/*** Updated 12/29/00 ***/ -void curvecoeff(EN_Project *p, int i, double q, double *h0, double *r) -/* -**------------------------------------------------------------------- -** Input: i = curve index -** q = flow rate -** Output: *h0 = head at zero flow (y-intercept) -** *r = dHead/dFlow (slope) -** Purpose: computes intercept and slope of head v. flow curve -** at current flow. -**------------------------------------------------------------------- -*/ -{ - int k1, k2, npts; - double *x, *y; - Scurve *curve; - - /* Remember that curve is stored in untransformed units */ - q *= p->Ucf[FLOW]; - curve = &p->network.Curve[i]; - x = curve->X; /* x = flow */ - y = curve->Y; /* y = head */ - npts = curve->Npts; - - /* Find linear segment of curve that brackets flow q */ - k2 = 0; - while (k2 < npts && x[k2] < q) - k2++; - - if (k2 == 0) - k2++; - - else if (k2 == npts) - k2--; - - k1 = k2 - 1; - - /* Compute slope and intercept of this segment */ - *r = (y[k2]-y[k1]) / (x[k2]-x[k1]); - *h0 = y[k1] - (*r)*x[k1]; - - /* Convert units */ - *h0 = (*h0) / p->Ucf[HEAD]; - *r = (*r)*p->Ucf[FLOW]/p->Ucf[HEAD]; -} /* End of curvecoeff */ - - -void gpvcoeff(EN_Project *p, int k) -/* -**-------------------------------------------------------------- -** Input: k = link index -** Output: none -** Purpose: computes P & Y coeffs. for general purpose valve -**-------------------------------------------------------------- -*/ -{ - double h0, /* Headloss curve intercept */ - q, /* Abs. value of flow */ - r; /* Flow resistance coeff. */ - - hydraulics_t *hyd = &p->hydraulics; - solver_t *s = &hyd->solver; - -/*** Updated 9/7/00 ***/ - /* Treat as a pipe if valve closed */ - if (hyd->LinkStatus[k] == CLOSED) { - valvecoeff(p,k); //pipecoeff(k); - } - /* Otherwise utilize headloss curve */ - /* whose index is stored in K */ - else { - /* Find slope & intercept of headloss curve. */ - q = ABS(hyd->LinkFlows[k]); - q = MAX(q,TINY); - -/*** Updated 10/25/00 ***/ -/*** Updated 12/29/00 ***/ - curvecoeff(p,(int)ROUND(hyd->LinkSetting[k]),q,&h0,&r); - - /* Compute inverse headloss gradient (P) */ - /* and flow correction factor (Y). */ - s->P[k] = 1.0 / MAX(r,hyd->RQtol); - s->Y[k] = s->P[k]*(h0 + r*q) * SGN(hyd->LinkFlows[k]); - } -} - - -void pbvcoeff(EN_Project *p, int k) -/* -**-------------------------------------------------------------- -** Input: k = link index -** Output: none -** Purpose: computes P & Y coeffs. for pressure breaker valve -**-------------------------------------------------------------- -*/ -{ - - Slink *link = &p->network.Link[k]; - hydraulics_t *hyd = &p->hydraulics; - solver_t *s = &hyd->solver; - - /* If valve fixed OPEN or CLOSED then treat as a pipe */ - if (hyd->LinkSetting[k] == MISSING || hyd->LinkSetting[k] == 0.0) { - valvecoeff(p,k); //pipecoeff(k); - } - - /* If valve is active */ - else { - /* Treat as a pipe if minor loss > valve setting */ - if (link->Km * SQR(hyd->LinkFlows[k]) > hyd->LinkSetting[k]) { - valvecoeff(p,k); //pipecoeff(k); - } - /* Otherwise force headloss across valve to be equal to setting */ - else { - s->P[k] = CBIG; - s->Y[k] = hyd->LinkSetting[k] * CBIG; - } - } -} /* End of pbvcoeff */ - - -void tcvcoeff(EN_Project *p, int k) -/* -**-------------------------------------------------------------- -** Input: k = link index -** Output: none -** Purpose: computes P & Y coeffs. for throttle control valve -**-------------------------------------------------------------- -*/ -{ - double km; - Slink *link = &p->network.Link[k]; - hydraulics_t *hyd = &p->hydraulics; - - /* Save original loss coeff. for open valve */ - km = link->Km; - - /* If valve not fixed OPEN or CLOSED, compute its loss coeff. */ - if (hyd->LinkSetting[k] != MISSING) { - link->Km = 0.02517 * hyd->LinkSetting[k] / (SQR(link->Diam)*SQR(link->Diam)); - } - - /* Then apply usual pipe formulas */ - valvecoeff(p,k); //pipecoeff(k); - - /* Restore original loss coeff. */ - link->Km = km; -} /* End of tcvcoeff */ - - -void prvcoeff(EN_Project *p, int k, int n1, int n2) -/* -**-------------------------------------------------------------- -** Input: k = link index -** n1 = upstream node of valve -** n2 = downstream node of valve -** Output: none -** Purpose: computes solution matrix coeffs. for pressure -** reducing valves -**-------------------------------------------------------------- -*/ -{ - int i,j; /* Rows of solution matrix */ - double hset; /* Valve head setting */ - hydraulics_t *hyd = &p->hydraulics; - solver_t *s = &hyd->solver; - - i = s->Row[n1]; /* Matrix rows of nodes */ - j = s->Row[n2]; - hset = p->network.Node[n2].El + hyd->LinkSetting[k]; /* Valve setting */ - - if (hyd->LinkStatus[k] == ACTIVE) - { - /* - Set coeffs. to force head at downstream - node equal to valve setting & force flow (when updated in - newflows()) equal to flow imbalance at downstream node. - */ - s->P[k] = 0.0; - s->Y[k] = hyd->LinkFlows[k] + hyd->X_tmp[n2]; /* Force flow balance */ - s->F[j] += (hset * CBIG); /* Force head = hset */ - s->Aii[j] += CBIG; /* at downstream node */ - if (hyd->X_tmp[n2] < 0.0) { - s->F[i] += hyd->X_tmp[n2]; - } - return; - } - - /* - For OPEN, CLOSED, or XPRESSURE valve - compute matrix coeffs. using the valvecoeff() function. - */ - valvecoeff(p,k); /*pipecoeff(k);*/ - s->Aij[s->Ndx[k]] -= s->P[k]; - s->Aii[i] += s->P[k]; - s->Aii[j] += s->P[k]; - s->F[i] += (s->Y[k] - hyd->LinkFlows[k]); - s->F[j] -= (s->Y[k] - hyd->LinkFlows[k]); -} /* End of prvcoeff */ - - -void psvcoeff(EN_Project *p, int k, int n1, int n2) -/* -**-------------------------------------------------------------- -** Input: k = link index -** n1 = upstream node of valve -** n2 = downstream node of valve -** Output: none -** Purpose: computes solution matrix coeffs. for pressure -** sustaining valve -**-------------------------------------------------------------- -*/ -{ - int i,j; /* Rows of solution matrix */ - double hset; /* Valve head setting */ - - hydraulics_t *hyd = &p->hydraulics; - solver_t *s = &hyd->solver; - - i = s->Row[n1]; /* Matrix rows of nodes */ - j = s->Row[n2]; - hset = p->network.Node[n1].El + hyd->LinkSetting[k]; /* Valve setting */ - - if (hyd->LinkStatus[k] == ACTIVE) - { - /* - Set coeffs. to force head at upstream - node equal to valve setting & force flow (when updated in - newflows()) equal to flow imbalance at upstream node. - */ - s->P[k] = 0.0; - s->Y[k] = hyd->LinkFlows[k] - hyd->X_tmp[n1]; /* Force flow balance */ - s->F[i] += (hset*CBIG); /* Force head = hset */ - s->Aii[i] += CBIG; /* at upstream node */ - if (hyd->X_tmp[n1] > 0.0) { - s->F[j] += hyd->X_tmp[n1]; - } - return; - } - - /* - For OPEN, CLOSED, or XPRESSURE valve - compute matrix coeffs. using the valvecoeff() function. - */ - valvecoeff(p,k); /*pipecoeff(k);*/ - s->Aij[s->Ndx[k]] -= s->P[k]; - s->Aii[i] += s->P[k]; - s->Aii[j] += s->P[k]; - s->F[i] += (s->Y[k] - hyd->LinkFlows[k]); - s->F[j] -= (s->Y[k] - hyd->LinkFlows[k]); -} /* End of psvcoeff */ - - -void fcvcoeff(EN_Project *p, int k, int n1, int n2) -/* -**-------------------------------------------------------------- -** Input: k = link index -** n1 = upstream node of valve -** n2 = downstream node of valve -** Output: none -** Purpose: computes solution matrix coeffs. for flow control -** valve -**-------------------------------------------------------------- -*/ -{ - int i,j; /* Rows in solution matrix */ - double q; /* Valve flow setting */ - - hydraulics_t *hyd = &p->hydraulics; - solver_t *s = &hyd->solver; - - q = hyd->LinkSetting[k]; - i = hyd->solver.Row[n1]; - j = hyd->solver.Row[n2]; - - /* - If valve active, break network at valve and treat - flow setting as external demand at upstream node - and external supply at downstream node. - */ - if (hyd->LinkStatus[k] == ACTIVE) - { - hyd->X_tmp[n1] -= q; - s->F[i] -= q; - hyd->X_tmp[n2] += q; - s->F[j] += q; - /*P[k] = 0.0;*/ - s->P[k] = 1.0/CBIG; - s->Aij[s->Ndx[k]] -= s->P[k]; - s->Aii[i] += s->P[k]; - s->Aii[j] += s->P[k]; - s->Y[k] = hyd->LinkFlows[k] - q; - } - /* - Otherwise treat valve as an open pipe - */ - else - { - valvecoeff(p,k); //pipecoeff(k); - s->Aij[s->Ndx[k]] -= s->P[k]; - s->Aii[i] += s->P[k]; - s->Aii[j] += s->P[k]; - s->F[i] += (s->Y[k] - hyd->LinkFlows[k]); - s->F[j] -= (s->Y[k] - hyd->LinkFlows[k]); - } -} /* End of fcvcoeff */ - - -/*** New function added. ***/ -void valvecoeff(EN_Project *pr, int k) -/* -**-------------------------------------------------------------- -** Input: k = link index -** Output: none -** Purpose: computes solution matrix coeffs. for a completely -** open, closed, or throttled control valve. -**-------------------------------------------------------------- -*/ -{ - double p; - - EN_Network *n = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - solver_t *s = &hyd->solver; - - Slink *link = &n->Link[k]; - - double flow = hyd->LinkFlows[k]; - - // Valve is closed. Use a very small matrix coeff. - if (hyd->LinkStatus[k] <= CLOSED) - { - s->P[k] = 1.0/CBIG; - s->Y[k] = flow; - return; - } - - // Account for any minor headloss through the valve - if (link->Km > 0.0) - { - p = 2.0 * link->Km * fabs(flow); - if ( p < hyd->RQtol ) { - p = hyd->RQtol; - } - s->P[k] = 1.0 / p; - s->Y[k] = flow / 2.0; - } - else - { - s->P[k] = 1.0 / hyd->RQtol; - s->Y[k] = flow; - } -} - /**************** END OF HYDRAUL.C ***************/ - diff --git a/src/hydsolver.c b/src/hydsolver.c new file mode 100644 index 0000000..0adf780 --- /dev/null +++ b/src/hydsolver.c @@ -0,0 +1,1049 @@ +/* +********************************************************************* + +HYDSOLVER.C -- Equilibrium hydraulic solver for the EPANET Program + +This module contains a pipe network hydraulic solver that computes +flows and pressures within the network at a specific point in time. + +The solver implements Todini's Global Gradient Algorithm. +******************************************************************* +*/ + +#include +#include +#ifndef __APPLE__ +#include +#else +#include +#endif +#include +#include "text.h" +#include "types.h" +#include "funcs.h" + +// Error in network being hydraulically balanced +typedef struct { + double maxheaderror; + double maxflowerror; + double maxflowchange; + int maxheadlink; + int maxflownode; + int maxflowlink; +} Hydbalance; + +// External functions +//int hydsolve(EN_Project *pr, int *iter, double *relerr); + +// Local functions +static int badvalve(EN_Project *pr, int); +static int valvestatus(EN_Project *pr); +static int linkstatus(EN_Project *pr); +static StatType cvstatus(EN_Project *pr, StatType, double, double); +static StatType pumpstatus(EN_Project *pr, int, double); +static StatType prvstatus(EN_Project *pr, int, StatType, double, double, double); +static StatType psvstatus(EN_Project *pr, int, StatType, double, double, double); +static StatType fcvstatus(EN_Project *pr, int, StatType, double, double); +static void tankstatus(EN_Project *pr, int, int, int); +static int pswitch(EN_Project *pr); + +static double newflows(EN_Project *pr, Hydbalance *hbal); +static double emitflowchange(EN_Project *pr, int i); +static void checkhydbalance(EN_Project *pr, Hydbalance *hbal); +static int hasconverged(EN_Project *pr, double *relerr, Hydbalance *hbal); +static void reporthydbal(EN_Project *pr, Hydbalance *hbal); + + +int hydsolve(EN_Project *pr, int *iter, double *relerr) +/* +**------------------------------------------------------------------- +** Input: none +** Output: *iter = # of iterations to reach solution +** *relerr = convergence error in solution +** returns error code +** Purpose: solves network nodal equations for heads and flows +** using Todini's Gradient algorithm +** +*** Updated 9/7/00 *** +*** Updated 2.00.11 *** +*** Updated 2.00.12 *** +** Notes: Status checks on CVs, pumps and pipes to tanks are made +** every hyd->CheckFreq iteration, up until MaxCheck iterations +** are reached. Status checks on control valves are made +** every iteration if DampLimit = 0 or only when the +** convergence error is at or below DampLimit. If DampLimit +** is > 0 then future computed flow changes are only 60% of +** their full value. A complete status check on all links +** is made when convergence is achieved. If convergence is +** not achieved in hyd->MaxIter trials and hyd->ExtraIter > 0 then +** another hyd->ExtraIter trials are made with no status changes +** made to any links and a warning message is generated. +** +** This procedure calls linsolve() which appears in SMATRIX.C. +**------------------------------------------------------------------- +*/ +{ + int i; /* Node index */ + int errcode = 0; /* Node causing solution error */ + int nextcheck; /* Next status check trial */ + int maxtrials; /* Max. trials for convergence */ + double newerr; /* New convergence error */ + int valveChange; /* Valve status change flag */ + int statChange; /* Non-valve status change flag */ + Hydbalance hydbal; /* Hydraulic balance errors */ + + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + report_options_t *rep = &pr->report; + + /* Initialize status checking & relaxation factor */ + nextcheck = hyd->CheckFreq; + hyd->RelaxFactor = 1.0; + + /* Repeat iterations until convergence or trial limit is exceeded. */ + /* (hyd->ExtraIter used to increase trials in case of status cycling.) */ + if (pr->report.Statflag == FULL) { + writerelerr(pr, 0, 0); + } + maxtrials = hyd->MaxIter; + if (hyd->ExtraIter > 0) { + maxtrials += hyd->ExtraIter; + } + *iter = 1; + while (*iter <= maxtrials) { + /* + ** Compute coefficient matrices A & F and solve A*H = F + ** where H = heads, A = Jacobian coeffs. derived from + ** head loss gradients, & F = flow correction terms. + ** Solution for H is returned in F from call to linsolve(). + */ + for (i = 1; i <= net->Nlinks; i++) { + hlosscoeff(pr, i); + } + matrixcoeffs(pr); + errcode = linsolve(&hyd->solver, net->Njuncs); + + /* Take action depending on error code */ + if (errcode < 0) { + break; /* Memory allocation problem */ + } + if (errcode > 0) { /* Ill-conditioning problem */ + /* If control valve causing problem, fix its status & continue, */ + /* otherwise end the iterations with no solution. */ + if (badvalve(pr, sol->Order[errcode])) { + continue; + } + else break; + } + + /* Update current solution. */ + /* (Row[i] = row of solution matrix corresponding to node i). */ + for (i = 1; i <= net->Njuncs; i++) { + hyd->NodeHead[i] = sol->F[sol->Row[i]]; /* Update heads */ + } + newerr = newflows(pr, &hydbal); /* Update flows */ + *relerr = newerr; + + /* Write convergence error to status report if called for */ + if (rep->Statflag == FULL) { + writerelerr(pr, *iter, *relerr); + } + + /* Apply solution damping & check for change in valve status */ + hyd->RelaxFactor = 1.0; + valveChange = FALSE; + if (hyd->DampLimit > 0.0) { + if (*relerr <= hyd->DampLimit) { + hyd->RelaxFactor = 0.6; + valveChange = valvestatus(pr); + } + } + else { + valveChange = valvestatus(pr); + } + + /* Check for convergence */ + if (hasconverged(pr, relerr, &hydbal)) { + /* We have convergence. Quit if we are into extra iterations. */ + if (*iter > hyd->MaxIter) { + break; + } + + /* Quit if no status changes occur. */ + statChange = FALSE; + if (valveChange) { + statChange = TRUE; + } + if (linkstatus(pr)) { + statChange = TRUE; + } + if (pswitch(pr)) { + statChange = TRUE; + } + if (!statChange) { + break; + } + + /* We have a status change so continue the iterations */ + nextcheck = *iter + hyd->CheckFreq; + } + + /* No convergence yet. See if its time for a periodic status */ + /* check on pumps, CV's, and pipes connected to tanks. */ + else if (*iter <= hyd->MaxCheck && *iter == nextcheck) + { + linkstatus(pr); + nextcheck += hyd->CheckFreq; + } + (*iter)++; + } + + /* Iterations ended. Report any errors. */ + if (errcode < 0) errcode = 101; /* Memory allocation error */ + else if (errcode > 0) + { + writehyderr(pr, sol->Order[errcode]); /* Ill-conditioned eqns. error */ + errcode = 110; + } + + /* Add any emitter flows to junction demands */ + for (i = 1; i <= net->Njuncs; i++) { + hyd->NodeDemand[i] += hyd->EmitterFlows[i]; + } + return(errcode); +} /* End of netsolve */ + + +int badvalve(EN_Project *pr, int n) +/* +**----------------------------------------------------------------- +** Input: n = node index +** Output: returns 1 if node n belongs to an active control valve, +** 0 otherwise +** Purpose: determines if a node belongs to an active control valve +** whose setting causes an inconsistent set of eqns. If so, +** the valve status is fixed open and a warning condition +** is generated. +**----------------------------------------------------------------- +*/ +{ + int i, k, n1, n2; + + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + report_options_t *rep = &pr->report; + time_options_t *time = &pr->time_options; + Slink *link; + + for (i = 1; i <= net->Nvalves; i++) { + k = net->Valve[i].Link; + link = &net->Link[k]; + n1 = link->N1; + n2 = link->N2; + if (n == n1 || n == n2) { + EN_LinkType t = link->Type; + if (t == EN_PRV || t == EN_PSV || t == EN_FCV) { + if (hyd->LinkStatus[k] == ACTIVE) { + if (rep->Statflag == FULL) { + sprintf(pr->Msg, FMT61, clocktime(rep->Atime, time->Htime), link->ID); + writeline(pr, pr->Msg); + } + if (link->Type == EN_FCV) { + hyd->LinkStatus[k] = XFCV; + } + else { + hyd->LinkStatus[k] = XPRESSURE; + } + return(1); + } + } + return(0); + } + } + return(0); +} + + +int valvestatus(EN_Project *pr) +/* +**----------------------------------------------------------------- +** Input: none +** Output: returns 1 if any pressure or flow control valve +** changes status, 0 otherwise +** Purpose: updates status for PRVs & PSVs whose status +** is not fixed to OPEN/CLOSED +**----------------------------------------------------------------- +*/ +{ + int change = FALSE, /* Status change flag */ + i, k, /* Valve & link indexes */ + n1, n2; /* Start & end nodes */ + StatType status; /* Valve status settings */ + double hset; /* Valve head setting */ + Slink *link; + + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + report_options_t *rep = &pr->report; + + for (i = 1; i <= net->Nvalves; i++) /* Examine each valve */ + { + k = net->Valve[i].Link; /* Link index of valve */ + link = &net->Link[k]; + if (hyd->LinkSetting[k] == MISSING) { + continue; /* Valve status fixed */ + } + n1 = link->N1; /* Start & end nodes */ + n2 = link->N2; + status = hyd->LinkStatus[k]; /* Save current status */ + +// if (s != CLOSED /* No change if flow is */ +// && ABS(hyd->LinkFlows[k]) < hyd->Qtol) continue; /* negligible. */ + + switch (link->Type) /* Evaluate new status: */ + { + case EN_PRV: + hset = net->Node[n2].El + hyd->LinkSetting[k]; + hyd->LinkStatus[k] = prvstatus(pr, k, status, hset, hyd->NodeHead[n1], hyd->NodeHead[n2]); + break; + case EN_PSV: + hset = net->Node[n1].El + hyd->LinkSetting[k]; + hyd->LinkStatus[k] = psvstatus(pr, k, status, hset, hyd->NodeHead[n1], hyd->NodeHead[n2]); + break; + +//// FCV status checks moved back into the linkstatus() function //// +// case FCV: S[k] = fcvstatus(k,s,hyd->NodeHead[n1],hyd->NodeHead[n2]); +// break; + + default: + continue; + } + + /*** Updated 9/7/00 ***/ + /* Do not reset flow in valve if its status changes. */ + /* This strategy improves convergence. */ + + /* Check for status change */ + if (status != hyd->LinkStatus[k]) + { + if (rep->Statflag == FULL) { + writestatchange(pr, k, status, hyd->LinkStatus[k]); + } + change = TRUE; + } + } + return(change); +} /* End of valvestatus() */ + + +/*** Updated 9/7/00 ***/ +int linkstatus(EN_Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: returns 1 if any link changes status, 0 otherwise +** Purpose: determines new status for pumps, CVs, FCVs & pipes +** to tanks. +**-------------------------------------------------------------- +*/ +{ + int change = FALSE, /* Status change flag */ + k, /* Link index */ + n1, /* Start node index */ + n2; /* End node index */ + double dh; /* Head difference */ + StatType status; /* Current status */ + + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + report_options_t *rep = &pr->report; + Slink *link; + + /* Examine each Slink */ + for (k = 1; k <= net->Nlinks; k++) + { + link = &net->Link[k]; + n1 = link->N1; + n2 = link->N2; + dh = hyd->NodeHead[n1] - hyd->NodeHead[n2]; + + /* Re-open temporarily closed links (status = XHEAD or TEMPCLOSED) */ + status = hyd->LinkStatus[k]; + if (status == XHEAD || status == TEMPCLOSED) { + hyd->LinkStatus[k] = OPEN; + } + + /* Check for status changes in CVs and pumps */ + if (link->Type == EN_CVPIPE) { + hyd->LinkStatus[k] = cvstatus(pr, hyd->LinkStatus[k], dh, hyd->LinkFlows[k]); + } + if (link->Type == EN_PUMP && hyd->LinkStatus[k] >= OPEN && hyd->LinkSetting[k] > 0.0) { + hyd->LinkStatus[k] = pumpstatus(pr, k, -dh); + } + + /* Check for status changes in non-fixed FCVs */ + if (link->Type == EN_FCV && hyd->LinkSetting[k] != MISSING) { + hyd->LinkStatus[k] = fcvstatus(pr, k, status, hyd->NodeHead[n1], hyd->NodeHead[n2]); // + } + + /* Check for flow into (out of) full (empty) tanks */ + if (n1 > net->Njuncs || n2 > net->Njuncs) { + tankstatus(pr, k, n1, n2); + } + + /* Note change in link status; do not revise link flow */ + if (status != hyd->LinkStatus[k]) { + change = TRUE; + if (rep->Statflag == FULL) { + writestatchange(pr, k, status, hyd->LinkStatus[k]); + } + + //if (S[k] <= CLOSED) hyd->LinkFlows[k] = QZERO; + //else setlinkflow(k, dh); + } + } + return(change); +} /* End of linkstatus */ + + +StatType cvstatus(EN_Project *pr, StatType s, double dh, double q) +/* +**-------------------------------------------------- +** Input: s = current status +** dh = headloss +** q = flow +** Output: returns new link status +** Purpose: updates status of a check valve. +**-------------------------------------------------- +*/ +{ + hydraulics_t *hyd = &pr->hydraulics; + + /* Prevent reverse flow through CVs */ + if (ABS(dh) > hyd->Htol) + { + if (dh < -hyd->Htol) { + return(CLOSED); + } + else if (q < -hyd->Qtol) { + return(CLOSED); + } + else { + return(OPEN); + } + } + else + { + if (q < -hyd->Qtol) { + return(CLOSED); + } + else { + return(s); + } + } +} + + +StatType pumpstatus(EN_Project *pr, int k, double dh) +/* +**-------------------------------------------------- +** Input: k = link index +** dh = head gain +** Output: returns new pump status +** Purpose: updates status of an open pump. +**-------------------------------------------------- +*/ +{ + int p; + double hmax; + hydraulics_t *hyd = &pr->hydraulics; + EN_Network *net = &pr->network; + + /* Prevent reverse flow through pump */ + p = findpump(net, k); + if (net->Pump[p].Ptype == CONST_HP) { + hmax = BIG; + } + else { + hmax = SQR(hyd->LinkSetting[k]) * net->Pump[p].Hmax; + } + if (dh > hmax + hyd->Htol) { + return(XHEAD); + } + + /*** Flow higher than pump curve no longer results in a status change ***/ + /* Check if pump cannot deliver flow */ + //if (hyd->LinkFlows[k] > K[k]*Pump[p].Qmax + hyd->Qtol) return(XFLOW); + return(OPEN); +} + + +StatType prvstatus(EN_Project *pr, int k, StatType s, double hset, double h1, double h2) +/* +**----------------------------------------------------------- +** Input: k = link index +** s = current status +** hset = valve head setting +** h1 = head at upstream node +** h2 = head at downstream node +** Output: returns new valve status +** Purpose: updates status of a pressure reducing valve. +**----------------------------------------------------------- +*/ +{ + StatType status; /* New valve status */ + double hml; /* Minor headloss */ + hydraulics_t *hyd = &pr->hydraulics; + + double htol = hyd->Htol; + Slink *link = &pr->network.Link[k]; + + status = s; + if (hyd->LinkSetting[k] == MISSING) + return(status); /* Status fixed by user */ + hml = link->Km*SQR(hyd->LinkFlows[k]); /* Head loss when open */ + + /*** Status rules below have changed. ***/ + + switch (s) + { + case ACTIVE: + if (hyd->LinkFlows[k] < -hyd->Qtol) { + status = CLOSED; + } + else if (h1 - hml < hset - htol) { + status = OPEN; + } + else { + status = ACTIVE; + } + break; + case OPEN: + if (hyd->LinkFlows[k] < -hyd->Qtol) { + status = CLOSED; + } + else if (h2 >= hset + htol) { + status = ACTIVE; + } + else { + status = OPEN; + } + break; + case CLOSED: + if (h1 >= hset + htol && h2 < hset - htol) { + status = ACTIVE; + } + else if (h1 < hset - htol && h1 > h2 + htol) { + status = OPEN; + } + else { + status = CLOSED; + } + break; + case XPRESSURE: + if (hyd->LinkFlows[k] < -hyd->Qtol) { + status = CLOSED; + } + break; + default: + break; + } + return(status); +} + + +StatType psvstatus(EN_Project *pr, int k, StatType s, double hset, double h1, double h2) +/* +**----------------------------------------------------------- +** Input: k = link index +** s = current status +** hset = valve head setting +** h1 = head at upstream node +** h2 = head at downstream node +** Output: returns new valve status +** Purpose: updates status of a pressure sustaining valve. +**----------------------------------------------------------- +*/ +{ + StatType status; /* New valve status */ + double hml; /* Minor headloss */ + hydraulics_t *hyd = &pr->hydraulics; + + double htol = hyd->Htol; + + Slink *link = &pr->network.Link[k]; + + status = s; + if (hyd->LinkSetting[k] == MISSING) { + return(status); /* Status fixed by user */ + } + hml = link->Km*SQR(hyd->LinkFlows[k]); /* Head loss when open */ + +/*** Status rules below have changed. ***/ + + switch (s) + { + case ACTIVE: + if (hyd->LinkFlows[k] < -hyd->Qtol) { + status = CLOSED; + } + else if (h2 + hml > hset + htol) { + status = OPEN; + } + else { + status = ACTIVE; + } + break; + case OPEN: + if (hyd->LinkFlows[k] < -hyd->Qtol) { + status = CLOSED; + } + else if (h1 < hset - htol) { + status = ACTIVE; + } + else { + status = OPEN; + } + break; + case CLOSED: + if (h2 > hset + htol && h1 > h2 + htol) { + status = OPEN; + } + else if (h1 >= hset + htol && h1 > h2 + htol) { + status = ACTIVE; + } + else { + status = CLOSED; + } + break; + case XPRESSURE: + if (hyd->LinkFlows[k] < -hyd->Qtol) { + status = CLOSED; + } + break; + default: + break; + } + return(status); +} + + +StatType fcvstatus(EN_Project *pr, int k, StatType s, double h1, double h2) +/* +**----------------------------------------------------------- +** Input: k = link index +** s = current status +** h1 = head at upstream node +** h2 = head at downstream node +** Output: returns new valve status +** Purpose: updates status of a flow control valve. +** +** Valve status changes to XFCV if flow reversal. +** If current status is XFCV and current flow is +** above setting, then valve becomes active. +** If current status is XFCV, and current flow +** positive but still below valve setting, then +** status remains same. +**----------------------------------------------------------- +*/ +{ + StatType status; /* New valve status */ + hydraulics_t *hyd = &pr->hydraulics; + + status = s; + if (h1 - h2 < -hyd->Htol) { + status = XFCV; + } + else if (hyd->LinkFlows[k] < -hyd->Qtol) { + status = XFCV; + } + else if (s == XFCV && hyd->LinkFlows[k] >= hyd->LinkSetting[k]) { + status = ACTIVE; + } + return(status); +} + + +/*** Updated 9/7/00 ***/ +/*** Updated 11/19/01 ***/ +void tankstatus(EN_Project *pr, int k, int n1, int n2) +/* +**---------------------------------------------------------------- +** Input: k = link index +** n1 = start node of link +** n2 = end node of link +** Output: none +** Purpose: closes link flowing into full or out of empty tank +**---------------------------------------------------------------- +*/ +{ + int i, n; + double h, q; + Stank *tank; + + hydraulics_t *hyd = &pr->hydraulics; + EN_Network *net = &pr->network; + Slink *link = &net->Link[k]; + + /* Make node n1 be the tank */ + q = hyd->LinkFlows[k]; + i = n1 - net->Njuncs; + if (i <= 0) { + i = n2 - net->Njuncs; + if (i <= 0) { + return; + } + n = n1; + n1 = n2; + n2 = n; + q = -q; + } + h = hyd->NodeHead[n1] - hyd->NodeHead[n2]; + tank = &net->Tank[i]; + /* Skip reservoirs & closed links */ + if (tank->A == 0.0 || hyd->LinkStatus[k] <= CLOSED) { + return; + } + + /* If tank full, then prevent flow into it */ + if (hyd->NodeHead[n1] >= tank->Hmax - hyd->Htol) + { + + /* Case 1: Link is a pump discharging into tank */ + if (link->Type == EN_PUMP) { + if (link->N2 == n1) + hyd->LinkStatus[k] = TEMPCLOSED; + } + + /* Case 2: Downstream head > tank head */ + /* (i.e., an open outflow check valve would close) */ + else if (cvstatus(pr, OPEN, h, q) == CLOSED) { + hyd->LinkStatus[k] = TEMPCLOSED; + } + } + + /* If tank empty, then prevent flow out of it */ + if (hyd->NodeHead[n1] <= tank->Hmin + hyd->Htol) { + + /* Case 1: Link is a pump discharging from tank */ + if (link->Type == EN_PUMP) { + if (link->N1 == n1) { + hyd->LinkStatus[k] = TEMPCLOSED; + } + } + + /* Case 2: Tank head > downstream head */ + /* (i.e., a closed outflow check valve would open) */ + else if (cvstatus(pr, CLOSED, h, q) == OPEN) { + hyd->LinkStatus[k] = TEMPCLOSED; + } + } +} /* End of tankstatus */ + + +int pswitch(EN_Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: returns 1 if status of any link changes, 0 if not +** Purpose: adjusts settings of links controlled by junction +** pressures after a hydraulic solution is found +**-------------------------------------------------------------- +*/ +{ + int i, /* Control statement index */ + k, /* Link being controlled */ + n, /* Node controlling Slink */ + reset, /* Flag on control conditions */ + change, /* Flag for status or setting change */ + anychange = 0; /* Flag for 1 or more changes */ + char s; /* Current link status */ + + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + report_options_t *rep = &pr->report; + Slink *link; + + /* Check each control statement */ + for (i = 1; i <= net->Ncontrols; i++) + { + reset = 0; + if ((k = net->Control[i].Link) <= 0) { + continue; + } + + /* Determine if control based on a junction, not a tank */ + if ((n = net->Control[i].Node) > 0 && n <= net->Njuncs) { + /* Determine if control conditions are satisfied */ + if (net->Control[i].Type == LOWLEVEL + && hyd->NodeHead[n] <= net->Control[i].Grade + hyd->Htol) { + reset = 1; + } + if (net->Control[i].Type == HILEVEL + && hyd->NodeHead[n] >= net->Control[i].Grade - hyd->Htol) { + reset = 1; + } + } + + /* Determine if control forces a status or setting change */ + if (reset == 1) + { + link = &net->Link[k]; + change = 0; + s = hyd->LinkStatus[k]; + if (link->Type == EN_PIPE) { + if (s != net->Control[i].Status) { + change = 1; + } + } + if (link->Type == EN_PUMP) { + if (hyd->LinkSetting[k] != net->Control[i].Setting) { + change = 1; + } + } + if (link->Type >= EN_PRV) { + if (hyd->LinkSetting[k] != net->Control[i].Setting) { + change = 1; + } + else if (hyd->LinkSetting[k] == MISSING && s != net->Control[i].Status) { + change = 1; + } + } + + /* If a change occurs, update status & setting */ + if (change) { + hyd->LinkStatus[k] = net->Control[i].Status; + if (link->Type > EN_PIPE) { + hyd->LinkSetting[k] = net->Control[i].Setting; + } + if (rep->Statflag == FULL) { + writestatchange(pr, k, s, hyd->LinkStatus[k]); + } + + /* Re-set flow if status has changed */ + // if (S[k] != s) initlinkflow(k, S[k], K[k]); + anychange = 1; + } + } + } + return(anychange); +} /* End of pswitch */ + + +double newflows(EN_Project *pr, Hydbalance *hbal) +/* +**---------------------------------------------------------------- +** Input: none +** Output: returns solution convergence error +** Purpose: updates link flows after new nodal heads computed +**---------------------------------------------------------------- +*/ +{ + double dh, /* Link head loss */ + dq; /* Link flow change */ + double dqsum, /* Network flow change */ + qsum; /* Network total flow */ + int k, n, n1, n2; + + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + Slink *link; + + /* Initialize net inflows (i.e., demands) at tanks */ + for (n = net->Njuncs + 1; n <= net->Nnodes; n++) { + hyd->NodeDemand[n] = 0.0; + } + + /* Initialize sum of flows & corrections */ + qsum = 0.0; + dqsum = 0.0; + + hbal->maxflowchange = 0.0; + hbal->maxflowlink = 1; + + /* Update flows in all links */ + for (k = 1; k <= net->Nlinks; k++) + { + link = &net->Link[k]; + /* + ** Apply flow update formula: + ** dq = Y - P*(new head loss) + ** P = 1/(dh/dq) + ** Y = P*(head loss based on current flow) + ** where P & Y were computed in newcoeffs(). + */ + + n1 = link->N1; + n2 = link->N2; + dh = hyd->NodeHead[n1] - hyd->NodeHead[n2]; + dq = sol->Y[k] - sol->P[k] * dh; + + /* Adjust flow change by the relaxation factor */ + dq *= hyd->RelaxFactor; + + /* Prevent flow in constant HP pumps from going negative */ + if (link->Type == EN_PUMP) { + n = findpump(net, k); + if (net->Pump[n].Ptype == CONST_HP && dq > hyd->LinkFlows[k]) { + dq = hyd->LinkFlows[k] / 2.0; + } + } + hyd->LinkFlows[k] -= dq; + + /* Update sum of absolute flows & flow corrections */ + qsum += ABS(hyd->LinkFlows[k]); + dqsum += ABS(dq); + + /* Update identity of link with max. flow change */ + if (ABS(dq) > hbal->maxflowchange) { + hbal->maxflowchange = ABS(dq); + hbal->maxflowlink = k; + } + + /* Update net flows to tanks */ + if (hyd->LinkStatus[k] > CLOSED) + { + if (n1 > net->Njuncs) { + hyd->NodeDemand[n1] -= hyd->LinkFlows[k]; + } + if (n2 > net->Njuncs) { + hyd->NodeDemand[n2] += hyd->LinkFlows[k]; + } + } + } + + /* Update emitter flows */ + for (k = 1; k <= net->Njuncs; k++) + { + if (net->Node[k].Ke == 0.0) { + continue; + } + dq = emitflowchange(pr, k); + hyd->EmitterFlows[k] -= dq; + qsum += ABS(hyd->EmitterFlows[k]); + dqsum += ABS(dq); + } + + /* Return ratio of total flow corrections to total flow */ + if (qsum > hyd->Hacc) { + return(dqsum / qsum); + } + else { + return(dqsum); + } +} /* End of newflows */ + + +double emitflowchange(EN_Project *pr, int i) +/* +**-------------------------------------------------------------- +** Input: i = node index +** Output: returns change in flow at an emitter node +** Purpose: computes flow change at an emitter node +**-------------------------------------------------------------- +*/ +{ + double ke, p; + hydraulics_t *hyd = &pr->hydraulics; + EN_Network *n = &pr->network; + Snode *node = &n->Node[i]; + + ke = MAX(CSMALL, node->Ke); + p = hyd->Qexp * ke * pow(ABS(hyd->EmitterFlows[i]), (hyd->Qexp - 1.0)); + if (p < hyd->RQtol) { + p = 1 / hyd->RQtol; + } + else { + p = 1.0 / p; + } + return(hyd->EmitterFlows[i] / hyd->Qexp - p * (hyd->NodeHead[i] - node->El)); +} + + +void checkhydbalance(EN_Project *pr, Hydbalance *hbal) +/* +**-------------------------------------------------------------- +** Input: hbal = hydraulic balance errors +** Output: none +** Purpose: finds the link with the largest head imbalance +**-------------------------------------------------------------- +*/ +{ + int k, n1, n2; + double dh, headerror, headloss; + + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + Slink *link; + hbal->maxheaderror = 0.0; + hbal->maxheadlink = 1; + for (k = 1; k <= net->Nlinks; k++) { + if (hyd->LinkStatus[k] <= CLOSED) continue; + hlosscoeff(pr, k); + if (sol->P[k] == 0.0) continue; + link = &net->Link[k]; + n1 = link->N1; + n2 = link->N2; + dh = hyd->NodeHead[n1] - hyd->NodeHead[n2]; + headloss = sol->Y[k] / sol->P[k]; + headerror = ABS(dh - headloss); + if (headerror > hbal->maxheaderror) { + hbal->maxheaderror = headerror; + hbal->maxheadlink = k; + } + } +} + + +int hasconverged(EN_Project *pr, double *relerr, Hydbalance *hbal) +/* +**-------------------------------------------------------------- +** Input: relerr = current total relative flow change +** hbal = current hydraulic balance errors +** Output: returns 1 if system has converged or 0 if not +** Purpose: checks various criteria to see if system has +** become hydraulically balanced +**-------------------------------------------------------------- +*/ +{ + hydraulics_t *hyd = &pr->hydraulics; + + if (*relerr > hyd->Hacc) return 0; + checkhydbalance(pr, hbal); + if (pr->report.Statflag == FULL) { + reporthydbal(pr, hbal); + } + if (hyd->HeadErrorLimit > 0.0 && + hbal->maxheaderror > hyd->HeadErrorLimit) return 0; + if (hyd->FlowChangeLimit > 0.0 && + hbal->maxflowchange > hyd->FlowChangeLimit) return 0; + return 1; +} + + +void reporthydbal(EN_Project *pr, Hydbalance *hbal) +/* +**-------------------------------------------------------------- +** Input: hbal = current hydraulic balance errors +** Output: none +** Purpose: identifies links with largest flow change and +** largest head loss error. +**-------------------------------------------------------------- +*/ +{ + double qchange = hbal->maxflowchange * pr->Ucf[FLOW]; + double herror = hbal->maxheaderror * pr->Ucf[HEAD]; + int qlink = hbal->maxflowlink; + int hlink = hbal->maxheadlink; + if (qlink >= 1) { + sprintf(pr->Msg, FMT66, qchange, pr->network.Link[qlink].ID); + writeline(pr, pr->Msg); + } + if (hlink >= 1) { + sprintf(pr->Msg, FMT67, herror, pr->network.Link[hlink].ID); + writeline(pr, pr->Msg); + } +} diff --git a/src/inpfile.c b/src/inpfile.c index 1a2b17d..8bf894d 100644 --- a/src/inpfile.c +++ b/src/inpfile.c @@ -618,6 +618,12 @@ int saveinpfile(EN_Project *pr, char *fname) fprintf(f, "\n CHECKFREQ %-d", hyd->CheckFreq); fprintf(f, "\n MAXCHECK %-d", hyd->MaxCheck); fprintf(f, "\n DAMPLIMIT %-.8f", hyd->DampLimit); + if (hyd->HeadErrorLimit > 0.0) { + fprintf(f, "\n HEADERROR %-.8f", hyd->HeadErrorLimit * pr->Ucf[HEAD]); + } + if (hyd->FlowChangeLimit > 0.0) { + fprintf(f, "\n FLOWCHANGE %-.8f", hyd->FlowChangeLimit * pr->Ucf[FLOW]); + } /* Write [REPORT] section */ diff --git a/src/input1.c b/src/input1.c index 8bb0658..b59096a 100644 --- a/src/input1.c +++ b/src/input1.c @@ -40,9 +40,9 @@ AUTHOR: L. Rossman --------------------- Module Global Variables ---------------------- */ -#define MAXITER 200 /* Default max. # hydraulic iterations */ -#define HACC 0.001 /* Default hydraulics convergence ratio */ -#define HTOL 0.0005 /* Default hydraulic head tolerance (ft) */ +#define MAXITER 200 /* Default max. # hydraulic iterations */ +#define HACC 0.001 /* Default hydraulics convergence ratio */ +#define HTOL 0.0005 /* Default hydraulic head tolerance (ft) */ /*** Updated 11/19/01 ***/ #define QTOL 0.0001 /* Default flow rate tolerance (cfs) */ @@ -117,19 +117,23 @@ void setdefaults(EN_Project *pr) strncpy(qu->ChemUnits, u_MGperL, MAXID); strncpy(par->DefPatID, DEFPATID, MAXID); out->Hydflag = SCRATCH; /* No external hydraulics file */ - qu->Qualflag = NONE; /* No quality simulation */ + qu->Qualflag = NONE; /* No quality simulation */ hyd->Formflag = HW; /* Use Hazen-Williams formula */ par->Unitsflag = US; /* US unit system */ par->Flowflag = GPM; /* Flow units are gpm */ par->Pressflag = PSI; /* Pressure units are psi */ rep->Tstatflag = SERIES; /* Generate time series output */ pr->Warnflag = FALSE; /* Warning flag is off */ - hyd->Htol = HTOL; /* Default head tolerance */ - hyd->Qtol = QTOL; /* Default flow tolerance */ - hyd->Hacc = HACC; /* Default hydraulic accuracy */ + hyd->Htol = HTOL; /* Default head tolerance */ + hyd->Qtol = QTOL; /* Default flow tolerance */ + hyd->Hacc = HACC; /* Default hydraulic accuracy */ + + hyd->FlowChangeLimit = 0.0; /* Default flow change limit */ + hyd->HeadErrorLimit = 0.0; /* Default head error limit */ + qu->Ctol = MISSING; /* No pre-set quality tolerance */ - hyd->MaxIter = MAXITER; /* Default max. hydraulic trials */ - hyd->ExtraIter = -1; /* Stop if network unbalanced */ + hyd->MaxIter = MAXITER; /* Default max. hydraulic trials */ + hyd->ExtraIter = -1; /* Stop if network unbalanced */ time->Dur = 0; /* 0 sec duration (steady state) */ time->Tstart = 0; /* Starting time of day */ time->Pstart = 0; /* Starting pattern period */ @@ -529,8 +533,8 @@ void initunits(EN_Project *pr) pr->Ucf[POWER] = wcf; pr->Ucf[VOLUME] = hcf * hcf * hcf; if (time->Hstep < 1800) /* Report time in mins. */ - { /* if hydraulic time step */ - pr->Ucf[TIME] = 1.0 / 60.0; /* is less than 1/2 hour. */ + { /* if hydraulic time step */ + pr->Ucf[TIME] = 1.0 / 60.0; /* is less than 1/2 hour. */ strcpy(rep->Field[TIME].Units, u_MINUTES); } else { pr->Ucf[TIME] = 1.0 / 3600.0; @@ -604,6 +608,10 @@ void convertunits(EN_Project *pr) tank->V1max *= tank->Vmax; } + /* Convert hydraulic convergence criteria */ + hyd->FlowChangeLimit /= pr->Ucf[FLOW]; + hyd->HeadErrorLimit /= pr->Ucf[HEAD]; + /* Convert WQ option concentration units */ qu->Climit /= pr->Ucf[QUALITY]; qu->Ctol /= pr->Ucf[QUALITY]; @@ -676,8 +684,9 @@ void convertunits(EN_Project *pr) } +////// Moved to inithyd() in hydraul.c /////// /* Compute flow resistances */ - resistance(pr, k); + //resistance(pr, k); } /* Convert units on control settings */ diff --git a/src/input3.c b/src/input3.c index a011af1..9c66f63 100644 --- a/src/input3.c +++ b/src/input3.c @@ -1913,6 +1913,10 @@ int optionvalue(EN_Project *pr, int n) ** SPECIFIC GRAVITY value ** TRIALS value ** ACCURACY value + +** HEADLIMIT value +** FLOWLIMIT value + ** TOLERANCE value ** SEGMENTS value (not used) ** ------ Undocumented Options ----- @@ -1928,18 +1932,20 @@ int optionvalue(EN_Project *pr, int n) hydraulics_t *hyd = &pr->hydraulics; quality_t *qu = &pr->quality; parser_data_t *par = &pr->parser; + char* tok0 = par->Tok[0]; int nvalue = 1; /* Index of token with numerical value */ double y; /* Check for obsolete SEGMENTS keyword */ - if (match(par->Tok[0], w_SEGMENTS)) - return (0); + //if (match(par->Tok[0], w_SEGMENTS)) + if (match(tok0, w_SEGMENTS)) + return (0); /* Check for missing value (which is permissible) */ - if (match(par->Tok[0], w_SPECGRAV) || match(par->Tok[0], w_EMITTER) || - match(par->Tok[0], w_DEMAND)) + if (match(tok0, w_SPECGRAV) || match(tok0, w_EMITTER) || + match(tok0, w_DEMAND)) nvalue = 2; if (n < nvalue) return (0); @@ -1949,7 +1955,7 @@ int optionvalue(EN_Project *pr, int n) return (213); /* Check for WQ tolerance option (which can be 0) */ - if (match(par->Tok[0], w_TOLERANCE)) { + if (match(tok0, w_TOLERANCE)) { if (y < 0.0) return (213); qu->Ctol = y; /* Quality tolerance*/ @@ -1957,7 +1963,7 @@ int optionvalue(EN_Project *pr, int n) } /* Check for Diffusivity option */ - if (match(par->Tok[0], w_DIFFUSIVITY)) { + if (match(tok0, w_DIFFUSIVITY)) { if (y < 0.0) return (213); qu->Diffus = y; @@ -1965,42 +1971,59 @@ int optionvalue(EN_Project *pr, int n) } /* Check for Damping Limit option */ - if (match(par->Tok[0], w_DAMPLIMIT)) { + if (match(tok0, w_DAMPLIMIT)) { hyd->DampLimit = y; return (0); } + + /* Check for flow change limit*/ + else if (match(tok0, w_FLOWCHANGE)) + { + if (y < 0.0) return 213; + hyd->FlowChangeLimit = y; + return 0; + } + + /* Check for head error limit*/ + else if (match(tok0, w_HEADERROR)) + { + if (y < 0.0) return 213; + hyd->HeadErrorLimit = y; + return 0; + } /* All other options must be > 0 */ if (y <= 0.0) return (213); /* Assign value to specified option */ - if (match(par->Tok[0], w_VISCOSITY)) + if (match(tok0, w_VISCOSITY)) hyd->Viscos = y; /* Viscosity */ - else if (match(par->Tok[0], w_SPECGRAV)) + else if (match(tok0, w_SPECGRAV)) hyd->SpGrav = y; /* Spec. gravity */ - else if (match(par->Tok[0], w_TRIALS)) + else if (match(tok0, w_TRIALS)) hyd->MaxIter = (int)y; /* Max. trials */ - else if (match(par->Tok[0], w_ACCURACY)) /* Accuracy */ + else if (match(tok0, w_ACCURACY)) /* Accuracy */ { y = MAX(y, 1.e-5); y = MIN(y, 1.e-1); hyd->Hacc = y; - } else if (match(par->Tok[0], w_HTOL)) + } + else if (match(tok0, w_HTOL)) hyd->Htol = y; - else if (match(par->Tok[0], w_QTOL)) + else if (match(tok0, w_QTOL)) hyd->Qtol = y; - else if (match(par->Tok[0], w_RQTOL)) { + else if (match(tok0, w_RQTOL)) { if (y >= 1.0) return (213); hyd->RQtol = y; - } else if (match(par->Tok[0], w_CHECKFREQ)) + } else if (match(tok0, w_CHECKFREQ)) hyd->CheckFreq = (int)y; - else if (match(par->Tok[0], w_MAXCHECK)) + else if (match(tok0, w_MAXCHECK)) hyd->MaxCheck = (int)y; - else if (match(par->Tok[0], w_EMITTER)) + else if (match(tok0, w_EMITTER)) hyd->Qexp = 1.0 / y; - else if (match(par->Tok[0], w_DEMAND)) + else if (match(tok0, w_DEMAND)) hyd->Dmult = y; else return (201); diff --git a/src/output.c b/src/output.c index 91c54ea..39687cf 100644 --- a/src/output.c +++ b/src/output.c @@ -252,60 +252,75 @@ int saveenergy(EN_Project *pr) */ { - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - out_file_t *out = &pr->out_files; - parser_data_t *par = &pr->parser; + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + out_file_t *out = &pr->out_files; + parser_data_t *par = &pr->parser; time_options_t *time = &pr->time_options; - FILE *outFile = out->OutFile; + FILE *outFile = out->OutFile; + Spump *pump; int i, j; INT4 index; - REAL4 x[6]; /* work array */ - double hdur, /* total time->Duration in hours */ - t; /* pumping time->Duration */ + REAL4 x[MAX_ENERGY_STATS]; // work array + double hdur, // total simulation duration in hours + t; // total pumping time duration hdur = time->Dur / 3600.0; for (i = 1; i <= net->Npumps; i++) { - Spump *pump = &net->Pump[i]; + pump = &net->Pump[i]; if (hdur == 0.0) { - for (j = 0; j < 5; j++) - x[j] = (REAL4)pump->Energy[j]; - x[5] = (REAL4)(pump->Energy[5] * 24.0); - } else { - t = pump->Energy[0]; - x[0] = (REAL4)(t / hdur); - x[1] = 0.0f; - x[2] = 0.0f; - x[3] = 0.0f; - x[4] = 0.0f; - if (t > 0.0) { - x[1] = (REAL4)(pump->Energy[1] / t); - x[2] = (REAL4)(pump->Energy[2] / t); - x[3] = (REAL4)(pump->Energy[3] / t); - } - x[4] = (REAL4)pump->Energy[4]; - x[5] = (REAL4)(pump->Energy[5] * 24.0 / hdur); + pump->Energy[TOTAL_COST] *= 24.0; } - x[0] *= 100.0f; - x[1] *= 100.0f; - /* Compute Kw-hr per MilGal (or per cubic meter) */ - if (par->Unitsflag == SI) - x[2] *= (REAL4)(1000.0 / LPSperCFS / 3600.0); - else - x[2] *= (REAL4)(1.0e6 / GPMperCFS / 60.0); - for (j = 0; j < 6; j++) - pump->Energy[j] = x[j]; + else { + // ... convert total hrs. online to fraction of total time online + t = pump->Energy[PCNT_ONLINE]; //currently holds total hrs. online + pump->Energy[PCNT_ONLINE] = t / hdur; + + // ... convert cumulative values to time-averaged ones + if (t > 0.0) { + pump->Energy[PCNT_EFFIC] /= t; + pump->Energy[KWH_PER_FLOW] /= t; + pump->Energy[TOTAL_KWH] /= t; + } + + // ... convert total cost to cost per day + pump->Energy[TOTAL_COST] *= 24.0 / hdur; + } + + // ... express time online and avg. efficiency as percentages + pump->Energy[PCNT_ONLINE] *= 100.0; + pump->Energy[PCNT_EFFIC] *= 100.0; + + // ... compute KWH per Million Gallons or per Cubic Meter + if (par->Unitsflag == SI) { + pump->Energy[KWH_PER_FLOW] *= (1000. / LPSperCFS / 3600.); + } + else { + pump->Energy[KWH_PER_FLOW] *= (1.0e6 / GPMperCFS / 60.); + } + + // ... save energy stats to REAL4 work array + for (j = 0; j < MAX_ENERGY_STATS; j++) { + x[j] = (REAL4)pump->Energy[j]; + } + + // ... save energy results to output file index = pump->Link; - if (fwrite(&index, sizeof(INT4), 1, outFile) < 1) + if (fwrite(&index, sizeof(INT4), 1, outFile) < 1) { return (308); - if (fwrite(x, sizeof(REAL4), 6, outFile) < 6) + } + if (fwrite(x, sizeof(REAL4), MAX_ENERGY_STATS, outFile) < MAX_ENERGY_STATS) { return (308); + } } + + // ... compute and save demand charge hyd->Emax = hyd->Emax * hyd->Dcost; x[0] = (REAL4)hyd->Emax; - if (fwrite(&x[0], sizeof(REAL4), 1, outFile) < 1) + if (fwrite(&x[0], sizeof(REAL4), 1, outFile) < 1) { return (308); + } return (0); } diff --git a/src/report.c b/src/report.c index 9620f1f..8d2101e 100644 --- a/src/report.c +++ b/src/report.c @@ -229,6 +229,15 @@ void writesummary(EN_Project *pr) sprintf(s, FMT27, hyd->Hacc); writeline(pr, s); + if (hyd->HeadErrorLimit > 0.0) { + sprintf(s, FMT27d, hyd->HeadErrorLimit*pr->Ucf[HEAD], rep->Field[HEAD].Units); + writeline(pr, s); + } + if (hyd->FlowChangeLimit > 0.0) { + sprintf(s, FMT27e, hyd->FlowChangeLimit*pr->Ucf[FLOW], rep->Field[FLOW].Units); + writeline(pr, s); + } + sprintf(s, FMT27a, hyd->CheckFreq); writeline(pr, s); sprintf(s, FMT27b, hyd->MaxCheck); @@ -378,6 +387,7 @@ void writeenergy(EN_Project *pr) EN_Network *net = &pr->network; hydraulics_t *hyd = &pr->hydraulics; report_options_t *rep = &pr->report; + Spump *pump; int j; double csum; @@ -388,14 +398,14 @@ void writeenergy(EN_Project *pr) writeheader(pr,ENERHDR, 0); csum = 0.0; for (j = 1; j <= net->Npumps; j++) { - Spump *pump = &net->Pump[j]; - csum += pump->Energy[5]; + pump = &net->Pump[j]; + csum += pump->Energy[TOTAL_COST]; if (rep->LineNum == (long)rep->PageSize) writeheader(pr, ENERHDR, 1); sprintf(s, "%-8s %6.2f %6.2f %9.2f %9.2f %9.2f %9.2f", - net->Link[pump->Link].ID, pump->Energy[0], pump->Energy[1], - pump->Energy[2], pump->Energy[3], pump->Energy[4], - pump->Energy[5]); + net->Link[pump->Link].ID, pump->Energy[PCNT_ONLINE], pump->Energy[PCNT_EFFIC], + pump->Energy[KWH_PER_FLOW], pump->Energy[TOTAL_KWH], pump->Energy[MAX_KW], + pump->Energy[TOTAL_COST]); writeline(pr, s); } fillstr(s, '-', 63); diff --git a/src/text.h b/src/text.h index be47e39..85159a9 100755 --- a/src/text.h +++ b/src/text.h @@ -140,7 +140,10 @@ AUTHOR: L. Rossman #define w_RQTOL "RQTOL" #define w_CHECKFREQ "CHECKFREQ" #define w_MAXCHECK "MAXCHECK" -#define w_DAMPLIMIT "DAMPLIMIT" +#define w_DAMPLIMIT "DAMPLIMIT" + +#define w_FLOWCHANGE "FLOWCHANGE" +#define w_HEADERROR "HEADERROR" #define w_SECONDS "SEC" #define w_MINUTES "MIN" @@ -364,7 +367,10 @@ AUTHOR: L. Rossman #define FMT27a " Status Check Frequency ............ %-d" #define FMT27b " Maximum Trials Checked ............ %-d" -#define FMT27c " Damping Limit Threshold ........... %-.6f" +#define FMT27c " Damping Limit Threshold ........... %-.6f" + +#define FMT27d " Headloss Error Limit .............. %-.6f %s" +#define FMT27e " Flow Change Limit ................. %-.6f %s" #define FMT28 " Maximum Trials .................... %-d" #define FMT29 " Quality Analysis .................. None" @@ -409,10 +415,13 @@ AUTHOR: L. Rossman #define FMT61 "%10s: Valve %s caused ill-conditioning" #define FMT62 "%10s: System ill-conditioned at node %s" #define FMT63 "%10s: %s %s changed by rule %s" -#define FMT64 "%10s: Balancing the network:" +#define FMT64 "%10s: Balancing the network:\n" #define FMT65 " Trial %2d: relative flow change = %-.6f" /*** End of update ***/ +#define FMT66 " maximum flow change = %.4f for Link %s" +#define FMT67 " maximum head error = %.4f for Link %s\n" + /* -------------------- Energy Report Table ------------------- */ #define FMT71 "Energy Usage:" #define FMT72 \ diff --git a/src/types.h b/src/types.h index a3aa962..3c26393 100755 --- a/src/types.h +++ b/src/types.h @@ -47,7 +47,7 @@ typedef int INT4; #define EOFMARK 0x1A /* Use 0x04 for UNIX systems */ #define MAXTITLE 3 /* Max. # title lines */ #define MAXID 31 /* Max. # characters in ID name */ -#define MAXMSG 79 /* Max. # characters in message text */ +#define MAXMSG 79 /* Max. # characters in message text */ #define MAXLINE 255 /* Max. # characters read from input line */ #define MAXFNAME 259 /* Max. # characters in file name */ #define MAXTOKS 40 /* Max. items per line of input */ @@ -59,6 +59,9 @@ typedef int INT4; #define TINY 1.E-6 #define MISSING -1.E10 +#define CBIG 1.e8 /* Big coefficient */ +#define CSMALL 1.e-6 /* Small coefficient */ + #ifdef M_PI #define PI M_PI #else @@ -135,21 +138,21 @@ typedef enum { USE, /* use from previous run */ SAVE, /* save after current run */ SCRATCH -} Hydtype; /* use temporary file */ +} Hydtype; /* use temporary file */ typedef enum { NONE, /* no quality analysis */ CHEM, /* analyze a chemical */ AGE, /* analyze water age */ TRACE -} QualType; /* trace % of flow from a source */ +} QualType; /* trace % of flow from a source */ typedef enum { - V_CURVE, /* volume curve */ - P_CURVE, /* pump curve */ - E_CURVE, /* efficiency curve */ + V_CURVE, /* volume curve */ + P_CURVE, /* pump curve */ + E_CURVE, /* efficiency curve */ H_CURVE -} CurveType; /* head loss curve */ +} CurveType; /* head loss curve */ typedef enum { CONST_HP, /* constant horsepower */ @@ -169,8 +172,8 @@ typedef enum { LOWLEVEL, /* act when grade below set level */ HILEVEL, /* act when grade above set level */ TIMER, /* act when set time reached */ - TIMEOFDAY -} ControlType; /* act when time of day occurs */ + TIMEOFDAY /* act when time of day occurs */ +} ControlType; typedef enum { XHEAD, /* pump cannot deliver head (closed) */ @@ -183,18 +186,18 @@ typedef enum { XPRESSURE, /* valve cannot supply pressure */ FILLING, /* tank filling */ EMPTYING -} StatType; /* tank emptying */ +} StatType; /* tank emptying */ typedef enum { HW, /* Hazen-Williams */ DW, /* Darcy-Weisbach */ - CM -} FormType; /* Chezy-Manning */ + CM /* Chezy-Manning */ +} FormType; typedef enum { US, /* US */ - SI -} UnitsType; /* SI (metric) */ + SI /* SI (metric) */ +} UnitsType; typedef enum { CFS, /* cubic feet per second */ @@ -206,40 +209,41 @@ typedef enum { LPM, /* liters per minute */ MLD, /* megaliters per day */ CMH, /* cubic meters per hour */ - CMD -} FlowUnitsType; /* cubic meters per day */ + CMD /* cubic meters per day */ +} FlowUnitsType; typedef enum { PSI, /* pounds per square inch */ KPA, /* kiloPascals */ - METERS -} PressUnitsType; /* meters */ + METERS /* meters */ +} PressUnitsType; typedef enum { LOW, /* lower limit */ HI, /* upper limit */ - PREC -} RangeType; /* precision */ + PREC /* precision */ +} RangeType; typedef enum { MIX1, /* 1-compartment model */ MIX2, /* 2-compartment model */ FIFO, /* First in, first out model */ - LIFO -} MixType; /* Last in, first out model */ + LIFO /* Last in, first out model */ +} MixType; typedef enum { SERIES, /* none */ AVG, /* time-averages */ MIN, /* minimum values */ MAX, /* maximum values */ - RANGE -} TstatType; /* max - min values */ + RANGE /* max - min values */ +} TstatType; -#define MAXVAR 21 /* Max. # types of network variables */ -/* (equals # items enumed below) */ + +#define MAXVAR 21 /* Max. # types of network variables */ + /* (equals # items enumed below) */ typedef enum { - ELEV = 0, /* nodal elevation */ + ELEV = 0, /* nodal elevation */ DEMAND, /* nodal demand flow */ HEAD, /* nodal hydraulic head */ PRESSURE, /* nodal pressure */ @@ -273,17 +277,27 @@ typedef enum { } SectType; typedef enum { - STATHDR, /* Hydraulic Status */ - ENERHDR, /* Energy Usage */ - NODEHDR, /* Node Results */ - LINKHDR -} HdrType; /* Link Results */ + STATHDR, /* Hydraulic Status header */ + ENERHDR, /* Energy Usage header */ + NODEHDR, /* Node Results header */ + LINKHDR /* Link Results header */ +} HdrType; typedef enum { POSITIVE, NEGATIVE } FlowDirection; +typedef enum { + PCNT_ONLINE, + PCNT_EFFIC, + KWH_PER_FLOW, + TOTAL_KWH, + MAX_KW, + TOTAL_COST, + MAX_ENERGY_STATS +} EnergyStats; + /* ------------------------------------------------------ Global Data Structures @@ -344,6 +358,16 @@ struct Sdemand /* DEMAND CATEGORY OBJECT */ }; typedef struct Sdemand *Pdemand; /* Pointer to demand object */ +typedef struct +{ + double hrsOnLine; // hours pump is online + double efficiency; // total time wtd. efficiency + double kwHrsPerCFS; // total kw-hrs per cfs of flow + double kwHrs; // total kw-hrs consumed + double maxKwatts; // max. kw consumed + double totalCost; // total pumping cost +} Senergy; + struct Ssource /* WQ SOURCE OBJECT */ { /*int Node;*/ /* Node index of source */ @@ -422,13 +446,7 @@ typedef struct /* PUMP OBJECT */ int Upat; /* Utilization pattern index */ int Epat; /* Energy cost pattern index */ double Ecost; /* Unit energy cost */ - double Energy[6]; /* Energy usage statistics: */ - /* 0 = pump utilization */ - /* 1 = avg. efficiency */ - /* 2 = avg. kW/flow */ - /* 3 = avg. kwatts */ - /* 4 = peak kwatts */ - /* 5 = cost/day */ + double Energy[MAX_ENERGY_STATS]; /* Energy usage statistics */ } Spump; typedef struct /* VALVE OBJECT */ @@ -474,7 +492,7 @@ typedef struct /* FIELD OBJECT of report table */ double RptLim[2]; /* Lower/upper report limits */ } SField; -typedef struct s_Premise /* Rule Premise Clause */ +typedef struct s_Premise /* Rule Premise Clause */ { int logop; /* Logical operator */ int object; /* Node or link */ @@ -486,7 +504,7 @@ typedef struct s_Premise /* Rule Premise Clause */ struct s_Premise *next; } Premise; -typedef struct s_Action /* Rule Action Clause */ +typedef struct s_Action /* Rule Action Clause */ { int link; /* Link index */ int status; /* Link's status */ @@ -494,19 +512,19 @@ typedef struct s_Action /* Rule Action Clause */ struct s_Action *next; } Action; -typedef struct s_aRule /* Control Rule Structure */ +typedef struct s_aRule /* Control Rule Structure */ { - char label[MAXID+1]; /* Rule character label */ - double priority; /* Priority level */ - Premise *Pchain; /* Linked list of premises */ - Action *Tchain; /* Linked list of actions if true */ - Action *Fchain; /* Linked list of actions if false */ + char label[MAXID+1]; /* Rule character label */ + double priority; /* Priority level */ + Premise *Pchain; /* Linked list of premises */ + Action *Tchain; /* Linked list of actions if true */ + Action *Fchain; /* Linked list of actions if false */ struct s_aRule *next; } aRule; -typedef struct s_ActItem /* Action list item */ +typedef struct s_ActItem /* Action list item */ { - int ruleindex; /* Index of rule action belongs to */ + int ruleindex; /* Index of rule action belongs to */ struct s_Action *action; /* An action structure */ struct s_ActItem *next; } ActItem; @@ -614,15 +632,15 @@ typedef struct { *Patlist, /* Temporary time pattern list */ *Curvelist; /* Temporary list of curves */ - double *X; // temporary array for curve data + double *X; // temporary array for curve data int - Ntokens, /* Number of tokens in input line */ - Ntitle; /* Number of title lines */ + Ntokens, /* Number of tokens in input line */ + Ntitle; /* Number of title lines */ - char *Tok[MAXTOKS]; /* Array of token strings */ + char *Tok[MAXTOKS]; /* Array of token strings */ char Comment[MAXMSG+1]; - STmplist *PrevPat; /* Pointer to pattern list element */ - STmplist *PrevCurve; /* Pointer to curve list element */ + STmplist *PrevPat; /* Pointer to pattern list element */ + STmplist *PrevCurve; /* Pointer to curve list element */ } parser_data_t; @@ -647,7 +665,7 @@ typedef struct { Rpt1Fname[MAXFNAME+1], /* Primary report file name */ Rpt2Fname[MAXFNAME+1]; /* Secondary report file name */ - SField Field[MAXVAR]; /* Output reporting fields */ + SField Field[MAXVAR]; /* Output reporting fields */ long LineNum; /* Current line number */ long PageNum; /* Current page number */ @@ -694,12 +712,11 @@ typedef struct { aRule *Rule; /* Array of rules */ ActItem *ActList; /* Linked list of action items */ - int RuleState; /* State of rule interpreter */ - long Time1; /* Start of rule evaluation time interval (sec) */ + int RuleState; /* State of rule interpreter */ + long Time1; /* Start of rule evaluation time interval (sec) */ Premise *Plast; /* Previous premise clause */ Action *Tlast; /* Previous true action */ Action *Flast; /* Previous false action */ - } rules_t; /* @@ -754,6 +771,10 @@ typedef struct { Qexp, /* Exponent in orifice formula */ Dmult, /* Demand multiplier */ Hacc, /* Hydraulics solution accuracy */ + + FlowChangeLimit, /* Hydraulics flow change limit */ + HeadErrorLimit, /* Hydraulics head error limit */ + DampLimit, /* Solution damping threshold */ Viscos, /* Kin. viscosity (sq ft/sec) */ SpGrav, /* Specific gravity */ @@ -797,7 +818,7 @@ typedef struct { } hydraulics_t; typedef struct { - int Nnodes, /* Number of network nodes */ + int Nnodes, /* Number of network nodes */ Ntanks, /* Number of tanks */ Njuncs, /* Number of junction nodes */ Nlinks, /* Number of network links */ @@ -810,17 +831,18 @@ typedef struct { Ncurves, /* Number of data curves */ Ncoords; /* Number of Coords */ - Snode *Node; /* Node data */ - Slink *Link; /* Link data */ - Stank *Tank; /* Tank data */ - Spump *Pump; /* Pump data */ - Svalve *Valve; /* Valve data */ - Spattern *Pattern; /* Time patterns */ - Scurve *Curve; /* Curve data */ - Scoord *Coord; /* Coordinate data */ - Scontrol *Control; /* Control data */ - ENHashTable *NodeHashTable, *LinkHashTable; /* Hash tables for ID labels */ - Padjlist *Adjlist; /* Node adjacency lists */ + Snode *Node; /* Node data */ + Slink *Link; /* Link data */ + Stank *Tank; /* Tank data */ + Spump *Pump; /* Pump data */ + Svalve *Valve; /* Valve data */ + Spattern *Pattern; /* Time patterns */ + Scurve *Curve; /* Curve data */ + Scoord *Coord; /* Coordinate data */ + Scontrol *Control; /* Control data */ + ENHashTable *NodeHashTable, + *LinkHashTable; /* Hash tables for ID labels */ + Padjlist *Adjlist; /* Node adjacency lists */ } EN_Network; @@ -852,7 +874,4 @@ struct EN_Project { }; - - - #endif diff --git a/tools/before-test.cmd b/tools/before-test.cmd index 0e8b41a..8871b44 100644 --- a/tools/before-test.cmd +++ b/tools/before-test.cmd @@ -25,8 +25,8 @@ set SCRIPT_HOME=%~dp0 set TEST_HOME=%~1 -set EXAMPLES_VER=1.0.1 -set BENCHMARK_VER=2012vs10 +set EXAMPLES_VER=1.0.2-dev +set BENCHMARK_VER=220dev-vs17 set TESTFILES_URL=https://github.com/OpenWaterAnalytics/epanet-example-networks/archive/v%EXAMPLES_VER%.zip diff --git a/tools/run-nrtest.cmd b/tools/run-nrtest.cmd index 0b7c17f..f5d4c8d 100644 --- a/tools/run-nrtest.cmd +++ b/tools/run-nrtest.cmd @@ -18,7 +18,7 @@ setlocal set NRTEST_SCRIPT_PATH=%~1 set TEST_SUITE_PATH=%~2 -set BENCHMARK_VER=2012vs10 +set BENCHMARK_VER=220dev-vs17 set NRTEST_EXECUTE_CMD=python %NRTEST_SCRIPT_PATH%\nrtest execute diff --git a/win_build/WinSDK/Makefile.bat b/win_build/WinSDK/Makefile.bat index eeb967e..d0a00c3 100644 --- a/win_build/WinSDK/Makefile.bat +++ b/win_build/WinSDK/Makefile.bat @@ -21,9 +21,9 @@ Find /i "x86" < checkOS.tmp > StringCheck.tmp If %ERRORLEVEL% == 1 ( CALL "%SDK_PATH%bin\"SetEnv.cmd /x64 /release rem : create EPANET2.DLL - cl -o epanet2.dll epanet.c hash.c hydraul.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /link /DLL + cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /link /DLL rem : create EPANET2.EXE - cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /I ..\src /link + cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /I ..\src /link md "%Build_PATH%"\64bit move /y "%SRC_PATH%"\*.dll "%Build_PATH%"\64bit move /y "%SRC_PATH%"\*.exe "%Build_PATH%"\64bit @@ -35,9 +35,9 @@ rem : 32 bit with DEF CALL "%SDK_PATH%bin\"SetEnv.cmd /x86 /release echo "32 bit with epanet2.def mapping" rem : create EPANET2.DLL -cl -o epanet2.dll epanet.c hash.c hydraul.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /link /DLL /def:..\win_build\WinSDK\epanet2.def /MAP +cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /link /DLL /def:..\win_build\WinSDK\epanet2.def /MAP rem : create EPANET2.EXE -cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /I ..\src /link +cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /I ..\src /link md "%Build_PATH%"\32bit move /y "%SRC_PATH%"\*.dll "%Build_PATH%"\32bit move /y "%SRC_PATH%"\*.exe "%Build_PATH%"\32bit