dvstep performs one step of the integration of an initial value problem for a system of ordinary differential equations.
dvstep calls subroutine dvnlsd for the solution of the nonlinear system arising in the time step. thus it is independent of the problem jacobian structure and the type of nonlinear system solution method. dvstep returns a completion flag kflag (in common). a return with kflag = -1 or -2 means either abs(h) = hmin or 10 consecutive failures occurred. on a return with kflag negative, the values of tn and the yh array are as of the beginning of the last step, and h is the last step size attempted.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| class(dvode_t), | intent(inout) | :: | me | |||
| real(kind=wp), | intent(inout) | :: | y(*) |
an array of length n used for the dependent variable vector. |
||
| real(kind=wp), | intent(inout) | :: | yh(ldyh,*) |
an ldyh by lmax array containing the dependent variables and their approximate scaled derivatives, where lmax = maxord + 1. yh(i,j+1) contains the approximate j-th derivative of y(i), scaled by h**j/factorial(j) (j = 0,1,...,nq). on entry for the first step, the first two columns of yh must be set from the initial values. |
||
| integer, | intent(in) | :: | ldyh |
a constant integer >= n, the first dimension of yh. n is the number of odes in the system. |
||
| real(kind=wp), | intent(inout) | :: | yh1(*) |
a one-dimensional array occupying the same space as yh. |
||
| real(kind=wp), | intent(in) | :: | ewt(*) |
an array of length n containing multiplicative weights for local error measurements. local errors in y(i) are compared to 1.0/ewt(i) in various error tests. |
||
| real(kind=wp) | :: | savf(*) |
an array of working storage, of length n. also used for input of yh(*,maxord+2) when jstart = -1 and maxord < the current order nq. |
|||
| real(kind=wp) | :: | vsav(*) |
a work array of length n passed to subroutine dvnlsd. |
|||
| real(kind=wp) | :: | acor(*) |
a work array of length n, used for the accumulated corrections. on a successful return, acor(i) contains the estimated one-step local error in y(i). |
|||
| real(kind=wp) | :: | wm(*) |
real work array associated with matrix operations in dvnlsd. |
|||
| integer | :: | iwm(*) |
integer work array associated with matrix operations in dvnlsd. |
subroutine dvstep(me,y,yh,ldyh,yh1,ewt,savf,vsav,acor,wm,iwm) class(dvode_t),intent(inout) :: me integer,intent(in) :: ldyh !! a constant integer >= n, the first dimension of yh. !! n is the number of odes in the system. real(wp),intent(inout) :: y(*) !! an array of length n used for the dependent variable vector. real(wp),intent(inout) :: yh(ldyh,*) !! an ldyh by lmax array containing the dependent variables !! and their approximate scaled derivatives, where !! lmax = maxord + 1. yh(i,j+1) contains the approximate !! j-th derivative of y(i), scaled by h**j/factorial(j) !! (j = 0,1,...,nq). on entry for the first step, the first !! two columns of yh must be set from the initial values. real(wp),intent(inout) :: yh1(*) !! a one-dimensional array occupying the same space as yh. real(wp),intent(in) :: ewt(*) !! an array of length n containing multiplicative weights !! for local error measurements. local errors in y(i) are !! compared to 1.0/ewt(i) in various error tests. real(wp) :: savf(*) !! an array of working storage, of length n. !! also used for input of yh(*,maxord+2) when jstart = -1 !! and maxord < the current order nq. real(wp) :: vsav(*) !! a work array of length n passed to subroutine [[dvnlsd]]. real(wp) :: acor(*) !! a work array of length n, used for the accumulated !! corrections. on a successful return, acor(i) contains !! the estimated one-step local error in y(i). real(wp) :: wm(*) !! real work array associated with matrix !! operations in [[dvnlsd]]. integer :: iwm(*) !! integer work array associated with matrix !! operations in [[dvnlsd]]. real(wp) :: cnquot , ddn , dsm , dup , etaqp1 , & flotl , r , told integer :: i , i1 , i2 , iback , j , jb , ncf , nflag integer,parameter :: kfc = -3 integer,parameter :: kfh = -7 integer,parameter :: mxncf = 10 real(wp),parameter :: one = 1.0_wp real(wp),parameter :: addon = 1.0e-6_wp real(wp),parameter :: bias1 = 6.0_wp real(wp),parameter :: bias2 = 6.0_wp real(wp),parameter :: bias3 = 10.0_wp real(wp),parameter :: etacf = 0.25_wp real(wp),parameter :: etamin = 0.1_wp real(wp),parameter :: etamxf = 0.2_wp real(wp),parameter :: etamx1 = 10000.0_wp real(wp),parameter :: etamx2 = 10.0_wp real(wp),parameter :: etamx3 = 10.0_wp real(wp),parameter :: onepsm = 1.00001_wp real(wp),parameter :: thresh = 1.5_wp me%dat%kflag = 0 told = me%dat%tn ncf = 0 me%dat%jcur = 0 nflag = 0 if ( me%dat%jstart<=0 ) then if ( me%dat%jstart==-1 ) then !----------------------------------------------------------------------- ! the following block handles preliminaries needed when jstart = -1. ! if n was reduced, zero out part of yh to avoid undefined references. ! if maxord was reduced to a value less than the tentative order newq, ! then nq is set to maxord, and a new h ratio eta is chosen. ! otherwise, we take the same preliminary actions as for jstart > 0. ! in any case, nqwait is reset to l = nq + 1 to prevent further ! changes in order for that many steps. ! the new h ratio eta is limited by the input h if kuth = 1, ! by hmin if kuth = 0, and by hmxi in any case. ! finally, the history array yh is rescaled. !----------------------------------------------------------------------- me%dat%lmax = me%dat%maxord + 1 if ( me%dat%n/=ldyh ) then i1 = 1 + (me%dat%newq+1)*ldyh i2 = (me%dat%maxord+1)*ldyh if ( i1<=i2 ) then do i = i1 , i2 yh1(i) = zero enddo endif endif if ( me%dat%newq>me%dat%maxord ) then flotl = real(me%dat%lmax,wp) if ( me%dat%maxord<me%dat%nq-1 ) then ddn = me%dvnorm(me%dat%n,savf(1:me%dat%n),ewt(1:me%dat%n))/me%dat%tq(1) me%dat%eta = one/((bias1*ddn)**(one/flotl)+addon) endif if ( me%dat%maxord==me%dat%nq .and. me%dat%newq==me%dat%nq+1 ) me%dat%eta = me%etaq if ( me%dat%maxord==me%dat%nq-1 .and. me%dat%newq==me%dat%nq+1 ) then me%dat%eta = me%etaqm1 call me%dvjust(yh,ldyh,-1) endif if ( me%dat%maxord==me%dat%nq-1 .and. me%dat%newq==me%dat%nq ) then ddn = me%dvnorm(me%dat%n,savf(1:me%dat%n),ewt(1:me%dat%n))/me%dat%tq(1) me%dat%eta = one/((bias1*ddn)**(one/flotl)+addon) call me%dvjust(yh,ldyh,-1) endif me%dat%eta = min(me%dat%eta,one) me%dat%nq = me%dat%maxord me%dat%l = me%dat%lmax endif if ( me%dat%kuth==1 ) me%dat%eta = min(me%dat%eta,abs(me%dat%h/me%dat%hscal)) if ( me%dat%kuth==0 ) me%dat%eta = max(me%dat%eta,me%dat%hmin/abs(me%dat%hscal)) me%dat%eta = me%dat%eta/max(one,abs(me%dat%hscal)*me%dat%hmxi*me%dat%eta) me%dat%newh = 1 me%dat%nqwait = me%dat%l if ( me%dat%newq>me%dat%maxord ) goto 300 else !----------------------------------------------------------------------- ! on the first call, the order is set to 1, and other variables are ! initialized. etamax is the maximum ratio by which h can be increased ! in a single step. it is normally 10, but is larger during the ! first step to compensate for the small initial h. if a failure ! occurs (in corrector convergence or error test), etamax is set to 1 ! for the next increase. !----------------------------------------------------------------------- me%dat%lmax = me%dat%maxord + 1 me%dat%nq = 1 me%dat%l = 2 me%dat%nqnyh = me%dat%nq*ldyh me%dat%tau(1) = me%dat%h me%dat%prl1 = one me%dat%rc = zero me%dat%etamax = etamx1 me%dat%nqwait = 2 me%dat%hscal = me%dat%h goto 400 end if else if ( me%dat%kuth==1 ) then !----------------------------------------------------------------------- ! take preliminary actions on a normal continuation step (jstart>0). ! if the driver changed h, then eta must be reset and newh set to 1. ! if a change of order was dictated on the previous step, then ! it is done here and appropriate adjustments in the history are made. ! on an order decrease, the history array is adjusted by dvjust. ! on an order increase, the history array is augmented by a column. ! on a change of step size h, the history array yh is rescaled. !----------------------------------------------------------------------- me%dat%eta = min(me%dat%eta,me%dat%h/me%dat%hscal) me%dat%newh = 1 endif if ( me%dat%newh==0 ) goto 400 if ( me%dat%newq<me%dat%nq ) then call me%dvjust(yh,ldyh,-1) me%dat%nq = me%dat%newq me%dat%l = me%dat%nq + 1 me%dat%nqwait = me%dat%l else if ( me%dat%newq>me%dat%nq ) then call me%dvjust(yh,ldyh,1) me%dat%nq = me%dat%newq me%dat%l = me%dat%nq + 1 me%dat%nqwait = me%dat%l endif ! rescale the history array for a change in h by a factor of eta. ------ 300 r = one do j = 2 , me%dat%l r = r*me%dat%eta call dscal(me%dat%n,r,yh(1,j),1) enddo me%dat%h = me%dat%hscal*me%dat%eta me%dat%hscal = me%dat%h me%dat%rc = me%dat%rc*me%dat%eta me%dat%nqnyh = me%dat%nq*ldyh !----------------------------------------------------------------------- ! this section computes the predicted values by effectively ! multiplying the yh array by the pascal triangle matrix. ! dvset is called to calculate all integration coefficients. ! rc is the ratio of new to old values of the coefficient h/el(2)=h/l1. !----------------------------------------------------------------------- 400 me%dat%tn = me%dat%tn + me%dat%h i1 = me%dat%nqnyh + 1 do jb = 1 , me%dat%nq i1 = i1 - ldyh do i = i1 , me%dat%nqnyh yh1(i) = yh1(i) + yh1(i+ldyh) enddo enddo call me%dvset() me%dat%rl1 = one/me%dat%el(2) me%dat%rc = me%dat%rc*(me%dat%rl1/me%dat%prl1) me%dat%prl1 = me%dat%rl1 ! call the nonlinear system solver. ------------------------------------ call me%dvnlsd(y,yh,ldyh,vsav,savf,ewt,acor,iwm,wm,nflag) if ( nflag==0 ) then !----------------------------------------------------------------------- ! the corrector has converged (nflag = 0). the local error test is ! made and control passes to statement 500 if it fails. !----------------------------------------------------------------------- dsm = me%dat%acnrm/me%dat%tq(2) if ( dsm>one ) then !----------------------------------------------------------------------- ! the error test failed. kflag keeps track of multiple failures. ! restore tn and the yh array to their previous values, and prepare ! to try the step again. compute the optimum step size for the ! same order. after repeated failures, h is forced to decrease ! more rapidly. !----------------------------------------------------------------------- me%dat%kflag = me%dat%kflag - 1 me%dat%netf = me%dat%netf + 1 nflag = -2 me%dat%tn = told i1 = me%dat%nqnyh + 1 do jb = 1 , me%dat%nq i1 = i1 - ldyh do i = i1 , me%dat%nqnyh yh1(i) = yh1(i) - yh1(i+ldyh) enddo enddo if ( abs(me%dat%h)<=me%dat%hmin*onepsm ) then !----------------------------------------------------------------------- ! all returns are made through this section. ! on a successful return, etamax is reset and acor is scaled. !----------------------------------------------------------------------- me%dat%kflag = -1 me%dat%jstart = 1 return else me%dat%etamax = one if ( me%dat%kflag>kfc ) then ! compute ratio of new h to current h at the current order. ------------ flotl = real(me%dat%l,wp) me%dat%eta = one/((bias2*dsm)**(one/flotl)+addon) me%dat%eta = max(me%dat%eta,me%dat%hmin/abs(me%dat%h),etamin) if ( (me%dat%kflag<=-2) .and. (me%dat%eta>etamxf) ) me%dat%eta = etamxf goto 300 !----------------------------------------------------------------------- ! control reaches this section if 3 or more consecutive failures ! have occurred. it is assumed that the elements of the yh array ! have accumulated errors of the wrong order. the order is reduced ! by one, if possible. then h is reduced by a factor of 0.1 and ! the step is retried. after a total of 7 consecutive failures, ! an exit is taken with kflag = -1. !----------------------------------------------------------------------- else if ( me%dat%kflag==kfh ) then me%dat%kflag = -1 me%dat%jstart = 1 return else if ( me%dat%nq==1 ) then me%dat%eta = max(etamin,me%dat%hmin/abs(me%dat%h)) me%dat%h = me%dat%h*me%dat%eta me%dat%hscal = me%dat%h me%dat%tau(1) = me%dat%h call me%f(me%dat%n,me%dat%tn,y(1:me%dat%n),savf(1:me%dat%n)) me%dat%nfe = me%dat%nfe + 1 do i = 1 , me%dat%n yh(i,2) = me%dat%h*savf(i) enddo me%dat%nqwait = 10 goto 400 else me%dat%eta = max(etamin,me%dat%hmin/abs(me%dat%h)) call me%dvjust(yh,ldyh,-1) me%dat%l = me%dat%nq me%dat%nq = me%dat%nq - 1 me%dat%nqwait = me%dat%l goto 300 endif endif else !----------------------------------------------------------------------- ! after a successful step, update the yh and tau arrays and decrement ! nqwait. if nqwait is then 1 and nq < maxord, then acor is saved ! for use in a possible order increase on the next step. ! if etamax = 1 (a failure occurred this step), keep nqwait >= 2. !----------------------------------------------------------------------- me%dat%kflag = 0 me%dat%nst = me%dat%nst + 1 me%dat%hu = me%dat%h me%dat%nqu = me%dat%nq do iback = 1 , me%dat%nq i = me%dat%l - iback me%dat%tau(i+1) = me%dat%tau(i) enddo me%dat%tau(1) = me%dat%h do j = 1 , me%dat%l call daxpy(me%dat%n,me%dat%el(j),acor,1,yh(1,j),1) enddo me%dat%nqwait = me%dat%nqwait - 1 if ( (me%dat%l/=me%dat%lmax) .and. (me%dat%nqwait==1) ) then call dcopy(me%dat%n,acor,1,yh(1,me%dat%lmax),1) me%dat%conp = me%dat%tq(5) endif if ( me%dat%etamax/=one ) then !----------------------------------------------------------------------- ! if nqwait = 0, an increase or decrease in order by one is considered. ! factors etaq, etaqm1, etaqp1 are computed by which h could ! be multiplied at order q, q-1, or q+1, respectively. ! the largest of these is determined, and the new order and ! step size set accordingly. ! a change of h or nq is made only if h increases by at least a ! factor of thresh. if an order change is considered and rejected, ! then nqwait is set to 2 (reconsider it after 2 steps). !----------------------------------------------------------------------- ! compute ratio of new h to current h at the current order. ------------ flotl = real(me%dat%l,wp) me%etaq = one/((bias2*dsm)**(one/flotl)+addon) if ( me%dat%nqwait==0 ) then me%dat%nqwait = 2 me%etaqm1 = zero if ( me%dat%nq/=1 ) then ! compute ratio of new h to current h at the current order less one. --- ddn = me%dvnorm(me%dat%n,yh(1,me%dat%l),ewt(1:me%dat%n))/me%dat%tq(1) me%etaqm1 = one/((bias1*ddn)**(one/(flotl-one))+addon) endif etaqp1 = zero if ( me%dat%l/=me%dat%lmax ) then ! compute ratio of new h to current h at current order plus one. ------- cnquot = (me%dat%tq(5)/me%dat%conp)*(me%dat%h/me%dat%tau(2))**me%dat%l do i = 1 , me%dat%n savf(i) = acor(i) - cnquot*yh(i,me%dat%lmax) enddo dup = me%dvnorm(me%dat%n,savf(1:me%dat%n),ewt(1:me%dat%n))/me%dat%tq(3) etaqp1 = one/((bias3*dup)**(one/(flotl+one))+addon) endif if ( me%etaq<etaqp1 ) then if ( etaqp1<=me%etaqm1 ) goto 420 me%dat%eta = etaqp1 me%dat%newq = me%dat%nq + 1 call dcopy(me%dat%n,acor,1,yh(1,me%dat%lmax),1) goto 450 else if ( me%etaq<me%etaqm1 ) then goto 420 endif endif me%dat%eta = me%etaq me%dat%newq = me%dat%nq goto 450 else if ( me%dat%nqwait<2 ) me%dat%nqwait = 2 me%dat%newq = me%dat%nq me%dat%newh = 0 me%dat%eta = one me%dat%hnew = me%dat%h goto 500 endif 420 me%dat%eta = me%etaqm1 me%dat%newq = me%dat%nq - 1 endif ! test tentative new h against thresh, etamax, and hmxi, then exit. ---- 450 if ( me%dat%eta<thresh .or. me%dat%etamax==one ) then me%dat%newq = me%dat%nq me%dat%newh = 0 me%dat%eta = one me%dat%hnew = me%dat%h else me%dat%eta = min(me%dat%eta,me%dat%etamax) me%dat%eta = me%dat%eta/max(one,abs(me%dat%h)*me%dat%hmxi*me%dat%eta) me%dat%newh = 1 me%dat%hnew = me%dat%h*me%dat%eta endif else !----------------------------------------------------------------------- ! the vnls routine failed to achieve convergence (nflag /= 0). ! the yh array is retracted to its values before prediction. ! the step size h is reduced and the step is retried, if possible. ! otherwise, an error exit is taken. !----------------------------------------------------------------------- ncf = ncf + 1 me%dat%ncfn = me%dat%ncfn + 1 me%dat%etamax = one me%dat%tn = told i1 = me%dat%nqnyh + 1 do jb = 1 , me%dat%nq i1 = i1 - ldyh do i = i1 , me%dat%nqnyh yh1(i) = yh1(i) - yh1(i+ldyh) enddo enddo if ( nflag<-1 ) then if ( nflag==-2 ) me%dat%kflag = -3 if ( nflag==-3 ) me%dat%kflag = -4 me%dat%jstart = 1 return else if ( abs(me%dat%h)<=me%dat%hmin*onepsm ) then me%dat%kflag = -2 me%dat%jstart = 1 return else if ( ncf==mxncf ) then me%dat%kflag = -2 me%dat%jstart = 1 return else me%dat%eta = etacf me%dat%eta = max(me%dat%eta,me%dat%hmin/abs(me%dat%h)) nflag = -1 goto 300 endif endif 500 me%dat%etamax = etamx3 if ( me%dat%nst<=10 ) me%dat%etamax = etamx2 r = one/me%dat%tq(2) call dscal(me%dat%n,r,acor,1) me%dat%jstart = 1 end subroutine dvstep