Error in IDASolve

[Yang]

Hey everyone,
I would like to ask if everyone can give me some suggestions to solve this error. I am using Mollerup 2006/2008 HIC Isotherm in Unified HIC Isotherm, but the following error occurred:

CompletedProcess(
    args=[PosixPath('/usr/local/bin/cadet-cli'), 'model.h5'],
    returncode=3,
    stdout=b"[Error: idasErrorHandler::200] In function 'IDASolve' of module 'IDAS', error code 'IDA_CONV_FAIL':\nAt t = 0 and h = 9.53674e-13, the corrector convergence failed repeatedly or with |h| = hmin.\n[Error: integrate::1367] IDASolve returned IDA_CONV_FAIL at t = 0\n",
    stderr=b'SOLVER ERROR: Error in IDASolve: IDA_CONV_FAIL at t = 0.000000\n'
)

The original code is also attached here.

model = Cadet()

model.root.input.model.nunits = 3

#  
n_comp = 2

model.root.input.model.unit_000.unit_type = 'INLET'
model.root.input.model.unit_000.ncomp = n_comp
model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'

# Column
model.root.input.model.unit_001.unit_type = 'LUMPED_RATE_MODEL_WITH_PORES'
model.root.input.model.unit_001.ncomp = n_comp

model.root.input.model.unit_001.col_length = 0.1  
model.root.input.model.unit_001.cross_section_area = 4.657e-5  
model.root.input.model.unit_001.col_porosity = 0.364  
model.root.input.model.unit_001.par_porosity = 0.799 
model.root.input.model.unit_001.par_radius = 3.75e-5  

model.root.input.model.unit_001.col_dispersion = 1.25e-7  
model.root.input.model.unit_001.film_diffusion = [1.25e-5, 1.9e-6]  


model.root.input.model.unit_001.adsorption_model = 'HIC_UNIFIED'
model.root.input.model.unit_001.adsorption.is_kinetic = True    # Kinetic binding
model.root.input.model.unit_001.adsorption.hicuni_ka = [0, 0.010459588]     
model.root.input.model.unit_001.adsorption.hicuni_kd = [0,0.02881429]     
model.root.input.model.unit_001.adsorption.hicuni_kp = [0,0.0]      
model.root.input.model.unit_001.adsorption.hicuni_ks = [0,4.279]    
model.root.input.model.unit_001.adsorption.hicuni_epsilon = [0,0.0]
model.root.input.model.unit_001.adsorption.hicuni_beta0 = 0
model.root.input.model.unit_001.adsorption.hicuni_beta1 = 0
model.root.input.model.unit_001.adsorption.hicuni_nu = [0,0.0]
model.root.input.model.unit_001.adsorption.hicuni_qmax = [0,5.0]

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

### Grid cells
model.root.input.model.unit_001.discretization.ncol = 20
model.root.input.model.unit_001.discretization.npar = 5

### Bound states
model.root.input.model.unit_001.discretization.nbound = [0, 1]

### Other options
model.root.input.model.unit_001.discretization.par_disc_type = 'EQUIDISTANT_PAR'
model.root.input.model.unit_001.discretization.use_analytic_jacobian = 1
model.root.input.model.unit_001.discretization.reconstruction = 'WENO'
model.root.input.model.unit_001.discretization.gs_type = 1
model.root.input.model.unit_001.discretization.max_krylov = 0
model.root.input.model.unit_001.discretization.max_restarts = 10
model.root.input.model.unit_001.discretization.schur_safety = 1.0e-8

model.root.input.model.unit_001.discretization.weno.boundary_model = 0
model.root.input.model.unit_001.discretization.weno.weno_eps = 1e-10
model.root.input.model.unit_001.discretization.weno.weno_order = 3

model.root.input.model.unit_002.unit_type = 'OUTLET'
model.root.input.model.unit_002.ncomp = n_comp

model.root.input.solver.sections.nsec = 3
model.root.input.solver.sections.section_times = [0.0, 1200, 4200, 5400]   
# model.root.input.solver.sections.section_continuity = [0, 0, 0]
model.root.input.model.unit_000.sec_000.const_coeff = [1575.0, 0.1] 
model.root.input.model.unit_000.sec_001.const_coeff = [1575.0, 0.0] 
model.root.input.model.unit_000.sec_001.lin_coeff = [-0.5,0.0]
model.root.input.model.unit_000.sec_002.const_coeff = [75,0.0] 


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

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

# Number of cores for parallel simulation
model.root.input.solver.nthreads = 6

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

# Return data
model.root.input['return'].split_components_data = 0
model.root.input['return'].split_ports_data = 0
model.root.input['return'].unit_000.write_solution_bulk = 1
model.root.input['return'].unit_000.write_solution_inlet = 1
model.root.input['return'].unit_000.write_solution_outlet = 1

# Copy settings to the other unit operations
model.root.input['return'].unit_001 = model.root.input['return'].unit_000
model.root.input['return'].unit_002 = model.root.input['return'].unit_000

# Solution times
model.root.input.solver.user_solution_times = np.linspace(0, 5400, 5401)

model.filename = 'model.h5'
model.save()

data = model.run()

if data.returncode == 0:
    print("Simulation completed successfully")
    model.load()
else:
    print(data)
    raise Exception("Simulation failed")

plt.figure()

time = model.root.output.solution.solution_times
c = model.root.output.solution.unit_001.solution_outlet
plt.plot(time/60, c)
plt.xlabel('$time~/~min$')
plt.ylabel('$Outlet~concentration~/~mol \cdot m^{-3} $')
plt.show()

If anyone could give me some advice, I would greatly appreciate it!
Thanks,
Yang

[ronald.jaepel]

A Ks of 4.3 is really high. If that was taken from another publication, it might have been used with concentrations of mol/L, not mol/m³ as is the case in CADET. Dividing it by 1000 results in a nice chromatogram.

(A small side note: your code to generate the figure used c = model.root.output.solution.unit_001.solution_outlet, which returns a 2D array of the salt concentration and the protein concentration. Plotting that on one y-axis leads to a very squished protein concentration curve, I’ve appended the code to generate the figure I made)

fig, ax = plt.subplots(1, 1)
ax2 = ax.twinx()

ax.plot(time / 60, c[:, 1])
ax2.plot(time / 60, c[:, 0], ":", color="grey")

ax.set_xlabel('$time~/~min$')
ax.set_ylabel('$Outlet~concentration~/~mol \cdot m^{-3} $')
ax2.set_ylabel('$Outlet~Salt~concentration~/~mol \cdot m^{-3} $')
fig.tight_layout()
plt.show()

[Yang]

Thanks for your reply! The simulation runs now. I am very grateful for your help and advice.