Skip to content

Commit

Permalink
Revised solver with fewer status checks
Browse files Browse the repository at this point in the history
  • Loading branch information
@LRossman
LRossman committed Apr 23, 2024
1 parent ef234a1 commit 970bbf4
Show file tree
Hide file tree
Showing 4 changed files with 43 additions and 151 deletions.
35 changes: 27 additions & 8 deletions src/hydcoeffs.c
View file
Original file line number Diff line number Diff line change
@@ -1,13 +1,13 @@
/*
******************************************************************************
Project: OWA EPANET
Version: 2.2
Version: 2.3
Module: hydcoeffs.c
Description: computes coefficients for a hydraulic solution matrix
Authors: see AUTHORS
Copyright: see AUTHORS
License: see LICENSE
Last Updated: 03/29/2023
Last Updated: 04/23/2024
******************************************************************************
*/

Expand All @@ -31,8 +31,8 @@ const double AB = 3.28895476345399058690e-03; // 5.74/(4000^.9)
const double AC = -5.14214965799093883760e-03; // AA*AB

// Definitions of very small and very big coefficients
const double CSMALL = 1.e-6;
const double CBIG = 1.e8;
////const double CSMALL = 1.e-6;
////const double CBIG = 1.e8;

// Exported functions
//void resistcoeff(Project *, int );
Expand Down Expand Up @@ -276,7 +276,7 @@ void headlosscoeffs(Project *pr)
case GPV:
gpvcoeff(pr, k);
break;
case FCV:
case FCV:
case PRV:
case PSV:
if (hyd->LinkSetting[k] == MISSING) valvecoeff(pr, k);
Expand Down Expand Up @@ -648,6 +648,14 @@ void pipecoeff(Project *pr, int k)
return;
}

// Do likewise for CV pipes with negative flow
if (pr->network.Link[k].Type == CVPIPE && hyd->LinkFlow[k] < 0.0)
{
hyd->P[k] = 1.0 / CBIG;
hyd->Y[k] = hyd->LinkFlow[k];
return;
}

// Use custom function for Darcy-Weisbach formula
if (hyd->Formflag == DW)
{
Expand Down Expand Up @@ -681,7 +689,7 @@ void pipecoeff(Project *pr, int k)

// Adjust head loss sign for flow direction
hloss *= SGN(hyd->LinkFlow[k]);

// P and Y coeffs.
hyd->P[k] = 1.0 / hgrad;
hyd->Y[k] = hloss / hgrad;
Expand Down Expand Up @@ -812,12 +820,23 @@ void pumpcoeff(Project *pr, int k)
hyd->Y[k] = hyd->LinkFlow[k];
return;
}

// Obtain reference to pump object
q = ABS(hyd->LinkFlow[k]);
p = findpump(&pr->network, k);
pump = &pr->network.Pump[p];

// Adjust head loss and gradient to prevent negative flow
if (hyd->LinkFlow[k] < 0.0)
{
double hmax = SQR(hyd->LinkSetting[k]) * pump->Hmax;
hloss = -hmax + CBIG * hyd->LinkFlow[k];
hgrad = CBIG;
hyd->P[k] = 1.0 / hgrad;
hyd->Y[k] = hloss / hgrad;
return;
}

// If no pump curve treat pump as an open valve
if (pump->Ptype == NOCURVE)
{
Expand Down Expand Up @@ -1221,7 +1240,7 @@ void fcvcoeff(Project *pr, int k, int n1, int n2)
// flow setting as external demand at upstream node
// and external supply at downstream node.

if (hyd->LinkStatus[k] == ACTIVE)
if (hyd->LinkFlow[k] > q)
{
hyd->Xflow[n1] -= q;
hyd->Xflow[n2] += q;
Expand Down
9 changes: 5 additions & 4 deletions src/hydsolver.c
View file
Original file line number Diff line number Diff line change
@@ -1,14 +1,14 @@
/*
******************************************************************************
Project: OWA EPANET
Version: 2.2
Version: 2.3
Module: hydsolver.c
Description: computes flows and pressures throughout a pipe network using
Todini's Global Gradient Algorithm
Authors: see AUTHORS
Copyright: see AUTHORS
License: see LICENSE
Last Updated: 02/14/2022
Last Updated: 04/23/2024
******************************************************************************
*/

Expand Down Expand Up @@ -595,8 +595,10 @@ void checkhydbalance(Project *pr, Hydbalance *hbal)
headlosscoeffs(pr);
for (k = 1; k <= net->Nlinks; k++)
{
// Skip closed links and active valves
if (hyd->LinkStatus[k] <= CLOSED) continue;
if (hyd->P[k] == 0.0) continue;
if (hyd->P[k] <= 1.0/CBIG) continue;

link = &net->Link[k];
n1 = link->N1;
n2 = link->N2;
Expand All @@ -611,7 +613,6 @@ void checkhydbalance(Project *pr, Hydbalance *hbal)
}
}


int hasconverged(Project *pr, double *relerr, Hydbalance *hbal)
/*
**--------------------------------------------------------------
Expand Down
140 changes: 2 additions & 138 deletions src/hydstatus.c
View file
Original file line number Diff line number Diff line change
@@ -1,13 +1,13 @@
/*
******************************************************************************
Project: OWA EPANET
Version: 2.2
Version: 2.3
Module: hydstatus.c
Description: updates hydraulic status of network elements
Authors: see AUTHORS
Copyright: see AUTHORS
License: see LICENSE
Last Updated: 02/03/2023
Last Updated: 04/23/2024
******************************************************************************
*/

Expand All @@ -22,11 +22,8 @@ int valvestatus(Project *);
int linkstatus(Project *);

// Local functions
static StatusType cvstatus(Project *, StatusType, double, double);
static StatusType pumpstatus(Project *, int, double);
static StatusType prvstatus(Project *, int, StatusType, double, double, double);
static StatusType psvstatus(Project *, int, StatusType, double, double, double);
static StatusType fcvstatus(Project *, int, StatusType, double, double);
static void tankstatus(Project *, int, int, double);


Expand Down Expand Up @@ -135,25 +132,6 @@ int linkstatus(Project *pr)
hyd->LinkStatus[k] = OPEN;
}

// Check for status changes in CVs and pumps
if (link->Type == CVPIPE)
{
hyd->LinkStatus[k] = cvstatus(pr, hyd->LinkStatus[k], dh,
hyd->LinkFlow[k]);
}
if (link->Type == PUMP && hyd->LinkStatus[k] >= OPEN &&
hyd->LinkSetting[k] > 0.0)
{
hyd->LinkStatus[k] = pumpstatus(pr, k, -dh);
}

// Check for status changes in non-fixed FCVs
if (link->Type == FCV && hyd->LinkSetting[k] != MISSING)
{
hyd->LinkStatus[k] = fcvstatus(pr, k, status, hyd->NodeHead[n1],
hyd->NodeHead[n2]);
}

// Check for flow into (out of) full (empty) tanks
if (n1 > net->Njuncs) tankstatus(pr, k, n1, hyd->LinkFlow[k]);
if (n2 > net->Njuncs) tankstatus(pr, k, n2, -hyd->LinkFlow[k]);
Expand All @@ -172,72 +150,6 @@ int linkstatus(Project *pr)
}


StatusType cvstatus(Project *pr, StatusType s, double dh, double q)
/*
**--------------------------------------------------
** Input: s = current link status
** dh = head loss across link
** q = link flow
** Output: returns new link status
** Purpose: updates status of a check valve link.
**--------------------------------------------------
*/
{
Hydraul *hyd = &pr->hydraul;

// Prevent reverse flow through CVs
if (ABS(dh) > hyd->Htol)
{
if (dh < -hyd->Htol) return CLOSED;
else if (q < -hyd->Qtol) return CLOSED;
else return OPEN;
}
else
{
if (q < -hyd->Qtol) return CLOSED;
else return s;
}
}


StatusType pumpstatus(Project *pr, int k, double dh)
/*
**--------------------------------------------------
** Input: k = link index
** dh = head gain across link
** Output: returns new pump status
** Purpose: updates status of an open pump.
**--------------------------------------------------
*/
{
Hydraul *hyd = &pr->hydraul;
Network *net = &pr->network;

int p;
double hmax;

// Find maximum head (hmax) pump can deliver
p = findpump(net, k);
if (net->Pump[p].Ptype == CONST_HP)
{
// Use huge value for constant HP pump
hmax = BIG;
if (hyd->LinkFlow[k] < TINY) return TEMPCLOSED;
}
else
{
// Use speed-adjusted shut-off head for other pumps
hmax = SQR(hyd->LinkSetting[k]) * net->Pump[p].Hmax;
}

// Check if currrent head gain exceeds pump's max. head
if (dh > hmax + hyd->Htol) return XHEAD;

// No check is made to see if flow exceeds pump's max. flow
return OPEN;
}


StatusType prvstatus(Project *pr, int k, StatusType s, double hset,
double h1, double h2)
/*
Expand Down Expand Up @@ -358,54 +270,6 @@ StatusType psvstatus(Project *pr, int k, StatusType s, double hset,
}


StatusType fcvstatus(Project *pr, int k, StatusType s, double h1, double h2)
/*
**-----------------------------------------------------------
** Input: k = link index
** s = current status
** h1 = head at upstream node
** h2 = head at downstream node
** Output: returns new valve status
** Purpose: updates status of a flow control valve.
**
** Valve status changes to XFCV if flow reversal.
** If current status is XFCV and current flow is
** above setting, then valve becomes active.
** If current status is XFCV, and current flow
** positive but still below valve setting, then
** status remains same.
**-----------------------------------------------------------
*/
{
Hydraul *hyd = &pr->hydraul;
StatusType status; // New valve status

status = s;
if (h1 - h2 < -hyd->Htol)
{
status = XFCV;
}
else if (hyd->LinkFlow[k] < -hyd->Qtol)
{
status = XFCV;
}
else if (s == XFCV && hyd->LinkFlow[k] >= hyd->LinkSetting[k])
{
status = ACTIVE;
}

// Active valve's loss coeff. can't be < fully open loss coeff.
else if (status == ACTIVE)
{
if ((h1 - h2) / SQR(hyd->LinkFlow[k]) < pr->network.Link[k].Km)
{
status = XFCV;
}
}
return status;
}


void tankstatus(Project *pr, int k, int n, double q)
/*
**----------------------------------------------------------------
Expand Down
10 changes: 9 additions & 1 deletion src/types.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
Authors: see AUTHORS
Copyright: see AUTHORS
License: see LICENSE
Last Updated: 07/17/2023
Last Updated: 04/23/2024
******************************************************************************
*/

Expand Down Expand Up @@ -63,6 +63,14 @@ typedef int INT4;
#define PI 3.141592654
#endif

/*
----------------------------------------------
Very small and very big coefficients
----------------------------------------------
*/
#define CSMALL 1.e-6
#define CBIG 1.e8

/*
----------------------------------------------
Flow units conversion factors
Expand Down

0 comments on commit 970bbf4

Please sign in to comment.