Support mixed pressure units

README.md

@@ -39,3 +39,7 @@ EN_getnodevalue(ph, index, EN_ELEVATION, &retrieved_elevation);
 - PR: [#594](https://github.com/OpenWaterAnalytics/EPANET/pull/609)
 - Issue: [#593](https://github.com/OpenWaterAnalytics/EPANET/issues/608)
 - Patch: [5b0a3f6](https://github.com/modelcreate/EPANET/commit/5b0a3f678e8c29fa47e8af70c8c1ea27f392273d)
+
+## Support for mixed pressure units
+
+- PR: [#864](https://github.com/OpenWaterAnalytics/EPANET/pull/862)

src/epanet.c

@@ -1376,7 +1376,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, xfactor, yfactor;
     double *Ucf = p->Ucf;
 
@@ -1388,6 +1388,7 @@ int DLLEXPORT EN_setflowunits(EN_Project p, int units)
     hfactor = Ucf[HEAD];
     efactor = Ucf[ELEV];
 
+    oldUnitFlag = p->parser.Unitsflag;
     p->parser.Flowflag = units;
     switch (units)
     {
@@ -1404,8 +1405,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 curves

src/input1.c

@@ -40,6 +40,7 @@ Last Updated: 07/08/2019
 // Defined in ENUMSTXT.H
 extern char *Fldname[];
 extern char *RptFlowUnitsTxt[];
+extern char *PressUnitsTxt[];
 
 int getdata(Project *pr)
 /*
@@ -96,7 +97,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
     parser->DefPat = 0;         // Default demand pattern index
     out->Hydflag = SCRATCH;     // No external hydraulics file
     rpt->Tstatflag = SERIES;    // Generate time series output
@@ -265,8 +266,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;
@@ -430,8 +434,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]);
@@ -448,16 +450,13 @@ void initunits(Project *pr)
         if (parser->Flowflag == CMD) qcf = CMDperCFS;
 
         hcf = MperFT;
-        if (parser->Pressflag == METERS) pcf = MperFT * hyd->SpGrav;
-        else pcf = KPAperPSI * PSIperFT * hyd->SpGrav;
         wcf = KWperHP;
     }
     else  // US units
     {
         strcpy(rpt->Field[DEMAND].Units, RptFlowUnitsTxt[parser->Flowflag]);
         strcpy(rpt->Field[ELEV].Units, u_FEET);
         strcpy(rpt->Field[HEAD].Units, u_FEET);
-        strcpy(rpt->Field[PRESSURE].Units, u_PSI);
         strcpy(rpt->Field[LENGTH].Units, u_FEET);
         strcpy(rpt->Field[DIAM].Units, u_INCHES);
         strcpy(rpt->Field[FLOW].Units, RptFlowUnitsTxt[parser->Flowflag]);
@@ -473,10 +472,14 @@ 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;
+
     strcpy(rpt->Field[QUALITY].Units, "");
     ccf = 1.0;
     if (qual->Qualflag == CHEM)
@@ -533,7 +536,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;
@@ -565,7 +568,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];

src/types.h

@@ -231,7 +231,8 @@ typedef enum {
 typedef enum {
   PSI,           // pounds per square inch
   KPA,           // kiloPascals
-  METERS         // meters
+  METERS,        // meters
+  DEFAULTUNIT    // default based on unit system (SI or US)
 } PressureUnitsType;
 
 typedef enum {