diff --git a/ReleaseNotes2_2.md b/ReleaseNotes2_2.md
index dd14ca3..618499b 100644
--- a/ReleaseNotes2_2.md
+++ b/ReleaseNotes2_2.md
@@ -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.
- *Pmin* is allowed to equal to *Preq*. 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 *Pmin*.
+## 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`
diff --git a/include/epanet2.bas b/include/epanet2.bas
index 7b738ab..cb62b96 100644
--- a/include/epanet2.bas
+++ b/include/epanet2.bas
@@ -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
diff --git a/include/epanet2.h b/include/epanet2.h
index 58a935e..9147912 100644
--- a/include/epanet2.h
+++ b/include/epanet2.h
@@ -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 {
diff --git a/include/epanet2.vb b/include/epanet2.vb
index d751762..0af9f70 100644
--- a/include/epanet2.vb
+++ b/include/epanet2.vb
@@ -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
diff --git a/run/main.c b/run/main.c
index ed79edf..5d811ee 100644
--- a/run/main.c
+++ b/run/main.c
@@ -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 []\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);
- }
+ 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);
+ // error //
+ printf(FMT11);
return(errcode);
}
else {
- // warning //
- writeConsole(FMT10);
+ // warning //
+ 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);
-}
-
-
diff --git a/src/epanet.c b/src/epanet.c
index 135e000..d021313 100644
--- a/src/epanet.c
+++ b/src/epanet.c
@@ -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;
diff --git a/src/funcs.h b/src/funcs.h
index 10fe575..8d267f9 100755
--- a/src/funcs.h
+++ b/src/funcs.h
@@ -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 */
diff --git a/src/hydraul.c b/src/hydraul.c
index 3984ec0..4f9ce46 100755
--- a/src/hydraul.c
+++ b/src/hydraul.c
@@ -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);
diff --git a/src/input1.c b/src/input1.c
index 46da66a..3eb5585 100644
--- a/src/input1.c
+++ b/src/input1.c
@@ -167,7 +167,8 @@ void setdefaults(EN_Project *pr)
hyd->RQtol = RQTOL; /* Default hydraulics parameters */
hyd->CheckFreq = CHECKFREQ;
hyd->MaxCheck = MAXCHECK;
- hyd->DampLimit = DAMPLIMIT;
+ hyd->DampLimit = DAMPLIMIT;
+ qu->massbalance.ratio = 0.0;
} /* End of setdefaults */
void initreport(report_options_t *r)
diff --git a/src/output.c b/src/output.c
index d9ec660..b91ab83 100644
--- a/src/output.c
+++ b/src/output.c
@@ -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 */
diff --git a/src/quality.c b/src/quality.c
index 1a70beb..b58acbb 100644
--- a/src/quality.c
+++ b/src/quality.c
@@ -1,50 +1,13 @@
/*
-*******************************************************************************
+*********************************************************************
-QUALITY.C -- Water Quality Simulator for EPANET Program
+QUALITY.C -- water quality engine module for the EPANET program
-VERSION: 2.00
-DATE: 5/29/00
- 9/7/00
- 10/25/00
- 8/15/07 (2.00.11)
- 2/14/08 (2.00.12)
-AUTHOR: L. Rossman
- US EPA - NRMRL
+This module works together with QUALROUTE.C and QUALREACT.C to
+compute transport and reaction of a water quality constituent within
+a pipe network over a series of time steps.
- This module contains the network water quality simulator.
-
- For each time period, hydraulic results are read in from the
- binary file HydFile, hydraulic and water quality results are
- written to the binary output file OutFile (if the current period
- is a reporting period), and the water quality is transported
- and reacted over the duration of the time period.
-
- The entry points for this module are:
- openqual() -- called from ENopenQ() in EPANET.C
- initqual() -- called from ENinitQ() in EPANET.C
- runqual() -- called from ENrunQ() in EPANET.C
- nextqual() -- called from ENnextQ() in EPANET.C
- stepqual() -- called from ENstepQ() in EPANET.C
- closequal() -- called from ENcloseQ() in EPANET.C
-
- Calls are made to:
- mempool_create()
- mempool_alloc()
- mempool_delete()
- in MEMPOOL.C to utilize a memory pool to prevent excessive malloc'ing
- when constantly creating and destroying pipe sub-segments during
- the water quality transport calculations.
-
- Calls are also made to:
- readhyd()
- readhydstep()
- savenetdata()
- saveoutput()
- savefinaloutput()
- in OUTPUT.C to retrieve hydraulic results and save all results.
-
-*******************************************************************************
+*********************************************************************
*/
#include
@@ -54,1797 +17,684 @@ AUTHOR: L. Rossman
#else
#include
#endif
+#include
-#include "hash.h"
-#include "text.h"
+#include "mempool.h"
#include "types.h"
#include "funcs.h"
-#include
-#include "mempool.h"
-/*
-** Macros to identify upstream & downstream nodes of a link
-** under the current flow and to compute link volume
-*/
-#define UP_NODE(x) \
- ((qu->FlowDir[(x)] == POSITIVE) ? net->Link[(x)].N1 : net->Link[(x)].N2)
-#define DOWN_NODE(x) \
- ((qu->FlowDir[(x)] == POSITIVE) ? net->Link[(x)].N2 : net->Link[(x)].N1)
-#define LINKVOL(k) (0.785398 * net->Link[(k)].Len * SQR(net->Link[(k)].Diam))
+// Stagnant flow tolerance
+const double QZERO = 0.005 / GPMperCFS; // 0.005 gpm = 1.114e-5 cfs
+
+// Exported Functions (declared in FUNCS.H)
+//int openqual(EN_Project *pr);
+//void initqual(EN_Project *pr);
+//int runqual(EN_Project *pr, long *);
+//int nextqual(EN_Project *pr, long *);
+//int stepqual(EN_Project *pr, long *);
+//int closequal(EN_Project *pr);
+//double avgqual(EN_Project *pr, int);
+double findsourcequal(EN_Project *pr, int, double, double, long);
+
+// Imported Functions
+extern char setreactflag(EN_Project *pr);
+extern double getucf(double);
+extern void ratecoeffs(EN_Project *pr);
+extern int buildilists(EN_Project *pr);
+extern void initsegs(EN_Project *pr);
+extern void reversesegs(EN_Project *pr, int);
+extern int sortnodes(EN_Project *pr);
+extern void transport(EN_Project *pr, long);
+
+// Local Functions
+static double sourcequal(EN_Project *pr, Psource);
+static void evalmassbalance(EN_Project *pr);
+static double findstoredmass(EN_Project *pr);
+static int flowdirchanged(EN_Project *pr);
+
int openqual(EN_Project *pr)
/*
**--------------------------------------------------------------
** Input: none
** Output: returns error code
-** Purpose: opens WQ solver system
+** Purpose: opens water quality solver
**--------------------------------------------------------------
*/
{
- int errcode = 0;
- int n;
+ int errcode = 0;
+ int n;
- quality_t *qu = &pr->quality;
- EN_Network *net = &pr->network;
+ quality_t *qual = &pr->quality;
+ EN_Network *net = &pr->network;
- /* Allocate memory pool for WQ segments */
- qu->OutOfMemory = FALSE;
- qu->SegPool = mempool_create();
- if (qu->SegPool == NULL) {
- errcode = 101;
- }
+ // Create a memory pool for water quality segments
+ qual->OutOfMemory = FALSE;
+ qual->SegPool = mempool_create();
+ if (qual->SegPool == NULL) errcode = 101;;
- /* Allocate scratch array & reaction rate array*/
- qu->TempQual = (double *)calloc(MAX((net->Nnodes + 1), (net->Nlinks + 1)),
- sizeof(double));
- qu->PipeRateCoeff = (double *)calloc((net->Nlinks + 1), sizeof(double));
- ERRCODE(MEMCHECK(qu->TempQual));
- ERRCODE(MEMCHECK(qu->PipeRateCoeff));
+ // Allocate arrays for link flow direction & reaction rates
+ n = net->Nlinks + 1;
+ qual->FlowDir = (FlowDirection *)calloc(n, sizeof(FlowDirection));
+ qual->PipeRateCoeff = (double *)calloc(n, sizeof(double));
- /* Allocate memory for WQ solver */
- n = net->Nlinks + net->Ntanks + 1;
- qu->FirstSeg = (Pseg *)calloc(n, sizeof(Pseg));
- qu->LastSeg = (Pseg *)calloc(n, sizeof(Pseg));
- qu->FlowDir = (FlowDirection *)calloc(n, sizeof(FlowDirection));
- n = net->Nnodes + 1;
- qu->VolIn = (double *)calloc(n, sizeof(double));
- qu->MassIn = (double *)calloc(n, sizeof(double));
- ERRCODE(MEMCHECK(qu->FirstSeg));
- ERRCODE(MEMCHECK(qu->LastSeg));
- ERRCODE(MEMCHECK(qu->FlowDir));
- ERRCODE(MEMCHECK(qu->VolIn));
- ERRCODE(MEMCHECK(qu->MassIn));
- return (errcode);
+ // Allocate arrays used for volume segments in links & tanks
+ n = net->Nlinks + net->Ntanks + 1;
+ qual->FirstSeg = (Pseg *)calloc(n, sizeof(Pseg));
+ qual->LastSeg = (Pseg *)calloc(n, sizeof(Pseg));
+
+ // Allocate memory for sorted nodes and link incidence lists
+ // ... Ilist contains the list of link indexes that are incident
+ // on each node in compact format
+ n = 2 * net->Nlinks + 2;
+ qual->Ilist = (int *)calloc(n, sizeof(int));
+ // ... IlistPtr holds the position in Ilist where the
+ // incidence list for each node begins
+ n = net->Nnodes + 1;
+ qual->IlistPtr = (int *)calloc(n + 2, sizeof(int));
+ // ... SortedNodes contains the list of node indexes in topological
+ // order
+ qual->SortedNodes = (int *)calloc(n, sizeof(int));
+
+ ERRCODE(MEMCHECK(qual->FlowDir));
+ ERRCODE(MEMCHECK(qual->PipeRateCoeff));
+ ERRCODE(MEMCHECK(qual->FirstSeg));
+ ERRCODE(MEMCHECK(qual->LastSeg));
+ ERRCODE(MEMCHECK(qual->Ilist));
+ ERRCODE(MEMCHECK(qual->IlistPtr));
+ ERRCODE(MEMCHECK(qual->SortedNodes));
+
+ // Build link incidence lists
+ if (!errcode) errcode = buildilists(pr);
+ return errcode;
}
-/* Local function to compute unit conversion factor for bulk reaction rates */
-double getucf(double order) {
- if (order < 0.0) {
- order = 0.0;
- }
- if (order == 1.0) {
- return (1.0);
- } else {
- return (1. / pow(LperFT3, (order - 1.0)));
- }
-}
-void initqual(EN_Project *pr)
+int initqual(EN_Project *pr)
/*
**--------------------------------------------------------------
** Input: none
** Output: none
-** Purpose: re-initializes WQ solver system
+** Purpose: re-initializes water quality solver
**--------------------------------------------------------------
*/
{
- int i;
- quality_t *qu = &pr->quality;
- EN_Network *net = &pr->network;
- hydraulics_t *hyd = &pr->hydraulics;
- time_options_t *time = &pr->time_options;
+ int i;
+ int errcode = 0;
+ EN_Network *net = &pr->network;
+ hydraulics_t *hyd = &pr->hydraulics;
+ quality_t *qual = &pr->quality;
- /* Initialize quality, tank volumes, & source mass flows */
- for (i = 1; i <= net->Nnodes; i++) {
- qu->NodeQual[i] = net->Node[i].C0;
- }
- for (i = 1; i <= net->Ntanks; i++) {
- net->Tank[i].C = net->Node[net->Tank[i].Node].C0;
- }
- for (i = 1; i <= net->Ntanks; i++) {
- net->Tank[i].V = net->Tank[i].V0;
- }
- for (i = 1; i <= net->Nnodes; i++) {
- if (net->Node[i].S != NULL)
- net->Node[i].S->Smass = 0.0;
- }
-
- qu->QTempVolumes =
- calloc(net->Ntanks,
- sizeof(double)); // keep track of next tank volumes.
- qu->QTankVolumes =
- calloc(net->Ntanks,
- sizeof(double)); // keep track of previous step's tank volumes.
- qu->QLinkFlow = calloc(
- net->Nlinks, sizeof(double)); // keep track of previous step's link flows.
-
- /* Set WQ parameters */
- qu->Bucf = 1.0;
- qu->Tucf = 1.0;
- qu->Reactflag = 0;
- if (qu->Qualflag != NONE) {
- /* Initialize WQ at trace node (if applicable) */
- if (qu->Qualflag == TRACE) {
- qu->NodeQual[qu->TraceNode] = 100.0;
+ // Re-position hydraulics file
+ if (!hyd->OpenHflag)
+ {
+ fseek(pr->out_files.HydFile, pr->out_files.HydOffset, SEEK_SET);
}
- /* Compute Schmidt number */
- if (qu->Diffus > 0.0)
- qu->Sc = hyd->Viscos / qu->Diffus;
- else
- qu->Sc = 0.0;
+ // Set elapsed times to zero
+ qual->Qtime = 0;
+ pr->time_options.Htime = 0;
+ pr->time_options.Rtime = pr->time_options.Rstart;
+ pr->report.Nperiods = 0;
- /* Compute unit conversion factor for bulk react. coeff. */
- qu->Bucf = getucf(qu->BulkOrder);
- qu->Tucf = getucf(qu->TankOrder);
+ // Initialize node quality
+ for (i = 1; i <= net->Nnodes; i++)
+ {
+ if (qual->Qualflag == TRACE) qual->NodeQual[i] = 0.0;
+ else qual->NodeQual[i] = net->Node[i].C0;
+ if (net->Node[i].S != NULL) net->Node[i].S->Smass = 0.0;
+ }
+ if (qual->Qualflag == NONE) return errcode;
- /* Check if modeling a reactive substance */
- qu->Reactflag = setReactflag(pr);
+ // Initialize tank quality
+ for (i = 1; i <= net->Ntanks; i++)
+ {
+ net->Tank[i].C = qual->NodeQual[net->Tank[i].Node];
+ }
- /* Reset memory pool */
- qu->FreeSeg = NULL;
- mempool_reset(qu->SegPool);
- }
+ // Initialize quality at trace node (if applicable)
+ if (qual->Qualflag == TRACE) qual->NodeQual[qual->TraceNode] = 100.0;
+
+ // Compute Schmidt number
+ if (qual->Diffus > 0.0) qual->Sc = hyd->Viscos / qual->Diffus;
+ else qual->Sc = 0.0;
- /* Initialize avg. reaction rates */
- qu->Wbulk = 0.0;
- qu->Wwall = 0.0;
- qu->Wtank = 0.0;
- qu->Wsource = 0.0;
+ // Compute unit conversion factor for bulk react. coeff.
+ qual->Bucf = getucf(qual->BulkOrder);
+ qual->Tucf = getucf(qual->TankOrder);
- /* Re-position hydraulics file */
- if (!hyd->OpenHflag) {
- fseek(pr->out_files.HydFile, pr->out_files.HydOffset, SEEK_SET);
- }
+ // Check if modeling a reactive substance
+ qual->Reactflag = setreactflag(pr);
- /* Set elapsed times to zero */
- time->Htime = 0;
- qu->Qtime = 0;
- pr->time_options.Rtime = pr->time_options.Rstart;
- pr->report.Nperiods = 0;
+ // Reset memory pool used for pipe & tank segments
+ qual->FreeSeg = NULL;
+ mempool_reset(qual->SegPool);
- initsegs(pr);
+ // Create initial set of pipe & tank segments
+ initsegs(pr);
+
+ // Initialize link flow direction indicator
+ for (i = 1; i <= net->Nlinks; i++) qual->FlowDir[i] = ZERO_FLOW;
+
+ // Initialize avg. reaction rates
+ qual->Wbulk = 0.0;
+ qual->Wwall = 0.0;
+ qual->Wtank = 0.0;
+ qual->Wsource = 0.0;
+
+ // Initialize mass balance components
+ qual->massbalance.initial = findstoredmass(pr);
+ qual->massbalance.inflow = 0.0;
+ qual->massbalance.outflow = 0.0;
+ qual->massbalance.reacted = 0.0;
+ qual->massbalance.final = 0.0;
+ qual->massbalance.ratio = 0.0;
+ return errcode;
}
+
int runqual(EN_Project *pr, long *t)
/*
**--------------------------------------------------------------
** Input: none
-** Output: t = pointer to current simulation time (sec)
+** Output: t = current simulation time (sec)
** Returns: error code
** Purpose: retrieves hydraulics for next hydraulic time step
-** (at time *t) and saves current results to file
+** (at time t) and saves current results to file
**--------------------------------------------------------------
*/
{
- long hydtime; /* Hydraulic solution time */
- long hydstep; /* Hydraulic time step */
- int errcode = 0;
- int i;
+ long hydtime; // Hydraulic solution time
+ long hydstep; // Hydraulic time step
+ int errcode = 0;
- EN_Network *net = &pr->network;
- hydraulics_t *hyd = &pr->hydraulics;
- quality_t *qu = &pr->quality;
- time_options_t *time = &pr->time_options;
+ hydraulics_t *hyd = &pr->hydraulics;
+ quality_t *qual = &pr->quality;
+ time_options_t *time = &pr->time_options;
- /* Update reported simulation time */
- *t = qu->Qtime;
+ // Update reported simulation time
+ *t = qual->Qtime;
- /* Read hydraulic solution from hydraulics file */
- if (qu->Qtime == time->Htime) {
- errcode = gethyd(pr, &hydtime, &hydstep);
- if (!hyd->OpenHflag) { // test for sequential vs stepwise
- // sequential
- time->Htime = hydtime + hydstep;
- } else {
- // stepwise calculation - hydraulic results are already in memory
- for (i = 1; i <= net->Ntanks; ++i) {
- qu->QTankVolumes[i - 1] = net->Tank[i].V;
- }
-
- for (i = 1; i <= net->Nlinks; ++i) {
- if (hyd->LinkStatus[i] <= CLOSED) {
- qu->QLinkFlow[i - 1] = hyd->LinkFlows[i];
+ // Read hydraulic solution from hydraulics file
+ if (qual->Qtime == time->Htime)
+ {
+ // Read hydraulic results from file
+ if (!hyd->OpenHflag)
+ {
+ if (!readhyd(pr, &hydtime)) return 307;
+ if (!readhydstep(pr, &hydstep)) return 307;
+ time->Htime = hydtime;
}
- }
- }
- } else {
- // stepwise calculation
- for (i = 1; i <= net->Ntanks; ++i) {
- qu->QTankVolumes[i - 1] = net->Tank[i].V;
- }
- for (i = 1; i <= net->Nlinks; ++i) {
- if (hyd->LinkStatus[i] <= CLOSED) {
- qu->QLinkFlow[i - 1] = hyd->LinkFlows[i];
- }
- }
- }
+ // Save current results to output file
+ if (time->Htime >= time->Rtime)
+ {
+ if (pr->save_options.Saveflag)
+ {
+ errcode = saveoutput(pr);
+ pr->report.Nperiods++;
+ }
+ time->Rtime += time->Rstep;
+ }
+ if (errcode) return errcode;
- return (errcode);
+ // If simulating water quality
+ if (qual->Qualflag != NONE && qual->Qtime < time->Dur)
+ {
+ // ... compute reaction rate coeffs.
+ if (qual->Reactflag && qual->Qualflag != AGE) ratecoeffs(pr);
+
+ // ... topologically sort network nodes if flow directions change
+ if (flowdirchanged(pr) == TRUE)
+ {
+ errcode = sortnodes(pr);
+ }
+ }
+ if (!hyd->OpenHflag) time->Htime = hydtime + hydstep;
+ }
+ return errcode;
}
+
int nextqual(EN_Project *pr, long *tstep)
/*
**--------------------------------------------------------------
** Input: none
-** Output: tstep = pointer to time step (sec)
+** Output: tstep = time step (sec) over which quality was updated
** Returns: error code
-** Purpose: updates WQ conditions until next hydraulic
-** solution occurs (after *tstep secs.)
+** Purpose: updates water quality in network until next hydraulic
+** event occurs (after tstep secs.)
**--------------------------------------------------------------
*/
{
- long hydstep; /* Hydraulic solution time step */
- int errcode = 0;
- //double *tankVolumes = NULL;
- int i;
- EN_Network *net;
- hydraulics_t *hyd;
- quality_t *qu;
- time_options_t *time;
- save_options_t *sav;
-
- /* Determine time step */
- *tstep = 0;
+ long hydstep; // Time step until next hydraulic event
+ long dt, qtime;
+ int errcode = 0;
- net = &pr->network;
- hyd = &pr->hydraulics;
- qu = &pr->quality;
- time = &pr->time_options;
- sav = &pr->save_options;
+ quality_t *qual = &pr->quality;
+ time_options_t *time = &pr->time_options;
- // hydstep = time->Htime - qu->Qtime;
-
- if (time->Htime <= time->Dur) {
- hydstep = time->Htime - qu->Qtime;
- } else {
+ // Find time step till next hydraulic event
+ *tstep = 0;
hydstep = 0;
- }
+ if (time->Htime <= time->Dur) hydstep = time->Htime - qual->Qtime;
- // if we're operating in stepwise mode, capture the tank levels so we can
- // restore them later.
- if (hyd->OpenHflag) {
- //tankVolumes = calloc(net->Ntanks, sizeof(double));
- for (i = 1; i <= net->Ntanks; ++i) {
- if (net->Tank[i].A != 0) { // skip reservoirs
- qu->QTempVolumes[i - 1] = net->Tank[i].V;
- }
+ // Perform water quality routing over this time step
+ if (qual->Qualflag != NONE && hydstep > 0)
+ {
+ // Repeat over each quality time step until tstep is reached
+ qtime = 0;
+ while (!qual->OutOfMemory && qtime < hydstep)
+ {
+ dt = MIN(qual->Qstep, hydstep - qtime);
+ qtime += dt;
+ transport(pr, dt);
+ }
+ if (qual->OutOfMemory) errcode = 101;
}
- // restore the previous step's tank volumes
- for (i = 1; i <= net->Ntanks; i++) {
- if (net->Tank[i].A != 0) { // skip reservoirs again
- int n = net->Tank[i].Node;
- net->Tank[i].V = qu->QTankVolumes[i - 1];
- hyd->NodeHead[n] = tankgrade(pr, i, net->Tank[i].V);
- }
+ // Update mass balance ratio
+ evalmassbalance(pr);
+
+ // Update current time
+ if (!errcode) *tstep = hydstep;
+ qual->Qtime += hydstep;
+
+ // If no more time steps remain
+ if (!errcode && *tstep == 0)
+ {
+ // ... report overall mass balance
+ if (qual->Qualflag != NONE && pr->report.Statflag)
+ {
+ writemassbalance(pr);
+ }
+
+ // ... write the final portion of the binary output file
+ if (pr->save_options.Saveflag) errcode = savefinaloutput(pr);
}
-
- // restore the previous step's pipe link flows
- for (i = 1; i <= net->Nlinks; i++) {
- if (hyd->LinkStatus[i] <= CLOSED) {
- hyd->LinkFlows[i] = 0.0;
- }
- }
- }
-
- /* Perform water quality routing over this time step */
- if (qu->Qualflag != NONE && hydstep > 0)
- transport(pr,hydstep);
-
- /* Update current time */
- if (qu->OutOfMemory)
- errcode = 101;
- if (!errcode)
- *tstep = hydstep;
- qu->Qtime += hydstep;
-
- /* Save final output if no more time steps */
- if (!errcode && sav->Saveflag && *tstep == 0)
- errcode = savefinaloutput(pr);
-
- // restore tank levels to post-runH state, if needed.
- if (hyd->OpenHflag) {
- for (i = 1; i <= net->Ntanks; i++) {
- if (net->Tank[i].A != 0) { // skip reservoirs again
- int n = net->Tank[i].Node;
- net->Tank[i].V = qu->QTempVolumes[i - 1];
- hyd->NodeHead[n] = tankgrade(pr, i, net->Tank[i].V);
- }
- }
-
- for (i = 1; i <= net->Nlinks; ++i) {
- if (hyd->LinkStatus[i] <= CLOSED) {
- hyd->LinkFlows[i] = qu->QLinkFlow[i - 1];
- }
- }
-
- //free(tankVolumes);
- }
-
- return (errcode);
+ return errcode;
}
+
int stepqual(EN_Project *pr, long *tleft)
/*
**--------------------------------------------------------------
** Input: none
-** Output: tleft = pointer to time left in simulation
+** Output: tleft = time left in simulation
** Returns: error code
-** Purpose: updates WQ conditions over a single WQ time step
+** Purpose: updates quality conditions over a single
+** quality time step
**--------------------------------------------------------------
*/
{
- quality_t *qu = &pr->quality;
- time_options_t *time = &pr->time_options;
- save_options_t *sav = &pr->save_options;
+ long dt, hstep, t, tstep;
+ int errcode = 0;
- long dt, hstep, t, tstep;
- int errcode = 0;
- tstep = qu->Qstep;
- do {
- dt = tstep;
- hstep = time->Htime - qu->Qtime;
- if (hstep < dt) {
- dt = hstep;
- if (qu->Qualflag != NONE) {
- transport(pr,dt);
- }
- qu->Qtime += dt;
- errcode = runqual(pr,&t);
- qu->Qtime = t;
- } else {
- if (qu->Qualflag != NONE)
- transport(pr,dt);
- qu->Qtime += dt;
- }
- tstep -= dt;
- if (qu->OutOfMemory) {
- errcode = 101;
- }
- } while (!errcode && tstep > 0);
- *tleft = time->Dur - qu->Qtime;
+ quality_t *qual = &pr->quality;
- if (!errcode && sav->Saveflag && *tleft == 0) {
- errcode = savefinaloutput(pr);
- }
- return (errcode);
+ tstep = qual->Qstep;
+ do
+ {
+ // Set local time step to quality time step
+ dt = tstep;
+
+ // Find time step until next hydraulic event
+ hstep = pr->time_options.Htime - qual->Qtime;
+
+ // If next hydraulic event occurs before end of local time step
+ if (hstep < dt)
+ {
+ // ... adjust local time step to next hydraulic event
+ dt = hstep;
+
+ // ... transport quality over local time step
+ if (qual->Qualflag != NONE) transport(pr, dt);
+ qual->Qtime += dt;
+
+ // ... quit if running quality concurrently with hydraulics
+ if (pr->hydraulics.OpenHflag) break;
+
+ // ... otherwise call runqual() to update hydraulics
+ errcode = runqual(pr, &t);
+ qual->Qtime = t;
+ }
+
+ // Otherwise transport quality over current local time step
+ else
+ {
+ if (qual->Qualflag != NONE) transport(pr, dt);
+ qual->Qtime += dt;
+ }
+
+ // Reduce quality time step by local time step
+ tstep -= dt;
+ if (qual->OutOfMemory) errcode = 101;
+
+ } while (!errcode && tstep > 0);
+
+ // Update mass balance ratio
+ evalmassbalance(pr);
+
+ // Update total simulation time left
+ *tleft = pr->time_options.Dur - qual->Qtime;
+
+ // If no more time steps remain
+ if (!errcode && *tleft == 0)
+ {
+ // ... report overall mass balance
+ if (qual->Qualflag != NONE && pr->report.Statflag)
+ {
+ writemassbalance(pr);
+ }
+
+ // ... write the final portion of the binary output file
+ if (pr->save_options.Saveflag) errcode = savefinaloutput(pr);
+ }
+ return errcode;
}
+
int closequal(EN_Project *pr)
/*
**--------------------------------------------------------------
** Input: none
** Output: returns error code
-** Purpose: closes WQ solver system
+** Purpose: closes water quality solver
**--------------------------------------------------------------
*/
{
- quality_t *qu = &pr->quality;
- int errcode = 0;
+ quality_t *qual = &pr->quality;
+ int errcode = 0;
- /* Free memory pool */
- if (qu->SegPool)
- {
- mempool_delete(qu->SegPool);
- }
-
- free(qu->FirstSeg);
- free(qu->LastSeg);
- free(qu->FlowDir);
- free(qu->VolIn);
- free(qu->MassIn);
- free(qu->PipeRateCoeff);
- free(qu->TempQual);
- free(qu->QTempVolumes);
- free(qu->QTankVolumes);
- free(qu->QLinkFlow);
- return (errcode);
+ if (qual->SegPool) mempool_delete(qual->SegPool);
+ FREE(qual->FirstSeg);
+ FREE(qual->LastSeg);
+ FREE(qual->PipeRateCoeff);
+ FREE(qual->FlowDir);
+ FREE(qual->Ilist);
+ FREE(qual->IlistPtr);
+ FREE(qual->SortedNodes);
+ return errcode;
}
-int gethyd(EN_Project *pr, long *hydtime, long *hydstep)
-/*
-**-----------------------------------------------------------
-** Input: none
-** Output: hydtime = pointer to hydraulic solution time
-** hydstep = pointer to hydraulic time step
-** Returns: error code
-** Purpose: retrieves hydraulic solution and hydraulic
-** time step for next hydraulic event
-**
-** NOTE: when this function is called, WQ results have
-** already been updated to the point in time when
-** the next hydraulic event occurs.
-**-----------------------------------------------------------
-*/
-{
- int errcode = 0;
-
- hydraulics_t *hyd = &pr->hydraulics;
- report_options_t *rep = &pr->report;
- quality_t *qu = &pr->quality;
- time_options_t *time = &pr->time_options;
- save_options_t *sav = &pr->save_options;
-
- // if hydraulics are not open, then we're operating in sequential mode.
- // else hydraulics are open, so use the hydraulic results in memory rather
- // than reading from the temp file.
- if (!hyd->OpenHflag) {
- /* Read hydraulic results from file */
- if (!readhyd(pr,hydtime)) {
- return (307);
- }
- if (!readhydstep(pr->out_files.HydFile, hydstep)) {
- return (307);
- }
- time->Htime = *hydtime;
- }
-
- /* Save current results to output file */
- if (time->Htime >= time->Rtime) {
- if (sav->Saveflag) {
- errcode = saveoutput(pr);
- rep->Nperiods++;
- }
- time->Rtime += time->Rstep;
- }
-
- /* If simulating WQ: */
- if (qu->Qualflag != NONE && qu->Qtime < time->Dur) {
-
- /* Compute reaction rate coeffs. */
- if (qu->Reactflag && qu->Qualflag != AGE) {
- ratecoeffs(pr);
- }
-
- /* Initialize pipe segments (at time 0) or */
- /* else re-orient segments if flow reverses.*/
- // if (qu->Qtime == 0)
- // initsegs();
- // else
- // if hydraulics are open, or if we're in sequential mode (where qtime can
- // increase)
- if (hyd->OpenHflag || qu->Qtime != 0) {
- reorientsegs(pr);
- }
- }
- return (errcode);
-}
-
-char setReactflag(EN_Project *pr)
-/*
-**-----------------------------------------------------------
-** Input: none
-** Output: returns 1 for reactive WQ constituent, 0 otherwise
-** Purpose: checks if reactive chemical being simulated
-**-----------------------------------------------------------
-*/
-{
- EN_Network *net = &pr->network;
- quality_t *qu = &pr->quality;
-
- int i;
- if (qu->Qualflag == TRACE) {
- return (0);
- } else if (qu->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);
-}
-
-void transport(EN_Project *pr, long tstep)
-/*
-**--------------------------------------------------------------
-** Input: tstep = length of current time step
-** Output: none
-** Purpose: transports constituent mass through pipe network
-** under a period of constant hydraulic conditions.
-**--------------------------------------------------------------
-*/
-{
- long qtime, dt;
- quality_t *qu = &pr->quality;
-
- /* Repeat until elapsed time equals hydraulic time step */
-
- qtime = 0;
- while (!qu->OutOfMemory && qtime < tstep) { /* Qstep is quality time step */
- dt = MIN(qu->Qstep, tstep - qtime); /* Current time step */
- qtime += dt; /* Update elapsed time */
- if (qu->Reactflag) {
- updatesegs(pr, dt); /* Update quality in inner link segs */
- }
- accumulate(pr, dt); /* Accumulate flow at nodes */
- updatenodes(pr, dt); /* Update nodal quality */
- sourceinput(pr, dt); /* Compute inputs from sources */
- release(pr, dt); /* Release new nodal flows */
- }
- updatesourcenodes(pr, tstep); /* Update quality at source nodes */
-}
-
-void initsegs(EN_Project *pr)
-/*
-**--------------------------------------------------------------
-** Input: none
-** Output: none
-** Purpose: initializes water quality segments
-**--------------------------------------------------------------
-*/
-{
- int j, k;
- double c, v;
- EN_Network *net = &pr->network;
- hydraulics_t *hyd = &pr->hydraulics;
- quality_t *qu = &pr->quality;
-
- /* Examine each link */
- for (k = 1; k <= net->Nlinks; k++) {
-
- /* Establish flow direction */
- qu->FlowDir[k] = POSITIVE;
- if (hyd->LinkFlows[k] < 0.) {
- qu->FlowDir[k] = NEGATIVE;
- }
-
- /* Set segs to zero */
- qu->LastSeg[k] = NULL;
- qu->FirstSeg[k] = NULL;
-
- /* Find quality of downstream node */
- j = DOWN_NODE(k);
- if (j <= net->Njuncs) {
- c = qu->NodeQual[j];
- } else {
- c = net->Tank[j - net->Njuncs].C;
- }
-
- /* Fill link with single segment with this quality */
- addseg(pr, k, LINKVOL(k), c);
- }
-
- /* Initialize segments in tanks that use them */
- for (j = 1; j <= net->Ntanks; j++) {
-
- /* Skip reservoirs & complete mix tanks */
- if (net->Tank[j].A == 0.0 || net->Tank[j].MixModel == MIX1)
- continue;
-
- /* Tank segment pointers are stored after those for links */
- k = net->Nlinks + j;
- c = net->Tank[j].C;
- qu->LastSeg[k] = NULL;
- qu->FirstSeg[k] = NULL;
-
- /* Add 2 segments for 2-compartment model */
- if (net->Tank[j].MixModel == MIX2) {
- v = MAX(0, net->Tank[j].V - net->Tank[j].V1max);
- addseg(pr, k, v, c);
- v = net->Tank[j].V - v;
- addseg(pr, k, v, c);
- }
-
- /* Add one segment for FIFO & LIFO models */
- else {
- v = net->Tank[j].V;
- addseg(pr, k, v, c);
- }
- }
-}
-
-void reorientsegs(EN_Project *pr)
-/*
-**--------------------------------------------------------------
-** Input: none
-** Output: none
-** Purpose: re-orients segments (if flow reverses)
-**--------------------------------------------------------------
-*/
-{
- Pseg seg, nseg, pseg;
- int k;
- FlowDirection newdir;
- EN_Network *net = &pr->network;
- hydraulics_t *hyd = &pr->hydraulics;
- quality_t *qu = &pr->quality;
-
- /* Examine each link */
- for (k = 1; k <= net->Nlinks; k++) {
-
- /* Find new flow direction */
- newdir = POSITIVE;
- if (hyd->LinkFlows[k] == 0.0) {
- newdir = qu->FlowDir[k];
- } else if (hyd->LinkFlows[k] < 0.0) {
- newdir = NEGATIVE;
- }
-
- /* If direction changes, then reverse order of segments */
- /* (first to last) and save new direction */
- if (newdir != qu->FlowDir[k]) {
- seg = qu->FirstSeg[k];
- qu->FirstSeg[k] = qu->LastSeg[k];
- qu->LastSeg[k] = seg;
- pseg = NULL;
- while (seg != NULL) {
- nseg = seg->prev;
- seg->prev = pseg;
- pseg = seg;
- seg = nseg;
- }
- qu->FlowDir[k] = newdir;
- }
- }
-}
-
-void updatesegs(EN_Project *pr, long dt)
-/*
-**-------------------------------------------------------------
-** Input: t = time from last WQ segment update
-** Output: none
-** Purpose: reacts material in pipe segments up to time t
-**-------------------------------------------------------------
-*/
-{
- int k;
- Pseg seg;
- double cseg, rsum, vsum;
- EN_Network *net = &pr->network;
- quality_t *qu = &pr->quality;
-
- /* Examine each link in network */
- for (k = 1; k <= net->Nlinks; k++) {
-
- /* Skip zero-length links (pumps & valves) */
- rsum = 0.0;
- vsum = 0.0;
- if (net->Link[k].Len == 0.0) {
- continue;
- }
-
- /* Examine each segment of the link */
- seg = qu->FirstSeg[k];
- while (seg != NULL) {
-
- /* React segment over time dt */
- cseg = seg->c;
- seg->c = pipereact(pr, k, seg->c, seg->v, dt);
-
- /* Accumulate volume-weighted reaction rate */
- if (qu->Qualflag == CHEM) {
- rsum += ABS((seg->c - cseg)) * seg->v;
- vsum += seg->v;
- }
- seg = seg->prev;
- }
-
- /* Normalize volume-weighted reaction rate */
- if (vsum > 0.0) {
- qu->PipeRateCoeff[k] = rsum / vsum / dt * SECperDAY;
- } else
- qu->PipeRateCoeff[k] = 0.0;
- }
-}
-
-void removesegs(EN_Project *pr, int k)
-/*
-**-------------------------------------------------------------
-** Input: k = link index
-** Output: none
-** Purpose: removes all segments in link k
-**-------------------------------------------------------------
-*/
-{
- Pseg seg;
- quality_t *qu = &pr->quality;
-
- seg = qu->FirstSeg[k];
- while (seg != NULL) {
- qu->FirstSeg[k] = seg->prev;
- seg->prev = qu->FreeSeg;
- qu->FreeSeg = seg;
- seg = qu->FirstSeg[k];
- }
- qu->LastSeg[k] = NULL;
-}
-
-void addseg(EN_Project *pr, int k, double v, double c)
-/*
-**-------------------------------------------------------------
-** Input: k = link segment
-** v = segment volume
-** c = segment quality
-** Output: none
-** Purpose: adds a segment to start of link k (i.e., upstream
-** of current last segment).
-**-------------------------------------------------------------
-*/
-{
- Pseg seg;
- quality_t *qu = &pr->quality;
-
- if (qu->FreeSeg != NULL) {
- seg = qu->FreeSeg;
- qu->FreeSeg = seg->prev;
- } else {
- seg = (struct Sseg *) mempool_alloc(qu->SegPool, sizeof(struct Sseg));
- if (seg == NULL) {
- qu->OutOfMemory = TRUE;
- return;
- }
- }
- seg->v = v;
- seg->c = c;
- seg->prev = NULL;
- if (qu->FirstSeg[k] == NULL) {
- qu->FirstSeg[k] = seg;
- }
- if (qu->LastSeg[k] != NULL) {
- qu->LastSeg[k]->prev = seg;
- }
- qu->LastSeg[k] = seg;
-}
-
-void accumulate(EN_Project *pr, long dt)
-/*
-**-------------------------------------------------------------
-** Input: dt = current WQ time step
-** Output: none
-** Purpose: accumulates mass flow at nodes and updates nodal
-** quality
-**-------------------------------------------------------------
-*/
-{
- int i, j, k;
- double cseg, v, vseg;
- Pseg seg;
-
- EN_Network *net = &pr->network;
- hydraulics_t *hyd = &pr->hydraulics;
- quality_t *qu = &pr->quality;
-
- /* Re-set memory used to accumulate mass & volume */
- memset(qu->VolIn, 0, (net->Nnodes + 1) * sizeof(double));
- memset(qu->MassIn, 0, (net->Nnodes + 1) * sizeof(double));
- memset(qu->TempQual, 0, (net->Nnodes + 1) * sizeof(double));
-
- /* Compute average conc. of segments adjacent to each node */
- /* (For use if there is no transport through the node) */
- for (k = 1; k <= net->Nlinks; k++) {
- j = DOWN_NODE(k); /* Downstream node */
- if (qu->FirstSeg[k] != NULL) /* Accumulate concentrations */
- {
- qu->MassIn[j] += qu->FirstSeg[k]->c;
- qu->VolIn[j]++;
- }
- j = UP_NODE(k); /* Upstream node */
- if (qu->LastSeg[k] != NULL) /* Accumulate concentrations */
- {
- qu->MassIn[j] += qu->LastSeg[k]->c;
- qu->VolIn[j]++;
- }
- }
-
- for (k = 1; k <= net->Nnodes; k++) {
- if (qu->VolIn[k] > 0.0) {
- qu->TempQual[k] = qu->MassIn[k] / qu->VolIn[k];
- }
- }
-
- /* Move mass from first segment of each pipe into downstream node */
- memset(qu->VolIn, 0, (net->Nnodes + 1) * sizeof(double));
- memset(qu->MassIn, 0, (net->Nnodes + 1) * sizeof(double));
- for (k = 1; k <= net->Nlinks; k++) {
- i = UP_NODE(k); /* Upstream node */
- j = DOWN_NODE(k); /* Downstream node */
- v = ABS(hyd->LinkFlows[k]) * dt; /* Flow volume */
-
- while (v > 0.0)
- {
- /* Identify leading segment in pipe */
- seg = qu->FirstSeg[k];
- if (seg == NULL)
- break;
-
- /* Volume transported from this segment is */
- /* minimum of flow volume & segment volume */
- /* (unless leading segment is also last segment) */
- vseg = seg->v;
- vseg = MIN(vseg, v);
- if (seg == qu->LastSeg[k])
- vseg = v;
-
- /* Update volume & mass entering downstream node */
- cseg = seg->c;
- qu->VolIn[j] += vseg;
- qu->MassIn[j] += vseg * cseg;
-
- /* Reduce flow volume by amount transported */
- v -= vseg;
-
- /* If all of segment's volume was transferred, then */
- /* replace leading segment with the one behind it */
- /* (Note that the current seg is recycled for later use.) */
- if (v >= 0.0 && vseg >= seg->v) {
- qu->FirstSeg[k] = seg->prev;
- if (qu->FirstSeg[k] == NULL)
- qu->LastSeg[k] = NULL;
- seg->prev = qu->FreeSeg;
- qu->FreeSeg = seg;
- }
-
- /* Otherwise reduce segment's volume */
- else {
- seg->v -= vseg;
- }
- } /* End while */
- } /* Next link */
-}
-
-void updatenodes(EN_Project *pr, long dt)
-/*
-**---------------------------------------------------------------------------
-** Input: dt = current WQ time step
-** Output: none
-** Purpose: updates concentration at all nodes to mixture of accumulated
-** inflow from connecting pipes.
-**
-** Note: Does not account for source flow effects. TempQual[i] contains
-** average concen. of segments adjacent to node i, used in case
-** there was no inflow into i.
-**---------------------------------------------------------------------------
-*/
-{
- int i;
- EN_Network *net = &pr->network;
- hydraulics_t *hyd = &pr->hydraulics;
- quality_t *qu = &pr->quality;
-
- /* Update junction quality */
- for (i = 1; i <= net->Njuncs; i++) {
- if (hyd->NodeDemand[i] < 0.0) {
- qu->VolIn[i] -= hyd->NodeDemand[i] * dt;
- }
- if (qu->VolIn[i] > 0.0) {
- qu->NodeQual[i] = qu->MassIn[i] / qu->VolIn[i];
- } else {
- qu->NodeQual[i] = qu->TempQual[i];
- }
- }
-
- /* Update tank quality */
- updatetanks(pr, dt);
-
- /* For flow tracing, set source node concen. to 100. */
- if (qu->Qualflag == TRACE)
- qu->NodeQual[qu->TraceNode] = 100.0;
-}
-
-void sourceinput(EN_Project *pr, long dt)
-/*
-**---------------------------------------------------------------------
-** Input: dt = current WQ time step
-** Output: none
-** Purpose: computes contribution (if any) of mass additions from WQ
-** sources at each node.
-**---------------------------------------------------------------------
-*/
-{
- int j, n;
- double massadded = 0.0, s, volout;
- double qout, qcutoff;
- Psource source;
-
- EN_Network *net = &pr->network;
- hydraulics_t *hyd = &pr->hydraulics;
- quality_t *qu = &pr->quality;
- time_options_t *time = &pr->time_options;
-
- /* Establish a flow cutoff which indicates no outflow from a node */
- qcutoff = 10.0 * TINY;
-
- /* Zero-out the work array TempQual */
- memset(qu->TempQual, 0, (net->Nnodes + 1) * sizeof(double));
- if (qu->Qualflag != CHEM)
- return;
-
- /* Consider each node */
- for (n = 1; n <= net->Nnodes; n++) {
- double thisDemand = hyd->NodeDemand[n];
- /* Skip node if no WQ source */
- source = net->Node[n].S;
- if (source == NULL)
- continue;
- if (source->C0 == 0.0)
- continue;
-
- /* Find total flow volume leaving node */
- if (n <= net->Njuncs) {
- volout = qu->VolIn[n]; /* Junctions */
- } else {
- volout = qu->VolIn[n] - (thisDemand * dt); /* Tanks */
- }
- qout = volout / (double)dt;
-
- /* Evaluate source input only if node outflow > cutoff flow */
- if (qout > qcutoff) {
-
- /* Mass added depends on type of source */
- s = sourcequal(pr,source);
- switch (source->Type) {
- /* Concen. Source: */
- /* Mass added = source concen. * -(demand) */
- case CONCEN:
-
- /* Only add source mass if demand is negative */
- if (thisDemand < 0.0) {
- massadded = -s * thisDemand * dt;
-
- /* If node is a tank then set concen. to 0. */
- /* (It will be re-set to true value in updatesourcenodes()) */
- if (n > net->Njuncs)
- qu->NodeQual[n] = 0.0;
- } else
- massadded = 0.0;
- break;
-
- /* Mass Inflow Booster Source: */
- case MASS:
- massadded = s * dt;
- break;
-
- /* Setpoint Booster Source: */
- /* Mass added is difference between source */
- /* & node concen. times outflow volume */
- case SETPOINT:
- if (s > qu->NodeQual[n]) {
- massadded = (s - qu->NodeQual[n]) * volout;
- } else {
- massadded = 0.0;
- }
- break;
-
- /* Flow-Paced Booster Source: */
- /* Mass added = source concen. times outflow volume */
- case FLOWPACED:
- massadded = s * volout;
- break;
- }
-
- /* Source concen. contribution = (mass added / outflow volume) */
- qu->TempQual[n] = massadded / volout;
-
- /* Update total mass added for time period & simulation */
- source->Smass += massadded;
- if (time->Htime >= time->Rstart) {
- qu->Wsource += massadded;
- }
- }
- }
-
- /* Add mass inflows from reservoirs to Wsource*/
- if (time->Htime >= time->Rstart) {
- for (j = 1; j <= net->Ntanks; j++) {
- if (net->Tank[j].A == 0.0) {
- n = net->Njuncs + j;
- volout = qu->VolIn[n] - hyd->NodeDemand[n] * dt;
- if (volout > 0.0)
- qu->Wsource += volout * qu->NodeQual[n];
- }
- }
- }
-}
-
-void release(EN_Project *pr, long dt)
-/*
-**---------------------------------------------------------
-** Input: dt = current WQ time step
-** Output: none
-** Purpose: creates new segments in outflow links from nodes.
-**---------------------------------------------------------
-*/
-{
- EN_Network *net = &pr->network;
- hydraulics_t *hyd = &pr->hydraulics;
- quality_t *qu = &pr->quality;
-
- int k, n;
- double c, q, v;
- Pseg seg;
-
- /* Examine each link */
- for (k = 1; k <= net->Nlinks; k++) {
-
- /* Ignore links with no flow */
- if (hyd->LinkFlows[k] == 0.0)
- continue;
-
- /* Find flow volume released to link from upstream node */
- /* (NOTE: Flow volume is allowed to be > link volume.) */
- n = UP_NODE(k);
- q = ABS(hyd->LinkFlows[k]);
- v = q * dt;
-
- /* Include source contribution in quality released from node. */
- c = qu->NodeQual[n] + qu->TempQual[n];
-
- /* If link has a last seg, check if its quality */
- /* differs from that of the flow released from node.*/
- if ((seg = qu->LastSeg[k]) != NULL) {
- /* Quality of seg close to that of node */
- if (ABS(seg->c - c) < qu->Ctol) {
- seg->c = (seg->c * seg->v + c * v) / (seg->v + v);
- seg->v += v;
- }
-
- /* Otherwise add a new seg to end of link */
- else
- addseg(pr, k, v, c);
- }
-
- /* If link has no segs then add a new one. */
- else
- addseg(pr, k, LINKVOL(k), c);
- }
-}
-
-void updatesourcenodes(EN_Project *pr, long dt)
-/*
-**---------------------------------------------------
-** Input: dt = current WQ time step
-** Output: none
-** Purpose: updates quality at source nodes.
-** (TempQual[n] = concen. added by source at node n)
-**---------------------------------------------------
-*/
-{
- EN_Network *net = &pr->network;
- quality_t *qu = &pr->quality;
-
- int i, n;
- Psource source;
-
- if (qu->Qualflag != CHEM)
- return;
-
- /* Examine each WQ source node */
- for (n = 1; n <= net->Nnodes; n++) {
- source = net->Node[n].S;
- if (source == NULL)
- continue;
-
- /* Add source to current node concen. */
- qu->NodeQual[n] += qu->TempQual[n];
-
- /* For tanks, node concen. = internal concen. */
- if (n > net->Njuncs) {
- i = n - net->Njuncs;
- if (net->Tank[i].A > 0.0)
- qu->NodeQual[n] = net->Tank[i].C;
- }
-
- /* Normalize mass added at source to time step */
- source->Smass /= (double)dt;
- }
-}
-
-void updatetanks(EN_Project *pr, long dt)
-/*
-**---------------------------------------------------
-** Input: dt = current WQ time step
-** Output: none
-** Purpose: updates tank volumes & concentrations
-**---------------------------------------------------
-*/
-{
- EN_Network *net = &pr->network;
- quality_t *qu = &pr->quality;
-
- int i, n;
-
- /* Examine each reservoir & tank */
- for (i = 1; i <= net->Ntanks; i++) {
- n = net->Tank[i].Node;
- /* Use initial quality for reservoirs */
- if (net->Tank[i].A == 0.0) {
- qu->NodeQual[n] = net->Node[n].C0;
- }
- /* Update tank WQ based on mixing model */
- else {
- switch (net->Tank[i].MixModel) {
- case MIX2:
- tankmix2(pr, i, dt);
- break;
- case FIFO:
- tankmix3(pr, i, dt);
- break;
- case LIFO:
- tankmix4(pr, i, dt);
- break;
- default:
- tankmix1(pr, i, dt);
- break;
- }
- }
- }
-}
-
-//// New version of tankmix1 ////
-void tankmix1(EN_Project *pr, int i, long dt)
-/*
-**---------------------------------------------
-** Input: i = tank index
-** dt = current WQ time step
-** Output: none
-** Purpose: complete mix tank model
-**---------------------------------------------
-*/
-{
- EN_Network *net = &pr->network;
- hydraulics_t *hyd = &pr->hydraulics;
- quality_t *qu = &pr->quality;
-
- int n;
- double cin;
- double c, cmax, vold, vin;
-
- Stank *tank = &net->Tank[i];
-
- /* React contents of tank */
- c = tankreact(pr, tank->C, tank->V, tank->Kb, dt);
-
- /* Determine tank & volumes */
- vold = tank->V;
- n = tank->Node;
- tank->V += hyd->NodeDemand[n] * dt;
- vin = qu->VolIn[n];
-
- /* Compute inflow concen. */
- if (vin > 0.0)
- cin = qu->MassIn[n] / vin;
- else
- cin = 0.0;
- cmax = MAX(c, cin);
-
- /* Mix inflow with tank contents */
- if (vin > 0.0)
- c = (c * vold + cin * vin) / (vold + vin);
- c = MIN(c, cmax);
- c = MAX(c, 0.0);
- tank->C = c;
- qu->NodeQual[n] = tank->C;
-}
-
-/*** Updated 10/25/00 ***/
-//// New version of tankmix2 ////
-void tankmix2(EN_Project *pr, int i, long dt)
-/*
-**------------------------------------------------
-** Input: i = tank index
-** dt = current WQ time step
-** Output: none
-** Purpose: 2-compartment tank model
-** (seg1 = mixing zone,
-** seg2 = ambient zone)
-**------------------------------------------------
-*/
-{
- EN_Network *net = &pr->network;
- hydraulics_t *hyd = &pr->hydraulics;
- quality_t *qu = &pr->quality;
-
- int k, n;
- double cin, /* Inflow quality */
- vin, /* Inflow volume */
- vt, /* Transferred volume */
- vnet, /* Net volume change */
- v1max; /* Full mixing zone volume */
- Pseg seg1, seg2; /* Compartment segments */
-
- Stank *tank = &pr->network.Tank[i];
-
- /* Identify segments for each compartment */
- k = net->Nlinks + i;
- seg1 = qu->LastSeg[k];
- seg2 = qu->FirstSeg[k];
- if (seg1 == NULL || seg2 == NULL)
- return;
-
- /* React contents of each compartment */
- seg1->c = tankreact(pr, seg1->c, seg1->v, tank->Kb, dt);
- seg2->c = tankreact(pr, seg2->c, seg2->v, tank->Kb, dt);
-
- /* Find inflows & outflows */
- n = tank->Node;
- vnet = hyd->NodeDemand[n] * dt;
- vin = qu->VolIn[n];
- if (vin > 0.0)
- cin = qu->MassIn[n] / vin;
- else
- cin = 0.0;
- v1max = tank->V1max;
-
- /* Tank is filling */
- vt = 0.0;
- if (vnet > 0.0) {
- vt = MAX(0.0, (seg1->v + vnet - v1max));
- if (vin > 0.0) {
- seg1->c = ((seg1->c) * (seg1->v) + cin * vin) / (seg1->v + vin);
- }
- if (vt > 0.0) {
- seg2->c = ((seg2->c) * (seg2->v) + (seg1->c) * vt) / (seg2->v + vt);
- }
- }
-
- /* Tank is emptying */
- if (vnet < 0.0) {
- if (seg2->v > 0.0) {
- vt = MIN(seg2->v, (-vnet));
- }
- if (vin + vt > 0.0) {
- seg1->c = ((seg1->c) * (seg1->v) + cin * vin + (seg2->c) * vt) /
- (seg1->v + vin + vt);
- }
- }
-
- /* Update segment volumes */
- if (vt > 0.0) {
- seg1->v = v1max;
- if (vnet > 0.0)
- seg2->v += vt;
- else
- seg2->v = MAX(0.0, ((seg2->v) - vt));
- } else {
- seg1->v += vnet;
- seg1->v = MIN(seg1->v, v1max);
- seg1->v = MAX(0.0, seg1->v);
- seg2->v = 0.0;
- }
- tank->V += vnet;
- tank->V = MAX(0.0, tank->V);
-
- /* Use quality of mixed compartment (seg1) to */
- /* represent quality of tank since this is where */
- /* outflow begins to flow from */
- tank->C = seg1->c;
- qu->NodeQual[n] = tank->C;
-}
-
-void tankmix3(EN_Project *pr, int i, long dt)
-/*
-**----------------------------------------------------------
-** Input: i = tank index
-** dt = current WQ time step
-** Output: none
-** Purpose: First-In-First-Out (FIFO) tank model
-**----------------------------------------------------------
-*/
-{
- EN_Network *net = &pr->network;
- hydraulics_t *hyd = &pr->hydraulics;
- quality_t *qu = &pr->quality;
-
- int k, n;
- double vin, vnet, vout, vseg;
- double cin, vsum, csum;
- Pseg seg;
- Stank *tank = &pr->network.Tank[i];
- k = net->Nlinks + i;
- if (qu->LastSeg[k] == NULL || qu->FirstSeg[k] == NULL)
- return;
-
- /* React contents of each compartment */
- if (qu->Reactflag) {
- seg = qu->FirstSeg[k];
- while (seg != NULL) {
- seg->c = tankreact(pr, seg->c, seg->v, tank->Kb, dt);
- seg = seg->prev;
- }
- }
-
- /* Find inflows & outflows */
- n = tank->Node;
- vnet = hyd->NodeDemand[n] * dt;
- vin = qu->VolIn[n];
- vout = vin - vnet;
- if (vin > 0.0)
- cin = qu->MassIn[n] / qu->VolIn[n];
- else
- cin = 0.0;
- tank->V += vnet;
- tank->V = MAX(0.0, tank->V);
-
- /* Withdraw flow from first segment */
- vsum = 0.0;
- csum = 0.0;
- while (vout > 0.0) {
- seg = qu->FirstSeg[k];
- if (seg == NULL)
- break;
- vseg = seg->v; /* Flow volume from leading seg */
- vseg = MIN(vseg, vout);
- if (seg == qu->LastSeg[k])
- vseg = vout;
- vsum += vseg;
- csum += (seg->c) * vseg;
- vout -= vseg; /* Remaining flow volume */
- if (vout >= 0.0 && vseg >= seg->v) /* Seg used up */
- {
- if (seg->prev)
- {
- qu->FirstSeg[k] = seg->prev;
- seg->prev = qu->FreeSeg;
- qu->FreeSeg = seg;
- }
- } else /* Remaining volume in segment */
- {
- seg->v -= vseg;
- }
- }
-
- /* Use quality withdrawn from 1st segment */
- /* to represent overall quality of tank */
- if (vsum > 0.0)
- tank->C = csum / vsum;
- else
- tank->C = qu->FirstSeg[k]->c;
- qu->NodeQual[n] = tank->C;
-
- /* Add new last segment for new flow entering tank */
- if (vin > 0.0) {
- if ((seg = qu->LastSeg[k]) != NULL) {
- /* Quality is the same, so just add flow volume to last seg */
- if (ABS(seg->c - cin) < qu->Ctol)
- seg->v += vin;
-
- /* Otherwise add a new seg to tank */
- else
- addseg(pr, k, vin, cin);
- }
-
- /* If no segs left then add a new one. */
- else
- addseg(pr, k, vin, cin);
- }
-}
-
-void tankmix4(EN_Project *pr, int i, long dt)
-/*
-**----------------------------------------------------------
-** Input: i = tank index
-** dt = current WQ time step
-** Output: none
-** Purpose: Last In-First Out (LIFO) tank model
-**----------------------------------------------------------
-*/
-{
- EN_Network *net = &pr->network;
- hydraulics_t *hyd = &pr->hydraulics;
- quality_t *qu = &pr->quality;
-
- int k, n;
- double vin, vnet, cin, vsum, csum, vseg;
- Pseg seg, tmpseg;
- Stank *tank = &pr->network.Tank[i];
- k = net->Nlinks + i;
- if (qu->LastSeg[k] == NULL || qu->FirstSeg[k] == NULL)
- return;
-
- /* React contents of each compartment */
- if (qu->Reactflag) {
- seg = qu->LastSeg[k];
- while (seg != NULL) {
- seg->c = tankreact(pr, seg->c, seg->v, tank->Kb, dt);
- seg = seg->prev;
- }
- }
-
- /* Find inflows & outflows */
- n = tank->Node;
- vnet = hyd->NodeDemand[n] * dt;
- vin = qu->VolIn[n];
- if (vin > 0.0)
- cin = qu->MassIn[n] / qu->VolIn[n];
- else
- cin = 0.0;
- tank->V += vnet;
- tank->V = MAX(0.0, tank->V);
- tank->C = qu->LastSeg[k]->c;
-
- /* If tank filling, then create new last seg */
- if (vnet > 0.0) {
- if ((seg = qu->LastSeg[k]) != NULL) {
- /* Quality is the same, so just add flow volume to last seg */
- if (ABS(seg->c - cin) < qu->Ctol)
- seg->v += vnet;
-
- /* Otherwise add a new last seg to tank */
- /* which points to old last seg */
- else {
- tmpseg = seg;
- qu->LastSeg[k] = NULL;
- addseg(pr, k, vnet, cin);
- qu->LastSeg[k]->prev = tmpseg;
- }
- }
-
- /* If no segs left then add a new one. */
- else
- addseg(pr, k, vnet, cin);
-
- /* Update reported tank quality */
- tank->C = qu->LastSeg[k]->c;
- }
-
- /* If net emptying then remove last segments until vnet consumed */
- else if (vnet < 0.0) {
- vsum = 0.0;
- csum = 0.0;
- vnet = -vnet;
- while (vnet > 0.0) {
- seg = qu->LastSeg[k];
- if (seg == NULL)
- break;
- vseg = seg->v;
- vseg = MIN(vseg, vnet);
- if (seg == qu->FirstSeg[k])
- vseg = vnet;
- vsum += vseg;
- csum += (seg->c) * vseg;
- vnet -= vseg;
- if (vnet >= 0.0 && vseg >= seg->v) /* Seg used up */
- {
- if (seg->prev)
- {
- qu->LastSeg[k] = seg->prev;
- seg->prev = qu->FreeSeg;
- qu->FreeSeg = seg;
- }
- } else /* Remaining volume in segment */
- {
- seg->v -= vseg;
- }
- }
- /* Reported tank quality is mixture of flow released and any inflow */
- tank->C = (csum + qu->MassIn[n]) / (vsum + vin);
- }
- qu->NodeQual[n] = tank->C;
-}
-
-double sourcequal(EN_Project *pr, Psource source)
-/*
-**--------------------------------------------------------------
-** Input: j = source index
-** Output: returns source WQ value
-** Purpose: determines source concentration in current time period
-**--------------------------------------------------------------
-*/
-{
- int i;
- long k;
- double c;
- EN_Network *net = &pr->network;
- quality_t *qu = &pr->quality;
- time_options_t *time = &pr->time_options;
-
- /* Get source concentration (or mass flow) in original units */
- c = source->C0;
-
- /* Convert mass flow rate from min. to sec. */
- /* and convert concen. from liters to cubic feet */
- if (source->Type == MASS)
- c /= 60.0;
- else
- c /= pr->Ucf[QUALITY];
-
- /* Apply time pattern if assigned */
- i = source->Pat;
- if (i == 0)
- return (c);
- k = ((qu->Qtime + time->Pstart) / time->Pstep) % (long)net->Pattern[i].Length;
- return (c * net->Pattern[i].F[k]);
-}
double avgqual(EN_Project *pr, int k)
/*
**--------------------------------------------------------------
** Input: k = link index
-** Output: returns WQ value
-** Purpose: computes average quality in link k
+** Output: returns quality concentration
+** Purpose: computes current average quality in link k
**--------------------------------------------------------------
*/
{
- double vsum = 0.0, msum = 0.0;
- Pseg seg;
- EN_Network *net = &pr->network;
- quality_t *qu = &pr->quality;
+ double vsum = 0.0, msum = 0.0;
+ Pseg seg;
- if (qu->Qualflag == NONE)
- return (0.);
- seg = qu->FirstSeg[k];
- while (seg != NULL) {
- vsum += seg->v;
- msum += (seg->c) * (seg->v);
- seg = seg->prev;
- }
- if (vsum > 0.0 && qu->Qtime > 0)
- return (msum / vsum);
- else
- return ((qu->NodeQual[net->Link[k].N1] + qu->NodeQual[net->Link[k].N2]) /
- 2.);
+ EN_Network *net = &pr->network;
+ quality_t *qual = &pr->quality;
+
+ if (qual->Qualflag == NONE) return 0.0;
+
+ // Sum up the quality and volume in each segment of the link
+ if (qual->FirstSeg != NULL)
+ {
+ seg = qual->FirstSeg[k];
+ while (seg != NULL)
+ {
+ vsum += seg->v;
+ msum += (seg->c) * (seg->v);
+ seg = seg->prev;
+ }
+ }
+
+ // Compute average quality if link has volume
+ if (vsum > 0.0) return (msum / vsum);
+
+ // Otherwise use the average quality of the link's end nodes
+ else
+ {
+ return ((qual->NodeQual[net->Link[k].N1] +
+ qual->NodeQual[net->Link[k].N2]) / 2.);
+ }
}
-void ratecoeffs(EN_Project *pr)
+
+double findsourcequal(EN_Project *pr, int n, double volin, double volout, long tstep)
+/*
+**---------------------------------------------------------------------
+** Input: n = node index
+** volin = volume of node inflow over time step
+** volout = volume of node outflow over time step
+** tstep = current quality time step
+** Output: returns concentration added by an external quality source.
+** Purpose: computes contribution (if any) of mass addition from an
+** external quality source at a node.
+**---------------------------------------------------------------------
+*/
+{
+ double massadded = 0.0, c;
+ Psource source;
+
+ EN_Network *net = &pr->network;
+ hydraulics_t *hyd = &pr->hydraulics;
+ quality_t *qual = &pr->quality;
+ time_options_t *time = &pr->time_options;
+
+ // Sources only apply to CHEMICAL analyses
+ if (qual->Qualflag != CHEM) return 0.0;
+
+ // Return 0 if node is not a quality source or has no outflow
+ source = net->Node[n].S;
+ if (source == NULL) return 0.0;
+ if (source->C0 == 0.0) return 0.0;
+ if (volout / tstep <= QZERO) return 0.0;
+
+ // Added source concentration depends on source type
+ c = sourcequal(pr, source);
+ switch (source->Type)
+ {
+ // Concentration Source:
+ case CONCEN:
+ if (net->Node[n].Type == EN_JUNCTION)
+ {
+ // ... source requires a negative demand at the node
+ if (hyd->NodeDemand[n] < 0.0)
+ {
+ c = -c * hyd->NodeDemand[n] * tstep / volout;
+ }
+ else c = 0.0;
+ }
+ break;
+
+ // Mass Inflow Booster Source:
+ case MASS:
+ // ... convert source input from mass/sec to concentration
+ c = c * tstep / volout;
+ break;
+
+ // Setpoint Booster Source:
+ // Source quality is difference between source strength
+ // & node quality
+ case SETPOINT:
+ c = MAX(c - qual->NodeQual[n], 0.0);
+ break;
+
+ // Flow-Paced Booster Source:
+ // Source quality equals source strength
+ case FLOWPACED:
+ break;
+ }
+
+ // Source mass added over time step = source concen. * outflow volume
+ massadded = c * volout;
+
+ // Update source's total mass added
+ source->Smass += massadded;
+
+ // Update Wsource
+ if (time->Htime >= time->Rstart)
+ {
+ qual->Wsource += massadded;
+ }
+ return c;
+}
+
+
+double sourcequal(EN_Project *pr, Psource source)
+/*
+**--------------------------------------------------------------
+** Input: source = a water quality source object
+** Output: returns strength of quality source
+** Purpose: determines source strength in current time period
+**--------------------------------------------------------------
+*/
+{
+ int i;
+ long k;
+ double c;
+
+ EN_Network *net = &pr->network;
+ quality_t *qual = &pr->quality;
+ time_options_t *time = &pr->time_options;
+
+ // Get source concentration (or mass flow) in original units
+ c = source->C0;
+
+ // Convert mass flow rate from min. to sec.
+ // and convert concen. from liters to cubic feet
+ if (source->Type == MASS) c /= 60.0;
+ else c /= pr->Ucf[QUALITY];
+
+ // Apply time pattern if assigned
+ i = source->Pat;
+ if (i == 0) return c;
+ k = ((qual->Qtime + time->Pstart) / time->Pstep) %
+ (long)net->Pattern[i].Length;
+ return (c * net->Pattern[i].F[k]);
+}
+
+
+void evalmassbalance(EN_Project *pr)
/*
**--------------------------------------------------------------
** Input: none
** Output: none
-** Purpose: determines wall reaction coeff. for each pipe
+** Purpose: computes the overall mass balance ratio of a
+** quality constituent.
**--------------------------------------------------------------
*/
{
- int k;
- double kw;
- EN_Network *net = &pr->network;
- quality_t *qu = &pr->quality;
+ double massin;
+ double massout;
+ double massreacted;
+ 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;
- qu->PipeRateCoeff[k] = 0.0;
- }
-} /* End of ratecoeffs */
-
-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, kf, kw, y, Re, Sh;
-
- EN_Network *net = &pr->network;
- hydraulics_t *hyd = &pr->hydraulics;
- quality_t *qu = &pr->quality;
-
- d = net->Link[k].Diam; /* Pipe diameter, ft */
-
- /* Ignore mass transfer if Schmidt No. is 0 */
- if (qu->Sc == 0.0) {
- if (qu->WallOrder == 0.0)
- return (BIG);
+ if (qual->Qualflag == NONE) qual->massbalance.ratio = 1.0;
else
- return (net->Link[k].Kw * (4.0 / d) / pr->Ucf[ELEV]);
- }
-
- /* Compute Reynolds No. */
- a = PI * d * d / 4.0;
- u = ABS(hyd->LinkFlows[k]) / a;
- Re = u * d / hyd->Viscos;
-
- /* 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(qu->Sc, 0.333);
-
- /* Compute Sherwood No. for laminar flow */
- /* using Graetz solution formula. */
- else {
- y = d / net->Link[k].Len * Re * qu->Sc;
- Sh = 3.65 + 0.0668 * y / (1.0 + 0.04 * pow(y, 0.667));
- }
-
- /* Compute mass transfer coeff. (in ft/sec) */
- kf = Sh * qu->Diffus / d;
-
- /* For zero-order reaction, return mass transfer coeff. */
- if (qu->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 + ABS(kw)); /* Wall coeff, 1/sec */
- return (kw);
-} /* End of piperate */
-
-double pipereact(EN_Project *pr, int k, double c, double v, long dt)
-/*
-**------------------------------------------------------------
-** Input: k = link index
-** c = current WQ 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 *qu = &pr->quality;
- time_options_t *time = &pr->time_options;
-
- /* For water age (hrs), update concentration by timestep */
- if (qu->Qualflag == AGE)
- return (c + (double)dt / 3600.0);
-
- /* Otherwise find bulk & wall reaction rates */
- rbulk = bulkrate(pr, c, net->Link[k].Kb, qu->BulkOrder) * qu->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 (time->Htime >= time->Rstart) {
- qu->Wbulk += ABS(dcbulk) * v;
- qu->Wwall += ABS(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 WQ 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 *qu = &pr->quality;
- time_options_t *time = &pr->time_options;
-
- /*** Updated 9/7/00 ***/
- /* If no reaction then return current WQ */
- if (!qu->Reactflag) {
- return (c);
- }
-
- /* For water age, update concentration by timestep */
- if (qu->Qualflag == AGE) {
- return (c + (double)dt / 3600.0);
- }
-
- /* Find bulk reaction rate */
- rbulk = bulkrate(pr, c, kb, qu->TankOrder) * qu->Tucf;
-
- /* Find concentration change & update quality */
- dc = rbulk * (double)dt;
- if (time->Htime >= time->Rstart) {
- qu->Wtank += ABS(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 *qu = &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 = qu->Climit + SGN(kb) * c;
- if (ABS(c1) < TINY)
- c1 = SGN(c1) * TINY;
- c = c / c1;
- }
-
- /* N-th order kinetics: */
- else {
- /* Account for limiting potential */
- if (qu->Climit == 0.0)
- c1 = c;
- else
- c1 = MAX(0.0, SGN(kb) * (qu->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 *qu = &pr->quality;
- if (kw == 0.0 || d == 0.0) {
- return (0.0);
- }
- if (qu->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 (ABS(kf) < ABS(kw)) { /* Reaction mass transfer limited */
- kw = kf;
+ {
+ qual->massbalance.final = findstoredmass(pr);
+ massin = qual->massbalance.initial + qual->massbalance.inflow;
+ massout = qual->massbalance.outflow + qual->massbalance.final;
+ massreacted = qual->massbalance.reacted;
+ if (massreacted > 0.0) massout += massreacted;
+ else massin -= massreacted;
+ if (massin == 0.0) qual->massbalance.ratio = 1.0;
+ else qual->massbalance.ratio = massout / massin;
}
- return (kw * 4.0 / d); /* Reaction rate (mass/ft3/sec) */
- } else
- return (c * kf); /* 1st-order reaction */
+}
+
+
+double findstoredmass(EN_Project *pr)
+/*
+**--------------------------------------------------------------
+** Input: none
+** Output: returns total constituent mass stored in the network
+** Purpose: finds the current mass of a constituent stored in
+** all pipes and tanks.
+**--------------------------------------------------------------
+*/
+{
+ int i, k;
+ double totalmass = 0.0;
+ Pseg seg;
+
+ EN_Network *net = &pr->network;
+ quality_t *qual = &pr->quality;
+
+ // Mass residing in each pipe
+ for (k = 1; k <= net->Nlinks; k++)
+ {
+ // Sum up the quality and volume in each segment of the link
+ seg = qual->FirstSeg[k];
+ while (seg != NULL)
+ {
+ totalmass += (seg->c) * (seg->v);
+ seg = seg->prev;
+ }
+ }
+
+ // Mass residing in each tank
+ for (i = 1; i <= net->Ntanks; i++)
+ {
+ // ... skip reservoirs
+ if (net->Tank[i].A == 0.0) continue;
+
+ // ... add up mass in each volume segment
+ else
+ {
+ k = net->Nlinks + i;
+ seg = qual->FirstSeg[k];
+ while (seg != NULL)
+ {
+ totalmass += seg->c * seg->v;
+ seg = seg->prev;
+ }
+ }
+ }
+ return totalmass;
+}
+
+int flowdirchanged(EN_Project *pr)
+/*
+**--------------------------------------------------------------
+** Input: none
+** Output: returns TRUE if flow direction changes in any link
+** Purpose: finds new flow directions for each network link.
+**--------------------------------------------------------------
+*/
+{
+ int k;
+ int result = FALSE;
+ int newdir;
+ int olddir;
+ double q;
+
+ hydraulics_t *hyd = &pr->hydraulics;
+ quality_t *qual = &pr->quality;
+
+ // Examine each network link
+ for (k = 1; k <= pr->network.Nlinks; k++)
+ {
+ // Determine sign (+1 or -1) of new flow rate
+ olddir = qual->FlowDir[k];
+ q = (hyd->LinkStatus[k] <= CLOSED) ? 0.0 : hyd->LinkFlows[k];
+ newdir = SGN(q);
+
+ // Indicate if flow is negligible
+ if (fabs(q) < QZERO) newdir = 0;
+
+ // Reverse link's volume segments if flow direction changes sign
+ if (newdir * olddir < 0) reversesegs(pr, k);
+
+ // If flow direction changes either sign or magnitude then set
+ // result to true (e.g., if a link's positive flow becomes
+ // negligible then the network still needs to be re-sorted)
+ if (newdir != olddir) result = TRUE;
+
+ // ... replace old flow direction with the new direction
+ qual->FlowDir[k] = newdir;
+ }
+ return result;
}
/************************* End of QUALITY.C ***************************/
diff --git a/src/qualreact.c b/src/qualreact.c
new file mode 100644
index 0000000..0778ed4
--- /dev/null
+++ b/src/qualreact.c
@@ -0,0 +1,721 @@
+/*
+*********************************************************************
+
+QUALREACT.C -- water quality reaction module for the EPANET program
+
+*********************************************************************
+*/
+
+#include
+#include
+#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);
+ }
+}
diff --git a/src/qualroute.c b/src/qualroute.c
new file mode 100644
index 0000000..9099f47
--- /dev/null
+++ b/src/qualroute.c
@@ -0,0 +1,776 @@
+/*
+*********************************************************************
+
+QUALROUTE.C -- water quality routing module for the EPANET program
+
+*********************************************************************
+*/
+
+#include
+#include
+#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;
+}
diff --git a/src/report.c b/src/report.c
index 80eb45a..058ab22 100644
--- a/src/report.c
+++ b/src/report.c
@@ -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)
/*
**-------------------------------------------------------------
diff --git a/src/types.h b/src/types.h
index 278d98c..c2350d0 100755
--- a/src/types.h
+++ b/src/types.h
@@ -23,6 +23,7 @@ AUTHOR: L. Rossman
#include "epanet2.h"
#include "hash.h"
#include "util/errormanager.h"
+#include
/*********************************************************/
@@ -34,7 +35,7 @@ AUTHOR: L. Rossman
-------------------------------------------
*/
typedef float REAL4;
-typedef int INT4;
+typedef int INT4;
/*
-----------------------------
@@ -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 {
@@ -361,7 +363,7 @@ struct Sdemand /* DEMAND CATEGORY OBJECT */
{
double Base; /* Baseline demand */
int Pat; /* Pattern index */
- char Name[MAXMSG+1]; /* Demand category name */
+ char Name[MAXMSG+1]; /* Demand category name */
struct Sdemand *next; /* Next record */
};
typedef struct Sdemand *Pdemand; /* Pointer to demand object */
@@ -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 {
diff --git a/src/vars.h b/src/vars.h
deleted file mode 100755
index 9c4c08f..0000000
--- a/src/vars.h
+++ /dev/null
@@ -1,16 +0,0 @@
-/*
-************************************************************************
- Global Variables for EPANET Program
-************************************************************************
-*/
-#ifndef VARS_H
-#define VARS_H
-
-#include
-#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
diff --git a/win_build/WinSDK/Makefile.bat b/win_build/WinSDK/Makefile.bat
index 91ffc95..a175d7e 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 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