Difference in colloidal isotherm behavior between v4.4 and v5.0

a.moser · October 30, 2024, 11:12am

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')

ronald.jaepel · October 30, 2024, 11:13am

Good catch, thanks for the report. I’ll look into it.

ronald.jaepel · October 30, 2024, 12:13pm

Ah, yes, now I remember

We did change it and didn’t list the change because the isotherm was a “new” addition in v5. Sorry about that. I hope it didn’t eat too much of your time.

Here are the changes, you can see them in the src/libcadet/model/binding/ColloidalBinding.cpp file in this PullRequest

Renamed COL_RADIUS to COL_PROTEIN_RADIUS
Fixed kappa factors not adhering to SI units by changing it from \kappa = \frac{10^9}{\kappa_f c_0^{\kappa_e} + \kappa_c} to \kappa = \frac{1}{\kappa_f c_0^{\kappa_e} + \kappa_c}. So you’ll need to adapt your COL_KAPPA_FACT and COL_KAPPA_CONST by 10^9.
Fixed a bug in the analytical Jacobian, so it works now and can speed up the simulation time if enabled.

a.moser · October 30, 2024, 3:34pm

Good to know, thanks. Got the expected results after changing the kappa units.