77 improve matrix class

src/simulator_core/solver/matrix/equation_object.py

@@ -38,6 +38,10 @@ def __init__(self) -> None:
         self.coefficients = np.array([], dtype=float)
         self.rhs = 0.0
 
+    def __len__(self) -> int:
+        """Return the number of coefficients in this equation."""
+        return len(self.coefficients)
+
     def to_list(self, length: int) -> list[float]:
         """Method to change the equation object to list which can be stored in the matrix.
 

src/simulator_core/solver/matrix/matrix.py

@@ -14,23 +14,25 @@
 #  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 """Module containing a matrix class to store the matrix and solve it using numpy."""
 import numpy as np
+import numpy.typing as npt
+import scipy as sp
 import csv
 from simulator_core.solver.matrix.equation_object import EquationObject
-from simulator_core.solver.matrix.utility import absolute_difference, relative_difference
+from simulator_core.solver.matrix.core_enum import NUMBER_CORE_QUANTITIES
 
 
 class Matrix:
     """Class which stores the matrix and can be used to solve it."""
 
+    num_unknowns: int = 0
+    sol_new: npt.NDArray = np.array([], dtype=float)
+    sol_old: npt.NDArray = np.array([], dtype=float)
+    relative_convergence: float = 1e-6
+    absolute_convergence: float = 1e-6
+
     def __init__(self) -> None:
         """Constructor of matrix class."""
-        self.num_unknowns: int = 0
-        self.mat: list[list[float]] = []
-        self.rhs: list[float] = []
-        self.sol_new: list[float] = []
-        self.sol_old: list[float] = []
-        self.relative_convergence: float = 1e-6
-        self.absolute_convergence: float = 1e-6
+        pass
 
     def add_unknowns(self, number_unknowns: int) -> int:
         """Method to add unknowns to the matrix.
@@ -42,58 +44,63 @@ def add_unknowns(self, number_unknowns: int) -> int:
         if number_unknowns < 1:
             raise ValueError("Number of unknowns should be at least 1.")
         self.num_unknowns += number_unknowns
-        self.sol_new += [1.0] * number_unknowns
-        self.sol_old += [0.0] * number_unknowns
-        return self.num_unknowns - number_unknowns
+        self.sol_new = np.concatenate([self.sol_new, np.ones(number_unknowns)])
+        self.sol_old = np.concatenate([self.sol_old, np.zeros(number_unknowns)])
 
-    def add_equation(self, equation_object: EquationObject) -> None:
-        """Method to add an equation to the matrix.
-
-        Add equation to the matrix it returns a unique id to the equation for easy access.
-        :param EquationObject equation_object: Object containing all information of the equation
-        :return:
-        """
-        row = len(self.rhs)
-        self.rhs.append(0.0)
-        self.mat.append([])
-        self.mat[row] = equation_object.to_list(self.num_unknowns)
-        self.rhs[row] = equation_object.rhs
+        return self.num_unknowns - number_unknowns
 
-    def solve(self, equations: list[EquationObject], dumb: bool = False) -> list[float]:
-        """Method to solve the system of equation given in the matrix.
+    def solve(self, equations: list[EquationObject], dump: bool = False) -> list[float]:
+        """Method to solve the system of equation given in the matrix using sparse matrix solver.
 
-        Solves the system of equations and returns the solution. The numpy linalg
-        library solve method is used. For this the matrix is converted to np arrays.
+        A sparse matrix solver is used, for this the coefficients, indices in the matrix are
+        converted to numpy arrays. These arrays are then used to create a csc_matrix which can
+        be solved by the sparse matrix solver of scipy.
+        :param dump: if true it will dump the matrix to a csv file
+        :param equations: list with the equations to solve.
         :return: list containing the solution of the system of equations.
         """
-        # TODO add checks if enough equations have been supplied.
-        # TODO check if matrix is solvable.
-        self.rhs = []
-        self.mat = []
-        for equation in equations:
-            self.add_equation(equation)
-        if dumb:
-            self.dump_matrix()
+        self.verify_equations(equations)
         self.sol_old = self.sol_new
-        a = np.array(self.mat)
-        b = np.array(self.rhs)
-        self.sol_new = np.linalg.solve(a, b).tolist()
-        return self.sol_new
+        coefficient_array = np.concatenate([equation.coefficients for equation in equations])
+        column_index_array = np.concatenate([equation.indices for equation in equations])
+        row_index_array = np.concatenate([np.full((len(equations[i])), i)
+                                          for i in range(len(equations))])
+        matrix = sp.sparse.csc_matrix((coefficient_array,
+                                       (row_index_array , column_index_array)),
+                                      shape=(self.num_unknowns, self.num_unknowns))
+        rhs = sp.sparse.csc_matrix([[equation.rhs] for equation in equations])
+        if dump:
+            self.dump_matrix(matrix=matrix, rhs_array=rhs)
+        self.sol_new = sp.sparse.linalg.spsolve(matrix, rhs)
+        if np.isnan(self.sol_new).any():
+            raise RuntimeError("Matrix is singular")
+        result: list[float] = self.sol_new.tolist()
+        return result
+
+    def verify_equations(self, equations: list[EquationObject]) -> None:
+        """Method to verify if the system of equations can be solved.
+
+        This method checks if the number of equations supplied is equal to the number of unknowns.
+        :param equations: list with the equations to verify.
+        :return: None
+        """
+        if len(equations) > self.num_unknowns:
+            raise RuntimeError(f"Too many equations supplied. Got {len(equations)} equations, "
+                               f"but number of unknowns is {self.num_unknowns}")
+        if len(equations) < self.num_unknowns:
+            raise RuntimeError(f"Not enough equation supplied. Got {len(equations)} equations, "
+                               f"but number of unknowns is {self.num_unknowns}")
 
     def is_converged(self) -> bool:
         """Returns true when the solution has converged and false when not.
 
-        This method uses both the absolute difference and the relative difference to determine if
-        the solution has converged. If one of them is converged true is returned.
+        This method uses the np.allclose method to calculate if the solution is converged
         :return: Bool whether the solution has converged based on the given convergence criteria.
         """
         if self.num_unknowns == 0:
             raise ValueError("No unknowns have been added to the matrix.")
-        rel_dif = max(
-            [relative_difference(old, new) for old, new in zip(self.sol_old, self.sol_new)])
-        abs_dif = max(
-            [absolute_difference(old, new) for old, new in zip(self.sol_old, self.sol_new)])
-        return rel_dif < self.relative_convergence or abs_dif < self.absolute_convergence
+        return np.allclose(self.sol_new, self.sol_old,
+                           atol=self.absolute_convergence, rtol=self.relative_convergence)
 
     def get_solution(self, index: int, number_of_unknowns: int) -> list[float]:
         """Method to get the solution of an asset.
@@ -106,20 +113,32 @@ def get_solution(self, index: int, number_of_unknowns: int) -> list[float]:
         be provided
         :return: list with teh solution
         """
-        #  TODO add checks if the index and number of unknowns are within the range of the solution.
-        return self.sol_new[index:index + number_of_unknowns]
-
-    def dump_matrix(self, file_name: str = 'dump.csv') -> None:
+        if index > self.num_unknowns:
+            raise IndexError(f"Index ({index}) greater than number "
+                             f"of unknowns ({self.num_unknowns})")
+        if (index + number_of_unknowns) > self.num_unknowns:
+            raise IndexError(f"Index ({index}) plus request number of elements "
+                             f"({number_of_unknowns}) greater than "
+                             f"number of unknowns {self.num_unknowns}")
+        result: list[float] = self.sol_new[index:index + number_of_unknowns].tolist()
+        return result
+
+    def dump_matrix(self, matrix: sp.sparse.csc_matrix, rhs_array: sp.sparse.csc_matrix,
+                    file_name: str = 'dump.csv') -> None:
         """Method to dump the matrix to a csv file.
 
-        :param str file_name: File name to dump the matrix in default=dump.csv
+        :param matrix: Matrix to be printed to the file
+        :param rhs_array: Right hand side to be printed to the file
+        :param file_name: File name to dump the matrix in default=dump.csv
         :return:
         """
         with open(file_name, 'w', newline='') as f:
             write = csv.writer(f)
-            for row, rhs in zip(self.mat, self.rhs):
-                write.writerow(row + [rhs])
+            write.writerow(
+                ['m', 'P', 'u'] * int(self.num_unknowns / NUMBER_CORE_QUANTITIES) + ['rhs'])
+            for row, rhs in zip(matrix.todense(), rhs_array.todense()):
+                write.writerow(row.tolist()[0] + rhs.tolist()[0])
 
     def reset_solution(self) -> None:
         """Method to reset the solution to 1, so the new iteration can start."""
-        self.sol_new = [1.0] * len(self.sol_new)
+        self.sol_new = np.ones(self.num_unknowns)

src/simulator_core/solver/network/network.py

@@ -14,6 +14,7 @@
 #  along with this program.  If not, see <https://www.gnu.org/licenses/>.
 """Module containing the network class."""
 import uuid
+import numpy.typing as npt
 
 from simulator_core.solver.network.assets.base_asset import BaseAsset
 from simulator_core.solver.network.assets.boundary import BaseBoundary
@@ -33,14 +34,15 @@ class Network:
         "Production": ProductionAsset,
         "Pipe": SolverPipe,
     }
+    assets: dict[uuid.UUID, BaseAsset]
+    nodes: dict[uuid.UUID, Node]
 
     def __init__(self) -> None:
         """Constructor of the network class.
 
         Initializes the class properties and loads the fluid properties.
         """
-        self.assets: dict[uuid.UUID, BaseAsset] = {}
-        self.nodes: dict[uuid.UUID, Node] = {}
+        self.assets = {}
 
     def add_asset(self, asset_type: str, name: uuid.UUID | None = None) -> uuid.UUID:
         """Method to add an asset to the network.
@@ -369,7 +371,7 @@ def check_connectivity(self) -> bool:
         """
         return self.check_connectivity_assets() and self.check_connectivity_nodes()
 
-    def set_result_asset(self, solution: list[float]) -> None:
+    def set_result_asset(self, solution: npt.NDArray) -> None:
         """Method to transfer the solution to the asset in the network.
 
         :param list[float] solution:Solution to be transferred to the assets.
@@ -378,9 +380,9 @@ def set_result_asset(self, solution: list[float]) -> None:
         for asset in self.assets:
             index = self.get_asset(asset_id=asset).matrix_index
             nou = self.get_asset(asset_id=asset).number_of_unknowns
-            self.get_asset(asset_id=asset).prev_sol = solution[index: index + nou]
+            self.get_asset(asset_id=asset).prev_sol = solution[index: index + nou].tolist()
 
-    def set_result_node(self, solution: list[float]) -> None:
+    def set_result_node(self, solution: npt.NDArray) -> None:
         """Method to transfer the solution to the nodes in the network.
 
         :param list[float] solution:Solution to be transferred to the nodes.
@@ -389,7 +391,7 @@ def set_result_node(self, solution: list[float]) -> None:
         for node in self.nodes:
             index = self.get_node(node_id=node).matrix_index
             nou = self.get_node(node_id=node).number_of_unknowns
-            self.get_node(node_id=node).prev_sol = solution[index: index + nou]
+            self.get_node(node_id=node).prev_sol = solution[index: index + nou].tolist()
 
     def print_result(self) -> None:
         """Method to print the result of the network."""

src/simulator_core/solver/solver.py

@@ -65,7 +65,7 @@ def solve(self) -> None:
         while not self.matrix.is_converged():
             iteration += 1
             equations = self.get_equations()
-            self.matrix.solve(equations, dumb=False)
+            self.matrix.solve(equations, dump=False)
             self.results_to_assets()
             if iteration > self._iteration_limit:
                 print("No converged solution reached")