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

Merge pull request #350 from LRossman/lrossman-dev

Browse Source 
Code cleanup
This commit is contained in:
Michael Tryby
2018-11-28 17:26:37 -05:00
committed by GitHub
parent 2988800448 20bc6358ff
commit 9ff0056b7c
36 changed files with 15193 additions and 15419 deletions
Download Patch File Download Diff File Expand all files Collapse all files
include/epanet2.h
+1100 -992
View File
File diff suppressed because it is too large Load Diff
src/enumstxt.h
+10 -10
View File
@@ -1,14 +1,14 @@
/*
***********************************************************************
ENUMSTXT.H -- Text strings for enumerated data types in EPANET
VERSION: 2.00
DATE: 5/8/00
AUTHOR: L. Rossman
US EPA - NRMRL
**********************************************************************
*****************************************************************************
Project: OWA EPANET
Version: 2.2
Module: enumstxt.h
Description: text strings for enumerated data types
Authors: see AUTHORS
Copyright: see AUTHORS
License: see LICENSE
Last Updated: 11/27/2018
******************************************************************************
*/
#ifndef ENUMSTXT_H
src/epanet.c
+4637 -5273
View File
File diff suppressed because it is too large Load Diff
src/errors.dat
+1 -4
View File
@@ -57,7 +57,7 @@ DAT(257,ENERR_NO_RULE,"function applied to nonexistent rule")
DAT(258,ENERR_NO_CONDITION_ACTION,"function applied to nonexistent rule clause")
DAT(260,ENERR_DEL_TRACE_NODE,"cannot delete node assigned as a Trace Node")
DAT(261,ENERR_DEL_NODE_LINK, "cannot delete a node or link contained in a control or rule")
DAT(261,ENERR_DEL_NODE_LINK, "cannot delete a node or link contained in a control")
DAT(301,ENERR_FILES_ARE_SAME,"identical file names")
DAT(302,ENERR_CANT_OPEN_INP,"cannot open input file")
@@ -69,9 +69,6 @@ DAT(307,ENERR_CANT_READ_HYD,"cannot read hydraulics file")
DAT(308,ENERR_CANT_SAVE_RES,"cannot save results to file")
DAT(309,ENERR_CANT_SAVE_RPT,"cannot save results to report file")
DAT(401,ENERR_QSTEP_NOT_DIVISIBLE,"Qstep is not dividable by Hstep")
DAT(411,ENERR_NO_MEM_ALLOC,"no memory allocated for results")
DAT(412,ENERR_NO_RESULTS_NO_FILE,"no results; binary file hasn't been opened")
DAT(435,ENERR_RUN_TERMINATED,"run terminated; no results in binary file")
src/funcs.h
+155 -207
View File
@@ -1,227 +1,175 @@
/*
**************************************************************************
FUNCS.H -- Function Prototypes for EPANET Program
VERSION: 2.00
DATE: 5/8/00
9/25/00
10/25/00
12/29/00
3/1/01
2/14/08 (2.00.12)
AUTHOR: L. Rossman
US EPA - NRMRL
**************************************************************************
******************************************************************************
Project: OWA EPANET
Version: 2.2
Module: funcs.h
Description: prototypes of external functions called by various modules
Authors: see AUTHORS
Copyright: see AUTHORS
License: see LICENSE
Last Updated: 11/27/2018
******************************************************************************
*/
/*****************************************************************/
/* Most float arguments have been changed to double - 7/3/07 */
/*****************************************************************/
/* ------- EPANET.C --------------------*/
/*
** NOTE: The exportable functions that can be called
** via the DLL are prototyped in TOOLKIT.H.
*/
#ifndef FUNCS_H
#define FUNCS_H
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, const char *,
const char *,const char *); /* Opens input & report files */
int openhydfile(EN_Project *pr); /* Opens hydraulics file */
int openoutfile(EN_Project *pr); /* Opens binary output file */
int strcomp(const char *, const char *); /* Compares two strings */
char* getTmpName(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 writewin(void (*vp)(char *), char *); /* Passes text to calling app */
// ------- PROJECT.C ------------
/* ------- INPUT1.C --------------------*/
int getdata(EN_Project *pr); /* Gets network data */
void setdefaults(EN_Project *pr); /* Sets default values */
void initreport(report_options_t *r); /* Initializes report options */
void adjustdata(EN_Project *pr); /* Adjusts input data */
int inittanks(EN_Project *pr); /* Initializes tank levels */
void initunits(EN_Project *pr); /* Determines reporting units */
void convertunits(EN_Project *pr); /* Converts data to std. units*/
void initpointers(Project *);
int allocdata(Project *);
void freeTmplist(STmplist *);
void freeFloatlist(SFloatlist *);
void freedata(Project *);
/* -------- INPUT2.C -------------------*/
int netsize(EN_Project *pr); /* Determines network size */
int readdata(EN_Project *pr); /* Reads in network data */
int newline(EN_Project *pr, int, char *); /* Processes new line of data */
int addnodeID(EN_Network *n, int, char *); /* Adds node ID to data base */
int addlinkID(EN_Network *n, int, char *); /* Adds link ID to data base */
int addpattern(parser_data_t *par, char *); /* Adds pattern to data base */
int addcurve(parser_data_t *par, char *); /* Adds curve to data base */
STmplist *findID(char *, STmplist *); /* Locates ID on linked list */
int unlinked(EN_Project *pr); /* Checks for unlinked nodes */
int getpumpparams(EN_Project *pr); /* Computes pump curve coeffs.*/
int updatepumpparams(EN_Project *pr, int); // Updates pump curve coeffs.
int getpatterns(EN_Project *pr); /* Gets pattern data from list*/
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 */
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 */
int openfiles(Project *, const char *, const char *,const char *);
int openhydfile(Project *);
int openoutfile(Project *);
/* ---------- INPUT3.C -----------------*/
int juncdata(EN_Project *pr); /* Processes junction data */
int tankdata(EN_Project *pr); /* Processes tank data */
int pipedata(EN_Project *pr); /* Processes pipe data */
int pumpdata(EN_Project *pr); /* Processes pump data */
int valvedata(EN_Project *pr); /* Processes valve data */
int patterndata(EN_Project *pr); /* Processes pattern data */
int curvedata(EN_Project *pr); /* Processes curve data */
int coordata(EN_Project *pr); /* Processes coordinate data */
int demanddata(EN_Project *pr); /* Processes demand data */
int controldata(EN_Project *pr); /* Processes simple controls */
int energydata(EN_Project *pr); /* Processes energy data */
int sourcedata(EN_Project *pr); /* Processes source data */
int emitterdata(EN_Project *pr); /* Processes emitter data */
int qualdata(EN_Project *pr); /* Processes quality data */
int reactdata(EN_Project *pr); /* Processes reaction data */
int mixingdata(EN_Project *pr); /* Processes tank mixing data */
int statusdata(EN_Project *pr); /* Processes link status data */
int reportdata(EN_Project *pr); /* Processes report options */
int timedata(EN_Project *pr); /* Processes time options */
int optiondata(EN_Project *pr); /* Processes analysis options */
int optionchoice(EN_Project *pr, int); /* Processes option choices */
int optionvalue(EN_Project *pr, int); /* Processes option values */
int getpumpcurve(EN_Project *pr, int); /* Constructs a pump curve */
int powercurve(double, double, double, /* Coeffs. of power pump curve*/
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 */
int buildadjlists(Network *);
void freeadjlists(Network *);
/* -------------- RULES.C --------------*/
void initrules(EN_Project *pr); /* Initializes rule base */
void addrule(parser_data_t *par, char *); /* Adds rule to rule base */
int allocrules(EN_Project *pr); /* Allocates memory for rule */
void adjustrules(EN_Project *pr, int, int); // Shifts object indices down
void adjusttankrules(EN_Project *pr); // Shifts tank indices up
int ruledata(EN_Project *pr); /* Processes rule input data */
int checkrules(EN_Project *pr, long); /* Checks all rules */
void deleterule(EN_Project *pr, int); // Deletes a rule
void freerules(EN_Project *pr); /* Frees rule base memory */
int writerule(EN_Project *pr, FILE *, int); /* Writes rule to an INP file */
void ruleerrmsg(EN_Project *pr); /* Reports rule parser error */
Spremise *getpremise(Spremise *, int); // Retrieves a rule's premise
Saction *getaction(Saction *, int); // Retrieves a rule's action
int incontrols(Project *, int, int);
int findnode(Network *, char *);
int findlink(Network *, char *);
int findtank(Network *, int);
int findvalve(Network *, int);
int findpump(Network *, int);
/* ------------- REPORT.C --------------*/
int writereport(EN_Project *pr); /* Writes formatted report */
void writelogo(EN_Project *pr); /* Writes program logo */
void writesummary(EN_Project *pr); /* Writes network summary */
void writehydstat(EN_Project *pr, int,double); /* Writes hydraulic status */
void writeenergy(EN_Project *pr); /* Writes energy usage */
int writeresults(EN_Project *pr); /* Writes node/link results */
void writeheader(EN_Project *pr, int,int); /* Writes heading on report */
void writeline(EN_Project *pr, char *); /* Writes line to report file */
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 */
void writehyderr(EN_Project *pr, int); /* Writes hydraulic error msg.*/
void writemassbalance(EN_Project *pr); // Writes mass balance ratio
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 */
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*/
int getnodetype(EN_Network *net, int); /* Determines node type */
char *getTmpName(char *);
int strcomp(const char *, const char *);
double interp(int, double [], double [], double);
char *geterrmsg(int, char *);
void errmsg(Project *, int);
void writewin(void (*vp)(char *), char *);
/* --------- HYDRAUL.C -----------------*/
int openhyd(EN_Project *pr); /* Opens hydraulics solver */
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 */
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 */
int tanktimestep(EN_Project *pr, long *); /* Time till tanks fill/drain */
void getenergy(EN_Project *pr, int, double *,
double *); /* Computes link energy use */
double tankvolume(EN_Project *pr, int, double); /* Finds tank vol. from grade */
double tankgrade(EN_Project *pr, int, double); /* Finds tank grade from vol. */
// ------- INPUT1.C ----------------
/* ----------- HYDSOLVER.C - ----------*/
int hydsolve(EN_Project *pr, int *,double *); /* Solves network equations */
int getdata(Project *);
void setdefaults(Project *);
void initreport(Report *);
void adjustdata(Project *);
int inittanks(Project *);
void initunits(Project *);
void convertunits(Project *);
/* ----------- HYDCOEFFS.C --------------*/
void resistcoeff(EN_Project *pr, int k); /* Finds pipe flow resistance */
void headlosscoeffs(EN_Project *pr); // Finds link head loss coeffs.
void matrixcoeffs(EN_Project *pr); /* Finds hyd. matrix coeffs. */
void emitheadloss(EN_Project *pr, int, // Finds emitter head loss
double *, double *);
double demandflowchange(EN_Project *pr, int, // Change in demand outflow
double, double);
void demandparams(EN_Project *pr, double *, // PDA function parameters
double *);
//-------- INPUT2.C -----------------
/* ----------- SMATRIX.C ---------------*/
int createsparse(EN_Project *pr); /* Creates sparse matrix */
void freesparse(EN_Project *pr); /* Frees matrix memory */
int linsolve(EN_Project *pr, int); /* Solves set of linear eqns. */
int netsize(Project *);
int readdata(Project *);
int updatepumpparams(Project *, int);
int getpatterns(Project *);
int getcurves(Project *);
int findmatch(char *, char *[]);
int match(const char *, const char *);
int gettokens(char *, char **, int, char *);
int getfloat(char *, double *);
double hour(char *, char *);
int setreport(Project *, char *);
/* ----------- QUALITY.C ---------------*/
int openqual(EN_Project *pr); /* Opens WQ solver system */
int initqual(EN_Project *pr); /* Initializes WQ solver */
int runqual(EN_Project *pr, long *); /* Gets current WQ results */
int nextqual(EN_Project *pr, long *); /* Updates WQ by hyd.timestep */
int stepqual(EN_Project *pr, long *); /* Updates WQ by WQ time step */
int closequal(EN_Project *pr); /* Closes WQ solver system */
double avgqual(EN_Project *pr, int); /* Finds avg. quality in pipe */
// ------- INPUT3.C -----------------
/* ------------ OUTPUT.C ---------------*/
int savenetdata(EN_Project *pr); /* Saves basic data to file */
int savehyd(EN_Project *pr, long *); /* Saves hydraulic solution */
int savehydstep(EN_Project *pr, long *); /* Saves hydraulic timestep */
int saveenergy(EN_Project *pr); /* Saves energy usage */
int readhyd(EN_Project *pr, long *); /* Reads hydraulics from file */
int readhydstep(EN_Project *pr, long *); /* Reads time step from file */
int saveoutput(EN_Project *pr); /* Saves results to file */
int nodeoutput(EN_Project *pr, int, REAL4 *,
double); /* Saves node results to file */
int linkoutput(EN_Project *pr, int, REAL4 *,
double); /* Saves link results to file */
int savefinaloutput(EN_Project *pr); /* Finishes saving output */
int savetimestat(EN_Project *pr, REAL4 *,
HdrType); /* Saves time stats to file */
int savenetreacts(EN_Project *pr, double,
double,double, double); /* Saves react. rates to file */
int saveepilog(EN_Project *pr); /* Saves output file epilog */
int juncdata(Project *);
int tankdata(Project *);
int pipedata(Project *);
int pumpdata(Project *);
int valvedata(Project *);
int patterndata(Project *);
int curvedata(Project *);
int coordata(Project *);
int demanddata(Project *);
int controldata(Project *);
int energydata(Project *);
int sourcedata(Project *);
int emitterdata(Project *);
int qualdata(Project *);
int reactdata(Project *);
int mixingdata(Project *);
int statusdata(Project *);
int reportdata(Project *);
int timedata(Project *);
int optiondata(Project *);
int valvecheck(Project *, int, int, int);
// ------- RULES.C ------------------
/* ------------ INPFILE.C --------------*/
int saveinpfile(EN_Project *pr, const char *); /* Saves network to text file */
void initrules(Project *);
void addrule(Parser *, char *);
void deleterule(Project *, int);
int allocrules(Project *);
void freerules(Project *);
int ruledata(Project *);
void ruleerrmsg(Project *);
void adjustrules(Project *, int, int);
void adjusttankrules(Project *);
Spremise *getpremise(Spremise *, int);
Saction *getaction(Saction *, int);
int writerule(Project *, FILE *, int);
int checkrules(Project *, long);
// ------- REPORT.C -----------------
int writereport(Project *);
void writelogo(Project *);
void writesummary(Project *);
void writehydstat(Project *, int, double);
void writeheader(Project *, int,int);
void writeline(Project *, char *);
void writerelerr(Project *, int, double);
void writestatchange(Project *, int,char,char);
void writecontrolaction(Project *, int, int);
void writeruleaction(Project *, int, char *);
int writehydwarn(Project *, int,double);
void writehyderr(Project *, int);
void writemassbalance(Project *);
void writetime(Project *, char *);
char *clocktime(char *, long);
// ------- HYDRAUL.C -----------------
int openhyd(Project *);
void inithyd(Project *, int initFlags);
int runhyd(Project *, long *);
int nexthyd(Project *, long *);
void closehyd(Project *);
void setlinkstatus(Project *, int, char, StatusType *, double *);
void setlinksetting(Project *, int, double, StatusType *, double *);
int tanktimestep(Project *, long *);
void getenergy(Project *, int, double *, double *);
double tankvolume(Project *, int, double);
double tankgrade(Project *, int, double);
// ------- HYDCOEFFS.C -----------------
void resistcoeff(Project *, int);
void headlosscoeffs(Project *);
void matrixcoeffs(Project *);
void emitheadloss(Project *, int, double *, double *);
double demandflowchange(Project *, int, double, double);
void demandparams(Project *, double *, double *);
// ------- QUALITY.C --------------------
int openqual(Project *);
int initqual(Project *);
int runqual(Project *, long *);
int nextqual(Project *, long *);
int stepqual(Project *, long *);
int closequal(Project *);
double avgqual(Project *, int);
// ------- OUTPUT.C ---------------------
int savenetdata(Project *);
int savehyd(Project *, long *);
int savehydstep(Project *, long *);
int saveenergy(Project *);
int readhyd(Project *, long *);
int readhydstep(Project *, long *);
int saveoutput(Project *);
int savefinaloutput(Project *);
// ------- INPFILE.C --------------------
int saveinpfile(Project *, const char *);
#endif
src/genmmd.c
-1
View File
@@ -4,7 +4,6 @@
*
*/
#include <math.h>
int genmmd(int* neqns, int* xadj, int* adjncy, int* invp, int* perm,
src/hash.c
+21 -17
View File
@@ -1,20 +1,14 @@
/*-----------------------------------------------------------------------------
** hash.c
**
** Implementation of a simple Hash Table that uses a string as a key
** and an associated integer as data.
**
** Written by L. Rossman
** Last Updated on 10/19/18
**
** Interface Functions:
** hashtable_create - creates a hash table
** hashtable_insert - inserts a string & its data value into a table
** hashtable_find - retrieves the data value associated with a string
** hashtable_findkey - retrieves the key associated with a data value
** hashtable_delete - deletes an entry from a table
** hashtable_free - frees a hash table
**
/*
******************************************************************************
Project: OWA EPANET
Version: 2.2
Module: hash.c
Description: implementation of a simple hash table
Authors: see AUTHORS
Copyright: see AUTHORS
License: see LICENSE
Last Updated: 11/27/2018
******************************************************************************
*/
#ifndef __APPLE__
@@ -27,6 +21,7 @@
#define HASHTABLEMAXSIZE 128000
// An entry in the hash table
typedef struct DataEntryStruct
{
char *key;
@@ -34,6 +29,7 @@ typedef struct DataEntryStruct
struct DataEntryStruct *next;
} DataEntry;
// Hash a string to an integer
unsigned int gethash(char *str)
{
unsigned int hash = 5381;
@@ -47,6 +43,7 @@ unsigned int gethash(char *str)
return retHash;
}
// Produce a duplicate string
char *dupstr(const char *s)
{
size_t size = strlen(s) + 1;
@@ -55,6 +52,7 @@ char *dupstr(const char *s)
return p;
}
// Create a hash table
HashTable *hashtable_create()
{
int i;
@@ -66,6 +64,7 @@ HashTable *hashtable_create()
return ht;
}
// Insert an entry into the hash table
int hashtable_insert(HashTable *ht, char *key, int data)
{
unsigned int i = gethash(key);
@@ -80,6 +79,7 @@ int hashtable_insert(HashTable *ht, char *key, int data)
return 1;
}
// Change the hash table's data entry for a particular key
int hashtable_update(HashTable *ht, char *key, int new_data)
{
unsigned int i = gethash(key);
@@ -99,6 +99,7 @@ int hashtable_update(HashTable *ht, char *key, int new_data)
return NOTFOUND;
}
// Delete an entry in the hash table
int hashtable_delete(HashTable *ht, char *key)
{
unsigned int i = gethash(key);
@@ -124,6 +125,7 @@ int hashtable_delete(HashTable *ht, char *key)
return NOTFOUND;
}
// Find the data for a particular key
int hashtable_find(HashTable *ht, char *key)
{
unsigned int i = gethash(key);
@@ -142,6 +144,7 @@ int hashtable_find(HashTable *ht, char *key)
return NOTFOUND;
}
// Find a particular key in the hash table
char *hashtable_findkey(HashTable *ht, char *key)
{
unsigned int i = gethash(key);
@@ -156,6 +159,7 @@ char *hashtable_findkey(HashTable *ht, char *key)
return NULL;
}
// Delete a hash table and free all of its memory
void hashtable_free(HashTable *ht)
{
DataEntry *entry, *nextentry;
src/hash.h
+11 -5
View File
@@ -1,9 +1,15 @@
/* HASH.H
**
** Header file for Hash Table module HASH.C
**
/*
******************************************************************************
Project: OWA EPANET
Version: 2.2
Module: hash.h
Description: header for a simple hash table
Authors: see AUTHORS
Copyright: see AUTHORS
License: see LICENSE
Last Updated: 11/27/2018
******************************************************************************
*/
#ifndef HASH_H
#define HASH_H
src/hydcoeffs.c
+226 -218
View File
@@ -1,9 +1,14 @@
/*
*********************************************************************
HYDCOEFFS.C -- hydraulic coefficients for the EPANET Program
*******************************************************************
******************************************************************************
Project: OWA EPANET
Version: 2.2
Module: hydcoeffs.c
Description: computes coefficients for a hydraulic solution matrix
Authors: see AUTHORS
Copyright: see AUTHORS
License: see LICENSE
Last Updated: 11/27/2018
******************************************************************************
*/
#include <stdio.h>
@@ -29,40 +34,44 @@ const double AA = -1.5634601348517065795e+00; // -2*.9*2/ln(10)
const double AB = 3.28895476345399058690e-03; // 5.74/(4000^.9)
const double AC = -5.14214965799093883760e-03; // AA*AB
// External functions
//void resistcoeff(EN_Project *pr, int k);
//void headlosscoeffs(EN_Project *pr);
//void matrixcoeffs(EN_Project *pr);
//void emitheadloss(EN_Project *pr, int i, double *hloss, double *dhdq);
//double demandflowchange(EN_Project *pr, int i, double dp, double n);
//void demandparams(EN_Project *pr, double *dp, double *n);
// Definitions of very small and very big coefficients
const double CSMALL = 1.e-6;
const double CBIG = 1.e8;
// Exported functions
//void resistcoeff(Project *, int );
//void headlosscoeffs(Project *);
//void matrixcoeffs(Project *);
//void emitheadloss(Project *, int, double *, double *);
//double demandflowchange(Project *, int, double, double);
//void demandparams(Project *, double *, double *);
// 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 linkcoeffs(Project *pr);
static void nodecoeffs(Project *pr);
static void valvecoeffs(Project *pr);
static void emittercoeffs(Project *pr);
static void demandcoeffs(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);
static void pipecoeff(Project *pr, int k);
static void DWpipecoeff(Project *pr, int k);
static double frictionFactor(double q, double e, double s, 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 pumpcoeff(Project *pr, int k);
static void curvecoeff(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);
static void valvecoeff(Project *pr, int k);
static void gpvcoeff(Project *pr, int k);
static void pbvcoeff(Project *pr, int k);
static void tcvcoeff(Project *pr, int k);
static void prvcoeff(Project *pr, int k, int n1, int n2);
static void psvcoeff(Project *pr, int k, int n1, int n2);
static void fcvcoeff(Project *pr, int k, int n1, int n2);
void resistcoeff(EN_Project *pr, int k)
void resistcoeff(Project *pr, int k)
/*
**--------------------------------------------------------------------
** Input: k = link index
@@ -71,10 +80,10 @@ void resistcoeff(EN_Project *pr, int k)
**--------------------------------------------------------------------
*/
{
double e, d, L;
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
double e, d, L;
Slink *link = &net->Link[k];
link->Qa = 0.0;
@@ -121,7 +130,7 @@ void resistcoeff(EN_Project *pr, int k)
}
void headlosscoeffs(EN_Project *pr)
void headlosscoeffs(Project *pr)
/*
**--------------------------------------------------------------
** Input: none
@@ -131,9 +140,10 @@ void headlosscoeffs(EN_Project *pr)
**--------------------------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
int k;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
for (k = 1; k <= net->Nlinks; k++)
{
@@ -159,13 +169,13 @@ void headlosscoeffs(EN_Project *pr)
case PRV:
case PSV:
if (hyd->LinkSetting[k] == MISSING) valvecoeff(pr, k);
else hyd->solver.P[k] = 0.0;
else hyd->P[k] = 0.0;
}
}
}
void matrixcoeffs(EN_Project *pr)
void matrixcoeffs(Project *pr)
/*
**--------------------------------------------------------------
** Input: none
@@ -174,16 +184,16 @@ void matrixcoeffs(EN_Project *pr)
**--------------------------------------------------------------
*/
{
hydraulics_t *hyd = &pr->hydraulics;
solver_t *sol = &hyd->solver;
EN_Network *net = &pr->network;
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Smatrix *sm = &hyd->smatrix;
// 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));
// coeffs. (Aij), r.h.s. coeffs. (F) and node excess flow (Xflow)
memset(sm->Aii, 0, (net->Nnodes + 1) * sizeof(double));
memset(sm->Aij, 0, (sm->Ncoeffs + 1) * sizeof(double));
memset(sm->F, 0, (net->Nnodes + 1) * sizeof(double));
memset(hyd->Xflow, 0, (net->Nnodes + 1) * sizeof(double));
// Compute matrix coeffs. from links, emitters, and nodal demands
linkcoeffs(pr);
@@ -199,7 +209,7 @@ void matrixcoeffs(EN_Project *pr)
}
void linkcoeffs(EN_Project *pr)
void linkcoeffs(Project *pr)
/*
**--------------------------------------------------------------
** Input: none
@@ -209,80 +219,81 @@ void linkcoeffs(EN_Project *pr)
**--------------------------------------------------------------
*/
{
int k, n1, n2;
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Smatrix *sm = &hyd->smatrix;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
solver_t *sol = &hyd->solver;
int k, n1, n2;
Slink *link;
// Examine each link of network
for (k = 1; k <= net->Nlinks; k++)
{
if (sol->P[k] == 0.0) continue;
if (hyd->P[k] == 0.0) continue;
link = &net->Link[k];
n1 = link->N1; // Start node of link
n2 = link->N2; // End node of link
// Update nodal flow balance (X_tmp)
// Update nodal flow excess (Xflow)
// (Flow out of node is (-), flow into node is (+))
hyd->X_tmp[n1] -= hyd->LinkFlows[k];
hyd->X_tmp[n2] += hyd->LinkFlows[k];
hyd->Xflow[n1] -= hyd->LinkFlow[k];
hyd->Xflow[n2] += hyd->LinkFlow[k];
// Add to off-diagonal coeff. of linear system matrix
sol->Aij[sol->Ndx[k]] -= sol->P[k];
sm->Aij[sm->Ndx[k]] -= hyd->P[k];
// Update linear system coeffs. associated with start node n1
// ... 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.
sm->Aii[sm->Row[n1]] += hyd->P[k]; // Diagonal coeff.
sm->F[sm->Row[n1]] += hyd->Y[k]; // RHS coeff.
}
// ... node n1 is a tank/reservoir
else sol->F[sol->Row[n2]] += (sol->P[k] * hyd->NodeHead[n1]);
else sm->F[sm->Row[n2]] += (hyd->P[k] * hyd->NodeHead[n1]);
// Update linear system coeffs. associated with end node n2
// ... 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.
sm->Aii[sm->Row[n2]] += hyd->P[k]; // Diagonal coeff.
sm->F[sm->Row[n2]] -= hyd->Y[k]; // RHS coeff.
}
// ... node n2 is a tank/reservoir
else sol->F[sol->Row[n1]] += (sol->P[k] * hyd->NodeHead[n2]);
else sm->F[sm->Row[n1]] += (hyd->P[k] * hyd->NodeHead[n2]);
}
}
void nodecoeffs(EN_Project *pr)
void nodecoeffs(Project *pr)
/*
**----------------------------------------------------------------
** Input: none
** Output: none
** Purpose: completes calculation of nodal flow balance array
** (X_tmp) & r.h.s. (F) of linearized hydraulic eqns.
** (Xflow) & r.h.s. (F) of linearized hydraulic eqns.
**----------------------------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Smatrix *sm = &hyd->smatrix;
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 balance & add flow balance to RHS array F
// flow excess & add flow excess to RHS array F
for (i = 1; i <= net->Njuncs; i++)
{
hyd->X_tmp[i] -= hyd->DemandFlows[i];
sol->F[sol->Row[i]] += hyd->X_tmp[i];
hyd->Xflow[i] -= hyd->DemandFlow[i];
sm->F[sm->Row[i]] += hyd->Xflow[i];
}
}
void valvecoeffs(EN_Project *pr)
void valvecoeffs(Project *pr)
/*
**--------------------------------------------------------------
** Input: none
@@ -293,10 +304,10 @@ void valvecoeffs(EN_Project *pr)
**--------------------------------------------------------------
*/
{
int i, k, n1, n2;
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
hydraulics_t *hyd = &pr->hydraulics;
EN_Network *net = &pr->network;
int i, k, n1, n2;
Slink *link;
Svalve *valve;
@@ -333,7 +344,7 @@ void valvecoeffs(EN_Project *pr)
}
void emittercoeffs(EN_Project *pr)
void emittercoeffs(Project *pr)
/*
**--------------------------------------------------------------
** Input: none
@@ -348,13 +359,13 @@ void emittercoeffs(EN_Project *pr)
**--------------------------------------------------------------
*/
{
int i, row;
double hloss, hgrad;
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Smatrix *sm = &hyd->smatrix;
hydraulics_t *hyd = &pr->hydraulics;
solver_t *sol = &hyd->solver;
EN_Network *net = &pr->network;
Snode *node;
int i, row;
double hloss, hgrad;
Snode *node;
for (i = 1; i <= net->Njuncs; i++)
{
@@ -366,19 +377,19 @@ void emittercoeffs(EN_Project *pr)
emitheadloss(pr, i, &hloss, &hgrad);
// Row of solution matrix
row = sol->Row[i];
row = sm->Row[i];
// Addition to matrix diagonal & r.h.s
sol->Aii[row] += 1.0 / hgrad;
sol->F[row] += (hloss + node->El) / hgrad;
sm->Aii[row] += 1.0 / hgrad;
sm->F[row] += (hloss + node->El) / hgrad;
// Update to node flow balance
hyd->X_tmp[i] -= hyd->EmitterFlows[i];
// Update to node flow excess
hyd->Xflow[i] -= hyd->EmitterFlow[i];
}
}
void emitheadloss(EN_Project *pr, int i, double *hloss, double *hgrad)
void emitheadloss(Project *pr, int i, double *hloss, double *hgrad)
/*
**-------------------------------------------------------------
** Input: i = node index
@@ -388,10 +399,11 @@ void emitheadloss(EN_Project *pr, int i, double *hloss, double *hgrad)
**-------------------------------------------------------------
*/
{
Hydraul *hyd = &pr->hydraul;
double ke;
double q;
double qa;
hydraulics_t *hyd = &pr->hydraulics;
// Set adjusted emitter coeff.
ke = MAX(CSMALL, pr->network.Node[i].Ke);
@@ -400,7 +412,7 @@ void emitheadloss(EN_Project *pr, int i, double *hloss, double *hgrad)
qa = pow(hyd->RQtol / ke / hyd->Qexp, 1.0 / (hyd->Qexp - 1.0));
// Use linear head loss relation for small flow
q = hyd->EmitterFlows[i];
q = hyd->EmitterFlow[i];
if (fabs(q) <= qa)
{
*hgrad = hyd->RQtol;
@@ -416,7 +428,7 @@ void emitheadloss(EN_Project *pr, int i, double *hloss, double *hgrad)
}
void demandparams(EN_Project *pr, double *dp, double *n)
void demandparams(Project *pr, double *dp, double *n)
/*
**--------------------------------------------------------------
** Input: none
@@ -427,7 +439,7 @@ void demandparams(EN_Project *pr, double *dp, double *n)
**--------------------------------------------------------------
*/
{
hydraulics_t *hyd = &pr->hydraulics;
Hydraul *hyd = &pr->hydraul;
// If required pressure equals minimum pressure, use a linear demand
// curve with a 0.01 PSI pressure range to approximate an all or
@@ -447,7 +459,7 @@ void demandparams(EN_Project *pr, double *dp, double *n)
}
void demandcoeffs(EN_Project *pr)
void demandcoeffs(Project *pr)
/*
**--------------------------------------------------------------
** Input: none
@@ -462,16 +474,16 @@ void demandcoeffs(EN_Project *pr)
**--------------------------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Smatrix *sm = &hyd->smatrix;
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);
@@ -483,18 +495,18 @@ void demandcoeffs(EN_Project *pr)
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,
demandheadloss(hyd->DemandFlow[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;
row = sm->Row[i];
sm->Aii[row] += 1.0 / hgrad;
sm->F[row] += (hloss + net->Node[i].El + hyd->Pmin) / hgrad;
}
}
double demandflowchange(EN_Project *pr, int i, double dp, double n)
double demandflowchange(Project *pr, int i, double dp, double n)
/*
**--------------------------------------------------------------
** Input: i = node index
@@ -506,16 +518,17 @@ double demandflowchange(EN_Project *pr, int i, double dp, double n)
**--------------------------------------------------------------
*/
{
double hloss, hgrad;
hydraulics_t *hyd = &pr->hydraulics;
Hydraul *hyd = &pr->hydraul;
demandheadloss(hyd->DemandFlows[i], hyd->NodeDemand[i], dp, n, &hloss, &hgrad);
double hloss, hgrad;
demandheadloss(hyd->DemandFlow[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)
double *hloss, double *hgrad)
/*
**--------------------------------------------------------------
** Input: d = actual junction demand (cfs)
@@ -563,7 +576,7 @@ void demandheadloss(double d, double dfull, double dp, double n,
}
void pipecoeff(EN_Project *pr, int k)
void pipecoeff(Project *pr, int k)
/*
**--------------------------------------------------------------
** Input: k = link index
@@ -575,20 +588,19 @@ void pipecoeff(EN_Project *pr, int k)
**--------------------------------------------------------------
*/
{
Hydraul *hyd = &pr->hydraul;
double hloss, // Head loss
hgrad, // Head loss gradient
ml, // Minor loss coeff.
q, // Abs. value of flow
r; // Resistance coeff.
hydraulics_t *hyd = &pr->hydraulics;
solver_t *sol = &hyd->solver;
// For closed pipe use headloss formula: hloss = CBIG*q
if (hyd->LinkStatus[k] <= CLOSED)
{
sol->P[k] = 1.0 / CBIG;
sol->Y[k] = hyd->LinkFlows[k];
hyd->P[k] = 1.0 / CBIG;
hyd->Y[k] = hyd->LinkFlow[k];
return;
}
@@ -599,7 +611,7 @@ void pipecoeff(EN_Project *pr, int k)
return;
}
q = ABS(hyd->LinkFlows[k]);
q = ABS(hyd->LinkFlow[k]);
ml = pr->network.Link[k].Km;
r = pr->network.Link[k].R;
@@ -625,15 +637,15 @@ void pipecoeff(EN_Project *pr, int k)
}
// Adjust head loss sign for flow direction
hloss *= SGN(hyd->LinkFlows[k]);
hloss *= SGN(hyd->LinkFlow[k]);
// P and Y coeffs.
sol->P[k] = 1.0 / hgrad;
sol->Y[k] = hloss / hgrad;
hyd->P[k] = 1.0 / hgrad;
hyd->Y[k] = hloss / hgrad;
}
void DWpipecoeff(EN_Project *pr, int k)
void DWpipecoeff(Project *pr, int k)
/*
**--------------------------------------------------------------
** Input: k = link index
@@ -643,16 +655,14 @@ void DWpipecoeff(EN_Project *pr, int k)
**--------------------------------------------------------------
*/
{
hydraulics_t *hyd = &pr->hydraulics;
solver_t *sol = &hyd->solver;
Slink *link = &pr->network.Link[k];
Hydraul *hyd = &pr->hydraul;
Slink *link = &pr->network.Link[k];
double q = ABS(hyd->LinkFlows[k]);
double q = ABS(hyd->LinkFlow[k]);
double r = link->R; // Resistance coeff.
double ml = link->Km; // Minor loss coeff.
double e = link->Kc / link->Diam; // Relative roughness
double s = hyd->Viscos * link->Diam; // Viscosity / diameter
double hloss, hgrad, f, dfdq, r1;
// Compute head loss and its derivative
@@ -660,7 +670,7 @@ void DWpipecoeff(EN_Project *pr, int k)
if (q <= A2 * s)
{
r = 16.0 * PI * s * r;
hloss = hyd->LinkFlows[k] * (r + ml * q);
hloss = hyd->LinkFlow[k] * (r + ml * q);
hgrad = r + 2.0 * ml * q;
}
@@ -670,13 +680,13 @@ void DWpipecoeff(EN_Project *pr, int k)
dfdq = 0.0;
f = frictionFactor(q, e, s, &dfdq);
r1 = f * r + ml;
hloss = r1 * q * hyd->LinkFlows[k];
hloss = r1 * q * hyd->LinkFlow[k];
hgrad = (2.0 * r1 * q) + (dfdq * r * q * q);
}
// Compute P and Y coefficients
sol->P[k] = 1.0 / hgrad;
sol->Y[k] = hloss / hgrad;
hyd->P[k] = 1.0 / hgrad;
hyd->Y[k] = hloss / hgrad;
}
@@ -730,7 +740,7 @@ double frictionFactor(double q, double e, double s, double *dfdq)
}
void pumpcoeff(EN_Project *pr, int k)
void pumpcoeff(Project *pr, int k)
/*
**--------------------------------------------------------------
** Input: k = link index
@@ -739,6 +749,8 @@ void pumpcoeff(EN_Project *pr, int k)
**--------------------------------------------------------------
*/
{
Hydraul *hyd = &pr->hydraul;
int p; // Pump index
double h0, // Shutoff head
q, // Abs. value of flow
@@ -748,22 +760,19 @@ void pumpcoeff(EN_Project *pr, int k)
qa, // Flow limit for linear head loss
hloss, // Head loss across pump
hgrad; // Head loss gradient
hydraulics_t *hyd = &pr->hydraulics;
solver_t *sol = &hyd->solver;
Spump *pump;
Spump *pump;
// Use high resistance pipe if pump closed or cannot deliver head
setting = hyd->LinkSetting[k];
if (hyd->LinkStatus[k] <= CLOSED || setting == 0.0)
{
sol->P[k] = 1.0 / CBIG;
sol->Y[k] = hyd->LinkFlows[k];
hyd->P[k] = 1.0 / CBIG;
hyd->Y[k] = hyd->LinkFlow[k];
return;
}
// Obtain reference to pump object & its speed setting
q = ABS(hyd->LinkFlows[k]);
q = ABS(hyd->LinkFlow[k]);
p = findpump(&pr->network, k);
pump = &pr->network.Pump[p];
@@ -783,7 +792,7 @@ void pumpcoeff(EN_Project *pr, int k)
// Compute head loss and its gradient
hgrad = pump->R * setting ;
hloss = pump->H0 * SQR(setting) + hgrad * hyd->LinkFlows[k];
hloss = pump->H0 * SQR(setting) + hgrad * hyd->LinkFlow[k];
}
else
{
@@ -798,23 +807,23 @@ void pumpcoeff(EN_Project *pr, int k)
if (q <= qa)
{
hgrad = hyd->RQtol;
hloss = h0 + hgrad * hyd->LinkFlows[k];
hloss = h0 + hgrad * hyd->LinkFlow[k];
}
// ... use original pump curve for normal flows
else
{
hgrad = n * r * pow(q, n - 1.0);
hloss = h0 + hgrad * hyd->LinkFlows[k] / n;
hloss = h0 + hgrad * hyd->LinkFlow[k] / n;
}
}
// P and Y coeffs.
sol->P[k] = 1.0 / hgrad;
sol->Y[k] = hloss / hgrad;
hyd->P[k] = 1.0 / hgrad;
hyd->Y[k] = hloss / hgrad;
}
void curvecoeff(EN_Project *pr, int i, double q, double *h0, double *r)
void curvecoeff(Project *pr, int i, double q, double *h0, double *r)
/*
**-------------------------------------------------------------------
** Input: i = curve index
@@ -854,7 +863,7 @@ void curvecoeff(EN_Project *pr, int i, double q, double *h0, double *r)
}
void gpvcoeff(EN_Project *pr, int k)
void gpvcoeff(Project *pr, int k)
/*
**--------------------------------------------------------------
** Input: k = link index
@@ -868,8 +877,7 @@ void gpvcoeff(EN_Project *pr, int k)
r, // Slope of head loss curve segment
q; // Abs. value of flow
hydraulics_t *hyd = &pr->hydraulics;
solver_t *sol = &hyd->solver;
Hydraul *hyd = &pr->hydraul;
// Treat as a pipe if valve closed
if (hyd->LinkStatus[k] == CLOSED) valvecoeff(pr, k);
@@ -883,7 +891,7 @@ void gpvcoeff(EN_Project *pr, int k)
i = (int)ROUND(hyd->LinkSetting[k]);
// Adjusted flow rate
q = ABS(hyd->LinkFlows[k]);
q = ABS(hyd->LinkFlow[k]);
q = MAX(q, TINY);
// Intercept and slope of curve segment containing q
@@ -891,13 +899,13 @@ void gpvcoeff(EN_Project *pr, int k)
r = MAX(r, TINY);
// Resulting P and Y coeffs.
sol->P[k] = 1.0 / r;
sol->Y[k] = (h0 / r + q) * SGN(hyd->LinkFlows[k]);
hyd->P[k] = 1.0 / r;
hyd->Y[k] = (h0 / r + q) * SGN(hyd->LinkFlow[k]);
}
}
void pbvcoeff(EN_Project *pr, int k)
void pbvcoeff(Project *pr, int k)
/*
**--------------------------------------------------------------
** Input: k = link index
@@ -906,8 +914,7 @@ void pbvcoeff(EN_Project *pr, int k)
**--------------------------------------------------------------
*/
{
hydraulics_t *hyd = &pr->hydraulics;
solver_t *sol = &hyd->solver;
Hydraul *hyd = &pr->hydraul;
Slink *link = &pr->network.Link[k];
// If valve fixed OPEN or CLOSED then treat as a pipe
@@ -920,21 +927,21 @@ void pbvcoeff(EN_Project *pr, int k)
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->LinkFlow[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;
hyd->P[k] = CBIG;
hyd->Y[k] = hyd->LinkSetting[k] * CBIG;
}
}
}
void tcvcoeff(EN_Project *pr, int k)
void tcvcoeff(Project *pr, int k)
/*
**--------------------------------------------------------------
** Input: k = link index
@@ -944,7 +951,7 @@ void tcvcoeff(EN_Project *pr, int k)
*/
{
double km;
hydraulics_t *hyd = &pr->hydraulics;
Hydraul *hyd = &pr->hydraul;
Slink *link = &pr->network.Link[k];
// Save original loss coeff. for open valve
@@ -964,7 +971,7 @@ void tcvcoeff(EN_Project *pr, int k)
}
void prvcoeff(EN_Project *pr, int k, int n1, int n2)
void prvcoeff(Project *pr, int k, int n1, int n2)
/*
**--------------------------------------------------------------
** Input: k = link index
@@ -976,13 +983,14 @@ void prvcoeff(EN_Project *pr, int k, int n1, int n2)
**--------------------------------------------------------------
*/
{
Hydraul *hyd = &pr->hydraul;
Smatrix *sm = &hyd->smatrix;
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];
i = sm->Row[n1]; // Matrix rows of nodes
j = sm->Row[n2];
hset = pr->network.Node[n2].El +
hyd->LinkSetting[k]; // Valve setting
@@ -991,15 +999,15 @@ void prvcoeff(EN_Project *pr, int k, int n1, int n2)
// Set coeffs. to force head at downstream
// node equal to valve setting & force flow
// to equal to flow imbalance at downstream node.
// to equal to flow excess 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)
hyd->P[k] = 0.0;
hyd->Y[k] = hyd->LinkFlow[k] + hyd->Xflow[n2]; // Force flow balance
sm->F[j] += (hset * CBIG); // Force head = hset
sm->Aii[j] += CBIG; // at downstream node
if (hyd->Xflow[n2] < 0.0)
{
sol->F[i] += hyd->X_tmp[n2];
sm->F[i] += hyd->Xflow[n2];
}
return;
}
@@ -1008,15 +1016,15 @@ void prvcoeff(EN_Project *pr, int k, int n1, int n2)
// 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]);
sm->Aij[sm->Ndx[k]] -= hyd->P[k];
sm->Aii[i] += hyd->P[k];
sm->Aii[j] += hyd->P[k];
sm->F[i] += (hyd->Y[k] - hyd->LinkFlow[k]);
sm->F[j] -= (hyd->Y[k] - hyd->LinkFlow[k]);
}
void psvcoeff(EN_Project *pr, int k, int n1, int n2)
void psvcoeff(Project *pr, int k, int n1, int n2)
/*
**--------------------------------------------------------------
** Input: k = link index
@@ -1028,13 +1036,14 @@ void psvcoeff(EN_Project *pr, int k, int n1, int n2)
**--------------------------------------------------------------
*/
{
Hydraul *hyd = &pr->hydraul;
Smatrix *sm = &hyd->smatrix;
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];
i = sm->Row[n1]; // Matrix rows of nodes
j = sm->Row[n2];
hset = pr->network.Node[n1].El +
hyd->LinkSetting[k]; // Valve setting
@@ -1042,15 +1051,15 @@ void psvcoeff(EN_Project *pr, int k, int n1, int n2)
{
// Set coeffs. to force head at upstream
// node equal to valve setting & force flow
// equal to flow imbalance at upstream node.
// equal to flow excess 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)
hyd->P[k] = 0.0;
hyd->Y[k] = hyd->LinkFlow[k] - hyd->Xflow[n1]; // Force flow balance
sm->F[i] += (hset * CBIG); // Force head = hset
sm->Aii[i] += CBIG; // at upstream node
if (hyd->Xflow[n1] > 0.0)
{
sol->F[j] += hyd->X_tmp[n1];
sm->F[j] += hyd->Xflow[n1];
}
return;
}
@@ -1059,15 +1068,15 @@ void psvcoeff(EN_Project *pr, int k, int n1, int n2)
// 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]);
sm->Aij[sm->Ndx[k]] -= hyd->P[k];
sm->Aii[i] += hyd->P[k];
sm->Aii[j] += hyd->P[k];
sm->F[i] += (hyd->Y[k] - hyd->LinkFlow[k]);
sm->F[j] -= (hyd->Y[k] - hyd->LinkFlow[k]);
}
void fcvcoeff(EN_Project *pr, int k, int n1, int n2)
void fcvcoeff(Project *pr, int k, int n1, int n2)
/*
**--------------------------------------------------------------
** Input: k = link index
@@ -1079,14 +1088,15 @@ void fcvcoeff(EN_Project *pr, int k, int n1, int n2)
**--------------------------------------------------------------
*/
{
Hydraul *hyd = &pr->hydraul;
Smatrix *sm = &hyd->smatrix;
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];
i = sm->Row[n1];
j = sm->Row[n2];
// If valve active, break network at valve and treat
// flow setting as external demand at upstream node
@@ -1094,15 +1104,15 @@ void fcvcoeff(EN_Project *pr, int k, int n1, int n2)
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;
hyd->Xflow[n1] -= q;
hyd->Xflow[n2] += q;
hyd->Y[k] = hyd->LinkFlow[k] - q;
sm->F[i] -= q;
sm->F[j] += q;
hyd->P[k] = 1.0 / CBIG;
sm->Aij[sm->Ndx[k]] -= hyd->P[k];
sm->Aii[i] += hyd->P[k];
sm->Aii[j] += hyd->P[k];
}
// Otherwise treat valve as an open pipe
@@ -1110,16 +1120,16 @@ void fcvcoeff(EN_Project *pr, int k, int n1, int n2)
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]);
sm->Aij[sm->Ndx[k]] -= hyd->P[k];
sm->Aii[i] += hyd->P[k];
sm->Aii[j] += hyd->P[k];
sm->F[i] += (hyd->Y[k] - hyd->LinkFlow[k]);
sm->F[j] -= (hyd->Y[k] - hyd->LinkFlow[k]);
}
}
void valvecoeff(EN_Project *pr, int k)
void valvecoeff(Project *pr, int k)
/*
**--------------------------------------------------------------
** Input: k = link index
@@ -1129,20 +1139,18 @@ void valvecoeff(EN_Project *pr, int k)
**--------------------------------------------------------------
*/
{
Hydraul *hyd = &pr->hydraul;
Slink *link = &pr->network.Link[k];
double flow, q, y, qa, hgrad;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
solver_t *sol = &hyd->solver;
Slink *link = &net->Link[k];
flow = hyd->LinkFlows[k];
flow = hyd->LinkFlow[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;
hyd->P[k] = 1.0 / CBIG;
hyd->Y[k] = flow;
return;
}
@@ -1164,15 +1172,15 @@ void valvecoeff(EN_Project *pr, int k)
}
// P and Y coeffs.
sol->P[k] = 1.0 / hgrad;
sol->Y[k] = y;
hyd->P[k] = 1.0 / hgrad;
hyd->Y[k] = y;
}
// If no minor loss coeff. specified use a
// low resistance linear head loss relation
else
{
sol->P[k] = 1.0 / CSMALL;
sol->Y[k] = flow;
hyd->P[k] = 1.0 / CSMALL;
hyd->Y[k] = flow;
}
}
src/hydraul.c
+768 -850
View File
File diff suppressed because it is too large Load Diff
src/hydsolver.c
+105 -124
View File
@@ -1,13 +1,15 @@
/*
*********************************************************************
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.
*******************************************************************
******************************************************************************
Project: OWA EPANET
Version: 2.2
Module: hydsolver.c
Description: computes flows and pressures throughout a pipe network using
Todini's Global Gradient Algorithm
Authors: see AUTHORS
Copyright: see AUTHORS
License: see LICENSE
Last Updated: 11/27/2018
******************************************************************************
*/
#include <stdio.h>
@@ -33,32 +35,29 @@ typedef struct {
int maxflowlink;
} Hydbalance;
// External functions
//int hydsolve(EN_Project *pr, int *iter, double *relerr);
//void headlosscoeffs(EN_Project *pr);
//void matrixcoeffs(EN_Project *pr);
// Exported functions
int hydsolve(Project *, int *, double *);
extern int valvestatus(EN_Project *pr); //(see HYDSTATUS.C)
extern int linkstatus(EN_Project *pr); //(see HYDSTATUS.C)
// Imported functions
extern int linsolve(Smatrix *, int); //(see SMATRIX.C)
extern int valvestatus(Project *); //(see HYDSTATUS.C)
extern int linkstatus(Project *); //(see HYDSTATUS.C)
// Local functions
static int badvalve(EN_Project *pr, int);
static int pswitch(EN_Project *pr);
static int badvalve(Project *, int);
static int pswitch(Project *);
static double newflows(EN_Project *pr, Hydbalance *hbal);
static void newlinkflows(EN_Project *pr, Hydbalance *hbal, double *qsum,
double *dqsum);
static void newemitterflows(EN_Project *pr, Hydbalance *hbal, double *qsum,
double *dqsum);
static void newdemandflows(EN_Project *pr, Hydbalance *hbal, double *qsum,
double *dqsum);
static double newflows(Project *, Hydbalance *);
static void newlinkflows(Project *, Hydbalance *, double *, double *);
static void newemitterflows(Project *, Hydbalance *, double *, double *);
static void newdemandflows(Project *, Hydbalance *, double *, double *);
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);
static void checkhydbalance(Project *, Hydbalance *);
static int hasconverged(Project *, double *, Hydbalance *);
static void reporthydbal(Project *, Hydbalance *);
int hydsolve(EN_Project *pr, int *iter, double *relerr)
int hydsolve(Project *pr, int *iter, double *relerr)
/*
**-------------------------------------------------------------------
** Input: none
@@ -68,11 +67,8 @@ int hydsolve(EN_Project *pr, int *iter, double *relerr)
** 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
** every 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
@@ -87,6 +83,11 @@ int hydsolve(EN_Project *pr, int *iter, double *relerr)
**-------------------------------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Smatrix *sm = &hyd->smatrix;
Report *rpt = &pr->report;
int i; // Node index
int errcode = 0; // Node causing solution error
int nextcheck; // Next status check trial
@@ -96,47 +97,32 @@ int hydsolve(EN_Project *pr, int *iter, double *relerr)
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.
// (ExtraIter used to increase trials in case of status cycling.)
if (pr->report.Statflag == FULL)
{
writerelerr(pr, 0, 0);
}
if (rpt->Statflag == FULL) writerelerr(pr, 0, 0);
maxtrials = hyd->MaxIter;
if (hyd->ExtraIter > 0)
{
maxtrials += hyd->ExtraIter;
}
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().
*/
// 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().
headlosscoeffs(pr);
matrixcoeffs(pr);
errcode = linsolve(pr, net->Njuncs);
// Quit if memory allocation error
if (errcode < 0) break;
errcode = linsolve(sm, net->Njuncs);
// Matrix ill-conditioning problem - if control valve causing problem,
// fix its status & continue, otherwise quit with no solution.
if (errcode > 0)
{
if (badvalve(pr, sol->Order[errcode])) continue;
if (badvalve(pr, sm->Order[errcode])) continue;
else break;
}
@@ -144,13 +130,13 @@ int hydsolve(EN_Project *pr, int *iter, double *relerr)
// (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
hyd->NodeHead[i] = sm->F[sm->Row[i]]; // Update heads
}
newerr = newflows(pr, &hydbal); // Update flows
newerr = newflows(pr, &hydbal); // Update flows
*relerr = newerr;
// Write convergence error to status report if called for
if (rep->Statflag == FULL)
if (rpt->Statflag == FULL)
{
writerelerr(pr, *iter, *relerr);
}
@@ -199,17 +185,16 @@ int hydsolve(EN_Project *pr, int *iter, double *relerr)
}
// Iterations ended - report any errors.
if (errcode < 0) errcode = 101; // Memory allocation error
else if (errcode > 0)
if (errcode > 0)
{
writehyderr(pr, sol->Order[errcode]); // Ill-conditioned matrix error
writehyderr(pr, sm->Order[errcode]); // Ill-conditioned matrix error
errcode = 110;
}
// Replace junction demands with total outflow values
for (i = 1; i <= net->Njuncs; i++)
{
hyd->NodeDemand[i] = hyd->DemandFlows[i] + hyd->EmitterFlows[i];
hyd->NodeDemand[i] = hyd->DemandFlow[i] + hyd->EmitterFlow[i];
}
// Save convergence info
@@ -218,11 +203,11 @@ int hydsolve(EN_Project *pr, int *iter, double *relerr)
hyd->MaxFlowChange = hydbal.maxflowchange;
hyd->Iterations = *iter;
return(errcode);
return errcode;
}
int badvalve(EN_Project *pr, int n)
int badvalve(Project *pr, int n)
/*
**-----------------------------------------------------------------
** Input: n = node index
@@ -235,12 +220,12 @@ int badvalve(EN_Project *pr, int n)
**-----------------------------------------------------------------
*/
{
int i, k, n1, n2;
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Report *rpt = &pr->report;
Times *time = &pr->times;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
report_options_t *rep = &pr->report;
time_options_t *time = &pr->time_options;
int i, k, n1, n2;
Slink *link;
LinkType t;
@@ -257,19 +242,14 @@ int badvalve(EN_Project *pr, int n)
{
if (hyd->LinkStatus[k] == ACTIVE)
{
if (rep->Statflag == FULL)
if (rpt->Statflag == FULL)
{
sprintf(pr->Msg, FMT61, clocktime(rep->Atime, time->Htime), link->ID);
sprintf(pr->Msg, FMT61,
clocktime(rpt->Atime, time->Htime), link->ID);
writeline(pr, pr->Msg);
}
if (link->Type == FCV)
{
hyd->LinkStatus[k] = XFCV;
}
else
{
hyd->LinkStatus[k] = XPRESSURE;
}
if (link->Type == FCV) hyd->LinkStatus[k] = XFCV;
else hyd->LinkStatus[k] = XPRESSURE;
return 1;
}
}
@@ -280,7 +260,7 @@ int badvalve(EN_Project *pr, int n)
}
int pswitch(EN_Project *pr)
int pswitch(Project *pr)
/*
**--------------------------------------------------------------
** Input: none
@@ -290,6 +270,10 @@ int pswitch(EN_Project *pr)
**--------------------------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Report *rpt = &pr->report;
int i, // Control statement index
k, // Index of link being controlled
n, // Node controlling link k
@@ -297,10 +281,6 @@ int pswitch(EN_Project *pr)
change, // Flag for status or setting change
anychange = 0; // Flag for 1 or more control actions
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
@@ -358,7 +338,7 @@ int pswitch(EN_Project *pr)
{
hyd->LinkSetting[k] = net->Control[i].Setting;
}
if (rep->Statflag == FULL)
if (rpt->Statflag == FULL)
{
writestatchange(pr, k, s, hyd->LinkStatus[k]);
}
@@ -366,11 +346,11 @@ int pswitch(EN_Project *pr)
}
}
}
return(anychange);
return anychange;
}
double newflows(EN_Project *pr, Hydbalance *hbal)
double newflows(Project *pr, Hydbalance *hbal)
/*
**----------------------------------------------------------------
** Input: hbal = ptr. to hydraulic balance information
@@ -380,12 +360,12 @@ double newflows(EN_Project *pr, Hydbalance *hbal)
**----------------------------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
double dqsum, // Network flow change
qsum; // Network total flow
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
// Initialize sum of flows & corrections
qsum = 0.0;
dqsum = 0.0;
@@ -404,7 +384,7 @@ double newflows(EN_Project *pr, Hydbalance *hbal)
}
void newlinkflows(EN_Project *pr, Hydbalance *hbal, double *qsum, double *dqsum)
void newlinkflows(Project *pr, Hydbalance *hbal, double *qsum, double *dqsum)
/*
**----------------------------------------------------------------
** Input: hbal = ptr. to hydraulic balance information
@@ -415,14 +395,13 @@ void newlinkflows(EN_Project *pr, Hydbalance *hbal, double *qsum, double *dqsum
**----------------------------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
double dh, /* Link head loss */
dq; /* Link flow change */
int k, n, n1, n2;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
solver_t *sol = &hyd->solver;
Slink *link;
Slink *link;
// Initialize net inflows (i.e., demands) at fixed grade nodes
for (n = net->Njuncs + 1; n <= net->Nnodes; n++)
@@ -445,7 +424,7 @@ void newlinkflows(EN_Project *pr, Hydbalance *hbal, double *qsum, double *dqsum
// where P & Y were computed in hlosscoeff() in hydcoeffs.c
dh = hyd->NodeHead[n1] - hyd->NodeHead[n2];
dq = sol->Y[k] - sol->P[k] * dh;
dq = hyd->Y[k] - hyd->P[k] * dh;
// Adjust flow change by the relaxation factor
dq *= hyd->RelaxFactor;
@@ -454,15 +433,15 @@ void newlinkflows(EN_Project *pr, Hydbalance *hbal, double *qsum, double *dqsum
if (link->Type == PUMP)
{
n = findpump(net, k);
if (net->Pump[n].Ptype == CONST_HP && dq > hyd->LinkFlows[k])
if (net->Pump[n].Ptype == CONST_HP && dq > hyd->LinkFlow[k])
{
dq = hyd->LinkFlows[k] / 2.0;
dq = hyd->LinkFlow[k] / 2.0;
}
}
// Update link flow and system flow summation
hyd->LinkFlows[k] -= dq;
*qsum += ABS(hyd->LinkFlows[k]);
hyd->LinkFlow[k] -= dq;
*qsum += ABS(hyd->LinkFlow[k]);
*dqsum += ABS(dq);
// Update identity of element with max. flow change
@@ -476,14 +455,14 @@ void newlinkflows(EN_Project *pr, Hydbalance *hbal, double *qsum, double *dqsum
// Update net flows to fixed grade nodes
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];
if (n1 > net->Njuncs) hyd->NodeDemand[n1] -= hyd->LinkFlow[k];
if (n2 > net->Njuncs) hyd->NodeDemand[n2] += hyd->LinkFlow[k];
}
}
}
void newemitterflows(EN_Project *pr, Hydbalance *hbal, double *qsum,
void newemitterflows(Project *pr, Hydbalance *hbal, double *qsum,
double *dqsum)
/*
**----------------------------------------------------------------
@@ -495,10 +474,11 @@ void newemitterflows(EN_Project *pr, Hydbalance *hbal, double *qsum,
**----------------------------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
int i;
double hloss, hgrad, dh, dq;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
// Examine each network junction
for (i = 1; i <= net->Njuncs; i++)
@@ -512,10 +492,10 @@ void newemitterflows(EN_Project *pr, Hydbalance *hbal, double *qsum,
// Find emitter flow change
dh = hyd->NodeHead[i] - net->Node[i].El;
dq = (hloss - dh) / hgrad;
hyd->EmitterFlows[i] -= dq;
hyd->EmitterFlow[i] -= dq;
// Update system flow summation
*qsum += ABS(hyd->EmitterFlows[i]);
*qsum += ABS(hyd->EmitterFlow[i]);
*dqsum += ABS(dq);
// Update identity of element with max. flow change
@@ -529,7 +509,7 @@ void newemitterflows(EN_Project *pr, Hydbalance *hbal, double *qsum,
}
void newdemandflows(EN_Project *pr, Hydbalance *hbal, double *qsum, double *dqsum)
void newdemandflows(Project *pr, Hydbalance *hbal, double *qsum, double *dqsum)
/*
**----------------------------------------------------------------
** Input: hbal = ptr. to hydraulic balance information
@@ -541,12 +521,13 @@ void newdemandflows(EN_Project *pr, Hydbalance *hbal, double *qsum, double *dqsu
**----------------------------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
double dp, // pressure range over which demand can vary (ft)
dq, // change in demand flow (cfs)
n; // exponent in head loss v. demand function
int k;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
// Get demand function parameters
if (hyd->DemandModel == DDA) return;
@@ -560,10 +541,10 @@ void newdemandflows(EN_Project *pr, Hydbalance *hbal, double *qsum, double *dqsu
// Find change in demand flow (see hydcoeffs.c)
dq = demandflowchange(pr, k, dp, n);
hyd->DemandFlows[k] -= dq;
hyd->DemandFlow[k] -= dq;
// Update system flow summation
*qsum += ABS(hyd->DemandFlows[k]);
*qsum += ABS(hyd->DemandFlow[k]);
*dqsum += ABS(dq);
// Update identity of element with max. flow change
@@ -577,7 +558,7 @@ void newdemandflows(EN_Project *pr, Hydbalance *hbal, double *qsum, double *dqsu
}
void checkhydbalance(EN_Project *pr, Hydbalance *hbal)
void checkhydbalance(Project *pr, Hydbalance *hbal)
/*
**--------------------------------------------------------------
** Input: hbal = hydraulic balance errors
@@ -586,25 +567,25 @@ void checkhydbalance(EN_Project *pr, Hydbalance *hbal)
**--------------------------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
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;
headlosscoeffs(pr);
for (k = 1; k <= net->Nlinks; k++)
{
if (hyd->LinkStatus[k] <= CLOSED) continue;
if (sol->P[k] == 0.0) continue;
if (hyd->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];
headloss = hyd->Y[k] / hyd->P[k];
headerror = ABS(dh - headloss);
if (headerror > hbal->maxheaderror)
{
@@ -615,7 +596,7 @@ void checkhydbalance(EN_Project *pr, Hydbalance *hbal)
}
int hasconverged(EN_Project *pr, double *relerr, Hydbalance *hbal)
int hasconverged(Project *pr, double *relerr, Hydbalance *hbal)
/*
**--------------------------------------------------------------
** Input: relerr = current total relative flow change
@@ -626,7 +607,7 @@ int hasconverged(EN_Project *pr, double *relerr, Hydbalance *hbal)
**--------------------------------------------------------------
*/
{
hydraulics_t *hyd = &pr->hydraulics;
Hydraul *hyd = &pr->hydraul;
if (*relerr > hyd->Hacc) return 0;
checkhydbalance(pr, hbal);
@@ -642,7 +623,7 @@ int hasconverged(EN_Project *pr, double *relerr, Hydbalance *hbal)
}
void reporthydbal(EN_Project *pr, Hydbalance *hbal)
void reporthydbal(Project *pr, Hydbalance *hbal)
/*
**--------------------------------------------------------------
** Input: hbal = current hydraulic balance errors
src/hydstatus.c
+81 -69
View File
@@ -1,29 +1,34 @@
/*
*********************************************************************
HYDSTATUS.C -- Hydraulic status updating for the EPANET Program
*******************************************************************
******************************************************************************
Project: OWA EPANET
Version: 2.2
Module: hydstatus.c
Description: updates hydraulic status of network elements
Authors: see AUTHORS
Copyright: see AUTHORS
License: see LICENSE
Last Updated: 11/27/2018
******************************************************************************
*/
#include <stdio.h>
#include "types.h"
#include "funcs.h"
// External functions
int valvestatus(EN_Project *pr);
int linkstatus(EN_Project *pr);
// Exported functions
int valvestatus(Project *);
int linkstatus(Project *);
// 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);
static StatusType cvstatus(Project *, StatusType, double, double);
static StatusType pumpstatus(Project *, int, double);
static StatusType prvstatus(Project *, int, StatusType, double, double, double);
static StatusType psvstatus(Project *, int, StatusType, double, double, double);
static StatusType fcvstatus(Project *, int, StatusType, double, double);
static void tankstatus(Project *, int, int, int);
int valvestatus(EN_Project *pr)
int valvestatus(Project *pr)
/*
**-----------------------------------------------------------------
** Input: none
@@ -34,16 +39,16 @@ int valvestatus(EN_Project *pr)
**-----------------------------------------------------------------
*/
{
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;
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Report *rpt = &pr->report;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
report_options_t *rep = &pr->report;
int change = FALSE, // Status change flag
i, k, // Valve & link indexes
n1, n2; // Start & end nodes
double hset; // Valve head setting
StatusType status; // Valve status settings
Slink *link;
// Examine each valve
for (i = 1; i <= net->Nvalves; i++)
@@ -80,7 +85,7 @@ int valvestatus(EN_Project *pr)
// Check for a status change
if (status != hyd->LinkStatus[k])
{
if (rep->Statflag == FULL)
if (rpt->Statflag == FULL)
{
writestatchange(pr, k, status, hyd->LinkStatus[k]);
}
@@ -88,10 +93,10 @@ int valvestatus(EN_Project *pr)
}
}
return change;
} /* End of valvestatus() */
}
int linkstatus(EN_Project *pr)
int linkstatus(Project *pr)
/*
**--------------------------------------------------------------
** Input: none
@@ -101,16 +106,16 @@ int linkstatus(EN_Project *pr)
**--------------------------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Report *rpt = &pr->report;
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;
StatusType status; // Current status
Slink *link;
// Examine each link
@@ -132,7 +137,7 @@ int linkstatus(EN_Project *pr)
if (link->Type == CVPIPE)
{
hyd->LinkStatus[k] = cvstatus(pr, hyd->LinkStatus[k], dh,
hyd->LinkFlows[k]);
hyd->LinkFlow[k]);
}
if (link->Type == PUMP && hyd->LinkStatus[k] >= OPEN &&
hyd->LinkSetting[k] > 0.0)
@@ -157,7 +162,7 @@ int linkstatus(EN_Project *pr)
if (status != hyd->LinkStatus[k])
{
change = TRUE;
if (rep->Statflag == FULL)
if (rpt->Statflag == FULL)
{
writestatchange(pr, k, status, hyd->LinkStatus[k]);
}
@@ -167,7 +172,7 @@ int linkstatus(EN_Project *pr)
}
StatType cvstatus(EN_Project *pr, StatType s, double dh, double q)
StatusType cvstatus(Project *pr, StatusType s, double dh, double q)
/*
**--------------------------------------------------
** Input: s = current link status
@@ -178,7 +183,7 @@ StatType cvstatus(EN_Project *pr, StatType s, double dh, double q)
**--------------------------------------------------
*/
{
hydraulics_t *hyd = &pr->hydraulics;
Hydraul *hyd = &pr->hydraul;
// Prevent reverse flow through CVs
if (ABS(dh) > hyd->Htol)
@@ -195,7 +200,7 @@ StatType cvstatus(EN_Project *pr, StatType s, double dh, double q)
}
StatType pumpstatus(EN_Project *pr, int k, double dh)
StatusType pumpstatus(Project *pr, int k, double dh)
/*
**--------------------------------------------------
** Input: k = link index
@@ -205,10 +210,11 @@ StatType pumpstatus(EN_Project *pr, int k, double dh)
**--------------------------------------------------
*/
{
Hydraul *hyd = &pr->hydraul;
Network *net = &pr->network;
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);
@@ -231,7 +237,7 @@ StatType pumpstatus(EN_Project *pr, int k, double dh)
}
StatType prvstatus(EN_Project *pr, int k, StatType s, double hset,
StatusType prvstatus(Project *pr, int k, StatusType s, double hset,
double h1, double h2)
/*
**-----------------------------------------------------------
@@ -245,28 +251,31 @@ StatType prvstatus(EN_Project *pr, int k, StatType s, double hset,
**-----------------------------------------------------------
*/
{
StatType status; // Valve's new status
double hml; // Head loss when fully opened
hydraulics_t *hyd = &pr->hydraulics;
Hydraul *hyd = &pr->hydraul;
double htol = hyd->Htol;
Slink *link = &pr->network.Link[k];
StatusType status; // Valve's new status
double hml; // Head loss when fully opened
double htol;
Slink *link;
htol = hyd->Htol;
link = &pr->network.Link[k];
// Head loss when fully open
hml = link->Km * SQR(hyd->LinkFlows[k]);
hml = link->Km * SQR(hyd->LinkFlow[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;
if (hyd->LinkFlow[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;
if (hyd->LinkFlow[k] < -hyd->Qtol) status = CLOSED;
else if (h2 >= hset + htol) status = ACTIVE;
else status = OPEN;
break;
@@ -278,7 +287,7 @@ StatType prvstatus(EN_Project *pr, int k, StatType s, double hset,
break;
case XPRESSURE:
if (hyd->LinkFlows[k] < -hyd->Qtol) status = CLOSED;
if (hyd->LinkFlow[k] < -hyd->Qtol) status = CLOSED;
break;
default:
@@ -288,7 +297,7 @@ StatType prvstatus(EN_Project *pr, int k, StatType s, double hset,
}
StatType psvstatus(EN_Project *pr, int k, StatType s, double hset,
StatusType psvstatus(Project *pr, int k, StatusType s, double hset,
double h1, double h2)
/*
**-----------------------------------------------------------
@@ -302,28 +311,31 @@ StatType psvstatus(EN_Project *pr, int k, StatType s, double hset,
**-----------------------------------------------------------
*/
{
StatType status; // Valve's new status
double hml; // Head loss when fully opened
hydraulics_t *hyd = &pr->hydraulics;
Hydraul *hyd = &pr->hydraul;
double htol = hyd->Htol;
Slink *link = &pr->network.Link[k];
StatusType status; // Valve's new status
double hml; // Head loss when fully opened
double htol;
Slink *link;
htol = hyd->Htol;
link = &pr->network.Link[k];
// Head loss when fully open
hml = link->Km * SQR(hyd->LinkFlows[k]);
hml = link->Km * SQR(hyd->LinkFlow[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;
if (hyd->LinkFlow[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;
if (hyd->LinkFlow[k] < -hyd->Qtol) status = CLOSED;
else if (h1 < hset - htol) status = ACTIVE;
else status = OPEN;
break;
@@ -335,7 +347,7 @@ StatType psvstatus(EN_Project *pr, int k, StatType s, double hset,
break;
case XPRESSURE:
if (hyd->LinkFlows[k] < -hyd->Qtol) status = CLOSED;
if (hyd->LinkFlow[k] < -hyd->Qtol) status = CLOSED;
break;
default:
@@ -345,7 +357,7 @@ StatType psvstatus(EN_Project *pr, int k, StatType s, double hset,
}
StatType fcvstatus(EN_Project *pr, int k, StatType s, double h1, double h2)
StatusType fcvstatus(Project *pr, int k, StatusType s, double h1, double h2)
/*
**-----------------------------------------------------------
** Input: k = link index
@@ -364,19 +376,19 @@ StatType fcvstatus(EN_Project *pr, int k, StatType s, double h1, double h2)
**-----------------------------------------------------------
*/
{
StatType status; // New valve status
hydraulics_t *hyd = &pr->hydraulics;
Hydraul *hyd = &pr->hydraul;
StatusType status; // New valve status
status = s;
if (h1 - h2 < -hyd->Htol)
{
status = XFCV;
}
else if (hyd->LinkFlows[k] < -hyd->Qtol)
else if (hyd->LinkFlow[k] < -hyd->Qtol)
{
status = XFCV;
}
else if (s == XFCV && hyd->LinkFlows[k] >= hyd->LinkSetting[k])
else if (s == XFCV && hyd->LinkFlow[k] >= hyd->LinkSetting[k])
{
status = ACTIVE;
}
@@ -384,7 +396,7 @@ StatType fcvstatus(EN_Project *pr, int k, StatType s, double h1, double h2)
}
void tankstatus(EN_Project *pr, int k, int n1, int n2)
void tankstatus(Project *pr, int k, int n1, int n2)
/*
**----------------------------------------------------------------
** Input: k = link index
@@ -395,19 +407,19 @@ void tankstatus(EN_Project *pr, int k, int n1, int n2)
**----------------------------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
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];
q = hyd->LinkFlow[k];
i = n1 - net->Njuncs;
if (i <= 0)
{
src/inpfile.c
+697 -646
View File
File diff suppressed because it is too large Load Diff
src/input1.c
+591 -612
View File
File diff suppressed because it is too large Load Diff
src/input2.c
+588 -651
View File
File diff suppressed because it is too large Load Diff
src/input3.c
+1570 -1684
View File
File diff suppressed because it is too large Load Diff
src/mempool.c
+14 -7
View File
@@ -1,10 +1,17 @@
/* mempool.c
**
** A simple fast pooled memory allocation package.
**
** Based on code by Steve Hill in Graphics Gems III,
** David Kirk (ed.), Academic Press, Boston, MA, 1992
**
/*
******************************************************************************
Project: OWA EPANET
Version: 2.2
Module: mempool.c
Description: a simple fast poooled memory allocation package
Authors: see AUTHORS
Copyright: see AUTHORS
License: see LICENSE
Last Updated: 11/27/2018
This module is based code by Steve Hill in Graphics Gems III,
David Kirk (ed.), Academic Press, Boston, MA, 1992
******************************************************************************
*/
#include <stdlib.h>
src/mempool.h
+10 -5
View File
@@ -1,9 +1,14 @@
/*
** mempool.h
**
** Header for mempool.c
**
** A simple pooled memory allocator
******************************************************************************
Project: OWA EPANET
Version: 2.2
Module: mempool.h
Description: header for a simple pooled memory allocator
Authors: see AUTHORS
Copyright: see AUTHORS
License: see LICENSE
Last Updated: 11/27/2018
******************************************************************************
*/
#ifndef MEMPOOL_H
src/output.c
+647 -644
View File
File diff suppressed because it is too large Load Diff
src/project.c
+838
View File
@@ -0,0 +1,838 @@
/*
******************************************************************************
Project: OWA EPANET
Version: 2.2
Module: project.c
Description: project data management routines
Authors: see AUTHORS
Copyright: see AUTHORS
License: see LICENSE
Last Updated: 11/27/2018
******************************************************************************
*/
#include <stdio.h>
#include <string.h>
#ifndef __APPLE__
#include <malloc.h>
#else
#include <stdlib.h>
#endif
//*** Need to define WINDOWS to use the getTmpName function ***//
// --- define WINDOWS
#undef WINDOWS
#ifdef _WIN32
#define WINDOWS
#endif
#ifdef __WIN32__
#define WINDOWS
#endif
#ifdef WINDOWS
#include <windows.h>
#endif
#include "types.h"
#include "funcs.h"
int openfiles(Project *pr, const char *f1, const char *f2, const char *f3)
/*----------------------------------------------------------------
** Input: f1 = pointer to name of input file
** f2 = pointer to name of report file
** f3 = pointer to name of binary output file
** Output: none
** Returns: error code
** Purpose: opens input & report files
**----------------------------------------------------------------
*/
{
// Initialize file pointers to NULL
pr->parser.InFile = NULL;
pr->report.RptFile = NULL;
pr->outfile.OutFile = NULL;
pr->outfile.HydFile = NULL;
// Save file names
strncpy(pr->parser.InpFname, f1, MAXFNAME);
strncpy(pr->report.Rpt1Fname, f2, MAXFNAME);
strncpy(pr->outfile.OutFname, f3, MAXFNAME);
if (strlen(f3) > 0) pr->outfile.Outflag = SAVE;
else pr->outfile.Outflag = SCRATCH;
// Check that file names are not identical
if (strcomp(f1, f2) || strcomp(f1, f3) ||
(strcomp(f2, f3) && (strlen(f2) > 0 || strlen(f3) > 0))) return 301;
// Attempt to open input and report files
if (strlen(f1) > 0)
{
if ((pr->parser.InFile = fopen(f1, "rt")) == NULL) return 302;
}
if (strlen(f2) == 0) pr->report.RptFile = stdout;
else if ((pr->report.RptFile = fopen(f2, "wt")) == NULL) return 303;
return 0;
}
int openhydfile(Project *pr)
/*----------------------------------------------------------------
** Input: none
** Output: none
** Returns: error code
** Purpose: opens file that saves hydraulics solution
**----------------------------------------------------------------
*/
{
const int Nnodes = pr->network.Nnodes;
const int Ntanks = pr->network.Ntanks;
const int Nlinks = pr->network.Nlinks;
const int Nvalves = pr->network.Nvalves;
const int Npumps = pr->network.Npumps;
INT4 nsize[6]; // Temporary array
INT4 magic;
INT4 version;
int errcode = 0;
// If HydFile currently open, then close it if its not a scratch file
if (pr->outfile.HydFile != NULL)
{
if (pr->outfile.Hydflag == SCRATCH) return 0;
fclose(pr->outfile.HydFile);
}
// Use Hydflag to determine the type of hydraulics file to use.
// Write error message if the file cannot be opened.
pr->outfile.HydFile = NULL;
switch (pr->outfile.Hydflag)
{
case SCRATCH:
strcpy(pr->outfile.HydFname, pr->TmpHydFname);
pr->outfile.HydFile = fopen(pr->outfile.HydFname, "w+b");
break;
case SAVE:
pr->outfile.HydFile = fopen(pr->outfile.HydFname, "w+b");
break;
case USE:
pr->outfile.HydFile = fopen(pr->outfile.HydFname, "rb");
break;
}
if (pr->outfile.HydFile == NULL) return 305;
// If a previous hydraulics solution is not being used, then
// save the current network size parameters to the file.
if (pr->outfile.Hydflag != USE)
{
magic = MAGICNUMBER;
version = ENGINE_VERSION;
nsize[0] = Nnodes;
nsize[1] = Nlinks;
nsize[2] = Ntanks;
nsize[3] = Npumps;
nsize[4] = Nvalves;
nsize[5] = (int)pr->times.Dur;
fwrite(&magic, sizeof(INT4), 1, pr->outfile.HydFile);
fwrite(&version, sizeof(INT4), 1, pr->outfile.HydFile);
fwrite(nsize, sizeof(INT4), 6, pr->outfile.HydFile);
}
// If a previous hydraulics solution is being used, then
// make sure its network size parameters match those of
// the current network
if (pr->outfile.Hydflag == USE)
{
fread(&magic, sizeof(INT4), 1, pr->outfile.HydFile);
if (magic != MAGICNUMBER) return 306;
fread(&version, sizeof(INT4), 1, pr->outfile.HydFile);
if (version != ENGINE_VERSION) return 306;
if (fread(nsize, sizeof(INT4), 6, pr->outfile.HydFile) < 6) return 306;
if (nsize[0] != Nnodes || nsize[1] != Nlinks || nsize[2] != Ntanks ||
nsize[3] != Npumps || nsize[4] != Nvalves ||
nsize[5] != pr->times.Dur
) return 306;
pr->outfile.SaveHflag = TRUE;
}
// Save current position in hydraulics file
// where storage of hydraulic results begins
pr->outfile.HydOffset = ftell(pr->outfile.HydFile);
return errcode;
}
int openoutfile(Project *pr)
/*----------------------------------------------------------------
** Input: none
** Output: none
** Returns: error code
** Purpose: opens binary output file.
**----------------------------------------------------------------
*/
{
int errcode = 0;
// Close output file if already opened
if (pr->outfile.OutFile != NULL)
{
fclose(pr->outfile.OutFile);
pr->outfile.OutFile = NULL;
}
if (pr->outfile.TmpOutFile != NULL)
{
fclose(pr->outfile.TmpOutFile);
pr->outfile.TmpOutFile = NULL;
}
// If output file name was supplied, then attempt to
// open it. Otherwise open a temporary output file.
if (pr->outfile.Outflag == SAVE)
{
pr->outfile.OutFile = fopen(pr->outfile.OutFname, "w+b");
if (pr->outfile.OutFile == NULL) errcode = 304;
}
else
{
strcpy(pr->outfile.OutFname, pr->TmpOutFname);
pr->outfile.OutFile = fopen(pr->outfile.OutFname, "w+b");
if (pr->outfile.OutFile == NULL) errcode = 304;
}
// Save basic network data & energy usage results
ERRCODE(savenetdata(pr));
pr->outfile.OutOffset1 = ftell(pr->outfile.OutFile);
ERRCODE(saveenergy(pr));
pr->outfile.OutOffset2 = ftell(pr->outfile.OutFile);
// Open temporary file if computing time series statistic
if (!errcode)
{
if (pr->report.Tstatflag != SERIES)
{
pr->outfile.TmpOutFile = fopen(pr->TmpStatFname, "w+b");
if (pr->outfile.TmpOutFile == NULL) errcode = 304;
}
else pr->outfile.TmpOutFile = pr->outfile.OutFile;
}
return errcode;
}
void initpointers(Project *pr)
/*----------------------------------------------------------------
** Input: none
** Output: none
** Purpose: initializes data array pointers to NULL
**----------------------------------------------------------------
*/
{
pr->hydraul.NodeDemand = NULL;
pr->hydraul.NodeHead = NULL;
pr->hydraul.LinkFlow = NULL;
pr->hydraul.LinkStatus = NULL;
pr->hydraul.LinkSetting = NULL;
pr->hydraul.OldStatus = NULL;
pr->hydraul.P = NULL;
pr->hydraul.Y = NULL;
pr->hydraul.Xflow = NULL;
pr->quality.NodeQual = NULL;
pr->quality.PipeRateCoeff = NULL;
pr->network.Node = NULL;
pr->network.Link = NULL;
pr->network.Tank = NULL;
pr->network.Pump = NULL;
pr->network.Valve = NULL;
pr->network.Pattern = NULL;
pr->network.Curve = NULL;
pr->network.Control = NULL;
pr->network.Adjlist = NULL;
pr->network.NodeHashTable = NULL;
pr->network.LinkHashTable = NULL;
pr->parser.Patlist = NULL;
pr->parser.Curvelist = NULL;
pr->hydraul.smatrix.Aii = NULL;
pr->hydraul.smatrix.Aij = NULL;
pr->hydraul.smatrix.F = NULL;
pr->hydraul.smatrix.Order = NULL;
pr->hydraul.smatrix.Row = NULL;
pr->hydraul.smatrix.Ndx = NULL;
pr->hydraul.smatrix.XLNZ = NULL;
pr->hydraul.smatrix.NZSUB = NULL;
pr->hydraul.smatrix.LNZ = NULL;
initrules(pr);
}
int allocdata(Project *pr)
/*----------------------------------------------------------------
** Input: none
** Output: none
** Returns: error code
** Purpose: allocates memory for network data structures
**----------------------------------------------------------------
*/
{
int n;
int errcode = 0;
// Allocate node & link ID hash tables
pr->network.NodeHashTable = hashtable_create();
pr->network.LinkHashTable = hashtable_create();
ERRCODE(MEMCHECK(pr->network.NodeHashTable));
ERRCODE(MEMCHECK(pr->network.LinkHashTable));
// Allocate memory for network nodes
//*************************************************************
// NOTE: Because network components of a given type are indexed
// starting from 1, their arrays must be sized 1
// element larger than the number of components.
//*************************************************************
if (!errcode)
{
n = pr->parser.MaxNodes + 1;
pr->network.Node = (Snode *)calloc(n, sizeof(Snode));
pr->hydraul.NodeDemand = (double *)calloc(n, sizeof(double));
pr->hydraul.NodeHead = (double *)calloc(n, sizeof(double));
pr->quality.NodeQual = (double *)calloc(n, sizeof(double));
ERRCODE(MEMCHECK(pr->network.Node));
ERRCODE(MEMCHECK(pr->hydraul.NodeDemand));
ERRCODE(MEMCHECK(pr->hydraul.NodeHead));
ERRCODE(MEMCHECK(pr->quality.NodeQual));
}
// Allocate memory for network links
if (!errcode)
{
n = pr->parser.MaxLinks + 1;
pr->network.Link = (Slink *)calloc(n, sizeof(Slink));
pr->hydraul.LinkFlow = (double *)calloc(n, sizeof(double));
pr->hydraul.LinkSetting = (double *)calloc(n, sizeof(double));
pr->hydraul.LinkStatus = (StatusType *)calloc(n, sizeof(StatusType));
ERRCODE(MEMCHECK(pr->network.Link));
ERRCODE(MEMCHECK(pr->hydraul.LinkFlow));
ERRCODE(MEMCHECK(pr->hydraul.LinkSetting));
ERRCODE(MEMCHECK(pr->hydraul.LinkStatus));
}
// Allocate memory for tanks, sources, pumps, valves,
// controls, demands, time patterns, & operating curves
if (!errcode)
{
pr->network.Tank =
(Stank *)calloc(pr->parser.MaxTanks + 1, sizeof(Stank));
pr->network.Pump =
(Spump *)calloc(pr->parser.MaxPumps + 1, sizeof(Spump));
pr->network.Valve =
(Svalve *)calloc(pr->parser.MaxValves + 1, sizeof(Svalve));
pr->network.Control =
(Scontrol *)calloc(pr->parser.MaxControls + 1, sizeof(Scontrol));
pr->network.Pattern =
(Spattern *)calloc(pr->parser.MaxPats + 1, sizeof(Spattern));
pr->network.Curve =
(Scurve *)calloc(pr->parser.MaxCurves + 1, sizeof(Scurve));
ERRCODE(MEMCHECK(pr->network.Tank));
ERRCODE(MEMCHECK(pr->network.Pump));
ERRCODE(MEMCHECK(pr->network.Valve));
ERRCODE(MEMCHECK(pr->network.Control));
ERRCODE(MEMCHECK(pr->network.Pattern));
ERRCODE(MEMCHECK(pr->network.Curve));
}
// Initialize pointers used in patterns, curves, and demand category lists
if (!errcode)
{
for (n = 0; n <= pr->parser.MaxPats; n++)
{
pr->network.Pattern[n].Length = 0;
pr->network.Pattern[n].F = NULL;
}
for (n = 0; n <= pr->parser.MaxCurves; n++)
{
pr->network.Curve[n].Npts = 0;
pr->network.Curve[n].Type = G_CURVE;
pr->network.Curve[n].X = NULL;
pr->network.Curve[n].Y = NULL;
}
for (n = 0; n <= pr->parser.MaxNodes; n++)
{
pr->network.Node[n].D = NULL; // node demand
}
}
// Allocate memory for rule base (see RULES.C)
if (!errcode) errcode = allocrules(pr);
return errcode;
}
void freeTmplist(STmplist *t)
/*----------------------------------------------------------------
** Input: t = pointer to start of a temporary list
** Output: none
** Purpose: frees memory used for temporary storage
** of pattern & curve data
**----------------------------------------------------------------
*/
{
STmplist *tnext;
while (t != NULL)
{
tnext = t->next;
freeFloatlist(t->x);
freeFloatlist(t->y);
free(t);
t = tnext;
}
}
void freeFloatlist(SFloatlist *f)
/*----------------------------------------------------------------
** Input: f = pointer to start of list of floats
** Output: none
** Purpose: frees memory used for storing list of floats
**----------------------------------------------------------------
*/
{
SFloatlist *fnext;
while (f != NULL)
{
fnext = f->next;
free(f);
f = fnext;
}
}
void freedata(Project *pr)
/*----------------------------------------------------------------
** Input: none
** Output: none
** Purpose: frees memory allocated for network data structures.
**----------------------------------------------------------------
*/
{
int j;
Pdemand demand, nextdemand;
Psource source;
// Free memory for computed results
free(pr->hydraul.NodeDemand);
free(pr->hydraul.NodeHead);
free(pr->hydraul.LinkFlow);
free(pr->hydraul.LinkSetting);
free(pr->hydraul.LinkStatus);
free(pr->quality.NodeQual);
// Free memory used for nodal adjacency lists
freeadjlists(&pr->network);
// Free memory for node data
if (pr->network.Node != NULL)
{
for (j = 0; j <= pr->parser.MaxNodes; j++)
{
// Free memory used for demand category list
demand = pr->network.Node[j].D;
while (demand != NULL)
{
nextdemand = demand->next;
free(demand);
demand = nextdemand;
}
// Free memory used for WQ source data
source = pr->network.Node[j].S;
if (source != NULL) free(source);
}
free(pr->network.Node);
}
// Free memory for other network objects
free(pr->network.Link);
free(pr->network.Tank);
free(pr->network.Pump);
free(pr->network.Valve);
free(pr->network.Control);
// Free memory for time patterns
if (pr->network.Pattern != NULL)
{
for (j = 0; j <= pr->parser.MaxPats; j++) free(pr->network.Pattern[j].F);
free(pr->network.Pattern);
}
// Free memory for curves
if (pr->network.Curve != NULL)
{
for (j = 0; j <= pr->parser.MaxCurves; j++)
{
free(pr->network.Curve[j].X);
free(pr->network.Curve[j].Y);
}
free(pr->network.Curve);
}
// Free memory for rule base (see RULES.C)
freerules(pr);
// Free hash table memory
if (pr->network.NodeHashTable != NULL)
{
hashtable_free(pr->network.NodeHashTable);
}
if (pr->network.LinkHashTable != NULL)
{
hashtable_free(pr->network.LinkHashTable);
}
}
int buildadjlists(Network *net)
/*
**--------------------------------------------------------------
** Input: none
** Output: returns error code
** Purpose: builds linked list of links adjacent to each node
**--------------------------------------------------------------
*/
{
int i, j, k;
int errcode = 0;
Padjlist alink;
// Create an array of adjacency lists
freeadjlists(net);
net->Adjlist = (Padjlist *)calloc(net->Nnodes + 1, sizeof(Padjlist));
if (net->Adjlist == NULL) return 101;
for (i = 0; i <= net->Nnodes; i++) net->Adjlist[i] = NULL;
// For each link, update adjacency lists of its end nodes
for (k = 1; k <= net->Nlinks; k++)
{
i = net->Link[k].N1;
j = net->Link[k].N2;
// Include link in start node i's list
alink = (struct Sadjlist *) malloc(sizeof(struct Sadjlist));
if (alink == NULL)
{
errcode = 101;
break;
}
alink->node = j;
alink->link = k;
alink->next = net->Adjlist[i];
net->Adjlist[i] = alink;
// Include link in end node j's list
alink = (struct Sadjlist *) malloc(sizeof(struct Sadjlist));
if (alink == NULL)
{
errcode = 101;
break;
}
alink->node = i;
alink->link = k;
alink->next = net->Adjlist[j];
net->Adjlist[j] = alink;
}
if (errcode) freeadjlists(net);
return errcode;
}
void freeadjlists(Network *net)
/*
**--------------------------------------------------------------
** Input: none
** Output: none
** Purpose: frees memory used for nodal adjacency lists
**--------------------------------------------------------------
*/
{
int i;
Padjlist alink;
if (net->Adjlist == NULL) return;
for (i = 0; i <= net->Nnodes; i++)
{
for (alink = net->Adjlist[i]; alink != NULL; alink = net->Adjlist[i])
{
net->Adjlist[i] = alink->next;
free(alink);
}
}
FREE(net->Adjlist);
}
int incontrols(Project *pr, int objType, int index)
/*----------------------------------------------------------------
** Input: objType = type of object (either NODE or LINK)
** index = the object's index
** Output: none
** Returns: 1 if any controls contain the object; 0 if not
** Purpose: determines if any simple or rule-based controls
** contain a particular node or link.
**----------------------------------------------------------------
*/
{
Network *net = &pr->network;
Rules *rules = &pr->rules;
int i, ruleObject;
Spremise *premise;
Saction *action;
// Check simple controls
for (i = 1; i <= net->Ncontrols; i++)
{
if (objType == NODE && net->Control[i].Node == index) return 1;
if (objType == LINK && net->Control[i].Link == index) return 1;
}
// Check rule-based controls
for (i = 1; i <= net->Nrules; i++)
{
// Convert objType to a rule object type
if (objType == NODE) ruleObject = 6;
else ruleObject = 7;
// Check rule's premises
premise = net->Rule[i].Premises;
while (premise != NULL)
{
if (ruleObject == premise->object && premise->index == index) return 1;
premise = premise->next;
}
// Rule actions only need to be checked for link objects
if (objType == LINK)
{
// Check rule's THEN actions
action = net->Rule[i].ThenActions;
while (action != NULL)
{
if (action->link == index) return 1;
action = action->next;
}
// Check rule's ELSE actions
action = net->Rule[i].ElseActions;
while (action != NULL)
{
if (action->link == index) return 1;
action = action->next;
}
}
}
return 0;
}
int findnode(Network *network, char *id)
/*----------------------------------------------------------------
** Input: id = node ID
** Output: none
** Returns: index of node with given ID, or 0 if ID not found
** Purpose: uses hash table to find index of node with given ID
**----------------------------------------------------------------
*/
{
return (hashtable_find(network->NodeHashTable, id));
}
int findlink(Network *network, char *id)
/*----------------------------------------------------------------
** Input: id = link ID
** Output: none
** Returns: index of link with given ID, or 0 if ID not found
** Purpose: uses hash table to find index of link with given ID
**----------------------------------------------------------------
*/
{
return (hashtable_find(network->LinkHashTable, id));
}
int findtank(Network *network, int index)
/*----------------------------------------------------------------
** Input: index = node index
** Output: none
** Returns: index of tank with given node id, or NOTFOUND if tank not found
** Purpose: for use in the deletenode function
**----------------------------------------------------------------
*/
{
int i;
for (i = 1; i <= network->Ntanks; i++)
{
if (network->Tank[i].Node == index) return i;
}
return NOTFOUND;
}
int findpump(Network *network, int index)
/*----------------------------------------------------------------
** Input: index = link ID
** Output: none
** Returns: index of pump with given link id, or NOTFOUND if pump not found
** Purpose: for use in the deletelink function
**----------------------------------------------------------------
*/
{
int i;
for (i = 1; i <= network->Npumps; i++)
{
if (network->Pump[i].Link == index) return i;
}
return NOTFOUND;
}
int findvalve(Network *network, int index)
/*----------------------------------------------------------------
** Input: index = link ID
** Output: none
** Returns: index of valve with given link id, or NOTFOUND if valve not found
** Purpose: for use in the deletelink function
**----------------------------------------------------------------
*/
{
int i;
for (i = 1; i <= network->Nvalves; i++)
{
if (network->Valve[i].Link == index) return i;
}
return NOTFOUND;
}
char *getTmpName(char *fname)
//----------------------------------------------------------------
// Input: fname = file name string
// Output: returns pointer to file name
// Purpose: creates a temporary file name with path prepended to it.
//----------------------------------------------------------------
{
#ifdef _WIN32
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;
// --- 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);
// --- for non-Windows systems:
#else
// --- use system function mkstemp() to create a temporary file name
strcpy(fname, "enXXXXXX");
mkstemp(fname);
#endif
return fname;
}
int strcomp(const char *s1, const char *s2)
/*---------------------------------------------------------------
** Input: s1 = character string
** s2 = character string
** Output: none
** Returns: 1 if s1 is same as s2, 0 otherwise
** Purpose: case insensitive comparison of strings s1 & s2
**---------------------------------------------------------------
*/
{
int i;
for (i = 0; UCHAR(s1[i]) == UCHAR(s2[i]); i++)
{
if (!s1[i + 1] && !s2[i + 1]) return 1;
}
return 0;
}
double interp(int n, double x[], double y[], double xx)
/*----------------------------------------------------------------
** Input: n = number of data pairs defining a curve
** x = x-data values of curve
** y = y-data values of curve
** xx = specified x-value
** Output: none
** Returns: y-value on curve at x = xx
** Purpose: uses linear interpolation to find y-value on a
** data curve corresponding to specified x-value.
** NOTE: does not extrapolate beyond endpoints of curve.
**----------------------------------------------------------------
*/
{
int k, m;
double dx, dy;
m = n - 1; // Highest data index
if (xx <= x[0]) return (y[0]); // xx off low end of curve
for (k = 1; k <= m; k++) // Bracket xx on curve
{
if (x[k] >= xx) // Interp. over interval
{
dx = x[k] - x[k - 1];
dy = y[k] - y[k - 1];
if (ABS(dx) < TINY) return (y[k]);
else return (y[k] - (x[k] - xx) * dy / dx);
}
}
return (y[m]); // xx off high end of curve
}
char *geterrmsg(int errcode, char *msg)
/*----------------------------------------------------------------
** Input: errcode = error code
** Output: none
** Returns: pointer to string with error message
** Purpose: retrieves text of error message
**----------------------------------------------------------------
*/
{
switch (errcode) { /* Warnings */
#define DAT(code,enumer,string) case code: strcpy(msg, string); break;
//#define DAT(code,enumer,string) case code: sprintf(msg, "Error %d: %s", code, string); break;
#include "errors.dat"
#undef DAT
default:
strcpy(msg, "");
}
return (msg);
}
void errmsg(Project *pr, int errcode)
/*----------------------------------------------------------------
** Input: errcode = error code
** Output: none
** Purpose: writes error message to report file
**----------------------------------------------------------------
*/
{
if (errcode == 309) /* Report file write error - */
{ /* Do not write msg to file. */
}
else if (pr->report.RptFile != NULL && pr->report.Messageflag)
{
writeline(pr, geterrmsg(errcode, pr->Msg));
}
}
void writewin(void(*vp)(char *), char *s)
/*----------------------------------------------------------------
** Input: text string
** Output: none
** Purpose: passes character string to viewprog() in
** application which calls the EPANET DLL
**----------------------------------------------------------------
*/
{
char progmsg[MAXMSG + 1];
if (vp != NULL)
{
strncpy(progmsg, s, MAXMSG);
vp(progmsg);
}
}
src/quality.c
+128 -136
View File
@@ -1,13 +1,14 @@
/*
*********************************************************************
QUALITY.C -- water quality engine module for the EPANET program
This module works together with QUALROUTE.C and QUALREACT.C to
compute transport and reaction of a water quality constituent within
a pipe network over a series of time steps.
*********************************************************************
******************************************************************************
Project: OWA EPANET
Version: 2.2
Module: quality.c
Description: implements EPANET's water quality engine
Authors: see AUTHORS
Copyright: see AUTHORS
License: see LICENSE
Last Updated: 11/27/2018
******************************************************************************
*/
#include <stdio.h>
@@ -24,36 +25,35 @@ a pipe network over a series of time steps.
#include "funcs.h"
// Stagnant flow tolerance
const double QZERO = 0.005 / GPMperCFS; // 0.005 gpm = 1.114e-5 cfs
const double Q_STAGNANT = 0.005 / GPMperCFS; // 0.005 gpm = 1.114e-5 cfs
// Exported Functions (declared in FUNCS.H)
//int openqual(EN_Project *pr);
//void initqual(EN_Project *pr);
//int runqual(EN_Project *pr, long *);
//int nextqual(EN_Project *pr, long *);
//int stepqual(EN_Project *pr, long *);
//int closequal(EN_Project *pr);
//double avgqual(EN_Project *pr, int);
double findsourcequal(EN_Project *pr, int, double, double, long);
// Exported functions (declared in FUNCS.H)
//int openqual(Project *);
//void initqual(Project *);
//int runqual(Project *, long *);
//int nextqual(Project *, long *);
//int stepqual(Project *, long *);
//int closequal(Project *);
//double avgqual(Project *, int);
double findsourcequal(Project *, int, double, double, long);
// Imported Functions
extern char setreactflag(EN_Project *pr);
// Imported functions
extern char setreactflag(Project *);
extern double getucf(double);
extern void ratecoeffs(EN_Project *pr);
extern int buildilists(EN_Project *pr);
extern void initsegs(EN_Project *pr);
extern void reversesegs(EN_Project *pr, int);
extern int sortnodes(EN_Project *pr);
extern void transport(EN_Project *pr, long);
extern void ratecoeffs(Project *);
extern void initsegs(Project *);
extern void reversesegs(Project *, int);
extern int sortnodes(Project *);
extern void transport(Project *, long);
// Local Functions
static double sourcequal(EN_Project *pr, Psource);
static void evalmassbalance(EN_Project *pr);
static double findstoredmass(EN_Project *pr);
static int flowdirchanged(EN_Project *pr);
// Local functions
static double sourcequal(Project *, Psource);
static void evalmassbalance(Project *);
static double findstoredmass(Project *);
static int flowdirchanged(Project *);
int openqual(EN_Project *pr)
int openqual(Project *pr)
/*
**--------------------------------------------------------------
** Input: none
@@ -62,16 +62,23 @@ int openqual(EN_Project *pr)
**--------------------------------------------------------------
*/
{
Network *net = &pr->network;
Quality *qual = &pr->quality;
int errcode = 0;
int n;
quality_t *qual = &pr->quality;
EN_Network *net = &pr->network;
// Build nodal adjacency lists if they don't already exist
if (net->Adjlist == NULL)
{
errcode = buildadjlists(net);
if (errcode ) return errcode;
}
// Create a memory pool for water quality segments
qual->OutOfMemory = FALSE;
qual->SegPool = mempool_create();
if (qual->SegPool == NULL) errcode = 101;;
if (qual->SegPool == NULL) errcode = 101;
// Allocate arrays for link flow direction & reaction rates
n = net->Nlinks + 1;
@@ -83,34 +90,19 @@ int openqual(EN_Project *pr)
qual->FirstSeg = (Pseg *)calloc(n, sizeof(Pseg));
qual->LastSeg = (Pseg *)calloc(n, sizeof(Pseg));
// Allocate memory for sorted nodes and link incidence lists
// ... Ilist contains the list of link indexes that are incident
// on each node in compact format
n = 2 * net->Nlinks + 2;
qual->Ilist = (int *)calloc(n, sizeof(int));
// ... IlistPtr holds the position in Ilist where the
// incidence list for each node begins
n = net->Nnodes + 1;
qual->IlistPtr = (int *)calloc(n + 2, sizeof(int));
// ... SortedNodes contains the list of node indexes in topological
// order
// Allocate memory for topologically sorted nodes
qual->SortedNodes = (int *)calloc(n, sizeof(int));
ERRCODE(MEMCHECK(qual->FlowDir));
ERRCODE(MEMCHECK(qual->PipeRateCoeff));
ERRCODE(MEMCHECK(qual->FirstSeg));
ERRCODE(MEMCHECK(qual->LastSeg));
ERRCODE(MEMCHECK(qual->Ilist));
ERRCODE(MEMCHECK(qual->IlistPtr));
ERRCODE(MEMCHECK(qual->SortedNodes));
// Build link incidence lists
if (!errcode) errcode = buildilists(pr);
return errcode;
}
int initqual(EN_Project *pr)
int initqual(Project *pr)
/*
**--------------------------------------------------------------
** Input: none
@@ -119,22 +111,24 @@ int initqual(EN_Project *pr)
**--------------------------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Quality *qual = &pr->quality;
Times *time = &pr->times;
int i;
int errcode = 0;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
// Re-position hydraulics file
if (!hyd->OpenHflag)
{
fseek(pr->out_files.HydFile, pr->out_files.HydOffset, SEEK_SET);
fseek(pr->outfile.HydFile, pr->outfile.HydOffset, SEEK_SET);
}
// Set elapsed times to zero
qual->Qtime = 0;
pr->time_options.Htime = 0;
pr->time_options.Rtime = pr->time_options.Rstart;
time->Qtime = 0;
time->Htime = 0;
time->Rtime = time->Rstart;
pr->report.Nperiods = 0;
// Initialize node quality
@@ -183,17 +177,17 @@ int initqual(EN_Project *pr)
qual->Wsource = 0.0;
// Initialize mass balance components
qual->massbalance.initial = findstoredmass(pr);
qual->massbalance.inflow = 0.0;
qual->massbalance.outflow = 0.0;
qual->massbalance.reacted = 0.0;
qual->massbalance.final = 0.0;
qual->massbalance.ratio = 0.0;
qual->MassBalance.initial = findstoredmass(pr);
qual->MassBalance.inflow = 0.0;
qual->MassBalance.outflow = 0.0;
qual->MassBalance.reacted = 0.0;
qual->MassBalance.final = 0.0;
qual->MassBalance.ratio = 0.0;
return errcode;
}
int runqual(EN_Project *pr, long *t)
int runqual(Project *pr, long *t)
/*
**--------------------------------------------------------------
** Input: none
@@ -204,19 +198,19 @@ int runqual(EN_Project *pr, long *t)
**--------------------------------------------------------------
*/
{
Hydraul *hyd = &pr->hydraul;
Quality *qual = &pr->quality;
Times *time = &pr->times;
long hydtime; // Hydraulic solution time
long hydstep; // Hydraulic time step
int errcode = 0;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
time_options_t *time = &pr->time_options;
// Update reported simulation time
*t = qual->Qtime;
*t = time->Qtime;
// Read hydraulic solution from hydraulics file
if (qual->Qtime == time->Htime)
if (time->Qtime == time->Htime)
{
// Read hydraulic results from file
if (!hyd->OpenHflag)
@@ -229,7 +223,7 @@ int runqual(EN_Project *pr, long *t)
// Save current results to output file
if (time->Htime >= time->Rtime)
{
if (pr->save_options.Saveflag)
if (pr->outfile.Saveflag)
{
errcode = saveoutput(pr);
pr->report.Nperiods++;
@@ -239,7 +233,7 @@ int runqual(EN_Project *pr, long *t)
if (errcode) return errcode;
// If simulating water quality
if (qual->Qualflag != NONE && qual->Qtime < time->Dur)
if (qual->Qualflag != NONE && time->Qtime < time->Dur)
{
// ... compute reaction rate coeffs.
if (qual->Reactflag && qual->Qualflag != AGE) ratecoeffs(pr);
@@ -256,7 +250,7 @@ int runqual(EN_Project *pr, long *t)
}
int nextqual(EN_Project *pr, long *tstep)
int nextqual(Project *pr, long *tstep)
/*
**--------------------------------------------------------------
** Input: none
@@ -267,17 +261,17 @@ int nextqual(EN_Project *pr, long *tstep)
**--------------------------------------------------------------
*/
{
Quality *qual = &pr->quality;
Times *time = &pr->times;
long hydstep; // Time step until next hydraulic event
long dt, qtime;
int errcode = 0;
quality_t *qual = &pr->quality;
time_options_t *time = &pr->time_options;
// Find time step till next hydraulic event
*tstep = 0;
hydstep = 0;
if (time->Htime <= time->Dur) hydstep = time->Htime - qual->Qtime;
if (time->Htime <= time->Dur) hydstep = time->Htime - time->Qtime;
// Perform water quality routing over this time step
if (qual->Qualflag != NONE && hydstep > 0)
@@ -286,7 +280,7 @@ int nextqual(EN_Project *pr, long *tstep)
qtime = 0;
while (!qual->OutOfMemory && qtime < hydstep)
{
dt = MIN(qual->Qstep, hydstep - qtime);
dt = MIN(time->Qstep, hydstep - qtime);
qtime += dt;
transport(pr, dt);
}
@@ -298,7 +292,7 @@ int nextqual(EN_Project *pr, long *tstep)
// Update current time
if (!errcode) *tstep = hydstep;
qual->Qtime += hydstep;
time->Qtime += hydstep;
// If no more time steps remain
if (!errcode && *tstep == 0)
@@ -310,13 +304,13 @@ int nextqual(EN_Project *pr, long *tstep)
}
// ... write the final portion of the binary output file
if (pr->save_options.Saveflag) errcode = savefinaloutput(pr);
if (pr->outfile.Saveflag) errcode = savefinaloutput(pr);
}
return errcode;
}
int stepqual(EN_Project *pr, long *tleft)
int stepqual(Project *pr, long *tleft)
/*
**--------------------------------------------------------------
** Input: none
@@ -327,19 +321,20 @@ int stepqual(EN_Project *pr, long *tleft)
**--------------------------------------------------------------
*/
{
Quality *qual = &pr->quality;
Times *time = &pr->times;
long dt, hstep, t, tstep;
int errcode = 0;
quality_t *qual = &pr->quality;
tstep = qual->Qstep;
tstep = time->Qstep;
do
{
// Set local time step to quality time step
dt = tstep;
// Find time step until next hydraulic event
hstep = pr->time_options.Htime - qual->Qtime;
hstep = time->Htime - time->Qtime;
// If next hydraulic event occurs before end of local time step
if (hstep < dt)
@@ -349,21 +344,21 @@ int stepqual(EN_Project *pr, long *tleft)
// ... transport quality over local time step
if (qual->Qualflag != NONE) transport(pr, dt);
qual->Qtime += dt;
time->Qtime += dt;
// ... quit if running quality concurrently with hydraulics
if (pr->hydraulics.OpenHflag) break;
if (pr->hydraul.OpenHflag) break;
// ... otherwise call runqual() to update hydraulics
errcode = runqual(pr, &t);
qual->Qtime = t;
time->Qtime = t;
}
// Otherwise transport quality over current local time step
else
{
if (qual->Qualflag != NONE) transport(pr, dt);
qual->Qtime += dt;
time->Qtime += dt;
}
// Reduce quality time step by local time step
@@ -376,7 +371,7 @@ int stepqual(EN_Project *pr, long *tleft)
evalmassbalance(pr);
// Update total simulation time left
*tleft = pr->time_options.Dur - qual->Qtime;
*tleft = time->Dur - time->Qtime;
// If no more time steps remain
if (!errcode && *tleft == 0)
@@ -388,13 +383,13 @@ int stepqual(EN_Project *pr, long *tleft)
}
// ... write the final portion of the binary output file
if (pr->save_options.Saveflag) errcode = savefinaloutput(pr);
if (pr->outfile.Saveflag) errcode = savefinaloutput(pr);
}
return errcode;
}
int closequal(EN_Project *pr)
int closequal(Project *pr)
/*
**--------------------------------------------------------------
** Input: none
@@ -403,7 +398,7 @@ int closequal(EN_Project *pr)
**--------------------------------------------------------------
*/
{
quality_t *qual = &pr->quality;
Quality *qual = &pr->quality;
int errcode = 0;
if (qual->SegPool) mempool_delete(qual->SegPool);
@@ -411,14 +406,12 @@ int closequal(EN_Project *pr)
FREE(qual->LastSeg);
FREE(qual->PipeRateCoeff);
FREE(qual->FlowDir);
FREE(qual->Ilist);
FREE(qual->IlistPtr);
FREE(qual->SortedNodes);
return errcode;
}
double avgqual(EN_Project *pr, int k)
double avgqual(Project *pr, int k)
/*
**--------------------------------------------------------------
** Input: k = link index
@@ -427,12 +420,12 @@ double avgqual(EN_Project *pr, int k)
**--------------------------------------------------------------
*/
{
Network *net = &pr->network;
Quality *qual = &pr->quality;
double vsum = 0.0, msum = 0.0;
Pseg seg;
EN_Network *net = &pr->network;
quality_t *qual = &pr->quality;
if (qual->Qualflag == NONE) return 0.0;
// Sum up the quality and volume in each segment of the link
@@ -459,7 +452,7 @@ double avgqual(EN_Project *pr, int k)
}
double findsourcequal(EN_Project *pr, int n, double volin, double volout, long tstep)
double findsourcequal(Project *pr, int n, double volin, double volout, long tstep)
/*
**---------------------------------------------------------------------
** Input: n = node index
@@ -472,14 +465,14 @@ double findsourcequal(EN_Project *pr, int n, double volin, double volout, long t
**---------------------------------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Quality *qual = &pr->quality;
Times *time = &pr->times;
double massadded = 0.0, c;
Psource source;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
time_options_t *time = &pr->time_options;
// Sources only apply to CHEMICAL analyses
if (qual->Qualflag != CHEM) return 0.0;
@@ -487,7 +480,7 @@ double findsourcequal(EN_Project *pr, int n, double volin, double volout, long t
source = net->Node[n].S;
if (source == NULL) return 0.0;
if (source->C0 == 0.0) return 0.0;
if (volout / tstep <= QZERO) return 0.0;
if (volout / tstep <= Q_STAGNANT) return 0.0;
// Added source concentration depends on source type
c = sourcequal(pr, source);
@@ -540,7 +533,7 @@ double findsourcequal(EN_Project *pr, int n, double volin, double volout, long t
}
double sourcequal(EN_Project *pr, Psource source)
double sourcequal(Project *pr, Psource source)
/*
**--------------------------------------------------------------
** Input: source = a water quality source object
@@ -549,14 +542,14 @@ double sourcequal(EN_Project *pr, Psource source)
**--------------------------------------------------------------
*/
{
Network *net = &pr->network;
Quality *qual = &pr->quality;
Times *time = &pr->times;
int i;
long k;
double c;
EN_Network *net = &pr->network;
quality_t *qual = &pr->quality;
time_options_t *time = &pr->time_options;
// Get source concentration (or mass flow) in original units
c = source->C0;
@@ -568,13 +561,13 @@ double sourcequal(EN_Project *pr, Psource source)
// Apply time pattern if assigned
i = source->Pat;
if (i == 0) return c;
k = ((qual->Qtime + time->Pstart) / time->Pstep) %
k = ((time->Qtime + time->Pstart) / time->Pstep) %
(long)net->Pattern[i].Length;
return (c * net->Pattern[i].F[k]);
}
void evalmassbalance(EN_Project *pr)
void evalmassbalance(Project *pr)
/*
**--------------------------------------------------------------
** Input: none
@@ -584,27 +577,28 @@ void evalmassbalance(EN_Project *pr)
**--------------------------------------------------------------
*/
{
Quality *qual = &pr->quality;
double massin;
double massout;
double massreacted;
quality_t *qual = &pr->quality;
if (qual->Qualflag == NONE) qual->massbalance.ratio = 1.0;
if (qual->Qualflag == NONE) qual->MassBalance.ratio = 1.0;
else
{
qual->massbalance.final = findstoredmass(pr);
massin = qual->massbalance.initial + qual->massbalance.inflow;
massout = qual->massbalance.outflow + qual->massbalance.final;
massreacted = qual->massbalance.reacted;
qual->MassBalance.final = findstoredmass(pr);
massin = qual->MassBalance.initial + qual->MassBalance.inflow;
massout = qual->MassBalance.outflow + qual->MassBalance.final;
massreacted = qual->MassBalance.reacted;
if (massreacted > 0.0) massout += massreacted;
else massin -= massreacted;
if (massin == 0.0) qual->massbalance.ratio = 1.0;
else qual->massbalance.ratio = massout / massin;
if (massin == 0.0) qual->MassBalance.ratio = 1.0;
else qual->MassBalance.ratio = massout / massin;
}
}
double findstoredmass(EN_Project *pr)
double findstoredmass(Project *pr)
/*
**--------------------------------------------------------------
** Input: none
@@ -614,13 +608,13 @@ double findstoredmass(EN_Project *pr)
**--------------------------------------------------------------
*/
{
Network *net = &pr->network;
Quality *qual = &pr->quality;
int i, k;
double totalmass = 0.0;
Pseg seg;
EN_Network *net = &pr->network;
quality_t *qual = &pr->quality;
// Mass residing in each pipe
for (k = 1; k <= net->Nlinks; k++)
{
@@ -654,7 +648,7 @@ double findstoredmass(EN_Project *pr)
return totalmass;
}
int flowdirchanged(EN_Project *pr)
int flowdirchanged(Project *pr)
/*
**--------------------------------------------------------------
** Input: none
@@ -663,25 +657,25 @@ int flowdirchanged(EN_Project *pr)
**--------------------------------------------------------------
*/
{
Hydraul *hyd = &pr->hydraul;
Quality *qual = &pr->quality;
int k;
int result = FALSE;
int newdir;
int olddir;
double q;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
// Examine each network link
for (k = 1; k <= pr->network.Nlinks; k++)
{
// Determine sign (+1 or -1) of new flow rate
olddir = qual->FlowDir[k];
q = (hyd->LinkStatus[k] <= CLOSED) ? 0.0 : hyd->LinkFlows[k];
q = (hyd->LinkStatus[k] <= CLOSED) ? 0.0 : hyd->LinkFlow[k];
newdir = SGN(q);
// Indicate if flow is negligible
if (fabs(q) < QZERO) newdir = 0;
if (fabs(q) < Q_STAGNANT) newdir = 0;
// Reverse link's volume segments if flow direction changes sign
if (newdir * olddir < 0) reversesegs(pr, k);
@@ -696,5 +690,3 @@ int flowdirchanged(EN_Project *pr)
}
return result;
}
/************************* End of QUALITY.C ***************************/
src/qualreact.c
+93 -85
View File
@@ -1,40 +1,45 @@
/*
*********************************************************************
QUALREACT.C -- water quality reaction module for the EPANET program
*********************************************************************
******************************************************************************
Project: OWA EPANET
Version: 2.2
Module: qualreact.c
Description: computes water quality reactions within pipes and tanks
Authors: see AUTHORS
Copyright: see AUTHORS
License: see LICENSE
Last Updated: 11/27/2018
******************************************************************************
*/
#include <stdio.h>
#include <math.h>
#include "types.h"
// Exported Functions
char setreactflag(EN_Project *pr);
// Exported functions
char setreactflag(Project *);
double getucf(double);
void ratecoeffs(EN_Project *pr);
void reactpipes(EN_Project *pr, long);
void reacttanks(EN_Project *pr, long);
double mixtank(EN_Project *pr, int, double, double ,double);
void ratecoeffs(Project *);
void reactpipes(Project *, long);
void reacttanks(Project *, long);
double mixtank(Project *, int, double, double ,double);
// Imported Functions
extern void addseg(EN_Project *pr, int, double, double);
// Imported functions
extern void addseg(Project *, int, double, double);
// Local Functions
static double piperate(EN_Project *pr, int);
static double pipereact(EN_Project *pr, int, double, double, long);
static double tankreact(EN_Project *pr, double, double, double, long);
static double bulkrate(EN_Project *pr, double, double, double);
static double wallrate(EN_Project *pr, double, double, double, double);
// Local functions
static double piperate(Project *, int);
static double pipereact(Project *, int, double, double, long);
static double tankreact(Project *, double, double, double, long);
static double bulkrate(Project *, double, double, double);
static double wallrate(Project *, double, double, double, double);
static void tankmix1(EN_Project *pr, int, double, double, double);
static void tankmix2(EN_Project *pr, int, double, double, double);
static void tankmix3(EN_Project *pr, int, double, double, double);
static void tankmix4(EN_Project *pr, int, double, double, double);
static void tankmix1(Project *, int, double, double, double);
static void tankmix2(Project *, int, double, double, double);
static void tankmix3(Project *, int, double, double, double);
static void tankmix4(Project *, int, double, double, double);
char setreactflag(EN_Project *pr)
char setreactflag(Project *pr)
/*
**-----------------------------------------------------------
** Input: none
@@ -43,8 +48,8 @@ char setreactflag(EN_Project *pr)
**-----------------------------------------------------------
*/
{
Network *net = &pr->network;
int i;
EN_Network *net = &pr->network;
if (pr->quality.Qualflag == TRACE) return 0;
else if (pr->quality.Qualflag == AGE) return 1;
@@ -82,7 +87,7 @@ double getucf(double order)
}
void ratecoeffs(EN_Project *pr)
void ratecoeffs(Project *pr)
/*
**--------------------------------------------------------------
** Input: none
@@ -91,12 +96,12 @@ void ratecoeffs(EN_Project *pr)
**--------------------------------------------------------------
*/
{
Network *net = &pr->network;
Quality *qual = &pr->quality;
int k;
double kw;
EN_Network *net = &pr->network;
quality_t *qual = &pr->quality;
for (k = 1; k <= net->Nlinks; k++)
{
kw = net->Link[k].Kw;
@@ -107,7 +112,7 @@ void ratecoeffs(EN_Project *pr)
}
void reactpipes(EN_Project *pr, long dt)
void reactpipes(Project *pr, long dt)
/*
**--------------------------------------------------------------
** Input: dt = time step
@@ -116,13 +121,13 @@ void reactpipes(EN_Project *pr, long dt)
**--------------------------------------------------------------
*/
{
Network *net = &pr->network;
Quality *qual = &pr->quality;
int k;
Pseg seg;
double cseg, rsum, vsum;
EN_Network *net = &pr->network;
quality_t *qual = &pr->quality;
// Examine each link in network
for (k = 1; k <= net->Nlinks; k++)
{
@@ -140,7 +145,7 @@ void reactpipes(EN_Project *pr, long dt)
seg->c = pipereact(pr, k, seg->c, seg->v, dt);
// Update reaction component of mass balance
qual->massbalance.reacted += (cseg - seg->c) * seg->v;
qual->MassBalance.reacted += (cseg - seg->c) * seg->v;
// Accumulate volume-weighted reaction rate
if (qual->Qualflag == CHEM)
@@ -158,7 +163,7 @@ void reactpipes(EN_Project *pr, long dt)
}
void reacttanks(EN_Project *pr, long dt)
void reacttanks(Project *pr, long dt)
/*
**--------------------------------------------------------------
** Input: dt = time step
@@ -167,13 +172,13 @@ void reacttanks(EN_Project *pr, long dt)
**--------------------------------------------------------------
*/
{
Network *net = &pr->network;
Quality *qual = &pr->quality;
int i, k;
double c;
Pseg seg;
EN_Network *net = &pr->network;
quality_t *qual = &pr->quality;
Stank *tank;
Stank *tank;
// Examine each tank in network
for (i = 1; i <= net->Ntanks; i++)
@@ -191,14 +196,14 @@ void reacttanks(EN_Project *pr, long dt)
{
c = seg->c;
seg->c = tankreact(pr, seg->c, seg->v, tank->Kb, dt);
qual->massbalance.reacted += (c - seg->c) * seg->v;
qual->MassBalance.reacted += (c - seg->c) * seg->v;
seg = seg->prev;
}
}
}
double piperate(EN_Project *pr, int k)
double piperate(Project *pr, int k)
/*
**--------------------------------------------------------------
** Input: k = link index
@@ -209,11 +214,11 @@ double piperate(EN_Project *pr, int k)
**--------------------------------------------------------------
*/
{
double a, d, u, q, kf, kw, y, Re, Sh;
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Quality *qual = &pr->quality;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
double a, d, u, q, kf, kw, y, Re, Sh;
d = net->Link[k].Diam; // Pipe diameter, ft
@@ -226,7 +231,7 @@ double piperate(EN_Project *pr, int k)
// Compute Reynolds No.
// Flow rate made consistent with how its saved to hydraulics file
q = (hyd->LinkStatus[k] <= CLOSED) ? 0.0 : hyd->LinkFlows[k];
q = (hyd->LinkStatus[k] <= CLOSED) ? 0.0 : hyd->LinkFlow[k];
a = PI * d * d / 4.0; // pipe area
u = fabs(q) / a; // flow velocity
Re = u * d / hyd->Viscos; // Reynolds number
@@ -258,7 +263,7 @@ double piperate(EN_Project *pr, int k)
}
double pipereact(EN_Project *pr, int k, double c, double v, long dt)
double pipereact(Project *pr, int k, double c, double v, long dt)
/*
**------------------------------------------------------------
** Input: k = link index
@@ -271,10 +276,10 @@ double pipereact(EN_Project *pr, int k, double c, double v, long dt)
**------------------------------------------------------------
*/
{
double cnew, dc, dcbulk, dcwall, rbulk, rwall;
Network *net = &pr->network;
Quality *qual = &pr->quality;
EN_Network *net = &pr->network;
quality_t *qual = &pr->quality;
double cnew, dc, dcbulk, dcwall, rbulk, rwall;
// For water age (hrs), update concentration by timestep
if (qual->Qualflag == AGE)
@@ -294,7 +299,7 @@ double pipereact(EN_Project *pr, int k, double c, double v, long dt)
dcwall = rwall * (double)dt;
// Update cumulative mass reacted
if (pr->time_options.Htime >= pr->time_options.Rstart)
if (pr->times.Htime >= pr->times.Rstart)
{
qual->Wbulk += fabs(dcbulk) * v;
qual->Wwall += fabs(dcwall) * v;
@@ -308,7 +313,7 @@ double pipereact(EN_Project *pr, int k, double c, double v, long dt)
}
double tankreact(EN_Project *pr, double c, double v, double kb, long dt)
double tankreact(Project *pr, double c, double v, double kb, long dt)
/*
**-------------------------------------------------------
** Input: c = current quality in tank
@@ -321,8 +326,9 @@ double tankreact(EN_Project *pr, double c, double v, double kb, long dt)
**-------------------------------------------------------
*/
{
Quality *qual = &pr->quality;
double cnew, dc, rbulk;
quality_t *qual = &pr->quality;
// For water age, update concentration by timestep
if (qual->Qualflag == AGE)
@@ -338,7 +344,7 @@ double tankreact(EN_Project *pr, double c, double v, double kb, long dt)
// Find concentration change & update quality
dc = rbulk * (double)dt;
if (pr->time_options.Htime >= pr->time_options.Rstart)
if (pr->times.Htime >= pr->times.Rstart)
{
qual->Wtank += fabs(dc) * v;
}
@@ -349,7 +355,7 @@ double tankreact(EN_Project *pr, double c, double v, double kb, long dt)
}
double bulkrate(EN_Project *pr, double c, double kb, double order)
double bulkrate(Project *pr, double c, double kb, double order)
/*
**-----------------------------------------------------------
** Input: c = current WQ concentration
@@ -360,8 +366,9 @@ double bulkrate(EN_Project *pr, double c, double kb, double order)
**-----------------------------------------------------------
*/
{
Quality *qual = &pr->quality;
double c1;
quality_t *qual = &pr->quality;
// Find bulk reaction potential taking into account
// limiting potential & reaction order.
@@ -396,7 +403,7 @@ double bulkrate(EN_Project *pr, double c, double kb, double order)
}
double wallrate(EN_Project *pr, double c, double d, double kw, double kf)
double wallrate(Project *pr, double c, double d, double kw, double kf)
/*
**------------------------------------------------------------
** Input: c = current WQ concentration
@@ -410,7 +417,8 @@ double wallrate(EN_Project *pr, double c, double d, double kw, double kf)
**------------------------------------------------------------
*/
{
quality_t *qual = &pr->quality;
Quality *qual = &pr->quality;
if (kw == 0.0 || d == 0.0) return (0.0);
if (qual->WallOrder == 0.0) // 0-order reaction */
@@ -424,7 +432,7 @@ double wallrate(EN_Project *pr, double c, double d, double kw, double kf)
}
double mixtank(EN_Project *pr, int n, double volin, double massin, double volout)
double mixtank(Project *pr, int n, double volin, double massin, double volout)
/*
**------------------------------------------------------------
** Input: n = node index
@@ -436,11 +444,11 @@ double mixtank(EN_Project *pr, int n, double volin, double massin, double volout
**------------------------------------------------------------
*/
{
Network *net = &pr->network;
Quality *qual = &pr->quality;
int i;
double vnet;
EN_Network *net = &pr->network;
quality_t *qual = &pr->quality;
i = n - net->Njuncs;
vnet = volin - volout;
switch (net->Tank[i].MixModel)
@@ -454,7 +462,7 @@ double mixtank(EN_Project *pr, int n, double volin, double massin, double volout
}
void tankmix1(EN_Project *pr, int i, double vin, double win, double vnet)
void tankmix1(Project *pr, int i, double vin, double win, double vnet)
/*
**---------------------------------------------
** Input: i = tank index
@@ -466,14 +474,14 @@ void tankmix1(EN_Project *pr, int i, double vin, double win, double vnet)
**---------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Quality *qual = &pr->quality;
int k;
double vnew;
Pseg seg;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
Stank *tank = &net->Tank[i];
Stank *tank = &net->Tank[i];
k = net->Nlinks + i;
seg = qual->FirstSeg[k];
@@ -488,7 +496,7 @@ void tankmix1(EN_Project *pr, int i, double vin, double win, double vnet)
}
void tankmix2(EN_Project *pr, int i, double vin, double win, double vnet)
void tankmix2(Project *pr, int i, double vin, double win, double vnet)
/*
**------------------------------------------------
** Input: i = tank index
@@ -500,15 +508,15 @@ void tankmix2(EN_Project *pr, int i, double vin, double win, double vnet)
**------------------------------------------------
*/
{
int k;
Network *net = &pr->network;
Quality *qual = &pr->quality;
int k;
double vt, // Transferred volume
vmz; // Full mixing zone volume
Pseg mixzone, // Mixing zone segment
stagzone; // Stagnant zone segment
EN_Network *net = &pr->network;
quality_t *qual = &pr->quality;
Stank *tank = &pr->network.Tank[i];
Stank *tank = &pr->network.Tank[i];
// Identify segments for each compartment
k = net->Nlinks + i;
@@ -568,7 +576,7 @@ void tankmix2(EN_Project *pr, int i, double vin, double win, double vnet)
}
void tankmix3(EN_Project *pr, int i, double vin, double win, double vnet)
void tankmix3(Project *pr, int i, double vin, double win, double vnet)
/*
**----------------------------------------------------------
** Input: i = tank index
@@ -580,15 +588,15 @@ void tankmix3(EN_Project *pr, int i, double vin, double win, double vnet)
**----------------------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Quality *qual = &pr->quality;
int k;
double vout, vseg;
double cin, vsum, wsum;
Pseg seg;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
Stank *tank = &pr->network.Tank[i];
Stank *tank = &pr->network.Tank[i];
k = net->Nlinks + i;
if (qual->LastSeg[k] == NULL || qual->FirstSeg[k] == NULL) return;
@@ -639,7 +647,7 @@ void tankmix3(EN_Project *pr, int i, double vin, double win, double vnet)
}
void tankmix4(EN_Project *pr, int i, double vin, double win, double vnet)
void tankmix4(Project *pr, int i, double vin, double win, double vnet)
/*
**----------------------------------------------------------
** Input: i = tank index
@@ -651,13 +659,13 @@ void tankmix4(EN_Project *pr, int i, double vin, double win, double vnet)
**----------------------------------------------------------
*/
{
Network *net = &pr->network;
Quality *qual = &pr->quality;
int k, n;
double cin, vsum, wsum, vseg;
Pseg seg;
EN_Network *net = &pr->network;
quality_t *qual = &pr->quality;
Stank *tank = &pr->network.Tank[i];
Stank *tank = &pr->network.Tank[i];
k = net->Nlinks + i;
if (qual->LastSeg[k] == NULL || qual->FirstSeg[k] == NULL) return;
src/qualroute.c
+143 -219
View File
@@ -1,46 +1,55 @@
/*
*********************************************************************
QUALROUTE.C -- water quality routing module for the EPANET program
*********************************************************************
******************************************************************************
Project: OWA EPANET
Version: 2.2
Module: qualroute.c
Description: computes water quality transport over a single time step
Authors: see AUTHORS
Copyright: see AUTHORS
License: see LICENSE
Last Updated: 11/27/2018
******************************************************************************
*/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <stdio.h>
#ifndef __APPLE__
#include <malloc.h>
#else
#include <stdlib.h>
#endif
#include <math.h>
#include "mempool.h"
#include "types.h"
// Macro to compute the volume of a link
#define LINKVOL(k) (0.785398 * net->Link[(k)].Len * SQR(net->Link[(k)].Diam))
// Macro to get link flow compatible with flow saved to hydraulics file
#define LINKFLOW(k) ((hyd->LinkStatus[k] <= CLOSED) ? 0.0 : hyd->LinkFlows[k])
#define LINKFLOW(k) ((hyd->LinkStatus[k] <= CLOSED) ? 0.0 : hyd->LinkFlow[k])
// Exported Functions
int buildilists(EN_Project *pr);
int sortnodes(EN_Project *pr);
void transport(EN_Project *pr, long);
void initsegs(EN_Project *pr);
void reversesegs(EN_Project *pr, int);
void addseg(EN_Project *pr, int, double, double);
// Exported functions
int sortnodes(Project *);
void transport(Project *, long);
void initsegs(Project *);
void reversesegs(Project *, int);
void addseg(Project *, int, double, double);
// Imported Functions
extern double findsourcequal(EN_Project *pr, int, double, double, long);
extern void reactpipes(EN_Project *pr, long);
extern void reacttanks(EN_Project *pr, long);
extern double mixtank(EN_Project *pr, int, double, double, double);
// Imported functions
extern double findsourcequal(Project *, int, double, double, long);
extern void reactpipes(Project *, long);
extern void reacttanks(Project *, long);
extern double mixtank(Project *, int, double, double, double);
// Local Functions
static void evalnodeinflow(EN_Project *pr, int, long, double *, double *);
static void evalnodeoutflow(EN_Project *pr, int, double, long);
static double findnodequal(EN_Project *pr, int, double, double, double, long);
static double noflowqual(EN_Project *pr, int);
static void updatemassbalance(EN_Project *pr, int, double, double, long);
static int selectnonstacknode(EN_Project *pr, int, int *);
// Local functions
static void evalnodeinflow(Project *, int, long, double *, double *);
static void evalnodeoutflow(Project *, int, double, long);
static double findnodequal(Project *, int, double, double, double, long);
static double noflowqual(Project *, int);
static void updatemassbalance(Project *, int, double, double, long);
static int selectnonstacknode(Project *, int, int *);
void transport(EN_Project *pr, long tstep)
void transport(Project *pr, long tstep)
/*
**--------------------------------------------------------------
** Input: tstep = length of current time step
@@ -50,12 +59,13 @@ void transport(EN_Project *pr, long tstep)
**--------------------------------------------------------------
*/
{
int i, j, k, m, n;
double volin, massin, volout, nodequal;
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Quality *qual = &pr->quality;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
EN_Network *net = &pr->network;
int j, k, m, n;
double volin, massin, volout, nodequal;
Padjlist alink;
// React contents of each pipe and tank
if (qual->Reactflag)
@@ -76,10 +86,10 @@ void transport(EN_Project *pr, long tstep)
volout = 0.0;
// ... examine each link with flow into the node
for (i = qual->IlistPtr[n]; i < qual->IlistPtr[n + 1]; i++)
for (alink = net->Adjlist[n]; alink != NULL; alink = alink->next)
{
// ... k is index of next link incident on node n
k = qual->Ilist[i];
k = alink->link;
// ... link has flow into node - add it to node's inflow
// (m is index of link's downstream node)
@@ -107,10 +117,10 @@ void transport(EN_Project *pr, long tstep)
nodequal = findnodequal(pr, n, volin, massin, volout, tstep);
// ... examine each link with flow out of the node
for (i = qual->IlistPtr[n]; i < qual->IlistPtr[n + 1]; i++)
for (alink = net->Adjlist[n]; alink != NULL; alink = alink->next)
{
// ... link k incident on node n has upstream node m equal to n
k = qual->Ilist[i];
k = alink->link;
m = net->Link[k].N1;
if (qual->FlowDir[k] < 0) m = net->Link[k].N2;
if (m == n)
@@ -123,27 +133,27 @@ void transport(EN_Project *pr, long tstep)
}
}
void evalnodeinflow(EN_Project *pr, int k, long tstep, double *volin,
double *massin)
/*
**--------------------------------------------------------------
** Input: k = link index
** tstep = quality routing time step
** Output: volin = flow volume entering a node
** massin = constituent mass entering a node
** Purpose: adds the contribution of a link's outflow volume
** and constituent mass to the total inflow into its
** downstream node over a time step.
**--------------------------------------------------------------
*/
void evalnodeinflow(Project *pr, int k, long tstep, double *volin,
double *massin)
/*
**--------------------------------------------------------------
** Input: k = link index
** tstep = quality routing time step
** Output: volin = flow volume entering a node
** massin = constituent mass entering a node
** Purpose: adds the contribution of a link's outflow volume
** and constituent mass to the total inflow into its
** downstream node over a time step.
**--------------------------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Quality *qual = &pr->quality;
double q, v, vseg;
Pseg seg;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
// Get flow rate (q) and flow volume (v) through link
q = LINKFLOW(k);
v = fabs(q) * tstep;
@@ -185,24 +195,24 @@ void evalnodeinflow(EN_Project *pr, int k, long tstep, double *volin,
}
double findnodequal(EN_Project *pr, int n, double volin,
double massin, double volout, long tstep)
/*
**--------------------------------------------------------------
** Input: n = node index
** volin = flow volume entering node
** massin = mass entering node
** volout = flow volume leaving node
** tstep = length of current time step
** Output: returns water quality in a node's outflow
** Purpose: computes a node's new quality from its inflow
** volume and mass, including any source contribution.
**--------------------------------------------------------------
*/
double findnodequal(Project *pr, int n, double volin,
double massin, double volout, long tstep)
/*
**--------------------------------------------------------------
** Input: n = node index
** volin = flow volume entering node
** massin = mass entering node
** volout = flow volume leaving node
** tstep = length of current time step
** Output: returns water quality in a node's outflow
** Purpose: computes a node's new quality from its inflow
** volume and mass, including any source contribution.
**--------------------------------------------------------------
*/
{
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Quality *qual = &pr->quality;
// Node is a junction - update its water quality
if (net->Node[n].Type == JUNCTION)
@@ -262,7 +272,7 @@ double findnodequal(EN_Project *pr, int n, double volin,
}
double noflowqual(EN_Project *pr, int n)
double noflowqual(Project *pr, int n)
/*
**--------------------------------------------------------------
** Input: n = node index
@@ -274,33 +284,34 @@ double noflowqual(EN_Project *pr, int n)
**--------------------------------------------------------------
*/
{
int i, k, inflow, kount = 0;
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Quality *qual = &pr->quality;
int k, inflow, kount = 0;
double c = 0.0;
FlowDirection dir;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
Padjlist alink;
// Examine each link incident on the node
for (i = qual->IlistPtr[n]; i < qual->IlistPtr[n + 1]; i++)
for (alink = net->Adjlist[n]; alink != NULL; alink = alink->next)
{
// ... index of an incident link
k = qual->Ilist[i];
k = alink->link;
dir = qual->FlowDir[k];
// Node n is link's downstream node - add quality
// Node n is link's downstream node - add quality
// of link's first segment to average
if (net->Link[k].N2 == n && dir >= 0) inflow = TRUE;
else if (net->Link[k].N1 == n && dir < 0) inflow = TRUE;
else inflow = FALSE;
else inflow = FALSE;
if (inflow == TRUE && qual->FirstSeg[k] != NULL)
{
c += qual->FirstSeg[k]->c;
kount++;
}
// Node n is link's upstream node - add quality
// Node n is link's upstream node - add quality
// of link's last segment to average
else if (inflow == FALSE && qual->LastSeg[k] != NULL)
{
@@ -313,7 +324,7 @@ double noflowqual(EN_Project *pr, int n)
}
void evalnodeoutflow(EN_Project *pr, int k, double c, long tstep)
void evalnodeoutflow(Project *pr, int k, double c, long tstep)
/*
**--------------------------------------------------------------
** Input: k = link index
@@ -325,12 +336,12 @@ void evalnodeoutflow(EN_Project *pr, int k, double c, long tstep)
**--------------------------------------------------------------
*/
{
Hydraul *hyd = &pr->hydraul;
Quality *qual = &pr->quality;
double v;
Pseg seg;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
// Find flow volume (v) released over time step
v = fabs(LINKFLOW(k)) * tstep;
if (v == 0.0) return;
@@ -358,27 +369,26 @@ void evalnodeoutflow(EN_Project *pr, int k, double c, long tstep)
}
void updatemassbalance(EN_Project *pr, int n, double massin,
double volout, long tstep)
/*
**--------------------------------------------------------------
** Input: n = node index
** massin = mass inflow to node
** volout = outflow volume from node
** Output: none
** Purpose: Adds a node's external mass inflow and outflow
** over the current time step to the network's
** overall mass balance.
**--------------------------------------------------------------
*/
void updatemassbalance(Project *pr, int n, double massin,
double volout, long tstep)
/*
**--------------------------------------------------------------
** Input: n = node index
** massin = mass inflow to node
** volout = outflow volume from node
** Output: none
** Purpose: Adds a node's external mass inflow and outflow
** over the current time step to the network's
** overall mass balance.
**--------------------------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Quality *qual = &pr->quality;
double masslost = 0.0,
massadded = 0.0;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
massadded = 0.0;
switch (net->Node[n].Type)
{
@@ -401,82 +411,12 @@ void updatemassbalance(EN_Project *pr, int n, double massin,
massadded = qual->SourceQual * volout;
break;
}
qual->massbalance.outflow += masslost;
qual->massbalance.inflow += massadded;
qual->MassBalance.outflow += masslost;
qual->MassBalance.inflow += massadded;
}
int buildilists(EN_Project *pr)
/*
**--------------------------------------------------------------
** Input: none
** Output: returns an error code
** Purpose: Creates link incidence lists for each node stored
** in compact form.
**--------------------------------------------------------------
*/
{
int i, j, k, n, n1, n2;
int *degree = NULL;
quality_t *qual = &pr->quality;
EN_Network *net = &pr->network;
// Allocate an array to count # links incident on each node
n = net->Nnodes + 1;
degree = (int *)calloc(n, sizeof(int));
if (degree == NULL) return 101;
// Count # links incident on each node
for (k = 1; k <= net->Nlinks; k++)
{
degree[net->Link[k].N1]++;
degree[net->Link[k].N2]++;
}
// Use incidence counts to determine start position of
// each node's incidence list in Xilist
qual->IlistPtr[1] = 1;
for (i = 1; i <= n; i++)
{
qual->IlistPtr[i + 1] = qual->IlistPtr[i] + degree[i];
}
// Add each link to the incidence lists of its start & end nodes
for (i = 1; i <= net->Nnodes; i++) degree[i] = 0;
for (k = 1; k <= net->Nlinks; k++)
{
// j is index of next unused location in link's start node list
n1 = net->Link[k].N1;
j = qual->IlistPtr[n1] + degree[n1];
qual->Ilist[j] = k;
degree[n1]++;
// Repeat same for end node
n2 = net->Link[k].N2;
j = qual->IlistPtr[n2] + degree[n2];
qual->Ilist[j] = k;
degree[n2]++;
}
free(degree);
/*//////// QA CHECK
for (i = 1; i <= net->Nnodes; i++)
{
printf("\nNode %s: ", net->Node[i].ID);
for (j = qual->IlistPtr[i]; j < qual->IlistPtr[i + 1]; j++)
{
printf(" %s,", net->Link[qual->Ilist[j]].ID);
}
}
*/
return 0;
}
int sortnodes(EN_Project *pr)
int sortnodes(Project *pr)
/*
**--------------------------------------------------------------
** Input: none
@@ -487,6 +427,10 @@ int sortnodes(EN_Project *pr)
**--------------------------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Quality *qual = &pr->quality;
int i, j, k, n;
int *indegree = NULL;
int *stack = NULL;
@@ -494,10 +438,7 @@ int sortnodes(EN_Project *pr)
int numsorted = 0;
int errcode = 0;
FlowDirection dir;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
Padjlist alink;
// Allocate an array to count # links with inflow to each node
// and for a stack to hold nodes waiting to be processed
@@ -548,10 +489,10 @@ int sortnodes(EN_Project *pr)
// ... for each outflow link from this node reduce the in-degree
// of its downstream node
for (j = qual->IlistPtr[i]; j < qual->IlistPtr[i + 1]; j++)
for (alink = net->Adjlist[i]; alink != NULL; alink = alink->next)
{
// ... k is the index of the next link incident on node i
k = qual->Ilist[j];
k = alink->link;
// ... skip link if flow is negligible
if (qual->FlowDir[k] == 0) continue;
@@ -579,25 +520,11 @@ int sortnodes(EN_Project *pr)
if (numsorted < net->Nnodes) errcode = 120;
FREE(indegree);
FREE(stack);
/*
/////////////////// QA CHECK
snprintf(pr->Msg, MAXMSG, "\n\nSorted Nodes:");
writeline(pr, pr->Msg);
for (i = 1; i <= net->Nnodes; i++)
{
j = qual->SortedNodes[i];
snprintf(pr->Msg, MAXMSG, "%s", net->Node[j].ID);
writeline(pr, pr->Msg);
}
//printf("\n");
//system("pause");
/////////////////
*/
return errcode;
}
int selectnonstacknode(EN_Project *pr, int numsorted, int *indegree)
int selectnonstacknode(Project *pr, int numsorted, int *indegree)
/*
**--------------------------------------------------------------
** Input: numsorted = number of nodes that have been sorted
@@ -607,24 +534,21 @@ int selectnonstacknode(EN_Project *pr, int numsorted, int *indegree)
**--------------------------------------------------------------
*/
{
int i, j, k, m, n;
Network *net = &pr->network;
Quality *qual = &pr->quality;
quality_t *qual = &pr->quality;
EN_Network *net = &pr->network;
int i, m, n;
Padjlist alink;
// Examine each sorted node in last in - first out order
for (i = numsorted; i > 0; i--)
{
// For each link connected to the sorted node
m = qual->SortedNodes[i];
for (j = qual->IlistPtr[m]; j < qual->IlistPtr[m + 1]; j++)
for (alink = net->Adjlist[m]; alink != NULL; alink = alink->next)
{
// ... k is index of next link incident on node m
k = qual->Ilist[j];
// ... n is the node of link k opposite to node m
n = net->Link[k].N2;
if (n == m) n = net->Link[k].N1;
n = alink->node;
// ... select node n if it still has inflow links
if (indegree[n] > 0) return n;
@@ -643,7 +567,7 @@ int selectnonstacknode(EN_Project *pr, int numsorted, int *indegree)
}
void initsegs(EN_Project *pr)
void initsegs(Project *pr)
/*
**--------------------------------------------------------------
** Input: none
@@ -653,13 +577,13 @@ void initsegs(EN_Project *pr)
**--------------------------------------------------------------
*/
{
Network *net = &pr->network;
Hydraul *hyd = &pr->hydraul;
Quality *qual = &pr->quality;
int j, k;
double c, v, v1;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
// Add one segment with assigned downstream node quality to each pipe
for (k = 1; k <= net->Nlinks; k++)
{
@@ -706,7 +630,7 @@ void initsegs(EN_Project *pr)
}
void reversesegs(EN_Project *pr, int k)
void reversesegs(Project *pr, int k)
/*
**--------------------------------------------------------------
** Input: k = link index
@@ -715,8 +639,8 @@ void reversesegs(EN_Project *pr, int k)
**--------------------------------------------------------------
*/
{
Pseg seg, nseg, pseg;
quality_t *qual = &pr->quality;
Quality *qual = &pr->quality;
Pseg seg, nseg, pseg;
seg = qual->FirstSeg[k];
qual->FirstSeg[k] = qual->LastSeg[k];
@@ -732,7 +656,7 @@ void reversesegs(EN_Project *pr, int k)
}
void addseg(EN_Project *pr, int k, double v, double c)
void addseg(Project *pr, int k, double v, double c)
/*
**-------------------------------------------------------------
** Input: k = segment chain index
@@ -744,8 +668,8 @@ void addseg(EN_Project *pr, int k, double v, double c)
**-------------------------------------------------------------
*/
{
Quality *qual = &pr->quality;
Pseg seg;
quality_t *qual = &pr->quality;
// Grab the next free segment from the segment pool if available
if (qual->FreeSeg != NULL)
src/report.c
+931 -936
View File
File diff suppressed because it is too large Load Diff
src/rules.c
+760 -872
View File
File diff suppressed because it is too large Load Diff
src/smatrix.c
+330 -373
View File
File diff suppressed because it is too large Load Diff
src/text.h
+39 -33
View File
@@ -1,22 +1,21 @@
/*
****************************************************
String Constants for EPANET Program
VERSION: 2.00
DATE: 5/8/00
10/25/00
8/15/07 (2.00.11)
2/14/08 (2.00.12)
AUTHOR: L. Rossman
US EPA - NRMRL
****************************************************
******************************************************************************
Project: OWA EPANET
Version: 2.2
Module: text.h
Description: string constants used throughout EPANET
Authors: see AUTHORS
Copyright: see AUTHORS
License: see LICENSE
Last Updated: 11/27/2018
******************************************************************************
*/
/* ------------ Keyword Dictionary ---------- */
#ifndef TEXT_H
#define TEXT_H
//------- Keyword Dictionary ------------------------------
#define w_USE "USE"
#define w_SAVE "SAVE"
@@ -194,7 +193,8 @@ AUTHOR: L. Rossman
#define w_ELSE "ELSE"
#define w_PRIORITY "PRIO"
/* ---------Input Section Names ---------- */
// ------ Input File Section Names ------------------------
#define s_TITLE "[TITLE]"
#define s_JUNCTIONS "[JUNCTIONS]"
#define s_RESERVOIRS "[RESERVOIRS]"
@@ -225,8 +225,8 @@ AUTHOR: L. Rossman
#define s_TAGS "[TAGS]"
#define s_END "[END]"
/* ---------------- Units ---------------- */
/*** Limit units to MAXID or less characters ***/
//------- Units -------------------------------------------
#define u_CFS "cfs"
#define u_GPM "gpm"
#define u_AFD "a-f/d"
@@ -257,13 +257,15 @@ AUTHOR: L. Rossman
#define u_per1000FT "/1000ft"
#define u_HP "hp"
/* -------------- Curve Types ----------------- */
//------- Curve Types -------------------------------------
#define c_HEADLOSS "HEADLOSS"
#define c_PUMP "PUMP"
#define c_EFFIC "EFFIC"
#define c_VOLUME "VOLUME"
/* ------------------ Text Phrases ------------------- */
//------- Text Phrases ------------------------------------
#define t_ABOVE "above"
#define t_BELOW "below"
#define t_HW "Hazen-Williams"
@@ -334,7 +336,8 @@ AUTHOR: L. Rossman
#define t_ORIFICE "Orifice Flow"
/* ------------------ Format Messages ------------------*/
//----- Summary Report Format Strings ---------------------
#define LOGO1 \
"******************************************************************"
#define LOGO2 \
@@ -399,10 +402,9 @@ AUTHOR: L. Rossman
#define FMT47 " with %s below %-.2f %s"
#define FMT48 " with %s above %-.2f %s"
/* ---------- Status Report Format Strings ------------ */
#define FMT49 "Hydraulic Status:"
//----- Status Report Format Strings ----------------------
/*** Updated 6/24/02 ***/
#define FMT49 "Hydraulic Status:"
#define FMT50 "%10s: Tank %s is %s at %-.2f %s"
#define FMT51 "%10s: Reservoir %s is %s"
#define FMT52 "%10s: %s %s %s"
@@ -422,13 +424,12 @@ AUTHOR: L. Rossman
#define FMT63 "%10s: %s %s changed by rule %s"
#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 flow change = %.4f for Node %s"
#define FMT68 " maximum head error = %.4f for Link %s\n"
/* -------------------- Energy Report Table ------------------- */
//----- Energy Report Table -------------------------------
#define FMT71 "Energy Usage:"
#define FMT72 \
" Usage Avg. Kw-hr Avg. Peak Cost"
@@ -437,18 +438,21 @@ AUTHOR: L. Rossman
#define FMT74 "%38s Demand Charge: %9.2f"
#define FMT75 "%38s Total Cost: %9.2f"
/* -------------------- Node Report Table --------------------- */
//----- Node Report Table ---------------------------------
#define FMT76 "%s Node Results:"
#define FMT77 "Node Results:"
#define FMT78 "Node Results at %s hrs:"
/* -------------------- Link Report Table --------------------- */
//----- Link Report Table ---------------------------------
#define FMT79 "%s Link Results:"
#define FMT80 "Link Results:"
#define FMT81 "Link Results at %s hrs:"
#define FMT82 "\n\f\n Page %-d %60.60s\n"
/* ------------------- Progress Messages ---------------------- */
//----- Progress Messages ---------------------------------
#define FMT100 " Retrieving network data ... "
#define FMT101 " Computing hydraulics at hour %-10s "
#define FMT102 " Computing water quality at hour %-10s "
@@ -457,16 +461,17 @@ AUTHOR: L. Rossman
#define FMT104 "Analysis begun %s"
#define FMT105 "Analysis ended %s"
/*------------------- Error Messages --------------------*/
//----- Rule Error Messages -------------------------------
#define R_ERR201 "Input Error 201: syntax error in following line of "
#define R_ERR202 "Input Error 202: illegal numeric value in following line of "
#define R_ERR203 "Input Error 203: undefined node in following line of "
#define R_ERR204 "Input Error 204: undefined link in following line of "
#define R_ERR207 "Input Error 207: attempt to control a CV in following line of "
#define R_ERR221 "Input Error 221: mis-placed clause in following line of "
/*-------------------- Specific Warning Messages -------------------------*/
//----- Specific Warning Messages -------------------------
#define WARN01 "WARNING: System unbalanced at %s hrs."
#define WARN02 \
"WARNING: Maximum trials exceeded at %s hrs. System may be unstable."
@@ -477,7 +482,8 @@ AUTHOR: L. Rossman
#define WARN05 "WARNING: %s %s %s at %s hrs."
#define WARN06 "WARNING: Negative pressures at %s hrs."
/*-------------------- General Warning Messages -------------------------*/
//----- General Warning Messages --------------------------
#define WARN1 "WARNING: System hydraulically unbalanced."
#define WARN2 "WARNING: System may be hydraulically unstable."
#define WARN3 "WARNING: System disconnected."
src/types.h
+663 -724
View File
File diff suppressed because it is too large Load Diff
tests/test_addrule.cpp
+1 -1
View File
@@ -45,7 +45,7 @@ int main(int argc, char *argv[])
int link113, node23, link22, pump9_before, pump9_after;
float priority;
EN_ProjectHandle ph = NULL;
EN_Project ph = NULL;
EN_createproject(&ph);
std::string path_inp = std::string(DATA_PATH_INP);
tests/test_demand_categories.cpp
+1 -1
View File
@@ -45,7 +45,7 @@ BOOST_AUTO_TEST_CASE(test_demand_categories)
int error = 0;
int Nindex, ndem;
EN_ProjectHandle ph = NULL;
EN_Project ph = NULL;
error = EN_createproject(&ph);
error = EN_open(ph, path_inp.c_str(), path_rpt.c_str(), path_out.c_str());
tests/test_net_builder.cpp
+1 -1
View File
@@ -37,7 +37,7 @@ BOOST_AUTO_TEST_CASE(test_net_builder)
float h_orig, h_build, h_build_loaded;
// first we load Net1.inp, run it and record the head in Tank 2 at the end of the simulation (h_orig)
EN_ProjectHandle ph = NULL;
EN_Project ph = NULL;
EN_createproject(&ph);
std::string path_inp = std::string(DATA_PATH_INP);
tests/test_reent.cpp
+1 -1
View File
@@ -23,7 +23,7 @@ using namespace std;
void epanet_thread(long i)
{
int errorcode = 0;
EN_ProjectHandle ph;
EN_Project ph;
string prefix = "example_";
string suffix = ".inp";
tests/test_setid.cpp
+22 -9
View File
@@ -5,10 +5,14 @@
This is a test for the API functions that change a node or link ID name.
A node and link name are changed, the network is saved, reopened and the new names are checked.
*/
//#define NO_BOOST
#ifndef NO_BOOST
#define BOOST_TEST_MODULE "toolkit"
#include <boost/test/included/unit_test.hpp>
#endif
#include <iostream>
#include <string>
#include "epanet2.h"
@@ -16,20 +20,28 @@ A node and link name are changed, the network is saved, reopened and the new nam
#define DATA_PATH_RPT "./test.rpt"
#define DATA_PATH_OUT "./test.out"
#ifdef NO_BOOST
#define BOOST_REQUIRE(x) (((x)) ? cout << "\nPassed at line " << __LINE__ : cout << "\nFailed at line " << __LINE__ )
#endif
using namespace std;
#ifndef NO_BOOST
BOOST_AUTO_TEST_SUITE (test_toolkit)
BOOST_AUTO_TEST_CASE(test_setid)
#else
int main(int argc, char *argv[])
#endif
{
string path_inp(DATA_PATH_INP);
string path_rpt(DATA_PATH_RPT);
string path_out(DATA_PATH_OUT);
string inp_save("net1_setid.inp");
int error = 0;
int index;
EN_ProjectHandle ph = NULL;
EN_Project ph = NULL;
EN_createproject(&ph);
error = EN_open(ph, path_inp.c_str(), path_rpt.c_str(), "");
@@ -56,32 +68,33 @@ BOOST_AUTO_TEST_CASE(test_setid)
BOOST_REQUIRE(error == 0);
// Save the project
error = EN_saveinpfile(ph, "net1_setid.inp");
error = EN_saveinpfile(ph, inp_save.c_str());
BOOST_REQUIRE(error == 0);
error = EN_close(ph);
BOOST_REQUIRE(error == 0);
EN_deleteproject(&ph);
// Re-open the saved project
// Re-open the saved project
EN_createproject(&ph);
error = EN_open(ph, "net1_setid.inp", path_rpt.c_str(), "");
error = EN_open(ph, inp_save.c_str(), path_rpt.c_str(), "");
BOOST_REQUIRE(error == 0);
// Check that 3rd node has its new name
error = EN_getnodeindex(ph, newid_2, &index);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(index == 3);
BOOST_REQUIRE(index == 3);
// Check that 3rd link has its new name
error = EN_getlinkindex(ph, newid_4, &index);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(index == 3);
BOOST_REQUIRE(index == 3);
error = EN_close(ph);
BOOST_REQUIRE(error == 0);
EN_deleteproject(&ph);
}
#ifndef NO_BOOST
BOOST_AUTO_TEST_SUITE_END()
#endif
tests/test_setlinktype.cpp
+1 -1
View File
@@ -29,7 +29,7 @@ BOOST_AUTO_TEST_CASE(test_setlinktype)
int p113, n31, p121, n113_1, n113_2;
float q113 = 0.0f, p31 = 0.0f, diam;
EN_ProjectHandle ph = NULL;
EN_Project ph = NULL;
EN_createproject(&ph);
std::string path_inp = std::string(DATA_PATH_INP);
tests/test_toolkit.cpp
+5 -4
View File
@@ -36,7 +36,7 @@ BOOST_AUTO_TEST_SUITE (test_toolkit)
BOOST_AUTO_TEST_CASE (test_alloc_free)
{
int error = 0;
EN_ProjectHandle ph = NULL;
EN_Project ph = NULL;
error = EN_createproject(&ph);
@@ -55,7 +55,7 @@ BOOST_AUTO_TEST_CASE (test_open_close)
string path_rpt(DATA_PATH_RPT);
string path_out(DATA_PATH_OUT);
EN_ProjectHandle ph = NULL;
EN_Project ph = NULL;
EN_createproject(&ph);
int error = EN_open(ph, path_inp.c_str(), path_rpt.c_str(), path_out.c_str());
@@ -74,7 +74,8 @@ BOOST_AUTO_TEST_CASE(test_save_reopen)
string path_rpt(DATA_PATH_RPT);
string path_out(DATA_PATH_OUT);
EN_ProjectHandle ph_save, ph_reopen;
EN_Project ph_save;
EN_Project ph_reopen;
EN_createproject(&ph_save);
@@ -126,7 +127,7 @@ struct Fixture{
string path_out;
int error;
EN_ProjectHandle ph;
EN_Project ph;
};
BOOST_AUTO_TEST_SUITE(test_epanet_fixture)
win_build/WinSDK/Makefile.bat
+4 -4
View File
@@ -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 util\errormanager.c hash.c hydraul.c hydcoeffs.c hydstatus.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c qualroute.c qualreact.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /link /DLL
cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydstatus.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c project.c quality.c qualroute.c qualreact.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /link /DLL
rem : create EPANET2.EXE
cl -o epanet2.exe epanet.c util\errormanager.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydstatus.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c qualroute.c qualreact.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /I ..\src /link
cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydstatus.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c project.c quality.c qualroute.c qualreact.c report.c rules.c smatrix.c genmmd.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 util\errormanager.c hash.c hydraul.c hydcoeffs.c hydstatus.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c qualroute.c qualreact.c report.c rules.c smatrix.c genmmd.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 hydstatus.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c project.c quality.c qualroute.c qualreact.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /link /DLL /def:..\win_build\WinSDK\epanet2.def /MAP
rem : create EPANET2.EXE
cl -o epanet2.exe epanet.c util\errormanager.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydstatus.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c qualroute.c qualreact.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /I ..\src /link
cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydstatus.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c project.c quality.c qualroute.c qualreact.c report.c rules.c smatrix.c genmmd.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