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}
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.
