Has anyone used DVODE and is it good?

Thank you all! @ivanpribec @certik @jacobwilliams @nicholaswogan

Yes @jacobwilliams,
I actually have both @prince_mahajan 's FLINT and your DOP853 and DDEABM at hand. The structure of both of your code are similar, both are modernized. I learned a lot from your code!
Perhaps according to FLINT’s benchmark, FLINT’s DOP853 is slightly faster. But I have not done through test myself.
My plan is to put FLINT, your DOP853 and DDEABM inside a package I am developing (which needs solving ODE to make prediction), and also add at least an ODE solver for stiff problem. Currently I may just add the f90 version DVODE.

@jacobwilliams ,
May I ask, what do you mean by modernized, Uhm, does that mean making Dvode object-oriented?
Will that make DVODE faster?
Thank you very much!

@nicholaswogan
Thank you very much for your suggestion of CVODE, which is a part of Sundials. It should be great.
Uhm, I currently may just try to use a Fortran ODE instead. It is relatively easy for me to use Sundials for my own research. But I am developing a package, adding Sundials inside my package may be a little bit not not too easy right now.

By the way, has anyone by any chance tried both DVODE and CVODE, if so, is CVODE faster than DVODE for the same ODEs?

What we have done and agreed upon for other packages such as FFTPACK, modernize means: compile with fpm, use free form, possibly add into modules, ensure all interfaces are present, remove implicit typing, add tests, documentation, examples and so on.

1 Like

Time difference between CVODE and DVODE will depend on the problem but I’m guessing it will be small, like within a factor of 2. The exception is if you want to do some really big problem. Then CVODE will win/allow a solution to be calculated because it has stuff like GMRES.

2 Likes

@certik @jacobwilliams
Alright. Uhm, do you think the current DVODE by Byrne and Thompson is modernized? link is below,
File Listing (radford.edu)

I have downloaded the DVODE.f90 distribution file there,
image

uhm, it seem it is somehow modernized a little bit.
It uses module and f90 format, and implicit none.
In windows it looks like below, and it can compile and run without problem,

Seem not too bad.

@CRquantum im not able to access the code via the link your originally posted. Do you know an alternative link?

Thank you very much for the message @nicholaswogan ,
true, strangely, the original link does not work any more. Fortunately, I have downloaded their whole website (the latest version) of dvode before.
Now I pushed it to a new gitlab repo, you could download it from the link below,

In the folder VODE_F90 Support Page , double click index.html,
then you should see the whole original website, and you can download all the corresponding files locally.
Hope that helps :slight_smile:

2 Likes

I have not test the speed of dvode, but I guess it should not be too bad.
dvode can handle both stiff and non-stiff problem.
The only issue seems to be that dvode does not natively support interpolation.

Thanks!! This looks great. I think the code does have interpolation. Section 5 of the documentation, which is all in dvode_f90_m.f90, talked about interpolation.

1 Like

Great! Thank you so much for letting me know! I did not read the file carefully enough LOL.
Right, from line 1554 there are instructions about how to do interpolation using subroutine DVINDY.
Then the dvode seems a quite versatile ODEs solver!

Hi @nicholaswogan ,
Have you by any chance tried the interpolation function in DVODE?
It is the subroutine DVINDY.

It is a little strange that, I solve the ODE ( just say Y(t), I want Y at given t ) within the time t range, say from T0 to Tf.
Then I can call DVINDY to interpolate Y for t>Tf. (But this is strange, it is more like extrapolation than interpolation)

However , if I just want to interpolate Y at t within T0 to Tf, then DVINDY give me some error, something like t is not within the range like below,

Error in DVINDY, T(=R1) is illegal. T is not
In the above message, R1 =   0.6000000000000D+01
 in interval TCUR - HU(= R1) to TCUR(=R2)
In the above message, R1 =   0.2364435159699D+02   R2 =   0.7499440289459D+02

Just curious, have you meet similar issues?

I think DVINDY can correctly do interpolation, but I just do not know which part did I do is wrong… LOL.

The interpolation routine DVINDY has the following requirement, which is stated in the source code:

For valid results, T must lie between TCUR - HU and TCUR.

In other words, if interpolation is desired, it should be performed after each successful integration step, which covers Tcur - Hu to Tcur (both evaluated with the current value of Tcur).

If, instead, you perform a number of successive integration steps and only then try to interpolate over the entire range of the independent variable, that will not be possible. The capability for doing that could have been provided, but that would have required a large amount of supporting information to be stored. At the time the software was published (1989), how much RAM did our PCs have?

You can verify for yourself the relation between T and Tcur by adding the following print statement at the beginning of subroutine DVINDY_CORE:

print '(A,3ES12.4)','T, TP, TN = ',T,TP,TN

and running the two example problems in your Gitlab repository.

1 Like

Thank you very much @mecej4 !

What I want is, call DVODE between time [0, 24], then interpolate points at [0,2,4,6,8,10,…22,24].

So, I mean if I call

CALL DVODE_F90(..., T=0,TOUT=24,...)

Where T=0 is the initial time, TOUT=24 is the final time.
So the integration is in the range of [0,24].
I would expect then I should be able to do things like

CALL DVINDY(Time,...)

where Time should be within [0,24], am I right?

But DVINDY does not seem to interpolate as I said above.
It seems after I do

CALL DVODE_F90(..., T,TOUT,...)

T becomes TOUT, then DVINDY can only interpolate at time between like [24-xxx, 24+yyy], such as [23, 75]. But Definitely not the original range [0,24] which is what I really want.

Overall, it seems DVINDY is not particularly useful. It can only interpolate at its own particular range TCUR - HU and TCUR. Instead of the original range [0,24].

I may just CALL DVODE_F90 between [0,2], then [2,4], then [4,6], … [22,24].

Are you saying adding the print you mentioned at the position below in dvode_f90_m.f90?


Line 6941 is where I added the print you mentioned.
But here TP seems has not been calculated yet?

Or, do I add the print at line 6947 below?

Do you know if there is quick way that I can add something in DVODE so that all those points are stored? Then I guess DVINDY will be more useful.

Edit: the follow-up post from @mecej4 addresses the actual problem of output at regular intervals, without the need to call DVINDY.


To use interpolation you should consult the instructions under itask in dvode (original or _f90):

! ITASK  = An index specifying the task to be performed.
!          Input only.  ITASK has the following values and meanings.
!          1  means normal computation of output values of y(t) at
!             t = TOUT (by overshooting and interpolating).
!          2  means take one step only and return.
!          3  means stop at the first internal mesh point at or
!             beyond t = TOUT and return.
!          4  means normal computation of output values of y(t) at
!             t = TOUT but without overshooting t = TCRIT.
!             TCRIT must be input as RWORK(1).  TCRIT may be equal to
!             or beyond TOUT, but not behind it in the direction of
!             integration.  This option is useful if the problem
!             has a singularity at or beyond t = TCRIT.
!          5  means take one step, without passing TCRIT, and return.
!             TCRIT must be input as RWORK(1).

I assume you are currently using itask = 1. Probably you’d want to use 2 (or 5 in case you have a known singularity) and integrate in “one-step” mode. This is a fairly common pattern in classic Fortran ODE solvers.

A similar mechanism is used also in RKC (Runge-Kutta-Chebyshev) with which I’m more familiar. The rkc subroutine has the interface:

      subroutine rkc(neqn,f,y,t,tend,rtol,atol,info,work,idid)

The array of integers info is used to control the procedure. Specifically info(1) fulfill the same role as itask:

c  INFO(1):  RKC integrates the initial value problem from T to TEND. 
c            This is done by computing approximate solutions at points 
c            chosen automatically throughout [T, TEND].  Ordinarily RKC
c            returns at each step with an approximate solution. These
c            approximations show how y behaves throughout the interval.
c            The subroutine RKCINT can be used to obtain answers anywhere
c            in the span of a step very inexpensively. This makes it 
c            possible to obtain answers at specific points in [T, TEND]
c            and to obtain many answers very cheaply when attempting to
c            locating where some function of the solution has a zero
c            (event location).  Sometimes you will be interested only in
c            a solution at TEND, so you can suppress the returns at each
c            step along the way if you wish.
c
c  INFO(1)  = 0 Return after each step on the way to TEND with a
c               solution Y(*) at the output value of T.
c
c          ` = 1 Compute a solution Y(*) at TEND only.
c

An example of integration with interpolation to produce output at regularly spaced time intervals can be found in rkc-exa.f.

1 Like

No, your expectations do not match the design of the VODE integrator package.

Look at the code in example1.f90 and see how it calls the integrator inside a loop. For your desired calculation you would have to set up a loop that is run 12 times, with Tout = 2, 4,…,24.

Note also that the driver program in example1.f90 does not itself call DVINDY. Such calls could be needed if, for example, you wanted to obtain many more values for plotting the solutions (which is sometimes labelled “dense output”).

1 Like

I think the purpose of DVINDY is actually to get derivatives of y as stated in this comment:

!  CALL DVINDY(,,,,,)         Provide derivatives of y, of various
!        (See below.)         orders, at a specified point T, if
!                             desired.  It may be called only after
!                             a successful return from DVODE.

and not interpolation per se.

While the solution with itask=2 might work, I agree with what @mecej4 just suggested.

1 Like

Note that the argument K of DVINDY, which is the order of the derivative, can be zero. In fact, subroutine DVODE makes calls with K = 0 to DVINDY_CORE. Such calls interpolate the function values. The user may make calls with other values of K if derivative values are wanted.

1 Like

Thank you very much @ivanpribec @mecej4 !
I really appreciate your insights!

I have tried ITASK = 2 and other values for DVODE, it seems no matter what number ITASK is, the current DVODE cannot do the interpolation I mentioned.
I guess like you suggested, DVINDY may be more for some dense output, and calculate some derivatives.

For now, I will just use DVODE and do integrals from 0 to 2, 2 to 4, 4 to 6, … 22 to 24, instead of do 2 to 24 then do interpolation at 2,4,6,8,…,22,24.

Here’s an example of using dvindy for dense output.

Example
program droplet

  implicit none

  integer, parameter :: dp = kind(1.0d0)
  real(dp), parameter :: pi = 4.0_dp*atan(1.0_dp)

  integer, parameter :: neq = 3

  integer ::  itol, itask, istate, iopt, mf, iflag, n, ngrid
  real(dp) :: rtol, s, send, ds
  real(dp), dimension(neq) :: atol, y

  real(dp), allocatable :: ygrid(:,:), sgrid(:)

  real(dp), allocatable :: rwork(:)
  integer, allocatable :: iwork(:)
  integer :: lrw, liw

  real(dp) :: rpar(2)
  integer :: idum(1)

  external :: dvode, dvindy

  rpar(1) = 3.3888_dp
  rpar(2) = 0.8414_dp

  s = 0.0_dp
  send = 2.0_dp

  ngrid = 41
  allocate(sgrid(ngrid),ygrid(neq,ngrid))
  sgrid(1) = 0.0_dp
  ygrid(:,1) = [0.0_dp, 0.0_dp, 0.0_dp]

  ds = (send - s)/real(ngrid - 1,dp)
  do n = 2, ngrid
    sgrid(n) = (n - 1)*ds
  end do
  n = 2

  ds = 0.00001_dp
  s = ds

  y(1) = ds
  y(2) = cos(y(1))*ds
  y(3) = sin(y(1))*ds


  itol = 2
  rtol = 1.0e-7_dp
  atol(1) = 1.0e-6_dp
  atol(2) = 1.0e-6_dp
  atol(3) = 1.0e-6_dp

  mf = 22
  lrw = len_rwork(neq,mf)
  liw = len_iwork(neq,mf)
  allocate(rwork(lrw), iwork(liw))

  istate = 1
  itask = 2
  iopt = 0

  integrate: do

    call dvode(young_laplace, neq, y, s, send, &
        itol, rtol, atol, itask, istate, iopt, &
        rwork, lrw, iwork, liw, &
        young_laplace_jacobian, &
        mf, rpar, idum)

    ! Automatic stepping
    write(*,'(4(ES12.4))') s, y

    if (istate < 0) then
      write(*,*) ' Failed at s = ',s,' with istate = ',istate
      error stop
    end if

    ! Interpolation
    interpolate: do while(s >= sgrid(n))
      call dvindy (sgrid(n), 0, rwork(21), neq, ygrid(:,n), iflag)
      n = n + 1
      if (n > ngrid) exit interpolate
    end do interpolate

    if (s >= send) then
      exit integrate
    end if

  end do integrate

  write(*,*)
  write(*,*)

  do n = 1, ngrid
    write(*,'(4(ES12.4))') sgrid(n), ygrid(:,n)
  end do

contains

  subroutine young_laplace(neq, t, y, ydot, rpar, ipar)
    integer, intent(in) :: neq
    real(dp), intent(in) :: t, y(neq)
    real(dp), intent(out) :: ydot(neq)
    real(dp), intent(in) :: rpar(*)
    integer, intent(in) :: ipar(*)

    associate(pL => rpar(1), drho => rpar(2))
      ydot(1) = pL - drho*y(3) - sin(y(1))/y(2)
      ydot(2) = cos(y(1))
      ydot(3) = sin(y(1))
    end associate

  end subroutine

  subroutine young_laplace_jacobian(neq, t, y, ml, mu, pd, nrpd, rpar, ipar)
    integer, intent(in) :: neq, nrpd, ml, mu
    real(dp), intent(out) :: pd(nrpd,neq)
    real(dp), intent(in) :: t, y(neq), rpar(*)
    integer, intent(in) :: ipar(*)

    ! dummy routine
    associate(pL => rpar(1), drho => rpar(2))
      pd(1,1) = -cos(y(1))/y(2)
      pd(2,1) = -sin(y(1))
      pd(3,1) =  cos(y(1))

      pd(1,2) = sin(y(1))/(y(2))**2
      pd(2,2) = 0.0_dp
      pd(3,2) = 0.0_dp

      pd(1,3) = -drho
      pd(2,3) = 0.0_dp
      pd(3,3) = 0.0_dp
    end associate

  end subroutine

  function len_rwork(neq, mf, mlu) result(lrw)
    integer, intent(in) :: neq, mf
    integer, intent(in), optional :: mlu(2)
      ! Upper and lower bandwidth in case of banded Jacobian

    integer :: lrw

! RWORK  = Real work array of length at least:
!             20 + 16*NEQ                      for MF = 10,
!             22 +  9*NEQ + 2*NEQ**2           for MF = 21 or 22,
!             22 + 11*NEQ + (3*ML + 2*MU)*NEQ  for MF = 24 or 25.

    select case(mf)
    case(10)
      lrw = 20 + 17*neq
    case(21,22)
      lrw = 22 + 9*neq + 2*neq**2
    case(24,25)
      if (present(mlu)) then
        lrw = 22 + 11*neq + (3*mlu(1) + 2*mlu(2))*neq
      else
        write(*,*) "len_rwork: Missing array mlu containing bandwidths [ml, mu]"
        error stop
      end if
    case default
      ! Wrong value for mf
      lrw = -1
    end select
    
  end function

  function len_iwork(neq, mf) result(liw)
    integer, intent(in) :: neq, mf

    integer :: liw

! IWORK  = Integer work array of length at least:
!             30        for MF = 10,
!             30 + NEQ  for MF = 21, 22, 24, or 25.
!          If MF = 24 or 25, input in IWORK(1),IWORK(2) the lower
!          and upper half-bandwidths ML,MU.

    select case (mf)
    case (10)
      liw = 30
    case (21, 22, 24, 25)
      liw = 30 + neq
    case default
      ! Wrong value for mf
      liw = -1
    end select

  end function

end program

The system solved are the Young-Laplace equations for a pendant drop in a cylindrical coordinate system:

\frac{d\phi}{ds} = \tilde{p}_L - \Delta\tilde{\rho}\,z - \frac{\sin{\phi}}{r} \\ \frac{dr}{ds} = \cos{\phi}\\ \frac{dz}{ds} = \sin{\phi}

An equivalent system of equations (but with a slightly different set of parameters) was first solved in 1883 by Francis Bashforth and John Couch Adams in An attempt to test the theories of capillary action by comparing the theoretical and measured forms of drops of fluid - GDZ. Notably, for this occasion Adams derived the family of explicit multi-step methods, now called Adams-Bashforth methods.

In VODE, when using mf = 10 the solver uses an implicit method from Adams. Here is plot of the interpolated (r,z) pairs super-imposed on the automatic step control:

The resulting profile traces half of a pendant drop. The interpolated points are spaced equidistantly in terms of the arc length s.

By fitting the theoretical drop shape to an experimental shape profile, the surface tension of a liquid can be determined. The method is called pendant drop tensiometry.

1 Like

That is awesome @ivanpribec . Thank you my man!
Cool, I see. So you set ITASK = 2, so all the automated s and y are saved. Then you interpolate y at corresponding s. Brilliant.

I think in this way as you did, we actually can achieve what a normal interpolation ODE would do for DVODE. Save all the time (in your case s) and the y array along the way between the initial and final time. Then do interpolation at corresponding s between the initial and final time!

Thanks a again! Very impressive!

PS. I used your droplet.f90 combined with the dvodeoriginal.f, I use intel fortran.
In the debug mode, it shows an error

droplet.f90(67): error #8284: If the actual argument is scalar, the dummy argument shall be scalar unless the actual argument is of type character or is an element of an array that is not assumed shape, pointer, or polymorphic.   [RTOL]

But in release mode it compile just fine and run perfectly.

But I checked, it does looks like your rtol is a scalar, do you know what is wrong? Perhaps it is a Intel compiler bug or somthing?

I was relying on an external interface. The issue originates here:

      SUBROUTINE DVODE (F, NEQ, Y, T, TOUT, ITOL, RTOL, ATOL, ITASK,
     1            ISTATE, IOPT, RWORK, LRW, IWORK, LIW, JAC, MF,
     2            RPAR, IPAR)
      EXTERNAL F, JAC
      DOUBLE PRECISION Y, T, TOUT, RTOL, ATOL, RWORK, RPAR
      INTEGER NEQ, ITOL, ITASK, ISTATE, IOPT, LRW, IWORK, LIW,
     1        MF, IPAR
      DIMENSION Y(*), RTOL(*), ATOL(*), RWORK(LRW), IWORK(LIW),
     1          RPAR(*), IPAR(*)

since RTOL(*) is an asummed-size array and not a scalar like in my calling code.

Easiest way to get rid of the error is to change rtol from a scalar to a one-element array in the calling code:

real(dp) :: rtol(1)  ! Array rtol
1 Like