diff --git a/wntr/epanet/io.py b/wntr/epanet/io.py index fea4e3651..af56b8bed 100644 --- a/wntr/epanet/io.py +++ b/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': ';'} diff --git a/wntr/epanet/util.py b/wntr/epanet/util.py index c06528573..c09e51832 100644 --- a/wntr/epanet/util.py +++ b/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 diff --git a/wntr/network/options.py b/wntr/network/options.py index 3d1d086a5..120cb39dc 100644 --- a/wntr/network/options.py +++ b/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) diff --git a/wntr/tests/test_epanet_io.py b/wntr/tests/test_epanet_io.py index 38d611734..8414ca127 100644 --- a/wntr/tests/test_epanet_io.py +++ b/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()