[Scilab-users] Sundials ode solvers

Wed Oct 27 16:18:25 CEST 2021

Hi Stéphane,

I tested the cvode_solve() function in my case study cited below and it actually ran 40% faster than the previous ode(‘stiff’,…), with easy one-to-one function’s parameters replacement (in my case including the list syntax, with function’s multiple parameters, gradient function, simple tolerance specifications).
As you advised, I need to investigate moving the differential equation code from Scilab to fortran external function to reduce further the simulation time.

Thank you for this new module, staying updated with Sundials !

Best,

David Chèze

tel. +33(0)479792130 – CEA ☼ INES – LITEN/DTS/SIRE/LELA –  Le Bourget du Lac (73), FRANCE – bât. HELIOS/3108 Liten.cea.fr

De : users <users-bounces at lists.scilab.org> De la part de Stéphane Mottelet
Envoyé : lundi 25 octobre 2021 16:35
À : Users mailing list for Scilab <users at lists.scilab.org>
Objet : [Scilab-users] Sundials ode solvers

Hi all,

As discussed in this thread last month

https://www.mail-archive.com/users@lists.scilab.org/msg10679.html

I am glad to annouce that a first version of the  sci-sundials toolbox (maybe part of Scilab in the future) is available on Atoms in the "Differential Equations" category (refresh the package list if you don't see it).  To have an idea of its features I pasted below the content of the README.md file on the gitlab project of sci-sundials (https://gitlab.com/mottelet/sci-sundials/).

If you appreciate the work and want to help, doc, demos, or even more if you know how to code in C and C++ you are welcome !

S.

What has been done and what is to be done

Until now only CVODE has been interfaced but many common features have been developped in the OdeManager class, hence interfacing IDA will be quite easy. ARKODE (various modern explicit, implicit and mixed explicit/implicit Runge-Kutta solvers) does not exist in the old 2.4.0 version of Sundials which is used in Scilab hence won't be interfaced unless an upgrade is done (Sundials is now at version 5.7.0). The support of Sparse Jacobians is also missing for the same reason.

Features

The CVODE gateway cvode_solve() implements the following features which were missing by the legacy LSODE/LSODA/LSODAR/... ode() gateway :

  *   full a posteriori access to solver continuous extension of solution at arbitrary points via a MList output when only one lhs is given:

sol = cvode_solve(f, [t0 tf], y0)

t = linspace(t0, tf, 1000)

plot(t, sol(t))

The sol MList fields gathers all information related to the obtained solution (time steps, solution at time steps, events, ...). The solver can be restarted by giving the MLlist as first argument to cvode_extend :

sol2 = cvode_extend(sol, tx, ...)

where tx is the time point to which solution has to be extended. The options of the call that yielded sol are used and can be changed as optional named parameters after tx.

  *   a user-friendly access to solver options via optional named parameters in cvode_solve call, i.e.

cvode_solve(f, tspan, y0, h0=0.01, rtol=1e-3)

  *   a really simpler way to give time span of integration allowing to choose between error driven solver internal time steps and user fixed time steps, i.e.

[t,y] = cvode_solve(f, [t0 tf], y0)

[t,y] = cvode_solve(f, [t0 t1 ... tf], y0)

[t,t] = cvode_solve(f, tspan, y0, t0=0)

the latter style being the closest to actual ode() behavior where solution is by default given at user time steps and t0 is not necessarily equal to tspan(1), which is the default in the two former calls above.

  *   a better and user-friendly specification of events via a variable number of outputs event function (giving value of event equations, wheter to stop integration for a given event and event direction selection), minimal style beeing a single output. Information about events is also simpler to get:

function [eq,term,dir] = evfun(t,y)

    eq = y(1)-1.7;

    term = %f;

    dir  = 1;

end

[t,y,te,ye,ie] = cvode_solve(f, tspan, y0, events = evfun)

sol = cvode_solve(f, tspan, y0, events = evfun)

in the latter call information about events is recovered in sol.te,sol.ye,sol.ie.

  *   support of complex solution with detection of complexity in y0 or f(t0,y0), e.g.

function out = crhs(t,y)

    out = 10*exp(2*%i*%pi*t)*y;

end

[t,y] = cvode_solve(crhs, [0,5], 1)

plot(t,real(y),t,imag(y))

  *   Support of a callback function called after each successfull step, giving access to current solver statistics and allowing to stop integration (e.g. by a "stop" button on a GUI), e.g.

function stop = scicallback(t,y,flag,stats)

    stop = %f

    if flag == "step"

        mprintf("%s : hlast=%g\n", flag, stats.hlast)

    end

end

[t,y] = cvode_solve(f, tspan, y0, intcb=scicallback);

  *   support of an arbitrary number of dimensions of ode state

Performance

cvode_solve() is roughly two times faster than ode() for both fixed methods (Adams and BDF). As ode() already did, compiled and dynamically linked C,C++ or Fortran externals are supported by cvode_solve(). When using such externals instead of Scilab functions cvode_solve() is generally an order of magnitude faster.

--

Stéphane Mottelet

Ingénieur de recherche

EA 4297 Transformations Intégrées de la Matière Renouvelable

Département Génie des Procédés Industriels

Sorbonne Universités - Université de Technologie de Compiègne

CS 60319, 60203 Compiègne cedex

Tel : +33(0)344234688

http://www.utc.fr/~mottelet
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20211027/af1d21ca/attachment.htm>