ddstp.f

Go to the documentation of this file.

C Work performed under the auspices of the U.S. Department of Energy
C by Lawrence Livermore National Laboratory under contract number
C W-7405-Eng-48.
C
      SUBROUTINE ddstp(X,Y,YPRIME,NEQ,RES,JAC,PSOL,H,WT,VT,
     *  jstart,idid,rpar,ipar,phi,savr,delta,e,wm,iwm,
     *  alpha,beta,gamma,psi,sigma,cj,cjold,hold,s,hmin,uround,
     *  epli,sqrtn,rsqrtn,epcon,iphase,jcalc,jflg,k,kold,ns,nonneg,
     *  ntype,nls)
C
C***BEGIN PROLOGUE  DDSTP
C***REFER TO  DDASPK
C***DATE WRITTEN   890101   (YYMMDD)
C***REVISION DATE  900926   (YYMMDD)
C***REVISION DATE  940909   (YYMMDD) (Reset PSI(1), PHI(*,2) at 690)
C
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C     DDSTP solves a system of differential/algebraic equations of
C     the form G(X,Y,YPRIME) = 0, for one step (normally from X to X+H).
C
C     The methods used are modified divided difference, fixed leading
C     coefficient forms of backward differentiation formulas.
C     The code adjusts the stepsize and order to control the local error
C     per step.
C
C
C     The parameters represent
C     X  --        Independent variable.
C     Y  --        Solution vector at X.
C     YPRIME --    Derivative of solution vector
C                  after successful step.
C     NEQ --       Number of equations to be integrated.
C     RES --       External user-supplied subroutine
C                  to evaluate the residual.  See RES description
C                  in DDASPK prologue.
C     JAC --       External user-supplied routine to update
C                  Jacobian or preconditioner information in the
C                  nonlinear solver.  See JAC description in DDASPK
C                  prologue.
C     PSOL --      External user-supplied routine to solve
C                  a linear system using preconditioning.
C                  (This is optional).  See PSOL in DDASPK prologue.
C     H --         Appropriate step size for next step.
C                  Normally determined by the code.
C     WT --        Vector of weights for error criterion used in Newton test.
C     VT --        Masked vector of weights used in error test.
C     JSTART --    Integer variable set 0 for
C                  first step, 1 otherwise.
C     IDID --      Completion code returned from the nonlinear solver.
C                  See IDID description in DDASPK prologue.
C     RPAR,IPAR -- Real and integer parameter arrays that
C                  are used for communication between the
C                  calling program and external user routines.
C                  They are not altered by DNSK
C     PHI --       Array of divided differences used by
C                  DDSTP. The length is NEQ*(K+1), where
C                  K is the maximum order.
C     SAVR --      Work vector for DDSTP of length NEQ.
C     DELTA,E --   Work vectors for DDSTP of length NEQ.
C     WM,IWM --    Real and integer arrays storing
C                  information required by the linear solver.
C
C     The other parameters are information
C     which is needed internally by DDSTP to
C     continue from step to step.
C
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C   NLS, DDWNRM, DDATRP
C
C***END PROLOGUE  DDSTP
C
C
      IMPLICIT DOUBLE PRECISION(a-h,o-z)
      dimension y(*),yprime(*),wt(*),vt(*)
      dimension phi(neq,*),savr(*),delta(*),e(*)
      dimension wm(*),iwm(*)
      dimension psi(*),alpha(*),beta(*),gamma(*),sigma(*)
      dimension rpar(*),ipar(*)
      EXTERNAL  res, jac, psol, nls
C
      parameter(lmxord=3)
      parameter(lnst=11, letf=14, lcfn=15)
C
C
C-----------------------------------------------------------------------
C     BLOCK 1.
C     Initialize.  On the first call, set
C     the order to 1 and initialize
C     other variables.
C-----------------------------------------------------------------------
C
C     Initializations for all calls
C
      xold=x
      ncf=0
      nef=0
      IF(jstart .NE. 0) go to 120
C
C     If this is the first step, perform
C     other initializations
C
      k=1
      kold=0
      hold=0.0d0
      psi(1)=h
      cj = 1.d0/h
      iphase = 0
      ns=0
120   CONTINUE
C
C
C
C
C
C-----------------------------------------------------------------------
C     BLOCK 2
C     Compute coefficients of formulas for
C     this step.
C-----------------------------------------------------------------------
200   CONTINUE
      kp1=k+1
      kp2=k+2
      km1=k-1
      IF(h.NE.hold.OR.k .NE. kold) ns = 0
      ns=min0(ns+1,kold+2)
      nsp1=ns+1
      IF(kp1 .LT. ns)go to 230
C
      beta(1)=1.0d0
      alpha(1)=1.0d0
      temp1=h
      gamma(1)=0.0d0
      sigma(1)=1.0d0
      DO 210 i=2,kp1
         temp2=psi(i-1)
         psi(i-1)=temp1
         beta(i)=beta(i-1)*psi(i-1)/temp2
         temp1=temp2+h
         alpha(i)=h/temp1
         sigma(i)=(i-1)*sigma(i-1)*alpha(i)
         gamma(i)=gamma(i-1)+alpha(i-1)/h
210      CONTINUE
      psi(kp1)=temp1
230   CONTINUE
C
C     Compute ALPHAS, ALPHA0
C
      alphas = 0.0d0
      alpha0 = 0.0d0
      DO 240 i = 1,k
        alphas = alphas - 1.0d0/i
        alpha0 = alpha0 - alpha(i)
240     CONTINUE
C
C     Compute leading coefficient CJ
C
      cjlast = cj
      cj = -alphas/h
C
C     Compute variable stepsize error coefficient CK
C
      ck = abs(alpha(kp1) + alphas - alpha0)
      ck = max(ck,alpha(kp1))
C
C     Change PHI to PHI STAR
C
      IF(kp1 .LT. nsp1) go to 280
      DO 270 j=nsp1,kp1
         DO 260 i=1,neq
260         phi(i,j)=beta(j)*phi(i,j)
270      CONTINUE
280   CONTINUE
C
C     Update time
C
      x=x+h
C
C     Initialize IDID to 1
C
      idid = 1
C
C
C
C
C
C-----------------------------------------------------------------------
C     BLOCK 3
C     Call the nonlinear system solver to obtain the solution and
C     derivative.
C-----------------------------------------------------------------------
C
      CALL nls(x,y,yprime,neq,
     *   res,jac,psol,h,wt,jstart,idid,rpar,ipar,phi,gamma,
     *   savr,delta,e,wm,iwm,cj,cjold,cjlast,s,
     *   uround,epli,sqrtn,rsqrtn,epcon,jcalc,jflg,kp1,
     *   nonneg,ntype,iernls)
C
      IF(iernls .NE. 0)go to 600
C
C
C
C
C
C-----------------------------------------------------------------------
C     BLOCK 4
C     Estimate the errors at orders K,K-1,K-2
C     as if constant stepsize was used. Estimate
C     the local error at order K and test
C     whether the current step is successful.
C-----------------------------------------------------------------------
C
C     Estimate errors at orders K,K-1,K-2
C
      enorm = ddwnrm(neq,e,vt,rpar,ipar)
      erk = sigma(k+1)*enorm
      terk = (k+1)*erk
      est = erk
      knew=k
      IF(k .EQ. 1)go to 430
      DO 405 i = 1,neq
405     delta(i) = phi(i,kp1) + e(i)
      erkm1=sigma(k)*ddwnrm(neq,delta,vt,rpar,ipar)
      terkm1 = k*erkm1
      IF(k .GT. 2)go to 410
      IF(terkm1 .LE. 0.5*terk)go to 420
      go to 430
410   CONTINUE
      DO 415 i = 1,neq
415     delta(i) = phi(i,k) + delta(i)
      erkm2=sigma(k-1)*ddwnrm(neq,delta,vt,rpar,ipar)
      terkm2 = (k-1)*erkm2
      IF(max(terkm1,terkm2).GT.terk)go to 430
C
C     Lower the order
C
420   CONTINUE
      knew=k-1
      est = erkm1
C
C
C     Calculate the local error for the current step
C     to see if the step was successful
C
430   CONTINUE
      err = ck * enorm
      IF(err .GT. 1.0d0)go to 600
C
C
C
C
C
C-----------------------------------------------------------------------
C     BLOCK 5
C     The step is successful. Determine
C     the best order and stepsize for
C     the next step. Update the differences
C     for the next step.
C-----------------------------------------------------------------------
      idid=1
      iwm(lnst)=iwm(lnst)+1
      kdiff=k-kold
      kold=k
      hold=h
C
C
C     Estimate the error at order K+1 unless
C        already decided to lower order, or
C        already using maximum order, or
C        stepsize not constant, or
C        order raised in previous step
C
      IF(knew.EQ.km1.OR.k.EQ.iwm(lmxord))iphase=1
      IF(iphase .EQ. 0)go to 545
      IF(knew.EQ.km1)go to 540
      IF(k.EQ.iwm(lmxord)) go to 550
      IF(kp1.GE.ns.OR.kdiff.EQ.1)go to 550
      DO 510 i=1,neq
510      delta(i)=e(i)-phi(i,kp2)
      erkp1 = (1.0d0/(k+2))*ddwnrm(neq,delta,vt,rpar,ipar)
      terkp1 = (k+2)*erkp1
      IF(k.GT.1)go to 520
      IF(terkp1.GE.0.5d0*terk)go to 550
      go to 530
520   IF(terkm1.LE.min(terk,terkp1))go to 540
      IF(terkp1.GE.terk.OR.k.EQ.iwm(lmxord))go to 550
C
C     Raise order
C
530   k=kp1
      est = erkp1
      go to 550
C
C     Lower order
C
540   k=km1
      est = erkm1
      go to 550
C
C     If IPHASE = 0, increase order by one and multiply stepsize by
C     factor two
C
545   k = kp1
      hnew = h*2.0d0
      h = hnew
      go to 575
C
C
C     Determine the appropriate stepsize for
C     the next step.
C
550   hnew=h
      temp2=k+1
      r=(2.0d0*est+0.0001d0)**(-1.0d0/temp2)
      IF(r .LT. 2.0d0) go to 555
      hnew = 2.0d0*h
      go to 560
555   IF(r .GT. 1.0d0) go to 560
      r = max(0.5d0,min(0.9d0,r))
      hnew = h*r
560   h=hnew
C
C
C     Update differences for next step
C
575   CONTINUE
      IF(kold.EQ.iwm(lmxord))go to 585
      DO 580 i=1,neq
580      phi(i,kp2)=e(i)
585   CONTINUE
      DO 590 i=1,neq
590      phi(i,kp1)=phi(i,kp1)+e(i)
      DO 595 j1=2,kp1
         j=kp1-j1+1
         DO 595 i=1,neq
595      phi(i,j)=phi(i,j)+phi(i,j+1)
      jstart = 1
      RETURN
C
C
C
C
C
C-----------------------------------------------------------------------
C     BLOCK 6
C     The step is unsuccessful. Restore X,PSI,PHI
C     Determine appropriate stepsize for
C     continuing the integration, or exit with
C     an error flag if there have been many
C     failures.
C-----------------------------------------------------------------------
600   iphase = 1
C
C     Restore X,PHI,PSI
C
      x=xold
      IF(kp1.LT.nsp1)go to 630
      DO 620 j=nsp1,kp1
         temp1=1.0d0/beta(j)
         DO 610 i=1,neq
610         phi(i,j)=temp1*phi(i,j)
620      CONTINUE
630   CONTINUE
      DO 640 i=2,kp1
640      psi(i-1)=psi(i)-h
C
C
C     Test whether failure is due to nonlinear solver
C     or error test
C
      IF(iernls .EQ. 0)go to 660
      iwm(lcfn)=iwm(lcfn)+1
C
C
C     The nonlinear solver failed to converge.
C     Determine the cause of the failure and take appropriate action.
C     If IERNLS .LT. 0, then return.  Otherwise, reduce the stepsize
C     and try again, unless too many failures have occurred.
C
      IF (iernls .LT. 0) go to 675
      ncf = ncf + 1
      r = 0.25d0
      h = h*r
      IF (ncf .LT. 10 .AND. abs(h) .GE. hmin) go to 690
      IF (idid .EQ. 1) idid = -7
      IF (nef .GE. 3) idid = -9
      go to 675
C
C
C     The nonlinear solver converged, and the cause
C     of the failure was the error estimate
C     exceeding the tolerance.
C
660   nef=nef+1
      iwm(letf)=iwm(letf)+1
      IF (nef .GT. 1) go to 665
C
C     On first error test failure, keep current order or lower
C     order by one.  Compute new stepsize based on differences
C     of the solution.
C
      k = knew
      temp2 = k + 1
      r = 0.90d0*(2.0d0*est+0.0001d0)**(-1.0d0/temp2)
      r = max(0.25d0,min(0.9d0,r))
      h = h*r
      IF (abs(h) .GE. hmin) go to 690
      idid = -6
      go to 675
C
C     On second error test failure, use the current order or
C     decrease order by one.  Reduce the stepsize by a factor of
C     one quarter.
C
665   IF (nef .GT. 2) go to 670
      k = knew
      r = 0.25d0
      h = r*h
      IF (abs(h) .GE. hmin) go to 690
      idid = -6
      go to 675
C
C     On third and subsequent error test failures, set the order to
C     one, and reduce the stepsize by a factor of one quarter.
C
670   k = 1
      r = 0.25d0
      h = r*h
      IF (abs(h) .GE. hmin) go to 690
      idid = -6
      go to 675
C
C
C
C
C     For all crashes, restore Y to its last value,
C     interpolate to find YPRIME at last X, and return.
C
C     Before returning, verify that the user has not set
C     IDID to a nonnegative value.  If the user has set IDID
C     to a nonnegative value, then reset IDID to be -7, indicating
C     a failure in the nonlinear system solver.
C
675   CONTINUE
      CALL ddatrp(x,x,y,yprime,neq,k,phi,psi)
      jstart = 1
      IF (idid .GE. 0) idid = -7
      RETURN
C
C
C     Go back and try this step again.
C     If this is the first step, reset PSI(1) and rescale PHI(*,2).
C
690   IF (kold .EQ. 0) THEN
        psi(1) = h
        DO 695 i = 1,neq
695       phi(i,2) = r*phi(i,2)
        ENDIF
      go to 200
C
C------END OF SUBROUTINE DDSTP------------------------------------------
      END