Bug fix in roughness unit conversion when using D-W (#450)

* Added D-W conversion for roughness coefficient for inp file i/o

* bug fix in roughness conversion

* Added a warning if the user changes headloss to/from D-W

* Added additional tests for changing headloss formula

* Updated test

* commented out graphics in test

wntr/epanet/io.py

@@ -674,6 +674,8 @@ def _write_tanks(self, f, wn, version=2.2):
         f.write('\n'.encode(sys_default_enc))
 
     def _read_pipes(self):
+        darcy_weisbach = self.wn.options.hydraulic.headloss == "D-W"
+
         for lnum, line in self.sections['[PIPES]']:
             line = line.split(';')[0]
             current = line.split()
@@ -701,7 +703,8 @@ def _read_pipes(self):
                         current[2],
                         to_si(self.flow_units, float(current[3]), HydParam.Length),
                         to_si(self.flow_units, float(current[4]), HydParam.PipeDiameter),
-                        float(current[5]),
+                        to_si(self.flow_units, float(current[5]), HydParam.RoughnessCoeff, 
+                              darcy_weisbach=darcy_weisbach),
                         minor_loss,
                         link_status,
                         check_valve)
@@ -711,6 +714,8 @@ def _read_pipes(self):
                 raise ENValueError(211, str(e.args[0]), line_num=lnum) from e
 
     def _write_pipes(self, f, wn):
+        darcy_weisbach = wn.options.hydraulic.headloss == "D-W"
+
         f.write('[PIPES]\n'.encode(sys_default_enc))
         f.write(_PIPE_LABEL.format(';ID', 'Node1', 'Node2', 'Length', 'Diameter',
                                    'Roughness', 'Minor Loss', 'Status').encode(sys_default_enc))
@@ -723,7 +728,9 @@ def _write_pipes(self, f, wn):
                  'node2': pipe.end_node_name,
                  'len': from_si(self.flow_units, pipe.length, HydParam.Length),
                  'diam': from_si(self.flow_units, pipe.diameter, HydParam.PipeDiameter),
-                 'rough': pipe.roughness,
+                 'rough': from_si(self.flow_units, pipe.roughness, 
+                                  HydParam.RoughnessCoeff, 
+                                  darcy_weisbach=darcy_weisbach),
                  'mloss': pipe.minor_loss,
                  'status': str(pipe.initial_status),
                  'com': ';'}

wntr/epanet/util.py

@@ -571,7 +571,7 @@ def _to_si(self, flow_units, data, darcy_weisbach=False):
 
         elif self in [HydParam.RoughnessCoeff] and darcy_weisbach:
             if flow_units.is_traditional:
-                data = data * (1000.0 * 0.3048)  # 1e-3 ft to m
+                data = data * (0.001 * 0.3048)  # 1e-3 ft to m
             elif flow_units.is_metric:
                 data = data * 0.001  # mm to m
 
@@ -668,7 +668,7 @@ def _from_si(self, flow_units, data, darcy_weisbach=False):
 
         elif self in [HydParam.RoughnessCoeff] and darcy_weisbach:
             if flow_units.is_traditional:
-                data = data / (1000.0 * 0.3048)  # 1e-3 ft from m
+                data = data / (0.001 * 0.3048)  # 1e-3 ft from m
             elif flow_units.is_metric:
                 data = data / 0.001  # mm from m
 

wntr/network/options.py

@@ -10,6 +10,7 @@
 """
 import re
 import logging
+import warnings
 import copy
 
 logger = logging.getLogger(__name__)
@@ -406,6 +407,19 @@ def __setattr__(self, name, value):
             value = str.upper(value)
             if value not in ['H-W', 'D-W', 'C-M']:
                 raise ValueError('headloss must be one of "H-W", "D-W", or "C-M"')
+            # If headloss is changed from ['H-W', 'C-M'] to/from 'D-W', print
+            # a warning, the units of the roughness coefficient cannot be 
+            # converted from unitless to length (0.001 ft or mm)
+            try:
+                orig_value = self.__dict__[name]
+            except:
+                orig_value = None
+            if orig_value is not None:
+                if (orig_value in ['H-W', 'C-M'] and value == 'D-W') or \
+                   (value in ['H-W', 'C-M'] and orig_value == 'D-W'):
+                   warnings.warn('Changing the headloss formula from ' +
+                                 orig_value + ' to ' + value +
+                                 ' will not change the units of the roughness coefficient.')
         elif name == 'hydraulics':
             if value is not None:
                 value = str.upper(value)

wntr/tests/test_epanet_io.py

@@ -625,7 +625,9 @@ def setUpClass(self):
     def tearDownClass(self):
         pass
 
-    def test_link_flowrate_units_convert(self):
+    def test_units_convert(self):
+        # Compares Net3 EpanetSimulator flowrate results using INP files saved 
+        # using GPM and CMH units
         for link_name, link in self.wn.links():
             for t in self.results2.link["flowrate"].index:
                 self.assertLessEqual(
@@ -636,7 +638,7 @@ def test_link_flowrate_units_convert(self):
                     0.00001,
                 )
 
-    def test_link_headloss_units_convert(self):
+    def test_link_headloss_units(self):
 
         # headloss = per unit length for pipes and CVs
         pipe_name = '123'
@@ -665,7 +667,78 @@ def test_link_headloss_units_convert(self):
             ),
             0.00001,
         )
-
 
+    def test_headloss_formula_roughness_units(self):
+        """
+        See Table 3.2 Roughness Coefficients, this test uses values for Plastic Pipes
+        https://epanet22.readthedocs.io/en/latest/3_network_model.html?highlight=roughness#id3
+        """
+        pressure = {}
+        for headloss in ['H-W', 'C-M', 'D-W']:
+            for units in ['GPM', 'LPS']:
+                file_prefix = 'temp_'+headloss+'_'+units
+
+                inp_file = join(ex_datadir, "Net3.inp")
+                wn = self.wntr.network.WaterNetworkModel(inp_file)
+                wn.options.hydraulic.demand_model = 'DDA'
+                wn.options.hydraulic.inpfile_units = units
+                wn.options.time.duration = 0 # steady state
+
+                wn.options.hydraulic.headloss = headloss
+                if headloss == 'H-W':
+                    roughness = 145 # unitless
+                elif headloss == 'D-W':
+                    roughness = 0.005 # 0.001*ft
+                    roughness = (roughness*0.001)*0.3048 #  1.524e-6 m == 0.001524 mm
+                else: # C-M
+                    roughness = 0.011 # unitless
+
+                for name, pipe in wn.pipes():
+                    pipe.roughness = roughness
+
+                sim = self.wntr.sim.EpanetSimulator(wn)
+                results = sim.run_sim(file_prefix=file_prefix)
+                temp = results.node['pressure'].loc[0,wn.junction_name_list]
+                pressure[headloss+', '+units] = temp
+
+        #import matplotlib.pylab as plt
+        #import pandas as pd
+        #plt.figure()
+        #pd.DataFrame(pressure).plot.bar()
+        #plt.legend() #loc='lower right')
+
+        ## Compares all results to H-W, GPM
+        threshold = 0.55
+        for key in pressure.keys():
+            MAE = (pressure['H-W, GPM'] - pressure[key]).abs().mean()
+            print(key, MAE)
+            self.assertLessEqual(MAE, threshold) # m
+
+        ## Changing the headloss formula from H-W to D-W, without changing roughness,
+        # will result in errors
+        inp_file = join(ex_datadir, "Net3.inp")
+        wn = self.wntr.network.WaterNetworkModel(inp_file)
+        assert wn.options.hydraulic.headloss == 'H-W'
+        assert wn.options.hydraulic.inpfile_units == 'GPM'
+        pressure_hw = pressure['H-W, GPM']
+        # Change to D-W without changing roughness
+        wn.options.hydraulic.headloss = 'D-W'
+        wn.options.hydraulic.demand_model = 'DDA'
+        wn.options.hydraulic.inpfile_units = units
+        wn.options.time.duration = 0 # steady state
+        sim = self.wntr.sim.EpanetSimulator(wn)
+        results = sim.run_sim()
+        pressure_dw = results.node['pressure'].loc[0,wn.junction_name_list]
+        # Compare results
+        MAE = (pressure_hw - pressure_dw).abs().mean()
+        with self.assertRaises(AssertionError):
+            self.assertLessEqual(MAE, threshold) # m
+
+        ## D-W is not supported by the WNTRSimulator
+        sim = self.wntr.sim.WNTRSimulator(wn)
+        with self.assertRaises(NotImplementedError):
+            results = sim.run_sim()
+
+
 if __name__ == "__main__":
     unittest.main()