diff --git a/include/epanet2.h b/include/epanet2.h index b216597..0543a15 100755 --- a/include/epanet2.h +++ b/include/epanet2.h @@ -9,8 +9,8 @@ 3/1/01 8/15/07 (2.00.11) 2/14/08 (2.00.12) - AUTHOR: L. Rossman - US EPA - NRMRL + AUTHORS: L. Rossman - US EPA - NRMRL + OpenWaterAnalytics members: see git stats for contributors ******************************************************************* */ @@ -248,7 +248,7 @@ extern "C" { int DLLEXPORT ENsetpattern(int index, EN_API_FLOAT_TYPE *f, int len); int DLLEXPORT ENsetpatternvalue(int index, int period, EN_API_FLOAT_TYPE value); int DLLEXPORT ENsettimeparam(int code, long value); - int DLLEXPORT ENsetoption(int, EN_API_FLOAT_TYPE); + int DLLEXPORT ENsetoption(int code, EN_API_FLOAT_TYPE v); int DLLEXPORT ENsetstatusreport(int code); int DLLEXPORT ENsetqualtype(int qualcode, char *chemname, char *chemunits, char *tracenode); int DLLEXPORT ENsetbasedemand(int nodeIndex, int demandIdx, EN_API_FLOAT_TYPE baseDemand); diff --git a/src/epanet.c b/src/epanet.c index 19f668e..2c4aa09 100755 --- a/src/epanet.c +++ b/src/epanet.c @@ -1443,19 +1443,19 @@ int DLLEXPORT ENgetnodevalue(int index, int code, EN_API_FLOAT_TYPE *value) break; //(2.00.11 - LR) case EN_DEMAND: - v = D[index]*Ucf[FLOW]; + v = NodeDemand[index]*Ucf[FLOW]; break; case EN_HEAD: - v = H[index]*Ucf[HEAD]; + v = NodeHead[index]*Ucf[HEAD]; break; case EN_PRESSURE: - v = (H[index] - Node[index].El)*Ucf[PRESSURE]; + v = (NodeHead[index] - Node[index].El)*Ucf[PRESSURE]; break; case EN_QUALITY: - v = C[index]*Ucf[QUALITY]; + v = NodeQual[index]*Ucf[QUALITY]; break; /*** New parameters added for retrieval begins here ***/ //(2.00.12 - LR) @@ -1518,7 +1518,7 @@ int DLLEXPORT ENgetnodevalue(int index, int code, EN_API_FLOAT_TYPE *value) case EN_TANKVOLUME: if (index <= Njuncs) return(251); - v = tankvolume(index-Njuncs, H[index])*Ucf[VOLUME]; + v = tankvolume(index-Njuncs, NodeHead[index])*Ucf[VOLUME]; break; default: return(251); @@ -1685,7 +1685,7 @@ int DLLEXPORT ENgetlinkvalue(int index, int code, EN_API_FLOAT_TYPE *value) case EN_FLOW: /*** Updated 10/25/00 ***/ - if (S[index] <= CLOSED) v = 0.0; + if (LinkStatus[index] <= CLOSED) v = 0.0; else v = Q[index]*Ucf[FLOW]; break; @@ -1693,7 +1693,7 @@ int DLLEXPORT ENgetlinkvalue(int index, int code, EN_API_FLOAT_TYPE *value) if (Link[index].Type == PUMP) v = 0.0; /*** Updated 11/19/01 ***/ - else if (S[index] <= CLOSED) v = 0.0; + else if (LinkStatus[index] <= CLOSED) v = 0.0; else { @@ -1706,26 +1706,31 @@ int DLLEXPORT ENgetlinkvalue(int index, int code, EN_API_FLOAT_TYPE *value) case EN_HEADLOSS: /*** Updated 11/19/01 ***/ - if (S[index] <= CLOSED) v = 0.0; + if (LinkStatus[index] <= CLOSED) v = 0.0; else { - h = H[Link[index].N1] - H[Link[index].N2]; + h = NodeHead[Link[index].N1] - NodeHead[Link[index].N2]; if (Link[index].Type != PUMP) h = ABS(h); v = h*Ucf[HEADLOSS]; } break; case EN_STATUS: - if (S[index] <= CLOSED) v = 0.0; + if (LinkStatus[index] <= CLOSED) v = 0.0; else v = 1.0; break; case EN_SETTING: - if (Link[index].Type == PIPE || Link[index].Type == CV) + if (Link[index].Type == PIPE || Link[index].Type == CV) { return(ENgetlinkvalue(index, EN_ROUGHNESS, value)); - if (K[index] == MISSING) v = 0.0; - else v = K[index]; + } + if (LinkSetting[index] == MISSING) { + v = 0.0; + } + else { + v = LinkSetting[index]; + } switch (Link[index].Type) { case PRV: @@ -1911,7 +1916,7 @@ int DLLEXPORT ENsetnodevalue(int index, int code, EN_API_FLOAT_TYPE v) Tank[j].Hmin += value; Tank[j].Hmax += value; Node[index].El += value; - H[index] += value; + NodeHead[index] += value; } break; @@ -1994,7 +1999,7 @@ int DLLEXPORT ENsetnodevalue(int index, int code, EN_API_FLOAT_TYPE v) Tank[j].Hmin = Tank[j].H0; Tank[j].Hmax = Tank[j].H0; Node[index].El = Tank[j].H0; - H[index] = Tank[j].H0; + NodeHead[index] = Tank[j].H0; } else { @@ -2005,7 +2010,7 @@ int DLLEXPORT ENsetnodevalue(int index, int code, EN_API_FLOAT_TYPE v) Tank[j].V0 = tankvolume(j, Tank[j].H0); // Resetting Volume in addition to initial volume Tank[j].V = Tank[j].V0; - H[index] = Tank[j].H0; + NodeHead[index] = Tank[j].H0; } break; @@ -2160,7 +2165,7 @@ int DLLEXPORT ENsetlinkvalue(int index, int code, EN_API_FLOAT_TYPE v) if (code == EN_INITSTATUS) setlinkstatus(index, s, &Link[index].Stat, &Link[index].Kc); else - setlinkstatus(index, s, &S[index], &K[index]); + setlinkstatus(index, s, &LinkStatus[index], &LinkSetting[index]); break; case EN_INITSETTING: @@ -2187,7 +2192,7 @@ int DLLEXPORT ENsetlinkvalue(int index, int code, EN_API_FLOAT_TYPE v) if (code == EN_INITSETTING) setlinksetting(index, value, &Link[index].Stat, &Link[index].Kc); else - setlinksetting(index, value, &S[index], &K[index]); + setlinksetting(index, value, &LinkStatus[index], &LinkSetting[index]); } break; @@ -2755,13 +2760,13 @@ void initpointers() **---------------------------------------------------------------- */ { - D = NULL; - C = NULL; - H = NULL; + NodeDemand = NULL; + NodeQual = NULL; + NodeHead = NULL; Q = NULL; PipeRateCoeff = NULL; - S = NULL; - K = NULL; + LinkStatus = NULL; + LinkSetting = NULL; OldStat = NULL; Node = NULL; @@ -2824,13 +2829,13 @@ int allocdata() { n = MaxNodes + 1; Node = (Snode *) calloc(n, sizeof(Snode)); - D = (double *) calloc(n, sizeof(double)); - C = (double *) calloc(n, sizeof(double)); - H = (double *) calloc(n, sizeof(double)); + NodeDemand = (double *) calloc(n, sizeof(double)); + NodeQual = (double *) calloc(n, sizeof(double)); + NodeHead = (double *) calloc(n, sizeof(double)); ERRCODE(MEMCHECK(Node)); - ERRCODE(MEMCHECK(D)); - ERRCODE(MEMCHECK(C)); - ERRCODE(MEMCHECK(H)); + ERRCODE(MEMCHECK(NodeDemand)); + ERRCODE(MEMCHECK(NodeQual)); + ERRCODE(MEMCHECK(NodeHead)); } /* Allocate memory for network links */ @@ -2839,12 +2844,12 @@ int allocdata() n = MaxLinks + 1; Link = (Slink *) calloc(n, sizeof(Slink)); Q = (double *) calloc(n, sizeof(double)); - K = (double *) calloc(n, sizeof(double)); - S = (char *) calloc(n, sizeof(char)); + LinkSetting = (double *) calloc(n, sizeof(double)); + LinkStatus = (char *) calloc(n, sizeof(char)); ERRCODE(MEMCHECK(Link)); ERRCODE(MEMCHECK(Q)); - ERRCODE(MEMCHECK(K)); - ERRCODE(MEMCHECK(S)); + ERRCODE(MEMCHECK(LinkSetting)); + ERRCODE(MEMCHECK(LinkStatus)); } /* Allocate memory for tanks, sources, pumps, valves, */ @@ -2955,12 +2960,12 @@ void freedata() Psource source; /* Free memory for computed results */ - free(D); - free(C); - free(H); + free(NodeDemand); + free(NodeQual); + free(NodeHead); free(Q); - free(K); - free(S); + free(LinkSetting); + free(LinkStatus); /* Free memory for node data */ if (Node != NULL) diff --git a/src/hydraul.c b/src/hydraul.c index 8ee6ff8..e9bcbb7 100755 --- a/src/hydraul.c +++ b/src/hydraul.c @@ -123,10 +123,10 @@ void inithyd(int initflag) for (i=1; i<=Ntanks; i++) { Tank[i].V = Tank[i].V0; - H[Tank[i].Node] = Tank[i].H0; + NodeHead[Tank[i].Node] = Tank[i].H0; /*** Updated 10/25/00 ***/ - D[Tank[i].Node] = 0.0; + NodeDemand[Tank[i].Node] = 0.0; OldStat[Nlinks+i] = TEMPCLOSED; } @@ -140,24 +140,24 @@ void inithyd(int initflag) for (i=1; i<=Nlinks; i++) { /* Initialize status and setting */ - S[i] = Link[i].Stat; - K[i] = Link[i].Kc; + LinkStatus[i] = Link[i].Stat; + LinkSetting[i] = Link[i].Kc; /* Start active control valves in ACTIVE position */ //(2.00.11 - LR) if ( (Link[i].Type == PRV || Link[i].Type == PSV || Link[i].Type == FCV) //(2.00.11 - LR) && (Link[i].Kc != MISSING) - ) S[i] = ACTIVE; //(2.00.11 - LR) + ) LinkStatus[i] = ACTIVE; //(2.00.11 - LR) /*** Updated 3/1/01 ***/ /* Initialize flows if necessary */ - if (S[i] <= CLOSED) Q[i] = QZERO; + if (LinkStatus[i] <= CLOSED) Q[i] = QZERO; else if (ABS(Q[i]) <= QZERO || initflag > 0) - initlinkflow(i, S[i], K[i]); + initlinkflow(i, LinkStatus[i], LinkSetting[i]); /* Save initial status */ - OldStat[i] = S[i]; + OldStat[i] = LinkStatus[i]; } /* Reset pump energy usage */ @@ -375,7 +375,7 @@ void setlinkflow(int k, double dh) /* use approx. inverse of formula. */ if (Formflag == DW) { - x = -log(K[k]/3.7/Link[k].Diam); + x = -log(LinkSetting[k]/3.7/Link[k].Diam); y = sqrt(ABS(dh)/Link[k].R/1.32547); Q[k] = x*y; } @@ -402,17 +402,17 @@ void setlinkflow(int k, double dh) /* For custom pump curve, interpolate from curve */ if (Pump[p].Ptype == CUSTOM) { - dh = -dh*Ucf[HEAD]/SQR(K[k]); + dh = -dh*Ucf[HEAD]/SQR(LinkSetting[k]); i = Pump[p].Hcurve; Q[k] = interp(Curve[i].Npts,Curve[i].Y,Curve[i].X, - dh)*K[k]/Ucf[FLOW]; + dh)*LinkSetting[k]/Ucf[FLOW]; } /* Otherwise use inverse of power curve */ else { - h0 = -SQR(K[k])*Pump[p].H0; - x = pow(K[k],2.0-Pump[p].N); + h0 = -SQR(LinkSetting[k])*Pump[p].H0; + x = pow(LinkSetting[k],2.0-Pump[p].N); x = ABS(h0-dh)/(Pump[p].R*x), y = 1.0/Pump[p].N; Q[k] = pow(x,y); @@ -606,7 +606,7 @@ void demands() if (djunc > 0.0) Dsystem += djunc; sum += djunc; } - D[i] = sum; + NodeDemand[i] = sum; } /* Update head at fixed grade nodes with time patterns. */ @@ -619,7 +619,7 @@ void demands() { k = p % (long) Pattern[j].Length; i = Tank[n].Node; - H[i] = Node[i].El*Pattern[j].F[k]; + NodeHead[i] = Node[i].El*Pattern[j].F[k]; } } } @@ -632,7 +632,7 @@ void demands() { i = Pump[n].Link; k = p % (long) Pattern[j].Length; - setlinksetting(i, Pattern[j].F[k], &S[i], &K[i]); + setlinksetting(i, Pattern[j].F[k], &LinkStatus[i], &LinkSetting[i]); } } } /* End of demands */ @@ -664,8 +664,8 @@ int controls() /* Link is controlled by tank level */ if ((n = Control[i].Node) > 0 && n > Njuncs) { - h = H[n]; - vplus = ABS(D[n]); + h = NodeHead[n]; + vplus = ABS(NodeDemand[n]); v1 = tankvolume(n-Njuncs,h); v2 = tankvolume(n-Njuncs,Control[i].Grade); if (Control[i].Type == LOWLEVEL && v1 <= v2 + vplus) @@ -689,16 +689,16 @@ int controls() /* Update link status & pump speed or valve setting */ if (reset == 1) { - if (S[k] <= CLOSED) s1 = CLOSED; + if (LinkStatus[k] <= CLOSED) s1 = CLOSED; else s1 = OPEN; s2 = Control[i].Status; - k1 = K[k]; + k1 = LinkSetting[k]; k2 = k1; if (Link[k].Type > PIPE) k2 = Control[i].Setting; if (s1 != s2 || k1 != k2) { - S[k] = s2; - K[k] = k2; + LinkStatus[k] = s2; + LinkSetting[k] = k2; if (Statflag) writecontrolaction(k,i); // if (s1 != s2) initlinkflow(k, S[k], K[k]); setsum++; @@ -764,8 +764,8 @@ void tanktimestep(long *tstep) { if (Tank[i].A == 0.0) continue; /* Skip reservoirs */ n = Tank[i].Node; - h = H[n]; /* Current tank grade */ - q = D[n]; /* Flow into tank */ + h = NodeHead[n]; /* Current tank grade */ + q = NodeDemand[n]; /* Flow into tank */ if (ABS(q) <= QZERO) continue; if (q > 0.0 && h < Tank[i].Hmax) { @@ -802,8 +802,8 @@ void controltimestep(long *tstep) if ( (n = Control[i].Node) > 0) /* Node control: */ { if ((j = n-Njuncs) <= 0) continue; /* Node is a tank */ - h = H[n]; /* Current tank grade */ - q = D[n]; /* Flow into tank */ + h = NodeHead[n]; /* Current tank grade */ + q = NodeDemand[n]; /* Flow into tank */ if (ABS(q) <= QZERO) continue; if ( (h < Control[i].Grade && @@ -838,8 +838,8 @@ void controltimestep(long *tstep) /* Check if rule actually changes link status or setting */ k = Control[i].Link; if ( - (Link[k].Type > PIPE && K[k] != Control[i].Setting) || - (S[k] != Control[i].Status) + (Link[k].Type > PIPE && LinkSetting[k] != Control[i].Setting) || + (LinkStatus[k] != Control[i].Status) ) *tstep = t; } @@ -954,7 +954,7 @@ void addenergy(long hstep) { /* Skip closed pumps */ k = Pump[j].Link; - if (S[k] <= CLOSED) continue; + if (LinkStatus[k] <= CLOSED) continue; q = MAX(QZERO, ABS(Q[k])); /* Find pump-specific energy cost */ @@ -1001,7 +1001,7 @@ void getenergy(int k, double *kw, double *eff) /*** Updated 6/24/02 ***/ /* No energy if link is closed */ - if (S[k] <= CLOSED) + if (LinkStatus[k] <= CLOSED) { *kw = 0.0; *eff = 0.0; @@ -1011,7 +1011,7 @@ void getenergy(int k, double *kw, double *eff) /* Determine flow and head difference */ q = ABS(Q[k]); - dh = ABS(H[Link[k].N1] - H[Link[k].N2]); + dh = ABS(NodeHead[Link[k].N1] - NodeHead[Link[k].N2]); /* For pumps, find effic. at current flow */ if (Link[k].Type == PUMP) @@ -1053,18 +1053,18 @@ void tanklevels(long tstep) /* Update the tank's volume & water elevation */ n = Tank[i].Node; - dv = D[n]*tstep; + dv = NodeDemand[n]*tstep; Tank[i].V += dv; /*** Updated 6/24/02 ***/ /* Check if tank full/empty within next second */ - if (Tank[i].V + D[n] >= Tank[i].Vmax) { + if (Tank[i].V + NodeDemand[n] >= Tank[i].Vmax) { Tank[i].V = Tank[i].Vmax; } - else if (Tank[i].V - D[n] <= Tank[i].Vmin) { + else if (Tank[i].V - NodeDemand[n] <= Tank[i].Vmin) { Tank[i].V = Tank[i].Vmin; } - H[n] = tankgrade(i,Tank[i].V); + NodeHead[n] = tankgrade(i,Tank[i].V); } } /* End of tanklevels */ @@ -1188,7 +1188,7 @@ int netsolve(int *iter, double *relerr) /* Update current solution. */ /* (Row[i] = row of solution matrix corresponding to node i). */ - for (i=1; i<=Njuncs; i++) H[i] = F[Row[i]]; /* Update heads */ + for (i=1; i<=Njuncs; i++) NodeHead[i] = F[Row[i]]; /* Update heads */ newerr = newflows(); /* Update flows */ *relerr = newerr; @@ -1244,7 +1244,7 @@ int netsolve(int *iter, double *relerr) } /* Add any emitter flows to junction demands */ - for (i=1; i<=Njuncs; i++) D[i] += E[i]; + for (i=1; i<=Njuncs; i++) NodeDemand[i] += E[i]; return(errcode); } /* End of netsolve */ @@ -1274,15 +1274,15 @@ int badvalve(int n) Link[k].Type == PSV || Link[k].Type == FCV) { - if (S[k] == ACTIVE) + if (LinkStatus[k] == ACTIVE) { if (Statflag == FULL) { sprintf(Msg,FMT61,clocktime(Atime,Htime),Link[k].ID); writeline(Msg); } - if (Link[k].Type == FCV) S[k] = XFCV; - else S[k] = XPRESSURE; + if (Link[k].Type == FCV) LinkStatus[k] = XFCV; + else LinkStatus[k] = XPRESSURE; return(1); } } @@ -1313,25 +1313,25 @@ int valvestatus() for (i=1; i<=Nvalves; i++) /* Examine each valve */ { k = Valve[i].Link; /* Link index of valve */ - if (K[k] == MISSING) continue; /* Valve status fixed */ + if (LinkSetting[k] == MISSING) continue; /* Valve status fixed */ n1 = Link[k].N1; /* Start & end nodes */ n2 = Link[k].N2; - s = S[k]; /* Save current status */ + s = LinkStatus[k]; /* Save current status */ // if (s != CLOSED /* No change if flow is */ //(2.00.11 - LR) // && ABS(Q[k]) < Qtol) continue; /* negligible. */ //(2.00.11 - LR) switch (Link[k].Type) /* Evaluate new status: */ { - case PRV: hset = Node[n2].El + K[k]; - S[k] = prvstatus(k,s,hset,H[n1],H[n2]); + case PRV: hset = Node[n2].El + LinkSetting[k]; + LinkStatus[k] = prvstatus(k,s,hset,NodeHead[n1],NodeHead[n2]); break; - case PSV: hset = Node[n1].El + K[k]; - S[k] = psvstatus(k,s,hset,H[n1],H[n2]); + case PSV: hset = Node[n1].El + LinkSetting[k]; + LinkStatus[k] = psvstatus(k,s,hset,NodeHead[n1],NodeHead[n2]); break; //// FCV status checks moved back into the linkstatus() function //// //(2.00.12 - LR) -// case FCV: S[k] = fcvstatus(k,s,H[n1],H[n2]); //(2.00.12 - LR) +// case FCV: S[k] = fcvstatus(k,s,NodeHead[n1],NodeHead[n2]); //(2.00.12 - LR) // break; //(2.00.12 - LR) default: continue; @@ -1342,9 +1342,9 @@ int valvestatus() /* This strategy improves convergence. */ /* Check for status change */ - if (s != S[k]) + if (s != LinkStatus[k]) { - if (Statflag == FULL) writestatchange(k,s,S[k]); + if (Statflag == FULL) writestatchange(k,s,LinkStatus[k]); change = TRUE; } } @@ -1375,29 +1375,29 @@ int linkstatus() { n1 = Link[k].N1; n2 = Link[k].N2; - dh = H[n1] - H[n2]; + dh = NodeHead[n1] - NodeHead[n2]; /* Re-open temporarily closed links (status = XHEAD or TEMPCLOSED) */ - status = S[k]; - if (status == XHEAD || status == TEMPCLOSED) S[k] = OPEN; + status = LinkStatus[k]; + if (status == XHEAD || status == TEMPCLOSED) LinkStatus[k] = OPEN; /* Check for status changes in CVs and pumps */ - if (Link[k].Type == CV) S[k] = cvstatus(S[k],dh,Q[k]); - if (Link[k].Type == PUMP && S[k] >= OPEN && K[k] > 0.0) //(2.00.11 - LR) - S[k] = pumpstatus(k,-dh); + if (Link[k].Type == CV) LinkStatus[k] = cvstatus(LinkStatus[k],dh,Q[k]); + if (Link[k].Type == PUMP && LinkStatus[k] >= OPEN && LinkSetting[k] > 0.0) //(2.00.11 - LR) + LinkStatus[k] = pumpstatus(k,-dh); /* Check for status changes in non-fixed FCVs */ - if (Link[k].Type == FCV && K[k] != MISSING) //(2.00.12 - LR)// - S[k] = fcvstatus(k,status,H[n1],H[n2]); //(2.00.12 - LR)// + if (Link[k].Type == FCV && LinkSetting[k] != MISSING) //(2.00.12 - LR)// + LinkStatus[k] = fcvstatus(k,status,NodeHead[n1],NodeHead[n2]); //(2.00.12 - LR)// /* Check for flow into (out of) full (empty) tanks */ if (n1 > Njuncs || n2 > Njuncs) tankstatus(k,n1,n2); /* Note change in link status; do not revise link flow */ //(2.00.11 - LR) - if (status != S[k]) + if (status != LinkStatus[k]) { change = TRUE; - if (Statflag == FULL) writestatchange(k,status,S[k]); + if (Statflag == FULL) writestatchange(k,status,LinkStatus[k]); //if (S[k] <= CLOSED) Q[k] = QZERO; //(2.00.11 - LR) //else setlinkflow(k, dh); //(2.00.11 - LR) @@ -1449,7 +1449,7 @@ char pumpstatus(int k, double dh) /* Prevent reverse flow through pump */ p = PUMPINDEX(k); if (Pump[p].Ptype == CONST_HP) hmax = BIG; - else hmax = SQR(K[k])*Pump[p].Hmax; + else hmax = SQR(LinkSetting[k])*Pump[p].Hmax; if (dh > hmax + Htol) return(XHEAD); /*** Flow higher than pump curve no longer results in a status change ***/ //(2.00.11 - LR) @@ -1477,7 +1477,7 @@ char prvstatus(int k, char s, double hset, double h1, double h2) double htol = Htol; status = s; - if (K[k] == MISSING) return(status); /* Status fixed by user */ + if (LinkSetting[k] == MISSING) return(status); /* Status fixed by user */ hml = Link[k].Km*SQR(Q[k]); /* Head loss when open */ /*** Status rules below have changed. ***/ //(2.00.11 - LR) @@ -1527,7 +1527,7 @@ char psvstatus(int k, char s, double hset, double h1, double h2) double htol = Htol; status = s; - if (K[k] == MISSING) return(status); /* Status fixed by user */ + if (LinkSetting[k] == MISSING) return(status); /* Status fixed by user */ hml = Link[k].Km*SQR(Q[k]); /* Head loss when open */ /*** Status rules below have changed. ***/ //(2.00.11 - LR) @@ -1580,9 +1580,15 @@ char fcvstatus(int k, char s, double h1, double h2) { char status; /* New valve status */ status = s; - if (h1 - h2 < -Htol) status = XFCV; - else if ( Q[k] < -Qtol ) status = XFCV; //(2.00.11 - LR) - else if (s == XFCV && Q[k] >= K[k]) status = ACTIVE; + if (h1 - h2 < -Htol) { + status = XFCV; + } + else if ( Q[k] < -Qtol ) { + status = XFCV; //(2.00.11 - LR) + } + else if (s == XFCV && Q[k] >= LinkSetting[k]) { + status = ACTIVE; + } return(status); } @@ -1615,39 +1621,39 @@ void tankstatus(int k, int n1, int n2) n2 = n; q = -q; } - h = H[n1] - H[n2]; + h = NodeHead[n1] - NodeHead[n2]; /* Skip reservoirs & closed links */ - if (Tank[i].A == 0.0 || S[k] <= CLOSED) return; + if (Tank[i].A == 0.0 || LinkStatus[k] <= CLOSED) return; /* If tank full, then prevent flow into it */ - if (H[n1] >= Tank[i].Hmax - Htol) + if (NodeHead[n1] >= Tank[i].Hmax - Htol) { /* Case 1: Link is a pump discharging into tank */ if ( Link[k].Type == PUMP ) { - if (Link[k].N2 == n1) S[k] = TEMPCLOSED; + if (Link[k].N2 == n1) LinkStatus[k] = TEMPCLOSED; } /* Case 2: Downstream head > tank head */ /* (i.e., an open outflow check valve would close) */ - else if (cvstatus(OPEN, h, q) == CLOSED) S[k] = TEMPCLOSED; + else if (cvstatus(OPEN, h, q) == CLOSED) LinkStatus[k] = TEMPCLOSED; } /* If tank empty, then prevent flow out of it */ - if (H[n1] <= Tank[i].Hmin + Htol) + if (NodeHead[n1] <= Tank[i].Hmin + Htol) { /* Case 1: Link is a pump discharging from tank */ if ( Link[k].Type == PUMP) { - if (Link[k].N1 == n1) S[k] = TEMPCLOSED; + if (Link[k].N1 == n1) LinkStatus[k] = TEMPCLOSED; } /* Case 2: Tank head > downstream head */ /* (i.e., a closed outflow check valve would open) */ - else if (cvstatus(CLOSED, h, q) == OPEN) S[k] = TEMPCLOSED; + else if (cvstatus(CLOSED, h, q) == OPEN) LinkStatus[k] = TEMPCLOSED; } } /* End of tankstatus */ @@ -1681,10 +1687,10 @@ int pswitch() { /* Determine if control conditions are satisfied */ if (Control[i].Type == LOWLEVEL - && H[n] <= Control[i].Grade + Htol ) + && NodeHead[n] <= Control[i].Grade + Htol ) reset = 1; if (Control[i].Type == HILEVEL - && H[n] >= Control[i].Grade - Htol ) + && NodeHead[n] >= Control[i].Grade - Htol ) reset = 1; } @@ -1692,28 +1698,28 @@ int pswitch() if (reset == 1) { change = 0; - s = S[k]; + s = LinkStatus[k]; if (Link[k].Type == PIPE) { if (s != Control[i].Status) change = 1; } if (Link[k].Type == PUMP) { - if (K[k] != Control[i].Setting) change = 1; + if (LinkSetting[k] != Control[i].Setting) change = 1; } if (Link[k].Type >= PRV) { - if (K[k] != Control[i].Setting) change = 1; - else if (K[k] == MISSING && + if (LinkSetting[k] != Control[i].Setting) change = 1; + else if (LinkSetting[k] == MISSING && s != Control[i].Status) change = 1; } /* If a change occurs, update status & setting */ if (change) { - S[k] = Control[i].Status; - if (Link[k].Type > PIPE) K[k] = Control[i].Setting; - if (Statflag == FULL) writestatchange(k,s,S[k]); + LinkStatus[k] = Control[i].Status; + if (Link[k].Type > PIPE) LinkSetting[k] = Control[i].Setting; + if (Statflag == FULL) writestatchange(k,s,LinkStatus[k]); /* Re-set flow if status has changed */ // if (S[k] != s) initlinkflow(k, S[k], K[k]); @@ -1741,7 +1747,7 @@ double newflows() int k, n, n1, n2; /* Initialize net inflows (i.e., demands) at tanks */ - for (n=Njuncs+1; n <= Nnodes; n++) D[n] = 0.0; + for (n=Njuncs+1; n <= Nnodes; n++) NodeDemand[n] = 0.0; /* Initialize sum of flows & corrections */ qsum = 0.0; @@ -1761,7 +1767,7 @@ double newflows() n1 = Link[k].N1; n2 = Link[k].N2; - dh = H[n1] - H[n2]; + dh = NodeHead[n1] - NodeHead[n2]; dq = Y[k] - P[k]*dh; /* Adjust flow change by the relaxation factor */ //(2.00.11 - LR) @@ -1780,10 +1786,10 @@ double newflows() dqsum += ABS(dq); /* Update net flows to tanks */ - if ( S[k] > CLOSED ) //(2.00.12 - LR) + if ( LinkStatus[k] > CLOSED ) //(2.00.12 - LR) { - if (n1 > Njuncs) D[n1] -= Q[k]; - if (n2 > Njuncs) D[n2] += Q[k]; + if (n1 > Njuncs) NodeDemand[n1] -= Q[k]; + if (n2 > Njuncs) NodeDemand[n2] += Q[k]; } } @@ -1861,7 +1867,7 @@ void linkcoeffs() case PRV: case PSV: /* If valve status fixed then treat as pipe */ /* otherwise ignore the valve for now. */ - if (K[k] == MISSING) valvecoeff(k); //pipecoeff(k); //(2.00.11 - LR) + if (LinkSetting[k] == MISSING) valvecoeff(k); //pipecoeff(k); //(2.00.11 - LR) else continue; break; default: continue; @@ -1877,13 +1883,13 @@ void linkcoeffs() Aii[Row[n1]] += P[k]; /* Diagonal coeff. */ F[Row[n1]] += Y[k]; /* RHS coeff. */ } - else F[Row[n2]] += (P[k]*H[n1]); /* Node n1 is a tank */ + else F[Row[n2]] += (P[k]*NodeHead[n1]); /* Node n1 is a tank */ if (n2 <= Njuncs) /* Node n2 is junction */ { Aii[Row[n2]] += P[k]; /* Diagonal coeff. */ F[Row[n2]] -= Y[k]; /* RHS coeff. */ } - else F[Row[n1]] += (P[k]*H[n2]); /* Node n2 is a tank */ + else F[Row[n1]] += (P[k]*NodeHead[n2]); /* Node n2 is a tank */ } } /* End of linkcoeffs */ @@ -1904,7 +1910,7 @@ void nodecoeffs() /* flow imbalance & add imbalance to RHS array F. */ for (i=1; i<=Njuncs; i++) { - X[i] -= D[i]; + X[i] -= NodeDemand[i]; F[Row[i]] += X[i]; } } /* End of nodecoeffs */ @@ -1925,7 +1931,7 @@ void valvecoeffs() for (i=1; i<=Nvalves; i++) /* Examine each valve */ { k = Valve[i].Link; /* Link index of valve */ - if (K[k] == MISSING) continue; /* Valve status fixed */ + if (LinkSetting[k] == MISSING) continue; /* Valve status fixed */ n1 = Link[k].N1; /* Start & end nodes */ n2 = Link[k].N2; switch (Link[k].Type) /* Call valve-specific */ @@ -1992,7 +1998,7 @@ double emitflowchange(int i) p = 1/RQtol; else p = 1.0/p; - return(E[i]/Qexp - p*(H[i] - Node[i].El)); + return(E[i]/Qexp - p*(NodeHead[i] - Node[i].El)); } @@ -2019,7 +2025,7 @@ void pipecoeff(int k) dfdq; /* Derivative of fric. fact. */ /* For closed pipe use headloss formula: h = CBIG*q */ - if (S[k] <= CLOSED) + if (LinkStatus[k] <= CLOSED) { P[k] = 1.0/CBIG; Y[k] = Q[k]; @@ -2148,8 +2154,10 @@ void pumpcoeff(int k) r, /* Flow resistance coeff. */ n; /* Flow exponent coeff. */ + double setting = LinkSetting[k]; + /* Use high resistance pipe if pump closed or cannot deliver head */ - if (S[k] <= CLOSED || K[k] == 0.0) + if (LinkStatus[k] <= CLOSED || setting == 0.0) { //pipecoeff(k); //(2.00.11 - LR) P[k] = 1.0/CBIG; //(2.00.11 - LR) @@ -2166,7 +2174,7 @@ void pumpcoeff(int k) { /* Find intercept (h0) & slope (r) of pump curve */ /* line segment which contains speed-adjusted flow. */ - curvecoeff(Pump[p].Hcurve,q/K[k],&h0,&r); + curvecoeff(Pump[p].Hcurve, q/setting, &h0, &r); /* Determine head loss coefficients. */ Pump[p].H0 = -h0; @@ -2175,9 +2183,9 @@ void pumpcoeff(int k) } /* Adjust head loss coefficients for pump speed. */ - h0 = SQR(K[k])*Pump[p].H0; + h0 = SQR(setting)*Pump[p].H0; n = Pump[p].N; - r = Pump[p].R*pow(K[k],2.0-n); + r = Pump[p].R*pow(setting,2.0-n); if (n != 1.0) r = n*r*pow(q,n-1.0); /* Compute inverse headloss gradient (P) and flow correction factor (Y) */ @@ -2241,7 +2249,7 @@ void gpvcoeff(int k) /*** Updated 9/7/00 ***/ /* Treat as a pipe if valve closed */ - if (S[k] == CLOSED) valvecoeff(k); //pipecoeff(k); //(2.00.11 - LR) + if (LinkStatus[k] == CLOSED) valvecoeff(k); //pipecoeff(k); //(2.00.11 - LR) /* Otherwise utilize headloss curve */ /* whose index is stored in K */ @@ -2253,7 +2261,7 @@ void gpvcoeff(int k) /*** Updated 10/25/00 ***/ /*** Updated 12/29/00 ***/ - curvecoeff((int)ROUND(K[k]),q,&h0,&r); + curvecoeff((int)ROUND(LinkSetting[k]),q,&h0,&r); /* Compute inverse headloss gradient (P) */ /* and flow correction factor (Y). */ @@ -2273,19 +2281,19 @@ void pbvcoeff(int k) */ { /* If valve fixed OPEN or CLOSED then treat as a pipe */ - if (K[k] == MISSING || K[k] == 0.0) valvecoeff(k); //pipecoeff(k); //(2.00.11 - LR) + if (LinkSetting[k] == MISSING || LinkSetting[k] == 0.0) valvecoeff(k); //pipecoeff(k); //(2.00.11 - LR) /* If valve is active */ else { /* Treat as a pipe if minor loss > valve setting */ - if (Link[k].Km*SQR(Q[k]) > K[k]) valvecoeff(k); //pipecoeff(k); //(2.00.11 - LR) + if (Link[k].Km*SQR(Q[k]) > LinkSetting[k]) valvecoeff(k); //pipecoeff(k); //(2.00.11 - LR) /* Otherwise force headloss across valve to be equal to setting */ else { P[k] = CBIG; - Y[k] = K[k]*CBIG; + Y[k] = LinkSetting[k]*CBIG; } } } /* End of pbvcoeff */ @@ -2306,8 +2314,8 @@ void tcvcoeff(int k) km = Link[k].Km; /* If valve not fixed OPEN or CLOSED, compute its loss coeff. */ - if (K[k] != MISSING) - Link[k].Km = 0.02517*K[k]/(SQR(Link[k].Diam)*SQR(Link[k].Diam)); + if (LinkSetting[k] != MISSING) + Link[k].Km = 0.02517*LinkSetting[k]/(SQR(Link[k].Diam)*SQR(Link[k].Diam)); /* Then apply usual pipe formulas */ valvecoeff(k); //pipecoeff(k); //(2.00.11 - LR) @@ -2333,9 +2341,9 @@ void prvcoeff(int k, int n1, int n2) double hset; /* Valve head setting */ i = Row[n1]; /* Matrix rows of nodes */ j = Row[n2]; - hset = Node[n2].El + K[k]; /* Valve setting */ + hset = Node[n2].El + LinkSetting[k]; /* Valve setting */ - if (S[k] == ACTIVE) + if (LinkStatus[k] == ACTIVE) { /* Set coeffs. to force head at downstream @@ -2379,9 +2387,9 @@ void psvcoeff(int k, int n1, int n2) double hset; /* Valve head setting */ i = Row[n1]; /* Matrix rows of nodes */ j = Row[n2]; - hset = Node[n1].El + K[k]; /* Valve setting */ + hset = Node[n1].El + LinkSetting[k]; /* Valve setting */ - if (S[k] == ACTIVE) + if (LinkStatus[k] == ACTIVE) { /* Set coeffs. to force head at upstream @@ -2423,7 +2431,7 @@ void fcvcoeff(int k, int n1, int n2) { int i,j; /* Rows in solution matrix */ double q; /* Valve flow setting */ - q = K[k]; + q = LinkSetting[k]; i = Row[n1]; j = Row[n2]; @@ -2432,7 +2440,7 @@ void fcvcoeff(int k, int n1, int n2) flow setting as external demand at upstream node and external supply at downstream node. */ - if (S[k] == ACTIVE) + if (LinkStatus[k] == ACTIVE) { X[n1] -= q; F[i] -= q; @@ -2474,7 +2482,7 @@ void valvecoeff(int k) double p; // Valve is closed. Use a very small matrix coeff. - if (S[k] <= CLOSED) + if (LinkStatus[k] <= CLOSED) { P[k] = 1.0/CBIG; Y[k] = Q[k]; diff --git a/src/input3.c b/src/input3.c index a19bfdc..d0922f6 100755 --- a/src/input3.c +++ b/src/input3.c @@ -99,9 +99,9 @@ int juncdata() demand->Pat = p; demand->next = Node[Njuncs].D; Node[Njuncs].D = demand; - D[Njuncs] = y; + NodeDemand[Njuncs] = y; } - else D[Njuncs] = MISSING; + else NodeDemand[Njuncs] = MISSING; /*** end of update ***/ return(0); } /* end of juncdata */ @@ -682,11 +682,11 @@ int demanddata() /*** Updated 6/24/02 ***/ demand = Node[j].D; - if (demand && D[j] != MISSING) + if (demand && NodeDemand[j] != MISSING) { demand->Base = y; demand->Pat = p; - D[j] = MISSING; + NodeDemand[j] = MISSING; } /*** End of update ***/ diff --git a/src/output.c b/src/output.c index 9d04531..0e18bb8 100755 --- a/src/output.c +++ b/src/output.c @@ -152,30 +152,30 @@ int savehyd(long *htime) fwrite(&t,sizeof(INT4),1,HydFile); /* Save current nodal demands (D) */ - for (i=1; i<=Nnodes; i++) x[i] = (REAL4)D[i]; + for (i=1; i<=Nnodes; i++) x[i] = (REAL4)NodeDemand[i]; fwrite(x+1,sizeof(REAL4),Nnodes,HydFile); /* Copy heads (H) to buffer of floats (x) and save buffer */ - for (i=1; i<=Nnodes; i++) x[i] = (REAL4)H[i]; + for (i=1; i<=Nnodes; i++) x[i] = (REAL4)NodeHead[i]; fwrite(x+1,sizeof(REAL4),Nnodes,HydFile); /* Force flow in closed links to be zero then save flows */ for (i=1; i<=Nlinks; i++) { - if (S[i] <= CLOSED) x[i] = 0.0f; + if (LinkStatus[i] <= CLOSED) x[i] = 0.0f; else x[i] = (REAL4)Q[i]; } fwrite(x+1,sizeof(REAL4),Nlinks,HydFile); /* Copy link status to buffer of floats (x) & write buffer */ - for (i=1; i<=Nlinks; i++) x[i] = (REAL4)S[i]; + for (i=1; i<=Nlinks; i++) x[i] = (REAL4)LinkStatus[i]; fwrite(x+1,sizeof(REAL4),Nlinks,HydFile); /* Save link settings & check for successful write-to-disk */ /* (We assume that if any of the previous fwrites failed, */ /* then this one will also fail.) */ - for (i=1; i<=Nlinks; i++) x[i] = (REAL4)K[i]; + for (i=1; i<=Nlinks; i++) x[i] = (REAL4)LinkSetting[i]; if (fwrite(x+1,sizeof(REAL4),Nlinks,HydFile) < (unsigned)Nlinks) errcode = 308; free(x); @@ -286,19 +286,19 @@ int readhyd(long *hydtime) *hydtime = t; if (fread(x+1,sizeof(REAL4),Nnodes,HydFile) < (unsigned)Nnodes) result = 0; - else for (i=1; i<=Nnodes; i++) D[i] = x[i]; + else for (i=1; i<=Nnodes; i++) NodeDemand[i] = x[i]; if (fread(x+1,sizeof(REAL4),Nnodes,HydFile) < (unsigned)Nnodes) result = 0; - else for (i=1; i<=Nnodes; i++) H[i] = x[i]; + else for (i=1; i<=Nnodes; i++) NodeHead[i] = x[i]; if (fread(x+1,sizeof(REAL4),Nlinks,HydFile) < (unsigned)Nlinks) result = 0; else for (i=1; i<=Nlinks; i++) Q[i] = x[i]; if (fread(x+1,sizeof(REAL4),Nlinks,HydFile) < (unsigned)Nlinks) result = 0; - else for (i=1; i<=Nlinks; i++) S[i] = (char) x[i]; + else for (i=1; i<=Nlinks; i++) LinkStatus[i] = (char) x[i]; if (fread(x+1,sizeof(REAL4),Nlinks,HydFile) < (unsigned)Nlinks) result = 0; - else for (i=1; i<=Nlinks; i++) K[i] = x[i]; + else for (i=1; i<=Nlinks; i++) LinkSetting[i] = x[i]; free(x); return result; @@ -361,16 +361,16 @@ int nodeoutput(int j, REAL4 *x, double ucf) switch(j) { case DEMAND: for (i=1; i<=Nnodes; i++) - x[i] = (REAL4)(D[i]*ucf); + x[i] = (REAL4)(NodeDemand[i]*ucf); break; case HEAD: for (i=1; i<=Nnodes; i++) - x[i] = (REAL4)(H[i]*ucf); + x[i] = (REAL4)(NodeHead[i]*ucf); break; case PRESSURE: for (i=1; i<=Nnodes; i++) - x[i] = (REAL4)((H[i] - Node[i].El)*ucf); + x[i] = (REAL4)((NodeHead[i] - Node[i].El)*ucf); break; case QUALITY: for (i=1; i<=Nnodes; i++) - x[i] = (REAL4)(C[i]*ucf); + x[i] = (REAL4)(NodeQual[i]*ucf); } /* Write x[1] to x[Nnodes] to output file */ @@ -413,10 +413,10 @@ int linkoutput(int j, REAL4 *x, double ucf) break; case HEADLOSS: for (i=1; i<=Nlinks; i++) { - if (S[i] <= CLOSED) x[i] = 0.0f; + if (LinkStatus[i] <= CLOSED) x[i] = 0.0f; else { - h = H[Link[i].N1] - H[Link[i].N2]; + h = NodeHead[Link[i].N1] - NodeHead[Link[i].N2]; if (Link[i].Type != PUMP) h = ABS(h); if (Link[i].Type <= PIPE) x[i] = (REAL4)(1000.0*h/Link[i].Len); @@ -428,25 +428,26 @@ int linkoutput(int j, REAL4 *x, double ucf) x[i] = (REAL4)(avgqual(i)*ucf); break; case STATUS: for (i=1; i<=Nlinks; i++) - x[i] = (REAL4)S[i]; + x[i] = (REAL4)LinkStatus[i]; break; case SETTING: for (i=1; i<=Nlinks; i++) { - if (K[i] != MISSING) + double setting = LinkSetting[i]; + if (setting != MISSING) switch (Link[i].Type) { case CV: - case PIPE: x[i] = (REAL4)K[i]; + case PIPE: x[i] = (REAL4)setting; break; - case PUMP: x[i] = (REAL4)K[i]; + case PUMP: x[i] = (REAL4)setting; break; case PRV: case PSV: - case PBV: x[i] = (REAL4)(K[i]*Ucf[PRESSURE]); + case PBV: x[i] = (REAL4)(setting*Ucf[PRESSURE]); break; - case FCV: x[i] = (REAL4)(K[i]*Ucf[FLOW]); + case FCV: x[i] = (REAL4)(setting*Ucf[FLOW]); break; - case TCV: x[i] = (REAL4)K[i]; + case TCV: x[i] = (REAL4)setting; break; default: x[i] = 0.0f; } @@ -464,7 +465,7 @@ int linkoutput(int j, REAL4 *x, double ucf) { if (Link[i].Type <= PIPE && ABS(Q[i]) > TINY) { - h = ABS(H[Link[i].N1] - H[Link[i].N2]); + h = ABS(NodeHead[Link[i].N1] - NodeHead[Link[i].N2]); f = 39.725*h*pow(Link[i].Diam,5)/Link[i].Len/SQR(Q[i]); x[i] = (REAL4)f; } @@ -642,11 +643,11 @@ int savetimestat(REAL4 *x, char objtype) /* Update internal output variables where applicable */ if (objtype == NODEHDR) switch (j) { - case DEMAND: for (i=1; i<=n; i++) D[i] = x[i]/Ucf[DEMAND]; + case DEMAND: for (i=1; i<=n; i++) NodeDemand[i] = x[i]/Ucf[DEMAND]; break; - case HEAD: for (i=1; i<=n; i++) H[i] = x[i]/Ucf[HEAD]; + case HEAD: for (i=1; i<=n; i++) NodeHead[i] = x[i]/Ucf[HEAD]; break; - case QUALITY: for (i=1; i<=n; i++) C[i] = x[i]/Ucf[QUALITY]; + case QUALITY: for (i=1; i<=n; i++) NodeQual[i] = x[i]/Ucf[QUALITY]; break; } else if (j == FLOW) for (i=1; i<=n; i++) Q[i] = x[i]/Ucf[FLOW]; diff --git a/src/quality.c b/src/quality.c index 8a5b8dd..be9193d 100755 --- a/src/quality.c +++ b/src/quality.c @@ -106,9 +106,9 @@ int openqual() if (SegPool == NULL) errcode = 101; //(2.00.11 - LR) /* Allocate scratch array & reaction rate array*/ - XC = (double *) calloc(MAX((Nnodes+1),(Nlinks+1)),sizeof(double)); + TempQual = (double *) calloc(MAX((Nnodes+1),(Nlinks+1)),sizeof(double)); PipeRateCoeff = (double *) calloc((Nlinks+1), sizeof(double)); - ERRCODE(MEMCHECK(XC)); + ERRCODE(MEMCHECK(TempQual)); ERRCODE(MEMCHECK(PipeRateCoeff)); /* Allocate memory for WQ solver */ @@ -148,7 +148,7 @@ void initqual() int i; /* Initialize quality, tank volumes, & source mass flows */ - for (i=1; i<=Nnodes; i++) C[i] = Node[i].C0; + for (i=1; i<=Nnodes; i++) NodeQual[i] = Node[i].C0; for (i=1; i<=Ntanks; i++) Tank[i].C = Node[Tank[i].Node].C0; for (i=1; i<=Ntanks; i++) Tank[i].V = Tank[i].V0; for (i=1; i<=Nnodes; i++) { @@ -165,7 +165,7 @@ void initqual() if (Qualflag != NONE) { /* Initialize WQ at trace node (if applicable) */ - if (Qualflag == TRACE) C[TraceNode] = 100.0; + if (Qualflag == TRACE) NodeQual[TraceNode] = 100.0; /* Compute Schmidt number */ if (Diffus > 0.0) @@ -242,7 +242,7 @@ int runqual(long *t) for (int i=1; i<= Nlinks; ++i) { - if (S[i] <= CLOSED) { + if (LinkStatus[i] <= CLOSED) { QLinkFlow[i-1] = Q[i]; } } @@ -257,7 +257,7 @@ int runqual(long *t) for (int i=1; i<= Nlinks; ++i) { - if (S[i] <= CLOSED) { + if (LinkStatus[i] <= CLOSED) { QLinkFlow[i-1] = Q[i]; } } @@ -306,13 +306,13 @@ int nextqual(long *tstep) if (Tank[i].A != 0) { // skip reservoirs again int n = Tank[i].Node; Tank[i].V = QTankVolumes[i-1]; - H[n] = tankgrade(i,Tank[i].V); + NodeHead[n] = tankgrade(i,Tank[i].V); } } // restore the previous step's pipe link flows for (int i=1; i<=Nlinks; i++) { - if (S[i] <= CLOSED) { + if (LinkStatus[i] <= CLOSED) { Q[i] = 0.0; } } @@ -336,12 +336,12 @@ int nextqual(long *tstep) if (Tank[i].A != 0) { // skip reservoirs again int n = Tank[i].Node; Tank[i].V = tankVolumes[i-1]; - H[n] = tankgrade(i,Tank[i].V); + NodeHead[n] = tankgrade(i,Tank[i].V); } } for (int i=1; i<=Nlinks; ++i) { - if (S[i] <= CLOSED) { + if (LinkStatus[i] <= CLOSED) { Q[i] = QLinkFlow[i-1]; } } @@ -415,7 +415,7 @@ int closequal() free(VolIn); free(MassIn); free(PipeRateCoeff); - free(XC); + free(TempQual); free(QTankVolumes); free(QLinkFlow); return(errcode); @@ -571,7 +571,7 @@ void initsegs() /* Find quality of downstream node */ j = DOWN_NODE(k); - if (j <= Njuncs) c = C[j]; + if (j <= Njuncs) c = NodeQual[j]; else c = Tank[j-Njuncs].C; /* Fill link with single segment with this quality */ @@ -781,7 +781,7 @@ void accumulate(long dt) /* Re-set memory used to accumulate mass & volume */ memset(VolIn,0,(Nnodes+1)*sizeof(double)); memset(MassIn,0,(Nnodes+1)*sizeof(double)); - memset(XC,0,(Nnodes+1)*sizeof(double)); + memset(TempQual,0,(Nnodes+1)*sizeof(double)); /* Compute average conc. of segments adjacent to each node */ /* (For use if there is no transport through the node) */ @@ -803,7 +803,7 @@ void accumulate(long dt) for (k=1; k<=Nnodes; k++) { if (VolIn[k] > 0.0) { - XC[k] = MassIn[k]/VolIn[k]; + TempQual[k] = MassIn[k]/VolIn[k]; } } @@ -824,7 +824,7 @@ void accumulate(long dt) { VolIn[j] += v; seg = FirstSeg[k]; - cseg = C[i]; + cseg = NodeQuali]; if (seg != NULL) cseg = seg->c; MassIn[j] += v*cseg; removesegs(k); @@ -886,7 +886,7 @@ void updatenodes(long dt) ** Purpose: updates concentration at all nodes to mixture of accumulated ** inflow from connecting pipes. ** -** Note: Does not account for source flow effects. XC[i] contains +** 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. **--------------------------------------------------------------------------- @@ -897,14 +897,14 @@ void updatenodes(long dt) /* Update junction quality */ for (i=1; i<=Njuncs; i++) { - if (D[i] < 0.0) { - VolIn[i] -= D[i]*dt; + if (NodeDemand[i] < 0.0) { + VolIn[i] -= NodeDemand[i]*dt; } if (VolIn[i] > 0.0) { - C[i] = MassIn[i]/VolIn[i]; + NodeQual[i] = MassIn[i]/VolIn[i]; } else { - C[i] = XC[i]; + NodeQual[i] = TempQual[i]; } } @@ -912,7 +912,7 @@ void updatenodes(long dt) updatetanks(dt); /* For flow tracing, set source node concen. to 100. */ - if (Qualflag == TRACE) C[TraceNode] = 100.0; + if (Qualflag == TRACE) NodeQual[TraceNode] = 100.0; } @@ -934,14 +934,14 @@ void sourceinput(long dt) /* Establish a flow cutoff which indicates no outflow from a node */ qcutoff = 10.0*TINY; - /* Zero-out the work array XC */ - memset(XC,0,(Nnodes+1)*sizeof(double)); + /* Zero-out the work array TempQual */ + memset(TempQual,0,(Nnodes+1)*sizeof(double)); if (Qualflag != CHEM) return; /* Consider each node */ for (n=1; n<=Nnodes; n++) { - + double thisDemand = NodeDemand[n]; /* Skip node if no WQ source */ source = Node[n].S; if (source == NULL) continue; @@ -949,7 +949,7 @@ void sourceinput(long dt) /* Find total flow volume leaving node */ if (n <= Njuncs) volout = VolIn[n]; /* Junctions */ - else volout = VolIn[n] - D[n]*dt; /* Tanks */ + else volout = VolIn[n] - (thisDemand * dt); /* Tanks */ qout = volout / (double) dt; /* Evaluate source input only if node outflow > cutoff flow */ @@ -965,13 +965,13 @@ void sourceinput(long dt) case CONCEN: /* Only add source mass if demand is negative */ - if (D[n] < 0.0) + if (thisDemand < 0.0) { - massadded = -s*D[n]*dt; + 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 > Njuncs) C[n] = 0.0; + if (n > Njuncs) NodeQual[n] = 0.0; } else massadded = 0.0; break; @@ -985,8 +985,8 @@ void sourceinput(long dt) /* Mass added is difference between source */ /* & node concen. times outflow volume */ case SETPOINT: - if (s > C[n]) { - massadded = (s-C[n])*volout; + if (s > NodeQual[n]) { + massadded = (s-NodeQual[n])*volout; } else { massadded = 0.0; @@ -1001,7 +1001,7 @@ void sourceinput(long dt) } /* Source concen. contribution = (mass added / outflow volume) */ - XC[n] = massadded/volout; + TempQual[n] = massadded/volout; /* Update total mass added for time period & simulation */ source->Smass += massadded; @@ -1017,8 +1017,8 @@ void sourceinput(long dt) if (Tank[j].A == 0.0) { n = Njuncs + j; - volout = VolIn[n] - D[n]*dt; - if (volout > 0.0) Wsource += volout*C[n]; + volout = VolIn[n] - NodeDemand[n]*dt; + if (volout > 0.0) Wsource += volout*NodeQual[n]; } } } @@ -1052,7 +1052,7 @@ void release(long dt) v = q*dt; /* Include source contribution in quality released from node. */ - c = C[n] + XC[n]; + c = NodeQual[n] + TempQual[n]; /* If link has a last seg, check if its quality */ /* differs from that of the flow released from node.*/ @@ -1081,7 +1081,7 @@ void updatesourcenodes(long dt) ** Input: dt = current WQ time step ** Output: none ** Purpose: updates quality at source nodes. -** (XC[n] = concen. added by source at node n) +** (TempQual[n] = concen. added by source at node n) **--------------------------------------------------- */ { @@ -1097,13 +1097,13 @@ void updatesourcenodes(long dt) if (source == NULL) continue; /* Add source to current node concen. */ - C[n] += XC[n]; + NodeQual[n] += TempQual[n]; /* For tanks, node concen. = internal concen. */ if (n > Njuncs) { i = n - Njuncs; - if (Tank[i].A > 0.0) C[n] = Tank[i].C; + if (Tank[i].A > 0.0) NodeQual[n] = Tank[i].C; } /* Normalize mass added at source to time step */ @@ -1130,7 +1130,7 @@ void updatetanks(long dt) /* Use initial quality for reservoirs */ if (Tank[i].A == 0.0) { - C[n] = Node[n].C0; + NodeQual[n] = Node[n].C0; } /* Update tank WQ based on mixing model */ else { @@ -1173,7 +1173,7 @@ void updatetanks(long dt) // /* Update tank volume & nodal quality */ // Tank[i].V += D[n]*dt; -// C[n] = Tank[i].C; +// NodeQual[n] = Tank[i].C; //} @@ -1198,7 +1198,7 @@ void tankmix1(int i, long dt) /* Determine tank & volumes */ vold = Tank[i].V; n = Tank[i].Node; - Tank[i].V += D[n]*dt; + Tank[i].V += NodeDemand[n]*dt; vin = VolIn[n]; /* Compute inflow concen. */ @@ -1211,7 +1211,7 @@ void tankmix1(int i, long dt) c = MIN(c, cmax); c = MAX(c, 0.0); Tank[i].C = c; - C[n] = Tank[i].C; + NodeQual[n] = Tank[i].C; } /*** Updated 10/25/00 ***/ @@ -1248,7 +1248,7 @@ void tankmix2(int i, long dt) /* Find inflows & outflows */ n = Tank[i].Node; - vnet = D[n]*dt; + vnet = NodeDemand[n]*dt; vin = VolIn[n]; if (vin > 0.0) cin = MassIn[n]/vin; else cin = 0.0; @@ -1304,7 +1304,7 @@ void tankmix2(int i, long dt) /* represent quality of tank since this is where */ /* outflow begins to flow from */ Tank[i].C = seg1->c; - C[n] = Tank[i].C; + NodeQual[n] = Tank[i].C; } @@ -1339,7 +1339,7 @@ void tankmix3(int i, long dt) /* Find inflows & outflows */ n = Tank[i].Node; - vnet = D[n]*dt; + vnet = NodeDemand[n]*dt; vin = VolIn[n]; vout = vin - vnet; if (vin > 0.0) cin = MassIn[n]/VolIn[n]; @@ -1380,7 +1380,7 @@ void tankmix3(int i, long dt) /* to represent overall quality of tank */ if (vsum > 0.0) Tank[i].C = csum/vsum; else Tank[i].C = FirstSeg[k]->c; - C[n] = Tank[i].C; + NodeQual[n] = Tank[i].C; /* Add new last segment for new flow entering tank */ if (vin > 0.0) @@ -1430,7 +1430,7 @@ void tankmix4(int i, long dt) /* Find inflows & outflows */ n = Tank[i].Node; - vnet = D[n]*dt; + vnet = NodeDemand[n]*dt; vin = VolIn[n]; if (vin > 0.0) cin = MassIn[n]/VolIn[n]; else cin = 0.0; @@ -1498,7 +1498,7 @@ void tankmix4(int i, long dt) /* Reported tank quality is mixture of flow released and any inflow */ Tank[i].C = (csum + MassIn[n])/(vsum + vin); } - C[n] = Tank[i].C; + NodeQual[n] = Tank[i].C; } @@ -1553,7 +1553,7 @@ double avgqual(int k) seg = seg->prev; } if (vsum > 0.0) return(msum/vsum); - else return( (C[Link[k].N1] + C[Link[k].N2])/2. ); + else return( (NodeQual[Link[k].N1] + NodeQual[Link[k].N2])/2. ); } diff --git a/src/report.c b/src/report.c index b021461..deffe03 100755 --- a/src/report.c +++ b/src/report.c @@ -296,15 +296,15 @@ void writehydstat(int iter, double relerr) for (i=1; i<=Ntanks; i++) { n = Tank[i].Node; - if (ABS(D[n]) < 0.001) newstat = CLOSED; - else if (D[n] > 0.0) newstat = FILLING; - else if (D[n] < 0.0) newstat = EMPTYING; + if (ABS(NodeDemand[n]) < 0.001) newstat = CLOSED; + else if (NodeDemand[n] > 0.0) newstat = FILLING; + else if (NodeDemand[n] < 0.0) newstat = EMPTYING; else newstat = OldStat[Nlinks+i]; if (newstat != OldStat[Nlinks+i]) { if (Tank[i].A > 0.0) sprintf(s1,FMT50,atime,Node[n].ID,StatTxt[newstat], - (H[n]-Node[n].El)*Ucf[HEAD],Field[HEAD].Units); + (NodeHead[n]-Node[n].El)*Ucf[HEAD],Field[HEAD].Units); else sprintf(s1,FMT51,atime,Node[n].ID,StatTxt[newstat]); writeline(s1); OldStat[Nlinks+i] = newstat; @@ -314,15 +314,15 @@ void writehydstat(int iter, double relerr) /* Display status changes for links */ for (i=1; i<=Nlinks; i++) { - if (S[i] != OldStat[i]) + if (LinkStatus[i] != OldStat[i]) { if (Htime == 0) sprintf(s1,FMT52,atime,LinkTxt[Link[i].Type],Link[i].ID, - StatTxt[S[i]]); + StatTxt[LinkStatus[i]]); else sprintf(s1,FMT53,atime,LinkTxt[Link[i].Type],Link[i].ID, - StatTxt[OldStat[i]],StatTxt[S[i]]); + StatTxt[OldStat[i]],StatTxt[LinkStatus[i]]); writeline(s1); - OldStat[i] = S[i]; + OldStat[i] = LinkStatus[i]; } } writeline(" "); @@ -763,7 +763,7 @@ void writestatchange(int k, char s1, char s2) { /*** Updated 10/25/00 ***/ - setting = K[k]; //Link[k].Kc; + setting = LinkSetting[k]; //Link[k].Kc; switch (Link[k].Type) { @@ -871,7 +871,7 @@ int writehydwarn(int iter, double relerr) /* Check for negative pressures */ for (i=1; i<=Njuncs; i++) { - if (H[i] < Node[i].El && D[i] > 0.0) + if (NodeHead[i] < Node[i].El && NodeDemand[i] > 0.0) { sprintf(Msg,WARN06,clocktime(Atime,Htime)); if (Messageflag) writeline(Msg); @@ -884,10 +884,10 @@ int writehydwarn(int iter, double relerr) for (i=1; i<=Nvalves; i++) { j = Valve[i].Link; - if (S[j] >= XFCV) + if (LinkStatus[j] >= XFCV) { sprintf(Msg,WARN05,LinkTxt[Link[j].Type],Link[j].ID, - StatTxt[S[j]],clocktime(Atime,Htime)); + StatTxt[LinkStatus[j]],clocktime(Atime,Htime)); if (Messageflag) writeline(Msg); flag = 5; } @@ -897,10 +897,10 @@ int writehydwarn(int iter, double relerr) for (i=1; i<=Npumps; i++) { j = Pump[i].Link; - s = S[j]; //(2.00.11 - LR) - if (S[j] >= OPEN) //(2.00.11 - LR) + s = LinkStatus[j]; //(2.00.11 - LR) + if (LinkStatus[j] >= OPEN) //(2.00.11 - LR) { //(2.00.11 - LR) - if (Q[j] > K[j]*Pump[i].Qmax) s = XFLOW; //(2.00.11 - LR) + if (Q[j] > LinkSetting[j]*Pump[i].Qmax) s = XFLOW; //(2.00.11 - LR) if (Q[j] < 0.0) s = XHEAD; //(2.00.11 - LR) } //(2.00.11 - LR) if (s == XHEAD || s == XFLOW) //(2.00.11 - LR) @@ -984,7 +984,7 @@ int disconnected() mcount = Ntanks; for (i=1; i<=Njuncs; i++) { - if (D[i] < 0.0) + if (NodeDemand[i] < 0.0) { mcount++; nodelist[mcount] = i; @@ -999,7 +999,7 @@ int disconnected() count = 0; for (i=1; i<=Njuncs; i++) { - if (!marked[i] && D[i] != 0.0) /* Skip if no demand */ + if (!marked[i] && NodeDemand[i] != 0.0) /* Skip if no demand */ { count++; if (count <= MAXCOUNT && Messageflag) @@ -1068,7 +1068,7 @@ void marknodes(int m, int *nodelist, char *marked) } /* Mark connection node if link not closed */ - if (S[k] > CLOSED) + if (LinkStatus[k] > CLOSED) { marked[j] = 1; m++; diff --git a/src/rules.c b/src/rules.c index aa2e003..7ad038f 100755 --- a/src/rules.c +++ b/src/rules.c @@ -704,7 +704,7 @@ int checkstatus(struct Premise *p) case IS_OPEN: case IS_CLOSED: case IS_ACTIVE: - i = S[p->index]; + i = LinkStatus[p->index]; if (i <= CLOSED) j = IS_CLOSED; else if (i == ACTIVE) j = IS_ACTIVE; else j = IS_OPEN; @@ -736,20 +736,20 @@ int checkvalue(struct Premise *p) /*** Updated 10/25/00 ***/ case r_DEMAND: if (p->object == r_SYSTEM) x = Dsystem*Ucf[DEMAND]; - else x = D[i]*Ucf[DEMAND]; + else x = NodeDemand[i]*Ucf[DEMAND]; break; case r_HEAD: - case r_GRADE: x = H[i]*Ucf[HEAD]; + case r_GRADE: x = NodeHead[i]*Ucf[HEAD]; break; - case r_PRESSURE: x = (H[i]-Node[i].El)*Ucf[PRESSURE]; + case r_PRESSURE: x = (NodeHead[i]-Node[i].El)*Ucf[PRESSURE]; break; - case r_LEVEL: x = (H[i]-Node[i].El)*Ucf[HEAD]; + case r_LEVEL: x = (NodeHead[i]-Node[i].El)*Ucf[HEAD]; break; case r_FLOW: x = ABS(Q[i])*Ucf[FLOW]; break; - case r_SETTING: if (K[i] == MISSING) return(0); - x = K[i]; + case r_SETTING: if (LinkSetting[i] == MISSING) return(0); + x = LinkSetting[i]; switch (Link[i].Type) { case PRV: @@ -761,14 +761,14 @@ int checkvalue(struct Premise *p) case r_FILLTIME: if (i <= Njuncs) return(0); j = i-Njuncs; if (Tank[j].A == 0.0) return(0); - if (D[i] <= TINY) return(0); - x = (Tank[j].Vmax - Tank[j].V)/D[i]; + if (NodeDemand[i] <= TINY) return(0); + x = (Tank[j].Vmax - Tank[j].V)/NodeDemand[i]; break; case r_DRAINTIME: if (i <= Njuncs) return(0); j = i-Njuncs; if (Tank[j].A == 0.0) return(0); - if (D[i] >= -TINY) return(0); - x = (Tank[j].Vmin - Tank[j].V)/D[i]; + if (NodeDemand[i] >= -TINY) return(0); + x = (Tank[j].Vmin - Tank[j].V)/NodeDemand[i]; break; default: return(0); } @@ -875,21 +875,21 @@ int takeactions() flag = FALSE; a = item->action; k = a->link; - s = S[k]; - v = K[k]; + s = LinkStatus[k]; + v = LinkSetting[k]; x = a->setting; /* Switch link from closed to open */ if (a->status == IS_OPEN && s <= CLOSED) { - setlinkstatus(k, 1, &S[k], &K[k]); + setlinkstatus(k, 1, &LinkStatus[k], &LinkSetting[k]); flag = TRUE; } /* Switch link from not closed to closed */ else if (a->status == IS_CLOSED && s > CLOSED) { - setlinkstatus(k, 0, &S[k], &K[k]); + setlinkstatus(k, 0, &LinkStatus[k], &LinkSetting[k]); flag = TRUE; } @@ -905,7 +905,7 @@ int takeactions() } if (ABS(x-v) > tol) { - setlinksetting(k, x, &S[k], &K[k]); + setlinksetting(k, x, &LinkStatus[k], &LinkSetting[k]); flag = TRUE; } } diff --git a/src/vars.h b/src/vars.h index eee1cda..59b92ad 100755 --- a/src/vars.h +++ b/src/vars.h @@ -144,17 +144,17 @@ AUTHOR: L. Rossman SField Field[MAXVAR]; /* Output reporting fields */ /* Array pointers not allocated and freed in same routine */ - char *S, /* Link status */ + char *LinkStatus, /* Link status */ *OldStat; /* Previous link/tank status */ - double *D, /* Node actual demand */ - *C, /* Node actual quality */ + double *NodeDemand, /* Node actual demand */ + *NodeQual, /* Node actual quality */ *E, /* Emitter flows */ - *K, /* Link settings */ + *LinkSetting, /* Link settings */ *Q, /* Link flows */ - *PipeRateCoeff, /* Pipe reaction rate */ + *PipeRateCoeff, /* Pipe reaction rate */ *X, /* General purpose array */ - *XC; /* General purpose array for water quality */ -EXTERN double *H; /* Node heads */ + *TempQual; /* General purpose array for water quality */ +EXTERN double *NodeHead; /* Node heads */ EXTERN double *QTankVolumes; EXTERN double *QLinkFlow; EXTERN STmplist *Patlist; /* Temporary time pattern list */