Replicating example 2 in CADET python tutorial

Hi, I tried to simulate example 2 in Cadet Python but encountered the following error. Kindly help as I am new to CADET

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[63], line 153
    151 time = langmuir_model.root.output.solution.solution_times
    152 c = langmuir_model.root.output.solution.unit_001.solution_outlet
--> 153 plt.plot(time, c)
    154 plt.title('Column (Outlet)')
    155 plt.xlabel('$time~/~min$')

File /opt/anaconda3/envs/cadet_env/lib/python3.11/site-packages/matplotlib/pyplot.py:3794, in plot(scalex, scaley, data, *args, **kwargs)
   3786 @_copy_docstring_and_deprecators(Axes.plot)
   3787 def plot(
   3788     *args: float | ArrayLike | str,
   (...)
   3792     **kwargs,
   3793 ) -> list[Line2D]:
-> 3794     return gca().plot(
   3795         *args,
   3796         scalex=scalex,
   3797         scaley=scaley,
   3798         **({"data": data} if data is not None else {}),
   3799         **kwargs,
   3800     )

File /opt/anaconda3/envs/cadet_env/lib/python3.11/site-packages/matplotlib/axes/_axes.py:1779, in Axes.plot(self, scalex, scaley, data, *args, **kwargs)
...
-> 2387     return x.to_numpy()
   2388 if hasattr(x, 'values'):
   2389     xtmp = x.values

TypeError: 'Dict' object is not callable

Hi Kauthar Maalim,

welcome to the CADET Forum.

Based on the code you shared I guess you were trying to run exercise 2 of Section 5:

Is that correct?

If so, please note the deprecation warning on the landing page. We are no longer recommending CADET-Python to new users but instead the new and improved CADET-Process interface, which has it’s own tutorial.

Regarding your error message: I can not reproduce it, on my machine the code runs fine.

Hi,

I use CADET python because I want to model radial flow chromatography which is not available in CADET process

I managed to solve an issue, but now I have a problem after adapting the code to the radial flow chromatography ( load and elute)

def breakthrough_radial_langmuir(Q,  load_duration, c=10, wash_elute_duration=110, kDa=150):
    
    Q_ = Q.to_base_units().magnitude
    c_ = ConcentrationConverter(c, kDa=kDa).mM.to_base_units().magnitude
    load_duration = round(load_duration.to_base_units().magnitude)

    langmuir_model = get_cadet_template(n_units=3) # 3 = Inlet + Column + Outlet

        # INLET
    langmuir_model.root.input.model.unit_000.unit_type = 'INLET'
    langmuir_model.root.input.model.unit_000.ncomp = 2
    langmuir_model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
    langmuir_model.root.input.model.unit_000.sec_000.const_coeff = [c_,] # 10mg/ml
    langmuir_model.root.input.model.unit_000.sec_001.const_coeff = [0.0,] # Wash & Elute step

        # Column
    langmuir_model.root.input.model.unit_001.unit_type = 'RADIAL_LUMPED_RATE_MODEL_WITH_PORES'
    langmuir_model.root.input.model.unit_001.ncomp = 2

    langmuir_model.root.input.model.unit_001.col_length = 0.12
    langmuir_model.root.input.model.unit_001.col_radius_inner = 0.0103 * 0.2
    langmuir_model.root.input.model.unit_001.col_radius_outer = 0.0103
    langmuir_model.root.input.model.unit_001.col_porosity = 0.29
    langmuir_model.root.input.model.unit_001.par_porosity = 0.95
    langmuir_model.root.input.model.unit_001.par_radius = 0.00025
        
    langmuir_model.root.input.model.unit_001.col_dispersion = 1e-8
    langmuir_model.root.input.model.unit_001.film_diffusion = [6.9e-6,]
        
    langmuir_model.root.input.model.unit_001.adsorption_model = 'MULTI_COMPONENT_LANGMUIR'

    langmuir_model.root.input.model.unit_001.adsorption.is_kinetic = 1
    langmuir_model.root.input.model.unit_001.adsorption.mcl_ka = [1.,]
    langmuir_model.root.input.model.unit_001.adsorption.mcl_kd = [1.,]
    langmuir_model.root.input.model.unit_001.adsorption.mcl_qmax = [0.1,] ############??????

    langmuir_model.root.input.model.unit_001.init_c = [0.0,]
    langmuir_model.root.input.model.unit_001.init_q = [0.0,]

        ## Outlet
    langmuir_model.root.input.model.unit_002.ncomp = 2
    langmuir_model.root.input.model.unit_002.unit_type = 'OUTLET'
        
    set_discretization(langmuir_model, n_bound=[1,])

        # Sections and connections
    langmuir_model.root.input.solver.sections.nsec = 1
    langmuir_model.root.input.solver.sections.section_times = [0.0, 50]
        
        # LOAD + (WASH + ELUTE) step
    langmuir_model.root.input.solver.sections.nsec = 2
    langmuir_model.root.input.solver.sections.section_times = [0.0, load_duration, wash_elute_duration]
        
        ## Inlet Profile
    langmuir_model.root.input.model.unit_000.sec_000.const_coeff = [c_,0] # 10mg/ml
    langmuir_model.root.input.model.unit_000.sec_001.const_coeff = [0.0,0] # Wash & Elute step
    langmuir_model.root.input.solver.sections.section_continuity = [0, 0]
        
        ## Switches
    langmuir_model.root.input.model.connections.nswitches = 1
    langmuir_model.root.input.model.connections.switch_000.section = 0
    langmuir_model.root.input.model.connections.switch_000.connections = [
        0, 1, -1, -1, Q_, # unit_000, unit_001, all components, all components, Q/ m^3/s
        1, 2, -1, -1, Q_, # unit_001, unit_002, all components, all components, Q/ m^3/s
    ]  

        # set the times that the simulator writes out data for
    langmuir_model.root.input.solver.user_solution_times = np.linspace(0, int(load_duration), int(load_duration)+1, int(wash_elute_duration), int(wash_elute_duration)+1) 

    return langmuir_model

run_simulation(langmuir_model)
time = langmuir_model.root.output.solution.solution_times
c1 = langmuir_model.root.output.solution.unit_002['solution_outlet_comp_000']
c2 = langmuir_model.root.output.solution.unit_002['solution_outlet_comp_001']
plt.plot(time, c1)
plt.plot(time, c2)
plt.title('Column (Outlet)')
plt.xlabel('$time~/~min$')
plt.ylabel('$concentration~/~mM$')
plt.show()  

Hi Kauthar,

I can not reproduce the error because the code you shared is not complete. Please include all imports and all function/class definitions (such as ConcentrationConverter).

import os
import platform
from pathlib import Path
from cadet import Cadet

from cadet import Cadet
from matplotlib import pyplot as plt
import numpy as np

import os

import numpy as np

import json

# Temporary files for simulation objects
import tempfile
tempfile.tempdir = os.path.join('')

def breakthrough_radial_langmuir(Q,  load_duration, c=10, wash_elute_duration=110, kDa=150):
    
    Q_ = Q.to_base_units().magnitude
    c_ = ConcentrationConverter(c, kDa=kDa).mM.to_base_units().magnitude
    load_duration = round(load_duration.to_base_units().magnitude)

    langmuir_model = get_cadet_template(n_units=3) # 3 = Inlet + Column + Outlet

        # INLET
    langmuir_model.root.input.model.unit_000.unit_type = 'INLET'
    langmuir_model.root.input.model.unit_000.ncomp = 2
    langmuir_model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
    langmuir_model.root.input.model.unit_000.sec_000.const_coeff = [c_,] # 10mg/ml
    langmuir_model.root.input.model.unit_000.sec_001.const_coeff = [0.0,] # Wash & Elute step

        # Column
    langmuir_model.root.input.model.unit_001.unit_type = 'RADIAL_LUMPED_RATE_MODEL_WITH_PORES'
    langmuir_model.root.input.model.unit_001.ncomp = 2

    langmuir_model.root.input.model.unit_001.col_length = 0.12
    langmuir_model.root.input.model.unit_001.col_radius_inner = 0.0103 * 0.2
    langmuir_model.root.input.model.unit_001.col_radius_outer = 0.0103
    langmuir_model.root.input.model.unit_001.col_porosity = 0.29
    langmuir_model.root.input.model.unit_001.par_porosity = 0.95
    langmuir_model.root.input.model.unit_001.par_radius = 0.00025
        
    langmuir_model.root.input.model.unit_001.col_dispersion = 1e-8
    langmuir_model.root.input.model.unit_001.film_diffusion = [6.9e-6,]
        
    langmuir_model.root.input.model.unit_001.adsorption_model = 'MULTI_COMPONENT_LANGMUIR'

    langmuir_model.root.input.model.unit_001.adsorption.is_kinetic = 1
    langmuir_model.root.input.model.unit_001.adsorption.mcl_ka = [1.,]
    langmuir_model.root.input.model.unit_001.adsorption.mcl_kd = [1.,]
    langmuir_model.root.input.model.unit_001.adsorption.mcl_qmax = [0.1,] ############??????

    langmuir_model.root.input.model.unit_001.init_c = [0.0,]
    langmuir_model.root.input.model.unit_001.init_q = [0.0,]

        ## Outlet
    langmuir_model.root.input.model.unit_002.ncomp = 2
    langmuir_model.root.input.model.unit_002.unit_type = 'OUTLET'
        
    set_discretization(langmuir_model, n_bound=[1,])

        # Sections and connections
    langmuir_model.root.input.solver.sections.nsec = 1
    langmuir_model.root.input.solver.sections.section_times = [0.0, 50]
        
        # LOAD + (WASH + ELUTE) step
    langmuir_model.root.input.solver.sections.nsec = 2
    langmuir_model.root.input.solver.sections.section_times = [0.0, load_duration, wash_elute_duration]
        
        ## Inlet Profile
    langmuir_model.root.input.model.unit_000.sec_000.const_coeff = [c_,0] # 10mg/ml
    langmuir_model.root.input.model.unit_000.sec_001.const_coeff = [0.0,0] # Wash & Elute step
    langmuir_model.root.input.solver.sections.section_continuity = [0, 0]
        
        ## Switches
    langmuir_model.root.input.model.connections.nswitches = 1
    langmuir_model.root.input.model.connections.switch_000.section = 0
    langmuir_model.root.input.model.connections.switch_000.connections = [
        0, 1, -1, -1, Q_, # unit_000, unit_001, all components, all components, Q/ m^3/s
        1, 2, -1, -1, Q_, # unit_001, unit_002, all components, all components, Q/ m^3/s
    ]  

        # set the times that the simulator writes out data for
    langmuir_model.root.input.solver.user_solution_times = np.linspace(0, int(load_duration), int(load_duration)+1, int(wash_elute_duration), int(wash_elute_duration)+1) 

    return langmuir_model

run_simulation(langmuir_model)
time = langmuir_model.root.output.solution.solution_times
c1 = langmuir_model.root.output.solution.unit_002['solution_outlet_comp_000']
c2 = langmuir_model.root.output.solution.unit_002['solution_outlet_comp_001']
plt.plot(time, c1)
plt.plot(time, c2)
plt.title('Column (Outlet)')
plt.xlabel('$time~/~min$')
plt.ylabel('$concentration~/~mM$')
plt.show()

Hi please find the updated code

Hi,

that code also didn’t include the import for ConcentrationConverter that I had asked for, but I made it work anyways with some guesswork. Please ensure your code runs as an isolated minimal reproducible example in a fresh Python Kernel. Also, please include the error message you get next time.

With that said, here are the errors in your code:

  1. Missing call to breakthrough_radial_langmuir(): You define breakthrough_radial_langmuir() but never call it. Therefore your code does not know what langmuir_model is in the global scope. If the reason is unclear please check out an article on scopes, e.g. this. This can be fixed with langmuir_model = breakthrough_radial_langmuir(Q=2e-6, load_duration=15)
  2. Inconsistent number of components & parameters. You specified a lot of parameters for one component (mcl_ka = [1., ], init_c = [0.0, ], etc.) but told the unit operation to expect two components unit_001.ncomp = 2. To fix: set n_comp = 1 or give all parameters that require a parameter per component two values mcl_ka = [1.0, 2.0] for example.
  3. Inproper arguments for np.linspace: You gave 5 arguments to np.linspace: np.linspace(0, int(load_duration), int(load_duration)+1, int(wash_elute_duration), int(wash_elute_duration)+1). np.linspace commonly takes three: start, stop, num (number of entries on the linear scale). So np.linspace(start=0, stop=int(wash_elute_duration), num=300) is I think what you wanted.
  4. Accessing component solutions without splitting components: You try to access the solution with unit_002['solution_outlet_comp_000']. Because you didn’t include the definition of get_cadet_template() I had to guess that you used the default one from the CADET-Python-Tutorial. This function does not set the split_components_data flag to True by default. Therefore unit_002['solution_outlet_comp_000'] is going to be empty. You can access the solution with unit_002['solution_outlet'] or use langmuir_model = get_cadet_template(n_units=3, split_components_data=True) with unit_002['solution_outlet_comp_000'].
  5. Missing unit operation type in set_discretization(): Again, you didn’t include the definition of set_discretization() so I can only guess. But if it’s the function from the CADET-Python-Tutorial, then it only sets discretizations for columns = {'GENERAL_RATE_MODEL', 'LUMPED_RATE_MODEL_WITH_PORES', 'LUMPED_RATE_MODEL_WITHOUT_PORES'} and not RADIAL_LUMPED_RATE_MODEL_WITH_PORES. This unit operation type needs to be added to that function.
  6. Some code duplications These didn’t break anything but still made reading and debugging harder, as they set the same variable in multiple places. Such as: line 33: langmuir_model.root.input.model.unit_000.sec_000.const_coeff = [c_,] # 10mg/ml langmuir_model.root.input.model.unit_000.sec_001.const_coeff = [0.0,] # Wash & Elute step
    line 74:
    langmuir_model.root.input.model.unit_000.sec_000.const_coeff = [c_,0] # 10mg/ml langmuir_model.root.input.model.unit_000.sec_001.const_coeff = [0.0,0] # Wash & Elute step

Here is a complete and working code example:

import os
import tempfile

import numpy as np
from matplotlib import pyplot as plt

from cadet import Cadet

tempfile.tempdir = os.path.join('')

Cadet.cadet_path = SET_CADET_PATH_HERE


def get_cadet_template(n_units=3, split_components_data=False):
    cadet_template = Cadet()

    cadet_template.root.input.model.nunits = n_units

    # Store solution
    cadet_template.root.input['return'].split_components_data = split_components_data
    cadet_template.root.input['return'].split_ports_data = 0
    cadet_template.root.input['return'].unit_000.write_solution_inlet = 1
    cadet_template.root.input['return'].unit_000.write_solution_outlet = 1
    cadet_template.root.input['return'].unit_000.write_solution_bulk = 1
    cadet_template.root.input['return'].unit_000.write_solution_particle = 1
    cadet_template.root.input['return'].unit_000.write_solution_solid = 1
    cadet_template.root.input['return'].unit_000.write_solution_flux = 1
    cadet_template.root.input['return'].unit_000.write_solution_volume = 1
    cadet_template.root.input['return'].unit_000.write_coordinates = 1
    cadet_template.root.input['return'].unit_000.write_sens_outlet = 1

    for unit in range(n_units):
        cadet_template.root.input['return']['unit_{0:03d}'.format(unit)] = cadet_template.root.input['return'].unit_000

    # Tolerances for the time integrator
    cadet_template.root.input.solver.time_integrator.abstol = 1e-6
    cadet_template.root.input.solver.time_integrator.algtol = 1e-10
    cadet_template.root.input.solver.time_integrator.reltol = 1e-6
    cadet_template.root.input.solver.time_integrator.init_step_size = 1e-6
    cadet_template.root.input.solver.time_integrator.max_steps = 1000000

    # Solver settings
    cadet_template.root.input.model.solver.gs_type = 1
    cadet_template.root.input.model.solver.max_krylov = 0
    cadet_template.root.input.model.solver.max_restarts = 10
    cadet_template.root.input.model.solver.schur_safety = 1e-8

    # Run the simulation on single thread
    cadet_template.root.input.solver.nthreads = 1

    return cadet_template


def set_discretization(model, n_bound=None, n_col=20, n_par_types=1):
    columns = {'GENERAL_RATE_MODEL', 'LUMPED_RATE_MODEL_WITH_PORES', 'LUMPED_RATE_MODEL_WITHOUT_PORES',
               "RADIAL_LUMPED_RATE_MODEL_WITH_PORES"}

    for unit_name, unit in model.root.input.model.items():
        if 'unit_' in unit_name and unit.unit_type in columns:
            unit.discretization.ncol = n_col
            unit.discretization.npar = 5
            unit.discretization.npartype = n_par_types

            if n_bound is None:
                n_bound = unit.ncomp * [0]
            unit.discretization.nbound = n_bound

            unit.discretization.par_disc_type = 'EQUIDISTANT_PAR'
            unit.discretization.use_analytic_jacobian = 1
            unit.discretization.reconstruction = 'WENO'
            unit.discretization.WENO_ORDER = 1
            unit.discretization.gs_type = 1
            unit.discretization.max_krylov = 0
            unit.discretization.max_restarts = 10
            unit.discretization.schur_safety = 1.0e-8

            unit.discretization.weno.boundary_model = 0
            unit.discretization.weno.weno_eps = 1e-10
            unit.discretization.weno.weno_order = 3


def run_simulation(cadet, file_name=None):
    if file_name is None:
        f = next(tempfile._get_candidate_names())
        cadet.filename = os.path.join(tempfile.tempdir, f + '.h5')
    else:
        cadet.filename = file_name
    # save the simulation
    cadet.save()

    # run the simulation and load results
    data = cadet.run(timeout=10)
    print(data)
    cadet.load()

    # Remove files
    if file_name is None:
        os.remove(os.path.join(tempfile.tempdir, f + '.h5'))


def breakthrough_radial_langmuir(Q: float, load_duration: float, c=10, wash_elute_duration=110, kDa=150):
    langmuir_model = get_cadet_template(n_units=3, split_components_data=True)  # 3 = Inlet + Column + Outlet

    # INLET
    langmuir_model.root.input.model.unit_000.unit_type = 'INLET'
    langmuir_model.root.input.model.unit_000.ncomp = 2
    langmuir_model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'

    # Column
    langmuir_model.root.input.model.unit_001.unit_type = 'RADIAL_LUMPED_RATE_MODEL_WITH_PORES'
    langmuir_model.root.input.model.unit_001.ncomp = 2

    langmuir_model.root.input.model.unit_001.col_length = 0.12
    langmuir_model.root.input.model.unit_001.col_radius_inner = 0.0103 * 0.2
    langmuir_model.root.input.model.unit_001.col_radius_outer = 0.0103
    langmuir_model.root.input.model.unit_001.col_porosity = 0.29
    langmuir_model.root.input.model.unit_001.par_porosity = 0.95
    langmuir_model.root.input.model.unit_001.par_radius = 0.00025
    langmuir_model.root.input.model.unit_001.VELOCITY_COEFF = 0.1
    langmuir_model.root.input.model.unit_001.CROSS_SECTION_AREA = 0.001

    langmuir_model.root.input.model.unit_001.col_dispersion = 1e-8
    langmuir_model.root.input.model.unit_001.film_diffusion = [6.9e-6, 6.9e-6]

    langmuir_model.root.input.model.unit_001.adsorption_model = 'MULTI_COMPONENT_LANGMUIR'

    langmuir_model.root.input.model.unit_001.adsorption.is_kinetic = 1
    langmuir_model.root.input.model.unit_001.adsorption.mcl_ka = [1., 2, ]
    langmuir_model.root.input.model.unit_001.adsorption.mcl_kd = [1., 1, ]
    langmuir_model.root.input.model.unit_001.adsorption.mcl_qmax = [100., 50., ]

    langmuir_model.root.input.model.unit_001.init_c = [0.0, 0.0]
    langmuir_model.root.input.model.unit_001.init_q = [0.0, 0.0]

    ## Outlet
    langmuir_model.root.input.model.unit_002.ncomp = 2
    langmuir_model.root.input.model.unit_002.unit_type = 'OUTLET'

    set_discretization(langmuir_model, n_bound=[0, 1], n_col=50)

    # LOAD + (WASH + ELUTE) step
    langmuir_model.root.input.solver.sections.nsec = 2
    langmuir_model.root.input.solver.sections.section_times = [0.0, load_duration, wash_elute_duration]

    ## Inlet Profile
    langmuir_model.root.input.model.unit_000.sec_000.const_coeff = [c, c]  # 10mg/ml
    langmuir_model.root.input.model.unit_000.sec_001.const_coeff = [0.0, 0.0]  # Wash & Elute step
    langmuir_model.root.input.solver.sections.section_continuity = [0, 0]

    ## Switches
    langmuir_model.root.input.model.connections.nswitches = 1
    langmuir_model.root.input.model.connections.switch_000.section = 0
    langmuir_model.root.input.model.connections.switch_000.connections = [
        0, 1, -1, -1, Q,  # unit_000, unit_001, all components, all components, Q/ m^3/s
        1, 2, -1, -1, Q,  # unit_001, unit_002, all components, all components, Q/ m^3/s
    ]

    # set the times that the simulator writes out data for
    langmuir_model.root.input.solver.user_solution_times = np.linspace(start=0, stop=int(wash_elute_duration), num=300)

    return langmuir_model


langmuir_model = breakthrough_radial_langmuir(Q=2e-6, load_duration=15)
run_simulation(langmuir_model)
time = langmuir_model.root.output.solution.solution_times
c1 = langmuir_model.root.output.solution.unit_002['solution_outlet_comp_000']
c2 = langmuir_model.root.output.solution.unit_002['solution_outlet_comp_001']
plt.plot(time, c1)
plt.plot(time, c2)
plt.title('Column (Outlet)')
plt.xlabel('$time~/~min$')
plt.ylabel('$concentration~/~mM$')
plt.show()
3 Likes

Hi Ronald,

Thanks for the feedback and corrections. I have another question: In radial flow chromatography, the velocity changes with column radius, and the column radius is not constant (the inlet radius is larger than the outlet). How does CADET consider these two features?

Hi Kauthar,

The required input parameters for radial flow models are described here, modeling is considered here.

The interstitial column velocity is calculated from the flow rate (that you specify for the column inlet) and the cylindrical geometry. This computation additionally relies on the column height. The computed velocity decreases towards the outlet and is defined as

u(r) := \frac{v}{r}

with v being some constant and r is the radial position. The constant is calculated as described above via

v = F_{in} / (2.0 * \pi * L * \varepsilon)

where L is the column height, F_{in} is the inlet flow rate and \varepsilon the column porosity.

Alternatively, you can specify the velocity ceofficient v directly and instead not specify the column height. If height is provided, the input of v will be ignored.

The derivation of the above velocity for cylinder coordinates and solely in radial direction can be found here

Best,
Jan

1 Like