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

Resolving merge conflicts

Browse Source
This commit is contained in:
Michael Tryby
2018-08-27 09:38:46 -04:00
parent 40cd25748c c4b7c90634
commit bdc05058ca
26 changed files with 1572 additions and 1059 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/enumstxt.h
View File

@@ -74,6 +74,10 @@ char *PressUnitsTxt[] = {w_PSI,
w_KPA,
w_METERS};
char *DemandModelTxt[] = { w_DDA,
w_PDA,
NULL };
char *QualTxt[] = {w_NONE,
w_CHEM,
w_AGE,

147
src/epanet.c
View File

@@ -129,9 +129,10 @@ execute function x and set the error code equal to its return value.
#include "funcs.h"
#include "text.h"
#include "types.h"
#define EXTERN
////////////////////////////////////////////#include "epanet2.h"
#include "vars.h"
// This single global variable is used only when the library is called
// in "legacy mode" with the 2.1-style API.
EN_Project *_defaultModel;
// Local functions
@@ -164,18 +165,26 @@ int DLLEXPORT ENepanet(const char *f1, const char *f2, const char *f3, void (*pv
{
int errcode = 0;
ERRCODE(EN_alloc(&_defaultModel));
ERRCODE(EN_epanet(_defaultModel, f1, f2, f3, pviewprog));
ERRCODE(EN_free(&_defaultModel));
ERRCODE(EN_createproject(&_defaultModel));
ERRCODE(EN_open(_defaultModel, f1, f2, f3));
_defaultModel->viewprog = pviewprog;
return errcode;
if (_defaultModel->out_files.Hydflag != USE) {
ERRCODE(EN_solveH(_defaultModel));
}
ERRCODE(EN_solveQ(_defaultModel));
ERRCODE(EN_report(_defaultModel));
EN_close(_defaultModel);
EN_deleteproject(&_defaultModel);
return (errcode);
}
int DLLEXPORT ENopen(char *f1, char *f2, char *f3) {
int errcode = 0;
ERRCODE(EN_alloc(&_defaultModel));
ERRCODE(EN_createproject(&_defaultModel));
EN_open(_defaultModel, f1, f2, f3);
return (errcode);
}
@@ -263,6 +272,16 @@ int DLLEXPORT ENsetflowunits(int code) {
return EN_setflowunits(_defaultModel, code);
}
int DLLEXPORT ENgetdemandmodel(int *type, EN_API_FLOAT_TYPE *pmin, EN_API_FLOAT_TYPE *preq,
EN_API_FLOAT_TYPE *pexp) {
return EN_getdemandmodel(_defaultModel, type, pmin, preq, pexp);
}
int DLLEXPORT ENsetdemandmodel(int type, EN_API_FLOAT_TYPE pmin, EN_API_FLOAT_TYPE preq,
EN_API_FLOAT_TYPE pexp) {
return EN_setdemandmodel(_defaultModel, type, pmin, preq, pexp);
}
int DLLEXPORT ENgetpatternindex(char *id, int *index) {
return EN_getpatternindex(_defaultModel, id, index);
}
@@ -537,6 +556,12 @@ int DLLEXPORT ENdeletenode(int index) {
return EN_deletenode(_defaultModel, index);
}
// int DLLEXPORT EN_epanet(char *inpFile, char *rptFile, char *binOutFile, void(*callback) (char *))
// {
// return ENepanet(inpFile, rptFile, binOutFile, callback);
// }
/*
----------------------------------------------------------------
Functions for opening & closing the EPANET system
@@ -544,7 +569,7 @@ int DLLEXPORT ENdeletenode(int index) {
*/
/// allocate a project pointer
int DLLEXPORT EN_alloc(EN_ProjectHandle *ph)
int DLLEXPORT EN_createproject(EN_ProjectHandle *ph)
{
int errorcode = 0;
EN_Project *project = calloc(1, sizeof(EN_Project));
@@ -559,7 +584,7 @@ int DLLEXPORT EN_alloc(EN_ProjectHandle *ph)
return errorcode;
}
int DLLEXPORT EN_free(EN_ProjectHandle *ph)
int DLLEXPORT EN_deleteproject(EN_ProjectHandle *ph)
{
int errorcode = 0;
EN_Project *p = (EN_Project*)(*ph);
@@ -577,25 +602,43 @@ int DLLEXPORT EN_free(EN_ProjectHandle *ph)
return 0;
}
int DLLEXPORT EN_epanet(EN_ProjectHandle ph, const char *f1, const char *f2,
const char *f3, void(*pviewprog)(char *))
{
int errcode = 0;
EN_Project *_p = (EN_Project*)ph;
// int DLLEXPORT EN_epanet(EN_ProjectHandle ph, const char *f1, const char *f2,
// const char *f3, void(*pviewprog)(char *))
// {
// int errcode = 0;
// EN_Project *_p = (EN_Project*)ph;
ERRCODE(EN_open(ph, f1, f2, f3));
// ERRCODE(EN_open(ph, f1, f2, f3));
_p->viewprog = pviewprog;
if (_p->out_files.Hydflag != USE) {
ERRCODE(EN_solveH(ph));
}
// _p->viewprog = pviewprog;
// if (_p->out_files.Hydflag != USE) {
// ERRCODE(EN_solveH(ph));
// }
ERRCODE(EN_solveQ(ph));
ERRCODE(EN_report(ph));
EN_close(ph);
// ERRCODE(EN_solveQ(ph));
// ERRCODE(EN_report(ph));
// EN_close(ph);
return set_error(_p->error_handle, errcode);
}
// return set_error(_p->error_handle, errcode);
// =======
/// Create an EPANET project
// int DLLEXPORT EN_createproject(EN_Project **p)
// {
// EN_Project *project = calloc(1, sizeof(EN_Project));
// if (project == NULL) return 101;
// *p = project;
// return 0;
// }
// /// Delete an EPANET project
// int DLLEXPORT EN_deleteproject(EN_Project *p)
// {
// if (p) free(p);
// p = NULL;
// return 0;
// >>>>>>> c4b7c90634c73c1f800b2a18ea01f4be5b6711c1
// }
int DLLEXPORT EN_init(EN_ProjectHandle *ph, char *f2, char *f3,
EN_FlowUnits UnitsType, EN_FormType HeadlossFormula)
@@ -603,9 +646,8 @@ int DLLEXPORT EN_init(EN_ProjectHandle *ph, char *f2, char *f3,
** Input:
** f2 = pointer to name of report file
** f3 = pointer to name of binary output file
UnitsType = flow units flag
HeadlossFormula = headloss formula flag
** UnitsType = flow units flag
** HeadlossFormula = headloss formula flag
** Output: none
** Returns: error code
** Purpose: opens EPANET
@@ -1622,7 +1664,39 @@ int DLLEXPORT EN_setflowunits(EN_ProjectHandle ph, int code) {
return(0);
}
// <<<<<<< HEAD
// int DLLEXPORT EN_getpatternindex(EN_ProjectHandle ph, char *id, int *index) {
// =======
int DLLEXPORT EN_getdemandmodel(EN_ProjectHandle ph, int *type, EN_API_FLOAT_TYPE *pmin,
EN_API_FLOAT_TYPE *preq, EN_API_FLOAT_TYPE *pexp)
{
EN_Project *p = (EN_Project*)ph;
*type = p->hydraulics.DemandModel;
*pmin = (EN_API_FLOAT_TYPE)(p->hydraulics.Pmin * p->Ucf[PRESSURE]);
*preq = (EN_API_FLOAT_TYPE)(p->hydraulics.Preq * p->Ucf[PRESSURE]);
*pexp = (EN_API_FLOAT_TYPE)(p->hydraulics.Pexp);
return 0;
}
int DLLEXPORT EN_setdemandmodel(EN_ProjectHandle ph, int type, EN_API_FLOAT_TYPE pmin,
EN_API_FLOAT_TYPE preq, EN_API_FLOAT_TYPE pexp)
{
EN_Project *p = (EN_Project*)ph;
if (type < 0 || type > EN_PDA) return 251;
if (pmin > preq || pexp <= 0.0) return 202;
p->hydraulics.DemandModel = type;
p->hydraulics.Pmin = pmin / p->Ucf[PRESSURE];
p->hydraulics.Preq = preq / p->Ucf[PRESSURE];
p->hydraulics.Pexp = pexp;
return 0;
}
int DLLEXPORT EN_getpatternindex(EN_ProjectHandle ph, char *id, int *index) {
//>>>>>>> c4b7c90634c73c1f800b2a18ea01f4be5b6711c1
int i;
EN_Project *p = (EN_Project*)ph;
@@ -1868,11 +1942,17 @@ int DLLEXPORT EN_getstatistic(EN_ProjectHandle ph, int code, EN_API_FLOAT_TYPE *
switch (code) {
case EN_ITERATIONS:
*value = (EN_API_FLOAT_TYPE)p->hydraulics.iterations;
*value = (EN_API_FLOAT_TYPE)p->hydraulics.Iterations;
break;
case EN_RELATIVEERROR:
*value = (EN_API_FLOAT_TYPE)p->hydraulics.relativeError;
*value = (EN_API_FLOAT_TYPE)p->hydraulics.RelativeError;
break;
case EN_MAXHEADERROR:
*value = (EN_API_FLOAT_TYPE)(p->hydraulics.MaxHeadError * p->Ucf[HEAD]);
break;
case EN_MAXFLOWCHANGE:
*value = (EN_API_FLOAT_TYPE)(p->hydraulics.MaxFlowChange * p->Ucf[FLOW]);
break;
default:
break;
}
@@ -1950,7 +2030,10 @@ int DLLEXPORT EN_getcoord(EN_ProjectHandle ph, int index, EN_API_FLOAT_TYPE *x,
*x = (EN_API_FLOAT_TYPE)p->network.Coord[index].X;
*y = (EN_API_FLOAT_TYPE)p->network.Coord[index].Y;
// <<<<<<< HEAD
// =======
// >>>>>>> c4b7c90634c73c1f800b2a18ea01f4be5b6711c1
return 0;
}

22
src/funcs.h
View File

@@ -161,28 +161,13 @@ void inithyd(EN_Project *pr, int initFlags); /* Re-sets initial condition
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 */
int allocmatrix(EN_Project *pr); /* Allocates matrix coeffs. */
void freematrix(EN_Project *pr); /* Frees matrix coeffs. */
void initlinkflow(EN_Project *pr, int, char,
double); /* Initializes link flow */
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 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. */
@@ -191,8 +176,13 @@ int hydsolve(EN_Project *pr, int *,double *); /* Solves network equations
/* ----------- 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 headlosscoeffs(EN_Project *pr); // Finds link head loss coeffs.
void matrixcoeffs(EN_Project *pr); /* Finds hyd. matrix coeffs. */
double emitflowchange(EN_Project *pr, int i); /* Change in emitter outflow */
double demandflowchange(EN_Project *pr, int i, // Change in demand outflow
double dp, double n);
void demandparams(EN_Project *pr, double *dp, // PDA function parameters
double *n);
/* ----------- SMATRIX.C ---------------*/
int createsparse(EN_Project *pr); /* Creates sparse matrix */

360
src/hydcoeffs.c
View File

@@ -29,15 +29,21 @@ 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);
//void resistcoeff(EN_Project *pr, int k);
//void headlosscoeffs(EN_Project *pr);
//void matrixcoeffs(EN_Project *pr);
//double emitflowchange(EN_Project *pr, int i);
//double demandflowchange(EN_Project *pr, int i, double dp, double n);
//void demandparams(EN_Project *pr, double *dp, double *n);
// 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 demandcoeffs(EN_Project *pr);
static void demandheadloss(double d, double dfull, double dp,
double n, double *hloss, double *hgrad);
static void pipecoeff(EN_Project *pr, int k);
static void DWpipecoeff(EN_Project *pr, int k);
@@ -55,7 +61,6 @@ 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)
/*
**--------------------------------------------------------------------
@@ -83,13 +88,14 @@ void resistcoeff(EN_Project *pr, int k)
switch (hyd->Formflag)
{
case HW:
link->R = 4.727*L / pow(e, hyd->Hexp) / pow(d, 4.871);
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);
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;
link->R = SQR(4.0 * e / (1.49 * PI * SQR(d))) *
pow((d / 4.0), -1.333) * L;
}
break;
@@ -106,45 +112,46 @@ void resistcoeff(EN_Project *pr, int k)
}
void hlosscoeff(EN_Project *pr, int k)
void headlosscoeffs(EN_Project *pr)
/*
**--------------------------------------------------------------
** Input: k = link index
** Input: none
** Output: none
** Purpose: computes P and Y coefficients for a link
** Purpose: computes coefficients P (1 / head loss gradient)
** and Y (head loss / gradient) for all links.
**--------------------------------------------------------------
*/
{
int k;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
solver_t *sol = &hyd->solver;
Slink *link = &net->Link[k];
switch (link->Type)
for (k = 1; k <= net->Nlinks; k++)
{
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);
switch (net->Link[k].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 hyd->solver.P[k] = 0.0;
}
else sol->P[k] = 0.0;
}
}
@@ -162,13 +169,23 @@ void matrixcoeffs(EN_Project *pr)
solver_t *sol = &hyd->solver;
EN_Network *net = &pr->network;
// Reset values of all diagonal coeffs. (Aii), off-diagonal
// coeffs. (Aij), r.h.s. coeffs. (F) and node flow balance (X_tmp)
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));
// Compute matrix coeffs. from links, emitters, and nodal demands
linkcoeffs(pr);
emittercoeffs(pr);
demandcoeffs(pr);
// Update nodal flow balances with demands and add onto r.h.s. coeffs.
nodecoeffs(pr);
// Finally, find coeffs. for PRV/PSV/FCV control valves whose
// status is not fixed to OPEN/CLOSED
valvecoeffs(pr);
}
@@ -189,7 +206,7 @@ void linkcoeffs(EN_Project *pr)
solver_t *sol = &hyd->solver;
Slink *link;
// Examine each link of network */
// Examine each link of network
for (k = 1; k <= net->Nlinks; k++)
{
if (sol->P[k] == 0.0) continue;
@@ -197,7 +214,7 @@ void linkcoeffs(EN_Project *pr)
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)
// Update net nodal inflows (X_tmp), 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];
@@ -212,21 +229,18 @@ void linkcoeffs(EN_Project *pr)
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 n1 is a tank/reservoir
else sol->F[sol->Row[n2]] += (sol->P[k] * hyd->NodeHead[n1]);
// Node n2 is junction
if (n2 <= net->Njuncs) {
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]);
}
// Node n2 is a tank/reservoir
else sol->F[sol->Row[n1]] += (sol->P[k] * hyd->NodeHead[n2]);
}
}
@@ -236,7 +250,7 @@ void nodecoeffs(EN_Project *pr)
**----------------------------------------------------------------
** Input: none
** Output: none
** Purpose: completes calculation of nodal flow imbalance (X)
** Purpose: completes calculation of nodal flow imbalance (X_tmp)
** & flow correction (F) arrays
**----------------------------------------------------------------
*/
@@ -250,7 +264,7 @@ void nodecoeffs(EN_Project *pr)
// flow imbalance & add imbalance to RHS array F.
for (i = 1; i <= net->Njuncs; i++)
{
hyd->X_tmp[i] -= hyd->NodeDemand[i];
hyd->X_tmp[i] -= hyd->DemandFlows[i];
sol->F[sol->Row[i]] += hyd->X_tmp[i];
}
}
@@ -276,15 +290,20 @@ void valvecoeffs(EN_Project *pr)
// Examine each valve
for (i = 1; i <= net->Nvalves; i++)
{
// Find valve's link index
valve = &net->Valve[i];
k = valve->Link; // Link index of valve
k = valve->Link;
// Coeffs. for fixed status valves have already been computed
if (hyd->LinkSetting[k] == MISSING) continue;
// Start & end nodes of valve's link
link = &net->Link[k];
if (hyd->LinkSetting[k] == MISSING) {
continue; // Valve status fixed
}
n1 = link->N1; // Start & end nodes
n1 = link->N1;
n2 = link->N2;
switch (link->Type) // Call valve-specific function
// Call valve-specific function
switch (link->Type)
{
case EN_PRV:
prvcoeff(pr, k, n1, n2);
@@ -330,17 +349,17 @@ void emittercoeffs(EN_Project *pr)
for (i = 1; i <= net->Njuncs; i++)
{
node = &net->Node[i];
if (node->Ke == 0.0) {
continue;
}
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) {
if (p < hyd->RQtol)
{
p = 1.0 / hyd->RQtol;
}
else {
else
{
p = 1.0 / p; // inverse head loss gradient
}
y = SGN(q)*z*p; // head loss / gradient
@@ -351,6 +370,173 @@ void emittercoeffs(EN_Project *pr)
}
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;
Snode *node = &pr->network.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 demandparams(EN_Project *pr, double *dp, double *n)
/*
**--------------------------------------------------------------
** Input: none
** Output: dp = pressure range over which demands can vary
** n = exponent in head loss v. demand function
** Purpose: retrieves parameters that define a pressure dependent
** demand function.
**--------------------------------------------------------------
*/
{
hydraulics_t *hyd = &pr->hydraulics;
// If required pressure equals minimum pressure, use a linear demand
// curve with a 0.01 PSI pressure range to approximate an all or
// nothing demand solution
if (hyd->Preq == hyd->Pmin)
{
*dp = 0.01 / PSIperFT;
*n = 1.0;
}
// Otherwise use the user-supplied demand curve parameters
else
{
*dp = hyd->Preq - hyd->Pmin;
*n = 1.0 / hyd->Pexp;
}
}
void demandcoeffs(EN_Project *pr)
/*
**--------------------------------------------------------------
** Input: none
** Output: none
** Purpose: computes matrix coeffs. for pressure dependent demands
**
** Note: Pressure dependent demands are modelled like emitters
** with Hloss = Pserv * (D / Dfull)^(1/Pexp)
** where D (actual demand) is zero for negative pressure
** and is Dfull above pressure Pserv.
**--------------------------------------------------------------
*/
{
int i, row;
double dp, // pressure range over which demand can vary (ft)
n, // exponent in head loss v. demand function
hloss, // head loss in supplying demand (ft)
hgrad; // gradient of demand head loss (ft/cfs)
hydraulics_t *hyd = &pr->hydraulics;
solver_t *sol = &hyd->solver;
EN_Network *net = &pr->network;
// Get demand function parameters
if (hyd->DemandModel == DDA) return;
demandparams(pr, &dp, &n);
// Examine each junction node
for (i = 1; i <= net->Njuncs; i++)
{
// Skip junctions with non-positive demands
if (hyd->NodeDemand[i] <= 0.0) continue;
// Find head loss for demand outflow at node's elevation
demandheadloss(hyd->DemandFlows[i], hyd->NodeDemand[i], dp, n,
&hloss, &hgrad);
// Update row of solution matrix A & its r.h.s. F
row = sol->Row[i];
sol->Aii[row] += 1.0 / hgrad;
sol->F[row] += (hloss + net->Node[i].El + hyd->Pmin) / hgrad;
}
}
double demandflowchange(EN_Project *pr, int i, double dp, double n)
/*
**--------------------------------------------------------------
** Input: i = node index
** dp = pressure range fro demand funtion (ft)
** n = exponent in head v. demand function
** Output: returns change in pressure dependent demand flow
** Purpose: computes change in outflow at at a node subject to
** pressure dependent demands
**--------------------------------------------------------------
*/
{
double hloss, hgrad;
hydraulics_t *hyd = &pr->hydraulics;
demandheadloss(hyd->DemandFlows[i], hyd->NodeDemand[i], dp, n, &hloss, &hgrad);
return (hloss - hyd->NodeHead[i] + pr->network.Node[i].El + hyd->Pmin) / hgrad;
}
void demandheadloss(double d, double dfull, double dp, double n,
double *hloss, double *hgrad)
/*
**--------------------------------------------------------------
** Input: d = actual junction demand (cfs)
** dfull = full junction demand required (cfs)
** dp = pressure range for demand function (ft)
** n = exponent in head v. demand function
** Output: hloss = head loss delivering demand d (ft)
** hgrad = gradient of head loss (ft/cfs)
** Purpose: computes head loss and its gradient for delivering
** a pressure dependent demand flow.
**--------------------------------------------------------------
*/
{
const double RB = 1.0e9;
const double EPS = 0.001;
double r = d / dfull;
// Use upper barrier function for excess demand above full value
if (r > 1.0)
{
*hgrad = RB;
*hloss = dp + RB * (d - dfull);
}
// Use lower barrier function for negative demand
else if (r < 0)
{
*hgrad = RB;
*hloss = RB * d;
}
// Use linear head loss function for near zero demand
else if (r < EPS)
{
*hgrad = dp * pow(EPS, n - 1.0) / dfull;
*hloss = (*hgrad) * d;
}
// Otherwise use power head loss function
else
{
*hgrad = n * dp * pow(r, n - 1.0) / dfull;
*hloss = (*hgrad) * d / n;
}
}
void pipecoeff(EN_Project *pr, int k)
/*
**--------------------------------------------------------------
@@ -382,8 +568,9 @@ void pipecoeff(EN_Project *pr, int k)
return;
}
// ... head loss formula is Darcy-Weisbach
if (hyd->Formflag == DW) {
// Use custom function for Darcy-Weisbach formula
if (hyd->Formflag == DW)
{
DWpipecoeff(pr, k);
return;
}
@@ -550,7 +737,8 @@ void pumpcoeff(EN_Project *pr, int k)
Spump *pump;
// Use high resistance pipe if pump closed or cannot deliver head
if (hyd->LinkStatus[k] <= CLOSED || setting == 0.0) {
if (hyd->LinkStatus[k] <= CLOSED || setting == 0.0)
{
sol->P[k] = 1.0 / CBIG;
sol->Y[k] = hyd->LinkFlows[k];
return;
@@ -580,9 +768,7 @@ void pumpcoeff(EN_Project *pr, int k)
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);
}
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);
@@ -615,15 +801,9 @@ void curvecoeff(EN_Project *pr, int i, double q, double *h0, double *r)
// 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--;
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
@@ -653,13 +833,12 @@ void gpvcoeff(EN_Project *pr, int k)
solver_t *sol = &hyd->solver;
// Treat as a pipe if valve closed
if (hyd->LinkStatus[k] == CLOSED) {
valvecoeff(pr, k);
}
if (hyd->LinkStatus[k] == CLOSED) valvecoeff(pr, k);
// Otherwise utilize headloss curve
// whose index is stored in K
else {
else
{
// Find slope & intercept of headloss curve.
q = ABS(hyd->LinkFlows[k]);
q = MAX(q, TINY);
@@ -687,18 +866,22 @@ void pbvcoeff(EN_Project *pr, int k)
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) {
if (hyd->LinkSetting[k] == MISSING || hyd->LinkSetting[k] == 0.0)
{
valvecoeff(pr, k);
}
// If valve is active
else {
else
{
// Treat as a pipe if minor loss > valve setting
if (link->Km * SQR(hyd->LinkFlows[k]) > hyd->LinkSetting[k]) {
if (link->Km * SQR(hyd->LinkFlows[k]) > hyd->LinkSetting[k])
{
valvecoeff(pr, k);
}
// Otherwise force headloss across valve to be equal to setting
else {
else
{
sol->P[k] = CBIG;
sol->Y[k] = hyd->LinkSetting[k] * CBIG;
}
@@ -723,7 +906,8 @@ void tcvcoeff(EN_Project *pr, int k)
km = link->Km;
// If valve not fixed OPEN or CLOSED, compute its loss coeff.
if (hyd->LinkSetting[k] != MISSING) {
if (hyd->LinkSetting[k] != MISSING)
{
link->Km = 0.02517 * hyd->LinkSetting[k] / (SQR(link->Diam)*SQR(link->Diam));
}
@@ -768,7 +952,8 @@ void prvcoeff(EN_Project *pr, int k, int n1, int n2)
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) {
if (hyd->X_tmp[n2] < 0.0)
{
sol->F[i] += hyd->X_tmp[n2];
}
return;
@@ -818,7 +1003,8 @@ void psvcoeff(EN_Project *pr, int k, int n1, int n2)
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) {
if (hyd->X_tmp[n1] > 0.0)
{
sol->F[j] += hyd->X_tmp[n1];
}
return;
@@ -919,9 +1105,7 @@ void valvecoeff(EN_Project *pr, int k)
if (link->Km > 0.0)
{
p = 2.0 * link->Km * fabs(flow);
if (p < hyd->RQtol) {
p = hyd->RQtol;
}
if (p < hyd->RQtol) p = hyd->RQtol;
sol->P[k] = 1.0 / p;
sol->Y[k] = flow / 2.0;
}

26
src/hydraul.c
View File

@@ -60,6 +60,20 @@ AUTHOR: L. Rossman
#define QZERO 1.e-6 /* Equivalent to zero flow */
// Local functions
int allocmatrix(EN_Project *pr);
void freematrix(EN_Project *pr);
void initlinkflow(EN_Project *pr, int, char, double);
void setlinkflow(EN_Project *pr, int, double);
void demands(EN_Project *pr);
int controls(EN_Project *pr);
long timestep(EN_Project *pr);
void controltimestep(EN_Project *pr, long *);
void ruletimestep(EN_Project *pr, long *);
void addenergy(EN_Project *pr, long);
void tanklevels(EN_Project *pr, long);
int openhyd(EN_Project *pr)
/*
*--------------------------------------------------------------
@@ -194,15 +208,12 @@ int runhyd(EN_Project *pr, long *t)
/* 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;
/*** Updated 3/1/01 ***/
/* If system unbalanced and no extra trials */
@@ -316,6 +327,7 @@ int allocmatrix(EN_Project *pr)
s->Aii = (double *) calloc(net->Nnodes+1,sizeof(double));
s->Aij = (double *) calloc(hyd->Ncoeffs+1,sizeof(double));
s->F = (double *) calloc(net->Nnodes+1,sizeof(double));
hyd->DemandFlows = (double *)calloc(net->Nnodes + 1, sizeof(double));
hyd->EmitterFlows = (double *) calloc(net->Nnodes+1,sizeof(double));
s->P = (double *) calloc(net->Nlinks+1,sizeof(double));
s->Y = (double *) calloc(net->Nlinks+1,sizeof(double));
@@ -324,6 +336,7 @@ int allocmatrix(EN_Project *pr)
ERRCODE(MEMCHECK(s->Aii));
ERRCODE(MEMCHECK(s->Aij));
ERRCODE(MEMCHECK(s->F));
ERRCODE(MEMCHECK(hyd->DemandFlows));
ERRCODE(MEMCHECK(hyd->EmitterFlows));
ERRCODE(MEMCHECK(s->P));
ERRCODE(MEMCHECK(s->Y));
@@ -348,6 +361,7 @@ void freematrix(EN_Project *pr)
free(s->Aii);
free(s->Aij);
free(s->F);
free(hyd->DemandFlows);
free(hyd->EmitterFlows);
free(s->P);
free(s->Y);
@@ -551,7 +565,6 @@ void setlinksetting(EN_Project *pr, int index, double value, StatType *s, doubl
}
void demands(EN_Project *pr)
/*
**--------------------------------------------------------------------
@@ -591,6 +604,9 @@ void demands(EN_Project *pr)
sum += djunc;
}
hyd->NodeDemand[i] = sum;
// Initialize pressure dependent demand
hyd->DemandFlows[i] = sum;
}
/* Update head at fixed grade nodes with time patterns. */

1007
src/hydsolver.c
View File

File diff suppressed because it is too large Load Diff

462
src/hydstatus.c Normal file
View File

@@ -0,0 +1,462 @@
/*
*********************************************************************
HYDSTATUS.C -- Hydraulic status updating for the EPANET Program
*******************************************************************
*/
#include <stdio.h>
#include "types.h"
#include "funcs.h"
// External functions
int valvestatus(EN_Project *pr);
int linkstatus(EN_Project *pr);
// Local functions
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);
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;
// Examine each valve
for (i = 1; i <= net->Nvalves; i++)
{
// Get valve's link and its index
k = net->Valve[i].Link;
link = &net->Link[k];
// Ignore valve if its status is fixed to OPEN/CLOSED
if (hyd->LinkSetting[k] == MISSING) continue;
// Get start/end node indexes & save current status
n1 = link->N1;
n2 = link->N2;
status = hyd->LinkStatus[k];
// Evaluate valve's new status
switch (link->Type)
{
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;
default:
continue;
}
// Check for a 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() */
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 across link
StatType status; // Current status
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
report_options_t *rep = &pr->report;
Slink *link;
// Examine each link
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 any 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]);
}
}
}
return change;
}
StatType cvstatus(EN_Project *pr, StatType s, double dh, double q)
/*
**--------------------------------------------------
** Input: s = current link status
** dh = head loss across link
** q = link flow
** Output: returns new link status
** Purpose: updates status of a check valve link.
**--------------------------------------------------
*/
{
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 across link
** 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;
// Find maximum head (hmax) pump can deliver
p = findpump(net, k);
if (net->Pump[p].Ptype == CONST_HP)
{
// Use huge value for constant HP pump
hmax = BIG;
}
else
{
// Use speed-adjusted shut-off head for other pumps
hmax = SQR(hyd->LinkSetting[k]) * net->Pump[p].Hmax;
}
// Check if currrent head gain exceeds pump's max. head
if (dh > hmax + hyd->Htol) return XHEAD;
// No check is made to see if flow exceeds pump's max. flow
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; // Valve's new status
double hml; // Head loss when fully opened
hydraulics_t *hyd = &pr->hydraulics;
double htol = hyd->Htol;
Slink *link = &pr->network.Link[k];
// Head loss when fully open
hml = link->Km * SQR(hyd->LinkFlows[k]);
// Rules for updating valve's status from current value s
status = s;
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; // Valve's new status
double hml; // Head loss when fully opened
hydraulics_t *hyd = &pr->hydraulics;
double htol = hyd->Htol;
Slink *link = &pr->network.Link[k];
// Head loss when fully open
hml = link->Km * SQR(hyd->LinkFlows[k]);
// Rules for updating valve's status from current value s
status = s;
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;
}
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];
// Return if link is closed
if (hyd->LinkStatus[k] <= CLOSED) return;
// Make node n1 be the tank, reversing flow (q) if need be
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;
}
// Ignore reservoirs
tank = &net->Tank[i];
if (tank->A == 0.0) return;
// Find head difference across link
h = hyd->NodeHead[n1] - hyd->NodeHead[n2];
// If tank is 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
// (e.g., an open outflow check valve would close)
else if (cvstatus(pr, OPEN, h, q) == CLOSED)
{
hyd->LinkStatus[k] = TEMPCLOSED;
}
}
// If tank is 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
// (e.g., a closed outflow check valve would open)
else if (cvstatus(pr, CLOSED, h, q) == OPEN)
{
hyd->LinkStatus[k] = TEMPCLOSED;
}
}
}

13
src/inpfile.c
View File

@@ -26,14 +26,14 @@ data describing a piping network to a file in EPANET's text format.
#else
#include <stdlib.h>
#endif
#include "hash.h"
#include "text.h"
#include "types.h"
#include "epanet2.h"
#include "funcs.h"
#include <math.h>
#define EXTERN extern
#include "vars.h"
//#define EXTERN extern
//#include "vars.h"
/* Defined in enumstxt.h in EPANET.C */
extern char *LinkTxt[];
@@ -624,6 +624,13 @@ int saveinpfile(EN_Project *pr, char *fname)
if (hyd->FlowChangeLimit > 0.0) {
fprintf(f, "\n FLOWCHANGE %-.8f", hyd->FlowChangeLimit * pr->Ucf[FLOW]);
}
if (hyd->DemandModel == PDA)
{
fprintf(f, "\n DEMAND MODEL PDA");
fprintf(f, "\n MINIMUM PRESSURE %-.4f", hyd->Pmin * pr->Ucf[PRESSURE]);
fprintf(f, "\n REQUIRED PRESSURE %-.4f", hyd->Preq * pr->Ucf[PRESSURE]);
fprintf(f, "\n PRESSURE EXPONENT %-.4f", hyd->Pexp);
}
/* Write [REPORT] section */

12
src/input1.c
View File

@@ -33,8 +33,6 @@ AUTHOR: L. Rossman
#include "epanet2.h"
#include "funcs.h"
#include <math.h>
#define EXTERN extern
#include "vars.h"
/*
--------------------- Module Global Variables ----------------------
@@ -128,8 +126,12 @@ void setdefaults(EN_Project *pr)
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 */
hyd->FlowChangeLimit = 0.0; // Default flow change limit
hyd->HeadErrorLimit = 0.0; // Default head error limit
hyd->DemandModel = DDA; // Demand driven analysis
hyd->Pmin = 0.0; // Minimum demand pressure (ft)
hyd->Preq = 0.0; // Required demand pressure (ft)
hyd->Pexp = 0.5; // Pressure function exponent
qu->Ctol = MISSING; /* No pre-set quality tolerance */
hyd->MaxIter = MAXITER; /* Default max. hydraulic trials */
@@ -581,6 +583,8 @@ void convertunits(EN_Project *pr)
demand->Base /= pr->Ucf[DEMAND];
}
}
hyd->Pmin /= pr->Ucf[PRESSURE];
hyd->Preq /= pr->Ucf[PRESSURE];
/* Convert emitter discharge coeffs. to head loss coeff. */
ucf = pow(pr->Ucf[FLOW], hyd->Qexp) / pr->Ucf[PRESSURE];

2
src/input2.c
View File

@@ -37,8 +37,6 @@ The following utility functions are all called from INPUT3.C
#include "epanet2.h"
#include "funcs.h"
#include <math.h>
#define EXTERN extern
#include "vars.h"
#define MAXERRS 10 /* Max. input errors reported */

301
src/input3.c
View File

@@ -32,12 +32,11 @@ All functions in this module are called from newline() in INPUT2.C.
#include "text.h"
#include "types.h"
#include <math.h>
#define EXTERN extern
#include "vars.h"
/* Defined in enumstxt.h in EPANET.C */
extern char *MixTxt[];
extern char *Fldname[];
extern char *DemandModelTxt[];
/* Defined in INPUT2.C */
@@ -1763,6 +1762,7 @@ int optionchoice(EN_Project *pr, int n)
** VERIFY filename
** UNBALANCED STOP/CONTINUE {Niter}
** PATTERN id
** DEMAND MODEL DDA/PDA/PPA
**--------------------------------------------------------------
*/
{
@@ -1771,131 +1771,125 @@ int optionchoice(EN_Project *pr, int n)
quality_t *qu = &pr->quality;
parser_data_t *par = &pr->parser;
out_file_t *out = &pr->out_files;
int choice;
/* Check if 1st token matches a parameter name and */
/* process the input for the matched parameter */
if (n < 0)
return (201);
if (match(par->Tok[0], w_UNITS)) {
if (n < 1)
return (0);
else if (match(par->Tok[1], w_CFS))
par->Flowflag = CFS;
else if (match(par->Tok[1], w_GPM))
par->Flowflag = GPM;
else if (match(par->Tok[1], w_AFD))
par->Flowflag = AFD;
else if (match(par->Tok[1], w_MGD))
par->Flowflag = MGD;
else if (match(par->Tok[1], w_IMGD))
par->Flowflag = IMGD;
else if (match(par->Tok[1], w_LPS))
par->Flowflag = LPS;
else if (match(par->Tok[1], w_LPM))
par->Flowflag = LPM;
else if (match(par->Tok[1], w_CMH))
par->Flowflag = CMH;
else if (match(par->Tok[1], w_CMD))
par->Flowflag = CMD;
else if (match(par->Tok[1], w_MLD))
par->Flowflag = MLD;
else if (match(par->Tok[1], w_SI))
par->Flowflag = LPS;
else
return (201);
} else if (match(par->Tok[0], w_PRESSURE)) {
if (n < 1)
return (0);
else if (match(par->Tok[1], w_PSI))
par->Pressflag = PSI;
else if (match(par->Tok[1], w_KPA))
par->Pressflag = KPA;
else if (match(par->Tok[1], w_METERS))
par->Pressflag = METERS;
else
return (201);
} else if (match(par->Tok[0], w_HEADLOSS)) {
if (n < 1)
return (0);
else if (match(par->Tok[1], w_HW))
hyd->Formflag = HW;
else if (match(par->Tok[1], w_DW))
hyd->Formflag = DW;
else if (match(par->Tok[1], w_CM))
hyd->Formflag = CM;
else
return (201);
} else if (match(par->Tok[0], w_HYDRAULIC)) {
if (n < 2)
return (0);
else if (match(par->Tok[1], w_USE))
out->Hydflag = USE;
else if (match(par->Tok[1], w_SAVE))
out->Hydflag = SAVE;
else
return (201);
if (n < 0) return (201);
if (match(par->Tok[0], w_UNITS))
{
if (n < 1) return (0);
else if (match(par->Tok[1], w_CFS)) par->Flowflag = CFS;
else if (match(par->Tok[1], w_GPM)) par->Flowflag = GPM;
else if (match(par->Tok[1], w_AFD)) par->Flowflag = AFD;
else if (match(par->Tok[1], w_MGD)) par->Flowflag = MGD;
else if (match(par->Tok[1], w_IMGD)) par->Flowflag = IMGD;
else if (match(par->Tok[1], w_LPS)) par->Flowflag = LPS;
else if (match(par->Tok[1], w_LPM)) par->Flowflag = LPM;
else if (match(par->Tok[1], w_CMH)) par->Flowflag = CMH;
else if (match(par->Tok[1], w_CMD)) par->Flowflag = CMD;
else if (match(par->Tok[1], w_MLD)) par->Flowflag = MLD;
else if (match(par->Tok[1], w_SI)) par->Flowflag = LPS;
else return (201);
}
else if (match(par->Tok[0], w_PRESSURE))
{
if (n < 1) return (0);
else if (match(par->Tok[1], w_EXPONENT)) return -1;
else if (match(par->Tok[1], w_PSI)) par->Pressflag = PSI;
else if (match(par->Tok[1], w_KPA)) par->Pressflag = KPA;
else if (match(par->Tok[1], w_METERS)) par->Pressflag = METERS;
else return (201);
}
else if (match(par->Tok[0], w_HEADLOSS))
{
if (n < 1) return (0);
else if (match(par->Tok[1], w_HW)) hyd->Formflag = HW;
else if (match(par->Tok[1], w_DW)) hyd->Formflag = DW;
else if (match(par->Tok[1], w_CM)) hyd->Formflag = CM;
else return (201);
}
else if (match(par->Tok[0], w_HYDRAULIC))
{
if (n < 2) return (0);
else if (match(par->Tok[1], w_USE)) out->Hydflag = USE;
else if (match(par->Tok[1], w_SAVE)) out->Hydflag = SAVE;
else return (201);
strncpy(out->HydFname, par->Tok[2], MAXFNAME);
} else if (match(par->Tok[0], w_QUALITY)) {
if (n < 1)
return (0);
else if (match(par->Tok[1], w_NONE))
qu->Qualflag = NONE;
else if (match(par->Tok[1], w_CHEM))
qu->Qualflag = CHEM;
else if (match(par->Tok[1], w_AGE))
qu->Qualflag = AGE;
else if (match(par->Tok[1], w_TRACE))
qu->Qualflag = TRACE;
else {
qu->Qualflag = CHEM;
strncpy(qu->ChemName, par->Tok[1], MAXID);
if (n >= 2)
strncpy(qu->ChemUnits, par->Tok[2], MAXID);
}
else if (match(par->Tok[0], w_QUALITY))
{
if (n < 1) return (0);
else if (match(par->Tok[1], w_NONE)) qu->Qualflag = NONE;
else if (match(par->Tok[1], w_CHEM)) qu->Qualflag = CHEM;
else if (match(par->Tok[1], w_AGE)) qu->Qualflag = AGE;
else if (match(par->Tok[1], w_TRACE)) qu->Qualflag = TRACE;
else
{
qu->Qualflag = CHEM;
strncpy(qu->ChemName, par->Tok[1], MAXID);
if (n >= 2) strncpy(qu->ChemUnits, par->Tok[2], MAXID);
}
if (qu->Qualflag == TRACE) /* Source tracing option */
{
/* Copy Trace Node ID to par->Tok[0] for error reporting */
strcpy(par->Tok[0], "");
if (n < 2)
return (212);
if (n < 2) return (212);
strcpy(par->Tok[0], par->Tok[2]);
qu->TraceNode = findnode(net,par->Tok[2]);
if (qu->TraceNode == 0)
return (212);
if (qu->TraceNode == 0) return (212);
strncpy(qu->ChemName, u_PERCENT, MAXID);
strncpy(qu->ChemUnits, par->Tok[2], MAXID);
}
if (qu->Qualflag == AGE) {
if (qu->Qualflag == AGE)
{
strncpy(qu->ChemName, w_AGE, MAXID);
strncpy(qu->ChemUnits, u_HOURS, MAXID);
}
} else if (match(par->Tok[0], w_MAP)) {
if (n < 1)
return (0);
}
else if (match(par->Tok[0], w_MAP))
{
if (n < 1) return (0);
strncpy(pr->MapFname, par->Tok[1], MAXFNAME); /* Map file name */
} else if (match(par->Tok[0], w_VERIFY)) {
}
else if (match(par->Tok[0], w_VERIFY))
{
/* Backward compatibility for verification file */
} else if (match(par->Tok[0], w_UNBALANCED)) /* Unbalanced option */
}
else if (match(par->Tok[0], w_UNBALANCED)) /* Unbalanced option */
{
if (n < 1)
return (0);
if (match(par->Tok[1], w_STOP))
hyd->ExtraIter = -1;
else if (match(par->Tok[1], w_CONTINUE)) {
if (n >= 2)
hyd->ExtraIter = atoi(par->Tok[2]);
else
hyd->ExtraIter = 0;
} else
return (201);
} else if (match(par->Tok[0], w_PATTERN)) /* Pattern option */
if (n < 1) return (0);
if (match(par->Tok[1], w_STOP)) hyd->ExtraIter = -1;
else if (match(par->Tok[1], w_CONTINUE))
{
if (n >= 2) hyd->ExtraIter = atoi(par->Tok[2]);
else hyd->ExtraIter = 0;
}
else return (201);
}
else if (match(par->Tok[0], w_PATTERN)) /* Pattern option */
{
if (n < 1)
return (0);
if (n < 1) return (0);
strncpy(par->DefPatID, par->Tok[1], MAXID);
} else
return (-1);
}
else if (match(par->Tok[0], w_DEMAND))
{
if (n < 2) return 0;
if (!match(par->Tok[1], w_MODEL)) return -1;
choice = findmatch(par->Tok[2], DemandModelTxt);
if (choice < 0) return 201;
hyd->DemandModel = choice;
}
else return (-1);
return (0);
} /* end of optionchoice */
@@ -1916,6 +1910,9 @@ int optionvalue(EN_Project *pr, int n)
** HEADLIMIT value
** FLOWLIMIT value
** MINIMUM PRESSURE value
** REQUIRED PRESSURE value
** PRESSURE EXPONENT value
** TOLERANCE value
** SEGMENTS value (not used)
@@ -1939,39 +1936,40 @@ int optionvalue(EN_Project *pr, int n)
double y;
/* Check for obsolete SEGMENTS keyword */
//if (match(par->Tok[0], w_SEGMENTS))
if (match(tok0, w_SEGMENTS))
return (0);
if (match(tok0, w_SEGMENTS)) return (0);
/* Check for missing value (which is permissible) */
if (match(tok0, w_SPECGRAV) || match(tok0, w_EMITTER) ||
match(tok0, w_DEMAND))
match(tok0, w_DEMAND) || match(tok0, w_MINIMUM) ||
match(tok0, w_REQUIRED) || match(tok0, w_PRESSURE) ||
match(tok0, w_PRECISION))
{
nvalue = 2;
if (n < nvalue)
return (0);
}
if (n < nvalue) return (0);
/* Check for valid numerical input */
if (!getfloat(par->Tok[nvalue], &y))
return (213);
if (!getfloat(par->Tok[nvalue], &y)) return (213);
/* Check for WQ tolerance option (which can be 0) */
if (match(tok0, w_TOLERANCE)) {
if (y < 0.0)
return (213);
if (match(tok0, w_TOLERANCE))
{
if (y < 0.0) return (213);
qu->Ctol = y; /* Quality tolerance*/
return (0);
}
/* Check for Diffusivity option */
if (match(tok0, w_DIFFUSIVITY)) {
if (y < 0.0)
return (213);
if (match(tok0, w_DIFFUSIVITY))
{
if (y < 0.0) return (213);
qu->Diffus = y;
return (0);
}
/* Check for Damping Limit option */
if (match(tok0, w_DAMPLIMIT)) {
if (match(tok0, w_DAMPLIMIT))
{
hyd->DampLimit = y;
return (0);
}
@@ -1992,41 +1990,51 @@ int optionvalue(EN_Project *pr, int n)
return 0;
}
/* Check for pressure dependent demand params */
else if (match(tok0, w_MINIMUM))
{
if (y < 0.0) return 213;
hyd->Pmin = y;
return 0;
}
else if (match(tok0, w_REQUIRED))
{
if (y < 0.0) return 213;
hyd->Preq = y;
return 0;
}
else if (match(tok0, w_PRESSURE))
{
if (y < 0.0) return 213;
hyd->Pexp = y;
return 0;
}
/* All other options must be > 0 */
if (y <= 0.0)
return (213);
if (y <= 0.0) return (213);
/* Assign value to specified option */
if (match(tok0, w_VISCOSITY))
hyd->Viscos = y; /* Viscosity */
else if (match(tok0, w_SPECGRAV))
hyd->SpGrav = y; /* Spec. gravity */
else if (match(tok0, w_TRIALS))
hyd->MaxIter = (int)y; /* Max. trials */
if (match(tok0, w_VISCOSITY)) hyd->Viscos = y; /* Viscosity */
else if (match(tok0, w_SPECGRAV)) hyd->SpGrav = y; /* Spec. gravity */
else if (match(tok0, w_TRIALS)) hyd->MaxIter = (int)y; /* Max. trials */
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(tok0, w_HTOL))
hyd->Htol = y;
else if (match(tok0, w_QTOL))
hyd->Qtol = y;
else if (match(tok0, w_RQTOL)) {
if (y >= 1.0)
return (213);
else if (match(tok0, w_HTOL)) hyd->Htol = y;
else if (match(tok0, w_QTOL)) hyd->Qtol = y;
else if (match(tok0, w_RQTOL))
{
if (y >= 1.0) return (213);
hyd->RQtol = y;
} else if (match(tok0, w_CHECKFREQ))
hyd->CheckFreq = (int)y;
else if (match(tok0, w_MAXCHECK))
hyd->MaxCheck = (int)y;
else if (match(tok0, w_EMITTER))
hyd->Qexp = 1.0 / y;
else if (match(tok0, w_DEMAND))
hyd->Dmult = y;
else
return (201);
}
else if (match(tok0, w_CHECKFREQ)) hyd->CheckFreq = (int)y;
else if (match(tok0, w_MAXCHECK)) hyd->MaxCheck = (int)y;
else if (match(tok0, w_EMITTER)) hyd->Qexp = 1.0 / y;
else if (match(tok0, w_DEMAND)) hyd->Dmult = y;
else return (201);
return (0);
} /* end of optionvalue */
@@ -2091,7 +2099,8 @@ int getpumpcurve(EN_Project *pr, int n)
return (0);
}
int powercurve(double h0, double h1, double h2, double q1, double q2, double *a, double *b, double *c)
int powercurve(double h0, double h1, double h2, double q1, double q2,
double *a, double *b, double *c)
/*
**---------------------------------------------------------
** Input: h0 = shutoff head

4
src/output.c
View File

@@ -19,14 +19,12 @@ AUTHOR: L. Rossman
#else
#include <stdlib.h>
#endif
#include "epanet2.h"
#include "funcs.h"
#include "text.h"
#include "types.h"
#include <math.h>
#define EXTERN extern
#include "hash.h"
#include "vars.h"
/* write x[1] to x[n] to file */
size_t f_save(REAL4 *x, int n, FILE *file) {

4
src/quality.c
View File

@@ -54,15 +54,13 @@ AUTHOR: L. Rossman
#else
#include <stdlib.h>
#endif
#include "hash.h"
#include "text.h"
#include "types.h"
#include "epanet2.h"
#include "funcs.h"
#include <math.h>
#define EXTERN extern
#include "mempool.h"
#include "vars.h"
/*
** Macros to identify upstream & downstream nodes of a link

9
src/report.c
View File

@@ -32,15 +32,13 @@ formatted string S to the report file.
#else
#include <stdlib.h>
#endif
#include "epanet2.h"
#include "funcs.h"
#include "hash.h"
#include "text.h"
#include "types.h"
#include <math.h>
#include <time.h>
#define EXTERN extern
#include "vars.h"
#undef WINDOWS
#ifdef _WIN32
@@ -50,12 +48,13 @@ formatted string S to the report file.
#define MAXCOUNT 10 /* Max. # of disconnected nodes listed */
/* Defined in enumstxt.h in EPANET.C */
/* Defined in enumstxt.h */
extern char *NodeTxt[];
extern char *LinkTxt[];
extern char *StatTxt[];
extern char *TstatTxt[];
extern char *RptFormTxt[];
extern char *DemandModelTxt[];
typedef REAL4 *Pfloat;
void writenodetable(EN_Project *pr, Pfloat *);
@@ -224,6 +223,8 @@ void writesummary(EN_Project *pr)
writeline(pr, s);
sprintf(s, FMT25, RptFormTxt[hyd->Formflag]);
writeline(pr, s);
sprintf(s, FMT25a, DemandModelTxt[hyd->DemandModel]);
writeline(pr, s);
sprintf(s, FMT26, time->Hstep * pr->Ucf[TIME], rep->Field[TIME].Units);
writeline(pr, s);
sprintf(s, FMT27, hyd->Hacc);

4
src/rules.c
View File

@@ -29,13 +29,11 @@ AUTHOR: L. Rossman
#else
#include <stdlib.h>
#endif
#include "epanet2.h"
#include "funcs.h"
#include "hash.h"
#include "text.h"
#include "types.h"
#define EXTERN extern
#include "vars.h"
enum Rulewords {
r_RULE,

13
src/text.h
View File

@@ -145,6 +145,12 @@ AUTHOR: L. Rossman
#define w_FLOWCHANGE "FLOWCHANGE"
#define w_HEADERROR "HEADERROR"
#define w_MODEL "MODEL"
#define w_DDA "DDA"
#define w_PDA "PDA"
#define w_REQUIRED "REQ"
#define w_EXPONENT "EXP"
#define w_SECONDS "SEC"
#define w_MINUTES "MIN"
#define w_HOURS "HOU"
@@ -323,6 +329,9 @@ AUTHOR: L. Rossman
#define t_perM3 " /m3"
#define t_perMGAL "/Mgal"
#define t_DIFFER "DIFFERENTIAL"
#define t_FIXED "Fixed Demands"
#define t_POWER "Power Function"
#define t_ORIFICE "Orifice Flow"
/* ------------------ Format Messages ------------------*/
@@ -362,6 +371,7 @@ AUTHOR: L. Rossman
#define FMT23 " Number of Pumps ................... %-d"
#define FMT24 " Number of Valves .................. %-d"
#define FMT25 " Headloss Formula .................. %s"
#define FMT25a " Nodal Demand Model ................ %s"
#define FMT26 " Hydraulic Timestep ................ %-.2f %s"
#define FMT27 " Hydraulic Accuracy ................ %-.6f"
@@ -420,7 +430,8 @@ AUTHOR: L. Rossman
/*** End of update ***/
#define FMT66 " maximum flow change = %.4f for Link %s"
#define FMT67 " maximum head error = %.4f for Link %s\n"
#define FMT67 " maximum flow change = %.4f for Node %s"
#define FMT68 " maximum head error = %.4f for Link %s\n"
/* -------------------- Energy Report Table ------------------- */
#define FMT71 "Energy Usage:"

24
src/types.h
View File

@@ -301,6 +301,11 @@ typedef enum {
MAX_ENERGY_STATS
} EnergyStats;
typedef enum {
DDA, // Demand Driven Analysis
PDA // Pressure Driven Analysis
} DemandModelType;
/*
------------------------------------------------------
Global Data Structures
@@ -762,7 +767,8 @@ typedef struct {
typedef struct {
double
*NodeDemand, /* Node actual demand */
*NodeDemand, // Node actual total outflow
*DemandFlows, // Demand outflows
*EmitterFlows, /* Emitter flows */
*LinkSetting, /* Link settings */
*LinkFlows, /* Link flows */
@@ -771,10 +777,13 @@ typedef struct {
Qtol, /* Flow rate tolerance */
RQtol, /* Flow resistance tolerance */
Hexp, /* Exponent in headloss formula */
Qexp, /* Exponent in orifice formula */
Qexp, /* Exponent in emitter formula */
Pexp, // Exponent in demand formula
Pmin, // Pressure needed for any demand
Preq, // Pressure needed for full demand
Dmult, /* Demand multiplier */
Hacc, /* Hydraulics solution accuracy */
Hacc, /* Hydraulics solution accuracy */
FlowChangeLimit, /* Hydraulics flow change limit */
HeadErrorLimit, /* Hydraulics head error limit */
@@ -790,7 +799,8 @@ typedef struct {
int
DefPat, /* Default demand pattern */
Epat; /* Energy cost time pattern */
Epat, /* Energy cost time pattern */
DemandModel; // Fixed or pressure dependent
StatType
*LinkStatus, /* Link status */
@@ -808,8 +818,10 @@ typedef struct {
Formflag; /* Hydraulic formula flag */
/* Info about hydraulic solution */
double relativeError;
int iterations;
double RelativeError;
double MaxHeadError;
double MaxFlowChange;
int Iterations;
/* Flag used to halt taking further time steps */
int Haltflag;