ddasic.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 ddasic (X, Y, YPRIME, NEQ, ICOPT, ID, RES, JAC, PSOL,
     *   H, WT, NIC, IDID, RPAR, IPAR, PHI, SAVR, DELTA, E, YIC, YPIC,
     *   PWK, WM, IWM, HMIN, UROUND, EPLI, SQRTN, RSQRTN, EPCONI,
     *   STPTOL, JFLG, ICNFLG, ICNSTR, NLSIC)
C
C***BEGIN PROLOGUE  DDASIC
C***REFER TO  DDASPK
C***DATE WRITTEN   940628   (YYMMDD)
C***REVISION DATE  941206   (YYMMDD)
C***REVISION DATE  950714   (YYMMDD)
C
C-----------------------------------------------------------------------
C***DESCRIPTION
C
C     DDASIC is a driver routine to compute consistent initial values
C     for Y and YPRIME.  There are two different options:
C     Denoting the differential variables in Y by Y_d, and
C     the algebraic variables by Y_a, the problem solved is either:
C     1.  Given Y_d, calculate Y_a and Y_d', or
C     2.  Given Y', calculate Y.
C     In either case, initial values for the given components
C     are input, and initial guesses for the unknown components
C     must also be provided as input.
C
C     The external routine NLSIC solves the resulting nonlinear system.
C
C     The parameters represent
C
C     X  --        Independent variable.
C     Y  --        Solution vector at X.
C     YPRIME --    Derivative of solution vector.
C     NEQ --       Number of equations to be integrated.
C     ICOPT     -- Flag indicating initial condition option chosen.
C                    ICOPT = 1 for option 1 above.
C                    ICOPT = 2 for option 2.
C     ID        -- Array of dimension NEQ, which must be initialized
C                  if option 1 is chosen.
C                    ID(i) = +1 if Y_i is a differential variable,
C                    ID(i) = -1 if Y_i is an algebraic variable.
C     RES --       External user-supplied subroutine to evaluate the
C                  residual.  See RES description in DDASPK prologue.
C     JAC --       External user-supplied routine to update Jacobian
C                  or preconditioner information in the nonlinear solver
C                  (optional).  See JAC description in DDASPK prologue.
C     PSOL --      External user-supplied routine to solve
C                  a linear system using preconditioning.
C                  See PSOL in DDASPK prologue.
C     H --         Scaling factor in iteration matrix.  DDASIC may
C                  reduce H to achieve convergence.
C     WT --        Vector of weights for error criterion.
C     NIC --       Input number of initial condition calculation call
C                  (= 1 or 2).
C     IDID --      Completion code.  See IDID 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 --       Work space for DDASIC of length at least 2*NEQ.
C     SAVR --      Work vector for DDASIC of length NEQ.
C     DELTA --     Work vector for DDASIC of length NEQ.
C     E --         Work vector for DDASIC of length NEQ.
C     YIC,YPIC --  Work vectors for DDASIC, each of length NEQ.
C     PWK --       Work vector for DDASIC of length NEQ.
C     WM,IWM --    Real and integer arrays storing
C                  information required by the linear solver.
C     EPCONI --    Test constant for Newton iteration convergence.
C     ICNFLG --    Flag showing whether constraints on Y are to apply.
C     ICNSTR --    Integer array of length NEQ with constraint types.
C
C     The other parameters are for use internally by DDASIC.
C
C-----------------------------------------------------------------------
C***ROUTINES CALLED
C   DCOPY, NLSIC
C
C***END PROLOGUE  DDASIC
C
C
      IMPLICIT DOUBLE PRECISION(a-h,o-z)
      dimension y(*),yprime(*),id(*),wt(*),phi(neq,*)
      dimension savr(*),delta(*),e(*),yic(*),ypic(*),pwk(*)
      dimension wm(*),iwm(*), rpar(*),ipar(*), icnstr(*)
      EXTERNAL res, jac, psol, nlsic
C
      parameter(lcfn=15)
      parameter(lmxnh=34)
C
C The following parameters are data-loaded here:
C     RHCUT  = factor by which H is reduced on retry of Newton solve.
C     RATEMX = maximum convergence rate for which Newton iteration
C              is considered converging.
C
      SAVE rhcut, ratemx
      DATA rhcut/0.1d0/, ratemx/0.8d0/
C
C
C-----------------------------------------------------------------------
C     BLOCK 1.
C     Initializations.
C     JSKIP is a flag set to 1 when NIC = 2 and NH = 1, to signal that
C     the initial call to the JAC routine is to be skipped then.
C     Save Y and YPRIME in PHI.  Initialize IDID, NH, and CJ.
C-----------------------------------------------------------------------
C
      mxnh = iwm(lmxnh)
      idid = 1
      nh = 1
      jskip = 0
      IF (nic .EQ. 2) jskip = 1
      CALL dcopy (neq, y, 1, phi(1,1), 1)
      CALL dcopy (neq, yprime, 1, phi(1,2), 1)
C
      IF (icopt .EQ. 2) THEN
        cj = 0.0d0
      ELSE
        cj = 1.0d0/h
      ENDIF
C
C-----------------------------------------------------------------------
C     BLOCK 2
C     Call the nonlinear system solver to obtain
C     consistent initial values for Y and YPRIME.
C-----------------------------------------------------------------------
C
 200  CONTINUE
      CALL nlsic(x,y,yprime,neq,icopt,id,res,jac,psol,h,wt,jskip,
     *   rpar,ipar,savr,delta,e,yic,ypic,pwk,wm,iwm,cj,uround,
     *   epli,sqrtn,rsqrtn,epconi,ratemx,stptol,jflg,icnflg,icnstr,
     *   iernls)
C
      IF (iernls .EQ. 0) RETURN
C
C-----------------------------------------------------------------------
C     BLOCK 3
C     The nonlinear solver was unsuccessful.  Increment NCFN.
C     Return with IDID = -12 if either
C       IERNLS = -1: error is considered unrecoverable,
C       ICOPT = 2: we are doing initialization problem type 2, or
C       NH = MXNH: the maximum number of H values has been tried.
C     Otherwise (problem 1 with IERNLS .GE. 1), reduce H and try again.
C     If IERNLS > 1, restore Y and YPRIME to their original values.
C-----------------------------------------------------------------------
C
      iwm(lcfn) = iwm(lcfn) + 1
      jskip = 0
C
      IF (iernls .EQ. -1) GO TO 350
      IF (icopt .EQ. 2) GO TO 350
      IF (nh .EQ. mxnh) GO TO 350
C
      nh = nh + 1
      h = h*rhcut
      cj = 1.0d0/h
C
      IF (iernls .EQ. 1) GO TO 200
C
      CALL dcopy (neq, phi(1,1), 1, y, 1)
      CALL dcopy (neq, phi(1,2), 1, yprime, 1)
      GO TO 200
C
 350  idid = -12
      RETURN
C
C------END OF SUBROUTINE DDASIC-----------------------------------------
      END