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
Thank you @ivanpribec ! That is great! My bad, did not look carefully enough that RTOL is assumed-size, LOL.
It seems for the Jacobian matrix, it is by default 0 in DVODE, so there is no need to set 0 elements to PD, pd(2,2), pd(2,3), pd(3,2), pd(3,3) can be commented.
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
By the way, for stiff problem, providing a Jacobian matrix .vs. without Jacobian matrix,
in your experience, will that cause a huge performance difference?
I found that with Jacobian matrix, there is basically no particular performance gain, but it does improve the DVODEās stability and ability of reaching convergence.
Thanks!
PS.
Besides DVINDY, every time I can learn something new from your code 
I learned function/subroutine can contain function/subroutine inside them before. I apply it in combining DVODE and FLINT ODE solver, such that the function FEX and the Jacobian JEX can be defined and used in FLINT ODE class.
With interpolation the DVODE in FLINT can be,
subroutine dvode_interpolate(me, T0, TOUT, Y0, Y, K, T_array, Y_array, Params)
class(PM_ODE_Sys), intent(inout) :: me !< Differential Equation object
real(kind = r8), intent(in), dimension(:), optional :: Params
real(kind =r8) :: T0, T, TOUT, RTOL
real(kind=r8) :: Y0(:), Y(:)
integer :: NEQ, ITASK, ISTATE, IFLAG
real(kind=r8) :: T_array(:), Y_array(:,:)
integer :: K
integer(kind=i8) :: i, ngrid
NEQ = me%n
Y = Y0
T = T0
ITASK = 2
ISTATE = 1
ngrid = size(T_array,kind=i8)
i = 1
integrate : do while (T <= TOUT)
CALL DVODE_F90(FEX,NEQ,Y,T,TOUT,ITASK,ISTATE,OPTIONS,J_FCN=JEX)
IF (ISTATE<0) THEN
write(6,*) 'DVODE Failed at ',T,' with istate = ',istate
error stop
ELSE ! do interpolation
interpolate: do while(T >= T_array(i))
call DVINDY (T_array(i), K, Y_array(:,i),IFLAG)
i = i + 1
if (i > ngrid) exit interpolate
end do interpolate
END IF
enddo integrate
return
contains
SUBROUTINE FEX(NEQ,T,Y,YDOT)
INTEGER, INTENT (IN) :: NEQ
DOUBLE PRECISION, INTENT (IN) :: T
DOUBLE PRECISION, INTENT (IN) :: Y(NEQ)
DOUBLE PRECISION, INTENT (OUT) :: YDOT(NEQ)
YDOT = me%F(T,Y,Params)
RETURN
END SUBROUTINE FEX
SUBROUTINE JEX(NEQ,T,Y,ML,MU,PD,NRPD)
INTEGER, INTENT (IN) :: NEQ, ML, MU, NRPD
DOUBLE PRECISION, INTENT (IN) :: T
DOUBLE PRECISION, INTENT (IN) :: Y(NEQ)
DOUBLE PRECISION, INTENT (OUT) :: PD(NRPD,NEQ)
PD = me%J(NEQ,T,Y,ML,MU,NRPD)
RETURN
END SUBROUTINE JEX
end subroutine dvode_interpolate
where me%F is the differential equations and me% is the Jacobian matrix defined in FLINT ODE class me.
This time I learned associate,
associate(pL => rpar(1), drho => rpar(2))
Cool!
Thank you @ivanpribec !
I just learned from Dr. Fortran @sblionel the [ā¦] trick,
The below code also work, just add [ ] around atol and rtol to make array interface when call dvode (I see this trick in FLINT ode solver by @prince_mahajan too),
call dvode(young_laplace, neq, y, s, send, &
itol, [rtol], [atol], itask, istate, iopt, &
rwork, lrw, iwork, liw, &
young_laplace_jacobian, &
mf, rpar, idum)
By the way, have you by any chance used DVODPK? You know like LSODPK, etc, those end with PK or K, means using some thing called Preconditioned Krylov methods?
LLNLās ODEPACK webpage mentioned DVODPK as a variant for DVODE,
https://computing.llnl.gov/projects/odepack
It looks like CVODE in SUNDIALS is basically VODE + VODPK written in C and added some parallelization stuff.
Do you know if those Preconditioned Krylov methods efficient?
Thanks a lot!
1 Like
PLEASE NOTE:
VODE has been widely used for decades, so there are several diverging versions and descendants that a new user may try, without being aware of the other versions. Here are some comments, mostly from my recollections, and to be verified before being acted upon. In particular, if someone plans to use this solver for large problems (high order coupled differential equations in many variables), they should spend some time reading reports and papers before adopting a particular versionā¦
The version from Netlib uses the venerable Linpack solvers DGEFA/DGESV (for dense matrix representation) and DGBFA/DGBSV (for band matrix representation). Consider replacing these by the following Lapack routines: DGETRF/DGETRS and DGBTRF/DGBTRS, especially if you are using a library such as Intelās MKL.
The F90 version from Byrne and Thompson, of which @CRQuantum has a Gitlab backup, allows the use of a modified version of MA28 from the HSL library. MA28 is itself quite old, and there is a successor, MA48, available (licensed, but free for academic use).
Look at Sundials from LLNL, which is a modern successor to VODE. It is written in C, but has hooks for calling from Fortran, and allows for parallelization.
A large number of other modern sparse solvers are available, including Pardiso, SuperLU and SuiteSparse, which are direct solvers. Other iterative solvers based on Krylov-Arnoldi are available, as well. Thus, a prospective user of VODE should devote some time to planning the work before committing to one combination, trying it out on a small set of small problems, and then scaling it up for testing on large systems.
2 Likes
Thank you very much @mecej4 ! You comments are very helpful!
RIght. Actually the F90 version from Byrne and Thompson natively support LAPACK with DGETRF and DGBTRF, as well as MA48. This may be the best DVODE version.
There are some changing point in Byrne and Thompson code, just make some changes according to the instructions in the code, e.g.,
I have tested, the lapack or intel msl lapack works fine 
The MA48, seems need to be download from HSL,
https://www.hsl.rl.ac.uk/catalogue/hsl_ma48.html
DOVDE seems also have some preconditioned Krylov method variants, such as krysi and vodpk,
https://computing.llnl.gov/projects/odepack
But I have not figured out how to use those preconditioned Krylov variants, LOL.
Are those preconditioned Krylov method good?
I have used the rkc solver you mentioned. It seems very fast if setting atol = rtol as its illustration code did! Thanks!
Thank you @mecej4 !
About MA48, I just got the f90 version of HSL MA48 files which is the one the F90 version DVODE from Byrne and Thompson can use.
I follow the instructions in the F90 DVODE code, but finally I have a problem of getting the Jacobian routines for MA48, as indicated at the end of the F90 DVODE code below,
! *****MA48 build change point. Insert the MA48 Jacobian related
! routines DVNLSS48, DVSOLS48, DVPREPS48, and DVJACS48 here.
! filename = jacobian_for_ma48.f90. Insert the following line
! after the first line of this file:
! USE hsl_ma48_double
I do see DVODE have MA28 version of the Jacobian routines such as DVNLSS28, DVJACS28, DVSOLS28.
But there is no subroutines DVNLSS48, DVSOLS48, DVPREPS48, and DVJACS48.
Do you or anyone by any chance know what are those Jacobian subroutines DVNLSS48, DVSOLS48, DVPREPS48, and DVJACS48?
Thank you very much in advance!
PS.
In the F90 DVODE, there may be a typo, around line 1882, perhaps it should be ZD11_TYPE, I did not seem to find ZD01_TYPE in the MA48 files.
Look near line 12950 of DVODE_m.f90, where it says āIf you have licensed access to the HSL library, an alternate version of DVODE_F90 based on the successor to MA28, MA48, is available from the authors.ā
1 Like
Thank you @mecej4 !
I actually have emailed S. Thompson before, but he has not replied yet. Perhaps he was retired?
Byrneās email does not work.
May I ask, do you by any chance have a MA48 version of DVODE? You have it, right? LOL.
The dvode f90 seems actually good, and with ma48 should be great for sparse Jacobian system
Thanks @ivanpribec !
I briefly tested some ODEs by using DVODE and LSODA in the ODEPACK. It seems LSODA is about 2 times faster than DVODE (DVODE is set in stiff ode mode since it cannot automatically switch between stiff and non-stiff) and LSODAR can automatically switch between stiff and non-stiff problem which is also very good.
At least for the ODEs I tested, the performance of LSODA is very similar with RKC solver you mentioned before.
For DVODE which is basically VODE, it seems in terms of speed, it is not among the most efficient, I found some benchmark
Choosing a Solver ā Odes 2.6.2 documentation (scikits-odes.readthedocs.io)
Although it is in python, I believe the underline code is Fortran or C/C++. So the benchmark is roughly kind of accurate.
Some of these Fortran subroutines are wrapped in the R package deSolve. CRAN - Package deSolve
1 Like
Of course it doesā¦
A common programmerās mistake when using any kind of numerical method that asks for a Jacobian is to omit it, and let the algorithm calculate the Jacobian āautomaticallyā (=numerically). Never do that. Numerical differentiation may seem an easy task but it is not; it is actually one of the hardest numerical methods and it is prone to erroneous results. In well-behaving āniceā functions, numerical differentiation usually works well but do not rely on that. Even when it works well, it does so with a specific numerical integration method, which is not necessarily the one implemented in DVODE or whatever else. For example, numerical differentiation via splines may give accurate results with a ānot-a-knotā constraint and not so accurate with the popular ānaturalā spline.
Providing a numerical Jacobian to an algorithm that uses it extensively means looking for trouble. It means you give potentially inaccurate information, and that not once or twice, but literally in every single step of the algorithm. You are at the mercy of the specific numerical integration used and the functions involved in the problem at hand: Inaccuracies in the Jacobian (even slight ones) may or may propagate in the next steps of the DVODE (or whenever else algorithm), then you will wonder why it converges slowly or doesnāt converge at all.
Whenever you can provide the Jacobian analytically, do so.
1 Like
Thank you @Pap ! Nice, now I see the importance of providing Jacobian!
Hmmm, talking about Jacobian, if both DVODE and LSODA are given Jacobian, in your opinion which one is faster (if both are solving say, non-stiff problem)? Do you think as @mecej4 pointed, it will be useful to use the MA48 version of dvode?
Or, does LSODA have something like MA48 already? Thanks!
I was looking for a MA48 version of DVODE before but cannot find such a version. I also emailed S. Thompson but perhaps he no longer use his email anymore so no reply 
But since LSODA is very decent I guess I just mainly use LOSDA.
One more very naĆÆve question @Pap , hmmm, it seems there are some ode solver does not require jacobian, does those solver reliable?
Like RKC, RKC: An explicit solver for parabolic PDEs - ScienceDirect
which can be found as rkc.f in ode
or the modern ones like FLINT by @prince_mahajan which mainly based on DOP853,
GitHub - princemahajan/FLINT: Fortran Library for numerical INTegration of differential equations
For non-stiff problems, any decent adaptive step size algorithm should do, and I do not expect a significant performance difference between DVODE and LSODA. Howeverā¦ are you sure your problem is not stiff? In many cases it is easy to tell if a problem is stiff (just a glimpse in the differential equation may smell stiffness ahead). However, it is not so easy to tell if a problem is not stiff. You will be surprised how many problems seem innocent and stiffness-free while they arenāt. Experience in the topic definitely helps in this case.
Using an algorithm for stiff problems to solve a non-stiff problem is a bad idea, but if the algorithm is flexible enough, able to switch methods when necessary, then it is worth trying that one first. See what kind of a solution you get; plot the results; do you see any bottlenecks? Any small regions where many integration points are taken? (like the nearly vertical segment in the second figure above). If you donāt see such things, it is relatively safe to assume the problem is not stiff, and you can use DVODE or RKF45 or any other decent method for similar problems; they should perform equally well. And I say ārelatively safeā because even the term stiffness doesnāt have a unique, universally adopted definition.
Alternatively, use DVODE and see what results you get. If you see nearly horizontal segments with many integration points taken there (like the first figure above), then you are using the wrong method.
I must add something here, though. You mention āperformanceā and āspeedā very often, way more often than you mention āconvergenceā. Sure, it is always nice to have speed, and tons of books dedicate huge sections in algorithm improvements for better performance. But they do so after they assured convergence. if performance is essential (and in many cases is), I would prefer a safe algorithm and will look for ways to parallelize the code, instead of a fast but potentially unsafe algorithm.
Solving a āreal worldā numerical problem is like a gunfight in old wild West: Your first and foremost priority is to make sure you shoot your opponent then to do that fast. You want to shoot him fast, sure, but you want to make sure your bullet goes to the right place first. A slight delay to aim well is better than a fast āintuitiveā shoot. Listen to what an old gunslinger says 
: Fast, inattentive gunslingers wonāt see many summers.
2 Likes
Note that modern solvers that use dual number based differentiation (rather than finite difference) donāt require the user to pass in jacobians (and can sometimes be faster by more aggressively combining the code to calculate the primal and jacobian values).
1 Like