From 4872eff732aad1f389dc2e316e4aec116a931696 Mon Sep 17 00:00:00 2001 From: "Marios S. Kyriakou" Date: Sun, 21 Apr 2024 22:26:32 +0300 Subject: [PATCH] fixes in getMSXOptions, setMSXOptions... - add testMSXoptions add example inp, msx --- epyt/__init__.py | 2 +- epyt/epanet.py | 189 +++++++++++-------------- epyt/networks/msx-examples/example.inp | 36 +++++ epyt/networks/msx-examples/example.msx | 58 ++++++++ epyt/tests/testMSXoptions.py | 73 ++++++++++ 5 files changed, 247 insertions(+), 111 deletions(-) create mode 100644 epyt/networks/msx-examples/example.inp create mode 100644 epyt/networks/msx-examples/example.msx create mode 100644 epyt/tests/testMSXoptions.py diff --git a/epyt/__init__.py b/epyt/__init__.py index 60ff64f..fb7613e 100644 --- a/epyt/__init__.py +++ b/epyt/__init__.py @@ -1,7 +1,7 @@ # -*- coding: utf-8 -*- __author__ = """Marios S. Kyriakou""" __email__ = "kiriakou.marios@ucy.ac.cy" -__version__ = "1.1.1" +__version__ = "1.1.2" __msxversion__ = "2.0.0" __copyright__ = """Copyright 2022, KIOS Research and Innovation Center of Excellence (KIOS CoE), University of Cyprus (www.kios.org.cy).""" diff --git a/epyt/epanet.py b/epyt/epanet.py index 5e964fd..7155e24 100644 --- a/epyt/epanet.py +++ b/epyt/epanet.py @@ -62,6 +62,7 @@ from inspect import getmembers, isfunction, currentframe, getframeinfo from ctypes import cdll, byref, create_string_buffer, c_uint64, c_uint32, c_void_p, c_int, c_double, c_float, c_long, \ c_char_p +from types import SimpleNamespace import matplotlib.pyplot as plt from datetime import datetime from epyt import __version__, __msxversion__ @@ -11196,6 +11197,17 @@ def loadMSXFile(self, msxname, customMSXlib=None, ignore_properties=False): d = epanet(inpname, msx=True,customlib=epanetlib) d.loadMSXFile(msxname,customMSXlib=msxlib)""" + if not os.path.exists(msxname): + for root, dirs, files in os.walk(resource_filename("epyt", "")): + for name in files: + if name.lower().endswith(".msx"): + if name == msxname: + msxname = os.path.join(root, msxname) + break + else: + continue + break + self.msxname = msxname[:-4] + '_temp.msx' copyfile(msxname, self.msxname) self.msx = epanetmsxapi(self.msxname, customMSXlib=customMSXlib) @@ -11422,66 +11434,58 @@ def getMSXError(self, code): d.getMSXError(510)""" self.msx.MSXgeterror(code) - def getMSXOptions(self, param="", getall=False): + def getMSXOptions(self): """ Retrieves all the options. Example: d=epanet('net2-cl2.inp') d.loadMSXFile('net2-cl2.msx') d.getMSXOptions()""" - msxname = self.msxname - options = {} - options["AREA_UNITS"] = "FT2" - options["RATE_UNITS"] = "HR" - options["SOLVER"] = "EUL" - options["TIMESTEP"] = 300 - options["ATOL"] = 0.01 - options["RTOL"] = 0.001 - options["COUPLING"] = "NONE" - options["COMPILER"] = "NONE" + # AREA_UNITS FT2/M2/CM2 + # RATE_UNITS SEC/MIN/HR/DAY + # SOLVER EUL/RK5/ROS2 + # COUPLING FULL/NONE + # TIMESTEP seconds + # ATOL value + # RTOL value + # COMPILER NONE/VC/GC + # SEGMENTS value + # PECLET value try: - with open(msxname, 'r') as f: - sect = 0 - for line in f: - if line.startswith("[OPTIONS]"): - sect = 1 - continue - elif line.startswith("[END") or line.startswith("[REPORTS]"): - break - elif sect == 1: - if not line.strip(): - continue - if line.startswith("["): - return options - key, value = line.split(None, 1) - if key == param or not param: - if key == "TIMESTEP": - options["TIMESTEP"] = float(value) - elif key == "AREA_UNITS": - options["AREA_UNITS"] = value.strip() - elif key == "RATE_UNITS": - options["RATE_UNITS"] = value.strip() - elif key == "SOLVER": - options["SOLVER"] = value.strip() - elif key == "RTOL": - options["RTOL"] = float(value) - elif key == "ATOL": - options["ATOL"] = float(value) - elif key == "COUPLING": - options["COUPLING"] = value.strip() - elif key == "COMPILER": - options["COMPILER"] = value.strip() + # Key-value pairs to search for + keys = ["AREA_UNITS", "RATE_UNITS", "SOLVER", "COUPLING", "TIMESTEP", "ATOL", "RTOL", "COMPILER", "SEGMENTS", \ + "PECLET"] + float_values = ["TIMESTEP", "ATOL", "RTOL", "SEGMENTS", "PECLET"] + values = {key: None for key in keys} + + # Flag to determine if we're in the [OPTIONS] section + in_options = False + + # Open and read the file + with open(self.msxname, 'r') as file: + for line in file: + # Check for [OPTIONS] section + if "[OPTIONS]" in line: + in_options = True + elif "[" in line and "]" in line: + in_options = False # We've reached a new section + + if in_options: + # Pattern to match the keys and extract values, ignoring comments and whitespace + pattern = re.compile(r'^\s*(' + '|'.join(keys) + r')\s+(.*?)\s*(?:;.*)?$') + match = pattern.search(line) + if match: + key, value = match.groups() + if key in float_values: + values[key] = float(value) else: - options[key] = value - - if not getall and param: - return options[param] + values[key] = value + return SimpleNamespace(**values) except FileNotFoundError: warnings.warn("Please load MSX File.") return {} - return options def getMSXTimeStep(self): """ Retrieves the time step. @@ -11492,7 +11496,7 @@ def getMSXTimeStep(self): d.getMSXTimeStep() See also setMSXTimeStep.""" - return self.getMSXOptions("TIMESTEP") + return self.getMSXOptions().TIMESTEP def getMSXRateUnits(self): """ Retrieves rate units. @@ -11502,7 +11506,7 @@ def getMSXRateUnits(self): d.getMSXRateUnits() See also setMSXRateUnits.""" - return self.getMSXOptions("RATE_UNITS") + return self.getMSXOptions().RATE_UNITS def getMSXAreaUnits(self): """ Retrieves Are units. @@ -11512,7 +11516,7 @@ def getMSXAreaUnits(self): d.getMSXAreaUnits() See also setMSXAreaUnits.""" - return self.getMSXOptions("AREA_UNITS") + return self.getMSXOptions().AREA_UNITS def getMSXCompiler(self): """ Retrieves the chemistry function compiler code. @@ -11529,7 +11533,7 @@ def getMSXCompiler(self): See also setMSXCompilerNONE, setMSXCompilerVC, setMSXCompilerGC.""" - return self.getMSXOptions("COMPILER") + return self.getMSXOptions().COMPILER def getMSXCoupling(self): """ Retrieves the degree of coupling for solving DAE's. @@ -11547,7 +11551,7 @@ def getMSXCoupling(self): d.getMSXCoupling() See also setMSXCouplingFULL, setMSXCouplingNONE.""" - return self.getMSXOptions("COUPLING") + return self.getMSXOptions().COUPLING def getMSXSolver(self): """ Retrieves the solver method. @@ -11563,7 +11567,7 @@ def getMSXSolver(self): d.getMSXSolver() See also setMSXSolverEUL, setMSXSolverRK5, setMSXSolverROS2.""" - return self.getMSXOptions("SOLVER") + return self.getMSXOptions().SOLVER def getMSXAtol(self): """ Retrieves the absolute tolerance. @@ -11574,7 +11578,7 @@ def getMSXAtol(self): d.getMSXAtol() See also getMSXRtol.""" - return self.getMSXOptions("ATOL") + return self.getMSXOptions().ATOL def getMSXRtol(self): """ Retrieves the relative accuracy level. @@ -11585,7 +11589,7 @@ def getMSXRtol(self): d.getMSXRtol() See also getMSXAtol.""" - return self.getMSXOptions("RTOL") + return self.getMSXOptions().RTOL def getMSXConstantsNameID(self, varagin=None): """ Retrieves the constant's ID. @@ -12340,55 +12344,22 @@ def getMSXSourceNodeNameID(self): nodes.append(i) return nodes - def setMSXOptions(self, *args): - - for i in range(len(args) // 2): - argument = args[2 * i].lower() - if argument == 'areaunits': - self.areaunits = args[2 * i + 1] - self.changeMSXOptions("AREA_UNITS", args[2 * i + 1]) - elif argument == 'rateunits': - self.rateunits = args[2 * i + 1] - self.changeMSXOptions("RATE_UNITS", args[2 * i + 1]) - elif argument == 'solver': - self.solver = args[2 * i + 1] - self.changeMSXOptions("SOLVER", args[2 * i + 1]) - elif argument == 'timestep': - self.timestep = args[2 * i + 1] - self.changeMSXOptions("TIMESTEP", args[2 * i + 1]) - elif argument == 'atol': - self.atol = args[2 * i + 1] - self.changeMSXOptions("ATOL", args[2 * i + 1]) - elif argument == 'rtol': - self.rtol = args[2 * i + 1] - self.changeMSXOptions("RTOL", args[2 * i + 1]) - elif argument == 'coupling': - self.coupling = args[2 * i + 1] - self.changeMSXOptions("COUPLING", args[2 * i + 1]) - elif argument == 'compiler': - self.compiler = args[2 * i + 1] - self.changeMSXOptions("COMPILER", args[2 * i + 1]) - else: - print('Invalid property found.') - return - def changeMSXOptions(self, param, change): - msxname = self.msxname - f = open(msxname, 'r+') - lines = f.readlines() - flag = 0 - for i, line in enumerate(lines): - if line.strip() == '[OPTIONS]': - options_index = i - if line.startswith(param): - lines[i] = param + "\t" + str(change) + "\n" - flag = 1 - if flag == 0: - lines = list(lines) - lines.insert(options_index + 1, param + "\t" + str(change) + "\n") - f.seek(0) - f.writelines(lines) - f.close() + with open(self.msxname, 'r+') as f: + lines = f.readlines() + options_index = -1 # Default to -1 in case the [OPTIONS] section does not exist + flag = 0 + for i, line in enumerate(lines): + if line.strip() == '[OPTIONS]': + options_index = i + elif line.strip().startswith(param): + lines[i] = param + "\t" + str(change) + "\n" + flag = 1 + if flag == 0 and options_index != -1: + lines.insert(options_index + 1, param + "\t" + str(change) + "\n") + f.seek(0) + f.writelines(lines) + f.truncate() def setMSXAreaUnitsCM2(self): """ Sets the area units to square centimeters. @@ -12784,7 +12755,7 @@ def getMSXComputedQualitySpecie(self, species=None): species_index_name = self.getMSXSpeciesIndex(species) node_count = self.getNodeCount() - link_count = self.getNodeCount() + link_count = self.getLinkCount() node_indices = list(range(1, node_count + 1)) link_indices = list(range(1, link_count + 1)) # Initialized quality and time @@ -13988,8 +13959,6 @@ def ENgetnodevalue(self, index, code_p): OWA-EPANET Toolkit: http://wateranalytics.org/EPANET/group___nodes.html """ - fValue = None - if self._ph is not None: fValue = c_double() self.errcode = self._lib.EN_getnodevalue(self._ph, int(index), code_p, byref(fValue)) @@ -15980,10 +15949,10 @@ def MSXgetinitqual(self, obj_type, index, species): Returns: value: the initial concetration of the species at the node or link of interest.""" - # obj_type = c_int(obj_type) value = c_double() - # species = c_int(species) - # index = c_int(index) + obj_type = c_int(obj_type) + species = c_int(species) + index = c_int(index) err = self.msx_lib.MSXgetinitqual(obj_type, index, species, byref(value)) if err: Warning(self.MSXerror(err)) @@ -16311,4 +16280,4 @@ def MSXgeterror(self, err): if e: Warning(errmsg.value.decode()) - return errmsg.value.decode() + return errmsg.value.decode() \ No newline at end of file diff --git a/epyt/networks/msx-examples/example.inp b/epyt/networks/msx-examples/example.inp new file mode 100644 index 0000000..0887e00 --- /dev/null +++ b/epyt/networks/msx-examples/example.inp @@ -0,0 +1,36 @@ +[TITLE] +EPANET-MSX Example Network + +[JUNCTIONS] +;ID Elev Demand Pattern + A 0 4.1 + B 0 3.4 + C 0 5.5 + D 0 2.3 + +[RESERVOIRS] +;ID Head Pattern + Source 100 + +[PIPES] +;ID Node1 Node2 Length Diameter Roughness + 1 Source A 1000 200 100 + 2 A B 800 150 100 + 3 A C 1200 200 100 + 4 B C 1000 150 100 + 5 C D 2000 150 100 + +[TIMES] + Duration 48 + Hydraulic Timestep 1:00 + Quality Timestep 0:05 + Report Timestep 2 + Report Start 0 + Statistic NONE + +[OPTIONS] + Units CMH + Headloss H-W + Quality NONE + +[END] diff --git a/epyt/networks/msx-examples/example.msx b/epyt/networks/msx-examples/example.msx new file mode 100644 index 0000000..b20e085 --- /dev/null +++ b/epyt/networks/msx-examples/example.msx @@ -0,0 +1,58 @@ +[TITLE] +Arsenic Oxidation/Adsorption Example + +[OPTIONS] + AREA_UNITS M2 ;Surface concentration is mass/m2 + RATE_UNITS HR ;Reaction rates are concentration/hour + SOLVER RK5 ;5-th order Runge-Kutta integrator + TIMESTEP 360 ;360 sec (5 min) solution time step + RTOL 0.001 ;Relative concentration tolerance + ATOL 0.0001 ;Absolute concentration tolerance + +[SPECIES] + BULK AS3 UG ;Dissolved arsenite + BULK AS5 UG ;Dissolved arsenate + BULK AStot UG ;Total dissolved arsenic + WALL AS5s UG ;Adsorbed arsenate + BULK NH2CL MG ;Monochloramine + +[COEFFICIENTS] + CONSTANT Ka 10.0 ;Arsenite oxidation rate coefficient + CONSTANT Kb 0.1 ;Monochloramine decay rate coefficient + CONSTANT K1 5.0 ;Arsenate adsorption coefficient + CONSTANT K2 1.0 ;Arsenate desorption coefficient + CONSTANT Smax 50 ;Arsenate adsorption saturation limit + +[TERMS] + Ks K1/K2 ;Equil. adsorption coeff. + +[PIPES] + ;Arsenite oxidation + RATE AS3 -Ka*AS3*NH2CL + ;Arsenate production + RATE AS5 Ka*AS3*NH2CL - Av*(K1*(Smax-AS5s)*AS5 - K2*AS5s) + ;Monochloramine decay + RATE NH2CL -Kb*NH2CL + ;Arsenate adsorption + EQUIL AS5s Ks*Smax*AS5/(1+Ks*AS5) - AS5s + ;Total bulk arsenic + FORMULA AStot AS3 + AS5 + +[TANKS] + RATE AS3 -Ka*AS3*NH2CL + RATE AS5 Ka*AS3*NH2CL + RATE NH2CL -Kb*NH2CL + FORMULA AStot AS3 + AS5 + +[QUALITY] + ;Initial conditions (= 0 if not specified here) + NODE Source AS3 10.0 + NODE Source NH2CL 2.5 + +[REPORT] + NODES C D ;Report results for nodes C and D + LINKS 5 ;Report results for pipe 5 + SPECIES AStot YES ;Report results for each specie + SPECIES AS5 YES + SPECIES AS5s YES + SPECIES NH2CL YES diff --git a/epyt/tests/testMSXoptions.py b/epyt/tests/testMSXoptions.py new file mode 100644 index 0000000..01a739e --- /dev/null +++ b/epyt/tests/testMSXoptions.py @@ -0,0 +1,73 @@ +from epyt import epanet, networks +import os + +d = epanet('example.inp') +d.loadMSXFile('example.msx') + +a = d.getMSXOptions() +# AREA_UNITS FT2/M2/CM2 +# RATE_UNITS SEC/MIN/HR/DAY +# SOLVER EUL/RK5/ROS2 +# COUPLING FULL/NONE +# TIMESTEP seconds +# ATOL value +# RTOL value +# COMPILER NONE/VC/GC +# SEGMENTS value +# PECLET value +d.printv(d.getMSXTimeStep()) +d.printv(d.getMSXRateUnits()) +d.printv(d.getMSXAreaUnits()) +d.printv(d.getMSXCompiler()) +d.printv(d.getMSXCoupling()) +d.printv(d.getMSXSolver()) +d.printv(d.getMSXAtol()) +d.printv(d.getMSXRtol()) +# d.printv(d.getMSXSegments()) +# d.printv(d.getMSXPeclet()) + +print('after set') +d.setMSXAreaUnitsCM2() +d.printv(d.getMSXAreaUnits()) + +d.setMSXAreaUnitsFT2() +d.printv(d.getMSXAreaUnits()) + +d.setMSXAreaUnitsM2() +d.printv(d.getMSXAreaUnits()) + +d.setMSXAtol(0.1) +d.printv(d.getMSXAtol()) + +d.setMSXRtol(0.2) +d.printv(d.getMSXRtol()) + +d.setMSXCompilerGC() +d.printv(d.getMSXCompiler()) + +d.setMSXCompilerVC() +d.printv(d.getMSXCompiler()) + +d.setMSXCompilerNONE() +d.printv(d.getMSXCompiler()) + +d.setMSXCouplingFULL() +d.printv(d.getMSXCoupling()) + +d.setMSXCouplingNONE() +d.printv(d.getMSXCoupling()) + +d.setMSXTimeStep(100) +d.printv(d.getMSXTimeStep()) + +d.setMSXRateUnitsSEC() +d.printv(d.getMSXRateUnits()) + +d.setMSXRateUnitsMIN() +d.printv(d.getMSXRateUnits()) + +d.setMSXRateUnitsHR() +d.printv(d.getMSXRateUnits()) + +d.setMSXRateUnitsDAY() +d.printv(d.getMSXRateUnits()) \ No newline at end of file