Help debugging GIEX

[moritz.imend]

Expected behavior.

Hi,

I am trying to extend an example from SMA to GIEX using cadet process. I am running it on ubuntu 20.04, with cadet 3.04 (using conda) and Cadet Process as 0.11.1. I would expect the GIEX example to run without major issues for this simple example.

Actual behavior

It seems that the solver gets stuck if I slightly change the linear dependence of the characteristic charge. e.g. from -0.2 to -0.3, although the function should return positive values over the chosen pH range. If I use a positive slope, the solver seems much more robust.

without the pH dependence of the characteristic charge the solver works well with the standard settings.

Can you please help me to debug this issue? E.g. find suitable paramters for the solver?

How to produce bug (including a minimal reproducible example)

import numpy as np
from CADETProcess.processModel import (
    ComponentSystem,
    Inlet,
    GeneralRateModel,
    GeneralizedIonExchange,
    Outlet,
    FlowSheet,
    Process,
)
from CADETProcess.simulator import Cadet  


import time
def calculate_column_volume(diameter_cm: float, length_cm: float) -> float:
    """Calculate column volume in mL."""
    return np.pi * diameter_cm**2/4 * length_cm

"""Create CADET Process for Load-Wash-Elute operation."""

inputs = {
    "temperature": 20,
    "load_ph": 3.0,
    "load_salt": 2.0,  #mM
    "wash_ph": 4.0,
    "wash_salt": 5.0,
    "elution_ph": 8.0,
    "elution_salt": 1000.0,
    "flow_rate": 2.0,  # mL/min
    "equilibration_volume": 3.0, #CV
    "load_volume": 3.0,  #CV
    "sample_concentration": 1.0,  # g/L
    "required_purity": 0.7,
    "wash_volume": 3,
    "elution_volume": 3.0,
    "column_length": 20,  # cm
    "column_diameter": 10.0,  # cm
    "particle_size": 45,  # µm
    "strip_volume": 3.0
}

    
    
# === Component System ===
component_system = ComponentSystem()
component_system.add_component("Salt")
component_system.add_component("pH Modifier")
component_system.add_component("IgG") 

# Binding Model
binding_model = GeneralizedIonExchange(component_system, name="GIEX")
binding_model.is_kinetic = True
binding_model.adsorption_rate = [0.0, 0.0, 5]
binding_model.adsorption_rate_linear = [0.0, 0.0, 0]
binding_model.adsorption_rate_quadratic = [0.0, 0.0, 0]
binding_model.desorption_rate = [0.0, 0.0, 1]
pcoeff = [-0.3, 6.0]  #  +0.3 works well
print(f"pH between 2 and 10: {np.polyval(pcoeff, np.arange(2,10))}")
# https://forum.cadet-web.de/t/giex-setup-in-cadet-process/867/3
binding_model.characteristic_charge_breaks = [3, 8, 3, 8, 3, 8]
binding_model.characteristic_charge = [0.0, 0.0, pcoeff[-1]]
binding_model.characteristic_charge_linear = [0, 0, pcoeff[0]]
#binding_model.characteristic_charge_quadratic = [0, 0, 0.14]
binding_model.steric_factor = [0, 0.0, 10]
binding_model.capacity = 1200.0
binding_model.reference_liquid_phase_conc = 1000
binding_model.reference_solid_phase_conc = 1200.0

# Create Units
# Inlet
inlet = Inlet(component_system, name="inlet")
flow_rate_m3_s = inputs.get("flow_rate", 1.0) * 1e-6 / 60.0  # mL/min to m³/s
inlet.flow_rate = flow_rate_m3_s

# Column
column = GeneralRateModel(component_system, name="column")
column.binding_model = binding_model

# Column geometry
column.length = inputs.get("column_length") / 100.0  # cm to m
column.diameter = inputs.get("column_diameter") / 100.0  # cm to m
column.bed_porosity = 0.37
column.particle_radius = inputs.get("particle_size") * 1e-6 / 2.0  # μm diameter to m radius
column.particle_porosity = 0.75

# Transport parameters (adjusted for temperature)
column.axial_dispersion = 5.75e-8
column.film_diffusion = [6.9e-6] * column.n_comp
column.pore_diffusion = [1e-9, 1e-9, 6e-10]
column.surface_diffusion = [0.0] * column.n_bound_states

# Initial conditions (equilibrated with load buffer)
load_salt = inputs.get("load_salt")
load_ph = inputs.get("load_ph")

column.c = [load_salt, load_ph, 0]
column.cp = [load_salt, load_ph, 0]

column.q = [binding_model.capacity, 0]
# Outlet
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(outlet, product_outlet=True)

# Connections
flow_sheet.add_connection(inlet, column)
flow_sheet.add_connection(column, outlet)

# Process

process = Process(flow_sheet, name="LWE")

# Calculate step durations
column_volume = calculate_column_volume(
    column.diameter * 100,  # m to cm
    column.length * 100  # m to cm
)  # mL
flow_rate_mL_min = inputs.get("flow_rate")

# Step volumes
load_volume_mL = inputs.get("load_volume") * column_volume
wash_volume_mL = inputs.get("wash_volume") * column_volume
elution_volume_mL = inputs.get("elution_volume") * column_volume
strip_volume_mL = inputs.get("strip_volume") * column_volume
equilibration_volume_mL = inputs.get("equilibration_volume") * column_volume

# Step durations in seconds
equilibration_duration = (equilibration_volume_mL / flow_rate_mL_min) * 60
load_duration = (load_volume_mL / flow_rate_mL_min) * 60
wash_duration = (wash_volume_mL / flow_rate_mL_min) * 60
elution_duration = (elution_volume_mL / flow_rate_mL_min) * 60
strip_duration = (strip_volume_mL / flow_rate_mL_min) * 60

# Total cycle time
process.cycle_time = (equilibration_duration + load_duration + wash_duration + 
                        elution_duration + strip_duration)

# Define Event Times
event_times = {
    'Equilibration': 0.0,
    'Load': equilibration_duration,
    'Wash': equilibration_duration + load_duration,
    'Elution': equilibration_duration + load_duration + wash_duration,
    'Strip': equilibration_duration + load_duration + wash_duration + elution_duration
}


# Get buffer parameters
load_salt_conc = inputs.get("load_salt")
wash_salt_conc = inputs.get("wash_salt")
elution_salt_conc = inputs.get("elution_salt")

sample_conc = inputs.get("sample_concentration")

load_ph = inputs.get("load_ph")
wash_ph = inputs.get("wash_ph")
elution_ph = inputs.get("elution_ph")

c_equilibration = [load_salt_conc, load_ph, 0.0]
c_load = [
    load_salt_conc,
    load_ph,
    sample_conc
]


c_wash = [wash_salt_conc, wash_ph, 0.0]
c_elution = [elution_salt_conc, elution_ph, 0.0]


# gradient slope from wash conditions to end of elution conditions
gradient_slope = (np.array(c_elution) - np.array(c_wash)) / elution_duration

# pH will change within X seconds from wash to elution pH
ph_gradient_duration = 60  # seconds
gradient_slope_pH = (np.array(c_elution)[1] - np.array(c_wash)[1]) / ph_gradient_duration

# first, gradient with accelerated pH conditions
c_gradient_poly_ph = np.array(list(zip(c_wash, gradient_slope))) # end of pH gradient, set pH to elution pH
c_gradient_poly_ph[1, 1] = gradient_slope_pH

# second, gradient only for salt
c_gradient_poly = np.array(list(zip(c_wash, gradient_slope)))
c_gradient_poly[1, 1] = 0
c_gradient_poly[1, 0] = c_elution[1]

# Set inlet concentration for each step
process.add_event('Equilibration', 'flow_sheet.inlet.c', c_equilibration)
process.add_event('Load', 'flow_sheet.inlet.c', c_load, time=event_times['Load'])
process.add_event('Wash', 'flow_sheet.inlet.c', c_wash, time=event_times['Wash'])
process.add_event('grad_start_pH', 'flow_sheet.inlet.c', c_gradient_poly_ph, event_times['Elution'])
process.add_event('grad_start_salt', 'flow_sheet.inlet.c', c_gradient_poly, event_times['Elution']+ph_gradient_duration)
process.add_event('grad_stop', 'flow_sheet.inlet.c', c_elution, time=event_times['Strip'])

#process.add_event('Strip', 'flow_sheet.inlet.c', c_strip, time=event_times['Strip'])
cadet = Cadet(install_path="/home/moritz/miniconda3/envs/cadet/bin/cadet-cli")
cadet.save_to_h5(process, "test.hdf5")

  
start = time.time()
results = cadet.simulate(process)
print(f"Elapsed time: {time.time() - start}")

from CADETProcess.plotting import SecondaryAxis
from matplotlib import pyplot as plt


sec = SecondaryAxis()
sec.components = ['Salt', "pH"]
sec.y_label = '$c_{salt}, pH$'

fig = results.solution.column.outlet.plot(secondary_axis=sec)[0]
plt.savefig("sol_simple1.png") 

fig = results.solution.column.inlet.plot(secondary_axis=sec)[0]
plt.savefig("inlet_simple.png")

File produced by conda env export > environment.yml

environment.yml (6.5 KB)

[alters]

Some suggestions:

You are using reference concentrations but your values for kA have not been adjusted. @a.moser made a nice post recently about this, I recommend looking at it.

You are not using reference pHs - this is always advised when using pH dependent adsorption models. For example if your ref pH is 5.5, load pH of 5 would be: 5 - 5.5 = -0.5. This really helps with scaling the problem, making it easier for the solver and also easier to interpret parameter values. The ref pH is usually the centerpoint of the pH range or one of the edges, depending on what works best. It is helpful to do some example calculations of the kA and charge to define the ideal parameter ranges.

Also, I think in general pH gradients are difficult to model. If the data is at a constant salt concentration but varying pH, you might be better off using a purely pH dependent model like MPM langmuir or colloidal. There have been some forum posts recently on doing this for proA, but the functional forms could also work for IEC. If you have multiple salt and pH, then GIEX is a good option. I think the SMA type models work much better if there is a dependence on salt concentration in the data.