Large commit, sry for mess reworked def unit operation, initial condi…

…tion for non circular system

CADETPythonSimulator/componentsystem.py

@@ -26,6 +26,8 @@ class CPSSpecies(Structure):
         Density of the species.
     molecular_volume : float
         The molecular volume of the species
+    specific_cake_resistance : float
+        The specific cake resistance of the species
 
     """
 
@@ -35,6 +37,8 @@ class CPSSpecies(Structure):
     density = UnsignedFloat()
     molecular_volume = UnsignedFloat()
     viscosity = UnsignedFloat()
+    pure_density = UnsignedFloat()
+    specific_cake_resistance = UnsignedFloat()
 
 
 class CPSComponent(Component):
@@ -65,6 +69,10 @@ class CPSComponent(Component):
         Molecular volume of component (including species).
     viscosity : list
         Viscosity of component(including species)
+    pure_density : list
+        pure_density of component(including species)
+    specific_cake_resistance : list
+        specific_cake_resistance of component(including species)
 
     See Also
     --------
@@ -80,7 +88,9 @@ def __init__(self,
                  molecular_weight=None,
                  density=None,
                  molecular_volume=None,
-                 viscosity=None):
+                 viscosity=None,
+                 pure_density=None,
+                 specific_cake_resistance=None):
         """Construct CPSComponent."""
         self.name = name
         self._species = []
@@ -92,7 +102,9 @@ def __init__(self,
                 molecular_weight,
                 density,
                 molecular_volume,
-                viscosity
+                viscosity,
+                pure_density,
+                specific_cake_resistance
             )
         elif isinstance(species, str):
             self.add_species(
@@ -101,7 +113,9 @@ def __init__(self,
                 molecular_weight,
                 density,
                 molecular_volume,
-                viscosity
+                viscosity,
+                pure_density,
+                specific_cake_resistance
             )
         elif isinstance(species, list):
             if charge is None:
@@ -114,14 +128,20 @@ def __init__(self,
                 molecular_volume = len(species) * [None]
             if viscosity is None:
                 viscosity = len(species) * [None]
+            if pure_density is None:
+                pure_density = len(species) * [None]
+            if specific_cake_resistance is None:
+                specific_cake_resistance = len(species) * [None]
             for i, spec in enumerate(species):
                 self.add_species(
                     spec,
                     charge[i],
                     molecular_weight[i],
                     density[i],
                     molecular_volume[i],
-                    viscosity[i]
+                    viscosity[i],
+                    pure_density[i],
+                    specific_cake_resistance[i]
                 )
         else:
             raise CADETPythonSimError("Could not determine number of species")
@@ -152,6 +172,16 @@ def viscosity(self):
         """List of float or None: The viscosity of the subspecies."""
         return [spec.viscosity for spec in self.species]
 
+    @property
+    def pure_density(self):
+        """List of float or None: Density of the subspecies."""
+        return [spec.pure_density for spec in self.species]
+
+    @property
+    def specific_cake_resistance(self):
+        """List of float or None: specific cake resistance of the subspecies."""
+        return [spec.specific_cake_resistance for spec in self.species]
+
 class CPSComponentSystem(ComponentSystem):
     """
     Component System Class.
@@ -187,6 +217,8 @@ class CPSComponentSystem(ComponentSystem):
         Molecular volume of all component species.
     viscosity : list
         Viscosity of all component species.
+    pure_density : list
+        Pure density of all component species.
 
     See Also
     --------
@@ -203,7 +235,9 @@ def __init__(
             molecular_weights=None,
             densities=None,
             molecular_volume=None,
-            viscosities=None
+            viscosities=None,
+            pure_densities=None,
+            specific_cake_resistances=None
         ):
         """
         Initialize the ComponentSystem object.
@@ -225,6 +259,10 @@ def __init__(
             The molecular volume of each component.
         viscosities : list, None
             viscosity of each component.
+        pure_densities : list, None
+            pure densities of each component.
+        specific_cake_resistances : list, None.
+            specific cake resistance of each component
 
         Raises
         ------
@@ -257,6 +295,11 @@ def __init__(
             molecular_volume = n_comp * [None]
         if viscosities is None:
             viscosities = n_comp * [None]
+        if pure_densities is None:
+            pure_densities = n_comp * [None]
+        if specific_cake_resistances is None:
+            specific_cake_resistances = n_comp * [None]
+
 
         for i, comp in enumerate(components):
             self.add_component(
@@ -265,7 +308,9 @@ def __init__(
                 molecular_weight=molecular_weights[i],
                 density=densities[i],
                 molecular_volume=molecular_volume[i],
-                viscosity=viscosities[i]
+                viscosity=viscosities[i],
+                pure_density=pure_densities[i],
+                specific_cake_resistance=specific_cake_resistances[i]
             )
 
     @wraps(CPSComponent.__init__)
@@ -320,3 +365,21 @@ def viscosities(self):
             viscosities +=comp.viscosity
 
         return viscosities
+
+    @property
+    def pure_densities(self):
+        """list: List of species density if pure."""
+        pure_densities = []
+        for comp in self.components:
+            pure_densities += comp.pure_density
+
+        return pure_densities
+
+    @property
+    def specific_cake_resistances(self):
+        """list: List of species cake resistance if pure."""
+        specific_cake_resistances = []
+        for comp in self.components:
+            specific_cake_resistances += comp.specific_cake_resistance
+
+        return specific_cake_resistances

CADETPythonSimulator/residual.py

@@ -63,7 +63,7 @@ def calculate_residual_concentration_cstr(
     Q_out : float
         Volume leaving the Unit
     c_in : np.ndarray
-        Initial concentration
+        Incomming concentration
 
     """
     if V < 0:
@@ -78,30 +78,33 @@ def calculate_residual_visc_cstr():
     return 0
 
 def calculate_residual_cake_vol_def(
-        V_dot_f: float,
+        V: float,
         rejection: np.ndarray,
-        molar_volume: np.ndarray,
+        densities: np.ndarray,
+        molecular_weights: np.ndarray,
         c_in: np.ndarray,
-        V_dot_C: float
+        V_C: float
         ) -> float:
     """
     Residual equation for the Cake Volume.
 
     Parameters
     ----------
-    V_dot_f : float
-        Flowrate of incoming feed
+    V : float
+        Volume of liquid pressured into the filter
     rejection : np.ndarray
         Rejection of the filter
-    molar_volume : np.ndarray
-        Volume of suspended material
+    densities : np.ndarray
+        densities of suspended material
+    molecular_weights : np.ndarray
+        molecular weights of the components
     c_in : np.array
         Incoming Concentration
-    V_dot_C : float
-        Change of Cake Volume
+    V_C : float
+        Cake Volume
 
     """
-    return -V_dot_C + np.sum(rejection * molar_volume * c_in * V_dot_f)
+    return -V_C + V * np.sum(rejection * c_in * molecular_weights / densities)
 
 def calculate_residual_press_easy_def(
         V_dot_Perm: float,
@@ -138,7 +141,7 @@ def calculate_residual_press_easy_def(
     """
     hyd_resistance = (Rm + alpha*V_C/A) * mu
 
-    return -V_dot_Perm + deltap * A *hyd_resistance
+    return deltap * A - V_dot_Perm* hyd_resistance
 
 def calculate_residual_visc_def():
     """Calculate the residual of the Viscosity equation of the CSTR."""

CADETPythonSimulator/solver.py

@@ -87,9 +87,7 @@ def _compute_residual(
             Array to save the calculated residual
 
         """
-        self._system.y = y
-        self._system.y_dot = y_dot
-        self._system.compute_residual(t)
+        self._system.compute_residual(t, y, y_dot)
         r[...] = self._system.r
 
     def initialize_solution_recorder(self) -> NoReturn:
@@ -157,7 +155,7 @@ def write_solution(
 
     def solve(self) -> NoReturn:
         """Simulate the system."""
-        self.initialize_system()
+        # self.initialize_system()
         self.initialize_solution_recorder()
 
         for section in self.sections:

CADETPythonSimulator/system.py

@@ -152,6 +152,22 @@ def y(self, y: np.ndarray) -> NoReturn:
             unit_operation.y = y[start_index:end_index]
             start_index = end_index
 
+    @property
+    def y_init(self) -> np.ndarray:
+        """np.ndarray: State array flattened into one dimension."""
+        return np.concatenate([
+            unit_operation.y_init for unit_operation in self.unit_operations.values()
+        ])
+
+    @y_init.setter
+    def y_init(self, y_init: np.ndarray) -> NoReturn:
+        """Setter for state in From of y_init array."""
+        start_index = 0
+        for unit_operation in self.unit_operations.values():
+            end_index = start_index + unit_operation.n_dof
+            unit_operation.y_init = y_init[start_index:end_index]
+            start_index = end_index
+
     @property
     def state_derivatives(self) -> dict[str, dict[str, State]]:
         """dict: State derivative array block of the system, indexed by name."""
@@ -237,7 +253,6 @@ def update_system_connectivity(self, connections: list) -> NoReturn:
         """
         self._compute_connectivity_matrix(connections)
         self.set_rates()
-        self.couple_unit_operations()
 
     def set_rates(self) -> NoReturn:
         """
@@ -251,6 +266,8 @@ def set_rates(self) -> NoReturn:
     def compute_residual(
             self,
             t: float,
+            y: np.ndarray,
+            y_dot: np.ndarray
             ) -> NoReturn:
         """
         Compute the residual for the differential-algebraic equations system.
@@ -259,8 +276,15 @@ def compute_residual(
         ----------
         t : float
             Current time point.
+        y : np.ndarray
+            initial state
+        y_dot : np.ndarray
+            initial state derivative
 
         """
+        self.y = y
+        self.y_dot = y_dot
+        self.r = y
         self.couple_unit_operations()
         for unit_operation in self.unit_operations.values():
             unit_operation.compute_residual(t)
@@ -382,20 +406,6 @@ def _compute_connectivity_matrix(self, connections: list) -> np.ndarray:
 
         self._connectivity = connections_matrix
 
-    def compute_residual_for_initial_values(self, y_dot: np.ndarray) -> np.ndarray:
-        """
-        Compute the residual for calculating the initial values.
-
-        Parameters
-        ----------
-        y_dot: np.ndarray
-            derivative array to calculate
-
-        """
-        self.y_dot = y_dot
-        self.compute_residual(self.t_zero)
-        return self.r
-
     def initialize_initial_values(self, t0: float):
         """
         Calculate initial values and set them up.
@@ -405,14 +415,12 @@ def initialize_initial_values(self, t0: float):
 
         """
         self.t_zero = t0
+        self.set_rates()
+        self.couple_unit_operations()
         for unit_operation in self.unit_operations.values():
             unit_operation.initialize_initial_values(t0)
+            self.couple_unit_operations()
 
-        y_dot = self.y_dot
-
-        y_dot = scipyopt.fsolve(self.compute_residual_for_initial_values, y_dot)
-
-        self.y_dot = y_dot
 
 class FlowSystem(SystemBase):
     """