Bi-Langmuir q-profiles

Hello,

I am trying to capture the data on the solid phase concentrations for a case with A and B under various inlet conditions. I can produce plots for outlet concentration but I cannot work out how to get the data of q. My working code is:

#%% (A)

import numpy as np
import matplotlib.pyplot as plt

from CADETProcess.processModel import ComponentSystem
from CADETProcess.processModel import BiLangmuir
from CADETProcess.processModel import (
    Inlet, LumpedRateModelWithoutPores, Outlet
)
from CADETProcess.processModel import FlowSheet, Process

C_in = [[7.5, 7.5], [10, 10], [15, 15]]

c_store = []
q_store = []

for i in range(len(C_in)):

    component_system = ComponentSystem()
    component_system.add_component('A')
    component_system.add_component('B')

    # --------------------------------------------------------
    # Bi-Langmuir binding model
    # --------------------------------------------------------

    binding_model = BiLangmuir(component_system, n_binding_sites=2, name='langmuir')
    binding_model.is_kinetic = True
        
    binding_model.adsorption_rate = [0.015, 0.03, 0.02, 0.035]
    binding_model.desorption_rate = [1, 1, 1, 1]
    binding_model.capacity = [100, 100, 100, 100]

    print('hi 1')
    # --------------------------------------------------------
    # Column system
    # --------------------------------------------------------

    feed = Inlet(component_system, name='feed')
    feed.c = C_in[i]

    eluent = Inlet(component_system, name='eluent')
    eluent.c = [0, 0]

    column = LumpedRateModelWithoutPores(component_system, name='column')
    column.binding_model = binding_model
    column.length = 0.6
    column.diameter = 0.024
    column.axial_dispersion = 4.7e-6
    column.total_porosity = 0.65

    column.solution_recorder.write_solution_bulk = True

    outlet = Outlet(component_system, name='outlet')

    flow_sheet = FlowSheet(component_system)

    flow_sheet.add_unit(feed, feed_inlet=True)
    flow_sheet.add_unit(eluent, eluent_inlet=True)
    flow_sheet.add_unit(column)
    flow_sheet.add_unit(outlet, product_outlet=True)

    flow_sheet.add_connection(feed, column)
    flow_sheet.add_connection(eluent, column)
    flow_sheet.add_connection(column, outlet)

    process = Process(flow_sheet, 'batch elution')
    print('hi 2')
    # --------------------------------------------------------
    # Events
    # --------------------------------------------------------

    Q = 60 / (60 * 1e6)

    process.add_event('feed_on', 'flow_sheet.feed.flow_rate', Q)
    process.add_event('feed_off', 'flow_sheet.feed.flow_rate', 0.0)

    process.add_event('eluent_on', 'flow_sheet.eluent.flow_rate', Q)
    process.add_event('eluent_off', 'flow_sheet.eluent.flow_rate', 0.0)

    process.add_duration('feed_duration')

    process.add_event_dependency('eluent_on', ['feed_off'])
    process.add_event_dependency('eluent_off', ['feed_on'])
    process.add_event_dependency('feed_off', ['feed_on','feed_duration'], [1,1])

    process.cycle_time = 750
    process.feed_duration.time = 60
    print('hi 3')
    # --------------------------------------------------------
    # Run CADET
    # --------------------------------------------------------
    print(binding_model.n_binding_sites)
    print(binding_model.n_bound_states)
    if __name__ == "__main__":
        from CADETProcess.simulator import Cadet
        print('hi 4')
        process_simulator = Cadet()
        print('hi 5')
        simulation_results = process_simulator.simulate(process)
        print('hi 6')
    # --------------------------------------------------------
    # Extract outlet solution
    # --------------------------------------------------------

    outlet_solution = simulation_results.solution.column.outlet
    solid_solution = simulation_results.solution.column.solid

    t_full = outlet_solution.time
    c = outlet_solution.solution

    c_store.append(c)

    solid_solution = simulation_results.solution.column.solid
    q_store.append(solid_solution.solution)

c_store = np.array(c_store)
#q_store = np.array(q_store)

print("DATA GENERATED\n")

plt.figure(figsize=(8,6))

for i in range(len(C_in)):

    c = c_store[i]

    cA = c[:,0]
    cB = c[:,1]

    plt.plot(t_full, cA, label=f"A dataset {i}")
    plt.plot(t_full, cB, '--', label=f"B dataset {i}")

plt.xlabel("Time (s)")
plt.ylabel("Outlet concentration")
plt.title("Outlet Chromatograms (Bi-Langmuir)")
plt.legend()
plt.grid()

plt.show()

When I add the line:

column.solution_recorder.write_solution_solid = True

I get the error:

ValueError: Expected size (751, 100, 2)

Which comes after (during the first iteration of the loop):

hi 1
hi 2
hi 3
2
4
hi 4
hi 5

I have also tried using (to the same result):

column.solution_recorder.write_solution_bound

I put all the print(‘hi n’) in to indicate where the code reaches. Note: hi 6 is never printed when i get the error.

Any ideas?

Thanks!