From b3ab8ea2c71100b472c930d423505aa087ed6ff7 Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Sat, 16 Jun 2018 11:02:18 -0400 Subject: [PATCH 1/9] Addresses issue #161 Adds new options HEADERROR and FLOWCHANGE to provide more rigorous criteria for hydraulic convergence. Also breaks HYDRAUL.C into 3 separate files to improve code readability. --- src/epanet.c | 59 +- src/epanet2.h | 1236 +++++++++++++++++++++++++++++ src/funcs.h | 111 +-- src/hydcoeffs.c | 919 +++++++++++++++++++++ src/hydraul.c | 2017 ++++------------------------------------------- src/hydsolver.c | 1052 ++++++++++++++++++++++++ src/inpfile.c | 2 + src/input1.c | 33 +- src/input3.c | 59 +- src/main.c | 135 ++++ src/output.c | 91 ++- src/report.c | 20 +- src/text.h | 15 +- src/types.h | 169 ++-- 14 files changed, 3800 insertions(+), 2118 deletions(-) create mode 100644 src/epanet2.h create mode 100644 src/hydcoeffs.c create mode 100644 src/hydsolver.c create mode 100644 src/main.c diff --git a/src/epanet.c b/src/epanet.c index 602e0ed..f0f0ff8 100644 --- a/src/epanet.c +++ b/src/epanet.c @@ -71,11 +71,12 @@ This module calls the following functions that reside in other modules: runhyd() nexthyd() closehyd() - resistance() tankvolume() getenergy() setlinkstatus() setlinksetting() + HYDCOEFFS + resistcoeff() QUALITY.C openqual() initqual() @@ -129,6 +130,7 @@ execute function x and set the error code equal to its return value. #include "text.h" #include "types.h" #define EXTERN +////////////////////////////////////////////#include "epanet2.h" #include "vars.h" /**************************************************************** @@ -1341,6 +1343,14 @@ int DLLEXPORT EN_getoption(EN_Project *pr, EN_Option code, case EN_DEMANDMULT: v = hyd->Dmult; break; + + case EN_HEADERROR: + v = hyd->HeadErrorLimit * Ucf[HEAD]; + break; + case EN_FLOWCHANGE: + v = hyd->FlowChangeLimit * Ucf[FLOW]; + break; + default: return (251); } @@ -2648,7 +2658,7 @@ int DLLEXPORT EN_setlinkvalue(EN_Project *p, int index, int code, r = Link[index].Diam / value; /* Ratio of old to new diam */ Link[index].Km *= SQR(r) * SQR(r); /* Adjust minor loss factor */ Link[index].Diam = value; /* Update diameter */ - resistance(p, index); /* Update resistance factor */ + resistcoeff(p, index); /* Update resistance coeff. */ } break; @@ -2657,7 +2667,7 @@ int DLLEXPORT EN_setlinkvalue(EN_Project *p, int index, int code, if (value <= 0.0) return (202); Link[index].Len = value / Ucf[ELEV]; - resistance(p, index); + resistcoeff(p, index); } break; @@ -2668,7 +2678,7 @@ int DLLEXPORT EN_setlinkvalue(EN_Project *p, int index, int code, Link[index].Kc = value; if (hyd->Formflag == DW) Link[index].Kc /= (1000.0 * Ucf[ELEV]); - resistance(p, index); + resistcoeff(p, index); } break; @@ -3171,6 +3181,18 @@ int DLLEXPORT EN_setoption(EN_Project *p, int code, EN_API_FLOAT_TYPE v) return (202); hyd->Dmult = value; break; + + case EN_HEADERROR: + if (value < 0.0) + return (202); + hyd->HeadErrorLimit = value / Ucf[HEAD]; + break; + case EN_FLOWCHANGE: + if (value < 0.0) + return (202); + hyd->FlowChangeLimit = value / Ucf[FLOW]; + break; + default: return (251); } @@ -3419,7 +3441,7 @@ int openhydfile(EN_Project *p) out->HydFile = NULL; switch (out->Hydflag) { case SCRATCH: - getTmpName(out->HydFname); + getTmpName(p, out->HydFname); out->HydFile = fopen(out->HydFname, "w+b"); break; case SAVE: @@ -3512,7 +3534,7 @@ int openoutfile(EN_Project *p) // else if ( (OutFile = tmpfile()) == NULL) else { - getTmpName(out->OutFname); + getTmpName(p, out->OutFname); if ((out->OutFile = fopen(out->OutFname, "w+b")) == NULL) { writecon(FMT08); @@ -3530,7 +3552,7 @@ int openoutfile(EN_Project *p) if (!errcode) { if (rep->Tstatflag != SERIES) { // if ( (TmpOutFile = tmpfile()) == NULL) errcode = 304; - getTmpName(out->TmpFname); + getTmpName(p, out->TmpFname); out->TmpOutFile = fopen(out->TmpFname, "w+b"); if (out->TmpOutFile == NULL) errcode = 304; @@ -3837,7 +3859,7 @@ void freedata(EN_Project *p) ---------------------------------------------------------------- */ -char *getTmpName(char *fname) +char *getTmpName(EN_Project *p, char *fname) // // Input: fname = file name string // Output: returns pointer to file name @@ -3847,25 +3869,24 @@ char *getTmpName(char *fname) #ifdef _WIN32 - char* name = NULL; + char* name = NULL; - // --- use Windows _tempnam function to get a pointer to an - // unused file name that begins with "en" - name = _tempnam(NULL, "en"); - if (name == NULL) return NULL; + // --- use Windows _tempnam function to get a pointer to an + // unused file name that begins with "en" + name = _tempnam(NULL, "en"); + if (name == NULL) return NULL; - // --- copy the file name to fname - if (strlen(name) < MAXFNAME) strncpy(fname, name, MAXFNAME); - else fname = NULL; + // --- copy the file name to fname + if (strlen(name) < MAXFNAME) strncpy(fname, name, MAXFNAME); + else fname = NULL; - // --- free the pointer returned by _tempnam - if (name) free(name); + // --- free the pointer returned by _tempnam + if (name) free(name); /* /////////////////// DEPRECATED ///////////////////////////////////// // --- for Windows systems: #ifdef WINDOWS - out_file_t *out = &p->out_files; // --- use system function tmpnam() to create a temporary file name char name[MAXFNAME + 1]; int n; diff --git a/src/epanet2.h b/src/epanet2.h new file mode 100644 index 0000000..81b469e --- /dev/null +++ b/src/epanet2.h @@ -0,0 +1,1236 @@ +/** @file epanet2.h + @see http://github.com/openwateranalytics/epanet + + */ + +/* + ******************************************************************* + + EPANET2.H - Prototypes for EPANET Functions Exported to DLL Toolkit + + VERSION: 2.00 + DATE: 5/8/00 + 10/25/00 + 3/1/01 + 8/15/07 (2.00.11) + 2/14/08 (2.00.12) + AUTHORS: L. Rossman - US EPA - NRMRL + OpenWaterAnalytics members: see git stats for contributors + + ******************************************************************* + */ + + +#ifndef EPANET2_H +#define EPANET2_H + +// the toolkit can be compiled with support for double-precision as well. +// just make sure that you use the correct #define in your client code. +#ifndef EN_API_FLOAT_TYPE + #define EN_API_FLOAT_TYPE float +#endif + +// --- define WINDOWS +#undef WINDOWS +#ifdef _WIN32 + #define WINDOWS +#endif +#ifdef __WIN32__ + #define WINDOWS +#endif + +// --- define DLLEXPORT +#ifndef DLLEXPORT + #ifdef WINDOWS + #ifdef __cplusplus + #define DLLEXPORT extern "C" __declspec(dllexport) + #else + #define DLLEXPORT __declspec(dllexport) __stdcall + #endif // __cplusplus + #elif defined(CYGWIN) + #define DLLEXPORT __stdcall + #elif defined(__APPLE__) + #ifdef __cplusplus + #define DLLEXPORT + #else + #define DLLEXPORT + #endif + #else + #define DLLEXPORT + #endif +#endif + +//#include "epanet_export.h" + +// --- Define the EPANET toolkit constants + +/// Node property codes +typedef enum { + EN_ELEVATION = 0, /**< Node Elevation */ + EN_BASEDEMAND = 1, /**< Node Base Demand, from last demand category */ + EN_PATTERN = 2, /**< Node Demand Pattern */ + EN_EMITTER = 3, /**< Node Emitter Coefficient */ + EN_INITQUAL = 4, /**< Node initial quality */ + EN_SOURCEQUAL = 5, /**< Node source quality */ + EN_SOURCEPAT = 6, /**< Node source pattern index */ + EN_SOURCETYPE = 7, /**< Node source type */ + EN_TANKLEVEL = 8, /**< Tank Level */ + EN_DEMAND = 9, /**< Node current simulated demand */ + EN_HEAD = 10, /**< Node Head value */ + EN_PRESSURE = 11, /**< Node pressure value */ + EN_QUALITY = 12, /**< Node quality value */ + EN_SOURCEMASS = 13, /**< Node source mass value */ + EN_INITVOLUME = 14, /**< Tank or Reservoir initial volume */ + EN_MIXMODEL = 15, /**< Tank mixing model */ + EN_MIXZONEVOL = 16, /**< Tank mixing zone volume */ + EN_TANKDIAM = 17, /**< Tank diameter */ + EN_MINVOLUME = 18, /**< Tank minimum volume */ + EN_VOLCURVE = 19, /**< Tank volume curve index */ + EN_MINLEVEL = 20, /**< Tank minimum level */ + EN_MAXLEVEL = 21, /**< Tank maximum level */ + EN_MIXFRACTION = 22, /**< Tank mixing fraction */ + EN_TANK_KBULK = 23, /**< Tank bulk decay coefficient */ + EN_TANKVOLUME = 24, /**< Tank current volume */ + EN_MAXVOLUME = 25 /**< Tank maximum volume */ +} EN_NodeProperty; + +/// Link property codes +typedef enum { + EN_DIAMETER = 0, + EN_LENGTH = 1, + EN_ROUGHNESS = 2, + EN_MINORLOSS = 3, + EN_INITSTATUS = 4, + EN_INITSETTING = 5, + EN_KBULK = 6, + EN_KWALL = 7, + EN_FLOW = 8, + EN_VELOCITY = 9, + EN_HEADLOSS = 10, + EN_STATUS = 11, + EN_SETTING = 12, + EN_ENERGY = 13, + EN_LINKQUAL = 14, + EN_LINKPATTERN = 15, + EN_EFFICIENCY = 16, + EN_HEADCURVE = 17, + EN_EFFICIENCYCURVE = 18, + EN_PRICEPATTERN = 19 +} EN_LinkProperty; + +/// Time parameter codes +typedef enum { + EN_DURATION = 0, + EN_HYDSTEP = 1, + EN_QUALSTEP = 2, + EN_PATTERNSTEP = 3, + EN_PATTERNSTART = 4, + EN_REPORTSTEP = 5, + EN_REPORTSTART = 6, + EN_RULESTEP = 7, + EN_STATISTIC = 8, + EN_PERIODS = 9, + EN_STARTTIME = 10, + EN_HTIME = 11, + EN_QTIME = 12, + EN_HALTFLAG = 13, + EN_NEXTEVENT = 14, + EN_NEXTEVENTIDX = 15 +} EN_TimeProperty; + + +typedef enum { + EN_ITERATIONS = 0, + EN_RELATIVEERROR = 1 +} EN_AnalysisStatistic; + +typedef enum { + EN_NODECOUNT = 0, /**< Number of Nodes (Juntions + Tanks + Reservoirs) */ + EN_TANKCOUNT = 1, /**< Number of Tanks and Reservoirs */ + EN_LINKCOUNT = 2, /**< Number of Links (Pipes + Pumps + Valves) */ + EN_PATCOUNT = 3, /**< Number of Time Patterns */ + EN_CURVECOUNT = 4, /**< Number of Curves */ + EN_CONTROLCOUNT = 5, /**< Number of Control Statements */ + EN_RULECOUNT = 6 /**< Number of Rule-based Control Statements */ +} EN_CountType; + +typedef enum { + EN_JUNCTION = 0, + EN_RESERVOIR = 1, + EN_TANK = 2 +} EN_NodeType; + +typedef enum { + EN_CVPIPE = 0, /* Link types. */ + EN_PIPE = 1, /* See LinkType in TYPES.H */ + EN_PUMP = 2, + EN_PRV = 3, + EN_PSV = 4, + EN_PBV = 5, + EN_FCV = 6, + EN_TCV = 7, + EN_GPV = 8 +} EN_LinkType; + +typedef enum { + EN_NONE = 0, /* Quality analysis types. */ + EN_CHEM = 1, /* See QualType in TYPES.H */ + EN_AGE = 2, + EN_TRACE = 3 +} EN_QualityType; + +typedef enum { + EN_CONCEN = 0, /* Source quality types. */ + EN_MASS = 1, /* See SourceType in TYPES.H. */ + EN_SETPOINT = 2, + EN_FLOWPACED = 3 +} EN_SourceType; + +typedef enum { /* Head loss formula: */ + EN_HW = 0, /* Hazen-Williams */ + EN_DW = 1, /* Darcy-Weisbach */ + EN_CM = 2 /* Chezy-Manning */ +} EN_FormType; /* See FormType in TYPES.H */ + +typedef enum { + EN_CFS = 0, /* Flow units types. */ + EN_GPM = 1, /* See FlowUnitsType */ + EN_MGD = 2, /* in TYPES.H. */ + EN_IMGD = 3, + EN_AFD = 4, + EN_LPS = 5, + EN_LPM = 6, + EN_MLD = 7, + EN_CMH = 8, + EN_CMD = 9 +} EN_FlowUnits; + + +/// Simulation Option codes +typedef enum { + EN_TRIALS = 0, + EN_ACCURACY = 1, + EN_TOLERANCE = 2, + EN_EMITEXPON = 3, + EN_DEMANDMULT = 4, + EN_HEADERROR = 5, + EN_FLOWCHANGE = 6 +} EN_Option; + +typedef enum { + EN_LOWLEVEL = 0, /* Control types. */ + EN_HILEVEL = 1, /* See ControlType */ + EN_TIMER = 2, /* in TYPES.H. */ + EN_TIMEOFDAY = 3 +} EN_ControlType; + + + +typedef enum { + EN_AVERAGE = 1, /* Time statistic types. */ + EN_MINIMUM = 2, /* See TstatType in TYPES.H */ + EN_MAXIMUM = 3, + EN_RANGE = 4 +} EN_StatisticType; + + + +typedef enum { + EN_MIX1 = 0, /* Tank mixing models */ + EN_MIX2 = 1, + EN_FIFO = 2, + EN_LIFO = 3 +} EN_MixingModel; + + + +typedef enum { + EN_NOSAVE = 0, + EN_SAVE = 1, + EN_INITFLOW = 10, + EN_SAVE_AND_INIT = 11 +} EN_SaveOption; + + + +typedef enum { + EN_CONST_HP = 0, /* constant horsepower */ + EN_POWER_FUNC = 1, /* power function */ + EN_CUSTOM = 2 /* user-defined custom curve */ +} EN_CurveType; + + + +// --- Declare the EPANET toolkit functions +#if defined(__cplusplus) +extern "C" { +#endif + + /** + @brief The EPANET Project wrapper object + */ + typedef struct EN_Project EN_Project; + typedef struct EN_Pattern EN_Pattern; + typedef struct EN_Curve EN_Curve; + + /** + @brief runs a complete EPANET simulation + @param inpFile pointer to name of input file (must exist) + @param rptFile pointer to name of report file (to be created) + @param binOutFile pointer to name of binary output file (to be created) + @param callback a callback function that takes a character string (char *) as its only parameter. + @return error code + + The callback function should reside in and be used by the calling + code to display the progress messages that EPANET generates + as it carries out its computations. If this feature is not + needed then the argument should be NULL. + */ + int DLLEXPORT ENepanet(char *inpFile, char *rptFile, char *binOutFile, void (*callback) (char *)); + + /** + @brief Initializes an EPANET session + @param rptFile pointer to name of report file (to be created) + @param binOutFile pointer to name of binary output file (to be created) + @param UnitsType flow units flag + @param HeadlossFormula headloss formula flag + @return error code + */ + int DLLEXPORT ENinit(char *rptFile, char *binOutFile, int UnitsType, int HeadlossFormula); + + /** + @brief Opens EPANET input file & reads in network data + @param inpFile pointer to name of input file (must exist) + @param rptFile pointer to name of report file (to be created) + @param binOutFile pointer to name of binary output file (to be created) + @return error code + */ + int DLLEXPORT ENopen(char *inpFile, char *rptFile, char *binOutFile); + + /** + @brief Saves current data to "INP" formatted text file. + @param filename The file path to create + @return Error code + */ + int DLLEXPORT ENsaveinpfile(char *filename); + + /** + @brief Frees all memory and files used by EPANET + @return Error code + */ + int DLLEXPORT ENclose(); + + /** + @brief Solves the network hydraulics for all time periods + @return Error code + */ + int DLLEXPORT ENsolveH(); + + /** + @brief Saves hydraulic results to binary file + @return Error code + + Must be called before ENreport() if no WQ simulation has been made. + Should not be called if ENsolveQ() will be used. + */ + int DLLEXPORT ENsaveH(); + + /** + @brief Sets up data structures for hydraulic analysis + @return Error code + */ + int DLLEXPORT ENopenH(); + + /** + @brief Initializes hydraulic analysis + @param initFlag 2-digit flag where 1st (left) digit indicates if link flows should be re-initialized (1) or not (0), and 2nd digit indicates if hydraulic results should be saved to file (1) or not (0). + @return Error code + */ + int DLLEXPORT ENinitH(int initFlag); + + /** + @brief Run a hydraulic solution period + @param[out] currentTime The current simulation time in seconds + @return Error or warning code + @see ENsolveH + + This function is used in a loop with ENnextH() to run + an extended period hydraulic simulation. + See ENsolveH() for an example. + */ + int DLLEXPORT ENrunH(long *currentTime); + + /** + @brief Determine time (in seconds) until next hydraulic event + @param[out] tStep Time (seconds) until next hydraulic event. 0 marks end of simulation period. + @return Error code + + This function is used in a loop with ENrunH() to run an extended period hydraulic simulation. + See ENsolveH() for an example. + */ + int DLLEXPORT ENnextH(long *tStep); + + + /** + @brief Frees data allocated by hydraulics solver + @return Error code + */ + int DLLEXPORT ENcloseH(); + + /** + @brief Copies binary hydraulics file to disk + @param filename Name of file to be created + @return Error code + */ + int DLLEXPORT ENsavehydfile(char *filename); + + /** + @brief Opens previously saved binary hydraulics file + @param filename Name of file to be used + @return Error code + */ + int DLLEXPORT ENusehydfile(char *filename); + + /** + @brief Solves for network water quality in all time periods + @return Error code + */ + int DLLEXPORT ENsolveQ(); + + /** + @brief Sets up data structures for WQ analysis + @return Error code + */ + int DLLEXPORT ENopenQ(); + + /** + @brief Initializes water quality analysis + @param saveFlag EN_SAVE (1) if results saved to file, EN_NOSAVE (0) if not + @return Error code + */ + int DLLEXPORT ENinitQ(int saveFlag); + + /** + @brief Retrieves hydraulic & WQ results at time t. + @param[out] currentTime Current simulation time, in seconds. + @return Error code + + This function is used in a loop with ENnextQ() to run + an extended period WQ simulation. See ENsolveQ() for + an example. + */ + int DLLEXPORT ENrunQ(long *currentTime); + + /** + @brief Advances WQ simulation to next hydraulic event. + @param[out] tStep Time in seconds until next hydraulic event. 0 marks end of simulation period. + @return Error code + + This function is used in a loop with ENrunQ() to run + an extended period WQ simulation. See ENsolveQ() for + an example. + */ + int DLLEXPORT ENnextQ(long *tStep); + + /** + @brief Advances WQ simulation by a single WQ time step + @param[out] timeLeft Time left in overall simulation (in seconds) + @return Error code + + This function is used in a loop with ENrunQ() to run + an extended period WQ simulation. + */ + int DLLEXPORT ENstepQ(long *timeLeft); + + /** + @brief Frees data allocated by water quality solver. + @return Error code. + */ + int DLLEXPORT ENcloseQ(); + + /** + @brief Writes line of text to the report file. + @param line Text string to write + @return Error code. + */ + int DLLEXPORT ENwriteline(char *line); + + /** + @brief Writes simulation report to the report file + @return Error code + */ + int DLLEXPORT ENreport(); + + /** + @brief Resets report options to default values + @return Error code + */ + int DLLEXPORT ENresetreport(); + + /** + @brief Processes a reporting format command + @return Error code + */ + int DLLEXPORT ENsetreport(char *reportFormat); + + /** + @brief Retrieves parameters that define a simple control + @param controlIndex Control index (position of control statement in the input file, starting from 1) + @param[out] controlType Control type code (see EPANET2.H) + @param[out] linkIndex Index of controlled link + @param[out] setting Control setting on link + @param[out] nodeIndex Index of controlling node (0 for TIMER or TIMEOFDAY control) + @param[out] level Control level (tank level, junction pressure, or time (seconds)) + @return Error code + */ + int DLLEXPORT ENgetcontrol(int controlIndex, int *controlType, int *linkIndex, EN_API_FLOAT_TYPE *setting, int *nodeIndex, EN_API_FLOAT_TYPE *level); + + /** + @brief Retrieves the number of components of a given type in the network. + @param code Component code (see EPANET2.H) + @param[out] count Number of components in network + @return Error code + */ + int DLLEXPORT ENgetcount(int code, int *count); + + /** + @brief Gets value for an analysis option + @param code Option code (see EPANET2.H) + @param[out] value Option value + @return Error code + */ + int DLLEXPORT ENgetoption(int code, EN_API_FLOAT_TYPE *value); + + /** + @brief Retrieves value of specific time parameter. + @param code Time parameter code + @param[out] value Value of time parameter. + @return Error code + */ + int DLLEXPORT ENgettimeparam(int code, long *value); + + /** + @brief Retrieves the flow units code + @param[out] code Code of flow units in use + @return Error code + */ + int DLLEXPORT ENgetflowunits(int *code); + + /** + @brief Sets the flow units + @param code Code of flow units to use + @return Error code + */ + int DLLEXPORT ENsetflowunits(int code); + + /** + @brief Retrieves the index of the time pattern with specified ID + @param id String ID of the time pattern + @param[out] index Index of the specified time pattern + @return Error code + @see ENgetpatternid + */ + int DLLEXPORT ENgetpatternindex(char *id, int *index); + + /** + @brief Retrieves ID of a time pattern with specific index. + @param index The index of a time pattern. + @param[out] id The string ID of the time pattern. + @return Error code + @see ENgetpatternindex + */ + int DLLEXPORT ENgetpatternid(int index, char *id); + + /** + @brief Retrieves the number of multipliers in a time pattern. + @param index The index of a time pattern. + @param[out] len The length of the time pattern. + @return Error code + */ + int DLLEXPORT ENgetpatternlen(int index, int *len); + + /** + @brief Retrive a multiplier from a pattern for a specific time period. + @param index The index of a time pattern + @param period The pattern time period. First time period is 1. + @param[out] value Pattern multiplier + @return Error code + */ + int DLLEXPORT ENgetpatternvalue(int index, int period, EN_API_FLOAT_TYPE *value); + + /** + @brief Retrieve the average multiplier value in a time pattern + @param index The index of a time pattern + @param[out] value The average of all of this time pattern's values + @return Error code + */ + int DLLEXPORT ENgetaveragepatternvalue(int index, EN_API_FLOAT_TYPE *value); + + /** + @brief Retrieve the type of quality analytis to be run. + @param[out] qualcode The quality analysis code number. + @param[out] tracenode The index of node being traced, if qualcode == trace + @return Error code + @see ENsetqualtype + */ + int DLLEXPORT ENgetqualtype(int *qualcode, int *tracenode); + + /** + @brief Get the text of an error code. + @param errcode The error code + @param[out] errmsg The error string represented by the code + @param maxLen The maximum number of characters to copy into the char pointer errmsg + @return Error code + */ + int DLLEXPORT ENgeterror(int errcode, char *errmsg, int maxLen); + + /** + @brief Get hydraulic simulation statistic + @param code The type of statistic to get + @param[out] value The value of the statistic + @return Error code + */ + int DLLEXPORT ENgetstatistic(int code, EN_API_FLOAT_TYPE* value); + + /** + @brief Get index of node with specified ID + @param id The string ID of the node + @param[out] index The node's index (first node is index 1) + @return Error code + @see ENgetnodeid + */ + int DLLEXPORT ENgetnodeindex(char *id, int *index); + + /** + @brief Get the string ID of the specified node. + @param index The index of the node (first node is index 1) + @param[out] id The string ID of the specified node. Up to MAXID characters will be copied, so id must be pre-allocated by the calling code to hold at least that many characters. + @return Error code + @see ENgetnodeindex + */ + int DLLEXPORT ENgetnodeid(int index, char *id); + + /** + @brief Get the type of node with specified index. + @param index The index of a node (first node is index 1) + @param[out] code The type code for the node. + @return Error code + */ + int DLLEXPORT ENgetnodetype(int index, int *code); + + /** + @brief Get a property value for specified node + @param index The index of a node (first node is index 1). + @param code The property type code + @param[out] value The value of the node's property. + @return Error code + @see EN_NodeProperty + */ + int DLLEXPORT ENgetnodevalue(int index, int code, EN_API_FLOAT_TYPE *value); + + /** + @brief Get coordinates (x,y) for a node. + @param index The index of a node (first node is index 1). + @param[out] x X-value of node's coordinate + @param[out] y Y-value of node's coordinate + @return Error code + @see ENsetcoord + */ + int DLLEXPORT ENgetcoord(int index, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y); + + /** + @brief Set coordinates (x,y) for a node. + @param index The index of a node (first node is index 1) + @param x X-value of node's coordinate + @param y Y-value of node's coordinate + @return Error code + @see ENgetcoord + */ + int DLLEXPORT ENsetcoord(int index, EN_API_FLOAT_TYPE x, EN_API_FLOAT_TYPE y); + + /** + @brief Get the number of demand categories for a node. + @param nodeIndex The index of a node (first node is index 1) + @param[out] numDemands The number of demand categories + @return Error code + */ + int DLLEXPORT ENgetnumdemands(int nodeIndex, int *numDemands); + + /** + @brief Get a node's base demand for a specified category. + @param nodeIndex The index of a node (first node is index 1) + @param demandIndex The index of the demand category (starting at 1) + @return Error code + */ + int DLLEXPORT ENgetbasedemand(int nodeIndex, int demandIndex, EN_API_FLOAT_TYPE *baseDemand); + + /** + @brief Get the index of the demand pattern assigned to a node for a category index. + @param nodeIndex The index of a node (first node is index 1). + @param demandIndex The index of a category (first category is index 1). + @param[out] pattIndex The index of the pattern for this node and category. + @return Error code + */ + int DLLEXPORT ENgetdemandpattern(int nodeIndex, int demandIndex, int *pattIndex); + + /** + @brief Get the index of a Link with specified ID. + @param id The string ID of a link. + @param[out] index The index of the named link (first link is index 1) + @return Error code + @see ENgetlinkid + */ + int DLLEXPORT ENgetlinkindex(char *id, int *index); + + /** + @brief Get the string ID of a link with specified index + @param index The index of a link (first link is index 1) + @param[out] id The ID of the link. Up to MAXID characters will be copied, so id must be pre-allocated by the calling code to hold at least that many characters. + @return Error code + @see ENgetlinkindex + */ + int DLLEXPORT ENgetlinkid(int index, char *id); + + /** + @brief Get the link type code for a specified link + @param index The index of a link (first link is index 1) + @param[out] code The type code of the link. + @return Error code + @see EN_LinkType + */ + int DLLEXPORT ENgetlinktype(int index, EN_LinkType *code); + + /** + @brief Set the link type code for a specified link + @param id The id of a link + @param type The type code of the link. + @return Error code + @see EN_LinkType + */ + int DLLEXPORT ENsetlinktype(char *id, EN_LinkType type); + + /** + @brief Get the indexes of a link's start- and end-nodes. + @param index The index of a link (first link is index 1) + @param[out] node1 The index of the link's start node (first node is index 1). + @param[out] node2 The index of the link's end node (first node is index 1). + @return Error code + @see ENgetnodeid, ENgetlinkid + */ + int DLLEXPORT ENgetlinknodes(int index, int *node1, int *node2); + + /** + @brief Get a property value for specified link. + @param index The index of a node (first node is index 1). + @param code The parameter desired. + @param[out] value The value of the link's specified property. + @return Error code + @see ENgetnodevalue, EN_LinkProperty + */ + int DLLEXPORT ENgetlinkvalue(int index, int code, EN_API_FLOAT_TYPE *value); + + /** + @brief Get a curve's properties. + @param curveIndex The index of a curve (first curve is index 1). + @param[out] id The curve's string ID. Client code must preallocate at least MAXID characters. + @param[out] nValues The number of values in the curve's (x,y) list. + @param[out] xValues The curve's x-values. Pointer must be freed by client. + @param[out] yValues The curve's y-values. Pointer must be freed by client. + @return Error code. + */ + int DLLEXPORT ENgetcurve(int curveIndex, char* id, int *nValues, EN_API_FLOAT_TYPE **xValues, EN_API_FLOAT_TYPE **yValues); + + /** + @brief Retrieves the curve index for a specified pump index. + @param pumpIndex The index of a pump + @param[out] curveIndex The index of the curve used by the pump. + @return Error code. + */ + int DLLEXPORT ENgetheadcurveindex(int pumpIndex, int *curveIndex); + + /** + @brief Sets the curve id for a specified pump index. + @param pumpIndex The index of the pump + @param curveIndex The index of the curve used by the pump + @return Error code. + */ + int DLLEXPORT ENsetheadcurveindex(int pumpIndex, int curveIndex); + + /** + @brief Get the type of pump + @param linkIndex The index of the pump element + @param[out] outType The integer-typed pump curve type signifier (output parameter) + @return Error code + @see EN_CurveType + */ + int DLLEXPORT ENgetpumptype(int linkIndex, int *outType); + + /** + @brief Get the version number. This number is to be interpreted with implied decimals, i.e., "20100" == "2(.)01(.)00" + @param[out] version The version of EPANET + @return Error code. + */ + int DLLEXPORT ENgetversion(int *version); + + /** + @brief Specify parameters to define a simple control + @param cindex The index of the control to edit. First control is index 1. + @param ctype The type code to set for this control. + @param lindex The index of a link to control. + @param setting The control setting applied to the link. + @param nindex The index of a node used to control the link, or 0 for TIMER / TIMEOFDAY control. + @param level control point (tank level, junction pressure, or time in seconds). + @return Error code. + */ + int DLLEXPORT ENsetcontrol(int cindex, int ctype, int lindex, EN_API_FLOAT_TYPE setting, int nindex, EN_API_FLOAT_TYPE level); + + /** + @brief Set a property value for a node. + @param index The index of a node. First node is index 1. + @param code The code for the proprty to set. + @param v The value to set for this node and property. + @return Error code. + @see EN_NodeProperty + */ + int DLLEXPORT ENsetnodevalue(int index, int code, EN_API_FLOAT_TYPE v); + + /** + @brief Set a proprty value for a link. + @param index The index of a link. First link is index 1. + @param code The code for the property to set. + @param v The value to set for this link and property. + @return Error code. + @see EN_LinkProperty + */ + int DLLEXPORT ENsetlinkvalue(int index, int code, EN_API_FLOAT_TYPE v); + + /** + @brief Add a new time pattern. + @param id The string ID of the pattern to add. + @return Error code. + @see ENgetpatternindex + */ + int DLLEXPORT ENaddpattern(char *id); + + /** + @brief Set multipliers for a specific pattern + @param index The index of a pattern. First pattern is index 1. + @param f An array of multipliers + @param len The length of array f. + @return Error code. + @see ENgetpatternindex + */ + int DLLEXPORT ENsetpattern(int index, EN_API_FLOAT_TYPE *f, int len); + + /** + @brief Set the multiplier for a specific pattern at a specific period. + @param index The index of a pattern. First pattern is index 1. + @param period The period of the pattern to set. + @param value The value of the multiplier to set. + @return Error code. + */ + int DLLEXPORT ENsetpatternvalue(int index, int period, EN_API_FLOAT_TYPE value); + + /** + @brief Set the value for a time parameter. + @param code The code for the parameter to set. + @param value The desired value of the parameter. + @return Error code. + @see EN_TimeProperty + */ + int DLLEXPORT ENsettimeparam(int code, long value); + + /** + @brief Set a value for an anlysis option. + @param code The code for the desired option. + @param v The desired value for the option specified. + @return Error code. + @see EN_Option + */ + int DLLEXPORT ENsetoption(int code, EN_API_FLOAT_TYPE v); + + /** + @brief Sets the level of hydraulic status reporting. + @param code Status reporting code. + @return Error code. + */ + int DLLEXPORT ENsetstatusreport(int code); + + /** + @brief Sets type of quality analysis called for + @param qualcode WQ parameter code, EN_QualityType + @param chemname Name of WQ constituent + @param chemunits Concentration units of WQ constituent + @param tracenode ID of node being traced (if applicable) + @return Error code. + @see EN_QualityType + + chemname and chemunits only apply when WQ analysis is for chemical. tracenode only applies when WQ analysis is source tracing. + */ + int DLLEXPORT ENsetqualtype(int qualcode, char *chemname, char *chemunits, char *tracenode); + + /** + @brief Get quality analysis information (type, chemical name, units, trace node ID) + @param[out] qualcode The EN_QualityType code being used. + @param[out] chemname The name of the WQ constituent. + @param[out] chemunits The cencentration units of the WQ constituent. + @param[out] tracenode The trace node ID. + @return Error code. + @see EN_QualityType + */ + int DLLEXPORT ENgetqualinfo(int *qualcode, char *chemname, char *chemunits, int *tracenode); + + /** + @brief Sets the node's base demand for a category. + @param nodeIndex The index of a node. + @param demandIdx The index of a demand category. + @param baseDemand The base demand multiplier for the selected category. + @return Error code. + @see ENgetbasedemand + */ + int DLLEXPORT ENsetbasedemand(int nodeIndex, int demandIdx, EN_API_FLOAT_TYPE baseDemand); + + /** + @brief Retrieves index of curve with specific ID. + @param id The ID of a curve. + @param[out] index The index of the named curve + @return Error code. + @see ENgetcurveid + */ + int DLLEXPORT ENgetcurveindex(char *id, int *index); + + /** + @brief Retrieves ID of a curve with specific index. + @param index The index of a curve. + @param[out] id The ID of the curve specified. + @return Error code. + @see ENsetcurveindex + + NOTE: 'id' must be able to hold MAXID characters + */ + int DLLEXPORT ENgetcurveid(int index, char *id); + + /** + @brief Retrieves number of points in a curve + @param index The index of a curve. + @param[out] len The length of the curve coordinate list + @return Error code. + @see ENgetcurvevalue + */ + int DLLEXPORT ENgetcurvelen(int index, int *len); + + /** + @brief retrieves x,y point for a specific point number and curve + @param curveIndex The index of a curve + @param pointIndex The index of a point in the curve + @param[out] x The x-value of the specified point. + @param[out] y The y-value of the specified point. + @return Error code. + @see ENgetcurvelen ENsetcurvevalue + */ + int DLLEXPORT ENgetcurvevalue(int curveIndex, int pointIndex, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y); + + /** + @brief Sets x,y point for a specific point and curve. + @param curveIndex The index of a curve. + @param pointIndex The index of a point in the curve. + @param x The x-value of the point. + @param y The y-value of the point. + @return Error code. + */ + int DLLEXPORT ENsetcurvevalue(int curveIndex, int pointIndex, EN_API_FLOAT_TYPE x, EN_API_FLOAT_TYPE y); + + /** + @brief Sets x,y values for a specified curve. + @param index The index of a curve. + @param x An array of x-values for the curve. + @param y An array of y-values for the curve. + @param len The length of the arrays x and y. + @return Error code. + */ + int DLLEXPORT ENsetcurve(int index, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y, int len); + + /** + @brief Adds a new curve appended to the end of the existing curves. + @param id The name of the curve to be added. + @return Error code. + @see ENgetcurveindex ENsetcurve + */ + int DLLEXPORT ENaddcurve(char *id); + + + /** + @brief Gets the number of premises, true actions, and false actions and the priority of an existing rule-based control. + @param index The index of a rule-based control. + @param nPremises The number of conditions in a rule-based control. + @param nTrueActions The number of actions that are executed when the conditions in the rule-based control are met. + @param nFalseActions The number of actions that are executed when the conditions in the rule-based control are not met. + @param priority The priority of a rule-based control. + @return Error code. + */ + int DLLEXPORT ENgetrule(int index, int *nPremises, int *nTrueActions, int *nFalseActions, EN_API_FLOAT_TYPE *priority); + + /** + @brief Sets the priority of the existing rule-based control. + @param index The index of a rule-based control. + @param priority The priority to be set in the rule-based control. + @return Error code. + */ + int DLLEXPORT ENsetrulepriority(int index, EN_API_FLOAT_TYPE priority); + + /** + @brief Gets the components of a premise/condition in an existing rule-based control. + @param indexRule The index of a rule-based control. + @param indexPremise The index of the premise. + @param logop The logiv operator (IF/AND/OR) in the premise + @param object The object (e.g. TANK) the premise is looking at. + @param indexObj The index of the object (e.g. the index of the tank). + @param variable The variable to be checked (e.g. level). + @param relop The relashionship operator (e.g. LARGER THAN) in the premise. + @param status The status of the object to be checked (e.g. CLOSED) + @param value The value of the variable to be checked (e.g. 5.5) + @return Error code. + */ + int DLLEXPORT ENgetpremise(int indexRule, int indexPremise, int *logop, int *object, int *indexObj, int *variable, int *relop, int *status, EN_API_FLOAT_TYPE *value); + + /** + @brief Sets the components of a premise/condition in an existing rule-based control. + @param indexRule The index of a rule-based control. + @param indexPremise The index of the premise. + @param logop The logiv operator (IF/AND/OR) in the premise + @param object The object (e.g. TANK) the premise is looking at. + @param indexObj The index of the object (e.g. the index of the tank). + @param variable The variable to be checked (e.g. level). + @param relop The relashionship operator (e.g. LARGER THAN) in the premise. + @param status The status of the object to be checked (e.g. CLOSED) + @param value The value of the variable to be checked (e.g. 5.5) + @return Error code. + */ + int DLLEXPORT ENsetpremise(int indexRule, int indexPremise, int logop, int object, int indexObj, int variable, int relop, int status, EN_API_FLOAT_TYPE value); + + /** + @brief Sets the index of an object in a premise of an existing rule-based control. + @param indexRule The index of a rule-based control. + @param indexPremise The index of the premise. + @param indexObj The index of the object (e.g. the index of the tank). + @return Error code. + */ + int DLLEXPORT ENsetpremiseindex(int indexRule, int indexPremise, int indexObj); + + /** + @brief Sets the status in a premise of an existing rule-based control. + @param indexRule The index of a rule-based control. + @param indexPremise The index of the premise. + @param status The status of the object to be checked (e.g. CLOSED) + @return Error code. + */ + int DLLEXPORT ENsetpremisestatus(int indexRule, int indexPremise, int status); + + /** + @brief Sets the value in a premise of an existing rule-based control. + @param indexRule The index of a rule-based control. + @param indexPremise The index of the premise. + @param value The value of the variable to be checked (e.g. 5.5) + @return Error code. + */ + int DLLEXPORT ENsetpremisevalue(int indexRule, int indexPremise, EN_API_FLOAT_TYPE value); + + /** + @brief Gets the components of a true-action in an existing rule-based control. + @param indexRule The index of a rule-based control. + @param indexAction The index of the action when the conditions in the rule are met. + @param indexLink The index of the link in the action (e.g. index of Pump 1) + @param status The status of the link (e.g. CLOSED) + @param setting The value of the link (e.g. pump speed 0.9) + @return Error code. + */ + int DLLEXPORT ENgettrueaction(int indexRule, int indexAction, int *indexLink, int *status, EN_API_FLOAT_TYPE *setting); + + /** + @brief Sets the components of a true-action in an existing rule-based control. + @param indexRule The index of a rule-based control. + @param indexAction The index of the action when the conditions in the rule are met. + @param indexLink The index of the link in the action (e.g. index of Pump 1) + @param status The status of the link (e.g. CLOSED) + @param setting The value of the link (e.g. pump speed 0.9) + @return Error code. + */ + int DLLEXPORT ENsettrueaction(int indexRule, int indexAction, int indexLink, int status, EN_API_FLOAT_TYPE setting); + + /** + @brief Gets the components of a false-action in an existing rule-based control. + @param indexRule The index of a rule-based control. + @param indexAction The index of the action when the conditions in the rule are not met. + @param indexLink The index of the link in the action (e.g. index of Pump 1) + @param status The status of the link (e.g. CLOSED) + @param setting The value of the link (e.g. pump speed 0.9) + @return Error code. + */ + int DLLEXPORT ENgetfalseaction(int indexRule, int indexAction, int *indexLink, int *status, EN_API_FLOAT_TYPE *setting); + + /** + @brief Sets the components of a false-action in an existing rule-based control. + @param indexRule The index of a rule-based control. + @param indexAction The index of the action when the conditions in the rule are not met. + @param indexLink The index of the link in the action (e.g. index of Pump 1) + @param status The status of the link (e.g. CLOSED) + @param setting The value of the link (e.g. pump speed 0.9) + @return Error code. + */ + int DLLEXPORT ENsetfalseaction(int indexRule, int indexAction, int indexLink, int status, EN_API_FLOAT_TYPE setting); + + /** + @brief Returns the ID of a rule. + @param indexRule The index of a rule-based control. + @param id The ID of the rule + @return Error code. + */ + int DLLEXPORT ENgetruleID(int indexRule, char* id); + + /** + @brief Adds a new node + @param id The name of the node to be added. + @param nodeType The node type code + @return Error code. + */ + int DLLEXPORT ENaddnode(char *id, EN_NodeType nodeType); + + /** + @brief Adds a new link + @param id The name of the link to be added. + @param linkType The link type code + @param fromNode The id of the from node + @param toNode The id of the to node + @return Error code. + */ + int DLLEXPORT ENaddlink(char *id, EN_LinkType linkType, char *fromNode, char *toNode); + + /** + @brief Deletes a node + @param nodeIndex The node index + @return Error code. + */ + int DLLEXPORT ENdeletenode(int nodeIndex); + + /** + @brief Deletes a link + @param linkIndex The link index + @return Error code. + */ + int DLLEXPORT ENdeletelink(int linkIndex); + + + + + /*************************************************** + + Threadsafe versions of all epanet functions + + ***************************************************/ + int DLLEXPORT EN_alloc(EN_Project **p); + int DLLEXPORT EN_free(EN_Project *p); + int DLLEXPORT EN_epanet(char *inpFile, char *rptFile, char *binOutFile, void (*callback) (char *)); + int DLLEXPORT EN_init(EN_Project *p, char *rptFile, char *binOutFile, EN_FlowUnits UnitsType, EN_FormType HeadlossFormula); + int DLLEXPORT EN_open(EN_Project *p, char *inpFile, char *rptFile, char *binOutFile); + int DLLEXPORT EN_saveinpfile(EN_Project *p, char *filename); + int DLLEXPORT EN_close(EN_Project *p); + int DLLEXPORT EN_solveH(EN_Project *p); + int DLLEXPORT EN_saveH(EN_Project *p); + int DLLEXPORT EN_openH(EN_Project *p); + int DLLEXPORT EN_initH(EN_Project *p, int EN_SaveOption); + int DLLEXPORT EN_runH(EN_Project *p, long *currentTime); + int DLLEXPORT EN_nextH(EN_Project *p, long *tStep); + int DLLEXPORT EN_closeH(EN_Project *p); + int DLLEXPORT EN_savehydfile(EN_Project *p, char *filename); + int DLLEXPORT EN_usehydfile(EN_Project *p, char *filename); + int DLLEXPORT EN_solveQ(EN_Project *p); + int DLLEXPORT EN_openQ(EN_Project *p); + int DLLEXPORT EN_initQ(EN_Project *p, int saveFlag); + int DLLEXPORT EN_runQ(EN_Project *p, long *currentTime); + int DLLEXPORT EN_nextQ(EN_Project *p, long *tStep); + int DLLEXPORT EN_stepQ(EN_Project *p, long *timeLeft); + int DLLEXPORT EN_closeQ(EN_Project *p); + int DLLEXPORT EN_writeline(EN_Project *p, char *line); + int DLLEXPORT EN_report(EN_Project *p); + int DLLEXPORT EN_resetreport(EN_Project *p); + int DLLEXPORT EN_setreport(EN_Project *p, char *reportFormat); + int DLLEXPORT EN_getcontrol(EN_Project *p, int controlIndex, int *controlType, int *linkIndex, EN_API_FLOAT_TYPE *setting, int *nodeIndex, EN_API_FLOAT_TYPE *level); + int DLLEXPORT EN_getcount(EN_Project *p, EN_CountType code, int *count); + int DLLEXPORT EN_getoption(EN_Project *p, EN_Option opt, EN_API_FLOAT_TYPE *value); + int DLLEXPORT EN_gettimeparam(EN_Project *p, int code, long *value); + int DLLEXPORT EN_getflowunits(EN_Project *p, int *code); + int DLLEXPORT EN_setflowunits(EN_Project *p, int code); + int DLLEXPORT EN_getpatternindex(EN_Project *p, char *id, int *index); + int DLLEXPORT EN_getpatternid(EN_Project *p, int index, char *id); + int DLLEXPORT EN_getpatternlen(EN_Project *p, int index, int *len); + int DLLEXPORT EN_getpatternvalue(EN_Project *p, int index, int period, EN_API_FLOAT_TYPE *value); + int DLLEXPORT EN_getaveragepatternvalue(EN_Project *p, int index, EN_API_FLOAT_TYPE *value); + int DLLEXPORT EN_getqualtype(EN_Project *p, int *qualcode, int *tracenode); + int DLLEXPORT EN_geterror(int errcode, char *errmsg, int maxLen); + int DLLEXPORT EN_getstatistic(EN_Project *p, int code, EN_API_FLOAT_TYPE* value); + int DLLEXPORT EN_getnodeindex(EN_Project *p, char *id, int *index); + int DLLEXPORT EN_getnodeid(EN_Project *p, int index, char *id); + int DLLEXPORT EN_getnodetype(EN_Project *p, int index, int *code); + int DLLEXPORT EN_getnodevalue(EN_Project *p, int index, int code, EN_API_FLOAT_TYPE *value); + int DLLEXPORT EN_getcoord(EN_Project *p, int index, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y); + int DLLEXPORT EN_setcoord(EN_Project *p, int index, EN_API_FLOAT_TYPE x, EN_API_FLOAT_TYPE y); + int DLLEXPORT EN_getnumdemands(EN_Project *p, int nodeIndex, int *numDemands); + int DLLEXPORT EN_getbasedemand(EN_Project *p, int nodeIndex, int demandIndex, EN_API_FLOAT_TYPE *baseDemand); + int DLLEXPORT EN_getdemandpattern(EN_Project *p, int nodeIndex, int demandIndex, int *pattIndex); + int DLLEXPORT EN_getlinkindex(EN_Project *p, char *id, int *index); + int DLLEXPORT EN_getlinkid(EN_Project *p, int index, char *id); + int DLLEXPORT EN_getlinktype(EN_Project *p, int index, EN_LinkType *code); + int DLLEXPORT EN_setlinktype(EN_Project *p, char *id, EN_LinkType type); + int DLLEXPORT EN_getlinknodes(EN_Project *p, int index, int *node1, int *node2); + int DLLEXPORT EN_getlinkvalue(EN_Project *p, int index, EN_LinkProperty code, EN_API_FLOAT_TYPE *value); + int DLLEXPORT EN_getcurve(EN_Project *p, int curveIndex, char* id, int *nValues, EN_API_FLOAT_TYPE **xValues, EN_API_FLOAT_TYPE **yValues); + int DLLEXPORT EN_getheadcurveindex(EN_Project *p, int pumpIndex, int *curveIndex); + int DLLEXPORT EN_setheadcurveindex(EN_Project *p, int pumpIndex, int curveIndex); + int DLLEXPORT EN_getpumptype(EN_Project *p, int linkIndex, int *outType); + int DLLEXPORT EN_getversion(int *version); + int DLLEXPORT EN_setcontrol(EN_Project *p, int cindex, int ctype, int lindex, EN_API_FLOAT_TYPE setting, int nindex, EN_API_FLOAT_TYPE level); + int DLLEXPORT EN_setnodevalue(EN_Project *p, int index, int code, EN_API_FLOAT_TYPE v); + int DLLEXPORT EN_setlinkvalue(EN_Project *p, int index, int code, EN_API_FLOAT_TYPE v); + int DLLEXPORT EN_addpattern(EN_Project *p, char *id); + int DLLEXPORT EN_setpattern(EN_Project *p, int index, EN_API_FLOAT_TYPE *f, int len); + int DLLEXPORT EN_setpatternvalue(EN_Project *p, int index, int period, EN_API_FLOAT_TYPE value); + int DLLEXPORT EN_settimeparam(EN_Project *p, int code, long value); + int DLLEXPORT EN_setoption(EN_Project *p, int code, EN_API_FLOAT_TYPE v); + int DLLEXPORT EN_setstatusreport(EN_Project *p, int code); + int DLLEXPORT EN_setqualtype(EN_Project *p, int qualcode, char *chemname, char *chemunits, char *tracenode); + int DLLEXPORT EN_getqualinfo(EN_Project *p, int *qualcode, char *chemname, char *chemunits, int *tracenode); + int DLLEXPORT EN_setbasedemand(EN_Project *p, int nodeIndex, int demandIdx, EN_API_FLOAT_TYPE baseDemand); + int DLLEXPORT EN_getcurveindex(EN_Project *p, char *id, int *index); + int DLLEXPORT EN_getcurveid(EN_Project *p, int index, char *id); + int DLLEXPORT EN_getcurvelen(EN_Project *p, int index, int *len); + int DLLEXPORT EN_getcurvevalue(EN_Project *p, int curveIndex, int pointIndex, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y); + int DLLEXPORT EN_setcurvevalue(EN_Project *p, int curveIndex, int pointIndex, EN_API_FLOAT_TYPE x, EN_API_FLOAT_TYPE y); + int DLLEXPORT EN_setcurve(EN_Project *p, int index, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y, int len); + int DLLEXPORT EN_addcurve(EN_Project *p, char *id); + int DLLEXPORT EN_getrule(EN_Project *p, int index, int *nPremises, int *nTrueActions, int *nFalseActions, EN_API_FLOAT_TYPE *priority); + int DLLEXPORT EN_setrulepriority(EN_Project *p, int index, EN_API_FLOAT_TYPE priority); + int DLLEXPORT EN_getpremise(EN_Project *p, int indexRule, int indexPremise, int *logop, int *object, int *indexObj, int *variable, int *relop, int *status, EN_API_FLOAT_TYPE *value); + int DLLEXPORT EN_setpremise(EN_Project *p, int indexRule, int indexPremise, int logop, int object, int indexObj, int variable, int relop, int status, EN_API_FLOAT_TYPE value); + int DLLEXPORT EN_setpremiseindex(EN_Project *p, int indexRule, int indexPremise, int indexObj); + int DLLEXPORT EN_setpremisestatus(EN_Project *p, int indexRule, int indexPremise, int status); + int DLLEXPORT EN_setpremisevalue(EN_Project *p, int indexRule, int indexPremise, EN_API_FLOAT_TYPE value); + int DLLEXPORT EN_gettrueaction(EN_Project *p, int indexRule, int indexAction, int *indexLink, int *status, EN_API_FLOAT_TYPE *setting); + int DLLEXPORT EN_settrueaction(EN_Project *p, int indexRule, int indexAction, int indexLink, int status, EN_API_FLOAT_TYPE setting); + int DLLEXPORT EN_getfalseaction(EN_Project *p, int indexRule, int indexAction, int *indexLink, int *status, EN_API_FLOAT_TYPE *setting); + int DLLEXPORT EN_setfalseaction(EN_Project *p, int indexRule, int indexAction, int indexLink, int status, EN_API_FLOAT_TYPE setting); + int DLLEXPORT EN_getruleID(EN_Project *p, int indexRule, char* id); + int DLLEXPORT EN_addnode(EN_Project *p, char *id, EN_NodeType nodeType); + int DLLEXPORT EN_addlink(EN_Project *p, char *id, EN_LinkType linkType, char *fromNode, char *toNode); + int DLLEXPORT EN_deletenode(EN_Project *p, int nodeIndex); + int DLLEXPORT EN_deletelink(EN_Project *p, int linkIndex); + + + + + + +#if defined(__cplusplus) +} +#endif + +#endif //EPANET2_H diff --git a/src/funcs.h b/src/funcs.h index 79490c6..b65f194 100755 --- a/src/funcs.h +++ b/src/funcs.h @@ -31,27 +31,28 @@ AUTHOR: L. Rossman #include "types.h" -void initpointers(EN_Project *p); /* Initializes pointers */ -int allocdata(EN_Project *p); /* Allocates memory */ -void freeTmplist(STmplist *); /* Frees items in linked list */ -void freeFloatlist(SFloatlist *); /* Frees list of floats */ -void freedata(EN_Project *p); /* Frees allocated memory */ -int openfiles(EN_Project *p, char *,char *,char *); /* Opens input & report files */ -int openhydfile(EN_Project *p); /* Opens hydraulics file */ -int openoutfile(EN_Project *p); /* Opens binary output file */ -int strcomp(char *, char *); /* Compares two strings */ -char* getTmpName(char* fname); /* Gets temporary file name */ -double interp(int, double *,double *, double); /* Interpolates a data curve */ +void initpointers(EN_Project *pr); /* Initializes pointers */ +int allocdata(EN_Project *pr); /* Allocates memory */ +void freeTmplist(STmplist *); /* Frees items in linked list */ +void freeFloatlist(SFloatlist *); /* Frees list of floats */ +void freedata(EN_Project *pr); /* Frees allocated memory */ +int openfiles(EN_Project *pr, char *,char *,char *); /* Opens input & report files */ +int openhydfile(EN_Project *pr); /* Opens hydraulics file */ +int openoutfile(EN_Project *pr); /* Opens binary output file */ +int strcomp(char *, char *); /* Compares two strings */ +char* getTmpName(EN_Project *p, char* fname); /* Gets temporary file name */ +double interp(int n, double x[], double y[], + double xx); /* Interpolates a data curve */ -int findnode(EN_Network *n, char *); /* Finds node's index from ID */ -int findlink(EN_Network *n, char *); /* Finds link's index from ID */ -int findtank(EN_Network *n, int); /* Find tank index from node index */ // (AH) -int findvalve(EN_Network *n, int); /* Find valve index from node index */ // (AH) -int findpump(EN_Network *n, int); /* Find pump index from node index */ // (AH) -char *geterrmsg(int errcode, char *msg); /* Gets text of error message */ -void errmsg(EN_Project *p, int); /* Reports program error */ -void writecon(char *); /* Writes text to console */ -void writewin(void (*vp)(char *), char *); /* Passes text to calling app */ +int findnode(EN_Network *n, char *); /* Finds node's index from ID */ +int findlink(EN_Network *n, char *); /* Finds link's index from ID */ +int findtank(EN_Network *n, int); /* Find tank index from node index */ // (AH) +int findvalve(EN_Network *n, int); /* Find valve index from node index */ // (AH) +int findpump(EN_Network *n, int); /* Find pump index from node index */ // (AH) +char *geterrmsg(int errcode, char *msg); /* Gets text of error message */ +void errmsg(EN_Project *p, int); /* Reports program error */ +void writecon(char *); /* Writes text to console */ +void writewin(void (*vp)(char *), char *); /* Passes text to calling app */ /* ------- INPUT1.C --------------------*/ int getdata(EN_Project *pr); /* Gets network data */ @@ -77,8 +78,9 @@ int getpatterns(EN_Project *pr); /* Gets pattern data from li int getcurves(EN_Project *pr); /* Gets curve data from list */ int findmatch(char *, char *[]); /* Finds keyword in line */ int match(const char *, const char *); /* Checks for word match */ -int gettokens(char *s, char** Tok, int maxToks, char *comment); /* Tokenizes input line */ -int getfloat(char *, double *); /* Converts string to double */ +int gettokens(char *s, char** Tok, int maxToks, + char *comment); /* Tokenizes input line */ +int getfloat(char *, double *); /* Converts string to double */ double hour(char *, char *); /* Converts time to hours */ int setreport(EN_Project *pr, char *); /* Processes reporting command*/ void inperrmsg(EN_Project *pr, int,int,char *); /* Input error message */ @@ -111,7 +113,8 @@ int powercurve(double, double, double, /* Coeffs. of power pump cur double, double, double *, double *, double *); int valvecheck(EN_Project *pr, int, int, int); /* Checks valve placement */ -void changestatus(EN_Network *net, int, StatType, double); /* Changes status of a link */ +void changestatus(EN_Network *net, int, StatType, + double); /* Changes status of a link */ /* -------------- RULES.C --------------*/ void initrules(rules_t *rules); /* Initializes rule base */ @@ -136,13 +139,14 @@ void writerelerr(EN_Project *pr, int, double); /* Writes convergence error void writestatchange(EN_Project *pr, int,char,char); /* Writes link status change */ void writecontrolaction(EN_Project *pr, int, int); /* Writes control action taken*/ void writeruleaction(EN_Project *pr, int, char *); /* Writes rule action taken */ -int writehydwarn(EN_Project *pr, int,double); /* Writes hydraulic warnings */ +int writehydwarn(EN_Project *pr, int,double); /* Writes hydraulic warnings */ void writehyderr(EN_Project *pr, int); /* Writes hydraulic error msg.*/ int disconnected(EN_Project *pr); /* Checks for disconnections */ void marknodes(EN_Project *pr, int, int *, char *); /* Identifies connected nodes */ void getclosedlink(EN_Project *pr, int, char *); /* Finds a disconnecting link */ void writelimits(EN_Project *pr, int,int); /* Writes reporting limits */ -int checklimits(report_options_t *rep, double *,int,int); /* Checks variable limits */ +int checklimits(report_options_t *rep, double *, + int,int); /* Checks variable limits */ void writetime(EN_Project *pr, char *); /* Writes current clock time */ char *clocktime(char *, long); /* Converts time to hrs:min */ char *fillstr(char *, char, int); /* Fills string with character*/ @@ -150,10 +154,7 @@ int getnodetype(EN_Network *net, int); /* Determines node typ /* --------- HYDRAUL.C -----------------*/ int openhyd(EN_Project *pr); /* Opens hydraulics solver */ - -/*** Updated 3/1/01 ***/ void inithyd(EN_Project *pr, int initFlags); /* Re-sets initial conditions */ - int runhyd(EN_Project *pr, long *); /* Solves 1-period hydraulics */ int nexthyd(EN_Project *pr, long *); /* Moves to next time period */ void closehyd(EN_Project *pr); /* Closes hydraulics solver */ @@ -164,64 +165,31 @@ void initlinkflow(EN_Project *pr, int, char, void setlinkflow(EN_Project *pr, int, double); /* Sets link flow via headloss*/ void setlinkstatus(EN_Project *pr, int, char, StatType *, double *); /* Sets link status */ - void setlinksetting(EN_Project *pr, int, double, StatType *, double *); /* Sets pump/valve setting */ -void resistance(EN_Project *pr, int); /* Computes resistance coeff. */ void demands(EN_Project *pr); /* Computes current demands */ int controls(EN_Project *pr); /* Controls link settings */ long timestep(EN_Project *pr); /* Computes new time step */ int tanktimestep(EN_Project *pr, long *); /* Time till tanks fill/drain */ void controltimestep(EN_Project *pr, long *); /* Time till control action */ void ruletimestep(EN_Project *pr, long *); /* Time till rule action */ + void addenergy(EN_Project *pr, long); /* Accumulates energy usage */ void getenergy(EN_Project *pr, int, double *, double *); /* Computes link energy use */ + void tanklevels(EN_Project *pr, long); /* Computes new tank levels */ double tankvolume(EN_Project *pr, int,double); /* Finds tank vol. from grade */ double tankgrade(EN_Project *pr, int,double); /* Finds tank grade from vol. */ -int netsolve(EN_Project *pr, int *,double *); /* Solves network equations */ -int badvalve(EN_Project *pr, int); /* Checks for bad valve */ -int valvestatus(EN_Project *pr); /* Updates valve status */ -int linkstatus(EN_Project *pr); /* Updates link status */ -StatType cvstatus(EN_Project *pr, StatType, - double,double); /* Updates CV status */ -StatType pumpstatus(EN_Project *pr, int,double); /* Updates pump status */ -StatType prvstatus(EN_Project *pr, int,StatType, - double,double,double); /* Updates PRV status */ - -StatType psvstatus(EN_Project *pr, int,StatType, - double,double,double); /* Updates PSV status */ - -StatType fcvstatus(EN_Project *pr, int,StatType, - double,double); /* Updates FCV status */ - -void tankstatus(EN_Project *pr, int,int,int); /* Checks if tank full/empty */ -int pswitch(EN_Project *pr); /* Pressure switch controls */ -double newflows(EN_Project *pr); /* Updates link flows */ -void newcoeffs(EN_Project *pr); /* Computes matrix coeffs. */ -void linkcoeffs(EN_Project *pr); /* Computes link coeffs. */ -void nodecoeffs(EN_Project *pr); /* Computes node coeffs. */ -void valvecoeffs(EN_Project *pr); /* Computes valve coeffs. */ -void pipecoeff(EN_Project *pr, int); /* Computes pipe coeff. */ -double DWcoeff(EN_Project *pr, int, double *); /* Computes D-W coeff. */ -void pumpcoeff(EN_Project *pr, int); /* Computes pump coeff. */ -/*** Updated 10/25/00 ***/ -/*** Updated 12/29/00 ***/ -void curvecoeff(EN_Project *pr, int,double, - double *,double *); /* Computes curve coeffs. */ - +/* ----------- HYDSOLVER.C - ----------*/ +int hydsolve(EN_Project *pr, int *,double *); /* Solves network equations */ -void gpvcoeff(EN_Project *pr, int iLink); /* Computes GPV coeff. */ -void pbvcoeff(EN_Project *pr, int iLink); /* Computes PBV coeff. */ -void tcvcoeff(EN_Project *pr, int iLink); /* Computes TCV coeff. */ -void prvcoeff(EN_Project *pr, int iLink, int n1, int n2); /* Computes PRV coeff. */ -void psvcoeff(EN_Project *pr, int iLink, int n1, int n2); /* Computes PSV coeff. */ -void fcvcoeff(EN_Project *pr, int iLink, int n1, int n2); /* Computes FCV coeff. */ -void emittercoeffs(EN_Project *pr); /* Computes emitter coeffs. */ -double emitflowchange(EN_Project *pr, int); /* Computes new emitter flow */ +/* ----------- HYDCOEFFS.C --------------*/ +void resistcoeff(EN_Project *pr, int k); /* Finds pipe flow resistance */ +void hlosscoeff(EN_Project *pr, int k); /* Finds link head loss coeff */ +void matrixcoeffs(EN_Project *pr); /* Finds hyd. matrix coeffs. */ /* ----------- SMATRIX.C ---------------*/ int createsparse(EN_Project *pr); /* Creates sparse matrix */ @@ -242,7 +210,7 @@ int storesparse(EN_Project *pr, int); /* Stores sparse matrix int ordersparse(hydraulics_t *h, int); /* Orders matrix storage */ void transpose(int,int *,int *, /* Transposes sparse matrix */ int *,int *,int *,int *,int *); -int linsolve(solver_t *s, int); /* Solves set of linear eqns. */ +int linsolve(solver_t *s, int); /* Solves set of inear eqns. */ /* ----------- QUALITY.C ---------------*/ int openqual(EN_Project *pr); /* Opens WQ solver system */ @@ -278,11 +246,10 @@ double pipereact(EN_Project *pr, int,double, double tankreact(EN_Project *pr, double,double, double,long); /* Reacts water in a tank */ double bulkrate(EN_Project *pr, double,double, - double); /* Finds bulk reaction rate */ + double); /* Finds bulk reaction rate */ double wallrate(EN_Project *pr, double,double, double,double); /* Finds wall reaction rate */ - /* ------------ OUTPUT.C ---------------*/ int savenetdata(EN_Project *pr); /* Saves basic data to file */ int savehyd(EN_Project *pr, long *); /* Saves hydraulic solution */ diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c new file mode 100644 index 0000000..a1bbca5 --- /dev/null +++ b/src/hydcoeffs.c @@ -0,0 +1,919 @@ +/* +********************************************************************* + +HYDCOEFFS.C -- hydraulic coefficients for the EPANET Program + +******************************************************************* +*/ + +#include +#include +#ifndef __APPLE__ +#include +#else +#include +#endif +#include +#include "types.h" +#include "funcs.h" + +// Constants used for computing Darcy-Weisbach friction factor +const double A1 = 0.314159265359e04; // 1000*PI +const double A2 = 0.157079632679e04; // 500*PI +const double A3 = 0.502654824574e02; // 16*PI +const double A4 = 6.283185307; // 2*PI +const double A8 = 4.61841319859; // 5.74*(PI/4)^.9 +const double A9 = -8.685889638e-01; // -2/ln(10) +const double AA = -1.5634601348; // -2*.9*2/ln(10) +const double AB = 3.28895476345e-03; // 5.74/(4000^.9) +const double AC = -5.14214965799e-03; // AA*AB + +// External functions +//void resistcoeff(EN_Project *pr, int k); +//void hlosscoeff(EN_Project *pr, int k); +//void matrixcoeffs(EN_Project *pr); + +// Local functions +static void linkcoeffs(EN_Project *pr); +static void nodecoeffs(EN_Project *pr); +static void valvecoeffs(EN_Project *pr); +static void emittercoeffs(EN_Project *pr); + +static void pipecoeff(EN_Project *pr, int k); +static void DWpipecoeff(EN_Project *pr, int k); +static double frictionFactor(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 valvecoeff(EN_Project *pr, int k); +static void gpvcoeff(EN_Project *pr, int k); +static void pbvcoeff(EN_Project *pr, int k); +static void tcvcoeff(EN_Project *pr, int k); +static void prvcoeff(EN_Project *pr, int k, int n1, int n2); +static void psvcoeff(EN_Project *pr, int k, int n1, int n2); +static void fcvcoeff(EN_Project *pr, int k, int n1, int n2); + + + +void resistcoeff(EN_Project *pr, int k) +/* +**-------------------------------------------------------------------- +** Input: k = link index +** Output: none +** Purpose: computes link flow resistance coefficient +**-------------------------------------------------------------------- +*/ +{ + double e, d, L; + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + Slink *link = &net->Link[k]; + + switch (link->Type) { + + // ... Link is a pipe. Compute resistance based on headloss formula. + // Friction factor for D-W formula gets included during head loss + // calculation. + case EN_CVPIPE: + case EN_PIPE: + e = link->Kc; // Roughness coeff. + d = link->Diam; // Diameter + L = link->Len; // Length + switch (hyd->Formflag) + { + case HW: + link->R = 4.727*L / pow(e, hyd->Hexp) / pow(d, 4.871); + break; + case DW: + link->R = L / 2.0 / 32.2 / d / SQR(PI*SQR(d) / 4.0); + break; + case CM: + link->R = SQR(4.0*e / (1.49*PI*d*d)) * pow((d / 4.0), -1.333)*L; + } + break; + + // ... Link is a pump. Use huge resistance. + case EN_PUMP: + link->R = CBIG; + break; + + // ... For all other links (e.g. valves) use a small resistance + default: + link->R = CSMALL; + break; + } +} + + +void hlosscoeff(EN_Project *pr, int k) +/* +**-------------------------------------------------------------- +** Input: k = link index +** Output: none +** Purpose: computes P and Y coefficients for a link +**-------------------------------------------------------------- +*/ +{ + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + Slink *link = &net->Link[k]; + + switch (link->Type) + { + case EN_CVPIPE: + case EN_PIPE: + pipecoeff(pr, k); + break; + case EN_PUMP: + pumpcoeff(pr, k); + break; + case EN_PBV: + pbvcoeff(pr, k); + break; + case EN_TCV: + tcvcoeff(pr, k); + break; + case EN_GPV: + gpvcoeff(pr, k); + break; + case EN_FCV: + case EN_PRV: + case EN_PSV: + if (hyd->LinkSetting[k] == MISSING) { + valvecoeff(pr, k); + } + else sol->P[k] = 0.0; + } +} + + +void matrixcoeffs(EN_Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: computes coefficients of linearized network eqns. +**-------------------------------------------------------------- +*/ +{ + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + EN_Network *net = &pr->network; + + memset(sol->Aii, 0, (net->Nnodes + 1) * sizeof(double)); + memset(sol->Aij, 0, (hyd->Ncoeffs + 1) * sizeof(double)); + memset(sol->F, 0, (net->Nnodes + 1) * sizeof(double)); + memset(hyd->X_tmp, 0, (net->Nnodes + 1) * sizeof(double)); + linkcoeffs(pr); + emittercoeffs(pr); + nodecoeffs(pr); + valvecoeffs(pr); +} + + +void linkcoeffs(EN_Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: computes matrix coefficients for links +**-------------------------------------------------------------- +*/ +{ + int k, n1, n2; + + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + Slink *link; + + // Examine each link of network */ + for (k = 1; k <= net->Nlinks; k++) + { + if (sol->P[k] == 0.0) continue; + link = &net->Link[k]; + n1 = link->N1; // Start node of link + n2 = link->N2; // End node of link + + + // Update net nodal inflows (X), solution matrix (A) and RHS array (F) + // (Use covention that flow out of node is (-), flow into node is (+)) + hyd->X_tmp[n1] -= hyd->LinkFlows[k]; + hyd->X_tmp[n2] += hyd->LinkFlows[k]; + + // Off-diagonal coeff. + sol->Aij[sol->Ndx[k]] -= sol->P[k]; + + // Node n1 is junction + if (n1 <= net->Njuncs) + { + sol->Aii[sol->Row[n1]] += sol->P[k]; // Diagonal coeff. + sol->F[sol->Row[n1]] += sol->Y[k]; // RHS coeff. + } + + // Node n1 is a tank + else { + sol->F[sol->Row[n2]] += (sol->P[k] * hyd->NodeHead[n1]); + } + + // Node n2 is junction + if (n2 <= net->Njuncs) { + sol->Aii[sol->Row[n2]] += sol->P[k]; // Diagonal coeff. + sol->F[sol->Row[n2]] -= sol->Y[k]; // RHS coeff. + } + + // Node n2 is a tank + else { + sol->F[sol->Row[n1]] += (sol->P[k] * hyd->NodeHead[n2]); + } + } +} + + +void nodecoeffs(EN_Project *pr) +/* +**---------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: completes calculation of nodal flow imbalance (X) +** & flow correction (F) arrays +**---------------------------------------------------------------- +*/ +{ + int i; + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + EN_Network *net = &pr->network; + + // For junction nodes, subtract demand flow from net + // flow imbalance & add imbalance to RHS array F. + for (i = 1; i <= net->Njuncs; i++) + { + hyd->X_tmp[i] -= hyd->NodeDemand[i]; + sol->F[sol->Row[i]] += hyd->X_tmp[i]; + } +} + + +void valvecoeffs(EN_Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: computes matrix coeffs. for PRVs, PSVs & FCVs +** whose status is not fixed to OPEN/CLOSED +**-------------------------------------------------------------- +*/ +{ + int i, k, n1, n2; + + hydraulics_t *hyd = &pr->hydraulics; + EN_Network *net = &pr->network; + Slink *link; + Svalve *valve; + + // Examine each valve + for (i = 1; i <= net->Nvalves; i++) + { + valve = &net->Valve[i]; + k = valve->Link; // Link index of valve + link = &net->Link[k]; + if (hyd->LinkSetting[k] == MISSING) { + continue; // Valve status fixed + } + n1 = link->N1; // Start & end nodes + n2 = link->N2; + switch (link->Type) // Call valve-specific function + { + case EN_PRV: + prvcoeff(pr, k, n1, n2); + break; + case EN_PSV: + psvcoeff(pr, k, n1, n2); + break; + case EN_FCV: + fcvcoeff(pr, k, n1, n2); + break; + default: continue; + } + } +} + + +void emittercoeffs(EN_Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: none +** Purpose: computes matrix coeffs. for emitters +** +** Note: Emitters consist of a fictitious pipe connected to +** a fictitious reservoir whose elevation equals that +** of the junction. The headloss through this pipe is +** Ke*(Flow)^hyd->Qexp, where Ke = emitter headloss coeff. +**-------------------------------------------------------------- +*/ +{ + int i; + double ke; + double p; + double q; + double y; + double z; + + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + EN_Network *net = &pr->network; + Snode *node; + + for (i = 1; i <= net->Njuncs; i++) + { + node = &net->Node[i]; + if (node->Ke == 0.0) { + continue; + } + ke = MAX(CSMALL, node->Ke); // emitter coeff. + q = hyd->EmitterFlows[i]; // emitter flow + z = ke * pow(ABS(q), hyd->Qexp); // emitter head loss + p = hyd->Qexp * z / ABS(q); // head loss gradient + if (p < hyd->RQtol) { + p = 1.0 / hyd->RQtol; + } + else { + p = 1.0 / p; // inverse head loss gradient + } + y = SGN(q)*z*p; // head loss / gradient + sol->Aii[sol->Row[i]] += p; // addition to main diagonal + sol->F[sol->Row[i]] += y + p * node->El; // addition to r.h.s. + hyd->X_tmp[i] -= q; // addition to net node inflow + } +} + + +void pipecoeff(EN_Project *pr, int k) +/* +**-------------------------------------------------------------- +** Input: k = link index +** Output: none +** Purpose: computes P & Y coefficients for pipe k +** +** P = inverse head loss gradient = 1/(dh/dQ) +** Y = flow correction term = h*P +**-------------------------------------------------------------- +*/ +{ + double hpipe, // Normal head loss + hml, // Minor head loss + ml, // Minor loss coeff. + p, // q*(dh/dq) + q, // Abs. value of flow + r; // Resistance coeff. + + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + Slink *link = &pr->network.Link[k]; + + // For closed pipe use headloss formula: h = CBIG*q + if (hyd->LinkStatus[k] <= CLOSED) + { + sol->P[k] = 1.0 / CBIG; + sol->Y[k] = hyd->LinkFlows[k]; + return; + } + + // ... head loss formula is Darcy-Weisbach + if (hyd->Formflag == DW) { + DWpipecoeff(pr, k); + return; + } + + // Evaluate headloss coefficients + q = ABS(hyd->LinkFlows[k]); // Absolute flow + ml = link->Km; // Minor loss coeff. + r = link->R; // Resistance coeff. + + // Use large P coefficient for small flow resistance product + if ( (r+ml)*q < hyd->RQtol) + { + sol->P[k] = 1.0 / hyd->RQtol; + sol->Y[k] = hyd->LinkFlows[k] / hyd->Hexp; + return; + } + + // Compute P and Y coefficients + hpipe = r*pow(q, hyd->Hexp); // Friction head loss + p = hyd->Hexp*hpipe; // Q*dh(friction)/dQ + if (ml > 0.0) + { + hml = ml*q*q; // Minor head loss + p += 2.0*hml; // Q*dh(Total)/dQ + } + else hml = 0.0; + p = hyd->LinkFlows[k] / p; // 1 / (dh/dQ) + sol->P[k] = ABS(p); + sol->Y[k] = p*(hpipe + hml); +} + + +void DWpipecoeff(EN_Project *pr, int k) +/* +**-------------------------------------------------------------- +** Input: k = link index +** Output: none +** Purpose: computes pipe head loss coeffs. for Darcy-Weisbach +** formula. +**-------------------------------------------------------------- +*/ +{ + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + Slink *link = &pr->network.Link[k]; + + double s = hyd->Viscos * link->Diam; + double q = ABS(hyd->LinkFlows[k]); + double dfdq = 0.0; + double e; // relative roughness height + double r; // total resistance coeff. + double f; // friction factor + double hloss; // head loss + double hgrad; // head loss gradient + + // ... use Hagen-Poiseuille formula for laminar flow (Re <= 2000) + if (q <= A2 * s) + { + r = 16.0 * PI * s * link->R; + hloss = q * (r + link->Km * q); + hgrad = r + 2.0 * link->Km * q; + } + + // ... use Colebrook formula for turbulent flow + else + { + e = link->Kc / link->Diam; + f = frictionFactor(q, e, s, &dfdq); + r = f * link->R + link->Km; + hloss = r * q * hyd->LinkFlows[k]; + hgrad = (2.0 * r * q) + (dfdq * link->R * q * q); + } + + // ... head loss has same sign as flow + hloss *= SGN(hyd->LinkFlows[k]); + + // ... compute P and Y coeffs. + sol->P[k] = 1.0 / hgrad; + sol->Y[k] = hloss / hgrad; +} + + +double frictionFactor(double q, double e, double s, double *dfdq) +/* +**-------------------------------------------------------------- +** Input: q = flow rate (cfs) +** e = pipe roughness / diameter +** s = viscosity * diameter +** Output: f = friction factor +** dfdq = derivative of f w.r.t. flow +** Purpose: computes Darcy-Weisbach friction factor and its +** derivative as a function of Reynolds Number (Re). +**-------------------------------------------------------------- +*/ +{ + double f; + double x1, x2, x3, x4, + y1, y2, y3, + fa, fb, r; + double w = q / s; // Re*Pi/4 + + // ... For Re >= 4000 use Colebrook Formula + if (w >= A1) + { + y1 = A8 / pow(w, 0.9); + y2 = e / 3.7 + y1; + y3 = A9 * log(y2); + f = 1.0 / (y3*y3); + *dfdq = 1.8 * f * y1 * A9 / y2 / y3 / q; + } + + // ... Use interpolating polynomials developed by + // E. Dunlop for transition flow from 2000 < Re < 4000. + else + { + y2 = e / 3.7 + AB; + y3 = A9 * log(y2); + fa = 1.0 / (y3*y3); + fb = (2.0 + AC / (y2*y3)) * fa; + r = w / A2; + x1 = 7.0 * fa - fb; + x2 = 0.128 - 17.0 * fa + 2.5 * fb; + x3 = -0.128 + 13.0 * fa - (fb + fb); + x4 = r * (0.032 - 3.0 * fa + 0.5 *fb); + f = x1 + r * (x2 + r * (x3 + x4)); + *dfdq = (x2 + 2.0 * r * (x3 + x4)) / s / A2; + } + return f; +} + + +void pumpcoeff(EN_Project *pr, int k) +/* +**-------------------------------------------------------------- +** Input: k = link index +** Output: none +** Purpose: computes P & Y coeffs. for pump in link k +**-------------------------------------------------------------- +*/ +{ + int p; // Pump index + double h0, // Shutoff head + q, // Abs. value of flow + r, // Flow resistance coeff. + n; // Flow exponent coeff. + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + double setting = hyd->LinkSetting[k]; + Spump *pump; + + // Use high resistance pipe if pump closed or cannot deliver head + if (hyd->LinkStatus[k] <= CLOSED || setting == 0.0) { + sol->P[k] = 1.0 / CBIG; + sol->Y[k] = hyd->LinkFlows[k]; + return; + } + + q = ABS(hyd->LinkFlows[k]); + q = MAX(q, TINY); + + // Obtain reference to pump object + p = findpump(&pr->network, k); + pump = &pr->network.Pump[p]; + + // Get pump curve coefficients for custom pump curve. + if (pump->Ptype == CUSTOM) + { + // Find intercept (h0) & slope (r) of pump curve + // line segment which contains speed-adjusted flow. + curvecoeff(pr, pump->Hcurve, q / setting, &h0, &r); + + // Determine head loss coefficients. + pump->H0 = -h0; + pump->R = -r; + pump->N = 1.0; + } + + // Adjust head loss coefficients for pump speed. + h0 = SQR(setting) * pump->H0; + n = pump->N; + r = pump->R * pow(setting, 2.0 - n); + if (n != 1.0) { + r = n * r * pow(q, n - 1.0); + } + + // Compute inverse headloss gradient (P) and flow correction factor (Y) + sol->P[k] = 1.0 / MAX(r, hyd->RQtol); + sol->Y[k] = hyd->LinkFlows[k] / n + sol->P[k] * h0; +} + + +void curvecoeff(EN_Project *pr, int i, double q, double *h0, double *r) +/* +**------------------------------------------------------------------- +** Input: i = curve index +** q = flow rate +** Output: *h0 = head at zero flow (y-intercept) +** *r = dHead/dFlow (slope) +** Purpose: computes intercept and slope of head v. flow curve +** at current flow. +**------------------------------------------------------------------- +*/ +{ + int k1, k2, npts; + double *x, *y; + Scurve *curve; + + // Remember that curve is stored in untransformed units + q *= pr->Ucf[FLOW]; + curve = &pr->network.Curve[i]; + x = curve->X; // x = flow + y = curve->Y; // y = head + npts = curve->Npts; + + // Find linear segment of curve that brackets flow q + k2 = 0; + while (k2 < npts && x[k2] < q) + k2++; + + if (k2 == 0) + k2++; + + else if (k2 == npts) + k2--; + + k1 = k2 - 1; + + // Compute slope and intercept of this segment + *r = (y[k2] - y[k1]) / (x[k2] - x[k1]); + *h0 = y[k1] - (*r)*x[k1]; + + // Convert units + *h0 = (*h0) / pr->Ucf[HEAD]; + *r = (*r) * pr->Ucf[FLOW] / pr->Ucf[HEAD]; +} + + +void gpvcoeff(EN_Project *pr, int k) +/* +**-------------------------------------------------------------- +** Input: k = link index +** Output: none +** Purpose: computes P & Y coeffs. for general purpose valve +**-------------------------------------------------------------- +*/ +{ + double h0, // Headloss curve intercept + q, // Abs. value of flow + r; // Flow resistance coeff. + + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + + // Treat as a pipe if valve closed + if (hyd->LinkStatus[k] == CLOSED) { + valvecoeff(pr, k); + } + + // Otherwise utilize headloss curve + // whose index is stored in K + else { + // Find slope & intercept of headloss curve. + q = ABS(hyd->LinkFlows[k]); + q = MAX(q, TINY); + curvecoeff(pr, (int)ROUND(hyd->LinkSetting[k]), q, &h0, &r); + + // Compute inverse headloss gradient (P) + // and flow correction factor (Y). + sol->P[k] = 1.0 / MAX(r, hyd->RQtol); + sol->Y[k] = sol->P[k] * (h0 + r*q) * SGN(hyd->LinkFlows[k]); + } +} + + +void pbvcoeff(EN_Project *pr, int k) +/* +**-------------------------------------------------------------- +** Input: k = link index +** Output: none +** Purpose: computes P & Y coeffs. for pressure breaker valve +**-------------------------------------------------------------- +*/ +{ + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + Slink *link = &pr->network.Link[k]; + + // If valve fixed OPEN or CLOSED then treat as a pipe + if (hyd->LinkSetting[k] == MISSING || hyd->LinkSetting[k] == 0.0) { + valvecoeff(pr, k); + } + + // If valve is active + else { + // Treat as a pipe if minor loss > valve setting + if (link->Km * SQR(hyd->LinkFlows[k]) > hyd->LinkSetting[k]) { + valvecoeff(pr, k); + } + // Otherwise force headloss across valve to be equal to setting + else { + sol->P[k] = CBIG; + sol->Y[k] = hyd->LinkSetting[k] * CBIG; + } + } +} + + +void tcvcoeff(EN_Project *pr, int k) +/* +**-------------------------------------------------------------- +** Input: k = link index +** Output: none +** Purpose: computes P & Y coeffs. for throttle control valve +**-------------------------------------------------------------- +*/ +{ + double km; + hydraulics_t *hyd = &pr->hydraulics; + Slink *link = &pr->network.Link[k]; + + // Save original loss coeff. for open valve + km = link->Km; + + // If valve not fixed OPEN or CLOSED, compute its loss coeff. + if (hyd->LinkSetting[k] != MISSING) { + link->Km = 0.02517 * hyd->LinkSetting[k] / (SQR(link->Diam)*SQR(link->Diam)); + } + + // Then apply usual valve formula + valvecoeff(pr, k); + + // Restore original loss coeff. + link->Km = km; +} + + +void prvcoeff(EN_Project *pr, int k, int n1, int n2) +/* +**-------------------------------------------------------------- +** Input: k = link index +** n1 = upstream node of valve +** n2 = downstream node of valve +** Output: none +** Purpose: computes solution matrix coeffs. for pressure +** reducing valves +**-------------------------------------------------------------- +*/ +{ + int i, j; // Rows of solution matrix + double hset; // Valve head setting + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + + i = sol->Row[n1]; // Matrix rows of nodes + j = sol->Row[n2]; + hset = pr->network.Node[n2].El + + hyd->LinkSetting[k]; // Valve setting + + if (hyd->LinkStatus[k] == ACTIVE) + { + + // Set coeffs. to force head at downstream + // node equal to valve setting & force flow + // to equal to flow imbalance at downstream node. + + sol->P[k] = 0.0; + sol->Y[k] = hyd->LinkFlows[k] + hyd->X_tmp[n2]; // Force flow balance + sol->F[j] += (hset * CBIG); // Force head = hset + sol->Aii[j] += CBIG; // at downstream node + if (hyd->X_tmp[n2] < 0.0) { + sol->F[i] += hyd->X_tmp[n2]; + } + return; + } + + // For OPEN, CLOSED, or XPRESSURE valve + // compute matrix coeffs. using the valvecoeff() function. + + valvecoeff(pr, k); + sol->Aij[sol->Ndx[k]] -= sol->P[k]; + sol->Aii[i] += sol->P[k]; + sol->Aii[j] += sol->P[k]; + sol->F[i] += (sol->Y[k] - hyd->LinkFlows[k]); + sol->F[j] -= (sol->Y[k] - hyd->LinkFlows[k]); +} + + +void psvcoeff(EN_Project *pr, int k, int n1, int n2) +/* +**-------------------------------------------------------------- +** Input: k = link index +** n1 = upstream node of valve +** n2 = downstream node of valve +** Output: none +** Purpose: computes solution matrix coeffs. for pressure +** sustaining valve +**-------------------------------------------------------------- +*/ +{ + int i, j; // Rows of solution matrix + double hset; // Valve head setting + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + + i = sol->Row[n1]; // Matrix rows of nodes + j = sol->Row[n2]; + hset = pr->network.Node[n1].El + + hyd->LinkSetting[k]; // Valve setting + + if (hyd->LinkStatus[k] == ACTIVE) + { + // Set coeffs. to force head at upstream + // node equal to valve setting & force flow + // equal to flow imbalance at upstream node. + + sol->P[k] = 0.0; + sol->Y[k] = hyd->LinkFlows[k] - hyd->X_tmp[n1]; // Force flow balance + sol->F[i] += (hset * CBIG); // Force head = hset + sol->Aii[i] += CBIG; // at upstream node + if (hyd->X_tmp[n1] > 0.0) { + sol->F[j] += hyd->X_tmp[n1]; + } + return; + } + + // For OPEN, CLOSED, or XPRESSURE valve + // compute matrix coeffs. using the valvecoeff() function. + + valvecoeff(pr, k); + sol->Aij[sol->Ndx[k]] -= sol->P[k]; + sol->Aii[i] += sol->P[k]; + sol->Aii[j] += sol->P[k]; + sol->F[i] += (sol->Y[k] - hyd->LinkFlows[k]); + sol->F[j] -= (sol->Y[k] - hyd->LinkFlows[k]); +} + + +void fcvcoeff(EN_Project *pr, int k, int n1, int n2) +/* +**-------------------------------------------------------------- +** Input: k = link index +** n1 = upstream node of valve +** n2 = downstream node of valve +** Output: none +** Purpose: computes solution matrix coeffs. for flow control +** valve +**-------------------------------------------------------------- +*/ +{ + int i, j; // Rows in solution matrix + double q; // Valve flow setting + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + + q = hyd->LinkSetting[k]; + i = hyd->solver.Row[n1]; + j = hyd->solver.Row[n2]; + + // If valve active, break network at valve and treat + // flow setting as external demand at upstream node + // and external supply at downstream node. + + if (hyd->LinkStatus[k] == ACTIVE) + { + hyd->X_tmp[n1] -= q; + sol->F[i] -= q; + hyd->X_tmp[n2] += q; + sol->F[j] += q; + sol->P[k] = 1.0 / CBIG; + sol->Aij[sol->Ndx[k]] -= sol->P[k]; + sol->Aii[i] += sol->P[k]; + sol->Aii[j] += sol->P[k]; + sol->Y[k] = hyd->LinkFlows[k] - q; + } + + // Otherwise treat valve as an open pipe + + else + { + valvecoeff(pr, k); + sol->Aij[sol->Ndx[k]] -= sol->P[k]; + sol->Aii[i] += sol->P[k]; + sol->Aii[j] += sol->P[k]; + sol->F[i] += (sol->Y[k] - hyd->LinkFlows[k]); + sol->F[j] -= (sol->Y[k] - hyd->LinkFlows[k]); + } +} + + +void valvecoeff(EN_Project *pr, int k) +/* +**-------------------------------------------------------------- +** Input: k = link index +** Output: none +** Purpose: computes solution matrix coeffs. for a completely +** open, closed, or throttled control valve. +**-------------------------------------------------------------- +*/ +{ + double p; + + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + Slink *link = &net->Link[k]; + + double flow = hyd->LinkFlows[k]; + + // Valve is closed. Use a very small matrix coeff. + if (hyd->LinkStatus[k] <= CLOSED) + { + sol->P[k] = 1.0 / CBIG; + sol->Y[k] = flow; + return; + } + + // Account for any minor headloss through the valve + if (link->Km > 0.0) + { + p = 2.0 * link->Km * fabs(flow); + if (p < hyd->RQtol) { + p = hyd->RQtol; + } + sol->P[k] = 1.0 / p; + sol->Y[k] = flow / 2.0; + } + else + { + sol->P[k] = 1.0 / hyd->RQtol; + sol->Y[k] = flow; + } +} diff --git a/src/hydraul.c b/src/hydraul.c index f4eb45c..37664a8 100755 --- a/src/hydraul.c +++ b/src/hydraul.c @@ -29,13 +29,13 @@ AUTHOR: L. Rossman closehyd() -- called from ENcloseH() in EPANET.C tankvolume() -- called from ENsetnodevalue() in EPANET.C setlinkstatus(), - setlinksetting(), - resistance()-- all called from ENsetlinkvalue() in EPANET.C + setlinksetting() -- all called from ENsetlinkvalue() in EPANET.C External functions called by this module are: createsparse() -- see SMATRIX.C freesparse() -- see SMATRIX.C - linsolve() -- see SMATRIX.C + resistcoeff() -- see HYDCOEFFS.C + hydsolve() -- see HYDSOLVER.C checkrules() -- see RULES.C interp() -- see EPANET.C savehyd() -- see OUTPUT.C @@ -54,34 +54,11 @@ AUTHOR: L. Rossman #include #endif #include -#include "hash.h" #include "text.h" #include "types.h" -#include "epanet2.h" #include "funcs.h" -#define EXTERN extern -#include "vars.h" #define QZERO 1.e-6 /* Equivalent to zero flow */ -#define CBIG 1.e8 /* Big coefficient */ -#define CSMALL 1.e-6 /* Small coefficient */ - -/* Constants used for computing Darcy-Weisbach friction factor */ -#define A1 0.314159265359e04 /* 1000*PI */ -#define A2 0.157079632679e04 /* 500*PI */ -#define A3 0.502654824574e02 /* 16*PI */ -#define A4 6.283185307 /* 2*PI */ -#define A8 4.61841319859 /* 5.74*(PI/4)^.9 */ -#define A9 -8.685889638e-01 /* -2/ln(10) */ -#define AA -1.5634601348 /* -2*.9*2/ln(10) */ -#define AB 3.28895476345e-03 /* 5.74/(4000^.9) */ -#define AC -5.14214965799e-03 /* AA*AB */ - -/*** Updated 3/1/01 ***/ - -/* Function to find flow coeffs. through open/closed valves */ -void valvecoeff(EN_Project *pr, int k); - int openhyd(EN_Project *pr) /* @@ -92,15 +69,15 @@ int openhyd(EN_Project *pr) *-------------------------------------------------------------- */ { - int i; - int errcode = 0; - ERRCODE(createsparse(pr)); /* See SMATRIX.C */ - ERRCODE(allocmatrix(pr)); /* Allocate solution matrices */ - for (i=1; i <= pr->network.Nlinks; i++) { /* Initialize flows */ - Slink *link = &pr->network.Link[i]; - initlinkflow(pr, i, link->Stat, link->Kc); - } - return(errcode); + int i; + int errcode = 0; + ERRCODE(createsparse(pr)); /* See SMATRIX.C */ + ERRCODE(allocmatrix(pr)); /* Allocate solution matrices */ + for (i=1; i <= pr->network.Nlinks; i++) { /* Initialize flows */ + Slink *link = &pr->network.Link[i]; + initlinkflow(pr, i, link->Stat, link->Kc); + } + return(errcode); } @@ -115,72 +92,80 @@ void inithyd(EN_Project *pr, int initflag) **-------------------------------------------------------------- */ { - int i,j; + int i,j; - time_options_t *time = &pr->time_options; - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - out_file_t *out = &pr->out_files; + time_options_t *time = &pr->time_options; + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + out_file_t *out = &pr->out_files; + Stank *tank; + Slink *link; + Spump *pump; - /* Initialize tanks */ - for (i=1; i <= net->Ntanks; i++) { - Stank *tank = &net->Tank[i]; - tank->V = tank->V0; - hyd->NodeHead[tank->Node] = tank->H0; - hyd->NodeDemand[tank->Node] = 0.0; - hyd->OldStat[net->Nlinks+i] = TEMPCLOSED; - } - - /* Initialize emitter flows */ - memset(hyd->EmitterFlows,0,(net->Nnodes+1)*sizeof(double)); - for (i=1; i <= net->Njuncs; i++) { - if (net->Node[i].Ke > 0.0) { - hyd->EmitterFlows[i] = 1.0; + /* Initialize tanks */ + for (i=1; i <= net->Ntanks; i++) { + tank = &net->Tank[i]; + tank->V = tank->V0; + hyd->NodeHead[tank->Node] = tank->H0; + hyd->NodeDemand[tank->Node] = 0.0; + hyd->OldStat[net->Nlinks+i] = TEMPCLOSED; } - } - /* Initialize links */ - for (i=1; i <= net->Nlinks; i++) { - Slink *link = &net->Link[i]; - /* Initialize status and setting */ - hyd->LinkStatus[i] = link->Stat; - hyd->LinkSetting[i] = link->Kc; - /* Start active control valves in ACTIVE position */ - if ( - (link->Type == EN_PRV || link->Type == EN_PSV - || link->Type == EN_FCV) - && (link->Kc != MISSING) - ) hyd->LinkStatus[i] = ACTIVE; + /* Initialize emitter flows */ + memset(hyd->EmitterFlows,0,(net->Nnodes+1)*sizeof(double)); + for (i=1; i <= net->Njuncs; i++) { + if (net->Node[i].Ke > 0.0) { + hyd->EmitterFlows[i] = 1.0; + } + } + + /* Initialize links */ + for (i=1; i <= net->Nlinks; i++) { + link = &net->Link[i]; + + /* Initialize status and setting */ + hyd->LinkStatus[i] = link->Stat; + hyd->LinkSetting[i] = link->Kc; + + /* Compute flow resistance */ + resistcoeff(pr, i); + + /* Start active control valves in ACTIVE position */ + if ( + (link->Type == EN_PRV || link->Type == EN_PSV + || link->Type == EN_FCV) && (link->Kc != MISSING) + ) hyd->LinkStatus[i] = ACTIVE; /*** Updated 3/1/01 ***/ - /* Initialize flows if necessary */ - if (hyd->LinkStatus[i] <= CLOSED) hyd->LinkFlows[i] = QZERO; - else if (ABS(hyd->LinkFlows[i]) <= QZERO || initflag > 0) - initlinkflow(pr, i, hyd->LinkStatus[i], hyd->LinkSetting[i]); + /* Initialize flows if necessary */ + if (hyd->LinkStatus[i] <= CLOSED) hyd->LinkFlows[i] = QZERO; + else if (ABS(hyd->LinkFlows[i]) <= QZERO || initflag > 0) + initlinkflow(pr, i, hyd->LinkStatus[i], hyd->LinkSetting[i]); - /* Save initial status */ - hyd->OldStat[i] = hyd->LinkStatus[i]; - } + /* Save initial status */ + hyd->OldStat[i] = hyd->LinkStatus[i]; + } - /* Reset pump energy usage */ - for (i=1; i <= net->Npumps; i++) - { - for (j=0; j<6; j++) { - net->Pump[i].Energy[j] = 0.0; - } - } + /* Initialize pump energy usage */ + for (i=1; i <= net->Npumps; i++) + { + pump = &net->Pump[i]; + for (j = 0; j < MAX_ENERGY_STATS; j++) { + pump->Energy[j] = 0.0; + } + } - /* Re-position hydraulics file */ - if (pr->save_options.Saveflag) { - fseek(out->HydFile,out->HydOffset,SEEK_SET); - } + /* Re-position hydraulics file */ + if (pr->save_options.Saveflag) { + fseek(out->HydFile,out->HydOffset,SEEK_SET); + } /*** Updated 3/1/01 ***/ - /* Initialize current time */ - hyd->Haltflag = 0; - time->Htime = 0; - time->Hydstep = 0; - time->Rtime = time->Rstep; + /* Initialize current time */ + hyd->Haltflag = 0; + time->Htime = 0; + time->Hydstep = 0; + time->Rtime = time->Rstep; } @@ -194,42 +179,42 @@ int runhyd(EN_Project *pr, long *t) **-------------------------------------------------------------- */ { - int iter; /* Iteration count */ - int errcode; /* Error code */ - double relerr; /* Solution accuracy */ + int iter; /* Iteration count */ + int errcode; /* Error code */ + double relerr; /* Solution accuracy */ - hydraulics_t *hyd = &pr->hydraulics; - time_options_t *time = &pr->time_options; - report_options_t *rep = &pr->report; + hydraulics_t *hyd = &pr->hydraulics; + time_options_t *time = &pr->time_options; + report_options_t *rep = &pr->report; - /* Find new demands & control actions */ - *t = time->Htime; - demands(pr); - controls(pr); + /* Find new demands & control actions */ + *t = time->Htime; + demands(pr); + controls(pr); - /* Solve network hydraulic equations */ - errcode = netsolve(pr,&iter,&relerr); - if (!errcode) { - /* Report new status & save results */ - if (rep->Statflag) { - writehydstat(pr,iter,relerr); - } + /* Solve network hydraulic equations */ + errcode = hydsolve(pr,&iter,&relerr); + if (!errcode) { + /* Report new status & save results */ + if (rep->Statflag) { + writehydstat(pr,iter,relerr); + } - /* solution info */ - hyd->relativeError = relerr; - hyd->iterations = iter; + /* solution info */ + hyd->relativeError = relerr; + hyd->iterations = iter; /*** Updated 3/1/01 ***/ - /* If system unbalanced and no extra trials */ - /* allowed, then activate the Haltflag. */ - if (relerr > hyd->Hacc && hyd->ExtraIter == -1) { - hyd->Haltflag = 1; - } + /* If system unbalanced and no extra trials */ + /* allowed, then activate the Haltflag. */ + if (relerr > hyd->Hacc && hyd->ExtraIter == -1) { + hyd->Haltflag = 1; + } - /* Report any warning conditions */ - if (!errcode) { - errcode = writehydwarn(pr,iter,relerr); - } + /* Report any warning conditions */ + if (!errcode) { + errcode = writehydwarn(pr,iter,relerr); + } } return(errcode); } /* end of runhyd */ @@ -459,7 +444,8 @@ void setlinkflow(EN_Project *pr, int k, double dh) dh = -dh * pr->Ucf[HEAD] / SQR(hyd->LinkSetting[k]); i = net->Pump[p].Hcurve; curve = &net->Curve[i]; - hyd->LinkFlows[k] = interp(curve->Npts,curve->Y,curve->X,dh) * hyd->LinkSetting[k] / pr->Ucf[FLOW]; + hyd->LinkFlows[k] = interp(curve->Npts,curve->Y,curve->X,dh) * + hyd->LinkSetting[k] / pr->Ucf[FLOW]; } /* Otherwise use inverse of power curve */ @@ -565,55 +551,6 @@ void setlinksetting(EN_Project *pr, int index, double value, StatType *s, doubl } -void resistance(EN_Project *pr, int k) -/* -**-------------------------------------------------------------------- -** Input: k = link index -** Output: none -** Purpose: computes link flow resistance -**-------------------------------------------------------------------- -*/ -{ - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - - double e,d,L; - - Slink *link = &net->Link[k]; - link->R = CSMALL; - - - switch (link->Type) { - /* Link is a pipe. Compute resistance based on headloss formula. */ - /* Friction factor for D-W formula gets included during solution */ - /* process in pipecoeff() function. */ - case EN_CVPIPE: - case EN_PIPE: - e = link->Kc; /* Roughness coeff. */ - d = link->Diam; /* Diameter */ - L = link->Len; /* Length */ - switch(hyd->Formflag) - { - case HW: - link->R = 4.727*L/pow(e,hyd->Hexp)/pow(d,4.871); - break; - case DW: - link->R = L/2.0/32.2/d/SQR(PI*SQR(d)/4.0); - break; - case CM: - link->R = SQR(4.0*e/(1.49*PI*d*d)) * pow((d/4.0),-1.333)*L; - } - break; - - /* Link is a pump. Use negligible resistance. */ - case EN_PUMP: - link->R = CBIG; //CSMALL; - break; - default: - break; - } -} - void demands(EN_Project *pr) /* @@ -1028,7 +965,7 @@ void addenergy(EN_Project *pr, long hstep) { int i,j,k; long m,n; - double c0,c, /* Energy cost (cost/kwh) */ + double c0,c, /* Energy cost (cost/kwh) */ f0, /* Energy cost factor */ dt, /* Time interval (hr) */ e, /* Pump efficiency (fraction) */ @@ -1038,6 +975,7 @@ void addenergy(EN_Project *pr, long hstep) EN_Network *net = &pr->network; hydraulics_t *hyd = &pr->hydraulics; + Spump *pump; /* Determine current time interval in hours */ if (pr->time_options.Dur == 0) { @@ -1053,7 +991,7 @@ void addenergy(EN_Project *pr, long hstep) if (dt == 0.0) { return; } - n = (pr->time_options.Htime + pr->time_options.Pstart) / pr->time_options.Pstep; + n = (pr->time_options.Htime + pr->time_options.Pstart) / pr->time_options.Pstep; /* Compute default energy cost at current time */ c0 = hyd->Ecost; @@ -1067,17 +1005,17 @@ void addenergy(EN_Project *pr, long hstep) /* Examine each pump */ for (j=1; j <= net->Npumps; j++) { - Spump *pump = &net->Pump[j]; /* Skip closed pumps */ + pump = &net->Pump[j]; k = pump->Link; - if (hyd->LinkStatus[k] <= CLOSED) { + if (hyd->LinkStatus[k] <= CLOSED) { continue; - } + } q = MAX(QZERO, ABS(hyd->LinkFlows[k])); /* Find pump-specific energy cost */ - if (pump->Ecost > 0.0) { - c = pump->Ecost; + if (pump->Ecost > 0.0) { + c = pump->Ecost; } else { c = c0; @@ -1096,12 +1034,12 @@ void addenergy(EN_Project *pr, long hstep) psum += p; /* Update pump's cumulative statistics */ - pump->Energy[0] += dt; /* Time on-line */ - pump->Energy[1] += e*dt; /* Effic.-hrs */ - pump->Energy[2] += p/q*dt; /* kw/cfs-hrs */ - pump->Energy[3] += p*dt; /* kw-hrs */ - pump->Energy[4] = MAX(pump->Energy[4],p); - pump->Energy[5] += c*p*dt; /* cost-hrs. */ + pump->Energy[PCNT_ONLINE] += dt; + pump->Energy[PCNT_EFFIC] += e*dt; + pump->Energy[KWH_PER_FLOW] += p/q*dt; + pump->Energy[TOTAL_KWH] += p*dt; + pump->Energy[MAX_KW] = MAX(pump->Energy[MAX_KW],p); + pump->Energy[TOTAL_COST] += c*p*dt; } /* Update maximum kw value */ @@ -1119,12 +1057,16 @@ void getenergy(EN_Project *pr, int k, double *kw, double *eff) **---------------------------------------------------------------- */ { - int i,j; - double dh, q, e; - double q4eff; //q4eff=flow at nominal speed to compute efficiency + int i, // efficiency curve index + j; // pump index + double dh, // head across pump (ft) + q, // flow through pump (cfs) + e; // pump efficiency + double q4eff; // flow at nominal pump speed of 1.0 + double speed; // current speed setting Scurve *curve; - EN_Network *net = &pr->network; + EN_Network *net = &pr->network; hydraulics_t *hyd = &pr->hydraulics; Slink *link = &net->Link[k]; @@ -1147,12 +1089,15 @@ void getenergy(EN_Project *pr, int k, double *kw, double *eff) { j = findpump(net,k); e = hyd->Epump; + speed = hyd->LinkSetting[k]; if ( (i = net->Pump[j].Ecurve) > 0) - { - q4eff = q / hyd->LinkSetting[k]; + { + q4eff = q / speed * pr->Ucf[FLOW]; curve = &net->Curve[i]; - e = interp(curve->Npts,curve->X, curve->Y, q4eff * pr->Ucf[FLOW]); - } + e = interp(curve->Npts,curve->X, curve->Y, q4eff); + /* Sarbu and Borza pump speed adjustment */ + e = 100.0 - ((100.0-e) * pow(1.0/speed, 0.1)); + } e = MIN(e, 100.0); e = MAX(e, 1.0); e /= 100.0; @@ -1219,6 +1164,8 @@ double tankvolume(EN_Project *pr, int i, double h) int j; EN_Network *net = &pr->network; Stank *tank = &net->Tank[i]; + Scurve *curve; + /* Use level*area if no volume curve */ j = tank->Vcurve; if (j == 0) { @@ -1228,8 +1175,10 @@ double tankvolume(EN_Project *pr, int i, double h) /* If curve exists, interpolate on h to find volume v */ /* remembering that volume curve is in original units.*/ else { - Scurve *curve = &net->Curve[j]; - return(interp(curve->Npts, curve->X, curve->Y, (h - net->Node[tank->Node].El) * pr->Ucf[HEAD]) / pr->Ucf[VOLUME]); + curve = &net->Curve[j]; + return(interp(curve->Npts, curve->X, curve->Y, + (h - net->Node[tank->Node].El) * + pr->Ucf[HEAD]) / pr->Ucf[VOLUME]); } } /* End of tankvolume */ @@ -1265,1688 +1214,4 @@ double tankgrade(EN_Project *pr, int i, double v) } /* End of tankgrade */ -int netsolve(EN_Project *pr, int *iter, double *relerr) -/* -**------------------------------------------------------------------- -** Input: none -** Output: *iter = # of iterations to reach solution -** *relerr = convergence error in solution -** returns error code -** Purpose: solves network nodal equations for heads and flows -** using Todini's Gradient algorithm -** -*** Updated 9/7/00 *** -*** Updated 2.00.11 *** -*** Updated 2.00.12 *** -** Notes: Status checks on CVs, pumps and pipes to tanks are made -** every hyd->CheckFreq iteration, up until MaxCheck iterations -** are reached. Status checks on control valves are made -** every iteration if DampLimit = 0 or only when the -** convergence error is at or below DampLimit. If DampLimit -** is > 0 then future computed flow changes are only 60% of -** their full value. A complete status check on all links -** is made when convergence is achieved. If convergence is -** not achieved in hyd->MaxIter trials and hyd->ExtraIter > 0 then -** another hyd->ExtraIter trials are made with no status changes -** made to any links and a warning message is generated. -** -** This procedure calls linsolve() which appears in SMATRIX.C. -**------------------------------------------------------------------- -*/ -{ - int i; /* Node index */ - int errcode = 0; /* Node causing solution error */ - int nextcheck; /* Next status check trial */ - int maxtrials; /* Max. trials for convergence */ - double newerr; /* New convergence error */ - int valveChange; /* Valve status change flag */ - int statChange; - - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - solver_t *sol = &hyd->solver; - report_options_t *rep = &pr->report; - - /* Initialize status checking & relaxation factor */ - nextcheck = hyd->CheckFreq; - hyd->RelaxFactor = 1.0; - - /* Repeat iterations until convergence or trial limit is exceeded. */ - /* (hyd->ExtraIter used to increase trials in case of status cycling.) */ - if (pr->report.Statflag == FULL) { - writerelerr(pr,0,0); - } - maxtrials = hyd->MaxIter; - if (hyd->ExtraIter > 0) { - maxtrials += hyd->ExtraIter; - } - *iter = 1; - while (*iter <= maxtrials) { - /* - ** Compute coefficient matrices A & F and solve A*H = F - ** where H = heads, A = Jacobian coeffs. derived from - ** head loss gradients, & F = flow correction terms. - ** Solution for H is returned in F from call to linsolve(). - */ - newcoeffs(pr); - errcode = linsolve(&hyd->solver,net->Njuncs); - - /* Take action depending on error code */ - if (errcode < 0) { - break; /* Memory allocation problem */ - } - if (errcode > 0) { /* Ill-conditioning problem */ - /* If control valve causing problem, fix its status & continue, */ - /* otherwise end the iterations with no solution. */ - if (badvalve(pr,sol->Order[errcode])) { - continue; - } - else break; - } - - /* Update current solution. */ - /* (Row[i] = row of solution matrix corresponding to node i). */ - for (i=1; i<=net->Njuncs; i++) { - hyd->NodeHead[i] = sol->F[sol->Row[i]]; /* Update heads */ - } - newerr = newflows(pr); /* Update flows */ - *relerr = newerr; - - /* Write convergence error to status report if called for */ - if (rep->Statflag == FULL) { - writerelerr(pr, *iter,*relerr); - } - - /* Apply solution damping & check for change in valve status */ - hyd->RelaxFactor = 1.0; - valveChange = FALSE; - if ( hyd->DampLimit > 0.0 ) { - if( *relerr <= hyd->DampLimit ) { - hyd->RelaxFactor = 0.6; - valveChange = valvestatus(pr); - } - } - else { - valveChange = valvestatus(pr); - } - - /* Check for convergence */ - if (*relerr <= hyd->Hacc) { - /* We have convergence. Quit if we are into extra iterations. */ - if (*iter > hyd->MaxIter) { - break; - } - - /* Quit if no status changes occur. */ - statChange = FALSE; - if (valveChange) { - statChange = TRUE; - } - if (linkstatus(pr)) { - statChange = TRUE; - } - if (pswitch(pr)) { - statChange = TRUE; - } - if (!statChange) { - break; - } - - /* We have a status change so continue the iterations */ - nextcheck = *iter + hyd->CheckFreq; - } - - /* No convergence yet. See if its time for a periodic status */ - /* check on pumps, CV's, and pipes connected to tanks. */ - else if (*iter <= hyd->MaxCheck && *iter == nextcheck) - { - linkstatus(pr); - nextcheck += hyd->CheckFreq; - } - (*iter)++; - } - - /* Iterations ended. Report any errors. */ - if (errcode < 0) errcode = 101; /* Memory allocation error */ - else if (errcode > 0) - { - writehyderr(pr, sol->Order[errcode]); /* Ill-conditioned eqns. error */ - errcode = 110; - } - - /* Add any emitter flows to junction demands */ - for (i=1; i<=net->Njuncs; i++) { - hyd->NodeDemand[i] += hyd->EmitterFlows[i]; - } - return(errcode); -} /* End of netsolve */ - - -int badvalve(EN_Project *pr, int n) -/* - **----------------------------------------------------------------- -** Input: n = node index -** Output: returns 1 if node n belongs to an active control valve, -** 0 otherwise -** Purpose: determines if a node belongs to an active control valve -** whose setting causes an inconsistent set of eqns. If so, -** the valve status is fixed open and a warning condition -** is generated. -**----------------------------------------------------------------- -*/ -{ - int i,k,n1,n2; - Slink *link; - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - report_options_t *rep = &pr->report; - time_options_t *time = &pr->time_options; - - for (i=1; i <= net->Nvalves; i++) { - k = net->Valve[i].Link; - link = &net->Link[k]; - n1 = link->N1; - n2 = link->N2; - if (n == n1 || n == n2) { - EN_LinkType t = link->Type; - if (t == EN_PRV || t == EN_PSV || t == EN_FCV) { - if (hyd->LinkStatus[k] == ACTIVE) { - if (rep->Statflag == FULL) { - sprintf(pr->Msg,FMT61,clocktime(rep->Atime,time->Htime),link->ID); - writeline(pr, pr->Msg); - } - if (link->Type == EN_FCV) { - hyd->LinkStatus[k] = XFCV; - } - else { - hyd->LinkStatus[k] = XPRESSURE; - } - return(1); - } - } - return(0); - } - } - return(0); -} - - -int valvestatus(EN_Project *pr) -/* -**----------------------------------------------------------------- -** Input: none -** Output: returns 1 if any pressure or flow control valve -** changes status, 0 otherwise -** Purpose: updates status for PRVs & PSVs whose status -** is not fixed to OPEN/CLOSED -**----------------------------------------------------------------- -*/ -{ - int change = FALSE, /* Status change flag */ - i,k, /* Valve & link indexes */ - n1,n2; /* Start & end nodes */ - StatType status; /* Valve status settings */ - double hset; /* Valve head setting */ - Slink *link; - - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - report_options_t *rep = &pr->report; - - for (i=1; i <= net->Nvalves; i++) /* Examine each valve */ - { - k = net->Valve[i].Link; /* Link index of valve */ - link = &net->Link[k]; - if (hyd->LinkSetting[k] == MISSING) { - continue; /* Valve status fixed */ - } - n1 = link->N1; /* Start & end nodes */ - n2 = link->N2; - status = hyd->LinkStatus[k]; /* Save current status */ - -// if (s != CLOSED /* No change if flow is */ -// && ABS(hyd->LinkFlows[k]) < hyd->Qtol) continue; /* negligible. */ - - switch (link->Type) /* Evaluate new status: */ - { - case EN_PRV: - hset = net->Node[n2].El + hyd->LinkSetting[k]; - hyd->LinkStatus[k] = prvstatus(pr,k,status,hset,hyd->NodeHead[n1],hyd->NodeHead[n2]); - break; - case EN_PSV: - hset = net->Node[n1].El + hyd->LinkSetting[k]; - hyd->LinkStatus[k] = psvstatus(pr,k,status,hset,hyd->NodeHead[n1],hyd->NodeHead[n2]); - break; - -//// FCV status checks moved back into the linkstatus() function //// -// case FCV: S[k] = fcvstatus(k,s,hyd->NodeHead[n1],hyd->NodeHead[n2]); -// break; - - default: - continue; - } - -/*** Updated 9/7/00 ***/ - /* Do not reset flow in valve if its status changes. */ - /* This strategy improves convergence. */ - - /* Check for status change */ - if (status != hyd->LinkStatus[k]) - { - if (rep->Statflag == FULL) { - writestatchange(pr, k,status,hyd->LinkStatus[k]); - } - change = TRUE; - } - } - return(change); -} /* End of valvestatus() */ - - -/*** Updated 9/7/00 ***/ -int linkstatus(EN_Project *pr) -/* -**-------------------------------------------------------------- -** Input: none -** Output: returns 1 if any link changes status, 0 otherwise -** Purpose: determines new status for pumps, CVs, FCVs & pipes -** to tanks. -**-------------------------------------------------------------- -*/ -{ - int change = FALSE, /* Status change flag */ - k, /* Link index */ - n1, /* Start node index */ - n2; /* End node index */ - double dh; /* Head difference */ - StatType status; /* Current status */ - - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - report_options_t *rep = &pr->report; - - /* Examine each Slink */ - for (k=1; k <= net->Nlinks; k++) - { - Slink *link = &net->Link[k]; - n1 = link->N1; - n2 = link->N2; - dh = hyd->NodeHead[n1] - hyd->NodeHead[n2]; - - /* Re-open temporarily closed links (status = XHEAD or TEMPCLOSED) */ - status = hyd->LinkStatus[k]; - if (status == XHEAD || status == TEMPCLOSED) { - hyd->LinkStatus[k] = OPEN; - } - - /* Check for status changes in CVs and pumps */ - if (link->Type == EN_CVPIPE) { - hyd->LinkStatus[k] = cvstatus(pr,hyd->LinkStatus[k],dh,hyd->LinkFlows[k]); - } - if (link->Type == EN_PUMP && hyd->LinkStatus[k] >= OPEN && hyd->LinkSetting[k] > 0.0) { - hyd->LinkStatus[k] = pumpstatus(pr,k,-dh); - } - - /* Check for status changes in non-fixed FCVs */ - if (link->Type == EN_FCV && hyd->LinkSetting[k] != MISSING) { // - hyd->LinkStatus[k] = fcvstatus(pr,k,status,hyd->NodeHead[n1],hyd->NodeHead[n2]); // - } - - /* Check for flow into (out of) full (empty) tanks */ - if (n1 > net->Njuncs || n2 > net->Njuncs) { - tankstatus(pr,k,n1,n2); - } - - /* Note change in link status; do not revise link flow */ - if (status != hyd->LinkStatus[k]) { - change = TRUE; - if (rep->Statflag == FULL) { - writestatchange(pr,k,status,hyd->LinkStatus[k]); - } - - //if (S[k] <= CLOSED) hyd->LinkFlows[k] = QZERO; - //else setlinkflow(k, dh); - } - } - return(change); -} /* End of linkstatus */ - - -StatType cvstatus(EN_Project *pr, StatType s, double dh, double q) -/* -**-------------------------------------------------- -** Input: s = current status -** dh = headloss -** q = flow -** Output: returns new link status -** Purpose: updates status of a check valve. -**-------------------------------------------------- -*/ -{ - hydraulics_t *hyd = &pr->hydraulics; - - /* Prevent reverse flow through CVs */ - if (ABS(dh) > hyd->Htol) - { - if (dh < -hyd->Htol) { - return(CLOSED); - } - else if (q < -hyd->Qtol) { - return(CLOSED); - } - else { - return(OPEN); - } - } - else - { - if (q < -hyd->Qtol) { - return(CLOSED); - } - else { - return(s); - } - } -} - - -StatType pumpstatus(EN_Project *pr, int k, double dh) -/* -**-------------------------------------------------- -** Input: k = link index -** dh = head gain -** Output: returns new pump status -** Purpose: updates status of an open pump. -**-------------------------------------------------- -*/ -{ - int p; - double hmax; - hydraulics_t *hyd = &pr->hydraulics; - EN_Network *net = &pr->network; - - /* Prevent reverse flow through pump */ - p = findpump(net,k); - if (net->Pump[p].Ptype == CONST_HP) { - hmax = BIG; - } - else { - hmax = SQR(hyd->LinkSetting[k]) * net->Pump[p].Hmax; - } - if (dh > hmax + hyd->Htol) { - return(XHEAD); - } - -/*** Flow higher than pump curve no longer results in a status change ***/ - /* Check if pump cannot deliver flow */ - //if (hyd->LinkFlows[k] > K[k]*Pump[p].Qmax + hyd->Qtol) return(XFLOW); - return(OPEN); -} - - -StatType prvstatus(EN_Project *pr, int k, StatType s, double hset, double h1, double h2) -/* -**----------------------------------------------------------- -** Input: k = link index -** s = current status -** hset = valve head setting -** h1 = head at upstream node -** h2 = head at downstream node -** Output: returns new valve status -** Purpose: updates status of a pressure reducing valve. -**----------------------------------------------------------- -*/ -{ - StatType status; /* New valve status */ - double hml; /* Minor headloss */ - hydraulics_t *hyd = &pr->hydraulics; - - double htol = hyd->Htol; - Slink *link = &pr->network.Link[k]; - - status = s; - if (hyd->LinkSetting[k] == MISSING) - return(status); /* Status fixed by user */ - hml = link->Km*SQR(hyd->LinkFlows[k]); /* Head loss when open */ - - /*** Status rules below have changed. ***/ - - switch (s) - { - case ACTIVE: - if (hyd->LinkFlows[k] < -hyd->Qtol) { - status = CLOSED; - } - else if (h1-hml < hset-htol) { - status = OPEN; - } - else { - status = ACTIVE; - } - break; - case OPEN: - if (hyd->LinkFlows[k] < -hyd->Qtol) { - status = CLOSED; - } - else if (h2 >= hset+htol) { - status = ACTIVE; - } - else { - status = OPEN; - } - break; - case CLOSED: - if ( h1 >= hset+htol && h2 < hset-htol) { - status = ACTIVE; - } - else if (h1 < hset-htol && h1 > h2+htol) { - status = OPEN; - } - else { - status = CLOSED; - } - break; - case XPRESSURE: - if (hyd->LinkFlows[k] < -hyd->Qtol) { - status = CLOSED; - } - break; - default: - break; - } - return(status); -} - - -StatType psvstatus(EN_Project *pr, int k, StatType s, double hset, double h1, double h2) -/* -**----------------------------------------------------------- -** Input: k = link index -** s = current status -** hset = valve head setting -** h1 = head at upstream node -** h2 = head at downstream node -** Output: returns new valve status -** Purpose: updates status of a pressure sustaining valve. -**----------------------------------------------------------- -*/ -{ - StatType status; /* New valve status */ - double hml; /* Minor headloss */ - hydraulics_t *hyd = &pr->hydraulics; - - double htol = hyd->Htol; - - Slink *link = &pr->network.Link[k]; - - status = s; - if (hyd->LinkSetting[k] == MISSING) { - return(status); /* Status fixed by user */ - } - hml = link->Km*SQR(hyd->LinkFlows[k]); /* Head loss when open */ - - /*** Status rules below have changed. ***/ - - switch (s) - { - case ACTIVE: - if (hyd->LinkFlows[k] < -hyd->Qtol) { - status = CLOSED; - } - else if (h2+hml > hset+htol) { - status = OPEN; - } - else { - status = ACTIVE; - } - break; - case OPEN: - if (hyd->LinkFlows[k] < -hyd->Qtol) { - status = CLOSED; - } - else if (h1 < hset-htol) { - status = ACTIVE; - } - else { - status = OPEN; - } - break; - case CLOSED: - if (h2 > hset+htol && h1 > h2+htol) { - status = OPEN; - } - else if (h1 >= hset+htol && h1 > h2+htol) { - status = ACTIVE; - } - else { - status = CLOSED; - } - break; - case XPRESSURE: - if (hyd->LinkFlows[k] < -hyd->Qtol) { - status = CLOSED; - } - break; - default: - break; - } - return(status); -} - - -StatType fcvstatus(EN_Project *pr, int k, StatType s, double h1, double h2) -/* -**----------------------------------------------------------- -** Input: k = link index -** s = current status -** h1 = head at upstream node -** h2 = head at downstream node -** Output: returns new valve status -** Purpose: updates status of a flow control valve. -** -** Valve status changes to XFCV if flow reversal. -** If current status is XFCV and current flow is -** above setting, then valve becomes active. -** If current status is XFCV, and current flow -** positive but still below valve setting, then -** status remains same. -**----------------------------------------------------------- -*/ -{ - StatType status; /* New valve status */ - hydraulics_t *hyd = &pr->hydraulics; - - status = s; - if (h1 - h2 < -hyd->Htol) { - status = XFCV; - } - else if ( hyd->LinkFlows[k] < -hyd->Qtol ) { - status = XFCV; - } - else if (s == XFCV && hyd->LinkFlows[k] >= hyd->LinkSetting[k]) { - status = ACTIVE; - } - return(status); -} - - -/*** Updated 9/7/00 ***/ -/*** Updated 11/19/01 ***/ -void tankstatus(EN_Project *pr, int k, int n1, int n2) -/* -**---------------------------------------------------------------- -** Input: k = link index -** n1 = start node of link -** n2 = end node of link -** Output: none -** Purpose: closes link flowing into full or out of empty tank -**---------------------------------------------------------------- -*/ -{ - int i,n; - double h,q; - Stank *tank; - - hydraulics_t *hyd = &pr->hydraulics; - EN_Network *net = &pr->network; - Slink *link = &net->Link[k]; - - /* Make node n1 be the tank */ - q = hyd->LinkFlows[k]; - i = n1 - net->Njuncs; - if (i <= 0) { - i = n2 - net->Njuncs; - if (i <= 0) { - return; - } - n = n1; - n1 = n2; - n2 = n; - q = -q; - } - h = hyd->NodeHead[n1] - hyd->NodeHead[n2]; - tank = &net->Tank[i]; - /* Skip reservoirs & closed links */ - if (tank->A == 0.0 || hyd->LinkStatus[k] <= CLOSED) { - return; - } - - /* If tank full, then prevent flow into it */ - if (hyd->NodeHead[n1] >= tank->Hmax - hyd->Htol) - { - - /* Case 1: Link is a pump discharging into tank */ - if ( link->Type == EN_PUMP ) { - if (link->N2 == n1) - hyd->LinkStatus[k] = TEMPCLOSED; - } - - /* Case 2: Downstream head > tank head */ - /* (i.e., an open outflow check valve would close) */ - else if (cvstatus(pr, OPEN, h, q) == CLOSED) { - hyd->LinkStatus[k] = TEMPCLOSED; - } - } - - /* If tank empty, then prevent flow out of it */ - if (hyd->NodeHead[n1] <= tank->Hmin + hyd->Htol) { - - /* Case 1: Link is a pump discharging from tank */ - if ( link->Type == EN_PUMP) { - if (link->N1 == n1) { - hyd->LinkStatus[k] = TEMPCLOSED; - } - } - - /* Case 2: Tank head > downstream head */ - /* (i.e., a closed outflow check valve would open) */ - else if (cvstatus(pr, CLOSED, h, q) == OPEN) { - hyd->LinkStatus[k] = TEMPCLOSED; - } - } -} /* End of tankstatus */ - - -int pswitch(EN_Project *pr) -/* -**-------------------------------------------------------------- -** Input: none -** Output: returns 1 if status of any link changes, 0 if not -** Purpose: adjusts settings of links controlled by junction -** pressures after a hydraulic solution is found -**-------------------------------------------------------------- -*/ -{ - int i, /* Control statement index */ - k, /* Link being controlled */ - n, /* Node controlling Slink */ - reset, /* Flag on control conditions */ - change, /* Flag for status or setting change */ - anychange = 0; /* Flag for 1 or more changes */ - char s; /* Current link status */ - - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - report_options_t *rep = &pr->report; - - /* Check each control statement */ - for (i=1; i <= net->Ncontrols; i++) - { - reset = 0; - if ( (k = net->Control[i].Link) <= 0) { - continue; - } - - /* Determine if control based on a junction, not a tank */ - if ( (n = net->Control[i].Node) > 0 && n <= net->Njuncs) { - /* Determine if control conditions are satisfied */ - if (net->Control[i].Type == LOWLEVEL - && hyd->NodeHead[n] <= net->Control[i].Grade + hyd->Htol ) { - reset = 1; - } - if (net->Control[i].Type == HILEVEL - && hyd->NodeHead[n] >= net->Control[i].Grade - hyd->Htol ) { - reset = 1; - } - } - - /* Determine if control forces a status or setting change */ - if (reset == 1) - { - Slink *link = &net->Link[k]; - change = 0; - s = hyd->LinkStatus[k]; - if (link->Type == EN_PIPE) { - if (s != net->Control[i].Status) { - change = 1; - } - } - if (link->Type == EN_PUMP) { - if (hyd->LinkSetting[k] != net->Control[i].Setting) { - change = 1; - } - } - if (link->Type >= EN_PRV) { - if (hyd->LinkSetting[k] != net->Control[i].Setting) { - change = 1; - } - else if (hyd->LinkSetting[k] == MISSING && s != net->Control[i].Status) { - change = 1; - } - } - - /* If a change occurs, update status & setting */ - if (change) { - hyd->LinkStatus[k] = net->Control[i].Status; - if (link->Type > EN_PIPE) { - hyd->LinkSetting[k] = net->Control[i].Setting; - } - if (rep->Statflag == FULL) { - writestatchange(pr, k,s,hyd->LinkStatus[k]); - } - - /* Re-set flow if status has changed */ -// if (S[k] != s) initlinkflow(k, S[k], K[k]); - anychange = 1; - } - } - } - return(anychange); -} /* End of pswitch */ - - -double newflows(EN_Project *pr) -/* -**---------------------------------------------------------------- -** Input: none -** Output: returns solution convergence error -** Purpose: updates link flows after new nodal heads computed -**---------------------------------------------------------------- -*/ -{ - double dh, /* Link head loss */ - dq; /* Link flow change */ - double dqsum, /* Network flow change */ - qsum; /* Network total flow */ - int k, n, n1, n2; - - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - solver_t *s = &hyd->solver; - - /* Initialize net inflows (i.e., demands) at tanks */ - for (n = net->Njuncs + 1; n <= net->Nnodes; n++) { - hyd->NodeDemand[n] = 0.0; - } - - /* Initialize sum of flows & corrections */ - qsum = 0.0; - dqsum = 0.0; - - /* Update flows in all links */ - for (k=1; k<=net->Nlinks; k++) - { - Slink *link = &net->Link[k]; - /* - ** Apply flow update formula: - ** dq = Y - P*(new head loss) - ** P = 1/(dh/dq) - ** Y = P*(head loss based on current flow) - ** where P & Y were computed in newcoeffs(). - */ - - n1 = link->N1; - n2 = link->N2; - dh = hyd->NodeHead[n1] - hyd->NodeHead[n2]; - dq = s->Y[k] - s->P[k] * dh; - - /* Adjust flow change by the relaxation factor */ - dq *= hyd->RelaxFactor; - - /* Prevent flow in constant HP pumps from going negative */ - if (link->Type == EN_PUMP) { - n = findpump(net,k); - if (net->Pump[n].Ptype == CONST_HP && dq > hyd->LinkFlows[k]) { - dq = hyd->LinkFlows[k]/2.0; - } - } - hyd->LinkFlows[k] -= dq; - - /* Update sum of absolute flows & flow corrections */ - qsum += ABS(hyd->LinkFlows[k]); - dqsum += ABS(dq); - - /* Update net flows to tanks */ - if ( hyd->LinkStatus[k] > CLOSED ) - { - if (n1 > net->Njuncs) { - hyd->NodeDemand[n1] -= hyd->LinkFlows[k]; - } - if (n2 > net->Njuncs) { - hyd->NodeDemand[n2] += hyd->LinkFlows[k]; - } - } - - } - - /* Update emitter flows */ - for (k=1; k<=net->Njuncs; k++) - { - if (net->Node[k].Ke == 0.0) { - continue; - } - dq = emitflowchange(pr,k); - hyd->EmitterFlows[k] -= dq; - qsum += ABS(hyd->EmitterFlows[k]); - dqsum += ABS(dq); - } - - /* Return ratio of total flow corrections to total flow */ - if (qsum > hyd->Hacc) { - return(dqsum/qsum); - } - else { - return(dqsum); - } -} /* End of newflows */ - - -void newcoeffs(EN_Project *pr) -/* -**-------------------------------------------------------------- -** Input: none -** Output: none -** Purpose: computes coefficients of linearized network eqns. -**-------------------------------------------------------------- -*/ -{ - hydraulics_t *hyd = &pr->hydraulics; - solver_t *s = &hyd->solver; - EN_Network *n = &pr->network; - - memset(s->Aii,0,(n->Nnodes+1)*sizeof(double)); /* Reset coeffs. to 0 */ - memset(s->Aij,0,(hyd->Ncoeffs+1)*sizeof(double)); - memset(s->F,0,(n->Nnodes+1)*sizeof(double)); - memset(hyd->X_tmp,0,(n->Nnodes+1)*sizeof(double)); - memset(s->P,0,(n->Nlinks+1)*sizeof(double)); - memset(s->Y,0,(n->Nlinks+1)*sizeof(double)); - linkcoeffs(pr); /* Compute link coeffs. */ - emittercoeffs(pr); /* Compute emitter coeffs.*/ - nodecoeffs(pr); /* Compute node coeffs. */ - valvecoeffs(pr); /* Compute valve coeffs. */ -} /* End of newcoeffs */ - - -void linkcoeffs(EN_Project *pr) -/* -**-------------------------------------------------------------- -** Input: none -** Output: none -** Purpose: computes solution matrix coefficients for links -**-------------------------------------------------------------- -*/ -{ - int k,n1,n2; - - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - solver_t *s = &hyd->solver; - EN_Network *n = &pr->network; - - /* Examine each link of network */ - for (k=1; k<=net->Nlinks; k++) - { - Slink *link = &net->Link[k]; - n1 = link->N1; /* Start node of Slink */ - n2 = link->N2; /* End node of link */ - - /* Compute P[k] = 1 / (dh/dQ) and Y[k] = h * P[k] */ - /* for each link k (where h = link head loss). */ - /* FCVs, PRVs, and PSVs with non-fixed status */ - /* are analyzed later. */ - - switch (link->Type) - { - case EN_CVPIPE: - case EN_PIPE: - pipecoeff(pr, k); - break; - case EN_PUMP: - pumpcoeff(pr, k); - break; - case EN_PBV: - pbvcoeff(pr, k); - break; - case EN_TCV: - tcvcoeff(pr, k); - break; - case EN_GPV: - gpvcoeff(pr, k); - break; - case EN_FCV: - case EN_PRV: - case EN_PSV: /* If valve status fixed then treat as pipe, otherwise ignore the valve for now. */ - if (hyd->LinkSetting[k] == MISSING) { - valvecoeff(pr, k); //pipecoeff(k); - } - else { - continue; - } - break; - default: - continue; - } - - /* Update net nodal inflows (X), solution matrix (A) and RHS array (F) */ - /* (Use covention that flow out of node is (-), flow into node is (+)) */ - hyd->X_tmp[n1] -= hyd->LinkFlows[k]; - hyd->X_tmp[n2] += hyd->LinkFlows[k]; - s->Aij[s->Ndx[k]] -= s->P[k]; /* Off-diagonal coeff. */ - if (n1 <= n->Njuncs) /* Node n1 is junction */ - { - s->Aii[s->Row[n1]] += s->P[k]; /* Diagonal coeff. */ - s->F[s->Row[n1]] += s->Y[k]; /* RHS coeff. */ - } - else { - s->F[s->Row[n2]] += (s->P[k]*hyd->NodeHead[n1]); /* Node n1 is a tank */ - } - if (n2 <= n->Njuncs) { /* Node n2 is junction */ - s->Aii[s->Row[n2]] += s->P[k]; /* Diagonal coeff. */ - s->F[s->Row[n2]] -= s->Y[k]; /* RHS coeff. */ - } - else { - s->F[s->Row[n1]] += (s->P[k] * hyd->NodeHead[n2]); /* Node n2 is a tank */ - } - } -} /* End of linkcoeffs */ - - -void nodecoeffs(EN_Project *pr) -/* -**---------------------------------------------------------------- -** Input: none -** Output: none -** Purpose: completes calculation of nodal flow imbalance (X) -** & flow correction (F) arrays -**---------------------------------------------------------------- -*/ -{ - int i; - hydraulics_t *hyd = &pr->hydraulics; - solver_t *s = &hyd->solver; - EN_Network *n = &pr->network; - - /* For junction nodes, subtract demand flow from net */ - /* flow imbalance & add imbalance to RHS array F. */ - for (i=1; i <= n->Njuncs; i++) - { - hyd->X_tmp[i] -= hyd->NodeDemand[i]; - s->F[s->Row[i]] += hyd->X_tmp[i]; - } -} /* End of nodecoeffs */ - - -void valvecoeffs(EN_Project *pr) -/* -**-------------------------------------------------------------- -** Input: none -** Output: none -** Purpose: computes matrix coeffs. for PRVs, PSVs & FCVs -** whose status is not fixed to OPEN/CLOSED -**-------------------------------------------------------------- -*/ -{ - int i,k,n1,n2; - - hydraulics_t *hyd = &pr->hydraulics; - EN_Network *n = &pr->network; - Slink *link; - Svalve *valve; - - for (i=1; i<=n->Nvalves; i++) /* Examine each valve */ - { - valve = &n->Valve[i]; - k = valve->Link; /* Link index of valve */ - link = &n->Link[k]; - if (hyd->LinkSetting[k] == MISSING) { - continue; /* Valve status fixed */ - } - n1 = link->N1; /* Start & end nodes */ - n2 = link->N2; - switch (link->Type) /* Call valve-specific */ - { /* function */ - case EN_PRV: - prvcoeff(pr, k,n1,n2); - break; - case EN_PSV: - psvcoeff(pr, k,n1,n2); - break; - case EN_FCV: - fcvcoeff(pr, k,n1,n2); - break; - default: continue; - } - } -} /* End of valvecoeffs */ - - -void emittercoeffs(EN_Project *pr) -/* -**-------------------------------------------------------------- -** Input: none -** Output: none -** Purpose: computes matrix coeffs. for emitters -** -** Note: Emitters consist of a fictitious pipe connected to -** a fictitious reservoir whose elevation equals that -** of the junction. The headloss through this pipe is -** Ke*(Flow)^hyd->Qexp, where Ke = emitter headloss coeff. -**-------------------------------------------------------------- -*/ -{ - int i; - double ke; - double p; - double q; - double y; - double z; - - hydraulics_t *hyd = &pr->hydraulics; - solver_t *s = &hyd->solver; - EN_Network *n = &pr->network; - - for (i=1; i <= n->Njuncs; i++) - { - Snode *node = &n->Node[i]; - if (node->Ke == 0.0) { - continue; - } - ke = MAX(CSMALL, node->Ke); - q = hyd->EmitterFlows[i]; - z = ke * pow(ABS(q),hyd->Qexp); - p = hyd->Qexp * z / ABS(q); - if (p < hyd->RQtol) { - p = 1.0 / hyd->RQtol; - } - else { - p = 1.0 / p; - } - y = SGN(q)*z*p; - s->Aii[s->Row[i]] += p; - s->F[s->Row[i]] += y + p * node->El; - hyd->X_tmp[i] -= q; - } -} - - -double emitflowchange(EN_Project *pr, int i) -/* -**-------------------------------------------------------------- -** Input: i = node index -** Output: returns change in flow at an emitter node -** Purpose: computes flow change at an emitter node -**-------------------------------------------------------------- -*/ -{ - double ke, p; - hydraulics_t *hyd = &pr->hydraulics; - EN_Network *n = &pr->network; - Snode *node = &n->Node[i]; - - ke = MAX(CSMALL, node->Ke); - p = hyd->Qexp * ke * pow(ABS(hyd->EmitterFlows[i]),(hyd->Qexp - 1.0)); - if (p < hyd->RQtol) { - p = 1/hyd->RQtol; - } - else { - p = 1.0/p; - } - return(hyd->EmitterFlows[i] / hyd->Qexp - p * (hyd->NodeHead[i] - node->El)); -} - - -void pipecoeff(EN_Project *pr, int k) -/* -**-------------------------------------------------------------- -** Input: k = link index -** Output: none -** Purpose: computes P & Y coefficients for pipe k -** -** P = inverse head loss gradient = 1/(dh/dQ) -** Y = flow correction term = h*P -**-------------------------------------------------------------- -*/ -{ - double hpipe, /* Normal head loss */ - hml, /* Minor head loss */ - ml, /* Minor loss coeff. */ - p, /* q*(dh/dq) */ - q, /* Abs. value of flow */ - r, /* Resistance coeff. */ - r1, /* Total resistance factor */ - f, /* D-W friction factor */ - dfdq; /* Derivative of fric. fact. */ - - hydraulics_t *hyd = &pr->hydraulics; - solver_t *s = &hyd->solver; - Slink *link = &pr->network.Link[k]; - - /* For closed pipe use headloss formula: h = CBIG*q */ - if (hyd->LinkStatus[k] <= CLOSED) - { - s->P[k] = 1.0/CBIG; - s->Y[k] = hyd->LinkFlows[k]; - return; - } - - /* Evaluate headloss coefficients */ - q = ABS(hyd->LinkFlows[k]); /* Absolute flow */ - ml = link->Km; /* Minor loss coeff. */ - r = link->R; /* Resistance coeff. */ - f = 1.0; /* D-W friction factor */ - if (hyd->Formflag == DW) f = DWcoeff(pr,k,&dfdq); - r1 = f*r+ml; - - /* Use large P coefficient for small flow resistance product */ - if (r1*q < hyd->RQtol) - { - s->P[k] = 1.0/hyd->RQtol; - s->Y[k] = hyd->LinkFlows[k]/hyd->Hexp; - return; - } - - /* Compute P and Y coefficients */ - if (hyd->Formflag == DW) /* D-W eqn. */ - { - hpipe = r1*SQR(q); /* Total head loss */ - p = 2.0*r1*q; /* |dh/dQ| */ - /* + dfdq*r*q*q;*/ /* Ignore df/dQ term */ - p = 1.0/p; - s->P[k] = p; - s->Y[k] = SGN(hyd->LinkFlows[k])*hpipe*p; - } - else /* H-W or C-M eqn. */ - { - hpipe = r*pow(q,hyd->Hexp); /* Friction head loss */ - p = hyd->Hexp*hpipe; /* Q*dh(friction)/dQ */ - if (ml > 0.0) - { - hml = ml*q*q; /* Minor head loss */ - p += 2.0*hml; /* Q*dh(Total)/dQ */ - } - else hml = 0.0; - p = hyd->LinkFlows[k]/p; /* 1 / (dh/dQ) */ - s->P[k] = ABS(p); - s->Y[k] = p*(hpipe + hml); - } -} /* End of pipecoeff */ - - -double DWcoeff(EN_Project *pr, int k, double *dfdq) -/* -**-------------------------------------------------------------- -** Input: k = link index -** Output: returns Darcy-Weisbach friction factor -** Purpose: computes Darcy-Weisbach friction factor -** -** Uses interpolating polynomials developed by -** E. Dunlop for transition flow from 2000 < Re < 4000. -** -** df/dq term is ignored as it slows convergence rate. -**-------------------------------------------------------------- -*/ -{ - double q, /* Abs. value of flow */ - f; /* Friction factor */ - double x1,x2,x3,x4, - y1,y2,y3, - fa,fb,r; - double s,w; - - hydraulics_t *hyd = &pr->hydraulics; - Slink *link = &pr->network.Link[k]; - - *dfdq = 0.0; - if (link->Type > EN_PIPE) - return(1.0); /* Only apply to pipes */ - q = ABS(hyd->LinkFlows[k]); - s = hyd->Viscos * link->Diam; - w = q/s; /* w = Re(Pi/4) */ - if (w >= A1) /* Re >= 4000; Colebrook Formula */ - { - y1 = A8/pow(w,0.9); - y2 = link->Kc/(3.7*link->Diam) + y1; - y3 = A9*log(y2); - f = 1.0/SQR(y3); - /* *dfdq = (2.0+AA*y1/(y2*y3))*f; */ /* df/dq */ - } - else if (w > A2) /* Re > 2000; Interpolation formula */ - { - y2 = link->Kc/(3.7*link->Diam) + AB; - y3 = A9*log(y2); - fa = 1.0/SQR(y3); - fb = (2.0+AC/(y2*y3))*fa; - r = w/A2; - x1 = 7.0*fa - fb; - x2 = 0.128 - 17.0*fa + 2.5*fb; - x3 = -0.128 + 13.0*fa - (fb+fb); - x4 = r*(0.032 - 3.0*fa + 0.5*fb); - f = x1 + r*(x2 + r*(x3+x4)); - /* *dfdq = (x1 + x1 + r*(3.0*x2 + r*(4.0*x3 + 5.0*x4))); */ - } - else if (w > A4) /* Laminar flow: Hagen-Poiseuille Formula */ - { - f = A3*s/q; - /* *dfdq = A3*s; */ - } - else - { - f = 8.0; - *dfdq = 0.0; - } - return(f); -} /* End of DWcoeff */ - - -/*** Updated 10/25/00 ***/ -/*** Updated 12/29/00 ***/ -void pumpcoeff(EN_Project *pr, int k) -/* -**-------------------------------------------------------------- -** Input: k = link index -** Output: none -** Purpose: computes P & Y coeffs. for pump in link k -**-------------------------------------------------------------- -*/ -{ - int p; /* Pump index */ - double h0, /* Shutoff head */ - q, /* Abs. value of flow */ - r, /* Flow resistance coeff. */ - n; /* Flow exponent coeff. */ - hydraulics_t *hyd; - solver_t *s; - double setting; - Spump *pump; - - hyd = &pr->hydraulics; - s = &hyd->solver; - - setting = hyd->LinkSetting[k]; - - /* Use high resistance pipe if pump closed or cannot deliver head */ - if (hyd->LinkStatus[k] <= CLOSED || setting == 0.0) { - s->P[k] = 1.0/CBIG; - s->Y[k] = hyd->LinkFlows[k]; - return; - } - - q = ABS(hyd->LinkFlows[k]); - q = MAX(q,TINY); - p = findpump(&pr->network,k); - - pump = &pr->network.Pump[p]; - /* Get pump curve coefficients for custom pump curve. */ - if (pump->Ptype == CUSTOM) - { - /* Find intercept (h0) & slope (r) of pump curve */ - /* line segment which contains speed-adjusted flow. */ - curvecoeff(pr,pump->Hcurve, q/setting, &h0, &r); - - /* Determine head loss coefficients. */ - pump->H0 = -h0; - pump->R = -r; - pump->N = 1.0; - } - - /* Adjust head loss coefficients for pump speed. */ - h0 = SQR(setting) * pump->H0; - n = pump->N; - r = pump->R * pow(setting,2.0-n); - if (n != 1.0) { - r = n * r * pow(q,n-1.0); - } - - /* Compute inverse headloss gradient (P) and flow correction factor (Y) */ - s->P[k] = 1.0/MAX(r,hyd->RQtol); - s->Y[k] = hyd->LinkFlows[k]/n + s->P[k]*h0; -} /* End of pumpcoeff */ - - -/*** Updated 10/25/00 ***/ -/*** Updated 12/29/00 ***/ -void curvecoeff(EN_Project *p, int i, double q, double *h0, double *r) -/* -**------------------------------------------------------------------- -** Input: i = curve index -** q = flow rate -** Output: *h0 = head at zero flow (y-intercept) -** *r = dHead/dFlow (slope) -** Purpose: computes intercept and slope of head v. flow curve -** at current flow. -**------------------------------------------------------------------- -*/ -{ - int k1, k2, npts; - double *x, *y; - Scurve *curve; - - /* Remember that curve is stored in untransformed units */ - q *= p->Ucf[FLOW]; - curve = &p->network.Curve[i]; - x = curve->X; /* x = flow */ - y = curve->Y; /* y = head */ - npts = curve->Npts; - - /* Find linear segment of curve that brackets flow q */ - k2 = 0; - while (k2 < npts && x[k2] < q) - k2++; - - if (k2 == 0) - k2++; - - else if (k2 == npts) - k2--; - - k1 = k2 - 1; - - /* Compute slope and intercept of this segment */ - *r = (y[k2]-y[k1]) / (x[k2]-x[k1]); - *h0 = y[k1] - (*r)*x[k1]; - - /* Convert units */ - *h0 = (*h0) / p->Ucf[HEAD]; - *r = (*r)*p->Ucf[FLOW]/p->Ucf[HEAD]; -} /* End of curvecoeff */ - - -void gpvcoeff(EN_Project *p, int k) -/* -**-------------------------------------------------------------- -** Input: k = link index -** Output: none -** Purpose: computes P & Y coeffs. for general purpose valve -**-------------------------------------------------------------- -*/ -{ - double h0, /* Headloss curve intercept */ - q, /* Abs. value of flow */ - r; /* Flow resistance coeff. */ - - hydraulics_t *hyd = &p->hydraulics; - solver_t *s = &hyd->solver; - -/*** Updated 9/7/00 ***/ - /* Treat as a pipe if valve closed */ - if (hyd->LinkStatus[k] == CLOSED) { - valvecoeff(p,k); //pipecoeff(k); - } - /* Otherwise utilize headloss curve */ - /* whose index is stored in K */ - else { - /* Find slope & intercept of headloss curve. */ - q = ABS(hyd->LinkFlows[k]); - q = MAX(q,TINY); - -/*** Updated 10/25/00 ***/ -/*** Updated 12/29/00 ***/ - curvecoeff(p,(int)ROUND(hyd->LinkSetting[k]),q,&h0,&r); - - /* Compute inverse headloss gradient (P) */ - /* and flow correction factor (Y). */ - s->P[k] = 1.0 / MAX(r,hyd->RQtol); - s->Y[k] = s->P[k]*(h0 + r*q) * SGN(hyd->LinkFlows[k]); - } -} - - -void pbvcoeff(EN_Project *p, int k) -/* -**-------------------------------------------------------------- -** Input: k = link index -** Output: none -** Purpose: computes P & Y coeffs. for pressure breaker valve -**-------------------------------------------------------------- -*/ -{ - - Slink *link = &p->network.Link[k]; - hydraulics_t *hyd = &p->hydraulics; - solver_t *s = &hyd->solver; - - /* If valve fixed OPEN or CLOSED then treat as a pipe */ - if (hyd->LinkSetting[k] == MISSING || hyd->LinkSetting[k] == 0.0) { - valvecoeff(p,k); //pipecoeff(k); - } - - /* If valve is active */ - else { - /* Treat as a pipe if minor loss > valve setting */ - if (link->Km * SQR(hyd->LinkFlows[k]) > hyd->LinkSetting[k]) { - valvecoeff(p,k); //pipecoeff(k); - } - /* Otherwise force headloss across valve to be equal to setting */ - else { - s->P[k] = CBIG; - s->Y[k] = hyd->LinkSetting[k] * CBIG; - } - } -} /* End of pbvcoeff */ - - -void tcvcoeff(EN_Project *p, int k) -/* -**-------------------------------------------------------------- -** Input: k = link index -** Output: none -** Purpose: computes P & Y coeffs. for throttle control valve -**-------------------------------------------------------------- -*/ -{ - double km; - Slink *link = &p->network.Link[k]; - hydraulics_t *hyd = &p->hydraulics; - - /* Save original loss coeff. for open valve */ - km = link->Km; - - /* If valve not fixed OPEN or CLOSED, compute its loss coeff. */ - if (hyd->LinkSetting[k] != MISSING) { - link->Km = 0.02517 * hyd->LinkSetting[k] / (SQR(link->Diam)*SQR(link->Diam)); - } - - /* Then apply usual pipe formulas */ - valvecoeff(p,k); //pipecoeff(k); - - /* Restore original loss coeff. */ - link->Km = km; -} /* End of tcvcoeff */ - - -void prvcoeff(EN_Project *p, int k, int n1, int n2) -/* -**-------------------------------------------------------------- -** Input: k = link index -** n1 = upstream node of valve -** n2 = downstream node of valve -** Output: none -** Purpose: computes solution matrix coeffs. for pressure -** reducing valves -**-------------------------------------------------------------- -*/ -{ - int i,j; /* Rows of solution matrix */ - double hset; /* Valve head setting */ - hydraulics_t *hyd = &p->hydraulics; - solver_t *s = &hyd->solver; - - i = s->Row[n1]; /* Matrix rows of nodes */ - j = s->Row[n2]; - hset = p->network.Node[n2].El + hyd->LinkSetting[k]; /* Valve setting */ - - if (hyd->LinkStatus[k] == ACTIVE) - { - /* - Set coeffs. to force head at downstream - node equal to valve setting & force flow (when updated in - newflows()) equal to flow imbalance at downstream node. - */ - s->P[k] = 0.0; - s->Y[k] = hyd->LinkFlows[k] + hyd->X_tmp[n2]; /* Force flow balance */ - s->F[j] += (hset * CBIG); /* Force head = hset */ - s->Aii[j] += CBIG; /* at downstream node */ - if (hyd->X_tmp[n2] < 0.0) { - s->F[i] += hyd->X_tmp[n2]; - } - return; - } - - /* - For OPEN, CLOSED, or XPRESSURE valve - compute matrix coeffs. using the valvecoeff() function. - */ - valvecoeff(p,k); /*pipecoeff(k);*/ - s->Aij[s->Ndx[k]] -= s->P[k]; - s->Aii[i] += s->P[k]; - s->Aii[j] += s->P[k]; - s->F[i] += (s->Y[k] - hyd->LinkFlows[k]); - s->F[j] -= (s->Y[k] - hyd->LinkFlows[k]); -} /* End of prvcoeff */ - - -void psvcoeff(EN_Project *p, int k, int n1, int n2) -/* -**-------------------------------------------------------------- -** Input: k = link index -** n1 = upstream node of valve -** n2 = downstream node of valve -** Output: none -** Purpose: computes solution matrix coeffs. for pressure -** sustaining valve -**-------------------------------------------------------------- -*/ -{ - int i,j; /* Rows of solution matrix */ - double hset; /* Valve head setting */ - - hydraulics_t *hyd = &p->hydraulics; - solver_t *s = &hyd->solver; - - i = s->Row[n1]; /* Matrix rows of nodes */ - j = s->Row[n2]; - hset = p->network.Node[n1].El + hyd->LinkSetting[k]; /* Valve setting */ - - if (hyd->LinkStatus[k] == ACTIVE) - { - /* - Set coeffs. to force head at upstream - node equal to valve setting & force flow (when updated in - newflows()) equal to flow imbalance at upstream node. - */ - s->P[k] = 0.0; - s->Y[k] = hyd->LinkFlows[k] - hyd->X_tmp[n1]; /* Force flow balance */ - s->F[i] += (hset*CBIG); /* Force head = hset */ - s->Aii[i] += CBIG; /* at upstream node */ - if (hyd->X_tmp[n1] > 0.0) { - s->F[j] += hyd->X_tmp[n1]; - } - return; - } - - /* - For OPEN, CLOSED, or XPRESSURE valve - compute matrix coeffs. using the valvecoeff() function. - */ - valvecoeff(p,k); /*pipecoeff(k);*/ - s->Aij[s->Ndx[k]] -= s->P[k]; - s->Aii[i] += s->P[k]; - s->Aii[j] += s->P[k]; - s->F[i] += (s->Y[k] - hyd->LinkFlows[k]); - s->F[j] -= (s->Y[k] - hyd->LinkFlows[k]); -} /* End of psvcoeff */ - - -void fcvcoeff(EN_Project *p, int k, int n1, int n2) -/* -**-------------------------------------------------------------- -** Input: k = link index -** n1 = upstream node of valve -** n2 = downstream node of valve -** Output: none -** Purpose: computes solution matrix coeffs. for flow control -** valve -**-------------------------------------------------------------- -*/ -{ - int i,j; /* Rows in solution matrix */ - double q; /* Valve flow setting */ - - hydraulics_t *hyd = &p->hydraulics; - solver_t *s = &hyd->solver; - - q = hyd->LinkSetting[k]; - i = hyd->solver.Row[n1]; - j = hyd->solver.Row[n2]; - - /* - If valve active, break network at valve and treat - flow setting as external demand at upstream node - and external supply at downstream node. - */ - if (hyd->LinkStatus[k] == ACTIVE) - { - hyd->X_tmp[n1] -= q; - s->F[i] -= q; - hyd->X_tmp[n2] += q; - s->F[j] += q; - /*P[k] = 0.0;*/ - s->P[k] = 1.0/CBIG; - s->Aij[s->Ndx[k]] -= s->P[k]; - s->Aii[i] += s->P[k]; - s->Aii[j] += s->P[k]; - s->Y[k] = hyd->LinkFlows[k] - q; - } - /* - Otherwise treat valve as an open pipe - */ - else - { - valvecoeff(p,k); //pipecoeff(k); - s->Aij[s->Ndx[k]] -= s->P[k]; - s->Aii[i] += s->P[k]; - s->Aii[j] += s->P[k]; - s->F[i] += (s->Y[k] - hyd->LinkFlows[k]); - s->F[j] -= (s->Y[k] - hyd->LinkFlows[k]); - } -} /* End of fcvcoeff */ - - -/*** New function added. ***/ -void valvecoeff(EN_Project *pr, int k) -/* -**-------------------------------------------------------------- -** Input: k = link index -** Output: none -** Purpose: computes solution matrix coeffs. for a completely -** open, closed, or throttled control valve. -**-------------------------------------------------------------- -*/ -{ - double p; - - EN_Network *n = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - solver_t *s = &hyd->solver; - - Slink *link = &n->Link[k]; - - double flow = hyd->LinkFlows[k]; - - // Valve is closed. Use a very small matrix coeff. - if (hyd->LinkStatus[k] <= CLOSED) - { - s->P[k] = 1.0/CBIG; - s->Y[k] = flow; - return; - } - - // Account for any minor headloss through the valve - if (link->Km > 0.0) - { - p = 2.0 * link->Km * fabs(flow); - if ( p < hyd->RQtol ) { - p = hyd->RQtol; - } - s->P[k] = 1.0 / p; - s->Y[k] = flow / 2.0; - } - else - { - s->P[k] = 1.0 / hyd->RQtol; - s->Y[k] = flow; - } -} - /**************** END OF HYDRAUL.C ***************/ - diff --git a/src/hydsolver.c b/src/hydsolver.c new file mode 100644 index 0000000..9f789f1 --- /dev/null +++ b/src/hydsolver.c @@ -0,0 +1,1052 @@ +/* +********************************************************************* + +HYDSOLVER.C -- Equilibrium hydraulic solver for the EPANET Program + +This module contains a pipe network hydraulic solver that computes +flows and pressures within the network at a specific point in time. + +The solver implements Todini's Global Gradient Algorithm. +******************************************************************* +*/ + +#include +#include +#ifndef __APPLE__ +#include +#else +#include +#endif +#include +#include "text.h" +#include "types.h" +#include "funcs.h" + +// Error in network being hydraulically balanced +typedef struct { + double maxheaderror; + double maxflowerror; + double maxflowchange; + int maxheadlink; + int maxflownode; + int maxflowlink; +} Hydbalance; + +// External functions +//int hydsolve(EN_Project *pr, int *iter, double *relerr); + +// Local functions +static int badvalve(EN_Project *pr, int); +static int valvestatus(EN_Project *pr); +static int linkstatus(EN_Project *pr); +static StatType cvstatus(EN_Project *pr, StatType, double, double); +static StatType pumpstatus(EN_Project *pr, int, double); +static StatType prvstatus(EN_Project *pr, int, StatType, double, double, double); +static StatType psvstatus(EN_Project *pr, int, StatType, double, double, double); +static StatType fcvstatus(EN_Project *pr, int, StatType, double, double); +static void tankstatus(EN_Project *pr, int, int, int); +static int pswitch(EN_Project *pr); + +static double newflows(EN_Project *pr, Hydbalance *hbal); +static double emitflowchange(EN_Project *pr, int i); +static void checkhydbalance(EN_Project *pr, Hydbalance *hbal); +static int hasconverged(EN_Project *pr, double *relerr, Hydbalance *hbal); +static void reporthydbal(EN_Project *pr, Hydbalance *hbal); + + +int hydsolve(EN_Project *pr, int *iter, double *relerr) +/* +**------------------------------------------------------------------- +** Input: none +** Output: *iter = # of iterations to reach solution +** *relerr = convergence error in solution +** returns error code +** Purpose: solves network nodal equations for heads and flows +** using Todini's Gradient algorithm +** +*** Updated 9/7/00 *** +*** Updated 2.00.11 *** +*** Updated 2.00.12 *** +** Notes: Status checks on CVs, pumps and pipes to tanks are made +** every hyd->CheckFreq iteration, up until MaxCheck iterations +** are reached. Status checks on control valves are made +** every iteration if DampLimit = 0 or only when the +** convergence error is at or below DampLimit. If DampLimit +** is > 0 then future computed flow changes are only 60% of +** their full value. A complete status check on all links +** is made when convergence is achieved. If convergence is +** not achieved in hyd->MaxIter trials and hyd->ExtraIter > 0 then +** another hyd->ExtraIter trials are made with no status changes +** made to any links and a warning message is generated. +** +** This procedure calls linsolve() which appears in SMATRIX.C. +**------------------------------------------------------------------- +*/ +{ + int i; /* Node index */ + int errcode = 0; /* Node causing solution error */ + int nextcheck; /* Next status check trial */ + int maxtrials; /* Max. trials for convergence */ + double newerr; /* New convergence error */ + int valveChange; /* Valve status change flag */ + int statChange; /* Non-valve status change flag */ + Hydbalance hydbal; /* Hydraulic balance errors */ + + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + report_options_t *rep = &pr->report; + + /* Initialize status checking & relaxation factor */ + nextcheck = hyd->CheckFreq; + hyd->RelaxFactor = 1.0; + + /* Compute initial head loss coefficients*/ + for (i = 1; i <= net->Nlinks; i++) { + hlosscoeff(pr, i); + } + + /* Repeat iterations until convergence or trial limit is exceeded. */ + /* (hyd->ExtraIter used to increase trials in case of status cycling.) */ + if (pr->report.Statflag == FULL) { + writerelerr(pr, 0, 0); + } + maxtrials = hyd->MaxIter; + if (hyd->ExtraIter > 0) { + maxtrials += hyd->ExtraIter; + } + *iter = 1; + while (*iter <= maxtrials) { + /* + ** Compute coefficient matrices A & F and solve A*H = F + ** where H = heads, A = Jacobian coeffs. derived from + ** head loss gradients, & F = flow correction terms. + ** Solution for H is returned in F from call to linsolve(). + */ + matrixcoeffs(pr); + errcode = linsolve(&hyd->solver, net->Njuncs); + + /* Take action depending on error code */ + if (errcode < 0) { + break; /* Memory allocation problem */ + } + if (errcode > 0) { /* Ill-conditioning problem */ + /* If control valve causing problem, fix its status & continue, */ + /* otherwise end the iterations with no solution. */ + if (badvalve(pr, sol->Order[errcode])) { + continue; + } + else break; + } + + /* Update current solution. */ + /* (Row[i] = row of solution matrix corresponding to node i). */ + for (i = 1; i <= net->Njuncs; i++) { + hyd->NodeHead[i] = sol->F[sol->Row[i]]; /* Update heads */ + } + newerr = newflows(pr, &hydbal); /* Update flows */ + *relerr = newerr; + + /* Check hydraulic balance & re-compute head loss coefficients */ + checkhydbalance(pr, &hydbal); + for (i = 1; i <= net->Nlinks; i++) { + hlosscoeff(pr, i); + } + + /* Write convergence error to status report if called for */ + if (rep->Statflag == FULL) { + writerelerr(pr, *iter, *relerr); + reporthydbal(pr, &hydbal); + } + + /* Apply solution damping & check for change in valve status */ + hyd->RelaxFactor = 1.0; + valveChange = FALSE; + if (hyd->DampLimit > 0.0) { + if (*relerr <= hyd->DampLimit) { + hyd->RelaxFactor = 0.6; + valveChange = valvestatus(pr); + } + } + else { + valveChange = valvestatus(pr); + } + + /* Check for convergence */ + if (hasconverged(pr, relerr, &hydbal)) { + /* We have convergence. Quit if we are into extra iterations. */ + if (*iter > hyd->MaxIter) { + break; + } + + /* Quit if no status changes occur. */ + statChange = FALSE; + if (valveChange) { + statChange = TRUE; + } + if (linkstatus(pr)) { + statChange = TRUE; + } + if (pswitch(pr)) { + statChange = TRUE; + } + if (!statChange) { + break; + } + + /* We have a status change so continue the iterations */ + nextcheck = *iter + hyd->CheckFreq; + } + + /* No convergence yet. See if its time for a periodic status */ + /* check on pumps, CV's, and pipes connected to tanks. */ + else if (*iter <= hyd->MaxCheck && *iter == nextcheck) + { + linkstatus(pr); + nextcheck += hyd->CheckFreq; + } + (*iter)++; + } + + /* Iterations ended. Report any errors. */ + if (errcode < 0) errcode = 101; /* Memory allocation error */ + else if (errcode > 0) + { + writehyderr(pr, sol->Order[errcode]); /* Ill-conditioned eqns. error */ + errcode = 110; + } + + /* Add any emitter flows to junction demands */ + for (i = 1; i <= net->Njuncs; i++) { + hyd->NodeDemand[i] += hyd->EmitterFlows[i]; + } + return(errcode); +} /* End of netsolve */ + + +int badvalve(EN_Project *pr, int n) +/* +**----------------------------------------------------------------- +** Input: n = node index +** Output: returns 1 if node n belongs to an active control valve, +** 0 otherwise +** Purpose: determines if a node belongs to an active control valve +** whose setting causes an inconsistent set of eqns. If so, +** the valve status is fixed open and a warning condition +** is generated. +**----------------------------------------------------------------- +*/ +{ + int i, k, n1, n2; + + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + report_options_t *rep = &pr->report; + time_options_t *time = &pr->time_options; + Slink *link; + + for (i = 1; i <= net->Nvalves; i++) { + k = net->Valve[i].Link; + link = &net->Link[k]; + n1 = link->N1; + n2 = link->N2; + if (n == n1 || n == n2) { + EN_LinkType t = link->Type; + if (t == EN_PRV || t == EN_PSV || t == EN_FCV) { + if (hyd->LinkStatus[k] == ACTIVE) { + if (rep->Statflag == FULL) { + sprintf(pr->Msg, FMT61, clocktime(rep->Atime, time->Htime), link->ID); + writeline(pr, pr->Msg); + } + if (link->Type == EN_FCV) { + hyd->LinkStatus[k] = XFCV; + } + else { + hyd->LinkStatus[k] = XPRESSURE; + } + return(1); + } + } + return(0); + } + } + return(0); +} + + +int valvestatus(EN_Project *pr) +/* +**----------------------------------------------------------------- +** Input: none +** Output: returns 1 if any pressure or flow control valve +** changes status, 0 otherwise +** Purpose: updates status for PRVs & PSVs whose status +** is not fixed to OPEN/CLOSED +**----------------------------------------------------------------- +*/ +{ + int change = FALSE, /* Status change flag */ + i, k, /* Valve & link indexes */ + n1, n2; /* Start & end nodes */ + StatType status; /* Valve status settings */ + double hset; /* Valve head setting */ + Slink *link; + + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + report_options_t *rep = &pr->report; + + for (i = 1; i <= net->Nvalves; i++) /* Examine each valve */ + { + k = net->Valve[i].Link; /* Link index of valve */ + link = &net->Link[k]; + if (hyd->LinkSetting[k] == MISSING) { + continue; /* Valve status fixed */ + } + n1 = link->N1; /* Start & end nodes */ + n2 = link->N2; + status = hyd->LinkStatus[k]; /* Save current status */ + +// if (s != CLOSED /* No change if flow is */ +// && ABS(hyd->LinkFlows[k]) < hyd->Qtol) continue; /* negligible. */ + + switch (link->Type) /* Evaluate new status: */ + { + case EN_PRV: + hset = net->Node[n2].El + hyd->LinkSetting[k]; + hyd->LinkStatus[k] = prvstatus(pr, k, status, hset, hyd->NodeHead[n1], hyd->NodeHead[n2]); + break; + case EN_PSV: + hset = net->Node[n1].El + hyd->LinkSetting[k]; + hyd->LinkStatus[k] = psvstatus(pr, k, status, hset, hyd->NodeHead[n1], hyd->NodeHead[n2]); + break; + +//// FCV status checks moved back into the linkstatus() function //// +// case FCV: S[k] = fcvstatus(k,s,hyd->NodeHead[n1],hyd->NodeHead[n2]); +// break; + + default: + continue; + } + + /*** Updated 9/7/00 ***/ + /* Do not reset flow in valve if its status changes. */ + /* This strategy improves convergence. */ + + /* Check for status change */ + if (status != hyd->LinkStatus[k]) + { + if (rep->Statflag == FULL) { + writestatchange(pr, k, status, hyd->LinkStatus[k]); + } + change = TRUE; + } + } + return(change); +} /* End of valvestatus() */ + + +/*** Updated 9/7/00 ***/ +int linkstatus(EN_Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: returns 1 if any link changes status, 0 otherwise +** Purpose: determines new status for pumps, CVs, FCVs & pipes +** to tanks. +**-------------------------------------------------------------- +*/ +{ + int change = FALSE, /* Status change flag */ + k, /* Link index */ + n1, /* Start node index */ + n2; /* End node index */ + double dh; /* Head difference */ + StatType status; /* Current status */ + + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + report_options_t *rep = &pr->report; + Slink *link; + + /* Examine each Slink */ + for (k = 1; k <= net->Nlinks; k++) + { + link = &net->Link[k]; + n1 = link->N1; + n2 = link->N2; + dh = hyd->NodeHead[n1] - hyd->NodeHead[n2]; + + /* Re-open temporarily closed links (status = XHEAD or TEMPCLOSED) */ + status = hyd->LinkStatus[k]; + if (status == XHEAD || status == TEMPCLOSED) { + hyd->LinkStatus[k] = OPEN; + } + + /* Check for status changes in CVs and pumps */ + if (link->Type == EN_CVPIPE) { + hyd->LinkStatus[k] = cvstatus(pr, hyd->LinkStatus[k], dh, hyd->LinkFlows[k]); + } + if (link->Type == EN_PUMP && hyd->LinkStatus[k] >= OPEN && hyd->LinkSetting[k] > 0.0) { + hyd->LinkStatus[k] = pumpstatus(pr, k, -dh); + } + + /* Check for status changes in non-fixed FCVs */ + if (link->Type == EN_FCV && hyd->LinkSetting[k] != MISSING) { + hyd->LinkStatus[k] = fcvstatus(pr, k, status, hyd->NodeHead[n1], hyd->NodeHead[n2]); // + } + + /* Check for flow into (out of) full (empty) tanks */ + if (n1 > net->Njuncs || n2 > net->Njuncs) { + tankstatus(pr, k, n1, n2); + } + + /* Note change in link status; do not revise link flow */ + if (status != hyd->LinkStatus[k]) { + change = TRUE; + if (rep->Statflag == FULL) { + writestatchange(pr, k, status, hyd->LinkStatus[k]); + } + + //if (S[k] <= CLOSED) hyd->LinkFlows[k] = QZERO; + //else setlinkflow(k, dh); + } + } + return(change); +} /* End of linkstatus */ + + +StatType cvstatus(EN_Project *pr, StatType s, double dh, double q) +/* +**-------------------------------------------------- +** Input: s = current status +** dh = headloss +** q = flow +** Output: returns new link status +** Purpose: updates status of a check valve. +**-------------------------------------------------- +*/ +{ + hydraulics_t *hyd = &pr->hydraulics; + + /* Prevent reverse flow through CVs */ + if (ABS(dh) > hyd->Htol) + { + if (dh < -hyd->Htol) { + return(CLOSED); + } + else if (q < -hyd->Qtol) { + return(CLOSED); + } + else { + return(OPEN); + } + } + else + { + if (q < -hyd->Qtol) { + return(CLOSED); + } + else { + return(s); + } + } +} + + +StatType pumpstatus(EN_Project *pr, int k, double dh) +/* +**-------------------------------------------------- +** Input: k = link index +** dh = head gain +** Output: returns new pump status +** Purpose: updates status of an open pump. +**-------------------------------------------------- +*/ +{ + int p; + double hmax; + hydraulics_t *hyd = &pr->hydraulics; + EN_Network *net = &pr->network; + + /* Prevent reverse flow through pump */ + p = findpump(net, k); + if (net->Pump[p].Ptype == CONST_HP) { + hmax = BIG; + } + else { + hmax = SQR(hyd->LinkSetting[k]) * net->Pump[p].Hmax; + } + if (dh > hmax + hyd->Htol) { + return(XHEAD); + } + + /*** Flow higher than pump curve no longer results in a status change ***/ + /* Check if pump cannot deliver flow */ + //if (hyd->LinkFlows[k] > K[k]*Pump[p].Qmax + hyd->Qtol) return(XFLOW); + return(OPEN); +} + + +StatType prvstatus(EN_Project *pr, int k, StatType s, double hset, double h1, double h2) +/* +**----------------------------------------------------------- +** Input: k = link index +** s = current status +** hset = valve head setting +** h1 = head at upstream node +** h2 = head at downstream node +** Output: returns new valve status +** Purpose: updates status of a pressure reducing valve. +**----------------------------------------------------------- +*/ +{ + StatType status; /* New valve status */ + double hml; /* Minor headloss */ + hydraulics_t *hyd = &pr->hydraulics; + + double htol = hyd->Htol; + Slink *link = &pr->network.Link[k]; + + status = s; + if (hyd->LinkSetting[k] == MISSING) + return(status); /* Status fixed by user */ + hml = link->Km*SQR(hyd->LinkFlows[k]); /* Head loss when open */ + + /*** Status rules below have changed. ***/ + + switch (s) + { + case ACTIVE: + if (hyd->LinkFlows[k] < -hyd->Qtol) { + status = CLOSED; + } + else if (h1 - hml < hset - htol) { + status = OPEN; + } + else { + status = ACTIVE; + } + break; + case OPEN: + if (hyd->LinkFlows[k] < -hyd->Qtol) { + status = CLOSED; + } + else if (h2 >= hset + htol) { + status = ACTIVE; + } + else { + status = OPEN; + } + break; + case CLOSED: + if (h1 >= hset + htol && h2 < hset - htol) { + status = ACTIVE; + } + else if (h1 < hset - htol && h1 > h2 + htol) { + status = OPEN; + } + else { + status = CLOSED; + } + break; + case XPRESSURE: + if (hyd->LinkFlows[k] < -hyd->Qtol) { + status = CLOSED; + } + break; + default: + break; + } + return(status); +} + + +StatType psvstatus(EN_Project *pr, int k, StatType s, double hset, double h1, double h2) +/* +**----------------------------------------------------------- +** Input: k = link index +** s = current status +** hset = valve head setting +** h1 = head at upstream node +** h2 = head at downstream node +** Output: returns new valve status +** Purpose: updates status of a pressure sustaining valve. +**----------------------------------------------------------- +*/ +{ + StatType status; /* New valve status */ + double hml; /* Minor headloss */ + hydraulics_t *hyd = &pr->hydraulics; + + double htol = hyd->Htol; + + Slink *link = &pr->network.Link[k]; + + status = s; + if (hyd->LinkSetting[k] == MISSING) { + return(status); /* Status fixed by user */ + } + hml = link->Km*SQR(hyd->LinkFlows[k]); /* Head loss when open */ + +/*** Status rules below have changed. ***/ + + switch (s) + { + case ACTIVE: + if (hyd->LinkFlows[k] < -hyd->Qtol) { + status = CLOSED; + } + else if (h2 + hml > hset + htol) { + status = OPEN; + } + else { + status = ACTIVE; + } + break; + case OPEN: + if (hyd->LinkFlows[k] < -hyd->Qtol) { + status = CLOSED; + } + else if (h1 < hset - htol) { + status = ACTIVE; + } + else { + status = OPEN; + } + break; + case CLOSED: + if (h2 > hset + htol && h1 > h2 + htol) { + status = OPEN; + } + else if (h1 >= hset + htol && h1 > h2 + htol) { + status = ACTIVE; + } + else { + status = CLOSED; + } + break; + case XPRESSURE: + if (hyd->LinkFlows[k] < -hyd->Qtol) { + status = CLOSED; + } + break; + default: + break; + } + return(status); +} + + +StatType fcvstatus(EN_Project *pr, int k, StatType s, double h1, double h2) +/* +**----------------------------------------------------------- +** Input: k = link index +** s = current status +** h1 = head at upstream node +** h2 = head at downstream node +** Output: returns new valve status +** Purpose: updates status of a flow control valve. +** +** Valve status changes to XFCV if flow reversal. +** If current status is XFCV and current flow is +** above setting, then valve becomes active. +** If current status is XFCV, and current flow +** positive but still below valve setting, then +** status remains same. +**----------------------------------------------------------- +*/ +{ + StatType status; /* New valve status */ + hydraulics_t *hyd = &pr->hydraulics; + + status = s; + if (h1 - h2 < -hyd->Htol) { + status = XFCV; + } + else if (hyd->LinkFlows[k] < -hyd->Qtol) { + status = XFCV; + } + else if (s == XFCV && hyd->LinkFlows[k] >= hyd->LinkSetting[k]) { + status = ACTIVE; + } + return(status); +} + + +/*** Updated 9/7/00 ***/ +/*** Updated 11/19/01 ***/ +void tankstatus(EN_Project *pr, int k, int n1, int n2) +/* +**---------------------------------------------------------------- +** Input: k = link index +** n1 = start node of link +** n2 = end node of link +** Output: none +** Purpose: closes link flowing into full or out of empty tank +**---------------------------------------------------------------- +*/ +{ + int i, n; + double h, q; + Stank *tank; + + hydraulics_t *hyd = &pr->hydraulics; + EN_Network *net = &pr->network; + Slink *link = &net->Link[k]; + + /* Make node n1 be the tank */ + q = hyd->LinkFlows[k]; + i = n1 - net->Njuncs; + if (i <= 0) { + i = n2 - net->Njuncs; + if (i <= 0) { + return; + } + n = n1; + n1 = n2; + n2 = n; + q = -q; + } + h = hyd->NodeHead[n1] - hyd->NodeHead[n2]; + tank = &net->Tank[i]; + /* Skip reservoirs & closed links */ + if (tank->A == 0.0 || hyd->LinkStatus[k] <= CLOSED) { + return; + } + + /* If tank full, then prevent flow into it */ + if (hyd->NodeHead[n1] >= tank->Hmax - hyd->Htol) + { + + /* Case 1: Link is a pump discharging into tank */ + if (link->Type == EN_PUMP) { + if (link->N2 == n1) + hyd->LinkStatus[k] = TEMPCLOSED; + } + + /* Case 2: Downstream head > tank head */ + /* (i.e., an open outflow check valve would close) */ + else if (cvstatus(pr, OPEN, h, q) == CLOSED) { + hyd->LinkStatus[k] = TEMPCLOSED; + } + } + + /* If tank empty, then prevent flow out of it */ + if (hyd->NodeHead[n1] <= tank->Hmin + hyd->Htol) { + + /* Case 1: Link is a pump discharging from tank */ + if (link->Type == EN_PUMP) { + if (link->N1 == n1) { + hyd->LinkStatus[k] = TEMPCLOSED; + } + } + + /* Case 2: Tank head > downstream head */ + /* (i.e., a closed outflow check valve would open) */ + else if (cvstatus(pr, CLOSED, h, q) == OPEN) { + hyd->LinkStatus[k] = TEMPCLOSED; + } + } +} /* End of tankstatus */ + + +int pswitch(EN_Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: returns 1 if status of any link changes, 0 if not +** Purpose: adjusts settings of links controlled by junction +** pressures after a hydraulic solution is found +**-------------------------------------------------------------- +*/ +{ + int i, /* Control statement index */ + k, /* Link being controlled */ + n, /* Node controlling Slink */ + reset, /* Flag on control conditions */ + change, /* Flag for status or setting change */ + anychange = 0; /* Flag for 1 or more changes */ + char s; /* Current link status */ + + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + report_options_t *rep = &pr->report; + Slink *link; + + /* Check each control statement */ + for (i = 1; i <= net->Ncontrols; i++) + { + reset = 0; + if ((k = net->Control[i].Link) <= 0) { + continue; + } + + /* Determine if control based on a junction, not a tank */ + if ((n = net->Control[i].Node) > 0 && n <= net->Njuncs) { + /* Determine if control conditions are satisfied */ + if (net->Control[i].Type == LOWLEVEL + && hyd->NodeHead[n] <= net->Control[i].Grade + hyd->Htol) { + reset = 1; + } + if (net->Control[i].Type == HILEVEL + && hyd->NodeHead[n] >= net->Control[i].Grade - hyd->Htol) { + reset = 1; + } + } + + /* Determine if control forces a status or setting change */ + if (reset == 1) + { + link = &net->Link[k]; + change = 0; + s = hyd->LinkStatus[k]; + if (link->Type == EN_PIPE) { + if (s != net->Control[i].Status) { + change = 1; + } + } + if (link->Type == EN_PUMP) { + if (hyd->LinkSetting[k] != net->Control[i].Setting) { + change = 1; + } + } + if (link->Type >= EN_PRV) { + if (hyd->LinkSetting[k] != net->Control[i].Setting) { + change = 1; + } + else if (hyd->LinkSetting[k] == MISSING && s != net->Control[i].Status) { + change = 1; + } + } + + /* If a change occurs, update status & setting */ + if (change) { + hyd->LinkStatus[k] = net->Control[i].Status; + if (link->Type > EN_PIPE) { + hyd->LinkSetting[k] = net->Control[i].Setting; + } + if (rep->Statflag == FULL) { + writestatchange(pr, k, s, hyd->LinkStatus[k]); + } + + /* Re-set flow if status has changed */ + // if (S[k] != s) initlinkflow(k, S[k], K[k]); + anychange = 1; + } + } + } + return(anychange); +} /* End of pswitch */ + + +double newflows(EN_Project *pr, Hydbalance *hbal) +/* +**---------------------------------------------------------------- +** Input: none +** Output: returns solution convergence error +** Purpose: updates link flows after new nodal heads computed +**---------------------------------------------------------------- +*/ +{ + double dh, /* Link head loss */ + dq; /* Link flow change */ + double dqsum, /* Network flow change */ + qsum; /* Network total flow */ + int k, n, n1, n2; + + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + Slink *link; + + /* Initialize net inflows (i.e., demands) at tanks */ + for (n = net->Njuncs + 1; n <= net->Nnodes; n++) { + hyd->NodeDemand[n] = 0.0; + } + + /* Initialize sum of flows & corrections */ + qsum = 0.0; + dqsum = 0.0; + + hbal->maxflowchange = 0.0; + hbal->maxflowlink = -1; + + /* Update flows in all links */ + for (k = 1; k <= net->Nlinks; k++) + { + link = &net->Link[k]; + /* + ** Apply flow update formula: + ** dq = Y - P*(new head loss) + ** P = 1/(dh/dq) + ** Y = P*(head loss based on current flow) + ** where P & Y were computed in newcoeffs(). + */ + + n1 = link->N1; + n2 = link->N2; + dh = hyd->NodeHead[n1] - hyd->NodeHead[n2]; + dq = sol->Y[k] - sol->P[k] * dh; + + /* Adjust flow change by the relaxation factor */ + dq *= hyd->RelaxFactor; + + /* Prevent flow in constant HP pumps from going negative */ + if (link->Type == EN_PUMP) { + n = findpump(net, k); + if (net->Pump[n].Ptype == CONST_HP && dq > hyd->LinkFlows[k]) { + dq = hyd->LinkFlows[k] / 2.0; + } + } + hyd->LinkFlows[k] -= dq; + + /* Update sum of absolute flows & flow corrections */ + qsum += ABS(hyd->LinkFlows[k]); + dqsum += ABS(dq); + + /* Update identity of link with max. flow change */ + if (ABS(dq) > hbal->maxflowchange) { + hbal->maxflowchange = ABS(dq); + hbal->maxflowlink = k; + } + + /* Update net flows to tanks */ + if (hyd->LinkStatus[k] > CLOSED) + { + if (n1 > net->Njuncs) { + hyd->NodeDemand[n1] -= hyd->LinkFlows[k]; + } + if (n2 > net->Njuncs) { + hyd->NodeDemand[n2] += hyd->LinkFlows[k]; + } + } + } + + /* Update emitter flows */ + for (k = 1; k <= net->Njuncs; k++) + { + if (net->Node[k].Ke == 0.0) { + continue; + } + dq = emitflowchange(pr, k); + hyd->EmitterFlows[k] -= dq; + qsum += ABS(hyd->EmitterFlows[k]); + dqsum += ABS(dq); + } + + /* Return ratio of total flow corrections to total flow */ + if (qsum > hyd->Hacc) { + return(dqsum / qsum); + } + else { + return(dqsum); + } +} /* End of newflows */ + + +double emitflowchange(EN_Project *pr, int i) +/* +**-------------------------------------------------------------- +** Input: i = node index +** Output: returns change in flow at an emitter node +** Purpose: computes flow change at an emitter node +**-------------------------------------------------------------- +*/ +{ + double ke, p; + hydraulics_t *hyd = &pr->hydraulics; + EN_Network *n = &pr->network; + Snode *node = &n->Node[i]; + + ke = MAX(CSMALL, node->Ke); + p = hyd->Qexp * ke * pow(ABS(hyd->EmitterFlows[i]), (hyd->Qexp - 1.0)); + if (p < hyd->RQtol) { + p = 1 / hyd->RQtol; + } + else { + p = 1.0 / p; + } + return(hyd->EmitterFlows[i] / hyd->Qexp - p * (hyd->NodeHead[i] - node->El)); +} + + +void checkhydbalance(EN_Project *pr, Hydbalance *hbal) +/* +**-------------------------------------------------------------- +** Input: hbal = hydraulic balance errors +** Output: none +** Purpose: finds the link with the largest head imbalance +**-------------------------------------------------------------- +*/ +{ + int k, n1, n2; + double dh, headerror, headloss; + + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + solver_t *sol = &hyd->solver; + Slink *link; + hbal->maxheaderror = 0.0; + hbal->maxheadlink = -1; + for (k = 1; k <= net->Nlinks; k++) { + if (sol->P[k] == 0.0) continue; + link = &net->Link[k]; + n1 = link->N1; + n2 = link->N2; + dh = hyd->NodeHead[n1] - hyd->NodeHead[n2]; + headloss = sol->Y[k] / sol->P[k]; + headerror = ABS(dh - headloss); + if (headerror > hbal->maxheaderror) { + hbal->maxheaderror = headerror; + hbal->maxheadlink = k; + } + } +} + + +int hasconverged(EN_Project *pr, double *relerr, Hydbalance *hbal) +/* +**-------------------------------------------------------------- +** Input: relerr = current total relative flow change +** hbal = current hydraulic balance errors +** Output: returns 1 if system has converged or 0 if not +** Purpose: checks various criteria to see if system has +** become hydraulically balanced +**-------------------------------------------------------------- +*/ +{ + hydraulics_t *hyd = &pr->hydraulics; + + if (*relerr > hyd->Hacc) return 0; + if (hyd->HeadErrorLimit > 0.0 && + hbal->maxheaderror > hyd->HeadErrorLimit) return 0; + if (hyd->FlowChangeLimit > 0.0 && + hbal->maxflowchange > hyd->FlowChangeLimit) return 0; + return 1; +} + + +void reporthydbal(EN_Project *pr, Hydbalance *hbal) +/* +**-------------------------------------------------------------- +** Input: hbal = current hydraulic balance errors +** Output: none +** Purpose: identifies links with largest flow change and +** largest head loss error. +**-------------------------------------------------------------- +*/ +{ + double qchange = hbal->maxflowchange * pr->Ucf[FLOW]; + double herror = hbal->maxheaderror * pr->Ucf[HEAD]; + int qlink = hbal->maxflowlink; + int hlink = hbal->maxheadlink; + if (qlink >= 1) { + sprintf(pr->Msg, FMT66, qchange, pr->network.Link[qlink].ID); + writeline(pr, pr->Msg); + } + if (hlink >= 1) { + sprintf(pr->Msg, FMT67, herror, pr->network.Link[hlink].ID); + writeline(pr, pr->Msg); + } +} diff --git a/src/inpfile.c b/src/inpfile.c index 1a2b17d..b6a5fe3 100644 --- a/src/inpfile.c +++ b/src/inpfile.c @@ -618,6 +618,8 @@ int saveinpfile(EN_Project *pr, char *fname) fprintf(f, "\n CHECKFREQ %-d", hyd->CheckFreq); fprintf(f, "\n MAXCHECK %-d", hyd->MaxCheck); fprintf(f, "\n DAMPLIMIT %-.8f", hyd->DampLimit); + fprintf(f, "\n HEADERROR %-.8f", hyd->HeadErrorLimit * pr->Ucf[HEAD]); + fprintf(f, "\n FLOWCHANGE %-.8f", hyd->FlowChangeLimit * pr->Ucf[FLOW]); /* Write [REPORT] section */ diff --git a/src/input1.c b/src/input1.c index 8bb0658..b59096a 100644 --- a/src/input1.c +++ b/src/input1.c @@ -40,9 +40,9 @@ AUTHOR: L. Rossman --------------------- Module Global Variables ---------------------- */ -#define MAXITER 200 /* Default max. # hydraulic iterations */ -#define HACC 0.001 /* Default hydraulics convergence ratio */ -#define HTOL 0.0005 /* Default hydraulic head tolerance (ft) */ +#define MAXITER 200 /* Default max. # hydraulic iterations */ +#define HACC 0.001 /* Default hydraulics convergence ratio */ +#define HTOL 0.0005 /* Default hydraulic head tolerance (ft) */ /*** Updated 11/19/01 ***/ #define QTOL 0.0001 /* Default flow rate tolerance (cfs) */ @@ -117,19 +117,23 @@ void setdefaults(EN_Project *pr) strncpy(qu->ChemUnits, u_MGperL, MAXID); strncpy(par->DefPatID, DEFPATID, MAXID); out->Hydflag = SCRATCH; /* No external hydraulics file */ - qu->Qualflag = NONE; /* No quality simulation */ + qu->Qualflag = NONE; /* No quality simulation */ hyd->Formflag = HW; /* Use Hazen-Williams formula */ par->Unitsflag = US; /* US unit system */ par->Flowflag = GPM; /* Flow units are gpm */ par->Pressflag = PSI; /* Pressure units are psi */ rep->Tstatflag = SERIES; /* Generate time series output */ pr->Warnflag = FALSE; /* Warning flag is off */ - hyd->Htol = HTOL; /* Default head tolerance */ - hyd->Qtol = QTOL; /* Default flow tolerance */ - hyd->Hacc = HACC; /* Default hydraulic accuracy */ + hyd->Htol = HTOL; /* Default head tolerance */ + hyd->Qtol = QTOL; /* Default flow tolerance */ + hyd->Hacc = HACC; /* Default hydraulic accuracy */ + + hyd->FlowChangeLimit = 0.0; /* Default flow change limit */ + hyd->HeadErrorLimit = 0.0; /* Default head error limit */ + qu->Ctol = MISSING; /* No pre-set quality tolerance */ - hyd->MaxIter = MAXITER; /* Default max. hydraulic trials */ - hyd->ExtraIter = -1; /* Stop if network unbalanced */ + hyd->MaxIter = MAXITER; /* Default max. hydraulic trials */ + hyd->ExtraIter = -1; /* Stop if network unbalanced */ time->Dur = 0; /* 0 sec duration (steady state) */ time->Tstart = 0; /* Starting time of day */ time->Pstart = 0; /* Starting pattern period */ @@ -529,8 +533,8 @@ void initunits(EN_Project *pr) pr->Ucf[POWER] = wcf; pr->Ucf[VOLUME] = hcf * hcf * hcf; if (time->Hstep < 1800) /* Report time in mins. */ - { /* if hydraulic time step */ - pr->Ucf[TIME] = 1.0 / 60.0; /* is less than 1/2 hour. */ + { /* if hydraulic time step */ + pr->Ucf[TIME] = 1.0 / 60.0; /* is less than 1/2 hour. */ strcpy(rep->Field[TIME].Units, u_MINUTES); } else { pr->Ucf[TIME] = 1.0 / 3600.0; @@ -604,6 +608,10 @@ void convertunits(EN_Project *pr) tank->V1max *= tank->Vmax; } + /* Convert hydraulic convergence criteria */ + hyd->FlowChangeLimit /= pr->Ucf[FLOW]; + hyd->HeadErrorLimit /= pr->Ucf[HEAD]; + /* Convert WQ option concentration units */ qu->Climit /= pr->Ucf[QUALITY]; qu->Ctol /= pr->Ucf[QUALITY]; @@ -676,8 +684,9 @@ void convertunits(EN_Project *pr) } +////// Moved to inithyd() in hydraul.c /////// /* Compute flow resistances */ - resistance(pr, k); + //resistance(pr, k); } /* Convert units on control settings */ diff --git a/src/input3.c b/src/input3.c index a011af1..9c66f63 100644 --- a/src/input3.c +++ b/src/input3.c @@ -1913,6 +1913,10 @@ int optionvalue(EN_Project *pr, int n) ** SPECIFIC GRAVITY value ** TRIALS value ** ACCURACY value + +** HEADLIMIT value +** FLOWLIMIT value + ** TOLERANCE value ** SEGMENTS value (not used) ** ------ Undocumented Options ----- @@ -1928,18 +1932,20 @@ int optionvalue(EN_Project *pr, int n) hydraulics_t *hyd = &pr->hydraulics; quality_t *qu = &pr->quality; parser_data_t *par = &pr->parser; + char* tok0 = par->Tok[0]; int nvalue = 1; /* Index of token with numerical value */ double y; /* Check for obsolete SEGMENTS keyword */ - if (match(par->Tok[0], w_SEGMENTS)) - return (0); + //if (match(par->Tok[0], w_SEGMENTS)) + if (match(tok0, w_SEGMENTS)) + return (0); /* Check for missing value (which is permissible) */ - if (match(par->Tok[0], w_SPECGRAV) || match(par->Tok[0], w_EMITTER) || - match(par->Tok[0], w_DEMAND)) + if (match(tok0, w_SPECGRAV) || match(tok0, w_EMITTER) || + match(tok0, w_DEMAND)) nvalue = 2; if (n < nvalue) return (0); @@ -1949,7 +1955,7 @@ int optionvalue(EN_Project *pr, int n) return (213); /* Check for WQ tolerance option (which can be 0) */ - if (match(par->Tok[0], w_TOLERANCE)) { + if (match(tok0, w_TOLERANCE)) { if (y < 0.0) return (213); qu->Ctol = y; /* Quality tolerance*/ @@ -1957,7 +1963,7 @@ int optionvalue(EN_Project *pr, int n) } /* Check for Diffusivity option */ - if (match(par->Tok[0], w_DIFFUSIVITY)) { + if (match(tok0, w_DIFFUSIVITY)) { if (y < 0.0) return (213); qu->Diffus = y; @@ -1965,42 +1971,59 @@ int optionvalue(EN_Project *pr, int n) } /* Check for Damping Limit option */ - if (match(par->Tok[0], w_DAMPLIMIT)) { + if (match(tok0, w_DAMPLIMIT)) { hyd->DampLimit = y; return (0); } + + /* Check for flow change limit*/ + else if (match(tok0, w_FLOWCHANGE)) + { + if (y < 0.0) return 213; + hyd->FlowChangeLimit = y; + return 0; + } + + /* Check for head error limit*/ + else if (match(tok0, w_HEADERROR)) + { + if (y < 0.0) return 213; + hyd->HeadErrorLimit = y; + return 0; + } /* All other options must be > 0 */ if (y <= 0.0) return (213); /* Assign value to specified option */ - if (match(par->Tok[0], w_VISCOSITY)) + if (match(tok0, w_VISCOSITY)) hyd->Viscos = y; /* Viscosity */ - else if (match(par->Tok[0], w_SPECGRAV)) + else if (match(tok0, w_SPECGRAV)) hyd->SpGrav = y; /* Spec. gravity */ - else if (match(par->Tok[0], w_TRIALS)) + else if (match(tok0, w_TRIALS)) hyd->MaxIter = (int)y; /* Max. trials */ - else if (match(par->Tok[0], w_ACCURACY)) /* Accuracy */ + else if (match(tok0, w_ACCURACY)) /* Accuracy */ { y = MAX(y, 1.e-5); y = MIN(y, 1.e-1); hyd->Hacc = y; - } else if (match(par->Tok[0], w_HTOL)) + } + else if (match(tok0, w_HTOL)) hyd->Htol = y; - else if (match(par->Tok[0], w_QTOL)) + else if (match(tok0, w_QTOL)) hyd->Qtol = y; - else if (match(par->Tok[0], w_RQTOL)) { + else if (match(tok0, w_RQTOL)) { if (y >= 1.0) return (213); hyd->RQtol = y; - } else if (match(par->Tok[0], w_CHECKFREQ)) + } else if (match(tok0, w_CHECKFREQ)) hyd->CheckFreq = (int)y; - else if (match(par->Tok[0], w_MAXCHECK)) + else if (match(tok0, w_MAXCHECK)) hyd->MaxCheck = (int)y; - else if (match(par->Tok[0], w_EMITTER)) + else if (match(tok0, w_EMITTER)) hyd->Qexp = 1.0 / y; - else if (match(par->Tok[0], w_DEMAND)) + else if (match(tok0, w_DEMAND)) hyd->Dmult = y; else return (201); diff --git a/src/main.c b/src/main.c new file mode 100644 index 0000000..bda132f --- /dev/null +++ b/src/main.c @@ -0,0 +1,135 @@ +#include +#include +#include "epanet2.h" + +#define MAXMSG 255 /* Max. # characters in message text */ +#define MAXWARNCODE 99 +/* text copied here, no more need of include "text.h" */ +#define FMT01 "\nEPANET Version %d.%d.%d" +#define FMT03 "\n Correct syntax is:\n %s \n" +#define FMT09 "\n\nEPANET completed." +#define FMT10 "\nEPANET completed. There are warnings." +#define FMT11 "\nEPANET completed. There are errors." + + +void writeConsole(char *s); + + +/* +---------------------------------------------------------------- + Entry point used to compile a stand-alone executable. +---------------------------------------------------------------- +*/ + + +int main(int argc, char *argv[]) +/*-------------------------------------------------------------- + ** Input: argc = number of command line arguments + ** *argv[] = array of command line arguments + ** Output: none + ** Purpose: main program segment + ** + ** Command line for stand-alone operation is: + ** progname f1 f2 f3 + ** where progname = name of executable this code was compiled to, + ** f1 = name of input file, + ** f2 = name of report file (optional, stdout if left blank) + ** f3 = name of binary output file (optional, nullfile if left blank). + **-------------------------------------------------------------- + */ +{ + char *f1,*f2,*f3; + char blank[] = ""; + char errmsg[MAXMSG+1]=""; + int errcode; + int version; + char s[256]; + int major; + int minor; + int patch; + + /* get version from DLL and trasform in Major.Minor.Patch format + instead of hardcoded version */ + ENgetversion(&version); + major= version/10000; + minor= (version%10000)/100; + patch= version%100; + sprintf(s,FMT01, major , minor, patch); + writeConsole(s); + + + /* Check for proper number of command line arguments */ + if (argc < 2) { + sprintf(s, FMT03, argv[0]); + writeConsole(s); + return(1); + } + + /* set inputfile name */ + f1 = argv[1]; + if (argc > 2) { + /* set rptfile name */ + f2 = argv[2]; + } + else { + /* use stdout for rptfile */ + f2 = blank; + } + if (argc > 3) { + /* set binary output file name */ + f3 = argv[3]; + } + else { + /* NO binary output*/ + f3 = blank; + } + + /* Call the main control function */ + if (strlen(f2)> 0) { + /* use stdout for progress messages */ + //errcode = ENepanet(f1,f2,f3,writeConsole); + errcode = ENepanet(f1, f2, f3, NULL); + } + else { + /* use stdout for reporting, no progress messages */ + errcode = ENepanet(f1,f2,f3,NULL); + } + + /* Error/Warning check */ + if (errcode == 0) { + /* no errors */ + writeConsole(FMT09); + return(0); + } + else { + if (errcode > MAXWARNCODE) printf("\n Fatal Error: "); + ENgeterror(errcode, errmsg, MAXMSG); + writeConsole(errmsg); + if (errcode > MAXWARNCODE) { + // error // + writeConsole(FMT11); + return(errcode); + } + else { + // warning // + writeConsole(FMT10); + return(0); + } + } + +} /* End of main */ + + +void writeConsole(char *s) +/*---------------------------------------------------------------- + ** Input: text string + ** Output: none + ** Purpose: writes string of characters to console + **---------------------------------------------------------------- + */ +{ + fprintf(stdout,"%s\n",s); + fflush(stdout); +} + + diff --git a/src/output.c b/src/output.c index 91c54ea..39687cf 100644 --- a/src/output.c +++ b/src/output.c @@ -252,60 +252,75 @@ int saveenergy(EN_Project *pr) */ { - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - out_file_t *out = &pr->out_files; - parser_data_t *par = &pr->parser; + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + out_file_t *out = &pr->out_files; + parser_data_t *par = &pr->parser; time_options_t *time = &pr->time_options; - FILE *outFile = out->OutFile; + FILE *outFile = out->OutFile; + Spump *pump; int i, j; INT4 index; - REAL4 x[6]; /* work array */ - double hdur, /* total time->Duration in hours */ - t; /* pumping time->Duration */ + REAL4 x[MAX_ENERGY_STATS]; // work array + double hdur, // total simulation duration in hours + t; // total pumping time duration hdur = time->Dur / 3600.0; for (i = 1; i <= net->Npumps; i++) { - Spump *pump = &net->Pump[i]; + pump = &net->Pump[i]; if (hdur == 0.0) { - for (j = 0; j < 5; j++) - x[j] = (REAL4)pump->Energy[j]; - x[5] = (REAL4)(pump->Energy[5] * 24.0); - } else { - t = pump->Energy[0]; - x[0] = (REAL4)(t / hdur); - x[1] = 0.0f; - x[2] = 0.0f; - x[3] = 0.0f; - x[4] = 0.0f; - if (t > 0.0) { - x[1] = (REAL4)(pump->Energy[1] / t); - x[2] = (REAL4)(pump->Energy[2] / t); - x[3] = (REAL4)(pump->Energy[3] / t); - } - x[4] = (REAL4)pump->Energy[4]; - x[5] = (REAL4)(pump->Energy[5] * 24.0 / hdur); + pump->Energy[TOTAL_COST] *= 24.0; } - x[0] *= 100.0f; - x[1] *= 100.0f; - /* Compute Kw-hr per MilGal (or per cubic meter) */ - if (par->Unitsflag == SI) - x[2] *= (REAL4)(1000.0 / LPSperCFS / 3600.0); - else - x[2] *= (REAL4)(1.0e6 / GPMperCFS / 60.0); - for (j = 0; j < 6; j++) - pump->Energy[j] = x[j]; + else { + // ... convert total hrs. online to fraction of total time online + t = pump->Energy[PCNT_ONLINE]; //currently holds total hrs. online + pump->Energy[PCNT_ONLINE] = t / hdur; + + // ... convert cumulative values to time-averaged ones + if (t > 0.0) { + pump->Energy[PCNT_EFFIC] /= t; + pump->Energy[KWH_PER_FLOW] /= t; + pump->Energy[TOTAL_KWH] /= t; + } + + // ... convert total cost to cost per day + pump->Energy[TOTAL_COST] *= 24.0 / hdur; + } + + // ... express time online and avg. efficiency as percentages + pump->Energy[PCNT_ONLINE] *= 100.0; + pump->Energy[PCNT_EFFIC] *= 100.0; + + // ... compute KWH per Million Gallons or per Cubic Meter + if (par->Unitsflag == SI) { + pump->Energy[KWH_PER_FLOW] *= (1000. / LPSperCFS / 3600.); + } + else { + pump->Energy[KWH_PER_FLOW] *= (1.0e6 / GPMperCFS / 60.); + } + + // ... save energy stats to REAL4 work array + for (j = 0; j < MAX_ENERGY_STATS; j++) { + x[j] = (REAL4)pump->Energy[j]; + } + + // ... save energy results to output file index = pump->Link; - if (fwrite(&index, sizeof(INT4), 1, outFile) < 1) + if (fwrite(&index, sizeof(INT4), 1, outFile) < 1) { return (308); - if (fwrite(x, sizeof(REAL4), 6, outFile) < 6) + } + if (fwrite(x, sizeof(REAL4), MAX_ENERGY_STATS, outFile) < MAX_ENERGY_STATS) { return (308); + } } + + // ... compute and save demand charge hyd->Emax = hyd->Emax * hyd->Dcost; x[0] = (REAL4)hyd->Emax; - if (fwrite(&x[0], sizeof(REAL4), 1, outFile) < 1) + if (fwrite(&x[0], sizeof(REAL4), 1, outFile) < 1) { return (308); + } return (0); } diff --git a/src/report.c b/src/report.c index 9620f1f..8d2101e 100644 --- a/src/report.c +++ b/src/report.c @@ -229,6 +229,15 @@ void writesummary(EN_Project *pr) sprintf(s, FMT27, hyd->Hacc); writeline(pr, s); + if (hyd->HeadErrorLimit > 0.0) { + sprintf(s, FMT27d, hyd->HeadErrorLimit*pr->Ucf[HEAD], rep->Field[HEAD].Units); + writeline(pr, s); + } + if (hyd->FlowChangeLimit > 0.0) { + sprintf(s, FMT27e, hyd->FlowChangeLimit*pr->Ucf[FLOW], rep->Field[FLOW].Units); + writeline(pr, s); + } + sprintf(s, FMT27a, hyd->CheckFreq); writeline(pr, s); sprintf(s, FMT27b, hyd->MaxCheck); @@ -378,6 +387,7 @@ void writeenergy(EN_Project *pr) EN_Network *net = &pr->network; hydraulics_t *hyd = &pr->hydraulics; report_options_t *rep = &pr->report; + Spump *pump; int j; double csum; @@ -388,14 +398,14 @@ void writeenergy(EN_Project *pr) writeheader(pr,ENERHDR, 0); csum = 0.0; for (j = 1; j <= net->Npumps; j++) { - Spump *pump = &net->Pump[j]; - csum += pump->Energy[5]; + pump = &net->Pump[j]; + csum += pump->Energy[TOTAL_COST]; if (rep->LineNum == (long)rep->PageSize) writeheader(pr, ENERHDR, 1); sprintf(s, "%-8s %6.2f %6.2f %9.2f %9.2f %9.2f %9.2f", - net->Link[pump->Link].ID, pump->Energy[0], pump->Energy[1], - pump->Energy[2], pump->Energy[3], pump->Energy[4], - pump->Energy[5]); + net->Link[pump->Link].ID, pump->Energy[PCNT_ONLINE], pump->Energy[PCNT_EFFIC], + pump->Energy[KWH_PER_FLOW], pump->Energy[TOTAL_KWH], pump->Energy[MAX_KW], + pump->Energy[TOTAL_COST]); writeline(pr, s); } fillstr(s, '-', 63); diff --git a/src/text.h b/src/text.h index be47e39..85159a9 100755 --- a/src/text.h +++ b/src/text.h @@ -140,7 +140,10 @@ AUTHOR: L. Rossman #define w_RQTOL "RQTOL" #define w_CHECKFREQ "CHECKFREQ" #define w_MAXCHECK "MAXCHECK" -#define w_DAMPLIMIT "DAMPLIMIT" +#define w_DAMPLIMIT "DAMPLIMIT" + +#define w_FLOWCHANGE "FLOWCHANGE" +#define w_HEADERROR "HEADERROR" #define w_SECONDS "SEC" #define w_MINUTES "MIN" @@ -364,7 +367,10 @@ AUTHOR: L. Rossman #define FMT27a " Status Check Frequency ............ %-d" #define FMT27b " Maximum Trials Checked ............ %-d" -#define FMT27c " Damping Limit Threshold ........... %-.6f" +#define FMT27c " Damping Limit Threshold ........... %-.6f" + +#define FMT27d " Headloss Error Limit .............. %-.6f %s" +#define FMT27e " Flow Change Limit ................. %-.6f %s" #define FMT28 " Maximum Trials .................... %-d" #define FMT29 " Quality Analysis .................. None" @@ -409,10 +415,13 @@ AUTHOR: L. Rossman #define FMT61 "%10s: Valve %s caused ill-conditioning" #define FMT62 "%10s: System ill-conditioned at node %s" #define FMT63 "%10s: %s %s changed by rule %s" -#define FMT64 "%10s: Balancing the network:" +#define FMT64 "%10s: Balancing the network:\n" #define FMT65 " Trial %2d: relative flow change = %-.6f" /*** End of update ***/ +#define FMT66 " maximum flow change = %.4f for Link %s" +#define FMT67 " maximum head error = %.4f for Link %s\n" + /* -------------------- Energy Report Table ------------------- */ #define FMT71 "Energy Usage:" #define FMT72 \ diff --git a/src/types.h b/src/types.h index 927eac9..83441aa 100755 --- a/src/types.h +++ b/src/types.h @@ -47,9 +47,9 @@ typedef int INT4; #define EOFMARK 0x1A /* Use 0x04 for UNIX systems */ #define MAXTITLE 3 /* Max. # title lines */ #define MAXID 31 /* Max. # characters in ID name */ -#define MAXMSG 79 /* Max. # characters in message text */ +#define MAXMSG 79 /* Max. # characters in message text */ #define MAXLINE 255 /* Max. # characters read from input line */ -#define MAXFNAME 260 /* Max. # characters in file name */ +#define MAXFNAME 259 /* Max. # characters in file name */ #define MAXTOKS 40 /* Max. items per line of input */ #define TZERO 1.E-4 /* Zero time tolerance */ #define TRUE 1 @@ -59,6 +59,9 @@ typedef int INT4; #define TINY 1.E-6 #define MISSING -1.E10 +#define CBIG 1.e8 /* Big coefficient */ +#define CSMALL 1.e-6 /* Small coefficient */ + #ifdef M_PI #define PI M_PI #else @@ -135,21 +138,21 @@ typedef enum { USE, /* use from previous run */ SAVE, /* save after current run */ SCRATCH -} Hydtype; /* use temporary file */ +} Hydtype; /* use temporary file */ typedef enum { NONE, /* no quality analysis */ CHEM, /* analyze a chemical */ AGE, /* analyze water age */ TRACE -} QualType; /* trace % of flow from a source */ +} QualType; /* trace % of flow from a source */ typedef enum { - V_CURVE, /* volume curve */ - P_CURVE, /* pump curve */ - E_CURVE, /* efficiency curve */ + V_CURVE, /* volume curve */ + P_CURVE, /* pump curve */ + E_CURVE, /* efficiency curve */ H_CURVE -} CurveType; /* head loss curve */ +} CurveType; /* head loss curve */ typedef enum { CONST_HP, /* constant horsepower */ @@ -169,8 +172,8 @@ typedef enum { LOWLEVEL, /* act when grade below set level */ HILEVEL, /* act when grade above set level */ TIMER, /* act when set time reached */ - TIMEOFDAY -} ControlType; /* act when time of day occurs */ + TIMEOFDAY /* act when time of day occurs */ +} ControlType; typedef enum { XHEAD, /* pump cannot deliver head (closed) */ @@ -183,18 +186,18 @@ typedef enum { XPRESSURE, /* valve cannot supply pressure */ FILLING, /* tank filling */ EMPTYING -} StatType; /* tank emptying */ +} StatType; /* tank emptying */ typedef enum { HW, /* Hazen-Williams */ DW, /* Darcy-Weisbach */ - CM -} FormType; /* Chezy-Manning */ + CM /* Chezy-Manning */ +} FormType; typedef enum { US, /* US */ - SI -} UnitsType; /* SI (metric) */ + SI /* SI (metric) */ +} UnitsType; typedef enum { CFS, /* cubic feet per second */ @@ -206,40 +209,41 @@ typedef enum { LPM, /* liters per minute */ MLD, /* megaliters per day */ CMH, /* cubic meters per hour */ - CMD -} FlowUnitsType; /* cubic meters per day */ + CMD /* cubic meters per day */ +} FlowUnitsType; typedef enum { PSI, /* pounds per square inch */ KPA, /* kiloPascals */ - METERS -} PressUnitsType; /* meters */ + METERS /* meters */ +} PressUnitsType; typedef enum { LOW, /* lower limit */ HI, /* upper limit */ - PREC -} RangeType; /* precision */ + PREC /* precision */ +} RangeType; typedef enum { MIX1, /* 1-compartment model */ MIX2, /* 2-compartment model */ FIFO, /* First in, first out model */ - LIFO -} MixType; /* Last in, first out model */ + LIFO /* Last in, first out model */ +} MixType; typedef enum { SERIES, /* none */ AVG, /* time-averages */ MIN, /* minimum values */ MAX, /* maximum values */ - RANGE -} TstatType; /* max - min values */ + RANGE /* max - min values */ +} TstatType; -#define MAXVAR 21 /* Max. # types of network variables */ -/* (equals # items enumed below) */ + +#define MAXVAR 21 /* Max. # types of network variables */ + /* (equals # items enumed below) */ typedef enum { - ELEV = 0, /* nodal elevation */ + ELEV = 0, /* nodal elevation */ DEMAND, /* nodal demand flow */ HEAD, /* nodal hydraulic head */ PRESSURE, /* nodal pressure */ @@ -273,17 +277,27 @@ typedef enum { } SectType; typedef enum { - STATHDR, /* Hydraulic Status */ - ENERHDR, /* Energy Usage */ - NODEHDR, /* Node Results */ - LINKHDR -} HdrType; /* Link Results */ + STATHDR, /* Hydraulic Status header */ + ENERHDR, /* Energy Usage header */ + NODEHDR, /* Node Results header */ + LINKHDR /* Link Results header */ +} HdrType; typedef enum { POSITIVE, NEGATIVE } FlowDirection; +typedef enum { + PCNT_ONLINE, + PCNT_EFFIC, + KWH_PER_FLOW, + TOTAL_KWH, + MAX_KW, + TOTAL_COST, + MAX_ENERGY_STATS +} EnergyStats; + /* ------------------------------------------------------ Global Data Structures @@ -344,6 +358,16 @@ struct Sdemand /* DEMAND CATEGORY OBJECT */ }; typedef struct Sdemand *Pdemand; /* Pointer to demand object */ +typedef struct +{ + double hrsOnLine; // hours pump is online + double efficiency; // total time wtd. efficiency + double kwHrsPerCFS; // total kw-hrs per cfs of flow + double kwHrs; // total kw-hrs consumed + double maxKwatts; // max. kw consumed + double totalCost; // total pumping cost +} Senergy; + struct Ssource /* WQ SOURCE OBJECT */ { /*int Node;*/ /* Node index of source */ @@ -422,13 +446,7 @@ typedef struct /* PUMP OBJECT */ int Upat; /* Utilization pattern index */ int Epat; /* Energy cost pattern index */ double Ecost; /* Unit energy cost */ - double Energy[6]; /* Energy usage statistics: */ - /* 0 = pump utilization */ - /* 1 = avg. efficiency */ - /* 2 = avg. kW/flow */ - /* 3 = avg. kwatts */ - /* 4 = peak kwatts */ - /* 5 = cost/day */ + double Energy[MAX_ENERGY_STATS]; /* Energy usage statistics */ } Spump; typedef struct /* VALVE OBJECT */ @@ -474,7 +492,7 @@ typedef struct /* FIELD OBJECT of report table */ double RptLim[2]; /* Lower/upper report limits */ } SField; -typedef struct s_Premise /* Rule Premise Clause */ +typedef struct s_Premise /* Rule Premise Clause */ { int logop; /* Logical operator */ int object; /* Node or link */ @@ -486,7 +504,7 @@ typedef struct s_Premise /* Rule Premise Clause */ struct s_Premise *next; } Premise; -typedef struct s_Action /* Rule Action Clause */ +typedef struct s_Action /* Rule Action Clause */ { int link; /* Link index */ int status; /* Link's status */ @@ -494,19 +512,19 @@ typedef struct s_Action /* Rule Action Clause */ struct s_Action *next; } Action; -typedef struct s_aRule /* Control Rule Structure */ +typedef struct s_aRule /* Control Rule Structure */ { - char label[MAXID+1]; /* Rule character label */ - double priority; /* Priority level */ - Premise *Pchain; /* Linked list of premises */ - Action *Tchain; /* Linked list of actions if true */ - Action *Fchain; /* Linked list of actions if false */ + char label[MAXID+1]; /* Rule character label */ + double priority; /* Priority level */ + Premise *Pchain; /* Linked list of premises */ + Action *Tchain; /* Linked list of actions if true */ + Action *Fchain; /* Linked list of actions if false */ struct s_aRule *next; } aRule; -typedef struct s_ActItem /* Action list item */ +typedef struct s_ActItem /* Action list item */ { - int ruleindex; /* Index of rule action belongs to */ + int ruleindex; /* Index of rule action belongs to */ struct s_Action *action; /* An action structure */ struct s_ActItem *next; } ActItem; @@ -614,15 +632,15 @@ typedef struct { *Patlist, /* Temporary time pattern list */ *Curvelist; /* Temporary list of curves */ - double *X; // temporary array for curve data + double *X; // temporary array for curve data int - Ntokens, /* Number of tokens in input line */ - Ntitle; /* Number of title lines */ + Ntokens, /* Number of tokens in input line */ + Ntitle; /* Number of title lines */ - char *Tok[MAXTOKS]; /* Array of token strings */ + char *Tok[MAXTOKS]; /* Array of token strings */ char Comment[MAXMSG+1]; - STmplist *PrevPat; /* Pointer to pattern list element */ - STmplist *PrevCurve; /* Pointer to curve list element */ + STmplist *PrevPat; /* Pointer to pattern list element */ + STmplist *PrevCurve; /* Pointer to curve list element */ } parser_data_t; @@ -647,7 +665,7 @@ typedef struct { Rpt1Fname[MAXFNAME+1], /* Primary report file name */ Rpt2Fname[MAXFNAME+1]; /* Secondary report file name */ - SField Field[MAXVAR]; /* Output reporting fields */ + SField Field[MAXVAR]; /* Output reporting fields */ long LineNum; /* Current line number */ long PageNum; /* Current page number */ @@ -694,12 +712,11 @@ typedef struct { aRule *Rule; /* Array of rules */ ActItem *ActList; /* Linked list of action items */ - int RuleState; /* State of rule interpreter */ - long Time1; /* Start of rule evaluation time interval (sec) */ + int RuleState; /* State of rule interpreter */ + long Time1; /* Start of rule evaluation time interval (sec) */ Premise *Plast; /* Previous premise clause */ Action *Tlast; /* Previous true action */ Action *Flast; /* Previous false action */ - } rules_t; /* @@ -754,6 +771,10 @@ typedef struct { Qexp, /* Exponent in orifice formula */ Dmult, /* Demand multiplier */ Hacc, /* Hydraulics solution accuracy */ + + FlowChangeLimit, /* Hydraulics flow change limit */ + HeadErrorLimit, /* Hydraulics head error limit */ + DampLimit, /* Solution damping threshold */ Viscos, /* Kin. viscosity (sq ft/sec) */ SpGrav, /* Specific gravity */ @@ -797,7 +818,7 @@ typedef struct { } hydraulics_t; typedef struct { - int Nnodes, /* Number of network nodes */ + int Nnodes, /* Number of network nodes */ Ntanks, /* Number of tanks */ Njuncs, /* Number of junction nodes */ Nlinks, /* Number of network links */ @@ -810,17 +831,18 @@ typedef struct { Ncurves, /* Number of data curves */ Ncoords; /* Number of Coords */ - Snode *Node; /* Node data */ - Slink *Link; /* Link data */ - Stank *Tank; /* Tank data */ - Spump *Pump; /* Pump data */ - Svalve *Valve; /* Valve data */ - Spattern *Pattern; /* Time patterns */ - Scurve *Curve; /* Curve data */ - Scoord *Coord; /* Coordinate data */ - Scontrol *Control; /* Control data */ - ENHashTable *NodeHashTable, *LinkHashTable; /* Hash tables for ID labels */ - Padjlist *Adjlist; /* Node adjacency lists */ + Snode *Node; /* Node data */ + Slink *Link; /* Link data */ + Stank *Tank; /* Tank data */ + Spump *Pump; /* Pump data */ + Svalve *Valve; /* Valve data */ + Spattern *Pattern; /* Time patterns */ + Scurve *Curve; /* Curve data */ + Scoord *Coord; /* Coordinate data */ + Scontrol *Control; /* Control data */ + ENHashTable *NodeHashTable, + *LinkHashTable; /* Hash tables for ID labels */ + Padjlist *Adjlist; /* Node adjacency lists */ } EN_Network; @@ -852,7 +874,4 @@ struct EN_Project { }; - - - #endif From 9afb6f0b28912dbfe08248a7e7c2b6aa11e94cb0 Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Sat, 16 Jun 2018 11:10:07 -0400 Subject: [PATCH 2/9] Merge branch 'contributor-lr' of https://github.com/LRossman/EPANET into contributor-lr --- .travis.yml | 2 +- CMakeLists.txt | 53 +++++++++++++++++++++++++------------------------- README.md | 8 +++++--- appveyor.yml | 3 ++- 4 files changed, 35 insertions(+), 31 deletions(-) diff --git a/.travis.yml b/.travis.yml index f3d61f2..3cf3c7f 100644 --- a/.travis.yml +++ b/.travis.yml @@ -16,7 +16,7 @@ before_install: before_script: - mkdir -p $BUILD_HOME - cd $BUILD_HOME - - cmake .. + - cmake -DBUILD_TESTS=1 .. script: - cmake --build . diff --git a/CMakeLists.txt b/CMakeLists.txt index 808b836..59e1356 100755 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,27 +1,26 @@ -# CMakeLists.txt - CMake configuration file for EPANET 2.0 -# -# CMake is a cross-platform build tool. CMake generates platform native -# makefiles that can be used with your compiler of choice. CMake uses a -# generator concept to represent different build tooling. CMake automatically +# CMakeLists.txt - CMake configuration file for EPANET 2.0 +# +# CMake is a cross-platform build tool. CMake generates platform native +# makefiles that can be used with your compiler of choice. CMake uses a +# generator concept to represent different build tooling. CMake automatically # detects the platform it is running on and generates the appropriate makefiles -# for the platform default compiler. Different generators can also be specified. -# -# Note: CMake requires that your platform build system and compiler are -# properly installed. Build using Visual Studio requires msbuild shell. -# -# Example Usage: -# cd build/cmake +# for the platform default compiler. Different generators can also be specified. +# +# Note: CMake requires that your platform build system and compiler are +# properly installed. Build using Visual Studio requires msbuild shell. +# +# Example Usage: # mkdir buildproducts # cd buildproducts # cmake .. -# make -# +# make +# # Building MSYS on Windows: # ... # cmake -G "MSYS Makefiles" .. -# make -# -# Building Visual Studio on Windows: +# make +# +# Building Visual Studio on Windows: # ... # cmake -G "Visual Studio 10 2010" .. # msbuild /p:Configuration=Release ALL_BUILD.vcxproj @@ -29,23 +28,25 @@ # Generic Invocation: # cmake -E make_directory buildprod # cd build -# cmake -G GENERATOR -DCMAKE_BUILD_TYPE=Release .. +# cmake -G GENERATOR -DCMAKE_BUILD_TYPE=Release -DBUILD_TESTS=1 .. # cmake --build . --target SOME_TARGET --config Release # # More information: # cmake --help -# +# # CMake is available at https://cmake.org/download/ -# +# cmake_minimum_required (VERSION 2.8.8) project(EPANET) add_subdirectory(tools/epanet-output) -add_subdirectory(tests) +IF (BUILD_TESTS) + add_subdirectory(tests) +ENDIF (BUILD_TESTS) -# Sets for output directory for executables and libraries. +# Sets for output directory for executables and libraries. set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin) set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib) set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/lib) @@ -55,8 +56,8 @@ SET(CMAKE_POSITION_INDEPENDENT_CODE ON) IF (APPLE) - SET(CMAKE_INSTALL_NAME_DIR @executable_path) - SET(CMAKE_BUILD_WITH_INSTALL_RPATH ON) +SET(INSTALL_NAME_DIR @executable_path/../lib) +SET(CMAKE_MACOSX_RPATH 1) ENDIF (APPLE) IF (MSVC) @@ -89,8 +90,8 @@ target_include_directories(epanet PUBLIC ${PROJECT_SOURCE_DIR}/include) # EXPORT_MACRO_NAME DLLEXPORT # EXPORT_FILE_NAME epanet_export.h # STATIC_DEFINE SHARED_EXPORTS_BUILT_AS_STATIC) -# -#file(COPY ${CMAKE_CURRENT_BINARY_DIR}/epanet_export.h +# +#file(COPY ${CMAKE_CURRENT_BINARY_DIR}/epanet_export.h # DESTINATION ${CMAKE_CURRENT_SOURCE_DIR}/include) diff --git a/README.md b/README.md index 4a697d0..dd6f5d7 100755 --- a/README.md +++ b/README.md @@ -5,9 +5,11 @@ EPANET {#epanet-readme} [![Build status](https://ci.appveyor.com/api/projects/status/19wpg4g2cmj3oihl?svg=true)](https://ci.appveyor.com/project/OpenWaterAnalytics/epanet) [![Build Status](https://travis-ci.org/OpenWaterAnalytics/EPANET.svg?branch=master)](https://travis-ci.org/OpenWaterAnalytics/EPANET) -The EPANET Library is a pressurized pipe network hydraulic and water quality analysis toolkit written in C. +## For EPANET-related questions and discussion +For community discussion, FAQ, and roadmapping of the project, go to the [Community Forum](http://community.wateranalytics.org/category/epanet). -If you are interested in using/extending EPANET for academic, personal, or commercial use, then you've come to the right place. For community discussion, FAQ, and roadmapping of the project, go to the [Community Forum](http://community.wateranalytics.org/category/epanet). +## What is on this Repository? +The EPANET Library is a pressurized pipe network hydraulic and water quality analysis toolkit written in C. If you are interested in using/extending EPANET for academic, personal, or commercial use, then you've come to the right place. Please see the [`version 2.1` Release Notes](https://github.com/OpenWaterAnalytics/EPANET/blob/master/ReleaseNotes2_1.md) for information relevant to users of the previous official version (2.00.12). If you would like to contribute by addressing any of the outstanding [Issues](https://github.com/OpenWaterAnalytics/EPANET/issues), then please comment on the Issue, then Fork this repo to your own account and base your commits on the [`dev` branch](https://github.com/OpenWaterAnalytics/EPANET/tree/dev). Once you are finished, you can open a Pull Request to test the code and discuss merging your changes back into the community respository. @@ -15,4 +17,4 @@ A step-by-step tutorial on how to contribute to EPANET using GitHub is also [ava __Note:__ This repository is not affiliated with, or endorsed by, the USEPA. For the last "official" release of EPANET (2.00.12 UI and Toolkit) please go to the [EPA's GitHub repo](https://github.com/USEPA/Water-Distribution-Network-Model) or [the USEPA website](http://www2.epa.gov/water-research/epanet). It is also not the graphical user interface version. This is the hydraulic and water quality solver engine. -However, if you are interested in extending EPANET for academic, personal, or commercial use, then you've come to the right place. For community discussion, FAQ, and roadmapping of the project, go to the [Community Forum](http://community.wateranalytics.org/category/epanet). \ No newline at end of file +However, if you are interested in extending EPANET for academic, personal, or commercial use, then you've come to the right place. For community discussion, FAQ, and roadmapping of the project, go to the [Community Forum](http://community.wateranalytics.org/category/epanet). diff --git a/appveyor.yml b/appveyor.yml index c7e97ea..d9ef393 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -23,7 +23,7 @@ environment: BOOST_ROOT: "C:/Libraries/boost" # New build on Visual Studio 15 2017 - APPVEYOR_BUILD_WORKER_IMAGE: Visual Studio 2017 - GENERATOR: "Visual Studio 15 2017 Win64" + GENERATOR: "Visual Studio 15 2017" GROUP: "EXPERIMENTAL" BOOST_ROOT: "C:/Libraries/boost_1_67_0" @@ -48,6 +48,7 @@ before_build: - mkdir %BUILD_HOME% - cd %BUILD_HOME% - cmake -G "%GENERATOR%" + -DBUILD_TESTS=1 -DBOOST_ROOT="%BOOST_ROOT%" -DBoost_USE_STATIC_LIBS="ON" .. From 42f70f579f92594b6b739ba0bf70572f115c0418 Mon Sep 17 00:00:00 2001 From: Elad Salomons Date: Sun, 17 Jun 2018 07:48:25 +0300 Subject: [PATCH 3/9] Remove strange string from types.h and fix win_build make file --- src/types.h | 2 +- win_build/WinSDK/Makefile.bat | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/types.h b/src/types.h index 89a5467..6babe9a 100755 --- a/src/types.h +++ b/src/types.h @@ -50,7 +50,7 @@ typedef int INT4; #define MAXMSG 79 /* Max. # characters in message text */ #define MAXLINE 255 /* Max. # characters read from input line */ #define MAXFNAME 259 /* Max. # characters in file name */ ->>>>>>> f11308fc72eccc8b74b18d5ddab0eb2ae0d36587 +//>>>>>>> f11308fc72eccc8b74b18d5ddab0eb2ae0d36587 #define MAXTOKS 40 /* Max. items per line of input */ #define TZERO 1.E-4 /* Zero time tolerance */ #define TRUE 1 diff --git a/win_build/WinSDK/Makefile.bat b/win_build/WinSDK/Makefile.bat index eeb967e..d0a00c3 100644 --- a/win_build/WinSDK/Makefile.bat +++ b/win_build/WinSDK/Makefile.bat @@ -21,9 +21,9 @@ Find /i "x86" < checkOS.tmp > StringCheck.tmp If %ERRORLEVEL% == 1 ( CALL "%SDK_PATH%bin\"SetEnv.cmd /x64 /release rem : create EPANET2.DLL - cl -o epanet2.dll epanet.c hash.c hydraul.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /link /DLL + cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /link /DLL rem : create EPANET2.EXE - cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /I ..\src /link + cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /I ..\src /link md "%Build_PATH%"\64bit move /y "%SRC_PATH%"\*.dll "%Build_PATH%"\64bit move /y "%SRC_PATH%"\*.exe "%Build_PATH%"\64bit @@ -35,9 +35,9 @@ rem : 32 bit with DEF CALL "%SDK_PATH%bin\"SetEnv.cmd /x86 /release echo "32 bit with epanet2.def mapping" rem : create EPANET2.DLL -cl -o epanet2.dll epanet.c hash.c hydraul.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /link /DLL /def:..\win_build\WinSDK\epanet2.def /MAP +cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /link /DLL /def:..\win_build\WinSDK\epanet2.def /MAP rem : create EPANET2.EXE -cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /I ..\src /link +cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /I ..\src /link md "%Build_PATH%"\32bit move /y "%SRC_PATH%"\*.dll "%Build_PATH%"\32bit move /y "%SRC_PATH%"\*.exe "%Build_PATH%"\32bit From c41e140b6aebc5bc2d41463c67a3db0c69756917 Mon Sep 17 00:00:00 2001 From: Elad Salomons Date: Tue, 19 Jun 2018 11:16:32 +0300 Subject: [PATCH 4/9] Update VB header files --- include/epanet2.bas | 2 ++ include/epanet2.vb | 2 ++ 2 files changed, 4 insertions(+) diff --git a/include/epanet2.bas b/include/epanet2.bas index 789a441..0826a3b 100644 --- a/include/epanet2.bas +++ b/include/epanet2.bas @@ -128,6 +128,8 @@ Public Const EN_ACCURACY = 1 Public Const EN_TOLERANCE = 2 Public Const EN_EMITEXPON = 3 Public Const EN_DEMANDMULT = 4 +Public Const EN_HEADERROR = 5 +Public Const EN_FLOWCHANGE = 6 Public Const EN_LOWLEVEL = 0 ' Control types Public Const EN_HILEVEL = 1 diff --git a/include/epanet2.vb b/include/epanet2.vb index 44a47a1..fcb71af 100644 --- a/include/epanet2.vb +++ b/include/epanet2.vb @@ -121,6 +121,8 @@ Public Const EN_ACCURACY = 1 Public Const EN_TOLERANCE = 2 Public Const EN_EMITEXPON = 3 Public Const EN_DEMANDMULT = 4 +Public Const EN_HEADERROR = 5 +Public Const EN_FLOWCHANGE = 6 Public Const EN_LOWLEVEL = 0 ' Control types Public Const EN_HILEVEL = 1 From e9303de078bb1a3cd910ecdec6bd86a570580776 Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Tue, 19 Jun 2018 10:30:02 -0400 Subject: [PATCH 5/9] Fixes to implement issue #161 --- src/epanet2.h | 1236 ----------------------------------------------- src/hydcoeffs.c | 10 +- src/hydraul.c | 13 +- src/hydsolver.c | 24 +- src/main.c | 135 ------ src/types.h | 1 - 6 files changed, 25 insertions(+), 1394 deletions(-) delete mode 100644 src/epanet2.h delete mode 100644 src/main.c diff --git a/src/epanet2.h b/src/epanet2.h deleted file mode 100644 index 81b469e..0000000 --- a/src/epanet2.h +++ /dev/null @@ -1,1236 +0,0 @@ -/** @file epanet2.h - @see http://github.com/openwateranalytics/epanet - - */ - -/* - ******************************************************************* - - EPANET2.H - Prototypes for EPANET Functions Exported to DLL Toolkit - - VERSION: 2.00 - DATE: 5/8/00 - 10/25/00 - 3/1/01 - 8/15/07 (2.00.11) - 2/14/08 (2.00.12) - AUTHORS: L. Rossman - US EPA - NRMRL - OpenWaterAnalytics members: see git stats for contributors - - ******************************************************************* - */ - - -#ifndef EPANET2_H -#define EPANET2_H - -// the toolkit can be compiled with support for double-precision as well. -// just make sure that you use the correct #define in your client code. -#ifndef EN_API_FLOAT_TYPE - #define EN_API_FLOAT_TYPE float -#endif - -// --- define WINDOWS -#undef WINDOWS -#ifdef _WIN32 - #define WINDOWS -#endif -#ifdef __WIN32__ - #define WINDOWS -#endif - -// --- define DLLEXPORT -#ifndef DLLEXPORT - #ifdef WINDOWS - #ifdef __cplusplus - #define DLLEXPORT extern "C" __declspec(dllexport) - #else - #define DLLEXPORT __declspec(dllexport) __stdcall - #endif // __cplusplus - #elif defined(CYGWIN) - #define DLLEXPORT __stdcall - #elif defined(__APPLE__) - #ifdef __cplusplus - #define DLLEXPORT - #else - #define DLLEXPORT - #endif - #else - #define DLLEXPORT - #endif -#endif - -//#include "epanet_export.h" - -// --- Define the EPANET toolkit constants - -/// Node property codes -typedef enum { - EN_ELEVATION = 0, /**< Node Elevation */ - EN_BASEDEMAND = 1, /**< Node Base Demand, from last demand category */ - EN_PATTERN = 2, /**< Node Demand Pattern */ - EN_EMITTER = 3, /**< Node Emitter Coefficient */ - EN_INITQUAL = 4, /**< Node initial quality */ - EN_SOURCEQUAL = 5, /**< Node source quality */ - EN_SOURCEPAT = 6, /**< Node source pattern index */ - EN_SOURCETYPE = 7, /**< Node source type */ - EN_TANKLEVEL = 8, /**< Tank Level */ - EN_DEMAND = 9, /**< Node current simulated demand */ - EN_HEAD = 10, /**< Node Head value */ - EN_PRESSURE = 11, /**< Node pressure value */ - EN_QUALITY = 12, /**< Node quality value */ - EN_SOURCEMASS = 13, /**< Node source mass value */ - EN_INITVOLUME = 14, /**< Tank or Reservoir initial volume */ - EN_MIXMODEL = 15, /**< Tank mixing model */ - EN_MIXZONEVOL = 16, /**< Tank mixing zone volume */ - EN_TANKDIAM = 17, /**< Tank diameter */ - EN_MINVOLUME = 18, /**< Tank minimum volume */ - EN_VOLCURVE = 19, /**< Tank volume curve index */ - EN_MINLEVEL = 20, /**< Tank minimum level */ - EN_MAXLEVEL = 21, /**< Tank maximum level */ - EN_MIXFRACTION = 22, /**< Tank mixing fraction */ - EN_TANK_KBULK = 23, /**< Tank bulk decay coefficient */ - EN_TANKVOLUME = 24, /**< Tank current volume */ - EN_MAXVOLUME = 25 /**< Tank maximum volume */ -} EN_NodeProperty; - -/// Link property codes -typedef enum { - EN_DIAMETER = 0, - EN_LENGTH = 1, - EN_ROUGHNESS = 2, - EN_MINORLOSS = 3, - EN_INITSTATUS = 4, - EN_INITSETTING = 5, - EN_KBULK = 6, - EN_KWALL = 7, - EN_FLOW = 8, - EN_VELOCITY = 9, - EN_HEADLOSS = 10, - EN_STATUS = 11, - EN_SETTING = 12, - EN_ENERGY = 13, - EN_LINKQUAL = 14, - EN_LINKPATTERN = 15, - EN_EFFICIENCY = 16, - EN_HEADCURVE = 17, - EN_EFFICIENCYCURVE = 18, - EN_PRICEPATTERN = 19 -} EN_LinkProperty; - -/// Time parameter codes -typedef enum { - EN_DURATION = 0, - EN_HYDSTEP = 1, - EN_QUALSTEP = 2, - EN_PATTERNSTEP = 3, - EN_PATTERNSTART = 4, - EN_REPORTSTEP = 5, - EN_REPORTSTART = 6, - EN_RULESTEP = 7, - EN_STATISTIC = 8, - EN_PERIODS = 9, - EN_STARTTIME = 10, - EN_HTIME = 11, - EN_QTIME = 12, - EN_HALTFLAG = 13, - EN_NEXTEVENT = 14, - EN_NEXTEVENTIDX = 15 -} EN_TimeProperty; - - -typedef enum { - EN_ITERATIONS = 0, - EN_RELATIVEERROR = 1 -} EN_AnalysisStatistic; - -typedef enum { - EN_NODECOUNT = 0, /**< Number of Nodes (Juntions + Tanks + Reservoirs) */ - EN_TANKCOUNT = 1, /**< Number of Tanks and Reservoirs */ - EN_LINKCOUNT = 2, /**< Number of Links (Pipes + Pumps + Valves) */ - EN_PATCOUNT = 3, /**< Number of Time Patterns */ - EN_CURVECOUNT = 4, /**< Number of Curves */ - EN_CONTROLCOUNT = 5, /**< Number of Control Statements */ - EN_RULECOUNT = 6 /**< Number of Rule-based Control Statements */ -} EN_CountType; - -typedef enum { - EN_JUNCTION = 0, - EN_RESERVOIR = 1, - EN_TANK = 2 -} EN_NodeType; - -typedef enum { - EN_CVPIPE = 0, /* Link types. */ - EN_PIPE = 1, /* See LinkType in TYPES.H */ - EN_PUMP = 2, - EN_PRV = 3, - EN_PSV = 4, - EN_PBV = 5, - EN_FCV = 6, - EN_TCV = 7, - EN_GPV = 8 -} EN_LinkType; - -typedef enum { - EN_NONE = 0, /* Quality analysis types. */ - EN_CHEM = 1, /* See QualType in TYPES.H */ - EN_AGE = 2, - EN_TRACE = 3 -} EN_QualityType; - -typedef enum { - EN_CONCEN = 0, /* Source quality types. */ - EN_MASS = 1, /* See SourceType in TYPES.H. */ - EN_SETPOINT = 2, - EN_FLOWPACED = 3 -} EN_SourceType; - -typedef enum { /* Head loss formula: */ - EN_HW = 0, /* Hazen-Williams */ - EN_DW = 1, /* Darcy-Weisbach */ - EN_CM = 2 /* Chezy-Manning */ -} EN_FormType; /* See FormType in TYPES.H */ - -typedef enum { - EN_CFS = 0, /* Flow units types. */ - EN_GPM = 1, /* See FlowUnitsType */ - EN_MGD = 2, /* in TYPES.H. */ - EN_IMGD = 3, - EN_AFD = 4, - EN_LPS = 5, - EN_LPM = 6, - EN_MLD = 7, - EN_CMH = 8, - EN_CMD = 9 -} EN_FlowUnits; - - -/// Simulation Option codes -typedef enum { - EN_TRIALS = 0, - EN_ACCURACY = 1, - EN_TOLERANCE = 2, - EN_EMITEXPON = 3, - EN_DEMANDMULT = 4, - EN_HEADERROR = 5, - EN_FLOWCHANGE = 6 -} EN_Option; - -typedef enum { - EN_LOWLEVEL = 0, /* Control types. */ - EN_HILEVEL = 1, /* See ControlType */ - EN_TIMER = 2, /* in TYPES.H. */ - EN_TIMEOFDAY = 3 -} EN_ControlType; - - - -typedef enum { - EN_AVERAGE = 1, /* Time statistic types. */ - EN_MINIMUM = 2, /* See TstatType in TYPES.H */ - EN_MAXIMUM = 3, - EN_RANGE = 4 -} EN_StatisticType; - - - -typedef enum { - EN_MIX1 = 0, /* Tank mixing models */ - EN_MIX2 = 1, - EN_FIFO = 2, - EN_LIFO = 3 -} EN_MixingModel; - - - -typedef enum { - EN_NOSAVE = 0, - EN_SAVE = 1, - EN_INITFLOW = 10, - EN_SAVE_AND_INIT = 11 -} EN_SaveOption; - - - -typedef enum { - EN_CONST_HP = 0, /* constant horsepower */ - EN_POWER_FUNC = 1, /* power function */ - EN_CUSTOM = 2 /* user-defined custom curve */ -} EN_CurveType; - - - -// --- Declare the EPANET toolkit functions -#if defined(__cplusplus) -extern "C" { -#endif - - /** - @brief The EPANET Project wrapper object - */ - typedef struct EN_Project EN_Project; - typedef struct EN_Pattern EN_Pattern; - typedef struct EN_Curve EN_Curve; - - /** - @brief runs a complete EPANET simulation - @param inpFile pointer to name of input file (must exist) - @param rptFile pointer to name of report file (to be created) - @param binOutFile pointer to name of binary output file (to be created) - @param callback a callback function that takes a character string (char *) as its only parameter. - @return error code - - The callback function should reside in and be used by the calling - code to display the progress messages that EPANET generates - as it carries out its computations. If this feature is not - needed then the argument should be NULL. - */ - int DLLEXPORT ENepanet(char *inpFile, char *rptFile, char *binOutFile, void (*callback) (char *)); - - /** - @brief Initializes an EPANET session - @param rptFile pointer to name of report file (to be created) - @param binOutFile pointer to name of binary output file (to be created) - @param UnitsType flow units flag - @param HeadlossFormula headloss formula flag - @return error code - */ - int DLLEXPORT ENinit(char *rptFile, char *binOutFile, int UnitsType, int HeadlossFormula); - - /** - @brief Opens EPANET input file & reads in network data - @param inpFile pointer to name of input file (must exist) - @param rptFile pointer to name of report file (to be created) - @param binOutFile pointer to name of binary output file (to be created) - @return error code - */ - int DLLEXPORT ENopen(char *inpFile, char *rptFile, char *binOutFile); - - /** - @brief Saves current data to "INP" formatted text file. - @param filename The file path to create - @return Error code - */ - int DLLEXPORT ENsaveinpfile(char *filename); - - /** - @brief Frees all memory and files used by EPANET - @return Error code - */ - int DLLEXPORT ENclose(); - - /** - @brief Solves the network hydraulics for all time periods - @return Error code - */ - int DLLEXPORT ENsolveH(); - - /** - @brief Saves hydraulic results to binary file - @return Error code - - Must be called before ENreport() if no WQ simulation has been made. - Should not be called if ENsolveQ() will be used. - */ - int DLLEXPORT ENsaveH(); - - /** - @brief Sets up data structures for hydraulic analysis - @return Error code - */ - int DLLEXPORT ENopenH(); - - /** - @brief Initializes hydraulic analysis - @param initFlag 2-digit flag where 1st (left) digit indicates if link flows should be re-initialized (1) or not (0), and 2nd digit indicates if hydraulic results should be saved to file (1) or not (0). - @return Error code - */ - int DLLEXPORT ENinitH(int initFlag); - - /** - @brief Run a hydraulic solution period - @param[out] currentTime The current simulation time in seconds - @return Error or warning code - @see ENsolveH - - This function is used in a loop with ENnextH() to run - an extended period hydraulic simulation. - See ENsolveH() for an example. - */ - int DLLEXPORT ENrunH(long *currentTime); - - /** - @brief Determine time (in seconds) until next hydraulic event - @param[out] tStep Time (seconds) until next hydraulic event. 0 marks end of simulation period. - @return Error code - - This function is used in a loop with ENrunH() to run an extended period hydraulic simulation. - See ENsolveH() for an example. - */ - int DLLEXPORT ENnextH(long *tStep); - - - /** - @brief Frees data allocated by hydraulics solver - @return Error code - */ - int DLLEXPORT ENcloseH(); - - /** - @brief Copies binary hydraulics file to disk - @param filename Name of file to be created - @return Error code - */ - int DLLEXPORT ENsavehydfile(char *filename); - - /** - @brief Opens previously saved binary hydraulics file - @param filename Name of file to be used - @return Error code - */ - int DLLEXPORT ENusehydfile(char *filename); - - /** - @brief Solves for network water quality in all time periods - @return Error code - */ - int DLLEXPORT ENsolveQ(); - - /** - @brief Sets up data structures for WQ analysis - @return Error code - */ - int DLLEXPORT ENopenQ(); - - /** - @brief Initializes water quality analysis - @param saveFlag EN_SAVE (1) if results saved to file, EN_NOSAVE (0) if not - @return Error code - */ - int DLLEXPORT ENinitQ(int saveFlag); - - /** - @brief Retrieves hydraulic & WQ results at time t. - @param[out] currentTime Current simulation time, in seconds. - @return Error code - - This function is used in a loop with ENnextQ() to run - an extended period WQ simulation. See ENsolveQ() for - an example. - */ - int DLLEXPORT ENrunQ(long *currentTime); - - /** - @brief Advances WQ simulation to next hydraulic event. - @param[out] tStep Time in seconds until next hydraulic event. 0 marks end of simulation period. - @return Error code - - This function is used in a loop with ENrunQ() to run - an extended period WQ simulation. See ENsolveQ() for - an example. - */ - int DLLEXPORT ENnextQ(long *tStep); - - /** - @brief Advances WQ simulation by a single WQ time step - @param[out] timeLeft Time left in overall simulation (in seconds) - @return Error code - - This function is used in a loop with ENrunQ() to run - an extended period WQ simulation. - */ - int DLLEXPORT ENstepQ(long *timeLeft); - - /** - @brief Frees data allocated by water quality solver. - @return Error code. - */ - int DLLEXPORT ENcloseQ(); - - /** - @brief Writes line of text to the report file. - @param line Text string to write - @return Error code. - */ - int DLLEXPORT ENwriteline(char *line); - - /** - @brief Writes simulation report to the report file - @return Error code - */ - int DLLEXPORT ENreport(); - - /** - @brief Resets report options to default values - @return Error code - */ - int DLLEXPORT ENresetreport(); - - /** - @brief Processes a reporting format command - @return Error code - */ - int DLLEXPORT ENsetreport(char *reportFormat); - - /** - @brief Retrieves parameters that define a simple control - @param controlIndex Control index (position of control statement in the input file, starting from 1) - @param[out] controlType Control type code (see EPANET2.H) - @param[out] linkIndex Index of controlled link - @param[out] setting Control setting on link - @param[out] nodeIndex Index of controlling node (0 for TIMER or TIMEOFDAY control) - @param[out] level Control level (tank level, junction pressure, or time (seconds)) - @return Error code - */ - int DLLEXPORT ENgetcontrol(int controlIndex, int *controlType, int *linkIndex, EN_API_FLOAT_TYPE *setting, int *nodeIndex, EN_API_FLOAT_TYPE *level); - - /** - @brief Retrieves the number of components of a given type in the network. - @param code Component code (see EPANET2.H) - @param[out] count Number of components in network - @return Error code - */ - int DLLEXPORT ENgetcount(int code, int *count); - - /** - @brief Gets value for an analysis option - @param code Option code (see EPANET2.H) - @param[out] value Option value - @return Error code - */ - int DLLEXPORT ENgetoption(int code, EN_API_FLOAT_TYPE *value); - - /** - @brief Retrieves value of specific time parameter. - @param code Time parameter code - @param[out] value Value of time parameter. - @return Error code - */ - int DLLEXPORT ENgettimeparam(int code, long *value); - - /** - @brief Retrieves the flow units code - @param[out] code Code of flow units in use - @return Error code - */ - int DLLEXPORT ENgetflowunits(int *code); - - /** - @brief Sets the flow units - @param code Code of flow units to use - @return Error code - */ - int DLLEXPORT ENsetflowunits(int code); - - /** - @brief Retrieves the index of the time pattern with specified ID - @param id String ID of the time pattern - @param[out] index Index of the specified time pattern - @return Error code - @see ENgetpatternid - */ - int DLLEXPORT ENgetpatternindex(char *id, int *index); - - /** - @brief Retrieves ID of a time pattern with specific index. - @param index The index of a time pattern. - @param[out] id The string ID of the time pattern. - @return Error code - @see ENgetpatternindex - */ - int DLLEXPORT ENgetpatternid(int index, char *id); - - /** - @brief Retrieves the number of multipliers in a time pattern. - @param index The index of a time pattern. - @param[out] len The length of the time pattern. - @return Error code - */ - int DLLEXPORT ENgetpatternlen(int index, int *len); - - /** - @brief Retrive a multiplier from a pattern for a specific time period. - @param index The index of a time pattern - @param period The pattern time period. First time period is 1. - @param[out] value Pattern multiplier - @return Error code - */ - int DLLEXPORT ENgetpatternvalue(int index, int period, EN_API_FLOAT_TYPE *value); - - /** - @brief Retrieve the average multiplier value in a time pattern - @param index The index of a time pattern - @param[out] value The average of all of this time pattern's values - @return Error code - */ - int DLLEXPORT ENgetaveragepatternvalue(int index, EN_API_FLOAT_TYPE *value); - - /** - @brief Retrieve the type of quality analytis to be run. - @param[out] qualcode The quality analysis code number. - @param[out] tracenode The index of node being traced, if qualcode == trace - @return Error code - @see ENsetqualtype - */ - int DLLEXPORT ENgetqualtype(int *qualcode, int *tracenode); - - /** - @brief Get the text of an error code. - @param errcode The error code - @param[out] errmsg The error string represented by the code - @param maxLen The maximum number of characters to copy into the char pointer errmsg - @return Error code - */ - int DLLEXPORT ENgeterror(int errcode, char *errmsg, int maxLen); - - /** - @brief Get hydraulic simulation statistic - @param code The type of statistic to get - @param[out] value The value of the statistic - @return Error code - */ - int DLLEXPORT ENgetstatistic(int code, EN_API_FLOAT_TYPE* value); - - /** - @brief Get index of node with specified ID - @param id The string ID of the node - @param[out] index The node's index (first node is index 1) - @return Error code - @see ENgetnodeid - */ - int DLLEXPORT ENgetnodeindex(char *id, int *index); - - /** - @brief Get the string ID of the specified node. - @param index The index of the node (first node is index 1) - @param[out] id The string ID of the specified node. Up to MAXID characters will be copied, so id must be pre-allocated by the calling code to hold at least that many characters. - @return Error code - @see ENgetnodeindex - */ - int DLLEXPORT ENgetnodeid(int index, char *id); - - /** - @brief Get the type of node with specified index. - @param index The index of a node (first node is index 1) - @param[out] code The type code for the node. - @return Error code - */ - int DLLEXPORT ENgetnodetype(int index, int *code); - - /** - @brief Get a property value for specified node - @param index The index of a node (first node is index 1). - @param code The property type code - @param[out] value The value of the node's property. - @return Error code - @see EN_NodeProperty - */ - int DLLEXPORT ENgetnodevalue(int index, int code, EN_API_FLOAT_TYPE *value); - - /** - @brief Get coordinates (x,y) for a node. - @param index The index of a node (first node is index 1). - @param[out] x X-value of node's coordinate - @param[out] y Y-value of node's coordinate - @return Error code - @see ENsetcoord - */ - int DLLEXPORT ENgetcoord(int index, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y); - - /** - @brief Set coordinates (x,y) for a node. - @param index The index of a node (first node is index 1) - @param x X-value of node's coordinate - @param y Y-value of node's coordinate - @return Error code - @see ENgetcoord - */ - int DLLEXPORT ENsetcoord(int index, EN_API_FLOAT_TYPE x, EN_API_FLOAT_TYPE y); - - /** - @brief Get the number of demand categories for a node. - @param nodeIndex The index of a node (first node is index 1) - @param[out] numDemands The number of demand categories - @return Error code - */ - int DLLEXPORT ENgetnumdemands(int nodeIndex, int *numDemands); - - /** - @brief Get a node's base demand for a specified category. - @param nodeIndex The index of a node (first node is index 1) - @param demandIndex The index of the demand category (starting at 1) - @return Error code - */ - int DLLEXPORT ENgetbasedemand(int nodeIndex, int demandIndex, EN_API_FLOAT_TYPE *baseDemand); - - /** - @brief Get the index of the demand pattern assigned to a node for a category index. - @param nodeIndex The index of a node (first node is index 1). - @param demandIndex The index of a category (first category is index 1). - @param[out] pattIndex The index of the pattern for this node and category. - @return Error code - */ - int DLLEXPORT ENgetdemandpattern(int nodeIndex, int demandIndex, int *pattIndex); - - /** - @brief Get the index of a Link with specified ID. - @param id The string ID of a link. - @param[out] index The index of the named link (first link is index 1) - @return Error code - @see ENgetlinkid - */ - int DLLEXPORT ENgetlinkindex(char *id, int *index); - - /** - @brief Get the string ID of a link with specified index - @param index The index of a link (first link is index 1) - @param[out] id The ID of the link. Up to MAXID characters will be copied, so id must be pre-allocated by the calling code to hold at least that many characters. - @return Error code - @see ENgetlinkindex - */ - int DLLEXPORT ENgetlinkid(int index, char *id); - - /** - @brief Get the link type code for a specified link - @param index The index of a link (first link is index 1) - @param[out] code The type code of the link. - @return Error code - @see EN_LinkType - */ - int DLLEXPORT ENgetlinktype(int index, EN_LinkType *code); - - /** - @brief Set the link type code for a specified link - @param id The id of a link - @param type The type code of the link. - @return Error code - @see EN_LinkType - */ - int DLLEXPORT ENsetlinktype(char *id, EN_LinkType type); - - /** - @brief Get the indexes of a link's start- and end-nodes. - @param index The index of a link (first link is index 1) - @param[out] node1 The index of the link's start node (first node is index 1). - @param[out] node2 The index of the link's end node (first node is index 1). - @return Error code - @see ENgetnodeid, ENgetlinkid - */ - int DLLEXPORT ENgetlinknodes(int index, int *node1, int *node2); - - /** - @brief Get a property value for specified link. - @param index The index of a node (first node is index 1). - @param code The parameter desired. - @param[out] value The value of the link's specified property. - @return Error code - @see ENgetnodevalue, EN_LinkProperty - */ - int DLLEXPORT ENgetlinkvalue(int index, int code, EN_API_FLOAT_TYPE *value); - - /** - @brief Get a curve's properties. - @param curveIndex The index of a curve (first curve is index 1). - @param[out] id The curve's string ID. Client code must preallocate at least MAXID characters. - @param[out] nValues The number of values in the curve's (x,y) list. - @param[out] xValues The curve's x-values. Pointer must be freed by client. - @param[out] yValues The curve's y-values. Pointer must be freed by client. - @return Error code. - */ - int DLLEXPORT ENgetcurve(int curveIndex, char* id, int *nValues, EN_API_FLOAT_TYPE **xValues, EN_API_FLOAT_TYPE **yValues); - - /** - @brief Retrieves the curve index for a specified pump index. - @param pumpIndex The index of a pump - @param[out] curveIndex The index of the curve used by the pump. - @return Error code. - */ - int DLLEXPORT ENgetheadcurveindex(int pumpIndex, int *curveIndex); - - /** - @brief Sets the curve id for a specified pump index. - @param pumpIndex The index of the pump - @param curveIndex The index of the curve used by the pump - @return Error code. - */ - int DLLEXPORT ENsetheadcurveindex(int pumpIndex, int curveIndex); - - /** - @brief Get the type of pump - @param linkIndex The index of the pump element - @param[out] outType The integer-typed pump curve type signifier (output parameter) - @return Error code - @see EN_CurveType - */ - int DLLEXPORT ENgetpumptype(int linkIndex, int *outType); - - /** - @brief Get the version number. This number is to be interpreted with implied decimals, i.e., "20100" == "2(.)01(.)00" - @param[out] version The version of EPANET - @return Error code. - */ - int DLLEXPORT ENgetversion(int *version); - - /** - @brief Specify parameters to define a simple control - @param cindex The index of the control to edit. First control is index 1. - @param ctype The type code to set for this control. - @param lindex The index of a link to control. - @param setting The control setting applied to the link. - @param nindex The index of a node used to control the link, or 0 for TIMER / TIMEOFDAY control. - @param level control point (tank level, junction pressure, or time in seconds). - @return Error code. - */ - int DLLEXPORT ENsetcontrol(int cindex, int ctype, int lindex, EN_API_FLOAT_TYPE setting, int nindex, EN_API_FLOAT_TYPE level); - - /** - @brief Set a property value for a node. - @param index The index of a node. First node is index 1. - @param code The code for the proprty to set. - @param v The value to set for this node and property. - @return Error code. - @see EN_NodeProperty - */ - int DLLEXPORT ENsetnodevalue(int index, int code, EN_API_FLOAT_TYPE v); - - /** - @brief Set a proprty value for a link. - @param index The index of a link. First link is index 1. - @param code The code for the property to set. - @param v The value to set for this link and property. - @return Error code. - @see EN_LinkProperty - */ - int DLLEXPORT ENsetlinkvalue(int index, int code, EN_API_FLOAT_TYPE v); - - /** - @brief Add a new time pattern. - @param id The string ID of the pattern to add. - @return Error code. - @see ENgetpatternindex - */ - int DLLEXPORT ENaddpattern(char *id); - - /** - @brief Set multipliers for a specific pattern - @param index The index of a pattern. First pattern is index 1. - @param f An array of multipliers - @param len The length of array f. - @return Error code. - @see ENgetpatternindex - */ - int DLLEXPORT ENsetpattern(int index, EN_API_FLOAT_TYPE *f, int len); - - /** - @brief Set the multiplier for a specific pattern at a specific period. - @param index The index of a pattern. First pattern is index 1. - @param period The period of the pattern to set. - @param value The value of the multiplier to set. - @return Error code. - */ - int DLLEXPORT ENsetpatternvalue(int index, int period, EN_API_FLOAT_TYPE value); - - /** - @brief Set the value for a time parameter. - @param code The code for the parameter to set. - @param value The desired value of the parameter. - @return Error code. - @see EN_TimeProperty - */ - int DLLEXPORT ENsettimeparam(int code, long value); - - /** - @brief Set a value for an anlysis option. - @param code The code for the desired option. - @param v The desired value for the option specified. - @return Error code. - @see EN_Option - */ - int DLLEXPORT ENsetoption(int code, EN_API_FLOAT_TYPE v); - - /** - @brief Sets the level of hydraulic status reporting. - @param code Status reporting code. - @return Error code. - */ - int DLLEXPORT ENsetstatusreport(int code); - - /** - @brief Sets type of quality analysis called for - @param qualcode WQ parameter code, EN_QualityType - @param chemname Name of WQ constituent - @param chemunits Concentration units of WQ constituent - @param tracenode ID of node being traced (if applicable) - @return Error code. - @see EN_QualityType - - chemname and chemunits only apply when WQ analysis is for chemical. tracenode only applies when WQ analysis is source tracing. - */ - int DLLEXPORT ENsetqualtype(int qualcode, char *chemname, char *chemunits, char *tracenode); - - /** - @brief Get quality analysis information (type, chemical name, units, trace node ID) - @param[out] qualcode The EN_QualityType code being used. - @param[out] chemname The name of the WQ constituent. - @param[out] chemunits The cencentration units of the WQ constituent. - @param[out] tracenode The trace node ID. - @return Error code. - @see EN_QualityType - */ - int DLLEXPORT ENgetqualinfo(int *qualcode, char *chemname, char *chemunits, int *tracenode); - - /** - @brief Sets the node's base demand for a category. - @param nodeIndex The index of a node. - @param demandIdx The index of a demand category. - @param baseDemand The base demand multiplier for the selected category. - @return Error code. - @see ENgetbasedemand - */ - int DLLEXPORT ENsetbasedemand(int nodeIndex, int demandIdx, EN_API_FLOAT_TYPE baseDemand); - - /** - @brief Retrieves index of curve with specific ID. - @param id The ID of a curve. - @param[out] index The index of the named curve - @return Error code. - @see ENgetcurveid - */ - int DLLEXPORT ENgetcurveindex(char *id, int *index); - - /** - @brief Retrieves ID of a curve with specific index. - @param index The index of a curve. - @param[out] id The ID of the curve specified. - @return Error code. - @see ENsetcurveindex - - NOTE: 'id' must be able to hold MAXID characters - */ - int DLLEXPORT ENgetcurveid(int index, char *id); - - /** - @brief Retrieves number of points in a curve - @param index The index of a curve. - @param[out] len The length of the curve coordinate list - @return Error code. - @see ENgetcurvevalue - */ - int DLLEXPORT ENgetcurvelen(int index, int *len); - - /** - @brief retrieves x,y point for a specific point number and curve - @param curveIndex The index of a curve - @param pointIndex The index of a point in the curve - @param[out] x The x-value of the specified point. - @param[out] y The y-value of the specified point. - @return Error code. - @see ENgetcurvelen ENsetcurvevalue - */ - int DLLEXPORT ENgetcurvevalue(int curveIndex, int pointIndex, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y); - - /** - @brief Sets x,y point for a specific point and curve. - @param curveIndex The index of a curve. - @param pointIndex The index of a point in the curve. - @param x The x-value of the point. - @param y The y-value of the point. - @return Error code. - */ - int DLLEXPORT ENsetcurvevalue(int curveIndex, int pointIndex, EN_API_FLOAT_TYPE x, EN_API_FLOAT_TYPE y); - - /** - @brief Sets x,y values for a specified curve. - @param index The index of a curve. - @param x An array of x-values for the curve. - @param y An array of y-values for the curve. - @param len The length of the arrays x and y. - @return Error code. - */ - int DLLEXPORT ENsetcurve(int index, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y, int len); - - /** - @brief Adds a new curve appended to the end of the existing curves. - @param id The name of the curve to be added. - @return Error code. - @see ENgetcurveindex ENsetcurve - */ - int DLLEXPORT ENaddcurve(char *id); - - - /** - @brief Gets the number of premises, true actions, and false actions and the priority of an existing rule-based control. - @param index The index of a rule-based control. - @param nPremises The number of conditions in a rule-based control. - @param nTrueActions The number of actions that are executed when the conditions in the rule-based control are met. - @param nFalseActions The number of actions that are executed when the conditions in the rule-based control are not met. - @param priority The priority of a rule-based control. - @return Error code. - */ - int DLLEXPORT ENgetrule(int index, int *nPremises, int *nTrueActions, int *nFalseActions, EN_API_FLOAT_TYPE *priority); - - /** - @brief Sets the priority of the existing rule-based control. - @param index The index of a rule-based control. - @param priority The priority to be set in the rule-based control. - @return Error code. - */ - int DLLEXPORT ENsetrulepriority(int index, EN_API_FLOAT_TYPE priority); - - /** - @brief Gets the components of a premise/condition in an existing rule-based control. - @param indexRule The index of a rule-based control. - @param indexPremise The index of the premise. - @param logop The logiv operator (IF/AND/OR) in the premise - @param object The object (e.g. TANK) the premise is looking at. - @param indexObj The index of the object (e.g. the index of the tank). - @param variable The variable to be checked (e.g. level). - @param relop The relashionship operator (e.g. LARGER THAN) in the premise. - @param status The status of the object to be checked (e.g. CLOSED) - @param value The value of the variable to be checked (e.g. 5.5) - @return Error code. - */ - int DLLEXPORT ENgetpremise(int indexRule, int indexPremise, int *logop, int *object, int *indexObj, int *variable, int *relop, int *status, EN_API_FLOAT_TYPE *value); - - /** - @brief Sets the components of a premise/condition in an existing rule-based control. - @param indexRule The index of a rule-based control. - @param indexPremise The index of the premise. - @param logop The logiv operator (IF/AND/OR) in the premise - @param object The object (e.g. TANK) the premise is looking at. - @param indexObj The index of the object (e.g. the index of the tank). - @param variable The variable to be checked (e.g. level). - @param relop The relashionship operator (e.g. LARGER THAN) in the premise. - @param status The status of the object to be checked (e.g. CLOSED) - @param value The value of the variable to be checked (e.g. 5.5) - @return Error code. - */ - int DLLEXPORT ENsetpremise(int indexRule, int indexPremise, int logop, int object, int indexObj, int variable, int relop, int status, EN_API_FLOAT_TYPE value); - - /** - @brief Sets the index of an object in a premise of an existing rule-based control. - @param indexRule The index of a rule-based control. - @param indexPremise The index of the premise. - @param indexObj The index of the object (e.g. the index of the tank). - @return Error code. - */ - int DLLEXPORT ENsetpremiseindex(int indexRule, int indexPremise, int indexObj); - - /** - @brief Sets the status in a premise of an existing rule-based control. - @param indexRule The index of a rule-based control. - @param indexPremise The index of the premise. - @param status The status of the object to be checked (e.g. CLOSED) - @return Error code. - */ - int DLLEXPORT ENsetpremisestatus(int indexRule, int indexPremise, int status); - - /** - @brief Sets the value in a premise of an existing rule-based control. - @param indexRule The index of a rule-based control. - @param indexPremise The index of the premise. - @param value The value of the variable to be checked (e.g. 5.5) - @return Error code. - */ - int DLLEXPORT ENsetpremisevalue(int indexRule, int indexPremise, EN_API_FLOAT_TYPE value); - - /** - @brief Gets the components of a true-action in an existing rule-based control. - @param indexRule The index of a rule-based control. - @param indexAction The index of the action when the conditions in the rule are met. - @param indexLink The index of the link in the action (e.g. index of Pump 1) - @param status The status of the link (e.g. CLOSED) - @param setting The value of the link (e.g. pump speed 0.9) - @return Error code. - */ - int DLLEXPORT ENgettrueaction(int indexRule, int indexAction, int *indexLink, int *status, EN_API_FLOAT_TYPE *setting); - - /** - @brief Sets the components of a true-action in an existing rule-based control. - @param indexRule The index of a rule-based control. - @param indexAction The index of the action when the conditions in the rule are met. - @param indexLink The index of the link in the action (e.g. index of Pump 1) - @param status The status of the link (e.g. CLOSED) - @param setting The value of the link (e.g. pump speed 0.9) - @return Error code. - */ - int DLLEXPORT ENsettrueaction(int indexRule, int indexAction, int indexLink, int status, EN_API_FLOAT_TYPE setting); - - /** - @brief Gets the components of a false-action in an existing rule-based control. - @param indexRule The index of a rule-based control. - @param indexAction The index of the action when the conditions in the rule are not met. - @param indexLink The index of the link in the action (e.g. index of Pump 1) - @param status The status of the link (e.g. CLOSED) - @param setting The value of the link (e.g. pump speed 0.9) - @return Error code. - */ - int DLLEXPORT ENgetfalseaction(int indexRule, int indexAction, int *indexLink, int *status, EN_API_FLOAT_TYPE *setting); - - /** - @brief Sets the components of a false-action in an existing rule-based control. - @param indexRule The index of a rule-based control. - @param indexAction The index of the action when the conditions in the rule are not met. - @param indexLink The index of the link in the action (e.g. index of Pump 1) - @param status The status of the link (e.g. CLOSED) - @param setting The value of the link (e.g. pump speed 0.9) - @return Error code. - */ - int DLLEXPORT ENsetfalseaction(int indexRule, int indexAction, int indexLink, int status, EN_API_FLOAT_TYPE setting); - - /** - @brief Returns the ID of a rule. - @param indexRule The index of a rule-based control. - @param id The ID of the rule - @return Error code. - */ - int DLLEXPORT ENgetruleID(int indexRule, char* id); - - /** - @brief Adds a new node - @param id The name of the node to be added. - @param nodeType The node type code - @return Error code. - */ - int DLLEXPORT ENaddnode(char *id, EN_NodeType nodeType); - - /** - @brief Adds a new link - @param id The name of the link to be added. - @param linkType The link type code - @param fromNode The id of the from node - @param toNode The id of the to node - @return Error code. - */ - int DLLEXPORT ENaddlink(char *id, EN_LinkType linkType, char *fromNode, char *toNode); - - /** - @brief Deletes a node - @param nodeIndex The node index - @return Error code. - */ - int DLLEXPORT ENdeletenode(int nodeIndex); - - /** - @brief Deletes a link - @param linkIndex The link index - @return Error code. - */ - int DLLEXPORT ENdeletelink(int linkIndex); - - - - - /*************************************************** - - Threadsafe versions of all epanet functions - - ***************************************************/ - int DLLEXPORT EN_alloc(EN_Project **p); - int DLLEXPORT EN_free(EN_Project *p); - int DLLEXPORT EN_epanet(char *inpFile, char *rptFile, char *binOutFile, void (*callback) (char *)); - int DLLEXPORT EN_init(EN_Project *p, char *rptFile, char *binOutFile, EN_FlowUnits UnitsType, EN_FormType HeadlossFormula); - int DLLEXPORT EN_open(EN_Project *p, char *inpFile, char *rptFile, char *binOutFile); - int DLLEXPORT EN_saveinpfile(EN_Project *p, char *filename); - int DLLEXPORT EN_close(EN_Project *p); - int DLLEXPORT EN_solveH(EN_Project *p); - int DLLEXPORT EN_saveH(EN_Project *p); - int DLLEXPORT EN_openH(EN_Project *p); - int DLLEXPORT EN_initH(EN_Project *p, int EN_SaveOption); - int DLLEXPORT EN_runH(EN_Project *p, long *currentTime); - int DLLEXPORT EN_nextH(EN_Project *p, long *tStep); - int DLLEXPORT EN_closeH(EN_Project *p); - int DLLEXPORT EN_savehydfile(EN_Project *p, char *filename); - int DLLEXPORT EN_usehydfile(EN_Project *p, char *filename); - int DLLEXPORT EN_solveQ(EN_Project *p); - int DLLEXPORT EN_openQ(EN_Project *p); - int DLLEXPORT EN_initQ(EN_Project *p, int saveFlag); - int DLLEXPORT EN_runQ(EN_Project *p, long *currentTime); - int DLLEXPORT EN_nextQ(EN_Project *p, long *tStep); - int DLLEXPORT EN_stepQ(EN_Project *p, long *timeLeft); - int DLLEXPORT EN_closeQ(EN_Project *p); - int DLLEXPORT EN_writeline(EN_Project *p, char *line); - int DLLEXPORT EN_report(EN_Project *p); - int DLLEXPORT EN_resetreport(EN_Project *p); - int DLLEXPORT EN_setreport(EN_Project *p, char *reportFormat); - int DLLEXPORT EN_getcontrol(EN_Project *p, int controlIndex, int *controlType, int *linkIndex, EN_API_FLOAT_TYPE *setting, int *nodeIndex, EN_API_FLOAT_TYPE *level); - int DLLEXPORT EN_getcount(EN_Project *p, EN_CountType code, int *count); - int DLLEXPORT EN_getoption(EN_Project *p, EN_Option opt, EN_API_FLOAT_TYPE *value); - int DLLEXPORT EN_gettimeparam(EN_Project *p, int code, long *value); - int DLLEXPORT EN_getflowunits(EN_Project *p, int *code); - int DLLEXPORT EN_setflowunits(EN_Project *p, int code); - int DLLEXPORT EN_getpatternindex(EN_Project *p, char *id, int *index); - int DLLEXPORT EN_getpatternid(EN_Project *p, int index, char *id); - int DLLEXPORT EN_getpatternlen(EN_Project *p, int index, int *len); - int DLLEXPORT EN_getpatternvalue(EN_Project *p, int index, int period, EN_API_FLOAT_TYPE *value); - int DLLEXPORT EN_getaveragepatternvalue(EN_Project *p, int index, EN_API_FLOAT_TYPE *value); - int DLLEXPORT EN_getqualtype(EN_Project *p, int *qualcode, int *tracenode); - int DLLEXPORT EN_geterror(int errcode, char *errmsg, int maxLen); - int DLLEXPORT EN_getstatistic(EN_Project *p, int code, EN_API_FLOAT_TYPE* value); - int DLLEXPORT EN_getnodeindex(EN_Project *p, char *id, int *index); - int DLLEXPORT EN_getnodeid(EN_Project *p, int index, char *id); - int DLLEXPORT EN_getnodetype(EN_Project *p, int index, int *code); - int DLLEXPORT EN_getnodevalue(EN_Project *p, int index, int code, EN_API_FLOAT_TYPE *value); - int DLLEXPORT EN_getcoord(EN_Project *p, int index, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y); - int DLLEXPORT EN_setcoord(EN_Project *p, int index, EN_API_FLOAT_TYPE x, EN_API_FLOAT_TYPE y); - int DLLEXPORT EN_getnumdemands(EN_Project *p, int nodeIndex, int *numDemands); - int DLLEXPORT EN_getbasedemand(EN_Project *p, int nodeIndex, int demandIndex, EN_API_FLOAT_TYPE *baseDemand); - int DLLEXPORT EN_getdemandpattern(EN_Project *p, int nodeIndex, int demandIndex, int *pattIndex); - int DLLEXPORT EN_getlinkindex(EN_Project *p, char *id, int *index); - int DLLEXPORT EN_getlinkid(EN_Project *p, int index, char *id); - int DLLEXPORT EN_getlinktype(EN_Project *p, int index, EN_LinkType *code); - int DLLEXPORT EN_setlinktype(EN_Project *p, char *id, EN_LinkType type); - int DLLEXPORT EN_getlinknodes(EN_Project *p, int index, int *node1, int *node2); - int DLLEXPORT EN_getlinkvalue(EN_Project *p, int index, EN_LinkProperty code, EN_API_FLOAT_TYPE *value); - int DLLEXPORT EN_getcurve(EN_Project *p, int curveIndex, char* id, int *nValues, EN_API_FLOAT_TYPE **xValues, EN_API_FLOAT_TYPE **yValues); - int DLLEXPORT EN_getheadcurveindex(EN_Project *p, int pumpIndex, int *curveIndex); - int DLLEXPORT EN_setheadcurveindex(EN_Project *p, int pumpIndex, int curveIndex); - int DLLEXPORT EN_getpumptype(EN_Project *p, int linkIndex, int *outType); - int DLLEXPORT EN_getversion(int *version); - int DLLEXPORT EN_setcontrol(EN_Project *p, int cindex, int ctype, int lindex, EN_API_FLOAT_TYPE setting, int nindex, EN_API_FLOAT_TYPE level); - int DLLEXPORT EN_setnodevalue(EN_Project *p, int index, int code, EN_API_FLOAT_TYPE v); - int DLLEXPORT EN_setlinkvalue(EN_Project *p, int index, int code, EN_API_FLOAT_TYPE v); - int DLLEXPORT EN_addpattern(EN_Project *p, char *id); - int DLLEXPORT EN_setpattern(EN_Project *p, int index, EN_API_FLOAT_TYPE *f, int len); - int DLLEXPORT EN_setpatternvalue(EN_Project *p, int index, int period, EN_API_FLOAT_TYPE value); - int DLLEXPORT EN_settimeparam(EN_Project *p, int code, long value); - int DLLEXPORT EN_setoption(EN_Project *p, int code, EN_API_FLOAT_TYPE v); - int DLLEXPORT EN_setstatusreport(EN_Project *p, int code); - int DLLEXPORT EN_setqualtype(EN_Project *p, int qualcode, char *chemname, char *chemunits, char *tracenode); - int DLLEXPORT EN_getqualinfo(EN_Project *p, int *qualcode, char *chemname, char *chemunits, int *tracenode); - int DLLEXPORT EN_setbasedemand(EN_Project *p, int nodeIndex, int demandIdx, EN_API_FLOAT_TYPE baseDemand); - int DLLEXPORT EN_getcurveindex(EN_Project *p, char *id, int *index); - int DLLEXPORT EN_getcurveid(EN_Project *p, int index, char *id); - int DLLEXPORT EN_getcurvelen(EN_Project *p, int index, int *len); - int DLLEXPORT EN_getcurvevalue(EN_Project *p, int curveIndex, int pointIndex, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y); - int DLLEXPORT EN_setcurvevalue(EN_Project *p, int curveIndex, int pointIndex, EN_API_FLOAT_TYPE x, EN_API_FLOAT_TYPE y); - int DLLEXPORT EN_setcurve(EN_Project *p, int index, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y, int len); - int DLLEXPORT EN_addcurve(EN_Project *p, char *id); - int DLLEXPORT EN_getrule(EN_Project *p, int index, int *nPremises, int *nTrueActions, int *nFalseActions, EN_API_FLOAT_TYPE *priority); - int DLLEXPORT EN_setrulepriority(EN_Project *p, int index, EN_API_FLOAT_TYPE priority); - int DLLEXPORT EN_getpremise(EN_Project *p, int indexRule, int indexPremise, int *logop, int *object, int *indexObj, int *variable, int *relop, int *status, EN_API_FLOAT_TYPE *value); - int DLLEXPORT EN_setpremise(EN_Project *p, int indexRule, int indexPremise, int logop, int object, int indexObj, int variable, int relop, int status, EN_API_FLOAT_TYPE value); - int DLLEXPORT EN_setpremiseindex(EN_Project *p, int indexRule, int indexPremise, int indexObj); - int DLLEXPORT EN_setpremisestatus(EN_Project *p, int indexRule, int indexPremise, int status); - int DLLEXPORT EN_setpremisevalue(EN_Project *p, int indexRule, int indexPremise, EN_API_FLOAT_TYPE value); - int DLLEXPORT EN_gettrueaction(EN_Project *p, int indexRule, int indexAction, int *indexLink, int *status, EN_API_FLOAT_TYPE *setting); - int DLLEXPORT EN_settrueaction(EN_Project *p, int indexRule, int indexAction, int indexLink, int status, EN_API_FLOAT_TYPE setting); - int DLLEXPORT EN_getfalseaction(EN_Project *p, int indexRule, int indexAction, int *indexLink, int *status, EN_API_FLOAT_TYPE *setting); - int DLLEXPORT EN_setfalseaction(EN_Project *p, int indexRule, int indexAction, int indexLink, int status, EN_API_FLOAT_TYPE setting); - int DLLEXPORT EN_getruleID(EN_Project *p, int indexRule, char* id); - int DLLEXPORT EN_addnode(EN_Project *p, char *id, EN_NodeType nodeType); - int DLLEXPORT EN_addlink(EN_Project *p, char *id, EN_LinkType linkType, char *fromNode, char *toNode); - int DLLEXPORT EN_deletenode(EN_Project *p, int nodeIndex); - int DLLEXPORT EN_deletelink(EN_Project *p, int linkIndex); - - - - - - -#if defined(__cplusplus) -} -#endif - -#endif //EPANET2_H diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c index a1bbca5..5a889f3 100644 --- a/src/hydcoeffs.c +++ b/src/hydcoeffs.c @@ -192,12 +192,18 @@ void linkcoeffs(EN_Project *pr) // Examine each link of network */ for (k = 1; k <= net->Nlinks; k++) { - if (sol->P[k] == 0.0) continue; +// if (sol->P[k] == 0.0) continue; link = &net->Link[k]; + switch (link->Type) { + case EN_PRV: + case EN_PSV: + case EN_FCV: + if (hyd->LinkSetting[k] != MISSING) continue; + } + n1 = link->N1; // Start node of link n2 = link->N2; // End node of link - // Update net nodal inflows (X), solution matrix (A) and RHS array (F) // (Use covention that flow out of node is (-), flow into node is (+)) hyd->X_tmp[n1] -= hyd->LinkFlows[k]; diff --git a/src/hydraul.c b/src/hydraul.c index 37664a8..545ed37 100755 --- a/src/hydraul.c +++ b/src/hydraul.c @@ -636,7 +636,8 @@ int controls(EN_Project *pr) double k1, k2; char s1, s2; Slink *link; - + Scontrol *control; + EN_Network *net = &pr->network; time_options_t *top = &pr->time_options; hydraulics_t *hyd = &pr->hydraulics; @@ -645,13 +646,13 @@ int controls(EN_Project *pr) setsum = 0; for (i=1; i <= net->Ncontrols; i++) { - Scontrol *control = &net->Control[i]; + control = &net->Control[i]; /* Make sure that link is defined */ reset = 0; - if ( (k = control->Link) <= 0) { - continue; - } - link = &net->Link[k]; + if ( (k = control->Link) <= 0) { + continue; + } + link = &net->Link[k]; /* Link is controlled by tank level */ if ((n = control->Node) > 0 && n > net->Njuncs) { diff --git a/src/hydsolver.c b/src/hydsolver.c index 9f789f1..a9cd94f 100644 --- a/src/hydsolver.c +++ b/src/hydsolver.c @@ -101,11 +101,6 @@ int hydsolve(EN_Project *pr, int *iter, double *relerr) nextcheck = hyd->CheckFreq; hyd->RelaxFactor = 1.0; - /* Compute initial head loss coefficients*/ - for (i = 1; i <= net->Nlinks; i++) { - hlosscoeff(pr, i); - } - /* Repeat iterations until convergence or trial limit is exceeded. */ /* (hyd->ExtraIter used to increase trials in case of status cycling.) */ if (pr->report.Statflag == FULL) { @@ -123,6 +118,9 @@ int hydsolve(EN_Project *pr, int *iter, double *relerr) ** head loss gradients, & F = flow correction terms. ** Solution for H is returned in F from call to linsolve(). */ + for (i = 1; i <= net->Nlinks; i++) { + hlosscoeff(pr, i); + } matrixcoeffs(pr); errcode = linsolve(&hyd->solver, net->Njuncs); @@ -147,16 +145,9 @@ int hydsolve(EN_Project *pr, int *iter, double *relerr) newerr = newflows(pr, &hydbal); /* Update flows */ *relerr = newerr; - /* Check hydraulic balance & re-compute head loss coefficients */ - checkhydbalance(pr, &hydbal); - for (i = 1; i <= net->Nlinks; i++) { - hlosscoeff(pr, i); - } - /* Write convergence error to status report if called for */ if (rep->Statflag == FULL) { writerelerr(pr, *iter, *relerr); - reporthydbal(pr, &hydbal); } /* Apply solution damping & check for change in valve status */ @@ -870,7 +861,7 @@ double newflows(EN_Project *pr, Hydbalance *hbal) dqsum = 0.0; hbal->maxflowchange = 0.0; - hbal->maxflowlink = -1; + hbal->maxflowlink = 1; /* Update flows in all links */ for (k = 1; k <= net->Nlinks; k++) @@ -988,8 +979,9 @@ void checkhydbalance(EN_Project *pr, Hydbalance *hbal) solver_t *sol = &hyd->solver; Slink *link; hbal->maxheaderror = 0.0; - hbal->maxheadlink = -1; + hbal->maxheadlink = 1; for (k = 1; k <= net->Nlinks; k++) { + hlosscoeff(pr, k); if (sol->P[k] == 0.0) continue; link = &net->Link[k]; n1 = link->N1; @@ -1019,6 +1011,10 @@ int hasconverged(EN_Project *pr, double *relerr, Hydbalance *hbal) hydraulics_t *hyd = &pr->hydraulics; if (*relerr > hyd->Hacc) return 0; + checkhydbalance(pr, hbal); + if (pr->report.Statflag == FULL) { + reporthydbal(pr, hbal); + } if (hyd->HeadErrorLimit > 0.0 && hbal->maxheaderror > hyd->HeadErrorLimit) return 0; if (hyd->FlowChangeLimit > 0.0 && diff --git a/src/main.c b/src/main.c deleted file mode 100644 index bda132f..0000000 --- a/src/main.c +++ /dev/null @@ -1,135 +0,0 @@ -#include -#include -#include "epanet2.h" - -#define MAXMSG 255 /* Max. # characters in message text */ -#define MAXWARNCODE 99 -/* text copied here, no more need of include "text.h" */ -#define FMT01 "\nEPANET Version %d.%d.%d" -#define FMT03 "\n Correct syntax is:\n %s \n" -#define FMT09 "\n\nEPANET completed." -#define FMT10 "\nEPANET completed. There are warnings." -#define FMT11 "\nEPANET completed. There are errors." - - -void writeConsole(char *s); - - -/* ----------------------------------------------------------------- - Entry point used to compile a stand-alone executable. ----------------------------------------------------------------- -*/ - - -int main(int argc, char *argv[]) -/*-------------------------------------------------------------- - ** Input: argc = number of command line arguments - ** *argv[] = array of command line arguments - ** Output: none - ** Purpose: main program segment - ** - ** Command line for stand-alone operation is: - ** progname f1 f2 f3 - ** where progname = name of executable this code was compiled to, - ** f1 = name of input file, - ** f2 = name of report file (optional, stdout if left blank) - ** f3 = name of binary output file (optional, nullfile if left blank). - **-------------------------------------------------------------- - */ -{ - char *f1,*f2,*f3; - char blank[] = ""; - char errmsg[MAXMSG+1]=""; - int errcode; - int version; - char s[256]; - int major; - int minor; - int patch; - - /* get version from DLL and trasform in Major.Minor.Patch format - instead of hardcoded version */ - ENgetversion(&version); - major= version/10000; - minor= (version%10000)/100; - patch= version%100; - sprintf(s,FMT01, major , minor, patch); - writeConsole(s); - - - /* Check for proper number of command line arguments */ - if (argc < 2) { - sprintf(s, FMT03, argv[0]); - writeConsole(s); - return(1); - } - - /* set inputfile name */ - f1 = argv[1]; - if (argc > 2) { - /* set rptfile name */ - f2 = argv[2]; - } - else { - /* use stdout for rptfile */ - f2 = blank; - } - if (argc > 3) { - /* set binary output file name */ - f3 = argv[3]; - } - else { - /* NO binary output*/ - f3 = blank; - } - - /* Call the main control function */ - if (strlen(f2)> 0) { - /* use stdout for progress messages */ - //errcode = ENepanet(f1,f2,f3,writeConsole); - errcode = ENepanet(f1, f2, f3, NULL); - } - else { - /* use stdout for reporting, no progress messages */ - errcode = ENepanet(f1,f2,f3,NULL); - } - - /* Error/Warning check */ - if (errcode == 0) { - /* no errors */ - writeConsole(FMT09); - return(0); - } - else { - if (errcode > MAXWARNCODE) printf("\n Fatal Error: "); - ENgeterror(errcode, errmsg, MAXMSG); - writeConsole(errmsg); - if (errcode > MAXWARNCODE) { - // error // - writeConsole(FMT11); - return(errcode); - } - else { - // warning // - writeConsole(FMT10); - return(0); - } - } - -} /* End of main */ - - -void writeConsole(char *s) -/*---------------------------------------------------------------- - ** Input: text string - ** Output: none - ** Purpose: writes string of characters to console - **---------------------------------------------------------------- - */ -{ - fprintf(stdout,"%s\n",s); - fflush(stdout); -} - - diff --git a/src/types.h b/src/types.h index 89a5467..3c26393 100755 --- a/src/types.h +++ b/src/types.h @@ -50,7 +50,6 @@ typedef int INT4; #define MAXMSG 79 /* Max. # characters in message text */ #define MAXLINE 255 /* Max. # characters read from input line */ #define MAXFNAME 259 /* Max. # characters in file name */ ->>>>>>> f11308fc72eccc8b74b18d5ddab0eb2ae0d36587 #define MAXTOKS 40 /* Max. items per line of input */ #define TZERO 1.E-4 /* Zero time tolerance */ #define TRUE 1 From 92e10a18b06f96a8da0dba43f222cdaa1219a74b Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Tue, 19 Jun 2018 10:41:30 -0400 Subject: [PATCH 6/9] Added new convergence options to epanet2.h --- include/epanet2.h | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/include/epanet2.h b/include/epanet2.h index ebd367b..81b469e 100644 --- a/include/epanet2.h +++ b/include/epanet2.h @@ -212,7 +212,9 @@ typedef enum { EN_ACCURACY = 1, EN_TOLERANCE = 2, EN_EMITEXPON = 3, - EN_DEMANDMULT = 4 + EN_DEMANDMULT = 4, + EN_HEADERROR = 5, + EN_FLOWCHANGE = 6 } EN_Option; typedef enum { From f86df2291a12003c9ba05f56fbc574db5bbe8f3c Mon Sep 17 00:00:00 2001 From: Elad Salomons Date: Thu, 21 Jun 2018 10:20:16 +0300 Subject: [PATCH 7/9] Saving HeadErrorLimit and FlowChangeLimit only if they are used To help a bit with backward computability the two new parameters will be saved to the INP file only if values were set for them. --- src/inpfile.c | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/inpfile.c b/src/inpfile.c index b6a5fe3..8bf894d 100644 --- a/src/inpfile.c +++ b/src/inpfile.c @@ -618,8 +618,12 @@ int saveinpfile(EN_Project *pr, char *fname) fprintf(f, "\n CHECKFREQ %-d", hyd->CheckFreq); fprintf(f, "\n MAXCHECK %-d", hyd->MaxCheck); fprintf(f, "\n DAMPLIMIT %-.8f", hyd->DampLimit); - fprintf(f, "\n HEADERROR %-.8f", hyd->HeadErrorLimit * pr->Ucf[HEAD]); - fprintf(f, "\n FLOWCHANGE %-.8f", hyd->FlowChangeLimit * pr->Ucf[FLOW]); + if (hyd->HeadErrorLimit > 0.0) { + fprintf(f, "\n HEADERROR %-.8f", hyd->HeadErrorLimit * pr->Ucf[HEAD]); + } + if (hyd->FlowChangeLimit > 0.0) { + fprintf(f, "\n FLOWCHANGE %-.8f", hyd->FlowChangeLimit * pr->Ucf[FLOW]); + } /* Write [REPORT] section */ From f5a1b9e5183da8daf2377538700d18b2e206153e Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Thu, 21 Jun 2018 12:05:51 -0400 Subject: [PATCH 8/9] Reverted to legacy Darcy-Weisbach method --- src/hydcoeffs.c | 154 +++++++++++++++++++++++++----------------------- src/hydsolver.c | 1 + 2 files changed, 82 insertions(+), 73 deletions(-) diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c index 5a889f3..9b830d9 100644 --- a/src/hydcoeffs.c +++ b/src/hydcoeffs.c @@ -41,7 +41,7 @@ static void emittercoeffs(EN_Project *pr); static void pipecoeff(EN_Project *pr, int k); static void DWpipecoeff(EN_Project *pr, int k); -static double frictionFactor(double q, double e, double s, double *dfdq); +static double frictionFactor(EN_Project *pr, int k, double *dfdq); static void pumpcoeff(EN_Project *pr, int k); static void curvecoeff(EN_Project *pr, int i, double q, double *h0, double *r); @@ -192,15 +192,8 @@ void linkcoeffs(EN_Project *pr) // Examine each link of network */ for (k = 1; k <= net->Nlinks; k++) { -// if (sol->P[k] == 0.0) continue; + if (sol->P[k] == 0.0) continue; link = &net->Link[k]; - switch (link->Type) { - case EN_PRV: - case EN_PSV: - case EN_FCV: - if (hyd->LinkSetting[k] != MISSING) continue; - } - n1 = link->N1; // Start node of link n2 = link->N2; // End node of link @@ -437,88 +430,103 @@ void DWpipecoeff(EN_Project *pr, int k) solver_t *sol = &hyd->solver; Slink *link = &pr->network.Link[k]; - double s = hyd->Viscos * link->Diam; double q = ABS(hyd->LinkFlows[k]); double dfdq = 0.0; - double e; // relative roughness height - double r; // total resistance coeff. - double f; // friction factor - double hloss; // head loss - double hgrad; // head loss gradient + double r, r1, f, ml, p, hloss; - // ... use Hagen-Poiseuille formula for laminar flow (Re <= 2000) - if (q <= A2 * s) + ml = link->Km; // Minor loss coeff. + r = link->R; // Resistance coeff. + f = frictionFactor(pr,k,&dfdq); // D-W friction factor + r1 = f*r+ml; + + // Use large P coefficient for small flow resistance product + if (r1*q < hyd->RQtol) { - r = 16.0 * PI * s * link->R; - hloss = q * (r + link->Km * q); - hgrad = r + 2.0 * link->Km * q; + sol->P[k] = 1.0/hyd->RQtol; + sol->Y[k] = hyd->LinkFlows[k]/hyd->Hexp; + return; } - // ... use Colebrook formula for turbulent flow - else - { - e = link->Kc / link->Diam; - f = frictionFactor(q, e, s, &dfdq); - r = f * link->R + link->Km; - hloss = r * q * hyd->LinkFlows[k]; - hgrad = (2.0 * r * q) + (dfdq * link->R * q * q); - } - - // ... head loss has same sign as flow - hloss *= SGN(hyd->LinkFlows[k]); - - // ... compute P and Y coeffs. - sol->P[k] = 1.0 / hgrad; - sol->Y[k] = hloss / hgrad; + // Compute P and Y coefficients + hloss = r1*SQR(q); // Total head loss + p = 2.0*r1*q; // |dHloss/dQ| + // + dfdq*r*q*q; // Ignore df/dQ term + p = 1.0 / p; + sol->P[k] = p; + sol->Y[k] = SGN(hyd->LinkFlows[k]) * hloss * p; } -double frictionFactor(double q, double e, double s, double *dfdq) +double frictionFactor(EN_Project *pr, int k, double *dfdq) /* **-------------------------------------------------------------- -** Input: q = flow rate (cfs) -** e = pipe roughness / diameter -** s = viscosity * diameter -** Output: f = friction factor -** dfdq = derivative of f w.r.t. flow +** Input: k = link index +** Output: returns friction factor and +** replaces dfdq (derivative of f w.r.t. flow) ** Purpose: computes Darcy-Weisbach friction factor and its ** derivative as a function of Reynolds Number (Re). +** +** Note: Current formulas for dfdq need to be corrected +** so dfdq returned as 0. **-------------------------------------------------------------- */ { - double f; - double x1, x2, x3, x4, - y1, y2, y3, - fa, fb, r; - double w = q / s; // Re*Pi/4 + double q, // Abs. value of flow + f; // Friction factor + double x1,x2,x3,x4, + y1,y2,y3, + fa,fb,r; + double s,w; + + hydraulics_t *hyd = &pr->hydraulics; + Slink *link = &pr->network.Link[k]; - // ... For Re >= 4000 use Colebrook Formula - if (w >= A1) - { - y1 = A8 / pow(w, 0.9); - y2 = e / 3.7 + y1; - y3 = A9 * log(y2); - f = 1.0 / (y3*y3); - *dfdq = 1.8 * f * y1 * A9 / y2 / y3 / q; - } + *dfdq = 0.0; + if (link->Type > EN_PIPE) + return(1.0); // Only apply to pipes + q = ABS(hyd->LinkFlows[k]); + s = hyd->Viscos * link->Diam; + w = q/s; // w = Re(Pi/4) + + // For Re >= 4000 use Colebrook Formula + if (w >= A1) + { + y1 = A8/pow(w,0.9); + y2 = link->Kc/(3.7*link->Diam) + y1; + y3 = A9*log(y2); + f = 1.0/SQR(y3); + /* *dfdq = (2.0+AA*y1/(y2*y3))*f; */ /* df/dq */ + } + + // For Re > 2000 use Interpolation Formula + else if (w > A2) + { + y2 = link->Kc/(3.7*link->Diam) + AB; + y3 = A9*log(y2); + fa = 1.0/SQR(y3); + fb = (2.0+AC/(y2*y3))*fa; + r = w/A2; + x1 = 7.0*fa - fb; + x2 = 0.128 - 17.0*fa + 2.5*fb; + x3 = -0.128 + 13.0*fa - (fb+fb); + x4 = r*(0.032 - 3.0*fa + 0.5*fb); + f = x1 + r*(x2 + r*(x3+x4)); + /* *dfdq = (x1 + x1 + r*(3.0*x2 + r*(4.0*x3 + 5.0*x4))); */ + } + + // For Re > 8 (Laminar Flow) use Hagen-Poiseuille Formula + else if (w > A4) + { + f = A3*s/q; // 16 * PI * Viscos * Diam / Flow + /* *dfdq = A3*s; */ + } + else + { + f = 8.0; + *dfdq = 0.0; + } + return(f); - // ... Use interpolating polynomials developed by - // E. Dunlop for transition flow from 2000 < Re < 4000. - else - { - y2 = e / 3.7 + AB; - y3 = A9 * log(y2); - fa = 1.0 / (y3*y3); - fb = (2.0 + AC / (y2*y3)) * fa; - r = w / A2; - x1 = 7.0 * fa - fb; - x2 = 0.128 - 17.0 * fa + 2.5 * fb; - x3 = -0.128 + 13.0 * fa - (fb + fb); - x4 = r * (0.032 - 3.0 * fa + 0.5 *fb); - f = x1 + r * (x2 + r * (x3 + x4)); - *dfdq = (x2 + 2.0 * r * (x3 + x4)) / s / A2; - } - return f; } diff --git a/src/hydsolver.c b/src/hydsolver.c index a9cd94f..0adf780 100644 --- a/src/hydsolver.c +++ b/src/hydsolver.c @@ -981,6 +981,7 @@ void checkhydbalance(EN_Project *pr, Hydbalance *hbal) hbal->maxheaderror = 0.0; hbal->maxheadlink = 1; for (k = 1; k <= net->Nlinks; k++) { + if (hyd->LinkStatus[k] <= CLOSED) continue; hlosscoeff(pr, k); if (sol->P[k] == 0.0) continue; link = &net->Link[k]; From 39dc089b58fd064862bacc4bcf12aa2201b62f66 Mon Sep 17 00:00:00 2001 From: Michael Tryby Date: Tue, 26 Jun 2018 16:36:23 -0400 Subject: [PATCH 9/9] Updating benchmark --- tools/before-test.cmd | 4 ++-- tools/run-nrtest.cmd | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/tools/before-test.cmd b/tools/before-test.cmd index 0e8b41a..8871b44 100644 --- a/tools/before-test.cmd +++ b/tools/before-test.cmd @@ -25,8 +25,8 @@ set SCRIPT_HOME=%~dp0 set TEST_HOME=%~1 -set EXAMPLES_VER=1.0.1 -set BENCHMARK_VER=2012vs10 +set EXAMPLES_VER=1.0.2-dev +set BENCHMARK_VER=220dev-vs17 set TESTFILES_URL=https://github.com/OpenWaterAnalytics/epanet-example-networks/archive/v%EXAMPLES_VER%.zip diff --git a/tools/run-nrtest.cmd b/tools/run-nrtest.cmd index 0b7c17f..f5d4c8d 100644 --- a/tools/run-nrtest.cmd +++ b/tools/run-nrtest.cmd @@ -18,7 +18,7 @@ setlocal set NRTEST_SCRIPT_PATH=%~1 set TEST_SUITE_PATH=%~2 -set BENCHMARK_VER=2012vs10 +set BENCHMARK_VER=220dev-vs17 set NRTEST_EXECUTE_CMD=python %NRTEST_SCRIPT_PATH%\nrtest execute