search

sys-bio / roadrunner

libRoadRunner: A high-performance SBML simulator
http://libroadrunner.org/
Other
39 stars 24 forks source link

Onestep simulations result in CVODE Error: CV_CONV_FAILURE #1232

Open matthiaskoenig opened 1 month ago

matthiaskoenig commented 1 month ago

Hi all, I see the following issue a lot. If I run the identical simulation to a simulate via a onestep simulation the integrator fails with CV_CONV_FAILURE.

It is the identical simulation and model just once run via simulate, the other time via onestep. The only thing I can imagine is that roadrunner/CVODE is using different simulator/integrator settings for the two cases (e.g. retries, minimal step size based on end time, ...). I see this in many models and the iterative one step fails very often with convergence failures.

We use this a lot in model coupling within PDEs where we just continue the simulation for one step.

Example attached for the following model: droplets_iri_liver.zip

"""Example script for porous media coupling for IRI & Droplet model."""
from pathlib import Path
import pandas as pd
import roadrunner

# -------------------------------------------------------------------------------------
# load model
# -------------------------------------------------------------------------------------
model_path = Path(__file__).parent /  "droplets_iri_liver.xml"

# loading the model in roadrunner
r = roadrunner.RoadRunner(str(model_path))

# -------------------------------------------------------------------------------------
# set changes for given simulator
# -------------------------------------------------------------------------------------
vol_point = 1.0  # [l]
f_fluid = 0.3  # [dimensionless] fluid fraction
f_solid = 1 - f_fluid  # [dimensionless] solid fraction

# These are the changes which have to be applied
changes = {
    # setting volumes
    "Vext": vol_point * f_fluid,  # [liter] fluid phase volume
    "Vli": vol_point * f_solid,  # [liter] solid phase volume

    # concentrations of external concentrations in fluid phase
    "[glc_ext]": 5.0,  # [mM]
    "[lac_ext]": 1.0,  # [mM]
    "[o2_ext]": 1.2,  # [mM]
    "[alt_ext]": 0.0,  # [mM]
    "[ast_ext]": 0.0,  # [mM]
    "[ffa_ext]": 0.0,  # [mM]
    "[vldl_ext]": 0.0,  # [mM]

    "[glc]": 5.0,  # [mM]
    "[lac]": 1.0,  # [mM]
    "[o2]": 1.2,  # [mM]
    "[alt]": 1.0,  # [mM]
    "[ast]": 0.7,  # [mM]
    "[pyr]": 0.8,  # [mM]
    "[atp]": 3.0,  # [mM]
    "[adp]": 1.0,  # [mM]
    "[nadh]": 3.9,  # [mM]
    "[nad]": 0.1,  # [mM]
    "[ros]": 0.0,  # [mM]

    # set temperature
    "temperature": 37 + 273.15,  # [Kelvin] warm simulation
    # "temperature": 4 + 273.15,  # [Kelvin] cold simulation
}

# -------------------------------------------------------------------------------------
# Timecourse
# -------------------------------------------------------------------------------------
# run simulation
tend = 600  # [min]
steps = 100
dt = tend/steps  # [min]

# single simulation (singel simulation works)
print("--- simulate --- ")
r.resetAll()
for key, value in changes.items():
    r.setValue(key, value)
s1 = r.simulate(start=0.0, end=tend, steps=steps) #*60
df1 = pd.DataFrame(s1, columns=s1.colnames)
print(df1.head())

# FEBio version (onestep)
print("--- onestep --- ")
r.resetAll()
for key, value in changes.items():
    r.setValue(key, value)

t_model = 0.0
while t_model < tend:
    s = r.oneStep(t_model, dt)
    t_model += dt
df2 = pd.DataFrame(s, columns=s.colnames)
print(df2.head)

Results in 

/home/mkoenig/.virtualenvs/pkdb_models_new/bin/python /home/mkoenig/git/pkdb_models/pkdb_models/models/tpm_fat/simulations/iri_droplets_liver_rr.py 
--- simulate --- 
   time     [ins]     [ffa]  [vldl_ext]      [dg]  [protein_dgat2]  [protein_ldsyn]       [tg1]       [tg2]      [tg3]     [tg4]  [glc_ext]  ...  [ast_ext]     [pyr]     [atp]     [adp]    [nadh]     [nad]  [ros]  [alt]  [ast]  [necrosis]       [tgd]  tsim
0   0.0  0.000001  2.142857    0.000000  1.714286         1.428571         1.428571    0.000000    0.000000   0.000000  0.000000   5.000000  ...        0.0  0.800000  3.000000  1.000000  3.900000  0.100000    0.0    1.0    0.7         0.0    0.000000   0.0
1   6.0  0.000001  0.000449    0.020722  1.726787         1.428571         1.428571  177.988176   10.942289   0.261533  0.005179   4.908386  ...        0.0  0.729103  3.038536  0.961464  3.989563  0.010437    0.0    1.0    0.7         0.0  189.197177   6.0
2  12.0  0.000001  0.000449    0.022024  1.727565         1.428571         1.428571  270.007558   37.756799   1.758146  0.044349   4.820259  ...        0.0  0.752057  3.060534  0.939466  3.989774  0.010226    0.0    1.0    0.7         0.0  309.566852  12.0
3  18.0  0.000001  0.000449    0.023327  1.728343         1.428571         1.428571  309.703061   70.507305   6.067988  0.181956   4.734098  ...        0.0  0.774237  3.070634  0.929366  3.989982  0.010018    0.0    1.0    0.7         0.0  386.460310  18.0
4  24.0  0.000001  0.000449    0.024630  1.729120         1.428571         1.428571  317.444936  102.337022  14.189627  0.545982   4.648687  ...        0.0  0.796205  3.071470  0.928530  3.990204  0.009796    0.0    1.0    0.7         0.0  434.517567  24.0

[5 rows x 30 columns]
--- onestep --- 
Error: CVODE Error: CV_CONV_FAILURE, Module: CVODES, Function: CVode, Message: At t = 180.806 and h = 4.48854e-07, the corrector convergence test failed repeatedly or with |h| = hmin.
Traceback (most recent call last):
  File "/home/mkoenig/git/pkdb_models/pkdb_models/models/tpm_fat/simulations/iri_droplets_liver_rr.py", line 78, in <module>
    s = r.oneStep(t_model, dt)
        ^^^^^^^^^^^^^^^^^^^^^^
  File "/home/mkoenig/.virtualenvs/pkdb_models_new/lib/python3.11/site-packages/roadrunner/roadrunner.py", line 2968, in oneStep
    return _roadrunner.RoadRunner_oneStep(self, currentTime, stepSize, reset)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: CVODE Error: CV_CONV_FAILURE: Convergence test failures occurred too many times (= MXNCF = 10) during one internal timestep or occurred with |h| = hmin.; In virtual double rr::CVODEIntegrator::integrate(double, double)

Process finished with exit code 1
hsauro commented 1 month ago

It’s possible onestep and simulate are different. One thing I do know about onestep is that one step has to restart the simulation each time, this is a requirement of cvode whereas simulate doesn’t do this.

We can take a look.

Herbert Sauro, Professor Director: NIH Center for model reproducibility University of Washington, Bioengineering 206-685-2119, www.sys-bio.org, http://reproduciblebiomodels.org/ Mobile: 206-880-8093 @.*** Books: http://books.analogmachine.org/

On Thu, Jul 25, 2024 at 3:17 AM Matthias König @.***> wrote:

Hi all, I see the following issue a lot. If I run the identical simulation to a simulate via a onestep simulation the integrator fails with CV_CONV_FAILURE.

It is the identical simulation and model just once run via simulate, the other time via onestep. The only thing I can imagine is that roadrunner/CVODE is using different simulator/integrator settings for the two cases (e.g. retries, minimal step size based on end time, ...). I see this in many models and the iterative one step fails very often with convergence failures.

We use this a lot in model coupling within PDEs where we just continue the simulation for one step.

Example attached for the following model: droplets_iri_liver.zip https://urldefense.com/v3/__https://github.com/user-attachments/files/16374266/droplets_iri_liver.zip__;!!K-Hz7m0Vt54!lgl9DsGLm8y33lDhQBLlTnGLzByINfmDrHT1aGrgrF6BmBcMXYaa3oHRsd6Hs-OZuyVELUY-p8H4y1PsavEzsjxc0igebQ$

"""Example script for porous media coupling for IRI & Droplet model."""from pathlib import Pathimport pandas as pdimport roadrunner

-------------------------------------------------------------------------------------# load model# -------------------------------------------------------------------------------------model_path = Path(file).parent / "droplets_iri_liver.xml"

loading the model in roadrunnerr = roadrunner.RoadRunner(str(model_path))

-------------------------------------------------------------------------------------# set changes for given simulator# -------------------------------------------------------------------------------------vol_point = 1.0 # [l]f_fluid = 0.3 # [dimensionless] fluid fractionf_solid = 1 - f_fluid # [dimensionless] solid fraction

These are the changes which have to be appliedchanges = {

# setting volumes
"Vext": vol_point * f_fluid,  # [liter] fluid phase volume
"Vli": vol_point * f_solid,  # [liter] solid phase volume

# concentrations of external concentrations in fluid phase
"[glc_ext]": 5.0,  # [mM]
"[lac_ext]": 1.0,  # [mM]
"[o2_ext]": 1.2,  # [mM]
"[alt_ext]": 0.0,  # [mM]
"[ast_ext]": 0.0,  # [mM]
"[ffa_ext]": 0.0,  # [mM]
"[vldl_ext]": 0.0,  # [mM]

"[glc]": 5.0,  # [mM]
"[lac]": 1.0,  # [mM]
"[o2]": 1.2,  # [mM]
"[alt]": 1.0,  # [mM]
"[ast]": 0.7,  # [mM]
"[pyr]": 0.8,  # [mM]
"[atp]": 3.0,  # [mM]
"[adp]": 1.0,  # [mM]
"[nadh]": 3.9,  # [mM]
"[nad]": 0.1,  # [mM]
"[ros]": 0.0,  # [mM]

# set temperature
"temperature": 37 + 273.15,  # [Kelvin] warm simulation
# "temperature": 4 + 273.15,  # [Kelvin] cold simulation

}

-------------------------------------------------------------------------------------# Timecourse# -------------------------------------------------------------------------------------# run simulationtend = 600 # [min]steps = 100dt = tend/steps # [min]

single simulation (singel simulation works)print("--- simulate --- ")r.resetAll()for key, value in changes.items():

r.setValue(key, value)s1 = r.simulate(start=0.0, end=tend, steps=steps) #*60df1 = pd.DataFrame(s1, columns=s1.colnames)print(df1.head())

FEBio version (onestep)print("--- onestep --- ")r.resetAll()for key, value in changes.items():

r.setValue(key, value)

t_model = 0.0while t_model < tend: s = r.oneStep(t_model, dt) t_model += dtdf2 = pd.DataFrame(s, columns=s.colnames)print(df2.head)

Results in

/home/mkoenig/.virtualenvs/pkdb_models_new/bin/python /home/mkoenig/git/pkdb_models/pkdb_models/models/tpm_fat/simulations/iri_droplets_liver_rr.py --- simulate --- time [ins] [ffa] [vldl_ext] [dg] [protein_dgat2] [protein_ldsyn] [tg1] [tg2] [tg3] [tg4] [glc_ext] ... [ast_ext] [pyr] [atp] [adp] [nadh] [nad] [ros] [alt] [ast] [necrosis] [tgd] tsim 0 0.0 0.000001 2.142857 0.000000 1.714286 1.428571 1.428571 0.000000 0.000000 0.000000 0.000000 5.000000 ... 0.0 0.800000 3.000000 1.000000 3.900000 0.100000 0.0 1.0 0.7 0.0 0.000000 0.0 1 6.0 0.000001 0.000449 0.020722 1.726787 1.428571 1.428571 177.988176 10.942289 0.261533 0.005179 4.908386 ... 0.0 0.729103 3.038536 0.961464 3.989563 0.010437 0.0 1.0 0.7 0.0 189.197177 6.0 2 12.0 0.000001 0.000449 0.022024 1.727565 1.428571 1.428571 270.007558 37.756799 1.758146 0.044349 4.820259 ... 0.0 0.752057 3.060534 0.939466 3.989774 0.010226 0.0 1.0 0.7 0.0 309.566852 12.0 3 18.0 0.000001 0.000449 0.023327 1.728343 1.428571 1.428571 309.703061 70.507305 6.067988 0.181956 4.734098 ... 0.0 0.774237 3.070634 0.929366 3.989982 0.010018 0.0 1.0 0.7 0.0 386.460310 18.0 4 24.0 0.000001 0.000449 0.024630 1.729120 1.428571 1.428571 317.444936 102.337022 14.189627 0.545982 4.648687 ... 0.0 0.796205 3.071470 0.928530 3.990204 0.009796 0.0 1.0 0.7 0.0 434.517567 24.0

[5 rows x 30 columns] --- onestep --- Error: CVODE Error: CV_CONV_FAILURE, Module: CVODES, Function: CVode, Message: At t = 180.806 and h = 4.48854e-07, the corrector convergence test failed repeatedly or with |h| = hmin. Traceback (most recent call last): File "/home/mkoenig/git/pkdb_models/pkdb_models/models/tpm_fat/simulations/iri_droplets_liver_rr.py", line 78, in s = r.oneStep(t_model, dt) ^^^^^^^^^^^^^^^^^^^^^^ File "/home/mkoenig/.virtualenvs/pkdb_models_new/lib/python3.11/site-packages/roadrunner/roadrunner.py", line 2968, in oneStep return _roadrunner.RoadRunner_oneStep(self, currentTime, stepSize, reset) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ RuntimeError: CVODE Error: CV_CONV_FAILURE: Convergence test failures occurred too many times (= MXNCF = 10) during one internal timestep or occurred with |h| = hmin.; In virtual double rr::CVODEIntegrator::integrate(double, double)

Process finished with exit code 1

— Reply to this email directly, view it on GitHub https://urldefense.com/v3/__https://github.com/sys-bio/roadrunner/issues/1232__;!!K-Hz7m0Vt54!lgl9DsGLm8y33lDhQBLlTnGLzByINfmDrHT1aGrgrF6BmBcMXYaa3oHRsd6Hs-OZuyVELUY-p8H4y1PsavEzsjyMtwYgjg$, or unsubscribe https://urldefense.com/v3/__https://github.com/notifications/unsubscribe-auth/AAIBSDVHHGSFGQXQJYZLYJTZODGBZAVCNFSM6AAAAABLOHERXGVHI2DSMVQWIX3LMV43ASLTON2WKOZSGQZDSNJYGY4TOMQ__;!!K-Hz7m0Vt54!lgl9DsGLm8y33lDhQBLlTnGLzByINfmDrHT1aGrgrF6BmBcMXYaa3oHRsd6Hs-OZuyVELUY-p8H4y1PsavEzsjyTFdTBUQ$ . You are receiving this because you are subscribed to this thread.Message ID: @.***>

matthiaskoenig commented 3 weeks ago

Hi all, any tips/info on this issue. The same thing also happens with simulate for the model. I.e., the model integrates successfully if run as:

r.simulate(start=0.0, end=tend, steps=steps)

But fails with

t_model = 0.0
while t_model < tend:
    s = r.simulate(start=t_model, end=t_model+dt, steps=2)
    t_model += dt

For me this does not make any sense. Either the integrator can handle the model with the given settings and tolerances or not. It should not matter if the integration is done in one run or stepwise ?! I can only imagine that certain CVODE settings are calculated based on tend (e.g. minimal_step_size or maximum_step_size or similar things) which could affect the integration.

How can I make the integration work?

We are using this for multi-framework coupling. I.e. we have to run small integration steps which are then coupled to FEM/PDE simulations as source and sink terms.