Expected behavior.
Hi everyone,
I am trying to use HICConstantWaterActivity isotherm. I found an unexpect phenomenon, when I added a ‘CSTR’ the simulation never stops, but when I annotated “CSTR” code, the simulation succeeds very quickly. Would anyone know how to solve this problem? Thank you for your help.
Regards
Yang
Actual behavior
The simulation never stops
How to produce bug (including a minimal reproducible example)
from CADETProcess.processModel import ComponentSystem
from CADETProcess.processModel import Inlet
from CADETProcess.processModel import Cstr
from CADETProcess.processModel import HICConstantWaterActivity
from CADETProcess.processModel import GeneralRateModel
from CADETProcess.processModel import Outlet
from CADETProcess.processModel import FlowSheet
from CADETProcess.processModel import Process
from CADETProcess.simulator import Cadet
from matplotlib.lines import Line2D
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
def Simulation():
# Component System
component_system = ComponentSystem()
component_system.add_component('Salt')
component_system.add_component('A')
component_system.add_component('B')
# Binding Model
binding_model = HICConstantWaterActivity(component_system)
binding_model.is_kinetic = True
binding_model.adsorption_rate = [0.0, 0.7, 1]
binding_model.desorption_rate = [0.0, 1000, 1000]
binding_model.hic_characteristic = [0.0, 10, 13]
binding_model.capacity = [0.0, 10000000, 10000000]
binding_model.beta_0 = 10. ** -0.5
binding_model.beta_1 = 10. ** -3.65
binding_model.bound_states = [0, 1, 1]
Q = 1.94e-8
salt_gradient_start_concentration = 3000
salt_gradient_end_concentration = 50
# Unit Operations
inlet = Inlet(component_system, name='inlet')
inlet.flow_rate = Q
inlet.c = [salt_gradient_start_concentration, 0, 0]
column = GeneralRateModel(component_system, name='column')
column.binding_model = binding_model
column.length = 0.014
column.diameter = 0.02
column.bed_porosity = 0.37
column.particle_radius = 4.5e-5
column.particle_porosity = 0.75
column.axial_dispersion = 5.75e-8
column.film_diffusion = column.n_comp * [6.9e-6]
column.pore_diffusion = [7e-10, 6.07e-11, 6.07e-11]
column.surface_diffusion = column.n_bound_states * [0.0]
column.c = [salt_gradient_start_concentration, 0, 0]
column.cp = [salt_gradient_start_concentration, 0, 0]
column.q = [0, 0]
# CSTR1
cstr1 = Cstr(component_system, 'cstr1')
cstr1.V = 1840.556e-9 # m^3
cstr1.c = [salt_gradient_start_concentration, 0, 0]
outlet = Outlet(component_system, name='outlet')
# Flow Sheet
flow_sheet = FlowSheet(component_system)
flow_sheet.add_unit(inlet)
flow_sheet.add_unit(column)
flow_sheet.add_unit(cstr1)
flow_sheet.add_unit(outlet, product_outlet=True)
flow_sheet.add_connection(inlet, cstr1)
flow_sheet.add_connection(cstr1, column)
flow_sheet.add_connection(column, outlet)
# flow_sheet.add_connection(inlet, column)
# flow_sheet.add_connection(column, outlet)
# Process
process = Process(flow_sheet, 'lwe')
process.cycle_time = 2000.0
load_duration = 9
t_gradient_start = 90.0
gradient_duration = process.cycle_time - t_gradient_start
c_load = np.array([salt_gradient_start_concentration, 1, 1])
c_wash = np.array([salt_gradient_start_concentration, 0.0, 0.0])
c_elute = np.array([salt_gradient_end_concentration, 0.0, 0.0])
gradient_slope = (c_elute - c_wash) / gradient_duration
c_gradient_poly = np.array(list(zip(c_wash, gradient_slope)))
process.add_event('load', 'flow_sheet.inlet.c', c_load)
process.add_event('wash', 'flow_sheet.inlet.c', c_wash, load_duration)
process.add_event('grad_start', 'flow_sheet.inlet.c', c_gradient_poly, t_gradient_start)
# simulate
simulator = Cadet()
results = simulator.simulate(process)
t_SIM = results.solution.outlet.outlet.time
c_A_SIM = results.solution.outlet.outlet.solution[:, 1]
c_B_SIM = results.solution.outlet.outlet.solution[:, 2]
t_salt_SIM = results.solution.outlet.outlet.time
c_salt_SIM = results.solution.outlet.outlet.solution[:, 0]
fig, ax1 = plt.subplots(figsize=(10, 6))
coloura = np.array([247, 172, 83]) / 255
colourcyto = np.array([236, 110, 102]) / 255
coloursalt = np.array([120, 149, 193]) / 255
ax1.set_xlabel('time s')
ax1.set_ylabel('protein mg/mL', color='tab:blue')
ax1.plot(t_SIM, c_A_SIM, '-', color=coloura)
ax1.plot(t_SIM, c_B_SIM, '-', color=colourcyto)
ax2 = ax1.twinx()
ax2.set_ylabel('salt mM', color='tab:green')
ax2.plot(t_salt_SIM, c_salt_SIM, '-', color=coloursalt)
ax1.set_xlim([0, max(t_salt_SIM)])
legend_elements_1 = [
Line2D([0], [0], color=coloura, marker='o', linestyle='None', label='a'),
Line2D([0], [0], color=colourcyto, marker='o', linestyle='None', label='cyto'),
Line2D([0], [0], color=coloursalt, marker='o', linestyle='None', label='salt'),
]
legend_elements_2 = [
Line2D([0], [0], color=coloura, linestyle='-', label='SIM'),
]
ax1.legend(handles=legend_elements_1, loc='upper left', fontsize=16)
ax2.legend(handles=legend_elements_2, loc='upper right', fontsize=16)
plt.tight_layout()
plt.show()
plt.savefig('simulation_results.png')
Simulation()