JIANG/EPANET
1
0
Fork 0
Code Issues Pull Requests Actions 4 Packages Projects Releases Wiki Activity

Fixes for PDA

Browse Source
This commit is contained in:
Lew Rossman
2024-06-26 17:50:24 -04:00
parent 30a3fcb4b9
commit 68b73a14f1
3 changed files with 31 additions and 23 deletions
Download Patch File Download Diff File Expand all files Collapse all files

25
src/flowbalance.c
View File

@@ -73,9 +73,7 @@ void updateflowbalance(Project *pr, long hstep)
flowBalance.storageDemand = 0.0; flowBalance.storageDemand = 0.0;
fullDemand = 0.0; fullDemand = 0.0;
// Initialize demand deficiency & leakage loss // Initialize leakage loss
hyd->DeficientNodes = 0;
hyd->DemandReduction = 0.0;
hyd->LeakageLoss = 0.0; hyd->LeakageLoss = 0.0;
// Examine each junction node // Examine each junction node
@@ -104,11 +102,8 @@ void updateflowbalance(Project *pr, long hstep)
if (hyd->DemandModel == PDA && hyd->FullDemand[i] > 0.0) if (hyd->DemandModel == PDA && hyd->FullDemand[i] > 0.0)
{ {
deficit = hyd->FullDemand[i] - hyd->DemandFlow[i]; deficit = hyd->FullDemand[i] - hyd->DemandFlow[i];
if (deficit > TINY) if (deficit > 0.0)
{
hyd->DeficientNodes++;
flowBalance.deficitDemand += deficit; flowBalance.deficitDemand += deficit;
}
} }
} }
@@ -118,8 +113,8 @@ void updateflowbalance(Project *pr, long hstep)
i = net->Tank[j].Node; i = net->Tank[j].Node;
v = hyd->NodeDemand[i]; v = hyd->NodeDemand[i];
// For a snapshot analysis or a reservoir node // For a reservoir node
if (time->Dur == 0 || net->Tank[j].A == 0.0) if (net->Tank[j].A == 0.0)
{ {
if (v >= 0.0) if (v >= 0.0)
flowBalance.totalOutflow += v; flowBalance.totalOutflow += v;
@@ -127,16 +122,16 @@ void updateflowbalance(Project *pr, long hstep)
flowBalance.totalInflow += (-v); flowBalance.totalInflow += (-v);
} }
// For tank under extended period analysis // For tank
else else
flowBalance.storageDemand += v; flowBalance.storageDemand += v;
} }
// Find % demand reduction & % leakage for current period // Find % leakage for current period
if (fullDemand > 0.0) v = flowBalance.totalInflow;
hyd->DemandReduction = flowBalance.deficitDemand / fullDemand * 100.0; if (flowBalance.storageDemand < 0.0) v += (-flowBalance.storageDemand);
if (flowBalance.totalInflow > 0.0) if (v > 0.0)
hyd->LeakageLoss = flowBalance.leakageDemand / flowBalance.totalInflow * 100.0; hyd->LeakageLoss = flowBalance.leakageDemand / v * 100.0;
// Update flow balance for entire run // Update flow balance for entire run
hyd->FlowBalance.totalInflow += flowBalance.totalInflow * dt; hyd->FlowBalance.totalInflow += flowBalance.totalInflow * dt;

26
src/hydsolver.c
View File

@@ -8,7 +8,7 @@
Authors: see AUTHORS Authors: see AUTHORS
Copyright: see AUTHORS Copyright: see AUTHORS
License: see LICENSE License: see LICENSE
Last Updated: 06/15/2024 Last Updated: 06/26/2024
****************************************************************************** ******************************************************************************
*/ */
@@ -703,10 +703,15 @@ int pdaconverged(Project *pr)
Hydraul *hyd = &pr->hydraul; Hydraul *hyd = &pr->hydraul;
const double QTOL = 0.0001; // 0.0001 cfs ~= 0.005 gpm ~= 0.2 lpm) const double QTOL = 0.0001; // 0.0001 cfs ~= 0.005 gpm ~= 0.2 lpm)
int i; int i, converged = 1;
double totalDemand = 0.0, totalReduction = 0.0;
double dp = hyd->Preq - hyd->Pmin; double dp = hyd->Preq - hyd->Pmin;
double p, q, r; double p, q, r;
hyd->DeficientNodes = 0;
hyd->DemandReduction = 0.0;
// Examine each network junction // Examine each network junction
for (i = 1; i <= pr->network.Njuncs; i++) for (i = 1; i <= pr->network.Njuncs; i++)
{ {
@@ -726,9 +731,20 @@ int pdaconverged(Project *pr)
} }
// Check if demand has not converged // Check if demand has not converged
if (fabs(q - hyd->DemandFlow[i]) > QTOL) return 0; if (fabs(q - hyd->DemandFlow[i]) > QTOL)
converged = 0;
// Accumulate demand deficient node count and demand deficit
if (hyd->DemandFlow[i] + QTOL < hyd->FullDemand[i])
{
hyd->DeficientNodes++;
totalDemand += hyd->FullDemand[i];
totalReduction += hyd->FullDemand[i] - hyd->DemandFlow[i];
}
} }
return 1; if (totalDemand > 0.0)
hyd->DemandReduction = totalReduction / totalDemand * 100.0;
return converged;
} }

3
tests/test_pda.cpp
View File

@@ -79,9 +79,6 @@ BOOST_AUTO_TEST_CASE(test_pda_model)
error = EN_getnodeindex(ph, (char *)"21", &index); error = EN_getnodeindex(ph, (char *)"21", &index);
BOOST_REQUIRE(error == 0); BOOST_REQUIRE(error == 0);
error = EN_getnodevalue(ph, index, EN_DEMANDDEFICIT, &reduction); error = EN_getnodevalue(ph, index, EN_DEMANDDEFICIT, &reduction);
printf("\nreduction = %f", reduction);
BOOST_REQUIRE(error == 0); BOOST_REQUIRE(error == 0);
BOOST_REQUIRE(abs(reduction - 413.67) < 0.01); BOOST_REQUIRE(abs(reduction - 413.67) < 0.01);