modelon-community / Assimulo

Assimulo is a simulation package for solving ordinary differential equations.
https://jmodelica.org/assimulo/index.html
GNU Lesser General Public License v3.0
70 stars 17 forks source link

ncp_list timepoint evaluations do not match output #78

Closed PHvanLent closed 11 months ago

PHvanLent commented 11 months ago

I have noticed that for some simulations for this code, the list of timepoints passed at ncp_list does not lead to the same number of evaluations in the output sol. Is this something that is to be expected?

mod = Explicit_Problem(self.func, self.y0, np.min(t))
sim=CVode(mod)

sim.rtol=self.rtol
sim.atol=self.atol
sim.verbosity=0
sim.maxsteps=10000
sim.time_limit=600
sim.linear_solver = 'SPGMR'
sol=sim.simulate(tfinal=np.max(t),ncp=0,ncp_list=list(t))

E.g. 20 timepoints as input, but only 17 have output. From what I saw, these values are typically some in-between values.

Kind regards,

Paul

PeterMeisrimelModelon commented 11 months ago

Hi Paul,

can you provide a complete example or minimal reproducer for this?

PHvanLent commented 11 months ago

Hello Peter, thanks for the fast reply. Here is a small example:

from assimulo.solvers import CVode
from assimulo.problem import Explicit_Problem
import numpy as np
import matplotlib.pyplot as plt

t=[ 0.          ,2.63157895,  5.26315789,  7.89473684 ,10.52631579, 13.15789474,
 15.78947368, 18.42105263, 21.05263158, 23.68421053, 26.31578947, 28.94736842,
 31.57894737, 34.21052632, 36.84210526, 39.47368421, 42.10526316, 44.73684211,
 47.36842105, 50.        ]

def Batch(t,y):

    qsmax=-0.4288179560674244 #CmolS.Cmolx-1.h_1
    Ks= 0.0231995436988264 #Cmol.S.L-1
    a=  -1.7781239845815182 #S requirement for biomass
    ms=-0.0003873332954076 #CmolS.Cmolx-1.h-1

    Cs,Cx=y
    if Cs<0:
        Cs=0
    #Kinetic equations
    qs=qsmax*(Cs/(Ks+Cs))
    qx=(qs-ms)/a
    #Rate equations
    Rs=qs*Cx*1
    Rx=qx*Cx*1
    #Rc=-(Rs+Rx)  #Rco2 and Ro2
    R=[Rs,Rx]
    return R

# {'qsmax': -0.4288179560674244, 'Ks': 0.0231995436988264, 'a': -1.7781239845815182, 'ms': -0.0003873332954076}
y0=[3, 0.01]
mod = Explicit_Problem(Batch, y0, t[0])
sim = CVode(mod)
sim.linear_solver='SPGMR'
sim.verbosity=40
print("# input evaluations ",len(t))
t,y=sim.simulate(tfinal=np.max(t),ncp=0,ncp_list=list(t))
print("# output evaluations ",len(t))

This outputs: input evaluations 20 output evaluations 16

I was previously using an older version assimulo (not sure which one), but recently updated to version 3.4.3. In this case it only happens when I use the linear solver SPGMR, but also when using the dense solver I have seen this issue.

Thank you very much.

Kind regards,

Paul

PeterMeisrimelModelon commented 11 months ago

Hi Paul,

I cannot reproduce the issue, input & output evaluations match up for me when running this example. Can you provide some more information on your setup, e.g., Python & numpy versions.

/Peter

PHvanLent commented 11 months ago

Hi Peter,

I am using python 3.8.18, numpy 1.24.4, assimulo 3.4.3, sundials 6.6.2.

I have also added a list with all other packages in the environment in the attachment.

Kind regards,

Paul

all_packages_env.txt

PeterMeisrimelModelon commented 11 months ago

Is this on Windows/Linux/Mac?

Can you provide the incorrect t list you get as output?

/Peter

PHvanLent commented 11 months ago

Hello Peter,

I am working on linux. It seems to miss values before the last one:

[0.0, 2.63157895, 5.26315789, 7.89473684, 10.52631579, 13.15789474, 15.78947368, 18.42105263, 21.05263158, 23.68421053, 26.31578947, 28.94736842, 31.57894737, 34.21052632, 36.84210526, 50.0]

PeterMeisrimelModelon commented 11 months ago

I can reproduce the issue, it appears to be an issue with higher Sundials versions. There are currently no plan to enhance support for higher Sundials versions and we recommend to use Sundials 2.7.0 instead.

Please let us know if you experience the same or similar issues with Sundials 2.7.0.