fixes in getMSXOptions and setMSXOptions .- add test_msx_options.

add example.inp and example.msx

epyt/__init__.py

@@ -1,7 +1,7 @@
 # -*- coding: utf-8 -*-
 __author__ = """Marios S. Kyriakou"""
 __email__ = "[email protected]"
-__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)."""

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__
@@ -11422,66 +11423,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 +11485,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 +11495,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 +11505,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 +11522,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 +11540,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 +11556,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 +11567,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 +11578,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 +12333,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 +12744,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

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]

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