diff --git a/CADETProcess/modelBuilder/carouselBuilder.py b/CADETProcess/modelBuilder/carouselBuilder.py index a4d677db..1ad176fe 100644 --- a/CADETProcess/modelBuilder/carouselBuilder.py +++ b/CADETProcess/modelBuilder/carouselBuilder.py @@ -283,12 +283,12 @@ def _add_inter_zone_connections(self, flow_sheet: FlowSheet) -> NoReturn: origin = unit.outlet_unit else: origin = unit + if connections.destinations: + for destination in connections.destinations[None]: + if isinstance(destination, ZoneBaseClass): + destination = destination.inlet_unit - for destination in connections.destinations: - if isinstance(destination, ZoneBaseClass): - destination = destination.inlet_unit - - flow_sheet.add_connection(origin, destination) + flow_sheet.add_connection(origin, destination) for zone in self.zones: output_state = self.flow_sheet.output_states[zone] diff --git a/CADETProcess/modelBuilder/compartmentBuilder.py b/CADETProcess/modelBuilder/compartmentBuilder.py index 9a99dde8..3a488e28 100644 --- a/CADETProcess/modelBuilder/compartmentBuilder.py +++ b/CADETProcess/modelBuilder/compartmentBuilder.py @@ -290,14 +290,19 @@ def validate_flow_rates(self): flow_rates = self.flow_sheet.get_flow_rates() for comp in self._real_compartments: - if not np.all( - np.isclose( - flow_rates[comp.name].total_in, - flow_rates[comp.name].total_out - )): - raise CADETProcessError( - f"Unbalanced flow rate for compartment '{comp.name}'." - ) + for port in flow_rates[comp.name].total_in: + if not np.all( + np.isclose( + flow_rates[comp.name].total_in[port], + flow_rates[comp.name].total_out[port] + )): + if comp.n_ports == 1: + msg = comp.name + else: + msg = f"{comp.name} at Port {port}" + raise CADETProcessError( + f"Unbalanced flow rate for compartment '{msg}'." + ) class CompartmentModel(Cstr): diff --git a/CADETProcess/processModel/discretization.py b/CADETProcess/processModel/discretization.py index a050d471..60b92c56 100644 --- a/CADETProcess/processModel/discretization.py +++ b/CADETProcess/processModel/discretization.py @@ -26,7 +26,8 @@ class for all other classes in this module and defines some common parameters. 'LRMPDiscretizationFV', 'LRMPDiscretizationDG', 'GRMDiscretizationFV', 'GRMDiscretizationDG', 'WenoParameters', 'ConsistencySolverParameters', - 'DGMixin' + 'DGMixin', + 'MCTDiscretizationFV', ] @@ -571,3 +572,30 @@ class ConsistencySolverParameters(Structure): 'solver_name', 'init_damping', 'min_damping', 'max_iterations', 'subsolvers' ] + + +class MCTDiscretizationFV(DiscretizationParametersBase): + """Discretization parameters of the FV version of the MCT. + + Attributes + ---------- + ncol : UnsignedInteger, optional + Number of axial column discretization cells. Default is 100. + use_analytic_jacobian : Bool, optional + If True, use analytically computed Jacobian matrix (faster). + If False, use Jacobians generated by algorithmic differentiation (slower). + Default is True. + reconstruction : Switch, optional + Method for spatial reconstruction. Valid values are 'WENO' (Weighted + Essentially Non-Oscillatory). Default is 'WENO'. + + """ + + ncol = UnsignedInteger(default=100) + use_analytic_jacobian = Bool(default=True) + reconstruction = Switch(default='WENO', valid=['WENO']) + + _parameters = [ + 'ncol', 'use_analytic_jacobian', 'reconstruction', + ] + _dimensionality = ['ncol'] diff --git a/CADETProcess/processModel/flowSheet.py b/CADETProcess/processModel/flowSheet.py index fa5985a3..2badae1d 100644 --- a/CADETProcess/processModel/flowSheet.py +++ b/CADETProcess/processModel/flowSheet.py @@ -1,5 +1,6 @@ from functools import wraps from warnings import warn +from collections import defaultdict import numpy as np from addict import Dict @@ -176,6 +177,32 @@ def get_unit_index(self, unit): return self.units.index(unit) + def get_port_index(self, unit, port): + """Return the port index of a unit + + Parameters + ---------- + unit : UnitBaseClass + UnitBaseClass object of wich port index is to be returned + port : string + Name of port which index is to be returned + Raises + ------ + CADETProcessError + If unit or port is not in the current flow sheet. + Returns + ------- + port_index : int + Returns the port index of the port of the unit_operation. + """ + if unit not in self.units: + raise CADETProcessError('Unit not in flow sheet') + + + port_index = unit.ports.index(port) + + return port_index + @property def inlets(self): """list: All Inlets in the system.""" @@ -237,11 +264,24 @@ def add_unit( raise CADETProcessError('Component systems do not match.') self._units.append(unit) - self._connections[unit] = Dict({ - 'origins': [], - 'destinations': [], - }) - self._output_states[unit] = [] + + + for port in unit.ports: + + if isinstance(unit, Inlet): + self._connections[unit]['origins'] = None + self._connections[unit]['destinations'][port] = defaultdict(list) + + elif isinstance(unit, Outlet): + self._connections[unit]['origins'][port] = defaultdict(list) + self._connections[unit]['destinations'] = None + + else: + self._connections[unit]['origins'][port] = defaultdict(list) + self._connections[unit]['destinations'][port] = defaultdict(list) + + + self._output_states[unit] = Dict() self._flow_rates[unit] = [] super().__setattr__(unit.name, unit) @@ -292,13 +332,24 @@ def remove_unit(self, unit): if unit is self.product_outlets: self.remove_product_outlet(unit) - origins = self.connections[unit].origins.copy() + origins = [] + destinations = [] + + if self._connections[unit]['origins'] is not None: + origins = [origin for ports in self._connections[unit]['origins'] for origin in self._connections[unit]['origins'][ports]] + + if self._connections[unit]['destinations'] is not None: + destinations = [destination for ports in self._connections[unit]['destinations'] for destination in self._connections[unit]['destinations'][ports]].copy() + for origin in origins: - self.remove_connection(origin, unit) + for origin_port in self._connections[unit]['origins']: + for unit_port in self._connections[unit]['origins'][origin_port][origin]: + self.remove_connection(origin, unit, origin_port, unit_port) - destinations = self.connections[unit].destinations.copy() for destination in destinations: - self.remove_connection(unit, destination) + for destination_port in self._connections[unit]['destinations']: + for unit_port in self._connections[unit]['destinations'][destination_port][destination]: + self.remove_connection(unit, destination, unit_port, destination_port) self._units.remove(unit) self._connections.pop(unit) @@ -319,7 +370,8 @@ def connections(self): @origin_destination_name_decorator @update_parameters_decorator - def add_connection(self, origin, destination): + def add_connection( + self, origin, destination, origin_port=None, destination_port=None): """Add connection between units 'origin' and 'destination'. Parameters @@ -328,6 +380,10 @@ def add_connection(self, origin, destination): UnitBaseClass from which the connection originates. destination : UnitBaseClass UnitBaseClass where the connection terminates. + origin_port : str + Port from which connection originates. + destination_port : str + Port where connection terminates. Raises ------ @@ -352,17 +408,43 @@ def add_connection(self, origin, destination): if isinstance(destination, Inlet): raise CADETProcessError("Inlet unit cannot have ingoing stream.") - if destination in self.connections[origin].destinations: - raise CADETProcessError('Connection already exists') - self._connections[origin].destinations.append(destination) - self._connections[destination].origins.append(origin) + if origin.has_ports and origin_port is None: + raise CADETProcessError('Missing `origin_port`') + if not origin.has_ports and origin_port is not None: + raise CADETProcessError(f'Origin unit does not support ports.') + if origin.has_ports and origin_port not in origin.ports: + raise CADETProcessError(f'Origin port "{origin_port}" not found in ports: {origin.ports}.') + if origin_port in self._connections[destination]['origins'][destination_port][origin]: + raise CADETProcessError("Connection already exists") + + + + if destination.has_ports and destination_port is None: + raise CADETProcessError('Missing `destination_port`') + if not destination.has_ports and destination_port is not None: + raise CADETProcessError('Destination unit does not support ports.') + if destination.has_ports and destination_port not in destination.ports: + raise CADETProcessError(f'destination port "{destination_port}" not found in ports: {destination.ports}.') + + if destination_port in self._connections[origin]['destinations'][origin_port][destination]: + raise CADETProcessError("Connection already exists") + + + if not destination.has_ports: + destination_port = destination.ports[0] + if not origin.has_ports: + origin_port = origin.ports[0] + + self._connections[destination]['origins'][destination_port][origin].append(origin_port) + self._connections[origin]['destinations'][origin_port][destination].append(destination_port) + + self.set_output_state(origin, 0, origin_port) - self.set_output_state(origin, 0) @origin_destination_name_decorator @update_parameters_decorator - def remove_connection(self, origin, destination): + def remove_connection( self, origin, destination, origin_port=None, destination_port=None): """Remove connection between units 'origin' and 'destination'. Parameters @@ -371,6 +453,10 @@ def remove_connection(self, origin, destination): UnitBaseClass from which the connection originates. destination : UnitBaseClass UnitBaseClass where the connection terminates. + origin_port : int + Port from which connection originates. + destination_port : int + Port where connection terminates. Raises ------ @@ -386,17 +472,31 @@ def remove_connection(self, origin, destination): """ if origin not in self._units: raise CADETProcessError('Origin not in flow sheet') + if origin.has_ports and origin_port is None: + raise CADETProcessError('Missing `origin_port`') + + if origin_port not in origin.ports: + raise CADETProcessError(f'Origin port {origin_port} not in Unit {origin.name}.') + if destination not in self._units: raise CADETProcessError('Destination not in flow sheet') + if destination.has_ports and origin_port is None: + raise CADETProcessError('Missing `destination_port`') + + if destination_port not in destination.ports: + raise CADETProcessError(f'Origin port {destination_port} not in Unit {destination.name}.') try: - self._connections[origin].destinations.remove(destination) - self._connections[destination].origins.remove(origin) + + self._connections[destination]['origins'][destination_port].pop(origin) + self._connections[origin]['destinations'][origin_port].pop(destination) except KeyError: raise CADETProcessError('Connection does not exist.') + + @origin_destination_name_decorator - def connection_exists(self, origin, destination): + def connection_exists(self, origin, destination, origin_port=None, destination_port=None): """bool: check if connection exists in flow sheet. Parameters @@ -407,8 +507,22 @@ def connection_exists(self, origin, destination): UnitBaseClass where the connection terminates. """ - if destination in self._connections[origin].destinations \ - and origin in self._connections[destination].origins: + if origin.has_ports and not origin_port: + raise CADETProcessError(f'Origin port needs to be specified for unit operation with ports {origin.name}.') + + if destination.has_ports and not destination_port: + raise CADETProcessError(f'Destination port needs to be specified for unit operation with ports {destination.name}.') + + if origin_port not in origin.ports: + raise CADETProcessError(f'{origin.name} does not have port {origin_port}') + + if destination_port not in destination.ports: + raise CADETProcessError(f'{destination.name} does not have port {destination_port}') + + if destination in self._connections[origin].destinations[origin_port]\ + and destination_port in self._connections[origin].destinations[origin_port][destination]\ + and origin in self._connections[destination].origins[destination_port]\ + and origin_port in self._connections[destination].origins[destination_port][origin]: return True return False @@ -432,14 +546,14 @@ def check_connections(self): flag = True for unit, connections in self.connections.items(): if isinstance(unit, Inlet): - if len(connections.origins) != 0: + if connections.origins != None: flag = False warn("Inlet unit cannot have ingoing stream.") if len(connections.destinations) == 0: flag = False warn(f" Unit '{unit.name}' does not have outgoing stream.") elif isinstance(unit, Outlet): - if len(connections.destinations) != 0: + if connections.destinations != None: flag = False warn("Outlet unit cannot have outgoing stream.") if len(connections.origins) == 0: @@ -450,10 +564,10 @@ def check_connections(self): flag = False warn("Cstr cannot have flow rate without outgoing stream.") else: - if len(connections.origins) == 0: + if all(len(port_list) == 0 for port_list in connections.origins.values()): flag = False warn(f"Unit '{unit.name}' does not have ingoing stream.") - if len(connections.destinations) == 0: + if all(len(port_list) == 0 for port_list in connections.destinations.values()): flag = False warn(f" Unit '{unit.name}' does not have outgoing stream.") @@ -487,11 +601,20 @@ def check_units_config(self): @property def output_states(self): - return self._output_states + + output_states_dict = self._output_states.copy() + + for unit, ports in output_states_dict.items(): + for port in ports: + if port == None: + output_states_dict[unit] = output_states_dict[unit][port] + + + return output_states_dict @unit_name_decorator @update_parameters_decorator - def set_output_state(self, unit, state): + def set_output_state(self, unit, state, port=None): """Set split ratio of outgoing streams for UnitOperation. Parameters @@ -500,6 +623,8 @@ def set_output_state(self, unit, state): UnitOperation of flowsheet. state : int or list of floats or dict new output state of the unit. + port : int + Port for which to set the output state. Raises ------ @@ -508,12 +633,40 @@ def set_output_state(self, unit, state): If state is integer and the state >= the state_length. If the length of the states is unequal the state_length. If the sum of the states is not equal to 1. - + If port exceeds the number of ports """ + def get_port_index(unit_connection_dict, destination, destination_port): + """helper function to classify the index of a connection for your outputstate + + Parameters + ---------- + + unit_connection_dict : Defaultdict + contains dict with connected units and their respective ports + destination : UnitBaseClass + destination object + destination_port : int + destination Port + + """ + + ret_index = 0 + for unit_destination in unit_connection_dict: + if unit_destination is destination: + ret_index+=unit_connection_dict[unit_destination].index(destination_port) + break + ret_index+=len(unit_connection_dict[unit_destination]) + return ret_index + + if unit not in self._units: raise CADETProcessError('Unit not in flow sheet') - state_length = len(self.connections[unit].destinations) + if port not in unit.ports: + raise CADETProcessError(f'Port {port} is not a port of Unit {unit.name}') + + + state_length = sum([len(self._connections[unit].destinations[port][unit_name]) for unit_name in self._connections[unit].destinations[port]]) if state_length == 0: output_state = [] @@ -529,18 +682,34 @@ def set_output_state(self, unit, state): output_state = [0.0] * state_length for dest, value in state.items(): try: - assert self.connection_exists(unit, dest) - except AssertionError: - raise CADETProcessError(f'{unit} does not connect to {dest}.') - dest = self[dest] - index = self.connections[unit].destinations.index(dest) - output_state[index] = value - - elif isinstance(state, list): + for destination_port, output_value in value.items(): + try: + assert self.connection_exists(unit, dest, port, destination_port) + except AssertionError: + raise CADETProcessError(f'{unit} on port {port} does not connect to {dest} on port {destination_port}.') + inner_dest = self[dest] + index = get_port_index(self._connections[unit].destinations[port], inner_dest, destination_port) + output_state[index] = output_value + except AttributeError: + destination_port = None + try: + assert self.connection_exists(unit, dest, port, destination_port) + except AssertionError: + raise CADETProcessError(f'{unit} on port {port} does not connect to {dest} on port {destination_port}.') + dest = self[dest] + index = get_port_index(self.connections[unit].destinations[port], dest, destination_port) + output_state[index] = value + + elif isinstance(state, (list)): if len(state) != state_length: raise CADETProcessError(f'Expected length {state_length}.') output_state = state + elif isinstance(state, np.ndarray): + if len(state) != state_length: + raise CADETProcessError(f'Expected length {state_length}.') + + output_state = list(state) else: raise TypeError("Output state must be integer, list or dict.") @@ -548,18 +717,21 @@ def set_output_state(self, unit, state): if state_length != 0 and not np.isclose(sum(output_state), 1): raise CADETProcessError('Sum of fractions must be 1') - self._output_states[unit] = output_state + self._output_states[unit][port] = output_state - def get_flow_rates(self, state=None): - """ - Calculate the volumetric flow rate for all connections in the process. + def get_flow_rates(self, state=None, eps=5.9e16): + """Calculate flow rate for all connections. + + Optionally, an additional output state can be passed to update the + current output states. Parameters ---------- state : Dict, optional Updated flow rates and output states for process sections. Default is None. - + eps : float, optional + eps as an upper boarder for condition of flow_rate calculation Returns ------- Dict @@ -590,43 +762,74 @@ def get_flow_rates(self, state=None): Forum discussion on flow rate calculation: https://forum.cadet-web.de/t/improving-the-flowrate-calculation/795 """ + port_index_list = [] + port_number = 0 + + for unit in self.units: + for port in unit.ports: + port_index_list.append((unit, port)) + port_number += 1 + flow_rates = { unit.name: unit.flow_rate for unit in (self.inlets + self.cstrs) if unit.flow_rate is not None } - output_states = self.output_states + output_states = self._output_states if state is not None: for param, value in state.items(): param = param.split('.') - unit_name = param[1] param_name = param[-1] if param_name == 'flow_rate': + unit_name = param[1] flow_rates[unit_name] = value[0] - elif unit_name == 'output_states': - unit = self.units_dict[param_name] - output_states[unit] = list(value.ravel()) - - n_units = self.number_of_units + elif param[1] == 'output_states': + unit_name = param[2] + unit = self.units_dict[unit_name] + if unit.has_ports: + port = param[2] + else: + port = None + output_states[unit][port] = list(value.ravel()) # Setup matrix with output states. - w_out = np.zeros((n_units, n_units)) + w_out = np.zeros((port_number, port_number)) for unit in self.units: - unit_index = self.get_unit_index(unit) + if unit.name in flow_rates: + unit_index = port_index_list.index((unit, None)) w_out[unit_index, unit_index] = 1 else: - for origin in self.connections[unit]['origins']: - o_index = self.get_unit_index(origin) - local_d_index = self.connections[origin].destinations.index(unit) - w_out[unit_index, o_index] = output_states[origin][local_d_index] - w_out[unit_index, unit_index] += -1 + + for port in self._connections[unit]['origins']: + + port_index = port_index_list.index((unit, port)) + + for origin_unit in self._connections[unit]['origins'][port]: + + for origin_port in self._connections[unit]['origins'][port][origin_unit]: + + o_index = port_index_list.index((origin_unit, origin_port)) + + local_d_index = 0 + + for inner_unit in self._connections[origin_unit]['destinations'][origin_port]: + if inner_unit == unit: + break + local_d_index += len(list(self._connections[origin_unit]['destinations'][origin_port][inner_unit])) + + local_d_index += self._connections[origin_unit]['destinations'][origin_port][unit].index(port) + + if output_states[origin_unit][origin_port][local_d_index]: + w_out[port_index, o_index] = output_states[origin_unit][origin_port][local_d_index] + + w_out[port_index, port_index] += -1 # Check for a singular matrix before the loop - if np.linalg.cond(w_out) == np.inf: + if np.linalg.cond(w_out) > eps: raise CADETProcessError( "Flow sheet connectivity matrix is singular, which may be due to " "unconnected units or missing flow rates. Please ensure all units are " @@ -634,7 +837,7 @@ def get_flow_rates(self, state=None): ) # Solve system of equations for each polynomial coefficient - total_flow_rate_coefficents = np.zeros((4, n_units)) + total_flow_rate_coefficents = np.zeros((4, port_number)) for i in range(4): if len(flow_rates) == 0: continue @@ -643,10 +846,10 @@ def get_flow_rates(self, state=None): if not np.any(coeffs): continue - Q_vec = np.zeros(n_units) + Q_vec = np.zeros(port_number) for unit_name in flow_rates: - unit_index = self.get_unit_index(self.units_dict[unit_name]) - Q_vec[unit_index] = flow_rates[unit_name][i] + port_index = port_index_list.index((self.units_dict[unit_name], None)) + Q_vec[port_index] = flow_rates[unit_name][i] try: total_flow_rate_coefficents[i, :] = np.linalg.solve(w_out, Q_vec) except np.linalg.LinAlgError: @@ -656,50 +859,103 @@ def get_flow_rates(self, state=None): ) # w_out_help is the same as w_out but it contains the origin flow for every unit - w_out_help = np.zeros((n_units, n_units)) + w_out_help = np.zeros((port_number, port_number)) + + + for unit in self.units: + if self._connections[unit]['origins']: + for port in self._connections[unit]['origins']: + + port_index = port_index_list.index((unit, port)) - for unit in self.connections: - unit_index = self.get_unit_index(unit) - for origin in self.connections[unit]['origins']: - o_index = self.get_unit_index(origin) - local_d_index = self.connections[origin].destinations.index(unit) - w_out_help[unit_index, o_index] = output_states[origin][local_d_index] + for origin_unit in self._connections[unit]['origins'][port]: + + for origin_port in self._connections[unit]['origins'][port][origin_unit]: + o_index = port_index_list.index((origin_unit, origin_port)) + + local_d_index = 0 + + for inner_unit in self._connections[origin_unit]['destinations'][origin_port]: + if inner_unit == unit: + break + local_d_index += len(list(self._connections[origin_unit]['destinations'][origin_port][inner_unit])) + + local_d_index += self._connections[origin_unit]['destinations'][origin_port][unit].index(port) + + if not output_states[origin_unit][origin_port][local_d_index]: + w_out_help[port_index, o_index] = 0 + else : + w_out_help[port_index, o_index] = output_states[origin_unit][origin_port][local_d_index] # Calculate total_in as a matrix in "one" step rather than iterating manually. total_in_matrix = w_out_help @ total_flow_rate_coefficents.T # Generate output dict return_flow_rates = Dict() - for index, unit in enumerate(self.units): + for unit in self.units: unit_solution_dict = Dict() if not isinstance(unit, Inlet): - unit_solution_dict['total_in'] = list(total_in_matrix[index]) + + unit_solution_dict['total_in'] = Dict() + + for unit_port in unit.ports: + + index = port_index_list.index((unit, unit_port)) + unit_solution_dict['total_in'][unit_port] = list(total_in_matrix[index]) if not isinstance(unit, Outlet): - unit_solution_dict['total_out'] = list(total_flow_rate_coefficents[:, index]) + + unit_solution_dict['total_out'] = Dict() + + for unit_port in unit.ports: + + index = port_index_list.index((unit, unit_port)) + unit_solution_dict['total_out'][unit_port] = list(total_flow_rate_coefficents[:, index]) + if not isinstance(unit, Inlet): - unit_solution_dict['origins'] = Dict( - { - origin.name: list( - total_flow_rate_coefficents[:, self.get_unit_index(origin)] - * w_out_help[index, self.get_unit_index(origin)] - ) - for origin in self.connections[unit].origins - } - ) + unit_solution_dict['origins'] = Dict() + + for unit_port in self._connections[unit]['origins']: + + if self._connections[unit]['origins'][unit_port]: + + unit_solution_dict['origins'][unit_port] = Dict() + + for origin_unit in self._connections[unit]['origins'][unit_port]: + + unit_solution_dict['origins'][unit_port][origin_unit.name] = Dict() + + for origin_port in self._connections[unit]['origins'][unit_port][origin_unit]: + origin_port_index = port_index_list.index((origin_unit, origin_port)) + unit_port_index = port_index_list.index((unit, unit_port)) + flow_list = list( + total_flow_rate_coefficents[:, origin_port_index] + * w_out_help[unit_port_index, origin_port_index]) + unit_solution_dict['origins'][unit_port][origin_unit.name][origin_port] = flow_list if not isinstance(unit, Outlet): - unit_solution_dict['destinations'] = Dict( - { - destination.name: list( - total_flow_rate_coefficents[:, index] - * w_out_help[self.get_unit_index(destination), index] - ) - for destination in self.connections[unit].destinations - } - ) + + unit_solution_dict['destinations'] = Dict() + + for unit_port in self._connections[unit]['destinations']: + + if self._connections[unit]['destinations'][unit_port]: + + unit_solution_dict['destinations'][unit_port] = Dict() + + for destination_unit in self._connections[unit]['destinations'][unit_port]: + + unit_solution_dict['destinations'][unit_port][destination_unit.name] = Dict() + + for destination_port in self._connections[unit]['destinations'][unit_port][destination_unit]: + destination_port_index = port_index_list.index((destination_unit, destination_port)) + unit_port_index = port_index_list.index((unit, unit_port)) + flow_list = list( + total_flow_rate_coefficents[:, unit_port_index] + * w_out_help[destination_port_index, unit_port_index]) + unit_solution_dict['destinations'][unit_port][destination_unit.name][destination_port] = flow_list return_flow_rates[unit.name] = unit_solution_dict @@ -871,7 +1127,12 @@ def parameters(self, parameters): output_states = parameters.pop('output_states') for unit, state in output_states.items(): unit = self.units_dict[unit] - self.set_output_state(unit, state) + # Hier, if unit.has_ports: iterate. + if not unit.has_ports: + self.set_output_state(unit, state) + else: + for port_i, state_i in state.items(): + self.set_output_state(unit, state_i, port_i) except KeyError: pass diff --git a/CADETProcess/processModel/process.py b/CADETProcess/processModel/process.py index 04db23dc..78fec48d 100644 --- a/CADETProcess/processModel/process.py +++ b/CADETProcess/processModel/process.py @@ -90,7 +90,7 @@ def m_feed(self): feed_all = np.zeros((self.n_comp,)) for feed in self.flow_sheet.feed_inlets: - feed_flow_rate_time_line = flow_rate_timelines[feed.name].total_out + feed_flow_rate_time_line = flow_rate_timelines[feed.name].total_out[None] feed_signal_param = f'flow_sheet.{feed.name}.c' if feed_signal_param in self.parameter_timelines: tl = self.parameter_timelines[feed_signal_param] @@ -122,7 +122,7 @@ def V_eluent(self): V_all = 0 for eluent in self.flow_sheet.eluent_inlets: - eluent_time_line = flow_rate_timelines[eluent.name]['total_out'] + eluent_time_line = flow_rate_timelines[eluent.name]['total_out'][None] V_eluent = eluent_time_line.integral().squeeze() V_all += V_eluent @@ -140,10 +140,10 @@ def flow_rate_timelines(self): """dict: TimeLine of flow_rate for all unit_operations.""" flow_rate_timelines = { unit.name: { - 'total_in': TimeLine(), - 'origins': defaultdict(TimeLine), - 'total_out': TimeLine(), - 'destinations': defaultdict(TimeLine) + 'total_in': defaultdict(TimeLine), + 'origins': defaultdict( lambda: defaultdict( lambda: defaultdict(TimeLine))), + 'total_out': defaultdict(TimeLine), + 'destinations': defaultdict( lambda: defaultdict( lambda: defaultdict(TimeLine))), } for unit in self.flow_sheet.units } @@ -160,40 +160,55 @@ def flow_rate_timelines(self): flow_rates = self.flow_sheet.get_flow_rates(state) - for unit, flow_rate in flow_rates.items(): + for unit, flow_rate_dict in flow_rates.items(): unit_flow_rates = flow_rate_timelines[unit] # If inlet, also use outlet for total_in + if isinstance(self.flow_sheet[unit], Inlet): - section = Section( - start, end, flow_rate.total_out, is_polynomial=True - ) + for port in flow_rate_dict['total_out']: + section = Section( + start, end, flow_rate_dict.total_out[port], is_polynomial=True + ) + unit_flow_rates['total_in'][port].add_section(section) else: - section = Section( - start, end, flow_rate.total_in, is_polynomial=True - ) - unit_flow_rates['total_in'].add_section(section) - for orig, flow_rate_orig in flow_rate.origins.items(): - section = Section( - start, end, flow_rate_orig, is_polynomial=True - ) - unit_flow_rates['origins'][orig].add_section(section) + for port in flow_rate_dict['total_in']: + section = Section( + start, end, flow_rate_dict.total_in[port], is_polynomial=True + ) + unit_flow_rates['total_in'][port].add_section(section) + + for port in flow_rate_dict.origins: + + for orig, origin_port_dict in flow_rate_dict.origins[port].items(): + for orig_port, flow_rate_orig in origin_port_dict.items(): + section = Section( + start, end, flow_rate_orig, is_polynomial=True + ) + unit_flow_rates['origins'][port][orig][orig_port].add_section(section) # If outlet, also use inlet for total_out + if isinstance(self.flow_sheet[unit], Outlet): - section = Section( - start, end, flow_rate.total_in, is_polynomial=True - ) + for port in flow_rate_dict['total_in']: + section = Section( + start, end, flow_rate_dict.total_in[port], is_polynomial=True + ) + unit_flow_rates['total_out'][port].add_section(section) else: - section = Section( - start, end, flow_rate.total_out, is_polynomial=True - ) - unit_flow_rates['total_out'].add_section(section) - for dest, flow_rate_dest in flow_rate.destinations.items(): - section = Section( - start, end, flow_rate_dest, is_polynomial=True - ) - unit_flow_rates['destinations'][dest].add_section(section) + for port in flow_rate_dict['total_out']: + section = Section( + start, end, flow_rate_dict.total_out[port], is_polynomial=True + ) + unit_flow_rates['total_out'][port].add_section(section) + + for port in flow_rate_dict.destinations: + for dest, dest_port_dict in flow_rate_dict.destinations[port].items(): + for dest_port, flow_rate_dest in dest_port_dict.items(): + section = Section( + start, end, flow_rate_dest, is_polynomial=True + ) + unit_flow_rates['destinations'][port][dest][dest_port].add_section(section) return Dict(flow_rate_timelines) @@ -203,10 +218,10 @@ def flow_rate_section_states(self): section_states = { time: { unit.name: { - 'total_in': [], - 'origins': defaultdict(dict), - 'total_out': [], - 'destinations': defaultdict(dict), + 'total_in': defaultdict(list), + 'origins': defaultdict( lambda: defaultdict( lambda: defaultdict(list))), + 'total_out': defaultdict(list), + 'destinations': defaultdict( lambda: defaultdict( lambda: defaultdict(list))), } for unit in self.flow_sheet.units } for time in self.section_times[0:-1] } @@ -214,26 +229,35 @@ def flow_rate_section_states(self): for sec_time in self.section_times[0:-1]: for unit, unit_flow_rates in self.flow_rate_timelines.items(): if isinstance(self.flow_sheet[unit], Inlet): - section_states[sec_time][unit]['total_in'] \ - = unit_flow_rates['total_out'].coefficients(sec_time) + for port in unit_flow_rates['total_out']: + section_states[sec_time][unit]['total_in'][port] \ + = unit_flow_rates['total_out'][port].coefficients(sec_time) else: - section_states[sec_time][unit]['total_in'] \ - = unit_flow_rates['total_in'].coefficients(sec_time) + for port in unit_flow_rates['total_in']: + section_states[sec_time][unit]['total_in'][port] \ + = unit_flow_rates['total_in'][port].coefficients(sec_time) + + for port, orig_dict in unit_flow_rates.origins.items(): + for origin in orig_dict: + for origin_port, tl in orig_dict[origin].items(): + section_states[sec_time][unit]['origins'][port][origin][origin_port]\ + = tl.coefficients(sec_time) - for orig, tl in unit_flow_rates.origins.items(): - section_states[sec_time][unit]['origins'][orig] \ - = tl.coefficients(sec_time) if isinstance(self.flow_sheet[unit], Outlet): - section_states[sec_time][unit]['total_out'] \ - = unit_flow_rates['total_in'].coefficients(sec_time) + for port in unit_flow_rates['total_in']: + section_states[sec_time][unit]['total_out'][port] \ + = unit_flow_rates['total_in'][port].coefficients(sec_time) else: - section_states[sec_time][unit]['total_out'] \ - = unit_flow_rates['total_out'].coefficients(sec_time) + for port in unit_flow_rates['total_out']: + section_states[sec_time][unit]['total_out'][port] \ + = unit_flow_rates['total_out'][port].coefficients(sec_time) - for dest, tl in unit_flow_rates.destinations.items(): - section_states[sec_time][unit]['destinations'][dest] \ - = tl.coefficients(sec_time) + for port, dest_dict in unit_flow_rates.destinations.items(): + for dest in dest_dict: + for dest_port, tl in dest_dict[dest].items(): + section_states[sec_time][unit]['destinations'][port][dest][dest_port] \ + = tl.coefficients(sec_time) return Dict(section_states) @@ -685,11 +709,13 @@ def check_cstr_volume(self): if cstr.flow_rate is None: continue V_0 = cstr.V - V_in = self.flow_rate_timelines[cstr.name].total_in.integral() - V_out = self.flow_rate_timelines[cstr.name].total_out.integral() - if V_0 + V_in - V_out < 0: - flag = False - warn(f'CSTR {cstr.name} runs empty during process.') + unit_index = self.flow_sheet.get_unit_index(cstr) + for port in self.flow_sheet.units[unit_index].ports: + V_in = self.flow_rate_timelines[cstr.name].total_in[port].integral() + V_out = self.flow_rate_timelines[cstr.name].total_out[port].integral() + if V_0 + V_in - V_out < 0: + flag = False + warn(f'CSTR {cstr.name} runs empty on port {port} during process.') return flag diff --git a/CADETProcess/processModel/solutionRecorder.py b/CADETProcess/processModel/solutionRecorder.py index cd57a8d4..7e98b3de 100644 --- a/CADETProcess/processModel/solutionRecorder.py +++ b/CADETProcess/processModel/solutionRecorder.py @@ -363,3 +363,22 @@ class CSTRRecorder( """ pass + + +class MCTRecorder( + SolutionRecorderBase, + BaseMixin, + IOMixin, + BulkMixin): + """Recorder for TubularReactor. + + See Also + -------- + BaseMixin + IOMixin + BulkMixin + CADETProcess.processModel.MCT + + """ + + pass diff --git a/CADETProcess/processModel/unitOperation.py b/CADETProcess/processModel/unitOperation.py index e1563029..fd76db4b 100644 --- a/CADETProcess/processModel/unitOperation.py +++ b/CADETProcess/processModel/unitOperation.py @@ -7,10 +7,10 @@ from CADETProcess.dataStructure import frozen_attributes from CADETProcess.dataStructure import Structure from CADETProcess.dataStructure import ( - Constant, UnsignedFloat, + Constant, UnsignedFloat, UnsignedInteger, String, Switch, SizedUnsignedList, - Polynomial, NdPolynomial, SizedList + Polynomial, NdPolynomial, SizedList, SizedNdArray ) from .componentSystem import ComponentSystem @@ -20,12 +20,14 @@ DiscretizationParametersBase, NoDiscretization, LRMDiscretizationFV, LRMDiscretizationDG, LRMPDiscretizationFV, LRMPDiscretizationDG, - GRMDiscretizationFV, GRMDiscretizationDG + GRMDiscretizationFV, GRMDiscretizationDG, + MCTDiscretizationFV, ) from .solutionRecorder import ( IORecorder, - TubularReactorRecorder, LRMRecorder, LRMPRecorder, GRMRecorder, CSTRRecorder + TubularReactorRecorder, LRMRecorder, LRMPRecorder, GRMRecorder, CSTRRecorder, + MCTRecorder, ) @@ -41,6 +43,7 @@ 'LumpedRateModelWithoutPores', 'LumpedRateModelWithPores', 'GeneralRateModel', + 'MCT' ] @@ -61,6 +64,8 @@ class UnitBaseClass(Structure): list of parameter names. name : String name of the unit operation. + has_ports : bool + flag if unit has ports. The default is False. binding_model : BindingBaseClass binding behavior of the unit. Defaults to NoBinding. solution_recorder : IORecorder @@ -79,23 +84,48 @@ class UnitBaseClass(Structure): _section_dependent_parameters = [] _initial_state = [] + has_ports = False supports_binding = False supports_bulk_reaction = False supports_particle_reaction = False discretization_schemes = () - def __init__(self, component_system, name, *args, **kwargs): + def __init__( + self, + component_system, + name, + *args, + binding_model=None, + bulk_reaction_model=None, + particle_reaction_model=None, + discretization=None, + solution_recorder=None, + **kwargs + ): self.name = name self.component_system = component_system - self.binding_model = NoBinding() + if binding_model is None: + binding_model = NoBinding() + self.binding_model = binding_model - self.bulk_reaction_model = NoReaction() - self.particle_reaction_model = NoReaction() + if bulk_reaction_model is None: + bulk_reaction_model = NoReaction() + self.bulk_reaction_model = bulk_reaction_model + + if particle_reaction_model is None: + particle_reaction_model = NoReaction() + self.particle_reaction_model = particle_reaction_model + + if discretization is None: + discretization = NoDiscretization() + self.discretization = discretization + + if solution_recorder is None: + solution_recorder = IORecorder() + self.solution_recorder = solution_recorder - self.discretization = NoDiscretization() - self.solution_recorder = IORecorder() super().__init__(*args, **kwargs) @@ -134,6 +164,14 @@ def discretization(self, discretization): def n_comp(self): return self.component_system.n_comp + @property + def ports(self): + return [None] + + @property + def n_ports(self): + return 1 + @property def parameters(self): """dict: Dictionary with parameter values.""" @@ -714,14 +752,17 @@ class TubularReactor(TubularReactorBase): _parameters = ['c'] def __init__(self, *args, discretization_scheme='FV', **kwargs): - super().__init__(*args, **kwargs) - if discretization_scheme == 'FV': - self.discretization = LRMDiscretizationFV() + discretization = LRMDiscretizationFV() elif discretization_scheme == 'DG': - self.discretization = LRMDiscretizationDG() + discretization = LRMDiscretizationDG() - self.solution_recorder = TubularReactorRecorder() + super().__init__( + *args, + discretization=discretization, + solution_recorder=TubularReactorRecorder(), + **kwargs + ) class LumpedRateModelWithoutPores(TubularReactorBase): @@ -1161,3 +1202,105 @@ def q(self, q): self._q = q self.parameters['q'] = q + + +class MCT(UnitBaseClass): + """Parameters for multi-channel transportmodel. + + Parameters + ---------- + length : UnsignedFloat + Length of column. + channel_cross_section_areas : List of unsinged floats. Lenght depends on nchannel. + Diameter of column. + axial_dispersion : UnsignedFloat + Dispersion rate of components in axial direction. + flow_direction : Switch + If 1: Forward flow. + If -1: Backwards flow. + c : List of unsigned floats. Length depends n_comp or n_comp*nchannel. + Initial concentration for components. + exchange_matrix : List of unsigned floats. Lenght depends on nchannel. + solution_recorder : MCTRecorder + Solution recorder for the unit operation. + n_channel : int number of channels + """ + has_ports = True + supports_bulk_reaction = True + + discretization_schemes = (MCTDiscretizationFV) + + length = UnsignedFloat() + channel_cross_section_areas = SizedList(size='nchannel') + axial_dispersion = UnsignedFloat() + flow_direction = Switch(valid=[-1, 1], default=1) + nchannel = UnsignedInteger() + + exchange_matrix = SizedNdArray(size=('nchannel', 'nchannel','n_comp')) + + _parameters = [ + 'length', + 'channel_cross_section_areas', + 'axial_dispersion', + 'flow_direction', + 'exchange_matrix', + 'nchannel' + ] + _section_dependent_parameters = \ + UnitBaseClass._section_dependent_parameters + \ + ['axial_dispersion', 'flow_direction'] + + _section_dependent_parameters = \ + UnitBaseClass._section_dependent_parameters + \ + [] + + c = SizedNdArray(size=('n_comp', 'nchannel'), default=0) + _initial_state = ['c'] + _parameters = _parameters + _initial_state + + def __init__(self, *args, nchannel, **kwargs): + discretization = MCTDiscretizationFV() + self._nchannel = nchannel + super().__init__( + *args, + discretization=discretization, + solution_recorder=MCTRecorder(), + **kwargs + ) + + @property + def nchannel(self): + return self._nchannel + + @nchannel.setter + def nchannel(self, nchannel): + self._nchannel = nchannel + + @property + def ports(self): + return [f"channel_{i}" for i in range(self.nchannel)] + + @property + def n_ports(self): + return self.nchannel + + @property + def volume(self): + """float: Combined Volumes of all channels. + + See Also + -------- + channel_cross_section_areas + + """ + return sum(self.channel_cross_section_areas) * self.length + + @property + def volume_liquid(self): + """float: Volume of the liquid phase. Equals the volume, since there is no solid phase.""" + return self.volume + + @property + def volume_solid(self): + """float: Volume of the solid phase. Equals zero, since there is no solid phase.""" + return 0 diff --git a/CADETProcess/simulationResults.py b/CADETProcess/simulationResults.py index 3088be13..96d82377 100644 --- a/CADETProcess/simulationResults.py +++ b/CADETProcess/simulationResults.py @@ -142,16 +142,31 @@ def solution(self): solution = Dict() for unit, solutions in self.solution_cycles.items(): - for sol, cycles in solutions.items(): - solution[unit][sol] = copy.deepcopy(cycles[0]) - solution_complete = cycles[0].solution_original - for i in range(1, self.n_cycles): - solution_complete = np.vstack(( - solution_complete, cycles[i].solution_original[1:] - )) - solution[unit][sol].time_original = time_complete - solution[unit][sol].solution_original = solution_complete - solution[unit][sol].reset() + + for sol, ports_cycles in solutions.items(): + if isinstance(ports_cycles, Dict): + ports = ports_cycles + for port, cycles in ports.items(): + solution[unit][sol][port] = copy.deepcopy(cycles[0]) + solution_complete = cycles[0].solution_original + for i in range(1, self.n_cycles): + solution_complete = np.vstack(( + solution_complete, cycles[i].solution_original[1:] + )) + solution[unit][sol][port].time_original = time_complete + solution[unit][sol][port].solution_original = solution_complete + solution[unit][sol][port].reset() + else: + cycles = ports_cycles + solution[unit][sol] = copy.deepcopy(cycles[0]) + solution_complete = cycles[0].solution_original + for i in range(1, self.n_cycles): + solution_complete = np.vstack(( + solution_complete, cycles[i].solution_original[1:] + )) + solution[unit][sol].time_original = time_complete + solution[unit][sol].solution_original = solution_complete + solution[unit][sol].reset() self._solution = solution @@ -166,18 +181,36 @@ def sensitivity(self): time_complete = self.time_complete sensitivity = Dict() - for unit, sensitivities in self.sensitivity_cycles.items(): - for sens_name, sensitivities in sensitivities.items(): - for sens, cycles in sensitivities.items(): - sensitivity[unit][sens_name][sens] = copy.deepcopy(cycles[0]) - sensitivity_complete = cycles[0].solution_original - for i in range(1, self.n_cycles): - sensitivity_complete = np.vstack(( - sensitivity_complete, cycles[i].solution_original[1:] - )) - sensitivity[unit][sens_name][sens].time_original = time_complete - sensitivity[unit][sens_name][sens].solution_original = sensitivity_complete - sensitivity[unit][sens_name][sens].reset() + for sens_name, sensitivities in self.sensitivity_cycles.items(): + + for unit, sensitivities in sensitivities.items(): + + for flow, ports_cycles in sensitivities.items(): + if isinstance(ports_cycles, Dict): + ports = ports_cycles + for port, cycles in ports.items(): + sensitivity[sens_name][unit][flow][port] = copy.deepcopy( + cycles[0]) + sensitivity_complete = cycles[0].solution_original + for i in range(1, self.n_cycles): + sensitivity_complete = np.vstack(( + sensitivity_complete, cycles[i].solution_original[1:] + )) + sensitivity[sens_name][unit][flow][port].time_original = time_complete + sensitivity[sens_name][unit][flow][port].solution_original = sensitivity_complete + sensitivity[sens_name][unit][flow][port].reset() + + else: + cycles = ports_cycles + sensitivity[sens_name][unit][flow] = copy.deepcopy(cycles[0]) + sensitivity_complete = cycles[0].solution_original + for i in range(1, self.n_cycles): + sensitivity_complete = np.vstack(( + sensitivity_complete, cycles[i].solution_original[1:] + )) + sensitivity[sens_name][unit][flow].time_original = time_complete + sensitivity[sens_name][unit][flow].solution_original = sensitivity_complete + sensitivity[sens_name][unit][flow].reset() self._sensitivity = sensitivity diff --git a/CADETProcess/simulator/cadetAdapter.py b/CADETProcess/simulator/cadetAdapter.py index 3af7ec3a..c6105d9c 100644 --- a/CADETProcess/simulator/cadetAdapter.py +++ b/CADETProcess/simulator/cadetAdapter.py @@ -1,3 +1,4 @@ +from CADETProcess.dataStructure import Structure, ParameterWrapper from collections import defaultdict from functools import wraps import os @@ -9,6 +10,7 @@ import time import tempfile import warnings +import re from addict import Dict import numpy as np @@ -447,6 +449,45 @@ def run(self, process, cadet=None, file_path=None): return results + def get_cadet_version(self) -> tuple[str, str]: + """ + Get version and branch name of the currently instanced CADET build. + + Returns + ------- + tuple[str, str] + The CADET version as x.x.x and the branch name. + + Raises + ------ + ValueError + If version and branch name cannot be found in the output string. + RuntimeError + If any unhandled event during running the subprocess occurs. + """ + try: + result = subprocess.run( + [self.cadet_path, '--version'], + check=True, + stdout=subprocess.PIPE, + stderr=subprocess.PIPE, + text=True + ) + version_output = result.stdout.strip() + + version_match = re.search( + r'cadet-cli version ([\d.]+) \(([^)]+)\)', version_output + ) + + if version_match: + cadet_version = version_match.group(1) + branch_name = version_match.group(2) + return cadet_version, branch_name + else: + raise ValueError("CADET version or branch name missing from output.") + except subprocess.CalledProcessError as e: + raise RuntimeError(f"Command execution failed: {e}") + def get_new_cadet_instance(self): cadet = CadetAPI() # Because the initialization in __init__ isn't guaranteed to be called in multiprocessing @@ -552,74 +593,105 @@ def get_simulation_results( try: solution = Dict() for unit in process.flow_sheet.units: + solution[unit.name] = defaultdict(list) + + port_flag = unit.has_ports + + if port_flag: + solution[unit.name]['inlet'] = defaultdict(list) + solution[unit.name]['outlet'] = defaultdict(list) + unit_index = self.get_unit_index(process, unit) unit_solution = cadet.root.output.solution[unit_index] + unit_coordinates = \ cadet.root.output.coordinates[unit_index].copy() particle_coordinates = \ unit_coordinates.pop('particle_coordinates_000', None) - flow_in = process.flow_rate_timelines[unit.name].total_in - flow_out = process.flow_rate_timelines[unit.name].total_out + for port in unit.ports: - start = 0 - for cycle in range(self.n_cycles): - end = start + len(time) + port_index = self.get_port_index(process.flow_sheet, unit, port) - if 'solution_inlet' in unit_solution.keys(): - sol_inlet = unit_solution.solution_inlet[start:end, :] - solution[unit.name]['inlet'].append( - SolutionIO( - unit.name, - unit.component_system, time, sol_inlet, - flow_in - ) - ) + flow_in = process.flow_rate_timelines[unit.name].total_in[port] + flow_out = process.flow_rate_timelines[unit.name].total_out[port] - if 'solution_outlet' in unit_solution.keys(): - sol_outlet = unit_solution.solution_outlet[start:end, :] - solution[unit.name]['outlet'].append( - SolutionIO( - unit.name, - unit.component_system, time, sol_outlet, - flow_out - ) - ) + start = 0 + for cycle in range(self.n_cycles): + end = start + len(time) + + if f'solution_inlet_port_{port_index:03d}' in unit_solution.keys(): + sol_inlet = unit_solution[f'solution_inlet_port_{port_index:03d}'][start:end,] + if port_flag: + solution[unit.name]['inlet'][port].append( + SolutionIO( + unit.name, + unit.component_system, time, sol_inlet, + flow_in + ) + ) + else: + solution[unit.name]['inlet'].append( + SolutionIO( + unit.name, + unit.component_system, time, sol_inlet, + flow_in + ) + ) - if 'solution_bulk' in unit_solution.keys(): - sol_bulk = unit_solution.solution_bulk[start:end, :] - solution[unit.name]['bulk'].append( - SolutionBulk( - unit.name, - unit.component_system, time, sol_bulk, - **unit_coordinates + if f'solution_outlet_port_{port_index:03d}' in unit_solution.keys(): + sol_outlet = unit_solution[f'solution_outlet_port_{port_index:03d}'][start:end, :] + if port_flag: + solution[unit.name]['outlet'][port].append( + SolutionIO( + unit.name, + unit.component_system, time, sol_outlet, + flow_out + ) + ) + else: + solution[unit.name]['outlet'].append( + SolutionIO( + unit.name, + unit.component_system, time, sol_outlet, + flow_out + ) + ) + + if 'solution_bulk' in unit_solution.keys(): + sol_bulk = unit_solution.solution_bulk[start:end, :] + solution[unit.name]['bulk'].append( + SolutionBulk( + unit.name, + unit.component_system, time, sol_bulk, + **unit_coordinates + ) ) - ) - if 'solution_particle' in unit_solution.keys(): - sol_particle = unit_solution.solution_particle[start:end, :] - solution[unit.name]['particle'].append( - SolutionParticle( - unit.name, - unit.component_system, time, sol_particle, - **unit_coordinates, - particle_coordinates=particle_coordinates + if 'solution_particle' in unit_solution.keys(): + sol_particle = unit_solution.solution_particle[start:end, :] + solution[unit.name]['particle'].append( + SolutionParticle( + unit.name, + unit.component_system, time, sol_particle, + **unit_coordinates, + particle_coordinates=particle_coordinates + ) ) - ) - if 'solution_solid' in unit_solution.keys(): - sol_solid = unit_solution.solution_solid[start:end, :] - solution[unit.name]['solid'].append( - SolutionSolid( - unit.name, - unit.component_system, - unit.binding_model.bound_states, - time, sol_solid, - **unit_coordinates, - particle_coordinates=particle_coordinates + if 'solution_solid' in unit_solution.keys(): + sol_solid = unit_solution.solution_solid[start:end, :] + solution[unit.name]['solid'].append( + SolutionSolid( + unit.name, + unit.component_system, + unit.binding_model.bound_states, + time, sol_solid, + **unit_coordinates, + particle_coordinates=particle_coordinates + ) ) - ) if 'solution_volume' in unit_solution.keys(): sol_volume = unit_solution.solution_volume[start:end, :] @@ -638,8 +710,17 @@ def get_simulation_results( sensitivity = Dict() for i, sens in enumerate(process.parameter_sensitivities): sens_index = f'param_{i:03d}' + for unit in process.flow_sheet.units: + sensitivity[sens.name][unit.name] = defaultdict(list) + + port_flag = unit.has_ports + + if port_flag: + sensitivity[sens.name][unit.name]['inlet'] = defaultdict(list) + sensitivity[sens.name][unit.name]['outlet'] = defaultdict(list) + unit_index = self.get_unit_index(process, unit) unit_sensitivity = cadet.root.output.sensitivity[sens_index][unit_index] unit_coordinates = \ @@ -647,77 +728,102 @@ def get_simulation_results( particle_coordinates = \ unit_coordinates.pop('particle_coordinates_000', None) - flow_in = process.flow_rate_timelines[unit.name].total_in - flow_out = process.flow_rate_timelines[unit.name].total_out - - for cycle in range(self.n_cycles): - start = cycle * len(time) - end = (cycle + 1) * len(time) - - if 'sens_inlet' in unit_sensitivity.keys(): - sens_inlet = unit_sensitivity.sens_inlet[start:end, :] - sensitivity[sens.name][unit.name]['inlet'].append( - SolutionIO( - unit.name, - unit.component_system, time, sens_inlet, - flow_in + for port in unit.ports: + + port_index = self.get_port_index(process.flow_sheet, unit, port) + + flow_in = process.flow_rate_timelines[unit.name].total_in[port] + flow_out = process.flow_rate_timelines[unit.name].total_out[port] + + start = 0 + for cycle in range(self.n_cycles): + end = start + len(time) + + if f'sens_inlet_port_{port_index:03d}' in unit_sensitivity.keys(): + sens_inlet = unit_sensitivity[f'sens_inlet_port_{port_index:03d}'][start:end, :] + if port_flag: + sensitivity[sens.name][unit.name]['inlet'][port].append( + SolutionIO( + unit.name, + unit.component_system, time, sens_inlet, + flow_in + ) + ) + + else: + sensitivity[sens.name][unit.name]['inlet'].append( + SolutionIO( + unit.name, + unit.component_system, time, sens_inlet, + flow_in + ) + ) + + if f'sens_outlet_port_{port_index:03d}' in unit_sensitivity.keys(): + sens_outlet = unit_sensitivity[f'sens_outlet_port_{port_index:03d}'][start:end, :] + if port_flag: + sensitivity[sens.name][unit.name]['outlet'][port].append( + SolutionIO( + unit.name, + unit.component_system, time, sens_outlet, + flow_out + ) + ) + else: + sensitivity[sens.name][unit.name]['outlet'].append( + SolutionIO( + unit.name, + unit.component_system, time, sens_outlet, + flow_out + ) + ) + + if 'sens_bulk' in unit_sensitivity.keys(): + sens_bulk = unit_sensitivity.sens_bulk[start:end, :] + sensitivity[sens.name][unit.name]['bulk'].append( + SolutionBulk( + unit.name, + unit.component_system, time, sens_bulk, + **unit_coordinates + ) ) - ) - if 'sens_outlet' in unit_sensitivity.keys(): - sens_outlet = unit_sensitivity.sens_outlet[start:end, :] - sensitivity[sens.name][unit.name]['outlet'].append( - SolutionIO( - unit.name, - unit.component_system, time, sens_outlet, - flow_out + if 'sens_particle' in unit_sensitivity.keys(): + sens_particle = unit_sensitivity.sens_particle[start:end, :] + sensitivity[sens.name][unit.name]['particle'].append( + SolutionParticle( + unit.name, + unit.component_system, time, sens_particle, + **unit_coordinates, + particle_coordinates=particle_coordinates + ) ) - ) - if 'sens_bulk' in unit_sensitivity.keys(): - sens_bulk = unit_sensitivity.sens_bulk[start:end, :] - sensitivity[sens.name][unit.name]['bulk'].append( - SolutionBulk( - unit.name, - unit.component_system, time, sens_bulk, - **unit_coordinates + if 'sens_solid' in unit_sensitivity.keys(): + sens_solid = unit_sensitivity.sens_solid[start:end, :] + sensitivity[sens.name][unit.name]['solid'].append( + SolutionSolid( + unit.name, + unit.component_system, + unit.binding_model.bound_states, + time, sens_solid, + **unit_coordinates, + particle_coordinates=particle_coordinates + ) ) - ) - if 'sens_particle' in unit_sensitivity.keys(): - sens_particle = unit_sensitivity.sens_particle[start:end, :] - sensitivity[sens.name][unit.name]['particle'].append( - SolutionParticle( - unit.name, - unit.component_system, time, sens_particle, - **unit_coordinates, - particle_coordinates=particle_coordinates + if 'sens_volume' in unit_sensitivity.keys(): + sens_volume = unit_sensitivity.sens_volume[start:end, :] + sensitivity[sens.name][unit.name]['volume'].append( + SolutionVolume( + unit.name, + unit.component_system, + time, + sens_volume + ) ) - ) - if 'sens_solid' in unit_sensitivity.keys(): - sens_solid = unit_sensitivity.sens_solid[start:end, :] - sensitivity[sens.name][unit.name]['solid'].append( - SolutionSolid( - unit.name, - unit.component_system, - unit.binding_model.bound_states, - time, sens_solid, - **unit_coordinates, - particle_coordinates=particle_coordinates - ) - ) - - if 'sens_volume' in unit_sensitivity.keys(): - sens_volume = unit_sensitivity.sens_volume[start:end, :] - sensitivity[sens.name][unit.name]['volume'].append( - SolutionVolume( - unit.name, - unit.component_system, - time, - sens_volume - ) - ) + start = end - 1 sensitivity = Dict(sensitivity) @@ -786,6 +892,8 @@ def get_model_connections(self, process): else: model_connections['CONNECTIONS_INCLUDE_DYNAMIC_FLOW'] = 1 + model_connections['CONNECTIONS_INCLUDE_PORTS'] = 1 + index = 0 section_states = process.flow_rate_section_states @@ -828,21 +936,31 @@ def cadet_connections(self, flow_rates, flow_sheet): for origin, unit_flow_rates in flow_rates.items(): origin = flow_sheet[origin] origin_index = flow_sheet.get_unit_index(origin) - for dest, flow_rate in unit_flow_rates.destinations.items(): - destination = flow_sheet[dest] - destination_index = flow_sheet.get_unit_index(destination) - if np.any(flow_rate): - table[enum] = [] - table[enum].append(int(origin_index)) - table[enum].append(int(destination_index)) - table[enum].append(-1) - table[enum].append(-1) - Q = flow_rate.tolist() - if self._force_constant_flow_rate: - table[enum] += [Q[0]] - else: - table[enum] += Q - enum += 1 + for origin_port, dest_dict in unit_flow_rates.destinations.items(): + for dest in dest_dict: + for dest_port, flow_rate in dest_dict[dest].items(): + destination = flow_sheet[dest] + destination_index = flow_sheet.get_unit_index(destination) + if np.any(flow_rate): + + origin_port_red = flow_sheet.get_port_index( + origin, origin_port) + dest_port_red = flow_sheet.get_port_index( + destination, dest_port) + + table[enum] = [] + table[enum].append(int(origin_index)) + table[enum].append(int(destination_index)) + table[enum].append(int(origin_port_red)) + table[enum].append(int(dest_port_red)) + table[enum].append(-1) + table[enum].append(-1) + Q = flow_rate.tolist() + if self._force_constant_flow_rate: + table[enum] += [Q[0]] + else: + table[enum] += Q + enum += 1 ls = [] for connection in table.values(): @@ -869,6 +987,24 @@ def get_unit_index(self, process, unit): index = process.flow_sheet.get_unit_index(unit) return f'unit_{index:03d}' + def get_port_index(self, flow_sheet, unit, port): + """Helper function for getting port index in CADET format xxx. + + Parameters + ---------- + port : string + Indexed port + unit : UnitOperation + port of unit + Returns + ------- + port_index : index + Return the port_index in CADET format xxx + + """ + + return flow_sheet.get_port_index(unit, port) + def get_model_units(self, process): """Config branches for all units /input/model/unit_000 ... unit_xxx. @@ -1198,7 +1334,6 @@ def __str__(self): return 'CADET' -from CADETProcess.dataStructure import Structure, ParameterWrapper class ModelSolverParameters(Structure): """Converter for model solver parameters from CADETProcess to CADET. @@ -1312,7 +1447,7 @@ class ModelSolverParameters(Structure): 'COL_LENGTH': 'length', 'CROSS_SECTION_AREA': 'cross_section_area', 'VELOCITY': 'flow_direction', - }, + }, 'fixed': { 'TOTAL_POROSITY': 1, }, @@ -1328,6 +1463,19 @@ class ModelSolverParameters(Structure): 'FLOWRATE_FILTER': 'flow_rate_filter', }, }, + 'MCT': { + 'name': 'MULTI_CHANNEL_TRANSPORT', + 'parameters': { + 'NCOMP': 'n_comp', + 'INIT_C': 'c', + 'COL_DISPERSION': 'axial_dispersion', + 'COL_LENGTH': 'length', + 'NCHANNEL': 'nchannel', + 'CHANNEL_CROSS_SECTION_AREAS': 'channel_cross_section_areas', + 'EXCHANGE_MATRIX': 'exchange_matrix', + 'VELOCITY': 'flow_direction', + }, + }, 'Inlet': { 'name': 'INLET', 'parameters': { @@ -1701,7 +1849,7 @@ class AdsorptionParameters(ParameterWrapper): 'mal_exponents_bulk_bwd': 'exponents_bwd', 'mal_kfwd_bulk': 'k_fwd', 'mal_kbwd_bulk': 'k_bwd', - } + } }, 'MassActionLawParticle': { 'name': 'MASS_ACTION_LAW', @@ -1881,11 +2029,12 @@ class ReturnParameters(Structure): write_solution_last = Bool(default=True) write_sens_last = Bool(default=True) split_components_data = Bool(default=False) - split_ports_data = Bool(default=False) + split_ports_data = Bool(default=True) + single_as_multi_port = Bool(default=True) _parameters = [ 'write_solution_times', 'write_solution_last', 'write_sens_last', - 'split_components_data', 'split_ports_data' + 'split_components_data', 'split_ports_data', 'single_as_multi_port', ] diff --git a/CADETProcess/solution.py b/CADETProcess/solution.py index 1ec57122..77dffa44 100755 --- a/CADETProcess/solution.py +++ b/CADETProcess/solution.py @@ -786,10 +786,11 @@ def __init__( self.component_system_original = component_system self.time_original = time + if axial_coordinates is not None and len(axial_coordinates) == 1: + axial_coordinates = None + self.axial_coordinates = axial_coordinates - # Account for dimension reduction in case of only one cell (e.g. LRMP) - if radial_coordinates is not None and len(radial_coordinates) == 1: - radial_coordinates = None + self.radial_coordinates = radial_coordinates self.solution_original = solution @@ -1058,14 +1059,21 @@ def __init__( particle_coordinates=None ): + if axial_coordinates is not None and len(axial_coordinates) == 1: + axial_coordinates = None + self.axial_coordinates = axial_coordinates # Account for dimension reduction in case of only one cell (e.g. LRMP) + if radial_coordinates is not None and len(radial_coordinates) == 1: radial_coordinates = None + self.radial_coordinates = radial_coordinates # Account for dimension reduction in case of only one cell (e.g. CSTR) + if particle_coordinates is not None and len(particle_coordinates) == 1: particle_coordinates = None + self.particle_coordinates = particle_coordinates super().__init__(name, component_system, time, solution) @@ -1146,7 +1154,6 @@ def _plot_1D( return ax - def _plot_2D(self, t, comp, vmax, ax=None): x = self.axial_coordinates y = self.particle_coordinates @@ -1228,6 +1235,8 @@ def __init__( self.bound_states = bound_states + if axial_coordinates is not None and len(axial_coordinates) == 1: + axial_coordinates = None self.axial_coordinates = axial_coordinates # Account for dimension reduction in case of only one cell (e.g. LRMP) if radial_coordinates is not None and len(radial_coordinates) == 1: @@ -1679,9 +1688,9 @@ def __init__(self, time, signal): if len(signal.shape) == 1: signal = np.array(signal, ndmin=2).transpose() self._solutions = [ - PchipInterpolator(time, signal[:, comp]) - for comp in range(signal.shape[1]) - ] + PchipInterpolator(time, signal[:, comp]) + for comp in range(signal.shape[1]) + ] self._derivatives = [signal.derivative() for signal in self._solutions] self._antiderivatives = [signal.antiderivative() for signal in self._solutions] diff --git a/tests/create_LWE.py b/tests/create_LWE.py new file mode 100644 index 00000000..4a875cbd --- /dev/null +++ b/tests/create_LWE.py @@ -0,0 +1,414 @@ +import numpy as np + +from CADETProcess.processModel import ( + ComponentSystem, FlowSheet, Process, + Inlet, Outlet, Cstr, GeneralRateModel, + TubularReactor, LumpedRateModelWithoutPores, + LumpedRateModelWithPores, MCT, StericMassAction +) + + +def create_lwe(unit_type: str = 'GeneralRateModel', **kwargs) -> Process: + """ + Create a process with the specified unit type and configuration. + + Parameters + ---------- + unit_type : str, optional + The type of unit operation, by default 'GeneralRateModel'. + **kwargs : dict + Additional parameters for configuring the unit operation. + + Returns + ------- + Process + The configured process. + """ + n_comp: int = kwargs.get('n_comp', 4) + component_system = ComponentSystem(n_comp) + + if unit_type == 'Cstr': + unit = configure_cstr(component_system, **kwargs) + elif unit_type == 'GeneralRateModel': + unit = configure_general_rate_model(component_system, **kwargs) + elif unit_type == 'TubularReactor': + unit = configure_tubular_reactor(component_system, **kwargs) + elif unit_type == 'LumpedRateModelWithoutPores': + unit = configure_lumped_rate_model_without_pores(component_system, **kwargs) + elif unit_type == 'LumpedRateModelWithPores': + unit = configure_lumped_rate_model_with_pores(component_system, **kwargs) + elif unit_type == 'MCT': + unit = configure_multichannel_transport_model(component_system, **kwargs) + else: + raise ValueError(f'Unknown unit operation type {unit_type}') + + flow_sheet = setup_flow_sheet(unit, component_system) + + process = Process(flow_sheet, 'process') + process.cycle_time = 120 * 60 + + c1_lwe = [[50.0], [0.0], [[100.0, 0.2]]] + cx_lwe = [[1.0], [0.0], [0.0]] + + process.add_event( + 'load', 'flow_sheet.inlet.c', + c1_lwe[0] + cx_lwe[0] * (n_comp - 1), 0 + ) + process.add_event( + 'wash', 'flow_sheet.inlet.c', + c1_lwe[1] + cx_lwe[1] * (n_comp - 1), 10 + ) + process.add_event( + 'elute', 'flow_sheet.inlet.c', + c1_lwe[2] + cx_lwe[2] * (n_comp - 1), 90 + ) + + return process + + +def setup_flow_sheet(unit, component_system: ComponentSystem) -> FlowSheet: + """ + Set up the flow sheet for the process. + + Parameters + ---------- + unit : UnitOperation + The unit operation to be added to the flow sheet. + component_system : ComponentSystem + The component system of the process. + + Returns + ------- + FlowSheet + The configured flow sheet. + """ + flow_sheet = FlowSheet(component_system) + inlet = Inlet(component_system, name='inlet') + inlet.flow_rate = 1.2e-3 + outlet = Outlet(component_system, name='outlet') + + flow_sheet.add_unit(inlet) + flow_sheet.add_unit(unit) + flow_sheet.add_unit(outlet) + + if unit.has_ports: + flow_sheet.add_connection(inlet, unit, destination_port='channel_0') + flow_sheet.add_connection(unit, outlet, origin_port='channel_0') + else: + flow_sheet.add_connection(inlet, unit) + flow_sheet.add_connection(unit, outlet) + + return flow_sheet + + +def configure_cstr(component_system: ComponentSystem, **kwargs) -> Cstr: + """ + Configure a continuous stirred-tank reactor (CSTR). + + Parameters + ---------- + component_system : ComponentSystem + The component system of the process. + **kwargs : dict + Additional parameters for configuring the CSTR. + + Returns + ------- + Cstr + The configured CSTR. + """ + cstr = Cstr(component_system, name='Cstr') + cstr.V = 1e-3 + cstr.porosity = 0.37 + (1.0 - 0.37) * 0.75 + + configure_solution_recorder(cstr, **kwargs) + configure_steric_mass_action(cstr, component_system, **kwargs) + + return cstr + + +def configure_general_rate_model(component_system: ComponentSystem, **kwargs) -> GeneralRateModel: + """ + Configure a general rate model. + + Parameters + ---------- + component_system : ComponentSystem + The component system of the process. + **kwargs : dict + Additional parameters for configuring the general rate model. + + Returns + ------- + GeneralRateModel + The configured general rate model. + """ + grm = GeneralRateModel(component_system, name='GeneralRateModel') + + grm.length = 0.014 + grm.diameter = 0.01 * 2 + grm.bed_porosity = 0.37 + grm.axial_dispersion = 5.75e-8 + grm.pore_diffusion = [7e-10, 6.07e-11, 6.07e-11, 6.07e-11] + + configure_solution_recorder(grm, **kwargs) + configure_discretization(grm, **kwargs) + configure_particles(grm, **kwargs) + configure_steric_mass_action(grm, component_system, **kwargs) + configure_film_diffusion(grm, component_system.n_comp) + configure_flow_direction(grm, **kwargs) + + return grm + + +def configure_tubular_reactor(component_system: ComponentSystem, **kwargs) -> TubularReactor: + """ + Configure a tubular reactor. + + Parameters + ---------- + component_system : ComponentSystem + The component system of the process. + **kwargs : dict + Additional parameters for configuring the tubular reactor. + + Returns + ------- + TubularReactor + The configured tubular reactor. + """ + tr = TubularReactor(component_system, name='TubularReactor') + + tr.length = 0.014 + tr.diameter = 0.01 * 2 + tr.axial_dispersion = 5.75e-8 + + configure_solution_recorder(tr, **kwargs) + configure_discretization(tr, **kwargs) + configure_flow_direction(tr, **kwargs) + + return tr + + +def configure_lumped_rate_model_without_pores(component_system: ComponentSystem, **kwargs) -> LumpedRateModelWithoutPores: + """ + Configure a lumped rate model without pores. + + Parameters + ---------- + component_system : ComponentSystem + The component system of the process. + **kwargs : dict + Additional parameters for configuring the lumped rate model. + + Returns + ------- + LumpedRateModelWithoutPores + The configured lumped rate model. + """ + lrm = LumpedRateModelWithoutPores( + component_system, name='LumpedRateModelWithoutPores' + ) + + lrm.length = 0.014 + lrm.diameter = 0.01 * 2 + lrm.total_porosity = 0.37 + (1.0 - 0.37) * 0.75 + lrm.axial_dispersion = 5.75e-8 + + configure_solution_recorder(lrm, **kwargs) + configure_discretization(lrm, **kwargs) + configure_steric_mass_action(lrm, component_system, **kwargs) + configure_flow_direction(lrm, **kwargs) + + return lrm + + +def configure_lumped_rate_model_with_pores(component_system: ComponentSystem, **kwargs) -> LumpedRateModelWithPores: + """ + Configure a lumped rate model with pores. + + Parameters + ---------- + component_system : ComponentSystem + The component system of the process. + **kwargs : dict + Additional parameters for configuring the lumped rate model. + + Returns + ------- + LumpedRateModelWithPores + The configured lumped rate model. + """ + lrmp = LumpedRateModelWithPores( + component_system, name='LumpedRateModelWithPores' + ) + + lrmp.length = 0.014 + lrmp.diameter = 0.01 * 2 + lrmp.bed_porosity = 0.37 + lrmp.axial_dispersion = 5.75e-8 + + configure_solution_recorder(lrmp, **kwargs) + configure_discretization(lrmp, **kwargs) + configure_particles(lrmp, **kwargs) + configure_steric_mass_action(lrmp, component_system, **kwargs) + configure_film_diffusion(lrmp, component_system.n_comp) + configure_flow_direction(lrmp, **kwargs) + + return lrmp + + +def configure_multichannel_transport_model(component_system: ComponentSystem, **kwargs) -> MCT: + """ + Configure a multichannel transport model. + + Parameters + ---------- + component_system : ComponentSystem + The component system of the process. + **kwargs : dict + Additional parameters for configuring the multichannel transport model. + + Returns + ------- + MCT + The configured multichannel transport model. + """ + mct = MCT(component_system, nchannel=3, name='MCT') + + mct.length = 0.014 + mct.channel_cross_section_areas = 3 * [2 * np.pi * (0.01 ** 2)] + mct.axial_dispersion = 5.75e-8 + + n_comp: int = component_system.n_comp + + mct.exchange_matrix = np.array([ + [n_comp * [0.0], n_comp * [0.001], n_comp * [0.0]], + [n_comp * [0.002], n_comp * [0.0], n_comp * [0.003]], + [n_comp * [0.0], n_comp * [0.0], n_comp * [0.0]] + ]) + + configure_solution_recorder(mct, **kwargs) + configure_discretization(mct, **kwargs) + configure_flow_direction(mct, **kwargs) + + return mct + + +def configure_discretization(unit_operation, **kwargs) -> None: + """ + Configure the discretization settings for a unit operation. + + Parameters + ---------- + unit_operation : UnitOperation + The unit operation to be configured. + **kwargs : dict + Additional parameters for configuring the discretization. + """ + n_col: int = kwargs.get('n_col', 100) + n_par: int = kwargs.get('n_par', 2) + ad_jacobian: bool = kwargs.get('ad_jacobian', False) + + if 'npar' in unit_operation.discretization.parameters: + unit_operation.discretization.npar = n_par + unit_operation.discretization.par_disc_type = 'EQUIDISTANT_PAR' + + unit_operation.discretization.ncol = n_col + unit_operation.discretization.use_analytic_jacobian = not ad_jacobian + + +def configure_particles(unit_operation) -> None: + """ + Configure the particle settings for a unit operation. + + Parameters + ---------- + unit_operation : UnitOperation + The unit operation to be configured. + """ + par_radius = 4.5e-5 + par_porosity = 0.75 + + unit_operation.particle_radius = par_radius + unit_operation.particle_porosity = par_porosity + unit_operation.discretization.par_geom = 'SPHERE' + + +def configure_steric_mass_action(unit_operation, component_system, **kwargs) -> None: + """ + Configure the steric mass action binding model for a unit operation. + + Parameters + ---------- + unit_operation : UnitOperation + The unit operation to be configured. + component_system : ComponentSystem + The component system of the process. + **kwargs : dict + Additional parameters for configuring the steric mass action binding model. + """ + is_kinetic = kwargs.get('is_kinetic', True) + + kA = 35.5 + kD = 1000.0 + nu = 4.7 + sigma = 11.83 + sma_lambda = 1.2e3 + + binding_model = StericMassAction(component_system) + + binding_model.is_kinetic = is_kinetic + binding_model.n_binding_sites = 1 + binding_model.adsorption_rate = kA + binding_model.desorption_rate = kD + binding_model.characteristic_charge = nu + binding_model.steric_factor = sigma + binding_model.capacity = sma_lambda + + unit_operation.binding_model = binding_model + + +def configure_film_diffusion(unit_operation, n_comp) -> None: + """ + Configure the film diffusion settings for a unit operation. + + Parameters + ---------- + unit_operation : UnitOperation + The unit operation to be configured. + n_comp : int + The number of components in the process. + """ + unit_operation.film_diffusion = [6.9e-6] * n_comp + + +def configure_flow_direction(unit_operation, **kwargs) -> None: + """ + Configure the flow direction for a unit operation. + + Parameters + ---------- + unit_operation : UnitOperation + The unit operation to be configured. + **kwargs : dict + Additional parameters for configuring the flow direction. + """ + reverse_flow = kwargs.get('reverse_flow', False) + unit_operation.flow_direction = -1 if reverse_flow else 1 + + +def configure_solution_recorder(unit_operation, **kwargs) -> None: + """ + Configure the solution recorder for a unit operation. + + Parameters + ---------- + unit_operation : UnitOperation + The unit operation to be configured. + **kwargs : dict + Additional parameters for configuring the solution recorder. + """ + for write_solution, value in unit_operation.solution_recorder.parameters.items(): + if value is False: + unit_operation.solution_recorder.parameters[write_solution] = True diff --git a/tests/test_cadet_adapter.py b/tests/test_cadet_adapter.py index b1ac69b1..17559e9b 100644 --- a/tests/test_cadet_adapter.py +++ b/tests/test_cadet_adapter.py @@ -2,8 +2,19 @@ import platform import shutil import unittest +import warnings +from typing import Optional +import pytest +import numpy as np +import numpy.testing as npt +from tests.create_LWE import create_lwe + +from CADETProcess import CADETProcessError +from CADETProcess.processModel import Process +from CADETProcess import SimulationResults from CADETProcess.simulator import Cadet +from CADETProcess.processModel.discretization import NoDiscretization def detect_cadet(): @@ -70,5 +81,596 @@ def tearDown(self): shutil.rmtree('./tmp', ignore_errors=True) -if __name__ == '__main__': +# TODO: Update when MCT is included in CADET-Python +# Include path to latest CADET-Core version here to test the MCT +cadet_install_path = None +process_simulator = Cadet(install_path=cadet_install_path) +cadet_version, branch_name = process_simulator.get_cadet_version() + +if cadet_version < '5.0.0' and branch_name == 'GITDIR-NOTFOUND branch': + unit_types = [ + 'Cstr', 'GeneralRateModel', 'TubularReactor', + 'LumpedRateModelWithoutPores', 'LumpedRateModelWithPores' + ] +else: + unit_types = [ + 'Cstr', 'GeneralRateModel', 'TubularReactor', + 'LumpedRateModelWithoutPores', 'LumpedRateModelWithPores', 'MCT' + ] + + +def run_simulation( + process: Process, + install_path: Optional[str] = None + ) -> SimulationResults: + """ + Run the CADET simulation for the given process and handle potential issues. + + Parameters + ---------- + process : Process + The process to simulate. + cadet_install_path : str, optional + The path to the CADET installation. + + Returns + ------- + SimulationResults + The results of the simulation. + + Raises + ------ + CADETProcessError + If the simulation fails with an error. + """ + try: + process_simulator = Cadet(install_path) + + cadet_version, branch_name = process_simulator.get_cadet_version() + if cadet_version < '5.0.0' and branch_name == 'GITDIR-NOTFOUND branch': + warnings.warn( + f"Your current CADET-Core version ({cadet_version}) does not " + "support all unit operations of this CADET-Process version. " + "Please update to the latest CADET-Core version from " + "https://github.com/cadet/CADET-Core for full compatibility.", + UserWarning + ) + + simulation_results = process_simulator.simulate(process) + + if not simulation_results.exit_flag == 0: + raise CADETProcessError( + f"LWE simulation failed with {simulation_results.exit_message}." + ) + + return simulation_results + + except Exception as e: + raise CADETProcessError(f"CADET simulation failed: {e}.") from e + + +@pytest.fixture(scope="class", params=unit_types) +def process(request: pytest.FixtureRequest): + """ + Fixture to set up the process for each unit type without running the simulation. + """ + unit_type = request.param + process = create_lwe(unit_type) + return process + + +@pytest.fixture(scope="class", params=unit_types) +def simulation_results(request: pytest.FixtureRequest): + """ + Fixture to set up the simulation for each unit type. + """ + unit_type = request.param + process = create_lwe(unit_type) + simulation_results = run_simulation(process, install_path=cadet_install_path) + return simulation_results + + +@pytest.mark.parametrize("process", unit_types, indirect=True) +class TestProcessWithLWE: + + def return_process_config(self, process: Process) -> dict: + """ + Returns the process configuration. + + Parameters + ---------- + process : Process + The process object. + + Returns + ------- + dict + The configuration of the process. + """ + process_simulator = Cadet(install_path=cadet_install_path) + process_config = process_simulator.get_process_config(process).input + return process_config + + def test_model_config(self, process: Process): + """ + Test the model configuration for various unit types in the process. + + Parameters + ---------- + process : Process + The process object. + """ + process_config = self.return_process_config(process) + + n_comp = process.component_system.n_comp + unit = process.flow_sheet.units[1] + + model_config = process_config.model + input_config = model_config.unit_000 + output_config = model_config.unit_002 + unit_config = model_config.unit_001 + + # ASSERT INPUT CONFIGURATION + c1_lwe = [[50.0], [0.0], [[100.0, 0.2]]] + cx_lwe = [[1.0], [0.0], [0.0]] + + expected_input_config = { + 'UNIT_TYPE': 'INLET', + 'NCOMP': n_comp, + 'INLET_TYPE': 'PIECEWISE_CUBIC_POLY', + 'discretization': { + 'nbound': n_comp * [0] + }, + 'sec_000': { + 'const_coeff': np.array(c1_lwe[0] + cx_lwe[0] * (n_comp - 1)), + 'lin_coeff': np.array([0.] * n_comp), + 'quad_coeff': np.array([0.] * n_comp), + 'cube_coeff': np.array([0.] * n_comp) + }, + 'sec_001': { + 'const_coeff': np.array(c1_lwe[1] + cx_lwe[1] * (n_comp - 1)), + 'lin_coeff': np.array([0.] * n_comp), + 'quad_coeff': np.array([0.] * n_comp), + 'cube_coeff': np.array([0.] * n_comp) + }, + 'sec_002': { + 'const_coeff': np.array([c1_lwe[2][0][0]] + cx_lwe[2] * (n_comp - 1)), + 'lin_coeff': np.array([c1_lwe[2][0][1]] + cx_lwe[2] * (n_comp - 1)), + 'quad_coeff': np.array([0.] * n_comp), + 'cube_coeff': np.array([0.] * n_comp) + } + } + + npt.assert_equal(input_config, expected_input_config) + + # ASSERT OUTPUT CONFIGURATION + expected_output_config = { + 'UNIT_TYPE': 'OUTLET', + 'NCOMP': n_comp, + 'discretization': { + 'nbound': n_comp * [0] + } + } + + npt.assert_equal(output_config, expected_output_config) + + # ASSERT MODEL CONFIGURATION + assert unit_config.NCOMP == n_comp + assert model_config.nunits == 3 + + if unit.name == 'Cstr': + self.check_cstr(unit, unit_config) + elif unit.name == 'GeneralRateModel': + self.check_general_rate_model(unit, unit_config) + elif unit.name == 'TubularReactor': + self.check_tubular_reactor(unit, unit_config) + elif unit.name == 'LumpedRateModelWithoutPores': + self.check_lumped_rate_model_without_pores(unit, unit_config) + elif unit.name == 'LumpedRateModelWithPores': + self.check_lumped_rate_model_with_pores(unit, unit_config) + elif unit.name == 'MCT': + self.check_mct(unit, unit_config) + + def check_cstr(self, unit, unit_config): + """ + Check the configuration for a CSTR unit. + + Parameters + ---------- + unit : Unit + The unit object. + unit_config : dict + The configuration of the unit. + """ + n_comp = unit.component_system.n_comp + + assert unit_config.UNIT_TYPE == 'CSTR' + assert unit_config.INIT_Q == n_comp * [0] + assert unit_config.INIT_C == n_comp * [0] + assert unit_config.INIT_VOLUME == 0.001 + assert unit_config.POROSITY == 0.8425 + assert unit_config.FLOWRATE_FILTER == 0.0 + assert unit_config.nbound == [1, 1, 1, 1] + + self.check_adsorption_config(unit, unit_config) + + def check_general_rate_model(self, unit, unit_config): + """ + Check the configuration for a General Rate Model unit. + + Parameters + ---------- + unit : Unit + The unit object. + unit_config : dict + The configuration of the unit. + """ + n_comp = unit.component_system.n_comp + + assert unit_config.UNIT_TYPE == 'GENERAL_RATE_MODEL' + assert unit_config.INIT_Q == n_comp * [0] + assert unit_config.INIT_C == n_comp * [0] + assert unit_config.INIT_CP == n_comp * [0] + assert unit_config.VELOCITY == unit.flow_direction + assert unit_config.COL_DISPERSION == 5.75e-08 + assert unit_config.CROSS_SECTION_AREA == np.pi * 0.01 ** 2 + assert unit_config.COL_LENGTH == 0.014 + assert unit_config.COL_POROSITY == 0.37 + assert unit_config.FILM_DIFFUSION == [6.9e-6] * n_comp + + self.check_particle_config(unit_config) + self.check_adsorption_config(unit, unit_config) + self.check_discretization(unit, unit_config) + + def check_tubular_reactor(self, unit, unit_config): + """ + Check the configuration for a Tubular Reactor unit. + + Parameters + ---------- + unit : Unit + The unit object. + unit_config : dict + The configuration of the unit. + """ + n_comp = unit.component_system.n_comp + + assert unit_config.UNIT_TYPE == 'LUMPED_RATE_MODEL_WITHOUT_PORES' + assert unit_config.INIT_C == n_comp * [0] + assert unit_config.VELOCITY == unit.flow_direction + assert unit_config.COL_DISPERSION == 5.75e-08 + assert unit_config.CROSS_SECTION_AREA == np.pi * 0.01 ** 2 + assert unit_config.COL_LENGTH == 0.014 + assert unit_config.TOTAL_POROSITY == 1 + + self.check_discretization(unit, unit_config) + + def check_lumped_rate_model_without_pores(self, unit, unit_config): + """ + Check the configuration for a Lumped Rate Model Without Pores unit. + + Parameters + ---------- + unit : Unit + The unit object. + unit_config : dict + The configuration of the unit. + """ + n_comp = unit.component_system.n_comp + + assert unit_config.UNIT_TYPE == 'LUMPED_RATE_MODEL_WITHOUT_PORES' + assert unit_config.INIT_C == n_comp * [0] + assert unit_config.VELOCITY == unit.flow_direction + assert unit_config.COL_DISPERSION == 5.75e-08 + assert unit_config.CROSS_SECTION_AREA == np.pi * 0.01 ** 2 + assert unit_config.COL_LENGTH == 0.014 + assert unit_config.TOTAL_POROSITY == 0.8425 + + self.check_adsorption_config(unit, unit_config) + self.check_discretization(unit, unit_config) + + def check_lumped_rate_model_with_pores(self, unit, unit_config): + """ + Check the configuration for a Lumped Rate Model With Pores unit. + + Parameters + ---------- + unit : Unit + The unit object. + unit_config : dict + The configuration of the unit. + """ + n_comp = unit.component_system.n_comp + + assert unit_config.UNIT_TYPE == 'LUMPED_RATE_MODEL_WITH_PORES' + assert unit_config.INIT_C == n_comp * [0] + assert unit_config.INIT_Q == n_comp * [0] + assert unit_config.INIT_CP == n_comp * [0] + assert unit_config.VELOCITY == unit.flow_direction + assert unit_config.COL_DISPERSION == 5.75e-08 + assert unit_config.CROSS_SECTION_AREA == np.pi * 0.01 ** 2 + assert unit_config.COL_LENGTH == 0.014 + assert unit_config.COL_POROSITY == 0.37 + assert unit_config.FILM_DIFFUSION == [6.9e-6] * n_comp + + self.check_adsorption_config(unit, unit_config) + self.check_discretization(unit, unit_config) + self.check_particle_config(unit_config) + + def check_mct(self, unit, unit_config): + """ + Check the configuration for a Multi-Channel Transport unit. + + Parameters + ---------- + unit : Unit + The unit object. + unit_config : dict + The configuration of the unit. + """ + n_comp = unit.component_system.n_comp + n_channel = unit.nchannel + + assert unit_config.UNIT_TYPE == 'MULTI_CHANNEL_TRANSPORT' + npt.assert_equal(unit_config.INIT_C, n_comp * [(n_channel * [0])]) + assert unit_config.VELOCITY == unit.flow_direction + npt.assert_equal( + unit_config.EXCHANGE_MATRIX, + np.array([ + [n_comp * [0.0], n_comp * [0.001], n_comp * [0.0]], + [n_comp * [0.002], n_comp * [0.0], n_comp * [0.003]], + [n_comp * [0.0], n_comp * [0.0], n_comp * [0.0]] + ]) + ) + assert unit_config.COL_DISPERSION == 5.75e-08 + assert unit_config.NCHANNEL == 3 + assert unit_config.CHANNEL_CROSS_SECTION_AREAS == 3 * [2 * np.pi * (0.01 ** 2)] + + self.check_discretization(unit, unit_config) + + def check_adsorption_config(self, unit, unit_config): + """ + Check the adsorption configuration. + + Parameters + ---------- + unit : Unit + The unit object. + unit_config : dict + The configuration of the unit. + """ + n_comp = unit.component_system.n_comp + + expected_adsorption_config = { + 'ADSORPTION_MODEL': 'STERIC_MASS_ACTION', + 'IS_KINETIC': True, + 'SMA_KA': n_comp * [35.5], + 'SMA_KD': n_comp * [1000.0], + 'SMA_LAMBDA': 1200.0, + 'SMA_NU': n_comp * [4.7], + 'SMA_SIGMA': n_comp * [11.83], + 'SMA_REFC0': 1.0, + 'SMA_REFQ': 1.0 + } + + assert unit_config.adsorption_model == 'STERIC_MASS_ACTION' + npt.assert_equal(unit_config.adsorption, expected_adsorption_config) + + def check_particle_config(self, unit_config): + """ + Check the particle configuration. + + Parameters + ---------- + unit_config : dict + The configuration of the unit. + """ + assert unit_config.PAR_POROSITY == 0.75 + assert unit_config.PAR_RADIUS == 4.5e-05 + assert unit_config.discretization.par_geom == 'SPHERE' + + def check_discretization(self, unit, unit_config): + """ + Check the discretization configuration. + + Parameters + ---------- + unit : Unit + The unit object. + unit_config : dict + The configuration of the unit. + """ + assert unit_config.discretization.ncol == unit.discretization.ncol + assert unit_config.discretization.use_analytic_jacobian == unit.discretization.use_analytic_jacobian + assert unit_config.discretization.reconstruction == 'WENO' + + npt.assert_equal( + unit_config.discretization.weno, { + 'boundary_model': 0, + 'weno_eps': 1e-10, + 'weno_order': 3 + } + ) + + npt.assert_equal( + unit_config.discretization.consistency_solver, { + 'solver_name': 'LEVMAR', + 'init_damping': 0.01, + 'min_damping': 0.0001, + 'max_iterations': 50, + 'subsolvers': 'LEVMAR' + } + ) + + if 'spatial_method' in unit.discretization.parameters: + assert unit_config.discretization.spatial_method == unit.discretization.spatial_method + + if 'npar' in unit.discretization.parameters: + assert unit_config.discretization.npar == unit.discretization.npar + assert unit_config.discretization.par_disc_type == 'EQUIDISTANT_PAR' + + def test_solver_config(self, process: Process): + """ + Test the solver configuration for the process. + + Parameters + ---------- + process : Process + The process object. + """ + process_config = self.return_process_config(process) + solver_config = process_config.solver + + expected_solver_config = { + 'nthreads': 1, + 'consistent_init_mode': 1, + 'consistent_init_mode_sens': 1, + 'user_solution_times': np.arange(0.0, 120 * 60 + 1), + 'sections': { + 'nsec': 3, + 'section_times': [0.0, 10.0, 90.0, 7200.0], + 'section_continuity': [0, 0] + }, + 'time_integrator': { + 'abstol': 1e-08, + 'algtol': 1e-12, + 'reltol': 1e-06, + 'reltol_sens': 1e-12, + 'init_step_size': 1e-06, + 'max_steps': 1000000, + 'max_step_size': 0.0, + 'errortest_sens': False, + 'max_newton_iter': 1000000, + 'max_errtest_fail': 1000000, + 'max_convtest_fail': 1000000, + 'max_newton_iter_sens': 1000000 + } + } + + npt.assert_equal(solver_config, expected_solver_config) + + def test_return_config(self, process: Process): + """ + Test the return configuration for the process. + + Parameters + ---------- + process : Process + The process object. + """ + process_config = self.return_process_config(process) + return_config = process_config['return'] + + # Assert that all values in return_config.unit_001 (model unit operation) are True + for key, value in return_config['unit_001'].items(): + assert value, f"The value for key '{key}' is not True. Found: {value}" + + def test_sensitivity_config(self, process: Process): + """ + Test the sensitivity configuration for the process. + + Parameters + ---------- + process : Process + The process object. + """ + process_config = self.return_process_config(process) + sensitivity_config = process_config.sensitivity + + expected_sensitivity_config = {'sens_method': 'ad1', 'nsens': 0} + npt.assert_equal(sensitivity_config, expected_sensitivity_config) + + +@pytest.mark.parametrize("simulation_results", unit_types, indirect=True) +class TestResultsWithLWE: + def test_trigger_simulation(self, simulation_results): + """ + Test to trigger the simulation. + """ + simulation_results = simulation_results + assert simulation_results is not None + + def test_compare_solution_shape(self, simulation_results): + """ + Compare the dimensions of the solution object against the expected solution shape. + """ + simulation_results = simulation_results + process = simulation_results.process + unit = process.flow_sheet.units[1] + + # for units without ports + if not unit.has_ports: + # assert solution inlet has shape (t, n_comp) + assert simulation_results.solution[unit.name].inlet.solution_shape == ( + int(process.cycle_time+1), process.component_system.n_comp + ) + # assert solution outlet has shape (t, n_comp) + assert simulation_results.solution[unit.name].outlet.solution_shape == ( + int(process.cycle_time+1), process.component_system.n_comp + ) + # assert solution bulk has shape (t, n_col, n_comp) + if not isinstance(unit.discretization, NoDiscretization): + assert simulation_results.solution[unit.name].bulk.solution_shape == ( + int(process.cycle_time+1), + unit.discretization.ncol, + process.component_system.n_comp + ) + + # for units with ports + else: + # assert solution inlet is given for each port + assert len(simulation_results.solution[unit.name].inlet) == unit.n_ports + # assert solution for channel 0 has shape (t, n_comp) + assert simulation_results.solution[unit.name].inlet.channel_0.solution_shape == ( + int(process.cycle_time+1), process.component_system.n_comp + ) + # assert solution bulk has shape (t, n_col, n_ports, n_comp) + assert simulation_results.solution[unit.name].bulk.solution_shape == ( + int(process.cycle_time+1), + unit.discretization.ncol, + unit.n_ports, + process.component_system.n_comp + ) + + # for units with particles + if unit.supports_binding and not isinstance(unit.discretization, NoDiscretization): + # for units with solid phase and particle discretization + if 'npar' in unit.discretization.parameters: + # assert solution solid has shape (t, n_col, n_par, n_comp) + assert simulation_results.solution[unit.name].solid.solution_shape == ( + int(process.cycle_time+1), + unit.discretization.ncol, + unit.discretization.npar, + process.component_system.n_comp + ) + # for units with solid phase and without particle discretization + else: + # assert solution solid has shape (t, ncol, n_comp) + assert simulation_results.solution[unit.name].solid.solution_shape == ( + int(process.cycle_time+1), + unit.discretization.ncol, + process.component_system.n_comp + ) + + # for units with particle mobile phase and particle discretization + if unit.supports_particle_reaction and unit.name != "LumpedRateModelWithoutPores": + # assert soluction particle has shape (t, n_col, n_par, n_comp) + if 'npar' in unit.discretization.parameters: + assert simulation_results.solution[unit.name].particle.solution_shape == ( + int(process.cycle_time+1), + unit.discretization.ncol, + unit.discretization.npar, + process.component_system.n_comp + ) + # for units with particle mobiles phase and particle discretization + else: + # assert solution particle has shape (t, n_col, n_par, n_comp) + assert simulation_results.solution[unit.name].particle.solution_shape == ( + int(process.cycle_time+1), + unit.discretization.ncol, + process.component_system.n_comp + ) + + +if __name__ == "__main__": unittest.main() diff --git a/tests/test_cadet_reactions.py b/tests/test_cadet_reactions.py index 241016b1..8ae00175 100644 --- a/tests/test_cadet_reactions.py +++ b/tests/test_cadet_reactions.py @@ -16,7 +16,7 @@ from CADETProcess.simulator import Cadet -from test_cadet_adapter import found_cadet +from tests.test_cadet_adapter import found_cadet def setup_process(unit_type): diff --git a/tests/test_carousel.py b/tests/test_carousel.py index 2a6ad922..6e643ad9 100644 --- a/tests/test_carousel.py +++ b/tests/test_carousel.py @@ -526,19 +526,19 @@ def test_flow_rates(self): column_0 = process.flow_rate_timelines['column_0'] column_1 = process.flow_rate_timelines['column_1'] - flow_rate = serial_inlet.total_in.value(0) + flow_rate = serial_inlet.total_in[None].value(0) flow_rate_expected = 2e-7 np.testing.assert_almost_equal(flow_rate, flow_rate_expected) - flow_rate = serial_outlet.total_in.value(0) + flow_rate = serial_outlet.total_in[None].value(0) flow_rate_expected = 2e-7 np.testing.assert_almost_equal(flow_rate, flow_rate_expected) - flow_rate = column_0.total_in.value(0) + flow_rate = column_0.total_in[None].value(0) flow_rate_expected = 2e-7 np.testing.assert_almost_equal(flow_rate, flow_rate_expected) - flow_rate = column_1.total_in.value(0) + flow_rate = column_1.total_in[None].value(0) flow_rate_expected = 2e-7 np.testing.assert_almost_equal(flow_rate, flow_rate_expected) @@ -550,19 +550,19 @@ def test_flow_rates(self): column_0 = process.flow_rate_timelines['column_0'] column_1 = process.flow_rate_timelines['column_1'] - flow_rate = parallel_inlet.total_in.value(0) + flow_rate = parallel_inlet.total_in[None].value(0) flow_rate_expected = 2e-7 np.testing.assert_almost_equal(flow_rate, flow_rate_expected) - flow_rate = parallel_inlet.total_in.value(0) + flow_rate = parallel_inlet.total_in[None].value(0) flow_rate_expected = 2e-7 np.testing.assert_almost_equal(flow_rate, flow_rate_expected) - flow_rate = column_0 .total_in.value(0) + flow_rate = column_0 .total_in[None].value(0) flow_rate_expected = 1e-7 np.testing.assert_almost_equal(flow_rate, flow_rate_expected) - flow_rate = column_1.total_in.value(0) + flow_rate = column_1.total_in[None].value(0) flow_rate_expected = 1e-7 np.testing.assert_almost_equal(flow_rate, flow_rate_expected) @@ -575,31 +575,31 @@ def test_flow_rates(self): column_0 = process.flow_rate_timelines['column_0'] column_2 = process.flow_rate_timelines['column_2'] - flow_rate = serial_inlet.total_in.value(0) + flow_rate = serial_inlet.total_in[None].value(0) flow_rate_expected = 2e-7 np.testing.assert_almost_equal(flow_rate, flow_rate_expected) - flow_rate = parallel_inlet.total_in.value(0) + flow_rate = parallel_inlet.total_in[None].value(0) flow_rate_expected = 3e-7 np.testing.assert_almost_equal(flow_rate, flow_rate_expected) # Initial state t = 0 - flow_rate = column_0.total_in.value(t) + flow_rate = column_0.total_in[None].value(t) flow_rate_expected = 2e-7 np.testing.assert_almost_equal(flow_rate, flow_rate_expected) - flow_rate = column_2.total_in.value(t) + flow_rate = column_2.total_in[None].value(t) flow_rate_expected = 1.5e-7 np.testing.assert_almost_equal(flow_rate, flow_rate_expected) # First position t = builder.switch_time - flow_rate = column_0.total_in.value(t) + flow_rate = column_0.total_in[None].value(t) flow_rate_expected = 1.5e-7 np.testing.assert_almost_equal(flow_rate, flow_rate_expected) - flow_rate = column_2.total_in.value(t) + flow_rate = column_2.total_in[None].value(t) flow_rate_expected = 2e-7 np.testing.assert_almost_equal(flow_rate, flow_rate_expected) diff --git a/tests/test_compartment.py b/tests/test_compartment.py index de5e6120..62a4074e 100644 --- a/tests/test_compartment.py +++ b/tests/test_compartment.py @@ -66,44 +66,132 @@ def test_complex(self): def test_connections(self): flow_rates_expected = { 'compartment_0': { - 'total_in': np.array([1., 0., 0., 0.]), - 'total_out': np.array([1., 0., 0., 0.]), + 'total_in': { + None: np.array([1., 0., 0., 0.]) + }, + 'total_out': { + None: np.array([1., 0., 0., 0.]) + }, 'origins': { - 'compartment_1': np.array([0.1, 0., 0., 0.]), - 'compartment_2': np.array([0.2, 0., 0., 0.]), - 'compartment_3': np.array([0.3, 0., 0., 0.]), - 'compartment_4': np.array([0.4, 0., 0., 0.]) + None: { + 'compartment_1': { + None: np.array([0.1, 0., 0., 0.]) + }, + 'compartment_2': { + None: np.array([0.2, 0., 0., 0.]) + }, + 'compartment_3': { + None: np.array([0.3, 0., 0., 0.]) + }, + 'compartment_4': { + None: np.array([0.4, 0., 0., 0.]) + } + } }, 'destinations': { - 'compartment_1': np.array([0.1, 0., 0., 0.]), - 'compartment_2': np.array([0.2, 0., 0., 0.]), - 'compartment_3': np.array([0.3, 0., 0., 0.]), - 'compartment_4': np.array([0.4, 0., 0., 0.]) + None: { + 'compartment_1': { + None: np.array([0.1, 0., 0., 0.]) + }, + 'compartment_2': { + None: np.array([0.2, 0., 0., 0.]) + }, + 'compartment_3': { + None: np.array([0.3, 0., 0., 0.]) + }, + 'compartment_4': { + None: np.array([0.4, 0., 0., 0.]) + } + } }, }, 'compartment_1': { - 'total_in': np.array([0.1, 0., 0., 0.]), - 'total_out': np.array([0.1, 0., 0., 0.]), - 'origins': {'compartment_0': np.array([0.1, 0., 0., 0.])}, - 'destinations': {'compartment_0': np.array([0.1, 0., 0., 0.])}, + 'total_in': { + None: np.array([0.1, 0., 0., 0.]) + }, + 'total_out': { + None: np.array([0.1, 0., 0., 0.]) + }, + 'origins': { + None: { + 'compartment_0': { + None: np.array([0.1, 0., 0., 0.]) + } + } + }, + 'destinations': { + None: { + 'compartment_0': { + None: np.array([0.1, 0., 0., 0.]) + } + } + }, }, 'compartment_2': { - 'total_in': np.array([0.2, 0., 0., 0.]), - 'total_out': np.array([0.2, 0., 0., 0.]), - 'origins': {'compartment_0': np.array([0.2, 0., 0., 0.])}, - 'destinations': {'compartment_0': np.array([0.2, 0., 0., 0.])}, + 'total_in': { + None: np.array([0.2, 0., 0., 0.]) + }, + 'total_out': { + None: np.array([0.2, 0., 0., 0.]) + }, + 'origins': { + None: { + 'compartment_0': { + None: np.array([0.2, 0., 0., 0.]) + } + } + }, + 'destinations': { + None: { + 'compartment_0': { + None: np.array([0.2, 0., 0., 0.]) + } + } + } }, 'compartment_3': { - 'total_in': np.array([0.3, 0., 0., 0.]), - 'total_out': np.array([0.3, 0., 0., 0.]), - 'origins': {'compartment_0': np.array([0.3, 0., 0., 0.])}, - 'destinations': {'compartment_0': np.array([0.3, 0., 0., 0.])}, + 'total_in': { + None: np.array([0.3, 0., 0., 0.]) + }, + 'total_out': { + None: np.array([0.3, 0., 0., 0.]) + }, + 'origins': { + None: { + 'compartment_0': { + None: np.array([0.3, 0., 0., 0.]) + } + } + }, + 'destinations': { + None: { + 'compartment_0': { + None: np.array([0.3, 0., 0., 0.]) + } + } + } }, 'compartment_4': { - 'total_in': np.array([0.4, 0., 0., 0.]), - 'total_out': np.array([0.4, 0., 0., 0.]), - 'origins': {'compartment_0': np.array([0.4, 0., 0., 0.])}, - 'destinations': {'compartment_0': np.array([0.4, 0., 0., 0.])}, + 'total_in': { + None: np.array([0.4, 0., 0., 0.]) + }, + 'total_out': { + None: np.array([0.4, 0., 0., 0.]) + }, + 'origins': { + None: { + 'compartment_0': { + None: np.array([0.4, 0., 0., 0.]) + } + } + }, + 'destinations': { + None: { + 'compartment_0': { + None: np.array([0.4, 0., 0., 0.]) + } + } + } } } flow_rates = self.builder_simple.flow_sheet.get_flow_rates().to_dict() diff --git a/tests/test_flow_sheet.py b/tests/test_flow_sheet.py index 74d009de..1e19563f 100644 --- a/tests/test_flow_sheet.py +++ b/tests/test_flow_sheet.py @@ -5,11 +5,46 @@ from CADETProcess import CADETProcessError from CADETProcess.processModel import ComponentSystem from CADETProcess.processModel import ( - Inlet, Cstr, LumpedRateModelWithoutPores, Outlet + Inlet, Cstr, MCT, LumpedRateModelWithoutPores, Outlet ) from CADETProcess.processModel import FlowSheet +def assert_almost_equal_dict( + dict_actual, dict_expected, decimal=7, verbose=True): + """Helper function to assert nested dicts are (almost) equal. + + Because of floating point calculations, it is necessary to use + `np.assert_almost_equal` to check the flow rates. However, this does not + work well with nested dicts which is why this helper function was written. + + Parameters + ---------- + dict_actual : dict + The object to check. + dict_expected : dict + The expected object. + decimal : int, optional + Desired precision, default is 7. + err_msg : str, optional + The error message to be printed in case of failure. + verbose : bool, optional + If True, the conflicting values are appended to the error message. + + """ + for key in dict_actual: + if isinstance(dict_actual[key], dict): + assert_almost_equal_dict(dict_actual[key], dict_expected[key]) + + else: + np.testing.assert_almost_equal( + dict_actual[key], dict_expected[key], + decimal=decimal, + err_msg=f'Dicts are not equal in key {key}.', + verbose=verbose + ) + + def setup_single_cstr_flow_sheet(component_system=None): """ Set up a simple `FlowSheet` with a single Continuous Stirred Tank Reactor (CSTR). @@ -172,24 +207,59 @@ def test_connections(self): outlet = self.ssr_flow_sheet['outlet'] expected_connections = { feed: { - 'origins': [], - 'destinations': [cstr], + 'origins': None, + 'destinations': { + None: { + cstr:[None], + }, + }, }, eluent: { - 'origins': [], - 'destinations': [column], + 'origins': None, + 'destinations': { + None: { + column:[None], + }, + }, }, + cstr: { - 'origins': [feed, column], - 'destinations': [column], + 'origins': { + None: { + feed:[None], + column:[None], + }, + }, + 'destinations': { + None: { + column:[None], + }, + }, }, + column: { - 'origins': [cstr, eluent], - 'destinations': [cstr, outlet], + 'origins': { + None: { + cstr:[None], + eluent:[None], + }, + }, + 'destinations': { + None: { + cstr:[None], + outlet:[None], + }, + }, }, + outlet: { - 'origins': [column], - 'destinations': [], + 'origins':{ + None: { + column:[None], + }, + }, + + 'destinations': None, }, } @@ -225,30 +295,154 @@ def test_name_decorator(self): flow_sheet.add_connection(eluent, 'column') flow_sheet.add_connection(column, cstr) flow_sheet.add_connection('column', outlet) - expected_connections = { feed: { - 'origins': [], - 'destinations': [cstr], + 'origins': None, + 'destinations': { + 0: { + cstr:[0], + }, + }, }, eluent: { - 'origins': [], - 'destinations': [column], + 'origins': None, + 'destinations': { + 0: { + column:[0], + }, + }, }, + cstr: { - 'origins': [feed, column], - 'destinations': [column], + 'origins': { + 0: { + feed:[0], + column:[0], + }, + }, + 'destinations': { + 0: { + column:[0], + }, + }, }, + column: { - 'origins': [cstr, eluent], - 'destinations': [cstr, outlet], + 'origins': { + 0: { + cstr:[0], + eluent:[0], + }, + }, + 'destinations': { + 0: { + cstr:[0], + outlet:[0], + }, + }, }, + outlet: { - 'origins': [column], - 'destinations': [], + 'origins':{ + 0: { + column:[0], + }, + }, + + 'destinations': None, }, } + self.assertDictEqual( + self.ssr_flow_sheet.connections, expected_connections + ) + + self.assertTrue(self.ssr_flow_sheet.connection_exists(feed, cstr)) + self.assertTrue(self.ssr_flow_sheet.connection_exists(eluent, column)) + self.assertTrue(self.ssr_flow_sheet.connection_exists(column, outlet)) + + self.assertFalse(self.ssr_flow_sheet.connection_exists(feed, eluent)) + + def test_name_decorator(self): + feed = Inlet(self.component_system, name='feed') + eluent = Inlet(self.component_system, name='eluent') + cstr = Cstr(self.component_system, name='cstr') + column = LumpedRateModelWithoutPores( + self.component_system, name='column' + ) + outlet = Outlet(self.component_system, name='outlet') + + flow_sheet = FlowSheet(self.component_system) + + flow_sheet.add_unit(feed) + flow_sheet.add_unit(eluent) + flow_sheet.add_unit(cstr) + flow_sheet.add_unit(column) + flow_sheet.add_unit(outlet) + + flow_sheet.add_connection('feed', 'cstr') + flow_sheet.add_connection(cstr, column) + flow_sheet.add_connection(eluent, 'column') + flow_sheet.add_connection(column, cstr) + flow_sheet.add_connection('column', outlet) + + expected_connections = { + feed: { + 'origins': None, + 'destinations': { + None: { + cstr:[None], + }, + }, + }, + eluent: { + 'origins': None, + 'destinations': { + None: { + column:[None], + }, + }, + }, + + cstr: { + 'origins': { + None: { + feed:[None], + column:[None], + }, + }, + 'destinations': { + None: { + column:[None], + }, + }, + }, + + column: { + 'origins': { + None: { + cstr:[None], + eluent:[None], + }, + }, + 'destinations': { + None: { + outlet:[None], + cstr:[None], + }, + }, + }, + + outlet: { + 'origins':{ + None: { + column:[None], + }, + }, + + 'destinations': None, + }, + } self.assertDictEqual(flow_sheet.connections, expected_connections) # Connection already exists @@ -272,45 +466,94 @@ def test_flow_rates(self): expected_flow_rates = { 'feed': { - 'total_out': (0, 0, 0, 0), + 'total_out':{ + None: (0, 0, 0, 0), + }, + 'destinations': { - 'cstr': (0, 0, 0, 0), + None: { + 'cstr': { + None: (0, 0, 0, 0), + }, + }, }, }, 'eluent': { - 'total_out': (0, 0, 0, 0), + 'total_out': { + None: (0, 0, 0, 0), + }, 'destinations': { - 'column': (0, 0, 0, 0), + None: { + 'column': { + None: (0, 0, 0, 0), + }, + }, }, }, 'cstr': { - 'total_in': (0.0, 0, 0, 0), - 'total_out': (1.0, 0, 0, 0), + 'total_in': { + None: (0, 0, 0, 0), + }, + 'total_out': { + None: (1.0, 0, 0, 0), + }, 'origins': { - 'feed': (0, 0, 0, 0), - 'column': (0, 0, 0, 0), + None: { + 'feed': { + None: (0, 0, 0, 0), + }, + 'column': { + None: (0, 0, 0, 0), + }, + }, }, 'destinations': { - 'column': (1.0, 0, 0, 0), + None: { + 'column': { + None: (1.0, 0, 0, 0), + }, + } }, }, 'column': { - 'total_in': (1.0, 0, 0, 0), - 'total_out': (1.0, 0, 0, 0), + 'total_in': { + None: (1.0, 0, 0, 0), + }, + 'total_out': { + None: (1.0, 0, 0, 0), + }, 'origins': { - 'cstr': (1.0, 0, 0, 0), - 'eluent': (0, 0, 0, 0), + None: { + 'cstr': { + None: (1.0, 0, 0, 0), + }, + 'eluent': { + None: (0, 0, 0, 0), + }, + }, }, 'destinations': { - 'cstr': (0, 0, 0, 0), - 'outlet': (1.0, 0, 0, 0), + None: { + 'cstr': { + None: (0, 0, 0, 0), + }, + 'outlet': { + None: (1.0, 0, 0, 0), + }, + }, }, }, 'outlet': { 'origins': { - 'column': (1.0, 0, 0, 0), + None: { + 'column': { + None: (1.0, 0, 0, 0), + }, + }, }, - 'total_in': (1.0, 0, 0, 0), + 'total_in': { + None: (1.0, 0, 0, 0), + }, }, } @@ -326,46 +569,94 @@ def test_flow_rates(self): expected_flow_rates = { 'feed': { - 'total_out': (1, 0, 0, 0), + 'total_out': { + None: (1, 0, 0, 0), + }, 'destinations': { - 'cstr': (1, 0, 0, 0), + None: { + 'cstr': { + None: (1, 0, 0, 0), + }, + } }, }, 'eluent': { - 'total_out': (1, 0, 0, 0), + 'total_out': { + None: (1, 0, 0, 0) + }, 'destinations': { - 'column': (1, 0, 0, 0), + None: { + 'column': { + None: (1, 0, 0, 0), + }, + }, }, }, 'cstr': { - 'total_in': (1.0, 0, 0, 0), - 'total_out': (0.0, 0, 0, 0), + 'total_in': { + None: (1.0, 0, 0, 0), + }, + 'total_out': { + None: (0.0, 0, 0, 0), + }, 'origins': { - 'feed': (1, 0, 0, 0), - 'column': (0, 0, 0, 0), + None: { + 'feed': { + None: (1, 0, 0, 0), + }, + 'column': { + None: (0, 0, 0, 0), + }, + }, }, 'destinations': { - 'column': (0.0, 0, 0, 0), + None: { + 'column': { + None: (0.0, 0, 0, 0), + }, + }, }, }, 'column': { - 'total_in': (1.0, 0, 0, 0), - 'total_out': (1.0, 0, 0, 0), + 'total_in': { + None: (1.0, 0, 0, 0), + }, + 'total_out': { + None: (1.0, 0, 0, 0), + }, 'origins': { - 'cstr': (0, 0, 0, 0), - 'eluent': (1, 0, 0, 0), + None: { + 'cstr': { + None: (0, 0, 0, 0), + }, + 'eluent': { + None: (1, 0, 0, 0), + }, + }, }, 'destinations': { - 'cstr': (0, 0, 0, 0), - 'outlet': (1.0, 0, 0, 0), + None: { + 'cstr': { + None: (0, 0, 0, 0), + }, + 'outlet': { + None: (1.0, 0, 0, 0), + }, + }, }, }, 'outlet': { 'origins': { - 'column': (1.0, 0, 0, 0), + None: { + 'column':{ + None: (1.0, 0, 0, 0), + }, + }, }, - 'total_in': (1.0, 0, 0, 0), + 'total_in': { + None: (1.0, 0, 0, 0), }, + }, } np.testing.assert_equal( self.ssr_flow_sheet.get_flow_rates(), expected_flow_rates @@ -379,45 +670,93 @@ def test_flow_rates(self): expected_flow_rates = { 'feed': { - 'total_out': (0, 0, 0, 0), + 'total_out': { + None: (0, 0, 0, 0), + }, 'destinations': { - 'cstr': (0, 0, 0, 0), + None: { + 'cstr': { + None: (0, 0, 0, 0), + }, + }, }, }, 'eluent': { - 'total_out': (1, 0, 0, 0), + 'total_out': { + None: (1, 0, 0, 0), + }, 'destinations': { - 'column': (1, 0, 0, 0), + None: { + 'column': { + None: (1, 0, 0, 0), + }, + }, }, }, 'cstr': { - 'total_in': (0.0, 0, 0, 0), - 'total_out': (0.0, 0, 0, 0), + 'total_in': { + None: (0, 0, 0, 0), + }, + 'total_out': { + None: (0, 0, 0, 0), + }, 'origins': { - 'feed': (0, 0, 0, 0), - 'column': (0, 0, 0, 0), + None: { + 'feed': { + None: (0, 0, 0, 0), + }, + 'column': { + None: (0, 0, 0, 0), + }, + }, }, 'destinations': { - 'column': (0.0, 0, 0, 0), + None: { + 'column': { + None: (0.0, 0, 0, 0), + }, + }, }, }, 'column': { - 'total_in': (1.0, 0, 0, 0), - 'total_out': (1.0, 0, 0, 0), + 'total_in': { + None: (1.0, 0, 0, 0), + }, + 'total_out': { + None: (1.0, 0, 0, 0), + }, 'origins': { - 'cstr': (0, 0, 0, 0), - 'eluent': (1, 0, 0, 0), + None: { + 'cstr': { + None: (0, 0, 0, 0), + }, + 'eluent': { + None: (1, 0, 0, 0), + }, + }, }, 'destinations': { - 'cstr': (0, 0, 0, 0), - 'outlet': (1.0, 0, 0, 0), + None: { + 'cstr': { + None: (0, 0, 0, 0), + }, + 'outlet': { + None: (1.0, 0, 0, 0), + }, + }, }, }, 'outlet': { 'origins': { - 'column': (1.0, 0, 0, 0), + None: { + 'column': { + None: (1.0, 0, 0, 0), + }, + }, }, - 'total_in': (1.0, 0, 0, 0), + 'total_in': { + None: (1.0, 0, 0, 0), + }, }, } @@ -433,46 +772,94 @@ def test_flow_rates(self): expected_flow_rates = { 'feed': { - 'total_out': (0, 0, 0, 0), + 'total_out': { + None: (0, 0, 0, 0), + }, 'destinations': { - 'cstr': (0, 0, 0, 0), + None: { + 'cstr': { + None: (0, 0, 0, 0), + }, + }, }, }, 'eluent': { - 'total_out': (1, 0, 0, 0), + 'total_out': { + None: (1, 0, 0, 0), + }, 'destinations': { - 'column': (1, 0, 0, 0), + None: { + 'column': { + None: (1, 0, 0, 0), + }, + }, }, }, 'cstr': { - 'total_in': (1.0, 0, 0, 0), - 'total_out': (0.0, 0, 0, 0), + 'total_in': { + None: (1.0, 0, 0, 0), + }, + 'total_out': { + None: (0.0, 0, 0, 0), + }, 'origins': { - 'feed': (0, 0, 0, 0), - 'column': (1, 0, 0, 0), + None: { + 'feed': { + None: (0, 0, 0, 0), + }, + 'column': { + None: (1, 0, 0, 0), + }, + }, }, 'destinations': { - 'column': (0.0, 0, 0, 0), + None: { + 'column': { + None: (0.0, 0, 0, 0), + }, + }, }, }, 'column': { - 'total_in': (1.0, 0, 0, 0), - 'total_out': (1.0, 0, 0, 0), + 'total_in': { + None: (1.0, 0, 0, 0), + }, + 'total_out': { + None: (1.0, 0, 0, 0), + }, 'origins': { - 'cstr': (0, 0, 0, 0), - 'eluent': (1, 0, 0, 0), + None: { + 'cstr': { + None: (0, 0, 0, 0), + }, + 'eluent': { + None: (1, 0, 0, 0), + }, + }, }, 'destinations': { - 'cstr': (1, 0, 0, 0), - 'outlet': (0.0, 0, 0, 0), + None: { + 'cstr': { + None: (1, 0, 0, 0), + }, + 'outlet': { + None: (0.0, 0, 0, 0), + }, + }, }, }, 'outlet': { 'origins': { - 'column': (0, 0, 0, 0), + None: { + 'column': { + None: (0, 0, 0, 0), + }, + }, }, - 'total_in': (0, 0, 0, 0), + 'total_in': { + None: (0.0, 0, 0, 0), }, + }, } np.testing.assert_equal( @@ -482,8 +869,12 @@ def test_flow_rates(self): # Single Cstr expected_flow_rates = { 'cstr': { - 'total_in': [0.0, 0.0, 0.0, 0.0], - 'total_out': [0.0, 0.0, 0.0, 0.0], + 'total_in': { + None:[0.0, 0.0, 0.0, 0.0], + }, + 'total_out': { + None:[0.0, 0.0, 0.0, 0.0] + }, 'origins': {}, 'destinations': {} } @@ -507,34 +898,40 @@ def test_check_connectivity(self): def test_output_state(self): column = self.ssr_flow_sheet.column + output_state_expected = {None: [1, 0]} + output_state = self.ssr_flow_sheet._output_states[column] + np.testing.assert_equal(output_state, output_state_expected) + output_state_expected = [1, 0] output_state = self.ssr_flow_sheet.output_states[column] np.testing.assert_equal(output_state, output_state_expected) self.ssr_flow_sheet.set_output_state(column, [0, 1]) - output_state_expected = [0, 1] - output_state = self.ssr_flow_sheet.output_states[column] + output_state_expected = {None: [0, 1]} + output_state = self.ssr_flow_sheet._output_states[column] np.testing.assert_equal(output_state, output_state_expected) self.ssr_flow_sheet.set_output_state(column, 0) - output_state_expected = [1, 0] - output_state = self.ssr_flow_sheet.output_states[column] + output_state_expected = {None: [1, 0]} + output_state = self.ssr_flow_sheet._output_states[column] np.testing.assert_equal(output_state, output_state_expected) self.ssr_flow_sheet.set_output_state(column, [0.5, 0.5]) - output_state_expected = [0.5, 0.5] - output_state = self.ssr_flow_sheet.output_states[column] + output_state_expected = {None: [0.5, 0.5]} + output_state = self.ssr_flow_sheet._output_states[column] np.testing.assert_equal(output_state, output_state_expected) self.ssr_flow_sheet.set_output_state( column, { 'cstr': 0.1, - 'outlet': 0.9, + 'outlet': { + None: 0.9, + } } ) - output_state_expected = [0.1, 0.9] - output_state = self.ssr_flow_sheet.output_states[column] + output_state_expected = {None: [0.1, 0.9]} + output_state = self.ssr_flow_sheet._output_states[column] np.testing.assert_equal(output_state, output_state_expected) with self.assertRaises(TypeError): @@ -548,7 +945,9 @@ def test_output_state(self): column, { 'column': 0.1, - 'outlet': 0.9, + 'outlet': { + 0: 0.9, + } } ) @@ -581,6 +980,9 @@ def test_add_connection_error(self): with self.assertRaises(CADETProcessError): self.ssr_flow_sheet.add_connection(inlet, column) + with self.assertRaisesRegex(CADETProcessError, "not a port of"): + self.ssr_flow_sheet.set_output_state(column, [0.5,0.5], 'channel_5') + class TestCstrFlowRate(unittest.TestCase): """ @@ -619,11 +1021,11 @@ def test_continuous_flow(self): flow_rates = self.flow_sheet.get_flow_rates() - cstr_in = flow_rates['cstr']['total_in'] + cstr_in = flow_rates['cstr']['total_in'][None] cstr_in_expected = [1., 0., 0., 0.] np.testing.assert_almost_equal(cstr_in, cstr_in_expected) - cstr_out = flow_rates['cstr']['total_out'] + cstr_out = flow_rates['cstr']['total_out'][None] cstr_out_expected = [1., 0., 0., 0.] np.testing.assert_almost_equal(cstr_out, cstr_out_expected) @@ -632,11 +1034,11 @@ def test_continuous_flow(self): flow_rates = self.flow_sheet.get_flow_rates() - cstr_in = flow_rates['cstr']['total_in'] + cstr_in = flow_rates['cstr']['total_in'][None] cstr_in_expected = [1., 1., 0., 0.] np.testing.assert_almost_equal(cstr_in, cstr_in_expected) - cstr_out = flow_rates['cstr']['total_out'] + cstr_out = flow_rates['cstr']['total_out'][None] cstr_out_expected = [1., 1., 0., 0.] np.testing.assert_almost_equal(cstr_out, cstr_out_expected) @@ -645,11 +1047,11 @@ def test_no_flow(self): flow_rates = self.flow_sheet.get_flow_rates() - cstr_in = flow_rates['cstr']['total_in'] + cstr_in = flow_rates['cstr']['total_in'][None] cstr_in_expected = [0., 0., 0., 0.] np.testing.assert_almost_equal(cstr_in, cstr_in_expected) - cstr_out = flow_rates['cstr']['total_out'] + cstr_out = flow_rates['cstr']['total_out'][None] cstr_out_expected = [0., 0., 0., 0.] np.testing.assert_almost_equal(cstr_out, cstr_out_expected) @@ -658,11 +1060,11 @@ def test_no_flow(self): flow_rates = self.flow_sheet.get_flow_rates() - cstr_in = flow_rates['cstr']['total_in'] + cstr_in = flow_rates['cstr']['total_in'][None] cstr_in_expected = [0., 0., 0., 0.] np.testing.assert_almost_equal(cstr_in, cstr_in_expected) - cstr_out = flow_rates['cstr']['total_out'] + cstr_out = flow_rates['cstr']['total_out'][None] cstr_out_expected = [0., 0., 0., 0.] np.testing.assert_almost_equal(cstr_out, cstr_out_expected) @@ -672,11 +1074,11 @@ def test_holdup(self): flow_rates = self.flow_sheet.get_flow_rates() - cstr_in = flow_rates['cstr']['total_in'] + cstr_in = flow_rates['cstr']['total_in'][None] cstr_in_expected = [1., 0., 0., 0.] np.testing.assert_almost_equal(cstr_in, cstr_in_expected) - cstr_out = flow_rates['cstr']['total_out'] + cstr_out = flow_rates['cstr']['total_out'][None] cstr_out_expected = [0., 0., 0., 0.] np.testing.assert_almost_equal(cstr_out, cstr_out_expected) @@ -689,14 +1091,472 @@ def test_state_update(self): flow_rates = self.flow_sheet.get_flow_rates(state) - cstr_in = flow_rates['cstr']['total_in'] + cstr_in = flow_rates['cstr']['total_in'][None] cstr_in_expected = [1., 1., 0., 0.] np.testing.assert_almost_equal(cstr_in, cstr_in_expected) - cstr_out = flow_rates['cstr']['total_out'] + cstr_out = flow_rates['cstr']['total_out'][None] cstr_out_expected = [2., 2., 0., 0.] np.testing.assert_almost_equal(cstr_out, cstr_out_expected) +class TestPorts(unittest.TestCase): + def __init__(self, methodName='runTest'): + super().__init__(methodName) + + def setUp(self): + self.setup_mct_flow_sheet() + + self.setup_ccc_flow_sheet() + + def setup_mct_flow_sheet(self): + self.component_system = ComponentSystem(1) + + mct_flow_sheet = FlowSheet(self.component_system) + + inlet = Inlet(self.component_system, name='inlet') + mct_3c = MCT(self.component_system,nchannel=3, name='mct_3c') + mct_2c1 = MCT(self.component_system,nchannel=2, name='mct_2c1') + mct_2c2 = MCT(self.component_system,nchannel=2, name='mct_2c2') + outlet1 = Outlet(self.component_system, name='outlet1') + outlet2 = Outlet(self.component_system, name='outlet2') + + mct_flow_sheet.add_unit(inlet) + mct_flow_sheet.add_unit(mct_3c) + mct_flow_sheet.add_unit(mct_2c1) + mct_flow_sheet.add_unit(mct_2c2) + mct_flow_sheet.add_unit(outlet1) + mct_flow_sheet.add_unit(outlet2) + + mct_flow_sheet.add_connection(inlet, mct_3c, destination_port='channel_0') + mct_flow_sheet.add_connection(mct_3c, mct_2c1, origin_port='channel_0', destination_port='channel_0') + mct_flow_sheet.add_connection(mct_3c, mct_2c1, origin_port='channel_0', destination_port='channel_1') + mct_flow_sheet.add_connection(mct_3c, mct_2c2, origin_port='channel_1', destination_port='channel_0') + mct_flow_sheet.add_connection(mct_2c1, outlet1, origin_port='channel_0') + mct_flow_sheet.add_connection(mct_2c1, outlet1, origin_port='channel_1') + mct_flow_sheet.add_connection(mct_2c2, outlet2, origin_port='channel_0') + + + self.mct_flow_sheet = mct_flow_sheet + + def setup_ccc_flow_sheet(self): + self.component_system = ComponentSystem(1) + + ccc_flow_sheet = FlowSheet(self.component_system) + + inlet = Inlet(self.component_system, name='inlet') + ccc1 = MCT(self.component_system,nchannel=2, name='ccc1') + ccc2 = MCT(self.component_system,nchannel=2, name='ccc2') + ccc3 = MCT(self.component_system,nchannel=2, name='ccc3') + outlet = Outlet(self.component_system, name='outlet') + + ccc_flow_sheet.add_unit(inlet) + ccc_flow_sheet.add_unit(ccc1) + ccc_flow_sheet.add_unit(ccc2) + ccc_flow_sheet.add_unit(ccc3) + ccc_flow_sheet.add_unit(outlet) + + ccc_flow_sheet.add_connection(inlet, ccc1, destination_port='channel_0') + ccc_flow_sheet.add_connection(ccc1, ccc2, origin_port='channel_0', destination_port='channel_0') + ccc_flow_sheet.add_connection(ccc2, ccc3, origin_port='channel_0', destination_port='channel_0') + ccc_flow_sheet.add_connection(ccc3, outlet, origin_port='channel_0') + + + self.ccc_flow_sheet = ccc_flow_sheet + + def test_mct_connections(self): + inlet = self.mct_flow_sheet['inlet'] + mct_3c = self.mct_flow_sheet['mct_3c'] + mct_2c1 = self.mct_flow_sheet['mct_2c1'] + mct_2c2 = self.mct_flow_sheet['mct_2c2'] + outlet1 = self.mct_flow_sheet['outlet1'] + outlet2 = self.mct_flow_sheet['outlet2'] + + expected_connections = { + inlet: { + 'origins': None, + 'destinations': { + None: { + mct_3c: ['channel_0'], + }, + }, + }, + mct_3c: { + 'origins': { + 'channel_0': { + inlet: [None], + }, + 'channel_1': { + + }, + 'channel_2': { + + }, + }, + 'destinations': { + 'channel_0': { + mct_2c1: ['channel_0', 'channel_1'], + }, + 'channel_1': { + mct_2c2: ['channel_0'], + }, + 'channel_2': { + + }, + }, + }, + + mct_2c1: { + 'origins': { + 'channel_0': { + mct_3c:['channel_0'], + }, + + 'channel_1': { + mct_3c:['channel_0'], + }, + }, + 'destinations': { + 'channel_0': { + outlet1: [None], + }, + 'channel_1': { + outlet1: [None], + }, + }, + }, + + mct_2c2: { + 'origins': { + 'channel_0': { + mct_3c:['channel_1'], + }, + 'channel_1': { + + }, + }, + + 'destinations': { + 'channel_0': { + outlet2: [None], + }, + 'channel_1': { + + }, + }, + }, + + outlet1: { + 'origins':{ + None: { + mct_2c1:['channel_0','channel_1'], + }, + }, + + 'destinations': None, + }, + + outlet2: { + 'origins':{ + None: { + mct_2c2:['channel_0'], + }, + }, + + 'destinations': None, + }, + } + + self.assertDictEqual( + self.mct_flow_sheet.connections, expected_connections + ) + + self.assertTrue(self.mct_flow_sheet.connection_exists(inlet, mct_3c, destination_port='channel_0')) + self.assertTrue(self.mct_flow_sheet.connection_exists(mct_3c, mct_2c1, 'channel_0', 'channel_0')) + self.assertTrue(self.mct_flow_sheet.connection_exists(mct_3c, mct_2c1, 'channel_0', 'channel_1')) + self.assertTrue(self.mct_flow_sheet.connection_exists(mct_3c, mct_2c2, 'channel_1', 'channel_0')) + self.assertTrue(self.mct_flow_sheet.connection_exists(mct_2c1, outlet1, 'channel_0')) + self.assertTrue(self.mct_flow_sheet.connection_exists(mct_2c1, outlet1, 'channel_1')) + self.assertTrue(self.mct_flow_sheet.connection_exists(mct_2c2, outlet2, 'channel_0')) + + self.assertFalse(self.mct_flow_sheet.connection_exists(mct_2c2, outlet2, 'channel_1')) + self.assertFalse(self.mct_flow_sheet.connection_exists(inlet, mct_2c1, destination_port='channel_0')) + + def test_ccc_connections(self): + inlet = self.ccc_flow_sheet['inlet'] + ccc1 = self.ccc_flow_sheet['ccc1'] + ccc2 = self.ccc_flow_sheet['ccc2'] + ccc3 = self.ccc_flow_sheet['ccc3'] + outlet = self.ccc_flow_sheet['outlet'] + + expected_connections = { + inlet: { + 'origins': None, + 'destinations': { + None: { + ccc1: ['channel_0'], + }, + }, + }, + ccc1: { + 'origins': { + 'channel_0': { + inlet: [None], + }, + 'channel_1': { + + } + }, + 'destinations': { + 'channel_0': { + ccc2: ['channel_0'], + }, + 'channel_1': { + + } + }, + }, + + ccc2: { + 'origins': { + 'channel_0': { + ccc1: ['channel_0'], + }, + 'channel_1': { + + } + }, + 'destinations': { + 'channel_0': { + ccc3: ['channel_0'], + }, + 'channel_1': { + + } + }, + }, + + ccc3: { + 'origins': { + 'channel_0': { + ccc2: ['channel_0'], + }, + 'channel_1': { + + } + }, + 'destinations': { + 'channel_0': { + outlet: [None], + }, + 'channel_1': { + + } + }, + }, + + outlet: { + 'origins':{ + None: { + ccc3:['channel_0'], + }, + }, + + 'destinations': None, + }, + } + + self.assertDictEqual( + self.ccc_flow_sheet.connections, expected_connections + ) + + self.assertTrue(self.ccc_flow_sheet.connection_exists(inlet, ccc1, destination_port='channel_0')) + self.assertTrue(self.ccc_flow_sheet.connection_exists(ccc1, ccc2, 'channel_0', 'channel_0')) + self.assertTrue(self.ccc_flow_sheet.connection_exists(ccc2, ccc3, 'channel_0', 'channel_0')) + self.assertTrue(self.ccc_flow_sheet.connection_exists(ccc3, outlet, 'channel_0')) + + self.assertFalse(self.ccc_flow_sheet.connection_exists(inlet, outlet)) + + def test_mct_flow_rate_calculation(self): + + mct_flow_sheet = self.mct_flow_sheet + + mct_flow_sheet.inlet.flow_rate = 1 + + inlet = self.mct_flow_sheet['inlet'] + mct_3c = self.mct_flow_sheet['mct_3c'] + mct_2c1 = self.mct_flow_sheet['mct_2c1'] + mct_2c2 = self.mct_flow_sheet['mct_2c2'] + +# mct_flow_sheet.set_output_state(mct_3c, [0.5, 0.5], 0) + + expected_flow = { + 'inlet': { + 'total_out': { + None: (1,0,0,0) + }, + 'destinations': { + None: { + 'mct_3c': { + 'channel_0': (1,0,0,0) + } + } + } + }, + 'mct_3c': { + 'total_in': { + 'channel_0': (1,0,0,0), + 'channel_1': (0,0,0,0), + 'channel_2': (0,0,0,0) + }, + 'origins': { + 'channel_0': { + 'inlet':{ + None: (1,0,0,0) + } + } + }, + 'total_out': { + 'channel_0': (1,0,0,0), + 'channel_1': (0,0,0,0), + 'channel_2': (0,0,0,0) + }, + 'destinations': { + 'channel_0': { + 'mct_2c1': { + 'channel_0': (1,0,0,0), + 'channel_1': (0,0,0,0) + } + }, + 'channel_1': { + 'mct_2c2': { + 'channel_0': (0,0,0,0) + } + } + } + }, + 'mct_2c1': { + 'total_in': { + 'channel_0': (1,0,0,0), + 'channel_1': (0,0,0,0) + }, + 'origins': { + 'channel_0': { + 'mct_3c': { + 'channel_0': (1,0,0,0) + } + }, + 'channel_1': { + 'mct_3c': { + 'channel_0': (0,0,0,0) + } + } + }, + 'total_out': { + 'channel_0': (1,0,0,0), + 'channel_1': (0,0,0,0) + }, + 'destinations': { + 'channel_0': { + 'outlet1': { + None: (1,0,0,0) + } + }, + 'channel_1': { + 'outlet1': { + None: (0,0,0,0) + } + } + } + }, + + 'mct_2c2': { + 'total_in': { + 'channel_0': (0,0,0,0), + 'channel_1': (0,0,0,0) + }, + 'origins': { + 'channel_0': { + 'mct_3c': { + 'channel_1': (0,0,0,0) + } + }, + }, + 'total_out': { + 'channel_0': (0,0,0,0), + 'channel_1': (0,0,0,0) + }, + 'destinations': { + 'channel_0': { + 'outlet2': { + None: (0,0,0,0) + } + } + } + }, + 'outlet1': { + 'total_in': { + None: (1,0,0,0) + }, + 'origins': { + None: { + 'mct_2c1': { + 'channel_0': (1,0,0,0), + 'channel_1': (0,0,0,0) + } + } + } + }, + 'outlet2': { + 'total_in': { + None: (0,0,0,0) + }, + 'origins': { + None: { + 'mct_2c2': { + 'channel_0': (0,0,0,0) + } + } + } + } + } + + + flow_rates = mct_flow_sheet.get_flow_rates() + + + assert_almost_equal_dict(flow_rates, expected_flow) + + + + def test_port_add_connection(self): + + inlet = self.mct_flow_sheet['inlet'] + mct_3c = self.mct_flow_sheet['mct_3c'] + mct_2c2 = self.mct_flow_sheet['mct_2c2'] + outlet1 = self.mct_flow_sheet['outlet1'] + + with self.assertRaises(CADETProcessError): + self.mct_flow_sheet.add_connection(inlet, mct_3c, origin_port=0, destination_port=5) + + with self.assertRaises(CADETProcessError): + self.mct_flow_sheet.add_connection(inlet, mct_3c, origin_port=5, destination_port=0) + + with self.assertRaises(CADETProcessError): + self.mct_flow_sheet.add_connection(mct_2c2, outlet1, origin_port=0, destination_port=5) + + def test_set_output_state(self): + + mct_3c = self.mct_flow_sheet['mct_3c'] + + with self.assertRaises(CADETProcessError): + self.mct_flow_sheet.set_output_state(mct_3c, [0.5,0.5]) + + with self.assertRaises(CADETProcessError): + self.mct_flow_sheet.set_output_state(mct_3c, [0.5,0.5], 'channel_5') + + with self.assertRaises(CADETProcessError): + self.mct_flow_sheet.set_output_state(mct_3c, {'mct_2c1': {'channel_0': 0.5 , 'channel_5': 0.5}}, 'channel_0') + class TestFlowRateMatrix(unittest.TestCase): """Test calculation of flow rates with another simple testcase by @daklauss""" @@ -740,84 +1600,99 @@ def test_matrix_example(self): expected_flow_rates = { 'inlet': { - 'total_out': (1, 0, 0, 0), + 'total_out': { + None: (1, 0, 0, 0), + }, 'destinations': { - 'cstr1': (0.3, 0, 0, 0), - 'cstr2': (0.7, 0, 0, 0), + None: { + 'cstr1': { + None: (0.3, 0, 0, 0), + }, + 'cstr2': { + None: (0.7, 0, 0, 0), + }, + }, }, }, 'cstr1': { - 'total_in': (0.65, 0, 0, 0), - 'total_out': (0.65, 0, 0, 0), + 'total_in': { + None: (0.65, 0, 0, 0), + }, + 'total_out': { + None: (0.65, 0, 0, 0), + }, 'origins': { - 'inlet': (0.3, 0, 0, 0), - 'cstr2': (0.35, 0, 0, 0), + None: { + 'inlet': { + None: (0.3, 0, 0, 0), + }, + 'cstr2': { + None: (0.35, 0, 0, 0), + }, + } }, 'destinations': { - 'outlet1': (0.65, 0, 0, 0), + None: { + 'outlet1': { + None: (0.65, 0, 0, 0), + }, + }, }, }, 'cstr2': { - 'total_in': (0.7, 0, 0, 0), - 'total_out': (0.7, 0, 0, 0), + 'total_in': { + None: (0.7, 0, 0, 0), + }, + 'total_out': { + None: (0.7, 0, 0, 0), + }, 'origins': { - 'inlet': (0.7, 0, 0, 0), + None: { + 'inlet': { + None: (0.7, 0, 0, 0), + }, + } }, 'destinations': { - 'cstr1': (0.35, 0, 0, 0), - 'outlet2': (0.35, 0, 0, 0), + None: { + 'cstr1': { + None: (0.35, 0, 0, 0), + }, + 'outlet2': { + None: (0.35, 0, 0, 0), + }, + }, }, }, 'outlet1': { 'origins': { - 'cstr1': (0.65, 0, 0, 0), + None: { + 'cstr1': { + None: (0.65, 0, 0, 0), + }, + }, + }, + 'total_in': { + None: (0.65, 0, 0, 0), }, - 'total_in': (0.65, 0, 0, 0), }, 'outlet2': { 'origins': { - 'cstr2': (0.35, 0, 0, 0), + None: { + 'cstr2': { + None: (0.35, 0, 0, 0), + }, + } }, - 'total_in': (0.35, 0, 0, 0), - } + 'total_in': { + None: (0.35, 0, 0, 0), + }, + } } calc_flow_rate = self.flow_sheet.get_flow_rates() - def assert_almost_equal_dict( - dict_actual, dict_expected, decimal=7, verbose=True): - """Helper function to assert nested dicts are (almost) equal. - - Because of floating point calculations, it is necessary to use - `np.assert_almost_equal` to check the flow rates. However, this does not - work well with nested dicts which is why this helper function was written. - - Parameters - ---------- - dict_actual : dict - The object to check. - dict_expected : dict - The expected object. - decimal : int, optional - Desired precision, default is 7. - err_msg : str, optional - The error message to be printed in case of failure. - verbose : bool, optional - If True, the conflicting values are appended to the error message. - - """ - for key in dict_actual: - if isinstance(dict_actual[key], dict): - assert_almost_equal_dict(dict_actual[key], dict_expected[key]) - else: - np.testing.assert_almost_equal( - dict_actual[key], dict_expected[key], - decimal=decimal, - err_msg=f'Dicts are not equal in key {key}.', - verbose=verbose - ) - assert_almost_equal_dict(calc_flow_rate, expected_flow_rates) @@ -853,29 +1728,57 @@ def test_matrix_self_example(self): expected_flow_rates = { 'inlet': { - 'total_out': (1, 0, 0, 0), + 'total_out': { + None: (1, 0, 0, 0), + }, 'destinations': { - 'cstr': (1, 0, 0, 0) + None: { + 'cstr': { + None: (1, 0, 0, 0), + }, + }, }, }, 'cstr': { - 'total_in': (2, 0, 0, 0), - 'total_out': (2, 0, 0, 0), + 'total_in': { + None: (2, 0, 0, 0), + }, + 'total_out': { + None: (2, 0, 0, 0), + }, 'origins': { - 'inlet': (1, 0, 0, 0), - 'cstr': (1, 0, 0, 0), + None: { + 'inlet': { + None: (1, 0, 0, 0), + }, + 'cstr': { + None: (1, 0, 0, 0), + }, + }, }, 'destinations': { - 'outlet': (1, 0, 0, 0), - 'cstr': (1, 0, 0, 0) + None: { + 'outlet': { + None: (1, 0, 0, 0), + }, + 'cstr': { + None: (1, 0, 0, 0), + }, + }, }, }, 'outlet': { - 'total_in': (1, 0, 0, 0), + 'total_in': { + None: (1, 0, 0, 0), + }, 'origins': { - 'cstr': (1, 0, 0, 0) - } - } + None: { + 'cstr': { + None: (1, 0, 0, 0), + }, + }, + }, + }, } calc_flow_rate = self.flow_sheet.get_flow_rates() np.testing.assert_equal(calc_flow_rate, expected_flow_rates) @@ -933,37 +1836,73 @@ def test_expelled_circuit_with_flow(self): # Solvable because both disconnected circles have their own flow rates expected_flow_rates = { 'inlet': { - 'total_out': (1, 0, 0, 0), + 'total_out': { + None: (1, 0, 0, 0), + }, 'destinations': { - 'outlet': (1, 0, 0, 0) + None: { + 'outlet': { + None: (1, 0, 0, 0), + }, + }, }, }, 'cstr1': { - 'total_in': (1, 0, 0, 0), - 'total_out': (1, 0, 0, 0), + 'total_in': { + None: (1, 0, 0, 0), + }, + 'total_out': { + None: (1, 0, 0, 0), + }, 'origins': { - 'cstr2': (1, 0, 0, 0), + None: { + 'cstr2': { + None: (1, 0, 0, 0), + }, + }, }, 'destinations': { - 'cstr2': (1, 0, 0, 0) + None: { + 'cstr2': { + None: (1, 0, 0, 0) + }, + }, }, }, 'cstr2': { - 'total_in': (1, 0, 0, 0), - 'total_out': (1, 0, 0, 0), + 'total_in': { + None: (1, 0, 0, 0), + }, + 'total_out': { + None: (1, 0, 0, 0), + }, 'origins': { - 'cstr1': (1, 0, 0, 0), + None: { + 'cstr1': { + None: (1, 0, 0, 0), + }, + }, }, 'destinations': { - 'cstr1': (1, 0, 0, 0) + None: { + 'cstr1': { + None: (1, 0, 0, 0), + }, + }, }, }, 'outlet': { - 'total_in': (1, 0, 0, 0), + 'total_in': { + None: (1, 0, 0, 0), + }, 'origins': { - 'inlet': (1, 0, 0, 0) - } - } + None: { + 'inlet': { + None: (1, 0, 0, 0), + }, + }, + }, + }, } flow_sheet = self.flow_sheet diff --git a/tests/test_optimization_problem.py b/tests/test_optimization_problem.py index c1c9d5d1..ee471728 100644 --- a/tests/test_optimization_problem.py +++ b/tests/test_optimization_problem.py @@ -11,7 +11,7 @@ Structure, Float, List, SizedList, SizedNdArray, Polynomial, NdPolynomial ) from CADETProcess.optimization import OptimizationProblem -from optimization_problem_fixtures import ( +from tests.optimization_problem_fixtures import ( LinearConstraintsSooTestProblem2, LinearEqualityConstraintsSooTestProblem ) @@ -330,7 +330,7 @@ def test_transform(self): ) -from test_events import TestHandler +from tests.test_events import TestHandler class Test_OptimizationVariableEvents(unittest.TestCase): def __init__(self, methodName='runTest'): diff --git a/tests/test_optimization_results.py b/tests/test_optimization_results.py index 9f27a19f..48ba9418 100644 --- a/tests/test_optimization_results.py +++ b/tests/test_optimization_results.py @@ -5,8 +5,8 @@ from CADETProcess.optimization import OptimizationResults from CADETProcess.optimization import U_NSGA3 -from test_population import setup_population -from test_optimization_problem import setup_optimization_problem +from tests.test_population import setup_population +from tests.test_optimization_problem import setup_optimization_problem def setup_optimization_results( diff --git a/tests/test_parallelization_adapter.py b/tests/test_parallelization_adapter.py index 9d441451..ba8c91b6 100644 --- a/tests/test_parallelization_adapter.py +++ b/tests/test_parallelization_adapter.py @@ -10,8 +10,8 @@ from CADETProcess.simulator import Cadet from CADETProcess.optimization import SequentialBackend -from test_cadet_adapter import detect_cadet -from test_optimization_problem import setup_optimization_problem +from tests.test_cadet_adapter import detect_cadet +from tests.test_optimization_problem import setup_optimization_problem parallel_backends_module = importlib.import_module( diff --git a/tests/test_population.py b/tests/test_population.py index 277f252a..45fa7e66 100644 --- a/tests/test_population.py +++ b/tests/test_population.py @@ -4,7 +4,7 @@ from CADETProcess import CADETProcessError from CADETProcess.optimization import Individual, Population, ParetoFront -from test_individual import setup_individual +from tests.test_individual import setup_individual enable_plot = False diff --git a/tests/test_pymoo.py b/tests/test_pymoo.py index be0a5210..a27a08cf 100644 --- a/tests/test_pymoo.py +++ b/tests/test_pymoo.py @@ -2,7 +2,7 @@ from CADETProcess.optimization import U_NSGA3 -from test_optimization_problem import setup_optimization_problem +from tests.test_optimization_problem import setup_optimization_problem class Test_OptimizationProblemSimple(unittest.TestCase): diff --git a/tests/test_unit_operation.py b/tests/test_unit_operation.py index 8dbfc986..7a04b1fe 100644 --- a/tests/test_unit_operation.py +++ b/tests/test_unit_operation.py @@ -8,10 +8,11 @@ import numpy as np +from CADETProcess import CADETProcessError from CADETProcess.processModel import ComponentSystem from CADETProcess.processModel import ( Inlet, Cstr, - TubularReactor, LumpedRateModelWithPores, LumpedRateModelWithoutPores + TubularReactor, LumpedRateModelWithPores, LumpedRateModelWithoutPores, MCT ) length = 0.6 @@ -27,6 +28,14 @@ axial_dispersion = 4.7e-7 +channel_cross_section_areas = [0.1,0.1,0.1] +exchange_matrix = np.array([ + [[0.0],[0.01],[0.0]], + [[0.02],[0.0],[0.03]], + [[0.0],[0.0],[0.0]] + ]) +flow_direction = 1 + class Test_Unit_Operation(unittest.TestCase): @@ -60,6 +69,16 @@ def create_tubular_reactor(self): return tube + def create_MCT(self, components): + mct = MCT(ComponentSystem(components), nchannel=3, name='test') + + mct.length = length + mct.channel_cross_section_areas = channel_cross_section_areas + mct.axial_dispersion = 0 + mct.flow_direction = flow_direction + + return mct + def create_lrmwop(self): lrmwop = LumpedRateModelWithoutPores( self.component_system, name='test' @@ -204,6 +223,7 @@ def test_parameters(self): 'q': [], 'V': volume, } + np.testing.assert_equal(parameters_expected, cstr.parameters) sec_dep_parameters_expected = { @@ -224,5 +244,65 @@ def test_parameters(self): self.assertEqual(cstr.required_parameters, ['V']) + def test_MCT(self): + """ + Notes + ----- + Tests Parameters, Volumes and Attributes depending on nchannel. Should be later integrated into general testing workflow. + """ + total_porosity = 1 + + mct = self.create_MCT(1) + + mct.exchange_matrix = exchange_matrix + + parameters_expected = { + 'c': np.array([[0., 0., 0.]]), + 'axial_dispersion' : 0, + 'channel_cross_section_areas' : channel_cross_section_areas, + 'length' : length, + 'exchange_matrix': exchange_matrix, + 'flow_direction' : 1, + 'nchannel' : 3 + } + np.testing.assert_equal(parameters_expected, {key: value for key, value in mct.parameters.items() if key != 'discretization'}) + + volume = length*sum(channel_cross_section_areas) + volume_liquid = volume*total_porosity + volume_solid = (total_porosity-1)*volume + + self.assertAlmostEqual(mct.volume_liquid, volume_liquid) + self.assertAlmostEqual(mct.volume_solid, volume_solid) + + with self.assertRaises(ValueError): + mct.exchange_matrix = np.array([[ + [0.0, 0.01, 0.0], + [0.02, 0.0, 0.03], + [0.0, 0.0, 0.0] + ]]) + + mct.nchannel = 2 + with self.assertRaises(ValueError): + mct.exchange_matrix + mct.channel_cross_section_areas + + self.assertTrue(mct.nchannel*mct.component_system.n_comp == mct.c.size) + + mct2 = self.create_MCT(2) + + with self.assertRaises(ValueError): + mct2.exchange_matrix = np.array([[ + [0.0, 0.01, 0.0], + [0.02, 0.0, 0.03], + [0.0, 0.0, 0.0] + ], + + [ + [0.0, 0.01, 0.0], + [0.02, 0.0, 0.03], + [0.0, 0.0, 0.0] + ]]) + + if __name__ == '__main__': unittest.main()