Has anyone used DVODE and is it good?

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