Cyclic steady state in CADET

CADET Troubleshooting
1

Hi there,

I found this topic where you dicuss the cyclic steady state for SMBs in CADET but I have encountered a problem ((Achieving cyclic steady state in CADET - #4 by j.schmoelder).
I use the last_state_y and last_state_ydot as initial states for a new simulation as stated in the discussion above. However, I have encountered that when a new simulation is run, there is an issue with some of the concentrations. In the attached simulation files and results as images, it is shown that when the new simulation starts, the concentration increases. However, this happens only for some concentrations (not shown).
sim_langmuir_file
When I run the simulation for the same time but ‘in one piece’, these increases are not seen.
sim_langmuir_file_2

I have not been able to figure out why this happens, and any help is appreciated.
I am not able to upload my code because I am a new user.

Thanks in advance,
Jesper

2

Hey Jesper and welcome to the forum!

Could you please provide us with the code? Otherwise it’s hard to reproduce the error. You can just paste everything here (for nice formatting, put it between backticks `).
Alternatively, you could also create a Github Gist.

Btw, did you check that the init_c/init_state parameters are lower case? That was also the issue in another recent discussion. We’re aware of it and try to improve the interface here.

3

Hi Johannes,

My code is lower case so that does not solve it unfortunately.
For the run with sequences of SMB runs, the code is:

	#Import packages and lead way to cadet cli
	import numpy as np
	import matplotlib.pyplot as plt
	import pandas as pd
	import timeit
	from scipy.optimize import  least_squares


	from cadet import Cadet
	Cadet.cadet_path = r'C:\Users\jespfra\Anaconda3\bin\cadet-cli' #working pc

	#%% Functions


	# A function to set the discretization for the units instead of writting in manully for every single unit
	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'}
		
		
		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.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
				
	# General model options
	def SMBmodel(Q1,Q2,Q3,Q4,QD,QE,QF,QR,ts,n_col,j):
		#Setting up the model
		smb_model = Cadet()
		
		
		#Speciy number of unit operations: inputs, CSTR, column and output
		smb_model.root.input.model.nunits = 8
		
		# number of components
		n_comp  = 2
		
		#First unit operation: inlet
		## Feed
		smb_model.root.input.model.unit_000.unit_type = 'INLET'
		smb_model.root.input.model.unit_000.ncomp = n_comp
		smb_model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
		
		
		## Eluent
		smb_model.root.input.model.unit_001.unit_type = 'INLET'
		smb_model.root.input.model.unit_001.ncomp = n_comp
		smb_model.root.input.model.unit_001.inlet_type = 'PIECEWISE_CUBIC_POLY'
		
		## Extract 
		smb_model.root.input.model.unit_002.ncomp = n_comp
		smb_model.root.input.model.unit_002.unit_type = 'OUTLET'
		
		## Raffinate
		smb_model.root.input.model.unit_003.ncomp = n_comp
		smb_model.root.input.model.unit_003.unit_type = 'OUTLET'
		
		
		## Columns
		smb_model.root.input.model.unit_004.unit_type = 'LUMPED_RATE_MODEL_WITHOUT_PORES'
		smb_model.root.input.model.unit_004.ncomp = n_comp
		
		smb_model.root.input.model.unit_004.col_dispersion = 3.8148E-20
		smb_model.root.input.model.unit_004.col_length = 0.25
		smb_model.root.input.model.unit_004.total_porosity = 0.83
		smb_model.root.input.model.unit_004.velocity = 1
		smb_model.root.input.model.unit_004.cross_section_area = 3.141592653589793E-4
		
		smb_model.root.input.model.unit_004.adsorption_model = 'MULTI_COMPONENT_LANGMUIR'
		smb_model.root.input.model.unit_004.adsorption.is_kinetic = 0
		smb_model.root.input.model.unit_004.adsorption.mcl_ka = [ 0.11, 0.148]
		smb_model.root.input.model.unit_004.adsorption.mcl_kd = [ 1.0, 1.0]
		smb_model.root.input.model.unit_004.adsorption.mcl_qmax = [ 5.72/0.11, 7.7/0.148]
		
		smb_model.root.input.model.unit_004.init_c = [0,0.0]
		smb_model.root.input.model.unit_004.init_q = [0.0,0.0]

		### Copy column models
		smb_model.root.input.model.unit_005 = smb_model.root.input.model.unit_004
		smb_model.root.input.model.unit_006 = smb_model.root.input.model.unit_004
		smb_model.root.input.model.unit_007 = smb_model.root.input.model.unit_004
		
		#To write out last output to check for steady state
		smb_model.root.input['return'].write_solution_last = True
		
		
		#% Input and connections
		n_cycles = 5
		switch_time = ts #s
		
		#Sections
		smb_model.root.input.solver.sections.nsec = 4*n_cycles
		smb_model.root.input.solver.sections.section_times = [0]
		for i in range(n_cycles):
			smb_model.root.input.solver.sections.section_times.append((4*i+1)*switch_time)
			smb_model.root.input.solver.sections.section_times.append((4*i+2)*switch_time)
			smb_model.root.input.solver.sections.section_times.append((4*i+3)*switch_time)
			smb_model.root.input.solver.sections.section_times.append((4*i+4)*switch_time)    
		
		## Feed and Eluent concentration
		smb_model.root.input.model.unit_000.sec_000.const_coeff = [CfA,CfB] #Inlet flowrate concentration
		smb_model.root.input.model.unit_001.sec_000.const_coeff = [ 0.0, 0.0]
		
		
		#Connections
		smb_model.root.input.model.connections.nswitches = 4
		
		smb_model.root.input.model.connections.switch_000.section = 0
		smb_model.root.input.model.connections.switch_000.connections = [
			4, 5, -1, -1, Q4,#flowrates, Q, m3/s
			5, 6, -1, -1, Q4,
			6, 7, -1, -1, Q2,
			7, 4, -1, -1, Q2,
			0, 4, -1, -1, QF,
			1, 6, -1, -1, QD,
			4, 2, -1, -1, QR,
			6, 3, -1, -1, QE
		]
		
		smb_model.root.input.model.connections.switch_001.section = 1
		smb_model.root.input.model.connections.switch_001.connections = [
			4,	5,	-1,	-1,	Q2,#flowrates, Q, m3/s
			5,	6,	-1,	-1,	Q4,
			6,	7,	-1,	-1,	Q4,
			7,	4,	-1,	-1,	Q2,
			0,	5,	-1,	-1,	QF,
			1,	7,	-1,	-1,	QD,
			5,	2,	-1,	-1,	QR,
			7,	3,	-1,	-1,	QE
		]
		
		smb_model.root.input.model.connections.switch_002.section = 2
		smb_model.root.input.model.connections.switch_002.connections = [
			4,	5,	-1,	-1,	Q2,#flowrates, Q, m3/s
			5,	6,	-1,	-1,	Q2,
			6,	7,	-1,	-1,	Q4,
			7,	4,	-1,	-1,	Q4,
			0,	6,	-1,	-1,	QF,
			1,	4,	-1,	-1,	QD,
			6,	2,	-1,	-1,	QR,
			4,	3,	-1,	-1,	QE
		]
		
		smb_model.root.input.model.connections.switch_003.section = 3
		smb_model.root.input.model.connections.switch_003.connections = [
			4,	5,	-1,	-1,	Q4,#flowrates, Q, m3/s
			5,	6,	-1,	-1,	Q2,
			6,	7,	-1,	-1,	Q2,
			7,	4,	-1,	-1,	Q4,
			0,	7,	-1,	-1,	QF,
			1,	5,	-1,	-1,	QD,
			7,	2,	-1,	-1,	QR,
			5,	3,	-1,	-1,	QE
		]
		
		#solution times
		smb_model.root.input.solver.user_solution_times = np.linspace(0, n_cycles*4*switch_time, int(n_cycles*4*switch_time))
		
		#% Solver options: spatial and time
		
		#Spatial for the units
		set_discretization(smb_model, n_col=n_col, n_bound = list(np.ones(n_comp)))
	 
		
		
		#Time 
		# Tolerances for the time integrator
		smb_model.root.input.solver.time_integrator.abstol = 1e-6 #absolute tolerance
		smb_model.root.input.solver.time_integrator.algtol = 1e-10
		smb_model.root.input.solver.time_integrator.reltol = 1e-6 #Relative tolerance
		smb_model.root.input.solver.time_integrator.init_step_size = 1e-6
		smb_model.root.input.solver.time_integrator.max_steps = 1000000
		
		
		
		#Solver options in general (not only for column although the same)
		smb_model.root.input.model.solver.gs_type = 1
		smb_model.root.input.model.solver.max_krylov = 0
		smb_model.root.input.model.solver.max_restarts = 10
		smb_model.root.input.model.solver.schur_safety = 1e-8
		
		# Number of cores for parallel simulation
		smb_model.root.input.solver.nthreads = 1
		
		
		#Specify which results we want to return
		# Return data
		smb_model.root.input['return'].split_components_data = 0
		smb_model.root.input['return'].split_ports_data = 0
		smb_model.root.input['return'].unit_000.write_solution_bulk = 1
		smb_model.root.input['return'].unit_000.write_solution_inlet = 1
		smb_model.root.input['return'].unit_000.write_solution_outlet = 1
		
		
		# Copy settings to the other unit operations
		smb_model.root.input['return'].unit_001 = smb_model.root.input['return'].unit_000
		smb_model.root.input['return'].unit_002 = smb_model.root.input['return'].unit_000
		smb_model.root.input['return'].unit_003 = smb_model.root.input['return'].unit_000
		smb_model.root.input['return'].unit_004 = smb_model.root.input['return'].unit_000
		smb_model.root.input['return'].unit_005 = smb_model.root.input['return'].unit_000
		smb_model.root.input['return'].unit_006 = smb_model.root.input['return'].unit_000
		smb_model.root.input['return'].unit_007 = smb_model.root.input['return'].unit_000
		
		
		#Saving data
		smb_model.filename = 'Langmuir_sim_gradient_multi_'+str(j)+'.h5'
		smb_model.save()
		
		#run model
		data = smb_model.run()
		
		if data.returncode == 0:
			print("Simulation completed successfully")
			smb_model.load()   
		else:
			print(data)
			raise Exception("Simulation failed")
		
		df = pd.DataFrame({'Time': smb_model.root.output.solution.solution_times,
						   'Ca_R': smb_model.root.output.solution.unit_002.solution_outlet[:,0],
						   'Cb_R': smb_model.root.output.solution.unit_002.solution_outlet[:,1],
						   'Ca_E': smb_model.root.output.solution.unit_003.solution_outlet[:,0],
						   'Cb_E': smb_model.root.output.solution.unit_003.solution_outlet[:,1] })
		
		
		#If cyclic steady state has not been reached
		state_y = smb_model.root.output.last_state_y
		state_ydot = smb_model.root.output.last_state_ydot

		
		#Try to run SMB_model_SS which is just a new simulation that stores the name differently
		for k in range(0,3):
			df, state_y, state_ydot = SMBmodel_ss(Q1, Q2, Q3, Q4, QD, QE, QF, QR, ts, n_col, j,state_y, state_ydot,k)
			df, state_y, state_ydot = df, state_y, state_ydot
		
		
		return df



	def SMBmodel_ss(Q1,Q2,Q3,Q4,QD,QE,QF,QR,ts,n_col,j,state_y, state_ydot,k):
		#Setting up the model
		smb_model = Cadet()
		
		
		#Speciy number of unit operations: inputs, CSTR, column and output
		smb_model.root.input.model.nunits = 8
		
		# number of components
		n_comp  = 2
		
		#First unit operation: inlet
		## Feed
		smb_model.root.input.model.unit_000.unit_type = 'INLET'
		smb_model.root.input.model.unit_000.ncomp = n_comp
		smb_model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
		
		
		## Eluent
		smb_model.root.input.model.unit_001.unit_type = 'INLET'
		smb_model.root.input.model.unit_001.ncomp = n_comp
		smb_model.root.input.model.unit_001.inlet_type = 'PIECEWISE_CUBIC_POLY'
		
		## Extract 
		smb_model.root.input.model.unit_002.ncomp = n_comp
		smb_model.root.input.model.unit_002.unit_type = 'OUTLET'
		
		## Raffinate
		smb_model.root.input.model.unit_003.ncomp = n_comp
		smb_model.root.input.model.unit_003.unit_type = 'OUTLET'
		
		
		## Columns
		
		smb_model.root.input.model.unit_004.unit_type = 'LUMPED_RATE_MODEL_WITHOUT_PORES'
		smb_model.root.input.model.unit_004.ncomp = n_comp
		
		smb_model.root.input.model.unit_004.col_dispersion = 3.8148E-20
		smb_model.root.input.model.unit_004.col_length = 0.25
		smb_model.root.input.model.unit_004.total_porosity = 0.83
		smb_model.root.input.model.unit_004.velocity = 1
		smb_model.root.input.model.unit_004.cross_section_area = 3.141592653589793E-4
		
		smb_model.root.input.model.unit_004.adsorption_model = 'MULTI_COMPONENT_LANGMUIR'
		smb_model.root.input.model.unit_004.adsorption.is_kinetic = 0
		smb_model.root.input.model.unit_004.adsorption.mcl_ka = [0.11, 0.148]
		smb_model.root.input.model.unit_004.adsorption.mcl_kd = [ 1.0, 1.0]
		smb_model.root.input.model.unit_004.adsorption.mcl_qmax = [ 5.72/0.11, 7.7/0.148]
		
		# smb_model.root.input.model.unit_004.init_c = [0.0,0.0]
		# smb_model.root.input.model.unit_004.init_q = [0.0,0.0]

		### Copy column models
		smb_model.root.input.model.unit_005 = smb_model.root.input.model.unit_004
		smb_model.root.input.model.unit_006 = smb_model.root.input.model.unit_004
		smb_model.root.input.model.unit_007 = smb_model.root.input.model.unit_004
		
		#To write out last output to check for steady state
		smb_model.root.input['return'].write_solution_last = True
		
		#Load initial conditions from last state
		smb_model.root.input.model.init_state_y = state_y
		smb_model.root.input.model.init_state_ydot = state_ydot
		
		
		# smb_model.root.input.model.unit_004.init_c = [200, 0.0, 0.0, 0] #HERE
		# smb_model.root.input.model.unit_004.init_q = [Lambda_, 0.0,0.0, 0] #HERE


		### Copy column models
		smb_model.root.input.model.unit_005 = smb_model.root.input.model.unit_004
		smb_model.root.input.model.unit_006 = smb_model.root.input.model.unit_004
		smb_model.root.input.model.unit_007 = smb_model.root.input.model.unit_004
		
		
		
		#% Input and connections
		n_cycles = 3
		switch_time = ts #s
		
		#Sections
		smb_model.root.input.solver.sections.nsec = 4*n_cycles
		smb_model.root.input.solver.sections.section_times = [0]
		for i in range(n_cycles):
			smb_model.root.input.solver.sections.section_times.append((4*i+1)*switch_time)
			smb_model.root.input.solver.sections.section_times.append((4*i+2)*switch_time)
			smb_model.root.input.solver.sections.section_times.append((4*i+3)*switch_time)
			smb_model.root.input.solver.sections.section_times.append((4*i+4)*switch_time)    
		
		## Feed and Eluent concentration
		smb_model.root.input.model.unit_000.sec_000.const_coeff = [CfA,CfB] #Inlet flowrate concentration
		smb_model.root.input.model.unit_001.sec_000.const_coeff = [0.0, 0.0]
		
		
		#Connections
		smb_model.root.input.model.connections.nswitches = 4
		
		smb_model.root.input.model.connections.switch_000.section = 0
		smb_model.root.input.model.connections.switch_000.connections = [
			4, 5, -1, -1, Q4,#flowrates, Q, m3/s
			5, 6, -1, -1, Q4,
			6, 7, -1, -1, Q2,
			7, 4, -1, -1, Q2,
			0, 4, -1, -1, QF,
			1, 6, -1, -1, QD,
			4, 2, -1, -1, QR,
			6, 3, -1, -1, QE
		]
		
		smb_model.root.input.model.connections.switch_001.section = 1
		smb_model.root.input.model.connections.switch_001.connections = [
			4,	5,	-1,	-1,	Q2,#flowrates, Q, m3/s
			5,	6,	-1,	-1,	Q4,
			6,	7,	-1,	-1,	Q4,
			7,	4,	-1,	-1,	Q2,
			0,	5,	-1,	-1,	QF,
			1,	7,	-1,	-1,	QD,
			5,	2,	-1,	-1,	QR,
			7,	3,	-1,	-1,	QE
		]
		
		smb_model.root.input.model.connections.switch_002.section = 2
		smb_model.root.input.model.connections.switch_002.connections = [
			4,	5,	-1,	-1,	Q2,#flowrates, Q, m3/s
			5,	6,	-1,	-1,	Q2,
			6,	7,	-1,	-1,	Q4,
			7,	4,	-1,	-1,	Q4,
			0,	6,	-1,	-1,	QF,
			1,	4,	-1,	-1,	QD,
			6,	2,	-1,	-1,	QR,
			4,	3,	-1,	-1,	QE
		]
		
		smb_model.root.input.model.connections.switch_003.section = 3
		smb_model.root.input.model.connections.switch_003.connections = [
			4,	5,	-1,	-1,	Q4,#flowrates, Q, m3/s
			5,	6,	-1,	-1,	Q2,
			6,	7,	-1,	-1,	Q2,
			7,	4,	-1,	-1,	Q4,
			0,	7,	-1,	-1,	QF,
			1,	5,	-1,	-1,	QD,
			7,	2,	-1,	-1,	QR,
			5,	3,	-1,	-1,	QE
		]
		
		#solution times
		smb_model.root.input.solver.user_solution_times = np.linspace(0, n_cycles*4*switch_time, int(n_cycles*4*switch_time))
		
		#% Solver options: spatial and time
		
		#Spatial for the units
		set_discretization(smb_model, n_col=n_col, n_bound = list(np.ones(n_comp)))
	 
		
		
		#Time 
		# Tolerances for the time integrator
		smb_model.root.input.solver.time_integrator.abstol = 1e-6 #absolute tolerance
		smb_model.root.input.solver.time_integrator.algtol = 1e-10
		smb_model.root.input.solver.time_integrator.reltol = 1e-6 #Relative tolerance
		smb_model.root.input.solver.time_integrator.init_step_size = 1e-6
		smb_model.root.input.solver.time_integrator.max_steps = 1000000
		
		
		
		#Solver options in general (not only for column although the same)
		smb_model.root.input.model.solver.gs_type = 1
		smb_model.root.input.model.solver.max_krylov = 0
		smb_model.root.input.model.solver.max_restarts = 10
		smb_model.root.input.model.solver.schur_safety = 1e-8
		
		# Number of cores for parallel simulation
		smb_model.root.input.solver.nthreads = 1
		
		
		#Specify which results we want to return
		# Return data
		smb_model.root.input['return'].split_components_data = 0
		smb_model.root.input['return'].split_ports_data = 0
		smb_model.root.input['return'].unit_000.write_solution_bulk = 1
		smb_model.root.input['return'].unit_000.write_solution_inlet = 1
		smb_model.root.input['return'].unit_000.write_solution_outlet = 1
		
		
		# Copy settings to the other unit operations
		smb_model.root.input['return'].unit_001 = smb_model.root.input['return'].unit_000
		smb_model.root.input['return'].unit_002 = smb_model.root.input['return'].unit_000
		smb_model.root.input['return'].unit_003 = smb_model.root.input['return'].unit_000
		smb_model.root.input['return'].unit_004 = smb_model.root.input['return'].unit_000
		smb_model.root.input['return'].unit_005 = smb_model.root.input['return'].unit_000
		smb_model.root.input['return'].unit_006 = smb_model.root.input['return'].unit_000
		smb_model.root.input['return'].unit_007 = smb_model.root.input['return'].unit_000
		
		
		#Saving data
		smb_model.filename = 'Langmuir_sim_gradient_multi_S_'+str(j)+'_'+str(k)+'.h5'
		smb_model.save()
		
		#run model
		data = smb_model.run()
		
		if data.returncode == 0:
			print("Simulation completed successfully")
			smb_model.load()   
		else:
			print(data)
			raise Exception("Simulation failed")
		
		df = pd.DataFrame({'Time': smb_model.root.output.solution.solution_times,
						   'Ca_R': smb_model.root.output.solution.unit_002.solution_outlet[:,0],
						   'Cb_R': smb_model.root.output.solution.unit_002.solution_outlet[:,1],
						   'Ca_E': smb_model.root.output.solution.unit_003.solution_outlet[:,0],
						   'Cb_E': smb_model.root.output.solution.unit_003.solution_outlet[:,1] })
		
		state_y = smb_model.root.output.last_state_y
		state_ydot = smb_model.root.output.last_state_ydot
		
		return df, state_y, state_ydot



	#Determining isotherm and function for determining the region of separation

	#%% binding parameters
	#We set the binding parameters
	Ka = 0.148 #l/g
	Kb = 0.11
	gamma_a = 7.7 
	gamma_b = 5.72
	Na = gamma_a/Ka #g/L max capacity
	Nb = gamma_b/Kb



	#Inlet concentrations
	CfA = 0.55   #g/L
	CfB = 0.55


	eps = 0.83
	VC = (0.25*0.02**2/4)*np.pi #m3
	ts = 180 #s

	data = []

	n_col = 80


	m2_c = 6
	m3_c = 9
	m1_min = gamma_a
	m4_max = 1/2*(gamma_b + m3_c + Kb*CfB*(m3_c-m2_c)-((gamma_b + m3_c + Kb*CfB*(m3_c-m2_c))**2-4*gamma_b*m3_c)**0.5)
	m1_c = 8.2
	m4_c = m4_max*0.5

	Q1 = -(m1_c*eps - m1_c - eps)*VC/ts
	Q2 = -(m2_c*eps - m2_c - eps)*VC/ts
	Q3 = -(m3_c*eps - m3_c - eps)*VC/ts
	Q4 = -(m4_c*eps - m4_c - eps)*VC/ts
	QD = -VC*(m1_c*eps - m4_c*eps - m1_c + m4_c)/ts
	QE = -VC*(m1_c*eps - m2_c*eps - m1_c + m2_c)/ts
	QF = VC*(m2_c*eps - m3_c*eps - m2_c + m3_c)/ts
	QR = -VC*(m3_c*eps - m4_c*eps - m3_c + m4_c)/ts

	j=0

	#Running the SMB model 
	df = SMBmodel(Q1,Q2,Q3,Q4,QD,QE,QF,QR,ts,n_col,j)

	#%% Load data and plot the issue
	ts = 180
	data = []
		
	filename = 'Langmuir_sim_gradient_multi_0.h5'
	smb_model = Cadet()
	smb_model.filename = filename
	smb_model.load()
	df = pd.DataFrame({'Time': smb_model.root.output.solution.solution_times,
					   'Ca_E': smb_model.root.output.solution.unit_002.solution_outlet[:,0],
					   'Cb_E': smb_model.root.output.solution.unit_002.solution_outlet[:,1],
					   'Ca_R': smb_model.root.output.solution.unit_003.solution_outlet[:,0],
					   'Cb_R': smb_model.root.output.solution.unit_003.solution_outlet[:,1] })
	data.append(df)
		


	filename = 'Langmuir_sim_gradient_multi_S_0_0.h5'
	smb_model = Cadet()
	smb_model.filename = filename
	smb_model.load()
	df = pd.DataFrame({'Time': smb_model.root.output.solution.solution_times,
					   'Ca_E': smb_model.root.output.solution.unit_002.solution_outlet[:,0],
					   'Cb_E': smb_model.root.output.solution.unit_002.solution_outlet[:,1],
					   'Ca_R': smb_model.root.output.solution.unit_003.solution_outlet[:,0],
					   'Cb_R': smb_model.root.output.solution.unit_003.solution_outlet[:,1] })
	data.append(df)

	filename = 'Langmuir_sim_gradient_multi_S_0_1.h5'
	smb_model = Cadet()
	smb_model.filename = filename
	smb_model.load()
	df = pd.DataFrame({'Time': smb_model.root.output.solution.solution_times,
					   'Ca_E': smb_model.root.output.solution.unit_002.solution_outlet[:,0],
					   'Cb_E': smb_model.root.output.solution.unit_002.solution_outlet[:,1],
					   'Ca_R': smb_model.root.output.solution.unit_003.solution_outlet[:,0],
					   'Cb_R': smb_model.root.output.solution.unit_003.solution_outlet[:,1] })
	data.append(df)

	filename = 'Langmuir_sim_gradient_multi_S_0_2.h5'
	smb_model = Cadet()
	smb_model.filename = filename
	smb_model.load()
	df = pd.DataFrame({'Time': smb_model.root.output.solution.solution_times,
					   'Ca_E': smb_model.root.output.solution.unit_002.solution_outlet[:,0],
					   'Cb_E': smb_model.root.output.solution.unit_002.solution_outlet[:,1],
					   'Ca_R': smb_model.root.output.solution.unit_003.solution_outlet[:,0],
					   'Cb_R': smb_model.root.output.solution.unit_003.solution_outlet[:,1] })
	data.append(df)


	#plotting the problem
	plt.plot(np.concatenate((np.array(data[0]['Ca_R'][:]),np.array(data[1]['Ca_R'][:]),np.array(data[2]['Ca_R'][:]),np.array(data[3]['Ca_R'][:]))),label='Ca in Raffinate')
	plt.plot([5*4*ts,5*4*ts],[0,0.005],'--',label='First run')
	plt.plot([5*4*ts+3*4*ts,5*4*ts+3*4*ts],[0,0.005],'--',label='second run')
	plt.plot([5*4*ts+3*4*ts,5*4*ts+3*4*ts],[0,0.005],'--',label='third run')
	plt.plot([5*4*ts+3*4*ts*2,5*4*ts+3*4*ts*2],[0,0.005],'--',label='fourth run')
	plt.legend()
	plt.xlabel('Time(s)')
	plt.ylabel('Concentration')
4

For the non-sequential simulation, the code is:

	#Import packages and lead way to cadet cli
	import numpy as np
	import matplotlib.pyplot as plt
	import pandas as pd
	from scipy.optimize import  least_squares


	from cadet import Cadet
	Cadet.cadet_path = r'C:\Users\jespfra\Anaconda3\bin\cadet-cli' #working pc

	#%% Functions


	# A function to set the discretization for the units instead of writting in manully for every single unit
	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'}
		
		
		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.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
				
	# General model options
	def SMBmodel(Q1,Q2,Q3,Q4,QD,QE,QF,QR,ts,n_col,j):
		#Setting up the model
		smb_model = Cadet()
		
		
		#Speciy number of unit operations: inputs, CSTR, column and output
		smb_model.root.input.model.nunits = 8
		
		# number of components
		n_comp  = 2
		
		#First unit operation: inlet
		## Feed
		smb_model.root.input.model.unit_000.unit_type = 'INLET'
		smb_model.root.input.model.unit_000.ncomp = n_comp
		smb_model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
		
		
		## Eluent
		smb_model.root.input.model.unit_001.unit_type = 'INLET'
		smb_model.root.input.model.unit_001.ncomp = n_comp
		smb_model.root.input.model.unit_001.inlet_type = 'PIECEWISE_CUBIC_POLY'
		
		## Extract 
		smb_model.root.input.model.unit_002.ncomp = n_comp
		smb_model.root.input.model.unit_002.unit_type = 'OUTLET'
		
		## Raffinate
		smb_model.root.input.model.unit_003.ncomp = n_comp
		smb_model.root.input.model.unit_003.unit_type = 'OUTLET'
		
		
		## Columns
		smb_model.root.input.model.unit_004.unit_type = 'LUMPED_RATE_MODEL_WITHOUT_PORES'
		smb_model.root.input.model.unit_004.ncomp = n_comp
		
		smb_model.root.input.model.unit_004.col_dispersion = 3.8148E-20
		smb_model.root.input.model.unit_004.col_length = 0.25
		smb_model.root.input.model.unit_004.total_porosity = 0.83
		smb_model.root.input.model.unit_004.velocity = 1
		smb_model.root.input.model.unit_004.cross_section_area = 3.141592653589793E-4
		
		smb_model.root.input.model.unit_004.adsorption_model = 'MULTI_COMPONENT_LANGMUIR'
		smb_model.root.input.model.unit_004.adsorption.is_kinetic = 0
		smb_model.root.input.model.unit_004.adsorption.mcl_ka = [ 0.11, 0.148]
		smb_model.root.input.model.unit_004.adsorption.mcl_kd = [ 1.0, 1.0]
		smb_model.root.input.model.unit_004.adsorption.mcl_qmax = [ 5.72/0.11, 7.7/0.148]
		
		smb_model.root.input.model.unit_004.init_c = [0,0.0]
		smb_model.root.input.model.unit_004.init_q = [0.0,0.0]

		### Copy column models
		smb_model.root.input.model.unit_005 = smb_model.root.input.model.unit_004
		smb_model.root.input.model.unit_006 = smb_model.root.input.model.unit_004
		smb_model.root.input.model.unit_007 = smb_model.root.input.model.unit_004
		
		#To write out last output to check for steady state
		smb_model.root.input['return'].WRITE_SOLUTION_LAST = True
		
		
		#% Input and connections
		n_cycles = 14
		switch_time = ts #s
		
		#Sections
		smb_model.root.input.solver.sections.nsec = 4*n_cycles
		smb_model.root.input.solver.sections.section_times = [0]
		for i in range(n_cycles):
			smb_model.root.input.solver.sections.section_times.append((4*i+1)*switch_time)
			smb_model.root.input.solver.sections.section_times.append((4*i+2)*switch_time)
			smb_model.root.input.solver.sections.section_times.append((4*i+3)*switch_time)
			smb_model.root.input.solver.sections.section_times.append((4*i+4)*switch_time)    
		
		## Feed and Eluent concentration
		smb_model.root.input.model.unit_000.sec_000.const_coeff = [CfA,CfB] #Inlet flowrate concentration
		smb_model.root.input.model.unit_001.sec_000.const_coeff = [ 0.0, 0.0]
		
		
		#Connections
		smb_model.root.input.model.connections.nswitches = 4
		
		smb_model.root.input.model.connections.switch_000.section = 0
		smb_model.root.input.model.connections.switch_000.connections = [
			4, 5, -1, -1, Q4,#flowrates, Q, m3/s
			5, 6, -1, -1, Q4,
			6, 7, -1, -1, Q2,
			7, 4, -1, -1, Q2,
			0, 4, -1, -1, QF,
			1, 6, -1, -1, QD,
			4, 2, -1, -1, QR,
			6, 3, -1, -1, QE
		]
		
		smb_model.root.input.model.connections.switch_001.section = 1
		smb_model.root.input.model.connections.switch_001.connections = [
			4,	5,	-1,	-1,	Q2,#flowrates, Q, m3/s
			5,	6,	-1,	-1,	Q4,
			6,	7,	-1,	-1,	Q4,
			7,	4,	-1,	-1,	Q2,
			0,	5,	-1,	-1,	QF,
			1,	7,	-1,	-1,	QD,
			5,	2,	-1,	-1,	QR,
			7,	3,	-1,	-1,	QE
		]
		
		smb_model.root.input.model.connections.switch_002.section = 2
		smb_model.root.input.model.connections.switch_002.connections = [
			4,	5,	-1,	-1,	Q2,#flowrates, Q, m3/s
			5,	6,	-1,	-1,	Q2,
			6,	7,	-1,	-1,	Q4,
			7,	4,	-1,	-1,	Q4,
			0,	6,	-1,	-1,	QF,
			1,	4,	-1,	-1,	QD,
			6,	2,	-1,	-1,	QR,
			4,	3,	-1,	-1,	QE
		]
		
		smb_model.root.input.model.connections.switch_003.section = 3
		smb_model.root.input.model.connections.switch_003.connections = [
			4,	5,	-1,	-1,	Q4,#flowrates, Q, m3/s
			5,	6,	-1,	-1,	Q2,
			6,	7,	-1,	-1,	Q2,
			7,	4,	-1,	-1,	Q4,
			0,	7,	-1,	-1,	QF,
			1,	5,	-1,	-1,	QD,
			7,	2,	-1,	-1,	QR,
			5,	3,	-1,	-1,	QE
		]
		
		#solution times
		smb_model.root.input.solver.user_solution_times = np.linspace(0, n_cycles*4*switch_time, int(n_cycles*4*switch_time))
		
		#% Solver options: spatial and time
		
		#Spatial for the units
		set_discretization(smb_model, n_col=n_col, n_bound = list(np.ones(n_comp)))
	 
		
		
		#Time 
		# Tolerances for the time integrator
		smb_model.root.input.solver.time_integrator.abstol = 1e-6 #absolute tolerance
		smb_model.root.input.solver.time_integrator.algtol = 1e-10
		smb_model.root.input.solver.time_integrator.reltol = 1e-6 #Relative tolerance
		smb_model.root.input.solver.time_integrator.init_step_size = 1e-6
		smb_model.root.input.solver.time_integrator.max_steps = 1000000
		
		
		
		#Solver options in general (not only for column although the same)
		smb_model.root.input.model.solver.gs_type = 1
		smb_model.root.input.model.solver.max_krylov = 0
		smb_model.root.input.model.solver.max_restarts = 10
		smb_model.root.input.model.solver.schur_safety = 1e-8
		
		# Number of cores for parallel simulation
		smb_model.root.input.solver.nthreads = 1
		
		
		#Specify which results we want to return
		# Return data
		smb_model.root.input['return'].split_components_data = 0
		smb_model.root.input['return'].split_ports_data = 0
		smb_model.root.input['return'].unit_000.write_solution_bulk = 1
		smb_model.root.input['return'].unit_000.write_solution_inlet = 1
		smb_model.root.input['return'].unit_000.write_solution_outlet = 1
		
		
		# Copy settings to the other unit operations
		smb_model.root.input['return'].unit_001 = smb_model.root.input['return'].unit_000
		smb_model.root.input['return'].unit_002 = smb_model.root.input['return'].unit_000
		smb_model.root.input['return'].unit_003 = smb_model.root.input['return'].unit_000
		smb_model.root.input['return'].unit_004 = smb_model.root.input['return'].unit_000
		smb_model.root.input['return'].unit_005 = smb_model.root.input['return'].unit_000
		smb_model.root.input['return'].unit_006 = smb_model.root.input['return'].unit_000
		smb_model.root.input['return'].unit_007 = smb_model.root.input['return'].unit_000
		
		
		#Saving data
		smb_model.filename = 'Langmuir_sim_gradient_multi2_'+str(j)+'.h5'
		smb_model.save()
		
		#run model
		data = smb_model.run()
		
		if data.returncode == 0:
			print("Simulation completed successfully")
			smb_model.load()   
		else:
			print(data)
			raise Exception("Simulation failed")
		
		df = pd.DataFrame({'Time': smb_model.root.output.solution.solution_times,
						   'Ca_E': smb_model.root.output.solution.unit_002.solution_outlet[:,0],
						   'Cb_E': smb_model.root.output.solution.unit_002.solution_outlet[:,1],
						   'Ca_R': smb_model.root.output.solution.unit_003.solution_outlet[:,0],
						   'Cb_R': smb_model.root.output.solution.unit_003.solution_outlet[:,1] })
		   
		
		return df




	#%% binding parameters
	Ka = 0.148 #l/g
	Kb = 0.11
	gamma_a = 7.7 
	gamma_b = 5.72
	Na = gamma_a/Ka #g/L max capacity
	Nb = gamma_b/Kb



	#Inlet concentrations
	CfA = 0.55   #g/L
	CfB = 0.55


	eps = 0.83
	VC = (0.25*0.02**2/4)*np.pi #m3
	ts = 180 #s

	data = []

	n_col = 80


	m2_c = 6
	m3_c = 9
	m1_min = gamma_a
	m4_max = 1/2*(gamma_b + m3_c + Kb*CfB*(m3_c-m2_c)-((gamma_b + m3_c + Kb*CfB*(m3_c-m2_c))**2-4*gamma_b*m3_c)**0.5)
	m1_c = 8.2
	m4_c = m4_max*0.5

	Q1 = -(m1_c*eps - m1_c - eps)*VC/ts
	Q2 = -(m2_c*eps - m2_c - eps)*VC/ts
	Q3 = -(m3_c*eps - m3_c - eps)*VC/ts
	Q4 = -(m4_c*eps - m4_c - eps)*VC/ts
	QD = -VC*(m1_c*eps - m4_c*eps - m1_c + m4_c)/ts
	QE = -VC*(m1_c*eps - m2_c*eps - m1_c + m2_c)/ts
	QF = VC*(m2_c*eps - m3_c*eps - m2_c + m3_c)/ts
	QR = -VC*(m3_c*eps - m4_c*eps - m3_c + m4_c)/ts

	j=0

	#Running the SMB model 
	df = SMBmodel(Q1,Q2,Q3,Q4,QD,QE,QF,QR,ts,n_col,j)


	plt.plot(df['Ca_R'],label='Ca in raffinate')
	plt.plot([5*4*ts,5*4*ts],[0,0.005],'--',label='First run')
	plt.plot([5*4*ts+3*4*ts,5*4*ts+3*4*ts],[0,0.005],'--',label='second run')
	plt.plot([5*4*ts+3*4*ts,5*4*ts+3*4*ts],[0,0.005],'--',label='third run')
	plt.plot([5*4*ts+3*4*ts*2,5*4*ts+3*4*ts*2],[0,0.005],'--',label='fourth run')
	plt.legend()
	plt.xlabel('Time(s)')
	plt.ylabel('Concentration')
5

Please try disabling consistent initialization for all but the first simulation:

smb_model.root.input.solver.consistent_init_mode = 6 # or 7 or 0

The first simulation needs to have consistent initialization enabled!

6

When I implement the consistent_init_mode =0 or 6 or 7 for all but the first simulation, the error I get from the second simulation is:

CompletedProcess(args=[‘C:\Users\jespfra\Anaconda3\bin\cadet-cli’, ‘Langmuir_sim_gradient_multi_S_0_0.h5’], returncode=3, stdout=b"[Error: idasErrorHandler::200] In function ‘IDASolve’ of module ‘IDAS’, error code ‘IDA_ERR_FAIL’:\r\nAt t = 0 and h = 2.44141e-10, the error test failed repeatedly or with |h| = hmin.\r\n[Error: integrate::1367] IDASolve returned IDA_ERR_FAIL at t = 0\r\n", stderr=b’SOLVER ERROR: Error in IDASolve: IDA_ERR_FAIL at t = 0.000000\r\n’)

If I try to set the consistent_init_mode =2 or 4, the same error occurs, just at switching time (IDA_ERR_FAIL at t = 0\180s in this case).

7

Hey Jesper,

I’m having a difficult time trying to understand your code. Maybe we can have a call at some point. Unfortunately, we’re quite busy with preparations for our workshop next week so I hope you can wait a bit longer until I can spend more time on this.

8

Hey Johannes,

That is totally fine, I’ll attend the workshop myself, then maybe we find time to go through it.
It is not nicely coded, I know :smiley:

1 Like
9

Hey Jesper,

I’m curious, did you get everything up and running?

Cheers

Jo

10

Hi Johannes,

It works when using CADET-process but it does not work when using CADET.
A simpler code that demonstrates it is shown below.

#Import packages and lead way to cadet cli
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
# from scipy.optimize import  least_squares


from cadet import Cadet
Cadet.cadet_path = r'C:\Users\jespfra\Anaconda3\bin\cadet-cli' #working pc

#%% Functions


# A function to set the discretization for the units instead of writting in manully for every single unit
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'}
    
    
    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.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
            
# General model options
def SMBmodel(Q1,Q2,Q3,Q4,QD,QE,QF,QR,ts,n_col,j,state_y,state_ydot):
    #Setting up the model
    smb_model = Cadet()
    
    
    #Speciy number of unit operations: inputs, CSTR, column and output
    smb_model.root.input.model.nunits = 8
    
    # number of components
    n_comp  = 2
    
    #First unit operation: inlet
    ## Feed
    smb_model.root.input.model.unit_000.unit_type = 'INLET'
    smb_model.root.input.model.unit_000.ncomp = n_comp
    smb_model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
    
    
    ## Eluent
    smb_model.root.input.model.unit_001.unit_type = 'INLET'
    smb_model.root.input.model.unit_001.ncomp = n_comp
    smb_model.root.input.model.unit_001.inlet_type = 'PIECEWISE_CUBIC_POLY'
    
    ## Extract 
    smb_model.root.input.model.unit_002.ncomp = n_comp
    smb_model.root.input.model.unit_002.unit_type = 'OUTLET'
    
    ## Raffinate
    smb_model.root.input.model.unit_003.ncomp = n_comp
    smb_model.root.input.model.unit_003.unit_type = 'OUTLET'
    
    
    ## Columns
    smb_model.root.input.model.unit_004.unit_type = 'LUMPED_RATE_MODEL_WITHOUT_PORES'
    smb_model.root.input.model.unit_004.ncomp = n_comp
    
    smb_model.root.input.model.unit_004.col_dispersion = 0
    smb_model.root.input.model.unit_004.col_length = 0.25
    smb_model.root.input.model.unit_004.total_porosity = 0.83
    smb_model.root.input.model.unit_004.velocity = 1
    smb_model.root.input.model.unit_004.cross_section_area = 3.141592653589793E-4
    
    smb_model.root.input.model.unit_004.adsorption_model = 'MULTI_COMPONENT_LANGMUIR'
    smb_model.root.input.model.unit_004.adsorption.is_kinetic = 0
    smb_model.root.input.model.unit_004.adsorption.mcl_ka = [ 0.11, 0.148]
    smb_model.root.input.model.unit_004.adsorption.mcl_kd = [ 1.0, 1.0]
    smb_model.root.input.model.unit_004.adsorption.mcl_qmax = [ 5.72/0.11, 7.7/0.148]
    
    if j == 0:
        smb_model.root.input.model.unit_004.init_c = [0,0.0]
        smb_model.root.input.model.unit_004.init_q = [0.0,0.0]
        #Disabling consistent initialization mode - Leweke
        # smb_model.root.input.solver.consistent_init_mode = 1
    else:
        smb_model.root.input.model.init_state_y = state_y
        smb_model.root.input.model.init_state_ydot = state_ydot
        #Disabling consistent initialization mode - Leweke
        # smb_model.root.input.solver.consistent_init_mode = 1

    ### Copy column models
    smb_model.root.input.model.unit_005 = smb_model.root.input.model.unit_004
    smb_model.root.input.model.unit_006 = smb_model.root.input.model.unit_004
    smb_model.root.input.model.unit_007 = smb_model.root.input.model.unit_004
    
    #To write out last output to check for steady state
    smb_model.root.input['return'].write_solution_last = True

    
    #% Input and connections
    n_cycles = 4
    switch_time = ts #s
    
    #Sections
    smb_model.root.input.solver.sections.nsec = 4*n_cycles
    smb_model.root.input.solver.sections.section_times = [0]
    for i in range(n_cycles):
        smb_model.root.input.solver.sections.section_times.append((4*i+1)*switch_time)
        smb_model.root.input.solver.sections.section_times.append((4*i+2)*switch_time)
        smb_model.root.input.solver.sections.section_times.append((4*i+3)*switch_time)
        smb_model.root.input.solver.sections.section_times.append((4*i+4)*switch_time)    
    
    ## Feed and Eluent concentration
    smb_model.root.input.model.unit_000.sec_000.const_coeff = [CfA,CfB] #Inlet flowrate concentration
    smb_model.root.input.model.unit_001.sec_000.const_coeff = [ 0.0, 0.0]
    
    
    #Connections
    smb_model.root.input.model.connections.nswitches = 4
    
    smb_model.root.input.model.connections.switch_000.section = 0
    smb_model.root.input.model.connections.switch_000.connections = [
        4, 5, -1, -1, Q4,#flowrates, Q, m3/s
        5, 6, -1, -1, Q4,
        6, 7, -1, -1, Q2,
        7, 4, -1, -1, Q2,
        0, 4, -1, -1, QF,
        1, 6, -1, -1, QD,
        4, 2, -1, -1, QR,
        6, 3, -1, -1, QE
    ]
    
    smb_model.root.input.model.connections.switch_001.section = 1
    smb_model.root.input.model.connections.switch_001.connections = [
        4,	5,	-1,	-1,	Q2,#flowrates, Q, m3/s
        5,	6,	-1,	-1,	Q4,
        6,	7,	-1,	-1,	Q4,
        7,	4,	-1,	-1,	Q2,
        0,	5,	-1,	-1,	QF,
        1,	7,	-1,	-1,	QD,
        5,	2,	-1,	-1,	QR,
        7,	3,	-1,	-1,	QE
    ]
    
    smb_model.root.input.model.connections.switch_002.section = 2
    smb_model.root.input.model.connections.switch_002.connections = [
        4,	5,	-1,	-1,	Q2,#flowrates, Q, m3/s
        5,	6,	-1,	-1,	Q2,
        6,	7,	-1,	-1,	Q4,
        7,	4,	-1,	-1,	Q4,
        0,	6,	-1,	-1,	QF,
        1,	4,	-1,	-1,	QD,
        6,	2,	-1,	-1,	QR,
        4,	3,	-1,	-1,	QE
    ]
    
    smb_model.root.input.model.connections.switch_003.section = 3
    smb_model.root.input.model.connections.switch_003.connections = [
        4,	5,	-1,	-1,	Q4,#flowrates, Q, m3/s
        5,	6,	-1,	-1,	Q2,
        6,	7,	-1,	-1,	Q2,
        7,	4,	-1,	-1,	Q4,
        0,	7,	-1,	-1,	QF,
        1,	5,	-1,	-1,	QD,
        7,	2,	-1,	-1,	QR,
        5,	3,	-1,	-1,	QE
    ]
    
    #solution times
    smb_model.root.input.solver.user_solution_times = np.linspace(0, n_cycles*4*switch_time, int(n_cycles*4*switch_time))
    
    #% Solver options: spatial and time
    
    #Spatial for the units
    set_discretization(smb_model, n_col=n_col, n_bound = list(np.ones(n_comp)))
 
    
    
    #Time 
    # Tolerances for the time integrator
    smb_model.root.input.solver.time_integrator.abstol = 1e-6 #absolute tolerance
    smb_model.root.input.solver.time_integrator.algtol = 1e-10
    smb_model.root.input.solver.time_integrator.reltol = 1e-6 #Relative tolerance
    smb_model.root.input.solver.time_integrator.init_step_size = 1e-6
    smb_model.root.input.solver.time_integrator.max_steps = 1000000
    
    
    
    #Solver options in general (not only for column although the same)
    smb_model.root.input.model.solver.gs_type = 1
    smb_model.root.input.model.solver.max_krylov = 0
    smb_model.root.input.model.solver.max_restarts = 10
    smb_model.root.input.model.solver.schur_safety = 1e-8
    
    # Number of cores for parallel simulation
    smb_model.root.input.solver.nthreads = 1
    
    
    #Specify which results we want to return
    # Return data
    smb_model.root.input['return'].split_components_data = 0
    smb_model.root.input['return'].split_ports_data = 0
    smb_model.root.input['return'].unit_000.write_solution_bulk = 1
    smb_model.root.input['return'].unit_000.write_solution_inlet = 1
    smb_model.root.input['return'].unit_000.write_solution_outlet = 1
    
    
    # Copy settings to the other unit operations
    smb_model.root.input['return'].unit_001 = smb_model.root.input['return'].unit_000
    smb_model.root.input['return'].unit_002 = smb_model.root.input['return'].unit_000
    smb_model.root.input['return'].unit_003 = smb_model.root.input['return'].unit_000
    smb_model.root.input['return'].unit_004 = smb_model.root.input['return'].unit_000
    smb_model.root.input['return'].unit_005 = smb_model.root.input['return'].unit_000
    smb_model.root.input['return'].unit_006 = smb_model.root.input['return'].unit_000
    smb_model.root.input['return'].unit_007 = smb_model.root.input['return'].unit_000
    
    
    #Saving data
    smb_model.filename = 'Langmuir_sim_gradient_multi_'+str(j)+'.h5'
    smb_model.save()
    
    #run model
    data = smb_model.run()
    
    if data.returncode == 0:
        print("Simulation completed successfully")
        smb_model.load()

    else:
        print(data)
        raise Exception("Simulation failed")
    
    df = pd.DataFrame({'Time': smb_model.root.output.solution.solution_times,
                       'Ca_R': smb_model.root.output.solution.unit_002.solution_outlet[:,0],
                       'Cb_R': smb_model.root.output.solution.unit_002.solution_outlet[:,1],
                       'Ca_E': smb_model.root.output.solution.unit_003.solution_outlet[:,0],
                       'Cb_E': smb_model.root.output.solution.unit_003.solution_outlet[:,1] })
    
    #If cyclic steady state has not been reached
    state_y = smb_model.root.output.last_state_y
    state_ydot = smb_model.root.output.last_state_ydot
    
    return df,state_y,state_ydot


#%% binding parameters
#We set the binding parameters
Ka = 0.148 #l/g
Kb = 0.11
gamma_a = 7.7 
gamma_b = 5.72


#Inlet concentrations
CfA = 0.55   #g/L
CfB = 0.55


eps = 0.83
VC = (0.25*0.02**2/4)*np.pi #m3
ts = 180 #s

n_col = 80

m2_c = 6
m3_c = 9
m1_min = gamma_a
m4_max = 1/2*(gamma_b + m3_c + Kb*CfB*(m3_c-m2_c)-((gamma_b + m3_c + Kb*CfB*(m3_c-m2_c))**2-4*gamma_b*m3_c)**0.5)
m1_c = 8.2
m4_c = m4_max*0.5

Q1 = -(m1_c*eps - m1_c - eps)*VC/ts
Q2 = -(m2_c*eps - m2_c - eps)*VC/ts
Q3 = -(m3_c*eps - m3_c - eps)*VC/ts
Q4 = -(m4_c*eps - m4_c - eps)*VC/ts
QD = -VC*(m1_c*eps - m4_c*eps - m1_c + m4_c)/ts
QE = -VC*(m1_c*eps - m2_c*eps - m1_c + m2_c)/ts
QF = VC*(m2_c*eps - m3_c*eps - m2_c + m3_c)/ts
QR = -VC*(m3_c*eps - m4_c*eps - m3_c + m4_c)/ts

#The first run is run without initialization
state_y = []
state_ydot = []
#Running the SMB model 
for j in range(0,4):
    df,state_y,state_ydot = SMBmodel(Q1,Q2,Q3,Q4,QD,QE,QF,QR,ts,n_col,j,state_y,state_ydot)
    state_y,state_ydot = state_y,state_ydot

#%% Load data and plot the issue
ts = 180
data = []


filename = 'Langmuir_sim_gradient_multi_0.h5'
smb_model = Cadet()
smb_model.filename = filename
smb_model.load()
df = pd.DataFrame({'Time': smb_model.root.output.solution.solution_times,
                   'Ca_E': smb_model.root.output.solution.unit_002.solution_outlet[:,0],
                   'Cb_E': smb_model.root.output.solution.unit_002.solution_outlet[:,1],
                   'Ca_R': smb_model.root.output.solution.unit_003.solution_outlet[:,0],
                   'Cb_R': smb_model.root.output.solution.unit_003.solution_outlet[:,1] })
data.append(df)

filename = 'Langmuir_sim_gradient_multi_1.h5'
smb_model = Cadet()
smb_model.filename = filename
smb_model.load()
df = pd.DataFrame({'Time': smb_model.root.output.solution.solution_times,
                   'Ca_E': smb_model.root.output.solution.unit_002.solution_outlet[:,0],
                   'Cb_E': smb_model.root.output.solution.unit_002.solution_outlet[:,1],
                   'Ca_R': smb_model.root.output.solution.unit_003.solution_outlet[:,0],
                   'Cb_R': smb_model.root.output.solution.unit_003.solution_outlet[:,1] })
data.append(df)
filename = 'Langmuir_sim_gradient_multi_2.h5'
smb_model = Cadet()
smb_model.filename = filename
smb_model.load()
df = pd.DataFrame({'Time': smb_model.root.output.solution.solution_times,
                   'Ca_E': smb_model.root.output.solution.unit_002.solution_outlet[:,0],
                   'Cb_E': smb_model.root.output.solution.unit_002.solution_outlet[:,1],
                   'Ca_R': smb_model.root.output.solution.unit_003.solution_outlet[:,0],
                   'Cb_R': smb_model.root.output.solution.unit_003.solution_outlet[:,1] })    
data.append(df)
filename = 'Langmuir_sim_gradient_multi_3.h5'
smb_model = Cadet()
smb_model.filename = filename
smb_model.load()
df = pd.DataFrame({'Time': smb_model.root.output.solution.solution_times,
                   'Ca_E': smb_model.root.output.solution.unit_002.solution_outlet[:,0],
                   'Cb_E': smb_model.root.output.solution.unit_002.solution_outlet[:,1],
                   'Ca_R': smb_model.root.output.solution.unit_003.solution_outlet[:,0],
                   'Cb_R': smb_model.root.output.solution.unit_003.solution_outlet[:,1] })
data.append(df)


#plotting the problem
plt.plot(np.concatenate((np.array(data[0]['Ca_R'][:]),np.array(data[1]['Ca_R'][:]),np.array(data[2]['Ca_R'][:]))))
plt.plot([4*4*ts,4*4*ts],[0,0.005],'--',label='First run')
plt.plot([4*4*ts+4*4*ts,4*4*ts+4*4*ts],[0,0.005],'--',label='second run')
plt.plot([4*4*ts+2*4*4*ts,4*4*ts+2*4*4*ts],[0,0.005],'--',label='third run')
plt.legend()
plt.xlabel('Time(s)')
plt.ylabel('Concentration')
11

Hey Jesper,
I’m not familiar with data frames, so I modified the post processing a bit.

#%% Load data and plot the issue
ts = 180
data = []


filename = 'Langmuir_sim_gradient_multi_0.h5'
smb_model = Cadet()
smb_model.filename = filename
smb_model.load()

Ca_E = smb_model.root.output.solution.unit_002.solution_outlet[:,0]
Cb_E = smb_model.root.output.solution.unit_002.solution_outlet[:,1]
Ca_R = smb_model.root.output.solution.unit_003.solution_outlet[:,0]
Cb_R = smb_model.root.output.solution.unit_003.solution_outlet[:,1]

filename = 'Langmuir_sim_gradient_multi_1.h5'
smb_model = Cadet()
smb_model.filename = filename
smb_model.load()
Ca_E = np.hstack((Ca_E, smb_model.root.output.solution.unit_002.solution_outlet[:,0]))
Cb_E = np.hstack((Cb_E, smb_model.root.output.solution.unit_002.solution_outlet[:,1]))
Ca_R = np.hstack((Ca_R, smb_model.root.output.solution.unit_003.solution_outlet[:,0]))
Cb_R = np.hstack((Cb_R, smb_model.root.output.solution.unit_003.solution_outlet[:,1]))

filename = 'Langmuir_sim_gradient_multi_2.h5'
smb_model = Cadet()
smb_model.filename = filename
smb_model.load()
Ca_E = np.hstack((Ca_E, smb_model.root.output.solution.unit_002.solution_outlet[:,0]))
Cb_E = np.hstack((Cb_E, smb_model.root.output.solution.unit_002.solution_outlet[:,1]))
Ca_R = np.hstack((Ca_R, smb_model.root.output.solution.unit_003.solution_outlet[:,0]))
Cb_R = np.hstack((Cb_R, smb_model.root.output.solution.unit_003.solution_outlet[:,1]))

filename = 'Langmuir_sim_gradient_multi_3.h5'
smb_model = Cadet()
smb_model.filename = filename
smb_model.load()
Ca_E = np.hstack((Ca_E, smb_model.root.output.solution.unit_002.solution_outlet[:,0]))
Cb_E = np.hstack((Cb_E, smb_model.root.output.solution.unit_002.solution_outlet[:,1]))
Ca_R = np.hstack((Ca_R, smb_model.root.output.solution.unit_003.solution_outlet[:,0]))
Cb_R = np.hstack((Cb_R, smb_model.root.output.solution.unit_003.solution_outlet[:,1]))

fig, ax = plt.subplots()
ax.plot(Ca_E)

The Raffinate still looks the same. So at least we know there is no issue there.
image

However, if we examine the Extract stream, I actually don’t think there is an issue with the initial values. Otherwise, this would look very different. Also in your Raffinate stream, the first cycle looks different than the others. So do you think it’s possible that SMB is not correctly configured yet? Where did you you get the parameters from? Do you account for non-idealities like binding kinetics, dispersion etc. in you process design?

image

12

Hey Johannes,

I dont think the parameters are the reason.
I have tried to run the example from the CADET-tutorials for SMB ( CADET-Tutorial/01_advanced_chromatography_demonstration_solution.ipynb at master · modsim/CADET-Tutorial · GitHub) with the same functions as used in the tutorials and using linear isotherm.
When I run the simulation, store initial values and run another round, it shows the same result.
The simple code is shown below (run, store values, run).

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'}
    
    
    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.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()
    cadet.load()
    
    # Remove files 
    if file_name is None:
        os.remove(os.path.join(tempfile.tempdir, f + '.h5'))

    # Raise error if simulation fails
    if data.returncode == 0:
        print("Simulation completed successfully")
    else:
        print(data)
        raise Exception("Simulation failed")
#%%

import matplotlib.pyplot as plt
import numpy as np

from cadet import Cadet
Cadet.cadet_path = r'C:\Users\jespfra\Anaconda3\bin\cadet-cli' #working pc



smb_model = get_cadet_template(n_units=8)

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

## Eluent
smb_model.root.input.model.unit_001.unit_type = 'INLET'
smb_model.root.input.model.unit_001.ncomp = 2
smb_model.root.input.model.unit_001.inlet_type = 'PIECEWISE_CUBIC_POLY'

## Extract 
smb_model.root.input.model.unit_002.ncomp = 2
smb_model.root.input.model.unit_002.unit_type = 'OUTLET'

## Raffinate
smb_model.root.input.model.unit_003.ncomp = 2
smb_model.root.input.model.unit_003.unit_type = 'OUTLET'

## Columns
smb_model.root.input.model.unit_004.unit_type = 'LUMPED_RATE_MODEL_WITHOUT_PORES'
smb_model.root.input.model.unit_004.ncomp = 2

smb_model.root.input.model.unit_004.col_dispersion = 3.8148E-20
smb_model.root.input.model.unit_004.col_length = 0.25
smb_model.root.input.model.unit_004.total_porosity = 0.83
smb_model.root.input.model.unit_004.velocity = 1
smb_model.root.input.model.unit_004.cross_section_area = 3.141592653589793E-4

smb_model.root.input.model.unit_004.adsorption_model = 'LINEAR'
smb_model.root.input.model.unit_004.adsorption.is_kinetic = 0
smb_model.root.input.model.unit_004.adsorption.lin_ka = [5.72, 7.7]
smb_model.root.input.model.unit_004.adsorption.lin_kd = [1.0, 1.0]

smb_model.root.input.model.unit_004.init_c = [0.0,0.0]
smb_model.root.input.model.unit_004.init_q = [0.0,0.0]

### Copy column models
smb_model.root.input.model.unit_005 = smb_model.root.input.model.unit_004
smb_model.root.input.model.unit_006 = smb_model.root.input.model.unit_004
smb_model.root.input.model.unit_007 = smb_model.root.input.model.unit_004

### Discretization column settings
set_discretization(smb_model, n_col=40, n_bound=[1,1])

#To write out last output to check for steady state
smb_model.root.input['return'].write_solution_last = True

n_cycles = 10
switch_time = 180

# Sections
smb_model.root.input.solver.sections.nsec = 4*n_cycles
smb_model.root.input.solver.sections.section_times = [0]
for i in range(n_cycles):
    smb_model.root.input.solver.sections.section_times.append((4*i+1)*switch_time)
    smb_model.root.input.solver.sections.section_times.append((4*i+2)*switch_time)
    smb_model.root.input.solver.sections.section_times.append((4*i+3)*switch_time)
    smb_model.root.input.solver.sections.section_times.append((4*i+4)*switch_time)    

## Feed and Eluent concentration
smb_model.root.input.model.unit_000.sec_000.const_coeff = [0.003,0.003]
smb_model.root.input.model.unit_001.sec_000.const_coeff = [0.0,0.0]

# Connections
smb_model.root.input.model.connections.nswitches = 4

smb_model.root.input.model.connections.switch_000.section = 0
smb_model.root.input.model.connections.switch_000.connections = [
    4, 5, -1, -1, 7.66e-7,
    5, 6, -1, -1, 7.66e-7,
    6, 7, -1, -1, 8.08e-7,
    7, 4, -1, -1, 8.08e-7,
    0, 4, -1, -1, 0.98e-7,
    1, 6, -1, -1, 1.96e-7,
    4, 2, -1, -1, 1.4e-7,
    6, 3, -1, -1, 1.54e-7
]

smb_model.root.input.model.connections.switch_001.section = 1
smb_model.root.input.model.connections.switch_001.connections = [
    4, 5, -1, -1, 8.08e-7,
    5, 6, -1, -1, 7.66e-7,
    6, 7, -1, -1, 7.66e-7,
    7, 4, -1, -1, 8.08e-7,
    0, 5, -1, -1, 0.98e-7,
    1, 7, -1, -1, 1.96e-7,
    5, 2, -1, -1, 1.4e-7,
    7, 3, -1, -1, 1.54e-7
]

smb_model.root.input.model.connections.switch_002.section = 2
smb_model.root.input.model.connections.switch_002.connections = [
    4, 5, -1, -1, 8.08e-7,
    5, 6, -1, -1, 8.08e-7,
    6, 7, -1, -1, 7.66e-7,
    7, 4, -1, -1, 7.66e-7,
    0, 6, -1, -1, 0.98e-7,
    1, 4, -1, -1, 1.96e-7,
    6, 2, -1, -1, 1.4e-7,
    4, 3, -1, -1, 1.54e-7
]

smb_model.root.input.model.connections.switch_003.section = 3
smb_model.root.input.model.connections.switch_003.connections = [
    4, 5, -1, -1, 7.66e-7,
    5, 6, -1, -1, 8.08e-7,
    6, 7, -1, -1, 8.08e-7,
    7, 4, -1, -1, 7.66e-7,
    0, 7, -1, -1, 0.98e-7,
    1, 5, -1, -1, 1.96e-7,
    7, 2, -1, -1, 1.4e-7,
    5, 3, -1, -1, 1.54e-7
]

#set the times that the simulator writes out data for
smb_model.root.input.solver.user_solution_times = np.linspace(0, n_cycles*4*switch_time, int(n_cycles*4*switch_time))





# Run simulation and plot results
smb_model.filename = 'Johannes.h5'
smb_model.save()

#run model
data = smb_model.run()

if data.returncode == 0:
    print("Simulation completed successfully")
    smb_model.load()

else:
    print(data)
    raise Exception("Simulation failed")
    
#If cyclic steady state has not been reached
state_y = smb_model.root.output.last_state_y
state_ydot = smb_model.root.output.last_state_ydot

#Run the same model the 2nd time 

smb_model1 = get_cadet_template(n_units=8)

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

## Eluent
smb_model1.root.input.model.unit_001.unit_type = 'INLET'
smb_model1.root.input.model.unit_001.ncomp = 2
smb_model1.root.input.model.unit_001.inlet_type = 'PIECEWISE_CUBIC_POLY'

## Extract 
smb_model1.root.input.model.unit_002.ncomp = 2
smb_model1.root.input.model.unit_002.unit_type = 'OUTLET'

## Raffinate
smb_model1.root.input.model.unit_003.ncomp = 2
smb_model1.root.input.model.unit_003.unit_type = 'OUTLET'

## Columns
smb_model1.root.input.model.unit_004.unit_type = 'LUMPED_RATE_MODEL_WITHOUT_PORES'
smb_model1.root.input.model.unit_004.ncomp = 2

smb_model1.root.input.model.unit_004.col_dispersion = 3.8148E-20
smb_model1.root.input.model.unit_004.col_length = 0.25
smb_model1.root.input.model.unit_004.total_porosity = 0.83
smb_model1.root.input.model.unit_004.velocity = 1
smb_model1.root.input.model.unit_004.cross_section_area = 3.141592653589793E-4

smb_model1.root.input.model.unit_004.adsorption_model = 'LINEAR'
smb_model1.root.input.model.unit_004.adsorption.is_kinetic = 0
smb_model1.root.input.model.unit_004.adsorption.lin_ka = [5.72, 7.7]
smb_model1.root.input.model.unit_004.adsorption.lin_kd = [1.0, 1.0]

# smb_model1.root.input.model.unit_004.init_c = [0.0,0.0]
# smb_model1.root.input.model.unit_004.init_q = [0.0,0.0]
smb_model1.root.input.model.init_state_y = state_y
smb_model1.root.input.model.init_state_ydot = state_ydot


### Copy column models
smb_model1.root.input.model.unit_005 = smb_model1.root.input.model.unit_004
smb_model1.root.input.model.unit_006 = smb_model1.root.input.model.unit_004
smb_model1.root.input.model.unit_007 = smb_model1.root.input.model.unit_004

### Discretization column settings
set_discretization(smb_model1, n_col=40, n_bound=[1,1])

#To write out last output to check for steady state
smb_model1.root.input['return'].write_solution_last = True

n_cycles = 10
switch_time = 180

# Sections
smb_model1.root.input.solver.sections.nsec = 4*n_cycles
smb_model1.root.input.solver.sections.section_times = [0]
for i in range(n_cycles):
    smb_model1.root.input.solver.sections.section_times.append((4*i+1)*switch_time)
    smb_model1.root.input.solver.sections.section_times.append((4*i+2)*switch_time)
    smb_model1.root.input.solver.sections.section_times.append((4*i+3)*switch_time)
    smb_model1.root.input.solver.sections.section_times.append((4*i+4)*switch_time)    

## Feed and Eluent concentration
smb_model1.root.input.model.unit_000.sec_000.const_coeff = [0.003,0.003]
smb_model1.root.input.model.unit_001.sec_000.const_coeff = [0.0,0.0]

# Connections
smb_model1.root.input.model.connections.nswitches = 4

smb_model1.root.input.model.connections.switch_000.section = 0
smb_model1.root.input.model.connections.switch_000.connections = [
    4, 5, -1, -1, 7.66e-7,
    5, 6, -1, -1, 7.66e-7,
    6, 7, -1, -1, 8.08e-7,
    7, 4, -1, -1, 8.08e-7,
    0, 4, -1, -1, 0.98e-7,
    1, 6, -1, -1, 1.96e-7,
    4, 2, -1, -1, 1.4e-7,
    6, 3, -1, -1, 1.54e-7
]

smb_model1.root.input.model.connections.switch_001.section = 1
smb_model1.root.input.model.connections.switch_001.connections = [
    4, 5, -1, -1, 8.08e-7,
    5, 6, -1, -1, 7.66e-7,
    6, 7, -1, -1, 7.66e-7,
    7, 4, -1, -1, 8.08e-7,
    0, 5, -1, -1, 0.98e-7,
    1, 7, -1, -1, 1.96e-7,
    5, 2, -1, -1, 1.4e-7,
    7, 3, -1, -1, 1.54e-7
]

smb_model1.root.input.model.connections.switch_002.section = 2
smb_model1.root.input.model.connections.switch_002.connections = [
    4, 5, -1, -1, 8.08e-7,
    5, 6, -1, -1, 8.08e-7,
    6, 7, -1, -1, 7.66e-7,
    7, 4, -1, -1, 7.66e-7,
    0, 6, -1, -1, 0.98e-7,
    1, 4, -1, -1, 1.96e-7,
    6, 2, -1, -1, 1.4e-7,
    4, 3, -1, -1, 1.54e-7
]

smb_model1.root.input.model.connections.switch_003.section = 3
smb_model1.root.input.model.connections.switch_003.connections = [
    4, 5, -1, -1, 7.66e-7,
    5, 6, -1, -1, 8.08e-7,
    6, 7, -1, -1, 8.08e-7,
    7, 4, -1, -1, 7.66e-7,
    0, 7, -1, -1, 0.98e-7,
    1, 5, -1, -1, 1.96e-7,
    7, 2, -1, -1, 1.4e-7,
    5, 3, -1, -1, 1.54e-7
]

#set the times that the simulator writes out data for
smb_model1.root.input.solver.user_solution_times = np.linspace(0, n_cycles*4*switch_time, int(n_cycles*4*switch_time))



# Run simulation and plot results
smb_model1.filename = 'Johannes1.h5'
smb_model1.save()

#run model
data = smb_model1.run()

if data.returncode == 0:
    print("Simulation completed successfully")
    smb_model1.load()

else:
    print(data)
    raise Exception("Simulation failed")

plt.figure()
plt.plot(np.concatenate((np.array(smb_model.root.output.solution.unit_002.solution_outlet[:,1]),np.array(smb_model1.root.output.solution.unit_002.solution_outlet[:,1])),axis=0))
plt.xlabel('Time(s)')
plt.ylabel('Concentration (mM)')
plt.title('Extract')

plt.figure()
plt.plot(np.concatenate((np.array(smb_model.root.output.solution.unit_003.solution_outlet[:,0]),np.array(smb_model1.root.output.solution.unit_003.solution_outlet[:,0])),axis=0))
# plt.plot(np.concatenate((np.array(smb_model.root.output.solution.unit_003.solution_outlet[:,1]),np.array(smb_model1.root.output.solution.unit_003.solution_outlet[:,1])),axis=0))
plt.xlabel('Time(s)')
plt.ylabel('Concentration (mM)')
plt.title('Raffinate')
13

For those interested.
I found a mistake in the SMB setup.
I had switched around the extract and raffinate in the flows.
Maybe it should be corrected in the tutorial for the SMB as well?

Below is the corrected code for the tutorial example above.

However, it does not solve the issue.

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'}
    
    
    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.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()
    cadet.load()
    
    # Remove files 
    if file_name is None:
        os.remove(os.path.join(tempfile.tempdir, f + '.h5'))

    # Raise error if simulation fails
    if data.returncode == 0:
        print("Simulation completed successfully")
    else:
        print(data)
        raise Exception("Simulation failed")
#%%

import matplotlib.pyplot as plt
import numpy as np

from cadet import Cadet
Cadet.cadet_path = r'C:\Users\jespfra\Anaconda3\bin\cadet-cli' #working pc



smb_model = get_cadet_template(n_units=8)

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

## Eluent
smb_model.root.input.model.unit_001.unit_type = 'INLET'
smb_model.root.input.model.unit_001.ncomp = 2
smb_model.root.input.model.unit_001.inlet_type = 'PIECEWISE_CUBIC_POLY'

## Extract 
smb_model.root.input.model.unit_002.ncomp = 2
smb_model.root.input.model.unit_002.unit_type = 'OUTLET'

## Raffinate
smb_model.root.input.model.unit_003.ncomp = 2
smb_model.root.input.model.unit_003.unit_type = 'OUTLET'

## Columns
smb_model.root.input.model.unit_004.unit_type = 'LUMPED_RATE_MODEL_WITHOUT_PORES'
smb_model.root.input.model.unit_004.ncomp = 2

smb_model.root.input.model.unit_004.col_dispersion = 3.8148E-20
smb_model.root.input.model.unit_004.col_length = 0.25
smb_model.root.input.model.unit_004.total_porosity = 0.83
smb_model.root.input.model.unit_004.velocity = 1
smb_model.root.input.model.unit_004.cross_section_area = 3.141592653589793E-4

smb_model.root.input.model.unit_004.adsorption_model = 'LINEAR'
smb_model.root.input.model.unit_004.adsorption.is_kinetic = 0
smb_model.root.input.model.unit_004.adsorption.lin_ka = [5.72, 7.7]
smb_model.root.input.model.unit_004.adsorption.lin_kd = [1.0, 1.0]

smb_model.root.input.model.unit_004.init_c = [0.0,0.0]
smb_model.root.input.model.unit_004.init_q = [0.0,0.0]

### Copy column models
smb_model.root.input.model.unit_005 = smb_model.root.input.model.unit_004
smb_model.root.input.model.unit_006 = smb_model.root.input.model.unit_004
smb_model.root.input.model.unit_007 = smb_model.root.input.model.unit_004

### Discretization column settings
set_discretization(smb_model, n_col=40, n_bound=[1,1])

#To write out last output to check for steady state
smb_model.root.input['return'].write_solution_last = True

n_cycles = 10
switch_time = 180

# Sections
smb_model.root.input.solver.sections.nsec = 4*n_cycles
smb_model.root.input.solver.sections.section_times = [0]
for i in range(n_cycles):
    smb_model.root.input.solver.sections.section_times.append((4*i+1)*switch_time)
    smb_model.root.input.solver.sections.section_times.append((4*i+2)*switch_time)
    smb_model.root.input.solver.sections.section_times.append((4*i+3)*switch_time)
    smb_model.root.input.solver.sections.section_times.append((4*i+4)*switch_time)    

## Feed and Eluent concentration
smb_model.root.input.model.unit_000.sec_000.const_coeff = [0.003,0.003]
smb_model.root.input.model.unit_001.sec_000.const_coeff = [0.0,0.0]

# Connections
smb_model.root.input.model.connections.nswitches = 4

smb_model.root.input.model.connections.switch_000.section = 0
smb_model.root.input.model.connections.switch_000.connections = [
    4, 5, -1, -1, 7.66e-7,
    5, 6, -1, -1, 7.66e-7,
    6, 7, -1, -1, 8.08e-7,
    7, 4, -1, -1, 8.08e-7,
    0, 4, -1, -1, 0.98e-7,
    1, 6, -1, -1, 1.96e-7,
    4, 3, -1, -1, 1.4e-7,
    6, 2, -1, -1, 1.54e-7
]

smb_model.root.input.model.connections.switch_001.section = 1
smb_model.root.input.model.connections.switch_001.connections = [
    4, 5, -1, -1, 8.08e-7,
    5, 6, -1, -1, 7.66e-7,
    6, 7, -1, -1, 7.66e-7,
    7, 4, -1, -1, 8.08e-7,
    0, 5, -1, -1, 0.98e-7,
    1, 7, -1, -1, 1.96e-7,
    5, 3, -1, -1, 1.4e-7,
    7, 2, -1, -1, 1.54e-7
]

smb_model.root.input.model.connections.switch_002.section = 2
smb_model.root.input.model.connections.switch_002.connections = [
    4, 5, -1, -1, 8.08e-7,
    5, 6, -1, -1, 8.08e-7,
    6, 7, -1, -1, 7.66e-7,
    7, 4, -1, -1, 7.66e-7,
    0, 6, -1, -1, 0.98e-7,
    1, 4, -1, -1, 1.96e-7,
    6, 3, -1, -1, 1.4e-7,
    4, 2, -1, -1, 1.54e-7
]

smb_model.root.input.model.connections.switch_003.section = 3
smb_model.root.input.model.connections.switch_003.connections = [
    4, 5, -1, -1, 7.66e-7,
    5, 6, -1, -1, 8.08e-7,
    6, 7, -1, -1, 8.08e-7,
    7, 4, -1, -1, 7.66e-7,
    0, 7, -1, -1, 0.98e-7,
    1, 5, -1, -1, 1.96e-7,
    7, 3, -1, -1, 1.4e-7,
    5, 2, -1, -1, 1.54e-7
]

#set the times that the simulator writes out data for
smb_model.root.input.solver.user_solution_times = np.linspace(0, n_cycles*4*switch_time, int(n_cycles*4*switch_time))





# Run simulation and plot results
smb_model.filename = 'Johannes.h5'
smb_model.save()

#run model
data = smb_model.run()

if data.returncode == 0:
    print("Simulation completed successfully")
    smb_model.load()

else:
    print(data)
    raise Exception("Simulation failed")
    
#If cyclic steady state has not been reached
state_y = smb_model.root.output.last_state_y
state_ydot = smb_model.root.output.last_state_ydot

#Run the same model the 2nd time 

smb_model1 = get_cadet_template(n_units=8)

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

## Eluent
smb_model1.root.input.model.unit_001.unit_type = 'INLET'
smb_model1.root.input.model.unit_001.ncomp = 2
smb_model1.root.input.model.unit_001.inlet_type = 'PIECEWISE_CUBIC_POLY'

## Extract 
smb_model1.root.input.model.unit_002.ncomp = 2
smb_model1.root.input.model.unit_002.unit_type = 'OUTLET'

## Raffinate
smb_model1.root.input.model.unit_003.ncomp = 2
smb_model1.root.input.model.unit_003.unit_type = 'OUTLET'

## Columns
smb_model1.root.input.model.unit_004.unit_type = 'LUMPED_RATE_MODEL_WITHOUT_PORES'
smb_model1.root.input.model.unit_004.ncomp = 2

smb_model1.root.input.model.unit_004.col_dispersion = 3.8148E-20
smb_model1.root.input.model.unit_004.col_length = 0.25
smb_model1.root.input.model.unit_004.total_porosity = 0.83
smb_model1.root.input.model.unit_004.velocity = 1
smb_model1.root.input.model.unit_004.cross_section_area = 3.141592653589793E-4

smb_model1.root.input.model.unit_004.adsorption_model = 'LINEAR'
smb_model1.root.input.model.unit_004.adsorption.is_kinetic = 0
smb_model1.root.input.model.unit_004.adsorption.lin_ka = [5.72, 7.7]
smb_model1.root.input.model.unit_004.adsorption.lin_kd = [1.0, 1.0]

# smb_model1.root.input.model.unit_004.init_c = [0.0,0.0]
# smb_model1.root.input.model.unit_004.init_q = [0.0,0.0]
smb_model1.root.input.model.init_state_y = state_y
smb_model1.root.input.model.init_state_ydot = state_ydot


### Copy column models
smb_model1.root.input.model.unit_005 = smb_model1.root.input.model.unit_004
smb_model1.root.input.model.unit_006 = smb_model1.root.input.model.unit_004
smb_model1.root.input.model.unit_007 = smb_model1.root.input.model.unit_004

### Discretization column settings
set_discretization(smb_model1, n_col=40, n_bound=[1,1])

#To write out last output to check for steady state
smb_model1.root.input['return'].write_solution_last = True

n_cycles = 10
switch_time = 180

# Sections
smb_model1.root.input.solver.sections.nsec = 4*n_cycles
smb_model1.root.input.solver.sections.section_times = [0]
for i in range(n_cycles):
    smb_model1.root.input.solver.sections.section_times.append((4*i+1)*switch_time)
    smb_model1.root.input.solver.sections.section_times.append((4*i+2)*switch_time)
    smb_model1.root.input.solver.sections.section_times.append((4*i+3)*switch_time)
    smb_model1.root.input.solver.sections.section_times.append((4*i+4)*switch_time)    

## Feed and Eluent concentration
smb_model1.root.input.model.unit_000.sec_000.const_coeff = [0.003,0.003]
smb_model1.root.input.model.unit_001.sec_000.const_coeff = [0.0,0.0]

# Connections
smb_model1.root.input.model.connections.nswitches = 4

smb_model1.root.input.model.connections.switch_000.section = 0
smb_model1.root.input.model.connections.switch_000.connections = [
    4, 5, -1, -1, 7.66e-7,
    5, 6, -1, -1, 7.66e-7,
    6, 7, -1, -1, 8.08e-7,
    7, 4, -1, -1, 8.08e-7,
    0, 4, -1, -1, 0.98e-7,
    1, 6, -1, -1, 1.96e-7,
    4, 3, -1, -1, 1.4e-7,
    6, 2, -1, -1, 1.54e-7
]

smb_model1.root.input.model.connections.switch_001.section = 1
smb_model1.root.input.model.connections.switch_001.connections = [
    4, 5, -1, -1, 8.08e-7,
    5, 6, -1, -1, 7.66e-7,
    6, 7, -1, -1, 7.66e-7,
    7, 4, -1, -1, 8.08e-7,
    0, 5, -1, -1, 0.98e-7,
    1, 7, -1, -1, 1.96e-7,
    5, 3, -1, -1, 1.4e-7,
    7, 2, -1, -1, 1.54e-7
]

smb_model1.root.input.model.connections.switch_002.section = 2
smb_model1.root.input.model.connections.switch_002.connections = [
    4, 5, -1, -1, 8.08e-7,
    5, 6, -1, -1, 8.08e-7,
    6, 7, -1, -1, 7.66e-7,
    7, 4, -1, -1, 7.66e-7,
    0, 6, -1, -1, 0.98e-7,
    1, 4, -1, -1, 1.96e-7,
    6, 3, -1, -1, 1.4e-7,
    4, 2, -1, -1, 1.54e-7
]

smb_model1.root.input.model.connections.switch_003.section = 3
smb_model1.root.input.model.connections.switch_003.connections = [
    4, 5, -1, -1, 7.66e-7,
    5, 6, -1, -1, 8.08e-7,
    6, 7, -1, -1, 8.08e-7,
    7, 4, -1, -1, 7.66e-7,
    0, 7, -1, -1, 0.98e-7,
    1, 5, -1, -1, 1.96e-7,
    7, 3, -1, -1, 1.4e-7,
    5, 2, -1, -1, 1.54e-7
]

#set the times that the simulator writes out data for
smb_model1.root.input.solver.user_solution_times = np.linspace(0, n_cycles*4*switch_time, int(n_cycles*4*switch_time))



# Run simulation and plot results
smb_model1.filename = 'Johannes1.h5'
smb_model1.save()

#run model
data = smb_model1.run()

if data.returncode == 0:
    print("Simulation completed successfully")
    smb_model1.load()

else:
    print(data)
    raise Exception("Simulation failed")


plt.figure()
plt.plot(np.concatenate((np.array(smb_model.root.output.solution.unit_002.solution_outlet[:,0]),np.array(smb_model1.root.output.solution.unit_002.solution_outlet[:,0])),axis=0))
plt.xlabel('Time(s)')
plt.ylabel('Concentration (mM)')
plt.title('Extract')

plt.figure()
plt.plot(np.concatenate((np.array(smb_model.root.output.solution.unit_003.solution_outlet[:,1]),np.array(smb_model1.root.output.solution.unit_003.solution_outlet[:,1])),axis=0))
# plt.plot(np.concatenate((np.array(smb_model.root.output.solution.unit_003.solution_outlet[:,1]),np.array(smb_model1.root.output.solution.unit_003.solution_outlet[:,1])),axis=0))
plt.xlabel('Time(s)')
plt.ylabel('Concentration (mM)')
plt.title('Raffinate')