Section dependent isotherms and alternative approaches

CADET Troubleshooting
1

I wonder if a time section dependent isotherm is possible in CADET? I wanted to simulate a situation where either the isotherm parameters change or the isotherm model changes during one single simulation (experiment).

Time section dependent isotherm

After looking at the CADET document, I found the inlet unit and transport parameters can be section dependent, but not adsorption models. This is not supported now.

Run simulations sequentially, use the last state of a previous simulation to initialize the next

This should be possible. Relevant fields to check are: write_solution_last_unit (or write_solution_last). This will write the state vectors at the last time point to last_state_y in the output group and time derivatives to last_state_ydot. Appending last_state_ydot to last_state_y gives the final vector that can be passed to the next simulation through the init_state in the input group.

Using external functions

A successful example is posted below. This can change isotherm parameters but not the model itself.

1 Like
2

Unfortunately the second approach did not work properly. See my test below.

original simulation

GMR + Multicomponent Langmuir. Single component, breakthrough.

image
image717Ă—527 28.7 KB

Break up this simulation to two parts: t \in (0,8000) and t \in (8000, 16000).

simulation 1

image
image726Ă—525 30.4 KB

set up init_state for simulation 2

Append last_state_y and last_state_ydot and pass it to init_state.

simulation 2

The same script in simulation 1 is used. Nothing changes except for the init_state. The results are:

image
image731Ă—532 24.7 KB

In theory I should get a horizontal line. This result would indicate something is wrong.

Troubleshooting

The last state of simulation 1:

image
image751Ă—623 23.4 KB

This plots q and c_p as well as c as a function of particle radius r at an axial position at a fixed time. Note that c is a constant at an axial position and is not dependent on r. To facilitate comparison I plotted it as a straight horizontal line here.

The initial state of simulation 2 (at the same axial position):

image
image710Ă—611 22.5 KB

The magnitude is already different. Upon a closer look at the outputted last_state_y vector in simulation 1:

image
image646Ă—107 3.32 KB

It already contains these large values that are not seen in the simulation 1 which I think caused the problem.

Is there any conversion between the output last_state_y and the input init_state to properly set up the simulations?

Thanks,
Flynn.

1 Like
3

In your particular example (breaking up the simulation in parts with same isotherm), I suggest to disable consistent initialization for the second part:

model.root.input.solver.consistent_init_mode = 0 # or 3

My theory is that consistent initialization somehow fails and produces those unreasonable values.

Another way to accomplish your goal (change isotherm parameters during simulation) are external functions. Here, the isotherm parameters are taken as

k_{a,i} = k_{a,i}(T) = k_{a,i,0} + k_{a,i,1} T + k_{a,i,2} T^2 + k_{a,i,3} T^3,

where T = T(t) is some profile provided to the simulator. So each isotherm parameter is a third-degree polynomial of a given profile. You could set the constant coefficients of the isotherm parameters to your first configuration and the linear coefficients to

second configuration - first configuration.

For the first part of the simulation, we need to set T = 0 and for the second, we change to T=1. Here are the docs for specifying a profile, where you’d set VELOCITY to inf: A profile is usually moved along the column over time. In your case, we want the profile to instantly apply to all spatial positions (solely depend on time).

We have some posts in the forum about external functions.

3 Likes
4

Thanks for the reply!

  1. Disabling the consistent initialization unfortunately did not work. IDAS would fail at t=0, which is kind of strange since the last_state_y and last_state_ydot from simulation 1 should already be consistent? I think the problem may have something to do with the last_state_y vector from simulation 1, it already contained some unreal values which I was not able to find at any t, r and z.

  2. Using external functions seems to be working. My test case is simple: GRM, Multi-component Langmuir

    t \in (0, 4000): q_{max}=21, \; k_a=100
    t \in (4000, 8000): q_{max}=23, \; k_a=110.

    A runnable Python code is below:

from cadet import Cadet
import numpy as np
import matplotlib.pyplot as plt

# replace the path
Cadet.cadet_path = r"C:\Users\zwend\CADET\cadetL2\bin\cadet-cli.exe" 

## set up
c_feed = 0.0113                                            

def create_model():
    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 = 'GENERAL_RATE_MODEL'
    model.root.input.model.unit_001.ncomp = 1

    ## Geometry
    model.root.input.model.unit_001.col_length = 0.05                         # m              
    model.root.input.model.unit_001.cross_section_area = 5.025e-6/0.05      # m^2
    model.root.input.model.unit_001.col_porosity = 0.37                     # 1      
    model.root.input.model.unit_001.par_porosity = 0.97                     # 1      
    model.root.input.model.unit_001.par_radius = 0.045e-3                              # m     

    ## Transport
    model.root.input.model.unit_001.col_dispersion = 2.061e-8                          # m^2/s       
    model.root.input.model.unit_001.film_diffusion = [1.18e-5]                         # m/s 
    model.root.input.model.unit_001.par_diffusion = [0.8e-11,]                              # m^2/s  [0.55e-11]
    model.root.input.model.unit_001.par_surfdiffusion = [0.0, ]    

    ## external adsorption isotherm
    model.root.input.model.unit_001.adsorption_model = 'EXT_MULTI_COMPONENT_LANGMUIR'
    model.root.input.model.unit_001.adsorption.extfun = 0
    
    model.root.input.model.unit_001.adsorption.is_kinetic = 1
    
    model.root.input.model.unit_001.adsorption.ext_mcl_ka =     100.0         
    model.root.input.model.unit_001.adsorption.ext_mcl_ka_t =   10.0
    model.root.input.model.unit_001.adsorption.ext_mcl_ka_tt =  0
    model.root.input.model.unit_001.adsorption.ext_mcl_ka_ttt = 0
    
    model.root.input.model.unit_001.adsorption.ext_mcl_kd =     1.0                                           ## s^-1
    model.root.input.model.unit_001.adsorption.ext_mcl_kd_t =   0
    model.root.input.model.unit_001.adsorption.ext_mcl_kd_tt =  0
    model.root.input.model.unit_001.adsorption.ext_mcl_kd_ttt = 0
    
    model.root.input.model.unit_001.adsorption.ext_mcl_qmax =     21.0        
    model.root.input.model.unit_001.adsorption.ext_mcl_qmax_t =   2.0
    model.root.input.model.unit_001.adsorption.ext_mcl_qmax_tt =  0
    model.root.input.model.unit_001.adsorption.ext_mcl_qmax_ttt = 0
    
    ## external function
    model.root.input.model.external.source_000.extfun_type = 'PIECEWISE_CUBIC_POLY'
    model.root.input.model.external.source_000.section_times = [0.0, 4000.0, 8000.0]
    model.root.input.model.external.source_000.velocity = 1e5                                                 # set it to a large number
    model.root.input.model.external.source_000.const_coeff = [0.0, 1.0]
    model.root.input.model.external.source_000.lin_coeff =   [0.0, 0.0]
    model.root.input.model.external.source_000.quad_coeff =  [0.0, 0.0]
    model.root.input.model.external.source_000.cube_coeff =  [0.0, 0.0]
    
    ## init condition
    model.root.input.model.unit_001.init_c = [0.0, ]
    model.root.input.model.unit_001.init_q = [0.0, ] 

    ### Grid cells
    model.root.input.model.unit_001.discretization.ncol = 50                    ## increase as needed
    model.root.input.model.unit_001.discretization.npar = 12                    ## same as above

    ### Bound states
    model.root.input.model.unit_001.discretization.nbound = [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 = 2

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

    model.root.input.solver.sections.nsec = 1
    model.root.input.solver.sections.section_times = [0.0, 8000]   # s
    model.root.input.solver.sections.section_continuity = [0,0]

    model.root.input.model.unit_000.sec_000.const_coeff = [c_feed]                                                        # mol / m^3;  mg/ml = kg/m^3;  1 kda = 1 kg/mol 1.25/150
    model.root.input.model.unit_000.sec_000.lin_coeff = [0.0]
    model.root.input.model.unit_000.sec_000.quad_coeff = [0.0]
    model.root.input.model.unit_000.sec_000.cube_coeff = [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, 2.5 /(6e7),
        1, 2, -1, -1, 2.5 /(6e7)]

    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
    model.root.input['return'].split_ports_data = 0
    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, 8000, 1000+1)
    
    return model

model = create_model()
model.filename = 'model.h5'
model.save()
data = model.run()
print(data)
model.load()   

time = model.root.output.solution.solution_times
c = model.root.output.solution.unit_001.solution_outlet_comp_000                                        ## 1 mol/m^3 * 150 kg/mol=150 kg/m^3 = 150 mg/mL
    
fig = plt.figure()
ax = fig.add_subplot()
ax.plot(time, c*150, c="orange") 
ax.set_xlabel(r't/s')
ax.set_ylabel(r'conc./ mg/mL')
plt.show()

The warning would still appear: Warning: setExternalFunctions::113] Index 0 exceeds number of passed external functions (0), external dependence is ignored\r\n, but it should be fine. The following plot would appear:

image
image717Ă—535 32.4 KB

If there is no change in the isotherm parameters:
image
image718Ă—530 30 KB

The sharp decrease at t=4000 made sense because q_{max} and k_a increased, the column continued to absorb the feed so the outlet concentration decreased.

3 Likes
Implementing Section Dependent Isotherms
5

Please ignore this warning. This is a known “problem”.

Thanks for sharing the solution!

6

This is actually something @j.breuer and I were discussing recently. I need to follow this exact approach to solve a dilemma described in the earlier post (Protein at inlet not showing up at outlet - #10 by Flynn). We were not able to figure out a solution (Jan mentioned there might be a problem with consistent initialization). I still have those hdf5 files lying around, so we could look into actually solving the problem again. To be honest, it was a highly atypical application of CADET and using a workaround with running simulations sequentially is totally fine. The ultimate outlet solution already requires three separate simulations, before we implement the workaround.

I have another application on the docket that also requires using section-dependent isotherm parameters, so I am interested in the outcome of your investigation @Flynn.

1 Like
7

Hey Flynn,

How’s it going? Did you get the second approach (sequential simulations) working? I need to finish up a paper and use sequential sims (can’t think of any other way to solve my problem). So, I have a vested interest in getting this to work.

I think the problem you described in this post is much simpler than what I am trying to do, so it’s a good test case for this approach. If it’s not working, do we have any plans for fixing this? @Jan also suggested to me to use the second approach and I am curious to hear his thoughts on this.

8

Hey Scott, unfortunately I did try to debug this but could not find the root cause. Deeper debugging will probably take me some time. Maybe Sam @s.leweke has some ideas or suggestions on this?

9

Just bumping this thread to see if maybe we can take a stab at fixing it :slight_smile:

10

@j.schmoelder @Flynn This is the thread I was referring to about consistent initialization

11

For some reason the init_state defined at the unit operation level is still not working, maybe it has something to do with the state structure, but I found that init_state_y and init_state_ydot defined for the entire system (or model) actually worked.

To achieve this, one can first set write_solution_last to true in the first simulation. Then in the output group one should be able to find .root.output.last_state_y and .root.output.last_state_ydot. For the second simulation, one can pass both to .root.input.model.init_state_y and .root.input.model.init_state_ydot, separately. This should work.

In the following I break up a full breakthrough curve using GRM with a single-component Langmuir isotherm.
Part I:

image
image696Ă—531 22.2 KB

Part II:

image
image716Ă—523 27.6 KB

The full result is:

image
image691Ă—533 22.2 KB

One can also modify the isotherm parameters for the second simulation.

I can update the first wiki if admin gives me access.

I can also provide an example written in cadet-core.

2 Likes
12

Thanks for sharing, Flynn. I will definitely try this for my problem. Can you share your cadet core example?

Implementing Section Dependent Isotherms
13

Here is a working example

from cadet import Cadet
import numpy as np
import matplotlib.pyplot as plt

Cadet.cadet_path =   ## add your path

n_sites = 1

Q = 2.5 / (6e7)                                           # volumetric flow rate m^3/s

## kd, eta, pka
kd = np.array([1,])

keq = np.array([10,])
qmax = np.array([62,]) 

## set up
c_feed = 8.0 / 150                                       # mol/m^3 (devide mg/mL by 150)
simulation_end_time1 = 1200                              

column_length = 0.05  # m
column_volume = 5.025e-6 # m^3

column_porosity = 0.37
particle_porosity = 0.97

protein_MW = 150

total_porosity = column_porosity + (1.0-column_porosity)*particle_porosity


#### part I simulation
def create_model(qmax, keq, Dp):
    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 = 'GENERAL_RATE_MODEL'
    model.root.input.model.unit_001.ncomp = 1

    ## Geometry
    model.root.input.model.unit_001.col_length = column_length                         # m              
    model.root.input.model.unit_001.cross_section_area = column_volume/column_length   # m^2
    model.root.input.model.unit_001.col_porosity = column_porosity                     # 1      
    model.root.input.model.unit_001.par_porosity = particle_porosity                   # 1      
    model.root.input.model.unit_001.par_radius = 0.045e-3                              # m     

    ## Transport
    model.root.input.model.unit_001.col_dispersion = 2.061e-8                          # m^2/s       
    model.root.input.model.unit_001.film_diffusion = [1.18e-5]                         # m/s 
    model.root.input.model.unit_001.par_diffusion = [Dp,]                         # m^2/s  [0.55e-11]
    model.root.input.model.unit_001.par_surfdiffusion = n_sites * [0.0, ]    

    model.root.input.model.unit_001.adsorption_model = 'MULTI_COMPONENT_LANGMUIR'

    model.root.input.model.unit_001.adsorption.is_kinetic = 0
    model.root.input.model.unit_001.adsorption.mcl_ka = protein_MW*(1.0-total_porosity) * keq         ##  m^3 solid phase / mol protein / s   ml/mg=m^3/kg    1 kda = 1 kg/mol
    model.root.input.model.unit_001.adsorption.mcl_kd = kd                                          ## s^-1
    model.root.input.model.unit_001.adsorption.mcl_qmax = qmax/protein_MW/(1.0-total_porosity)   ##  mol/m^3 solid phase mg/ml=kg/m^3 / 150 kg/mol = 1/150 mol/m^3

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

    ### Grid cells
    model.root.input.model.unit_001.discretization.ncol = 50
    model.root.input.model.unit_001.discretization.npar = 12

    ### Bound states
    model.root.input.model.unit_001.discretization.nbound = [n_sites]

    ### 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 = 2

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

    model.root.input.solver.sections.nsec = 1
    model.root.input.solver.sections.section_times = [0.0, simulation_end_time1, simulation_end_time1+1]   # s
    model.root.input.solver.sections.section_continuity = [0,0]

    model.root.input.model.unit_000.sec_000.const_coeff = [c_feed]                                                        # mol / m^3;  mg/ml = kg/m^3;  1 kda = 1 kg/mol 1.25/150
    model.root.input.model.unit_000.sec_000.lin_coeff = [0.0]
    model.root.input.model.unit_000.sec_000.quad_coeff = [0.0]
    model.root.input.model.unit_000.sec_000.cube_coeff = [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, Q,
        1, 2, -1, -1, Q]

    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
    model.root.input['return'].split_ports_data = 0
    model.root.input['return'].unit_000.write_solution_outlet = 1
    model.root.input['return'].write_solution_last = 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, simulation_end_time1, 1*simulation_end_time1+1)
    
    return model

model = create_model(61.9, 473.8, 8.283e-12)
model.filename = 'model.h5'
model.save()
data = model.run()
model.load()   

time = model.root.output.solution.solution_times
c = model.root.output.solution.unit_001.solution_outlet_comp_000                                        ## 1 mol/m^3 * 150 kg/mol=150 kg/m^3 = 150 mg/mL
    
fig = plt.figure()
ax = fig.add_subplot()
ax.plot(time, c*150, c="orange", label="Simulated")
ax.set_xlabel(r't/s')
ax.set_ylabel(r'conc.')
plt.show()

init_s = model.root.output.last_state_y
init_sdot = model.root.output.last_state_ydot


#### part II simulation
def create_model2(qmax, keq, Dp):
    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 = 'GENERAL_RATE_MODEL'
    model.root.input.model.unit_001.ncomp = 1

    ## Geometry
    model.root.input.model.unit_001.col_length = column_length                         # m              
    model.root.input.model.unit_001.cross_section_area = column_volume/column_length   # m^2
    model.root.input.model.unit_001.col_porosity = column_porosity                     # 1      
    model.root.input.model.unit_001.par_porosity = particle_porosity                   # 1      
    model.root.input.model.unit_001.par_radius = 0.045e-3                              # m     

    ## Transport
    model.root.input.model.unit_001.col_dispersion = 2.061e-8                          # m^2/s       
    model.root.input.model.unit_001.film_diffusion = [1.18e-5]                         # m/s 
    model.root.input.model.unit_001.par_diffusion = [Dp,]                              # m^2/s
    model.root.input.model.unit_001.par_surfdiffusion = n_sites * [0.0, ]    

    model.root.input.model.unit_001.adsorption_model = 'MULTI_COMPONENT_LANGMUIR'

    model.root.input.model.unit_001.adsorption.is_kinetic = 0
    model.root.input.model.unit_001.adsorption.mcl_ka = protein_MW*(1.0-total_porosity) * keq       ##  m^3 solid phase / mol protein / s   ml/mg=m^3/kg    1 kda = 1 kg/mol
    model.root.input.model.unit_001.adsorption.mcl_kd = kd                                          ## s^-1
    model.root.input.model.unit_001.adsorption.mcl_qmax = qmax/protein_MW/(1.0-total_porosity)      ##  mol/m^3 solid phase mg/ml=kg/m^3 / 150 kg/mol = 1/150 mol/m^3

    model.root.input.model.init_state_y = init_s
    model.root.input.model.init_state_ydot = init_sdot

    ### Grid cells
    model.root.input.model.unit_001.discretization.ncol = 50
    model.root.input.model.unit_001.discretization.npar = 12

    ### Bound states
    model.root.input.model.unit_001.discretization.nbound = [n_sites]

    ### 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 = 2

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

    model.root.input.solver.sections.nsec = 1
    model.root.input.solver.sections.section_times = [0, 1.0*simulation_end_time1+1]   # s
    model.root.input.solver.sections.section_continuity = [0,0]

    model.root.input.model.unit_000.sec_000.const_coeff = [c_feed]                                                        # mol / m^3;  mg/ml = kg/m^3;  1 kda = 1 kg/mol 1.25/150
    model.root.input.model.unit_000.sec_000.lin_coeff = [0.0]
    model.root.input.model.unit_000.sec_000.quad_coeff = [0.0]
    model.root.input.model.unit_000.sec_000.cube_coeff = [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, Q,
        1, 2, -1, -1, Q]

    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
    model.root.input['return'].split_ports_data = 0
    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, 1.0*simulation_end_time1+1, simulation_end_time1+1)
    
    return model

model = create_model2(61.9, 473.8, 8.283e-12)
model.filename = 'practice1.h5'
model.save()
data = model.run()
model.load()

time = model.root.output.solution.solution_times
c = model.root.output.solution.unit_001.solution_outlet_comp_000
    
fig = plt.figure()
ax = fig.add_subplot()
ax.plot(time, c*150, c="orange", label="Simulated")
ax.set_xlabel(r't/s')
ax.set_ylabel(r'conc.')
plt.show()
2 Likes