From 970bbf483e8da7cb7a78a0707eb6986c3f564b35 Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Tue, 23 Apr 2024 17:33:00 -0400 Subject: [PATCH] Revised solver with fewer status checks --- src/hydcoeffs.c | 35 +++++++++--- src/hydsolver.c | 9 ++-- src/hydstatus.c | 140 +----------------------------------------------- src/types.h | 10 +++- 4 files changed, 43 insertions(+), 151 deletions(-) diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c index c00d2f9f..b9bb4954 100644 --- a/src/hydcoeffs.c +++ b/src/hydcoeffs.c @@ -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 ****************************************************************************** */ @@ -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 ); @@ -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); @@ -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) { @@ -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; @@ -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) { @@ -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; diff --git a/src/hydsolver.c b/src/hydsolver.c index 25c6752c..7841c31e 100644 --- a/src/hydsolver.c +++ b/src/hydsolver.c @@ -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 ****************************************************************************** */ @@ -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; @@ -611,7 +613,6 @@ void checkhydbalance(Project *pr, Hydbalance *hbal) } } - int hasconverged(Project *pr, double *relerr, Hydbalance *hbal) /* **-------------------------------------------------------------- diff --git a/src/hydstatus.c b/src/hydstatus.c index e746e586..b8fcb3b4 100644 --- a/src/hydstatus.c +++ b/src/hydstatus.c @@ -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 ****************************************************************************** */ @@ -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); @@ -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]); @@ -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) /* @@ -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) /* **---------------------------------------------------------------- diff --git a/src/types.h b/src/types.h index ff7c3599..3dac4fe5 100755 --- a/src/types.h +++ b/src/types.h @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 07/17/2023 + Last Updated: 04/23/2024 ****************************************************************************** */ @@ -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