Core integrator for dop853. parameters same as in dop853 with workspace added.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(dop853_class), | intent(inout) | :: | me | |||
| real(kind=wp), | intent(inout) | :: | x | |||
| real(kind=wp), | intent(inout), | dimension(:) | :: | y | ||
| real(kind=wp), | intent(in) | :: | xend | |||
| real(kind=wp), | intent(inout) | :: | hmax | |||
| real(kind=wp), | intent(inout) | :: | h | |||
| real(kind=wp), | intent(in), | dimension(:) | :: | rtol | ||
| real(kind=wp), | intent(in), | dimension(:) | :: | atol | ||
| integer, | intent(in) | :: | itol | |||
| integer, | intent(in) | :: | iprint | |||
| integer, | intent(in) | :: | iout | |||
| integer, | intent(out) | :: | idid | |||
| integer, | intent(in) | :: | nmax | |||
| integer, | intent(in) | :: | nstiff | |||
| real(kind=wp), | intent(in) | :: | safe | |||
| real(kind=wp), | intent(in) | :: | beta | |||
| real(kind=wp), | intent(in) | :: | fac1 | |||
| real(kind=wp), | intent(in) | :: | fac2 | |||
| integer, | intent(inout) | :: | nfcn | |||
| integer, | intent(inout) | :: | nstep | |||
| integer, | intent(inout) | :: | naccpt | |||
| integer, | intent(inout) | :: | nrejct |
subroutine dp86co(me,x,y,xend,hmax,h,rtol,atol,itol,iprint, &
iout,idid,nmax,nstiff,safe, &
beta,fac1,fac2, &
nfcn,nstep,naccpt,nrejct)
implicit none
class(dop853_class),intent(inout) :: me
real(wp),intent(inout) :: x
real(wp),dimension(:),intent(inout) :: y
real(wp),intent(in) :: xend
real(wp),intent(inout) :: hmax
real(wp),intent(inout) :: h
real(wp),dimension(:),intent(in) :: rtol
real(wp),dimension(:),intent(in) :: atol
integer,intent(in) :: itol
integer,intent(in) :: iprint
integer,intent(in) :: iout
integer,intent(out) :: idid
integer,intent(in) :: nmax
integer,intent(in) :: nstiff
real(wp),intent(in) :: safe
real(wp),intent(in) :: beta
real(wp),intent(in) :: fac1
real(wp),intent(in) :: fac2
integer,intent(inout) :: nfcn
integer,intent(inout) :: nstep
integer,intent(inout) :: naccpt
integer,intent(inout) :: nrejct
real(wp),dimension(me%n) :: y1,k1,k2,k3,k4,k5,k6,k7,k8,k9,k10
real(wp) :: atoli,bspl,deno,err,err2,erri,expo1,fac,fac11,&
facc1,facc2,facold,hlamb,hnew,posneg,rtoli,&
sk,stden,stnum,xout,xph,ydiff
integer :: i,iasti,iord,irtrn,j,nonsti,nrd
logical :: reject,last,event,abort
! initializations
nrd = me%nrdens
facold = 1.0e-4_wp
expo1 = 1.0_wp/8.0_wp - beta*0.2_wp
facc1 = 1.0_wp/fac1
facc2 = 1.0_wp/fac2
posneg = sign(1.0_wp,xend-x)
irtrn = 0 ! these were not initialized
nonsti = 0 ! in the original code
xout = 0.0_wp !
! initial preparations
atoli = atol(1)
rtoli = rtol(1)
last = .false.
hlamb = 0.0_wp
iasti = 0
call me%fcn(x,y,k1)
hmax = abs(hmax)
iord = 8
if ( h==0.0_wp ) then
h = me%hinit(x,y,posneg,k1,iord,hmax,atol,rtol,itol)
end if
nfcn = nfcn + 2
reject = .false.
me%xold = x
if ( iout/=0 ) then
irtrn = 1
me%hout = 1.0_wp
call me%solout(naccpt+1,me%xold,x,y,irtrn,xout)
abort = ( irtrn<0 )
else
abort = .false.
end if
if (.not. abort) then
do
! basic integration step
if ( nstep>nmax ) then
if ( iprint/=0 ) &
write (iprint,'(A,E18.4)') ' exit of dop853 at x=', x
if ( iprint/=0 ) &
write (iprint,*) ' more than nmax =' , nmax , 'steps are needed'
idid = -2
return
elseif ( 0.1_wp*abs(h)<=abs(x)*uround ) then
if ( iprint/=0 ) &
write (iprint,'(A,E18.4)') ' exit of dop853 at x=', x
if ( iprint/=0 ) &
write (iprint,*) ' step size too small, h=' , h
idid = -3
return
else
if ( (x+1.01_wp*h-xend)*posneg>0.0_wp ) then
h = xend - x
last = .true.
end if
nstep = nstep + 1
! the twelve stages
if ( irtrn>=2 ) call me%fcn(x,y,k1)
y1 = y + h*a21*k1
call me%fcn(x+c2*h,y1,k2)
y1 = y + h*(a31*k1+a32*k2)
call me%fcn(x+c3*h,y1,k3)
y1 = y + h*(a41*k1+a43*k3)
call me%fcn(x+c4*h,y1,k4)
y1 = y + h*(a51*k1+a53*k3+a54*k4)
call me%fcn(x+c5*h,y1,k5)
y1 = y + h*(a61*k1+a64*k4+a65*k5)
call me%fcn(x+c6*h,y1,k6)
y1 = y + h*(a71*k1+a74*k4+a75*k5+a76*k6)
call me%fcn(x+c7*h,y1,k7)
y1 = y + h*(a81*k1+a84*k4+a85*k5+a86*k6+a87*k7)
call me%fcn(x+c8*h,y1,k8)
y1 = y + h*(a91*k1+a94*k4+a95*k5+a96*k6+a97*k7+a98*k8)
call me%fcn(x+c9*h,y1,k9)
y1 = y + h*(a101*k1+a104*k4+a105*k5+a106*k6+a107*k7+&
a108*k8+a109*k9)
call me%fcn(x+c10*h,y1,k10)
y1 = y + h*(a111*k1+a114*k4+a115*k5+a116*k6+a117*k7+&
a118*k8+a119*k9+a1110*k10)
call me%fcn(x+c11*h,y1,k2)
xph = x + h
y1 = y + h*(a121*k1+a124*k4+a125*k5+a126*k6+a127*k7+&
a128*k8+a129*k9+a1210*k10+a1211*k2)
call me%fcn(xph,y1,k3)
nfcn = nfcn + 11
k4 = b1*k1+b6*k6+b7*k7+b8*k8+b9*k9+b10*k10+b11*k2+b12*k3
k5 = y + h*k4
! error estimation
err = 0.0_wp
err2 = 0.0_wp
if ( itol==0 ) then
do i = 1 , me%n
sk = atoli + rtoli*max(abs(y(i)),abs(k5(i)))
erri = k4(i) - bhh1*k1(i) - bhh2*k9(i) - bhh3*k3(i)
err2 = err2 + (erri/sk)**2
erri = er1*k1(i) + er6*k6(i) + er7*k7(i) + er8*k8(i) &
+ er9*k9(i) + er10*k10(i) + er11*k2(i) &
+ er12*k3(i)
err = err + (erri/sk)**2
end do
else
do i = 1 , me%n
sk = atol(i) + rtol(i)*max(abs(y(i)),abs(k5(i)))
erri = k4(i) - bhh1*k1(i) - bhh2*k9(i) - bhh3*k3(i)
err2 = err2 + (erri/sk)**2
erri = er1*k1(i) + er6*k6(i) + er7*k7(i) + er8*k8(i) &
+ er9*k9(i) + er10*k10(i) + er11*k2(i) &
+ er12*k3(i)
err = err + (erri/sk)**2
end do
end if
deno = err + 0.01_wp*err2
if ( deno<=0.0_wp ) deno = 1.0_wp
err = abs(h)*err*sqrt(1.0_wp/(me%n*deno))
! computation of hnew
fac11 = err**expo1
! lund-stabilization
fac = fac11/facold**beta
! we require fac1 <= hnew/h <= fac2
fac = max(facc2,min(facc1,fac/safe))
hnew = h/fac
if ( err<=1.0_wp ) then
! step is accepted
facold = max(err,1.0e-4_wp)
naccpt = naccpt + 1
call me%fcn(xph,k5,k4)
nfcn = nfcn + 1
! stiffness detection
if ( mod(naccpt,nstiff)==0 .or. iasti>0 ) then
stnum = 0.0_wp
stden = 0.0_wp
do i = 1 , me%n
stnum = stnum + (k4(i)-k3(i))**2
stden = stden + (k5(i)-y1(i))**2
end do
if ( stden>0.0_wp ) hlamb = abs(h)*sqrt(stnum/stden)
if ( hlamb>6.1_wp ) then
nonsti = 0
iasti = iasti + 1
if ( iasti==15 ) then
if ( iprint/=0 ) &
write (iprint,*) &
' the problem seems to become stiff at x = ', x
if ( iprint==0 ) then
idid = -4 ! fail exit
return
end if
end if
else
nonsti = nonsti + 1
if ( nonsti==6 ) iasti = 0
end if
end if
! final preparation for dense output
event = (iout==3) .and. (xout<=xph)
if ( iout==2 .or. event ) then
! save the first function evaluations
do j = 1 , nrd
i = me%icomp(j)
me%cont(j) = y(i)
ydiff = k5(i) - y(i)
me%cont(j+nrd) = ydiff
bspl = h*k1(i) - ydiff
me%cont(j+nrd*2) = bspl
me%cont(j+nrd*3) = ydiff - h*k4(i) - bspl
me%cont(j+nrd*4) = d41*k1(i)+d46*k6(i)+d47*k7(i)+&
d48*k8(i)+d49*k9(i)+d410*k10(i)+&
d411*k2(i)+d412*k3(i)
me%cont(j+nrd*5) = d51*k1(i)+d56*k6(i)+d57*k7(i)+&
d58*k8(i)+d59*k9(i)+d510*k10(i)+&
d511*k2(i)+d512*k3(i)
me%cont(j+nrd*6) = d61*k1(i)+d66*k6(i)+d67*k7(i)+&
d68*k8(i)+d69*k9(i)+d610*k10(i)+&
d611*k2(i)+d612*k3(i)
me%cont(j+nrd*7) = d71*k1(i)+d76*k6(i)+d77*k7(i)+&
d78*k8(i)+d79*k9(i)+d710*k10(i)+&
d711*k2(i)+d712*k3(i)
end do
! the next three function evaluations
y1 = y + h*(a141*k1+a147*k7+a148*k8+a149*k9+&
a1410*k10+a1411*k2+a1412*k3+a1413*k4)
call me%fcn(x+c14*h,y1,k10)
y1 = y + h*(a151*k1+a156*k6+a157*k7+a158*k8+&
a1511*k2+a1512*k3+a1513*k4+a1514*k10)
call me%fcn(x+c15*h,y1,k2)
y1 = y + h*(a161*k1+a166*k6+a167*k7+a168*k8+a169*k9+&
a1613*k4+a1614*k10+a1615*k2)
call me%fcn(x+c16*h,y1,k3)
nfcn = nfcn + 3
! final preparation
do j = 1 , nrd
i = me%icomp(j)
me%cont(j+nrd*4) = h*(me%cont(j+nrd*4)+d413*k4(i)+&
d414*k10(i)+d415*k2(i)+d416*k3(i))
me%cont(j+nrd*5) = h*(me%cont(j+nrd*5)+d513*k4(i)+&
d514*k10(i)+d515*k2(i)+d516*k3(i))
me%cont(j+nrd*6) = h*(me%cont(j+nrd*6)+d613*k4(i)+&
d614*k10(i)+d615*k2(i)+d616*k3(i))
me%cont(j+nrd*7) = h*(me%cont(j+nrd*7)+d713*k4(i)+&
d714*k10(i)+d715*k2(i)+d716*k3(i))
end do
me%hout = h
end if
k1 = k4
y = k5
me%xold = x
x = xph
if ( iout==1 .or. iout==2 .or. event ) then
call me%solout(naccpt+1,me%xold,x,y,irtrn,xout)
if ( irtrn<0 ) exit !abort
end if
! normal exit
if ( last ) then
h = hnew
idid = 1
return
end if
if ( abs(hnew)>hmax ) hnew = posneg*hmax
if ( reject ) hnew = posneg*min(abs(hnew),abs(h))
reject = .false.
else
! step is rejected
hnew = h/min(facc1,fac11/safe)
reject = .true.
if ( naccpt>=1 ) nrejct = nrejct + 1
last = .false.
end if
h = hnew
end if
end do
end if
if ( iprint/=0 ) write (iprint,'(A,E18.4)') ' exit of dop853 at x=', x
idid = 2
end subroutine dp86co