chem/KPP/kpp/kpp-2.1/int/oldies/ros2_cts_adj.f

   1       SUBROUTINE ros2_cts_adj(N,T,Tnext,Hmin,Hmax,Hstart,
   2      +                   y,Lambda,Fix,Rconst,AbsTol,RelTol,
   3      +                   Info)
   4       INCLUDE 'KPP_ROOT_params.h'
   5       INCLUDE 'KPP_ROOT_global.h'
   6
   7 C  INPUT ARGUMENTS:
   8 C     y = Vector of (NVAR) concentrations, contains the
   9 C         initial values on input
  10 C     [T, Tnext] = the integration interval
  11 C     Hmin, Hmax = lower and upper bounds for the selected step-size.
  12 C          Note that for Step = Hmin the current computed
  13 C          solution is unconditionally accepted by the error
  14 C          control mechanism.
  15 C     AbsTol, RelTol = (NVAR) dimensional vectors of
  16 C          componentwise absolute and relative tolerances.
  17 C     FUN = name of routine of derivatives. KPP syntax.
  18 C          See the header below.
  19 C     JAC_SP = name of routine that computes the Jacobian, in
  20 C          sparse format. KPP syntax. See the header below.
  21 C     Info(1) = 1  for  autonomous   system
  22 C             = 0  for nonautonomous system
  23 C     Info(2) = 1  for third order embedded formula
  24 C             = 0  for first order embedded formula
  25 C
  26 C   Note: Stage 3 used to build strongly A-stable order 3 formula for error control
  27 C            Embed3 = (Info(2).EQ.1)
  28 C         IF Embed3 = .true. THEN the third order embedded formula is used
  29 C                     .false. THEN a first order embedded  formula is used
  30 C
  31 C
  32 C  OUTPUT ARGUMENTS:
  33 C     y = the values of concentrations at TEND.
  34 C     T = equals TEND on output.
  35 C     Info(2) = # of FUN CALLs.
  36 C     Info(3) = # of JAC_SP CALLs.
  37 C     Info(4) = # of accepted steps.
  38 C     Info(5) = # of rejected steps.
  39
  40       INTEGER max_no_steps
  41       PARAMETER (max_no_steps = 200)
  42       KPP_REAL  Trajectory(NVAR,max_no_steps)
  43       KPP_REAL  StepSize(max_no_steps)
  44
  45       KPP_REAL  K1(NVAR), K2(NVAR), K3(NVAR)
  46       KPP_REAL  F1(NVAR), JAC(LU_NONZERO)
  47       KPP_REAL  DFDT(NVAR)(NRAD)
  48       KPP_REAL  Fix(NFIX), Rconst(NREACT)
  49       KPP_REAL  Hmin,Hmax,Hstart,ghinv,uround
  50       KPP_REAL  y(NVAR), Ynew(NVAR)
  51       KPP_REAL  AbsTol(NVAR), RelTol(NVAR)
  52       KPP_REAL  T, Tnext, H, Hold, Tplus
  53       KPP_REAL  ERR, factor, facmax
  54       KPP_REAL  Lambda(NVAR), K11(NVAR), JAC1(LU_NONZERO)
  55       INTEGER    n,nfcn,njac,Naccept,Nreject,i,j
  56       INTEGER    Info(5)
  57       LOGICAL    IsReject, Autonomous, Embed3
  58       EXTERNAL    FUN, JAC_SP
  59
  60       KPP_REAL gamma, m1, m2, alpha, beta1, beta2, delta, w, e
  61       KPP_REAL ginv
  62 c     Initialization of counters, etc.
  63       Autonomous = Info(1) .EQ. 1
  64       Embed3  = Info(2) .EQ. 1
  65       uround  = 1.d-15
  66       dround  = dsqrt(uround)
  67       H = DMAX1(Hstart,DMAX1(1.d-8, Hmin))
  68       Tplus = T
  69       IsReject = .false.
  70       Naccept  = 0
  71       Nreject  = 0
  72       Nfcn     = 0
  73       Njac     = 0
  74
  75 C     Method Parameters
  76       gamma = 1.d0 + 1.d0/sqrt(2.d0)
  77       a21   = - 1.d0/gamma
  78       m1 = -3.d0/(2.d0*gamma)
  79       m2 = -1.d0/(2.d0*gamma)
  80       c31 = -1.0D0/gamma**2*(1.0D0-7.0D0*gamma+9.0D0*gamma**2)
  81      &         /(-1.0D0+2.0D0*gamma)
  82       c32 = -1.0D0/gamma**2*(1.0D0-6.0D0*gamma+6.0D0*gamma**2)
  83      &          /(-1.0D0+2.0D0*gamma)/2
  84       gamma3 = 0.5D0 - 2*gamma
  85       d1 = ((-9.0D0*gamma+8.0D0*gamma**2+2.0D0)/gamma**2/
  86      &          (-1.0D0+2*gamma))/6.0D0
  87       d2 = ((-1.0D0+3.0D0*gam)/gamma**2/
  88      &          (-1.0D0+2.0D0*gamma))/6.0D0
  89       d3 = -1.0D0/(3.0D0*gamma)
  90
  91       Trajectory(1:NVAR,1) = Ynew(1)
  92
  93 C === Starting the time loop ===
  94  10    CONTINUE
  95        Tplus = T + H
  96        IF ( Tplus .gt. Tnext ) THEN
  97           H = Tnext - T
  98           Tplus = Tnext
  99        END IF
 100
 101        CALL Jac_SP( Y, Fix, Rconst, JAC )
 102
 103        Njac = Njac+1
 104        ghinv = -1.0d0/(gamma*H)
 105        DO 20 j=1,NVAR
 106          JAC(LU_DIAG_V(j)) = JAC(LU_DIAG_V(j)) + ghinv
 107  20    CONTINUE
 108        CALL KppDecomp (NVAR, JAC, ier)
 109
 110        IF (ier.ne.0) THEN
 111          IF ( H.gt.Hmin) THEN
 112             H = 5.0d-1*H
 113             GO TO 10
 114          else
 115             PRINT *,'IER <> 0, H=',H
 116             STOP
 117          END IF
 118        END IF
 119
 120        CALL Fun( Y, Fix, Rconst, F1 )
 121
 122 C ====== NONAUTONOMOUS CASE ===============
 123        IF (.not. Autonomous) THEN
 124          tau = dsign(dround*dmax1( 1.0d-6, dabs(T) ), T)
 125          CALL Fun( Y, Fix, Rconst, K2 )
 126          nfcn=nfcn+1
 127          DO 30 j = 1,NVAR
 128            DFDT(j) = ( K2(j)-F1(j) )/tau
 129  30      CONTINUE
 130        END IF ! .NOT.Autonomous
 131
 132 C ----- STAGE 1 -----
 133        DO 40 j = 1,NVAR
 134           K1(j) =  F1(j)
 135  40    CONTINUE
 136        IF (.NOT.Autonomous) THEN
 137          delta = gamma*H
 138          DO 45 j = 1,NVAR
 139            K1(j) = K1(j) + delta*DFDT(j)
 140  45      CONTINUE
 141        END IF ! .NOT.Autonomous
 142        CALL KppSolve (JAC, K1)
 143
 144 C ----- STAGE 2  -----
 145        DO 50 j = 1,NVAR
 146          Ynew(j) = y(j) + a21*K1(j)
 147  50    CONTINUE
 148        CALL Fun( Ynew, Fix, Rconst, F1 )
 149        nfcn=nfcn+1
 150        beta = 2.d0/(gamma*H)
 151        delta = -gamma*H
 152        DO 55 j = 1,NVAR
 153          K2(j) = F1(j) + beta*K1(j)
 154  55    CONTINUE
 155        IF (.NOT.Autonomous) THEN
 156          DO 56 j = 1,NVAR
 157            K2(j) = K2(j) + delta*DFDT(j)
 158  56      CONTINUE
 159        END IF ! .NOT.Autonomous
 160        CALL KppSolve (JAC, K2)
 161
 162 C ----- STAGE 3  -----
 163        IF (Embed3) THEN
 164        beta1 = -c31/H
 165        beta2 = -c32/H
 166        delta = gamma3*H
 167        DO 57 j = 1,NVAR
 168          K3(j) = F1(j) + beta1*K1(j) + beta2*K2(j)
 169  57    CONTINUE
 170        IF (.NOT.Autonomous) THEN
 171          DO 58 j = 1,NVAR
 172            K3(j) = K3(j) + delta*DFDT(j)
 173  58      CONTINUE
 174        END IF ! .NOT.Autonomous
 175        CALL KppSolve (JAC, K3)
 176        END IF ! Embed3
 177
 178
 179 C ---- The Solution ---
 180        DO 120 j = 1,NVAR
 181          Ynew(j) = y(j) + m1*K1(j) + m2*K2(j)
 182  120   CONTINUE
 183
 184
 185 C ====== Error estimation ========
 186
 187         ERR=0.d0
 188         DO 130 i=1,NVAR
 189            w = AbsTol(i) + RelTol(i)*DMAX1(DABS(y(i)),DABS(Ynew(i)))
 190            IF ( Embed3 ) THEN
 191              e = d1*K1(i) + d2*K2(i) + d3*K3(i)
 192            ELSE
 193              e = 1.d0/(2.d0*gamma)*(K1(i)+K2(i))
 194            END IF  ! Embed3
 195            ERR = ERR + ( e/w )**2
 196  130    CONTINUE
 197         ERR = DMAX1( uround, DSQRT( ERR/NVAR ) )
 198
 199 C ======= Choose the stepsize ===============================
 200
 201         IF ( Embed3 ) THEN
 202             elo = 3.0D0 ! estimator local order
 203         ELSE
 204             elo = 2.0D0
 205         END IF
 206         factor = DMAX1(2.0D-1,DMIN1(6.0D0,ERR**(1.0D0/elo)/.9D0))
 207         Hnew   = DMIN1(Hmax,DMAX1(Hmin, H/factor))
 208         Hold = H
 209
 210 C ======= Rejected/Accepted Step ============================
 211
 212         IF ( (ERR.gt.1).and.(H.gt.Hmin) ) THEN
 213           IsReject = .true.
 214           H = DMIN1(H/10,Hnew)
 215           Nreject  = Nreject+1
 216         ELSE
 217           DO 140 i=1,NVAR
 218              y(i)  = Ynew(i)
 219  140      CONTINUE
 220           T = Tplus
 221           IF (.NOT.IsReject) THEN
 222               H = Hnew   ! Do not increase stepsize IF previous step was rejected
 223           END IF
 224           IsReject = .false.
 225           Naccept = Naccept+1
 226           IF (Naccept+1>max_no_steps) THEN
 227             PRINT*,'Error in Adjoint Ros2: more steps than allowed'
 228             STOP
 229           END IF
 230           Trajectory(1:NVAR,Naccept+1) = Ynew(1:NVAR)
 231           StepSize(Naccept) = Hold
 232 !         CALL TRAJISTORE(y,hold)
 233         END IF
 234 C ======= END of the time loop ===============================
 235       IF ( T .lt. Tnext ) GO TO 10
 236
 237 C ======= Output Information =================================
 238       Info(2) = Nfcn
 239       Info(3) = Njac
 240       Info(4) = Naccept
 241       Info(5) = Nreject
 242
 243       ginv = 1.d0/gamma
 244 C -- The backwards loop for the CONTINUOUS ADJOINT
 245       DO istep = Naccept,1,-1
 246
 247         h = StepSize(istep)
 248         y(1:NVAR) = Trajectory(1:NVAR,istep+1)
 249         gHinv = -ginv/H
 250
 251         CALL Jac_SP(Y, Fix, Rconst, JAC)
 252         JAC1(1:LU_NONZERO)=JAC(1:LU_NONZERO)
 253         DO j=1,NVAR
 254           JAC(lu_diag_v(j)) = JAC(lu_diag_v(j)) + gHinv
 255         END DO
 256         CALL KppDecomp (NVAR,JAC,ier)
 257 ccc equivalent to function evaluation in forward integration
 258 ccc is J^T*Lambda in backward integration
 259         CALL JacTR_SP_Vec ( JAC1, Lambda, F1)
 260
 261 C ----- STAGE 1 (AUTONOMOUS) -----
 262         K11(1:NVAR) =  F1(1:NVAR)
 263         CALL KppSolveTR (JAC,K11,K1)
 264 C ----- STAGE 2 (AUTONOMOUS) -----
 265         y(1:NVAR) = Trajectory(1:NVAR,istep)
 266         CALL Jac_SP(Y, Fix, Rconst, JAC1)
 267         Ynew(1:NVAR) = Lambda(1:NVAR) - ginv*K1(1:NVAR)
 268         CALL JacTR_SP_Vec ( JAC1, Ynew, F1)
 269         beta = -2.d0*ghinv
 270         K11(1:NVAR) = F1(1:NVAR) + beta*K1(1:NVAR)
 271         CALL KppSolveTR (JAC,K11,K2)
 272 c ---- The solution
 273         Lambda(1:NVAR) = Lambda(1:NVAR)+m1*K1(1:NVAR)+m2*K2(1:NVAR)
 274
 275       END DO ! istep
 276
 277
 278       RETURN
 279       END
 280
 281
 282       SUBROUTINE FUNC_CHEM(N, T, Y, P)
 283       INCLUDE 'KPP_ROOT_params.h'
 284       INCLUDE 'KPP_ROOT_global.h'
 285       INTEGER N
 286       KPP_REAL   T, Told
 287       KPP_REAL   Y(NVAR), P(NVAR)
 288       Told = TIME
 289       TIME = T
 290       CALL Update_SUN()
 291       CALL Update_RCONST()
 292       CALL Fun( Y,  FIX, RCONST, P )
 293       TIME = Told
 294       RETURN
 295       END
 296
 297
 298       SUBROUTINE JAC_CHEM(N, T, Y, J)
 299       INCLUDE 'KPP_ROOT_params.h'
 300       INCLUDE 'KPP_ROOT_global.h'
 301       INTEGER N
 302       KPP_REAL   Told, T
 303       KPP_REAL   Y(NVAR), J(LU_NONZERO)
 304       Told = TIME
 305       TIME = T
 306       CALL Update_SUN()
 307       CALL Update_RCONST()
 308       CALL Jac_SP( Y,  FIX, RCONST, J )
 309       TIME = Told
 310       RETURN
 311       END
 312