diff --git a/include/epanet2_2.h b/include/epanet2_2.h index 18ef127..4393073 100644 --- a/include/epanet2_2.h +++ b/include/epanet2_2.h @@ -638,6 +638,16 @@ typedef struct Project *EN_Project; */ int DLLEXPORT EN_getstatistic(EN_Project ph, int type, double* out_value); + + /** + @brief Get information about upcoming time step events, and what causes them. + @param ph an EPANET project handle. + @param[out] eventType the type of event that will occur. + @param[out] duration the amount of time in the future this event will occur + @param[out] elementIndex the index of the element causing the event. + **/ + int DLLEXPORT EN_timeToNextEvent(EN_Project ph, EN_TimestepEvent *eventType, long *duration, int *elementIndex); + /** @brief Retrieves the order in which a node or link appears in an @ref OutFile "output file". @param ph an EPANET project handle. diff --git a/include/epanet2_enums.h b/include/epanet2_enums.h index a000b0c..e0aa9ae 100644 --- a/include/epanet2_enums.h +++ b/include/epanet2_enums.h @@ -125,6 +125,18 @@ typedef enum { EN_NEXTEVENTTANK = 15 //!< Index of tank with shortest time to become empty or full (read only) } EN_TimeParameter; + +/** +These are the types of events that can cause a timestep to end. +**/ +typedef enum { + EN_STEP_REPORT = 0, + EN_STEP_HYD = 1, + EN_STEP_WQ = 2, + EN_STEP_TANKEVENT = 3, + EN_STEP_CONTROLEVENT = 4 +} EN_TimestepEvent; + /// Analysis convergence statistics /** These statistics report the convergence criteria for the most current hydraulic analysis diff --git a/src/epanet.c b/src/epanet.c index 9a3c194..9c5783c 100644 --- a/src/epanet.c +++ b/src/epanet.c @@ -1610,6 +1610,41 @@ int DLLEXPORT EN_settimeparam(EN_Project p, int param, long value) return 0; } + +/// get the time to next event, and give a reason for the time step truncation +int DLLEXPORT EN_timeToNextEvent(EN_Project p, EN_TimestepEvent *eventType, long *duration, int *elementIndex) +{ + Times *time = &p->times; + long hydStep, tankStep, controlStep; + + hydStep = time->Hstep; + tankStep = hydStep; + controlStep = hydStep; + + int iTank = tanktimestep(p, &tankStep); + int iControl = controltimestep(p, &controlStep); + + // return the lesser of the three step lengths + if (controlStep < tankStep) { + *eventType = EN_STEP_CONTROLEVENT; + *duration = controlStep; + *elementIndex = iControl; + } + else if (tankStep < hydStep) { + *eventType = EN_STEP_TANKEVENT; + *duration = tankStep; + *elementIndex = iTank; + } + else { + *eventType = EN_STEP_HYD; + *duration = hydStep; + *elementIndex = 0; + } + + return 0; +} + + int DLLEXPORT EN_getqualinfo(EN_Project p, int *qualType, char *chemName, char *chemUnits, int *traceNode) /*---------------------------------------------------------------- @@ -1774,9 +1809,9 @@ int DLLEXPORT EN_addnode(EN_Project p, char *id, int nodeType, int *index) // Check if a node with same id already exists if (EN_getnodeindex(p, id, &i) == 0) return 215; - + // Check for valid node type - if (nodeType < EN_JUNCTION || nodeType > EN_TANK) return 251; + if (nodeType < EN_JUNCTION || nodeType > EN_TANK) return 251; // Grow node-related arrays to accomodate the new node size = (net->Nnodes + 2) * sizeof(Snode); @@ -1797,7 +1832,7 @@ int DLLEXPORT EN_addnode(EN_Project p, char *id, int nodeType, int *index) hashtable_update(net->NodeHashTable, net->Node[i].ID, i + 1); net->Node[i + 1] = net->Node[i]; } - + // set index of new Junction node net->Njuncs++; nIdx = net->Njuncs; @@ -2254,20 +2289,20 @@ int DLLEXPORT EN_getnodevalue(EN_Project p, int index, int property, double *val if (Node[index].Type != TANK) return 0; v = Tank[index - nJuncs].CanOverflow; break; - + case EN_DEMANDDEFICIT: if (index > nJuncs) return 0; // After an analysis, DemandFlow contains node's required demand // while NodeDemand contains delivered demand + emitter flow if (hyd->DemandFlow[index] < 0.0) return 0; - v = (hyd->DemandFlow[index] - + v = (hyd->DemandFlow[index] - (hyd->NodeDemand[index] - hyd->EmitterFlow[index])) * Ucf[FLOW]; break; - + case EN_NODE_INCONTROL: v = (double)incontrols(p, NODE, index); break; - + default: return 251; } @@ -2483,7 +2518,7 @@ int DLLEXPORT EN_setnodevalue(EN_Project p, int index, int property, double valu Tank[j].Vmin = tankvolume(p, j, Tank[j].Hmin); // new min. volume Tank[j].V0 = tankvolume(p, j, Tank[j].H0); // new init. volume Tank[j].Vmax = tankvolume(p, j, Tank[j].Hmax); // new max. volume - Tank[j].A = (curve->Y[n] - curve->Y[0]) / // nominal area + Tank[j].A = (curve->Y[n] - curve->Y[0]) / // nominal area (curve->X[n] - curve->X[0]); break; @@ -3793,7 +3828,7 @@ int DLLEXPORT EN_getlinkvalue(EN_Project p, int index, int property, double *val v = (double)Pump[findpump(&p->network, index)].Epat; } break; - + case EN_GPV_CURVE: if (Link[index].Type == GPV) { @@ -3804,7 +3839,7 @@ int DLLEXPORT EN_getlinkvalue(EN_Project p, int index, int property, double *val case EN_LINK_INCONTROL: v = (double)incontrols(p, LINK, index); break; - + default: return 251; } @@ -4011,7 +4046,7 @@ int DLLEXPORT EN_setlinkvalue(EN_Project p, int index, int property, double valu net->Pump[pumpIndex].Epat = patIndex; } break; - + case EN_GPV_CURVE: if (Link[index].Type == GPV) { @@ -4077,20 +4112,20 @@ int DLLEXPORT EN_getvertexcount(EN_Project p, int index, int *count) */ { Network *net = &p->network; - + Slink *Link = net->Link; Pvertices vertices; - + // Check that link exists *count = 0; if (!p->Openflag) return 102; if (index <= 0 || index > net->Nlinks) return 204; - + // Set count to number of vertices vertices = Link[index].Vertices; if (vertices) *count = vertices->Npts; return 0; -} +} int DLLEXPORT EN_getvertex(EN_Project p, int index, int vertex, double *x, double *y) /*---------------------------------------------------------------- @@ -4104,22 +4139,22 @@ int DLLEXPORT EN_getvertex(EN_Project p, int index, int vertex, double *x, doubl */ { Network *net = &p->network; - + Slink *Link = net->Link; Pvertices vertices; - + // Check that link exists *x = MISSING; *y = MISSING; if (!p->Openflag) return 102; if (index <= 0 || index > net->Nlinks) return 204; - + // Check that vertex exists vertices = Link[index].Vertices; if (vertices == NULL) return 255; if (vertex <= 0 || vertex > vertices->Npts) return 255; *x = vertices->X[vertex - 1]; - *y = vertices->Y[vertex - 1]; + *y = vertices->Y[vertex - 1]; return 0; } @@ -4135,23 +4170,23 @@ int DLLEXPORT EN_setvertex(EN_Project p, int index, int vertex, double x, double */ { Network *net = &p->network; - + Slink *Link = net->Link; Pvertices vertices; - + // Check that link exists if (!p->Openflag) return 102; if (index <= 0 || index > net->Nlinks) return 204; - + // Check that vertex exists vertices = Link[index].Vertices; if (vertices == NULL) return 255; if (vertex <= 0 || vertex > vertices->Npts) return 255; vertices->X[vertex - 1] = x; - vertices->Y[vertex - 1] = y; + vertices->Y[vertex - 1] = y; return 0; } - + int DLLEXPORT EN_setvertices(EN_Project p, int index, double *x, double *y, int count) /*---------------------------------------------------------------- ** Input: index = link index @@ -4164,11 +4199,11 @@ int DLLEXPORT EN_setvertices(EN_Project p, int index, double *x, double *y, int */ { Network *net = &p->network; - + Slink *link; int i; int err = 0; - + // Check that link exists if (!p->Openflag) return 102; if (index <= 0 || index > net->Nlinks) return 204; @@ -4176,7 +4211,7 @@ int DLLEXPORT EN_setvertices(EN_Project p, int index, double *x, double *y, int // Delete existing set of vertices freelinkvertices(link); - + // Add each new vertex to the link for (i = 0; i < count; i++) { @@ -4185,7 +4220,7 @@ int DLLEXPORT EN_setvertices(EN_Project p, int index, double *x, double *y, int } if (err) freelinkvertices(link); return err; -} +} /******************************************************************** @@ -4269,16 +4304,16 @@ int DLLEXPORT EN_setheadcurveindex(EN_Project p, int linkIndex, int curveIndex) pump = &p->network.Pump[pumpIndex]; oldCurveIndex = pump->Hcurve; newCurveType = p->network.Curve[curveIndex].Type; - + // Assign the new curve to the pump pump->Ptype = NOCURVE; pump->Hcurve = curveIndex; if (curveIndex == 0) return 0; - + // Update the pump's head curve parameters (which also changes // the new curve's Type to PUMP_CURVE) err = updatepumpparams(p, pumpIndex); - + // If the parameter updating failed (new curve was not a valid pump curve) // restore the pump's original curve and its parameters if (err > 0) @@ -4288,8 +4323,8 @@ int DLLEXPORT EN_setheadcurveindex(EN_Project p, int linkIndex, int curveIndex) pump->Hcurve = oldCurveIndex; if (oldCurveIndex == 0) return err; updatepumpparams(p, pumpIndex); - } - + } + // Convert the units of the updated pump parameters to feet and cfs if (pump->Ptype == POWER_FUNC) { @@ -4778,7 +4813,7 @@ int DLLEXPORT EN_setcurvetype(EN_Project p, int index, int type) Network *net = &p->network; if (!p->Openflag) return 102; if (index < 1 || index > net->Ncurves) return 206; - if (type < 0 || type > EN_GENERIC_CURVE) return 251; + if (type < 0 || type > EN_GENERIC_CURVE) return 251; net->Curve[index].Type = type; return 0; } @@ -4852,7 +4887,7 @@ int DLLEXPORT EN_setcurvevalue(EN_Project p, int curveIndex, int pointIndex, // Insert new point into curve curve->X[n] = x; curve->Y[n] = y; - + // Adjust parameters for pumps using curve as a head curve return adjustpumpparams(p, curveIndex); } @@ -4926,7 +4961,7 @@ int DLLEXPORT EN_setcurve(EN_Project p, int index, double *xValues, curve->X[j] = xValues[j]; curve->Y[j] = yValues[j]; } - + // Adjust parameters for pumps using curve as a head curve return adjustpumpparams(p, index); } diff --git a/src/funcs.h b/src/funcs.h index 366e5c7..1e6be71 100755 --- a/src/funcs.h +++ b/src/funcs.h @@ -155,6 +155,7 @@ void closehyd(Project *); void setlinkstatus(Project *, int, char, StatusType *, double *); void setlinksetting(Project *, int, double, StatusType *, double *); int tanktimestep(Project *, long *); +int controltimestep(Project *, long *); void getenergy(Project *, int, double *, double *); double tankvolume(Project *, int, double); double tankgrade(Project *, int, double); @@ -164,7 +165,7 @@ double tankgrade(Project *, int, double); void resistcoeff(Project *, int); void headlosscoeffs(Project *); void matrixcoeffs(Project *); -void emitterheadloss(Project *, int, double *, double *); +void emitterheadloss(Project *, int, double *, double *); void demandheadloss(Project *, int, double, double, double *, double *); // ------- QUALITY.C -------------------- diff --git a/src/hydraul.c b/src/hydraul.c index 36cd67b..fecc9f8 100755 --- a/src/hydraul.c +++ b/src/hydraul.c @@ -34,7 +34,6 @@ void initlinkflow(Project *, int, char, double); void demands(Project *); int controls(Project *); long timestep(Project *); -void controltimestep(Project *, long *); void ruletimestep(Project *, long *); void addenergy(Project *, long); void tanklevels(Project *, long); @@ -192,7 +191,7 @@ int runhyd(Project *pr, long *t) int iter; // Iteration count int errcode; // Error code double relerr; // Solution accuracy - + // Find new demands & control actions *t = time->Htime; demands(pr); @@ -390,7 +389,7 @@ void setlinkstatus(Project *pr, int index, char value, StatusType *s, double *k if (t == PUMP) { *k = 1.0; - // Check if a re-opened pump needs its flow reset + // Check if a re-opened pump needs its flow reset if (*s == CLOSED) resetpumpflow(pr, index); } if (t > PUMP && t != GPV) *k = MISSING; @@ -596,11 +595,11 @@ int controls(Project *pr) k1 = hyd->LinkSetting[k]; k2 = k1; if (link->Type > PIPE) k2 = control->Setting; - + // Check if a re-opened pump needs its flow reset if (link->Type == PUMP && s1 == CLOSED && s2 == OPEN) resetpumpflow(pr, k); - + if (s1 != s2 || k1 != k2) { hyd->LinkStatus[k] = s2; @@ -704,7 +703,7 @@ int tanktimestep(Project *pr, long *tstep) } -void controltimestep(Project *pr, long *tstep) +int controltimestep(Project *pr, long *tstep) /* **------------------------------------------------------------------ ** Input: *tstep = current time step @@ -717,7 +716,7 @@ void controltimestep(Project *pr, long *tstep) Network *net = &pr->network; Hydraul *hyd = &pr->hydraul; - int i, j, k, n; + int i, j, k, n, controlIndex; double h, q, v; long t, t1, t2; Slink *link; @@ -774,9 +773,14 @@ void controltimestep(Project *pr, long *tstep) k = control->Link; link = &net->Link[k]; if ( (link->Type > PIPE && hyd->LinkSetting[k] != control->Setting) - || (hyd->LinkStatus[k] != control->Status) ) *tstep = t; + || (hyd->LinkStatus[k] != control->Status) ) + { + *tstep = t; + controlIndex = i; + } } } + return controlIndex; } @@ -1009,7 +1013,7 @@ void getallpumpsenergy(Project *pr) getenergy(pr, pump->Link, &(pump->Energy.CurrentPower), &(pump->Energy.CurrentEffic)); } -} +} void tanklevels(Project *pr, long tstep) @@ -1129,6 +1133,5 @@ void resetpumpflow(Project *pr, int i) Network *net = &pr->network; Spump *pump = &net->Pump[findpump(net, i)]; if (pump->Ptype == CONST_HP) - pr->hydraul.LinkFlow[i] = pump->Q0; + pr->hydraul.LinkFlow[i] = pump->Q0; } -