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

Fixed water quality mass balance issue (#160)

Browse Source
This commit is contained in:
Lew Rossman
2018-10-09 12:53:20 -04:00
parent 4848f692f6
commit 7c021cf533
17 changed files with 2262 additions and 1884 deletions
Download Patch File Download Diff File Expand all files Collapse all files

51
ReleaseNotes2_2.md
View File

@@ -6,30 +6,30 @@ This document describes the changes and updates that have been made to version 2
## Thread-Safe API Functions
A duplicate set of the version 2.1 API functions has been provided that allow multiple EPANET projects to be analyzed concurrently in a thread-safe manner. These functions maintain the same name as the original but use a `EN_` prefix instead of `EN`. In addition, the first argument to each of these functions is a pointer to an `EN_Project` structure that encapsulates the network data for the particular project being analyzed. For example, instead of writing:
A duplicate set of the version 2.1 API functions has been provided that allow multiple EPANET projects to be analyzed concurrently in a thread-safe manner. These functions maintain the same name as the original but use a `EN_` prefix instead of `EN`. In addition, the first argument to each of these functions is a handle that identifies the network data for the particular project being analyzed. For example, instead of writing:
`ENgetnodevalue(nodeIndex, EN_ELEVATION, &elev)`
one would use:
`EN_getnodevalue(pr, nodeIndex, EN_ELEVATION, &elev)`
`EN_getnodevalue(ph, nodeIndex, EN_ELEVATION, &elev)`
where `pr` is the pointer to an `EN_Project`.
where `ph` is the handle assigned to the project.
Two new functions have been added to the API to manage the creation and deletion of project pointers. `EN_createproject` creates a new project along with a pointer to it, while `EN_deleteproject` deletes a project. An example of using the thread-safe version of the API is shown below:
Two new functions have been added to the API to manage the creation and deletion of project handles. `EN_createproject` creates a new project along with its handle, while `EN_deleteproject` deletes a project. An example of using the thread-safe version of the API is shown below:
```
#include "epanet2.h"
int runEpanet(char *finp, char *frpt)
{
EN_Project *pr = NULL;
EN_ProjectHandle ph = 0;
int err;
err = EN_createproject(&pr);
err = EN_createproject(&ph);
if (err) return err;
err = EN_open(pr, finp, frpt, "");
if (!err) err = EN_solveH(pr);
if (!err) err = EN_report(pr);
EN_close(pr);
EN_deleteproject(pr);
err = EN_open(ph, finp, frpt, "");
if (!err) err = EN_solveH(ph);
if (!err) err = EN_report(ph);
EN_close(ph);
EN_deleteproject(ph);
return err;
}
```
@@ -49,6 +49,11 @@ EPANET's hydraulic solver requires solving a system of linear equations over a s
EPANET's original node re-ordering scheme has been replaced by the more powerful **Multiple Minimum Degree (MMD)** algorithm. On a series of eight networks ranging in size from 7,700 to 100,000 nodes **MMD** reduced the solution time for a single period (steady state) hydraulic analysis by an average of 58%.
## Improved Handling of Near-Zero Flows
EPANET's hydraulic solver can generate an ill-conditioned solution matrix when pipe flows approach zero unless some adjustment is made (i.e., as a pipe's flow approaches 0 its head loss gradient also approaches 0 causing its reciprocal, which is used to form the solution matrix's coefficients, to approach infinity). EPANET 2.0 used an arbitrary cutoff on head loss gradient to prevent it from becoming 0. This approach doesn't allow a pipe to follow any head loss v. flow relation in the region below the cutoff and can produce incorrect solutions for some networks (see [Estrada et al., 2009](https://ascelibrary.org/doi/full/10.1061/%28ASCE%29IR.1943-4774.0000100)).
The hydraulic solver has been modified to use a linear head loss v. flow relation for flows approaching zero. For the Darcy-Weisbach equation, the linear Hagen-Poiseuille formula is used for laminar flow where the Reynolds Number is <= 2000. For the Hazen-Williams and Chezy-Manning equations, a flow limit `Qa` is established for each pipe, equal to the flow that produces the EPANET 2 gradient cutoff. For flows below this a linear head loss relation is used between 0 and the head loss at `Qa` and the gradient always equals the cutoff. EPANET 2.2 is now able to correctly solve the examples presented in Estrada et al. (2009) as well as those in [Gorev et al., (2013)](https://ascelibrary.org/doi/10.1061/%28ASCE%29HY.1943-7900.0000694) and [Elhay and Simpson (2011)](https://ascelibrary.org/doi/10.1061/%28ASCE%29HY.1943-7900.0000411).
## Pressure Dependent Demands
EPANET has always employed a Demand Driven Analysis (**DDA**) when modeling network hydraulics. Under this approach nodal demands at a given point in time are fixed values that must be delivered no matter what nodal heads and link flows are produced by a hydraulic solution. This can result in situations where required demands are satisfied at nodes that have negative pressures - a physical impossibility.
@@ -81,9 +86,26 @@ int EN_getdemandmodel(EN_Project *pr, int *modelType, double *pMin, double *pReq
for the thread-safe API. Some additional points regarding the new **PDA** option are:
- If no DEMAND MODEL and its parameters are specified then the analysis defaults to being demand driven (**DDA**).
- This implementation of **PDA** assumes that the same parameters apply to all nodes in the network. Extending the framework to allow different parameters for specific nodes is straightforward to do but is left as a future feature to implement.
- This implementation of **PDA** assumes that the same parameters apply to all nodes in the network. Extending the framework to allow different parameters for specific nodes is left as a future feature to implement.
- *P<sub>min</sub>* is allowed to equal to *P<sub>req</sub>*. This condition can be used to find a solution that results in the smallest amount of demand reductions needed to insure that no node delivers positive demand at a pressure below *P<sub>min</min>*.
## Improved Water Quality Mass Balance
As described by [Davis et al. (2018)](https://www.drink-water-eng-sci.net/11/25/2018/dwes-11-25-2018.pdf) EPANET's water quality simulations can result in some significant mass balance errors when modeling short term mass injections (errors are much smaller for continuous source flows). The entire water quality engine has been re-written to eliminate these errors. It still uses the Lagrangian Time Driven transport method but now analyzes each network node in topologically sorted order rather than in arbitrary order.
A Mass Balance Report now appears the end of a simulation's Status Report that lists the various components (inflow, outflow, reaction) that comprise the network's overall mass balance. In addition `EN_MASSBALANCE` can be used as a parameter in the `ENgetstatistic` (or `EN_getstatistic`) function to retrieve the Mass Balance Ratio (Total Outflow Mass / Total Inflow Mass) at any point during a water quality simulation.
Mass balance ratio (MBR) results for two of the networks analyzed by Davis et al. (2018) are shown in the following table. MBR-2.0 is for EPANET 2.0.012 as reported by Davis et al. while MBR-2.2 is for the re-written quality engine.
| Network | Time Step (s) | MBR-2.0 | MBR-2.2 |
|--|--|--|--|
| N2 | 900 | 16.63 | 1.00 |
| | 300 | 23.45 | 1.00 |
| | 60 | 6.49 | 1.00 |
| N4 | 900 | 0.09 | 1.00 |
| | 300 | 0.70 | 1.00 |
| | 60 | 0.98 | 1.00 |
Both network files are available [here](https://doi.org/10.23719/1375314).
## Code Changes
@@ -96,6 +118,10 @@ for the thread-safe API. Some additional points regarding the new **PDA** option
- `hydcoeffs.c` computes values of the matrix coefficients (derived from link head losses and their gradients) used by the hydraulic solver.
- `hydstatus.c` checks for status changes in valves and pumps as requested by the hydraulic solver.
- The Multiple Minimum Degree re-ordering algorithm appears in a new file named `genmmd.c`. This is 1990's legacy code that is readily available on the web and can be found in several linear equation solver libraries.
- The original `quality.c` file has also been split into three separate files:
- `quality.c` initializes the quality solver and supervises the quality calculations over each simulation time step.
- `qualreact.c` reacts the quality constituent within each pipe and tank over a single time step and also implements the various tank mixing models.
- `qualroute.c` topologically sorts the network's nodes when flow directions change and implements the Lagrangian Time Driven transport algorithm over a single time step.
## General changes
- Read and write demand categories names
@@ -152,6 +178,7 @@ for the thread-safe API. Some additional points regarding the new **PDA** option
- `EN_HEADERROR`
- `EN_FLOWCHANGE`
- `EN_DEMANDDEFPAT`
- `EN_MASSBALANCE`
### Curve types:
- `EN_V_CURVE`
- `EN_P_CURVE`

1
include/epanet2.bas
View File

@@ -80,6 +80,7 @@ Public Const EN_ITERATIONS = 0
Public Const EN_RELATIVEERROR = 1
Public Const EN_MAXHEADERROR = 2
Public Const EN_MAXFLOWCHANGE = 3
Public Const EN_MASSBALANCE = 4
Public Const EN_NODECOUNT = 0 'Component counts
Public Const EN_TANKCOUNT = 1

3
include/epanet2.h
View File

@@ -148,7 +148,8 @@ typedef enum {
EN_ITERATIONS = 0,
EN_RELATIVEERROR = 1,
EN_MAXHEADERROR = 2,
EN_MAXFLOWCHANGE = 3
EN_MAXFLOWCHANGE = 3,
EN_MASSBALANCE = 4
} EN_AnalysisStatistic;
typedef enum {

10
include/epanet2.vb
View File

@@ -72,6 +72,16 @@ Public Const EN_RULESTEP = 7
Public Const EN_STATISTIC = 8
Public Const EN_PERIODS = 9
Public Const EN_STARTTIME = 10 'ES
Public Const EN_HTIME = 11
Public Const EN_QTIME = 12
Public Const EN_HALTFLAG = 13
Public Const EN_NEXTEVENT = 14
Public Const EN_ITERATIONS = 0
Public Const EN_RELATIVEERROR = 1
Public Const EN_MAXHEADERROR = 2
Public Const EN_MAXFLOWCHANGE = 3
Public Const EN_MASSBALANCE = 4
Public Const EN_NODECOUNT = 0 'Component counts
Public Const EN_TANKCOUNT = 1

46
run/main.c
View File

@@ -7,9 +7,9 @@
/* text copied here, no more need of include "text.h" */
#define FMT01 "\nEPANET Version %d.%d.%d\n"
#define FMT03 "\nUsage:\n %s <input_filename> <report_filename> [<binary_filename>]\n"
#define FMT09 "\n\nEPANET completed."
#define FMT10 "\nEPANET completed. There are warnings."
#define FMT11 "\nEPANET completed. There are errors."
#define FMT09 "\n\nEPANET completed.\n"
#define FMT10 "\nEPANET completed. There are warnings.\n"
#define FMT11 "\nEPANET completed. There are errors.\n"
void writeConsole(char *s);
@@ -43,7 +43,6 @@ int main(int argc, char *argv[])
char errmsg[MAXMSG+1]="";
int errcode;
int version;
char s[256];
int major;
int minor;
int patch;
@@ -54,14 +53,11 @@ int main(int argc, char *argv[])
major= version/10000;
minor= (version%10000)/100;
patch= version%100;
sprintf(s,FMT01, major , minor, patch);
writeConsole(s);
printf(FMT01, major, minor, patch);
/* Check for proper number of command line arguments */
if (argc < 2) {
sprintf(s, FMT03, argv[0]);
writeConsole(s);
printf(FMT03, argv[0]);
return(1);
}
@@ -85,51 +81,27 @@ int main(int argc, char *argv[])
}
/* 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);
printf(FMT09);
return(0);
}
else {
if (errcode > MAXWARNCODE) printf("\n Fatal Error: ");
ENgeterror(errcode, errmsg, MAXMSG);
writeConsole(errmsg);
printf("%s\n", errmsg);
if (errcode > MAXWARNCODE) {
// error //
writeConsole(FMT11);
printf(FMT11);
return(errcode);
}
else {
// warning //
writeConsole(FMT10);
printf(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);
}

11
src/epanet.c
View File

@@ -1997,6 +1997,9 @@ int DLLEXPORT EN_getstatistic(EN_ProjectHandle ph, int code, EN_API_FLOAT_TYPE *
case EN_MAXFLOWCHANGE:
*value = (EN_API_FLOAT_TYPE)(p->hydraulics.MaxFlowChange * p->Ucf[FLOW]);
break;
case EN_MASSBALANCE:
*value = (EN_API_FLOAT_TYPE)(p->quality.massbalance.ratio);
break;
default:
break;
}
@@ -3759,7 +3762,7 @@ int DLLEXPORT EN_setoption(EN_ProjectHandle ph, int code, EN_API_FLOAT_TYPE v)
}
else
{
error = EN_getpatternid(p, value, &*tmpId);
error = EN_getpatternid(p, (int)value, tmpId);
if (error != 0)
return set_error(p->error_handle, error);
}
@@ -3768,12 +3771,12 @@ int DLLEXPORT EN_setoption(EN_ProjectHandle ph, int code, EN_API_FLOAT_TYPE v)
Snode *node = &net->Node[i];
for (demand = node->D; demand != NULL; demand = demand->next) {
if (demand->Pat == tmpPat) {
demand->Pat = value;
demand->Pat = (int)value;
}
}
}
strncpy(p->parser.DefPatID, tmpId, MAXID);
hyd->DefPat = value;
hyd->DefPat = (int)value;
break;
default:
@@ -4952,7 +4955,7 @@ int DLLEXPORT EN_addnode(EN_ProjectHandle ph, char *id, EN_NodeType nodeType) {
Snode *node;
Scoord *coord;
Scontrol *control;
rules_t *rule;
// rules_t *rule;
Premise *pchain, *pnext;

36
src/funcs.h
View File

@@ -127,8 +127,6 @@ void freerules(EN_Project *pr); /* Frees rule base memory
int writeRuleinInp(EN_Project *pr, FILE *f, /* Writes rule to an INP file */
int RuleIdx);
int writeRuleinInp(EN_Project *pr, FILE *f, int RuleIdx);
/* ------------- REPORT.C --------------*/
int writereport(EN_Project *pr); /* Writes formatted report */
void writelogo(EN_Project *pr); /* Writes program logo */
@@ -144,6 +142,7 @@ void writecontrolaction(EN_Project *pr, int, int); /* Writes control acti
void writeruleaction(EN_Project *pr, int, char *); /* Writes rule action taken */
int writehydwarn(EN_Project *pr, int,double); /* Writes hydraulic warnings */
void writehyderr(EN_Project *pr, int); /* Writes hydraulic error msg.*/
void writemassbalance(EN_Project *pr); // Writes mass balance ratio
int disconnected(EN_Project *pr); /* Checks for disconnections */
void marknodes(EN_Project *pr, int, int *, char *); /* Identifies connected nodes */
void getclosedlink(EN_Project *pr, int, char *); /* Finds a disconnecting link */
@@ -192,41 +191,12 @@ int linsolve(EN_Project *pr, int); /* Solves set of linear eqns
/* ----------- QUALITY.C ---------------*/
int openqual(EN_Project *pr); /* Opens WQ solver system */
void initqual(EN_Project *pr); /* Initializes WQ solver */
int initqual(EN_Project *pr); /* Initializes WQ solver */
int runqual(EN_Project *pr, long *); /* Gets current WQ results */
int nextqual(EN_Project *pr, long *); /* Updates WQ by hyd.timestep */
int stepqual(EN_Project *pr, long *); /* Updates WQ by WQ time step */
int closequal(EN_Project *pr); /* Closes WQ solver system */
int gethyd(EN_Project *pr, long *, long *); /* Gets next hyd. results */
char setReactflag(EN_Project *pr); /* Checks for reactive chem. */
void transport(EN_Project *pr, long); /* Transports mass in network */
void initsegs(EN_Project *pr); /* Initializes WQ segments */
void reorientsegs(EN_Project *pr); /* Re-orients WQ segments */
void updatesegs(EN_Project *pr, long); /* Updates quality in segments*/
void removesegs(EN_Project *pr, int); /* Removes a WQ segment */
void addseg(EN_Project *pr, int,double,double); /* Adds a WQ segment to pipe */
void accumulate(EN_Project *pr, long); /* Sums mass flow into node */
void updatenodes(EN_Project *pr, long); /* Updates WQ at nodes */
void sourceinput(EN_Project *pr, long); /* Computes source inputs */
void release(EN_Project *pr, long); /* Releases mass from nodes */
void updatetanks(EN_Project *pr, long); /* Updates WQ in tanks */
void updatesourcenodes(EN_Project *pr, long); /* Updates WQ at source nodes */
void tankmix1(EN_Project *pr, int, long); /* Complete mix tank model */
void tankmix2(EN_Project *pr, int, long); /* 2-compartment tank model */
void tankmix3(EN_Project *pr, int, long); /* FIFO tank model */
void tankmix4(EN_Project *pr, int, long); /* LIFO tank model */
double sourcequal(EN_Project *pr, Psource); /* Finds WQ input from source */
double avgqual(EN_Project *pr, int); /* Finds avg. quality in pipe */
void ratecoeffs(EN_Project *pr); /* Finds wall react. coeffs. */
double piperate(EN_Project *pr, int); /* Finds wall react. coeff. */
double pipereact(EN_Project *pr, int,double,
double,long); /* Reacts water in a pipe */
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 wallrate(EN_Project *pr, double,double,
double,double); /* Finds wall reaction rate */
/* ------------ OUTPUT.C ---------------*/
int savenetdata(EN_Project *pr); /* Saves basic data to file */
@@ -234,7 +204,7 @@ int savehyd(EN_Project *pr, long *); /* Saves hydraulic solution
int savehydstep(EN_Project *pr, long *); /* Saves hydraulic timestep */
int saveenergy(EN_Project *pr); /* Saves energy usage */
int readhyd(EN_Project *pr, long *); /* Reads hydraulics from file */
int readhydstep(FILE *hydFile, long *); /* Reads time step from file */
int readhydstep(EN_Project *pr, long *); /* Reads time step from file */
int saveoutput(EN_Project *pr); /* Saves results to file */
int nodeoutput(EN_Project *pr, int, REAL4 *,
double); /* Saves node results to file */

17
src/hydraul.c
View File

@@ -276,20 +276,21 @@ int nexthyd(EN_Project *pr, long *tstep)
addenergy(pr,hydstep);
}
/* Update current time. */
if (top->Htime < top->Dur) /* More time remains */
/* More time remains - update current time. */
if (top->Htime < top->Dur)
{
top->Htime += hydstep;
if (top->Htime >= top->Rtime) {
top->Rtime += top->Rstep;
if (!pr->quality.OpenQflag)
{
if (top->Htime >= top->Rtime) top->Rtime += top->Rstep;
}
}
/* No more time remains - force completion of analysis. */
else
{
top->Htime++; /* Force completion of analysis */
if (pr->quality.OpenQflag) {
pr->quality.Qtime++; // force completion of wq analysis too
}
top->Htime++;
if (pr->quality.OpenQflag) pr->quality.Qtime++;
}
*tstep = hydstep;
return(errcode);

1
src/input1.c
View File

@@ -168,6 +168,7 @@ void setdefaults(EN_Project *pr)
hyd->CheckFreq = CHECKFREQ;
hyd->MaxCheck = MAXCHECK;
hyd->DampLimit = DAMPLIMIT;
qu->massbalance.ratio = 0.0;
} /* End of setdefaults */
void initreport(report_options_t *r)

6
src/output.c
View File

@@ -386,7 +386,7 @@ int readhyd(EN_Project *pr, long *hydtime)
return result;
} /* End of readhyd */
int readhydstep(FILE *hydFile, long *hydstep)
int readhydstep(EN_Project *pr, long *hydstep)
/*
**--------------------------------------------------------------
** Input: none
@@ -397,8 +397,8 @@ int readhydstep(FILE *hydFile, long *hydstep)
*/
{
INT4 t;
if (fread(&t, sizeof(INT4), 1, hydFile) < 1)
return (0);
FILE *hydFile = pr->out_files.HydFile;
if (fread(&t, sizeof(INT4), 1, hydFile) < 1) return (0);
*hydstep = t;
return (1);
} /* End of readhydstep */

2106
src/quality.c
View File

File diff suppressed because it is too large Load Diff

721
src/qualreact.c Normal file
View File

@@ -0,0 +1,721 @@
/*
*********************************************************************
QUALREACT.C -- water quality reaction module for the EPANET program
*********************************************************************
*/
#include <stdio.h>
#include <math.h>
#include "types.h"
// Exported Functions
char setreactflag(EN_Project *pr);
double getucf(double);
void ratecoeffs(EN_Project *pr);
void reactpipes(EN_Project *pr, long);
void reacttanks(EN_Project *pr, long);
double mixtank(EN_Project *pr, int, double, double ,double);
// Imported Functions
extern void addseg(EN_Project *pr, int, double, double);
// Local Functions
static double piperate(EN_Project *pr, int);
static double pipereact(EN_Project *pr, int, double, double, long);
static double tankreact(EN_Project *pr, double, double, double, long);
static double bulkrate(EN_Project *pr, double, double, double);
static double wallrate(EN_Project *pr, double, double, double, double);
static void tankmix1(EN_Project *pr, int, double, double, double);
static void tankmix2(EN_Project *pr, int, double, double, double);
static void tankmix3(EN_Project *pr, int, double, double, double);
static void tankmix4(EN_Project *pr, int, double, double, double);
char setreactflag(EN_Project *pr)
/*
**-----------------------------------------------------------
** Input: none
** Output: returns 1 for reactive WQ constituent, 0 otherwise
** Purpose: checks if reactive chemical being simulated
**-----------------------------------------------------------
*/
{
int i;
EN_Network *net = &pr->network;
if (pr->quality.Qualflag == TRACE) return 0;
else if (pr->quality.Qualflag == AGE) return 1;
else
{
for (i = 1; i <= net->Nlinks; i++)
{
if (net->Link[i].Type <= EN_PIPE)
{
if (net->Link[i].Kb != 0.0 || net->Link[i].Kw != 0.0) return 1;
}
}
for (i = 1; i <= net->Ntanks; i++)
{
if (net->Tank[i].Kb != 0.0) return 1;
}
}
return 0;
}
double getucf(double order)
/*
**--------------------------------------------------------------
** Input: order = bulk reaction order
** Output: returns a unit conversion factor
** Purpose: converts bulk reaction rates from per Liter to
** per FT3 basis
**--------------------------------------------------------------
*/
{
if (order < 0.0) order = 0.0;
if (order == 1.0) return (1.0);
else return (1. / pow(LperFT3, (order - 1.0)));
}
void ratecoeffs(EN_Project *pr)
/*
**--------------------------------------------------------------
** Input: none
** Output: none
** Purpose: determines wall reaction coeff. for each pipe
**--------------------------------------------------------------
*/
{
int k;
double kw;
EN_Network *net = &pr->network;
quality_t *qual = &pr->quality;
for (k = 1; k <= net->Nlinks; k++)
{
kw = net->Link[k].Kw;
if (kw != 0.0) kw = piperate(pr, k);
net->Link[k].Rc = kw;
qual->PipeRateCoeff[k] = 0.0;
}
}
void reactpipes(EN_Project *pr, long dt)
/*
**--------------------------------------------------------------
** Input: dt = time step
** Output: none
** Purpose: reacts water within each pipe over a time step.
**--------------------------------------------------------------
*/
{
int k;
Pseg seg;
double cseg, rsum, vsum;
EN_Network *net = &pr->network;
quality_t *qual = &pr->quality;
// Examine each link in network
for (k = 1; k <= net->Nlinks; k++)
{
// Skip non-pipe links (pumps & valves)
if (net->Link[k].Type != EN_PIPE) continue;
rsum = 0.0;
vsum = 0.0;
// Examine each segment of the pipe
seg = qual->FirstSeg[k];
while (seg != NULL)
{
// React segment over time dt
cseg = seg->c;
seg->c = pipereact(pr, k, seg->c, seg->v, dt);
// Update reaction component of mass balance
qual->massbalance.reacted += (cseg - seg->c) * seg->v;
// Accumulate volume-weighted reaction rate
if (qual->Qualflag == CHEM)
{
rsum += fabs(seg->c - cseg) * seg->v;
vsum += seg->v;
}
seg = seg->prev;
}
// Normalize volume-weighted reaction rate
if (vsum > 0.0) qual->PipeRateCoeff[k] = rsum / vsum / dt * SECperDAY;
else qual->PipeRateCoeff[k] = 0.0;
}
}
void reacttanks(EN_Project *pr, long dt)
/*
**--------------------------------------------------------------
** Input: dt = time step
** Output: none
** Purpose: reacts water within each tank over a time step.
**--------------------------------------------------------------
*/
{
int i, k;
double c;
Pseg seg;
EN_Network *net = &pr->network;
quality_t *qual = &pr->quality;
Stank *tank;
// Examine each tank in network
for (i = 1; i <= net->Ntanks; i++)
{
// Skip reservoirs
tank = &net->Tank[i];
if (tank->A == 0.0) continue;
// k is segment chain belonging to tank i
k = net->Nlinks + i;
// React each volume segment in the chain
seg = qual->FirstSeg[k];
while (seg != NULL)
{
c = seg->c;
seg->c = tankreact(pr, seg->c, seg->v, tank->Kb, dt);
qual->massbalance.reacted += (c - seg->c) * seg->v;
seg = seg->prev;
}
}
}
double piperate(EN_Project *pr, int k)
/*
**--------------------------------------------------------------
** Input: k = link index
** Output: returns reaction rate coeff. for 1st-order wall
** reactions or mass transfer rate coeff. for 0-order
** reactions
** Purpose: finds wall reaction rate coeffs.
**--------------------------------------------------------------
*/
{
double a, d, u, q, kf, kw, y, Re, Sh;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
d = net->Link[k].Diam; // Pipe diameter, ft
// Ignore mass transfer if Schmidt No. is 0
if (qual->Sc == 0.0)
{
if (qual->WallOrder == 0.0) return BIG;
else return (net->Link[k].Kw * (4.0 / d) / pr->Ucf[ELEV]);
}
// Compute Reynolds No.
// Flow rate made consistent with how its saved to hydraulics file
q = (hyd->LinkStatus[k] <= CLOSED) ? 0.0 : hyd->LinkFlows[k];
a = PI * d * d / 4.0; // pipe area
u = fabs(q) / a; // flow velocity
Re = u * d / hyd->Viscos; // Reynolds number
// Compute Sherwood No. for stagnant flow
// (mass transfer coeff. = Diffus./radius)
if (Re < 1.0) Sh = 2.0;
// Compute Sherwood No. for turbulent flow using the Notter-Sleicher formula.
else if (Re >= 2300.0) Sh = 0.0149 * pow(Re, 0.88) * pow(qual->Sc, 0.333);
// Compute Sherwood No. for laminar flow using Graetz solution formula.
else
{
y = d / net->Link[k].Len * Re * qual->Sc;
Sh = 3.65 + 0.0668 * y / (1.0 + 0.04 * pow(y, 0.667));
}
// Compute mass transfer coeff. (in ft/sec)
kf = Sh * qual->Diffus / d;
// For zero-order reaction, return mass transfer coeff.
if (qual->WallOrder == 0.0) return kf;
// For first-order reaction, return apparent wall coeff.
kw = net->Link[k].Kw / pr->Ucf[ELEV]; // Wall coeff, ft/sec
kw = (4.0 / d) * kw * kf / (kf + fabs(kw)); // Wall coeff, 1/sec
return kw;
}
double pipereact(EN_Project *pr, int k, double c, double v, long dt)
/*
**------------------------------------------------------------
** Input: k = link index
** c = current quality in segment
** v = segment volume
** dt = time step
** Output: returns new WQ value
** Purpose: computes new quality in a pipe segment after
** reaction occurs
**------------------------------------------------------------
*/
{
double cnew, dc, dcbulk, dcwall, rbulk, rwall;
EN_Network *net = &pr->network;
quality_t *qual = &pr->quality;
// For water age (hrs), update concentration by timestep
if (qual->Qualflag == AGE)
{
dc = (double)dt / 3600.0;
cnew = c + dc;
cnew = MAX(0.0, cnew);
return cnew;
}
// Otherwise find bulk & wall reaction rates
rbulk = bulkrate(pr, c, net->Link[k].Kb, qual->BulkOrder) * qual->Bucf;
rwall = wallrate(pr, c, net->Link[k].Diam, net->Link[k].Kw, net->Link[k].Rc);
// Find change in concentration over timestep
dcbulk = rbulk * (double)dt;
dcwall = rwall * (double)dt;
// Update cumulative mass reacted
if (pr->time_options.Htime >= pr->time_options.Rstart)
{
qual->Wbulk += fabs(dcbulk) * v;
qual->Wwall += fabs(dcwall) * v;
}
// Update concentration
dc = dcbulk + dcwall;
cnew = c + dc;
cnew = MAX(0.0, cnew);
return cnew;
}
double tankreact(EN_Project *pr, double c, double v, double kb, long dt)
/*
**-------------------------------------------------------
** Input: c = current quality in tank
** v = tank volume
** kb = reaction coeff.
** dt = time step
** Output: returns new WQ value
** Purpose: computes new quality in a tank after
** reaction occurs
**-------------------------------------------------------
*/
{
double cnew, dc, rbulk;
quality_t *qual = &pr->quality;
// For water age, update concentration by timestep
if (qual->Qualflag == AGE)
{
dc = (double)dt / 3600.0;
}
// For chemical analysis apply bulk reaction rate
else
{
// Find bulk reaction rate
rbulk = bulkrate(pr, c, kb, qual->TankOrder) * qual->Tucf;
// Find concentration change & update quality
dc = rbulk * (double)dt;
if (pr->time_options.Htime >= pr->time_options.Rstart)
{
qual->Wtank += fabs(dc) * v;
}
}
cnew = c + dc;
cnew = MAX(0.0, cnew);
return cnew;
}
double bulkrate(EN_Project *pr, double c, double kb, double order)
/*
**-----------------------------------------------------------
** Input: c = current WQ concentration
** kb = bulk reaction coeff.
** order = bulk reaction order
** Output: returns bulk reaction rate
** Purpose: computes bulk reaction rate (mass/volume/time)
**-----------------------------------------------------------
*/
{
double c1;
quality_t *qual = &pr->quality;
// Find bulk reaction potential taking into account
// limiting potential & reaction order.
// Zero-order kinetics:
if (order == 0.0) c = 1.0;
// Michaelis-Menton kinetics:
else if (order < 0.0)
{
c1 = qual->Climit + SGN(kb) * c;
if (fabs(c1) < TINY) c1 = SGN(c1) * TINY;
c = c / c1;
}
// N-th order kinetics:
else
{
// Account for limiting potential
if (qual->Climit == 0.0) c1 = c;
else c1 = MAX(0.0, SGN(kb) * (qual->Climit - c));
// Compute concentration potential
if (order == 1.0) c = c1;
else if (order == 2.0) c = c1 * c;
else c = c1 * pow(MAX(0.0, c), order - 1.0);
}
// Reaction rate = bulk coeff. * potential
if (c < 0) c = 0;
return kb * c;
}
double wallrate(EN_Project *pr, double c, double d, double kw, double kf)
/*
**------------------------------------------------------------
** Input: c = current WQ concentration
** d = pipe diameter
** kw = intrinsic wall reaction coeff.
** kf = mass transfer coeff. for 0-order reaction
** (ft/sec) or apparent wall reaction coeff.
** for 1-st order reaction (1/sec)
** Output: returns wall reaction rate in mass/ft3/sec
** Purpose: computes wall reaction rate
**------------------------------------------------------------
*/
{
quality_t *qual = &pr->quality;
if (kw == 0.0 || d == 0.0) return (0.0);
if (qual->WallOrder == 0.0) // 0-order reaction */
{
kf = SGN(kw) * c * kf; //* Mass transfer rate (mass/ft2/sec)
kw = kw * SQR(pr->Ucf[ELEV]); // Reaction rate (mass/ft2/sec)
if (fabs(kf) < fabs(kw)) kw = kf; // Reaction mass transfer limited
return (kw * 4.0 / d); // Reaction rate (mass/ft3/sec)
}
else return (c * kf); // 1st-order reaction
}
double mixtank(EN_Project *pr, int n, double volin, double massin, double volout)
/*
**------------------------------------------------------------
** Input: n = node index
** volin = inflow volume to tank over time step
** massin = mass inflow to tank over time step
** volout = outflow volume from tank over time step
** Output: returns new quality for tank
** Purpose: mixes inflow with tank's contents to update its quality.
**------------------------------------------------------------
*/
{
int i;
double vnet;
EN_Network *net = &pr->network;
quality_t *qual = &pr->quality;
i = n - net->Njuncs;
vnet = volin - volout;
switch (net->Tank[i].MixModel)
{
case MIX1: tankmix1(pr, i, volin, massin, vnet); break;
case MIX2: tankmix2(pr, i, volin, massin, vnet); break;
case FIFO: tankmix3(pr, i, volin, massin, vnet); break;
case LIFO: tankmix4(pr, i, volin, massin, vnet); break;
}
return net->Tank[i].C;
}
void tankmix1(EN_Project *pr, int i, double vin, double win, double vnet)
/*
**---------------------------------------------
** Input: i = tank index
** vin = inflow volume
** win = mass inflow
** vnet = inflow - outflow
** Output: none
** Purpose: updates quality in a complete mix tank model
**---------------------------------------------
*/
{
int k;
double vnew;
Pseg seg;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
Stank *tank = &net->Tank[i];
k = net->Nlinks + i;
seg = qual->FirstSeg[k];
if (seg)
{
vnew = seg->v + vin;
if (vnew > 0.0) seg->c = (seg->c * seg->v + win) / vnew;
seg->v += vnet;
seg->v = MAX(0.0, seg->v);
tank->C = seg->c;
}
}
void tankmix2(EN_Project *pr, int i, double vin, double win, double vnet)
/*
**------------------------------------------------
** Input: i = tank index
** vin = inflow volume
** win = mass inflow
** vnet = inflow - outflow
** Output: none
** Purpose: updates quality in a 2-compartment tank model
**------------------------------------------------
*/
{
int k;
double vt, // Transferred volume
vmz; // Full mixing zone volume
Pseg mixzone, // Mixing zone segment
stagzone; // Stagnant zone segment
EN_Network *net = &pr->network;
quality_t *qual = &pr->quality;
Stank *tank = &pr->network.Tank[i];
// Identify segments for each compartment
k = net->Nlinks + i;
mixzone = qual->LastSeg[k];
stagzone = qual->FirstSeg[k];
if (mixzone == NULL || stagzone == NULL) return;
// Full mixing zone volume
vmz = tank->V1max;
// Tank is filling
vt = 0.0;
if (vnet > 0.0)
{
vt = MAX(0.0, (mixzone->v + vnet - vmz));
if (vin > 0.0)
{
mixzone->c = ((stagzone->c) * (stagzone->v) + win) /
(mixzone->v + vin);
}
if (vt > 0.0)
{
stagzone->c = ((stagzone->c) * (stagzone->v) +
(mixzone->c) * vt) / (stagzone->v + vt);
}
}
// Tank is emptying
else if (vnet < 0.0)
{
if (stagzone->v > 0.0) vt = MIN(stagzone->v, (-vnet));
if (vin + vt > 0.0)
{
mixzone->c = ((mixzone->c) * (mixzone->v) + win +
(stagzone->c) * vt) / (mixzone->v + vin + vt);
}
}
// Update segment volumes
if (vt > 0.0)
{
mixzone->v = vmz;
if (vnet > 0.0) stagzone->v += vt;
else stagzone->v = MAX(0.0, ((stagzone->v) - vt));
}
else
{
mixzone->v += vnet;
mixzone->v = MIN(mixzone->v, vmz);
mixzone->v = MAX(0.0, mixzone->v);
stagzone->v = 0.0;
}
// Use quality of mixing zone to represent quality of
// tank since this is where outflow begins to flow from
tank->C = mixzone->c;
}
void tankmix3(EN_Project *pr, int i, double vin, double win, double vnet)
/*
**----------------------------------------------------------
** Input: i = tank index
** vin = inflow volume
** win = mass inflow
** vnet = inflow - outflow
** Output: none
** Purpose: Updates quality in a First-In-First-Out (FIFO) tank model.
**----------------------------------------------------------
*/
{
int k;
double vout, vseg;
double cin, vsum, wsum;
Pseg seg;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
Stank *tank = &pr->network.Tank[i];
k = net->Nlinks + i;
if (qual->LastSeg[k] == NULL || qual->FirstSeg[k] == NULL) return;
// Add new last segment for flow entering the tank
if (vin > 0.0)
{
// ... increase segment volume if inflow has same quality as segment
cin = win / vin;
seg = qual->LastSeg[k];
if (fabs(seg->c - cin) < qual->Ctol) seg->v += vin;
// ... otherwise add a new last segment to the tank
else addseg(pr, k, vin, cin);
}
// Withdraw flow from first segment
vsum = 0.0;
wsum = 0.0;
vout = vin - vnet;
while (vout > 0.0)
{
seg = qual->FirstSeg[k];
if (seg == NULL) break;
vseg = seg->v; // Flow volume from leading seg
vseg = MIN(vseg, vout);
if (seg == qual->LastSeg[k]) vseg = vout;
vsum += vseg;
wsum += (seg->c) * vseg;
vout -= vseg; // Remaining flow volume
if (vout >= 0.0 && vseg >= seg->v) // Seg used up
{
if (seg->prev)
{
qual->FirstSeg[k] = seg->prev;
seg->prev = qual->FreeSeg;
qual->FreeSeg = seg;
}
}
else seg->v -= vseg; // Remaining volume in segment
}
// Use quality withdrawn from 1st segment
// to represent overall quality of tank
if (vsum > 0.0) tank->C = wsum / vsum;
else if (qual->FirstSeg[k] == NULL) tank->C = 0.0;
else tank->C = qual->FirstSeg[k]->c;
}
void tankmix4(EN_Project *pr, int i, double vin, double win, double vnet)
/*
**----------------------------------------------------------
** Input: i = tank index
** vin = inflow volume
** win = mass inflow
** vnet = inflow - outflow
** Output: none
** Purpose: Updates quality in a Last In-First Out (LIFO) tank model.
**----------------------------------------------------------
*/
{
int k, n;
double cin, vsum, wsum, vseg;
Pseg seg;
EN_Network *net = &pr->network;
quality_t *qual = &pr->quality;
Stank *tank = &pr->network.Tank[i];
k = net->Nlinks + i;
if (qual->LastSeg[k] == NULL || qual->FirstSeg[k] == NULL) return;
// Find inflows & outflows
n = tank->Node;
if (vin > 0.0) cin = win / vin;
else cin = 0.0;
// If tank filling, then create new last seg
tank->C = qual->LastSeg[k]->c;
seg = qual->LastSeg[k];
if (vnet > 0.0)
{
// ... quality is the same, so just add flow volume to last seg
if (fabs(seg->c - cin) < qual->Ctol) seg->v += vnet;
// ... otherwise add a new last seg to tank which points to old last seg
else
{
qual->LastSeg[k] = NULL;
addseg(pr, k, vnet, cin);
qual->LastSeg[k]->prev = seg;
}
// ... update reported tank quality
tank->C = qual->LastSeg[k]->c;
}
// If tank emptying then remove last segments until vnet consumed
else if (vnet < 0.0)
{
vsum = 0.0;
wsum = 0.0;
vnet = -vnet;
while (vnet > 0.0)
{
seg = qual->LastSeg[k];
if (seg == NULL) break;
vseg = seg->v;
vseg = MIN(vseg, vnet);
if (seg == qual->FirstSeg[k]) vseg = vnet;
vsum += vseg;
wsum += (seg->c) * vseg;
vnet -= vseg;
if (vnet >= 0.0 && vseg >= seg->v) // Seg used up
{
if (seg->prev)
{
qual->LastSeg[k] = seg->prev;
seg->prev = qual->FreeSeg;
qual->FreeSeg = seg;
}
}
else seg->v -= vseg; // Remaining volume in segment
}
// Reported tank quality is mixture of flow released and any inflow
tank->C = (wsum + win) / (vsum + vin);
}
}

776
src/qualroute.c Normal file
View File

@@ -0,0 +1,776 @@
/*
*********************************************************************
QUALROUTE.C -- water quality routing module for the EPANET program
*********************************************************************
*/
#include <stdio.h>
#include <math.h>
#include "mempool.h"
#include "types.h"
// Macro to compute the volume of a link
#define LINKVOL(k) (0.785398 * net->Link[(k)].Len * SQR(net->Link[(k)].Diam))
// Macro to get link flow compatible with flow saved to hydraulics file
#define LINKFLOW(k) ((hyd->LinkStatus[k] <= CLOSED) ? 0.0 : hyd->LinkFlows[k])
// Exported Functions
int buildilists(EN_Project *pr);
int sortnodes(EN_Project *pr);
void transport(EN_Project *pr, long);
void initsegs(EN_Project *pr);
void reversesegs(EN_Project *pr, int);
void addseg(EN_Project *pr, int, double, double);
// Imported Functions
extern double findsourcequal(EN_Project *pr, int, double, double, long);
extern void reactpipes(EN_Project *pr, long);
extern void reacttanks(EN_Project *pr, long);
extern double mixtank(EN_Project *pr, int, double, double, double);
// Local Functions
static void evalnodeinflow(EN_Project *pr, int, long, double *, double *);
static void evalnodeoutflow(EN_Project *pr, int, double, long);
static double findnodequal(EN_Project *pr, int, double, double, double, long);
static double noflowqual(EN_Project *pr, int);
static void updatemassbalance(EN_Project *pr, int, double, double, long);
static int selectnonstacknode(EN_Project *pr, int, int *);
void transport(EN_Project *pr, long tstep)
/*
**--------------------------------------------------------------
** Input: tstep = length of current time step
** Output: none
** Purpose: transports constituent mass through the pipe network
** under a period of constant hydraulic conditions.
**--------------------------------------------------------------
*/
{
int i, j, k, m, n;
double volin, massin, volout, nodequal;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
EN_Network *net = &pr->network;
// React contents of each pipe and tank
if (qual->Reactflag)
{
reactpipes(pr, tstep);
reacttanks(pr, tstep);
}
// Analyze each node in topological order
for (j = 1; j <= net->Nnodes; j++)
{
// ... index of node to be processed
n = qual->SortedNodes[j];
// ... zero out mass & flow volumes for this node
volin = 0.0;
massin = 0.0;
volout = 0.0;
// ... examine each link with flow into the node
for (i = qual->IlistPtr[n]; i < qual->IlistPtr[n + 1]; i++)
{
// ... k is index of next link incident on node n
k = qual->Ilist[i];
// ... link has flow into node - add it to node's inflow
// (m is index of link's downstream node)
m = net->Link[k].N2;
if (qual->FlowDir[k] < 0) m = net->Link[k].N1;
if (m == n)
{
evalnodeinflow(pr, k, tstep, &volin, &massin);
}
// ... link has flow out of node - add it to node's outflow
else volout += fabs(LINKFLOW(k));
}
// ... if node is a junction, add on any external outflow (e.g., demands)
if (net->Node[n].Type == EN_JUNCTION)
{
volout += MAX(0.0, hyd->NodeDemand[n]);
}
// ... convert from outflow rate to volume
volout *= tstep;
// ... find the concentration of flow leaving the node
nodequal = findnodequal(pr, n, volin, massin, volout, tstep);
// ... examine each link with flow out of the node
for (i = qual->IlistPtr[n]; i < qual->IlistPtr[n + 1]; i++)
{
// ... link k incident on node n has upstream node m equal to n
k = qual->Ilist[i];
m = net->Link[k].N1;
if (qual->FlowDir[k] < 0) m = net->Link[k].N2;
if (m == n)
{
// ... send flow at new node concen. into link
evalnodeoutflow(pr, k, nodequal, tstep);
}
}
updatemassbalance(pr, n, massin, volout, tstep);
}
}
void evalnodeinflow(EN_Project *pr, int k, long tstep, double *volin,
double *massin)
/*
**--------------------------------------------------------------
** Input: k = link index
** tstep = quality routing time step
** Output: volin = flow volume entering a node
** massin = constituent mass entering a node
** Purpose: adds the contribution of a link's outflow volume
** and constituent mass to the total inflow into its
** downstream node over a time step.
**--------------------------------------------------------------
*/
{
double q, v, vseg;
Pseg seg;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
// Get flow rate (q) and flow volume (v) through link
q = LINKFLOW(k);
v = fabs(q) * tstep;
// Transport flow volume v from link's leading segments into downstream
// node, removing segments once their full volume is consumed
while (v > 0.0)
{
seg = qual->FirstSeg[k];
if (!seg) break;
// ... volume transported from first segment is smaller of
// remaining flow volume & segment volume
vseg = seg->v;
vseg = MIN(vseg, v);
// ... update total volume & mass entering downstream node
*volin += vseg;
*massin += vseg * seg->c;
// ... reduce remaining flow volume by amount transported
v -= vseg;
// ... if all of segment's volume was transferred
if (v >= 0.0 && vseg >= seg->v)
{
// ... replace this leading segment with the one behind it
qual->FirstSeg[k] = seg->prev;
if (qual->FirstSeg[k] == NULL) qual->LastSeg[k] = NULL;
// ... recycle the used up segment
seg->prev = qual->FreeSeg;
qual->FreeSeg = seg;
}
// ... otherwise just reduce this segment's volume
else seg->v -= vseg;
}
}
double findnodequal(EN_Project *pr, int n, double volin,
double massin, double volout, long tstep)
/*
**--------------------------------------------------------------
** Input: n = node index
** volin = flow volume entering node
** massin = mass entering node
** volout = flow volume leaving node
** tstep = length of current time step
** Output: returns water quality in a node's outflow
** Purpose: computes a node's new quality from its inflow
** volume and mass, including any source contribution.
**--------------------------------------------------------------
*/
{
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
// Node is a junction - update its water quality
if (net->Node[n].Type == EN_JUNCTION)
{
// ... dilute inflow with any external negative demand
volin -= MIN(0.0, hyd->NodeDemand[n]) * tstep;
// ... new concen. is mass inflow / volume inflow
if (volin > 0.0) qual->NodeQual[n] = massin / volin;
// ... if no inflow adjust quality for reaction in connecting pipes
else if (qual->Reactflag) qual->NodeQual[n] = noflowqual(pr, n);
}
// Node is a tank - use its mixing model to update its quality
else if (net->Node[n].Type == EN_TANK)
{
qual->NodeQual[n] = mixtank(pr, n, volin, massin, volout);
}
// Add any external quality source onto node's concen.
qual->SourceQual = 0.0;
// For source tracing analysis find tracer added at source node
if (qual->Qualflag == TRACE)
{
if (n == qual->TraceNode)
{
// ... quality added to network is difference between tracer
// concentration (100 mg/L) and current node quality
if (net->Node[n].Type == EN_RESERVOIR) qual->SourceQual = 100.0;
else qual->SourceQual = MAX(100.0 - qual->NodeQual[n], 0.0);
qual->NodeQual[n] = 100.0;
}
return qual->NodeQual[n];
}
// Find quality contribued by any external chemical source
else qual->SourceQual = findsourcequal(pr, n, volin, volout, tstep);
if (qual->SourceQual == 0.0) return qual->NodeQual[n];
// Combine source quality with node quality
switch (net->Node[n].Type)
{
case EN_JUNCTION:
qual->NodeQual[n] += qual->SourceQual;
return qual->NodeQual[n];
case EN_TANK:
return qual->NodeQual[n] + qual->SourceQual;
case EN_RESERVOIR:
qual->NodeQual[n] = qual->SourceQual;
return qual->SourceQual;
}
return qual->NodeQual[n];
}
double noflowqual(EN_Project *pr, int n)
/*
**--------------------------------------------------------------
** Input: n = node index
** Output: quality for node n
** Purpose: sets the quality for a junction node that has no
** inflow to the average of the quality in its
** adjoining link segments.
** Note: this function is only used for reactive substances.
**--------------------------------------------------------------
*/
{
int i, k, inflow, kount = 0;
double c = 0.0;
FlowDirection dir;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
// Examine each link incident on the node
for (i = qual->IlistPtr[n]; i < qual->IlistPtr[n + 1]; i++)
{
// ... index of an incident link
k = qual->Ilist[i];
dir = qual->FlowDir[k];
// Node n is link's downstream node - add quality
// of link's first segment to average
if (net->Link[k].N2 == n && dir >= 0) inflow = TRUE;
else if (net->Link[k].N1 == n && dir < 0) inflow = TRUE;
else inflow = FALSE;
if (inflow == TRUE && qual->FirstSeg[k] != NULL)
{
c += qual->FirstSeg[k]->c;
kount++;
}
// Node n is link's upstream node - add quality
// of link's last segment to average
else if (inflow == FALSE && qual->LastSeg[k] != NULL)
{
c += qual->LastSeg[k]->c;
kount++;
}
}
if (kount > 0) c = c / (double)kount;
return c;
}
void evalnodeoutflow(EN_Project *pr, int k, double c, long tstep)
/*
**--------------------------------------------------------------
** Input: k = link index
** c = quality from upstream node
** tstep = time step
** Output: none
** Purpose: releases flow volume and mass from the upstream
** node of a link over a time step.
**--------------------------------------------------------------
*/
{
double v;
Pseg seg;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
// Find flow volume (v) released over time step
v = fabs(LINKFLOW(k)) * tstep;
if (v == 0.0) return;
// Release flow and mass into upstream end of the link
// ... case where link has a last (most upstream) segment
seg = qual->LastSeg[k];
if (seg)
{
// ... if node quality close to segment quality then mix
// the nodal outflow volume with the segment's volume
if (fabs(seg->c - c) < qual->Ctol)
{
seg->c = (seg->c*seg->v + c*v) / (seg->v + v);
seg->v += v;
}
// ... otherwise add a new segment at upstream end of link
else addseg(pr, k, v, c);
}
// ... link has no segments so add one
else addseg(pr, k, v, c);
}
void updatemassbalance(EN_Project *pr, int n, double massin,
double volout, long tstep)
/*
**--------------------------------------------------------------
** Input: n = node index
** massin = mass inflow to node
** volout = outflow volume from node
** Output: none
** Purpose: Adds a node's external mass inflow and outflow
** over the current time step to the network's
** overall mass balance.
**--------------------------------------------------------------
*/
{
double masslost = 0.0,
massadded = 0.0;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
switch (net->Node[n].Type)
{
// Junctions lose mass from outflow demand & gain it from source inflow
case EN_JUNCTION:
masslost = MAX(0.0, hyd->NodeDemand[n]) * tstep * qual->NodeQual[n];
massadded = qual->SourceQual * volout;
break;
// Reservoirs add mass from quality source if specified or from a fixed
// initial quality
case EN_RESERVOIR:
masslost = massin;
if (qual->SourceQual > 0.0) massadded = qual->SourceQual * volout;
else massadded = qual->NodeQual[n] * volout;
break;
// Tanks add mass only from external source inflow
case EN_TANK:
massadded = qual->SourceQual * volout;
break;
}
qual->massbalance.outflow += masslost;
qual->massbalance.inflow += massadded;
}
int buildilists(EN_Project *pr)
/*
**--------------------------------------------------------------
** Input: none
** Output: returns an error code
** Purpose: Creates link incidence lists for each node stored
** in compact form.
**--------------------------------------------------------------
*/
{
int i, j, k, n, n1, n2;
int *degree = NULL;
quality_t *qual = &pr->quality;
EN_Network *net = &pr->network;
// Allocate an array to count # links incident on each node
n = net->Nnodes + 1;
degree = (int *)calloc(n, sizeof(int));
if (degree == NULL) return 101;
// Count # links incident on each node
for (k = 1; k <= net->Nlinks; k++)
{
degree[net->Link[k].N1]++;
degree[net->Link[k].N2]++;
}
// Use incidence counts to determine start position of
// each node's incidence list in Xilist
qual->IlistPtr[1] = 1;
for (i = 1; i <= n; i++)
{
qual->IlistPtr[i + 1] = qual->IlistPtr[i] + degree[i];
}
// Add each link to the incidence lists of its start & end nodes
for (i = 1; i <= net->Nnodes; i++) degree[i] = 0;
for (k = 1; k <= net->Nlinks; k++)
{
// j is index of next unused location in link's start node list
n1 = net->Link[k].N1;
j = qual->IlistPtr[n1] + degree[n1];
qual->Ilist[j] = k;
degree[n1]++;
// Repeat same for end node
n2 = net->Link[k].N2;
j = qual->IlistPtr[n2] + degree[n2];
qual->Ilist[j] = k;
degree[n2]++;
}
free(degree);
/*//////// QA CHECK
for (i = 1; i <= net->Nnodes; i++)
{
printf("\nNode %s: ", net->Node[i].ID);
for (j = qual->IlistPtr[i]; j < qual->IlistPtr[i + 1]; j++)
{
printf(" %s,", net->Link[qual->Ilist[j]].ID);
}
}
*/
return 0;
}
int sortnodes(EN_Project *pr)
/*
**--------------------------------------------------------------
** Input: none
** Output: returns an error code
** Purpose: topologically sorts nodes from upstream to downstream.
** Note: links with negligible flow are ignored since they can
** create spurious cycles that cause the sort to fail.
**--------------------------------------------------------------
*/
{
int i, j, k, n;
int *indegree = NULL;;
int *stack = NULL;
int stacksize = 0;
int numsorted = 0;
int errcode = 0;
FlowDirection dir;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
// Allocate an array to count # links with inflow to each node
// and for a stack to hold nodes waiting to be processed
indegree = (int *)calloc(net->Nnodes + 1, sizeof(int));
stack = (int *)calloc(net->Nnodes + 1, sizeof(int));
if (indegree && stack)
{
// Count links with "non-negligible" inflow to each node
for (k = 1; k <= net->Nlinks; k++)
{
dir = qual->FlowDir[k];
if (dir == POSITIVE) n = net->Link[k].N2;
else if (dir == NEGATIVE) n = net->Link[k].N1;
else continue;
indegree[n]++;
}
// Place nodes with no inflow onto a stack
for (i = 1; i <= net->Nnodes; i++)
{
if (indegree[i] == 0)
{
stacksize++;
stack[stacksize] = i;
}
}
// Examine each node on the stack until none are left
while (numsorted < net->Nnodes)
{
// ... if stack is empty then a cycle exists
if (stacksize == 0)
{
// ... add a non-sorted node connected to a sorted one to stack
j = selectnonstacknode(pr, numsorted, indegree);
if (j == 0) break; // This shouldn't happen.
indegree[j] = 0;
stacksize++;
stack[stacksize] = j;
}
// ... make the last node added to the stack the next
// in sorted order & remove it from the stack
i = stack[stacksize];
stacksize--;
numsorted++;
qual->SortedNodes[numsorted] = i;
// ... for each outflow link from this node reduce the in-degree
// of its downstream node
for (j = qual->IlistPtr[i]; j < qual->IlistPtr[i + 1]; j++)
{
// ... k is the index of the next link incident on node i
k = qual->Ilist[j];
// ... skip link if flow is negligible
if (qual->FlowDir[k] == 0) continue;
// ... link has flow out of node (downstream node n not equal to i)
n = net->Link[k].N2;
if (qual->FlowDir[k] < 0) n = net->Link[k].N1;
// ... reduce degree of node n
if (n != i && indegree[n] > 0)
{
indegree[n]--;
// ... no more degree left so add node n to stack
if (indegree[n] == 0)
{
stacksize++;
stack[stacksize] = n;
}
}
}
}
}
else errcode = 101;
if (numsorted < net->Nnodes) errcode = 120;
FREE(indegree);
FREE(stack);
/*
/////////////////// QA CHECK
snprintf(pr->Msg, MAXMSG, "\n\nSorted Nodes:");
writeline(pr, pr->Msg);
for (i = 1; i <= net->Nnodes; i++)
{
j = qual->SortedNodes[i];
snprintf(pr->Msg, MAXMSG, "%s", net->Node[j].ID);
writeline(pr, pr->Msg);
}
//printf("\n");
//system("pause");
/////////////////
*/
return errcode;
}
int selectnonstacknode(EN_Project *pr, int numsorted, int *indegree)
/*
**--------------------------------------------------------------
** Input: numsorted = number of nodes that have been sorted
** indegree = number of inflow links to each node
** Output: returns a node index
** Purpose: selects a next node for sorting when a cycle exists.
**--------------------------------------------------------------
*/
{
int i, j, k, m, n;
quality_t *qual = &pr->quality;
EN_Network *net = &pr->network;
// Examine each sorted node in last in - first out order
for (i = numsorted; i > 0; i--)
{
// For each link connected to the sorted node
m = qual->SortedNodes[i];
for (j = qual->IlistPtr[m]; j < qual->IlistPtr[m + 1]; j++)
{
// ... k is index of next link incident on node m
k = qual->Ilist[j];
// ... n is the node of link k opposite to node m
n = net->Link[k].N2;
if (n == m) n = net->Link[k].N1;
// ... select node n if it still has inflow links
if (indegree[n] > 0) return n;
}
}
// If no node was selected by the above process then return the
// first node that still has inflow links remaining
for (i = 1; i <= net->Nnodes; i++)
{
if (indegree[i] > 0) return i;
}
// If all else fails return 0 indicating that no node was selected
return 0;
}
void initsegs(EN_Project *pr)
/*
**--------------------------------------------------------------
** Input: none
** Output: none
** Purpose: initializes water quality volume segments in each
** pipe and tank.
**--------------------------------------------------------------
*/
{
int j, k;
double c, v, v1;
EN_Network *net = &pr->network;
hydraulics_t *hyd = &pr->hydraulics;
quality_t *qual = &pr->quality;
// Add one segment with assigned downstream node quality to each pipe
for (k = 1; k <= net->Nlinks; k++)
{
qual->FirstSeg[k] = NULL;
qual->LastSeg[k] = NULL;
if (net->Link[k].Type == EN_PIPE)
{
v = LINKVOL(k);
j = net->Link[k].N2;
c = qual->NodeQual[j];
addseg(pr, k, v, c);
}
}
// Initialize segments in tanks
for (j = 1; j <= net->Ntanks; j++)
{
// Skip reservoirs
if (net->Tank[j].A == 0.0) continue;
// Establish initial tank quality & volume
k = net->Tank[j].Node;
c = net->Node[k].C0;
v = net->Tank[j].V0;
// Create one volume segment for entire tank
k = net->Nlinks + j;
qual->FirstSeg[k] = NULL;
qual->LastSeg[k] = NULL;
addseg(pr, k, v, c);
// Create a 2nd segment for the 2-compartment tank model
if (net->Tank[j].MixModel == MIX2)
{
// ... mixing zone segment
v1 = MAX(0, v - net->Tank[j].V1max);
qual->FirstSeg[k]->v = v1;
// ... stagnant zone segment
v = v - v1;
addseg(pr, k, v, c);
}
}
}
void reversesegs(EN_Project *pr, int k)
/*
**--------------------------------------------------------------
** Input: k = link index
** Output: none
** Purpose: re-orients a link's segments when flow reverses.
**--------------------------------------------------------------
*/
{
Pseg seg, nseg, pseg;
quality_t *qual = &pr->quality;
seg = qual->FirstSeg[k];
qual->FirstSeg[k] = qual->LastSeg[k];
qual->LastSeg[k] = seg;
pseg = NULL;
while (seg != NULL)
{
nseg = seg->prev;
seg->prev = pseg;
pseg = seg;
seg = nseg;
}
}
void addseg(EN_Project *pr, int k, double v, double c)
/*
**-------------------------------------------------------------
** Input: k = segment chain index
** v = segment volume
** c = segment quality
** Output: none
** Purpose: adds a segment to the start of a link
** upstream of its current last segment.
**-------------------------------------------------------------
*/
{
Pseg seg;
quality_t *qual = &pr->quality;
// Grab the next free segment from the segment pool if available
if (qual->FreeSeg != NULL)
{
seg = qual->FreeSeg;
qual->FreeSeg = seg->prev;
}
// Otherwise allocate a new segment
else
{
seg = (struct Sseg *) mempool_alloc(qual->SegPool, sizeof(struct Sseg));
if (seg == NULL)
{
qual->OutOfMemory = TRUE;
return;
}
}
// Assign volume and quality to the segment
seg->v = v;
seg->c = c;
// Add the new segment to the end of the segment chain
seg->prev = NULL;
if (qual->FirstSeg[k] == NULL) qual->FirstSeg[k] = seg;
if (qual->LastSeg[k] != NULL) qual->LastSeg[k]->prev = seg;
qual->LastSeg[k] = seg;
}

46
src/report.c
View File

@@ -376,6 +376,52 @@ void writehydstat(EN_Project *pr, int iter, double relerr)
writeline(pr, " ");
} /* End of writehydstat */
void writemassbalance(EN_Project *pr)
/*
**-------------------------------------------------------------
** Input: none
** Output: none
** Purpose: writes water quality mass balance ratio
** (Outflow + Final Storage) / Inflow + Initial Storage)
** to report file.
**-------------------------------------------------------------
*/
{
char s1[MAXMSG+1];
char *units[] = {"", " (mg)", " (ug)", " (hrs)"};
int kunits = 0;
quality_t *qual = &pr->quality;
if (qual->Qualflag == TRACE) kunits = 1;
else if (qual->Qualflag == AGE) kunits = 3;
else
{
if (match(qual->ChemUnits, "mg")) kunits = 1;
else if (match(qual->ChemUnits, "ug")) kunits = 2;
}
snprintf(s1, MAXMSG, "Water Quality Mass Balance%s", units[kunits]);
writeline(pr, s1);
snprintf(s1, MAXMSG, "================================");
writeline(pr, s1);
snprintf(s1, MAXMSG, "Initial Mass: %12.5e", qual->massbalance.initial);
writeline(pr, s1);
snprintf(s1, MAXMSG, "Mass Inflow: %12.5e", qual->massbalance.inflow);
writeline(pr, s1);
snprintf(s1, MAXMSG, "Mass Outflow: %12.5e", qual->massbalance.outflow);
writeline(pr, s1);
snprintf(s1, MAXMSG, "Mass Reacted: %12.5e", qual->massbalance.reacted);
writeline(pr, s1);
snprintf(s1, MAXMSG, "Final Mass: %12.5e", qual->massbalance.final);
writeline(pr, s1);
snprintf(s1, MAXMSG, "Mass Ratio: %-0.5f", qual->massbalance.ratio);
writeline(pr, s1);
snprintf(s1, MAXMSG, "================================\n");
writeline(pr, s1);
}
void writeenergy(EN_Project *pr)
/*
**-------------------------------------------------------------

95
src/types.h
View File

@@ -23,6 +23,7 @@ AUTHOR: L. Rossman
#include "epanet2.h"
#include "hash.h"
#include "util/errormanager.h"
#include <stdio.h>
/*********************************************************/
@@ -286,8 +287,9 @@ typedef enum {
} HdrType;
typedef enum {
POSITIVE,
NEGATIVE
NEGATIVE = -1, // Flow in reverse of pre-assigned direction
ZERO_FLOW = 0, // Zero flow
POSITIVE = 1 // Flow in pre-assigned direction
} FlowDirection;
typedef enum {
@@ -538,61 +540,74 @@ typedef struct s_ActItem /* Action list item */
struct s_ActItem *next;
} ActItem;
typedef struct
{
double initial;
double inflow;
double outflow;
double reacted;
double final;
double ratio;
} MassBalance;
// Forward declaration of the Mempool structure defined in mempool.h
struct Mempool;
typedef struct {
char
Qualflag, /// Water quality flag
OpenQflag, /// Quality system opened flag
Reactflag; /// Reaction indicator
Qualflag, // Water quality flag
OpenQflag, // Quality system opened flag
Reactflag, // Reaction indicator
OutOfMemory; // Out of memory indicator
char
ChemName[MAXID+1], /* Name of chemical */
ChemUnits[MAXID+1]; /* Units of chemical */
ChemName[MAXID+1], // Name of chemical
ChemUnits[MAXID+1]; // Units of chemical
int
TraceNode; /// Source node for flow tracing
TraceNode, // Source node for flow tracing
*SortedNodes, // Topologically sorted node indexes
*Ilist, // Link incidence lists for all nodes
*IlistPtr; // Start index of each node in Ilist
double
Ctol, /// Water quality tolerance
Diffus, /// Diffusivity (sq ft/sec)
Wbulk, /// Avg. bulk reaction rate
Wwall, /// Avg. wall reaction rate
Wtank, /// Avg. tank reaction rate
Wsource, /// Avg. mass inflow
Rfactor, /// Roughness-reaction factor
BulkOrder, /// Bulk flow reaction order
WallOrder, /// Pipe wall reaction order
TankOrder, /// Tank reaction order
Kbulk, /// Global bulk reaction coeff.
Kwall, /// Global wall reaction coeff.
Climit, /// Limiting potential quality
*NodeQual, /// Node quality state
*TempQual, /// General purpose array for water quality
*QTempVolumes,
*QTankVolumes,
*QLinkFlow,
*PipeRateCoeff;
Ctol, // Water quality tolerance
Diffus, // Diffusivity (sq ft/sec)
Wbulk, // Avg. bulk reaction rate
Wwall, // Avg. wall reaction rate
Wtank, // Avg. tank reaction rate
Wsource, // Avg. mass inflow
Rfactor, // Roughness-reaction factor
Sc, // Schmidt Number
Bucf, // Bulk reaction units conversion factor
Tucf, // Tank reaction units conversion factor
BulkOrder, // Bulk flow reaction order
WallOrder, // Pipe wall reaction order
TankOrder, // Tank reaction order
Kbulk, // Global bulk reaction coeff.
Kwall, // Global wall reaction coeff.
Climit, // Limiting potential quality
SourceQual, // External source quality
*NodeQual, // Reported node quality state
*PipeRateCoeff; // Pipe reaction rate coeffs.
long
Qstep, /// Quality time step (sec)
Qtime; /// Current quality time (sec)
Qstep, // Quality time step (sec)
Qtime; // Current quality time (sec)
char OutOfMemory; /* Out of memory indicator */
struct Mempool *SegPool; // Memory pool for water quality segments
struct Mempool
*SegPool; // Memory pool for water quality segments
Pseg FreeSeg; /* Pointer to unused segment */
Pseg *FirstSeg, /* First (downstream) segment in each pipe */
*LastSeg; /* Last (upstream) segment in each pipe */
FlowDirection *FlowDir; /* Flow direction for each pipe */
double *VolIn; /* Total volume inflow to node */
double *MassIn; /* Total mass inflow to node */
double Sc; /* Schmidt Number */
double Bucf; /* Bulk reaction units conversion factor */
double Tucf; /* Tank reaction units conversion factor */
Pseg
FreeSeg, // Pointer to unused segment
*FirstSeg, // First (downstream) segment in each pipe
*LastSeg; // Last (upstream) segment in each pipe
FlowDirection
*FlowDir; // Flow direction for each pipe
MassBalance
massbalance; // Mass balance components
} quality_t;
typedef struct {

16
src/vars.h
View File

@@ -1,16 +0,0 @@
/*
************************************************************************
Global Variables for EPANET Program
************************************************************************
*/
#ifndef VARS_H
#define VARS_H
#include <stdio.h>
#include "hash.h"
// this single global variable is used only when the library is called in "legacy mode"
// with the 2.1-style API.
EXTERN void *_defaultModel;
#endif

8
win_build/WinSDK/Makefile.bat
View File

@@ -21,9 +21,9 @@ Find /i "x86" < checkOS.tmp > StringCheck.tmp
If %ERRORLEVEL% == 1 (
CALL "%SDK_PATH%bin\"SetEnv.cmd /x64 /release
rem : create EPANET2.DLL
cl -o epanet2.dll epanet.c util\errormanager.c hash.c hydraul.c hydcoeffs.c hydstatus.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /link /DLL
cl -o epanet2.dll epanet.c util\errormanager.c hash.c hydraul.c hydcoeffs.c hydstatus.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c qualroute.c qualreact.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /link /DLL
rem : create EPANET2.EXE
cl -o epanet2.exe epanet.c util\errormanager.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydstatus.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /I ..\src /link
cl -o epanet2.exe epanet.c util\errormanager.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydstatus.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c qualroute.c qualreact.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /I ..\src /link
md "%Build_PATH%"\64bit
move /y "%SRC_PATH%"\*.dll "%Build_PATH%"\64bit
move /y "%SRC_PATH%"\*.exe "%Build_PATH%"\64bit
@@ -35,9 +35,9 @@ rem : 32 bit with DEF
CALL "%SDK_PATH%bin\"SetEnv.cmd /x86 /release
echo "32 bit with epanet2.def mapping"
rem : create EPANET2.DLL
cl -o epanet2.dll epanet.c util\errormanager.c hash.c hydraul.c hydcoeffs.c hydstatus.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /link /DLL /def:..\win_build\WinSDK\epanet2.def /MAP
cl -o epanet2.dll epanet.c util\errormanager.c hash.c hydraul.c hydcoeffs.c hydstatus.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c qualroute.c qualreact.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /link /DLL /def:..\win_build\WinSDK\epanet2.def /MAP
rem : create EPANET2.EXE
cl -o epanet2.exe epanet.c util\errormanager.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydstatus.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /I ..\src /link
cl -o epanet2.exe epanet.c util\errormanager.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydstatus.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c qualroute.c qualreact.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /I ..\src /link
md "%Build_PATH%"\32bit
move /y "%SRC_PATH%"\*.dll "%Build_PATH%"\32bit
move /y "%SRC_PATH%"\*.exe "%Build_PATH%"\32bit