Move residual calculation to separate module

CADETPythonSimulator/exception.py

@@ -1,3 +1,8 @@
-class NotInitializedError(Exception):
-    """Exception raised when a unit operation is not yet initialized."""
+
+class CADETPythonSimError(Exception):
+    """Typical Exception for Error Handling"""
     pass
+
+class NotInitializedError(CADETPythonSimError):
+    """Exception raised when a unit operation is not yet initialized."""
+    pass

CADETPythonSimulator/residual.py

@@ -0,0 +1,84 @@
+import numpy as np
+from CADETPythonSimulator.exception import CADETPythonSimError
+import warnings
+
+def calculate_residual_volume_cstr(
+        V : float,
+        V_dot : float,
+        Q_in : float ,
+        Q_out : float
+        ) -> float:
+    """
+    Calculates the residual equations of the volume of a cstr.
+
+    Parameters
+    ----------
+    V : float
+        Volume within the CSTR
+    V_dot : float
+        Volume change rate of the CSTR
+    Q_in : float
+        Volume entering the Unit
+    Q_out : float
+        Volume leaving the Unit
+    Returns
+    -------
+    float 
+        Residual of the Flow equation of the CSTR with dimensions like the inpu
+    """
+
+    if V < 0:
+        raise CADETPythonSimError("V can't be less then zero")
+
+    return V_dot - Q_in + Q_out
+
+def calculate_residual_concentration_cstr(
+        c : np.ndarray,
+        c_dot : np.ndarray,
+        V : float,
+        V_dot : float,
+        Q_in : float,
+        Q_out : float,
+        c_in : np.ndarray
+        ) -> np.ndarray :
+    """
+    Calculates the residual equations of the concentration of a cstr
+    
+    Parameters
+    ----------
+    c : np.ndarray
+        Concentration
+    c_dot : np.ndarray
+        Changing of the concentration
+    V : float
+        Volume within the CSTR
+    V_dot : float
+        Volume change rate of the CSTR
+    Q_in : float
+        Volume entering the Unit
+    Q_out : float
+        Volume leaving the Unit
+    c_in : np.ndarray
+        Initial concentration
+    """
+    if V < 0:
+        raise CADETPythonSimError("V can't be less then zero")
+
+
+    return c_dot * V + V_dot * c - Q_in * c_in + Q_out * c
+
+
+def calculate_residuals_visc_cstr():
+    """
+    Calculates the residual of the Viscosity equation of the CSTR
+    """
+    warnings.warn("Viscosity of CSTR not yet implemented")
+
+    return 0
+
+
+def calculate_residual_def():
+    """
+    Calculates the residual equations fo a dead end filtration equation.
+    """
+    raise NotImplementedError

CADETPythonSimulator/unit_operation.py

@@ -11,10 +11,14 @@
 )
 from CADETProcess.dynamicEvents import Section
 
-from CADETPythonSimulator.exception import NotInitializedError
+from CADETPythonSimulator.exception import NotInitializedError, CADETPythonSimError
 from CADETPythonSimulator.state import State, state_factory
+from CADETPythonSimulator.residual import (
+    calculate_residual_volume_cstr, calculate_residual_concentration_cstr, calculate_residuals_visc_cstr
+    )
 from CADETPythonSimulator.rejection import RejectionBase
 from CADETPythonSimulator.cake_compressibility import CakeCompressibilityBase
+from CADETPythonSimulator.viscosity import LogarithmicMixingViscosity, ViscosityBase
 
 
 class UnitOperationBase(Structure):
@@ -475,9 +479,9 @@ def compute_residual(
 
         """
         # Inlet DOFs are simply copied to the residual.
-        for i in range(self.n_dof_coupling):
-            residual[i] = self.y[i]
 
+        self.residuals['outlet']['c_poly'] = self.states['outlet']['c_poly']
+        self.residuals['outlet']['viscosity'] = self.states['outlet']['viscosity']
 
 class Outlet(UnitOperationBase):
     """System outlet."""
@@ -532,6 +536,7 @@ class Cstr(UnitOperationBase):
     }
     _state_structures = ['inlet', 'bulk']
 
+
     def compute_residual(
             self,
             t: float,
@@ -565,12 +570,11 @@ def compute_residual(
         # for i in range(self.n_comp):
         #     self.residuals['bulk']['c'][i] = c_dot[i] * V + V_dot * c[i] - Q_in * c_in[i] + Q_out * c[i]
         # Alternative: Can we vectorize this?
-        self.residuals['bulk']['c'] = c_dot * V + V_dot * c - Q_in * c_in + Q_out * c
-
-        self.residuals['bulk']['Volume'] = V_dot - self.Q_in + Q_out
+        self.residuals['bulk']['c'] = calculate_residual_concentration_cstr(c, c_dot, V, V_dot,  Q_in, Q_out, c_in) 
 
-        # TODO: What about viscosities?
+        self.residuals['bulk']['Volume'] = calculate_residual_volume_cstr(V, V_dot, Q_in, Q_out)
 
+        self.residuals['inlet']['viscosity'] = calculate_residuals_visc_cstr()
 
 class DeadEndFiltration(UnitOperationBase):
     """
@@ -653,6 +657,10 @@ def compute_residual(
 
         residual[self.n_dof_coupling + 2] = ((self.c(t) * y_dot[0]) / (1-self.c(t)/self.density)) - y_dot[2]
 
+        self.residuals['retentate']
+        self.residuals['permeate']
+
+
 
 class CrossFlowFiltration(UnitOperationBase):
     """
@@ -778,4 +786,4 @@ def compute_residual(
             Residual of the unit operation.
 
         """
-        raise NotImplementedError()
+        raise NotImplementedError()

pyproject.toml

@@ -64,3 +64,8 @@ markers = [
 
 [tool.setuptools.dynamic]
 version = { attr = "CADETPythonSimulator.__version__" }
+
+[tool.ruff]
+# Same as Black.
+line-length = 88
+indent-width = 4

tests/test_residual.py

@@ -0,0 +1,103 @@
+from CADETPythonSimulator.residual import (
+    calculate_residual_volume_cstr, calculate_residual_concentration_cstr
+)
+from CADETPythonSimulator.exception import CADETPythonSimError
+import pytest
+import numpy as np
+
+
+TestCaseConc = {
+    "values" : {
+        "c" : np.array([1, 2, 3]),
+        "c_dot" : np.array([4, 5, 6]),
+        "V" : 1,
+        "V_dot" : 2,
+        "Q_in" : 3,
+        "Q_out" : 4,
+        "c_in" : np.array([7, 8, 9])
+    },
+    "expected" : np.array([-11,-7,-3])
+}
+
+
+
+@pytest.mark.parametrize(
+    "parameters",
+    [
+        TestCaseConc
+    ]
+)
+class TestResidualConcCSTR():
+    def test_calculation_concentration_cstr(self, parameters):
+
+        param_vec_conc = parameters["values"].values()
+
+        np.testing.assert_equal(calculate_residual_concentration_cstr(*param_vec_conc), parameters["expected"])
+
+
+
+
+
+TestCaseVol = {
+    "values" : {
+        "V" : 1,
+        "V_dot" : 2,
+        "Q_in" : 3,
+        "Q_out" : 4,
+    },
+    "expected" : 3
+}
+
+
+
+@pytest.mark.parametrize(
+    "parameters",
+    [
+        TestCaseVol
+    ]
+)
+
+class TestResidualVolCSTR():
+    def test_calculation_cstr(self, parameters):
+
+        param_vec_volume = parameters["values"].values()
+
+        np.testing.assert_equal(calculate_residual_volume_cstr(*param_vec_volume), parameters["expected"])
+
+TestCaseConcError = {
+    "values" : {
+        "c" : np.array([1, 2, 3]),
+        "c_dot" : np.array([4, 5, 6]),
+        "V" : -1,
+        "V_dot" : 2,
+        "Q_in" : 3,
+        "Q_out" : 4,
+        "c_in" : np.array([7, 8, 9])
+    },
+    "expected" : np.array([-11,-7,-3])
+}
+
+
+@pytest.mark.parametrize(
+    "parameters",
+    [
+        TestCaseConcError
+    ]
+)
+
+
+class TestResidualError():
+
+    def test_calculation_vol_cstr_error(self, parameters):
+
+        param_vec_volume = parameters["values"].values()
+
+        with pytest.raises(CADETPythonSimError):
+            calculate_residual_volume_cstr(*list(param_vec_volume)[2:6])
+
+    def test_calculation_concentration_cstr_error(self, parameters):
+
+        param_vec_volume = parameters["values"].values()
+
+        with pytest.raises(CADETPythonSimError):
+            calculate_residual_concentration_cstr(*param_vec_volume)

tests/test_unit_operation.py

@@ -70,6 +70,7 @@ def __init__(self, component_system=None, name='cstr', *args, **kwargs):
         super().__init__(component_system, name, *args, **kwargs)
 
 
+
 class DeadEndFiltrationFixture(UnitOperationFixture, DeadEndFiltration):
     def __init__(self,
                  component_system=None,
@@ -335,7 +336,7 @@ def test_initialize(self, unit_operation: UnitOperationBase, expected: dict):
 
 
 @pytest.mark.parametrize(
-    "unit_operation, case, expected",
+    "unit_operation, case, residualfunc, expected",
     [
         # (
         #     InletFixture(),
@@ -354,15 +355,43 @@ def test_initialize(self, unit_operation: UnitOperationBase, expected: dict):
         (
             CstrFixture(),
             {
-                'y': [0., 0., 0., 0., 0., 0., 0.],
-                'y_dot': [0., 0., 0., 0., 0., 0., 0.],
-                'Q_in': 0,
-                'Q_out': 0,
-                't': 0,
+                'states' : {
+                    'inlet' : {
+                        'c' : np.array([7, 8]),
+                        'viscosity' : [3]
+                    },
+                    'bulk' : {
+                        'c' : np.array([1, 2]),
+                        'Volume' : 1
+                    }
+                },
+                'state_derivatives' : {
+                    'inlet' : {
+                        'c' : [6, 7]
+                    },
+                    'bulk' : {
+                        'c' : np.array([4, 5]),
+                        'Volume' : 2
+                    }
+                },
+                'Q_in' : [3],
+                'Q_out' : [4]
             },
+            [
+                ("calculate_residual_concentration_cstr", lambda c, c_dot, V, V_dot,  Q_in, Q_out, c_in: c_dot * V + V_dot * c - Q_in * c_in + Q_out * c),
+                ("calculate_residuals_visc_cstr", lambda *args : 0),
+                ("calculate_residual_volume_cstr", lambda V, V_dot, Q_in, Q_out: V_dot - Q_in + Q_out)
+            ],
             {
-                'residual': [0., 0., 0., 0., 0., 0., 0.]
-            },
+                'inlet' : {
+                    'c' : np.array([7, 8]),
+                    'viscosity' : 0
+                },
+                'bulk' : {
+                    'c' :  np.array([-11,-7]),
+                    'Volume' : 3
+                }
+            }
         ),
         # (
         #     DeadEndFiltrationFixture(),
@@ -391,20 +420,33 @@ class TestUnitResidual():
 
     def test_unit_residual(
             self,
+            monkeypatch,
             unit_operation: UnitOperationBase,
             case: dict,
-            expected: dict,
+            residualfunc: dict,
+            expected: dict
             ) -> NoReturn:
         """Test the residual of unit operations."""
-        unit_operation.y = case['y']
-        unit_operation.y_dot = case['y_dot']
 
-        unit_operation.compute_residual(case['t'])
+        for funcname, func in residualfunc:
+            monkeypatch.setattr('CADETPythonSimulator.unit_operation.'+funcname, func )
+
+        for key, value in case['states'].items():
+            unit_operation.states[key] = value
+
+        for key, value in case['state_derivatives'].items():
+            unit_operation.state_derivatives[key] = value
+
+        unit_operation._Q_in = case['Q_in']
+        unit_operation._Q_out = case['Q_out']
 
-        np.testing.assert_almost_equal(unit_operation.r, expected['residual'])
+        unit_operation.compute_residual(3)
 
+        for unit_module, module_dict in expected.items():
+            for property, value in module_dict.items():
+                np.testing.assert_equal(value, unit_operation.residuals[unit_module][property])
 
 # %% Run tests
 
 if __name__ == "__main__":
-    pytest.main(["test_unit_operation.py"])
+    pytest.main(["test_unit_operation.py"])