Moved to single equation model and system uses 2D numpy array

splitfxm/model.py

@@ -3,15 +3,15 @@ def __init__(self):
         """
         Initialize a Model object.
         """
-        self._equations = []
+        self._equation = None
 
-    def equations(self):
+    def equation(self):
         """
-        Get the list of equations associated with this Model.
+        Get the equation associated with this model
 
         Returns
         -------
         list
             The list of equations
         """
-        return self._equations
+        return self._equation

splitfxm/system.py

@@ -50,18 +50,21 @@ def residuals(self, d: Domain):
         nb_left = d.nb(btype.LEFT)
         nb_right = d.nb(btype.RIGHT)
 
+        # Fetch equation for model
+        eq = self._model.equation()
+
         rhs_list = []
 
         for i in range(ilo, ihi + 1):
-            rhs = np.array([])
-            # Append residuals from each equation
-            for eq in self._model.equations():
-                # Send two-sided stencil
-                # Let model decide computation
-                cell_sub = [cells[i + offset]
-                            for offset in range(-nb_left, nb_right + 1)]
-                rhs = np.concatenate(
-                    (rhs, eq.residuals(cell_sub, self._scheme, self._scheme_opts)))
+            # Define the neighborhood and band around the current cell
+            cell_sub = [cells[i + offset]
+                        for offset in range(-nb_left, nb_right + 1)]
+            # Send two-sided stencil
+            # Let model decide computation
+            rhs = eq.residuals(cell_sub, self._scheme, self._scheme_opts)
             rhs_list.append(rhs)
 
-        return rhs_list
+        # Stack the residuals to form a 2D array
+        rhs_array = np.vstack(rhs_list)
+
+        return rhs_array