Hi all,
We’ve observed a difference in the behavior of our model with the colloidal isotherm when switching from cadet v4.4 to 5.0. To try to make sure it wasn’t coming from anything else, I set up the same model with the multicomponent Langmuir model and some dummy parameters, but the difference in the output chromatograms only showed up with the colloidal isotherm.
The script I used to generate the figures is copied below. The only thing I’ve changed is whether the Langmuir or colloidal isotherm is used and whether the cadet_bin_path references the environment where I have cadet 4.4.0 or 5.0.0 installed.
I realize the colloidal isotherm wasn’t official before, but does this mean that the previous version was in some way incorrect?
Code below:
# -*- coding: utf-8 -*-
"""
Inlet > column > outlet with GRM and colloidal isotherm.
"""
import platform
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
from cadet import Cadet
def cadet_config(cadet_bin_path, binding):
print('\nCADET simulation running!\n')
set_CADET_path(cadet_bin_path)
# Create model object
model = Cadet()
# Number of unit operations
model.root.input.model.nunits = 3
n_comp = 2
init_c = [1e3*10**(-7.4), 0]
# Inlet ####################################################################
model.root.input.model.unit_000.unit_type = 'INLET'
model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
model.root.input.model.unit_000.ncomp = n_comp
model.root.input.model.unit_000.sec_000.const_coeff = [1e3*10**(-7.4), 0.0338]
model.root.input.model.unit_000.sec_001.const_coeff = [1e3*10**(-7.4), 0]
model.root.input.model.unit_000.sec_002.const_coeff = [1e3*10**(-3.5), 0]
# General Rate Model #######################################################
model.root.input.model.unit_001.unit_type = 'GENERAL_RATE_MODEL'
model.root.input.model.unit_001.ncomp = n_comp
model.root.input.model.unit_001.init_c = init_c
model.root.input.model.unit_001.init_cp = init_c
model.root.input.model.unit_001.init_q = [0, 0]
## Geometry
model.root.input.model.unit_001.col_length = 0.018
model.root.input.model.unit_001.cross_section_area = 1.963e-05
model.root.input.model.unit_001.col_porosity = 0.4
model.root.input.model.unit_001.par_porosity = 0.92
model.root.input.model.unit_001.par_radius = 0.5*5.41e-05
## Transport
model.root.input.model.unit_001.col_dispersion = [1.80e-08, 1.80e-08, 1.80e-08, 1.80e-08]
model.root.input.model.unit_001.film_diffusion = [[0.00011, 0.00011, 0.00011, 0.00011], [1.32e-05, 1.32e-05, 1.32e-05, 1.32e-05]]
model.root.input.model.unit_001.par_diffusion = [8.85e-10, 1.54e-12]
model.root.input.model.unit_001.adsorption.is_kinetic = True
model.root.input.model.unit_001.discretization.ncol = 100
model.root.input.model.unit_001.discretization.npar = 40
## Adsorption
if binding=='Langmuir':
model.root.input.model.unit_001.adsorption_model = 'MULTI_COMPONENT_LANGMUIR'
model.root.input.model.unit_001.adsorption.is_kinetic = True # Kinetic binding
model.root.input.model.unit_001.adsorption.mcl_ka = [0, 0.2] # m^3 / (mol * s) (mobile phase)
model.root.input.model.unit_001.adsorption.mcl_kd = [0, 1.0] # 1 / s (desorption)
model.root.input.model.unit_001.adsorption.mcl_qmax = [0, 300.0,] # mol / m^3 (solid phase)
#### Colloidal
elif binding=='colloidal':
model.root.input.model.unit_001.adsorption_model = 'MULTI_COMPONENT_COLLOIDAL'
model.root.input.model.unit_001.adsorption.col_phi = 591000000.0
model.root.input.model.unit_001.adsorption.col_kappa_exp = 0
model.root.input.model.unit_001.adsorption.col_kappa_fact = 0
model.root.input.model.unit_001.adsorption.col_kappa_const = 4
model.root.input.model.unit_001.adsorption.col_cordnum = 6
model.root.input.model.unit_001.adsorption.col_logkeq_ph_exp = [0] * n_comp
model.root.input.model.unit_001.adsorption.col_logkeq_salt_powexp = [0] * n_comp
model.root.input.model.unit_001.adsorption.col_logkeq_salt_powfact = [0] * n_comp
model.root.input.model.unit_001.adsorption.col_logkeq_salt_expfact = [0, 30.6]
model.root.input.model.unit_001.adsorption.col_logkeq_salt_expargmult = [0, -20.9]
model.root.input.model.unit_001.adsorption.COL_BPP_PH_EXP = [0] * n_comp
model.root.input.model.unit_001.adsorption.COL_BPP_SALT_POWEXP = [0, -0.152]
model.root.input.model.unit_001.adsorption.COL_BPP_SALT_POWFACT = [0, 5.73]
model.root.input.model.unit_001.adsorption.COL_BPP_SALT_EXPFACT = [0] * n_comp
model.root.input.model.unit_001.adsorption.COL_BPP_SALT_EXPARGMULT = [0] * n_comp
model.root.input.model.unit_001.adsorption.col_radius = [1.78e-10, 4.488e-09]
model.root.input.model.unit_001.adsorption.col_kKin = [1e8] * n_comp
model.root.input.model.unit_001.adsorption.COL_LINEAR_THRESHOLD = 1e-8
model.root.input.model.unit_001.adsorption.COL_USE_PH = False
model.root.input.model.unit_001.discretization.nbound = [0, 1]
###########################################################################
model.root.input.model.unit_001.discretization.par_disc_type = 'EQUIDISTANT_PAR'
model.root.input.model.unit_001.discretization.use_analytic_jacobian = 0
model.root.input.model.unit_001.discretization.reconstruction = 'WENO'
model.root.input.model.unit_001.discretization.gs_type = 1
model.root.input.model.unit_001.discretization.max_krylov = 0
model.root.input.model.unit_001.discretization.max_restarts = 10
model.root.input.model.unit_001.discretization.schur_safety = 1e-8
model.root.input.model.unit_001.discretization.weno.boundary_model = 0
model.root.input.model.unit_001.discretization.weno.weno_eps = 1e-10
model.root.input.model.unit_001.discretization.weno.weno_order = 3
## Outlet Model ############################################################
model.root.input.model.unit_002.unit_type = 'OUTLET'
model.root.input.model.unit_002.ncomp = n_comp
############################################################################
# Sections
model.root.input.solver.sections.nsec = 3
model.root.input.solver.sections.section_times = [0.0, 1440.4, 2040.4, 2640.4] # s
model.root.input.solver.sections.section_continuity = [0, 0, 0]
# units: inlet > Column > Outlet
# Switches
model.root.input.model.connections.nswitches = 1
model.root.input.model.connections.switch_000.section = 0
model.root.input.model.connections.switch_000.connections = [
0, 1, -1, -1, 2.95e-09, # [unit_000, unit_002, all components, all components, Q/ m^3*s^-1
1, 2, -1, -1, 2.95e-09, # [unit_002, unit_004, all components, all components, Q/ m^3*s^-1
]
# Solver settings
model.root.input.model.solver.gs_type = 1
model.root.input.model.solver.max_krylov = 0
model.root.input.model.solver.max_restarts = 10
model.root.input.model.solver.schur_safety = 1e-8
# Number of cores for parallel simulation
model.root.input.solver.nthreads = 4
# Tolerances for the time integrator
model.root.input.solver.time_integrator.abstol = 1e-8
model.root.input.solver.time_integrator.algtol = 1e-12
model.root.input.solver.time_integrator.reltol = 1e-6
model.root.input.solver.time_integrator.init_step_size = 1e-6
model.root.input.solver.time_integrator.max_steps = 1000000
# Return data
# dont return data except for outlet unit
model.root.input['return'].split_components_data = 0
model.root.input['return'].split_ports_data = 0
model.root.input['return'].unit_000.write_solution_bulk = 0
model.root.input['return'].unit_000.write_solution_inlet = 1
model.root.input['return'].unit_000.write_solution_outlet = 0
model.root.input['return'].unit_001.write_solution_bulk = 0
model.root.input['return'].unit_001.write_solution_inlet = 0
model.root.input['return'].unit_001.write_solution_outlet = 0
# return outlet data for outlet unit
model.root.input['return'].unit_002.write_solution_bulk = 0
model.root.input['return'].unit_002.write_solution_inlet = 0
model.root.input['return'].unit_002.write_solution_outlet = 1
# Solution times
model.root.input.solver.user_solution_times = np.linspace(0, 2640.4, 5000)
# Save and run simulation
model.filename = 'cadet_version_check.h5'
model.save()
data = model.run()
if data.returncode == 0:
model.load()
else:
print(data)
raise Exception('Simulation failed')
sim_solution = model.root.output.solution
return sim_solution
def set_CADET_path(cadet_bin_path):
if platform.system() == 'Windows':
cadet_path = cadet_bin_path / 'cadet-cli.exe'
lwe_path = cadet_bin_path / 'createLWE.exe'
else:
cadet_path = cadet_bin_path / 'cadet-cli'
lwe_path = cadet_bin_path / 'createLWE'
if cadet_path.exists() and lwe_path.exists():
Cadet.cadet_path = cadet_path.as_posix()
elif cadet_path.exists() and not lwe_path.exists():
print('CADET was found but createLWE.exe was not found. Please make sure that none of the files have been moved.')
else:
print('CADET could not be found. Please check the bin path')
def compare(cadet_bin_path, binding):
CADET_solution = cadet_config(cadet_bin_path, binding)
CADET_time = CADET_solution.solution_times
inlet_array = CADET_solution.unit_000.solution_inlet
CADET_solution_array = CADET_solution.unit_002.solution_outlet
CADET_protein = CADET_solution_array.T[1]
fig, ax = plt.subplots()
ax.plot(CADET_time/60, CADET_solution_array, linewidth=1, color='k')
plt.show()
# ax.plot(CADET_time/60, CADET_protein, linewidth=1, color='k')
if __name__=='__main__':
cadet_bin_path = Path('/Users/angelamoser/miniconda3/envs/cadet4-env/bin')
compare(cadet_bin_path, binding='colloidal')