Merge pull request #862 from lbutler/dev

Decoupled pressure units from the flow unit system
This commit is contained in:
Lew Rossman
2025-07-13 15:27:42 -04:00
committed by GitHub
13 changed files with 266 additions and 30 deletions

View File

@@ -42,6 +42,8 @@ This document describes the changes and updates that have been made in version 2
- `EN_PRESS_UNITS` can now be used with `EN_getoption` and `EN_setoption` to get or set the pressure unit used in EPANET.
- Decoupled pressure units from the flow unit system, allowing them to be set independently to support mixed-unit conventions (e.g., using LPS for flow and PSI for pressure).
- The following constants can be used with EN_getnodevalue to retrieve the components of a node's total demand at a given point in time:
- `EN_FULLDEMAND` - the consumer demand requested
- `EN_DEMANDFLOW` - the consumer demand delivered

View File

@@ -166,6 +166,8 @@ Public Const EN_CMS = 10
Public Const EN_PSI = 0 ' Pressure units types
Public Const EN_KPA = 1
Public Const EN_METERS = 2
Public Const EN_BAR = 3
Public Const EN_FEET = 4
Public Const EN_DDA = 0 ' Demand driven analysis
Public Const EN_PDA = 1 ' Pressure driven analysis

View File

@@ -170,6 +170,8 @@ namespace EpanetCSharpLibrary
public const int EN_PSI = 0; //Pressure units types
public const int EN_KPA = 1;
public const int EN_METERS = 2;
public const int EN_BAR = 3;
public const int EN_FEET = 4;
public const int EN_DDA = 0; //Demand driven analysis
public const int EN_PDA = 1; //Pressure driven analysis

View File

@@ -174,6 +174,8 @@ const
EN_PSI = 0; { Pressure units types }
EN_KPA = 1;
EN_METERS = 2;
EN_BAR = 3;
EN_FEET = 4;
EN_DDA = 0; { Demand model types }
EN_PDA = 1;

View File

@@ -161,6 +161,8 @@ Public Const EN_CMS = 10
Public Const EN_PSI = 0 ' Pressure units types
Public Const EN_KPA = 1
Public Const EN_METERS = 2
Public Const EN_BAR = 3
Public Const EN_FEET = 4
Public Const EN_DDA = 0 ' Demand driven analysis
Public Const EN_PDA = 1 ' Pressure driven analysis

View File

@@ -301,14 +301,14 @@ typedef enum {
/// Pressure units
/**
The available choices for pressure units for the `EN_PRESS_UNITS` option in @ref EN_getoption
and @ref EN_setoption. For networks using US Customary units for flow (`EN_CFS` through
`EN_AFD`) pressure units can only be set as PSI. For network using metric units, you can
select either `EN_METERS` or `EN_KPA`.
and @ref EN_setoption.
*/
typedef enum {
EN_PSI = 0, //!< Pounds per square inch
EN_KPA = 1, //!< Kilopascals
EN_METERS = 2 //!< Meters
EN_METERS = 2, //!< Meters
EN_BAR = 3, //!< Bar
EN_FEET = 4 //!< Feet
} EN_PressUnits;
/// Demand models

View File

@@ -76,7 +76,9 @@ char *FlowUnitsTxt[] = {w_CFS,
char *PressUnitsTxt[] = {w_PSI,
w_KPA,
w_METERS};
w_METERS,
w_BAR,
w_FEET};
char *DemandModelTxt[] = { w_DDA,
w_PDA,

View File

@@ -1355,7 +1355,12 @@ int DLLEXPORT EN_setoption(EN_Project p, int option, double value)
case EN_SP_GRAVITY:
if (value <= 0.0) return 213;
Ucf[PRESSURE] *= (value / hyd->SpGrav);
if (p->parser.Pressflag == PSI ||
p->parser.Pressflag == KPA ||
p->parser.Pressflag == BAR)
{
Ucf[PRESSURE] *= (value / hyd->SpGrav);
}
hyd->SpGrav = value;
break;
@@ -1410,9 +1415,7 @@ int DLLEXPORT EN_setoption(EN_Project p, int option, double value)
case EN_PRESS_UNITS:
unit = ROUND(value);
if (unit < 0 || unit > METERS) return 205;
if (p->parser.Unitsflag == US && unit > PSI) return 0;
if (p->parser.Unitsflag == SI && unit == PSI) return 0;
if (unit < 0 || unit > FEET) return 205;
p->parser.Pressflag = unit;
dfactor = Ucf[DEMAND];
@@ -1467,7 +1470,7 @@ int DLLEXPORT EN_setflowunits(EN_Project p, int units)
{
Network *net = &p->network;
int i, j;
int i, j, oldUnitFlag;
double qfactor, vfactor, hfactor, efactor, pfactor, dfactor, xfactor, yfactor;
double dcf, pcf, hcf, qcf;
double *Ucf = p->Ucf;
@@ -1482,6 +1485,7 @@ int DLLEXPORT EN_setflowunits(EN_Project p, int units)
pfactor = Ucf[PRESSURE];
dfactor = Ucf[DEMAND];
oldUnitFlag = p->parser.Unitsflag;
p->parser.Flowflag = units;
switch (units)
{
@@ -1499,8 +1503,11 @@ int DLLEXPORT EN_setflowunits(EN_Project p, int units)
}
// Revise pressure units depending on flow units
if (p->parser.Unitsflag != SI) p->parser.Pressflag = PSI;
else if (p->parser.Pressflag == PSI) p->parser.Pressflag = METERS;
if (oldUnitFlag != p->parser.Unitsflag)
{
if (p->parser.Unitsflag == US) p->parser.Pressflag = PSI;
else p->parser.Pressflag = METERS;
}
initunits(p);
// Update pressure units in rules
@@ -2235,8 +2242,10 @@ int DLLEXPORT EN_getnodevalue(EN_Project p, int index, int property, double *val
Network *net = &p->network;
Hydraul *hyd = &p->hydraul;
Quality *qual = &p->quality;
Parser *parser = &p->parser;
double v = 0.0;
double ecfTmp; // Unit conversion factor for emitter pressure
Psource source;
Snode *Node = net->Node;
@@ -2281,7 +2290,8 @@ int DLLEXPORT EN_getnodevalue(EN_Project p, int index, int property, double *val
v = 0.0;
if (Node[index].Ke > 0.0)
{
v = Ucf[FLOW] / pow((Ucf[PRESSURE] * Node[index].Ke), (1.0 / hyd->Qexp));
ecfTmp = (parser->Unitsflag == US) ? (PSIperFT * hyd->SpGrav) : MperFT;
v = Ucf[FLOW] / pow((ecfTmp * Node[index].Ke), (1.0 / hyd->Qexp));
}
break;
@@ -2470,6 +2480,7 @@ int DLLEXPORT EN_setnodevalue(EN_Project p, int index, int property, double valu
Network *net = &p->network;
Hydraul *hyd = &p->hydraul;
Quality *qual = &p->quality;
Parser *parser = &p->parser;
Snode *Node = net->Node;
Stank *Tank = net->Tank;
@@ -2483,7 +2494,7 @@ int DLLEXPORT EN_setnodevalue(EN_Project p, int index, int property, double valu
int i, j, n;
Psource source;
double hTmp;
double hTmp, ecfTmp;
if (!p->Openflag) return 102;
if (index <= 0 || index > nNodes) return 203;
@@ -2525,7 +2536,11 @@ int DLLEXPORT EN_setnodevalue(EN_Project p, int index, int property, double valu
case EN_EMITTER:
if (index > nJuncs) return 0;
if (value < 0.0) return 209;
if (value > 0.0) value = pow((Ucf[FLOW] / value), hyd->Qexp) / Ucf[PRESSURE];
if (value > 0.0)
{
ecfTmp = (parser->Unitsflag == US) ? (PSIperFT * hyd->SpGrav) : MperFT;
value = pow((Ucf[FLOW] / value), hyd->Qexp) / ecfTmp;
}
Node[index].Ke = value;
if (hyd->EmitterFlow[index] == 0.0) hyd->EmitterFlow[index] = 1.0;
break;

View File

@@ -40,6 +40,7 @@ Last Updated: 04/19/2025
// Defined in ENUMSTXT.H
extern char *Fldname[];
extern char *RptFlowUnitsTxt[];
extern char *PressUnitsTxt[];
extern void reindextanks(Project *pr);
@@ -103,7 +104,7 @@ void setdefaults(Project *pr)
pr->Warnflag = FALSE; // Warning flag is off
parser->Unitsflag = US; // US unit system
parser->Flowflag = GPM; // Flow units are gpm
parser->Pressflag = PSI; // Pressure units are psi
parser->Pressflag = DEFAULTUNIT; // Pressure units set based on unit system
out->Hydflag = SCRATCH; // No external hydraulics file
rpt->Tstatflag = SERIES; // Generate time series output
@@ -269,8 +270,11 @@ void adjustdata(Project *pr)
}
// Revise pressure units depending on flow units
if (parser->Unitsflag != SI) parser->Pressflag = PSI;
else if (parser->Pressflag == PSI) parser->Pressflag = METERS;
if (parser->Pressflag == DEFAULTUNIT)
{
if (parser->Unitsflag == SI) parser->Pressflag = METERS;
else parser->Pressflag = PSI;
}
// Store value of viscosity & diffusivity
ucf = 1.0;
@@ -404,8 +408,6 @@ void initunits(Project *pr)
strcpy(rpt->Field[DEMAND].Units, RptFlowUnitsTxt[parser->Flowflag]);
strcpy(rpt->Field[ELEV].Units, u_METERS);
strcpy(rpt->Field[HEAD].Units, u_METERS);
if (parser->Pressflag == METERS) strcpy(rpt->Field[PRESSURE].Units, u_METERS);
else strcpy(rpt->Field[PRESSURE].Units, u_KPA);
strcpy(rpt->Field[LENGTH].Units, u_METERS);
strcpy(rpt->Field[DIAM].Units, u_MMETERS);
strcpy(rpt->Field[FLOW].Units, RptFlowUnitsTxt[parser->Flowflag]);
@@ -423,8 +425,6 @@ void initunits(Project *pr)
if (parser->Flowflag == CMS) qcf = CMSperCFS;
hcf = MperFT;
if (parser->Pressflag == METERS) pcf = MperFT * hyd->SpGrav;
else pcf = KPAperPSI * PSIperFT * hyd->SpGrav;
wcf = KWperHP;
}
else // US units
@@ -448,10 +448,16 @@ void initunits(Project *pr)
if (parser->Flowflag == IMGD) qcf = IMGDperCFS;
if (parser->Flowflag == AFD) qcf = AFDperCFS;
hcf = 1.0;
pcf = PSIperFT * hyd->SpGrav;
wcf = 1.0;
}
strcpy(rpt->Field[PRESSURE].Units, PressUnitsTxt[parser->Pressflag]);
pcf = PSIperFT * hyd->SpGrav; // Default to PSI
if (parser->Pressflag == METERS) pcf = MperFT;
if (parser->Pressflag == KPA) pcf = KPAperPSI * PSIperFT * hyd->SpGrav;
if (parser->Pressflag == BAR) pcf = BARperPSI * PSIperFT * hyd->SpGrav;
if (parser->Pressflag == FEET) pcf = 1.0;
strcpy(rpt->Field[QUALITY].Units, "");
ccf = 1.0;
if (qual->Qualflag == CHEM)
@@ -508,7 +514,7 @@ void convertunits(Project *pr)
Parser *parser = &pr->parser;
int i, j, k;
double ucf; // Unit conversion factor
double ucf, ecf; // Unit conversion factor
Pdemand demand; // Pointer to demand record
Snode *node;
Stank *tank;
@@ -539,7 +545,9 @@ void convertunits(Project *pr)
hyd->Preq /= pr->Ucf[PRESSURE];
// Convert emitter discharge coeffs. to head loss coeff.
ucf = pow(pr->Ucf[FLOW], hyd->Qexp) / pr->Ucf[PRESSURE];
ecf = (parser->Unitsflag == US) ? (PSIperFT * hyd->SpGrav) : (MperFT);
ucf = pow(pr->Ucf[FLOW], hyd->Qexp) / ecf;
for (i = 1; i <= net->Njuncs; i++)
{
node = &net->Node[i];

View File

@@ -1926,6 +1926,8 @@ int optionchoice(Project *pr, int n)
else if (match(parser->Tok[1], w_PSI)) parser->Pressflag = PSI;
else if (match(parser->Tok[1], w_KPA)) parser->Pressflag = KPA;
else if (match(parser->Tok[1], w_METERS)) parser->Pressflag = METERS;
else if (match(parser->Tok[1], w_BAR)) parser->Pressflag = BAR;
else if (match(parser->Tok[1], w_FEET)) parser->Pressflag = FEET;
else return setError(parser, 1, 213);
}

View File

@@ -90,6 +90,8 @@
#define w_PSI "PSI"
#define w_KPA "KPA"
#define w_METERS "METERS"
#define w_BAR "BAR"
#define w_FEET "FEET"
#define w_ELEV "ELEV"
#define w_DEMAND "DEMAND"

View File

@@ -83,6 +83,7 @@ typedef int INT4;
#define MperFT 0.3048
#define PSIperFT 0.4333
#define KPAperPSI 6.895
#define BARperPSI 0.068948
#define KWperHP 0.7457
#define SECperDAY 86400
@@ -238,7 +239,10 @@ typedef enum {
typedef enum {
PSI, // pounds per square inch
KPA, // kiloPascals
METERS // meters
METERS, // meters
BAR, // bar
FEET, // feet
DEFAULTUNIT // default based on unit system (SI or US)
} PressureUnitsType;
typedef enum {

View File

@@ -74,13 +74,13 @@ BOOST_FIXTURE_TEST_CASE(test_pressure_units, FixtureInitClose)
BOOST_REQUIRE(error == 0);
BOOST_CHECK(units == EN_PSI);
// Change to pressure from PSI to meters and check it's still PSI
// Change to pressure from PSI to meters and check it is meters
error = EN_setoption(ph, EN_PRESS_UNITS, EN_METERS);
BOOST_REQUIRE(error == 0);
error = EN_getoption(ph, EN_PRESS_UNITS, &units);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(units == EN_PSI);
BOOST_CHECK(units == EN_METERS);
// Change flow units to LPS to change to metric units and rerun simulation
error = EN_setflowunits(ph, EN_LPS);
@@ -108,12 +108,12 @@ BOOST_FIXTURE_TEST_CASE(test_pressure_units, FixtureInitClose)
BOOST_REQUIRE(error == 0);
BOOST_CHECK(abs(p - 298.76035) < 1.e-5);
// Set pressure to PSI and check that it remains in kPa
// Set pressure to PSI and check that it has changed to PSI
error = EN_setoption(ph, EN_PRESS_UNITS, EN_PSI);
BOOST_REQUIRE(error == 0);
error = EN_getoption(ph, EN_PRESS_UNITS, &units);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(units == EN_KPA);
BOOST_CHECK(units == EN_PSI);
error = EN_closeH(ph);
BOOST_REQUIRE(error == 0);
@@ -264,4 +264,197 @@ BOOST_FIXTURE_TEST_CASE(test_rule_unit_change, FixtureOpenClose)
}
BOOST_FIXTURE_TEST_CASE(test_decoupled_pressure_units, FixtureInitClose)
{
int index;
long t;
double p, units;
// Create basic network
error = EN_addnode(ph, "R1", EN_RESERVOIR, &index);
BOOST_REQUIRE(error == 0);
error = EN_setnodevalue(ph, index, EN_ELEVATION, 100);
BOOST_REQUIRE(error == 0);
error = EN_addnode(ph, "J1", EN_JUNCTION, &index);
BOOST_REQUIRE(error == 0);
error = EN_addlink(ph, "P1", EN_PIPE, "R1", "J1", &index);
BOOST_REQUIRE(error == 0);
// Test 1: Start with US flow units (GPM) and change to PSI
error = EN_setflowunits(ph, EN_GPM);
BOOST_REQUIRE(error == 0);
// Should succeed in setting PSI pressure units
error = EN_setoption(ph, EN_PRESS_UNITS, EN_PSI);
BOOST_REQUIRE(error == 0);
error = EN_getoption(ph, EN_PRESS_UNITS, &units);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(units == EN_PSI);
// Test 2: With US flow units, set pressure to meters (should now work)
error = EN_setoption(ph, EN_PRESS_UNITS, EN_METERS);
BOOST_REQUIRE(error == 0);
error = EN_getoption(ph, EN_PRESS_UNITS, &units);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(units == EN_METERS);
// Test 3: With US flow units, set pressure to kPa (should now work)
error = EN_setoption(ph, EN_PRESS_UNITS, EN_KPA);
BOOST_REQUIRE(error == 0);
error = EN_getoption(ph, EN_PRESS_UNITS, &units);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(units == EN_KPA);
// Test 4: Change to SI flow units (LPS) but keep kPa pressure
error = EN_setflowunits(ph, EN_LPS);
BOOST_REQUIRE(error == 0);
// Pressure units should change to metric default of meters
error = EN_getoption(ph, EN_PRESS_UNITS, &units);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(units == EN_METERS);
// Test 5: With SI flow units, set pressure to PSI (should now work)
error = EN_setoption(ph, EN_PRESS_UNITS, EN_PSI);
BOOST_REQUIRE(error == 0);
error = EN_getoption(ph, EN_PRESS_UNITS, &units);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(units == EN_PSI);
// Test 6: Run simulation and check pressure values are correctly converted
error = EN_openH(ph);
BOOST_REQUIRE(error == 0);
error = EN_initH(ph, EN_NOSAVE);
BOOST_REQUIRE(error == 0);
error = EN_runH(ph, &t);
BOOST_REQUIRE(error == 0);
// Get pressure in PSI (should be ~43.33 PSI for 100 ft head)
error = EN_getnodevalue(ph, 1, EN_PRESSURE, &p);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(abs(p - 43.33) < 1.e-5);
// Change pressure units to meters during simulation
error = EN_setoption(ph, EN_PRESS_UNITS, EN_METERS);
BOOST_REQUIRE(error == 0);
// Pressure should now be in meters (~30.48 m for 100 ft head)
error = EN_getnodevalue(ph, 1, EN_PRESSURE, &p);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(abs(p - 30.48) < 1.e-5);
error = EN_closeH(ph);
BOOST_REQUIRE(error == 0);
}
BOOST_FIXTURE_TEST_CASE(test_automatic_pressure_unit_switching, FixtureInitClose)
{
int index;
double pressure_units;
// Create basic network
error = EN_addnode(ph, "R1", EN_RESERVOIR, &index);
BOOST_REQUIRE(error == 0);
error = EN_setnodevalue(ph, index, EN_ELEVATION, 100);
BOOST_REQUIRE(error == 0);
error = EN_addnode(ph, "J1", EN_JUNCTION, &index);
BOOST_REQUIRE(error == 0);
error = EN_addlink(ph, "P1", EN_PIPE, "R1", "J1", &index);
BOOST_REQUIRE(error == 0);
// Test 1: Start with US flow units (CFS) - should have PSI pressure units
error = EN_setflowunits(ph, EN_CFS);
BOOST_REQUIRE(error == 0);
error = EN_getoption(ph, EN_PRESS_UNITS, &pressure_units);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(pressure_units == EN_PSI);
// Test 2: Change from US flow units (CFS) to metric flow units (LPS)
// Pressure units should automatically change from PSI to METERS
error = EN_setflowunits(ph, EN_LPS);
BOOST_REQUIRE(error == 0);
error = EN_getoption(ph, EN_PRESS_UNITS, &pressure_units);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(pressure_units == EN_METERS);
// Test 3: Change from metric flow units (LPS) back to US flow units (GPM)
// Pressure units should automatically change from METERS to PSI
error = EN_setflowunits(ph, EN_GPM);
BOOST_REQUIRE(error == 0);
error = EN_getoption(ph, EN_PRESS_UNITS, &pressure_units);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(pressure_units == EN_PSI);
// Test 4: Change from US flow units (GPM) to another metric flow unit (MLD)
// Pressure units should automatically change from PSI to METERS
error = EN_setflowunits(ph, EN_MLD);
BOOST_REQUIRE(error == 0);
error = EN_getoption(ph, EN_PRESS_UNITS, &pressure_units);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(pressure_units == EN_METERS);
// Test 5: Manually set pressure units to kPa while using metric flow units
error = EN_setoption(ph, EN_PRESS_UNITS, EN_KPA);
BOOST_REQUIRE(error == 0);
error = EN_getoption(ph, EN_PRESS_UNITS, &pressure_units);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(pressure_units == EN_KPA);
// Test 6: Change from metric flow units (MLD) to US flow units (MGD)
// Pressure units should automatically change from kPa to PSI
error = EN_setflowunits(ph, EN_MGD);
BOOST_REQUIRE(error == 0);
error = EN_getoption(ph, EN_PRESS_UNITS, &pressure_units);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(pressure_units == EN_PSI);
// Test 7: Change from US flow units (MGD) to metric flow units (CMH)
// Pressure units should automatically change from PSI to METERS
error = EN_setflowunits(ph, EN_CMH);
BOOST_REQUIRE(error == 0);
error = EN_getoption(ph, EN_PRESS_UNITS, &pressure_units);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(pressure_units == EN_METERS);
// Test 8: Set pressure to kPa again with metric flow units
error = EN_setoption(ph, EN_PRESS_UNITS, EN_KPA);
BOOST_REQUIRE(error == 0);
// Test 9: Change between metric flow units (CMH to CMD)
// Pressure units should remain kPa (not changed to METERS since not switching from PSI)
error = EN_setflowunits(ph, EN_CMD);
BOOST_REQUIRE(error == 0);
error = EN_getoption(ph, EN_PRESS_UNITS, &pressure_units);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(pressure_units == EN_KPA);
// Test 10: Change from metric flow units (CMD) to US flow units (AFD)
// Pressure units should automatically change from kPa to PSI
error = EN_setflowunits(ph, EN_AFD);
BOOST_REQUIRE(error == 0);
error = EN_getoption(ph, EN_PRESS_UNITS, &pressure_units);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(pressure_units == EN_PSI);
// Test 11: Change between US flow units (AFD to IMGD)
// Pressure units should remain PSI
error = EN_setflowunits(ph, EN_IMGD);
BOOST_REQUIRE(error == 0);
error = EN_getoption(ph, EN_PRESS_UNITS, &pressure_units);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(pressure_units == EN_PSI);
// Test 12: Final test - metric flow units (CMS) should change PSI to METERS
error = EN_setflowunits(ph, EN_CMS);
BOOST_REQUIRE(error == 0);
error = EN_getoption(ph, EN_PRESS_UNITS, &pressure_units);
BOOST_REQUIRE(error == 0);
BOOST_CHECK(pressure_units == EN_METERS);
}
BOOST_AUTO_TEST_SUITE_END()