Hi,
that code also didn’t include the import for ConcentrationConverter
that I had asked for, but I made it work anyways with some guesswork. Please ensure your code runs as an isolated minimal reproducible example in a fresh Python Kernel. Also, please include the error message you get next time.
With that said, here are the errors in your code:
- Missing call to breakthrough_radial_langmuir(): You define
breakthrough_radial_langmuir()
but never call it. Therefore your code does not know what langmuir_model
is in the global scope. If the reason is unclear please check out an article on scopes, e.g. this. This can be fixed with langmuir_model = breakthrough_radial_langmuir(Q=2e-6, load_duration=15)
- Inconsistent number of components & parameters. You specified a lot of parameters for one component (
mcl_ka = [1., ]
, init_c = [0.0, ]
, etc.) but told the unit operation to expect two components unit_001.ncomp = 2
. To fix: set n_comp = 1
or give all parameters that require a parameter per component two values mcl_ka = [1.0, 2.0]
for example.
- Inproper arguments for np.linspace: You gave 5 arguments to np.linspace:
np.linspace(0, int(load_duration), int(load_duration)+1, int(wash_elute_duration), int(wash_elute_duration)+1)
. np.linspace commonly takes three: start, stop, num (number of entries on the linear scale). So np.linspace(start=0, stop=int(wash_elute_duration), num=300)
is I think what you wanted.
- Accessing component solutions without splitting components: You try to access the solution with
unit_002['solution_outlet_comp_000']
. Because you didn’t include the definition of get_cadet_template()
I had to guess that you used the default one from the CADET-Python-Tutorial. This function does not set the split_components_data
flag to True
by default. Therefore unit_002['solution_outlet_comp_000']
is going to be empty. You can access the solution with unit_002['solution_outlet']
or use langmuir_model = get_cadet_template(n_units=3, split_components_data=True)
with unit_002['solution_outlet_comp_000']
.
- Missing unit operation type in set_discretization(): Again, you didn’t include the definition of
set_discretization()
so I can only guess. But if it’s the function from the CADET-Python-Tutorial, then it only sets discretizations for columns = {'GENERAL_RATE_MODEL', 'LUMPED_RATE_MODEL_WITH_PORES', 'LUMPED_RATE_MODEL_WITHOUT_PORES'}
and not RADIAL_LUMPED_RATE_MODEL_WITH_PORES
. This unit operation type needs to be added to that function.
- Some code duplications These didn’t break anything but still made reading and debugging harder, as they set the same variable in multiple places. Such as: line 33:
langmuir_model.root.input.model.unit_000.sec_000.const_coeff = [c_,] # 10mg/ml langmuir_model.root.input.model.unit_000.sec_001.const_coeff = [0.0,] # Wash & Elute step
line 74:
langmuir_model.root.input.model.unit_000.sec_000.const_coeff = [c_,0] # 10mg/ml langmuir_model.root.input.model.unit_000.sec_001.const_coeff = [0.0,0] # Wash & Elute step
Here is a complete and working code example:
import os
import tempfile
import numpy as np
from matplotlib import pyplot as plt
from cadet import Cadet
tempfile.tempdir = os.path.join('')
Cadet.cadet_path = SET_CADET_PATH_HERE
def get_cadet_template(n_units=3, split_components_data=False):
cadet_template = Cadet()
cadet_template.root.input.model.nunits = n_units
# Store solution
cadet_template.root.input['return'].split_components_data = split_components_data
cadet_template.root.input['return'].split_ports_data = 0
cadet_template.root.input['return'].unit_000.write_solution_inlet = 1
cadet_template.root.input['return'].unit_000.write_solution_outlet = 1
cadet_template.root.input['return'].unit_000.write_solution_bulk = 1
cadet_template.root.input['return'].unit_000.write_solution_particle = 1
cadet_template.root.input['return'].unit_000.write_solution_solid = 1
cadet_template.root.input['return'].unit_000.write_solution_flux = 1
cadet_template.root.input['return'].unit_000.write_solution_volume = 1
cadet_template.root.input['return'].unit_000.write_coordinates = 1
cadet_template.root.input['return'].unit_000.write_sens_outlet = 1
for unit in range(n_units):
cadet_template.root.input['return']['unit_{0:03d}'.format(unit)] = cadet_template.root.input['return'].unit_000
# Tolerances for the time integrator
cadet_template.root.input.solver.time_integrator.abstol = 1e-6
cadet_template.root.input.solver.time_integrator.algtol = 1e-10
cadet_template.root.input.solver.time_integrator.reltol = 1e-6
cadet_template.root.input.solver.time_integrator.init_step_size = 1e-6
cadet_template.root.input.solver.time_integrator.max_steps = 1000000
# Solver settings
cadet_template.root.input.model.solver.gs_type = 1
cadet_template.root.input.model.solver.max_krylov = 0
cadet_template.root.input.model.solver.max_restarts = 10
cadet_template.root.input.model.solver.schur_safety = 1e-8
# Run the simulation on single thread
cadet_template.root.input.solver.nthreads = 1
return cadet_template
def set_discretization(model, n_bound=None, n_col=20, n_par_types=1):
columns = {'GENERAL_RATE_MODEL', 'LUMPED_RATE_MODEL_WITH_PORES', 'LUMPED_RATE_MODEL_WITHOUT_PORES',
"RADIAL_LUMPED_RATE_MODEL_WITH_PORES"}
for unit_name, unit in model.root.input.model.items():
if 'unit_' in unit_name and unit.unit_type in columns:
unit.discretization.ncol = n_col
unit.discretization.npar = 5
unit.discretization.npartype = n_par_types
if n_bound is None:
n_bound = unit.ncomp * [0]
unit.discretization.nbound = n_bound
unit.discretization.par_disc_type = 'EQUIDISTANT_PAR'
unit.discretization.use_analytic_jacobian = 1
unit.discretization.reconstruction = 'WENO'
unit.discretization.WENO_ORDER = 1
unit.discretization.gs_type = 1
unit.discretization.max_krylov = 0
unit.discretization.max_restarts = 10
unit.discretization.schur_safety = 1.0e-8
unit.discretization.weno.boundary_model = 0
unit.discretization.weno.weno_eps = 1e-10
unit.discretization.weno.weno_order = 3
def run_simulation(cadet, file_name=None):
if file_name is None:
f = next(tempfile._get_candidate_names())
cadet.filename = os.path.join(tempfile.tempdir, f + '.h5')
else:
cadet.filename = file_name
# save the simulation
cadet.save()
# run the simulation and load results
data = cadet.run(timeout=10)
print(data)
cadet.load()
# Remove files
if file_name is None:
os.remove(os.path.join(tempfile.tempdir, f + '.h5'))
def breakthrough_radial_langmuir(Q: float, load_duration: float, c=10, wash_elute_duration=110, kDa=150):
langmuir_model = get_cadet_template(n_units=3, split_components_data=True) # 3 = Inlet + Column + Outlet
# INLET
langmuir_model.root.input.model.unit_000.unit_type = 'INLET'
langmuir_model.root.input.model.unit_000.ncomp = 2
langmuir_model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
# Column
langmuir_model.root.input.model.unit_001.unit_type = 'RADIAL_LUMPED_RATE_MODEL_WITH_PORES'
langmuir_model.root.input.model.unit_001.ncomp = 2
langmuir_model.root.input.model.unit_001.col_length = 0.12
langmuir_model.root.input.model.unit_001.col_radius_inner = 0.0103 * 0.2
langmuir_model.root.input.model.unit_001.col_radius_outer = 0.0103
langmuir_model.root.input.model.unit_001.col_porosity = 0.29
langmuir_model.root.input.model.unit_001.par_porosity = 0.95
langmuir_model.root.input.model.unit_001.par_radius = 0.00025
langmuir_model.root.input.model.unit_001.VELOCITY_COEFF = 0.1
langmuir_model.root.input.model.unit_001.CROSS_SECTION_AREA = 0.001
langmuir_model.root.input.model.unit_001.col_dispersion = 1e-8
langmuir_model.root.input.model.unit_001.film_diffusion = [6.9e-6, 6.9e-6]
langmuir_model.root.input.model.unit_001.adsorption_model = 'MULTI_COMPONENT_LANGMUIR'
langmuir_model.root.input.model.unit_001.adsorption.is_kinetic = 1
langmuir_model.root.input.model.unit_001.adsorption.mcl_ka = [1., 2, ]
langmuir_model.root.input.model.unit_001.adsorption.mcl_kd = [1., 1, ]
langmuir_model.root.input.model.unit_001.adsorption.mcl_qmax = [100., 50., ]
langmuir_model.root.input.model.unit_001.init_c = [0.0, 0.0]
langmuir_model.root.input.model.unit_001.init_q = [0.0, 0.0]
## Outlet
langmuir_model.root.input.model.unit_002.ncomp = 2
langmuir_model.root.input.model.unit_002.unit_type = 'OUTLET'
set_discretization(langmuir_model, n_bound=[0, 1], n_col=50)
# LOAD + (WASH + ELUTE) step
langmuir_model.root.input.solver.sections.nsec = 2
langmuir_model.root.input.solver.sections.section_times = [0.0, load_duration, wash_elute_duration]
## Inlet Profile
langmuir_model.root.input.model.unit_000.sec_000.const_coeff = [c, c] # 10mg/ml
langmuir_model.root.input.model.unit_000.sec_001.const_coeff = [0.0, 0.0] # Wash & Elute step
langmuir_model.root.input.solver.sections.section_continuity = [0, 0]
## Switches
langmuir_model.root.input.model.connections.nswitches = 1
langmuir_model.root.input.model.connections.switch_000.section = 0
langmuir_model.root.input.model.connections.switch_000.connections = [
0, 1, -1, -1, Q, # unit_000, unit_001, all components, all components, Q/ m^3/s
1, 2, -1, -1, Q, # unit_001, unit_002, all components, all components, Q/ m^3/s
]
# set the times that the simulator writes out data for
langmuir_model.root.input.solver.user_solution_times = np.linspace(start=0, stop=int(wash_elute_duration), num=300)
return langmuir_model
langmuir_model = breakthrough_radial_langmuir(Q=2e-6, load_duration=15)
run_simulation(langmuir_model)
time = langmuir_model.root.output.solution.solution_times
c1 = langmuir_model.root.output.solution.unit_002['solution_outlet_comp_000']
c2 = langmuir_model.root.output.solution.unit_002['solution_outlet_comp_001']
plt.plot(time, c1)
plt.plot(time, c2)
plt.title('Column (Outlet)')
plt.xlabel('$time~/~min$')
plt.ylabel('$concentration~/~mM$')
plt.show()