Error in dextrane puls simulation

CADET Troubleshooting
1

Hello All,
I was trying to simualte injection of dextrane to a column like Dextran pulse excise. I get an error like this.

“CompletedProcess(args=[‘c:/cadet/bin/cadet-cli.exe’, ‘model.h5’], returncode=3221225477, stdout=b’‘, stderr=b’')”

The errror message is simple. Anyone has some hints for me? Thanks !

import numpy as np
import matplotlib.pyplot as plt

import numpy as np
import matplotlib.pyplot as plt
from cadet import Cadet
Cadet.cadet_path = 'c:/cadet/bin/cadet-cli.exe'
model = Cadet()
model.root.input.model.nunits = 3
model.root.input.model.unit_000.unit_type = 'INLET'
model.root.input.model.unit_000.ncomp = 1
model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
model.root.input.model.unit_001.unit_type = 'LUMPED_RATE_MODEL_WITH_PORES'
model.root.input.model.unit_001.ncomp = 1
​
# Geometry
model.root.input.model.unit_001.col_length = 0.6                # col len in m
model.root.input.model.unit_001.cross_section_area = 1.e-4      # CX section area in m2
model.root.input.model.unit_001.col_porosity = 0.37             # col porosity 
model.root.input.model.unit_001.par_porosity = 0.33             # particle or bead porosity 
model.root.input.model.unit_001.par_radius = 4.5e-5             # particle or bead radius 
                                                                
## Transport
model.root.input.model.unit_001.col_dispersion = 2.0e-7          # m^2 / s (interstitial volume)
model.root.input.model.unit_001.film_diffusion = [0.0]           # m / s
​
model.root.input.model.unit_001.init_c = [0.0,]
model.root.input.model.unit_002.unit_type = 'OUTLET'
model.root.input.model.unit_002.ncomp = 1
### Grid cells
model.root.input.model.unit_001.discretization.ncol = 100
model.root.input.model.unit_001.discretization.npar = 5
model.root.input.model.unit_001.discretization.npartype = 1
​
### Bound states
model.root.input.model.unit_001.discretization.nbound = 3*[0]
​
### 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.solver.sections.nsec = 2
model.root.input.solver.sections.section_times = [0.0, 50, 1200,]   # s
model.root.input.solver.sections.section_continuity = []
## Inlet Profile
model.root.input.model.unit_000.sec_000.const_coeff = [1.0,] # mol / m^3
​
## Switches
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, 1e-7,  # [unit_000, unit_001, all components, all components, Q/ m^3*s^-1 
    1, 2, -1, -1, 1e-7]  # [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 = 1
​
# 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 = True
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
#set the times that the simulator writes out data for
model.root.input.solver.user_solution_times = np.linspace(0, 1200, 1001)
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")
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")
CompletedProcess(args=['c:/cadet/bin/cadet-cli.exe', 'model.h5'], returncode=3221225477, stdout=b'', stderr=b'')
---------------------------------------------------------------------------
Exception                                 Traceback (most recent call last)
~\AppData\Local\Temp\ipykernel_11896\2242269908.py in <module>
      9 else:
     10     print(data)
---> 11     raise Exception("Simulation failed")

Exception: Simulation failed
2

Could you please try again with:

model.root.input.solver.sections.section_continuity = [0, 0]

If that fixes the issue, we should definitely improve the error message.

3

@j.schmoelder Thank you so much! It fixed the run error.
I found another error that prevents returning data or plotting data in my way.

model.root.input['return'].split_components_data = True

It should be False or 0 to fix the problem in my way. If you have better way, please let us know.

Here is the working version. I know it is not the best version. I feel it is good for beginner.


Code
import numpy as np
import matplotlib.pyplot as plt
from cadet import Cadet
Cadet.cadet_path = 'c:/cadet/bin/cadet-cli.exe'
import numpy as np
import matplotlib.pyplot as plt
from cadet import Cadet
Cadet.cadet_path = 'c:/cadet/bin/cadet-cli.exe'
model = Cadet()
model.root.input.model.nunits = 3
model.root.input.model.unit_000.unit_type = 'INLET'
model.root.input.model.unit_000.ncomp = 1
model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
6
model.root.input.model.unit_001.unit_type = 'LUMPED_RATE_MODEL_WITH_PORES'
model.root.input.model.unit_001.ncomp = 1
​
# Geometry
model.root.input.model.unit_001.col_length = 0.3               # col len in m
model.root.input.model.unit_001.cross_section_area = 1.e-4      # CX section area in m2
model.root.input.model.unit_001.col_porosity = 0.37             # col porosity 
model.root.input.model.unit_001.par_porosity = 0.63             # particle or bead porosity 
model.root.input.model.unit_001.par_radius = 4.5e-5             # particle or bead radius 
                                                                
## Transport
model.root.input.model.unit_001.col_dispersion = 2.0e-7          # m^2 / s (interstitial volume)
model.root.input.model.unit_001.film_diffusion = [0.0]           # m / s
​
model.root.input.model.unit_001.init_c = 1*[0.0,]
0
### Grid cells
model.root.input.model.unit_001.discretization.ncol = 100
model.root.input.model.unit_001.discretization.npar = 5
model.root.input.model.unit_001.discretization.npartype = 1
​
### Bound states
model.root.input.model.unit_001.discretization.nbound = [0,]   # Number of bound states for each component. This is linked to n_comp. 
​
### 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 = 1
2
model.root.input.solver.sections.nsec = 2
model.root.input.solver.sections.section_times = [0.0, 20, 1200,]   # s
model.root.input.solver.sections.section_continuity = [0, 0]
## Inlet Profile
model.root.input.model.unit_000.sec_000.const_coeff = [1.0,] # mol / m^3
model.root.input.model.unit_000.sec_001.const_coeff = [0.0,]
## Switches
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, 1e-7,  # [unit_000, unit_001, all components, all components, Q/ m^3*s^-1 
    1, 2, -1, -1, 1e-7]  # [unit_001, unit_002, all components, all components, Q/ m^3*s^-1 
## Switches
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, 1e-7,  # [unit_000, unit_001, all components, all components, Q/ m^3*s^-1 
    1, 2, -1, -1, 1e-7]  # [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 = 1
​
# 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 # True (1) or false (0). Determines whether a joint dataset (matrix or tensor) for all components is created or if each component is put in a separate dataset
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
#set the times that the simulator writes out data for
model.root.input.solver.user_solution_times = np.linspace(0, 1200, 1001)
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")
Simulation completed successfully
plt.figure()

time = model.root.output.solution.solution_times
c = model.root.output.solution.unit_001.solution_outlet

plt.plot(time, c)
plt.title('Column (Outlet)')
plt.xlabel('$time~/~s$')
plt.ylabel('$concentration~/~mM$')
plt.show()
plt.figure()
​
time = model.root.output.solution.solution_times
c = model.root.output.solution.unit_001.solution_outlet
​
plt.plot(time, c)
plt.title('Column (Outlet)')
plt.xlabel('$time~/~s$')
plt.ylabel('$concentration~/~mM$')
plt.show()

Again, thanks a lot! And you can lose the thread.

4

That’s actually optional (and I personally also prefer setting it to False). If it’s True, it will create a separate key for each component in the output. E.g.

print(model.root.output.solution.unit_001.keys())
Out[13]: dict_keys(['solution_bulk', 'solution_inlet_comp_000', 'solution_outlet_comp_000'])

So you would also have to read the solution separately for each component. E.g.

time = model.root.output.solution.solution_times
solution_0 = model.root.output.solution.unit_001.solution_outlet_comp_000

fig, ax = plt.subplots()
ax.plot(time, solution_0, label='Component 0')
5

@j.schmoelder Thanks so much for the comments. You comfirmed my speculation. I suspect the data saved in a different format if I use TRUE. I have error to plot if I use my way. If I use True, I need use alternative way to pull out the data points for Y-axis as you pointed out.
I tested your way by using True and reading data points from separated componment.
Again, thanks for helping me to undrstand this.

Here is the alternative version with splitting component data true.

import numpy as np
import matplotlib.pyplot as plt
from cadet import Cadet
Cadet.cadet_path = 'c:/cadet/bin/cadet-cli.exe'
model = Cadet()
model.root.input.model.nunits = 3
model.root.input.model.unit_000.unit_type = 'INLET'
model.root.input.model.unit_000.ncomp = 1
model.root.input.model.unit_000.inlet_type = 'PIECEWISE_CUBIC_POLY'
6
model.root.input.model.unit_001.unit_type = 'LUMPED_RATE_MODEL_WITH_PORES'
model.root.input.model.unit_001.ncomp = 1
​
# Geometry
model.root.input.model.unit_001.col_length = 0.3               # col len in m
model.root.input.model.unit_001.cross_section_area = 1.e-4      # CX section area in m2
model.root.input.model.unit_001.col_porosity = 0.37             # col porosity 
model.root.input.model.unit_001.par_porosity = 0.63             # particle or bead porosity 
model.root.input.model.unit_001.par_radius = 4.5e-5             # particle or bead radius 
                                                                
## Transport
model.root.input.model.unit_001.col_dispersion = 2.0e-7          # m^2 / s (interstitial volume)
model.root.input.model.unit_001.film_diffusion = [0.0]           # m / s
​
model.root.input.model.unit_001.init_c = 1*[0.0,]
### Grid cells
model.root.input.model.unit_001.discretization.ncol = 100
model.root.input.model.unit_001.discretization.npar = 5
model.root.input.model.unit_001.discretization.npartype = 1

### Bound states
model.root.input.model.unit_001.discretization.nbound = [0,]   # Number of bound states for each component. This is linked to n_comp. 

### 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
### Grid cells
model.root.input.model.unit_001.discretization.ncol = 100
model.root.input.model.unit_001.discretization.npar = 5
model.root.input.model.unit_001.discretization.npartype = 1
​
### Bound states
model.root.input.model.unit_001.discretization.nbound = [0,]   # Number of bound states for each component. This is linked to n_comp. 
​
### 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 = 1
2
model.root.input.solver.sections.nsec = 2
model.root.input.solver.sections.section_times = [0.0, 20, 1200,]   # s
model.root.input.solver.sections.section_continuity = [0, 0]
## Inlet Profile
model.root.input.model.unit_000.sec_000.const_coeff = [1.0,] # mol / m^3
model.root.input.model.unit_000.sec_001.const_coeff = [0.0,]
## Switches
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, 1e-7,  # [unit_000, unit_001, all components, all components, Q/ m^3*s^-1 
    1, 2, -1, -1, 1e-7]  # [unit_001, unit_002, all components, all components, Q/ m^3*s^-1 
## Switches
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, 1e-7,  # [unit_000, unit_001, all components, all components, Q/ m^3*s^-1 
    1, 2, -1, -1, 1e-7]  # [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 = 1

# 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
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 = 1
​
# 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 = 1 # True (1) or false (0). Determines whether a joint dataset (matrix or tensor) for all components is created or if each component is put in a separate dataset
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
#set the times that the simulator writes out data for
model.root.input.solver.user_solution_times = np.linspace(0, 1200, 1001)
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")
Simulation completed successfully
plt.figure()

time = model.root.output.solution.solution_times
solution_0 = model.root.output.solution.unit_001.solution_outlet_comp_000

plt.plot(time, solution_0)
plt.title('Column (Outlet)')
plt.xlabel('$time~/~s$')
plt.ylabel('$concentration~/~mM$')
plt.show()
1 Like